(China Petroleum Exploration and Development Research Institute, Beijing 100083)

The purpose of seismic exploration is to obtain accurate imaging of underground structures. Due to human factors and environmental reasons, seismic data are often irregularly sampled or missing in the spatial direction, so it is often necessary to reconstruct the missing seismic data in the spatial direction. The minimum norm Fourier reconstruction method is based on estimating the Fourier coefficients of irregularly sampled seismic data. Once these coefficients are accurately obtained, seismic data can be reconstructed to any suitable spatial position by inverse Fourier transform. The main advantage of this method is that it can deal with both regular sampling data with empty tracks and irregular sampling data. The disadvantage of this method is that it is impossible to reconstruct seismic data with spatial aliasing and excessive gaps. Aiming at the problem of seismic data reconstruction with spatial aliasing, the minimum norm Fourier reconstruction method is combined with multi-step autoregressive method to overcome the shortcomings of the minimum norm Fourier reconstruction method. Through the verification of different theories and actual earthquake data examples, it shows that the reconstruction method is effective and practical.

Seismic data reconstruction; Minimum norm inversion; Fourier transform; Multistep autoregressive

Study on reconstruction method of prestack seismic data

Huo Zhiming

(China Petrochemical Exploration and Production Research Institute, Beijing 100083)

The goal of exploration seismology is to obtain accurate images of the underground. Due to the influence of human factors and environmental conditions, irregular sampling or missing sampling often occurs in the spatial direction of seismic data. Therefore, it is often necessary to reconstruct the missing seismic data along the spatial direction. Fourier reconstruction with minimum normalized transform is based on the estimation of Fourier coefficients describing irregularly sampled seismic data, and once these coefficients are obtained, seismic data can be reconstructed at any suitable spatial position by inverse Fourier transform. The main advantage of Fourier reconstruction is flexibility, because it can handle not only regular sampling data with gaps, but also irregular sampling data. The disadvantage of this method is that it can't deal with seismic data with spatial aliasing and seismic data with large gaps. In this paper, the minimum normalized Fourier reconstruction and multi-step autoregressive method are combined to solve the problem of spatial aliasing seismic data reconstruction. This method overcomes the shortcomings of Fourier reconstruction method. The multi-step autoregressive method is used to reconstruct several different theoretical and practical seismic data, which proves the effectiveness and practicability of this method.

Keywords seismic data reconstruction; Minimum norm inversion; Fourier transform; Multistep autoregressive

As we all know, the acquisition of seismic data seriously affects the final imaging results of seismic data, and a common problem in seismic data acquisition is that seismic data is irregularly sampled or diluted along the spatial direction. The reason of sparse sampling of seismic data in spatial direction is mainly due to economic factors. Sparse sampling is more economical, but it means less data is collected, and it will lead to spatial aliasing in seismic data, especially in 3D seismic exploration. The main reason for irregular sampling of seismic data in spatial direction is the existence of surface obstacles (buildings, roads, bridges, etc.). ) or topographic conditions (forbidden mining areas and mountainous areas, forests, river network areas, etc.). ), bad acquisition trajectory caused by instrument hardware (geophone, air gun, cable, etc.). ), and the feathering drift of the cable during marine seismic data acquisition. In the process of seismic data processing, irregular sampling and sparse sampling will not only cause human error, but also have a serious impact on the processing results of DMO, FK domain filtering, velocity analysis, multiple attenuation, spectrum estimation and wave equation migration imaging based on multi-channel technology. Therefore, by reconstructing the original seismic data, the geophysical information contained in it can reflect the geophysical characteristics of underground geological bodies more truly. Subsequent seismic data processing can better meet the requirements of detailed description of complex geological structures and provide more effective guidance and help for oil and gas exploration, which has important practical significance [1, 2].

The reconstruction method of seismic data based on Fourier transform does not need geological or geophysical assumptions, but only requires limited bandwidth of seismic data and high computational efficiency. Fourier reconstruction method estimates the Fourier coefficient of irregular sampling data through least square inversion, and how to estimate the Fourier coefficient better is the core of this method. Once the Fourier coefficients are correctly estimated, the data can be reconstructed on any sampling grid. Duijndam et al [3] applied Fourier reconstruction method to regularization of irregularly sampled seismic data, and successfully solved a series of problems such as parameter selection. Hindriks and Duijndam [4] extended this method to 3D seismic data reconstruction. Liu and Sachhi [5] proposed a Fourier reconstruction method with minimum weighted norm interpolation. This band-limited reconstruction method uses the regularization term of the adaptive spectral weighting norm to constrain the solution of the inverse equation, and takes the bandwidth of the data and the shape of the spectrum as the prior information of the band-limited seismic data reconstruction problem, so it obtains a better solution than the traditional Fourier reconstruction method of band-limited data, but does not give a good anti-aliasing method. Zwartjes and Sachhi [6] proposed a sparse constrained Fourier reconstruction method using non-quadratic regularization term to improve the reconstruction effect of wide channel seismic data and solve the reconstruction problem of spatial aliasing seismic data. Fourier reconstruction method can reconstruct not only regularly sampled seismic data, but also irregularly and randomly sampled seismic data, but it can't reconstruct spatially overlapped seismic data well.

In this paper, the reconstruction method of Fourier seismic data based on minimum norm solution is studied and analyzed, and the coefficients in Fourier domain are obtained by least square inversion method to reconstruct seismic data. In order to improve the shortcomings that the minimum norm Fourier reconstruction method can not reconstruct seismic data with too large spatial interval and has spatial aliasing, this paper adopts the idea of combining the minimum norm Fourier reconstruction method with the multi-step autoregressive method to reconstruct seismic data. This method can not only reconstruct seismic data with large spatial interval, but also reconstruct seismic data with spatial aliasing.

1 minimum norm Fourier reconstruction method

Fourier reconstruction is a method to recover signals from irregularly sampled data, which is based on sampling theorem, that is, band-limited continuous signals can be recovered from regularly sampled data. If the average sampling rate of the irregularly sampled signal exceeds the Nyquist sampling rate, the irregularly sampled signal can also be reconstructed. In the case of regular sampling, the discrete Fourier transform is an orthogonal transform. However, when the sampling is irregular, the basis functions of Fourier transform are no longer orthogonal, which means that directly using discrete Fourier transform to calculate Fourier coefficients will produce errors. It is a remedial measure to calculate Fourier coefficient by least square inversion [7].

Assuming that the data is irregularly sampled in the spatial direction, the position of each sampling point is [x0, …, xn, …, xn- 1] respectively. Using the midpoint rule of actual sampling position and sampling interval, the discrete Fourier transform of irregular sampling data can be expressed as the following discrete summation form:

Oil and gas accumulation theory and exploration and development technology (5)

The above formula is non-uniform discrete Fourier transform. Where the spatial sampling interval △xn is defined as:

Oil and gas accumulation theory and exploration and development technology (5)

Regular sampling in wavenumber domain means that data is periodic in spatial domain, so x is the length of irregularly sampled data. If NDFT (Non-uniform Discrete Fourier Transform) is directly used to calculate the wave number, the irregular sampling will cause great errors, so the least square inversion is usually used to calculate the wave number in practical calculation.

Firstly, the mathematical transformation of sampling data at any spatial position calculated by regular sampling wave number is defined and regarded as a forward model. Assuming that the bandwidth of band-limited data in wavenumber domain is [-m △k, m △k], the sampling in wavenumber domain is regular, and △ k is the sampling interval of spatial wavenumber, the inverse discrete Fourier transform of any spatial position xn reconstructed from wavenumber domain is as follows.

Oil and gas accumulation theory and exploration and development technology (5)

Remember that the coefficient matrix is irregular, the sampling data is DN = P (xn, ω), and the required regular wave number is

Oil and gas accumulation theory and exploration and development technology (5)

Then, the formula (3) is written in matrix form, which is used as the oil and gas accumulation theory and exploration and development technology (5).

In the actual seismic data processing, because the data may not be completely band-limited, some spatial wave number components will exceed the defined frequency band range, and these excess components constitute the error and noise of the forward model, so the noise term needs to be used in the above formula:

Oil and gas accumulation theory and exploration and development technology (5)

Duijndam et al. [3] get the spatial wave number of irregular sampling data d(xn, t) through least square inversion estimation, and calculate the unknown Fourier coefficient vector of regular sampling from the irregular sampling data vector D, which can be reduced to solving an ill-posed linear inversion problem, and it needs regularization and constructing an appropriate solution by using some prior information. Any desired parameter estimation technique can be used. First, we assume that noise n = n (0, Cn) and prior information.

Oil and gas accumulation theory and exploration and development technology (5)

The covariance matrix of noise is Cn, and its average value is zero. Oil and gas accumulation theory and exploration and development technology for finding posterior probability density function by Bayesian parameter inversion method (5)

Where is the likelihood function, representing the prior distribution of model vectors. Satisfy separately

Oil and gas accumulation theory and exploration and development technology (5)

Oil and gas accumulation theory and exploration and development technology (5)

The maximum posterior probability solution is transformed into the minimum solution of the subsequent objective function, and the objective function is established.

Oil and gas accumulation theory and exploration and development technology (5)

Minimize the objective function:

Oil and gas accumulation theory and exploration and development technology (5)

Here, in order to calculate the regular sampling wavenumber to be obtained, AH is the * * * yoke transpose matrix of matrix A and the covariance matrix of the prior model.

Let's simplify the following equation (9). First of all, for seismic data, there is usually no prior model information, so there is generally no reason to assume the correlation between spatial wave numbers, so it is a diagonal matrix, usually in the form of prior model variance. It is unrealistic to accurately express the covariance matrix Cn of noise, because the detailed information about noise is unknown. The noise covariance matrix given by Duijndam et al. [3] is cn = C2W- 1, and c is a constant; W is the diagonal matrix of weight coefficient, that is, w = diag (△ xn). According to the theory of discrete Fourier transform, we should choose △k≤2π/X, where x = ∑ n △ xn, which is the length of data, that is, x = xn- 1-x0, then equation (9) becomes.

Oil and gas accumulation theory and exploration and development technology (5)

It is called the damping factor here. λ can be determined by L curve or generalized cross-validation (GCV), and the best selection method is [4]:

Oil and gas accumulation theory and exploration and development technology (5)

Where: f is a constant given by the user, indicating the expected data signal-to-noise ratio. However, in the process of actual seismic data reconstruction, λ is generally 1% of the main diagonal element of AHWA matrix.

The solution of equation (10) is called the minimum norm solution, which is also called the damped least squares solution. This reconstruction method is called the minimum norm Fourier reconstruction (FRMN) with minimum norm inversion [8]. Generally, when sampling irregularly, the coefficient matrix AHWA of formula (10) is a ill-conditioned toplitz matrix. When the matrix W is unweighted, the ill-conditioned degree of Toeplitz matrix formed by AHA is controlled by the density between irregularly sampled data. The closer the seismic trace is in the irregularly sampled seismic data, the smaller the distance △x is, and the greater the condition number of Toplitz matrix is, the more difficult it is to solve. After adding the weight coefficient matrix W, the ill-conditioned degree of the Toplitz matrix formed by AHWA is controlled by the maximum gap △xa between data, △ xa = max (△ xn). The relationship between the condition number of the coefficient matrix AHWA and the maximum gap △xa is as follows [7]:

Oil and gas accumulation theory and exploration and development technology (5)

It can be seen from the above formula that the larger the maximum gap △xa is, the more ill-conditioned the matrix AHWA is, and the more difficult it is to converge when solving the equation. If the spatial Nyquist sampling interval is defined as

Oil and gas accumulation theory and exploration and development technology (5)

Then when △xa≥3△xNyq, the coefficient matrix AHWA can no longer guarantee iterative convergence [3]. That is to say, when the gap between irregularly sampled seismic data is too large, satisfactory reconstruction results cannot be obtained. This is the inherent disadvantage of Fourier reconstruction method.

Equation (10) is usually solved frequency by frequency in the frequency domain. When solving the equation, because the low-frequency part can completely reconstruct the data with only a small wave number bandwidth, the scale of solving the equation (10) is small and the solution is relatively easy. However, the high frequency part needs a large wavenumber bandwidth, so there are many unknowns in the solution formula (10), and the calculation time is long, so the solution is unstable. Therefore, the low-frequency part of seismic data reconstructed by minimum norm Fourier method has high accuracy.

2 Multi-step autoregressive method

Autoregressive model (predictive filter) is widely used in the field of signal processing, and it is a technology to simulate signal evolution [9]. Autoregressive model can be applied to signal prediction and denoising [10], seismic trace interpolation [1 1, 12] and parameter spectrum analysis [13]. The linear in-phase transformation from T-X domain to F-X domain is a complex sine function, which can be simulated by autoregressive operators. Spitz [1 1] and Porsani [12] proposed an autoregressive reconstruction method, which successfully solved the interpolation problem of regularly sampled seismic data with spatial aliasing. These methods use low-frequency information to recover the high-frequency part of data. However, this method is only applicable to the case that the original seismic data are sampled regularly in space, and can only be used for encryption interpolation.

Multi-step autoregressive (MSAR) [14] is an extension of Spitz's one-step prediction method, which extends its application range from only track encryption interpolation to interpolation reconstruction of irregularly missing seismic data. It is assumed that the seismic data contains a limited number of linear in-phase axes and consists of n equally spaced seismic traces, and some seismic traces are missing. Firstly, seismic data is converted from time domain to frequency domain. In the F-X domain, seismic data can be represented by vectors x(f), XT (f) = [x 1 (f), x2(f), x3(f), …, xn (f)], in which only the data of m channels are known. The known data are n = {n (1), n(2), n(3), …, n (m)} and m = {m (1), m(2), m(3), …, m (n-m)} respectively.

Seismic data consisting of l approximately linear in-phase axes can be expressed as follows in F-X domain.

Oil and gas accumulation theory and exploration and development technology (5)

Where: △x and △f respectively represent the sampling interval in the spatial domain and the frequency domain; Pj represents the slope of the j-th linear in-phase axis; Aj stands for amplitude. For each frequency component f, the above formula shows that each linear in-phase axis in the F-X domain can be represented by a complex harmonic function. When △ x ′ = α△ x and △ f ′ = △ f/α, we get:

Oil and gas accumulation theory and exploration and development technology (5)

In addition, in the form of autoregressive model, the superposition of L harmonic functions can be expressed as

Oil and gas accumulation theory and exploration and development technology (5)

Where P(j, n△f) represents the prediction filter coefficient. Similarly, for △x' and △f', there are

Oil and gas accumulation theory and exploration and development technology (5)

Comparing the expressions (15), (16) and (17), we can get:

Oil and gas accumulation theory and exploration and development technology (5)

This formula is the basis of multi-step autoregressive method. It shows that each component of the prediction filter is predictable on the frequency axis. This means that if the prediction filters of some frequencies are known, the prediction filters of other frequencies can be predicted. That is to say, we can extract the prediction filter of high frequency component from the low frequency component prediction filter without spatial aliasing reconstructed by Fourier method, and then reconstruct the high frequency component of missing seismic trace.

Assuming that the frequency range of low-frequency data reconstructed by the minimum norm Fourier method is f ∈ [fminr, fmaxr]], the multi-step prediction filter for linear in-phase axial forward and backward prediction in the F-X domain can be determined by the following equation:

Oil and gas accumulation theory and exploration and development technology (5)

In which: * stands for compound yoke; L represents the length of the prediction filter; Pj(f) stands for prediction filter. These equations correspond to a special type of autoregressive model. The forward autoregressive equation (19) and backward autoregressive equation (20) are realized by jumping forward and backward by α steps at a time. The prediction filter Pj(f) of step α can be calculated by autoregressive equations (19) and (20). Parameters α = 1, 2, ..., αmax are step factors, which are used to extract the prediction filter of frequency αf from frequency F. Since the step factor is a positive integer, it is obvious that the low frequency part provides important information for the data reconstruction algorithm. The upper limit of the step size αmax depends on the number n of seismic traces and the length l of the prediction filter, and this parameter is given by the following formula.

Oil and gas accumulation theory and exploration and development technology (5)

Here it is. ] represents the integer part.

When the prediction filter of high-frequency data x(f') is calculated from reconstructed low-frequency data x (f) by multi-step autoregressive method, similar to Spitz interpolation method, missing data can be reconstructed by known data and prediction filter. Forward and reverse autoregressive reconstruction equations are as follows

Oil and gas accumulation theory and exploration and development technology (5)

Assuming that seismic data contains L linear in-phase axes with different slopes, the effective frequency band range of seismic data is [[fmin, fmax]]. The realization steps of spatial aliasing irregular track missing seismic data reconstruction are as follows: (1) First, the original seismic data are transformed into the F-X domain, and the low-frequency seismic data without spatial aliasing are reconstructed by using the minimum norm Fourier method [[fminr=fminr]] to obtain the low-frequency earthquake. For limited bandwidth signals without spatial aliasing, the seismic data reconstructed by FRMN has high accuracy. (2) using equations (19) and (20), the prediction filter PJ(f') is used to extract high frequency components [[fminr, fmaxr]]；] from the low frequency band; (3) using the known data and the prediction filter Pj(f') to reconstruct the missing seismic data; (4) Finally, the reconstructed seismic data are converted back to T-X domain. When encountering complex seismic data, the in-phase axis may not meet the linear assumption, so the seismic data can be divided into several small time-space windows for reconstruction. To sum up, the prediction filter of missing data with high frequency component f ′ = α f ′ = α f is extracted from the low frequency [[fminr, fmaxr]] data without spatial aliasing, and then the high frequency component of missing data is calculated by using the known data and the prediction filter, and finally the multi-step autoregressive reconstruction is completed.

3 theoretical data examples

In order to verify the effectiveness of the multi-step autoregressive algorithm, in this section, we apply the algorithm to theoretical data, reconstruct the missing trajectory and encrypt the interpolation. The first theoretical data, as shown in figure 1(a), consists of seven linear in-phase axes with different slopes, and its F-K spectrum contains serious spatial aliasing (as shown in figure 1(c)). * * * There are 8 1 channels, the channel spacing is 5m, the sampling interval is 2ms, and the number of sampling points is 90 1. Figure 1(b) is the data obtained by randomly selecting 40% seismic traces from the original data. The number 1(d) is the F-K spectrum corresponding to the number 1(b). As can be seen from Figure 1(d), the lack of seismic traces leads to serious noise on the F-K spectrum.

Figure 1 theoretical example of multi-step autoregressive method

Fig. 2 Theoretical joint application of minimum norm Fourier reconstruction method and multi-step autoregressive method (1)

Fig. 2(a) shows low-frequency data reconstructed by FRMN method, and its F-K spectrum is shown in fig. 2(c). MSAR algorithm uses the reconstructed low-frequency data to extract the prediction filter to reconstruct the high-frequency part of the data. A prediction filter that estimates the low-frequency end of data by extrapolation of the prediction filter. The complete data reconstructed by FRMn+MSAR method is shown in Figure 2(b), and its corresponding F-K spectrum is shown in Figure 2(d). Compared with the F-K spectrum of the original data (figure 1(c)), the noise caused by sampling omission is eliminated. Compared with the original data (figure 1(a)), the missing seismic trace is filled, and the continuity of the linear event axis is better.

Fig. 3 Theoretical joint application of minimum norm Fourier reconstruction method and multi-step autoregressive method (2)

Fig. 4 corresponds to the f-k spectrum of the data in fig. 3.

Fig. 5 Practical application of minimum norm Fourier reconstruction method and multi-step autoregressive method.

In order to further verify the applicability of the algorithm in complex situations, we selected single shot data from Marmousi model data (Figure 3(a)). * * * There are 96 tracks, the track spacing is 25m, the sampling interval is 4ms, and there are 750 sampling points. Randomly remove the data of 27 channels (Figure 3(b)), and reconstruct the data by FRMN+MSAR method. Fig. 3(c) shows the low frequency data reconstructed by FRMn method, and fig. 3(d) shows the complete single shot data reconstructed by FRMn+MSAR method. Because the model is very complex, there is spatial aliasing in the F-K spectrum of the original single shot data (Figure 4(a)). Fig. 4(b) is the F-K spectrum corresponding to fig. 3(b), and it can be seen that it contains serious noise. Figs. 4(c) and 4(d) are F-K spectra corresponding to 3(c) and 3(d) respectively. After reconstruction, the noise in the F-K spectrum of the data is eliminated, the missing trace is filled, and the in-phase axis also maintains good continuity.

Fig. 6 corresponds to the f-k spectrum of the data in fig. 5.

4 Examples of actual data

In this section, we will reconstruct the actual data to verify the applicability of FRMN +MSAR method. Select some data of * * * offset seismic profile (Figure 5(a)), with a total of 20 1 trace, a trace distance of12.5m, and a sampling interval of 2ms ... randomly select 30% seismic traces for reconstruction (Figure 5(b)), and Figure 5(c) shows the low-frequency data reconstructed by FRMN method. Figs. 6(a), 6(b), 6(c) and 6(d) are F-K spectra corresponding to figs. 5(a), 5(d), 5(c) and 5(d) respectively. It can be seen that the F-K spectrum of the data changes little before and after reconstruction. After reconstruction, the missing trajectory of the data is restored, and the in-phase axis is continuous, and the reconstruction result is close to the original data.

5 conclusion

Based on the minimum norm Fourier reconstruction method and the multi-step autoregressive method, this paper reconstructs the seismic data with spatial aliasing. Multi-step autoregressive method is an extension of Spitz method and is also based on the assumption of approximate linear in-phase axis. Therefore, it is usually difficult to satisfy this assumption when dealing with complex seismic data. At this time, the method of small time-space window can be used to calculate, and it can be considered that the assumption of approximate linearity is satisfied in the small time-space window. However, if the time-space window is too small, the amount of data will be insufficient, resulting in poor reconstruction results or may not be reconstructed. As we all know, in order to solve most geophysical problems, it must be based on some assumptions. Usually, these assumptions will be partially violated when dealing with actual data. In fact, this method can also achieve ideal reconstruction results for moderately curved in-phase axes, which shows that the reconstruction method in this paper has good stability. In fact, this method has also achieved good reconstruction results for seismic data with large spacing. Theoretical data and practical data reconstruction experiments verify the effectiveness and practicability of the reconstruction method in this paper. In addition, the reconstruction effect of seismic data is also related to the complexity of original data, the nature of frequency spectrum, the number and location of missing seismic traces and the distance between missing traces. It is necessary to further study the influence of these factors on the reconstruction algorithm.

refer to

[1]Eiken O, Haugen G U, Schonewile M A, and Duijndam A J W. Mature method for collecting seismic data of towed cables with high repeatability [J]. Geophysics, 2003,68 (4):1303 ~1309.

[2]Wever A, Spetzler J. Criteria for locating sources and receivers in time-lapse seismic acquisition: 74th Ann. Internal. MTG., SEG, Extended Abstract, 2004:23 19~2322.

Duijndam AJ W, Schonewille M A, and Hindricks C O H. Reconstruction of band-limited signals with irregular sampling along a spatial direction [J]. Geophysics,1999,64 (2): 524 ~ 538.

Hindrix K. Duidam A.J.W. Reconstruction of 3D seismic signals irregularly sampled along two spatial coordinates [J]. Geophysics, 2000,65 (1): 253 ~ 263.

Liu B, Sacchi M. Minimum weighted norm interpolation of seismic records [J]. Geophysics, 2004,69 (6):1560 ~1568.

Zhang Zhiyong, Zhang Zhiyong. Non-uniform sampling and aliasing Fourier reconstruction of seismic data [J]. Journal of Seismogeology, 2002. Geophysics, 2007,72 (1): 21~ 32.

[7]Feichtinger H G, Grochenig K, Strohmer T. Effective numerical methods in non-uniform sampling theory [J].Numerische Mathematik,1995,69: 423 ~ 440.

Zwartjes P.M., Fourier Reconstruction of Sparse Inversion: Doctoral Thesis of Delft University of Technology, 2005.

Takalo, Hetty, Ihalainen. A course of univariate autoregressive spectrum analysis [J]. Journal of Clinical Monitoring and Computing, 2005,19 (6): 401~ 410.

[10] canals L. Random noise reduction: 54th ampere. Internal installation. , SEG, extended abstract,1984: session: s10.1.

[1 1]Spitz S. f-x seismic trace interpolation [J]. Geophysics, 199 1, 56(6):785 ~794.

[12]Porsani M J seismic trace interpolation using half-step prediction filter [J]. Geophysics,1999,64 (5):1461~1467.

[13] Marple S.L. Digital spectral analysis and its application. Cleaves, Englewood, New Jersey: prentiss Hall Company, 1987.

Multi-step autoregressive reconstruction of seismic records [J]. Geophysics, 2007,72 (6):11-18.