A method and system for enhancing the stability of mirror boundary

Abstract

The present invention provides a method and system for enhancing the stability of mirror boundary, which belongs to the field of seismic exploration data processing. The method for enhancing the stability of mirror boundary uses the characteristics of data edge signals to calculate its mirror signal, and self-learns the seismic data after edge expansion to form big data after edge expansion, and then removes the edge expansion after various numerical operations on the big data after edge expansion to obtain the processed result track without boundary effect. The present invention can obviously remove the boundary effect and enhance the stability of mirror boundary, which has obvious advantages over the traditional zero value and random noise boundary.

Description

Technical Field

[0001] The invention belongs to the field of seismic exploration data processing, and in particular relates to a method and a system for enhancing the stability of a mirror image boundary. Background Art

[0002] In geophysical signal processing, due to the limited nature of the computational domain, appropriate artificial boundary conditions must be used at the free boundaries. Artificial boundaries involve almost every aspect of seismic signal processing. Inappropriate artificial boundary conditions will reduce the calculation accuracy, disrupt the calculation results, and even cause calculation divergence. This is particularly important for numerical simulation wave fields and various data processing methods.

[0003] Generally speaking, non-reflecting boundary conditions are the most widely used, which can be roughly divided into absorbing boundaries, characteristic analysis boundaries and asymptotic solution boundary conditions. Among them, absorbing boundaries are widely used. Scientists have done a lot of work in this regard. The most typical ones are Newton damping absorbing boundary conditions and PML boundary conditions. Israeli & Orszag proposed the classic Newton damping absorbing boundary conditions, which are simple in form and widely used, but in actual use, their accuracy and stability are average, and they strongly depend on the experience of users. PML is a non-reflecting boundary technology developed in recent years. Hu extended the PML boundary conditions to the full Euler equation, and then Parrish & Hu developed PML to be applicable to the case of inviscid and oblique mean flow. Based on previous research, Hu et al. established a conservative nonlinear PML equation applicable to the full Euler and Navier-Stok Lines equations, considering the wide application of the non-conservative form of the Euler equation. Considering the wide application of the non-conservative form of the Euler equation, Lin et al. developed the corresponding non-conservative nonlinear PML equation, further expanding the scope of application of the PML boundary conditions. The PML boundary condition has good stability and accuracy, but due to its complex form, it needs to introduce additional auxiliary variables and equations, which will increase the computing memory and calculation amount in the application.

[0004] Chinese patent publication CN102023077A discloses a supersonic axisymmetric boundary layer wind tunnel, which includes: a front transition section, used to introduce airflow and perform a first-stage rectification on the airflow; a stabilization section, connected to the downstream of the transition section, used to perform a second-stage rectification on the airflow; and a nozzle test section, connected to the downstream of the stabilization section, and having a fan-shaped axisymmetric structure relative to a rotation axis, including an outer peripheral wall, an inner peripheral wall, and two side walls connecting the outer peripheral wall and the inner peripheral wall, wherein at least a transparent window is formed on the side wall, and the nozzle test section includes: a nozzle portion, which is configured as a single-sided expansion nozzle structure, and a test portion, which is located in the nozzle. Downstream of the nozzle section, the axisymmetric boundary layer flow field formed in the nozzle test section can be observed and optically non-contact tested to facilitate the study of the flow field characteristics of the supersonic boundary layer; Chinese patent publication CN102792186A discloses reverse time migration with absorbing boundaries and random boundaries, which explains that images related to the underground area can be generated at least in part based on backward-propagating source wave fields and receiver wave fields, and the source wave field can be forward-propagated from the initial wave field state to the final time state in time through the earth model associated with the underground area, and the backward-propagating source wave field can be forward-propagated from the final time state to the initial time state. The initial wavefield state is determined by propagating the source wavefield backward in time from the final time state to reconstruct the initial wavefield state, and the receiver wavefield can be propagated through the earth model from the final time state, and the earth model can include at least one boundary region that can be defined as having one or more of absorption characteristics, boosting characteristics, random perturbation characteristics, and / or other characteristics, so that the wavefield can be damped, amplified, randomly scattered, and / or otherwise changed at the at least one boundary region, and these wavefields can be used to construct an image of the underground area with an improved signal-to-noise ratio; Chinese patent publication CN101771789A discloses a scanning method for realizing boundary positioning The invention discloses a scanning device and a method for realizing boundary positioning thereof, which describes a method for realizing boundary positioning of a scanning device, and a scanning device which effectively combines optical, electronic and mechanical technologies under this method and can realize intelligent judgment of boundary positioning of scanned files. The scanning device comprises a regular color ribbon and a sensor module capable of detecting the color ribbon, and the file to be scanned is placed between the color ribbon and the sensor module; the sensor module first scans the color ribbon information and stores it in a storage module of the scanning device; the file to be scanned is placed between the sensor module and the color ribbon, and the sensor can obtain the boundary position of the file to be scanned by comparing and measuring the portion of the color ribbon covered by the file to be scanned.

[0005] In geophysical signal processing, due to the limited data, there are often boundary problems. If no boundary is added or the boundary is not added ideally, the signal processing result may be distorted or deformed, resulting in serious boundary effects. Directly adding a mirror boundary greatly improves the boundary effect, but if there is an abnormal signal at the boundary, there will also be an abnormal signal after the mirror, which will lead to poor processing results. Summary of the invention

[0006] The purpose of the present invention is to solve the problems existing in the above-mentioned prior art, and to provide a method and system for enhancing the stability of the mirror boundary, which utilizes the characteristics of the data edge signal to calculate its self-learning mirror signal for edge expansion, improves the traditional boundary mode, optimizes the mirror signal edge expansion, and achieves better processing effect.

[0007] The present invention is achieved through the following technical solutions:

[0008] The first aspect of the present invention provides a method for enhancing the stability of mirror boundaries. The method uses the characteristics of data edge signals to calculate their mirror signals, and performs self-learning on the expanded seismic data to form expanded big data. After performing various numerical operations on the expanded big data, the expanded edges are removed to obtain a processed result set without boundary effects.

[0009] A further improvement of the present invention is that the method comprises:

[0010] (1) Input the original seismic signal shot gather or gather data, that is, the original matrix S;

[0011] (2) Construct a large matrix S1;

[0012] (3) Obtain the value of each point in the large matrix S1;

[0013] (4) Numerical calculation is performed on the large matrix S1 to obtain the calculated result, and then the expanded edge is removed from the calculated result to obtain the result gather.

[0014] A further improvement of the present invention is that the size of the original matrix S in step (1) is m×n, wherein m is the number of sampling points and n is the number of seismic traces.

[0015] A further improvement of the present invention is that the operation of step (2) comprises:

[0016] Add H points outward from the top and bottom of the sampling point direction of the original matrix S to form a top edge expansion blank area and a bottom edge expansion blank area. The size of the top edge expansion blank area and the bottom edge expansion blank area are both H*n;

[0017] Then, K points are added outwards on the left and right sides of the seismic trace direction to form a left-side expansion blank area and a right-side expansion blank area. The sizes of the left-side expansion blank area and the right-side expansion blank area are both K*(m+2*H);

[0018] The size of the large matrix S1 obtained after edge expansion is (m+2*H)×(n+2*K).

[0019] A further improvement of the present invention is that the operation of step (3) comprises:

[0020] (31) Mirror the top edge expansion blank area and the bottom edge expansion blank area along the sampling point direction, and remove the distortion points through self-learning:

[0021] (32) The values ​​of the left and right expansion blank areas are obtained by mirroring along the seismic trace direction, and the distortion points are removed by self-learning.

[0022] A further improvement of the present invention is that the operation of step (31) includes:

[0023] Put S into the center of the large matrix S1, extract the i-th channel in S, and assume that S i,j S is a point on the top boundary of the i-th track; for each track, the values ​​of each point in the top expanded blank area are obtained by mirroring along the sampling point direction using formula (7):

[0024] S i,j-k =S i,j+k k≤H (7)

[0025] Similarly, obtain the value of each point in the bottom margin blank area;

[0026] A model channel is established through a sliding window, and the similarity coefficient between each channel and the model channel is calculated. When the similarity coefficient is less than a first threshold, the value of the channel is replaced by the value of the model channel, that is, the distortion point is removed;

[0027] The calculation of the boundaries of all channels in the direction of the sampling point is completed. At this time, the size of S is (m+2*H)×n.

[0028] A further improvement of the present invention is that the operation of step (32) includes:

[0029] Put S into the center of the large matrix S1, extract the value of each seismic trace at the tth moment in S, and assume that S t,p is a point on the leftmost boundary of S at time t; at each moment, the value of each point in the left expanded blank area is obtained by mirroring along the seismic trace direction using formula (8):

[0030] S t,p-q =S t,p+q q≤K (8)

[0031] Similarly, obtain the values ​​of each point in the right expanded blank area;

[0032] A model channel is established through a sliding window, and a similarity coefficient between each channel and the model channel is calculated. When the similarity coefficient is less than a second threshold, the value of the channel is replaced by the value of the model channel, that is, the distortion point is removed;

[0033] The calculation of the boundaries at all times in the direction of the seismic trace is completed, and the size of S is now (m+2*H)×(n+2*K).

[0034] A further improvement of the present invention is that the first threshold is 0.3;

[0035] The second threshold is 0.1.

[0036] A further improvement of the present invention is that the operation of removing the edge expansion from the result of the calculation in step (4) includes:

[0037] Taking the center point of the large matrix S1 as the center point of the result matrix, a matrix with the same number of traces and sampling points as the original matrix S is extracted from the expanded large matrix S1. The extracted matrix is ​​the result trace gather.

[0038] A second aspect of the present invention provides a system for enhancing the stability of mirror boundaries, the system comprising: a memory, a processor, and a computer program stored in the memory, the computer program executing the following steps when executed by the processor:

[0039] (1) Input the original seismic signal shot gather or gather data, that is, the original matrix S;

[0040] (2) Construct a large matrix S1;

[0041] (3) Obtain the value of each point in the large matrix S1;

[0042] (4) Numerical calculation is performed on the large matrix S1 to obtain the calculated result, and then the expanded edge is removed from the calculated result to obtain the result gather.

[0043] Compared with the prior art, the beneficial effects of the present invention are: the present invention can significantly remove the boundary effect, enhance the stability of the mirror boundary, and has obvious advantages over the traditional zero value and random noise boundaries. BRIEF DESCRIPTION OF THE DRAWINGS

[0044] Figure 1 Schematic diagram of building a large matrix

[0045] Figure 2 A flowchart of the steps of the method of the present invention;

[0046] Figure 3-1Original simulation gathers;

[0047] Figure 3-2 The gathers without edge denoising;

[0048] Figure 3-3 The gather after adding zero edge denoising;

[0049] Figure 3-4 Gather after adding mirror edge denoising;

[0050] Figure 3-5 Gathers with enhanced mirror boundary stability;

[0051] Figure 4-1 Actual data gathers;

[0052] Figure 4-2 The actual data plus the zero-edge denoising gather;

[0053] Figure 4-3 The actual data plus the gather after image boundary denoising;

[0054] Figure 4-4 Actual data plus enhanced mirror boundary stability denoising. DETAILED DESCRIPTION

[0055] The present invention is further described in detail below in conjunction with the accompanying drawings:

[0056] In geophysical signal processing, due to the limited nature of the computational domain, appropriate artificial boundary conditions must be used at the free boundaries. Artificial boundaries involve almost every aspect of seismic signal processing. Inappropriate artificial boundary conditions will reduce the calculation accuracy, disrupt the calculation results, and even cause calculation divergence. This is particularly important for numerical simulation wave fields and various data processing methods.

[0057] Existing boundary algorithms used in data processing often only alleviate the boundary effect, but do not completely remove the boundary effect. For example, the classic absorbing boundary technology is briefly described as follows:

[0058] The application of classical absorbing boundaries is mainly in wave field simulation and wave field description. The truncation boundary surface of the wave field is generally far away from the earthquake source. Therefore, near the truncation boundary surface, the wave field satisfies the wave equation. Let the operator L be as shown in formula (1):

[0059]

[0060] Where x, y and z are spatial coordinate variables, t is the time variable, v is the wave propagation velocity, and P(x,z,t) is the compression wave field (P wave).

[0061] Then near the truncation boundary surface:

[0062] LP(x,y,z,t)=0 (2)

[0063] Decompose the operator L:

[0064] L=L - L + (3)

[0065] in:

[0066]

[0067] Among them, B. Engquuist et al. have proved that if P is incident on the left outer traveling wave on the left interface x = 0 of the calculation area at any angle, then:

[0068] L - P(x,y,z,t)=0 (5)

[0069] The incident wave will be completely absorbed without causing any reflection, so equation (5) can be used as the accurate absorption boundary condition on the left boundary x = 0. Similarly, if P is a right outgoing wave incident at any angle on the right interface x = x0 of the calculation area

[0070] L + P(x,y,z,t)=0 (6)

[0071] It can be used as an accurate absorbing boundary condition on the right interface x = x0. However, L - and L + They cannot absorb right-traveling waves and left-traveling waves separately, so they are called single-phase wave operators.

[0072] The invention discloses a method for enhancing the stability of a mirror boundary, which is mainly used to solve the boundary problem in data processing and has an improving effect on data distortion caused by data calculation boundary.

[0073] The present invention provides a method for enhancing the stability of a mirror image boundary. The method uses the characteristics of a data edge signal to calculate its mirror image signal, and performs self-learning on the seismic data after edge expansion to form more optimized, stable and smooth edge expansion data, and then performs various numerical operations on the big data after edge expansion. Finally, the edge expansion is removed from the processed result to obtain a processed result gather without boundary effect.

[0074] like Figure 2 As shown, the method of the present invention comprises:

[0075] (1) Input the original seismic signal shot gather or trace gather data S, that is, the original matrix S, the size of S is m sampling points, n seismic traces, that is, the size of S is m×n;

[0076] (2) Construct a large matrix S1:

[0077] For the original matrix S, H points are added outward along the top and bottom of the sampling point direction to form the top expansion blank area and the bottom expansion blank area. The size of the top expansion blank area and the bottom expansion blank area are both H*n. K points are added outward along the left and right sides of the seismic trace direction to form the left expansion blank area and the right expansion blank area. The size of the left expansion blank area and the right expansion blank area are both K*(m+2*H). The top expansion blank area, the bottom expansion blank area, the left expansion blank area and the right expansion blank area are collectively referred to as the expansion blank area. The seismic data after expansion is obtained, that is, the large matrix S1. The size of S1 is (m+2*H)×(n+2*K), such as Figure 1 As shown;

[0078] (3) Mirror expansion is performed on the blank area of ​​expansion along the sampling point direction and the seismic trace direction respectively: the mirror value of the new expanded data in the blank area of ​​expansion is iteratively obtained point by point, and the expanded data is self-learned to obtain the large matrix S1;

[0079] The step (3) is achieved as follows:

[0080] (31) Calculate the values ​​of the top edge expansion blank area and the bottom edge expansion blank area by mirroring along the sampling point direction:

[0081] Put S into the center of the large matrix S1, extract the i-th track in S, and assume that the upper boundary point of the track is j, then S i,j is a point on the upper boundary of the i-th track; for S1, S i,j is a point inside it. For each trace, the value of the blank area of ​​the expanded edge is obtained by mirroring along the sampling point direction using formula (7):

[0082] S i,j-k =S i,j+k k≤H (7)

[0083] The value of each point on the expanded edge is calculated using formula (7) for each channel. The value of each point is actually the mirror image signal of the point near the expanded edge in the original matrix S.

[0084] Similarly, assuming S i,j S is a point on the bottom boundary of the i-th track, and the same method as above is used to mirror and obtain the value of each point in the bottom expanded blank area.

[0085] The expanded seismic data is subjected to self-learning. Specifically, since the seismic data has similarity in the direction of the trace, that is, the sync axis has continuity, a model trace is established through a sliding window (generally, a few traces are taken to calculate the average value as the model trace, and then another few traces are slid to calculate the average value as another model trace). The similarity coefficient between each trace and the model trace is obtained (obtained by a conventional algorithm, which will not be repeated here). According to experience, when the similarity coefficient is less than the first threshold (for example, 0.3) (if it is greater than or equal to 0.3, the original value is kept unchanged), the trace value is replaced by the model trace value, so that the distortion point can be removed.

[0086] Complete the calculation of the boundaries of all channels in the sampling point direction. At this time, the size of S is (m+2*H)×n;

[0087] (32) Calculate the value of the blank area by mirroring along the track direction:

[0088] Put S into the center of the large matrix S1, extract the value of each channel in S at the tth time, and assume that the left boundary point of the channel is p, then S t,p is a point on the leftmost boundary of S at time t; for S1, S t,p is a point inside it. At each moment, the value of the blank area is obtained by mirroring along the track direction using formula (8):

[0089] S t,p-q =S t,p+q q≤K (8)

[0090] The value of each point on the expanded edge is calculated using formula (8) for each trace. The value of each point is actually the mirror image signal of the point near the expanded edge in the original matrix S.

[0091] Similarly, assuming S t,p S is a point on the rightmost boundary of S at the tth moment, and the same method as above is used to mirror and obtain the value of each point in the right expanded blank area.

[0092] Similarly, the expanded seismic data is subjected to self-learning. Specifically, since the seismic data has similarity in the direction of the trace, that is, the posited axis has continuity, a model trace is established through a sliding window, and the similarity coefficient between each trace and the model trace is obtained. According to experience, when the similarity coefficient is less than the second threshold value (for example, 0.1) (if it is greater than or equal to 0.1, the original value remains unchanged), the trace value is replaced by the model trace value, so that the distortion point can be removed.

[0093] The calculation of the boundaries in the road direction at all times is completed, and the size of S is now (m+2*H)×(n+2*K).

[0094] (4) Numerical calculation is performed on the large matrix S1 to obtain the calculated result, and then the expanded edge is removed from the calculated result to obtain the result gather.

[0095] After the edge expansion is completed, the data will become a matrix of size (m+2*H)×(n+2*K). Numerical operations are performed on the matrix, such as solving linear and nonlinear equations, matrix operations, interpolation, numerical integration, fitting, finding extreme values, mathematical transformation and filtering, etc. Different numerical operation methods are used according to different purposes. However, no matter what numerical operation is performed, a matrix of (m+2*H)×(n+2*K) will be obtained after the operation is completed.

[0096] The step (4) of removing the expanded edge from the result after the operation is achieved as follows: taking the center point of the large matrix S1 as the center point of the result matrix, extracting a matrix with the same number of channels and sampling points as the original matrix S from the expanded large matrix S1, and the extracted matrix is ​​the result channel set that needs to be output.

[0097] In the following embodiments, median filtering is used to analyze the algorithm and its application effect.

[0098] Embodiments of the present invention are as follows:

[0099] [Example 1]

[0100] like Figure 3-1 to Figure 3-5 This is the denoising effect of simulated seismic data and enhanced mirror boundary stability. Figure 3-1 This is the original simulation gather, with 2 reflection layers, at 1000ms and 2000ms respectively, 3000 sample points in total, 27 channels, sampling rate of 1ms, and random noise added. Figure 3-2 This is the data set after median filtering denoising without adding boundaries. It can be seen that the signal-to-noise ratio of the denoised data set is improved, but there is signal distortion at the boundaries. Figure 3-3 This is the data set after adding zero boundary denoising. It can be seen that the signal-to-noise ratio after denoising is also improved. The signal at the boundary is still distorted, and the value at the boundary is almost close to zero. Figure 3-4 The denoised data set after adding the conventional mirror boundary. It can be seen that the signal-to-noise ratio of the denoised data set is improved, and the boundary effect at the boundary is lighter. Figure 3-5 For the gather after the mirror image boundary stability is enhanced by the method of the present invention, the boundary effect at the boundary of the gather display is almost gone.

[0101] [Example 2]

[0102] like Figure 4-1 to Figure 4-4 This is the denoising effect of actual seismic data plus enhanced mirror boundary stability. Figure 4-1 This is the original actual gather, and we can see that the random noise is strong. Figure 4-2This is the data set after adding zero boundary and performing median filtering to denoise. It can be seen that the signal-to-noise ratio of the denoised data set is improved, but there is signal distortion at the boundary (in the black frame). Figure 4-3 This is the data set after adding zero boundary denoising. It can be seen that the signal-to-noise ratio after denoising is also improved. The signal at the boundary (in the black box) is still distorted, and the value at the boundary is almost close to zero. Figure 4-4 For the gather after the method of the present invention is used to enhance the stability of the mirror boundary, it can be seen that the signal-to-noise ratio of the gather after denoising is improved, the boundary effect at the boundary of the gather display (in the black frame) is almost gone, and the similarity with the original signal is the best.

[0103] The present invention utilizes the characteristics of data edge signals to calculate the self-learning mirror signal for edge expansion, thereby enhancing stability, improving the traditional boundary mode, optimizing the mirror signal edge expansion, and achieving better processing effect.

[0104] Finally, it should be explained that the above technical solution is only one implementation method of the present invention. For those skilled in the art, it is easy to make various types of improvements or modifications based on the application methods and principles disclosed in the present invention, and it is not limited to the method described in the above specific implementation method of the present invention. Therefore, the method described above is only preferred and does not have a restrictive meaning.

Claims

1. A method for enhancing the stability of a mirror image boundary, characterized in that: The method uses the characteristics of data edge signals to calculate their mirror signals, and performs self-learning on the expanded seismic data to form expanded big data, performs various numerical operations on the expanded big data, and then removes the expanded data to obtain a processed result gather without boundary effects; The method comprises: (1) Input the original seismic signal shot gather or gather data, that is, the original matrix S; (2) Construct a large matrix S1; (3) Obtain the value of each point in the large matrix S1; (4) performing numerical operations on the large matrix S1 to obtain a result after the operation, and then removing the expanded edge from the result after the operation to obtain a result gather; The operation of step (2) includes: Add H points outward from the top and bottom of the sampling point direction of the original matrix S to form a top edge expansion blank area and a bottom edge expansion blank area. The size of the top edge expansion blank area and the bottom edge expansion blank area are both H*n; Then, K points are added outwards on the left and right sides of the seismic trace direction to form a left-side expansion blank area and a right-side expansion blank area. The sizes of the left-side expansion blank area and the right-side expansion blank area are both K*(m+2*H); The size of the large matrix S1 obtained after edge expansion is (m+2*H)×(n+2*K); The operation of step (3) includes: (31) Mirror the top edge expansion blank area and the bottom edge expansion blank area along the sampling point direction, and remove the distortion points through self-learning: (32) Calculate the values ​​of the left and right expansion blank areas in a mirrored manner along the seismic trace direction, and remove the distortion points through self-learning; The operation of step (31) includes: Put S into the center of the large matrix S1, extract the i-th channel in S, and assume that S i,j S is a point on the top boundary of the i-th track; for each track, the values ​​of each point in the top expanded blank area are obtained by mirroring along the sampling point direction using formula (7): S i,j-k =S i,j+k k≤H (7) Similarly, obtain the value of each point in the bottom margin blank area; A model channel is established through a sliding window, and the similarity coefficient between each channel and the model channel is calculated. When the similarity coefficient is less than a first threshold, the value of the channel is replaced by the value of the model channel, that is, the distortion point is removed; Complete the calculation of the boundaries of all channels in the sampling point direction. At this time, the size of S is (m+2*H)×n; The operation of step (32) includes: Put S into the center of the large matrix S1, extract the value of each seismic trace at the tth moment in S, and assume that S t,p is a point on the leftmost boundary of S at time t; at each moment, the value of each point in the left expanded blank area is obtained by mirroring along the seismic trace direction using formula (8): S t,p-q =S t,p+q q≤K (8) Similarly, obtain the values ​​of each point in the right expanded blank area; A model channel is established through a sliding window, and a similarity coefficient between each channel and the model channel is calculated. When the similarity coefficient is less than a second threshold, the value of the channel is replaced by the value of the model channel, that is, the distortion point is removed; The calculation of the boundaries at all times in the direction of the seismic trace is completed, and the size of S is now (m+2*H)×(n+2*K).

2. The method for enhancing the stability of mirror image boundaries according to claim 1, characterized in that: The size of the original matrix S in step (1) is m×n, where m is the number of sampling points and n is the number of seismic traces.

3. The method for enhancing the stability of mirror image boundaries according to claim 1, characterized in that: The first threshold is 0.3; The second threshold is 0.

1.

4. The method for enhancing the stability of mirror image boundaries according to claim 1, characterized in that: The operation of removing the edge expansion from the result of the calculation in step (4) includes: Taking the center point of the large matrix S1 as the center point of the result matrix, a matrix with the same number of traces and sampling points as the original matrix S is extracted from the expanded large matrix S1. The extracted matrix is ​​the result trace gather.

5. A system for enhancing the stability of mirror boundary, characterized in that: The system comprises: a memory, a processor, and a computer program stored in the memory, wherein the computer program executes the following steps when executed by the processor: (1) Input the original seismic signal shot gather or gather data, that is, the original matrix S; (2) Construct a large matrix S1; (3) Obtain the value of each point in the large matrix S1; (4) performing numerical operations on the large matrix S1 to obtain a result after the operation, and then removing the expanded edge from the result after the operation to obtain a result gather; The operation of step (2) includes: Add H points outward from the top and bottom of the sampling point direction of the original matrix S to form a top edge expansion blank area and a bottom edge expansion blank area. The size of the top edge expansion blank area and the bottom edge expansion blank area are both H*n; Then, K points are added outwards on the left and right sides of the seismic trace direction to form a left-side expansion blank area and a right-side expansion blank area. The sizes of the left-side expansion blank area and the right-side expansion blank area are both K*(m+2*H); The size of the large matrix S1 obtained after edge expansion is (m+2*H)×(n+2*K); The operation of step (3) includes: (31) Mirror the top edge expansion blank area and the bottom edge expansion blank area along the sampling point direction, and remove the distortion points through self-learning: (32) Calculate the values ​​of the left and right expansion blank areas in a mirrored manner along the seismic trace direction, and remove the distortion points through self-learning; The operation of step (31) includes: Put S into the center of the large matrix S1, extract the i-th channel in S, and assume that S i,j S is a point on the top boundary of the i-th track; for each track, the values ​​of each point in the top expanded blank area are obtained by mirroring along the sampling point direction using formula (7): S i,j-k =S i,j+k k≤H (7) Similarly, obtain the value of each point in the bottom margin blank area; A model channel is established through a sliding window, and the similarity coefficient between each channel and the model channel is calculated. When the similarity coefficient is less than a first threshold, the value of the channel is replaced by the value of the model channel, that is, the distortion point is removed; Complete the calculation of the boundaries of all channels in the sampling point direction. At this time, the size of S is (m+2*H)×n; The operation of step (32) includes: Put S into the center of the large matrix S1, extract the value of each seismic trace at the tth moment in S, and assume that S t,p is a point on the leftmost boundary of S at time t; at each moment, the value of each point in the left expanded blank area is obtained by mirroring along the seismic trace direction using formula (8): S t,p-q =S t,p+q q≤K (8) Similarly, obtain the values ​​of each point in the right expanded blank area; A model channel is established through a sliding window, and a similarity coefficient between each channel and the model channel is calculated. When the similarity coefficient is less than a second threshold, the value of the channel is replaced by the value of the model channel, that is, the distortion point is removed; The calculation of the boundaries at all times in the direction of the seismic trace is completed, and the size of S is now (m+2*H)×(n+2*K).

Images