Entries by Sergey Fomel

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:

• sfnmo
• sfpow
• sfpwd
• sfricker1
• sfhalfint
• sfbandpass
• sfkirmod
• sfiwarp
• sfpick
• sffft3

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 http://www.ahay.org/wiki/Austin_2013. The application deadline is July 1.

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 :

• Reinhart, Rogoff... and Herndon: The student who caught out the profs (BBS News)
• How much unemployment did Reinhart and Rogoff's arithmetic mistake cause? (Guardian)
• Is the evidence for austerity based on an Excel spreadsheet error? (Washington Post)

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.

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.

• 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.

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+\genfrac{}{}{0.1ex}{}{x^2}{v^2\left(t_0\right)}}$ where $t$ is time before the correction, $t_0$ is time after the correction, $x$ is offset, and $v\left(t_0\right)$ 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 (1-\genfrac{}{}{0.1ex}{}{1}{s\left(t_0\right)})t_0+\genfrac{}{}{0.1ex}{}{1}{s\left(t_0\right)}\sqrt{t_0^2+\genfrac{}{}{0.1ex}{}{s\left(t_0\right)x^2}{v^2\left(t_0\right)}}$. 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\approx \sqrt{t_0^2+\genfrac{}{}{0.1ex}{}{x^2}{v^2\left(t_0\right)+a\left(t_0\right)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:

• sfpow
• sfpwd
• sfricker1
• sfhalfint
• sfbandpass
• sfkirmod
• sfiwarp
• sfpick
• sffft3
• sfdip

A new paper is added to the collection of reproducible documents:

Automatic traveltime picking using the instantaneous traveltime

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.

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.

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.

sfpow multiplies the input data by a gain function of the form $$G\left(x_1,x_2,\dots ,x_n\right)=x_1^{p_1}{x}_2^{p_2}\cdots x_n^{p_n}$$ The powers $$p_1,p_2,\dots ,p_n$$ are given by pow1=, pow2=, etc. parameters.

For backward compatibility, sftpow tpow= is an alias for sfpow pow1=.

The following example from geo391/hw1/tpow shows an application of sfpow pow1=2 to a shot gather. The gain of pow1=2 for seismic data was advocated by Jon Claerbout.



To estimate an appropriate power gain from the data, you can try sffpow.

10 previous programs of the month:

• sfpwd
• sfricker1
• sfhalfint
• sfbandpass
• sfkirmod
• sfiwarp
• sfpick
• sffft3
• sfdip
• sfderiv

The official birthday of the Madagascar project is June 11, 2006, when the project was first publicly announced. However, the initial work on the project (called RSF at the time) started sometime in early 2003, 10 years ago.

The age of a car is usually measured in miles, not years. Similarly, the age of an open-source project is best measured in the number of revisions and commits to the repository. Madagascar revision number 10,000 happened yesterday! The author of the anniversary revision, Professor William Symes from Rice University will receive a commemorative Madagascar stein.

According to Ohloh, the number of commits to the Madagascar repository has increased by 63% in the last 12 months, and the number of contributors during the same period has increased from 21 to 27.

Another old paper is added to the collection of reproducible documents:

Wavefront construction using waverays

A new paper is added to the collection of reproducible documents:

Lowrank finite-differences and lowrank Fourier finite-differences for seismic wave extrapolation in the acoustic approximation

A Winter School on Reproducible Science And Modern Scientific Software, the first of its kind, took place in Norway on January 20-25, 2013, and was organized by SINTEF.

Abstract: A major problem with the computational science community today is that many publications are impossible to reproduce. Results published in a paper are seldom accompanied by the source code used to produce these results. Even when the source code is available the published results can only be reproduced if run the code is compiled with a specific compiler and run on a specific architecture using a specific set of parameters. Reproducibility aims to make the process of publishing reproducible science as simple as possible, and it has gained a lot of momentum as a desirable principle of the scientific method. Tightly coupled with reproducible science is modern software development. Tools and methodologies including version control, unit testing, verification and validation, and continuous integration make the process of publishing reproducible science much simpler.

The winter school will give an introduction to the state-of-the-art in reproducible science and modern scientific software development. The aim is that participants will be able to apply the learned techniques to their own research. Topics that will be covered include reproducible research, verification and validation, software testing, and continuous integration.

Lectures were given by Rasmus Benestad, Andre Brodtkorb, Fernando Perez, and Johan Seland. The school program with supplemental materials is now available on the web.