Three-dimensional seismic data reconstruction method based on weighted fast iterative shrinkage threshold algorithm

Abstract

The invention discloses a three-dimensional seismic data reconstruction method based on a weighted fast iterative shrinkage threshold algorithm, and the method aims at solving a problem that a conventional three-dimensional seismic data reconstruction method neglects the potential correlation between adjacent slices, and comprises the steps: firstly dividing three-dimensional seismic data into a series of time slices; and carrying out multi-scale multi-direction two-dimensional curvelet transformation to obtain a curvelet domain coefficient, then introducing a fast iterative shrinkage threshold algorithm, and carrying out reconstruction by adopting an index threshold parameter formula. In the process, according to the characteristic that curvelet coefficient support sets of two adjacent time slices have a large intersection, a weighting operator is constructed in a curvelet domain by using a prior support set provided by a previous reconstructed time slice, and then the weighting operator is fused into a fast iterative shrinkage threshold algorithm. According to the method, the signal-to-noise ratio of the reconstructed signal is increased, the calculation speed is increased, the requirement for a computer memory is reduced, weak effective wave signals are protected, and the reflected wave event is more continuous and clearer.

Description

Technical Field

[0001] This invention relates to the field of seismic data reconstruction methods with irregular missing traces, and specifically discloses a three-dimensional seismic data reconstruction method based on a weighted fast iterative shrinkage threshold algorithm. Background Technology

[0002] With the rapid development of high-density 3D seismic exploration and wide-azimuth seismic exploration methods, there is a desire to acquire complete and regular seismic data. However, due to the influence of complex terrain conditions and acquisition environment, the acquired seismic data often contains a large number of irregular missing seismic traces, which seriously affects the requirements of high-resolution seismic data processing and interpretation, and reduces the accuracy of seismic exploration (Chen et al., 2016; Liu et al., 2018; Zhang et al., 2020). Due to the constraints of project exploration costs, it is generally not possible to re-acquire missing seismic traces in the field. In this case, it is necessary to use effective data reconstruction methods to recover the missing seismic traces during the budget processing stage, providing more geophysical information for subsequent processing, reducing the ambiguity of inversion, and meeting the requirements of high-precision seismic exploration (Gou et al., 2014; Huang et al., 2016; Li et al., 2022). Currently, many scholars both domestically and internationally have proposed various effective seismic data reconstruction methods, such as wave equation reconstruction methods (Fomel, 2003), predictive filter reconstruction methods (Liu et al., 2022), rank reduction reconstruction methods (Cavalcante and Porsani, 2022), sparse transform reconstruction methods (Zhang et al., 2017a), and artificial intelligence reconstruction methods (Chen et al., 2023). Among sparse transform-based data reconstruction methods, the most widely used is the Fourier transform-based seismic data reconstruction method (Cui et al., 2003; Gao et al., 2010). Generally speaking, this reconstruction method has simple parameter settings, does not require subsurface velocity structure information, has relatively fast computational efficiency, relatively high accuracy, and is relatively mature in commercial applications. Meanwhile, many scholars have also developed reconstruction methods based on other sparse transforms, such as the Ladon transform (Xue et al., 2014), wavelet transform (Liu et al., 2024), Shearlet transform (Yang et al., 2020), Dreamlet transform (Wang et al., 2015), Seislet transform (Liu et al., 2013), and curvelet transform (Mansour et al., 2013), which have also achieved good reconstruction results. To further improve accuracy, many scholars both domestically and internationally have extended seismic data reconstruction methods based on sparse transforms to three-dimensional or five-dimensional data spaces, further improving reconstruction accuracy by utilizing data information from different spatial domains. However, for three-dimensional seismic data reconstruction, a common approach is to divide the seismic data into individual slices for separate processing, increasing computation time and saving memory space (Zhang et al., 2013; Zhang et al., 2017b; Zhang et al., 2019).However, in general, there is a certain correlation between the effective wave energy of adjacent slice data (Friedlander et al., 2011). Traditional methods, which reconstruct slice by slice independently, cannot effectively utilize the prior information hidden between slices, thus reducing the reconstruction accuracy. Moreover, in the process of reconstructing 3D seismic data, traditional algorithms are relatively slow, especially when using curvelet transform with high redundancy as a sparse basis for data reconstruction. Although the reconstruction accuracy can be guaranteed, the computing speed is far from meeting the needs of processing massive seismic data (Zhang et al., 2018).

[0003] To this end, this invention employs curvelet transform to transform each time slice of 3D seismic data into the curvelet domain, and introduces a fast iterative threshold shrinkage algorithm to reconstruct each time slice sequentially based on the iterative threshold shrinkage algorithm. Summary of the Invention

[0004] The purpose of this invention is to enable high-precision reconstruction of irregular and missing seismic data in the field, significantly improve reconstruction accuracy and calculation speed, save memory space, and protect weak effective wave signals, thereby making the phase axis of reflected waves more continuous and clear. Therefore, a three-dimensional seismic data reconstruction method based on a weighted fast iterative shrinkage threshold algorithm is proposed.

[0005] To achieve the above objectives, the technical solution adopted by the present invention is as follows:

[0006] A method for reconstructing 3D seismic data based on a weighted fast iterative shrinkage threshold algorithm includes the following steps:

[0007] S1. Divide the three-dimensional seismic data into a series of time slices and perform multi-scale, multi-directional two-dimensional curvelet transform on it to obtain curvelet domain coefficients;

[0008] S2. Introduce a fast iterative shrinkage threshold algorithm and use an exponential threshold parameter formula for reconstruction;

[0009] S3. Based on the large overlap between the support sets of the curve wave coefficients of two adjacent time slices, a weighted operator is constructed in the curve wave domain using the prior support set provided by the previously reconstructed time slice. Then, the weighted operator is integrated into the fast iterative shrinkage threshold algorithm to realize a three-dimensional seismic data reconstruction method based on the weighted fast iterative shrinkage threshold algorithm.

[0010] Furthermore, the relationship between irregularly missing field data and complete seismic data can be represented by the following linear equation:

[0011] y = Mf

[0012] Where y∈R n This represents incomplete seismic data; f∈RN And N≥n, representing the complete seismic data that needs to be reconstructed; M∈R n×N Let be the sampling matrix. Assuming the sparsity coefficients of the seismic data f in the curvelet domain C are x, then the equation can be written as:

[0013] y = Ax and A = MC H

[0014] The superscript H represents the conjugate transpose.

[0015] Furthermore, the fast iterative threshold shrinkage algorithm is as follows:

[0016] Solving the above equation is an underdetermined ill-conditioned problem, and it is not easy to obtain an exact solution. Therefore, it can be transformed into an l1 minimum norm model, that is:

[0017] sty=MC H x

[0018] For ease of expression, the solution is written in the following form:

[0019]

[0020] The choice of the threshold parameter λ is crucial for solving this equation, and it needs to be continuously modified until the final solution meets certain accuracy requirements. Therefore, the Fast Iterative Thresholding Algorithm (FISTA) can be used to solve it, and its iterative formula is as follows:

[0021]

[0022]

[0023] y l This represents the l-th time slice. This represents the curvelet coefficient after the nth iteration of the l-th time slice, with an initial value of zero. The threshold is λ n The threshold function is denoted by t, where t represents the update factor.

[0024] Furthermore, the formula for the exponential threshold parameter is as follows:

[0025]

[0026] In the formula, Max is the maximum value of |Cy|, that is, the maximum value of the absolute value of the curve coefficient.

[0027] Furthermore, based on the large overlap between the support sets of the curvelet coefficients of two adjacent time slices, a weighted operator is constructed in the curvelet domain using the prior support set provided by the previously reconstructed time slice; the weighted operator is expressed as follows:

[0028]

[0029] Where J represents the maximum number of curve coefficients, w i This represents the weighted value, and its possible values ​​are as follows:

[0030]

[0031] The prior support set provided by the curvelet domain coefficients of the (l-1)th time slice is represented by the number of curvelet coefficients corresponding to a certain energy value, or directly by the percentage of the maximum number of curvelet coefficients. Since this weighting strategy is similar to a pre-filtering process, the weighting value γ outside the prior support set is generally between 0 and 1, i.e., 0 < γ ≤ 1. If it is too small, it may damage the effective signal; if it is too large, it will not play a filtering role. In particular, when γ = 1, this method will degenerate into a conventional fast iterative shrinking threshold algorithm.

[0032] Furthermore, the Weighted Fast Iterative Shrinking Threshold Algorithm (WFISTA) is expressed as follows:

[0033]

[0034]

[0035] The solution of the equation is continuously converged until the accuracy requirement is met. Then, the coefficients after N iterations are subjected to inverse curvelet transform to obtain the reconstructed complete seismic data. This not only improves the signal-to-noise ratio of the reconstructed signal, but also speeds up the calculation.

[0036] To address the issue that traditional 3D seismic data reconstruction methods neglect the potential correlation between adjacent slices, this paper first divides the 3D seismic data into a series of time slices and performs multi-scale, multi-directional 2D curvelet transform to obtain curvelet domain coefficients. Then, a fast iterative shrinkage thresholding algorithm is introduced, employing an exponential threshold parameter formula for reconstruction. During this process, considering the significant overlap between the curvelet coefficient support sets of two adjacent time slices, a weighted operator is constructed in the curvelet domain using the prior support set provided by the previously reconstructed time slice. This weighted operator is then integrated into the fast iterative shrinkage thresholding algorithm, thus forming a 3D seismic data reconstruction method based on a weighted fast iterative shrinkage thresholding algorithm.

[0037] The beneficial effects of this invention are as follows: 1. Compared with the wave equation reconstruction method, the technology of this invention can directly reconstruct irregular and missing seismic data without changing the amplitude and phase, with better fidelity, and does not require prior information on underground structures, thus reducing the computational workload.

[0038] 2. Compared with the Fourier transform reconstruction method, the technology of this invention can process nonlinear phase axis seismic data and seismic data with large missing intervals with high accuracy, and the reconstruction accuracy is high.

[0039] 3. Compared with 3D seismic data reconstruction based on the fast iterative shrinkage threshold algorithm, the technology of this invention not only improves the signal-to-noise ratio of the reconstructed signal and speeds up the calculation, but also reduces the requirements for computer memory, protects weak effective wave signals, and makes the phase axis of the reflected wave more continuous and clear.

[0040] In summary, this invention overcomes the shortcomings of conventional reconstruction methods, such as their inability to quickly and accurately reconstruct irregularly missing seismic data, low reconstruction accuracy, and limited application. Furthermore, this invention can reconstruct not only irregularly missing seismic data and spurious-frequency seismic data with large gaps in the missing traces, but also regularly missing seismic data. It significantly improves the computation speed and signal-to-noise ratio of the reconstructed signal, saves computational space, and protects weak effective wave signals, thereby making the reflected wave phase axis more continuous and clear. This has significant practical value for guiding seismic data acquisition and missing trace reconstruction in complex areas. Attached Figure Description

[0041] To more clearly illustrate the technical solutions in the embodiments of this application, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0042] Figure 1 This is a flowchart of the reconstruction process for irregularly missing seismic data in an embodiment of the present invention;

[0043] Figure 2 This is a theoretical seismic data diagram of the present invention;

[0044] Figure 3 These are reconstruction results of the present invention under different algorithms;

[0045] Figure 4 These are reconstruction error maps using different iterative algorithms of this invention;

[0046] Figure 5 This is a graph showing the number of iterations and the signal-to-noise ratio after reconstruction in this invention.

[0047] Figure 6 This is a signal-to-noise ratio curve after reconstruction with different weighting values ​​according to the present invention;

[0048] Figure 7 This is a signal-to-noise ratio curve of different time slices according to the present invention;

[0049] Figure 8 This is a noisy seismic data map in an embodiment of the present invention;

[0050] Figure 9 This is a diagram showing the reconstruction result of noisy seismic data in an embodiment of the present invention. Detailed Implementation

[0051] The technical solution of the present invention will now be clearly and completely described with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0052] Reference Figure 1 A three-dimensional seismic data reconstruction method based on a weighted fast iterative shrinkage threshold algorithm is proposed. The main steps of implementing this method include: establishing a seismic data reconstruction model, using a conventional fast iterative shrinkage threshold algorithm, formulating an exponential threshold parameter, establishing a weighted operator, and implementing the weighted fast iterative shrinkage threshold algorithm.

[0053] The specific steps are as follows:

[0054] Step 1: Data reconstruction model establishment; the relationship between irregular missing data in the field and complete seismic data can be represented by the following linear equation:

[0055] y = Mf

[0056] Where y∈R n This represents incomplete seismic data; f∈R N And N≥n, representing the complete seismic data that needs to be reconstructed; M∈R n×N Let be the sampling matrix. Assuming the sparsity coefficients of the seismic data f in the curvelet domain C are x, then the equation can be written as:

[0057] y = Ax and A = MC H

[0058] The superscript H represents the conjugate transpose.

[0059] Step 2: Conventional fast iterative threshold shrinkage algorithm; Solving the above equation is an underdetermined ill-conditioned problem, which is not easy to obtain an exact solution. Therefore, it can be transformed into an l1 minimum norm model, that is:

[0060] sty=MC H x

[0061] For ease of expression and solution, it can be written in the following form:

[0062]

[0063] The choice of the threshold parameter λ is crucial for solving this equation, and it needs to be continuously modified until the final solution meets certain accuracy requirements. Therefore, the Fast Iterative Thresholding Algorithm (FISTA) can be used to solve it, and its iterative formula is as follows:

[0064]

[0065] y l This represents the l-th time slice. This represents the curvelet coefficient after the nth iteration of the l-th time slice, with an initial value of zero. The threshold is λ n The threshold function is denoted by t, where t represents the update factor.

[0066] Step 3: Formula for the exponential threshold parameter. In solving the above nonlinear problem, the threshold parameter λ... n Satisfy, ||Cy|| ∞ =λ1>λ2>…>ε, where ε is a small value close to zero. Commonly used threshold parameters include linear threshold parameters and exponential threshold parameters. To improve efficiency, this invention introduces a threshold parameter with exponential decay, resulting in a larger gradual decrease in the threshold parameter. This ensures both solution accuracy and faster convergence speed, meeting the requirements of processing massive seismic data. The formula for this threshold parameter is:

[0067]

[0068] In the formula, Max is the maximum value of |Cy|, that is, the maximum value of the absolute value of the curve coefficient.

[0069] Step 4: Establish the weighted operator formulation. Based on the large overlap between the support sets of curvelet coefficients of two adjacent time slices, a weighted operator is constructed in the curvelet domain using the prior support set provided by the previously reconstructed time slice. The weighted operator formulation is as follows:

[0070]

[0071] Where J represents the maximum number of curve coefficients, w i This represents the weighted value, and its possible values ​​are as follows:

[0072]

[0073] here, The prior support set provided by the curvelet domain coefficients of the (l-1)th time slice represents the number of curvelet coefficients corresponding to a certain energy value, or simply the percentage of the largest curvelet coefficient. Since this weighting strategy is similar to a pre-filtering process, the weighting value γ outside the prior support set is generally between 0 and 1, i.e., 0 < γ ≤ 1. If the value is too small, it may damage the effective signal; if it is too large, it will not perform the filtering process. Especially when γ = 1, this invention will degenerate into a conventional fast iterative thresholding algorithm.

[0074] Step 5: Weighted Fast Iterative Shrinking Threshold Algorithm. The weighting operator W is incorporated into the fast iterative shrinking threshold algorithm, resulting in the following weighted fast iterative shrinking threshold algorithm formula (WFISTA):

[0075]

[0076]

[0077] This formula can converge continuously until the accuracy requirement is met. Then, by performing inverse curvelet transform on the coefficients after N iterations, the reconstructed complete seismic data can be obtained. This not only improves the signal-to-noise ratio of the reconstructed signal, but also speeds up the calculation.

[0078] The specific steps to implement this method are as follows:

[0079] To measure the effectiveness of data reconstruction methods, the signal-to-noise ratio (SNR) formula is defined as SNR = 20log 10 The signal-to-noise ratio (SNR) is compared using ||x0||2 / ||x-x0||2, with the unit being dB. Here, x0 represents the original model data, and x represents the reconstruction result. The higher the SNR, the closer the reconstruction result is to the model data, and the more ideal the processing effect.

[0080] First, we arrange the seismic data collected along the two-dimensional survey line into a three-dimensional data volume of 256×256×256 along the receiver point, shot point and time coordinate. Figure 2 'a' represents theoretical seismic data, with trace and shot distances both at 12m and a sampling rate of 4ms. To simulate the irregular missing data in field seismic data, 50% of the time slices are randomly sampled in two dimensions, such as... Figure 2As shown in b. To save space and for a more intuitive display, subsequent reconstruction processing only shows the 110th time slice of this data. Then, four methods were used for reconstruction: Iterative Shrink Thresholding Algorithm (ISTA), Weighted Iterative Shrink Thresholding Algorithm (WISTA), Fast Iterative Shrink Thresholding Algorithm (FISTA), and Enhanced Fast Iterative Shrink Thresholding Algorithm (WFISTA). All methods involved 15 iterations, with 5 scales for the curvelet transform and 16 angles at the coarsest scale. The weighting value γ was set to 0.4, and the prior support set consisted of the top 2% of the largest curvelet coefficients. The reconstruction results are shown below. Figure 3-4 As shown, the reconstructed signal-to-noise ratios (SNRs) are 16.7 dB, 20.9 dB, 23.4 dB, and 24.6 dB, respectively. This demonstrates that all methods can recover missing seismic traces. However, overall, the WISTA and WFISTA methods show higher reconstruction accuracy and smaller errors than their corresponding ISTA and FISTA methods. This is because they employ weighted operators, which effectively utilize prior information from the previous time slice to provide support set information for the next adjacent slice to be reconstructed, delineating the main effective wave energy distribution range and suppressing low-energy random noise outside the range, thus enabling high-precision reconstruction of even weak effective signals. Furthermore, it is noted that the FISTA method achieves a higher SNR than the WISTA method after only 15 iterations, indicating that the fast iterative shrinkage threshold algorithm can achieve a higher SNR with fewer iterations. Moreover, by introducing weighted operators, the computational speed can be further improved, resulting in even better reconstruction performance.

[0081] To further compare the reconstruction performance of the method of this invention with that of traditional methods at different iteration numbers, the weighting value γ was set to 0.4, and the prior support set was set to the number of the top 2% of the largest curvelet coefficients. Then, the maximum number of iterations was varied from 2 to 60, and the theoretically missing data was reconstructed sequentially. The relationship between the maximum number of iterations and the signal-to-noise ratio after reconstruction was calculated, as shown in the graph. Figure 5As shown, before 20 iterations, the WFISTA method had the highest signal-to-noise ratio after reconstruction, followed by the FISTA and WISTA methods, and finally the ISTA method. However, after 20 iterations, the WISTA method gradually achieved higher reconstruction accuracy than the FISTA method. Furthermore, with the increase in the number of iterations, the final convergence results of the WFISTA and WISTA methods were similar, as were the final convergence results of the FISTA and ISTA methods, but overall, the FISTA and ISTA methods were lower than the WFISTA and WISTA methods. This demonstrates that utilizing prior information from adjacent time slices can effectively improve reconstruction accuracy. It can also be seen that to obtain a reconstruction result with a signal-to-noise ratio of 24dB, the ISTA method requires 40 iterations, the FISTA method requires 35 iterations, and the WISTA method requires 20 iterations, while the WFISTA method proposed in this invention only requires about 11 iterations, saving nearly four times the computation time compared to the ISTA method. This has a significant advantage in processing massive amounts of seismic data. Although the weighted operator operation requires some time, it is negligible compared to the overall computation time. In general, the signal-to-noise ratio (SNR) after reconstruction using all four iterative algorithms increases with the number of iterations. However, too many iterations waste computation time, requiring a trade-off. For the WFISTA method, the increase in SNR with increasing iteration count becomes relatively small after 15 iterations. Therefore, this invention uses 15 iterations for both theoretical simulations and practical data applications.

[0082] In fact, the choice of weighting value γ is particularly important during the reconstruction process; choosing too small or too large a value will not yield good reconstruction results. Therefore, we kept the weighting value varying from 0 to 1, and then calculated the signal-to-noise ratio curves after reconstruction using the WFISTA and WISTA methods. The results are as follows. Figure 6 As shown in the figure, the signal-to-noise ratio (SNR) is highest when γ is 0.4. The SNR decreases when the weighting value is greater than or less than 0.4, especially when the weighting value is 1, in which case the proposed method degenerates into a conventional iterative shrinkage threshold algorithm. Therefore, in the theoretical noise-free data simulation experiment, the weighting value was chosen to be 0.4.

[0083] Since the theoretical data has 256 time slices, in order to further compare the advantages of the method of this invention, Figure 7The image shows the signal-to-noise ratio (SNR) curves after reconstruction of the entire time slice using different iterative methods. Overall, in each time slice reconstruction, WFISTA and WISTA show higher SNRs than their corresponding FISTA and ISTA reconstructions. This indicates a significant overlap between the curvelet coefficient support sets of adjacent time slices, suggesting that the previous time slice can provide prior support set information for the next adjacent slice to be reconstructed. Previous methods, however, neglected the potential correlation between adjacent slices, thus affecting the reconstruction results. This invention utilizes prior information between different slices and adds constraints, resulting in better reconstruction performance. It should be noted that because this invention only uses 15 iterations, FISTA achieves better reconstruction results than WISTA, as... Figure 4 As shown, if the number of iterations is increased by 20, the reconstruction effect of the WISTA method will be better than that of the FISTA method. Meanwhile, as time increases, the signal-to-noise ratio curve oscillates and alternates, while the overall signal-to-noise ratio gradually decreases, indicating that the effective wave energy at depth decreases, leading to poor recovery of some weak energy signals.

[0084] Since field seismic data is generally quite noisy, it is necessary to test the reconstruction effect of the method of this invention under such conditions. Therefore, a certain proportion of Gaussian random noise was added to the original data to simulate noisy field seismic records. The noise-added result is as follows: Figure 8 As shown in a, Figure 8 b represents the 50% two-dimensional random undersampling result. Then, four iterative methods—ISTA, FISTA, WISTA, and WFISTA—were used to reconstruct the noisy missing channel data, and the reconstruction effect was tested under conditions of relatively high noise levels. Figure 9 The reconstruction results and errors of the method of the present invention show that all the information of the noisy and missing seismic traces has been well recovered with small errors. The phase axes of the reconstructed seismic waves are very continuous. On the other hand, it also shows that the method of the present invention has strong noise resistance in the reconstruction of noisy seismic data and can meet the needs of processing actual data.

[0085] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit them. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features therein. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for reconstructing three-dimensional seismic data based on a weighted fast iterative shrinkage threshold algorithm, characterized in that, Includes the following steps: S1. Divide the three-dimensional seismic data into a series of time slices and perform multi-scale, multi-directional two-dimensional curvelet transform on it to obtain curvelet domain coefficients; S2. Introduce a fast iterative shrinkage threshold algorithm and use an exponential threshold parameter formula for reconstruction; S3. Based on the large overlap between the support sets of the curve wave coefficients of two adjacent time slices, a weighted operator is constructed in the curve wave domain using the prior support set provided by the previously reconstructed time slice. Then, the weighted operator is integrated into the fast iterative shrinkage threshold algorithm to realize a three-dimensional seismic data reconstruction method based on the weighted fast iterative shrinkage threshold algorithm.

2. The three-dimensional seismic data reconstruction method based on a weighted fast iterative shrinkage threshold algorithm according to claim 1, characterized in that, The relationship between irregular, missing field data and complete seismic data can be represented by the following linear equation: y = Mf Where y∈R n This represents incomplete seismic data; f∈R N And N≥n, representing the complete seismic data that needs to be reconstructed; M∈R n×N Let be the sampling matrix. Assuming the sparsity coefficients of the seismic data f in the curvelet domain C are x, then the equation can be written as: y = Ax and A = MC H The superscript H represents the conjugate transpose.

3. The three-dimensional seismic data reconstruction method based on a weighted fast iterative shrinkage threshold algorithm according to claim 2, characterized in that, The fast iterative threshold shrinkage algorithm is as follows: Solving the above equation is an underdetermined ill-conditioned problem, and it is not easy to obtain an exact solution. Therefore, it can be transformed into an l1 minimum norm model, that is: s.t.y=MC H x For ease of expression, the solution is written in the following form: The choice of the threshold parameter λ is crucial for solving this equation, and it needs to be continuously modified until the final solution meets certain accuracy requirements. Therefore, the Fast Iterative Thresholding Algorithm (FISTA) can be used to solve it, and its iterative formula is as follows: t0＝1 y l This represents the l-th time slice. This represents the curvelet coefficient after the nth iteration of the l-th time slice, with an initial value of zero. The threshold is λ n The threshold function is given by t, where t represents the update factor.

4. A three-dimensional seismic data reconstruction method based on a weighted fast iterative shrinkage threshold algorithm according to claim 2, characterized in that, The formula for the exponential threshold parameter is as follows: In the formula, Max is the maximum value of |Cy|, that is, the maximum value of the absolute value of the curve coefficient.

5. A three-dimensional seismic data reconstruction method based on a weighted fast iterative shrinkage threshold algorithm according to claim 2, characterized in that, Based on the large overlap between the support sets of the curve wave coefficients of two adjacent time slices, a weighted operator is constructed in the curve wave domain using the prior support set provided by the previously reconstructed time slice. The weighted operator is expressed as follows: Where J represents the maximum number of curve coefficients, w i This represents the weighted value, and its possible values ​​are as follows: The prior support set provided by the curvelet domain coefficients of the (l-1)th time slice is represented by the number of curvelet coefficients corresponding to a certain energy value, or directly by the percentage of the maximum number of curvelet coefficients. Since this weighting strategy is similar to a pre-filtering process, the weighting value γ outside the prior support set is generally between 0 and 1, i.e., 0 < γ ≤ 1. If it is too small, it may damage the effective signal; if it is too large, it will not play a filtering role. In particular, when γ = 1, this method will degenerate into a conventional fast iterative shrinking threshold algorithm.

6. A three-dimensional seismic data reconstruction method based on a weighted fast iterative shrinkage threshold algorithm according to claim 5, characterized in that, The weighted fast iterative shrinkage threshold algorithm (WFISTA) is expressed as follows: t0＝1 The solution of the equation is continuously converged until the accuracy requirement is met. Then, the coefficients after N iterations are subjected to inverse curvelet transform to obtain the reconstructed complete seismic data, which improves the signal-to-noise ratio of the reconstructed signal and speeds up the calculation.

Images