<p>Do your stacking properly. You should not simply sum them, bt since geophysics is a science, please provide a reference for your choice of stacking process.</p>
<p><a href="https://i.sstatic.net/8zsN0.gif" rel="nofollow noreferrer"><img src="https://i.sstatic.net/8zsN0.gif" alt="take into account the offset between geophones!" /></a>
(image from here: <a href="http://www.xsgeo.com/acq.htm" rel="nofollow noreferrer">http://www.xsgeo.com/acq.htm</a> )</p>

CC BY-SA 4.0
2023-05-01T12:52:27.097Z
2023-05-01T12:52:27.097Z
EarlGrey
<p>I have some simple seismic traces as shown here. <a href="https://i.sstatic.net/xDEFP.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/xDEFP.png" alt="Seismic Gather" /></a></p>
<p>When attempting to perform a dispersion analysis using the technique described by this <a href="https://www.kgs.ku.edu/software/surfseis/Publications/ParkEtAl1998.pdf" rel="nofollow noreferrer">paper</a>, I get the following image.</p>
<p><a href="https://i.sstatic.net/L37Pz.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/L37Pz.png" alt="Dispersion Analysis" /></a></p>
<p>Unfortunately, I cannot for the life of me figure out what it is I am doing wrong here. So I feel that I must be missing something conceptually. Below is my python script for processing the dispersion image.</p>
<pre><code>## . . Perform dispersion analysis

## . . Define parameters
ng    = 20       # Number of geophones
nVp   = 200      # Number of Vp test points
nfreq = 200      # Number of frequency test points
ndata = 600      # Number of data points per trace
dt    = 1/16000  # Time_Delta

## . . Define gridspace
freq = np.linspace(0, 40, nfreq)  # Frequency spectrum (Hz)    
Vp   = np.linspace(1, 500, nVp)   # Phase velocity spectrum (m/s)
X, Y = np.meshgrid(freq, Vp)

## . . Setup gridspace for the slant stack function
S = np.zeros((ng, nfreq, nVp, ndata), dtype=complex)

## . . Calculate dispersion with varying phase velocity and frequency
for idx, trace in enumerate(seismic_gather[0].trace):
  ## Get the fourier transformation and normalization of the waveform
  U = np.fft.fft(trace.data)  # Fourier transform 
  N = U/np.abs(U)             # Normalize waveform

  ## Apply dynamic linear moveout
  P = 2*np.pi*X*trace.offset/Y             # Dispersion property
  S[idx] = np.exp(1j*P)[:,:,np.newaxis]*N

  
## . . Stack traces
S = np.sum(S, 0)

## . . Extract amplitudes
A = np.abs(S)/ng
</code></pre>
<p>To describe the code a little, I am defining a meshgrid space of 200 testing points for the Phase Velocity and Frequency. I then take my seismic data and transform it into the frequency domain and normalize it. Once normalized, I apply the dynamic linear moveout to each trace, and then stack the the traces for each pair of phase velocity and frequency points along the meshgrid.</p>
<p>Unfortunately, the resulting image doesn't appear to be what I expect, which I expect to be similar to the paper linked. Am I missing a step here?</p>

<geophysics><python>
I am trying to write a simple python script for dispersion analysis of surface waves, but I think I am missing something conceptually
TylerSingleton
# I am trying to write a simple python script for dispersion analysis of surface waves, but I think I am missing something conceptually

(Asked by: username_0 on 2023-04-08)

I have some simple seismic traces as shown here. 
When attempting to perform a dispersion analysis using the technique described by this paper, I get the following image.

Unfortunately, I cannot for the life of me figure out what it is I am doing wrong here. So I feel that I must be missing something conceptually. Below is my python script for processing the dispersion image.
## . . Perform dispersion analysis

## . . Define parameters
ng    = 20       # Number of geophones
nVp   = 200      # Number of Vp test points
nfreq = 200      # Number of frequency test points
ndata = 600      # Number of data points per trace
dt    = 1/16000  # Time_Delta

## . . Define gridspace
freq = np.linspace(0, 40, nfreq)  # Frequency spectrum (Hz)    
Vp   = np.linspace(1, 500, nVp)   # Phase velocity spectrum (m/s)
X, Y = np.meshgrid(freq, Vp)

## . . Setup gridspace for the slant stack function
S = np.zeros((ng, nfreq, nVp, ndata), dtype=complex)

## . . Calculate dispersion with varying phase velocity and frequency
for idx, trace in enumerate(seismic_gather[0].trace):
  ## Get the fourier transformation and normalization of the waveform
  U = np.fft.fft(trace.data)  # Fourier transform 
  N = U/np.abs(U)             # Normalize waveform

  ## Apply dynamic linear moveout
  P = 2*np.pi*X*trace.offset/Y             # Dispersion property
  S[idx] = np.exp(1j*P)[:,:,np.newaxis]*N

  
## . . Stack traces
S = np.sum(S, 0)

## . . Extract amplitudes
A = np.abs(S)/ng

To describe the code a little, I am defining a meshgrid space of 200 testing points for the Phase Velocity and Frequency. I then take my seismic data and transform it into the frequency domain and normalize it. Once normalized, I apply the dynamic linear moveout to each trace, and then stack the the traces for each pair of phase velocity and frequency points along the meshgrid.
Unfortunately, the resulting image doesn't appear to be what I expect, which I expect to be similar to the paper linked. Am I missing a step here?

**Tags:** geophysics, python


---

## 1 answers

### Answer by username_1 on 2023-05-01. Score: -1

Do your stacking properly. You should not simply sum them, bt since geophysics is a science, please provide a reference for your choice of stacking process.

[![take into account the offset between geophones!](https://i.sstatic.net/8zsN0.gif)](https://i.sstatic.net/8zsN0.gif)
(image from here: <http://www.xsgeo.com/acq.htm> )


---

