- Email: [email protected]
Attenuation of multiple in reflection seismic data using Kalman–Bucy filter
Applied Mathematics and Computation 189 (2007) 805–815 www.elsevier.com/locate/amc
Attenuation of multiple in reflection seismic data using Kalman–Bucy filter Marcus P.C. Rocha
b
a,*
, Lourenildo W.B. Leite b, Mauro de L. Santos a, Valcir J.da C. Farias a
a Department of Mathematics, Federal University of Para´, Brazil Department of Geophysics, Graduate Course in Geophysics, Federal University of Para´, Brazil
Abstract The main objective of this work is the study and the application of the Kalman–Bucy method in the deconvolution process with prediction, considering the observed data as non-stationary. We propose a new prediction deconvolution operator (attenuation of multiple). The data used in this work are synthetic and, based on this, this work has characteristics of a numerical research. Ó 2007 Elsevier Inc. All rights reserved. Keywords: Stochastic process; Kalman–Bucy filter; Deconvolution; State space
1. Introduction The main objective of this paper is the study and the application of the Kalman–Bucy method in the deconvolution process with prediction considering the observed data as non-stationary. The Kalman–Bucy method is considered an alternative process in time domain in addition to the Wiener–Hopf theory. The valorization and importance of the Kalman–Bucy filter derives from its applicability and versatility as well as the essence of its conceptualization, this is a basic condition to geophysical processes. The basic application references in seismic in this paper are Crump [1], Mendel et al. [6], Mendel [4,5] and Robinson [7]. The prediction operator (KBCP) is based on Crump [7] and Mendel et al. [6] theories. Its structure resembles that of Wiener–Hopf filter, where the coefficients of the operator (WHLP) are obtained through autocorrelation, in the case of (KBCP) are obtained from the function bi ðkÞ. The objective of the seismic section of reflection is to allow the interpretation of the registered data. The detailed representation of the seismic signal demands a relatively complex model, and its processing makes use of a set of techniques based on information stochastic properties. Non-stationary stochastic processes are the basic characteristics of geophysical data and it is necessary to the application of the method used here. *
Corresponding author. E-mail address: [email protected] (M.P.C. Rocha).
0096-3003/$ - see front matter Ó 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2006.11.186
806
M.P.C. Rocha et al. / Applied Mathematics and Computation 189 (2007) 805–815
There are different techniques that can be used to get the best structural image of the undersurface. Some methods are based on the assumption that the seismogram includes only primary reflections, however, the seismogram in reality also allows multiple reflections that can be strong, making the deep reflectors totally invisible. In general, in the acquisition of maritime data, the waves, on the water surface, are repeatedly reflected on the surface of the ocean without any significant occurrence of range of loss. In order to identify and locate a target reflector (oil reservoir), it is necessary to eliminate the multiples or at least to attenuate it. In this work we propose and discuss a new method to attenuate the surface multiples. This method does not identify the multiples and it does not require measures of stationarity. The KBCP method was applied only to synthetic data of common source section obtained from models with continuums interfaces and homogeneous layers. Four different models were selected to represent the geological structure, i.e., a horizontal-plane layer, two dipping planar interface, simulation of a segment of the Solimes Basin and curved layers. As previously mentioned, this work aims specifically to develop a prediction deconvolution operator. Therefore, in this manner we supply an alternative deconvolution technique to aid in seismic sections processing. In Section 2, we present the Wiener–Kolmogorov problem, specifying the Wiener–Hopf equation and later we describe the solution presented by Kalman–Bucy in the continuous and discrete form for processes nonstationary. We propose in Section 3, a new prediction deconvolution operator (attenuation of multiple). In this section, we define the fundamental theories based on the Kalman–Bucy theory and extended by Crump [1], in the deconvolution of the pulse, and by Mendel et al. [6] to generate synthetic seismogram. The identification of the problem, presentation of the solution, fluxogram, diagram of blocks, results obtained from the synthetic data are present in Sections 3 and 4. In Section 5, we present our conclusions on the Kalman–Bucy method to the attenuation of multiple. 2. Deconvolution of non-stationary processes 2.1. The Wiener–Kolmogorov problem We initiate with the description of the simplest Wiener–Kolmogorov problem, the stationary form, whose objective it is to obtain an optimum-filter function, time-invariant, hðtÞ, that operates on the measured signal, zðtÞ and minimizes the mean quadratic error between real output ^xðtÞ and the desired output xðtÞ. The objective function of minimizing is given by 2
EfhðtÞg ¼ Ef½^xðtÞ xðtÞ g;
ð1Þ
that results in the normal equations between deviations and observations, Ef½^xðtÞ xðtÞzðtÞg ¼ 0: The basic formulation of the operation is given by the integral Z 1 ^xðtÞ ¼ zðt; sÞhðt sÞds;
ð2Þ
ð3Þ
1
where ^xðtÞ is the estimated sign, zðtÞ is the input and hðtÞ the optimum filter. It is shown [4] that the filter response to the impulse, time-invariant, satisfies the Wiener–Hopf equation: Z 1 /xz ðtÞ ¼ hw ðsÞ/zz ðt sÞds; ð4Þ 1
where /xz ðtÞ and /zz ðtÞ are the cross correlation functions and stochastic autocorrelations, and it is taken into account that, xðtÞ and vðtÞ are stationary random processes, and together with zðtÞ, hðtÞ, /xz ðtÞ and /zz ðtÞ are real, continuous, time unbound, and convergent functions. The Wiener–Hopf filter (WHF) presents a disadvantage for the seismic deconvolution which is the assumption of the stationarity of the process, and which is
M.P.C. Rocha et al. / Applied Mathematics and Computation 189 (2007) 805–815
807
not enough to solve the problem proposed in this paper, which is non-stationary. The stationarity concept used here is restrict, therefore, it is another a priori condition, which means that it is defined over the expectation and covariance as being constants in the WH case. Problem generalizations: (1) extensions to non-stationarity, hðt; sÞ; (2) operation limitation to a moving window, t 6 s 6 T and t0 6 r 6 T ; and (3) matrix equation to multichannel. These conditions are not satisfied by the integral of the convolution, with the output real, ^xðtÞ, rewritten in the form of a moving average according Wiener–Kolmogorov theory which established the relation: ^xðtÞ ¼
Z
t
hðt; sÞzðsÞds:
ð5Þ
t0
Eq. (5) is related to the correlations by /xz ðt; rÞe/zz ðs; rÞ which may be written as /xz ðt; rÞ ¼
Z
t
hðt; sÞ/zz ðs; rÞds:
ð6Þ
t0
The integral equation (6) is of the first type and represents the inherent difficulties to the general solution of deconvolution problems and geophysical inversion. On the other hand, it is useful in representing the multidimensional stochastic and non-stochastic processes where finite observations are included as well as time-variant estimates. The Kalman and Bucy [2] have converted the integral equation in ordinary linear and nonlinear differential equations which are more suitable to numerical analysis. 2.2. Kalman–Bucy solution 2.2.1. Continuous and discrete forms The problem solution ((5) and (6)) is initiated with an ordinary differential equation of order N 1 that expresses the relation between the input and output in a system, given by N 1 X n¼0
an ðtÞ
dn yðtÞ ¼ wðtÞ ða1 ¼ 1Þ: dtn
ð7Þ
The transformation for state variables xn ðtÞ and x_ n ðtÞ is done by substituting the derivatives of yðtÞ according to the rule: x1 ¼ y; x2 ¼ y_ ; x3 ¼ €y ; . . . ; xN ¼ y N 1 x_ 1 ¼ y; x2 ¼ x3 ; x_ 3 ¼ x4 ; . . . ; x_ N 1 ¼ xn ; x_ n ¼ ðan x1 þ an1 x2 þ þ a1 xn Þ þ w:
ð8Þ ð9Þ ð10Þ
The change of variables results in the state dynamic equations. The most general case, Eqs. (11) and (12), are continuous and time-variant: x_ ðtÞ ¼ F ðtÞxðtÞ þ GðtÞwðtÞ ðsystemÞ;
ð11Þ
zðtÞ ¼ H ðtÞxðtÞ þ vðtÞ ðmeasurementÞ;
ð12Þ
where the vector xðtÞ is state variable and the matrices F ðtÞ, GðtÞ and H ðtÞ are dependent of t; the vector wðtÞ is state generating (signal); and the vector zðtÞ is a the selected output from H ðtÞ; vðtÞ is the added noise to the process. In Kalman–Bucy solution it is necessary to define the a priori stochastic properties to the processes zðtÞ; wðtÞ and vðtÞ related in the autocorrelation and the stochastic cross correlation. The a priori conditions is formed by the autocorrelations of white series, and by null cross correlations in the window of definition of the governing integral equation, and they are
808
M.P.C. Rocha et al. / Applied Mathematics and Computation 189 (2007) 805–815
Table 2.1 Dynamics equations of the Kalman–Bucy filter (FKB) in the continue form
Ricatti equation
x_ ðtÞ ¼ F ðtÞxðtÞ þ GðtÞwðtÞ zðtÞ ¼ HðtÞxðtÞ þ vðtÞ Efxð0Þg ¼ 0; Efðxð0Þ ^x0 Þðxð0Þ ^x0 ÞT g ¼ P 0 x_ ðtÞ ¼ F ðtÞ^xðtÞ þ KðtÞ½zðtÞ þ HðtÞ^xðtÞ P_ ðtÞ ¼ F ðtÞP ðtÞ þ P ðtÞF T ðtÞ þ GðtÞQðtÞGT ðtÞ KðtÞRðtÞK T ðtÞ
Gain matrix
KðtÞ ¼ P ðtÞH T ðtÞV 1 ðtÞ
System Measurement Initial conditions State estimate
Table 2.2 Dynamics equations of the Kalman–Bucy filter (FKB) in the discretized form System
x_ kþ1 ¼ Uk xk þ wk
Measurement Initial conditions State estimate extrapolation Error covariance extrapolation
zðtÞ ¼ HðtÞxðtÞ þ vðtÞ Efxð0Þg ¼ 0; Efðxð0Þ ^x0 Þðxð0Þ ^x0 ÞT g ¼ P 0 x_ kþ1 ðÞ ¼ Uk xk ðþÞ P kþ1 ðÞ ¼ Uk P k ðþÞUTk þ Qk
State estimate updated Error covariance updated
^xkþ1 ðþÞ ¼ ^xkþ1 ðÞ þ K kþ1 ½zkþ1 H kþ1 ^xkþ1 ðÞ P kþ1 ðþÞ ¼ ½I K kþ1 H kþ1 P kþ1 ðÞ
EfwðtÞg ¼ 0;
/ww ðt; sÞ ¼ EfwðtÞwT ðsÞg ¼ QðtÞdðt sÞ;
ð13Þ
EfvðtÞg ¼ 0;
/vv ðt; sÞ ¼ EfvðtÞvT ðsÞg ¼ RV ðtÞdðt sÞ;
ð14Þ
T
T
/wz ðt; sÞ ¼ EfwðtÞz ðtÞg ¼ 0;
/wv ðt; sÞ ¼ EfwðtÞv ðsÞg ¼ 0;
ð15Þ
/xv ðt; sÞ ¼ EfxðtÞvT ðtÞg ¼ 0;
/wx ðt; sÞ ¼ EfwðtÞxT ðsÞg ¼ 0:
ð16Þ
dðtÞ is the distribution delta of Dirac that multiplied for QðtÞ and RðtÞ define diagonal matrices to the autocorrelation matrices of the channels as priori conditions. Reorganizing Eqs. (13)–(16), we obtain linear and nonlinear differential equations. The solution in the continuous form is summarized in Table 2.1, and Table 2.2 contains the discretized form. 3. Predictive deconvolution The fundamental theory to multiple attenuation of seismic sections (1D and 2D), based on the Kalman– Bucy theory, is presented in this section. The Kalman–Bucy method was extended by Crump [1] to the deconvolution to the impulse, and by Mendel et al. [6] to generate synthetic seismogram. The KB method allows the solution of different problems, and the objective here is to obtain an operator of prediction to attenuation of multiple and this is done by means of the representation of the state variable. Crump [1] represents the state variable as a distribution of the reflection coefficients, xðtÞ ¼ rðtÞ and in this case the recursive process to the generation of the state vector is represented by the criterion: rðkÞ ¼
L X
bCi ðkÞrðk iÞ þ wðk 1Þ
ðsystemÞ;
ð17Þ
i¼1
gðkÞ ¼ sðkÞ rðkÞ þ vðkÞ
ðoutputÞ;
ð18Þ
where vðkÞ and wðkÞ are theoretically considered as white stochastic processes. Mendel et al. [6] define the state variable as upgoing and downgoing waves from the representation: xðtÞ ¼ AxðtÞ þ bmðtÞ T
gðtÞ ¼ c xðtÞ þ r0 mðtÞ
ðsystemÞ; ðoutputÞ;
ð19Þ ð20Þ
M.P.C. Rocha et al. / Applied Mathematics and Computation 189 (2007) 805–815
where mðtÞ is the function of generation of state (signal) and A 1 A2 A¼ ; A 3 A4 3 2 0 0 0 0 0 6 ð1 þ r1 Þ 0 0 0 07 7 6 7 6 6 0 ð1 þ rk Þ 0 0 07 7 6 ; A1 ¼ 6 0 0 ð1 þ rk Þ 0 07 7 6 7 6 .. 7 .. .. .. .. 6 4 .5 . . . . 0 0 0 ð1 þ rk1 Þ 0 3 2 0 0 0 0 ð1 r1 Þ 7 60 0 ð1 r2 Þ 0 0 7 6 7 6 7 60 0 0 ð1 r3 Þ 0 7 6 A4 ¼ 6 . 7; . . . . 7 6 .. .. . .. .. . 7 6 7 6 40 0 0 0 ð1 rk1 Þ 5 0 0 0 0 0 A2 ¼ dig½ r0 A3 ¼ dig½ r1
r1 r2
r2 r3
r3 r4
rk1 ;
809
ð21Þ
ð22Þ
ð23Þ
ð24Þ
rk ;
ð25Þ
T
b ¼ ½ 1 þ r0 0 0 0 0 ; T c ¼ ½ 0 1 r0 0 0 0 ;
ð26Þ ð27Þ
and rk are the reflection coefficients. The state variable that defines the seismogram is represented by the upgoing waves, dðtÞ, and downgoing, uðtÞ as xðtÞ ¼ ½ u1 ðtÞ
d 1 ðtÞ
u2 ðtÞ
d 2 ðtÞ
uk ðtÞ
T d k ðtÞ :
ð28Þ
3.1. Identification of the problem To write the pair of system-output equations it is necessary to make the identify of the state variable with the non-stationary model. The selection of the state vector x has not uniqueness and it is defined as xðkÞ ¼ ½ gðkÞ
gðk 1Þ
gðk 2Þ
gðk L þ 1Þ ;
ð29Þ
where gðkÞ is the seismogram with noise (Eq. (18)). The dynamic equation of the recursive process of the generation of state vector will be completed with the equation: L X gðkÞ ¼ bi ðk 1Þgðk iÞ þ uðk 1Þ ðsystemÞ; ð30Þ i¼1
where uðkÞ is considered theoretically as white stochastic process. The equation of state transition is xðkÞ ¼ Uðk; k 1Þxðk 1Þ þ Suðk 1Þ ðsystemÞ; where
2
b0 ðN 1Þ 6 1 6 6 6 0 Uðk; k 1Þ6 6 0 4 0
3 bL2 ðN 1Þ bL1 ðN 1Þ bL ðN 1Þ 7 0 0 0 7 7 7: 0 0 0 7 7 1 0 0 5 0 1 0
ð31Þ
ð32Þ
810
M.P.C. Rocha et al. / Applied Mathematics and Computation 189 (2007) 805–815
Eqs. (30) and (31) are the mathematical problem to be solved, where gðkÞ is the seismic trace and bi ðkÞ is the parameter to be estimated. 3.2. Solution of the problem To obtain the values of bi ðkÞ, the equation gðkÞ ¼
L X
bi ðk 1Þgðk iÞ þ uðk 1Þ;
ð33Þ
i¼1
it is analyzed from the matrix form 2 3 2 0 0 g1 6g 7 6 0 0 6 27 6 6 7 6 6 g3 7 6 0 b1 ð2Þ 6 7¼6 6 . 7 6 .. .. .. 6 . 7 6 4 . 5 4 . . . b1 ðN 1Þ
gN
bL2 ðN 1Þ
0 b1 ð1Þ b2 ð2Þ .. .
32
3 2 gL1 7 6 6 b2 ð1Þ 7 76 ... 7 6 76 7 6 7 6 6 b3 ð2Þ 7 76 g2 7 þ 6 7 7 6 6 .. 76 7 6 5 4 g1 5 4 .
b1 ð0Þ
bL1 ðN 1Þ bL N 1
g0
u0
3
u1 7 7 7 u2 7 7: .. 7 7 . 5
ð34Þ
uN 1
It is defined that L¼
P X
Ti
ðP ¼ 1; 2; . . . ; k ¼ 1; 2; 3; . . . ; N Þ;
ð35Þ
i¼1
where L is the length of the operator bi ðkÞ; N is the length of input signal; uðk 1Þ is considered theoretically as a white stochastic process, P represents the periodicity of events to be attenuated and T is the difference among the travel time of the primary and the multiple. Eq. (28) projects the seismic trace through the average of the previous points. It is hard to obtain the solution of Eq. (29) because it has an atypical structure. To solve Eq. (29) it was necessary to include a parameter, b, in the process recursive of generation of the function bi ðkÞ, as g For k ¼ 1 b1 ð0Þ ¼ 1 ; ð36Þ g0 g ½b2 ð1Þg0 for k ¼ 2 b2 ð1Þ ¼ bb1 ð0Þ; b1 ð1Þ ¼ 2 ; ð37Þ g1 for k ¼ N b2 ðN 1Þ ¼ bb1 ðN 2Þ; bN 1 ðN 1Þ ¼ bbN 2 ðN 2Þ; ð38Þ b1 ðN 1Þ ¼
gN ½b2 ðN 1ÞgN 2 þ þ bL1 ðN 1Þg1 þ bL ðN 1Þg0 : gN 1
ð39Þ
The term uðk 1Þ in Eq. (30) has average zero and is non-correlated with any previous information, then its contribution was neglected. In each line the parameter b was considered as constant in each trace in the interval 0 < b < 1. In this way, it is possible to determine of recursive form of the values of bi ðkÞ. The prediction operator KBCP 1D and 2D are defined in a similar way to WHLP (see [8]), where the coefficients of the operator WHLP are obtained by means of the autocorrelation. In this case KBCP are obtained from the function bi ðkÞ, and they have the form presented in Figs. 1 and 2. 1D
2D
hk ðiÞ ¼ ½1; 0; . . . ; 0; bT 1 ðkÞ; 0; . . . ; 0; bT 2 ðkÞ; . . . hk ðiÞ ¼ ½1; 0; . . . ; 0; bT 1 ðkÞ; bT iþ1 ðkÞ; . . . ; 0; bT n ðkÞ; . . .
hjk ðiÞ ¼ ½1; 0; . . . ; 0; bT 1 ðxÞ ðkÞ; 0; . . . ; 0; bT 2 ðxÞ ðkÞ; . . . hjk ðiÞ ¼ ½1; 0; . . . ; 0; bT 1 ðxÞ ðkÞ; bT iþ1 ðkÞ ; . . . ; 0; bT n ðxÞ ðkÞ; . . .
The KBCP operator makes your operation in a window that includes the primary and its multiple. The reason of this strategy is to avoid distorters generated by the events out of the window.
M.P.C. Rocha et al. / Applied Mathematics and Computation 189 (2007) 805–815
1,2,3,K , n
bi T1 (k )
bi T1 (k ) bi
bk (i)
T1 T
811
T2 (k )
T2 2T
bn
i
T1 T
Tn (k )
Tn T (n 1 )
Fig. 1. Structures of the operator KBCP, bi ðkÞ, case 1D. T constant.
1
1 bi T1 (k )
b i T1 (k ) bi
hk ( j , i)
T2 (k ) n Tn (k )
T1 T ( x j , t1)
T2 T ( x j , t2 )
i
1
j 1
Tn
T (n ) 1
Fig. 2. Structures of the operator KBCP, bij ðkÞ, case 2D. T ¼ T ðx; tÞ variable.
The algorithm, organized in eight stages, describes the recursive form to obtain the KBCP operator and how to establish its initial conditions is shown as follows: 1. Input information: To choose the seismic section. 2. To identify the state variable (system) to set up Eqs. (24) and (25). 3. To define the initial conditions to start the recursive process of generation of the function bi ðkÞ: the initial attenuation factor b0, the increment D, the travel time of the primary and of its multiple. 4. Calculate the operator KBCP and the real input, y ðkÞ, to each increment Dt; b0 6 bj ¼ bj1 þ Dt 6 bp , where j ¼ 1; 2; . . . ; p. 5. To calculate the minimum error among the wanted output, yðkÞ, and the real output, y ðkÞ, to the form: EðkÞ ¼ jyðkÞ y ðkÞj;
ð40Þ
where yðkÞ is the seismogram without multiple, and y ðkÞ is the output after the KBCP operator. The attenuation factor, b, and the length of the operator ðLÞ are modified to obtain the matrix Eij ðkÞ to each y ij ðkÞ, i being the number of variations of b and j the number of variations of the length of the operator. Eij ðkÞ is a matrix, and the global minimum error is obtained through an automatic search (sub-routine Matlab). 6. To select the attenuation factor b and the length of the operator to obtain the smallest error. In this form, we obtain the best value of bi ðkÞ and, consequently, the operator KBCP more appropriate to each trace of the seismic section. 7. To obtain the convolution of the prediction operator KBCP with the trace corresponding of the seismic section. 8. Steps 3, 4 and 5 are accomplished in all seismic sections.
4. Results In this section, the results of the application of the KBCP method in synthetic data are presented. The purpose of this application is to evaluate the resolution in the sections of source-common. The synthetic data models are considered with continuous interfaces and homogeneous layers. Four models with different geological structures were selected; that is, a plan-horizontal layer, a plan-tilted layer, two plan-tilted layers and the simulation of a space of the Basin of Solimes. By means of an algorithm of ray trace the primary and multiple
812
M.P.C. Rocha et al. / Applied Mathematics and Computation 189 (2007) 805–815
reflections to first interface were generated. The group of data consists of a section of source-common with 50 receivers with uniform interval of 25 m. The sign-source is of the type Gabor with a dominant frequency of 40 Hz, and interval of sampling of two bad. 4.1. Model 1 This model is composed by a homogeneous layer separated by plan-horizontal interfaces on a half space (Fig. 3a), and velocities of 1800 m/s and 2300 m/s, and thickness 1300 m and 1200 m, respectively. With the purpose of evaluating the stability of the KBCP operator, we added white aleatory noise (Fig. 4a). The relation signal/noise (S/R) is calculated by reason of the estimate of the variances. The results obtained without noise are presented in Fig. 3b, and Fig. 4b shows the results with noise. We can notice that the KBCP operator attenuated the selected multiple in the two cases (without or with noise). 4.2. Model 2 This model is composed by two homogeneous layers separated by two planar-tilted continuous interfaces on a semi-space. The velocities of 3000 m/s, 4000 m/s and 5500 m/s were considered to each medium (Fig. 5a). The results obtained without noise are shown in Fig. 5c. We can notice that the operator KBCP attenuates the selected multiples. 4.3. Model 3 This model is formed by seven homogeneous layers with continuous, plan-tilted and horizontal interfaces (Fig. 6a). To improve the visualization of the interfaces and of the multiples (Fig. 6b), we applied dynamic gain in the geological section. The results obtained from the Model 3 are illustrated in Fig. 6c. We noticed
Fig. 3. (a) Source-common seismic section without noise and with a multiple of surface. (b) Result of the application of the operator KBCP. We observe that the multiple is attenuated.
M.P.C. Rocha et al. / Applied Mathematics and Computation 189 (2007) 805–815
813
Fig. 4. (a) Seismic section with noise (S/R = 75) with a symmetrical multiple. (b) Section after the application of the operator KBCP. We observe that the multiple is attenuate.
Fig. 5. (a) Model of two plan-tilted layers on a semi-space. The velocities were of 3000 m/s, 4000 m/s and 5500 m/s. Fifty receivers with uniform spacing of 25 m were used to realize the measurements. The rays source-receiver are illustrated. (b) Source-common seismic section obtained using the program seis88. (c) Result of the application of the operator KBCP.
that the operator KBCP attenuates the selected multiple. However, due to this model being more complex than the previous models, the precision was less than the previous ones. 4.4. Model 4 This macro-model is defined as a longitudinal geological section with five homogeneous layers on a semispace. The layers are separated by continuous interfaces (Fig. 7a). The values of velocities and thickness of the layers are indicated in Fig. 7a. We applied dynamic gain in the section with the purpose to improve the visualization of the primary events and of the multiple presents in the section (Fig. 7b). The results obtained from this model are illustrated in
814
M.P.C. Rocha et al. / Applied Mathematics and Computation 189 (2007) 805–815
Fig. 6. (a) Simulation with a selected space of the geological section (c) showing the path of the rays in the section. Fifty receivers with an uniform spacing of 25 m were used to realize the measurements, (b) source-common seismic section for the geological section using the program seis88. (c) Seismic section after the application of the operator KBCP.
Fig. 7. (a) Geological section used for simulation of the seismic section evidencing the velocities of the model. (b) Source-common seismic section for the geological section using the program seis88 with application of gain in the section. (c) Seismic section after the application of the operator KBCP.
Fig. 7c. We noticed that the operator KBCP attenuates the selected multiple, however, the precision is less than previous model. 5. Conclusions The comparison with other parallel studies shows that the implementation of FKBD and KBCP can be simpler than the one of conventional FWH, and theoretically the application of FKB neglected the measure of stationary. The versatility of FKB is directly related to the capacity of generalization of the problem WH in the sense of the non-stationary in the window of the data, that is a natural condition of geophysical data. The results obtained here show that FKBD operator solves the problem of the compression of the pulsesource. The FKBD has an ascending lineal response with the frequency, which is a disadvantage of the process. This subject (FKBD) was published partially in the Brazilian Journal of Geophysics (see [3]).
M.P.C. Rocha et al. / Applied Mathematics and Computation 189 (2007) 805–815
815
The great advantage of the FKBD operator is the non-stationary and the matrix of the effective pulsesource to include any form of pulse (minimum, maximum and mixed) with any distribution of phase. The last advantage turns the method potentially more versatile. In the synthetic models the operator KBCP carried out a good attenuation of multiple of surface. Here, it was admitted that the travel times of the multiple was known, so that the KBCP operator does not allow the identification of multiple. To attenuate multiple with the operator proposed in this paper it is necessary to associate it with a method by which the prediction of multiple can be carried out. Acknowledgements The authors wish to thank the unknown reviewers and Prof. E.A. Robinson for their suggestions and comments that were important for improving this paper. The second author wishes also to acknowledge CAPES for his abroad fellowship that made possible an important partial conduction of this work, and to the WIT Consortium for the support in the sabbatical year. References [1] N. Crump, A Kalman filter approach to the deconvolution of seismic signals, Geophysics 39 (1) (1974) 1–13. [2] R.E. Kalman, R.E. Bucy, New results in linear filtering and prediction theory, Transactions of ASME, Series D, Journal of Basic Engineering 83 (1961) 95–107. [3] L.W.B. Leite, M.P.C. Rocha, Deconvoluo de Processo Ssmico No-Estacionrio, Brazilian Journal of Geophysics 18 (1) (2000) 75–89. [4] J.M. Mendel, Optimal Seismic Deconvolution, Academic Press, New York, 1983. [5] J.M. Mendel, Maximum-Likelihood Deconvolution, A Journey into Model-Based Signal Processing, Springer-Verlag, 1990. [6] J.M. Mendel, N.E. Nahi, M. Chan, Synthetic seismogram using the state-space approach, Geophysics 44 (5) (1979) 880–895. [7] E.A. Robinson, Seismic Inversion and Deconvolution. Part B: Dual-sensor Technology, Pergamon Press, Amsterdam, 1999. [8] E.A. Robinson, H. Wold, Structural properties of stationary stochastic processes with applications, in: Proceedings of the Symposium on Time Series Analysis at Brown University, 11, 1962, pp. 170–196.
Attenuation of multiple in reflection seismic data using Kalman–Bucy filter Marcus P.C. Rocha
b
a,*
, Lourenildo W.B. Leite b, Mauro de L. Santos a, Valcir J.da C. Farias a
a Department of Mathematics, Federal University of Para´, Brazil Department of Geophysics, Graduate Course in Geophysics, Federal University of Para´, Brazil
Abstract The main objective of this work is the study and the application of the Kalman–Bucy method in the deconvolution process with prediction, considering the observed data as non-stationary. We propose a new prediction deconvolution operator (attenuation of multiple). The data used in this work are synthetic and, based on this, this work has characteristics of a numerical research. Ó 2007 Elsevier Inc. All rights reserved. Keywords: Stochastic process; Kalman–Bucy filter; Deconvolution; State space
1. Introduction The main objective of this paper is the study and the application of the Kalman–Bucy method in the deconvolution process with prediction considering the observed data as non-stationary. The Kalman–Bucy method is considered an alternative process in time domain in addition to the Wiener–Hopf theory. The valorization and importance of the Kalman–Bucy filter derives from its applicability and versatility as well as the essence of its conceptualization, this is a basic condition to geophysical processes. The basic application references in seismic in this paper are Crump [1], Mendel et al. [6], Mendel [4,5] and Robinson [7]. The prediction operator (KBCP) is based on Crump [7] and Mendel et al. [6] theories. Its structure resembles that of Wiener–Hopf filter, where the coefficients of the operator (WHLP) are obtained through autocorrelation, in the case of (KBCP) are obtained from the function bi ðkÞ. The objective of the seismic section of reflection is to allow the interpretation of the registered data. The detailed representation of the seismic signal demands a relatively complex model, and its processing makes use of a set of techniques based on information stochastic properties. Non-stationary stochastic processes are the basic characteristics of geophysical data and it is necessary to the application of the method used here. *
Corresponding author. E-mail address: [email protected] (M.P.C. Rocha).
0096-3003/$ - see front matter Ó 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2006.11.186
806
M.P.C. Rocha et al. / Applied Mathematics and Computation 189 (2007) 805–815
There are different techniques that can be used to get the best structural image of the undersurface. Some methods are based on the assumption that the seismogram includes only primary reflections, however, the seismogram in reality also allows multiple reflections that can be strong, making the deep reflectors totally invisible. In general, in the acquisition of maritime data, the waves, on the water surface, are repeatedly reflected on the surface of the ocean without any significant occurrence of range of loss. In order to identify and locate a target reflector (oil reservoir), it is necessary to eliminate the multiples or at least to attenuate it. In this work we propose and discuss a new method to attenuate the surface multiples. This method does not identify the multiples and it does not require measures of stationarity. The KBCP method was applied only to synthetic data of common source section obtained from models with continuums interfaces and homogeneous layers. Four different models were selected to represent the geological structure, i.e., a horizontal-plane layer, two dipping planar interface, simulation of a segment of the Solimes Basin and curved layers. As previously mentioned, this work aims specifically to develop a prediction deconvolution operator. Therefore, in this manner we supply an alternative deconvolution technique to aid in seismic sections processing. In Section 2, we present the Wiener–Kolmogorov problem, specifying the Wiener–Hopf equation and later we describe the solution presented by Kalman–Bucy in the continuous and discrete form for processes nonstationary. We propose in Section 3, a new prediction deconvolution operator (attenuation of multiple). In this section, we define the fundamental theories based on the Kalman–Bucy theory and extended by Crump [1], in the deconvolution of the pulse, and by Mendel et al. [6] to generate synthetic seismogram. The identification of the problem, presentation of the solution, fluxogram, diagram of blocks, results obtained from the synthetic data are present in Sections 3 and 4. In Section 5, we present our conclusions on the Kalman–Bucy method to the attenuation of multiple. 2. Deconvolution of non-stationary processes 2.1. The Wiener–Kolmogorov problem We initiate with the description of the simplest Wiener–Kolmogorov problem, the stationary form, whose objective it is to obtain an optimum-filter function, time-invariant, hðtÞ, that operates on the measured signal, zðtÞ and minimizes the mean quadratic error between real output ^xðtÞ and the desired output xðtÞ. The objective function of minimizing is given by 2
EfhðtÞg ¼ Ef½^xðtÞ xðtÞ g;
ð1Þ
that results in the normal equations between deviations and observations, Ef½^xðtÞ xðtÞzðtÞg ¼ 0: The basic formulation of the operation is given by the integral Z 1 ^xðtÞ ¼ zðt; sÞhðt sÞds;
ð2Þ
ð3Þ
1
where ^xðtÞ is the estimated sign, zðtÞ is the input and hðtÞ the optimum filter. It is shown [4] that the filter response to the impulse, time-invariant, satisfies the Wiener–Hopf equation: Z 1 /xz ðtÞ ¼ hw ðsÞ/zz ðt sÞds; ð4Þ 1
where /xz ðtÞ and /zz ðtÞ are the cross correlation functions and stochastic autocorrelations, and it is taken into account that, xðtÞ and vðtÞ are stationary random processes, and together with zðtÞ, hðtÞ, /xz ðtÞ and /zz ðtÞ are real, continuous, time unbound, and convergent functions. The Wiener–Hopf filter (WHF) presents a disadvantage for the seismic deconvolution which is the assumption of the stationarity of the process, and which is
M.P.C. Rocha et al. / Applied Mathematics and Computation 189 (2007) 805–815
807
not enough to solve the problem proposed in this paper, which is non-stationary. The stationarity concept used here is restrict, therefore, it is another a priori condition, which means that it is defined over the expectation and covariance as being constants in the WH case. Problem generalizations: (1) extensions to non-stationarity, hðt; sÞ; (2) operation limitation to a moving window, t 6 s 6 T and t0 6 r 6 T ; and (3) matrix equation to multichannel. These conditions are not satisfied by the integral of the convolution, with the output real, ^xðtÞ, rewritten in the form of a moving average according Wiener–Kolmogorov theory which established the relation: ^xðtÞ ¼
Z
t
hðt; sÞzðsÞds:
ð5Þ
t0
Eq. (5) is related to the correlations by /xz ðt; rÞe/zz ðs; rÞ which may be written as /xz ðt; rÞ ¼
Z
t
hðt; sÞ/zz ðs; rÞds:
ð6Þ
t0
The integral equation (6) is of the first type and represents the inherent difficulties to the general solution of deconvolution problems and geophysical inversion. On the other hand, it is useful in representing the multidimensional stochastic and non-stochastic processes where finite observations are included as well as time-variant estimates. The Kalman and Bucy [2] have converted the integral equation in ordinary linear and nonlinear differential equations which are more suitable to numerical analysis. 2.2. Kalman–Bucy solution 2.2.1. Continuous and discrete forms The problem solution ((5) and (6)) is initiated with an ordinary differential equation of order N 1 that expresses the relation between the input and output in a system, given by N 1 X n¼0
an ðtÞ
dn yðtÞ ¼ wðtÞ ða1 ¼ 1Þ: dtn
ð7Þ
The transformation for state variables xn ðtÞ and x_ n ðtÞ is done by substituting the derivatives of yðtÞ according to the rule: x1 ¼ y; x2 ¼ y_ ; x3 ¼ €y ; . . . ; xN ¼ y N 1 x_ 1 ¼ y; x2 ¼ x3 ; x_ 3 ¼ x4 ; . . . ; x_ N 1 ¼ xn ; x_ n ¼ ðan x1 þ an1 x2 þ þ a1 xn Þ þ w:
ð8Þ ð9Þ ð10Þ
The change of variables results in the state dynamic equations. The most general case, Eqs. (11) and (12), are continuous and time-variant: x_ ðtÞ ¼ F ðtÞxðtÞ þ GðtÞwðtÞ ðsystemÞ;
ð11Þ
zðtÞ ¼ H ðtÞxðtÞ þ vðtÞ ðmeasurementÞ;
ð12Þ
where the vector xðtÞ is state variable and the matrices F ðtÞ, GðtÞ and H ðtÞ are dependent of t; the vector wðtÞ is state generating (signal); and the vector zðtÞ is a the selected output from H ðtÞ; vðtÞ is the added noise to the process. In Kalman–Bucy solution it is necessary to define the a priori stochastic properties to the processes zðtÞ; wðtÞ and vðtÞ related in the autocorrelation and the stochastic cross correlation. The a priori conditions is formed by the autocorrelations of white series, and by null cross correlations in the window of definition of the governing integral equation, and they are
808
M.P.C. Rocha et al. / Applied Mathematics and Computation 189 (2007) 805–815
Table 2.1 Dynamics equations of the Kalman–Bucy filter (FKB) in the continue form
Ricatti equation
x_ ðtÞ ¼ F ðtÞxðtÞ þ GðtÞwðtÞ zðtÞ ¼ HðtÞxðtÞ þ vðtÞ Efxð0Þg ¼ 0; Efðxð0Þ ^x0 Þðxð0Þ ^x0 ÞT g ¼ P 0 x_ ðtÞ ¼ F ðtÞ^xðtÞ þ KðtÞ½zðtÞ þ HðtÞ^xðtÞ P_ ðtÞ ¼ F ðtÞP ðtÞ þ P ðtÞF T ðtÞ þ GðtÞQðtÞGT ðtÞ KðtÞRðtÞK T ðtÞ
Gain matrix
KðtÞ ¼ P ðtÞH T ðtÞV 1 ðtÞ
System Measurement Initial conditions State estimate
Table 2.2 Dynamics equations of the Kalman–Bucy filter (FKB) in the discretized form System
x_ kþ1 ¼ Uk xk þ wk
Measurement Initial conditions State estimate extrapolation Error covariance extrapolation
zðtÞ ¼ HðtÞxðtÞ þ vðtÞ Efxð0Þg ¼ 0; Efðxð0Þ ^x0 Þðxð0Þ ^x0 ÞT g ¼ P 0 x_ kþ1 ðÞ ¼ Uk xk ðþÞ P kþ1 ðÞ ¼ Uk P k ðþÞUTk þ Qk
State estimate updated Error covariance updated
^xkþ1 ðþÞ ¼ ^xkþ1 ðÞ þ K kþ1 ½zkþ1 H kþ1 ^xkþ1 ðÞ P kþ1 ðþÞ ¼ ½I K kþ1 H kþ1 P kþ1 ðÞ
EfwðtÞg ¼ 0;
/ww ðt; sÞ ¼ EfwðtÞwT ðsÞg ¼ QðtÞdðt sÞ;
ð13Þ
EfvðtÞg ¼ 0;
/vv ðt; sÞ ¼ EfvðtÞvT ðsÞg ¼ RV ðtÞdðt sÞ;
ð14Þ
T
T
/wz ðt; sÞ ¼ EfwðtÞz ðtÞg ¼ 0;
/wv ðt; sÞ ¼ EfwðtÞv ðsÞg ¼ 0;
ð15Þ
/xv ðt; sÞ ¼ EfxðtÞvT ðtÞg ¼ 0;
/wx ðt; sÞ ¼ EfwðtÞxT ðsÞg ¼ 0:
ð16Þ
dðtÞ is the distribution delta of Dirac that multiplied for QðtÞ and RðtÞ define diagonal matrices to the autocorrelation matrices of the channels as priori conditions. Reorganizing Eqs. (13)–(16), we obtain linear and nonlinear differential equations. The solution in the continuous form is summarized in Table 2.1, and Table 2.2 contains the discretized form. 3. Predictive deconvolution The fundamental theory to multiple attenuation of seismic sections (1D and 2D), based on the Kalman– Bucy theory, is presented in this section. The Kalman–Bucy method was extended by Crump [1] to the deconvolution to the impulse, and by Mendel et al. [6] to generate synthetic seismogram. The KB method allows the solution of different problems, and the objective here is to obtain an operator of prediction to attenuation of multiple and this is done by means of the representation of the state variable. Crump [1] represents the state variable as a distribution of the reflection coefficients, xðtÞ ¼ rðtÞ and in this case the recursive process to the generation of the state vector is represented by the criterion: rðkÞ ¼
L X
bCi ðkÞrðk iÞ þ wðk 1Þ
ðsystemÞ;
ð17Þ
i¼1
gðkÞ ¼ sðkÞ rðkÞ þ vðkÞ
ðoutputÞ;
ð18Þ
where vðkÞ and wðkÞ are theoretically considered as white stochastic processes. Mendel et al. [6] define the state variable as upgoing and downgoing waves from the representation: xðtÞ ¼ AxðtÞ þ bmðtÞ T
gðtÞ ¼ c xðtÞ þ r0 mðtÞ
ðsystemÞ; ðoutputÞ;
ð19Þ ð20Þ
M.P.C. Rocha et al. / Applied Mathematics and Computation 189 (2007) 805–815
where mðtÞ is the function of generation of state (signal) and A 1 A2 A¼ ; A 3 A4 3 2 0 0 0 0 0 6 ð1 þ r1 Þ 0 0 0 07 7 6 7 6 6 0 ð1 þ rk Þ 0 0 07 7 6 ; A1 ¼ 6 0 0 ð1 þ rk Þ 0 07 7 6 7 6 .. 7 .. .. .. .. 6 4 .5 . . . . 0 0 0 ð1 þ rk1 Þ 0 3 2 0 0 0 0 ð1 r1 Þ 7 60 0 ð1 r2 Þ 0 0 7 6 7 6 7 60 0 0 ð1 r3 Þ 0 7 6 A4 ¼ 6 . 7; . . . . 7 6 .. .. . .. .. . 7 6 7 6 40 0 0 0 ð1 rk1 Þ 5 0 0 0 0 0 A2 ¼ dig½ r0 A3 ¼ dig½ r1
r1 r2
r2 r3
r3 r4
rk1 ;
809
ð21Þ
ð22Þ
ð23Þ
ð24Þ
rk ;
ð25Þ
T
b ¼ ½ 1 þ r0 0 0 0 0 ; T c ¼ ½ 0 1 r0 0 0 0 ;
ð26Þ ð27Þ
and rk are the reflection coefficients. The state variable that defines the seismogram is represented by the upgoing waves, dðtÞ, and downgoing, uðtÞ as xðtÞ ¼ ½ u1 ðtÞ
d 1 ðtÞ
u2 ðtÞ
d 2 ðtÞ
uk ðtÞ
T d k ðtÞ :
ð28Þ
3.1. Identification of the problem To write the pair of system-output equations it is necessary to make the identify of the state variable with the non-stationary model. The selection of the state vector x has not uniqueness and it is defined as xðkÞ ¼ ½ gðkÞ
gðk 1Þ
gðk 2Þ
gðk L þ 1Þ ;
ð29Þ
where gðkÞ is the seismogram with noise (Eq. (18)). The dynamic equation of the recursive process of the generation of state vector will be completed with the equation: L X gðkÞ ¼ bi ðk 1Þgðk iÞ þ uðk 1Þ ðsystemÞ; ð30Þ i¼1
where uðkÞ is considered theoretically as white stochastic process. The equation of state transition is xðkÞ ¼ Uðk; k 1Þxðk 1Þ þ Suðk 1Þ ðsystemÞ; where
2
b0 ðN 1Þ 6 1 6 6 6 0 Uðk; k 1Þ6 6 0 4 0
3 bL2 ðN 1Þ bL1 ðN 1Þ bL ðN 1Þ 7 0 0 0 7 7 7: 0 0 0 7 7 1 0 0 5 0 1 0
ð31Þ
ð32Þ
810
M.P.C. Rocha et al. / Applied Mathematics and Computation 189 (2007) 805–815
Eqs. (30) and (31) are the mathematical problem to be solved, where gðkÞ is the seismic trace and bi ðkÞ is the parameter to be estimated. 3.2. Solution of the problem To obtain the values of bi ðkÞ, the equation gðkÞ ¼
L X
bi ðk 1Þgðk iÞ þ uðk 1Þ;
ð33Þ
i¼1
it is analyzed from the matrix form 2 3 2 0 0 g1 6g 7 6 0 0 6 27 6 6 7 6 6 g3 7 6 0 b1 ð2Þ 6 7¼6 6 . 7 6 .. .. .. 6 . 7 6 4 . 5 4 . . . b1 ðN 1Þ
gN
bL2 ðN 1Þ
0 b1 ð1Þ b2 ð2Þ .. .
32
3 2 gL1 7 6 6 b2 ð1Þ 7 76 ... 7 6 76 7 6 7 6 6 b3 ð2Þ 7 76 g2 7 þ 6 7 7 6 6 .. 76 7 6 5 4 g1 5 4 .
b1 ð0Þ
bL1 ðN 1Þ bL N 1
g0
u0
3
u1 7 7 7 u2 7 7: .. 7 7 . 5
ð34Þ
uN 1
It is defined that L¼
P X
Ti
ðP ¼ 1; 2; . . . ; k ¼ 1; 2; 3; . . . ; N Þ;
ð35Þ
i¼1
where L is the length of the operator bi ðkÞ; N is the length of input signal; uðk 1Þ is considered theoretically as a white stochastic process, P represents the periodicity of events to be attenuated and T is the difference among the travel time of the primary and the multiple. Eq. (28) projects the seismic trace through the average of the previous points. It is hard to obtain the solution of Eq. (29) because it has an atypical structure. To solve Eq. (29) it was necessary to include a parameter, b, in the process recursive of generation of the function bi ðkÞ, as g For k ¼ 1 b1 ð0Þ ¼ 1 ; ð36Þ g0 g ½b2 ð1Þg0 for k ¼ 2 b2 ð1Þ ¼ bb1 ð0Þ; b1 ð1Þ ¼ 2 ; ð37Þ g1 for k ¼ N b2 ðN 1Þ ¼ bb1 ðN 2Þ; bN 1 ðN 1Þ ¼ bbN 2 ðN 2Þ; ð38Þ b1 ðN 1Þ ¼
gN ½b2 ðN 1ÞgN 2 þ þ bL1 ðN 1Þg1 þ bL ðN 1Þg0 : gN 1
ð39Þ
The term uðk 1Þ in Eq. (30) has average zero and is non-correlated with any previous information, then its contribution was neglected. In each line the parameter b was considered as constant in each trace in the interval 0 < b < 1. In this way, it is possible to determine of recursive form of the values of bi ðkÞ. The prediction operator KBCP 1D and 2D are defined in a similar way to WHLP (see [8]), where the coefficients of the operator WHLP are obtained by means of the autocorrelation. In this case KBCP are obtained from the function bi ðkÞ, and they have the form presented in Figs. 1 and 2. 1D
2D
hk ðiÞ ¼ ½1; 0; . . . ; 0; bT 1 ðkÞ; 0; . . . ; 0; bT 2 ðkÞ; . . . hk ðiÞ ¼ ½1; 0; . . . ; 0; bT 1 ðkÞ; bT iþ1 ðkÞ; . . . ; 0; bT n ðkÞ; . . .
hjk ðiÞ ¼ ½1; 0; . . . ; 0; bT 1 ðxÞ ðkÞ; 0; . . . ; 0; bT 2 ðxÞ ðkÞ; . . . hjk ðiÞ ¼ ½1; 0; . . . ; 0; bT 1 ðxÞ ðkÞ; bT iþ1 ðkÞ ; . . . ; 0; bT n ðxÞ ðkÞ; . . .
The KBCP operator makes your operation in a window that includes the primary and its multiple. The reason of this strategy is to avoid distorters generated by the events out of the window.
M.P.C. Rocha et al. / Applied Mathematics and Computation 189 (2007) 805–815
1,2,3,K , n
bi T1 (k )
bi T1 (k ) bi
bk (i)
T1 T
811
T2 (k )
T2 2T
bn
i
T1 T
Tn (k )
Tn T (n 1 )
Fig. 1. Structures of the operator KBCP, bi ðkÞ, case 1D. T constant.
1
1 bi T1 (k )
b i T1 (k ) bi
hk ( j , i)
T2 (k ) n Tn (k )
T1 T ( x j , t1)
T2 T ( x j , t2 )
i
1
j 1
Tn
T (n ) 1
Fig. 2. Structures of the operator KBCP, bij ðkÞ, case 2D. T ¼ T ðx; tÞ variable.
The algorithm, organized in eight stages, describes the recursive form to obtain the KBCP operator and how to establish its initial conditions is shown as follows: 1. Input information: To choose the seismic section. 2. To identify the state variable (system) to set up Eqs. (24) and (25). 3. To define the initial conditions to start the recursive process of generation of the function bi ðkÞ: the initial attenuation factor b0, the increment D, the travel time of the primary and of its multiple. 4. Calculate the operator KBCP and the real input, y ðkÞ, to each increment Dt; b0 6 bj ¼ bj1 þ Dt 6 bp , where j ¼ 1; 2; . . . ; p. 5. To calculate the minimum error among the wanted output, yðkÞ, and the real output, y ðkÞ, to the form: EðkÞ ¼ jyðkÞ y ðkÞj;
ð40Þ
where yðkÞ is the seismogram without multiple, and y ðkÞ is the output after the KBCP operator. The attenuation factor, b, and the length of the operator ðLÞ are modified to obtain the matrix Eij ðkÞ to each y ij ðkÞ, i being the number of variations of b and j the number of variations of the length of the operator. Eij ðkÞ is a matrix, and the global minimum error is obtained through an automatic search (sub-routine Matlab). 6. To select the attenuation factor b and the length of the operator to obtain the smallest error. In this form, we obtain the best value of bi ðkÞ and, consequently, the operator KBCP more appropriate to each trace of the seismic section. 7. To obtain the convolution of the prediction operator KBCP with the trace corresponding of the seismic section. 8. Steps 3, 4 and 5 are accomplished in all seismic sections.
4. Results In this section, the results of the application of the KBCP method in synthetic data are presented. The purpose of this application is to evaluate the resolution in the sections of source-common. The synthetic data models are considered with continuous interfaces and homogeneous layers. Four models with different geological structures were selected; that is, a plan-horizontal layer, a plan-tilted layer, two plan-tilted layers and the simulation of a space of the Basin of Solimes. By means of an algorithm of ray trace the primary and multiple
812
M.P.C. Rocha et al. / Applied Mathematics and Computation 189 (2007) 805–815
reflections to first interface were generated. The group of data consists of a section of source-common with 50 receivers with uniform interval of 25 m. The sign-source is of the type Gabor with a dominant frequency of 40 Hz, and interval of sampling of two bad. 4.1. Model 1 This model is composed by a homogeneous layer separated by plan-horizontal interfaces on a half space (Fig. 3a), and velocities of 1800 m/s and 2300 m/s, and thickness 1300 m and 1200 m, respectively. With the purpose of evaluating the stability of the KBCP operator, we added white aleatory noise (Fig. 4a). The relation signal/noise (S/R) is calculated by reason of the estimate of the variances. The results obtained without noise are presented in Fig. 3b, and Fig. 4b shows the results with noise. We can notice that the KBCP operator attenuated the selected multiple in the two cases (without or with noise). 4.2. Model 2 This model is composed by two homogeneous layers separated by two planar-tilted continuous interfaces on a semi-space. The velocities of 3000 m/s, 4000 m/s and 5500 m/s were considered to each medium (Fig. 5a). The results obtained without noise are shown in Fig. 5c. We can notice that the operator KBCP attenuates the selected multiples. 4.3. Model 3 This model is formed by seven homogeneous layers with continuous, plan-tilted and horizontal interfaces (Fig. 6a). To improve the visualization of the interfaces and of the multiples (Fig. 6b), we applied dynamic gain in the geological section. The results obtained from the Model 3 are illustrated in Fig. 6c. We noticed
Fig. 3. (a) Source-common seismic section without noise and with a multiple of surface. (b) Result of the application of the operator KBCP. We observe that the multiple is attenuated.
M.P.C. Rocha et al. / Applied Mathematics and Computation 189 (2007) 805–815
813
Fig. 4. (a) Seismic section with noise (S/R = 75) with a symmetrical multiple. (b) Section after the application of the operator KBCP. We observe that the multiple is attenuate.
Fig. 5. (a) Model of two plan-tilted layers on a semi-space. The velocities were of 3000 m/s, 4000 m/s and 5500 m/s. Fifty receivers with uniform spacing of 25 m were used to realize the measurements. The rays source-receiver are illustrated. (b) Source-common seismic section obtained using the program seis88. (c) Result of the application of the operator KBCP.
that the operator KBCP attenuates the selected multiple. However, due to this model being more complex than the previous models, the precision was less than the previous ones. 4.4. Model 4 This macro-model is defined as a longitudinal geological section with five homogeneous layers on a semispace. The layers are separated by continuous interfaces (Fig. 7a). The values of velocities and thickness of the layers are indicated in Fig. 7a. We applied dynamic gain in the section with the purpose to improve the visualization of the primary events and of the multiple presents in the section (Fig. 7b). The results obtained from this model are illustrated in
814
M.P.C. Rocha et al. / Applied Mathematics and Computation 189 (2007) 805–815
Fig. 6. (a) Simulation with a selected space of the geological section (c) showing the path of the rays in the section. Fifty receivers with an uniform spacing of 25 m were used to realize the measurements, (b) source-common seismic section for the geological section using the program seis88. (c) Seismic section after the application of the operator KBCP.
Fig. 7. (a) Geological section used for simulation of the seismic section evidencing the velocities of the model. (b) Source-common seismic section for the geological section using the program seis88 with application of gain in the section. (c) Seismic section after the application of the operator KBCP.
Fig. 7c. We noticed that the operator KBCP attenuates the selected multiple, however, the precision is less than previous model. 5. Conclusions The comparison with other parallel studies shows that the implementation of FKBD and KBCP can be simpler than the one of conventional FWH, and theoretically the application of FKB neglected the measure of stationary. The versatility of FKB is directly related to the capacity of generalization of the problem WH in the sense of the non-stationary in the window of the data, that is a natural condition of geophysical data. The results obtained here show that FKBD operator solves the problem of the compression of the pulsesource. The FKBD has an ascending lineal response with the frequency, which is a disadvantage of the process. This subject (FKBD) was published partially in the Brazilian Journal of Geophysics (see [3]).
M.P.C. Rocha et al. / Applied Mathematics and Computation 189 (2007) 805–815
815
The great advantage of the FKBD operator is the non-stationary and the matrix of the effective pulsesource to include any form of pulse (minimum, maximum and mixed) with any distribution of phase. The last advantage turns the method potentially more versatile. In the synthetic models the operator KBCP carried out a good attenuation of multiple of surface. Here, it was admitted that the travel times of the multiple was known, so that the KBCP operator does not allow the identification of multiple. To attenuate multiple with the operator proposed in this paper it is necessary to associate it with a method by which the prediction of multiple can be carried out. Acknowledgements The authors wish to thank the unknown reviewers and Prof. E.A. Robinson for their suggestions and comments that were important for improving this paper. The second author wishes also to acknowledge CAPES for his abroad fellowship that made possible an important partial conduction of this work, and to the WIT Consortium for the support in the sabbatical year. References [1] N. Crump, A Kalman filter approach to the deconvolution of seismic signals, Geophysics 39 (1) (1974) 1–13. [2] R.E. Kalman, R.E. Bucy, New results in linear filtering and prediction theory, Transactions of ASME, Series D, Journal of Basic Engineering 83 (1961) 95–107. [3] L.W.B. Leite, M.P.C. Rocha, Deconvoluo de Processo Ssmico No-Estacionrio, Brazilian Journal of Geophysics 18 (1) (2000) 75–89. [4] J.M. Mendel, Optimal Seismic Deconvolution, Academic Press, New York, 1983. [5] J.M. Mendel, Maximum-Likelihood Deconvolution, A Journey into Model-Based Signal Processing, Springer-Verlag, 1990. [6] J.M. Mendel, N.E. Nahi, M. Chan, Synthetic seismogram using the state-space approach, Geophysics 44 (5) (1979) 880–895. [7] E.A. Robinson, Seismic Inversion and Deconvolution. Part B: Dual-sensor Technology, Pergamon Press, Amsterdam, 1999. [8] E.A. Robinson, H. Wold, Structural properties of stationary stochastic processes with applications, in: Proceedings of the Symposium on Time Series Analysis at Brown University, 11, 1962, pp. 170–196.