Project topics and materials
Get Complete Project Material File(s) Now! »
Chapter 5 Computational analysis
In the present chapter we shall briefly outline certain aspects of the numerical and computational methods used. Furthermore we shall present the statistical approach used in the simultaneous extraction of the reflectivity sequence and seismic wavelet.
Single channel inversion
The Marchenko integral equation will be solved for the kernel K(x, y) on the triangular region −x ≤ y ≤ x in the xy – plane by applying the orthogonal collocation method. Given a value of x, K(x, y) can be expanded in terms of N localized functions. The building blocks of our expansion are the Hermite interpolation polynomials that allow us to fix the value of the function to be interpolated and its derivative at the beginning and the end of any interval. For the unit interval, these polynomials are defined in By subdividing the interval [−x, x] into n equidistant sub-intervals, we can expand the output kernel in terms of the above interpolation polynomials adapted to each of those sub-intervals. Due to the requirements of continuity, this is equivalent to expanding K in terms of functions gj (y) that are spliced together by joining interpolation poly-nomials on adjacent sub-intervals, as shown in Fig. (5.2). The splicing together of the interpolation polynomials can be efficiently handled by creating an index array ji,ν . This array will indicate to which gj , the local Hermite interpolation polynomial with index ν in the i − th sub-interval, belongs. Due to the boundary condition given in Eq. (4.9), we do not use the function which equals unity at the left boundary y = −x, instead we use only 2(n + 1) − 1 = 2n + 1 functions. Therefore where cj are coefficients and gj are functions. Orthogonal collocation is then effected by requiring Eq. (4.8) to be satisfied at N = 2n + 1 collocation points yi, chosen to be the Gauss-Legendre integration points of order 2 for the sub-intervals of the grid and yN = x.
A solution of the quantum inverse scattering problem is generally a two step process. The first step handles the calculation of the B−function from the reflection coefficient by employing the Fast Fourier transform. The second step is the actual inversion pro-cess which uses the B−function as the input kernel to the Marchenko integral equation. The results are the reconstructed potentials. This is illustrated in the following discus-sion. For a numerical treatment, the Marchenko integral equation Eq. (4.8), has to be discretized such that the following linear system of equations are obtained In matrix form this set of equations may be written as where the matrices U and K and the vector b are respectively defined by the following equations The integral required in Eq. (5.9) is evaluated via repeated Gauss-Legendre integration over the sub-intervals of our discretization. We solve for the coefficient vector c in Eq. (5.7) via LU-decomposition using the subroutines from LAPACK [63]. To achieve our goal, that is, to obtain the potential, we need to take a derivative according to Eq. (4.10). Thus, another numerical procedure is required to evaluate this derivative. To this end we calculate K(x, x) on the grid points suitable for applying interpolation via Chebyshev polynomials. Those grid points are with respect to the interval [xmin, xmax] given by xi = (xmin + xmax) + cos i − 0.5 π (xmax − xmin) , (5.11) where nch is the number of Chebyshev grid points used. In our case xmin = 0 and xmax = a where a is the range of the potential. Once the coefficients of this interpola-tion are known, we can easily obtain the required derivatives since the coefficients of the Chebyshev interpolation of the derivative can be recursively calculated from the interpolation coefficients for K(x, x).
Coupled channel inversion
The Marchenko fundamental equation Eq. (4.8), can be generalized to a system of coupled channel integral equations [49] where α, β and γ are channel indices of the matrices. It is therefore observed that the input and output kernels of the Marchenko integral equation have been replaced by the matrix kernels as seen in Eq. (5.12). The input kernel Bαβ (x + y) and the output kernel Kαβ (x, y) are nc × nc matrices. Thus, the potential matrix becomes in analogy to Eq. (4.10). The symmetrical input kernel of Eq. (4.6) now has the form where Rαβ (k) is the reflection coefficient matrix. The scattering data Rαβ (k) determine the input kernel Bαβ which are needed to solve the Marchenko integral matrix equation. For a constant value of x at a collocation point yi, Eq. (5.12) can be written as The Marchenko integral matrix equation should be satisfied at N points. Thus, the kernel can then be expanded in terms of gj,x where j numbers the splines and x denotes the dependence of cj and gj on it. Inserting this expansion in Eq. (5.15) and requiring the equation to be fulfilled for all yi we obtain where x has been dropped for the sake of simplicity and δβγ is the Kronecker delta.
Reconstruction of impedance from synthetic data
Knowing η(0) in Eq. (4.12) the seismic impedance profile η(ξ), for ξ > 0 can be recovered from the solution K(ξ, τ ) of the Marchenko integral equation. In order to reconstruct the seismic impedance profile from the output kernel K(ξ, τ ) in the Marchenko integral equation, the expansion in terms of N localized functions can be considered.
Blind deconvolution method
Deconvolution process
Before discussing the deconvolution process we present a brief description of the con-volution model. This model [41] can be described schematically as measured output = output + noise = wavelet ∗ x + noise , where x is the reflectivity sequence or reflector series. Mathematically, it can be written as where z is an observed seismic trace of length N + M − 1, h represents the seismic wavelet of length N , x stands for the white reflectivity sequence of the medium of length M and n is a zero-mean white noise of Gaussian type. The noise sequence is characterized by its variance σ2 [41]. Eq. (5.33) can be written in a convolutional form Our objective is to separate the reflectivity sequence and seismic wavelet from each other by applying the blind deconvolution procedure.
In the literature the system’s unit response is called the reflectivity sequence. In our model it will also include multiple reflections (only a finite number is needed) effected by the system provided the seismic wavelet is shorter than the travel-time distance between the consecutive interfaces. For our numerical computations we identify the reflectivity sequence up to a scaling factor with the unit response of the medium and call it the B−function.
Deconvolution of the seismic reflection data series z when the source wavelet h is known, is a well-understood problem. However, in some investigations such as ours, only the marine seismic reflection data have been provided from which both reflectivity sequence and seismic wavelet have to be retrieved. For this purpose, we apply the blind deconvolution method to estimate the reflectivity sequence and seismic wavelet.
We assume the sea floor to consist of several homogeneous layers that are separated by interfaces. Such an assumption makes it possible to express the reflectivity sequence in terms of a Bernoulli-Gaussian (BG) sequence [64, 65]. Thus, the reflectivity se-quence that defines the generalized BG sequence can be modeled by using two random sequences in the form [41]
where r = (rk ) denotes a zero-mean Gaussian white sequence with variance σr2 and
q = (qk ) stands for the Bernoulli sequence with the probability parameter λ being equal to its mean value [66]. For the probabilities associated with this sequence we have that is, the random variable qk is one with probability λ and zero with probability 1 − λ. The probability of the whole sequence q reads where n is the number of ones in the sequence.
Markov Chain Monte Carlo method
The Markov Chain Monte Carlo (MCMC) method is concerned with the simulation via a Markov chain. It is a Monte Carlo integration procedure in which the random samples are generated by evolving a Markov chain. We are concerned with the MCMC method in a pure Bayesian approach. Upon using the Bayes’ rule, we can write the probability distribution as [44] where θ stands for all parameters of the problem [67], that is, and h, x, λ, and σ2 have the same meaning as above. P (θ|z) is the posterior probability of the model conditional on the observed data z, P (θ), and P (z) describe the prior knowledge and data respectively while P (z|θ) describes the discrepancy between the model and observation.
The complete joint probability distribution is expressed in the form since P (z) is in this case a normalizing constant. The MCMC algorithm is iterative and may require a number of iterations before the first sample from the correct distribution can be provided. These initial iterations are called burn-in iterations and should not be used in the statistical analysis. Thus, the estimation of the reflectivity sequence and seismic wavelet is determined by the simulation of random variables via the MCMC algorithms based on the Gibbs sampler [46], which is regarded as the best known and most popular of the MCMC algorithms [47]. It is an algorithm in which the vector θ (k+1) is obtained from θ(k) by updating the vector elements one at a time. The prior distribution in Eq. (5.43) can be written as where P (h) = P (σ2) = P (λ) = 1 are flat prior probability densities and P (x|λ) is described by a Bernoulli-Gaussian distribution. Upon using Eq. (5.44), we can write Eq. (5.41) in the form For our purpose, we need to calculate the distributions P (h|z, x, σ2, λ), P (x|z, h, λ, σ2), P (σ2|z, x, h, λ), and P (λ|z, x, h, σ2). Thus, Eq. (5.45) can be employed to handle the relevant re-sampling processes.
(a) Re-sampling of the amplitude of the reflectivity sequence
The reflectivity sequence x contains information about the Earth’s structure. In order to statistically separate it from the seismic wavelet we use the BG white sequence model [65]
where N is a Gaussian distribution with specified mean and variance and where λ ∈ (0, 1) is the probability that xm = 1, and both λ and ̺2 are unknown.
It can be shown [65] that the posterior probabilities involving single components of the vector x remain Gaussian mixtures with a structure comparable to that of Eq. (5.48). More precisely, where x(m) is identical to x except for x(mm) = 0. Using Eqs. (5.49)-(5.52), the compo-nents xm of the reflectivity sequence can be re-sampled, one at a time.
(b) Re-sampling of the seismic wavelet
In order to re-sample the seismic wavelet h, we deduce from the Bayes’ rule that P (h|σ2, x, z) ∝ P (z|σ2, x, h), given P (h) = 1, where P (z|σ2, x, h) is given by Eq.
This allows us to conclude that the posterior probability of h is a multivariate Gaussian with mean vector and with covariance matrix R. The latter probability is easy to sample according to h = + Qε, where ε is a normalized Gaussian white noise and QT is a square root matrix of R (that is, such that R = QQT ), such as the one resulting from the Cholesky decomposition.
Contents
1 Introduction
2 Quantum scattering in 1-D
2.1 Single channel Schr¨odinger equation
2.2 Coupled channel problem
3 Seismic scattering in layered media
3.1 Two-layered medium
3.2 Three-layered medium
3.3 Multi-layered medium
3.4 A simulated reflected seismic signal
3.5 Seismic wave equation
3.6 Transformation of seismic wave equation
4 Marchenko inversion method
4.1 Marchenko integral equation
4.2 Obtaining impedance from the Marchenko kernel
5 Computational analysis
5.1 Single channel inversion
5.2 Coupled channel inversion
5.3 Reconstruction of impedance from synthetic data
5.4 Blind deconvolution method
6 Application I: Quantum inversion
6.1 Single channel inversion
6.2 Coupled channel inversion
7 Application II: Seismic inversion
7.1 Slab-type seismic impedance
8 Application III: Marine data
8.1 Estimation of seismic impedance
8.2 Scaling factor
9 Conclusions and directions for future work
9.1 Single and coupled channel inversions
9.2 Synthetic seismic reflection data
9.3 Marine seismic reflection data
9.4 Directions for future work
GET THE COMPLETE PROJECT
