Model fitting by least squares

FLATTENING 3-D SEISMIC DATA

Here we follow the doctoral dissertation of Jesse Lomask. In Figure  we have seen how to flatten 2-D seismic data. The 3-D process is much more interesting. To see why, consider this: Starting from the origin $(x , y ) = (0, 0)$ we move along the $x$-axis flattening until we come to $(10, 0)$. From there we move along the $y$-axis flattening until we get to $(x , y ) = (10, 10)$. Move backwards on the $x$-axis to $(0, 10)$ flattening as you go. Finally, return to the origin. Along our journey around this square we have integrated $p = d t /d x$ (and $d t /d y$) to find the total time shift. Upon returning to the starting point, we would like the total time shift to return to zero. Dealing with real data of less than perfect coherence this might not happen. Old time seismologists would say, ``The survey lines don't tie.'' As we push to the limits of our knowledge (which we normally do) this problem always arises. We would like a solution that gives the best fit of all the data in a volume. Given a volume of data $u (t , x , y )$ we seek the best $\tau (x , y )$ such that $w(t , x , y )= u (t-\tau (x , y ), x , y )$ is flattened. Let's get it.

Here is an expression that on first sight seems to say nothing

$$\nabla \tau \eq \left[ \begin{array}{c} \frac{\partial \tau... ...} \\ \\ \frac{\partial \tau}{\partial y} \end{array}\right]$$ (86)

Equation (86) looks like a tautology, a restatement of basic mathematical notation. This is so, however, only if $\tau (x , y )$ is known and the derivatives are derived from it. When $\tau (x , y )$ is not known but the partial derivatives are observed, then we have two measurements at each $(x , y )$ location for the one unknown $\tau$ at that location. The same is true at all locations, so we write it as a regression, a residual ${\bf r}$ that we will work to get small to find a best fitting $\tau (x , y )$ or maybe $\tau (x , y , t )$. Let $d$ be the measurements in the vector in equation (86), the measurements throughout the $(t , x , y )$-volume. Expressed as a regression equation (86) becomes

$${\bf0} \quad \approx\quad {\bf r} \eq \nabla \tau \ -\ {\bf d}$$ (87)

chev
Figure 8. Chevron data cube from the Gulf of Mexico. A salt dome (lower left corner in the top plane) has pushed upwards, dragging bedding planes (seen in the bottom two orthogonal planes) along with it.

Let us see how the coming 3-D illustrations were created. First we need code for vector gradient with its adjoint, negative vector divergence. Here it is: api/c/igrad2.c

    for (i2=0; i2 < n2-1; i2++) {  
	for (i1=0; i1 < n1-1; i1++) {
	    i = i1+i2*n1;
	    if (adj) {
		p[i+1]  += r[i]; 
		p[i+n1] += r[i+n12];
		p[i]    -= (r[i] + r[i+n12]);
	    } else {
		r[i]     += (p[i+1]  - p[i]); 
		r[i+n12] += (p[i+n1] - p[i]);
	    }
	}
    }

In a kind of magic, all we need to fit our regression (86) is to pass the igrad2 module to the Krylov subspace solver, simple solver using cgstep, but first we need to compute d by calculating $dt/dx$ and $d t /d y$ between all the mesh points.

do iy=1,ny { # Calculate x-direction dips: px 
             call puck2d(dat(:,:,iy),coh_x,px,res_x,boxsz,nt,nx) 
           } 
do ix=1,nx { # Calculate y-direction dips: py 
             call puck2d(dat(:,ix,:),coh_y,py,res_y,boxsz,nt,ny) 
           } 
do it=1,nt { # Integrate dips: tau 
             call dipinteg(px(it,:,:),py(it,:,:),tau,niter,verb,nx,ny) 
           }

The code says first to initialize the gradient operator. Convert the 2-D plane of $dt/dx$ to a vector. Likewise $d t /d y$. Concatenate these two vectors into a single column vector d like in equation (86). Tell the simple solver to make its steps with to use cgstep with the linear operator igrad2.

Model fitting by least squares