Erupting Mount St. Helens

May 27, 2014 Examples No comments

The example in rsf/tutorials/sthelens reproduces the analysis by Evan Bianco of the digital elevation data from Mount St. Helens before and after its catastrophic eruption.


Madagascar users are encouraged to try improving the results.

Light Bartlein color palette

May 15, 2014 Examples No comments

The orange-white-purple diverging color palette was suggested in the article

Light, A. and P.J. Bartlein (2004) The end of the rainbow? Color schemes for improved data graphics. EOS Transactions of the American Geophysical Union 85(40):385

Matteo Niccoli recommends it for seismic data as a replacement for the familiar red-white-blue pallete (color=g in Madagascar). Now the Light-Bartlein palette is available to Madagascar plotting programs, such as sfgrey, as color=lb. See the following example from rsf/rsf/sfgrey:

More information:

Program of the month: sfhelicon

May 13, 2014 Programs No comments

sfhelicon performs multidimensional convolution and inverse convolution (recursive filtering) using the helix transform.

The theory behind helical convolution is explained by Jon Claerbout in the paper

Claerbout, J., 1998, Multidimensional recursive filters via a helix: Geophysics, 63, 1532-1541

and in the corresponding chapter in his book Image Estimation by Example. The following example from gee/hlx/helicon illustrates inverse convolution on a helix.

sfhelicon works in N dimensions. The filter for is supplied by filt= parameter, which points to a real-valued RSF file with filter coefficents. The filter lags on a helix can be stored in a separate integer-value RSF file specified wih lag= parameter (which can be optionally contained inside the file with filter coefficients). The dimensions used for specifying the filter lags are not necessarily the same as the dimensions of the input data and can be specified with n= parameter (which can be optionally contained inside the file with filter lags).

The choice between forward and inverse convolution is controlled by div= parameter. The adjoint flag is supplied by adj= parameter. The following figure illustrates the act of the adjoint convolution and inverse convolution on a helix.

For an example of least-squares inversion with sfhelicon using the conjugate-gradient algorithm, see the documentation for sfconjgrad.

10 previous programs of the month:

madagascar-1.6 released

May 9, 2014 Celebration No comments

The 1.6 stable release features fifteen new reproducible papers and multiple other enhancements including the addition of the seismic migration gallery.

According to the SourceForge statistics, the previous 1.5 stable distribution has been downloaded more than 5,000 times. The record number of downloads in September 2013 is probably due to the fact that Madagascar is being used for teaching at different universities. The top country (with 40% of all downloads) was China, followed by the US, Brazil, Mexico, and Australia.

According to Ohloh.net, the year before the 1.6 release was the period of a high development activity, with 41 contributors (up 41% compared to the previous year) making 1,762 commits to the repository. Ohloh.net says that Madagascar “has a well established, mature codebase maintained by a very large development team with stable year-over-year commits” and estimated 201 man-years of effort.

IWAVE update

April 28, 2014 Documentation No comments

The IWAVE package included in Madagascar went through a significant redesign. The current version is explained in the paper IWAVE Structure and Basic Use Cases by Bill Symes.

The IWAVE control structure facilitates construction of wave simulators with flexible specification of input and output. This document describes synthesis of seismograms and wavefield movies from initial data and from single and multiple sources (right-hand sides), and linearized (“Born”) and linearized adjoint (reverse time migration) modeling. The choice of physical model and simulation method – constant density acoustics with Dirichlet boundary conditions and (2,2k) finite difference schemes – is the simplest possible, but the framework accommodates any regularly gridded stencil-based discretization of arbitrary wave physics in the same way.

GPU implementation of RTM

April 26, 2014 Documentation No comments

A new paper is added to the collection of reproducible documents:
RTM using effective boundary saving: A staggered grid GPU implementation

GPU has become a booming technology in reverse time migration (RTM) to perform the intensive computation. Compared with saving forward modeled wavefield on the disk, RTM via wavefield reconstruction using saved boundaries on device is a more efficient method because computation is much faster than CPU-GPU data transfer. In this paper, we introduce the effective boundary saving strategy in backward reconstruction for RTM. The minimum storage requirement for regular and staggered grid finite difference is determined for perfect reconstruction of the source wavefield. Particularly, we implement RTM using GPU programming, combining staggered finite difference scheme with convolutional perfectly matched layer (CPML) boundary condition. We demonstrate the validity of the proposed approach and CUDA codes with numerical example and imaging of benchmark models.

This paper is the first contribution from Xi’an Jiaotong University in China.

Random noise attenuation by EMD predictive filtering

April 23, 2014 Documentation No comments

A new paper is added to the collection of reproducible documents:
Random noise attenuation by f-x empirical mode decomposition predictive filtering

Random noise attenuation always plays an important role in seismic data processing. One of the most widely used methods for suppressing random noise is f-x predictive filtering. When the subsurface structure becomes complex, this method suffers from higher prediction errors owing to the large number of different dip components that need to be predicted. In this paper, we propose a novel denoising method termed f-x empirical mode decomposition predictive filtering (EMDPF). This new scheme solves the problem that makes f-x empirical mode decomposition (EMD) ineffective with complex seismic data. Also, by making the prediction more precise, the new scheme removes the limitation of conventional f-x predictive filtering when dealing with multi-dip seismic profiles. In this new method, we first apply EMD to each frequency slice in the f-x domain and obtain several intrinsic mode functions (IMF). Then an auto-regressive (AR) model is applied to the sum of the first few IMFs, which contain the high-dip-angle components, in order to predict the useful steeper events. Finally, the predicted events are added to the sum of the remaining IMFs. This process improves the prediction precision by utilizing an EMD based dip filter to reduce the dip components before f-x predictive filtering. Both synthetic and real data sets demonstrate the performance of our proposed method in preserving more useful energy.

Pseudo-pure wave modes in anisotropic media

April 20, 2014 Documentation No comments

A new paper is added to the collection of reproducible documents:
Simulating propagation of separated wave modes in general anisotropic media, Part I: qP-wave propagators

Wave propagation in an anisotropic medium is inherently described by elastic wave equations, with P- and S-wave modes intrinsically coupled. We present an approach to simulate propagation of separated wave modes for forward modeling, migration, waveform inversion and other applications in general anisotropic media. The proposed approach consists of two cascaded computational steps. First, we simulate equivalent elastic anisotropic wavefields with a minimal second-order coupled system (that we call here a pseudo-pure-mode wave equation), which describes propagation of all wave modes with a partial wave mode separation. Such a system for qP-wave is derived from the inverse Fourier transform of the Christoffel equation after a similarity transformation, which aims to project the original vector displacement wavefields onto isotropic references of the polarization directions of qP-waves. It accurately describes the kinematics of all wave modes and enhances qP-waves when the pseudo-pure-mode wavefield components are summed. The second step is a filtering to further project the pseudo-pure-mode wavefields onto the polarization directions of qP-waves so that residual qS-wave energy is removed and scalar qP-wave fields are accurately separated at each time step during wavefield extrapolation. As demonstrated in the numerical examples, pseudo-pure-mode wave equation plus correction of projection deviation provides a robust and flexible tool for simulating propagation of separated wave modes in transversely isotropic and orthorhombic media. The synthetic example of Hess VTI model shows that the pseudo-pure-mode qP-wave equation can be used in prestack reverse-time migration (RTM) applications.

This paper is the first contribution from Tongji University in Shanghai, China.

Program of the month: sfcostaper

April 2, 2014 Programs No comments

sfcostaper applies cosine tapering to the edges of the input data.

Cosine tapering amounts to multiplying the edges by
$$\sin^2\left(\frac{\pi\,x}{2}\right) = \frac{1+\cos{(\pi\,x)}}{2}$$,
where $x$ is the relative distance from the start of the taper.

sfcostaper works in N dimensions, the size of the taper in #th dimension in samples is specified by nw#=. In the following example from gee/ajt/igrad1 the images are tapered at both edges by specifying nw1=50 nw2=50.

10 previous programs of the month:

Optimal aperture in Kirchhoff migration

March 25, 2014 Documentation No comments

A new paper is added to the collection of reproducible documents:
Selecting an optimal aperture in Kirchhoff migration using dip-angle images

We present a method for selecting a migration aperture in Kirchhoff migration. We first split migrated data into constant-dip-angle partial images. Then, in every partial image, we estimate the consistency between each event and the constant dip of the analyzed section. We filter out events whose slope is far from the corresponding dip. Stacking of the filtered partial images corresponds to migration having an optimal aperture. Synthetic and real data examples demonstrate that the proposed approach to migration-aperture optimization is able to reduce migration noise while preserving diffraction energy, which characterizes small geological objects and brings additional resolution to the image.