An adaptive low-rank seismic data denoising method based on symplectic geometry decomposition
The low-rank seismic data denoising method based on adaptive symplectic geometric decomposition utilizes symplectic QR decomposition and symplectic geometric entropy to select the rank, thus solving the problems of high computational cost and poor preservation of effective signals in seismic data processing and achieving efficient denoising results.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- CHINA PETROLEUM & CHEMICAL CORP
- Filing Date
- 2024-12-20
- Publication Date
- 2026-06-23
AI Technical Summary
Existing seismic data processing methods have high computational costs in the denoising process, making them difficult to apply in actual production, especially when the data volume is large, and the preservation of effective signals is poor.
A low-rank seismic data denoising method using adaptive symplectic geometric decomposition is proposed. The seismic data is transformed into the frequency domain through Fourier transform, a Hankel matrix is constructed and symplectic QR decomposition is performed, the rank is adaptively selected using symplectic geometric entropy, and the data is restored by combining low-rank approximate reconstruction and inverse Fourier transform.
It improves computational efficiency, can adaptively retain more effective signals in complex structures, reduces computational costs, and enhances denoising performance.
Smart Images
Figure CN122260408A_ABST
Abstract
Description
Technical Field
[0001] This invention belongs to the field of geophysical exploration seismic data processing technology, specifically a method for denoising low-rank seismic data using adaptive symplectic geometric decomposition. Background Technology
[0002] In seismic data processing, random noise is frequently encountered, significantly impacting subsequent seismic data analysis. Therefore, random noise suppression is a crucial step before applying seismic data. Numerous methods have been proposed and applied to suppress random noise in real-world seismic data. In recent years, many scholars have conducted in-depth research on low-rank methods for signals, and it is generally believed that seismic signals, after undergoing specific reconstruction techniques, possess low-rank characteristics. These methods demonstrate excellent processing performance; for example, by stacking non-locally similar blocks in seismic signals to construct a matrix, this matrix exhibits low-rank characteristics in a low-dimensional subspace of a high-dimensional space. Oropeza utilizes the conversion of seismic signals into the frequency domain to construct a Hankel matrix with a lower rank. Returning the reconstructed Hamiltonian (Hankel) matrix from the lower rank to the frequency domain yields excellent processing results. The Fxy-mssa method was proposed in Geophysics 76, pp. 25-32 in 2011. However, as this method has been developed, when the data volume is large, the SVD decomposition of the resulting Hankel matrix significantly impacts the computational cost. Furthermore, the choice of rank has constrained the development of this method, making it difficult to apply in practical production.
[0003] The patent "Earthquake Denoising Method Based on Sparse Low-Rank Constraints of Lp Pseudo-norm and γ-norm" (application number 202210534230.0) uses low-rank constraints to denoise time-frequency domain seismic records. After converting the seismic record to the frequency domain, singular value distributions are obtained through SVD decomposition. Then, a γ-norm is added for low-rank constraints, and iterative calculations are performed using the alternating direction multiplier method to obtain the denoised seismic record. However, the low-rank method itself is inefficient due to the need for extensive SVD processing. Furthermore, the method adds alternating direction multipliers for iterative processing after each SVD calculation, significantly increasing the time cost and making it difficult to apply in practical production.
[0004] The patent "Earthquake Denoising and Interpolation Method Based on Weighted Frame Transform-Low-Dimensional Manifold Model" (application number 202210058763.6) also utilizes the low-rank characteristic of the Hankel matrix for denoising and interpolation. It constructs the weights of the wf-LDMM regularization term based on the transformation speed of singular values to constrain the data fidelity term. This concentrates the energy of the updated frame transform coefficients on the first few larger singular values and directly discards the smaller singular values. However, it is difficult to determine the boundary between effective and invalid singular values, which can easily lead to the leakage of effective signals.
[0005] The paper "Earthquake Rank Reduction and Denoising Based on Adaptive Rank Selection Method" published in the Journal of Chengdu University of Technology (Natural Science Edition) in 2020, Vol. 47, No. 06, pp. 742-755, calculates the singular value signal-to-noise ratio threshold by using the ratio of the median of the singular value sequence of the observation matrix to the number of rows and columns of the matrix when performing block calculations on seismic data. Adaptive selection is performed for each singular value decomposition. However, in the selection process, the AMSE framework is established and iterative optimization is performed, which increases the time cost of the paper's method by nearly three times, making it more difficult to apply the method to actual production. Summary of the Invention
[0006] In view of the above problems, the present invention is proposed to provide an adaptive symplectic geometric decomposition method for denoising low-rank seismic data, which overcomes or at least partially solves the above problems.
[0007] To achieve the above objectives, the present invention adopts the following technical solution:
[0008] A method for denoising low-rank seismic data using adaptive symplectic geometric decomposition, the method comprising:
[0009] S1. Use Fourier transform to convert 3D seismic data into the frequency domain;
[0010] S2. Slice the frequency domain seismic data into Hankel matrices according to frequency, and then form the Hankel matrices into Hamiltonian matrix form.
[0011] S3. Use symmetric QR decomposition to find the eigenvalues and eigenvectors of the Hamiltonian matrix;
[0012] S4. Calculate the symplectic geometric entropy of the eigenvalues, and obtain the number of selection ranks from the adaptive symplectic geometric entropy;
[0013] S5. Reconstruct the block Hankel matrix using a low-rank approximation.
[0014] S6. Convert the processed frequency domain to the time domain using inverse Fourier transform.
[0015] Optionally, in step S1, if the size of the three-dimensional seismic data volume is (nt*nx*ny), where nt represents the number of sampling points, i.e., the number of samples in the time dimension; nx represents the number of points in the X-axis direction of the data volume, i.e., the number of samples in the horizontal direction of the seismic data; and ny represents the number of points in the Y-axis direction, i.e., the number of samples in the other horizontal direction of the seismic data; the data is cut into multiple seismic blocks of local size (ntwin*nxwin*nywin) using a windowing function, where ntwin is the number of samples in the window in the nt direction, nxwin is the number of samples in the window in the nx direction, and nywin is the number of samples in the window in the ny direction; then, the time domain data D(t, x, y) is transformed into frequency domain data F(ω, x, y) using a Fourier transform, where t represents time, ω represents frequency, x is the horizontal direction in the plane, and y is the other horizontal direction in the plane.
[0016] Optionally, step S2 includes:
[0017] The frequency component slice F(ω) of each frequency component in the frequency domain data is represented as:
[0018]
[0019] In the formula: nx represents the number of data points along the X-axis, i.e., the number of seismic data samples in the horizontal direction; ny represents the number of data points along the Y-axis, i.e., the number of seismic data samples in the other horizontal direction; ω represents the frequency.
[0020] The Hankel matrix formed by the frequency blocks in the i-th row is shown below:
[0021]
[0022] Where dx represents the embedding dimension in the X-axis direction; H x It is the Hankel matrix constructed along the X-axis; nx represents the number of data volumes along the X-axis, i.e., the number of samples of seismic data in the horizontal direction; ω represents the frequency;
[0023] Then construct a block Hankel matrix along the Y-axis:
[0024]
[0025] In the above formula, dy represents the embedding dimension in the Y-axis direction, and the data size of matrix H(ω) is I*J, where i = dx*dy, J = (nx - dx + 1) * (ny - dy + 1); ny represents the number in the Y-axis direction; ω represents the frequency.
[0026] Construct the Hamiltonian matrix M from the block Hankel matrix H:
[0027]
[0028] Where A is a symmetric matrix and satisfies A = H T H, H T A is the transpose of matrix H. T It is the transpose of A.
[0029] Optionally, step S3 includes:
[0030] After constructing the Hamiltonian matrix, we use the square matrix N of the Hamiltonian matrix, i.e., N = M. 2 We can construct a symplectic orthogonal matrix G such that Where C is the upper Hessenberg matrix, and G represents a symplectic orthogonal matrix, which can be replaced by a given Household matrix Q. According to the symplectic similarity theorem of the symplectic QR algorithm, the matrix [Q 0; 0 Q] is also a symplectic orthogonal matrix, used to replace G. Then, the symplectic similarity theorem is used to ensure that the symplectic orthogonal matrix maintains the structure of the Hamiltonian system, and its eigenvalues follow... Where λ(A) are the eigenvalues of the symmetric matrix A, and λ(C) are the eigenvalues of the upper Hessenberg matrix C.
[0031] Optionally, step S4 includes:
[0032] If the vector formed by the Hamiltonian matrices is s, where the length of s is L;
[0033] The first step is to normalize the eigenvectors:
[0034]
[0035] v j It is the j-th component in the standardized eigenvalue vector; It is the square of the j-th eigenvalue in the original eigenvalue vector s; Let be the sum of the squares of all eigenvalues in the standardized original eigenvalue vector s, k = 1, 2, ..., L, where L is the length of the eigenvalue vector;
[0036] The normalized symplectic geometric entropy e(N) is defined as:
[0037]
[0038] To achieve adaptive selection, the contribution of the i-th feature to the entropy must first be calculated, defined as ce. i , satisfy ce i =e(N)-e(N) -i ), where N -i This represents the i-th feature that was removed;
[0039] Features are ranked based on their relative contribution to entropy; the average of all ce values is defined as c, and their standard deviation is defined as d, where ce i -cd>0, features with high ce values are considered relevant.
[0040] Optionally, step S5 includes:
[0041] The matrix is reconstructed using the eigenvalue order determined by the symplectic geometric entropy. First, the corresponding coefficient matrix is constructed. If the selected eigenrank is m, the low-rank block Hankel matrix is obtained by summing the first m reconstructed matrices, and then the data is reconstructed.
[0042]
[0043] in This refers to the reconstruction matrix of the block Hankel matrix H; Q represents the reconstructed orthogonal matrix after determining the eigenvalue order m through symplectic geometric entropy, and its columns are the eigenvectors of the original data matrix.
[0044] Optionally, step S6 includes:
[0045] The matrix after low-rank reconstruction of the components at each frequency. The new matrix, formed by the diagonal averaging method, is then subjected to inverse Fourier transform to obtain the low-rank time-domain result.
[0046] In summary, due to the adoption of the above technical solution, the beneficial effects of the present invention are:
[0047] 1. This invention utilizes symplectic QR instead of SVD, which not only improves the efficiency of the method, but also proposes an adaptive rank selection method in combination with symplectic geometric entropy, which can adaptively retain more effective signals in complex structures.
[0048] 2. This invention performs QR decomposition (symplectic QR method) within a symplectic geometric framework. Symplectic QR greatly improves computational efficiency by replacing SVD. Furthermore, it provides an optimal rank reduction method by combining symplectic geometric entropy. This invention not only greatly improves computational efficiency but also adaptively retains more effective signals in complex structures.
[0049] 3. This invention extracts the eigenvalues and eigenvectors of the Hankel matrix by constructing the Hankel matrix and applying symplectic QR decomposition. It then uses symplectic geometric entropy to adaptively determine the number of ranks, thereby achieving low-rank approximate reconstruction, effectively denoising and preserving the effective signals in the seismic data. Attached Figure Description
[0050] Figure 1 This is a schematic flowchart of an adaptive symplectic geometric decomposition method for denoising low-rank seismic data, provided in an embodiment of this application.
[0051] Figure 2 This is a schematic diagram of a noisy linear model in an embodiment of this application.
[0052] Figure 3 This is a schematic diagram of the fxy-sgmd method for denoising a linear model according to an embodiment of this application.
[0053] Figure 4 This is a schematic diagram of the denoising residual of the linear model using the fxy-sgmd method in an embodiment of this application.
[0054] Figure 5 This is a schematic diagram comparing the time of different methods under different volumes of the linear model in the embodiments of this application.
[0055] Figure 6 This is a schematic diagram of a noisy curve model in an embodiment of this application.
[0056] Figure 7 This is a schematic diagram of the fxy-sgmd (rank ratio method) noise reduction method for the curve model in the embodiment of this application.
[0057] Figure 8 This is a schematic diagram of the denoising residual of the curve model in the embodiment of this application using the fxy-sgmd (rank ratio method).
[0058] Figure 9 The noise reduction effect diagram after processing in the embodiments of this application is shown.
[0059] Figure 10 The residual diagram of the processed data and the original data in the embodiments of this application.
[0060] Figure 11 This is pre-stack time migration data for work area A1 in region Z.
[0061] Figure 12 This application provides an example of a schematic diagram illustrating an insufficient selection of eigenvalue rank.
[0062] Figure 13 This is a schematic diagram illustrating the normal selection of the rank of the feature value in the embodiments of this application.
[0063] Figure 14 This is a schematic diagram illustrating an excessive number of eigenvalue ranks in an embodiment of this application.
[0064] Figure 15 This is a schematic diagram of the fxy-mssa (rank ratio method) noise reduction method according to an embodiment of this application.
[0065] Figure 16 This is a schematic diagram of the denoising residual of the fxy-mssa (rank ratio method) method in an embodiment of this application.
[0066] Figure 17This is a schematic diagram of similar attribute bodies after processing by the fxy-mssa (rank ratio method) method in an embodiment of this application.
[0067] Figure 18 This is a denoising diagram of an embodiment of this application.
[0068] Figure 19 This is a schematic diagram of the denoised residual in an embodiment of this application.
[0069] Figure 20 Similar property diagrams for removing noise in embodiments of this application.
[0070] Figure 21 for Figure 11 A partial cross-sectional view within the box.
[0071] Figure 22 for Figure 15 A partial cross-sectional view within the box.
[0072] Figure 23 for Figure 16 A partial cross-sectional view within the box. Figure 24 This is a schematic diagram of a partial cross-section of a similar attribute body after processing by the fxy-mssa (rank ratio method) in an embodiment of this application. Figure 25 This is a partial cross-sectional schematic diagram of the noise reduction process in an embodiment of this application. Figure 26 This is a partial cross-sectional schematic diagram of the denoised residual in an embodiment of this application. Figure 27 This is a partial cross-sectional schematic diagram of a similar attribute body after processing according to an embodiment of this application. Detailed Implementation
[0073] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative effort are within the scope of protection of the present invention.
[0074] This embodiment is applied to pre-stack time-migrated 3D seismic data from region Z of the Junggar Basin. Please refer to [link / reference]. Figure 1 This embodiment provides a method for denoising low-rank seismic data using adaptive symplectic geometric decomposition, the method comprising:
[0075] S1. Use Fourier transform to convert 3D seismic data into the frequency domain;
[0076] In step S1, if the size of the 3D seismic data volume is (nt*nx*ny), where nt represents the number of sampling points, i.e., the number of samples in the time dimension; nx represents the number of points in the X-axis direction, i.e., the number of seismic data samples in the horizontal direction; and ny represents the number of points in the Y-axis direction, i.e., the number of seismic data samples in another horizontal direction, the data is cut into multiple seismic blocks of local size (ntwin*nxwin*nywin) using a windowing function, where ntwin is the number of samples in the window in the nt direction, nxwin is the number of samples in the window in the nx direction, and nywin is the number of samples in the window in the ny direction. Then, the time domain data D(y, x, y) is transformed into frequency domain data F(ω, x, y) using a Fourier transform, where t represents time, ω represents frequency, x is the horizontal direction in the plane, and y is another horizontal direction in the plane.
[0077] Specifically, the volume of work area A1 in region Z is 1500*300*240 mm in 3D seismic data. Figure 11 The data of the three-dimensional work area is displayed. The window size of the local division is 200*50*50, and the overlap in each direction is 50%. The time domain D(t, x, y) is transformed into the complex value F(w, x, y) in the frequency domain through Fourier transform.
[0078] S2. Slice the frequency domain seismic data into Hankel matrices according to frequency, and then form the Hankel matrices into Hamiltonian matrix form.
[0079] Step S2 includes:
[0080] The frequency component slice F(ω) of each frequency component in the frequency domain data is represented as:
[0081]
[0082] In the formula: nx represents the number of data points along the X-axis, i.e., the number of seismic data samples in the horizontal direction; ny represents the number of data points along the Y-axis, i.e., the number of seismic data samples in the other horizontal direction; ω represents the frequency.
[0083] The Hankel matrix formed by the frequency blocks in the i-th row is shown below:
[0084]
[0085] Where dx represents the embedding dimension in the X-axis direction; H x It is the Hankel matrix constructed along the X-axis; nx represents the number of data volumes along the X-axis, i.e., the number of samples of seismic data in the horizontal direction; ω represents the frequency;
[0086] Then construct a block Hankel matrix along the Y-axis:
[0087]
[0088] In the above formula, dy represents the embedding dimension in the Y-axis direction, and the data size of matrix H(ω) is I*J, where I = dx*dy, J = (nx - dx + 1) * (ny - dy + 1); ny represents the number in the Y-axis direction; ω represents the frequency.
[0089] Construct the Hamiltonian matrix M from the block Hankel matrix H:
[0090]
[0091] Where a is a symmetric matrix and satisfies A = H T H, H T A is the transpose of matrix H. T It is the transpose of A.
[0092] S3. Use symmetric QR decomposition to find the eigenvalues and eigenvectors of the Hamiltonian matrix;
[0093] Step S3 includes:
[0094] After constructing the Hamiltonian matrix, we use the square matrix N of the Hamiltonian matrix, i.e., N = M. 2 We can construct a symplectic orthogonal matrix G such that Where c is the upper Hessenberg matrix, and G represents a symplectic orthogonal matrix, which can be replaced by a given Household matrix Q. According to the symplectic similarity theorem of the symplectic QR algorithm, the matrix [Q 0; 0 Q] is also a symplectic orthogonal matrix, used to replace G. Then, the symplectic similarity theorem is used to ensure that the symplectic orthogonal matrix maintains the structure of the Hamiltonian system, and its eigenvalues follow... Where λ(A) is the eigenvalue of the symmetric matrix A, and λ(C) is the eigenvalue of the upper Hessenberg matrix c.
[0095] S4. Calculate the symplectic geometric entropy of the eigenvalues, and obtain the number of selection ranks from the adaptive symplectic geometric entropy;
[0096] Step S4 includes:
[0097] If the vector formed by the Hamiltonian matrices is s, where the length of s is L;
[0098] The first step is to normalize the eigenvectors:
[0099]
[0100] v j It is the j-th component in the standardized eigenvalue vector; It is the square of the j-th eigenvalue in the original eigenvalue vector s; Let be the sum of the squares of all eigenvalues in the standardized original eigenvalue vector s, k = 1, 2, ..., L, where L is the length of the eigenvalue vector;
[0101] The normalized symplectic geometric entropy e(N) is defined as:
[0102]
[0103] To achieve adaptive selection, the contribution of the i-th feature to the entropy must first be calculated, defined as ce. i , satisfy ce i =e(N)-e(N) -i ), where n -i This represents the i-th feature that was removed;
[0104] Features are ranked based on their relative contribution to entropy; the average of all ce values is defined as c, and their standard deviation is defined as d, where ce i -cd>0, features with high ce values are considered relevant.
[0105] S5. Reconstruct the block Hankel matrix using a low-rank approximation.
[0106] Step S5 includes:
[0107] The matrix is reconstructed using the eigenvalue order determined by the symplectic geometric entropy. First, the corresponding coefficient matrix is constructed. If the selected eigenrank is m, the low-rank block Hankel matrix is obtained by summing the first m reconstructed matrices, and then the data is reconstructed.
[0108]
[0109] in This refers to the reconstruction matrix of the block Hankel matrix H; Q represents the reconstructed orthogonal matrix after determining the eigenvalue order m through symplectic geometric entropy, and its columns are the eigenvectors of the original data matrix.
[0110] S6. Convert the processed frequency domain to the time domain using inverse Fourier transform.
[0111] Step S6 includes:
[0112] The matrix after low-rank reconstruction of the components at each frequency. The new matrix, formed by the diagonal averaging method, is then subjected to inverse Fourier transform to obtain the low-rank time-domain result.
[0113] In contrast. Figure 15 This is the denoising effect image after Fxy-mssa (rank ratio method) processing. Figure 16The corresponding residual plot is shown. This method takes 994 seconds. Since there is no pure signal used to calculate the SNR in this study, only local similarity metrics can be used to evaluate the denoising performance. The general standard is that if there is no signal leakage in the removed noise, the local similarity between the denoised data and the removed noise should be negligible. Figure 17 This is a similarity attribute map of the noise removed by Fxy-mssa (rank ratio method). Figure 18 The image shows the denoising effect after processing by the method of this invention. The program time for Fxy-sgmd (symplectic geometric entropy method) is 335s, which shows that the efficiency of the method of this invention is greatly improved. Figure 19 It is the corresponding residual plot. Figure 20 This is a similarity attribute diagram of the noise removed by the present invention. To more clearly demonstrate the effectiveness of the present invention, Figure 21 for Figure 11 The partial cross-sectional view within the box. Figure 22-27 for Figure 15-20 The partial cross-sectional view in the corresponding box shows that, when the structure is complex, the symplectic geometric entropy method can preserve the effective signal better than the rank ratio method.
[0114] This becomes clearer in the cross-sectional view within the red box area. The rank ratio method exhibits significant leakage of effective signals in similar attribute graphs, while the symplectic geometric entropy method more effectively preserves signals and prevents leakage. Therefore, this invention not only greatly improves computational efficiency but also adaptively preserves more effective signals in complex structures.
[0115] To verify the effectiveness and efficiency of the method proposed in this embodiment, tests were conducted on two types of synthetic seismic data: a noisy linear seismic model and a noisy curvilinear seismic model, as follows:
[0116] like Figure 2 The image shown is a linear, noisy synthetic earthquake used for testing, with a size of 500*40*40. Figure 3 This is the denoising result image after adaptive symplectic geometric decomposition (Fxy-sgmd) processing in this embodiment. Figure 4 This is a residual plot corresponding to the tested earthquakes, comparing Fxy-sgmd with other low-rank class methods in recent years. Figure 5 This is a graph comparing the efficiency of different methods for different volumes, while Table 1 shows the statistical results of each method:
[0117] Table 1 Comparison of processing effects and efficiency of different low-rank methods
[0118]
[0119] Table 1 compares the processing performance and efficiency of different low-rank methods. Considering the signal-to-noise ratio (SNR) of the processed results, it can be seen that Fxy-hosvd is only more efficient than the method presented in this paper when the volume is greater than 60*60, but the SNR after processing is low, and it is difficult to define the rank selection in each dimension. The Fxy-rmssa method is poor in both efficiency and SNR after processing. Although the Lanczos method is efficient, it requires knowing the number of ranks in all computational blocks in advance, which is impossible in reality. The method presented in this paper, Fxy-sgmd, can save about 50-70% of the computational cost compared to Fxy-mssa (rank ratio method). When the volume is less than 40*40, its efficiency is even better than that of Fxy-hosvd. However, when the window size increases, the time consumed by the algorithm increases sharply. Therefore, in actual operation, the window range should be set between 40-60. This avoids the Hankel matrix becoming too large, which would lead to a sharp drop in efficiency, and also avoids the window becoming too small, which would fail to consider the overall information.
[0120] This time, a noisy curve model was also used. Figure 6 The denoising effect was tested, and the size of the data was 500*100*100 when processing earthquake data. Figure 7 This is the denoising effect image after Fxy-mssa (rank ratio method) processing. Figure 8 It is the corresponding residual plot. Figure 9 This is a denoising effect image after processing by the method of the present invention. Figure 10 This is the corresponding residual plot. As can be seen from the local profile, the Fxy-mssa method struggles to preserve the complete signal of the effective layer, and the curve intersections are treated as noise, resulting in severe signal leakage. Its residuals contain significant effective signals. In contrast, the Fxy-sgmd method, after using symplectic geometric entropy for rank selection, is significantly more accurate than the Fxy-mssa method, and the effective signals at the curve intersections are well preserved.
[0121] In the method for selecting rank, the symplectic geometric entropy method is adopted. Compared with the commonly used singular value ratio method and energy ratio method for rank selection, the symplectic geometric entropy method has the characteristics of accurate noise localization and small error after denoising, respectively, and is effective against weak noise ( Figure 12 ), normal noise ( Figure 13 ), loud noise ( Figure 14In weak noise conditions, both the singular value ratio and energy ratio methods choose an optimal rank of 1, leading to the loss of signal features after denoising. The symplectic geometric entropy method avoids this phenomenon. Under normal noise conditions, it can be seen that the singular value ratio, energy ratio, and symplectic geometric entropy methods can accurately select feature signals. Under strong noise conditions, it can be seen that errors occur in the noise part due to the small denominator of the singular value ratio. Both the singular value ratio method and the energy ratio method incorrectly determine the information at the noise location, thus retaining a lot of noise. The symplectic geometric entropy method can effectively avoid this phenomenon.
[0122] In summary, this embodiment utilizes symplectic QR instead of SVD, which not only improves the efficiency of the method, but also proposes an adaptive rank selection method based on symplectic geometric entropy, which can adaptively retain more effective signals in complex structures.
[0123] This embodiment performs QR decomposition (symplectic QR method) within a symplectic geometric framework. Symplectic QR greatly improves computational efficiency by replacing SVD, and it also provides an optimal rank reduction method by combining symplectic geometric entropy. This invention not only greatly improves computational efficiency, but also adaptively retains more effective signals in complex structures.
[0124] This embodiment constructs a Hankel matrix and applies symplectic QR decomposition to extract the eigenvalues and eigenvectors of the Hankel matrix. It then uses symplectic geometric entropy to adaptively determine the number of ranks, thereby achieving low-rank approximate reconstruction, effectively denoising and preserving the effective signals in the seismic data.
[0125] The above description, in conjunction with specific preferred embodiments, provides a further detailed explanation of the present invention. It should not be construed that the specific implementation of the present invention is limited to these descriptions. For those skilled in the art, various simple deductions or substitutions can be made without departing from the concept of the present invention, and all such modifications and substitutions should be considered within the scope of protection of the present invention.
Claims
1. A method for denoising low-rank seismic data using adaptive symplectic geometric decomposition, characterized in that, The method includes: S1. Use Fourier transform to convert 3D seismic data into the frequency domain; S2. Slice the frequency domain seismic data into Hankel matrices according to frequency, and then form the Hankel matrices into Hamiltonian matrix form. S3. Use symmetric QR decomposition to find the eigenvalues and eigenvectors of the Hamiltonian matrix; S4. Calculate the symplectic geometric entropy of the eigenvalues, and obtain the number of selection ranks from the adaptive symplectic geometric entropy; S5. Reconstruct the block Hankel matrix using a low-rank approximation. S6. Convert the processed frequency domain to the time domain using inverse Fourier transform.
2. The method for denoising low-rank seismic data using adaptive symplectic geometric decomposition as described in claim 1, characterized in that, In step S1, if the size of the 3D seismic data volume is (nt*nx*ny), where nt represents the number of sampling points, i.e., the number of samples in the time dimension; nx represents the number of points in the X-axis direction, i.e., the number of seismic data samples in the horizontal direction; and ny represents the number of points in the Y-axis direction, i.e., the number of seismic data samples in another horizontal direction, the data is divided into multiple seismic blocks of local size (ntwin*nxwin*nywin) using a windowing function, where ntwin is the number of samples in the window in the nt direction, nxwin is the number of samples in the window in the nx direction, and nywin is the number of samples in the window in the ny direction. Then, the time domain data D(t,x,y) is transformed into frequency domain data F(ω,x,y) using a Fourier transform, where t represents time, ω represents frequency, x is the horizontal direction in the plane, and y is another horizontal direction in the plane.
3. The method for denoising low-rank seismic data using adaptive symplectic geometric decomposition as described in claim 1, characterized in that, Step S2 includes: The frequency component slice F(ω) of each frequency component in the frequency domain data is represented as: In the formula: nx represents the number of data points along the X-axis, i.e., the number of seismic data samples in the horizontal direction; ny represents the number of data points along the Y-axis, i.e., the number of seismic data samples in the other horizontal direction; ω represents the frequency. The Hankel matrix formed by the frequency blocks in the i-th row is shown below: Where dx represents the embedding dimension in the X-axis direction; H x It is the Hankel matrix constructed along the X-axis; nx represents the number of data volumes along the X-axis, i.e., the number of samples of seismic data in the horizontal direction; ω represents the frequency; Then construct a block Hankel matrix along the Y-axis: In the above formula, dy represents the embedding dimension in the Y-axis direction, and the data size of matrix H(ω) is I*J, where I = dx*dy, J = (nx - dx + 1) * (ny - dy + 1); ny represents the number in the Y-axis direction; ω represents the frequency. Construct the Hamiltonian matrix M from the block Hankel matrix H: Where A is a symmetric matrix and satisfies A = H T H, H T A is the transpose of matrix H. T It is the transpose of A.
4. The method for denoising low-rank seismic data using adaptive symplectic geometric decomposition as described in claim 1, characterized in that, Step S3 includes: After constructing the Hamiltonian matrix, we use the square matrix N of the Hamiltonian matrix, i.e., N = M. 2 We can construct a symplectic orthogonal matrix G such that Where C is the upper Hessenberg matrix, and G represents a symplectic orthogonal matrix, which can be replaced by a given Household matrix Q. According to the symplectic similarity theorem of the symplectic QR algorithm, the matrix [Q 0; 0 Q] is also a symplectic orthogonal matrix, used to replace G. Then, the symplectic similarity theorem is used to ensure that the symplectic orthogonal matrix maintains the structure of the Hamiltonian system, and its eigenvalues follow... Where λ(A) is the eigenvalue of the symmetric matrix A, and λ(C) is the eigenvalue of the upper Hessenberg matrix c.
5. The method for denoising low-rank seismic data using adaptive symplectic geometric decomposition as described in claim 1, characterized in that, Step S4 includes: If the vector formed by the Hamiltonian matrices is s, where the length of s is L; The first step is to normalize the eigenvectors: v j It is the j-th component in the standardized eigenvalue vector; It is the square of the j-th eigenvalue in the original eigenvalue vector s; Let be the sum of the squares of all eigenvalues in the standardized original eigenvalue vector s, k = 1, 2, ..., L, where L is the length of the eigenvalue vector; The normalized symplectic geometric entropy e(N) is defined as: To achieve adaptive selection, the contribution of the i-th feature to the entropy must first be calculated, defined as ce. i , satisfy ce i =e(N)-e(N) -i ), where N -i This represents the i-th feature that was removed; Features are ranked based on their relative contribution to entropy; the average of all ce values is defined as c, and their standard deviation is defined as d, where ce i -cd>0, features with high ce values are considered relevant.
6. The method for denoising low-rank seismic data using adaptive symplectic geometric decomposition as described in claim 1, characterized in that, Step S5 includes: The matrix is reconstructed using the eigenvalue order determined by the symplectic geometric entropy. First, the corresponding coefficient matrix is constructed. If the selected eigenrank is m, the low-rank block Hankel matrix is obtained by summing the first m reconstructed matrices, and then the data is reconstructed. in This refers to the reconstruction matrix of the block Hankel matrix H; Q represents the reconstructed orthogonal matrix after determining the eigenvalue order m through symplectic geometric entropy, and its columns are the eigenvectors of the original data matrix.
7. The method for denoising low-rank seismic data using adaptive symplectic geometric decomposition as described in claim 1, characterized in that, Step S6 includes: The matrix after low-rank reconstruction of the components at each frequency. The new matrix, formed by the diagonal averaging method, is then subjected to inverse Fourier transform to obtain the low-rank time-domain result.