CN105182419B - A Global Optimization-Based Event Leveling Method for Prestack Seismic Signals - Google Patents

Abstract

Method is evened up the invention discloses a kind of pre-stack seismic signal lineups based on global optimization, by the use of when window center point as seed point, continuous sliding window try to achieve similarity between road collection it is maximum when amount of movement as seed point amount of movement, global optimization is carried out to seed point amount of movement, the prestack road collection amplitude after being evened up as initial data interpolation using seed point amount of movement finally by cubic spline function.The present invention chooses benchmark seismic channel by neighbour's similar propagation, similarity road collection high can be automatically chosen as benchmark road, it is not necessary to manual intervention, compared to directly optional one degree of accuracy is higher in similarity road collection higher from similarity matrix.The present invention designs new object function and constraints after seed point amount of movement is obtained and carries out global optimization to seed point amount of movement, " kick " and wave distortion during solving the problems, such as to even up.

Description

A Global Optimization-Based Event Leveling Method for Prestack Seismic Signals

technical field

The invention belongs to the technical field of seismic exploration, and in particular relates to the design of a global optimization-based pre-stack seismic signal event leveling method.

Background technique

As seismic exploration becomes more and more precise, prestack migration and AVO inversion play a very important role in complex structural imaging and complex lithology reservoir prediction. In order to improve the accuracy of lithology imaging, the pre-stack elastic parameter inversion method combined with pre-stack migration technology is a good choice. As the basis of inversion, pre-stack gather data plays a vital role in the inversion effect. Due to the influence of the anisotropy of the formation medium, there is a large amount of residual moveout in the prestack gather, which leads to the unevenness of the event in the prestack gather. The unevenness of the event will lead to inaccurate imaging and affect the inversion effect.

The commonly used pre-stack gather leveling methods are mainly divided into two categories: the leveling method based on velocity adjustment and the leveling method based on statistical effect. The flattening method based on statistical effects refers to first establishing the objective function of the L2 norm, the objective function is mainly represented by AVO I or III, and then using the sliding time window in the time direction to generate a moving solution for each track, and finally minimizing the objective function, The corresponding solution is the optimal solution. The flattening method based on velocity adjustment assumes that the non-flattened seismic waves in the original prestack seismic gather are mainly caused by residual moveout (RMO), so high-precision velocity estimation using the second-order or fourth-order RMS velocity field can flatten the event of the gather .

1. Leveling method based on speed adjustment

The velocity-adjusted gather flattening method assumes that the unflattened gathers in the original pre-stack seismic trace after noise reduction are mainly caused by residual moveout (RMO), so high-precision velocity estimation using the second-order or fourth-order RMS velocity field can be flattened The event of the gather. The RMO correction is based on formula (1):

where τ is the dynamic correction, x is the offset, t is the time at zero offset, V _ref is the reference velocity function, and V is the update velocity. Then, the development of seismic exploration made the distance between the seismic emission point and the receiving point farther and farther. This reason makes it increasingly difficult to model velocity using RMO curves at far offsets. After improvement, the Alkhalifah time difference model is used, and the model formula is as formula (2):

For the high-order remaining time difference, the time difference δt is determined by the speed difference δV and the travel time δη, as shown in formula (3):

In order to improve the infinitesimal offset distance, set ζ=x/Vt ₀ and the infinitesimal velocity Δ=δV/V, and substitute into formula (3) to get:

δV and δη are obtained by minimizing the error between the input data and the time shift δt, and the number of iterations can be set by yourself.

2. Leveling method based on statistical effects

Hinkley proposed a dynamic gather flattening method (DGF) in 2004, which is a statistical gather flattening method. First of all, this method is one-to-one mapping in the process of processing, that is, each output sample point data is obtained by processing the input data at the same time point in the same track, and can be expressed more conveniently by formula (5):

D _a (t,x)=D _b {(t+m(t,x)),x} (5)

Among them, x is the offset distance, which can also be regarded as the number of gathers sorted from small to large in this gather leveling method; T is time, and a and b represent the data after leveling and before leveling, respectively. By opening the time window for the horizontal offset and the vertical time, moving track by track makes the 2-norm between the two tracks the smallest, that is, to solve formula (6):

The time shift τ _ij of any two channels can be obtained. In the offset direction, 5 tracks are randomly selected as a group, where T ₁ indicates the movement amount between the first track and the second track, T ₂ indicates the movement amount between the third track and the first track, and T ₃ indicates the movement amount between the first track and the second track. The amount of movement between the fourth track and the first track, T ₄ means the amount of movement between the fifth track and the first track. Ten moving amounts can be obtained from 5 traces of data, that is, there is a moving amount between any two traces, and the above 4 moving amounts can be obtained in the least square sense. The calculation formula is shown in formula (7):

T ₁ =(2T _1,2 +T _1,3 +T _1,4 +T _1,5 -T _2,3 -T _2,4 -T _2,5 )/5

T ₂ ＝(T _1,2 +2T _1,3 +T _1,4 +T _1,5 +T _2,3 -T _3,4 -T _3,5 )/5

T ₃ ＝(T _1,2 +T _1,3 +2T _1,4 +T _1,5 +T _2,4 +T _3,4 -T _4,5 )/5

T ₄ =(T _1,2 +T _1,3 +T _1,4 +2T _1,5 +T _2,5 -T _3,5 -T _4,5 )/5 (7)

In prestack gather optimization, although the event can be basically flattened through dynamic correction, due to some factors, such as the inaccurate dynamic correction caused by the fluctuation of the surface, the time-distance curve caused by the horizontal layered isotropic medium The error of the equation. Due to the existence of these effects, the pre-stack gather events are still uneven, and further fine-leveling is required.

Contents of the invention

The purpose of the present invention is to solve the problem that in the prior art, due to the influence of the anisotropy of the formation medium, there is a large amount of residual moveout in the pre-stack gather, which leads to uneven events in the pre-stack gather, which in turn leads to inaccurate imaging effects, thus affecting the inversion To solve the problem of the effect, a global optimization-based event leveling method for pre-stack seismic signals is proposed.

The technical solution of the present invention is: a method for leveling the event of pre-stack seismic signals based on global optimization, comprising the following steps:

S1. Initialize gather leveling parameters;

S2, select the reference track;

S3. Calculating the movement amount of the seed point;

S4. The events are leveled.

Further, the gather flattening parameters in step S1 include the time window size of the prestack gather, the window shift amount, the search radius and the quantile threshold of the similarity matrix.

Further, step S2 includes the following sub-steps:

S21. Calculate the similarity matrix C of any two gathers;

S22. Initialize the attractiveness matrix and the belongingness matrix based on the similarity matrix;

S23. Iteratively updating the attractiveness matrix and the belongingness matrix;

S24. Calculate the gather k that maximizes the sum of the attractiveness matrix and the belongingness matrix;

S25, judging whether the number of iterations reaches the specified number of times, if so, enter step S3, otherwise enter step S26;

S26. Judging whether the gather k is consistent with the result of the last iteration, if so, proceed to step S3, otherwise return to step S23.

Further, step S3 includes the following sub-steps:

S31. Calculate the maximum similarity matrix C _max of the two gathers and define the optimal movement matrix S;

S32. Calculate the value c _m of the matrix C _max corresponding to the quantile threshold;

_S33 . Count the number of gathers greater than the cm value, select the row with the largest number, and the movement amount corresponding to the similarity of this row is the movement amount of the current time window seed point;

S34, interpolating to obtain the movement amount of the seed point of the remaining time window;

S35. Perform global optimization on the movement amount of the seed point in each time window.

Further, step S35 includes the following sub-steps:

S351. Calculate the movement matrix with the largest similarity between gathers;

S352. Calculate the displacement difference matrix in the horizontal direction and the vertical direction;

S353. Determine whether the displacement difference matrix in the horizontal direction and the vertical direction satisfies the constraint condition, and if so, proceed to step S4; otherwise, select a movement amount with a suboptimal similarity to form a new movement amount matrix and return to step S352.

Further, in step S4, the method of interpolating the movement amount matrix by cubic spline interpolation is used to achieve event leveling.

The beneficial effects of the present invention are: the present invention uses the center point of the time window as the seed point, and continuously slides the time window to obtain the movement amount when the similarity between gathers is the largest as the movement amount of the seed point, and globally optimizes the movement amount of the seed point, Finally, the pre-stack gather amplitude after flattening is obtained by interpolating the movement of the seed point as the original data through the cubic spline function, which can achieve the following beneficial effects:

(1) Selecting the reference seismic trace through the similarity propagation of the nearest neighbor can automatically select the gather with high similarity as the reference trace without manual intervention. higher degree.

(2) After obtaining the movement amount of the seed point, design a new objective function and constraint conditions to optimize the movement amount of the seed point globally, and solve the problems of "sudden jump" and waveform distortion in the leveling process.

Description of drawings

Fig. 1 is a flowchart of a global optimization-based pre-stack seismic signal event leveling method provided by the present invention.

Fig. 2 is a sub-step flowchart of step S2 of the present invention.

Fig. 3 is a sub-step flowchart of step S3 of the present invention.

Figure 4 is a graph showing the relationship between similarity and movement.

FIG. 5 is a sub-step flowchart of step S35 of the present invention.

Figure 6 is a cross-sectional view of an unflattened gather.

Fig. 7 is a cross-sectional view of a flattened gather without optimization of seed point movement.

Fig. 8 is a flattened cross-sectional view of the gather after the global optimization of the movement amount of the seed point.

detailed description

Embodiments of the present invention will be further described below in conjunction with the accompanying drawings.

The present invention provides a method for leveling the event of pre-stack seismic signals based on global optimization, as shown in Figure 1, comprising the following steps:

S1. Initialize gather leveling parameters.

Gather flattening parameters include prestack gather time window size N _w , window shift N _s , search radius (less than window shift) N _r and similarity matrix quantile threshold Toi.

Usually, the number of data samples in the time direction of the pre-stack gather is thousands of samples, and the time window is taken in the time direction, and the window can move up and down. The number of data points contained in the time window is recorded as the time window size N _w . The size of N _w is taken from the two cycles of the gather waveform, and one of the gathers can be extracted arbitrarily for Fourier transform, and the main frequency is estimated to determine the number of data points in one cycle.

The range in which the window moves up and down does not exceed the window movement amount N _s . The window moving amount N _s usually takes the size of a cycle data amount.

The search radius (less than the window movement amount) N _r refers to the maximum movement amount when the time window slides up and down.

The similarity matrix quantile threshold Toi, before introducing how to set this parameter, briefly introduce the quantile. The quantile comes from every equal interval point in the data set. Before calculating the quantile, it is necessary to calculate the quantile in the data set. The numbers are sorted from largest to smallest. For example, the 2-quantile is to divide the data into two parts, and the data volume of the data points corresponding to less than and greater than the 2-quantile each accounts for half of the entire data set. Similarly, it can be seen that the 4-quantile corresponds to 3 data point, these 3 data points divide the data set into four equal parts, so that each part of the data represents a quarter of the entire data distribution. The quantile threshold Toi of the similarity matrix usually takes the third 4-quantile value of the similarity matrix, that is, the value is 0.75. For the pre-stack gathers with poor gather quality, Toi can be taken between 0.8-0.9.

S2. Select a reference track.

Before performing gather leveling, it is first necessary to determine the reference gather. In the pre-stack gather profile, a profile usually has dozens of traces, and one of these dozens of traces needs to be selected as the reference trace.

The selection of reference traces depends on the calculation of the similarity between gathers, and iteratively obtains reference traces by using the nearest neighbor similarity propagation, as shown in Figure 2, including the following steps:

S21. For each trace of the pre-stack gather profile starting from time 0 ms, take the time window size N _w. as the length to intercept, and calculate the similarity between any two gathers according to formula (8):

Among them, j ₁ and j ₂ represent the labels of different traces, and the value range is 1...N, where N is the trace number of the gather. The calculated matrix C is the similarity matrix.

S22. Based on the similarity matrix C, define the attractiveness matrix R(i,k) and the belongingness matrix A(i,k) to represent two types of information between gathers, where R(i,k) is from the Set i points to candidate reference track k, which reflects the suitability of gather k to be the reference gather of gather i; A(i, k) points from candidate reference track k to gather i, which reflects the Choose k as the suitability of its reference track.

In order to facilitate the calculation, each element in the matrix C is given a negative sign, and two objective functions, the attractiveness matrix R(i,k) and the belongingness matrix A(i,k), are defined, and R(i,k) and A are iteratively solved (i,k). The size of the attractiveness and belongingness matrix is the same as that of the similarity matrix C, both of which are N×N. Initialize these two matrices as a 0 matrix, define the attractiveness objective function and the belongingness objective function for iterative update, and define the function by the formula (9), (10), (11), (12) indicate:

R(k,k)=P(k)-max{A(k,i')+C(k,i')} (11)

Among them, P(k) takes the mean value of matrix C.

S23. Iteratively updating the attractiveness matrix and the belongingness matrix through the formula (13):

In the formula, λ(0<λ<1) is the convergence coefficient (damping coefficient), which is mainly used to adjust the convergence speed of the algorithm and the stability of the iterative process. The larger the λ, the less obvious the oscillation phenomenon becomes, which means the better the effect Good, but the disadvantage is that the convergence speed will be slower. On the contrary, if the λ is reduced, although the convergence speed is fast, the oscillation phenomenon is obvious. The subscripts old and new represent the final results of the last and current update messages respectively. The iteration termination condition is that the number of iterations is greater than a preset number of times or that R(i, k) and A(i, k) remain unchanged after a number of iterations.

S24. Calculate the gather k that maximizes the sum of the attractiveness matrix and the belongingness matrix by formula (14):

In the formula, if i=k, then gather i is the reference gather of the pre-stack gather profile; if i≠k, then gather k is the reference gather of gather i.

S25 , judging whether the number of iterations reaches the specified number, if so, go to step S3 , otherwise go to step S26 .

S26. Judging whether the gather k is consistent with the result of the last iteration, if so, proceed to step S3, otherwise return to step S23.

In order to make the attractiveness matrix and belongingness matrix converge quickly, a shrinkage factor ρ is introduced, and the definition of ρ is shown in formula (15):

in, for the constructor variable. Update the two matrices R(i,k) and A(i,k) to get:

In the embodiment of the present invention, The value is 4.1, so the value of ρ can be calculated to be 0.729.

S3. Calculate the moving amount of the seed point.

After the reference seismic trace is selected, some data points are selected as seed points in the prestack gather profile and the movement of the seed points is calculated. Take the time window size N _w as the length, open the time window in the time direction, slide the time window with the time window movement N _s as the size until the time window covers all the sample points in time, and use the center point of each time window as the seed point .

As shown in Figure 3, the calculation of the seed point movement mainly includes the following steps:

S31. Record the moving amount when the similarity between gathers is the largest, and define the optimal moving amount matrix S, and the calculation of S is shown in formula (18):

Let i _s =(k-1)N _m , define the set In the formula, g _l is a vector with defined length N _s +1, the i-th element of g _l g _l (i)=D(i _s +1+i-1,j ₂ ), D is the input prestack earthquake Gather data.

The time window of gather j ₁ is fixed, and the time window of gather j ₂ is moved up and down, each movement interval is one data point, and the moving boundary does not exceed the search radius N _r . During the search process, the movement amount of the time window with the highest similarity is recorded, and the similarity is calculated according to formula (8), and the matrix C _max recording the maximum similarity is obtained.

S(j ₁ , j ₂ ) is an asymmetric matrix, such as: S(2,1) and S(1,2), S(2,1) represents the 1 gather obtained by moving gather 1 from fixed gather 2 The optimal moving amount of , and S(1,2) represents the optimal moving amount obtained by fixing gather 1 and moving gather 2. The elements on the diagonal are all 0, because the elements on the diagonal represent the autocorrelation of the gather itself. Obviously, the similarity is the largest when the time window does not perform any movement. For a more vivid description, Figure 4 shows the relationship between the similarity and the amount of movement, where the abscissa represents the amount of movement, the amount of movement to the left is a negative number, the amount of movement to the right is a positive number, and the ordinate is the similarity coefficient . It can be seen from the figure that the correlation coefficient reaches the maximum when the gather 2 moves upward for 3 sampling intervals, and the recorded movement amount is the movement amount of the seed point.

S32. Calculating the value cm of the matrix C _max corresponding to the _quantile threshold. For example, assuming that there are 6 pre-stack gathers to obtain the C _max matrix as shown in the table below, the calculated _cm = 0.85 for this matrix.

Number of channels 1 2 3 4 5 6 Number 1 1 0.9 0.71 0.35 0.89 0.87 4 2 0.9 1 0.76 0.43 0.95 0.96 4

3 0.71 0.76 1 0.78 0.72 0.69 1 4 0.35 0.43 0.78 1 0.64 0.41 1 5 0.89 0.95 0.72 0.64 1 0.93 4 6 0.87 0.96 0.69 0.41 0.93 1 4

Count the number of _gathers greater than the value of cm, and select the row with the largest number, which means that this group of gathers has the strongest similarity, and the gathers are very similar, and the similarity between gathers can be enhanced by moving . For the case where the largest number of rows is not unique, one row can be selected from the largest row. According to the above table, the first row can be selected, and S(1,1), S(1, 2), S(1,5), S(1,6) are reserved, others are discarded. The movement amount corresponding to the similarity of this row is used as the movement amount of the seed point of the current time window, so the movement amount of the seed point 1, 2, 5, and 6 in the above table has been determined.

S34, obtain according to formula (19) interpolation for the seed point that does not determine moving amount:

Among them, j _b represents the subscript of the selected gather, and j _min and j _max respectively represent the track number with the largest moving amount and the smallest track number. Represents the selected gather set. Therefore, channels 3 and 4 in the above table can be obtained by interpolating the movement of the other four channels.

Through the above steps, the movement amount of the center point of each time window can be obtained, and the movement amount of the seed point on the entire pre-stack gather section is composed into a matrix m. The shifts obtained at this time are rough, and if these shifts are used to move the pre-stack gathers to obtain a "jump" phenomenon in the event after the flattening. There are two reasons for this phenomenon: on the one hand, the calculation of similarity is to measure the similarity of a whole segment of signals, and finding the movement with the strongest similarity does not mean that the event in this segment of the waveform Significant (larger amplitude) positional alignment. On the other hand, part of the movement of the seed points is obtained by interpolation, which cannot accurately reflect the degree of event offset between the gather and the reference trace, so these movements cannot be directly used as the movement of the seed points. Optimize the amount of movement to get more accurate movement.

S35. Perform global optimization on the movement amount of the seed point in each time window.

In addition to the above-mentioned "jump" phenomenon, in the time direction, if the difference in the movement amount of the seed point at the center of the adjacent time window exceeds the search radius, it will also cause the waveform to be compressed or stretched too much in the final implementation of gather leveling, resulting in waveform distortion at within the unacceptable range. Therefore, it is necessary to globally optimize the movement amount of the seed point in each time window, as shown in Figure 5, and the specific steps are as follows:

S351. Calculate the movement matrix with the largest similarity between gathers according to formula (20):

Among them, g _x and g _y are different pre-stack gather waveforms in the same time window, Δl is the movement amount, and the arg symbol indicates that the function value takes the dependent variable.

S352. Calculate the displacement difference matrix in the horizontal direction and the vertical direction and

S353, judging whether the displacement difference matrix in the horizontal direction and the vertical direction satisfies the constraint condition, the constraint condition is shown in formula (21):

In the formula, st represents the constraint condition, R represents the threshold value, and |||| _∞ is the matrix infinite norm symbol, which means to find the largest element of the matrix.

If so, go to step S4, otherwise select the movement amount with the second best similarity to form a new movement amount matrix and return to step S352.

The constraints ensure that the difference between the movement of adjacent seed points in the horizontal and vertical directions does not exceed the search radius N _r . In this way, a globally optimized spatial correction factor, that is, the amount of movement in the horizontal direction and the vertical direction, can be obtained.

Figure 6 shows the profile of the original pre-stack seismic gather. The event flatness of the entire profile is greatly improved after the gather is flattened using the gather flattening algorithm, but some events in the circle appear to "jump" Phenomenon, the original one event has become two events, as shown in Figure 7, after the global optimization of the movement of the seed point, Figure 8 is obtained, and the part in the circle can be seen that the event "jumps" while flattening Problem solved.

S4. The events are leveled.

After calculating the movement amount matrix m of the seed point, assuming that the center point coordinate of the qth time window is i _q , the corresponding movement amount is q=1,2,...,Q, define the set Define a two-dimensional array X ^new with the same size as the pre-stack gather data to store the movement of each data sample point in the pre-stack gather, which can be obtained by interpolation, as shown in formula (22):

X ^new stores the stretched sampling coordinates, the matrix D _new of the same size represents the pre-stack gather after event flattening, the matrix X represents the coordinates of the original pre-stack gather, and the matrix D represents the amplitude of the original pre-stack gather. D _new is obtained by performing cubic spline interpolation on D, X, and X ^new . The interpolation is done one by one. For any given one j, the jth column of X is used as the independent variable, and the jth column of D is used as the function value to construct the cubic spline interpolation function as shown in formula (23)

The cubic spline difference can get a smoother result. The interpolation principle is to divide the given n+1 data points into n intervals. The cubic spline equation satisfies the following conditions:

(1) In each segment interval [x _i , x _i+1 ] (i=0,1,...,n-1, x increases), S(x)=S _i (x) is a cubic polynomial.

(2) S( _xi )=y _i (i=0,1,...,n) is satisfied.

(3) The first-order derivative S'(x) and the second-order derivative S″(x) of S(x) are continuous in the [a,b] interval, that is, the S(x) curve is smooth. So n A cubic polynomial piecewise can be written as:

S _i (x)＝a _i +b _i (xx _i )+c _i (xx _i ) ² +d _i (xx _i ) ³ (24)

Among them, a _i , b _i , c _i , and d _i represent 4n unknown coefficients, and the specific expressions are as follows:

m _i is the curve coefficient of the spline.

Those skilled in the art will appreciate that the embodiments described here are to help readers understand the principles of the present invention, and it should be understood that the protection scope of the present invention is not limited to such specific statements and embodiments. Those skilled in the art can make various other specific modifications and combinations based on the technical revelations disclosed in the present invention without departing from the essence of the present invention, and these modifications and combinations are still within the protection scope of the present invention.

Claims (3)

1. A pre-stack seismic signal event leveling method based on global optimization, characterized in that, comprising the following steps: S1. Initialize gather leveling parameters; S2, select the reference track; The step S2 includes the following sub-steps: S21. Calculate the similarity matrix C of any two gathers; S22. Initialize the attractiveness matrix and the belongingness matrix based on the similarity matrix; S23. Iteratively updating the attractiveness matrix and the belongingness matrix; S24. Calculate the gather k that maximizes the sum of the attractiveness matrix and the belongingness matrix; S25, judging whether the number of iterations reaches the specified number of times, if so, enter step S3, otherwise enter step S26; S26. Determine whether the gather k is consistent with the result of the last iteration, if so, enter step S3, otherwise return to step S23; S3. Calculating the movement amount of the seed point; The step S3 includes the following sub-steps: S31. Calculate the maximum similarity matrix C _max of the two gathers and define the optimal movement matrix S; S32. Calculate the value c _m of the matrix C _max corresponding to the quantile threshold; _S33 . Count the number of gathers greater than the cm value, select the row with the largest number, and the movement amount corresponding to the similarity of this row is the movement amount of the current time window seed point; S34, interpolating to obtain the movement amount of the seed point of the remaining time window; S35. Globally optimize the movement amount of the seed point in each time window; Described step S35 comprises following sub-steps: S351. Calculate the movement matrix with the largest similarity between gathers; S352. Calculate the displacement difference matrix in the horizontal direction and the vertical direction; S353, judging whether the displacement difference matrix in the horizontal direction and the vertical direction satisfies the constraint condition, if so, enter step S4, otherwise select the movement amount with the suboptimal similarity to form a new movement amount matrix and return to step S352; S4. The events are leveled. 2. The pre-stack seismic signal event flattening method according to claim 1, characterized in that the gather flattening parameters in the step S1 include the pre-stack gather time window size, window movement, search radius and similarity matrix Quantile threshold. 3. The method for event leveling of pre-stack seismic signals according to claim 1, characterized in that, in the step S4, cubic spline interpolation is used to interpolate the movement matrix to achieve event leveling.