Shell scripts: pullshots.sh | sort2cmp.sh | supergather.sh | mute.sh

• pullshots.sh - extract shots from tape
• sort2cmp.sh - sort shot-gathered data to cmp-gathered data
• supergather.sh - combine adjacent cmp's to produce a supergather
• mute.sh - mute bad traces and filter if required

  Notes

• The tau-p step (step 2) and source estimation step (step 3) include the designation of a bandpass filter which is applied to the input data. For plotting purposes, it may be beneficial to filter the data with a similar minimum phase band-pass at this point. However, for consistency ensure that the data which is passed to the tau-p or source estimation step is unfilter ( ARG).

2. Tau-p TRANSFORM

Shell script: taup.sh

Prep Phase

• sort original wide-angle data by range (if original ranges are not well sorted) [segy_sort]
• extract data by time and range (e.g., from 1.0 s to 7.0 s and from 0 km to 8 km) [segy_cut]
• band-pass filter (e.g., cut frequencies: 7 - 120 Hz, pass frequencies 10 - 90 Hz, with cosine taper) [segy_ffil]
• decimate to desired sampling interval (e.g., 4 ms) [segy_decim]

Tau-p Phase

• run segy_ftaup with initial guess of scaling factor
  
• use segy_diff to decide an optimal scaling factor (read Korenaga et al. 1997 for this.)
  

Notes on control parameters for segy_ftaup

• Tapering data: only needed at far offsets because data is symmetric about x=0 in the transform.
• Alias removal: uses primary v_p model to interpolate traces
  • -A uses Satish Singh's method.
  • -As uses Jun Korenaga's method. With this option, the time domain convolution is done more carefully than with -A, but the differences are usually quite small (so says Jun Korenaga.)
  • -v - initial v_p model is entered with a pair of parameters per line: h - thickness of the layer in km and c - velocity in km/s [vel.dat].
• Scaling factor: use the equation in Korenaga's JGR paper.
  • -S = 5 makes sense, as this is close to 2pi, or one cycle for a dominant frequency of 45 Hz.
• Differentiation of seismogram :
  • -D - differentiate tau-p seismogram with respect to time. This must be done because the tau-p transform and it's inverse, as defined by Chapman, are not symmetric (i.e., changes in frequency content will occur. Differentiating with respect to time compensates for this. Also, tau-p seismograms generated in the waveform inversion are differentiated, so this maintains consistency.

Current problem with tau-p transform

• Zero or NaN traces: Sometimes, when a supergather is transformed, there is a problem with the resulting tau-p data. When plotting the data using sioseis routines, the data is seen to contain all zeros or NaN's. These same data can be plotted without any trouble using plotsec. However, if this phenomenon occurs then somehow the data is corrupted for the purposes of running the 1D Inversion code. The source of this problem is not yet understood.

3. SOURCE WAVELET ESTIMATE

Shell scripts: srcest.sh | gunsig.sh | rsf.sh

• srcest.sh - calculates a source signature
• gunsig.sh - windows the source signature to a shorter length to be used in following steps.
• rsf.sh - calculates the reflection co-efficient of the seafloor based on near vertical traces.

  Notes on use of srcest.sh

• Choose clean input traces. It is not necessary to use the nearest offset traces, but offset will have an affect so use discretion.
• Be sure to choose appropriate time windows in order to select the ocean bottom reflection and (if applicable) its multiple.
• Use segy_ffil within the script to filter the data to the same bandpass range as that used in the previous step for the tau-p transform.

4. WAVEFORM INVERSION PREPARATION

Shell scripts: vel.sh | src.sh | data.sh | scale.sh

• vel.sh - make velocity file for reflectivity
  • rsb: seabed reflection co-efficient (from vertical incident data).
    This parameter is determined by the shell sccript, rsf.sh in the previous section.
  • dzcomp: compensation interval for density.
    Hamilton's empirical V_p - density - V_s relationship is used to determine density in the model. However, because we can calculate the reflection coefficient of the seafloor, this determines density a priori at the seafloor, given our V_p model. Therefore, we need to set a depth (sub-seafloor) at which the density profile based on the seafloor information is brought into agreement with Hamilton's relationship. In the current example, the value of dzcomp has been set to 0.25 km.
  • files created:
    • fort.1: control parameters for creation of velocity model
    • fort.10: the starting model. The first line contains the number of layers and the succeeding lines contain two columns in a line, i.e., depth and velocity.
  • velprep2 is executed. This creates fort.20 (start.model) which is used as the starting velocity model for succeeding shell scripts.

• src.sh - make source wavelet file
  • transforms source wavelet from time domain to frequency domain.
  • resamples data in time if necessary
  • t0: time shift (preceding first break). Determined by an iterative process involving the comparison of first arrival synthetic seismogram with real data.
  • fact: scaling factor set by matching amplitudes of real and synthetic first arrivals.
  • this shell uses Minshull's program src, because the FFT used is identical to that used in the waveform inversion.

• data.sh - make data input file
  • takes the tau-p data (as generated in the shell taup.sh) and Fourier transforms it into the omega-p (wp) domain.
  • ntd: number of samples (as determined by segy_cut step (li> set t0: start time, and ntf: the number of samples, to appropriate values (the latter should be a power of two so that it can be used in a fast Fourier transform. These two values are input to the Cambridge program fourier which requires that the input data start at t=0.

• scale.sh - determine scaling factor

• input.param - a file containing input parameters

  Care must be taken in setting up the input parameter list that values agree with those used in the valious shell scripts used (especially vel.sh and data.sh.)

   $IN
   NL  =  651  , TH  =  0.002  , THO =  2.70 ,
   XMIN  = -1. , XMAX  = -1.   ,
   NT0   = 1024  , TL0 =  2.048  , TW0 =  2.4,
   NT  = 1024  , TL  =  2.048  ,
   NP0   = 49  , PMIN0 =  0.01 , PMAX0 =  0.49 ,
   NP  = 35  , NPI =  5  , NDP =  1  ,
   IRICK =  0  , RCDEP =  0.0  , VW0 =  1.520  , POL = 1  , 
   ITERMAX =  0  , COST0 =  1.E-10 , CD  =  0.0025 , CM  = 1. ,
   F1  =  0.08 , F2  =  0.14 , F3  =  0.40 , F4  = 0.60 ,
   NORM  = .F. , RADS  = .004  , PCG = .F.   ,
   DGEN  = .F. , VISU  = .F.   , NRUN  =  1  ,
   IOBH  =  0  ,
   $
  

  Parameter Definitions
  Note (WSH 4Oct00): Certain parameters are limited in size. For example, if NT > 1024, the program must be recompiled after editing the value for ntm in inv1d.dim
  • NL = Number of Layers, including half-space - can be determined by running vel.sh
  • TH = Thickness of the layers (should be constant) (km)
  • THO = Thickness of water layer (km) - must match that used in vel.sh
  • XMIN = Minimum offset (km) - not used if set to -1
  • XMAX = Maximum offset (km) - not used if set to -1
    Note: XMIN and XMAX are ordinarily used to determine what part of tau-p space to use in comparison to the synthetic data if the real data are limited in x-range.

  • NT0 = The total number of time points of the data - must be a power of 2
  • TL0 = The total length of the time series (s)
  • TW0 = The time delay (starting time of first sample) - must be the same as t0 in data.sh
  • NT = The total number of time points of the data used for inversion - must be a power of 2
  • TL = The total length of the time series (s) used for inversion

  • NP0 = The total number of traces in the data file
  • PMIN0 = The minimum value of the slowness (s/km)
  • PMAX0 = The maximum value of the slowness (s/km)
  • NP = The total number of traces used for inversion
  • NPI = The number of the first trace used for inversion
  • NDP = The the number of traces skipped between two traces for inversion - to invert all traces, set to 1

  • ITERMAX: The maximum no. of iteration for inversion - 1st test to stop the algorithm. If ITERMAX = 0 then only the forward calculation is performed; there will be no iteration (set this way for the /prep/ stage but not for the /inv/ stage.)
  • COST0 : Cutoff value under which the recursion will stop - second test to stop the algorithm.
  • CD : Covariance of the data. Constant element of the diagonal of the covariance of the data. 0.5 used by Jun Korenaga - determined empirically.
  • CM : Covariance of the model. Constant element of the diagonal of the covariance of the model. 1.0 used by Jun Korenaga - determined empirically.

  • F1-F4 : Definition of the band-pass filter in frequency as a fraction of Nyquist.

       F1  > Lower pass  (0)    1_ _ _  _______  _ _ _1
       F2  >  "  " (1)   /|   |\
       F3  > Upper pass  (1)    0_ ___/_|_ _ _|_\___ _0 ->Freq.
       F4  >  "  " (0) |  | |   | |  |
       as fraction of Fmax 0 F1 F2 F3 F4 Fmax=1
      No filtering if Fi=-1, i=1,2,3 or 4
    

  Other parameters:
  • NORM = .TRUE. Normalisation of the data to adapt to the energy of the synthetic seismogram - ignored in the current version of the program.
  • RADS : Applicable if NORM=.TRUE. Amplitude agreement between the real and synthetic data. Will be equal to unity in proximity to the real model. (RADS.GE.1)
  • PCG = .TRUE. Preconditioning of the gradient to enhance the amplitude of the perturbation of the deeper model.
  • DGEN = .TRUE. Creation of the Omega-P data corresponding to the initial model (written on unit 49) - ignored in the current version of the program
  • VISU = .FALSE.
  • NRUN : Number of the run of inversion with the same data set. (Used only for output information.)

  • IOBH = 1 - use OBH geometry (modification by Korenaga)

5. WAVEFORM INVERSION

Shell scripts: inv.sh | run?.sh

• inv.sh - do-all-the-stuff script. This calls the N run?.sh scripts, where ? is an integer from 1 to N. N is the number of runs desired for a particular inversion. In general, a good approach to the inversion is to start with a small window of slownesses and frequencies and then progress to a full range of slowness and frequencies for the given data set.
• run?.sh - run script for each step. These N files control the parameters for a given inversion run. The values here must mach those found in input.param.run?.

N.B. For each run of inv1d, the following input files are required:
• input.param.run?
• start.model
• wpfile - created by data.sh in ../../prep/run_x.
• src.fft - created by src.sh in ../../prep/run_x.

Notes on input.param.run?

• Be sure that the parameters here correspond to those used for ../../prep/run_x. A few additional concerns are required:
  • NP, NPI, and NDP will need to be redefined (for each run of the inversion).
  • ITERMAX = 9 (anything more than 5 or 6 is likely fine; this can be determined empirically.)
  • F1 ... F4 - use only low frequency part of the data for the first run and expand to higher frequencies for later runs.
  • NZMIN = 16 - this stops the inversion from varying anything abover the 16th layer to avoid near sea floor effects. Be sure to change this parameter if the depth sample interval is changed back in the tau-p transform step.

Calculating the results

Shell scripts: run.sh | misfit.sh

• run.sh - similar to run_x script in ../ but adds conversion of the results to segy format and prints them using plot.sh.
• misfit.sh - calculates the following information and places it in a file named chi.tmp. Requires a file named box.aroundbsr.xtp which windows the data to be used in the calculation.
  • csum_count - number of data
  • csum: - misfit (this is what you want)
  • corr: - correlation.

Plotting the results

Shell script: gmt-taup.wsyn.sh

• gmt-taup.wsyn.sh - GMT script for plotting.