- Email: [email protected]
Dynamic conversion of head waves in seismic data processing
Russian Geology and Geophysics 49 (2008) 780–789 www.elsevier.com/locate/rgg
Dynamic conversion of head waves in seismic data processing A.F. Emanov *, B.C. Seleznev, N.A. Korshik Geophysical Surveys of the Siberian Branch of the RAS, Altai-Sayan Department, 3 prosp. Akad. Koptyuga, Novosibirsk, 630090, Russia Received 21 August 2007; accepted 26 December 2007
Abstract We suggest a new technique for processing head-wave data with a multistage Wiener filter, which provides conversion and stacking of head waves from a receiver array into a user-specified point. Different conversion algorithms are compared and shown to be a combination of parallel and serial Wiener filters. The Wiener filter accuracy for a conversion iteration depends on the in-line fold and on the coherence spectrum. The latter describes the signal/noise ratio dependence where head waves are the signal and the noise includes all other components and microseisms. Equations for the accuracy of head-wave traces have been derived from a corresponding equation for a single conversion iteration using the multistage Wiener filter schemes and proceeding from the theory of coupling stochastic systems. The obtained equations are applicable to different field data processing algorithms. The discussed conversion algorithms differ in noise regularization and are fit for data of any quality. © 2008, IGM, Siberian Branch of the RAS. Published by Elsevier B.V. All rights reserved. Keywords: Head waves; coherent signals; Wiener filter; seismic exploration
Introduction
Conversion with a single Wiener filter
Processing algorithms for head-wave data are much less developed than for high-fold reflection profiling data. An intention to apply reflection-data processing techniques to head waves gave rise to a common depth surface (CDS) method (Epinat’eva et al., 1990) in which, however, the dynamic parameters of head waves are distorted, and there are some other drawbacks. New algorithms for head-wave data processing became possible on the basis of the dynamic spectral ratio for head waves recorded at four correlated points of a measurement system (Krylov and Sergeev, 1985). The conversion algorithms we suggest combine Wiener filters and dynamic conversion of head waves (Emanov et al., 1998; Seleznev and Emanov, 1998). Head-wave data are processed in three steps: constraining refractor coverage, compressing redundant data within each coverage in the time section, and inversion for subsurface parameters. The suggested conversion algorithms were designed for the case of already constrained coverage, and in this paper we focus on the second processing step, that of picking head waves from the measured wavefield.
Krylov and Sergeev (1985) derived equations for headwave fields measured at four correlated points (Fig. 1): Fi, j(ω)Fi + 1, j + 1(ω) = Fi + 1, j(ω)Fi, j + 1(ω),
(1)
where F(ω) are the spectra of head waves at stacking points, and i and j are the numbers of sources and receivers, respectively. According to (1), the wavefield can be converted from one stacking point to another using a linear system. The conversion from the point (i, j) into the point (i + 1, j) is governed by the spectral equation ^ ^ Fi + 1, j(ω) = Fi, j(ω)hii, i + 1(ω),
(2)
where hii, i + 1(ω) is the frequency response of the conversion ^ linear filter, Fi, j(ω) is the spectrum of the signal plus random ^ noise at the point (i, j), and Fi + 1, j(ω) is the spectrum of the seismogram converted into (i + 1, j). For head waves, using (1), one can write the condition of equality of filters for multiple correlations (Fig. 1)
* Corresponding author.
E-mail address: [email protected] (A.F. Emanov) hii, i + 1(ω) = hii,+i +1 1(ω) = hii,+i +2 1(ω) 1068-7971/$ - see front matter D 2008, IGM, Siberian Branch of the RAS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.rgg.2007.12.015
+ … + hii,+i +k 1(ω).
(3)
781
A.F. Emanov et al. / Russian Geology and Geophysics 49 (2008) 780–789
averaged autocorrelation function at the input of the linear system, we obtain R12(θ) − hopt(τ)∗R11(τ − θ) = 0
(5)
or, in the frequency domain, hopt(ω) =
Fig. 1. Observation system and ray scheme. 1 — i-th source; 2 — j-th receiver; 3 — point into which head waves are converted; 4 — point from which head waves are converted; 5 — points used for filter construction. R is refractor.
This equation is valid in the absence of noise in seismograms and implies that conversion from (i, j) into (i + 1, j) is possible with a filter obtained for any correlation. In the presence of noise, when the filters for different correlations are dissimilar, one can apply multiple conversion of the same signal with different filters and its equilibrium (Krylov and Sergeev, 1985) or weighted (Mitrofanov and Sergeev, 1986) stacking. There is an alternative way of conversion, assuming, according to (3), that a set of linear systems is actually a set of multiple input and output signal and noise realizations of the same linear system. Thus, we have a linear system with a known set of input-output signal pairs recorded on the background of random noise and a single realization of an input head-wave record, which is to be converted into an output record. We search an optimal conversion filter using the approach applicable to single Wiener filters (Bendat and Piersol, 1980; Gol’din, 1974; Seleznev and Emanov, 1998). Let the input signal of the frequency h(t) be FI, j(t) and the output signal be Fi, j+1(t) = Fi, j(t)*h(t), where (*) denotes convolution. One ^ can measure only Fi, j(t) = Fi, j(t) + Wi, j(t), where Wi,j(t) is the random noise realization. The optimal Wiener filter to convert ^ the input signal Fi, j(t) into the output signal Fi, j+1(t) is found by minimizing, along h(t), the mathematical expectation of the ^ square misfit between the filtered Fi, j(t) signal and the true output signal Fi, j+1(t): ^ (4) M |Fi, j + 1(t) − Fi, j (t)∗hopt(t)2| = min. Zeroing the derivative of (4) along h(t), after simple transformation, gives ^ M [Fi, j + 1(t)Fi, j(t − θ)] − ^ ^ hopt(τ)∗M [Fi, j (t − τ)Fi, j (t − θ)] = 0. ^ Taking into account that M [Fi, j + 1(t)Fi, j (t − θ)] = R12(θ) is the averaged cross-correlation function between the input ^ ^ and output signals and M [Fi, j (t)Fi, j (t − θ)] = R11(τ − θ) is the
^∗ M [Fi, j + 1(ω)Fi, j(ω)] = , ^ 2 R11(ω) M |Fi, j(ω)|
R12(ω)
(6)
where (*) is the spectrum complex conjugation. The main problem is how to calculate the average reciprocal and auto spectra in (6). We have a set of input and output signal realizations against random noise realizations. The spectra can be averaged over the realization set using records along parallel correlations. Then, ^ ^ R12(ω) = ∑ Fi + k, j + 1(ω)F∗i + k, j(ω), k
^ R11(ω) = ∑ |Fi + k, j (ω)|2, k
and the equation for the head-wave conversion from (i, j) into (i, j + 1) becomes ^
^
^ Fi, j + 1(ω) = Fi, j (ω)
∑ [Fi + k, j + 1(ω)F∗i + k, j(ω)] k
.
^
∑ |Fi + k, j (ω)|2
(7)
k
A conversion iteration in case (1) is from three seismograms while each iteration in (7) uses a number of seismogram pairs from parallel correlations. In case (1), for a measurement system as in Fig. 1, the same wave is converted and stacked in many iterations, while (7) requires a single conversion iteration with a more accurate filter based on many seismograms. It is easy to show that ∗
Mhopt(ω) =
∗
M {[Fi, j + 1(ω) + Wi, j + 1(ω)] [Fi, j(ω) + Wi, j(ω)]} 2
M |Fi, j (ω) + Wi, j (ω)|
=
h (ω) , N (ω) 1 + __ 2 |F (ω)|
__ where N(ω) is the noise power spectrum and N (ω)/|F(ω)|2 is the square noise/signal ratio averaged over the set of seismograms. The mathematical expectation of the optimal Wiener filter for conversion of head waves differs from the true filter in its regularization with respect to the average noise/signal ratio. The conversion procedure for an observation system of the fold n is as shown in Fig. 2. Conversion into O is required only for seismograms from the correlations that cross at this point (Seleznev and Emanov, 1998). Assuming that each signal from the selected correlations is converted into the point O with a single filter, we obtain a single conversion option, using the seismograms from the triangle only (Fig. 2, A). This scheme has a drawback that ever fewer parallel correlations are used for obtaining the filters as the distance between the starting and target points increases. As a result, the mathematical expectations of the reciprocal and auto spectra are estimated from different realization numbers, and the filters differ in accuracy.
782
A.F. Emanov et al. / Russian Geology and Geophysics 49 (2008) 780–789
Fig. 2. Schemes for conversion of head wavefield into a point. — triangle; B — parallelogram; C — trapezium. 1 — point into which wave filed is converted; 2 — point from which wave field is converted; 3 — seismogram points used for filter construction.
It is possible to provide a uniform accuracy of the conversion filters by using a parallelogram as in (Fig. 2, B) for a correlation instead of the single-step conversion. The seismogram from the top point is converted into the point next to the filter calculated from (2n – 3) correlations; then the sum of two seismograms is converted a step further and is summed with the third one, and thus the whole wavefield becomes stacked successively at a given point. The number of averaging iterations is (2n – 3) for any step. The back correlation points are likewise converted using a parallelogram as in Fig. 2, B, but with an opposite tilt, and the conversion of two correlations together covers a trapezium-shaped area (Fig. 2, C).
Conversion from two stacking points into one We convert signals into the point O(i, j) (Fig. 2, C) from two arbitrary points which cross at O, one of a forward (i, j + b) and the other of a back (i – a, j) correlation. Conversion using a two-stage Wiener filter (Emanov et al., 2001) follows ^ ^ F′i, j (ω) = Fi − a, j (ω)Li,i −j a, j(ω) + Fi, j + b (ω)Li,i, jj + b (ω),
(8)
where Li,i −j a, j (ω) and Li,i, jj + b (ω) are the frequency responses of the Wiener filter for conversion from two points into one; prime discriminates between the measured and computed records. The responses of the two-stage Wiener filter are found using a standard approach, by minimizing the mathematical expectation of the square misfit
^ M |Fi, j (ω) − Li,i −j a, j (ω)Fi − a, j (ω) − ^ Li,i, jj + b (ω)Fi, j + b (ω)|2 = min.
(9)
Differentiating along each filter and zeroing the derivatives gives j i, j i, j + b i, j Rii −− a, a, j (ω)Li − a, j (ω) + Ri − a, j (ω)Li, j + b (ω) −
Ri,i −j a, j (ω) = 0,
(10)
Ri,i, jj ++ bb (ω)Li,i, jj + b (ω) + Rii,−j +a,bj (ω)Li,i −j a, j (ω) − Ri,i, jj + b (ω) = 0, ^ 2 j where Rii −− a, a, j (ω) = M |Fi − a, j (ω)| is the mathematical expectation of the trace square spectrum at the stacking point (i – a, ^ j); Ri,i, jj ++ bb (ω) = M |Fi, j + b (ω)|2 is the same at (i, j + b); ^ ^ Ri,i −j +a,bj (mω) = Rii,−j +a,bj (ω) = M [Fi − a, j (ω)F∗i, j + b (ω)] is the reciprocal spectrum between traces at (i – a, j) and (i, j + b) averaged over the set of realizations; Ri,i −j a, j (ω) and Ri,i, jj + b (ω) are the mathematical expectations of the reciprocal spectra between the records at (i, j) and (i – a, j) and at (i, j) and (i, j + b); Ri,i −j a, j (ω) and Ri,i, jj + b (ω) are the same for the pairs (i – a, j) and (i, j) and (i, j + b) and (i, j), respectively. The frequency responses of the two-stage Wiener filter are
Li,i −j a, j (ω) =
i, j + b i, j Ri − a, j (ω)Ri, j + b (ω) i, j Ri − a, j (ω) 1 − i, j i, j + b Ri − a, j (ω)Ri, j + b (ω) i − a, j
Ri − a, j (ω) 1 − γab (ω) 2
,
783
A.F. Emanov et al. / Russian Geology and Geophysics 49 (2008) 780–789
Li,i, jj + b (ω) =
i − a, j i, j Ri, j + b (ω)Ri − a, j (ω) i, j Ri, j + b (ω) 1 − i, j i − a, j Ri, j + b (ω)Ri − a, j (ω) 2 i, j + b Ri, j + b (ω) 1 − γab(ω)
^
∑ Li,i, jj + b (ω)Fi, j + b (ω)|2 = min. ,
(11)
2
where γ2ab(ω) =
Ri, j + b (ω) i − a, j i − a, j i, j + b Ri − a, j (ω)Ri, j + b (ω)
is the square coherence spec-
trum between the traces at (i – a, j) and (i, j + b). These are the equations for a standard two-stage Wiener filter. The mathematical expectations of the reciprocal and auto spectra in (11) are found using (3), i.e., the spectra are averaged over parallel correlations. Then, ^ ^ Ri,i −j a, j(ω) = ∑ Fi − a + k, j (ω)Fi∗+ k, j(ω), k
where the subscript k marks passage to parallel correlations and is either positive or negative, but never zero, ^ 2 j Rii −− a, a, j (ω) = ∑ |Fi − a + k, j (ω)| , k
Ri,i, jj ++ bb (ω)
^ = ∑ |Fi, j + b + k (ω)|2, k
^ ^ Ri,i, jj + b (ω) = ∑ Fi, j + b + k (ω)Fi,∗j + k (ω), k
^ ^ Ri,i −j +a,bj (ω) = ∑ Fi − a + k, j (ω)Fi,∗j + b + k (ω). k
The spectral responses of the two-stage Wiener filter include single Wiener filters and a frequency-dependent term as a multiplier which provides weighing in the stacking of converted seismograms.
A multistage conversion Wiener filter Judging by the above result, we may expect to improve the conversion scheme by employing a multistage Wiener filter to convert head waves simultaneously from all points of the forward and back correlations into a stacking point, using the equation (Gol’din, 1974) ^ F′i, j (ω) = ∑ Fi − a, j (ω)Li,i −j a, j (ω) + a
^
∑ Fi, j + b (ω)Li,i, jj + b (ω),
(12)
b
where a and b are nonzero integers which provide passage from one seismogram to another along the correlations. Conversion of waves from a given set of points into a single point requires an algorithm for a multistage Wiener filter. First we minimize the mathematical expectation of the square misfit: ^ M |Fi, j (ω) − ∑ Li,i −j a, j (ω)Fi − a, j (ω) − a
(13)
b
Zeroing the derivatives along each filter gives a system of equations written as the matrix ~ − 1), j Ri, j (ω) Ri − a, j(ω) (ω) Rii −− (a a, j ~ i − a, j i − a, j i − a, j i − (a − 1), j Ri, j ~ i − (a − 1), j(ω) Ri − (a − 1), j(ω) Ri − (a − 1), j(ω) ⋅⋅ i − (a − 1), j Ri, j i − a, j i − (a − 2), j(ω) Ri − (a − 2), j(ω) Ri − (a − 2), j(ω) ⋅⋅ ⋅ ⋅ ⋅ ⋅⋅ ⋅⋅ ⋅ ~ ⋅ ⋅ ⋅⋅ i, j i − a, j − 1), j (ω) Rii −+ (a Ri + 1, j(ω) Ri + 1, j(ω) ⋅⋅ 1, j ~ i, j = i − a, j i − (a − 1), j ⋅⋅ (ω) (ω) (ω) R R R i, j + b i, j + b i, j + b ~ i, j i − a, j ⋅⋅ i − (a − 1), j Ri, j + (b − 1)(ω) Ri, j + (b − 1)(ω) Ri, j + (b − 1)(ω) ⋅⋅ ~ i, j i − a, j ⋅⋅ Ri, j + (b − 2)(ω) Ri, j + (b − 2)(ω) Rii,−j +(a(b− −1),2)j(ω) ⋅⋅ ⋅ ⋅ ⋅ ⋅⋅ ⋅ ⋅ ⋅ ~ Ri, j (ω) Ri − a, j(ω) Ri, j + 1i − (a − 1), j(ω) i, j + 1 i, j + 1 j Rii +− 1, a, j(ω)
Ri,i −j +a,bj(ω)
Rii,+j +1,bj(ω)
Ri,i, jj ++ bb(ω)
j i, j + b Rii +− 1, (a − 1), j(ω) Ri − (a − 1), j(ω) j i, j + b Rii +− 1, (a − 2), j(ω) Ri − (a − 2), j(ω) ⋅ ⋅ ⋅ ⋅ i + 1, j i, j + b Ri + 1, j(ω) Ri + 1, j(ω) j i, j + b Rii,+j +1,(b − 1)(ω) Ri, j + (b − 1)(ω) j i, j + b Rii,+j +1,(b − 2)(ω) Ri, j + (b − 2)(ω) ⋅ ⋅ ⋅ ⋅ Rii,+j +1,1j(ω) Ri,i, jj ++ b1(ω)
Ri,i −j +a,(bj − 1)(ω) ⋅ ⋅ ⋅⋅ Ri,i −j +(a(b− −1),1)j(ω) ⋅⋅ Ri,i −j +(a(b− −2),1)j(ω) ⋅ ⋅ ⋅⋅ ⋅ ⋅⋅ ⋅ i, j + (b − 1) (ω) ⋅ ⋅ Ri + 1, j ⋅⋅ − 1) (ω) ⋅ ⋅ Ri,i, jj ++ (b b − 1) ⋅⋅ Ri,i, jj ++ (b (b − 1)(ω) ⋅⋅ i, j + (b − 1) Ri, j + (b − 2)(ω) ⋅ ⋅ ⋅ ⋅⋅ ⋅ ⋅⋅ Ri,i, jj ++ 1(b − 1)(ω) ⋅ ⋅
Li, j i − 1, j i, j L i − (a − 1), j Li, j i − (a − 2), j ⋅ ⋅ ⋅ ⋅ i, j i, j + 1 Ri + 1, j(ω) Li + 1, j (14) × i, j , i, j + 1 Ri, j + b(ω) Li, j + b i, j Ri,i, jj ++ 1(b − 1)(ω) Li, j + (b − 1) i, j Li, j + (b − 2) Ri,i, jj ++ 1(b − 2)(ω) ⋅ ⋅ ⋅ ⋅ i, j + 1 Ri, j + 1(ω) i, j L 1 + i, j ^ i − a, j where Ri − a, j(ω) = M |Fi − a, j (ω)|2 is the mathematical expectation of the trace square spectrum at the stacking point (i – a, Ri,i −j +a,1j(ω) Ri,i −j +(a1− 1), j(ω) Ri,i −j +(a1− 2), j(ω)
784
A.F. Emanov et al. / Russian Geology and Geophysics 49 (2008) 780–789
^ j); Ri,i, jj ++ bb (ω) = M |Fi, j + b (ω)|2 is the same at (i, j + b); ^ ^ Ri,i −j +a,bj (ω) = Rii,−j +a,bj (ω) = M [Fi − a, j (ω)F∗i, j + b (ω)] is the reciprocal spectrum between traces at (i – a, j) and (i, j + b) averaged over the set of realizations; ^∗ ^ j i − a, j Rii +− 1, a, j (ω) = Ri + 1, j (ω) = M [Fi − a, j (ω)Fi + 1, j (ω)] is the auto spectrum between traces at (i – a, j) and (i + 1, j) averaged over the set of realizations; − 1) (ω) = Ri,i, jj ++ b(b − 1) (ω) = Ri,i, jj ++ (b b ^ ^ M [Fi, j + b (ω)F∗i, j + (b − 1) (ω)]
is the auto spectrum between traces at (i, j + b) and (i, j + (b – 1)) averaged over the set of realizations; ~ ~ ^ ^ Ri,i −j (a − 1), j (ω) = Rii,−j (a − 1), j (ω) = M [Fi − (a − 1), j (ω)F∗i, j (ω)] ~ and Ri,i, jj + b (ω) are the mathematical expectations of the reciprocal spectra between the pair (i – (a – 1), j) and (i, j) and the pair (i, j + b) and (I, j); they are the same for the other points. The mathematical expectations in auto and reciprocal spectra are substituted by averaging over parallel correlations. This gives a matrix for which the solution can be obtained from the column vector of the spectra of the multistage Wiener filter. Solving (14) and converting head waves from many points into one point with a multistage Wiener filter is a standard problem. The multistage Wiener filter is coupled with dynamic conversion of head waves at the step of averaging of the auto and reciprocal spectra (matrix elements) over parallel correlations.
The accuracy of head-wave conversion into a time section The idea of the accuracy of a head-wave seismic section stems from the meaning of the latter. In the reported processing technique, a head-wave time section for a given offset is meant as a set of traces containing head waves only. The algorithms are actually a tool for selecting parallel refractions, which are the head waves that arrive at a receiver. Any head-wave section is obviously approximate. Below we discuss how the accuracy depends on the in-line fold and on the wavefield features. The accuracy is estimated for a whole trace consisting of several samples within a given time window rather than for some individual parameter. The accuracy of a trace is meant as a difference between the computed and true values, the true value being a time function describing head waves at a given receiver. The accuracy of a filter obtained according to (7) is estimated as follows. The uncertainty in (7) is found as (Bendat and Piersol, 1980) M {[h′0i(ω) − h0i(ω)]2} = M {[h′0i(ω) − M [h′0i(ω)]]2} + M {[M [h′0i(ω)] − h0i(ω)]2},
where h′0i is the computed response, the first term is the scatter D[h′0i(ω)] and the second term is the square misfit b2[h′0i(ω)]. As it was shown before, the misfit leads to regularization of filter’s amplitude response, which improves the solution stability, while the uncertainty in filter construction is defined by the scatter. To investigate the random error in the conversion filter obtained by (7), we consider separately the phase and amplitude errors. According to (7), the phase response of the conversion filter is found as its mathematical expectation for the reciprocal spectrum between the filter input and output. The reciprocal spectrum consists of a real and an imaginary parts: |R′0i(ω)|2 = C20i (ω) + Q20i (ω), and tg θ0i (ω) = C0i (ω)/Q0i (ω),
(15)
where C0i(ω) is the real part, Q0i(ω) is the imaginary part, and θ0i(ω) is the phase spectrum. The formula for random error in the filter’s phase response is derived using the mathematical expectations of the square real and imaginary parts of the reciprocal spectrum 1 M [C20i (ω)] = [R0(ω)Ri (ω) + 3C20i (ω) − Q20i (ω)], 2
(16)
1 M [Q20i (ω)] = [R0(ω)Ri (ω) + 3Q20i (ω) − C20i (ω)].
(17)
2
Equations (16) and (17) are easy to derive if the real and imaginary parts of the reciprocal spectrum are expressed via the real and imaginary spectrum parts at the filter input and output, squared and defined as mathematical expectations. The scatter in the real part of the spectrum is D [C0i (ω)] = M [C20i (ω)] − C20i (ω) = 1 2
[R0(ω)Ri (ω) + C20i (ω) − Q20i (ω)],
(18)
and that for the imaginary part is D [Q0i (ω)] = M [Q20i (ω)] − Q20i (ω) = 1 2
[R0(ω)Ri (ω) + Q20i (ω) − C20i (ω)].
(19)
The expression for the covariance between the real and imaginary parts of the reciprocal spectrum can be obtained in the same way: cov[C′0i (ω), Q′0i (ω)] = M [C′0i (ω)Q′0i (ω)] − C0i (ω)Q0i (ω) = C0i (ω)Q0i (ω).
(20)
Differentiation along both terms of (15) and finite-difference approximation gives sec2θ0i (ω)⋅∆θ0i (ω) ≈
C0i (ω)∆Q0i (ω) _ Q0i (ω)∆C0i (ω) 2
C0i (ω)
.
Taking into account that sec2θ0i (ω) = |R0i (ω)|2/C20i (ω), we obtain
A.F. Emanov et al. / Russian Geology and Geophysics 49 (2008) 780–789
∆θ0i (ω) ≈
C0i (ω)∆Q0i (ω) − Q0i (ω)∆C0i (ω) 2
|R0i (ω)|
.
(21)
Thus, the random error scatter in the filter phase response, without averaging, is D [θ′0i (ω)] = M [(∆θ0i (ω))2] ≈ 1 4
M [(C0i (ω)∆Q0i (ω) − Q0i (ω)∆C0i (ω))2] =
4
{C20i (ω)D [Q′0i (ω)] −
|R0i (ω)| 1 |R0i (ω)|
2C0i (ω)Q0i (ω) cov [C′0i (ω), Q′0i (ω)] + Q20i (ω)D [C′0i (ω)]}. Substitution of (18), (19), and (20) into (21) gives D [θ′0i (ω)] ≈ 2 (1 − γ0i (ω)) 2 2γ0i (ω)
4 |R0i (ω)|2 |R0i (ω)| (ω) − (ω)R R = 0 i 4 2 2 |R0i (ω)|
1
,
(22)
where 2
γ (ω) =
|R0i (ω)|
r00(ω)Rii (ω)
,
(23)
where R00(ω) is the input spectrum of the autocorrelation function of traces and Rii(ω) is the same at the filter output. Using model (1) gives R0i (ω) = h0i (ω)M |F0(ω)|2, R00(ω) = M |F0(ω)|2 + M |W0(ω)|2, 2
2
2
Rii (ω) = |h0i (ω)| M |F0(ω)| + M |Wi (ω)| . These equations are valid assuming the independence of input and output noise realizations and the absence of their correlation with the signal. Substituting these equations into (23) gives γ2(ω) =
1 , (1 + α0(ω)) (1 + αi (ω)) 2
where α0(ω) =
M |W0(ω)|
2
M |F0(ω)|
(24) 2
, αi (ω) =
M |Wi (ω)|
2
M |Fi (ω)|
are the input/out-
put ratios of weighted mean square noise/signal ratios. Taking into account that the phase spectra are averaged over (2n – 3) independent records, the random error scatter may be expected to become (2n – 3) times as low. Therefore, the rms error in the filter phase response is σθ(ω) ≈
The conversion filter is constructed using a set of input and output signal and noise records. Obtaining the explicit formula for the amplitude error is hard for a model implying noise both at the input and at the output but is easy (Bendat and Piersol, 1980) for the case when the input noise is vanishing relative to the signal. Inasmuch as only the γ20i(ω) value matters for the filter uncertainty and it is of no importance whether input or output noise is involved, the corresponding expression for this specific case will be valid in the general case as well. The records of head waves from the same refractors at two different points are related as linear equation system (2). In fact, we have the records of head waves mixed with the noise ^ Fi, j (t) = Fi, j (t) + Wi, j (t), (26) Fi + 1, j (t) = Fi + 1, j (t) + Wi + 1, j (t) at the input and output of the system. Equations (2) and (26) constitute a model of signals which lies at the base of the conversion algorithms. As a consequence of (2), in linear systems of this kind, Ri + 1, j (ω) = |h (ω)|2Ri, j (ω),
2
2 1 − γ0i(ω) √
|γ0i(ω)| √ 2(2n − 3)
.
(25)
Using (25) we can find the frequency-dependent phase response in radians. The same estimation for amplitude is more complicated. One might calculate the confidence intervals of the filter amplitude response at each conversion iteration (Bendat and Piersol, 1980) but having an explicit formula is much more convenient.
785
(27)
where Ri + 1, j (ω) = M |Fi + 1, j (ω)|2 and Ri, j (ω) = M |Fi, j (ω)|2 are the auto spectra of head waves at the system input and output. Calculation of finite differences for both terms of (27) gives ∆Ri + 1, j (ω) ≈ |h (ω)|2∆Ri, j (ω) + Ri, j (ω)∆|h (ω)|2.
(28)
Having transformed (28) into Ri, j (ω)∆|h (ω)|2 ≈ ∆Ri + 1, j (ω) − |h (ω)|2∆Ri, j (ω) and having found the mathematical expectation for both terms of this equation, we obtain R2i, j (ω)D [|h (ω)|2] ≈ D [Ri + 1, j (ω)] − 2|h (ω)|2 × cov [Ri, j (ω), Ri + 1, j (ω)] + |h (ω)|4D [Ri, j (ω)].
(29)
In the further consideration we assume that noise is vanishing relative to the signal at the point (i, j) but increases correspondingly at (I + 1, j). This assumption allows obtaining per se input and output amplitude spectra of head waves from measured data. Assuming vanishing noise at (i, j), we apply the Fourier transform to find the head wave spectrum and then obtain their spectrum at the filter output with reference to the coherence spectrum: ^ Ri + 1, j (ω) = γ2(ω)Ri + 1, j (ω), (30) ^ 2 (ω) = − γ (ω), (ω)]R [1 RW i + 1, j i + 1, j ^ ^ where Ri + 1, j (ω) = M |Fi + 1, j (ω)|2 is the auto spectrum of the seismogram at (i + 1, j) and RW i + 1, j (ω) is the noise auto spectrum at the same point. Thus, with the above assumptions, one can identify the input and output amplitude spectra of head waves and the amplitude spectrum of noise.
786
A.F. Emanov et al. / Russian Geology and Geophysics 49 (2008) 780–789
The scatter in the estimate of the input auto spectrum is, by definition, D [Ri, j (ω)] = M [R2i, j (ω)] − R2i, j (ω), and the mathematical expectation of the square auto spectrum of a stochastic process of a duration T will be 2R2i, j (ω) (Bendat and Piersol, 1980); correspondingly, D [Ri, j (ω)] = R2i, j (ω). The random error in the auto spectrum is fairly large and its scatter is independent of the respective duration. The problem is that the resolution of the spectrum increases with the time window while the accuracy remains invariable. In our case, finding the mathematical expectation is replaced by averaging over parallel correlations, whereby the random error scatter diminishes by a factor of (2n – 3), to give D [Ri, j (ω)] = R2i, j (ω)/(2n − 3). ^ As it follows from (26), D [Ri + 1, j (ω)] = D [Ri + 1, j (ω)] + D [RW i + 1, j (ω)], and the random error scatter in the spectrum of head waves at (i + 1, j) is obtained using (30) as 2
D [Ri + 1, j (ω)] =
2
[2 − γ (ω)]Ri + 1, j (ω) 2
γ (ω) (2n − 3)
(31)
and cov [Ri, j (ω), Ri + 1, j (ω)] = M [∆Ri, j (ω)∆Ri + 1, j (ω)] = Ri, j (ω)Ri + 1, j (ω) 2n − 3
.
(32)
Substitution of (30)–(32) into (29) gives D [|h (ω)|2] ≈
2
4
2[1 − γ (ω)] |h (ω)| 2
γ (ω)[2n − 3]
.
(33)
Using the equation D [A2] = 4A2D [A] (Bendat and Piersol, 1980), we obtain
D [|h (ω)|] =
2
2
[1 − γ (ω)]⋅|h (ω)| 2
2γ (ω)⋅[2n − 3]
,
(34)
wherefrom there follows the expression for the relative filter amplitude error 2 1/2 ε[|h (ω)|] = σ[|h (ω)|] = [1 − γ (ω)] .
|h (ω)|
|γ(ω)|⋅√ 2(2n − 3)
(35)
According to (25) and (35), the accuracy of the conversion filter improves as the in-line fold increases. We have discussed the accuracy of a Wiener filter for conversion of head waves from one point to another but actually we concern about the accuracy of the trace obtained by conversion from many points. We showed before (Seleznev and Emanov, 1998) that conversion involves head waves from points of the forward and back correlations that cross at a given point. Then, we deal with an algorithm, which brings together many Wiener filters according to a certain scheme (Pugachev and Sinitsin, 2000), and study their properties as of a single whole proceeding from the properties of the constituent linear systems. Thus, the conversion algorithms in our case correspond to a multistage Wiener filter. The above conversion schemes are actually various schemes of a multistage Wiener filter with different possibilities for solving the same problem. See Fig. 3 for flow charts of several conversion schemes. They are a parallel scheme (Fig. 3, A) in the case of a single conversion (Fig. 2, A) and a serial scheme (Fig. 3, B) for successive conversion of the wavefield into a point (Fig. 2, B) when the computed signal is summed with the measured one after each conversion iteration. Inasmuch as the input signal enters the filter series after each element, the scheme is not purely serial. The scheme C in Fig. 3 is a combination of serial and parallel Wiener filters and is the same as B but summation is at the
Fig. 3. Multistage Wiener filter algorithms of head-wave data processing. — parallel conversion; B — multistage filter corresponding to conversion scheme of Fig. 2, B; C — mixed serial-parallel multistage filter equivalent to case B.
787
A.F. Emanov et al. / Russian Geology and Geophysics 49 (2008) 780–789
end of the complete conversion procedure. There can be other possibilities besides the schemes of Fig. 3 but the latter are good examples because they show the extreme cases. A multistage Wiener filter can be expressed as a single whole via a pulse response. In our case, the pulse response is obtained from measured seismograms and is, because of noise, an approximation to the true response. We write a multistage Wiener filter as a mathematical expectation of the pulse response and its scatter. The response for the parallel scheme (Fig. 3, ) is the sum of frequency responses of the constituent elements (Pugachev and Sinitsin, 2000). Since the frequency responses of the filter elements bear an error, the mathematical expectation of the multistage filter response will be the sum of the mathematical expectations of the individual filters. For the scheme of Fig. 3, A it is easy to show that 2n − 1
2n − 1
k=1
l=1
M [hi,Σ j (ω)] = ∑ M [hi,i −j k, j (ω)] + ∑ M [hi,i, jj − l (ω)].
1+
2 i, j M |Wi − k, j (ω)| 2 i, j M |Fi − k, j (ω)|
and that, according to the latter equation, the mathematical expectation of the phase-amplitude response is independent of noise, we obtain the amplitude response of the multistage filter 2n − 1
(ω)| = ∑
2
k=1
2n − 1
∑ 1 + M |W
2
1 + M |Wi − k, j (ω)| /M |Fi − k, j (ω)|
2
2
(ω)| /M |Fi, j − l (ω)|
2n − 1
2n − 1
k=1
l=1
2n − 1 2 i, j i, j [1 − (γi − k, j (ω)) ] |hi − k, j|
2n − 1 2 i, j i, j [1 − (γi, j − l (ω)) ] |hi, j − l (ω)| 2
i, j
2
i, j
2(γi − k, j (ω)) [2n − 3]
k=1
2(γi, j − l (ω)) [2n − 3]
2n − 1
(ω)| = ∑
k=1
i, j
2 1/2
i, j
|γi − k, j (ω)|⋅√ 2(2n − 3)
l=1
[1 − (γi, j − l (ω)) ]
.
2
i, j
2(γi, j − l (ω)) (2n − 3)
In the case of serial conversion (Fig. 2, B) performed using the equivalent schemes B and C of Fig. 3, we actually deal with a multistage filter of several parallel series of Wiener filters. First we consider conversion with a series of filters. For a series of (2n – 1) filters, the mathematical expectation of the frequency response (hP) is the product of the mathematical expectations of the frequency responses of the elementary filters: 2n − 1
M |hi,P j (ω)| = ∏
k=1
i − k, j
|hi − k, j (ω)| 2
2
1 + M |Wi − k, j (ω)| /M |Fi − k, j (ω)|
.
For a series of Wiener filters, it is easy to obtain the equation for the relative error in the amplitude response as a sum of relative errors of the responses in the product (ω)| = ∑
j ε|hii −− k, k, j
2n − 1
(ω)| = ∑
k=1
i − k, j
2
[1 − (γi − k, j (ω)) ] i − k − 1, j
|γi − k, j
2n − 1
2n − 1
k=1
k=1
D [θi,Σ j (ω)] = ∑ D [θii −− kk,−j 1, j (ω)] = ∑
2(2n − 3) (ω)|⋅√
.
i − k − 1, j
[1 − (γi − k, j i − k − 1, j
|γi − k, j
2
(ω)) ]
2(2n − 3) (ω)|⋅√
.
For a combination of serial and parallel multistage filters (C in Fig. 3), the mathematical expectation of its frequency response becomes
+
l=1
k=1
i, j
i, j
∏
g=1
|γi, j − l (ω)|⋅√ 2(2n − 3)
M [hi,i, jj − l − g (ω)] .
2n − 1
.
p=1
M [hi,i −j k − p, j (ω)] +
The frequency response is obtained in this case assuming that summation applies to the relative errors of elements in the series and to the response scatters of each parallel series. Following the rules from (Pugachev and Sinitsin, 2000), the equation for the scatter of the frequency response errors in the mixed serial-parallel multistage Wiener filter is easily obtained as
2 1/2
[1 − (γi, j − l (ω)) ]
∏
2n − 1 2n − l − 1 l=1
2n − 1
+∑
M [hi,Σ j (ω)] = ∑
∑
.
[1 − (γi − k, j (ω)) ]
+∑e
2
2(γi − k, j (ω)) (2n − 3)
2
i, j
−jω
2n − 1 2n − k − 1
With the latter equation it is easy to pass to the relative error of the multistage filter amplitude response: ε|hi,Σ j
2n − 1
[1 − (γi − k, j (ω)) ] i, j
i, j
+ ∑ e−jωD [θi, j − l (ω)] =
The scatter of the phase response of the series is the sum of the respective scatters for each filter:
.
Now we consider separately the scatter of the absolute frequency response of the multistage filter and the phase response scatter. For the former, with regard to (35),
∑
∑e
k=1
2n − 1 l=1
2
i, j
−jω
k=1
D [hi,Σ j (ω)] = ∑ D [hi,i −j k, j (ω)] + ∑ D [hi,i, jj − l (ω)].
l=1
2n − 1
ε|hi,P j
For the case of parallel conversion, the scatter of the frequency response equals the sum of the respective scatters of the filter elements:
D |hi,Σ j (ω)| = ∑
=∑e
i, j
−jωD [θi − k, j (ω)]
k=1
2n − 1
+
i, j
|hi, j − l(ω)|
i, j − l
l=1
e
i, j
|hi − k, j (ω)|
2n − 1
i, j
−jωD [θSIGMS (ω)]
i, j
hi − k, j (ω)
Taking into account that M [hi,i −j k, j (ω)] =
M |hi,Σ j
The total relative error in the multistage filter amplitude response is found as a sum of the respective errors in the constituent filters. For the phase response error,
D |θi,Σ j
(ω)|
2n − k − 1 j D |hîi −− kk −− 1, p− 1, j (ω)| i, j = |hi−k, j| j |hîi −− kk −− 1, p− 1, j (ω)| k=1 p=1
∑
∑
+
788 2n − 1
∑
l=1
A.F. Emanov et al. / Russian Geology and Geophysics 49 (2008) 780–789 2n − l − 1
|hi,i, jj − l|
∑
g=1
D |hî,i, jj −− ll −1 − g − 1 (ω)| |hi,i, jj −− ll −− 1g− 1 (ω)|
.
The multistage filter schemes differ in both mathematical expectation and scatter of random errors. In fact, the filters have different regularization degrees and different errors. The algorithm accuracy improves as regularization increases, i.e., at narrower frequency ranges. Regularization is adjusted automatically in each conversion scheme, which is useful to estimate the filter accuracy. The scatter of errors in a multistage filter depends on the in-line fold, on the conversion scheme, and on the coherence spectra for each conversion step. All these parameters are easy to find from the measured wave field, and, besides the section of head waves, one can obtain a set of values along the profile which characterize the conversion accuracy. Comparing the accuracy of different conversion schemes, one can choose the best scheme for processing specific field data. The accuracy of a head-waves trace in a section is defined by errors in the filters and by the signal/noise ratio at the input of the multistage filter; the accuracy of the filters, in turn, depends on the in-line fold and the coherence spectrum. When planning a study, one is often free to choose the fold for each specific wave required to achieve the wanted accuracy of the section. However, sometimes one may want to process head waves additionally to the processing of CMP data collected according to the reflection profiling practice. In this case the in-line fold is independent of that for each specific head wave. It often occurs that each wave has a fixed fold, low for low-coverage waves. Thus, there arises the question how to improve the processing accuracy at a fixed fold. The only possibility is to increase γ(ω), for which the are two ways. The coherence spectrum at a fixed frequency is known to be proportional to the signal/noise ratio (see above). The signal includes only output head waves (parallel refractions) while other wave types and random noise are the noise component. Head waves are especially intense at the first arrivals. For a time section of head waves, one specifies the trace length in
order to obtain both first and subsequent arrivals. The intensity of head waves being the greatest at the first arrivals and decreasing with time, the coherence likewise decreases with time. A way to improve the accuracy is to reduce the trace length. By reducing the time window, one simplifies the problem and improves the solution accuracy at a fixed in-line fold. Another way of accuracy improvement is to process the signal to improve the signal/noise ratio in the initial seismograms. Head waves have their apparent velocities different from those of other components. Thus, one can increase the signal/noise ratio, as well as the coherence spectrum (i.e., to improve the quality of seismic sections) by velocity filtering of the initial traces.
Examples of processing field head-wave data Dynamic conversion of head waves into time sections is actually inversion of coherent head waves along parallel correlations for a system of given points. A time section with a user-specified offset is a specific processing case. In fact, it can be both dynamic conversion of head waves and wavefields from a complexly configured system of stacking points. In this study we focus on the theory of head-wave processing, and the cited field examples highlight briefly the potentialities of the suggested processing technique. Application of the developed algorithms is an extensive issue to be discussed in a special publication. The conversion examples are for different data sets. We applied the algorithm to common-depth-surface (CDS) data (Fig. 4) from the Nepa dome (Siberian craton) and to CMP-DSS data. The seismograms were courtesy of people from the Irkutskgeofizika Association. The CDS system recorded head-wave first arrivals from the craton basement surface. The shotpoints were spaced at 200 m along the profile. The spread of receivers had an end-on geometry with offsets from 7200 to 11,800 m, and a spacing of 50 m between the receivers.
Fig. 4. Dynamic head-wave conversion: CDS data. A — time section, L= 10 km; B — wavefield converted into common shotpoint seismogram.
A.F. Emanov et al. / Russian Geology and Geophysics 49 (2008) 780–789
789
Fig. 5. Dynamic head-wave conversion: CMP data. A — Time section, L = 2.5 km, second refractor coverage; B — time section, L = 3.5 km, third refractor coverage.
These head-wave data were fit for conversion without additional analysis, and no coverage study was required. We obtained dynamic time sections and the respective converted wavefields. The apparent velocities varied from 5.4 to 6.7 km/s. Head waves from the basement surface were multiphase, with apparent frequencies of 9 to 14 Hz. The time-distance section (Fig. 4, A) showed prominent basement surface in first arrivals and the wavefield converted into common shotpoint seismogram (Fig. 4, B) allowed recognizing multiple reflected refractions in second arrivals. The second phase was the most intense at an offset of L = 10 km; the first phase was the most prominent in the central part of the section but was comparable with noise in the right-hand part. The wave in the right end of the profile was almost single-phase, the first and second phases dominated in the center, and the second and third phases dominated in the profile left end. Dynamic conversion of head waves is actually a technique providing selection of parallel refractions. The noise/signal ratio increases notably in low-quality fragments. Comparison of field and computed seismograms shows that conversion leaves only head waves in the records while other components disappear. Thus, processing of this kind simplifies the wave pattern and makes the interpretation much easier because wave-type identification is no longer needed. The CMP-DSS data were collected with a main and an additional systems. In the former, the in-line fold was 48, the shotpoints were spaced at 100 m, and the shot-receiver distance was from 0 to 4.8 km; in the latter the same parameters were, respectively, 24, 200 m, and 4.8–9.6 km. Unlike the CDS data, CMP-DSS data processing requires coverage constraining because of wave change. We selected three head waves from kinematic data: for the first wave the coverage was from L = 50 m to 2.5 km, the apparent velocities were 1.8–2.2 km/s, the signal was single-phase, and the visible period was ∼62−66 ms; for the second wave, the offset was L = 2.5 to 3.5 km, the apparent velocities were ∼3.8−4.2 km/s, the signal duration was ∼1.5−2 periods, and the period was
∼60−64 ms; for the third wave, the parameters were, respectively, 3.5 to 9 km, ∼5.6−5.8 km/s, and a lower frequency, i.e., a longer period of ∼46−54 ms. See the sections for the second and third head waves in Fig. 5, where the inphase axis highlights a single refractor in the former (Fig. 5, A) and inphase axes appearing in second arrivals in the latter (Fig. 5, B). Time sections can be obtained successively with each first-arrival head wave for different offsets.
References Bendat, J.S., Piersol, A.G., 1980. Engineering Applications of Correlation and Spectral Analysis. John Wiley & Sons, New York. Emanov, A.F., Seleznev, B.C., Korshik, N.A., 2001. Automated processing of multiple-coverage refracted wave data, in: Setting up a National Network of Reference Geophysical Profiles: Methods of Data Acquisition, Processing, and Interpretation [in Russian]. SNIIGGiMS, Novosibirsk, pp. 145–161. Emanov, A.F., Seleznev, B.C., Solov’ev, V.M., Larkin, G.V., Fateev, A.V., Korshik, N.A., Gritsenko, S.A., Ivanov, N.K., 1998. Automated processing of multiple-coverage head waves, in: Lithospheric Structure: Methods of Investigation and Monitoring [in Russian]. NIC OIGGM SO RAN, Novosibirsk, pp. 197–204. Epinat’eva, A.M., Goloshubin, G.M., Litvin, A.P., Pavlenkin, A.D., Petrashen’, G.I., Starobinets, A.E., Shneerson, M.B., 1990. Method of Head Waves [in Russian]. Nedra, Moscow. Gol’din, S.V., 1974. Linear Conversion of Seismic Signals [in Russian]. Nedra, Moscow. Krylov, S.V., Sergeev, V.N., 1985. Properties of head waves and new possibilities of automation of their processing. Geologiya i Geofizika (Soviet Geology and Geophysics) 26 (4), 92–102 (85–94). Mitrofanov, G.M., Sergeev, V.N., 1986. Investigation of linearized model for head wave with regard to processing CMRW data. Geologiya i Geofizika (Soviet Geology and Geophysics) 27 (8), 98–108 (92–101). Pugachev, B.C., Sinitsin, I.N., 2000. The Theory of Stochastic Systems [in Russian]. Logos, Moscow. Seleznev, B.C., Emanov, A.F., 1998. Conversion of head-wave fields by the Wiener filters. Geologiya i Geofizika (Russian Geology and Geophysics) 39 (4), 536–546 (548–559).
Editorial responsibility: M.I. Epov
Dynamic conversion of head waves in seismic data processing A.F. Emanov *, B.C. Seleznev, N.A. Korshik Geophysical Surveys of the Siberian Branch of the RAS, Altai-Sayan Department, 3 prosp. Akad. Koptyuga, Novosibirsk, 630090, Russia Received 21 August 2007; accepted 26 December 2007
Abstract We suggest a new technique for processing head-wave data with a multistage Wiener filter, which provides conversion and stacking of head waves from a receiver array into a user-specified point. Different conversion algorithms are compared and shown to be a combination of parallel and serial Wiener filters. The Wiener filter accuracy for a conversion iteration depends on the in-line fold and on the coherence spectrum. The latter describes the signal/noise ratio dependence where head waves are the signal and the noise includes all other components and microseisms. Equations for the accuracy of head-wave traces have been derived from a corresponding equation for a single conversion iteration using the multistage Wiener filter schemes and proceeding from the theory of coupling stochastic systems. The obtained equations are applicable to different field data processing algorithms. The discussed conversion algorithms differ in noise regularization and are fit for data of any quality. © 2008, IGM, Siberian Branch of the RAS. Published by Elsevier B.V. All rights reserved. Keywords: Head waves; coherent signals; Wiener filter; seismic exploration
Introduction
Conversion with a single Wiener filter
Processing algorithms for head-wave data are much less developed than for high-fold reflection profiling data. An intention to apply reflection-data processing techniques to head waves gave rise to a common depth surface (CDS) method (Epinat’eva et al., 1990) in which, however, the dynamic parameters of head waves are distorted, and there are some other drawbacks. New algorithms for head-wave data processing became possible on the basis of the dynamic spectral ratio for head waves recorded at four correlated points of a measurement system (Krylov and Sergeev, 1985). The conversion algorithms we suggest combine Wiener filters and dynamic conversion of head waves (Emanov et al., 1998; Seleznev and Emanov, 1998). Head-wave data are processed in three steps: constraining refractor coverage, compressing redundant data within each coverage in the time section, and inversion for subsurface parameters. The suggested conversion algorithms were designed for the case of already constrained coverage, and in this paper we focus on the second processing step, that of picking head waves from the measured wavefield.
Krylov and Sergeev (1985) derived equations for headwave fields measured at four correlated points (Fig. 1): Fi, j(ω)Fi + 1, j + 1(ω) = Fi + 1, j(ω)Fi, j + 1(ω),
(1)
where F(ω) are the spectra of head waves at stacking points, and i and j are the numbers of sources and receivers, respectively. According to (1), the wavefield can be converted from one stacking point to another using a linear system. The conversion from the point (i, j) into the point (i + 1, j) is governed by the spectral equation ^ ^ Fi + 1, j(ω) = Fi, j(ω)hii, i + 1(ω),
(2)
where hii, i + 1(ω) is the frequency response of the conversion ^ linear filter, Fi, j(ω) is the spectrum of the signal plus random ^ noise at the point (i, j), and Fi + 1, j(ω) is the spectrum of the seismogram converted into (i + 1, j). For head waves, using (1), one can write the condition of equality of filters for multiple correlations (Fig. 1)
* Corresponding author.
E-mail address: [email protected] (A.F. Emanov) hii, i + 1(ω) = hii,+i +1 1(ω) = hii,+i +2 1(ω) 1068-7971/$ - see front matter D 2008, IGM, Siberian Branch of the RAS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.rgg.2007.12.015
+ … + hii,+i +k 1(ω).
(3)
781
A.F. Emanov et al. / Russian Geology and Geophysics 49 (2008) 780–789
averaged autocorrelation function at the input of the linear system, we obtain R12(θ) − hopt(τ)∗R11(τ − θ) = 0
(5)
or, in the frequency domain, hopt(ω) =
Fig. 1. Observation system and ray scheme. 1 — i-th source; 2 — j-th receiver; 3 — point into which head waves are converted; 4 — point from which head waves are converted; 5 — points used for filter construction. R is refractor.
This equation is valid in the absence of noise in seismograms and implies that conversion from (i, j) into (i + 1, j) is possible with a filter obtained for any correlation. In the presence of noise, when the filters for different correlations are dissimilar, one can apply multiple conversion of the same signal with different filters and its equilibrium (Krylov and Sergeev, 1985) or weighted (Mitrofanov and Sergeev, 1986) stacking. There is an alternative way of conversion, assuming, according to (3), that a set of linear systems is actually a set of multiple input and output signal and noise realizations of the same linear system. Thus, we have a linear system with a known set of input-output signal pairs recorded on the background of random noise and a single realization of an input head-wave record, which is to be converted into an output record. We search an optimal conversion filter using the approach applicable to single Wiener filters (Bendat and Piersol, 1980; Gol’din, 1974; Seleznev and Emanov, 1998). Let the input signal of the frequency h(t) be FI, j(t) and the output signal be Fi, j+1(t) = Fi, j(t)*h(t), where (*) denotes convolution. One ^ can measure only Fi, j(t) = Fi, j(t) + Wi, j(t), where Wi,j(t) is the random noise realization. The optimal Wiener filter to convert ^ the input signal Fi, j(t) into the output signal Fi, j+1(t) is found by minimizing, along h(t), the mathematical expectation of the ^ square misfit between the filtered Fi, j(t) signal and the true output signal Fi, j+1(t): ^ (4) M |Fi, j + 1(t) − Fi, j (t)∗hopt(t)2| = min. Zeroing the derivative of (4) along h(t), after simple transformation, gives ^ M [Fi, j + 1(t)Fi, j(t − θ)] − ^ ^ hopt(τ)∗M [Fi, j (t − τ)Fi, j (t − θ)] = 0. ^ Taking into account that M [Fi, j + 1(t)Fi, j (t − θ)] = R12(θ) is the averaged cross-correlation function between the input ^ ^ and output signals and M [Fi, j (t)Fi, j (t − θ)] = R11(τ − θ) is the
^∗ M [Fi, j + 1(ω)Fi, j(ω)] = , ^ 2 R11(ω) M |Fi, j(ω)|
R12(ω)
(6)
where (*) is the spectrum complex conjugation. The main problem is how to calculate the average reciprocal and auto spectra in (6). We have a set of input and output signal realizations against random noise realizations. The spectra can be averaged over the realization set using records along parallel correlations. Then, ^ ^ R12(ω) = ∑ Fi + k, j + 1(ω)F∗i + k, j(ω), k
^ R11(ω) = ∑ |Fi + k, j (ω)|2, k
and the equation for the head-wave conversion from (i, j) into (i, j + 1) becomes ^
^
^ Fi, j + 1(ω) = Fi, j (ω)
∑ [Fi + k, j + 1(ω)F∗i + k, j(ω)] k
.
^
∑ |Fi + k, j (ω)|2
(7)
k
A conversion iteration in case (1) is from three seismograms while each iteration in (7) uses a number of seismogram pairs from parallel correlations. In case (1), for a measurement system as in Fig. 1, the same wave is converted and stacked in many iterations, while (7) requires a single conversion iteration with a more accurate filter based on many seismograms. It is easy to show that ∗
Mhopt(ω) =
∗
M {[Fi, j + 1(ω) + Wi, j + 1(ω)] [Fi, j(ω) + Wi, j(ω)]} 2
M |Fi, j (ω) + Wi, j (ω)|
=
h (ω) , N (ω) 1 + __ 2 |F (ω)|
__ where N(ω) is the noise power spectrum and N (ω)/|F(ω)|2 is the square noise/signal ratio averaged over the set of seismograms. The mathematical expectation of the optimal Wiener filter for conversion of head waves differs from the true filter in its regularization with respect to the average noise/signal ratio. The conversion procedure for an observation system of the fold n is as shown in Fig. 2. Conversion into O is required only for seismograms from the correlations that cross at this point (Seleznev and Emanov, 1998). Assuming that each signal from the selected correlations is converted into the point O with a single filter, we obtain a single conversion option, using the seismograms from the triangle only (Fig. 2, A). This scheme has a drawback that ever fewer parallel correlations are used for obtaining the filters as the distance between the starting and target points increases. As a result, the mathematical expectations of the reciprocal and auto spectra are estimated from different realization numbers, and the filters differ in accuracy.
782
A.F. Emanov et al. / Russian Geology and Geophysics 49 (2008) 780–789
Fig. 2. Schemes for conversion of head wavefield into a point. — triangle; B — parallelogram; C — trapezium. 1 — point into which wave filed is converted; 2 — point from which wave field is converted; 3 — seismogram points used for filter construction.
It is possible to provide a uniform accuracy of the conversion filters by using a parallelogram as in (Fig. 2, B) for a correlation instead of the single-step conversion. The seismogram from the top point is converted into the point next to the filter calculated from (2n – 3) correlations; then the sum of two seismograms is converted a step further and is summed with the third one, and thus the whole wavefield becomes stacked successively at a given point. The number of averaging iterations is (2n – 3) for any step. The back correlation points are likewise converted using a parallelogram as in Fig. 2, B, but with an opposite tilt, and the conversion of two correlations together covers a trapezium-shaped area (Fig. 2, C).
Conversion from two stacking points into one We convert signals into the point O(i, j) (Fig. 2, C) from two arbitrary points which cross at O, one of a forward (i, j + b) and the other of a back (i – a, j) correlation. Conversion using a two-stage Wiener filter (Emanov et al., 2001) follows ^ ^ F′i, j (ω) = Fi − a, j (ω)Li,i −j a, j(ω) + Fi, j + b (ω)Li,i, jj + b (ω),
(8)
where Li,i −j a, j (ω) and Li,i, jj + b (ω) are the frequency responses of the Wiener filter for conversion from two points into one; prime discriminates between the measured and computed records. The responses of the two-stage Wiener filter are found using a standard approach, by minimizing the mathematical expectation of the square misfit
^ M |Fi, j (ω) − Li,i −j a, j (ω)Fi − a, j (ω) − ^ Li,i, jj + b (ω)Fi, j + b (ω)|2 = min.
(9)
Differentiating along each filter and zeroing the derivatives gives j i, j i, j + b i, j Rii −− a, a, j (ω)Li − a, j (ω) + Ri − a, j (ω)Li, j + b (ω) −
Ri,i −j a, j (ω) = 0,
(10)
Ri,i, jj ++ bb (ω)Li,i, jj + b (ω) + Rii,−j +a,bj (ω)Li,i −j a, j (ω) − Ri,i, jj + b (ω) = 0, ^ 2 j where Rii −− a, a, j (ω) = M |Fi − a, j (ω)| is the mathematical expectation of the trace square spectrum at the stacking point (i – a, ^ j); Ri,i, jj ++ bb (ω) = M |Fi, j + b (ω)|2 is the same at (i, j + b); ^ ^ Ri,i −j +a,bj (mω) = Rii,−j +a,bj (ω) = M [Fi − a, j (ω)F∗i, j + b (ω)] is the reciprocal spectrum between traces at (i – a, j) and (i, j + b) averaged over the set of realizations; Ri,i −j a, j (ω) and Ri,i, jj + b (ω) are the mathematical expectations of the reciprocal spectra between the records at (i, j) and (i – a, j) and at (i, j) and (i, j + b); Ri,i −j a, j (ω) and Ri,i, jj + b (ω) are the same for the pairs (i – a, j) and (i, j) and (i, j + b) and (i, j), respectively. The frequency responses of the two-stage Wiener filter are
Li,i −j a, j (ω) =
i, j + b i, j Ri − a, j (ω)Ri, j + b (ω) i, j Ri − a, j (ω) 1 − i, j i, j + b Ri − a, j (ω)Ri, j + b (ω) i − a, j
Ri − a, j (ω) 1 − γab (ω) 2
,
783
A.F. Emanov et al. / Russian Geology and Geophysics 49 (2008) 780–789
Li,i, jj + b (ω) =
i − a, j i, j Ri, j + b (ω)Ri − a, j (ω) i, j Ri, j + b (ω) 1 − i, j i − a, j Ri, j + b (ω)Ri − a, j (ω) 2 i, j + b Ri, j + b (ω) 1 − γab(ω)
^
∑ Li,i, jj + b (ω)Fi, j + b (ω)|2 = min. ,
(11)
2
where γ2ab(ω) =
Ri, j + b (ω) i − a, j i − a, j i, j + b Ri − a, j (ω)Ri, j + b (ω)
is the square coherence spec-
trum between the traces at (i – a, j) and (i, j + b). These are the equations for a standard two-stage Wiener filter. The mathematical expectations of the reciprocal and auto spectra in (11) are found using (3), i.e., the spectra are averaged over parallel correlations. Then, ^ ^ Ri,i −j a, j(ω) = ∑ Fi − a + k, j (ω)Fi∗+ k, j(ω), k
where the subscript k marks passage to parallel correlations and is either positive or negative, but never zero, ^ 2 j Rii −− a, a, j (ω) = ∑ |Fi − a + k, j (ω)| , k
Ri,i, jj ++ bb (ω)
^ = ∑ |Fi, j + b + k (ω)|2, k
^ ^ Ri,i, jj + b (ω) = ∑ Fi, j + b + k (ω)Fi,∗j + k (ω), k
^ ^ Ri,i −j +a,bj (ω) = ∑ Fi − a + k, j (ω)Fi,∗j + b + k (ω). k
The spectral responses of the two-stage Wiener filter include single Wiener filters and a frequency-dependent term as a multiplier which provides weighing in the stacking of converted seismograms.
A multistage conversion Wiener filter Judging by the above result, we may expect to improve the conversion scheme by employing a multistage Wiener filter to convert head waves simultaneously from all points of the forward and back correlations into a stacking point, using the equation (Gol’din, 1974) ^ F′i, j (ω) = ∑ Fi − a, j (ω)Li,i −j a, j (ω) + a
^
∑ Fi, j + b (ω)Li,i, jj + b (ω),
(12)
b
where a and b are nonzero integers which provide passage from one seismogram to another along the correlations. Conversion of waves from a given set of points into a single point requires an algorithm for a multistage Wiener filter. First we minimize the mathematical expectation of the square misfit: ^ M |Fi, j (ω) − ∑ Li,i −j a, j (ω)Fi − a, j (ω) − a
(13)
b
Zeroing the derivatives along each filter gives a system of equations written as the matrix ~ − 1), j Ri, j (ω) Ri − a, j(ω) (ω) Rii −− (a a, j ~ i − a, j i − a, j i − a, j i − (a − 1), j Ri, j ~ i − (a − 1), j(ω) Ri − (a − 1), j(ω) Ri − (a − 1), j(ω) ⋅⋅ i − (a − 1), j Ri, j i − a, j i − (a − 2), j(ω) Ri − (a − 2), j(ω) Ri − (a − 2), j(ω) ⋅⋅ ⋅ ⋅ ⋅ ⋅⋅ ⋅⋅ ⋅ ~ ⋅ ⋅ ⋅⋅ i, j i − a, j − 1), j (ω) Rii −+ (a Ri + 1, j(ω) Ri + 1, j(ω) ⋅⋅ 1, j ~ i, j = i − a, j i − (a − 1), j ⋅⋅ (ω) (ω) (ω) R R R i, j + b i, j + b i, j + b ~ i, j i − a, j ⋅⋅ i − (a − 1), j Ri, j + (b − 1)(ω) Ri, j + (b − 1)(ω) Ri, j + (b − 1)(ω) ⋅⋅ ~ i, j i − a, j ⋅⋅ Ri, j + (b − 2)(ω) Ri, j + (b − 2)(ω) Rii,−j +(a(b− −1),2)j(ω) ⋅⋅ ⋅ ⋅ ⋅ ⋅⋅ ⋅ ⋅ ⋅ ~ Ri, j (ω) Ri − a, j(ω) Ri, j + 1i − (a − 1), j(ω) i, j + 1 i, j + 1 j Rii +− 1, a, j(ω)
Ri,i −j +a,bj(ω)
Rii,+j +1,bj(ω)
Ri,i, jj ++ bb(ω)
j i, j + b Rii +− 1, (a − 1), j(ω) Ri − (a − 1), j(ω) j i, j + b Rii +− 1, (a − 2), j(ω) Ri − (a − 2), j(ω) ⋅ ⋅ ⋅ ⋅ i + 1, j i, j + b Ri + 1, j(ω) Ri + 1, j(ω) j i, j + b Rii,+j +1,(b − 1)(ω) Ri, j + (b − 1)(ω) j i, j + b Rii,+j +1,(b − 2)(ω) Ri, j + (b − 2)(ω) ⋅ ⋅ ⋅ ⋅ Rii,+j +1,1j(ω) Ri,i, jj ++ b1(ω)
Ri,i −j +a,(bj − 1)(ω) ⋅ ⋅ ⋅⋅ Ri,i −j +(a(b− −1),1)j(ω) ⋅⋅ Ri,i −j +(a(b− −2),1)j(ω) ⋅ ⋅ ⋅⋅ ⋅ ⋅⋅ ⋅ i, j + (b − 1) (ω) ⋅ ⋅ Ri + 1, j ⋅⋅ − 1) (ω) ⋅ ⋅ Ri,i, jj ++ (b b − 1) ⋅⋅ Ri,i, jj ++ (b (b − 1)(ω) ⋅⋅ i, j + (b − 1) Ri, j + (b − 2)(ω) ⋅ ⋅ ⋅ ⋅⋅ ⋅ ⋅⋅ Ri,i, jj ++ 1(b − 1)(ω) ⋅ ⋅
Li, j i − 1, j i, j L i − (a − 1), j Li, j i − (a − 2), j ⋅ ⋅ ⋅ ⋅ i, j i, j + 1 Ri + 1, j(ω) Li + 1, j (14) × i, j , i, j + 1 Ri, j + b(ω) Li, j + b i, j Ri,i, jj ++ 1(b − 1)(ω) Li, j + (b − 1) i, j Li, j + (b − 2) Ri,i, jj ++ 1(b − 2)(ω) ⋅ ⋅ ⋅ ⋅ i, j + 1 Ri, j + 1(ω) i, j L 1 + i, j ^ i − a, j where Ri − a, j(ω) = M |Fi − a, j (ω)|2 is the mathematical expectation of the trace square spectrum at the stacking point (i – a, Ri,i −j +a,1j(ω) Ri,i −j +(a1− 1), j(ω) Ri,i −j +(a1− 2), j(ω)
784
A.F. Emanov et al. / Russian Geology and Geophysics 49 (2008) 780–789
^ j); Ri,i, jj ++ bb (ω) = M |Fi, j + b (ω)|2 is the same at (i, j + b); ^ ^ Ri,i −j +a,bj (ω) = Rii,−j +a,bj (ω) = M [Fi − a, j (ω)F∗i, j + b (ω)] is the reciprocal spectrum between traces at (i – a, j) and (i, j + b) averaged over the set of realizations; ^∗ ^ j i − a, j Rii +− 1, a, j (ω) = Ri + 1, j (ω) = M [Fi − a, j (ω)Fi + 1, j (ω)] is the auto spectrum between traces at (i – a, j) and (i + 1, j) averaged over the set of realizations; − 1) (ω) = Ri,i, jj ++ b(b − 1) (ω) = Ri,i, jj ++ (b b ^ ^ M [Fi, j + b (ω)F∗i, j + (b − 1) (ω)]
is the auto spectrum between traces at (i, j + b) and (i, j + (b – 1)) averaged over the set of realizations; ~ ~ ^ ^ Ri,i −j (a − 1), j (ω) = Rii,−j (a − 1), j (ω) = M [Fi − (a − 1), j (ω)F∗i, j (ω)] ~ and Ri,i, jj + b (ω) are the mathematical expectations of the reciprocal spectra between the pair (i – (a – 1), j) and (i, j) and the pair (i, j + b) and (I, j); they are the same for the other points. The mathematical expectations in auto and reciprocal spectra are substituted by averaging over parallel correlations. This gives a matrix for which the solution can be obtained from the column vector of the spectra of the multistage Wiener filter. Solving (14) and converting head waves from many points into one point with a multistage Wiener filter is a standard problem. The multistage Wiener filter is coupled with dynamic conversion of head waves at the step of averaging of the auto and reciprocal spectra (matrix elements) over parallel correlations.
The accuracy of head-wave conversion into a time section The idea of the accuracy of a head-wave seismic section stems from the meaning of the latter. In the reported processing technique, a head-wave time section for a given offset is meant as a set of traces containing head waves only. The algorithms are actually a tool for selecting parallel refractions, which are the head waves that arrive at a receiver. Any head-wave section is obviously approximate. Below we discuss how the accuracy depends on the in-line fold and on the wavefield features. The accuracy is estimated for a whole trace consisting of several samples within a given time window rather than for some individual parameter. The accuracy of a trace is meant as a difference between the computed and true values, the true value being a time function describing head waves at a given receiver. The accuracy of a filter obtained according to (7) is estimated as follows. The uncertainty in (7) is found as (Bendat and Piersol, 1980) M {[h′0i(ω) − h0i(ω)]2} = M {[h′0i(ω) − M [h′0i(ω)]]2} + M {[M [h′0i(ω)] − h0i(ω)]2},
where h′0i is the computed response, the first term is the scatter D[h′0i(ω)] and the second term is the square misfit b2[h′0i(ω)]. As it was shown before, the misfit leads to regularization of filter’s amplitude response, which improves the solution stability, while the uncertainty in filter construction is defined by the scatter. To investigate the random error in the conversion filter obtained by (7), we consider separately the phase and amplitude errors. According to (7), the phase response of the conversion filter is found as its mathematical expectation for the reciprocal spectrum between the filter input and output. The reciprocal spectrum consists of a real and an imaginary parts: |R′0i(ω)|2 = C20i (ω) + Q20i (ω), and tg θ0i (ω) = C0i (ω)/Q0i (ω),
(15)
where C0i(ω) is the real part, Q0i(ω) is the imaginary part, and θ0i(ω) is the phase spectrum. The formula for random error in the filter’s phase response is derived using the mathematical expectations of the square real and imaginary parts of the reciprocal spectrum 1 M [C20i (ω)] = [R0(ω)Ri (ω) + 3C20i (ω) − Q20i (ω)], 2
(16)
1 M [Q20i (ω)] = [R0(ω)Ri (ω) + 3Q20i (ω) − C20i (ω)].
(17)
2
Equations (16) and (17) are easy to derive if the real and imaginary parts of the reciprocal spectrum are expressed via the real and imaginary spectrum parts at the filter input and output, squared and defined as mathematical expectations. The scatter in the real part of the spectrum is D [C0i (ω)] = M [C20i (ω)] − C20i (ω) = 1 2
[R0(ω)Ri (ω) + C20i (ω) − Q20i (ω)],
(18)
and that for the imaginary part is D [Q0i (ω)] = M [Q20i (ω)] − Q20i (ω) = 1 2
[R0(ω)Ri (ω) + Q20i (ω) − C20i (ω)].
(19)
The expression for the covariance between the real and imaginary parts of the reciprocal spectrum can be obtained in the same way: cov[C′0i (ω), Q′0i (ω)] = M [C′0i (ω)Q′0i (ω)] − C0i (ω)Q0i (ω) = C0i (ω)Q0i (ω).
(20)
Differentiation along both terms of (15) and finite-difference approximation gives sec2θ0i (ω)⋅∆θ0i (ω) ≈
C0i (ω)∆Q0i (ω) _ Q0i (ω)∆C0i (ω) 2
C0i (ω)
.
Taking into account that sec2θ0i (ω) = |R0i (ω)|2/C20i (ω), we obtain
A.F. Emanov et al. / Russian Geology and Geophysics 49 (2008) 780–789
∆θ0i (ω) ≈
C0i (ω)∆Q0i (ω) − Q0i (ω)∆C0i (ω) 2
|R0i (ω)|
.
(21)
Thus, the random error scatter in the filter phase response, without averaging, is D [θ′0i (ω)] = M [(∆θ0i (ω))2] ≈ 1 4
M [(C0i (ω)∆Q0i (ω) − Q0i (ω)∆C0i (ω))2] =
4
{C20i (ω)D [Q′0i (ω)] −
|R0i (ω)| 1 |R0i (ω)|
2C0i (ω)Q0i (ω) cov [C′0i (ω), Q′0i (ω)] + Q20i (ω)D [C′0i (ω)]}. Substitution of (18), (19), and (20) into (21) gives D [θ′0i (ω)] ≈ 2 (1 − γ0i (ω)) 2 2γ0i (ω)
4 |R0i (ω)|2 |R0i (ω)| (ω) − (ω)R R = 0 i 4 2 2 |R0i (ω)|
1
,
(22)
where 2
γ (ω) =
|R0i (ω)|
r00(ω)Rii (ω)
,
(23)
where R00(ω) is the input spectrum of the autocorrelation function of traces and Rii(ω) is the same at the filter output. Using model (1) gives R0i (ω) = h0i (ω)M |F0(ω)|2, R00(ω) = M |F0(ω)|2 + M |W0(ω)|2, 2
2
2
Rii (ω) = |h0i (ω)| M |F0(ω)| + M |Wi (ω)| . These equations are valid assuming the independence of input and output noise realizations and the absence of their correlation with the signal. Substituting these equations into (23) gives γ2(ω) =
1 , (1 + α0(ω)) (1 + αi (ω)) 2
where α0(ω) =
M |W0(ω)|
2
M |F0(ω)|
(24) 2
, αi (ω) =
M |Wi (ω)|
2
M |Fi (ω)|
are the input/out-
put ratios of weighted mean square noise/signal ratios. Taking into account that the phase spectra are averaged over (2n – 3) independent records, the random error scatter may be expected to become (2n – 3) times as low. Therefore, the rms error in the filter phase response is σθ(ω) ≈
The conversion filter is constructed using a set of input and output signal and noise records. Obtaining the explicit formula for the amplitude error is hard for a model implying noise both at the input and at the output but is easy (Bendat and Piersol, 1980) for the case when the input noise is vanishing relative to the signal. Inasmuch as only the γ20i(ω) value matters for the filter uncertainty and it is of no importance whether input or output noise is involved, the corresponding expression for this specific case will be valid in the general case as well. The records of head waves from the same refractors at two different points are related as linear equation system (2). In fact, we have the records of head waves mixed with the noise ^ Fi, j (t) = Fi, j (t) + Wi, j (t), (26) Fi + 1, j (t) = Fi + 1, j (t) + Wi + 1, j (t) at the input and output of the system. Equations (2) and (26) constitute a model of signals which lies at the base of the conversion algorithms. As a consequence of (2), in linear systems of this kind, Ri + 1, j (ω) = |h (ω)|2Ri, j (ω),
2
2 1 − γ0i(ω) √
|γ0i(ω)| √ 2(2n − 3)
.
(25)
Using (25) we can find the frequency-dependent phase response in radians. The same estimation for amplitude is more complicated. One might calculate the confidence intervals of the filter amplitude response at each conversion iteration (Bendat and Piersol, 1980) but having an explicit formula is much more convenient.
785
(27)
where Ri + 1, j (ω) = M |Fi + 1, j (ω)|2 and Ri, j (ω) = M |Fi, j (ω)|2 are the auto spectra of head waves at the system input and output. Calculation of finite differences for both terms of (27) gives ∆Ri + 1, j (ω) ≈ |h (ω)|2∆Ri, j (ω) + Ri, j (ω)∆|h (ω)|2.
(28)
Having transformed (28) into Ri, j (ω)∆|h (ω)|2 ≈ ∆Ri + 1, j (ω) − |h (ω)|2∆Ri, j (ω) and having found the mathematical expectation for both terms of this equation, we obtain R2i, j (ω)D [|h (ω)|2] ≈ D [Ri + 1, j (ω)] − 2|h (ω)|2 × cov [Ri, j (ω), Ri + 1, j (ω)] + |h (ω)|4D [Ri, j (ω)].
(29)
In the further consideration we assume that noise is vanishing relative to the signal at the point (i, j) but increases correspondingly at (I + 1, j). This assumption allows obtaining per se input and output amplitude spectra of head waves from measured data. Assuming vanishing noise at (i, j), we apply the Fourier transform to find the head wave spectrum and then obtain their spectrum at the filter output with reference to the coherence spectrum: ^ Ri + 1, j (ω) = γ2(ω)Ri + 1, j (ω), (30) ^ 2 (ω) = − γ (ω), (ω)]R [1 RW i + 1, j i + 1, j ^ ^ where Ri + 1, j (ω) = M |Fi + 1, j (ω)|2 is the auto spectrum of the seismogram at (i + 1, j) and RW i + 1, j (ω) is the noise auto spectrum at the same point. Thus, with the above assumptions, one can identify the input and output amplitude spectra of head waves and the amplitude spectrum of noise.
786
A.F. Emanov et al. / Russian Geology and Geophysics 49 (2008) 780–789
The scatter in the estimate of the input auto spectrum is, by definition, D [Ri, j (ω)] = M [R2i, j (ω)] − R2i, j (ω), and the mathematical expectation of the square auto spectrum of a stochastic process of a duration T will be 2R2i, j (ω) (Bendat and Piersol, 1980); correspondingly, D [Ri, j (ω)] = R2i, j (ω). The random error in the auto spectrum is fairly large and its scatter is independent of the respective duration. The problem is that the resolution of the spectrum increases with the time window while the accuracy remains invariable. In our case, finding the mathematical expectation is replaced by averaging over parallel correlations, whereby the random error scatter diminishes by a factor of (2n – 3), to give D [Ri, j (ω)] = R2i, j (ω)/(2n − 3). ^ As it follows from (26), D [Ri + 1, j (ω)] = D [Ri + 1, j (ω)] + D [RW i + 1, j (ω)], and the random error scatter in the spectrum of head waves at (i + 1, j) is obtained using (30) as 2
D [Ri + 1, j (ω)] =
2
[2 − γ (ω)]Ri + 1, j (ω) 2
γ (ω) (2n − 3)
(31)
and cov [Ri, j (ω), Ri + 1, j (ω)] = M [∆Ri, j (ω)∆Ri + 1, j (ω)] = Ri, j (ω)Ri + 1, j (ω) 2n − 3
.
(32)
Substitution of (30)–(32) into (29) gives D [|h (ω)|2] ≈
2
4
2[1 − γ (ω)] |h (ω)| 2
γ (ω)[2n − 3]
.
(33)
Using the equation D [A2] = 4A2D [A] (Bendat and Piersol, 1980), we obtain
D [|h (ω)|] =
2
2
[1 − γ (ω)]⋅|h (ω)| 2
2γ (ω)⋅[2n − 3]
,
(34)
wherefrom there follows the expression for the relative filter amplitude error 2 1/2 ε[|h (ω)|] = σ[|h (ω)|] = [1 − γ (ω)] .
|h (ω)|
|γ(ω)|⋅√ 2(2n − 3)
(35)
According to (25) and (35), the accuracy of the conversion filter improves as the in-line fold increases. We have discussed the accuracy of a Wiener filter for conversion of head waves from one point to another but actually we concern about the accuracy of the trace obtained by conversion from many points. We showed before (Seleznev and Emanov, 1998) that conversion involves head waves from points of the forward and back correlations that cross at a given point. Then, we deal with an algorithm, which brings together many Wiener filters according to a certain scheme (Pugachev and Sinitsin, 2000), and study their properties as of a single whole proceeding from the properties of the constituent linear systems. Thus, the conversion algorithms in our case correspond to a multistage Wiener filter. The above conversion schemes are actually various schemes of a multistage Wiener filter with different possibilities for solving the same problem. See Fig. 3 for flow charts of several conversion schemes. They are a parallel scheme (Fig. 3, A) in the case of a single conversion (Fig. 2, A) and a serial scheme (Fig. 3, B) for successive conversion of the wavefield into a point (Fig. 2, B) when the computed signal is summed with the measured one after each conversion iteration. Inasmuch as the input signal enters the filter series after each element, the scheme is not purely serial. The scheme C in Fig. 3 is a combination of serial and parallel Wiener filters and is the same as B but summation is at the
Fig. 3. Multistage Wiener filter algorithms of head-wave data processing. — parallel conversion; B — multistage filter corresponding to conversion scheme of Fig. 2, B; C — mixed serial-parallel multistage filter equivalent to case B.
787
A.F. Emanov et al. / Russian Geology and Geophysics 49 (2008) 780–789
end of the complete conversion procedure. There can be other possibilities besides the schemes of Fig. 3 but the latter are good examples because they show the extreme cases. A multistage Wiener filter can be expressed as a single whole via a pulse response. In our case, the pulse response is obtained from measured seismograms and is, because of noise, an approximation to the true response. We write a multistage Wiener filter as a mathematical expectation of the pulse response and its scatter. The response for the parallel scheme (Fig. 3, ) is the sum of frequency responses of the constituent elements (Pugachev and Sinitsin, 2000). Since the frequency responses of the filter elements bear an error, the mathematical expectation of the multistage filter response will be the sum of the mathematical expectations of the individual filters. For the scheme of Fig. 3, A it is easy to show that 2n − 1
2n − 1
k=1
l=1
M [hi,Σ j (ω)] = ∑ M [hi,i −j k, j (ω)] + ∑ M [hi,i, jj − l (ω)].
1+
2 i, j M |Wi − k, j (ω)| 2 i, j M |Fi − k, j (ω)|
and that, according to the latter equation, the mathematical expectation of the phase-amplitude response is independent of noise, we obtain the amplitude response of the multistage filter 2n − 1
(ω)| = ∑
2
k=1
2n − 1
∑ 1 + M |W
2
1 + M |Wi − k, j (ω)| /M |Fi − k, j (ω)|
2
2
(ω)| /M |Fi, j − l (ω)|
2n − 1
2n − 1
k=1
l=1
2n − 1 2 i, j i, j [1 − (γi − k, j (ω)) ] |hi − k, j|
2n − 1 2 i, j i, j [1 − (γi, j − l (ω)) ] |hi, j − l (ω)| 2
i, j
2
i, j
2(γi − k, j (ω)) [2n − 3]
k=1
2(γi, j − l (ω)) [2n − 3]
2n − 1
(ω)| = ∑
k=1
i, j
2 1/2
i, j
|γi − k, j (ω)|⋅√ 2(2n − 3)
l=1
[1 − (γi, j − l (ω)) ]
.
2
i, j
2(γi, j − l (ω)) (2n − 3)
In the case of serial conversion (Fig. 2, B) performed using the equivalent schemes B and C of Fig. 3, we actually deal with a multistage filter of several parallel series of Wiener filters. First we consider conversion with a series of filters. For a series of (2n – 1) filters, the mathematical expectation of the frequency response (hP) is the product of the mathematical expectations of the frequency responses of the elementary filters: 2n − 1
M |hi,P j (ω)| = ∏
k=1
i − k, j
|hi − k, j (ω)| 2
2
1 + M |Wi − k, j (ω)| /M |Fi − k, j (ω)|
.
For a series of Wiener filters, it is easy to obtain the equation for the relative error in the amplitude response as a sum of relative errors of the responses in the product (ω)| = ∑
j ε|hii −− k, k, j
2n − 1
(ω)| = ∑
k=1
i − k, j
2
[1 − (γi − k, j (ω)) ] i − k − 1, j
|γi − k, j
2n − 1
2n − 1
k=1
k=1
D [θi,Σ j (ω)] = ∑ D [θii −− kk,−j 1, j (ω)] = ∑
2(2n − 3) (ω)|⋅√
.
i − k − 1, j
[1 − (γi − k, j i − k − 1, j
|γi − k, j
2
(ω)) ]
2(2n − 3) (ω)|⋅√
.
For a combination of serial and parallel multistage filters (C in Fig. 3), the mathematical expectation of its frequency response becomes
+
l=1
k=1
i, j
i, j
∏
g=1
|γi, j − l (ω)|⋅√ 2(2n − 3)
M [hi,i, jj − l − g (ω)] .
2n − 1
.
p=1
M [hi,i −j k − p, j (ω)] +
The frequency response is obtained in this case assuming that summation applies to the relative errors of elements in the series and to the response scatters of each parallel series. Following the rules from (Pugachev and Sinitsin, 2000), the equation for the scatter of the frequency response errors in the mixed serial-parallel multistage Wiener filter is easily obtained as
2 1/2
[1 − (γi, j − l (ω)) ]
∏
2n − 1 2n − l − 1 l=1
2n − 1
+∑
M [hi,Σ j (ω)] = ∑
∑
.
[1 − (γi − k, j (ω)) ]
+∑e
2
2(γi − k, j (ω)) (2n − 3)
2
i, j
−jω
2n − 1 2n − k − 1
With the latter equation it is easy to pass to the relative error of the multistage filter amplitude response: ε|hi,Σ j
2n − 1
[1 − (γi − k, j (ω)) ] i, j
i, j
+ ∑ e−jωD [θi, j − l (ω)] =
The scatter of the phase response of the series is the sum of the respective scatters for each filter:
.
Now we consider separately the scatter of the absolute frequency response of the multistage filter and the phase response scatter. For the former, with regard to (35),
∑
∑e
k=1
2n − 1 l=1
2
i, j
−jω
k=1
D [hi,Σ j (ω)] = ∑ D [hi,i −j k, j (ω)] + ∑ D [hi,i, jj − l (ω)].
l=1
2n − 1
ε|hi,P j
For the case of parallel conversion, the scatter of the frequency response equals the sum of the respective scatters of the filter elements:
D |hi,Σ j (ω)| = ∑
=∑e
i, j
−jωD [θi − k, j (ω)]
k=1
2n − 1
+
i, j
|hi, j − l(ω)|
i, j − l
l=1
e
i, j
|hi − k, j (ω)|
2n − 1
i, j
−jωD [θSIGMS (ω)]
i, j
hi − k, j (ω)
Taking into account that M [hi,i −j k, j (ω)] =
M |hi,Σ j
The total relative error in the multistage filter amplitude response is found as a sum of the respective errors in the constituent filters. For the phase response error,
D |θi,Σ j
(ω)|
2n − k − 1 j D |hîi −− kk −− 1, p− 1, j (ω)| i, j = |hi−k, j| j |hîi −− kk −− 1, p− 1, j (ω)| k=1 p=1
∑
∑
+
788 2n − 1
∑
l=1
A.F. Emanov et al. / Russian Geology and Geophysics 49 (2008) 780–789 2n − l − 1
|hi,i, jj − l|
∑
g=1
D |hî,i, jj −− ll −1 − g − 1 (ω)| |hi,i, jj −− ll −− 1g− 1 (ω)|
.
The multistage filter schemes differ in both mathematical expectation and scatter of random errors. In fact, the filters have different regularization degrees and different errors. The algorithm accuracy improves as regularization increases, i.e., at narrower frequency ranges. Regularization is adjusted automatically in each conversion scheme, which is useful to estimate the filter accuracy. The scatter of errors in a multistage filter depends on the in-line fold, on the conversion scheme, and on the coherence spectra for each conversion step. All these parameters are easy to find from the measured wave field, and, besides the section of head waves, one can obtain a set of values along the profile which characterize the conversion accuracy. Comparing the accuracy of different conversion schemes, one can choose the best scheme for processing specific field data. The accuracy of a head-waves trace in a section is defined by errors in the filters and by the signal/noise ratio at the input of the multistage filter; the accuracy of the filters, in turn, depends on the in-line fold and the coherence spectrum. When planning a study, one is often free to choose the fold for each specific wave required to achieve the wanted accuracy of the section. However, sometimes one may want to process head waves additionally to the processing of CMP data collected according to the reflection profiling practice. In this case the in-line fold is independent of that for each specific head wave. It often occurs that each wave has a fixed fold, low for low-coverage waves. Thus, there arises the question how to improve the processing accuracy at a fixed fold. The only possibility is to increase γ(ω), for which the are two ways. The coherence spectrum at a fixed frequency is known to be proportional to the signal/noise ratio (see above). The signal includes only output head waves (parallel refractions) while other wave types and random noise are the noise component. Head waves are especially intense at the first arrivals. For a time section of head waves, one specifies the trace length in
order to obtain both first and subsequent arrivals. The intensity of head waves being the greatest at the first arrivals and decreasing with time, the coherence likewise decreases with time. A way to improve the accuracy is to reduce the trace length. By reducing the time window, one simplifies the problem and improves the solution accuracy at a fixed in-line fold. Another way of accuracy improvement is to process the signal to improve the signal/noise ratio in the initial seismograms. Head waves have their apparent velocities different from those of other components. Thus, one can increase the signal/noise ratio, as well as the coherence spectrum (i.e., to improve the quality of seismic sections) by velocity filtering of the initial traces.
Examples of processing field head-wave data Dynamic conversion of head waves into time sections is actually inversion of coherent head waves along parallel correlations for a system of given points. A time section with a user-specified offset is a specific processing case. In fact, it can be both dynamic conversion of head waves and wavefields from a complexly configured system of stacking points. In this study we focus on the theory of head-wave processing, and the cited field examples highlight briefly the potentialities of the suggested processing technique. Application of the developed algorithms is an extensive issue to be discussed in a special publication. The conversion examples are for different data sets. We applied the algorithm to common-depth-surface (CDS) data (Fig. 4) from the Nepa dome (Siberian craton) and to CMP-DSS data. The seismograms were courtesy of people from the Irkutskgeofizika Association. The CDS system recorded head-wave first arrivals from the craton basement surface. The shotpoints were spaced at 200 m along the profile. The spread of receivers had an end-on geometry with offsets from 7200 to 11,800 m, and a spacing of 50 m between the receivers.
Fig. 4. Dynamic head-wave conversion: CDS data. A — time section, L= 10 km; B — wavefield converted into common shotpoint seismogram.
A.F. Emanov et al. / Russian Geology and Geophysics 49 (2008) 780–789
789
Fig. 5. Dynamic head-wave conversion: CMP data. A — Time section, L = 2.5 km, second refractor coverage; B — time section, L = 3.5 km, third refractor coverage.
These head-wave data were fit for conversion without additional analysis, and no coverage study was required. We obtained dynamic time sections and the respective converted wavefields. The apparent velocities varied from 5.4 to 6.7 km/s. Head waves from the basement surface were multiphase, with apparent frequencies of 9 to 14 Hz. The time-distance section (Fig. 4, A) showed prominent basement surface in first arrivals and the wavefield converted into common shotpoint seismogram (Fig. 4, B) allowed recognizing multiple reflected refractions in second arrivals. The second phase was the most intense at an offset of L = 10 km; the first phase was the most prominent in the central part of the section but was comparable with noise in the right-hand part. The wave in the right end of the profile was almost single-phase, the first and second phases dominated in the center, and the second and third phases dominated in the profile left end. Dynamic conversion of head waves is actually a technique providing selection of parallel refractions. The noise/signal ratio increases notably in low-quality fragments. Comparison of field and computed seismograms shows that conversion leaves only head waves in the records while other components disappear. Thus, processing of this kind simplifies the wave pattern and makes the interpretation much easier because wave-type identification is no longer needed. The CMP-DSS data were collected with a main and an additional systems. In the former, the in-line fold was 48, the shotpoints were spaced at 100 m, and the shot-receiver distance was from 0 to 4.8 km; in the latter the same parameters were, respectively, 24, 200 m, and 4.8–9.6 km. Unlike the CDS data, CMP-DSS data processing requires coverage constraining because of wave change. We selected three head waves from kinematic data: for the first wave the coverage was from L = 50 m to 2.5 km, the apparent velocities were 1.8–2.2 km/s, the signal was single-phase, and the visible period was ∼62−66 ms; for the second wave, the offset was L = 2.5 to 3.5 km, the apparent velocities were ∼3.8−4.2 km/s, the signal duration was ∼1.5−2 periods, and the period was
∼60−64 ms; for the third wave, the parameters were, respectively, 3.5 to 9 km, ∼5.6−5.8 km/s, and a lower frequency, i.e., a longer period of ∼46−54 ms. See the sections for the second and third head waves in Fig. 5, where the inphase axis highlights a single refractor in the former (Fig. 5, A) and inphase axes appearing in second arrivals in the latter (Fig. 5, B). Time sections can be obtained successively with each first-arrival head wave for different offsets.
References Bendat, J.S., Piersol, A.G., 1980. Engineering Applications of Correlation and Spectral Analysis. John Wiley & Sons, New York. Emanov, A.F., Seleznev, B.C., Korshik, N.A., 2001. Automated processing of multiple-coverage refracted wave data, in: Setting up a National Network of Reference Geophysical Profiles: Methods of Data Acquisition, Processing, and Interpretation [in Russian]. SNIIGGiMS, Novosibirsk, pp. 145–161. Emanov, A.F., Seleznev, B.C., Solov’ev, V.M., Larkin, G.V., Fateev, A.V., Korshik, N.A., Gritsenko, S.A., Ivanov, N.K., 1998. Automated processing of multiple-coverage head waves, in: Lithospheric Structure: Methods of Investigation and Monitoring [in Russian]. NIC OIGGM SO RAN, Novosibirsk, pp. 197–204. Epinat’eva, A.M., Goloshubin, G.M., Litvin, A.P., Pavlenkin, A.D., Petrashen’, G.I., Starobinets, A.E., Shneerson, M.B., 1990. Method of Head Waves [in Russian]. Nedra, Moscow. Gol’din, S.V., 1974. Linear Conversion of Seismic Signals [in Russian]. Nedra, Moscow. Krylov, S.V., Sergeev, V.N., 1985. Properties of head waves and new possibilities of automation of their processing. Geologiya i Geofizika (Soviet Geology and Geophysics) 26 (4), 92–102 (85–94). Mitrofanov, G.M., Sergeev, V.N., 1986. Investigation of linearized model for head wave with regard to processing CMRW data. Geologiya i Geofizika (Soviet Geology and Geophysics) 27 (8), 98–108 (92–101). Pugachev, B.C., Sinitsin, I.N., 2000. The Theory of Stochastic Systems [in Russian]. Logos, Moscow. Seleznev, B.C., Emanov, A.F., 1998. Conversion of head-wave fields by the Wiener filters. Geologiya i Geofizika (Russian Geology and Geophysics) 39 (4), 536–546 (548–559).
Editorial responsibility: M.I. Epov