KIRCHHOFF INTEGRAL AND PORTER-BOJARSKI INTEGRAL

The basis of the iterative methods is the non-iterative ``deterministic'' approach to seismic inversion. It ignores the random components of seismic data and is strictly speaking less deterministic than the iterative inversion. The reason it is sometimes considered as deterministic is that it is based on the forward scattering problem that can be solved exactly by the Kirchhoff integral. People took the forward solutions and simply changed the direction of the wave propagation. The resultant integrals converge only between finite integration limits that cut off evanescent energy. Therefore, this method constructs images by backpropagating only the homogeneous part of the registered wavefield. Iterative methods treat evanescent energy at least statistically so that they should lead to more accurate images.

For didactical reasons we will now consider the Kirchhoff integral in its simplest form: for an homogeneous acoustic medium. The following concepts can be extended to inhomogeneous elastic media.

The Kirchhoff integral takes a wavefield $p(R,\omega)$ and its normal derivative ${\frac{\partial}{\partial{\underline {n}}}}p(R,\omega)$ measured on an arbitrary closed surface S surrounding the sources and propagates the field forward, away from the sources to an arbitrary place R₀ outside the surface. The normal vector $\underline {n}$ points to the outside of V that is away from the sources. The volume V is enclosed by S and contains the source volume V_S. Inside of V the Kirchhoff integral is zero. A Green's function $G(R-R_0,\omega)$ extrapolates the registered wave field from R to R₀.

If R₀ is inside of V we have:

$$\oint_{S} [p(R,\omega) {\frac{\partial}{\partial{\underline ... ...}{\partial{\underline {n}}}}p(R,\omega) G(R-R_0,\omega)] dS = 0$$ (2)

If R₀ is outside of V we have:

$$\oint_{S} [p(R,\omega) {\frac{\partial}{\partial{\underline ... ...underline {n}}}}p(R,\omega) G(R-R_0,\omega)] dS = p(R_0,\omega)$$ (3)

The Kirchhoff integral has to be zero inside V because we don't know where the sources are. If we knew the exact position of the source volume V_S, we could change the sign of $\underline {n}$ and propagate until close to that region. But there we had to stop because the source volume contains singularities that would make our integral diverge. Mathematically, it makes only a quantitative difference whether the sources inside are primary or secondary sources. In any case, we will run into trouble if we try to backpropagate evanescent energy, that is, if we come too close to the sources. Numerically, the whole thing could work if only weak scatterers have to be considered.

The physical interpretation of the Kirchhoff integral is that of a forward propagation and superposition of Huygen's sources on S. We assume two kinds of sources: monopoles $p(R,\omega)$and dipoles ${\frac{\partial}{\partial{\underline {n}}}}p(R,\omega)$. In case of stiff boundary conditions (von Neumann conditions) on S, we can drop the term ${\frac{\partial}{\partial{\underline {n}}}}p(R,\omega)$ and for a weak boundary S only the gradient ${\frac{\partial}{\partial{\underline {n}}}}p(R,\omega)$ is of importance and we drop the term containing $p(R,\omega)$ (Dirichlet boundary condition).

So far we learned how to forward propagation. For backpropagation, i.e. if we propagate towards the sources, the Kirchhoff integral has to be modified. Two changes have to be made:

1.

change the sign of the normal vector $\underline {n}$; i.e., let $\underline {n}$ point towards the source region.

2.

change the propagation direction, i.e., the Green's function G.

To reverse the wave propagation direction, we can change the sign of time in the Green's function G taken in the time domain. Because we are working in the frequency domain, we have to take the complex conjugate Green's function G^*. G^* performs the propagation of an imploding wavefield while G propagates an exploding wavefield. The two substitutions yield to the Porter-Bojarski integral (Langenberg, K.L., 1986)

$$- \oint_{S} [p(R,\omega) {\frac{\partial}{\partial{\underlin... ...e {n}}}}p(R,\omega) G^*(R-R_0,\omega)] dS = p(R_0,\omega) \ \ .$$ (4)

The implosion will take place at the source location and at the time the measured field exploded.

G and G^* solve the inhomogeneous wave equation for a point source:

$$\Delta G(R-R_0,\omega) - {\frac{\partial{^2}}{\partial{t^2}}}G (R-R_0,\omega) = - \delta(R-R_0)$$ (5)

and

$$\Delta G^*(R-R_0,\omega) - {\frac{\partial{^2}}{\partial{t^2}}}G^* (R-R_0,\omega) = - \delta(R-R_0)$$ (6)

The difference between the two inhomogeneous wave equations describes the homogeneous wavefield

$$\Delta (G^* - G) - {\frac{\partial{^2}}{\partial{t^2}}}(G^*-G) = 0 \ \ .$$ (7)

In a three-dimensional medium, the new Green's function is

$$\begin{eqnarray} G^* - G = & \frac{\exp (-i \omega (t- \frac{\vert R-R_{0}\vert}... ... \frac{\vert R-R_{0}\vert}{v}))}}{2 \pi \vert R-R_{0}\vert} \ \ . \end{eqnarray}$$ (8)

The Green's function of the homogeneous wavefield that is backpropagated has an interesting interpretation that is more obvious in its 3D time domain representation

$$G^*(R-R_0,t) - G(R-R_0,t) = {\frac{\delta (t-{\frac{\vert R-... ...t+{\frac{\vert R-R_0\vert}{v}})}{4 \pi \vert R-R_0\vert}} \ \ .$$ (9)

Our Green's function propagates the wavefield of time symmetric point sources on the surface S. The time symmetric waves implode for negative times, run through a focus at t=0, and explode for t>0.

The homogeneous wavefield is used in standard seismic imaging where we have to superpose the forward extrapolated wavefield and the backward extrapolated waves to get a representation of the subsurface by homogeneous waves at the time they were scattered from a discontinuity. Usually we propagate with the same sign of $\underline {n}$ for both forward and backward propagating Green's functions. This explains why we subtracted the forward and backward extrapolated waves, whereas in prestack migration, they have to be added.

The surface integral of Porter and Bojarski also has a volume integral representation

$$\int \int_{V} \int q(R,\omega) [G^*(R-R_0,\omega) - G(R-R_0,\omega)] dV = p(R_0,\omega)$$ (10)

where $q(R,\omega)$ is a source at R inside of V. This tells us that the scattered wavefield at R₀ can be obtained by the superposition of volume sources located in the volume enclosed by S. G-G^* extrapolates only an homogeneous wavefield.