A seismic data reconstruction method and system based on a TV regularization constraint time-frequency joint self-attention diffusion model

Abstract

The application belongs to the field of seismic data processing, and discloses a seismic data reconstruction method and system based on a time-frequency joint self-attention diffusion model with TV regularization constraint. The method comprises the following steps: obtaining two-dimensional ideal seismic data, and constructing a seismic data set for diffusion model training based on the two-dimensional ideal seismic data; based on the seismic data set, a multi-channel feature map is constructed, and a diffusion model based on a self-attention mechanism is trained using the multi-channel feature map to obtain a seismic data reconstruction model; based on the multi-directional wavelet coefficients output by the seismic data reconstruction model, a preliminary prediction result is obtained through inverse wavelet transform, and then the preliminary prediction result is fused with the original observation data according to a mask matrix to obtain complete seismic data after final reconstruction. The application realizes high-precision recovery of missing seismic data, improves the reconstruction accuracy, and ensures the continuity and structural integrity of the seismic data.

Description

Technical Field

[0001] This invention belongs to the field of earthquake data processing technology, and particularly relates to an earthquake data reconstruction method and system based on a time-frequency joint self-attention diffusion model with TV regularization constraints. Background Technology

[0002] In actual seismic data acquisition, due to the complex and varied geological conditions, rugged terrain, diverse subsurface media properties, and limitations in instrument deployment and acquisition costs, seismic data often becomes missing or incomplete. This missing data not only reduces overall data coverage but also significantly lowers the signal-to-noise ratio (SNR) of seismic records, affecting the clarity of seismic events. Furthermore, missing data disrupts the lateral continuity of seismic reflection signals, making it difficult to accurately identify subsurface stratigraphic structures and increasing the difficulty of seismic data interpretation and processing. The existence of missing data adversely affects subsequent seismic data processing and analysis, such as inversion, imaging, and reservoir prediction, severely limiting the accuracy and reliability of seismic exploration. Therefore, effectively recovering missing data and improving data quality has become one of the key challenges in seismic data processing.

[0003] Seismic data reconstruction is a crucial step in seismic data processing. Its primary goal is to recover missing seismic traces caused by limitations in acquisition conditions, instrument deployment, or the complexity of the field geological environment. The reconstruction problem can essentially be formulated as an inverse problem: reconstructing a complete subsurface seismic wavefield based on limited observed seismic data. Because seismic signals exhibit high spatial and temporal correlation, recovering missing data requires not only maintaining the consistency of observed data but also reasonably inferring structural information about the missing areas.

[0004] Traditional seismic data reconstruction methods are typically based on certain physical models or prior assumptions about the seismic data. These methods mainly include wavefield extension operator-based reconstruction, transform domain-based reconstruction, matrix rank reduction-based reconstruction, and predictive filtering-based reconstruction. Ronen proposed a de-aliasing reconstruction method based on dip time difference operators, primarily targeting frequency aliasing caused by missing or damaged portions of seismic data, effectively reducing the impact of aliasing. Wavefield extension operator-based seismic data reconstruction methods often face high computational costs, with computational resource requirements frequently becoming a bottleneck, resulting in limited research in recent years. When facing reconstruction problems with missing rules, many sparse transforms are susceptible to frequency aliasing interference. Since aliasing generally occurs in the high-frequency components of the data, while low-frequency components typically have less aliasing, some studies utilize the less aliased low-frequency information to aid in the reconstruction of high-frequency information, achieving de-aliasing reconstruction. Naghizadeh and Sacchi proposed a curvelet transform-based de-aliasing reconstruction method, assuming that seismic data has similar local slopes at coarse and fine scales. They filter aliasing coefficients using a mask function and invert curvelet coefficients for de-aliasing reconstruction. Zhang et al. combined the classic rank-reduced seismic data reconstruction algorithm MSSA with the POCS iterative framework based on Fourier transform to simultaneously reconstruct and denoise 3D seismic data, obtaining high-precision reconstruction results. Spitz proposed an fx-domain reconstruction method that uses the low-frequency components in the seismic data to estimate a prediction error filter, and then applies this prediction error filter to the high-frequency components of the data, thereby de-aliasing the high-frequency parts and significantly improving the reconstruction accuracy.

[0005] Traditional methods constrain missing data by establishing mathematical models; however, these methods often face limitations in complex data environments. First, they are based on fixed assumptions and lack sufficient adaptability, especially when dealing with high-frequency signals or complex missing data, where their performance falls short of expectations. Second, traditional reconstruction methods typically rely on meticulous parameter selection, with reconstruction results heavily dependent on manual parameter tuning. This selection process is highly subjective and prone to uncertainty in the reconstruction results. With the continuous advancement of seismic exploration technology, especially the introduction of "two-wide-one-high" (wide azimuth, wide bandwidth, and high density) technology, oil and gas geophysical exploration has gradually entered a new stage of high-precision, large-scale data processing. Accompanying this transformation, the dimensionality and scale of seismic data are constantly increasing, with data volume rapidly jumping from GB to TB or even PB levels. Therefore, traditional reconstruction methods based on assumptions and parameter adjustment can no longer meet the demands of current seismic data reconstruction tasks. To address this challenge, developing more intelligent and efficient reconstruction methods, especially technologies that reduce reliance on manual parameter selection, has become an urgent priority to adapt to the intelligent development needs of the petroleum geophysical exploration field.

[0006] In recent years, with the development of deep learning and generative models, data-driven methods can learn the intrinsic patterns of seismic wavefields, providing more flexible and efficient solutions to reconstruction problems, thus showing significant advantages in improving data quality and maintaining the continuity of seismic signals. Li et al. introduced a coordinate attention module into the Unet network for the reconstruction of two-dimensional continuous missing traces, and used a hybrid loss function combining structural similarity and L1 norm as the loss function to further improve the reconstruction performance of the designed network; Wang et al. introduced a residual learning network for seismic data reconstruction, performing anti-spoofing reconstruction on seismic data with missing rules, enhancing the lateral continuity of seismic data; Oliveira et al. proposed an innovative method of applying conditional generative adversarial networks to seismic data reconstruction, solving the problem of missing data caused by seismic data acquisition issues.

[0007] Most deep learning-based seismic data reconstruction methods employ Convolutional Neural Networks (CNNs) to learn the mapping between missing and complete data. Due to their localized perception, CNNs can effectively learn the mapping from local regions in the input data to the central pixel. Specifically, in seismic data reconstruction, CNNs are a point-scale prediction task; the network predicts for each data point, maintaining a consistent input and output dimension. However, because CNN layers are almost entirely local operators, at each pixel in the output image, the network only connects to a specific receptive field (the region centered on that data point) in the input. Therefore, CNNs have limitations in capturing global dependencies and long-range relationships in the data, especially when dealing with spatially widely dependent seismic data, making it difficult to capture the complex relationships between distant seismic traces. Furthermore, due to the complexity and non-stationarity of seismic data, the reconstruction problem is often ill-conditioned and underdetermined, meaning there are multiple possible solutions, requiring the introduction of prior information or regularization constraints to stabilize the reconstruction process. Summary of the Invention

[0008] The purpose of this invention is to provide a seismic data reconstruction method and system based on a time-frequency joint self-attention diffusion model with TV regularization constraints. Low-frequency global structure and high-frequency detail information are extracted through multi-directional wavelet transform, and sub-bands are processed in parallel, thereby accelerating the training and convergence of the diffusion model. Simultaneously, a self-attention mechanism is introduced to achieve non-local feature interaction, overcoming the problem of limited receptive field in traditional convolutional networks. Combined with total variational regularization constraints (TV), high-precision recovery of missing seismic data is achieved, improving reconstruction accuracy while ensuring the continuity and structural integrity of the seismic data.

[0009] This invention provides the following technical solution:

[0010] A seismic data reconstruction method based on a time-frequency joint self-attention diffusion model with TV regularization constraints includes the following steps:

[0011] Two-dimensional ideal seismic data is acquired, and a seismic dataset for training a diffusion model is constructed based on the two-dimensional ideal seismic data;

[0012] Based on the earthquake dataset, a multi-channel feature map is constructed, and the multi-channel feature map is used to train a diffusion model based on a self-attention mechanism to obtain the earthquake data reconstruction model.

[0013] Based on the multi-directional wavelet coefficients output by the earthquake data reconstruction model, preliminary prediction results are obtained through inverse wavelet transform. Subsequently, the preliminary prediction results are fused with the original observation data according to the mask matrix to obtain the final reconstructed complete earthquake data.

[0014] Preferably, the earthquake dataset includes paired data, namely ideal earthquake data and training earthquake data;

[0015] The method for constructing an earthquake dataset for training a diffusion model based on the aforementioned two-dimensional ideal seismic data includes:

[0016] The ideal seismic data is then normalized.

[0017] The training seismic data is generated by sampling the normalized ideal seismic data and generating a corresponding mask matrix, and then performing a Hadamard product operation on the ideal seismic data and the mask matrix.

[0018] After preprocessing the ideal seismic data and the training seismic data, the final seismic dataset is formed.

[0019] Preferably, the seismic data reconstruction model consists of three components: multi-directional feature map hybrid modeling, diffusion model reconstruction network construction based on self-attention mechanism, and time-frequency joint loss function setting;

[0020] The multi-directional feature map hybrid modeling is used to construct multi-channel feature maps based on the earthquake dataset;

[0021] The self-attention mechanism-based diffusion model reconstruction network is used to introduce the self-attention mechanism into the diffusion model and is trained using the multi-channel feature map to establish a seismic data reconstruction model.

[0022] The time-frequency joint loss function is used to set the reconstruction loss function during the training process of the diffusion model.

[0023] Preferably, the method for constructing a multi-channel feature map based on the earthquake dataset includes:

[0024] Using the ideal seismic data as the original complete seismic data and the training seismic data as the missing seismic data, two-dimensional discrete wavelet transforms are performed on both to obtain multi-directional sub-bands, including a low-frequency sub-band. , and high-frequency subband , ;

[0025] The multi-directional sub-bands are spliced ​​together in multiple channels to form the multi-channel feature map.

[0026] Preferably, the method of introducing a self-attention mechanism into the diffusion model and training it using the multi-channel feature map includes:

[0027] Flatten the multi-channel feature map The sequence is used to input the feature map into the self-attention mechanism for processing. Here, B represents the batch size, i.e. the number of samples input into the model at one time, H represents the height of the feature map, W represents the width of the feature map, and C represents the number of channels.

[0028] The self-attention mechanism achieves non-local feature interaction at each position in the sequence by weighted summation of the value matrix using the attention matrix.

[0029] Preferably, during the forward diffusion process of the diffusion model, Gaussian noise is added to the multi-channel feature map at each time step t.

[0030] In the reverse diffusion process of the diffusion model, from the noise characteristics Initially, seismic data is gradually restored, noise is predicted using a trained diffusion model, and reverse recovery is guided by wavelet multi-directional features.

[0031] Preferably, the reconstruction loss function is:

[0032] ,

[0033] in, Losses due to reconstruction of missing areas; This represents the loss of consistency across the observation area. For total variational regularization loss, These are the weighting coefficients.

[0034] The present invention also provides a seismic data reconstruction system based on a time-frequency joint self-attention diffusion model with TV regularization constraints. The system applies the aforementioned method and includes a dataset construction module, a reconstruction model module, and a data reconstruction module.

[0035] The dataset construction module is used to acquire two-dimensional ideal seismic data and construct an earthquake dataset for training the diffusion model based on the two-dimensional ideal seismic data;

[0036] The reconstruction model module is used to construct a multi-channel feature map based on the earthquake dataset, and use the multi-channel feature map to train a diffusion model based on a self-attention mechanism to obtain the earthquake data reconstruction model.

[0037] The data reconstruction module is used to obtain preliminary prediction results based on the multi-directional wavelet coefficients output by the seismic data reconstruction model through inverse wavelet transform. Subsequently, the preliminary prediction results are fused with the original observation data according to the mask matrix to obtain the final reconstructed complete seismic data.

[0038] Preferably, the reconstruction model module includes a feature map unit, a model training unit, and a loss function unit;

[0039] The feature map unit is used to construct a multi-channel feature map based on the earthquake dataset;

[0040] The model training unit is used to introduce the self-attention mechanism into the diffusion model and to train it using the multi-channel feature map to establish a seismic data reconstruction model.

[0041] The loss function unit is used to set the reconstruction loss function during the training process of the diffusion model.

[0042] The beneficial effects of this invention are as follows:

[0043] This invention provides a seismic data reconstruction method and system based on a time-frequency joint self-attention diffusion model with TV regularization constraints, which has the following advantages compared with traditional methods:

[0044] 1. A stepwise denoising mechanism based on the diffusion model is proposed, which utilizes wavelet transform to extract multi-directional features and parallel processing to accelerate the training and convergence of the diffusion model and improve the reconstruction computation efficiency.

[0045] 2. By introducing a self-attention mechanism to achieve non-local feature interaction, structural guidance is provided during the gradual denoising / reconstruction process of diffusion, overcoming the problem of limited receptive field in traditional convolutional networks, thereby more accurately recovering missing seismic data;

[0046] 3. In the stepwise denoising reconstruction process of the diffusion model, a total variational regularization constraint is introduced and a gradient constraint is applied to enhance the smoothness and continuity of the reconstruction results and ensure the structural integrity of the seismic data.

[0047] 4. By integrating time and frequency domain information, high-precision seismic data reconstruction is achieved. Compared with traditional convolution methods, the reconstruction accuracy and signal-to-noise ratio are significantly improved, verifying the feasibility and superiority of the method in seismic data reconstruction and restoration.

[0048] 5. Robustness to multiple missing patterns: It exhibits stability in three scenarios: random missing, regular missing, and continuous missing.

[0049] This invention innovatively combines multiple advanced technologies to effectively restore the overall structure and local details of seismic data while ensuring reconstruction accuracy. It significantly improves the lateral continuity and reliability of the data, while also taking into account computational efficiency, providing a feasible and efficient solution for large-scale seismic data processing. Attached Figure Description

[0050] To more clearly illustrate the technical solution of the present invention, the drawings used in the embodiments are briefly described below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0051] Figure 1 This is a flowchart of a seismic data reconstruction method based on a time-frequency joint self-attention diffusion model with TV regularization constraints, according to an embodiment of the present invention.

[0052] Figure 2 This is a flowchart of the training method according to an embodiment of the present invention;

[0053] Figure 3 This is a schematic diagram of the network structure of the seismic data reconstruction model according to an embodiment of the present invention;

[0054] Figure 4 This is a schematic diagram of self-attention feature transformation and loss calculation in seismic data reconstruction according to an embodiment of the present invention;

[0055] Figure 5 The image shows a sample of test data with 30% random missing values, excluding the self-attention mechanism module and the TV regularization term, as well as the reconstruction result image, according to an embodiment of the present invention.

[0056] Figure 6 The image shown is a sample of test data with 30% random missing data, along with the reconstruction result image, selected from any embodiment of the present invention.

[0057] Figure 7 The above is a sample image of test data and reconstruction result image of an embodiment of the present invention with 50% missing rule of any one of the rules, which does not contain the self-attention mechanism module and the TV regularization term.

[0058] Figure 8The above are sample images of test data and reconstruction results for any one rule in an embodiment of the present invention that are 50% missing.

[0059] Figure 9 The above is a sample image of test data and reconstruction results for any one of the following embodiments of the present invention: 16 consecutive missing channels without the self-attention mechanism module and TV regularization term.

[0060] Figure 10 The image shows a sample test data image and a reconstruction result image of any one of the 16 consecutive missing channels in Example 1.

[0061] Figure 11 This is a comparison chart showing the changes in SNR of different ablation configurations with training rounds under the condition of 30% random missing values ​​in the validation set in this embodiment of the invention.

[0062] Figure 12 This is a comparison chart showing the changes in SNR of different ablation configurations with training rounds under the condition that 50% of the rules are missing in the validation set in this embodiment of the invention.

[0063] Figure 13 This is a comparison chart showing the SNR of different ablation configurations with training rounds under the condition of 16 consecutive missing channels in the validation set of this invention.

[0064] Figure 14 This is a comparison diagram of the inference time of reconstruction methods based on wavelet transform and non-wavelet transform in embodiments of the present invention. Detailed Implementation

[0065] 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, and not all embodiments. 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.

[0066] Example 1

[0067] This invention provides a seismic data reconstruction method based on a time-frequency joint self-attention diffusion model with TV regularization constraints. The overall process consists of three main parts, which are described below in conjunction with... Figure 1-2 The details of each part are explained in detail:

[0068] I. Obtain two-dimensional ideal seismic data and construct a seismic dataset for training the diffusion model based on the two-dimensional ideal seismic data.

[0069] In this embodiment, the acquired seismic data includes signals whose amplitude varies over time for each seismic trace. The seismic data is represented in the form of a digital matrix, where each row represents a time sampling point and each column represents a different seismic trace.

[0070] In this embodiment, the earthquake dataset includes pairs of data, namely ideal earthquake data Sei_data_label and training earthquake data Sei_data_train.

[0071] The ideal seismic data Sei_data_label sources in this embodiment include: ① Standard model forward modeling seismic dataset: high-quality synthetic seismic signals without missing values ​​generated based on typical geological structure models using numerical forward modeling methods such as finite difference; ② Actual acquired and processed seismic dataset: reconstruction processing of real acquired data (such as missing value imputation, signal reconstruction, etc.) to simulate missing or damaged conditions in seismic data, forming relatively complete seismic data. These data are used to enhance the model's ability to learn different types of missing data patterns, thereby improving the model's performance in handling missing data and seismic data reconstruction tasks.

[0072] When constructing the training set, the two-dimensional ideal seismic data is first normalized to map the amplitude to a uniform range of [-1,1], which helps to reduce the impact of differences in data amplitude under different acquisition conditions on model training.

[0073] The training seismic data Sei_data_train, i.e., defective data (or missing seismic data), is generated in this embodiment by sampling the ideal seismic data Sei_data_label to simulate the missing traces that occur during actual seismic data acquisition. The specific steps are as follows:

[0074] Ideal seismic data is sampled to generate a corresponding mask matrix. The mask matrix is ​​a binary matrix, where a value of 1 indicates that the original data at that location is preserved, and a value of 0 indicates that the data at that location is missing. Mask generation methods include random missing, regular missing, and continuous missing. Handling of missing locations is achieved through the Hadamard product, which is calculated by multiplying corresponding elements of two matrices. The ideal seismic data (Sei_data_label) and the mask matrix are then compared. The formula for performing the Hadamard product operation is expressed as follows: .in, The sample is from the ideal seismic data Sei_data_labe, i=1,2,3,...; Hadamard product symbol, mask matrix The shape is the same as that of the Sei_data_label training seismic data, with missing trace values ​​being 0 and all others being 1.

[0075] The missing type is determined by the `mask_type` parameter, and the missing strength is controlled by `keep_ratio`. By adjusting `keep_ratio`, you can control the proportion of data retained, thereby setting the strength of the missing data.

[0076] Furthermore, in the later training phase, to increase the diversity of training data, random perturbations are introduced into the mask matrix. Specifically, the mask matrix for each sample is randomly perturbed based on a baseline random seed and the training epoch index, ensuring that the missing locations of different samples within the same training epoch are different, thus enhancing the diversity of training data and improving the robustness of the model. During the validation phase, to ensure the reproducibility of the evaluation results, the mask matrix remains fixed, avoiding random interference during the validation process.

[0077] Finally, before constructing the training set, the training seismic data needs to be preprocessed. Specifically, firstly, both the training seismic data and the ideal complete seismic data are subjected to sliding window cropping, dividing each dataset into two-dimensional blocks of size N×T, where N and T are the preset block height and width, respectively. During this process, each data block is processed using a sliding window method to ensure that each block contains the corresponding regions from both the training and ideal datasets. Next, blocks in the ideal seismic data with less than 20% of their data volume (i.e., blocks with less than 20% of their positions containing valid data) are discarded, and corresponding missing data blocks in the training data are also discarded simultaneously. This method ensures that the training dataset contains only valid blocks with sufficient data volume, avoiding interference from excessive invalid data and improving the quality of the training dataset and the model training effect.

[0078] II. Constructing multi-channel feature maps based on earthquake datasets and using multi-channel feature maps A diffusion model based on a self-attention mechanism is trained to obtain a seismic data reconstruction model. The network structure of this seismic data reconstruction model is illustrated below. Figure 3 As shown. This section's technical content is further divided into three parts: multi-directional feature map hybrid modeling, construction of a diffusion model reconstruction network based on a self-attention mechanism, and setting of the time-frequency joint loss function. Among them, multi-directional feature map hybrid modeling is used to construct multi-channel feature maps based on seismic datasets. A self-attention-based diffusion model reconstruction network is constructed to introduce the self-attention mechanism into the diffusion model and uses multi-channel feature maps. Training is then performed. The time-frequency joint loss function is used to define the reconstruction loss function during the training process of the diffusion model. The following sections detail the three parts:

[0079] 1) Multi-directional feature map hybrid modeling

[0080] In this embodiment, a two-dimensional discrete wavelet transform (DWT) is first performed on the ideal seismic data Sei_data_label and the training seismic data Sei_data_train. Specifically, the ideal seismic data Sei_data_label is used as the original complete seismic data. The training seismic data Sei_data_train is used as the missing seismic data. , respectively and Perform two-dimensional discrete wavelet transform (DWT) to obtain the low-frequency subband. , and multiple high-frequency subbands , The low-frequency subband contains global structural information of the signal, while the high-frequency subband contains detailed information. DWT effectively separates seismic signals of different frequency components, providing multi-directional features for subsequent reconstruction tasks. Furthermore, the high-frequency subband contains local details of the signal in different directions, enhancing these local details.

[0081] Furthermore, in this embodiment, to achieve the fusion of low-frequency and high-frequency information, the multi-directional sub-bands obtained by wavelet transform are subjected to multi-channel splicing processing. Specifically, the low-frequency and high-frequency sub-band feature maps from different frequency sub-bands are spliced ​​as different channels to form a multi-channel input feature map. , Multi-channel feature maps can fully integrate information from different frequency levels, providing rich multi-directional representations for the reconstruction model. This not only enables effective fusion of low-frequency and high-frequency features but also significantly accelerates model convergence and reduces the number of training iterations during subsequent training, thereby effectively improving overall training efficiency while ensuring reconstruction accuracy.

[0082] During training, multi-channel feature maps corresponding to complete seismic data are used. As the initial sample for the forward diffusion process, multi-channel feature maps corresponding to missing seismic data were used. The model is trained using this as a conditional input. During the forward diffusion process, noise is gradually added to the multi-directional features. The noise characteristics at different time steps are obtained. Specifically, for multi-channel wavelet feature maps, Gaussian noise is added at each time step t: ,in Let represent the multidirectional noise characteristics at time step t, q be the conditional probability distribution during the forward diffusion process, N denote a Gaussian distribution, and I be the identity matrix. This represents the noise variance. During the backdiffusion process, the model utilizes information from both low-frequency and high-frequency subbands to gradually recover seismic data from the noise state, preserving the global structure in the low-frequency portion and progressively reconstructing local details in the high-frequency portion. Specifically, from the noise characteristics... Initially, earthquake data is gradually recovered, and a trained diffusion model is used to predict noise: ,in π represents the mathematical multiplication symbol, and s represents the s-th diffusion time step. This represents the noise prediction function learned by the model, which utilizes wavelet multi-directional features to guide inverse recovery. Through the above processing method, the diffusion model can simultaneously combine multi-directional features in the time and frequency domains to achieve smooth recovery in the time domain and detailed reconstruction in the frequency domain, effectively completing the reconstruction of missing seismic data, improving reconstruction accuracy, and maintaining the continuity of the overall structure.

[0083] (ii) Construction of a diffusion model based on self-attention mechanism for network reconstruction

[0084] To overcome the limitation of traditional convolutional layers that can only capture local receptive field features, this embodiment introduces a self-attention mechanism into the diffusion model to enhance the nonlocal interaction of multi-directional seismic data features. A schematic diagram of the self-attention feature transformation is shown below. Figure 4 As shown. Specifically, the wavelet transform multi-directional sub-band feature map, after multi-channel splicing, is flattened into... The sequence is given, where B represents the batch size (i.e., the number of samples input into the model at one time), H represents the height of the feature map, W represents the width of the feature map, and C represents the number of channels. This feature map is then input into the self-attention mechanism module for processing. The self-attention mechanism uses the formula... Calculate the attention matrix and obtain the query result through a linear transformation. and key vector The attention matrix A is used to adjust the value matrix. We perform weighted summation to achieve nonlocal feature interaction at each position in the sequence: .in, This indicates the transpose operation. Represents the key vector Feature dimensions, This indicates that the input sequence will be... Mapped to query vector The linear transformation function, This indicates that the input sequence will be... Mapped to key vector The linear transformation function, This indicates that the input sequence will be... Mapped to value vector The linear transformation function.

[0085] Furthermore, during the forward diffusion process, the self-attention mechanism calculates the similarity between different locations in the input data to maintain global consistency. During the backward diffusion process, the self-attention mechanism enhances the model's non-local dependencies on different locations, enabling the model to rely not only on local information but also fully utilize global information, thereby accurately recovering missing portions of the seismic data.

[0086] Furthermore, the sequence weighted by the self-attention mechanism... The shape is readjusted back to the format of a 2D feature map, input into the diffusion model for reconstruction of missing seismic data, and a complete reconstructed seismic data output is generated. Non-local feature interaction enhances the model's ability to model dependencies between distant seismic traces, thereby improving reconstruction accuracy and the continuity of the overall seismic structure.

[0087] (iii) Setting the joint time-frequency loss function

[0088] In this embodiment, the reconstruction loss function of the overall network is:

[0089]

[0090] The loss function is mainly divided into three parts: ① Loss of missing region reconstruction ② Loss of consistency in the observation area ③ Total variational regularization loss . These are the weighting coefficients.

[0091] The specific forms of each loss function are as follows:

[0092] ①The first step is to calculate the reconstruction loss of the missing seismic data for the missing area. This is used to measure the prediction error of the model in repairing missing data regions. Specifically, it is implemented as follows: Let... This represents the i-th earthquake data element predicted by the model. This represents the i-th element of the ideal seismic data. For the mask matrix, This indicates that the i-th position is missing. This indicates that there is observation data at the i-th position. (Using the formula...) Calculate the reconstruction error of the missing region to ensure that the loss calculation is performed only in the missing region, so that the model can focus on calculating the reconstruction error of the missing region;

[0093] ② Secondly, for the seismic data of the observation area, calculate the consistency loss of the observation area. This is used to ensure that the observed data is not incorrectly modified during training. Specifically, it involves setting a mask matrix. , This indicates that the i-th position is missing. This indicates that data is known at the i-th position. (Using the formula...) Calculate the observation area error to ensure that the loss is calculated only in the observation data area, thereby reducing the interference of the training process on the observation data.

[0094] ③ Finally, total variational regularization (TV) loss is introduced. This is used to constrain the smoothness of the reconstruction results and suppress noise and discontinuities. Specifically, it involves: processing the prediction results... The difference between adjacent data points is penalized to improve the continuity of the reconstruction. The calculation formula is as follows:

[0095] ,

[0096] Where N is the number of samples, and H and W represent the height and width of the data block, respectively. This represents taking the absolute value, used to calculate the gradient between adjacent data points, and penalizes the differences between adjacent data points in the horizontal and vertical directions of the predicted data. This represents the position index of the data block in the height direction, i.e., the [index]. OK; This indicates the position index of the data block in the width direction, i.e., the... List.

[0097] The above loss functions are combined to train the seismic data reconstruction model to ensure that the observed data remain unchanged while recovering the missing data, and to improve the continuity and smoothness of the reconstruction results, thereby improving the overall reconstruction accuracy of the seismic data.

[0098] Third, based on the multi-directional wavelet coefficients output by the seismic data reconstruction model, preliminary prediction results are obtained through inverse wavelet transform. Subsequently, the preliminary prediction results are fused with the original observation data (the actual observation records obtained during seismic acquisition, including missing input seismic data) according to the mask matrix to obtain the final reconstructed complete seismic data.

[0099] Specifically, after the seismic data reconstruction model training converges and the reconstruction accuracy stabilizes, the multi-directional wavelet coefficients output by the seismic data reconstruction model are converted back into time-domain seismic data through inverse wavelet transform (IDWT). This converts the original wavelet domain multi-directional wavelet coefficients into seismic profile data represented on the time axis, which is the preliminary prediction result. Subsequently, the preliminary prediction result is fused with the original observation data using a mask matrix. This means that the original data is retained in the observation area, and the preliminary prediction result is filled in the missing areas. The observation area represents the locations in the input data that have been actually acquired and are not missing; these are marked as "known and valid" in the mask matrix. The original data represents the corresponding actual seismic observation values ​​within the observation area—that is, the seismic data actually recorded by the acquisition equipment, not the model prediction, but the existing real values. The preliminary prediction result itself is an estimate of a complete dataset, including values ​​corresponding to both the observation area and the missing areas; however, during fusion, only the portion corresponding to the missing areas is actually used. Therefore, in practice, the prediction results corresponding to these traces are filled into the missing areas. This yields the final reconstructed complete seismic data, which serves as the model's output.

[0100] The prediction results of this technology can be evaluated by calculating the SNR. Specifically, the training dataset is divided into a training set and a test set in a ratio of 8:2. After training the model on the training set, the test set is used to test the model's reconstruction and restoration performance, and the experimental results are recorded and analyzed according to the evaluation metric. The evaluation metric used is SNR, and the SNR calculation formula is shown below:

[0101] ,

[0102] in, Representing ideal seismic data, This represents the data predicted by the model for recovery. This indicates the number of data points. By calculating the SNR index, the difference between the restored data and the real data can be quantitatively assessed. The higher the SNR value, the closer the model reconstruction result is to the real data, and the higher the restoration quality.

[0103] Example 2

[0104] The experimental platform configuration for this embodiment of the invention is as follows: Intel(R) Core(TM) CPU i7-10700K, 16 GB RAM, NVIDIA GeForce RTX 5070 Ti GPU, PyTorch 2.9 framework, and CUDA version 12.9. The network is updated using the Adam optimizer and undergoes 100 iterations. This embodiment includes the following steps:

[0105] 1. Read the Seg-Y file;

[0106] Read ideal seismic data in Seg-Y format using the ReadSegy function from the toolkit. The calling method is as follows:

[0107] data = ReadSegy(filename);

[0108] To simulate missing regions in seismic data for each training sample, this invention employs a mask matrix generation strategy. The mask matrix generation introduces random perturbations based on the training round index and a baseline random seed. In each training round, if the mask perturbation strategy is enabled (per_sample_jitter=True), an independent sub-random source is derived for each sample, preventing all samples from having identical mask positions within the same training round, thus increasing the diversity of training data and improving the model's generalization ability. For the case where perturbation is disabled (per_sample_jitter=False), a uniform random seed is used to ensure the mask matrix remains stable within the same round, providing reproducible results.

[0109] if self.per_sample_jitter:

[0110] local_rng = np.random.RandomState(self.rng.randint(0, 2**31 - 1))

[0111] else:

[0112] local_rng = self.rng

[0113] 2. Generate the mask matrix;

[0114] Specifically, to simulate missing regions in seismic data, the implementation uses the `make_mask_cols` function. The mask matrix generation considers the specific shape, missing ratio (`keep_ratio`), column gap (`line_gap`), and block width (`block_width`) of each sample. Random perturbation is introduced through `local_rng` to ensure that the mask for each sample is different within the same round. A mask matrix `mask` is generated for each sample; the mask is a binary matrix where 1 indicates that the original data is retained, and 0 indicates missing data. The generated mask matrix is ​​then subjected to a Hadamard product with the original seismic data `x_true` to obtain the missing input data `x_sparse`, which contains only a portion of the observed seismic signals.

[0115] mask=make_mask_cols(

[0116] x_true.shape,

[0117] mask_type=self.mask_type,

[0118] keep_ratio = self.keep_ratio,

[0119] line_gap = self.line_gap,

[0120] block_width = self.block_width,

[0121] rng=local_rng)

[0122] x_sparse=x_true*mask

[0123] Subsequently, the generated mask matrix was compared with the original complete seismic data. The missing seismic input data was obtained by performing Hadamard accumulation. ,in It only includes a portion of the observed seismic signals for subsequent model training.

[0124] 3. Multi-directional feature map hybrid modeling;

[0125] Two-dimensional discrete wavelet transform (DWT) was performed on complete and missing seismic data to obtain the low-frequency subband. With high-frequency subband The low-frequency subband retains global structural information, while the high-frequency subband contains local details in different directions. The DWT decomposition is implemented as follows:

[0126] ll_true,hh_true=self.dwt(x_true)

[0127] ll_sparse,hh_sparse=self.dwt(x_sparse)

[0128] To integrate low-frequency and high-frequency information, the sub-bands are concatenated into a multi-channel input. The concatenated feature maps provide multi-directional time-frequency information for the diffusion model, enhancing its reconstruction capabilities.

[0129] x_input=torch.cat([ll_sparse],hh_sparse,dim=1)

[0130] 4. Construction of the diffusion model reconstruction network based on self-attention mechanism, including the self-attention mechanism module;

[0131] The self-attention mechanism enhances the representational power of input features. It provides a way to establish correlations between different locations in the input tensor, thereby effectively capturing the non-local dependencies of the input data. The specific processing procedure is as follows:

[0132] First, before processing the input feature x_input, the data is normalized using Group Normalization, which helps improve the stability of training and accelerate convergence.

[0133] Secondly, the normalized feature h= Performing four nonlinear transformations, it can be expressed as:

[0134]

[0135] In the formula, NIN represents a module specifically designed for 1×1 convolution and non-linear activation, capable of mapping input channels to a specified dimension and performing feature transformation. This is equivalent to... , indicating that the input feature tensor x is processed by the NIN module and the output size is adjusted to the number of channels y, where y is the number of output channels. is a non-linear activation function, representing the *convolution operation.

[0136] The code for generating q (query), k (key), and v (value) from the normalized h is as follows:

[0137] h = self.GroupNorm_0(x)

[0138] q = self.NIN_0(h)

[0139] k = self.NIN_1(h)

[0140] v = self.NIN_2(h)

[0141] Each location For all locations Calculate the attention weight w, which represents the relevance of this location to all other locations:

[0142]

[0143] in, Indicates the first In each sample, spatial location Spatial location Attention weights Indicates the first In each sample, the location The query vector at that location. Indicates the first In each sample, the location The key vector at the location, where C represents the number of channels.

[0144] The corresponding code implementation is as follows:

[0145] w=torch.einsum('bchw,bcij->bhwij',q,k)*(int(C)**(-0.5))

[0146] w=torch.reshape(w,(B,H,W,H*W))

[0147] w = F.softmax(w, dim = -1)

[0148] w=torch.reshape(w,(B,H,W,H,W))

[0149] Furthermore, the attention matrix is ​​used to perform a weighted summation of the value matrix to achieve global information fusion:

[0150]

[0151] in, Indicates the first In each sample, the location The c-th dimension output feature obtained after attention aggregation Indicates the first In each sample, the location The value feature on the c-th channel.

[0152] The corresponding code is:

[0153] h=torch.einsum('bhwij,bcij->bchw',w,v)

[0154] This module overcomes the problem of limited receptive field in convolution, enabling feature interaction between distant pixels; it integrates global contextual information into local convolutional features to improve the accuracy of missing data recovery; when combined with a diffusion model, it enhances the guiding ability of multi-directional features in the forward and reverse diffusion process, so that low frequencies maintain the global structure and high frequencies accurately reconstruct local details.

[0155] 5. Forward diffusion process;

[0156] During the forward diffusion process, Gaussian noise is gradually added to the multi-channel feature map enhanced by self-attention. In this process, the low-frequency subband ensures the consistency of the global structure, while the high-frequency subband provides guidance for detail recovery, thus achieving multi-directional noise injection.

[0157] beta_t=compute_beta(t) noise=torch.randn_like(x_input)

[0158] x_t=torch.sqrt(1-beta_t)*x_input+torch.sqrt(beta_t)*noise

[0159] 6. Reverse process of diffusion model;

[0160] During the back-diffusion process, the model gradually recovers the seismic data from the noisy state. The trained diffusion model predicts noise. :

[0161] epsilon_theta=self.diffusion_model(x_recon,t)

[0162] x_recon=(x_recon-torch.sqrt(beta_t)*epsilon_theta) / torch.sqrt(1-beta_t)

[0163] 7. Loss function setting;

[0164] Furthermore, in this embodiment, to effectively guide the training of the diffusion model, a reconstruction loss function is defined to measure the model's recovery accuracy in missing regions, consistency of observation areas, and overall continuity. The specific steps are as follows:

[0165] Isotropic total variational regularization loss:

[0166] For the input feature map Calculate the difference between adjacent data points; the horizontal difference is divided into:

[0167] ,

[0168] The vertical difference is divided into:

[0169] ,

[0170] Where B is the batch size, C is the number of channels, and H and W are the height and width of the feature map, respectively;

[0171] Then perform the following calculations for TVs of the same gender: Align to the common area. Then, the gradient magnitude of each pixel is calculated, and the average is taken over all pixels:

[0172]

[0173] in To prevent the use of tiny constants with zero square roots and ensure numerical stability, this invention employs isotropic TV regularization for seismic data reconstruction. This regularization constrains the prediction results by simultaneously calculating the gradient magnitudes in both the horizontal and vertical directions. This allows the model to smooth the reconstruction results while preserving the continuity and local details of seismic bedding, avoiding directional biases that may be introduced by anisotropic processing, and ensuring the accuracy and structural consistency of the reconstructed seismic data. The loss function is:

[0174] deftv_loss(x,isotropic=True,eps=1e-6):

[0175] dx = x[...,1:,:] - x[...,:-1,:]

[0176] dy = x[...,:,1:] - x[...,:,:-1]

[0177] dx_c=dx[...,:,:-1]

[0178] dy_c=dy[...,:-1,:]

[0179] returntorch.sqrt(dx_c**2+dy_c**2+eps).mean()

[0180] The total training loss is formed by combining the missing region reconstruction loss, the observation region consistency loss, and TV regularization. The core code is as follows:

[0181] L_missing=((x_recon-x_true).abs()*(1-mask)).sum() / (1-mask).sum().clamp_min(1.0)

[0182] L_obs=((x_recon-x_sparse).abs()*mask).sum() / mask.sum().clamp_min(1.0)

[0183] L_tv = tv_loss(x_recon)

[0184] rec_loss=L_missing+lam_obs*L_obs+lam_tv*L_tv

[0185] rec_loss.backward()

[0186] optimizer.step()

[0187] During training, the loss function guides model learning through the combined effects of time domain (TV regularization) and frequency domain (wavelet multi-directional features).

[0188] 8. Inverse wavelet transform and output fusion;

[0189] The predicted wavelet coefficients are reconstructed from the time-domain seismic data using inverse wavelet transform (IDWT), and then fused with the original observation data according to the mask matrix to obtain the final fully reconstructed seismic data.

[0190] X_rec_idwt=idwt2(X_recon)

[0191] X_final=X_rec_idwt*(1-mask)+x_true*mask

[0192] 9. Performance evaluation;

[0193] The recovery results are quantified, and the signal-to-noise ratio (SNR) is used as the evaluation index: the higher the SNR, the higher the reconstruction accuracy. By comparing the SNR under different missing types and missing ratios, the reconstruction effectiveness and stability of this method can be verified.

[0194] signal_power=(X_true**2).mean()

[0195] noise_power=((X_final-X_true)**2).mean()

[0196] snr=10*torch.log10(signal_power / noise_power)

[0197] Through the above steps, this embodiment closely integrates multi-directional wavelet features, self-attention mechanism, and the forward / reverse processes of diffusion model to achieve high-precision reconstruction of missing seismic data while maintaining the continuity of global structure and local details.

[0198] Implementation results:

[0199] First, the reconstruction results of different methods are presented under three typical scenarios: random missing data, regular missing data, and continuous missing data. Using the same missing dataset, the reconstruction and recovery results in this embodiment are as follows: Figure 6 , 8 As shown in Figures 1 and 10, the reconstruction and restoration results of the comparison methods are as follows: Figure 5 , 7As shown in Figures 9 and 1, by comparing the reconstruction results of the two methods pairwise, it is clear that the method of this invention achieves a higher degree of similarity between the reconstructed seismic record and the actual data, with continuous event structure and uniform energy distribution. In contrast, the comparative method exhibits obvious waveform breaks and energy attenuation in the missing regions, resulting in poorer reconstruction performance. Taking a 30% random missing portion as an example... Figure 5 The comparison method exhibits significant lateral discontinuities and a lack of high-frequency details; Figure 6 As a result of the method of the present invention, the layered structure is continuous, the events are clear, and the recovered area is highly consistent with the real data.

[0200] Figure 11 – Figure 13 This figure shows a comparison of the SNR values ​​of the validation set under different missing scenarios in this embodiment. As can be seen from the figure, the reconstruction and restoration performance of the method of this invention (including wavelet transform, self-attention mechanism, and regularization constraints) is significantly better than the comparative methods in all missing scenarios. With the increase of training rounds, the SNR value of the method of this invention shows an upward trend with small fluctuations, demonstrating the stability of the training process and good convergence performance. Compared with the comparative methods, its SNR growth rate is significantly faster, and the final SNR value at convergence is also higher. This indicates that the method of this invention, under the constraints of the self-attention mechanism and TV regularization term, can effectively improve the signal-to-noise ratio in the early stages of training, accelerate model convergence, and maintain a stable performance advantage in the later stages. Especially in more challenging scenarios such as regular missing and continuous missing, the advantages of the method of this invention are more prominent, demonstrating that it can still maintain strong reconstruction ability and stability under complex missing scenarios.

[0201] Figure 14 The inference times of two reconstruction methods, one including wavelet transform and the other not, were compared at different batch sizes. The results show that while the times are similar for small batches, the wavelet transform strategy significantly reduces inference time as the batch size increases, especially at 256 batches where the difference is nearly double. The multi-directional features of wavelets reduce the spatial dimension of the diffusion model input, significantly improving inference efficiency and meeting the practical processing needs of large-scale seismic data.

[0202] Example 3

[0203] This invention also provides a seismic data reconstruction system based on a time-frequency joint self-attention diffusion model with TV regularization constraints, which mainly includes three components: a dataset construction module, a reconstruction model module, and a data reconstruction module. Specifically:

[0204] The dataset building module is used to acquire two-dimensional ideal seismic data and build a seismic dataset for training diffusion models based on the two-dimensional ideal seismic data.

[0205] The reconstruction model module is used to construct multi-channel feature maps based on earthquake datasets, and to train a diffusion model based on a self-attention mechanism using the multi-channel feature maps to obtain an earthquake data reconstruction model.

[0206] The data reconstruction module is used to obtain preliminary recovery results of complete seismic data based on the seismic data reconstruction model. Subsequently, the predicted results are fused with the original observation data according to the mask matrix to obtain the final reconstructed complete seismic data.

[0207] The reconstruction model module is further divided into a feature map unit, a model training unit, and a loss function unit, specifically:

[0208] The feature map unit is used to construct multi-channel feature maps based on seismic datasets.

[0209] The model training unit is used to introduce the self-attention mechanism into the diffusion model and is trained using multi-channel feature maps.

[0210] The loss function unit is used to set the reconstruction loss function during the training process of the diffusion model.

[0211] The embodiments described above are merely preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Various modifications and improvements made to the technical solutions of the present invention by those skilled in the art without departing from the spirit of the present invention should fall within the protection scope defined by the claims of the present invention.

Claims

1. A seismic data reconstruction method based on a time-frequency joint self-attention diffusion model with TV regularization constraints, characterized in that, Includes the following steps: Two-dimensional ideal seismic data is acquired, and a seismic dataset for training a diffusion model is constructed based on the two-dimensional ideal seismic data; the seismic dataset includes pairs of data, namely ideal seismic data and training seismic data; The method for constructing an earthquake dataset for diffusion model training based on the two-dimensional ideal earthquake data includes: normalizing the ideal earthquake data; sampling the normalized ideal earthquake data to generate a corresponding mask matrix, and then performing a Hadamard product operation on the ideal earthquake data and the mask matrix to generate the training earthquake data; and preprocessing the ideal earthquake data and the training earthquake data to form the final earthquake dataset. Based on the earthquake dataset, a multi-channel feature map is constructed, and a diffusion model based on a self-attention mechanism is trained using the multi-channel feature map to obtain the earthquake data reconstruction model. The method for constructing the multi-channel feature map based on the earthquake dataset includes: using the ideal earthquake data as the original complete earthquake data and the training earthquake data as the missing earthquake data, performing two-dimensional discrete wavelet transforms on each to obtain multi-directional sub-bands, including a low-frequency sub-band. , and high-frequency subband , The multi-directional sub-bands are spliced ​​together in multiple channels to form the multi-channel feature map. The method of introducing a self-attention mechanism into a diffusion model and training it using the multi-channel feature map includes: flattening the multi-channel feature map into... The sequence of features is processed by a self-attention mechanism. Here, B represents the batch size (number of samples input into the model at one time), H represents the height of the feature map, W represents the width of the feature map, and C represents the number of channels. The self-attention mechanism uses the attention matrix to perform a weighted summation of the value matrix, achieving non-local feature interaction at each position in the sequence. During the forward diffusion process of the diffusion model, Gaussian noise is added to the multi-channel feature map at each time step t. During the backward diffusion process of the diffusion model, Gaussian noise is added from the noise features... Initially, the seismic data is gradually restored, noise is predicted using a trained diffusion model, and the multi-directional wavelet features are used to guide the reverse recovery. Based on the multi-directional wavelet coefficients output by the earthquake data reconstruction model, preliminary prediction results are obtained through inverse wavelet transform. Subsequently, the preliminary prediction results are fused with the original observation data according to the mask matrix to obtain the final reconstructed complete earthquake data.

2. The method according to claim 1, characterized in that, The earthquake data reconstruction model consists of three parts: multi-directional feature map hybrid modeling, diffusion model reconstruction network construction based on self-attention mechanism, and time-frequency joint loss function setting. The multi-directional feature map hybrid modeling is used to construct multi-channel feature maps based on the earthquake dataset; The self-attention mechanism-based diffusion model reconstruction network is used to introduce the self-attention mechanism into the diffusion model and is trained using the multi-channel feature map to establish a seismic data reconstruction model. The time-frequency joint loss function is used to set the reconstruction loss function during the training process of the diffusion model.

3. The method according to claim 2, characterized in that, The reconstruction loss function is: ， in, Losses due to reconstruction of missing areas; This represents the loss of consistency across the observation area. For total variational regularization loss, These are the weighting coefficients.

4. A seismic data reconstruction system based on a time-frequency joint self-attention diffusion model with TV regularization constraints, wherein the system applies the method described in any one of claims 1-3, characterized in that, It includes a dataset construction module, a reconstruction model module, and a data reconstruction module; The dataset construction module is used to acquire two-dimensional ideal seismic data and construct an earthquake dataset for training the diffusion model based on the two-dimensional ideal seismic data; The reconstruction model module is used to construct a multi-channel feature map based on the earthquake dataset, and use the multi-channel feature map to train a diffusion model based on a self-attention mechanism to obtain the earthquake data reconstruction model. The data reconstruction module is used to obtain preliminary prediction results based on the multi-directional wavelet coefficients output by the seismic data reconstruction model through inverse wavelet transform. Subsequently, the preliminary prediction results are fused with the original observation data according to the mask matrix to obtain the final reconstructed complete seismic data.

5. The system according to claim 4, characterized in that, The reconstruction model module includes a feature map unit, a model training unit, and a loss function unit; The feature map unit is used to construct a multi-channel feature map based on the earthquake dataset; The model training unit is used to introduce the self-attention mechanism into the diffusion model and to train it using the multi-channel feature map to establish a seismic data reconstruction model. The loss function unit is used to set the reconstruction loss function during the training process of the diffusion model.

Images