Estimation of seismic local slope using Hilbert transform noise correction method

In seismology, the local slope is essential for various seismic forward and inversion methods, including seismic migration imaging, seismic tomography, structure prediction, and other core seismic data processing methods. However, high-frequency random noise is unavoidable in real seismic data. Therefore, a robust estimation of event-slope attributes becomes important under high-frequency random noise. This research explores a robust and efficient method for estimating local slopes. The relationship between the frequency response of the Hilbert transform and Fourier transform is analyzed. The derivative operator constructed by the Hilbert transform can effectively weaken the effect on the energy enhancement of the high-frequency random noise. Minimizing the quadratic function of the proportionality factors was used to obtain the proportionality factors of the noise correction. Finally, the Hilbert transform and noise correction derived a new linear plane-wave destruction filter operator. The linear operator can effectively estimate the seismic events’ local slope with high-frequency random noise. The calculation results of the numerical examples show that the linear plane-wave destruction filter estimation local slope method proposed in this study has high calculation accuracy and efficiency.


Introduction
It is a basic process to estimate the local slope attribute of the seismic events in the seismic data processing, which is widely used.These processing methods are: stereo tomography (Billette and Lambaré, 1998;Lambaré, 2008;Wang et al., 2016), migration-velocity analysis (Lambaré, 2008;Li et al., 2014), wave-field separation (Liu and Li, 2020;Sun et al., 2023), data-denoising and steering filtering (Liu et al., 2011;Liu et al., 2014), Seislet-transform (Fomel and Liu, 2010), automatic normal moveout (ANMO) and imaging (Cooke et al., 2009;Fomel, 2007b), common-reflection-surface (CRS) stack (Hertweck et al., 2007) and Gaussian beam modeling (Sun et al., 2021;Yue et al., 2021).Conventional normal moveout (NMO) requires a lot of manual work to pick up the stack velocity to level the event, while ANMO does not require the manual work to pick up, which can save a lot of resources.
There are many ways to estimate the local slopes (Liu et al., 2018).Ottolini (1983) used the local-slant stacking method to extract the local slopes based on the correlation or the similarity function by using the maximum energy value after stacking.The efficiency of this method is not high, however, and the resolution of the estimated value is low.Claerbout (1992) proposed the concept of the linear plane-wave destruction (PWD) filter and gradually improved upon the method.Based on the idea of a multidimensional prediction error filter and the local plane wave model of the seismic data, the local plane wave was estimated by solving the differential equation of a plane wave with a finite difference template.It could well deal with the aliasing problem and the plane wave whose slope varies with space and time.However, the use of this method was not popular at that time.One of the reasons was that many parameters needed to be selected artificially, such as the filter type, the filter coefficient, and the size, number, and the shape of the local patches estimated by filtering, iteration times and regularization degree, which are not easy to control.Clapp et al. (1998) proposed a control filter by applying the method of the plane wave destruction filter in the preconditioned tomographic inversion problem.Crawley et al. (1999) reduced the parameters that needed to be set in the Claerbout multi-dimensional prediction-error filter, which was widely used in the offshore data regularization, and multiple suppression.Fomel (2002) proposed a nonlinear PWD, using the maxflat filter to approximate the phase shift operator and derived a polynomial expression for the local slope.The PWD filter was further improved by reducing the required parameters to one, and estimating only a quantity with a clear physical meaning, i.e., the slope parameter of the local plane-wave field.The slope was estimated as a smooth continuous function of the coordinates of the data, which did not require a time window to restrict it.However, this method ignored the influence of outliers on the interpolation results, and there was a large error in the estimated slope when there was high-frequency noise in the seismic data.Zhang (2009) discussed a design problem of filter design infinite-impulse-response (IIR) filters with the maxflat frequency response in frequency domain.Schleicher et al. (2009) compared several different methods for finding the event slopes and derived an expression of noise correction for the local slope.Liu et al. (2015) introduced the Hilbert transform to solve the plane wave differential equation, which enhanced the anti-noise ability of the slope field, and also had higher computational efficiency.proposed a deep learning network for the estimation of the slope field, which was trained using only synthetic seismic data.It has good performance in some real-seismic data, but not suitable for all of them.
Due to the restriction of the environment and human factors, the data collected in seismic exploration are often affected by random noise.Random noise is characterized by irregularity, frequency, apparent velocity, and random distribution in time and space.There are many reasons for this noise (Sheriff, 1975), such as receiving and excitation conditions, problems of the instrument itself, and environmental effects.In fact, it is that how to accurately estimate the local slope of events in noisy data, and estimate a stable and smooth slope field becomes very important.
In this paper, the basic principle of linear PWD operator using the Hilbert transform and noise correction is proposed, and Fomel's PWD algorithm is reviewed.Then, the proposed method here, and other methods for estimating the slope field, is applied to ANMO, and the effects are compared and analyzed under the conditions of non-noise and noisy seismic data.A new linear plane-wave destruction filter operator with noise correction is constructed to estimate the local slope field of a plane wave through the derivative operator using the Hilbert transform.This method effectively overcame the influence of the random noise on the slope field estimation.The results of the numerical examples show that this method can efficiently estimate the stable and smooth slope field in noisy data.

Linear plane-wave destruction filter
A wave with a plane front is called a plane wave.Plane waves do not exist in the real world.Spherical waves can be decomposed into plane waves, and plane waves can be synthesized into spherical waves.When a spherical wave travels a certain distance, the wavefront of the spherical wave can be regarded as a plane wave when viewed locally.In this case, the propagation of the wave can be viewed as a plane wave.Plane waves have certain directionality.If the wave field T -X in the domain is U(x;t), the plane wave propagating along the horizontal direction can be characterized by the one-dimensional wave equation: The Equation 1can be considered to have the following form: ( ) Where p = sin θ v is the ray parameter or event slope; θ is the angle between the propagation direction of the wavefront and the perpendicular direction; represent the derivatives of the wavefield with respect to x direction and t direction, respectively.Under the condition of θ=0, it means that the plane wave is perpendicular to the x axis and propagates along the positive direction of the x axis; p specifies the propagation direction of the plane waves.Equation ( 2) is also called the local plane wave differential equation.
In this way, the value of p can be determined from the space and time derivatives of the wave-field through equation ( 2).In the synthetic data without noise, equation (2) can quickly and accurately get the results of p.The frequency response of the derivative operator for discrete data is: (3) However, there is high-frequency random noise in the seismic field.From equation (3), in calculating the derivative, the energy of the high-frequency component of the noise will be significantly strengthened, which, in turn, means that the accuracy of estimating p gets reduced.To weaken this effect, Liu et al. (2015) used Hilbert transform instead of Fourier transform to obtain the derivative of the seismic data by using the relationship between the Fourier transform and the Hilbert transform on the frequency response of derivative operator (Pei and Wang, 2001).The ideal frequency response of the Hilbert transform derivative operator is: The derivation of the discrete data by the Hilbert transform can effectively reduce the effect of energy enhancement of high-frequency random noise.
To further improve the accuracy of p estimation in the presence of noise, and correct the noise, the linear PWD method proposed by Claerbout (1992) is used here.For the stability of the calculation results, establish a small space-time window, W, around each point, (x 0 ,t 0 ) to be calculated in the seismic profile, is established in which the space-time window, (i,j) ∈ W, contains points (x i ,t j ).Note that ( ) represent the derivatives of the wave field in direction x and direction t respectively, and the derivative value obtained by the Hilbert transform.The ray parameter is estimated by using the equation (2), , and equation ( 2) is converted into a quadratic function as follows: (5) The quadratic function is minimized with respect to p, and let ∂Q| ∂p=0, and A more stable ray parameter value can be estimated by using equation ( 6).We used the sliding window normalized correlation (Claerbout, 1992) to determine whether the estimated p value matched the data well: (7) From equation ( 7) that 0 ≤ C ≤ 1, and when C is close to 0, it indicates a low correlation, and when C approaches 1, it indicates high correlation.Except for when p = 0 is C = 0, equation ( 7) is very close to the classical similarity function method.Kington (2015) summarized and analyzed the relationship betweed p and C. The expression for the similarity function (Stoffa et al., 1981) is as follows: Where N represents the number of events calculated in the time-space domain, and w is the length of the time window.Schleicher et al. (2009) proposed a correction method for equation (2), and pointed out that when p ≠ 0, there is another corresponding parameter, q = 1/p, namely: Correspondingly, the expression of the minimized quadratic function of q can be estimated from the following: (10) By comparing equations ( 6) and ( 10), the values of p and q are equivalent, and the amount of estimation is the same.It follows, then, from equation (7) that: ( ) The normalized correlation function can be estimated directly by calculating p and q.When there is no noise in the data, C = 1.However, there is always random noise in the real data, and so the estimated values of p and q are inaccurate.They can be expressed as the actual true value multiplied by a scale factor, ε, less than or equal to 1.It is further assumed here that the influence of the noise on the estimation of p and q is the same.Then the proportionality factors are considered to be approximately the same, and this assumption is reasonable, and which can be expressed as follows: (12) According to equation ( 12), the linear plane-wave destruction filter operator using the Hilbert transform, noise correction is obtained, and the true value expression of the slope field is: The implementation flow of the algorithm is shown in Figure 1.

Nonlinear plane-wave destruction filter
The PWD was considered to be one of the methods to estimate the value of p with higher accuracy, and which utilizes a Taylor-series expansion to fit the low-frequency response of a time-domain all-pass filter.The estimation of the local slope was used to predict the data of each event by PWD through adjacent seismic events.According to the continuity of the slope of adjacent seismic events, the objective function of the difference between the predicted data of each event and the original seismic event data is constructed.The objective function was used to reach the optimal iterative solution in the sense of the least squares, and was also used to estimate the local slope field of the events.The estimation of the local slope field of the events can be regarded as a nonlinear inversion problem, which can be solved by the Gauss-Newton iteration method.The disadvantage of estimation method using PWD is that it needs iteration, and so the estimation speed is slow.
According to the local plane-wave differential equation ( 2), if the slope, p, does not depend on the change of time, t, it can be transformed into the frequency domain and the differential equation becomes: The general solution of the differential equation in ( 14) is: Where Û is the Fourier transform of U. The exponential term in the above equation represents the time shift of the seismic trace with slope p and trace spacing x.In other words, a plane wave can be perfectly predicted by a two-stage frequency-space-domain prediction-error filter.And the amplitude spectrum of the complex exponential term in the frequency domain is 1.This maintains the important property of a plane wave, which is, the invariance of wave propagation energy.Fomel (2002) proposed a Z-transform version of the three-point finitedifference PWD filter that approximates the time-shift operator as follows: where It can derive filter forms of more orders.The phase accurately fit the frequency response of the phase-shift operator at low frequencies.A twodimensional PWD filter using the expression ( 16) was constructed as follows: If D(p) was used to represent the convolution operator between the two-dimensional PWD filter, D(Z t , Z x ), and seismic data, d, then p is the local slope field.Subsequently, the objective function, D(p)d, was constructed.The estimation of p was transformed into finding the minimum of the objective function: where || • || 2 2 denotes the 2-norm, Ɛ denotes the regularization perturbation factor, and S is the shaping regularization operator (Fomel, 2007a), which controls the smoothness of the estimated local slope field.Due of the continuity of the slope field, the stability and smoothness of its estimated values were also expected to achieve and conform to the actual situation.This is also illustrated by the following numerical examples.

Experiments and Analysis
To verify the method proposed in this paper, the slope field estimated by the different P calculation methods was used for ANMO of common-midpoint (CMP) gathers under the conditions of the non-noise and noise data.The different methods for estimating the P value were: method 1, which was the constraint ratio of Fourier transform in equation (2); method 2, which was the constraint ratio of Hilbert transform of equation (2) given by Liu et al. (2015); method 3 which was equation ( 19) using PWD given by Fomel (2002); method 4, with the Equation 13 estimation method using the Hilbert transform and noise correction, that is, the method proposed in this paper.
The classical hyperbolic reflection time-distance curve equation is: where t is the travel time relative to offset x, t 0 is the zero offset travel time, and v nmo (t 0 ) is the stack or RMS velocity.The expression of p is: Therefore, v nmo (t 0 ) can be replaced by mathematical transformation from equations ( 20) and ( 21) to obtain the NMO correction equation using the slope: The above equation gives an automatic NMO technique using the slope information of the data, called the ANMO correction.The method performs the NMO correction according to the slope field information of the data, and avoids the problems in that a lot of manual work is needed to pick up the stack velocity for correction, and that the wavelet stretching occurs at the far offset in the conventional NMO correction.At the same time, this correction technique is sensitive to the smoothness and accuracy of the slope field.
A synthetic noiseless CMP gather-record was selected, as shown in Figure 2a.There were 201 traces in total, with a trace spacing of 10 m, and a sampling interval of 4 ms. Figure 2c shows the local slope field estimated by Method 4. The CMP gather-record of Figure 2a was superimposed on Figure 2c.For each calculation, location point, (x 0 ,t 0 ) here, the size of the space-time window, N t = 10, N x = 10, was taken.The smooth and the stable slope field was estimated by calculating the ratio of the division in the equation with the shaping regularization operator.Finally, the ANMO was carried out with equation ( 22), and the results are shown in Figure 2d.The events of each layer are flat, which proved the correctness of method 4, and the smooth stability of the calculation results.For the non-noise CMP gather-records, each slope estimation method can accurately estimate the value of p, and then flatten the events.To test the stability and accuracy of the different methods to estimate the value of p under the condition of noise, Gaussian white noise was added to original gather records, whose results were shown in Figure 2b.The slope field (Figure 3a, c, e, g) estimated by the different methods were superimposed with the CMP gather records in Figure 2a.The results of the ANMO flattening events (Figure 3b, d, f, h) showed the effectiveness of the different methods used here.
Because the slope field estimated by method 3 was not smooth enough, it led to serious "event jitter" and "scratch" phenomena.The triangular smoothing function was selected with the size of the space-time window, N t = 10 and N x = 10, to smooth the slope field.In this paper, the triangular smoothing function was used to repeat smoothing 5 times, and the result was shown in Figure 4a.The slope field became smoother, and the mutation values were reduced.The effect after ANMO was shown in Figure 4b.
The numerical results verified the effectiveness of the proposed method.With the same smoothing time window and regularization constraints, different methods estimated the local smooth-slope field.However, it was also that methods 1 and 3 were very sensitive to the noise, and the slope field estimated by methods 1 and 3 fluctuated greatly in the global range, especially in the region where there was no effective signal.For the event part of the effective signal, the estimated slope field was easily affected by the noise, which resulted in an anomaly and an inaccurate value.The ANMO results (Figure 3b) using method 1 showed that the event was not flattened, which indicated that the slope at the effective signal was distorted by the noise and thus the estimate was inaccurate.Fourier transform enhanced the energy of high-frequency noise.There were many "scratches" in the figure because there were many abrupt changes in the estimated slope field, and the slope gradient was large.This phenomenon was also evident in Figure 3f.The ANMO results (Figure 3f) using method 3 showed that the event was flattened, which proved that the estimation of the slope field was accurate.However, the "event jitter" phenomenon, which occurred at the upper and lower edges of the event, indicated that the estimated slope field was distorted by the noise.It was affected by the noise at the far offset.Compared with the result of method 1, the slope field estimated by method 3 was reliable at the effective signal.For all that, the slope field was still affected by the influence of the noise.Both methods 2 and 4 were based on the Hilbert transform, which would not enhance the energy of the high-frequency noise in the process of operation.Figure 3c showed that the estimated slope field results were stable and smooth.The estimated slope field was accurate at the effective signal.The influence of the noise was suppressed at the non-signal.The ANMO results (Figure 3d) using method 2 showed that the event was flattened, and only the "scratch" phenomenon exists at the edge.The accuracy of the slope estimation was not good in the presence of the noise for shallow reflecting layers with large offset.The slope field (Figure 3g) estimated by method 4 was more stable and smoother than the results (Figure 3c) estimated by method 2. The event shown in Figure 3h was flattened except for a slight misalignment at the far offset, and was minimally affected by the noise, which, in turn, showed good anti-noise ability.
Figure 3 showed that under the same smoothing condition and regularization constraint solution, the method in this paper had stronger noise immunity, and the stability and smoothness of the slope field estimated were also better.At the same time, the calculation efficiency was the same with methods 1 and 2, but faster than that of method 3. Compared with the results shown in Figure 3f, the "event jitter" phenomenon was well suppressed, but there was still no good effect at the edge of the event due to the inaccurate slope field (Figure 4b).At the same time, the "scratch" phenomenon has not been well controlled, which was caused by the large gradient of the slope field at the scratch position.

Conclusions
Several methods for estimating the local slope of events are introduced in this paper.On their bases, a method is proposed for estimating the local slope of the linear plane-wave destruction filter operator using the Hilbert transform and noise correction.It is a stable, efficient, and anti-noise slope field estimation method.The slope field estimated by the ANMO under the conditions of the non-noise and noise data was used to carry out the flattening experiment on the CMP gather records.The results of the numerical examples showed that the slope field estimated by the proposed method had high stability, smoothness, and strong anti-noise ability.

Figure 1 .
Figure 1.Flow chart of the linear plane-wave destruction filter operator using the Hilbert transform and noise correction

Figure 2 .
Figure 2. Schematic diagram of trial calculation results for a synthetic noiseless CMP gather record.(a) CMP gather record; (b) CMP gather record with white noise added; (c) slope field (CMP gather records are superimposed in the figure); (d) ANMO flattens the events.

Figure 3 .
Figure 3.The slope fields estimated by different methods with noise as shown in Figure 2a and the ANMO results.(a) Slope field estimated by method 1 (the CMP gather record is superimposed in the figure, the same below); (b) ANMO results of slope field estimated by method 1; (c) slope field estimated by method 2; (d) ANMO results of slope field estimated by method 2; (e) slope field estimated by method 3; (f) ANMO results of slope field estimated by method 3; (g) slope field estimated by method 4; (h) ANMO results of slope field estimated by method 4.

Figure 4 .
Figure 4.The slope fields estimated by method 3 using the smooth function, (a) method 3 estimates the slope field using the smooth function (CMP gather records are superimposed in the figure); (b) ANMO results of slope field estimated by method 3.