Program of the month: sfvscan

May 4, 2013 Programs No comments

sfvscan implements seismic velocity analysis by scanning stacking velocities. This transformation is also known as the velocity transform or the hyperbolic Radon transform.

The following example from bei/vela/vscan shows an example for transforming a CMP (common midpoint) gather to a velocity (actually slowness) scan.

By default, sfvscan uses velocity as the horizontal axis. To change it to slowness, use slowness=y. It is also possible to use velocity or slowness squared by specifying squared=y. The range of velocities or slownesses to scan is given by v0=, dv=, nv=. In addition, it is possible to scan heterogeneity parameters in the shifted-hyperbola approximation using smax= and ns=.

The offset range in the input file is specified similarly to how it is done in sfnmo: The offset can be either regular (specified as the second axis in the input file) or irregular (specified in offset= file). By default, half-offset is used. To use full offset, specify half=n. Additionally, it is possible to specify a mask (with mask= file containing 1s and 0s) for skipping certain input traces during the scan.

The NMO correction involved in velocity scanning is associated with the phenomenon of “NMO stretch”, a non-linear stretching of time at large offsets. The maximum allowed relative stretch is controlled by str= parameter. The part of the data that is stretched more than the allowed stretch gets muted. The width of the mute zone is controlled by mute= parameter.

By default, sfvscan outputs simple stack. To output semblance instead, use semblance=y or type=semblance. It is also possible to output other kinds of semblance attributes: differential semblance (type=w), AB semblance (type=a), or weighted semblance (type=w). See

Li, J., & Symes, W. W. (2007). Interval velocity estimation via NMO-based differential semblance. Geophysics, 72(6), U75-U88.

Sarkar, D., Castagna, J. P., & Lamb, W. J. (2001). AVO and velocity analysis. Geophysics, 66(4), 1284-1293.

Sarkar, D., Baumel, R. T., & Larner, K. L. (2002). Velocity analysis in the presence of amplitude variation. Geophysics, 67(5), 1664-1672.

Fomel, S. (2009). Velocity analysis using AB semblance. Geophysical Prospecting, 57(3), 311-321.

Luo, S., & Hale, D. (2012). Velocity analysis using weighted semblance. Geophysics, 77(2), U15-U22.

The following example from jsg/avo/avo2 shows a comparison between velocity scans computed using the conventional semblance and AB semblance.

For robustness, semblance values are averaged in time. The length of the averaging window in samples is given by nb= (the default is 2 time samples). By default, the input samples are stacked along the hyperbola with an asymptotically pseudounitary weight equal to the absolute value of offset times velocity. To apply a uniform weight, use weight=n. For the justification of the pseudounitary weighting, see

Claerbout, J. F., 1995, Basic Earth Imaging: Stanford Exploration Project.

S. Fomel, 2003, Asymptotic pseudounitary stacking operators: Geophysics, v. 68, 1032-1042.

For the operation inverse or adjoint to the velocity scan, use sfveltran, sfvelmod or sfvelinv.

For a faster implementation of the velocity scan (using the buttefly algorithm), use sfradon2.

From the output of the semblance scan computed with sfvscan, the velocity trend can be picked automatically with sfpick.

10 previous programs of the month

Events of the year

May 2, 2013 Celebration No comments

A Madagascar School will take place on August 15, 2013, in Melbourne, Australia, during the ASEG convention and is organized by Jeff Shragge from the University of Western Australia.

In addition, the First Madagascar Working Workshop will take place on July 25-27, 2013, in Austin, Texas. The declared objectives of the workshop are

  1. To expand Madagascar’s collection of reproducible papers. New reproducible papers may include papers written using Madagascar programs as well as papers written using other open-source software tools.
  2. To create and expand a seismic migration gallery. Migration gallery is a matrix where rows are different benchmark datasets and columns are different seismic migration algorithms.

If you are interested in participating in this event, please apply for a registration by following the link at https://www.ahay.org/wiki/Austin_2013. The application deadline is July 1.

Reasons not to share your code

April 22, 2013 Links No comments

In Top Ten Reasons To Not Share Your Code (and why you should anyway) published by SIAM News, Randy LeVeque, Professor of Applied Mathematics from the University of Washington, elegantly destroys common excuses computational scientists and applied mathematicians come up with when they refuse to share their software codes.

Today, most mathematicians find the idea of publishing a theorem without its proof laughable, even though many great mathematicians of the past apparently found it quite natural. Mathematics has since matured in healthy ways, and it seems inevitable that computational mathematics will follow a similar path, no matter how inconvenient it may seem. I sense growing concern among young people in particular about the way we’ve been doing things and the difficulty of understanding or building on earlier work […] We can all help our field mature by making the effort to share the code that supports our research.

As if to illustrate LeVeque’s point, major news outlets report on a story about a reproducibility error (Excel bug) discovered in the famous politically-influential paper by economists Reinhart and Rogoff:

The validity of the Reinhart-Rogoff assertion “once debt exceeds 90 percent of GDP, economic growth drops off sharply” continues to be debated by economists, but it is clear now that their data were flawed and that the error would have been discovered much easier if the publication had followed the reproducible-research discipline.

Zooming in Vplot

April 17, 2013 Programs No comments

sfzoom is a simple Tkinter script that allows for interactive zooming when viewing Vplot files generated by sfgrey. The zoomed images are created by windowing the original RSF file and therefore have the highest possible resolution.

Joe Dellinger has a more comprehensive plan for adding interactivity to Vplot graphics.

How can I read and write RSF files in MATLAB?

April 13, 2013 FAQ No comments

  • The most straightforward way is to install the MATLAB interface to Madagascar. When installing Madagascar, run 
    ./configure API=matlab

    The configure script will try to find and test matlab and mex executibles on your system. If they are not in your PATH, you can specify them with 

    ./configure API=matlab MATLAB=/path/to/matlab MEX=/path/to/mex

    Install Madagascar as usual, set MATLAB path to $RSFROOT/lib, and you will able to read and write RSF files from MATLAB using rsf_read, rsf_write, and other functions from the Madagascar interface.

  • Alternatively, you can try reading binary data using MATLAB functions, as in the following example 

    
% get in=, n1=, and n2= parameters from file.rsf
    
    [stat,in] = unix(‘sfget in parform=n < file.rsf’)
    
in = strtrim(in)
    
[stat,n1] = unix(‘sfget n1 parform=n < file.rsf’)
    
n1 = str2num(n1)
    
[stat,n2] = unix(‘sfget n2 parform=n < file.rsf’)
    
n2 = str2num(n2)
    
% read binary data
    
    fid = fopen(in,‘rb’)
    
data = fread(fid,n1*n2,‘float32’);
    
% reshape to 2-D matrix
    
    data = reshape(data,n1,n2);
  • An even better alternative is to abandone MATLAB in favor of free software, such as GNU Octave, Python with NumPy, Sage, etc. A Python interface to Madagascar is installed by default.

Program of the month: sfnmo

April 8, 2013 Programs No comments

sfnmo implements normal moveout (NMO) correction, one of the most fundamental operations in seismic reflection data processing.

The following example from jsg/avo/avo shows synthetic data before and after NMO correction.

NMO transforms prestack seismic gathers by stretching each trace according to the equation
$$t = \sqrt{t_0^2 + \frac{x^2}{v^2(t_0)}}\;,$$
where $t$ is time before the correction, $t_0$ is time after the correction, $x$ is offset, and $v(t_0)$ is velocity specified in velocity= file. The NMO velocity can be picked (for example, with sfpick) from velocity scans (for example, produced with sfvscan).

If s= file is provided, sfnmo uses Malovichko’s shifted-hyperbola approximation
$$t \approx \left(1-\frac{1}{s(t_0)}\right) + \frac{1}{s(t_0)}\,\sqrt{t_0^2 + s(t_0)\,\frac{x^2}{v^2(t_0)}}$$
See

  • Malovichko, A. A., 1978, A new representation of the traveltime curve of reflected waves in horizontally layered media: Applied Geophysics (in Russian), 91, 47-53
  • Sword, C. H., 1987, A Soviet look at datum shift, in SEP-51: Stanford Exploration Project, 313-316.
  • de Bazelaire, E., 1988, Normal moveout revisited – Inhomogeneous media and curved interfaces: Geophysics, 53, 143-157.
  • Castle, R. J., 1994, Theory of normal moveout: Geophysics, 59, 983-999.

If a= file is provided, sfnmo uses Taner’s velocity-acceleration approximation
$$t = \sqrt{t_0^2 + \frac{x^2}{v^2(t_0)+a(t_0)\,x^2}}$$
See

  • Taner, M. T., S. Treitel, and M. Al-Chalabi, 2005, A new travel time estimation method for horizontal strata: 75th Ann. Internat. Mtg, Soc. of Expl. Geophys., 2273-2276.
  • Taner, M. T., S. Treitel, M. Al-Chalabi, and S. Fomel, 2007, An offset dependent NMO velocity model: EAGE 69th Conference and Exhibition, EAGE, P036.

The offset $x$ can be either regular (specified as the second axis in the input file) or irregular (specified in offset= file). By default, half-offset is used $h=x/2$. To use full offset, specify half=n.

Additionally, it is possible to specify a mask (with mask= file containing 1s and 0s) for skipping certain traces during the correction. See the following example from rsf/scons/rsf

The NMO correction is associated with the phenomenon of “NMO stretch”, a non-linear stretching of time at large offsets. The maximum allowed relative stretch is controlled by str= parameter. The part of the data that is stretched more than the allowed stretch gets muted. The width of the mute zone is controlled by mute= parameter.

NMO with constant velocity can be accomplished with sfnmostretch. For 3-D azimuthally-anisotropic moveout correction, try sfnmo3. For NMO in the tau-p domain, try sftaupmo. For an inverse NMO operation, try sfnimo.

10 previous programs of the month

Automatic traveltime picking

April 2, 2013 Documentation No comments

A new paper is added to the collection of reproducible documents:
Automatic traveltime picking using the instantaneous traveltime

IPython

March 30, 2013 Celebration No comments

On March 23, Fernando Perez, the creator of IPython, received the Award for the Advancement of Free Software. The award was presented by Richard Stallman, the president of the Free Software Foundation, at a ceremony in Cambridge, Massachusetts.

IPython plays well with Madagascar through the use of the Python interface. A particularly convenient way of using IPython’s interactive capabilities is the web-based notebook interface.

We can expect new enhancements to IPython’s notebook interface from Fernando Perez and his collaborators thanks to a generous grant from the Sloan foundation.

US government gets serious about reproducible research

March 23, 2013 Links No comments

The debate on open science and reproducible research has reached Washington, DC.
On February 22, John Holdren (Assistant to the President for Science and Technology and Director of the White House Office of Science and Technology Policy) issued a Memorandum on Increasing Access to the Results of Federally Funded Scientific Research to the heads of all federal agencies that sponsor research and development projects. The memo states

Access to digital data sets resulting from federally funded research allows companies to focus resources and efforts on understanding and exploiting discoveries. For example, open weather data underpins the forecasting industry, and making genome sequences publicly available has spawned many biotechnology innovations. In addition, wider availability of peer reviewed publications and scientific data in digital formats will create innovative economic markets for services related to curation, preservation, analysis, and visualization. Policies that mobilize these publications and data for re-use through preservation and broader public access also maximize the impact and accountability of the Federal research investment. These policies will accelerate scientific breakthroughs and innovation, promote entrepreneurship, and enhance economic growth and job creation.

The memo obliges every agency to come up with a strategy for making both scientific publications and digital scientific data resulting from Federally funded research publicly available.

On March 5, the Subcommittee on Research of the US House Committee on Science, Space, and Technology’s held a hearing on the issue of access to data from federally funded published research. In his opening statement, Dan Lipinski, a democratic U.S. Representative from Illinois, said:

..the more data are open, the faster we will validate new theories and overturn old ones, and the more efficiently we will transform new discoveries into innovations that will create jobs and make us healthier and more prosperous. The movement toward open data is not primarily about scientific integrity, it’s mostly about speeding up the process of scientific discovery and innovation.

Victoria Stodden, an Assistant Professor of Statistics at Columbia University and a famous advocate for reproducible research, testified:

Making research data and software conveniently available also has valuable corollary effects beyond validating the original associated published results. Other researchers can use them for new research, linking datasets and augmenting results in other areas, or applying the software and methods to new research applications. These powerful benefits will accelerate scientific discovery. Benefits can also accrue to private industry. Again, data and software availability permit business to apply these methods to their own research problems, link with their own datasets, and accelerate innovation and economic growth.

Color palettes

March 19, 2013 Examples No comments

In addition to the standard color palettes, it is now possible to use a custom palette by specifying it in a simple comma-separated text file (CSV). The following example from rsf/rsf/sfgrey uses the perceptional rainbow palette created by Matteo Niccoli. See Matteo’s blog posts for an explanation.