next up previous [pdf]

Next: Completing the assignment Up: Homework 5 Previous: Homework 5

Computational part

  1. The first example is the model from Homework 1 with a hyperbolic reflector under a constant velocity layer. The model is shown in Figure 6a. Figure 2 shows a shot gather modeled at the surface and extrapolated to a level of 1 km in depth using two extrapolation operators:
    1. The exact phase shift filter
      \begin{displaymath}
\hat{U}(z,k,\omega) = \hat{U}(0,k,\omega) e^{i \sqrt{S^2 \omega^2 - k^2} z}\;.
\end{displaymath} (1)

    2. Its approximation
      \begin{displaymath}
\hat{U}(z,k,\omega) \approx \hat{U}(0,k,\omega) 
e^{i ...
... \left(\cos{(k \Delta x)}-1\right) z}{2 (\Delta x)^2}}\;.
\end{displaymath} (2)

      Approximation (2) is suitable for an implementation in the space domain with a digital recursive filter. However, its accuracy is limited, which is evident both from Figure 2 and from Figure 3, which compares the phases of the exact and the approximate extrapolators. We can see that approximation (2) is accurate only for small angles from the vertical $\theta$, defined by

      \begin{displaymath}
\theta = \arcsin\left(\frac{k}{S\,\omega}\right)\;.
\end{displaymath}

    Your task: Design an approximation that would be more accurate than approximation (2). Your approximation should be suitable for a digital filter implementation in the space domain. Therefore, it can involve $k$ only through $\cos{(n\,k\,\Delta
x)}$ functions with integer $n$.

    1. Change directory
      cd hw5/hyper
      
    2. Run
      scons view
      
      to generate figures and display them on your screen.
    3. Edit the SConstruct file to change the approximate extrapolator.
    4. Run
      scons view
      
      again to observe the differences.

    model
    model
    Figure 1.
    Synthetic velocity model with a hyperbolic reflector.
    [pdf] [png] [scons]

    shot
    shot
    Figure 2.
    Synthetic shot gather. Left: Modeled for receivers at the surface. Middle: Receivers extrapolated to 1 km in depth with an exact phase-shift extrapolation operator. Right: Receivers extrapolated to 1 km in depth with an approximate extrapolation operator.
    [pdf] [png] [scons]

    phase
    Figure 3.
    Phase of the extrapolation operator at 5 Hz frequency as a function of the wave propagation angle. Solid line: exact extrapolator. Dashed line: approximate extrapolator.
    phase
    [pdf] [png] [scons]

    from rsf.proj import *
    from math import pi
    
    # Make a reflector model 
    Flow('refl',None,
         '''
         math n1=10001 o1=-4 d1=0.001 output="sqrt(1+x1*x1)" 
         ''')
    Flow('model','refl',
         '''
         unif2 n1=201 d1=0.01 v00=1,2 |
         put label1=Depth unit1=km label2=Lateral unit2=km
         label=Velocity unit=km/s | 
         smooth rect1=3
         ''')
    Result('model',
           '''
           window min2=-1 max2=2 j2=10 |
           grey allpos=y title=Model 
           scalebar=y barreverse=y 
           ''')
    
    # Reflector dip
    Flow('dip','refl','math output="x1/input" ')
    
    # Kirchoff modeling
    Flow('shot','refl dip',
         '''
         kirmod nt=1501 twod=y vel=1 freq=10  
         ns=1 s0=1 ds=0.01 nh=801 dh=0.01 h0=-4
         dip=${SOURCES[1]} verb=y |
         costaper nw2=200
         ''')
    Plot('shot',
         '''
         window min2=-2 max2=2 min1=1 max1=4 |
         grey title="0 km Depth" 
         label1=Time unit1=s label2=Offset unit2=km
         labelsz=10 titlesz=15 grid=y
         ''')
    
    # Double Fourier transform
    
    Flow('kw','shot','fft1 | fft3 axis=2')
    
    # Extrapolation filters
    
    dx = 0.01
    dz = 1
    
    dx2 = 2*pi*dx
    dz2 = 2*pi*dz
    a = 0.5*dz2/(dx2*dx2)
    
    # In the expressions below,
    # x1 refers to omega and x2 refers to k
    
    Flow('exact','kw',
         'math output="exp(I*sqrt(x1*x1-x2*x2)*%g)" ' % dz2)
    
    Flow('approximate','kw',
         '''
         window f1=1 |
         math output="exp(I*x1*%g)* 
         (x1+I*(cos(x2*%g)-1)*%g)/
         (x1-I*(cos(x2*%g)-1)*%g)" |
         pad beg1=1 
         ''' % (dz2,dx2,a,dx2,a))
    
    # Extrapolation
    
    for case in ('exact','approximate'):
        Flow('phase-'+case,case,
             '''
             window n1=1 min1=5 min2=-2 max2=2 | 
             math output="log(input)" | sfimag
             ''')
        Flow('angle-'+case,'phase-'+case,
             '''
             math output="%g*asin(x1/5)" |
             cmplx $SOURCE
             ''' % (180/pi))
    
        Flow('shot-'+case,['kw',case],
             '''
             mul ${SOURCES[1]} | 
             fft3 axis=2 inv=y | fft1 inv=y
             ''')
        Plot('shot-'+case,
             '''
             window min2=-2 max2=2 min1=1 max1=4 |
             grey title="%g km Depth (%s)" 
             label1=Time unit1=s label2=Offset unit2=km
             labelsz=10 titlesz=15 grid=y
             ''' % (dz, case.capitalize()))
    
    Result('shot','shot shot-exact shot-approximate',
           'SideBySideAniso')
    
    Result('phase','angle-exact angle-approximate',
           '''
           cat axis=2 ${SOURCES[1]} | 
           graph title=Phase dash=0,1 
           label1="Angle From Vertical" unit1="ô_"
           ''')
    
    End()
    

  2. The second example uses the Sigsbee synthetic velocity model which you encountered in Homework 4. The model is shown in Figure 4. We will select one slice of the model (at 15 kft depth) to analyze different one-way wavefield extrapolators. Figure 5 compares the phases of two one-way wave extrapolation operators:
    1. The exact non-stationary phase shift filter
      \begin{displaymath}
E(k,x) = e^{i \sqrt{S^2(x) \omega^2 - k^2} \Delta z}
\end{displaymath} (3)

    2. Its split-step approximation
      \begin{displaymath}
E(k,x) \approx e^{i S(x) \omega \Delta z} e^{i \sqrt{S_0^2 \omega^2 - k^2} \Delta z-i S_0 \omega \Delta z}
\end{displaymath} (4)

      Approximation (4) can be implemented efficiently, because it involves only diagonal matrix multiplications in space and wavenumber domains. However, its accuracy is limited.

    Your task: Design an approximation that would be more accurate than approximation (4). Your approximation should be suitable for an efficient implementation. It can involve products of functions of $x$ and functions of $k$ but not any functions that mix $x$ and $k$.

    1. Change directory
      cd hw5/sigsbee
      
    2. Run
      scons view
      
      to generate figures and display them on your screen.
    3. Edit the SConstruct file to change the approximate extrapolator.
    4. Run
      scons view
      
      again to observe the differences.

    vel
    vel
    Figure 4.
    Sigsbee velocity model. A slice of the model at 15 kft depth is selected for analyzing wavefield extrapolation operators.
    [pdf] [png] [scons]

    phase2
    phase2
    Figure 5.
    Phase of the wavefield extrapolation operators for a depth slice at 15 kft and frequency of 5 Hz as a function of space $x$ and wavenumber $k$. Top: Exact mixed-domain extrapolator. Bottom: Split-step approximation.
    [pdf] [png] [scons]

    from rsf.proj import *
    from math import pi
    
    # Download velocity model from the data server
    ##############################################
    vstr = 'sigsbee2a_stratigraphy.sgy'
    Fetch(vstr,'sigsbee')
    Flow('zvstr zvstr.asc zvstr.bin',vstr,
         '''
         segyread read=data
         hfile=${TARGETS[1]}
         bfile=${TARGETS[2]}
         ''')
    
    Flow('vel','zvstr',
         '''
         put d1=0.025 o2=10.025 d2=0.025 |
         window f1=200 | scale dscale=0.001
         ''')
    Plot('vel',
         '''
         grey wanttitle=n color=j allpos=y 
         label1=Depth unit1=kft label2=Lateral unit2=kft
         screenratio=0.3125 screenht=4 labelsz=4
         ''')
    
    # Take a slice at 15 kft
    ########################
    
    Flow('slice','vel','window n1=1 min1=15')
    
    Plot('line','slice',
         '''
         math output=15 | graph min2=5 max2=30 yreverse=y 
         pad=n wanttitle=n wantaxis=n plotcol=7 plotfat=10
         screenratio=0.3125 screenht=4
         ''')
    Result('vel','vel line','Overlay')
    
    # Compute extrapolation matrix
    ##############################
    
    w = 5      # frequency
    dz = 0.025 # depth step
    v0 = 9.38  # mean velocity
    
    # x-k plane
    Flow('xk','slice',
         '''
         spray axis=2 n=201 d=0.01 o=-1 | rtoc | 
         put label1=x unit1=kft label2="k_x" unit2=1/kft
         ''')
    
    Flow('Exact','xk',
         '''
         math output="exp(I*(sqrt((%g/input)^2-x2^2)*%g))" 
         ''' % (w,2*pi*dz))
    
    Flow('Split-step','xk',
         '''
         math output="exp(I*((%g/input+sqrt(%g^2-x2^2)-%g)*%g))" 
         ''' % (w,w/v0,w/v0,2*pi*dz))
    
    for case in ('Exact','Split-step'):
        Plot(case,
             '''
             math output="log(input)" | imag | 
             grey transp=n title=%s allpos=y
             labelsz=10 titlesz=12
             ''' % case)
    Result('phase2','Exact Split-step','OverUnderAniso')
    
    End()
    

  3. In the third example, we return to the synthetic model shown in Figure 6a to experiment with modeling and migration by the phase-shift method in a $V(z)$ medium.

    Figure 6b shows an idealized image (band-passed reflectivity) for the synthetic model. In ``exploding-reflector'' modeling, this image is assumed to be the seismic wavefield frozen at zero time. The modeling program extrapolates the wavefield to the surface and is implemented in phaseshift.c. The program operates in the frequency-wavenumber domain and establishes a linear relationship between reflectivity as a function of depth and zero-offset data as a function of frequency for one space wavenumber. If we think of this linear relationship as a matrix multiplication, we can associate forward modeling with multiplication

    \begin{displaymath}
\mathbf{d} = \mathbf{A}\,\mathbf{m}
\end{displaymath} (5)

    and migration with the adjoint multiplication
    \begin{displaymath}
\widehat{\mathbf{m}} = \mathbf{A}^T\,\mathbf{d}\;,
\end{displaymath} (6)

    where $\mathbf{A}^T$ is the conjugate transpose of $\mathbf{A}$.

    Figure 7 shows the result of forward modeling (multiplication by $\mathbf{A}$) after inverse Fourier transform to time and space. Your task is to implement the corresponding migration (multiplication by $\mathbf{A}^T$).

    Instead of simply applying the adjoint operator, we can also try to compute the least-squares inverse

    \begin{displaymath}
\widehat{\mathbf{m}} = \left(\mathbf{A}^T\,\mathbf{A}\right)^{-1}\,\mathbf{A}^T\,\mathbf{d}\;,
\end{displaymath} (7)

    which corresponds to the minimum of the least-squares misfit function $\vert\mathbf{A}\,\mathbf{m}-\mathbf{d}\vert^2$. In practice, inversion in equation (7) is implemented with an iterative conjugate-gradient algorithm, which applies $\mathbf{A}$ and $\mathbf{A}^T$ (modeling and migration) at each iteration without the need to form any matrices explicitly.

    model expl
    model,expl
    Figure 6.
    (a) Synthetic model: curved reflectors in a $V(z)$ velocity. (b) ``Exploding reflector'', an ideal image of band-passed reflectivity.
    [pdf] [pdf] [png] [png] [scons]

    modl
    modl
    Figure 7.
    Zero-offset data generated by exploding-reflector phase-shift modeling.
    [pdf] [png] [scons]

    1. Change directory
      cd hw5/lsmig
      
    2. Run
      scons view
      
      to generate the figures and display them on your screen.
    3. Modify the program in the phaseshift.c file to fill the missing part and to implement phase-shift migration as the adjoint of phase-shift modeling.
    4. Modify the SConstruct file to uncomment the part related to migration. Check your result by running
      scons migr.view
      
    5. Test if your migration code is truly the adjoint of modeling by running the dot-product test
      sfcdottest ./phaseshift.exe mod=kexpl.rsf dat=kmodl.rsf \
      vel=vofz.rsf nw=247 dw=0.16276
      
      On a machine with multiple CPUs, you can also try
      sfcdottest sfomp ./phaseshift.exe split=2 \
      mod=kexpl.rsf dat=kmodl.rsf vel=vofz.rsf nw=247 dw=0.16276
      
      If your adjoint is correct, you should see two identical sets of numbers, such as
      sfcdottest:  L[m]*d=(25264,-20273.9)
      sfcdottest: L'[d]*m=(25264,-20273.9)
      
      Your actual numbers will be different because of random input vectors but they should be the same between L[m]*d and L'[d]*m.
    6. Now you are ready for testing least-squares migration. Uncomment the corresponding lines in SConstruct and run
      scons invs.view
      
      Do you notice any difference with the previous result? Increase the number of iterations from 10 to 100 and repeat the experiment.
    7. Include your figures in this document by editing hw5/paper.tex.

    from rsf.proj import *
    
    # Generate a reflector model
    
    layers = (
         ((0,2),(3.5,2),(4.5,2.5),(5,2.25),
           (5.5,2),(6.5,2.5),(10,2.5)),
         ((0,2.5),(10,3.5)),
         ((0,3.2),(3.5,3.2),(5,3.7),
           (6.5,4.2),(10,4.2)),
         ((0,4.5),(10,4.5))
    )
    
    nlays = len(layers)
    for i in range(nlays):
         inp = 'inp%d' % i
         Flow(inp+'.asc',None,
              '''
              echo %s in=$TARGET
              data_format=ascii_float n1=2 n2=%d
              ''' %            (' '.join(map(lambda x: ' '.join(map(str,x)),
                               layers[i])),len(layers[i])))
    
    dim1 = 'o1=0 d1=0.05 n1=201'
    Flow('lay1','inp0.asc',
         'dd form=native | spline %s fp=0,0' % dim1)
    Flow('lay2',None  ,
         'math %s output="2.5+x1*0.1" '      % dim1)
    Flow('lay3','inp2.asc',
         'dd form=native | spline %s fp=0,0' % dim1)
    Flow('lay4',None  ,'math %s output=4.5'  % dim1)
    
    Flow('lays','lay1 lay2 lay3 lay4',
         'cat axis=2 ${SOURCES[1:4]}')
    
    graph = '''
    graph min1=2.5 max1=7.5 min2=0 max2=5
    yreverse=y wantaxis=n wanttitle=n screenratio=1
    '''
    Plot('lays0','lays',graph + ' plotfat=10 plotcol=0')
    Plot('lays1','lays',graph + ' plotfat=2 plotcol=7')
    Plot('lays2','lays',graph + ' plotfat=2')
    
    # Velocity
    
    Flow('vofz',None,
         '''
         math output="1.5+0.25*x1"
         d2=0.05 n2=201 d1=0.01 n1=501
         label1=Depth unit1=km
         label2=Distance unit2=km
         label=Velocity unit=km/s
         ''')
    Plot('vofz',
         '''
         window min2=2.75 max2=7.25 |
         grey color=j allpos=y bias=1.5
         title=Model screenratio=1
         ''')
    
    Result('model','vofz lays0 lays1','Overlay')
    
    # Exploding reflector
    
    Flow('expl','lays vofz',
         '''
         unif2 d1=0.01 n1=501 v00=1,2,3,4,5 |
         depth2time velocity=${SOURCES[1]} |
         ai2refl | ricker1 frequency=10 |
         time2depth velocity=${SOURCES[1]} |
         put label1=Depth unit1=km 
         label2=Distance unit2=km
         ''')
    
    Result('expl',
           '''
           window max1=5 min2=2.5 max2=7.5 | 
           grey title="Exploding Reflector" screenratio=1
           ''')
    
    # Modeling
    ##########
    
    proj = Project()
    prog = proj.Program('phaseshift.c')
    
    # Cosine Fourier transform
    Flow('kexpl','expl','cosft sign2=1 | rtoc')
    
    Flow('kmodl','kexpl %s vofz' % prog[0],
         '''
         ./${SOURCES[1]} vel=${SOURCES[2]} 
         nw=247 dw=0.16276
         ''',split=[2,'omp',[0]])
    
    Flow('modl','kmodl',
         'pad n1=769 | fft1 inv=y | cosft sign2=-1')
    
    Result('modl',
           '''
           grey title="Exploding Reflector Modeling" 
           label1=Time unit1=s
           ''')
    
    # Migration
    ###########
    
    Flow('kdata','modl',
         'cosft sign2=1 | fft1 | window n1=247')
    
    Flow('kmigr','kdata %s vofz' % prog[0],
         './${SOURCES[1]} vel=${SOURCES[2]} adj=1',
         split=[2,'omp',[0]])
    
    Flow('migr','kmigr','real | cosft sign2=-1')
    
    # !!! UNCOMMENT BELOW !!!
    
    #Result('migr',
    #       '''
    #       window max1=5 min2=2.5 max2=7.5 | 
    #       grey title="Exploding Reflector Migration" 
    #       screenratio=1 label1=Depth unit1=km
    #       ''')
    
    # Least-Squares Migration
    #########################
    
    Flow('kinvs','kdata %s vofz kmigr' % prog[0],
         '''
         cconjgrad omp ./${SOURCES[1]} split=2 
         nw=247 dw=0.16276
         vel=${SOURCES[2]} mod=${SOURCES[3]} niter=10
         ''')
    
    Flow('invs','kinvs','real | cosft sign2=-1')
    
    # !!! UNCOMMENT BELOW !!!
    
    #Result('invs',
    #       '''
    #       window max1=5 min2=2.5 max2=7.5 | 
    #       grey title="Exploding Reflector LS Migration" 
    #       screenratio=1 label1=Depth unit1=km
    #       ''')
    
    
    End()
    

    /* Phase-shift modeling and migration */
    
    #include <rsf.h>
    
    int main(int argc, char* argv[])
    {
        bool adj;
        int ik, iw, nk, nw, iz, nz;
        float k, dk, k0, dw, dz, z0, *v, eps;
        sf_complex *dat, *mod, w2, ps, wave;
        sf_file inp, out, vel;
    
        sf_init(argc,argv);
    
        if (!sf_getbool("adj",&adj)) adj=false;
        /* adjoint flag, 0: modeling, 1: migration */
    
        inp = sf_input("in");
        out = sf_output("out");
        vel = sf_input("vel"); /* velocity file */
    
        if (!sf_histint(vel,"n1",&nz)) 
    	sf_error("No n1= in vel");
        if (!sf_histfloat(vel,"d1",&dz)) 
    	sf_error("No d1= in vel");
        if (!sf_histfloat(vel,"o1",&z0)) z0=0.0;
        
        if (!sf_histint(inp,"n2",&nk)) 
    	sf_error("No n2= in input");
        if (!sf_histfloat(inp,"d2",&dk)) 
    	sf_error("No d2= in input");
        if (!sf_histfloat(inp,"o2",&k0)) 
    	sf_error("No o2= in input");
    
        dk *= 2*SF_PI;
        k0 *= 2*SF_PI;
    
        if (adj) { /* migration */
    	if (!sf_histint(inp,"n1",&nw)) 
    	    sf_error("No n1= in input");
    	if (!sf_histfloat(inp,"d1",&dw)) 
    	    sf_error("No d1= in input");
    
    	sf_putint(out,"n1",nz);
    	sf_putfloat(out,"d1",dz);
    	sf_putfloat(out,"o1",z0);
        } else {  /* modeling */
    	if (!sf_getint("nw",&nw)) sf_error("No nw=");
    	if (!sf_getfloat("dw",&dw)) sf_error("No dw=");
    
    	sf_putint(out,"n1",nw);
    	sf_putfloat(out,"d1",dw);
    	sf_putfloat(out,"o1",0.0);	
        }
    
        if (!sf_getfloat("eps",&eps)) eps = 1.0f;
        /* stabilization parameter */
    
        dw *= 2*SF_PI;
    
        dat = sf_complexalloc(nw);
        mod = sf_complexalloc(nz);
        
        /* read velocity, convert to slowness squared */
        v = sf_floatalloc(nz);
        sf_floatread(v,nz,vel);
        for (iz=0; iz < nz; iz++) { 
    	v[iz] = 2.0f/v[iz];
    	v[iz] *= v[iz];
        }
    
        for (ik=0; ik<nk; ik++) { /* wavenumber */
    	sf_warning("wavenumber %d of %d;",ik+1,nk);
    	k=k0+ik*dk;
    	k *= k;
    	
    	if (adj) {
    	    sf_complexread(dat,nw,inp);
    	    for (iz=0; iz < nz; iz++) { 
    		mod[iz]=sf_cmplx(0.0,0.0);
    	    }
    	} else {
    	    sf_complexread(mod,nz,inp);
    	}
    
    	for (iw=0; iw<nw; iw++) { /* frequency */
    	    w2 = sf_cmplx(eps*dw,iw*dw);
    	    w2 *= w2;
    
    	    if (adj) { /* migration */
    		wave = dat[iw];
    		for (iz=0; iz < nz; iz++) { 
    		    /* !!! FILL MISSING LINES !!! */
    		}
    	    } else {  /* modeling */
    		wave = mod[nz-1];
    		for (iz=nz-2; iz >=0; iz-) {
    		    ps = cexpf(-csqrt(w2*v[iz]+k)*dz);
    		    wave = wave*ps+mod[iz];
    		}
    		dat[iw] = wave;
    	    }
    	}
    	
    	if (adj) {
    	    sf_complexwrite(mod,nz,out);
    	} else {
    	    sf_complexwrite(dat,nw,out);
    	}
        } 
        sf_warning(".");
    
        exit(0);
    }
    


next up previous [pdf]

Next: Completing the assignment Up: Homework 5 Previous: Homework 5

2013-11-07