Abstract
In the field of complex underground geological structures and irregular topography, prestack seismic data often have a low signal-to-noise ratio (SNR), in which weakly reflected signals are buried beneath strong incoherent, and scattered noise. Stacking, such as beamforming along the moveout surfaces of coherent local events, can significantly improve seismic data quality. Accurate and efficient estimation of the moveout for an irregular acquisition geometry and uneven illumination is important in a complex environment. In this paper, a new optimal stacking approach for enhancing weak prestack reflection signals is presented. The proposed method mainly includes regional division and moveout estimation. Optimal stacking should be implemented within local time and space domains. Based on beam-ray theory, we designed a reasonable regional division of the common-shot (CS), common-receiver (CR) and common-middle-point (CMP) domains. Then, we proposed using the sparse radon transform and dynamic waveform matching method to estimate the moveout surfaces of local reflection events. The sparse radon transform was applied to obtain the linear moveout to ensure the correctness of the reflection wave direction. The residual nonlinear disturbance was estimated using the dynamic waveform matching method. Tests on synthetic and field data demonstrated the effectiveness of the proposed method, which can effectively improve the SNR of prestack seismic data and attenuate incoherent noise.
1. Introduction
Seismic data can be produced by the excitation of the surface of a source wavelet, which propagates in underground media and is then received by a geophone. The data are extremely complex due to the underground geological structure and irregular topography, and contained various strong noises, particularly in the piedmont zone data. The fusion of seismic signals and interference noise makes the signal-to-noise ratio (SNR) of prestack seismic data extremely low, and irregular source and receiver acquisition geometries and uneven common-middle-point (CMP) fold times also decrease the data quality. Specialized heavy preprocessing is therefore required to suppress noise and improve the SNR of the seismic data.
Simple source and receiver arrays cannot obtain reflected data for reliable processing and imaging. The conventional method uses smart multidimensional local stacking of neighboring traces to enhance the seismic data. The premise of the stacking is that the wavelet of the reflection signal is flattened and has the same travel time. However, severe time differences between seismic traces cause wavelet distortion and signal-resolution degradation by a simple supergroup. The supergrouping method after normal-moveout (NMO) correction can effectively enhance the seismic data (Bakulin et al.2016, 2018a, 2018b). For example, based on NMO-corrected CMP gathers, we can use the tau-p transform (Tatham et al.1983), F-X deconvolution (Canales 1984), FXY filter (Chase 1992), independent component analysis and local singular value decomposition (Bekara & Van der Baan 2006) to suppress noise and enhance prestack seismic data. However, NMO-corrected supergrouping relies on the global hyperbolicity of reflection events and is model-driven, and these methods are not suitable for complex near-surface or overburdened conditions.
Actual moveouts can be estimated by the data-driven multidimensional local stacking method, which does not require simplified assumptions regarding the global hyperbolicity of events, and this can further enhance seismic data. Baykulov & Gajewski (2009) proposed partial stacking along the common-reflection-surface (CRS) travel time in the common-offset (CO) domain to construct a new supergather. However, the use of global zero-offset CRS operators is only applicable to the hyperbolic travel-time behavior of seismic events, and only a partial summation can be made with the non-hyperbolic travel-time behavior. Li et al. (2011) used the multiparameter CO CRS travel-time formula to obtain a finite offset section to calculate partially stacked CRS supergathers. Xie (2017) proposed using the differential evolution algorithm to determine eight wavefront attributes of the travel-time operator of the 3D CO CRS method. Buzlukov et al. (2010), and Buzlukov & Landa (2013) proposed using the local parabolic second-order travel-time approximation of reflection events in the CO domain to search for locally coherent events and partial summation along the estimated trajectories, to enhance the prestack data. Neklyudov et al. (2015) proposed the use of smart supergrouping to enhance low-frequency data and suggested a CO summation and a weighted stack with complex coefficients in the poorly sampled direction to preserve higher frequencies. Khaidukov et al. (2016) applied this method with operator-oriented interpolation, to coherently stack 3D marine data using floating ice acquisition.
In addition, the multi-focusing (MF) method (Berkovitch et al.1994, 2008) can be applied to improve the prestack SNR by partially stacking coherent seismic events that do not necessarily belong to the same CMP gather. This method requires different moveout corrections. For example, the zero-offset MF (ZOMF) formula is based on a quasi-hyperbolic approximation of the travel-time surfaces. However, this operator is only valid for short offsets, and for long offsets or strong lateral velocity variations, the quasi-hyperbolic zero-offset operator is no longer applicable. Rauch-Davies et al. (2013) proposed the use of a local CO MF approximation for a travel-time stacking surface, to accurately approximate the non-hyperbolic travel times of seismic events. Multidimensional data-driven local stacking can be considered a delay-and-sum beamforming method. In complex geology, the time delay is a nonlinear function of the distance in the local region. Bakulin et al. (2018c) proposed a nonlinear beamforming (NLBF) approach in the CMP offset domain to attenuate noise and enhance prestack land seismic data, and Bakulin et al. (2019, 2020) also applied the NLBF method in the cross-spread domain to enhance the local coherent nonlinear events for the 3D prestack land seismic data. Sun et al. (2022) proposed an efficiency-improved genetic algorithm to estimate the operator coefficients of the NLBF and improve computational efficiency.
By combining model-driven robustness and data-driven accuracy, a new optimal stacking method was proposed in this study for enhancing prestack seismic data. The stacking method is implemented in a local region, and the time delay is divided into linear and nonlinear disturbances. We propose using the sparse radon transform method to estimate the linear moveout and the dynamic waveform matching method to correct the residual nonlinear disturbance.
2. Regional division
A prerequisite for enhancing seismic data is that it should be relevant and similar in the time and space domains. Based on the acquisition-propagation-receiver process of seismic data in subsurface media, reflection events can be regarded as being stacked constructively from a half-Fresnel band, where the dominant energy along the ray path of the reflection wave coming from the preferred direction is bounded into a beam-ray tube. Therefore, this process must be implemented within the local time and space domains. In a local region, the similarity of seismic data can be used to attenuate incoherent and random noise and enhance the SNR of the data.
Regional divisions require a deliberate design. If the region is too small, there are too few valid signals to enhance data and eliminate noise. When the region is excessively large, the effective signals originate from the nonlocal underground reflection area and have no strong similarity; therefore, coherent stacking fails.
The arrangement of prestack seismic data uses different methods, and the implementation of regional divisions also varies. For the common-shot (CS) and common-receiver (CR) gathers, due to the interchangeability of shot and receiver positions, the regional division in the CS domain is the same as in the CR domain. An example of the CS domain was shown in figure 1. The original shot and receiver positions are represented by the black and green circles, respectively. The blue boxes represent the selected local regions. As the offset increases, the apparent velocity of the reflected wavelet propagating to the surface increases, such that the wavelet becomes longer at a large offset. Therefore, we set a small region for a short offset and a large region for a long offset. The size of the local region is shown in equation (1).
where |${{{\bf x}}}_{max}$| is the maximum offset, |${{\bf x}}$| is the distance of source and receiver, and a and b are the minimum and maximum regions, respectively.
The regional division in the CS domain.
Figure 2 shows the results of the regional division in the CMP domain. The first and fourth rows show the source and receiver positions, respectively. The CMP gather was rearranged as shown in the black box in figure 2, and the regional size increases with the offset.
The regional division in the CMP domain.
3. The constructive superposition
When the regional division is determined, the difference in travel time in the local region must be estimated to realize constructive stacking. The time difference |$\xi ( {{{\bf x}},t} )$| can be divided into two parts, the linear phase change |$\Delta ( {{{{\bf x}}}_i,t} )$| and nonlinear disturbance |$\delta ( {{{{\bf x}}}_i,t} )$|, in which the linear part is estimated using the sparse radon transform and inverse transform and the nonlinear disturbance is estimated using the dynamic waveform matching method. The linear part is important as the effectively reflected wave propagation has a preferred directionality, and the travel-time difference can be seen as a linear change for the local horizontal surfaces. Therefore, the total time difference becomes equation (2):
where Ω is the local region.
3.1. Sparse radon transform estimate of the linear difference
For prestack seismic data, local linear reflection events can be represented by stacking beam rays from different directions. This method can be implemented using tau-p or radon transform (RT) (Trad et al.2003). The RT and inverse RT are expressed as shown in equations (3) and (4):
where p is the ray parameter, and |${{\bf x}}$| is the spatial position. By applying the Fourier transform in the τ direction, we obtained the frequency domain RT, as shown in equation (5).
where Ω is the circular frequency. Equation (5) can be expressed as a matrix form, as shown in equation (6):
where |${{\bf m}},\,{{\bf d}}$| is a complex vector and the matrix |${{{\bf G}}}^T$| is a RT operator. This leads to equation (7):
Similarly, the matrix expression of the inverse RT can be obtained, as shown in equations (8) and (9):
where the matrix |${{\bf G}}$| is an inverse RT operator.
As |${{{\bf G}}}^T$| is usually rank deficient, directly solving the equation (6) is ill-conditioned and unstable. Herein, we introduced additional regularized constraints, as shown in equation (10):
where |$q = 1$|, λ is the hyperparameter. A linear radon spectrum can be obtained by solving the inverse problem described before. For the local reflection event, the maximum value was extracted as the beam-ray parameter, as shown in equation (11):
Based on the sparse radon spectrum, the local reflection event can be estimated using the inverse RT, as shown in equation (12):
Therefore, the travel-time difference |$\Delta ( {{{\bf x}},t} )$| between the local seismic trace |$d( {{\bf x}} )$| and the seismic trace |$d( {{{{\bf x}}}_i} )$|, where |${{{\bf x}}}_i$| is the central position |${{\bf x}}$|, can be obtained using equation (13):
3.2. Dynamic waveform estimate of the nonlinear disturbance
The dynamic waveform matching (DWM) method can calculate the signal similarity by extending and shortening the signal sequences, and searching for the shortest path between the reference and test signals (Anderson & Gaby 1983). The DWM can be used to estimate the nonlinear disturbance between the seismic trace |$d( {{{\bf x}},t} )$| and reference trace |${d}_1( {{{\bf x}},t} )$|, which is the result obtained by the sparse RT. For a single trace of prestack data, the difference between the reference trace |${d}_1( t )$| and test trace |$d( t )$| at each sampling time is expressed as shown in equation (14):
where |$e( {i,j( i )} ) = e[ {i,j} ]$| is the time-warping function or time-warping matrix, where the size of the matrix is |$nt*nt$|, and |$nt$| is the number of time samples of the seismic trace. Figure 3 shows the two signals and the corresponding time-warping matrix.
Diagrammatic sketch for the two test signals of the DWM method.
The DWM has an optimization problem. Its purpose is to find a warping path |$W = {w}_1,....,{w}_{nt}$| from the time-warping matrix, as shown by the black dotted line in figure 3. When the reference and test signals matched perfectly, the cumulative differences were minimal. However, there are several paths to choose from the time-warping matrix. To determine the optimal path from a large array, the following constraints are considered:
Boundary condition: the initial and termination positions of the path are |$[1,1]$| and |$[nt,nt]$|, respectively.
Consistency condition: the two signals contain the same effective reflection events.
Monotonicity condition: the sequence of effective reflection events was the same for both signals.
Under the constraints of continuity and monotonicity, the path of each point has only three directions. For example, if the path has already passed through the one-point |$( {i-1,j-1} )$|, as shown in figure 4, the next passing point may only be one of the following three cases: |$( {i,j-1} )$|, |$( {i-1,j} )$| or |$( {i,j} )$|.
Schematic diagram of the searching path direction.
Starting from the initial position, |$[1,1]$|, the cumulative distances of all paths was calculated. The cumulative error is calculated by selecting the minimum value from the possible paths, as shown in equation (15).
where |$\gamma ( {1,1} ) = e( {1,1} )$|. The desired optimal path minimizes the errors that accumulate along the path. When in the terminated position |$[ {nt,nt} ]$|, the corresponding cumulative error is |$\gamma ( {nt,nt} )$|. Therefore, we can calculate the minimum distance at time |$i = nt$| using equation (16):
where r denotes the final cumulative error. Based on the shortest cumulative error r, we searched for a path in the reverse direction to obtain a regular path |$w( i )$|, as shown in equation (17).
where |$w( {nt} ) = \mathop {arg\min }\limits_j \,\,\,\,\gamma ( {nt,j} )$|.
The travel-time differences between the seismic trace |$d( {{{\bf x}},t} )$| and reference trace |${d}_1( {{{\bf x}},t} )$| can be calculated using the regular path |$w( {{{\bf x}},t} )$|, as demonstrated by equation (18):
where |$\delta ( {{{\bf x}},t} )$| denotes time difference. The corrected prestack data |${d}_c( {{{\bf x}},t} )$| can be obtained as follows in equation (19):
3.3. The optimal stack for enhancing the prestack seismic data
Based on linear and nonlinear disturbances, the total time differences were estimated using equation (2). The optimal trace was obtained by stacking all traces in the local region, as shown in equation (20).
where |$d( {{{\bf x}},t} )$| represents the original seismic data, |${\alpha }_i$| represents the weighting coefficient, |$\xi ( {{{{\bf x}}}_i,t} )$| represents the estimated time difference for the local trace and |$\tilde{d}( {{{\bf x}},t} )$| represents the stacked seismic data. The weighting coefficient |${\alpha }_i$| is estimated by the local correlation (Fomel 2007; Wu et al.2019).
4. Example
In this section, the synthetic model and field data were used to verify its validity for enhancing seismic data. Figure 5a shows the synthetic local seismic data with different dominant frequencies of the Ricker wavelet and moveouts. The random disturbances and white Gaussian noise (SNR = 2) shown in figure 5b are used to simulate the complex geological structure, irregular topography and acquisition noise. The black line shown in figure 5c represents ideal seismic data, where the result is the central trace of figure 5a. The supergrouping and stacking of the local seismic data are indicated by the blue line in figure 5c. The severe moveout between the seismic traces and the random noise suppressed the reflection signal, and direct stacking cannot distinguish between the effective seismic signals and noise, as the wavelet of the reflected signal has a different travel time and is non-flattened. Stacking seismic data along the time difference can effectively enhance the seismic data quality. The red line shown in figure 5c is the stacking result, which shows that the proposed method can effectively enhance the reflection signals and suppress noise.
The synthetic data test: (a) the ideal data; (b) the noisy and random disturbance data; (c) the stack result where the black line is the original wavelet. The blue line is the direct result, and the red line is the result from the proposed method.
Then, we tested the proposed method using field data. Figure 6a shows the original seismic data. As can be seen, in addition to effective reflection waves, a large amount of linear noise and random noise is included in the seismic data. Through estimating the nonlinear moveout surfaces by the sparse RT and DWM method, we can stack the local seismic data along moveout surfaces to enhance the reflection events, the stacking results are shown in figure 6b. With the optimal stacking method, the reflection coherent events are preserved, and undesired incoherent features and noise are suppressed. The zoomed-in original data and optimal stacking result in the blue and red boxes in figure 6 are shown in figure 7. For local seismic data, the moveout of the reflection wavelet is nonlinear and the shape of wavelets between traces changes drastically in figure 7a and c. After stacking with the local traces, the results are shown in figure 7b and d, respectively. As can be seen, the reflection events are enhanced and the noise is suppressed.
The field data test in the CS domain (shot number 2412): (a) the original seismic data and (b) the optimal stacking result using the proposed method.
Furthermore, the test results of the two CS seismic data are shown in figures 8 and 10. The zoomed-in data are shown in figures 9 and 11, respectively. With the proposed nonlinear moveout optimal stacking, the wavelet features with lateral variations are effectively eliminated. The quality of seismic data can be effectively improved. Figure 12 parts a and b show the imaging result by reverse time migration method using the original seismic data and the enhanced seismic data, respectively. The migration artifact and unwanted events is effective attenuated and the continuity of reflection events have a significant improvement in figure 12b.
The field data test in the CS domain (shot number 2213): (a) the original seismic data and (b) the optimal stacking result using the proposed method.
The field data test in the CS domain (shot number 2333): (a) the original seismic data and (b) the optimal stacking result using the proposed method.
The migration result of field data: (a) using the original seismic data and (b) using the seismic data by the proposed method.
5. Conclusions
In prestack seismic data, an optimal stacking method was developed along the nonlinear moveout surface of a local reflection event. The nonlinear moveout was estimated in the local region, where the regional division ensured that all seismic traces had similar reflection signals. The time delay of the local seismic data was divided into a dominant linear difference and a tiny nonlinear disturbance to improve the accuracy of the nonlinear moveout. The linear parameter was estimated using sparse RT to determine the preferred direction of propagation. DWM was proposed to estimate the nonlinear disturbances. Tests on synthetic and field data demonstrated that the proposed method can enhance reflection signals, attenuate incoherent noise and suppress near-surface scattered noise to effectively improve the SNR of prestack seismic data.
Acknowledgments
The authors thank the sponsors of WPI group for their financial support and helps. WPI's research work is financially supported by the National Natural Science Foundation of China (grant no. 42174135 and 42074143), the National Key R&D Program of China (grant no. 2018YFA0702503) and the SINOPEC Key Laboratory of Geophysics (grant no. 33550006–20-ZC0699-0011 and 33550006–22-FW0399-0019). The authors also thank the support of the Fund of the National Key Research and Development Program of China (grant no. 2019YFC0604902).
Conflict of interest statement. None declared.
