VSP shaft wave suppression method based on physical constraint improved U-Net

Abstract

The invention discloses a VSP wellbore wave suppression method for improving U-Net based on physical constraints, which comprises the following steps: introducing apparent velocity constraints of an F-K transform domain, constructing physical constraint conditions, simultaneously realizing time sequence feature extraction and prediction of wellbore waves in combination with a GRU network, and learning a mapping relation between a noise image and a clean image by using the network so as to improve the suppression efficiency of the wellbore waves. The details and the structure of a noise image are accurately restored through a training network, the difference between a real de-noised image and network output is minimized so as to achieve the purpose of image de-noising, the model retains effective signals through a residual learning mechanism, and meanwhile precise suppression of shaft wave noise is achieved. The method provided by the invention can accurately identify and effectively suppress wellbore waves and reflected waves thereof in a complex VSP wave field, significantly improve the signal-to-noise ratio, reduce residual noise, improve data processing efficiency while improving data quality, provide technical support for high-precision seismic imaging under a complex geological condition, and have a wide application prospect. And meanwhile, a feasible solution is provided for popularization and application in an actual work area.

Description

Technical Field

[0001] The present invention belongs to the technical field of oil and gas exploration, and particularly relates to a VSP wellbore wave suppression method based on physically constrained improved U-Net. Background Art

[0002] In the fine exploration of deep and complex exploration targets in China, the vertical seismic profile (VSP) technology collects seismic signals by using sensors in a well drilled deep underground, and realizes high-resolution underground structure imaging at a relatively low cost. With its rich wavefield information and high-quality data performance, the VSP technology is receiving more and more extensive attention. With the increase in the complexity of exploration targets and the improvement of the requirements for wavefield processing accuracy and efficiency, how to effectively suppress the strong noise interference in VSP data and achieve high-quality wavefield separation has become the research focus in the field of signal processing.

[0003] The wellbore wave is a kind of guided wave that propagates along the interface between the fluid (mud) and the solid (wellbore wall, etc.) in the wellbore, and has characteristics such as strong energy and slow attenuation. Especially during the reflection process from the bottom of the well to the liquid surface, the wellbore wave may be reflected multiple times and interfere with the effective signal. Conventional methods suppress the wellbore wave from two aspects: data acquisition and signal processing. In terms of data acquisition, measures such as avoiding small offset measurements, adding inflatable airbags around the geophones, and placing effective baffles between the geophones and the wellbore wall are usually adopted to reduce the generation of wellbore waves or enhance the isolation effect of interference signals. In terms of signal processing, the wellbore wave is usually suppressed through various filtering techniques. Among them, the filtering method based on apparent velocity has less damage to the effective signal, but has limited suppression effect on the wellbore wave and has a large residual energy, and is prone to erroneously filtering out shear wave signals with similar apparent velocities. Although the method based on wavelet analysis has certain applicability, the selection of window parameters is relatively complex and subjective in practical applications. Existing wellbore wave suppression methods all have certain limitations and deficiencies, and the effective suppression of wellbore waves is still a difficult point in VSP fidelity signal processing and a research direction that needs to be studied in the field of signal processing. Summary of the Invention

[0004] To solve the above technical problems, the present invention proposes a VSP wellbore wave suppression method based on physically constrained improved U-Net. By introducing the apparent velocity constraint in the F-K transform domain into the existing U-Net network architecture, physical constraint conditions are constructed based on the differences in apparent velocity and the frequency spectrum in the F-K domain. At the same time, dynamic mechanisms such as GRU are introduced, and the wellbore wave suppression without manual parameter selection is realized through few-shot learning.

[0005] The technical solution adopted by the present invention is as follows: a VSP wellbore wave suppression method based on physically constrained improved U-Net, and the specific steps are as follows:

[0006] S1. Establish a velocity model, perform forward simulation to solve the elastic wave equation by the finite difference method, then use a three-component geophone to record the three-component gather, generate pure and noise-free VSP data, and then obtain a synthetic dataset by adding background noise with different intensities and synthetic borehole wave data with different intensities.

[0007] Among them, the particle vibration directions received by the three-component geophone are horizontal and vertical directions, obtaining the vibration components of the wave field along the Z direction and the X direction, and generating pure and noise-free VSP data.

[0008] The pure and noise-free VSP data is added with background noise of different intensities, and the relevant parameters of the borehole wave are adjusted to simulate complex geological conditions and noise interference characteristics, generating borehole waves and borehole reflection wave data with different starting positions and energy intensities to construct a synthetic dataset.

[0009] Among them, the relevant parameters of the borehole wave include: the amplitude, waveform, and spectral characteristics of the borehole wave.

[0010] S2. Preprocess the actual VSP data, then use the F-K filtering method to extract the borehole wave noise data and the effective signal data without borehole waves, and then use the τ-P transform domain filtering to solve the dispersion phenomenon generated by the F-K filtering, obtaining the actual dataset, that is, more accurate borehole wave noise data and effective signal data without borehole waves.

[0011] Among them, the actual VSP data is the VSP data with borehole wave noise obtained from the actual work area; the preprocessing includes: removing bad traces and horizontal correction.

[0012] S3. Fuse the synthetic dataset in step S1 with the actual dataset in step S2 to construct a total dataset, and slice the total dataset using the sliding slice method, that is, enhance the full-wave field data through the sliding slice technique.

[0013] S4. Classify and organize the total dataset after slicing in step S3, complete the construction of the sample and label dataset, and then perform normalization processing on the sample data and the label data respectively.

[0014] S5. Construct an improved U-Net network model, expand the data in step S4 after normalization into a four-dimensional structure as the original input, and divide the normalized data into a training set and a validation set through K-Fold cross-validation. Use the training set to train the model and the validation set to evaluate the model. According to the evaluation results, repeat the training steps until a trained improved U-Net network model is obtained, that is, the target noise suppression and effective signal prediction model.

[0015] Among them, the improved U-Net network model includes: an input layer, an encoder part, a skip connection, a decoder part, an apparent velocity physical constraint layer, a GRU layer, and an output layer.

[0016] The improved U-Net network model takes the original noisy data and borehole wave noise data after slicing and normalization as inputs, combines the structure of U-Net, time series processing, and apparent velocity physical constraints, uses residual learning, and denoises the borehole wave noise by learning the borehole wave characteristics.

[0017] S6. Based on the trained improved U-Net network model obtained in step S5, use the newly collected actual VSP data in real time to verify the model generalization.

[0018] Furthermore, the specific steps of step S3 are as follows:

[0019] Fuse the synthetic data set in step S1 and the actual data set in step S2 to construct a total data set, and set the VSP data in the total data set to be a time series matrix D∈R composed of multiple receivers M×N .

[0020] Among them, M represents the number of time sampling points, N represents the number of seismic traces, and R represents the complex domain.

[0021] The sliding slicing method adopts a rectangular M×N sliding slicing method. Each slice includes M rows and N columns, the sliding step size in the X direction is S x =16, and the sliding step size in the Y direction is S y =128. Then the extraction formula for each small block is as follows:

[0022] P i，j =D[j:j+M, i:i+N] (1)

[0023] Among them, P i , j represents a small block extracted from the matrix D, the position is determined by the indexes i and j, and the local area in the matrix D has been extracted; D[j:j+M, i:i+N] represents extracting a sub-matrix from the j-th row to the j+M-1-th row and from the i-th column to the i+N-1-th column from the matrix D.

[0024] Then slice the total data set according to the sliding step sizes of S x =16 and S y =128. The overlap degree between each small slice is 50%, and the borehole wave slice data with a size of 256×32 and the original slice data containing borehole waves are obtained.

[0025] Furthermore, the specific steps of step S4 are as follows:

[0026] S41. Classify and organize the total data set after slicing in step S3 to complete the construction of the sample and label data sets;

[0027] Among them, the original noisy data is the sample, and the borehole wave noise data is the label.

[0028] S42. Based on step S41, normalize the maximum value of the data amplitude, that is, perform normalization processing on the sample data and the label data respectively;

[0029] Based on step S41, calculate the absolute maximum value of the amplitude in the input sample and label gather data respectively, and normalize the input sample and label data to the range of [-1, 1] with their respective maximum values.

[0030] The input data is D ∈ R M×N , calculate the absolute maximum value of the sample / label amplitude according to formula (2), and use formula (3) for normalization. The expressions are as follows:

[0031] T = max(|D i , l |) (2)

[0032]

[0033] Among them, |D i，l | represents the absolute value of the element in the i-th row and l-th column of matrix D, that is, the absolute value of the amplitude in the wave field. T represents its absolute maximum value. If T = 0, return the original wave field data. If T ≠ 0, normalize each data to obtain the normalized matrix D'.

[0034] Furthermore, the specific steps of step S5 are as follows:

[0035] The improved U-Net network model expands the sliced and normalized original noisy data and borehole wave noise data into a four-dimensional structure as the original input, and uses K-Fold cross-validation to divide the normalized data into a training set and a validation set.

[0036] Among them, during the model training process, different hyperparameter settings are selected to determine the best training effect. The hyperparameters include: learning rate, batch size, and the number of hidden layer units.

[0037] K-Fold cross-validation divides the data set into K non-overlapping subsets. In each cross-validation, K - 1 subsets are used as the training set, and the remaining one subset is used as the validation set for model evaluation. Assume that the data set has a total of L samples. Divide the data set Z into K subsets, and the size of each subset is L / K.

[0038] Among them, K = 5, and the division ratio of the training set to the validation set is 4:1.

[0039] In the K-th verification, the subset Z is used K denotes the validation set, and the remaining Z\Z K denotes the training set. Then the data training set data validation set has the following mathematical expressions:

[0040]

[0041] Repeat the verification process of formula (4) K times, each time using a different validation set, and finally summarize the evaluation results of all validation sets to obtain the overall performance of the model.

[0042] Let X denote the total dataset after slicing and normalization, X train denotes the training set, X val denotes the validation set, split idx denotes the partitioning index, and the expression is as follows:

[0043]

[0044] X train = X[split idx X val = X[split idx :] (6)

[0045] where Q represents the total number of samples in the dataset.

[0046] In the total dataset after slicing and normalization, the sample data and label data are paired one by one to form training pairs, and the DataLoader is used to load the dataset for training and validation, provided in batches of a fixed size Batch, while shuffling the data loading order. The batch data expression is as follows:

[0047]

[0048] where Batch q represents the data of the q-th batch, x high,i and x low,i represent the sample and label data pairs respectively, B represents the size of each batch, q represents the index of the batch, that is, the current is the q-th batch, represents the total number of batches and rounds up.

[0049] The random order of the data is achieved through the permutation index π, and the expression is as follows:

[0050] Shuffle((X) = Permute(X, π) π represents the random permutation index (8)

[0051] Further, in step S5, the improved U-Net network model is specifically as follows:

[0052] (1) Input layer;

[0053] The input layer is used to receive the sliced VSP original noisy data and the wellbore wave data. Taking the sliced original noisy data as samples and the wellbore wave noise data as labels as the training input, when inputting, for each two-dimensional data sample, it is expanded into a four-dimensional structure, and the samples and labels are in one-to-one correspondence, and the order of data loading is shuffled.

[0054] Among them, the dimensions of the four-dimensional structure respectively represent the batch size, the number of channels, the data height, and the width.

[0055] (2) Encoder part;

[0056] The encoder part of the network includes: feature extraction and downsampling of the input data.

[0057] Among them, each layer performs feature extraction and fusion through multi-resolution convolution, and in each level of feature extraction, global context information is extracted through global average pooling, and the channel weights are adjusted. Downsampling is to reduce the resolution of the feature map through operations while retaining the main features, and finally output the feature map for the decoding stage, which is implemented using max pooling MaxPool2d, that is, taking the maximum value within the local window of the feature map as the output, and the expression is as follows:

[0058]

[0059] Among them, y(i, j) represents the output feature element after pooling at position (i, j), window(i,j) represents the 2x2 window centered on (i, j), x(s, t) represents the pixel value within the window, and the final feature map size is 1 / 2 of the input, and the width and height change as Height in represents the height of the input feature map, Weight in represents the width of the input feature map, Height out represents the height of the output feature map after pooling, Weight out represents the width of the output feature map after pooling.

[0060] (3) Skip connection;

[0061] The skip connection is used to fuse the feature map of the encoder part with the feature map of the corresponding layer of the decoder part, and retain the low-level feature information, and the expression is as follows:

[0062]

[0063] Among them, Represents the input of the skip connection, Represents the output calculated through the intermediate layer, Represents that the final output contains the original input and the output after transformation.

[0064] (4) GRU layer;

[0065] In GRU, Update Z t and Reset R t Control the dynamic adjustment of the hidden state, through the current input x t and the hidden state H at the previous time step t-1 The weighted sum of, via W z and W R Passed to the Sigmoid activation function σ for calculation, the expression is as follows:

[0066] Z t = σ(W z [H t-1 , x t ) (11)

[0067] R t = σ(W R [H t-1 , x t ) (12)

[0068] Among them, W z Represents the weight matrix for calculating the update gate Z t , W R Represents the weight matrix for calculating the update gate R t ; Z t Determines how much of the current hidden state is inherited from the previous state. The closer the value is to 1, the more information from the previous state is retained. The closer it is to 0, the more it focuses on the current input; R t Controls the influence of the previous hidden state on the candidate state.

[0069] Candidate hidden state Combined with the current input x t and the previous hidden state H t-1 Calculated, S represents the weight matrix, the expression is as follows:

[0070]

[0071] Final hidden state H t Is obtained by the weighted sum of the candidate state and the previous hidden state H t-1 The expression is as follows:

[0072]

[0073] In a bidirectional GRU, the forward GRU and the backward GRU process sequence information simultaneously, generating the forward hidden state and the backward hidden state respectively. The two are combined through concatenation or weighted combination to obtain the final hidden state h t , and the expression is as follows:

[0074]

[0075]

[0076]

[0077] (5) Decoder part;

[0078] The decoder is used to gradually map the deep features extracted by the encoder back to the spatial resolution of the original input to generate the final output, specifically including: upsampling, feature map concatenation, and feature fusion.

[0079] The decoder performs bilinear interpolation upsampling multiple times through Up sample. After each upsampling, the number of channels is reduced through the Conv2d and BatchNorm2d modules. Combining with the corresponding feature maps of the encoder, feature fusion is achieved through torch.cat. On the feature fusion result after each upsampling, a multi-resolution convolution module is applied to further extract and reconstruct features. The specific calculation expression is as follows:

[0080] x concat = Concat(x up , x encoder ) (18)

[0081] x out = Activation(Conv2D(x concat ) (19)

[0082] x SE (c) = x(c).σ(W 2 * ReLU(W 1 * GAP(x))) (20)

[0083] Among them, x concat represents the feature map fused after upsampling and skip connection at a certain layer of the decoder, x up represents the feature map after upsampling, x encoder represents the feature map of the corresponding layer in the encoder, x out represents the output map of the corresponding layer in the encoder; x SE(c) represents the c-th channel of the feature map after SE module processing, x(c) represents the c-th channel of the input feature map, σ represents the Sigmoid activation function, W 2 represents the weight matrix of the fully connected layer, which is used to map the feature map generated by W 1 and GAP operation to the final channel weights, W 1 represents the weight matrix of another fully connected layer, which is used to transform the feature map obtained by GAP operation into an intermediate representation, and GAP(x) represents global average pooling.

[0084] (6) Apparent velocity physical constraint layer;

[0085] The apparent velocity physical constraint layer calculates the difference between the predicted noise and the label noise and their corresponding apparent velocity errors by transforming the predicted noise and the label noise into the frequency-wavenumber domain, so as to improve the accuracy of the model output with physical constraints.

[0086] First, the input data is transformed by F-K transform to obtain the spectral values spec in the frequency f and wavenumber k x domain, and then the apparent velocity v is calculated. The expression is as follows:

[0087]

[0088] By transforming the input predicted noise and label noise data from the time-space domain to the frequency-wavenumber domain, and then calculating the spectrum and the apparent velocity error, the loss at this time consists of two parts: the spectrum difference between the predicted noise and the label noise in the F-K domain and the apparent velocity difference between the predicted value and the label value. The expression is as follows:

[0089]

[0090] Among them, A represents the number of samples, B represents the batch size, C represents the number of channels, H i , W i represents the frequency dimension and wavenumber dimension in the F-K domain, represents the spectral value of the predicted noise, represents the spectral value of the label noise, f h represents the frequency component at the h-th time.

[0091] The apparent velocity is the ratio of the frequency f and the wavenumber k x , and the apparent velocity loss is the mean square error between the predicted value and the label value. The expression is as follows:

[0092]

[0093]

[0094] Among them, represents the apparent velocity of the predicted noise, Represents the apparent velocity of label noise. When k x , w = 0, the corresponding position can be ignored to avoid the error of division by zero.

[0095] Frequency - wavenumber domain loss Total FK-loss Calculated by weighted sum of two parts, FK loss and Velocity Loss The expression is as follows:

[0096] Total FK-loss = Velocity Loss + FK loss (25)

[0097] (7) Output layer;

[0098] The output layer is used to map the last feature map of the network to the desired output space. The output layer converts the multi - channel feature map into the target output noise through a 1x1 convolution. For the input feature map F, the conversion expression is as follows:

[0099]

[0100] Among them, the size of the feature map F is B×C×H×W, H represents the height of the data, and W represents the width of the data; W conv represents the convolution kernel weight, b conv represents the bias of the convolutional layer, represents the predicted output.

[0101] Finally, the predicted noise signal is converted into two - dimensional data, the predicted noise signal is output, and physical rationality constraints are imposed in the frequency domain through the apparent velocity physical constraint layer. At this time, the data distribution range is in the normalized range [-1, 1]. Then the predicted output data is restored to its original amplitude range, that is, the original amplitude of the data is restored by multiplying by the maximum amplitude T used during normalization. The expression is as follows:

[0102] D = D′×T (27)

[0103] Among them, D′ represents the predicted output after normalization, T represents the absolute maximum amplitude of the corresponding label data during the normalization process, and D represents the data after denormalization, that is, the target predicted output is restored to the original amplitude range.

[0104] (8) Loss function design;

[0105] Compare the predicted noise signal with the true signal, and comprehensively design MSE, MAE and R 2Optimize the model with a three - part loss function, and control their contributions to the total loss by defining the weights of each loss: Weight_mse, Weight_mae, Weight_r2. The expression is as follows:

[0106]

[0107]

[0108]

[0109] Among them, represents the predicted value, y i represents the true value, n represents the number of samples; NR and NT represent the number of geophones and the number of signal samplings respectively, m represents the m - th sampling point in the time series, E[S(x j , t m ) represents the average value of all data in a single shot gather record, S(x j , t m ) represents the labeled signal, that is, the true value of the seismic signal at position x j and time t m , S pre represents the predicted signal obtained through the model, and the magnitude of this value is close to 0, S pre and S represent the predicted and target shot gather records. R 2 reflects the goodness of fit of the prediction, and the value range is (-∞, 1].

[0110] Set two - part losses in the time - space domain Total time loss and the frequency - wavenumber domain Total FK-loss , and add them weighted according to the dynamically adjusted weight ratio FK weight to form the total loss Total loss , and the expression is as follows:

[0111] Total loss = Total time loss + FK weight * Total FK-loss (31)

[0112] Advantages of the present invention: By introducing apparent velocity constraints in the F-K transform domain, the method of the present invention constructs physical constraint conditions. At the same time, combined with the GRU network, it realizes the extraction and prediction of the temporal characteristics of borehole waves, and uses the network to learn the mapping relationship between the noise image and the clean image. By training the network, it accurately restores the details and structure of the noise image, minimizes the difference between the real denoised image and the network output, so as to achieve the purpose of image denoising. The model retains the effective signal through the residual learning mechanism, and at the same time realizes the precise suppression of borehole wave noise. By fully exploring the physical information in the borehole wave noise, the method of the present invention applies physical constraint conditions in combination with the differences in apparent velocity and the frequency spectrum in the F-K domain, constructs a constraint model applicable to various VSP wave fields, can accurately identify and effectively suppress borehole waves and their reflected waves in complex VSP wave fields, significantly improve the signal-to-noise ratio, reduce the residual noise, improve the data processing efficiency while improving the data quality, provides technical support for high-precision seismic imaging under complex geological conditions, and at the same time provides a feasible solution for the popularization and application in actual work areas. Brief Description of the Drawings

[0113] Figure 1 It is a flow chart of a method for suppressing VSP borehole waves based on physically constrained improved U-Net of the present invention.

[0114] Figure 2 It is an architecture diagram of the improved U-Net network in the embodiment of the present invention.

[0115] Figure 3 It is a schematic diagram of the velocity model in the embodiment of the present invention.

[0116] Figure 4 It is a schematic diagram of the synthetic data in the embodiment of the present invention.

[0117] Figure 5 It is a schematic diagram of adding noise samples in the embodiment of the present invention.

[0118] Figure 6 It is a schematic diagram of the samples and labels after data slicing in the embodiment of the present invention.

[0119] Figure 7 It is a schematic diagram of the skip connection in the embodiment of the present invention.

[0120] Figure 8 It is a schematic diagram of the GRU module in the embodiment of the present invention.

[0121] Figure 9 It is a schematic diagram of the FK Constraint Layer module in the embodiment of the present invention.

[0122] Figure 10 It is a schematic diagram of the refined processing samples, borehole wave noise, effective signal and energy statistics in the embodiment of the present invention.

[0123] Figure 11 It is a schematic diagram of the original noisy data, wellbore wave and effective signal, and energy statistics under the existing filtering method in the embodiment of the present invention.

[0124] Figure 12 It is a schematic diagram of the predicted residual data, predicted wellbore wave data, predicted effective signal, and energy statistics by the U-Net method in the embodiment of the present invention.

[0125] Figure 13 It is a schematic diagram of the predicted residual data, predicted wellbore wave data, predicted effective signal, and energy statistics by the method of the present invention in the embodiment of the present invention. Detailed implementation manners

[0126] The method of the present invention will be further described below in conjunction with the accompanying drawings and embodiments.

[0127] As Figure 1 shown, the flowchart of a VSP wellbore wave suppression method based on physically constrained improved U-Net of the present invention is as follows:

[0128] S1. Establish a velocity model, solve the elastic wave equation by forward simulation using the finite difference method, then use a three-component geophone to record the three-component gather, generate pure noise-free VSP data, and then obtain a synthetic dataset by adding background noise with different intensities and synthetic wellbore wave data with different intensities.

[0129] Among them, the particle vibration directions received by the three-component geophone are horizontal and vertical directions, and the vibration components of the wave field along the Z direction (vertically downward) and the X direction (horizontally to the right) are obtained to generate pure noise-free VSP data.

[0130] The pure noise-free VSP data is added with background noise of different intensities, and the relevant parameters of the wellbore wave are adjusted to simulate complex geological conditions and noise interference characteristics, and wellbore waves and wellbore reflection wave data with different starting positions and energy intensities are generated to construct a synthetic dataset.

[0131] Among them, the relevant parameters of the wellbore wave include: the amplitude, waveform, and spectral characteristics of the wellbore wave.

[0132] S2. Preprocess the actual VSP data, then use the F-K filtering method to extract the wellbore wave noise data and the effective signal data without wellbore waves, and then use the τ-P transform domain filtering to solve the dispersion phenomenon generated by the F-K filtering to obtain the actual dataset, that is, more accurate wellbore wave noise data and effective signal data without wellbore waves.

[0133] Among them, the actual VSP data is the VSP data containing borehole wave noise obtained from the actual work area; the preprocessing includes: removing bad channels and horizontal correction.

[0134] S3. Fuse the synthetic data set in step S1 with the actual data set in step S2 to construct a total data set, and slice the total data set using the sliding slice method, that is, enhance the full-wavefield data through the sliding slice technique;

[0135] S4. Classify and organize the total data set sliced in step S3 to complete the construction of the sample and label data sets, and then perform normalization processing on the sample data and label data respectively;

[0136] S5. Construct an improved U-Net network model, expand the data in step S4 after normalization into a four-dimensional structure as the original input, and divide the normalized data into a training set and a validation set through K-Fold cross-validation. Use the training set to train the model, use the validation set to evaluate the model, and according to the evaluation results, repeat the training steps until an improved U-Net network model that is trained is obtained, that is, the target noise suppression and effective signal prediction model;

[0137] Among them, as Figure 2 shown, the improved U-Net network model includes: an input layer, an encoder part (Encoder), skip connections (Skip Connections), a decoder part (Decoder), a FK constraint layer (FK ConstraintLayer), a GRU layer, and an output layer.

[0138] For the design and training of the improved U-Net model, this embodiment combines the encoder-decoder architecture with skip connections, which not only retains the image details but also extracts rich features. Aiming at the suppression problem of borehole waves and borehole reflection waves in the complex VSP wavefield, the model takes the sliced noisy data (samples) and borehole wave noise (labels) as inputs, introduces the apparent velocity constraint in the F-K transform domain, and establishes physical constraint conditions by combining the differences in apparent velocity and the frequency spectra in the F-K domain, so that it participates in the loss calculation.

[0139] At the same time, combine the GRU network to capture the temporal characteristics of borehole waves and borehole reflection waves, and further improve the modeling effect of temporal information. In order to more accurately suppress the borehole wave noise, the model learns the mapping relationship between the noise image and the clean image, and minimizes the difference between the network output and the real denoised image during training. Combining these improvements, the model realizes the accurate suppression of borehole wave noise while retaining the details and structure of the target signal.

[0140] S6. Based on the trained improved U-Net network model obtained in step S5, use the newly collected actual VSP data in real time to verify the model generalization ability.

[0141] In this embodiment, the specific steps of step S1 are as follows:

[0142] The diversity of training samples and the accuracy of label data play important roles in the performance and generalization ability of the network. To expand the training data set, in this embodiment, based on the velocity model, synthetic pure noise-free VSP data is generated by forward simulating the elastic wave equation.

[0143] As Figure 3 shown, a dipping layered medium containing 6 different velocity layers is constructed. When simulating VSP data with variable well-source distances, the wellbore is placed at the position of the red frame, the horizontal coordinate is 3500 meters, the depth range is from 400 meters to 3500 meters, and the seismic source is placed at the 5000-meter position. A P-wave source with a dominant frequency of 30 Hz Ricker wavelet is used, and 161 geophones are arranged at 20-meter intervals underground.

[0144] The elastic wave simulation is carried out in a two-dimensional underground medium to obtain the vibration components of the wave field in the Z direction (vertically downward) and the X direction (horizontally to the right). Figure 4 Show the synthetic single-component data. Further, different intensities of noise are added to the synthetic data, and relevant parameters are adjusted to simulate complex geological conditions and noise interference characteristics, as Figure 5 shown.

[0145] Among them, Figure 4 (a) and (d) are the X and Z components of synthetic data 1 respectively, Figure 4 (b) and (e) are the X and Z components of synthetic data 2 respectively, Figure 4 (c) and (f) are the X and Z components of synthetic data 3 respectively; Figure 5 (a) is a sample of synthetic data containing wellbore wave noise (weak), Figure 5 (b) is a sample of synthetic data containing wellbore wave noise (strong), Figure 5 (c) is a sample of synthetic data containing wellbore and reflected wave noise, Figure 5 (d) is the label of wellbore wave noise (weak), Figure 5 (e) is the label of wellbore wave noise (strong), Figure 5 (f) is the label of wellbore and reflected wave noise.

[0146] In this embodiment, the specific steps of step S3 are as follows:

[0147] Fuse the synthetic data set in step S1 and the actual data set in step S2 to construct a total data set. It is set that the VSP data in the total data set is represented as a time series matrix D∈R composed of multiple receivers M×N. The time dimension is usually greater than the number of seismic traces, and the data usually appears in the form of a spatio-temporal two-dimensional matrix. The time axis is called the long axis, and the spatial axis is the short axis, and there are obvious differences in their lengths. The overall data presents a narrow and long strip data structure.

[0148] Among them, M represents the number of time sampling points, N represents the number of seismic traces, and R represents the complex domain.

[0149] The sliding slice method adopts a rectangular M×N sliding slice method. Each slice includes M rows and N columns, and the sliding step size in the X direction is S x = 16, and the sliding step size in the Y direction is S y = 128. Then the extraction formula for each small block is as follows:

[0150] P i，j = D[j:j+M, i:i+N] (1)

[0151] Among them, P i , j represents a small block extracted from the matrix D, and its position is determined by the indexes i and j. The local area in the matrix D has been extracted; D[j:j+M, i:i+N] represents extracting a submatrix from the j-th row to the j+M-1-th row and from the i-th column to the i+N-1-th column from the matrix D.

[0152] By adjusting the sliding step size, more training samples can be obtained, thereby increasing the diversity of the data. In this embodiment, the total data set is sliced according to the sliding step sizes of S x = 16 and S y = 128. The overlap degree between each small slice is 50%, and wellbore wave slice data with a size of 256×32 and original slice data containing wellbore waves are obtained. The sliced data is shown as Figure 6 shown.

[0153] Among them, Figure 6 (a) and (f) are respectively the sample and label schematic diagrams of slice data 1, Figure 6 (b) and (g) are respectively the sample and label schematic diagrams of slice data 2, Figure 6 (c) and (h) are respectively the sample and label schematic diagrams of slice data 3, Figure 6 (d) and (i) are respectively the sample and label schematic diagrams of slice data 4, Figure 6 (e) and (j) are respectively the sample and label schematic diagrams of slice data 5.

[0154] In this embodiment, the specific steps of step S4 are as follows:

[0155] S41. Classify and sort the total data set sliced in step S3 to complete the construction of the sample and label data sets;

[0156] Among them, the original noisy data is the sample, and the borehole wave noise data is the label.

[0157] S42. Based on step S41, normalize the maximum value of the data amplitude, that is, perform normalization processing on the sample data and the label data respectively;

[0158] To help the model better learn and train, normalization is used to scale the data into a standard interval to ensure that different features have the same dimension, thereby improving the stability and efficiency of model training. Seismic data includes positive and negative amplitudes, indicating the direction of wave propagation (such as the forward and backward propagation of waves). If normalized to [0, 1], the positive and negative information of the data will be lost. Normalizing to the range [-1, 1] can retain the positive and negative information of the amplitude, enabling the model to more accurately capture the physical characteristics of the data and pay more attention to the relative changes in the data.

[0159] Based on step S41, calculate the absolute maximum value of the amplitude in the input sample and the label gather data respectively, and normalize the input sample and the label data to the range [-1, 1] with their respective maximum values, avoiding uneven gradients caused by differences in dimensions between different features, and thus accelerating the convergence speed. At this time, ensure that the scales of all input data are consistent, and the model can fairly learn the contribution of each feature.

[0160] The input data is D ∈ R M×N , calculate the absolute maximum value of the sample / label amplitude according to Equation (2), and perform normalization using Equation (3). The expressions are as follows:

[0161] T = max(|D i , l |) (2)

[0162]

[0163] Among them, |D i，l | represents the absolute value of the element in the i-th row and the l-th column of the matrix D, that is, the absolute value of the amplitude in the wave field. T is its absolute maximum value. If T = 0, return the original wave field data. If T ≠ 0, then normalize each data to obtain the normalized matrix D'. This normalization method can ensure that the relative changes in the wave field data are evenly distributed in the value range [-1, 1], thereby enhancing the model's ability to capture relative differences.

[0164] In this embodiment, the specific steps of step S5 are as follows:

[0165] The improved U-Net network model expands the sliced and normalized original noisy data and borehole wave noise data into a four-dimensional structure as the original input, enabling the data to be effectively input into the convolutional neural network model for training and prediction. K-Fold cross-validation is adopted to divide the normalized data into a training set and a validation set.

[0166] During the model training process, different hyperparameter settings are selected to determine the best training effect. The hyperparameters include: learning rate, batch size, and the number of hidden layer units.

[0167] K-Fold cross-validation is a commonly used technique for validating the model effect. It divides the dataset into K non-overlapping subsets. In each cross-validation, K - 1 subsets are used as the training set, and the remaining one subset is used as the validation set for model evaluation. Suppose the dataset has L samples. The dataset Z is divided into K subsets, and the size of each subset is L / K.

[0168] Among them, K = 5, and the division ratio of the training set to the validation set is 4:1.

[0169] In the Kth validation, the subset Z K represents the validation set, and the remaining Z\Z K represents the training set. Then the data training set data validation set The mathematical expressions are as follows:

[0170]

[0171] Repeat the validation process in Equation (4) K times, each time using a different validation set, and finally summarize all the validation set evaluation results to obtain the overall performance of the model.

[0172] Suppose X represents the total dataset after slicing and normalization, X train represents the training set, X val represents the validation set, and split idx represents the division index. The expressions are as follows:

[0173]

[0174] X train = X[split idx X val = X[split idx :] (6)

[0175] Among them, Q represents the total number of samples in the dataset.

[0176] In the total dataset after slicing and normalization, the sample data and label data are paired one by one to form training pairs. The dataset is loaded using DataLoader for training and validation, provided in batches of a fixed size (Batch), and at the same time, the data loading order is shuffled to enhance the generalization of the model. The batch data expression is as follows:

[0177]

[0178] Among them, Batch q represents the data of the q-th batch, x high,i and x low,i represent the sample and label data pairs respectively, B represents the size of each batch (Batch size), q represents the index of the batch, that is, the current is the q-th batch, represents the total number of batches and is rounded up.

[0179] The random order of the data is achieved through the permutation index π, and the expression is as follows:

[0180] Shuffle((X)=Permute(X, π) π represents the random permutation index (8)

[0181] In this embodiment, in the step S5, the improved U-Net network model is specifically as follows:

[0182] (1) Input layer;

[0183] The input layer is used to receive the sliced VSP original noisy data and wellbore wave data. In this embodiment, the network slices the sampling time and the number of seismic traces with a sliding step of 50% each, and the overlap degree between each small slice is 50%. After slicing, the original noisy data after slicing is used as the sample, and the wellbore wave noise data is used as the label as the training input. When inputting, for each two-dimensional data sample, it is expanded into a four-dimensional structure. By corresponding the sample and the label one by one, and at the same time, to ensure the randomness of training, the data loading order is shuffled, so that the data can be effectively input into the convolutional neural network for training and prediction.

[0184] Among them, the dimensions represent the batch size (Batch Size), the number of channels (Channel), the data height (Height), and the width (Width) respectively.

[0185] (2) Encoder part (Encoder);

[0186] The encoder part of the network includes: feature extraction and downsampling of the input data.

[0187] Among them, each layer performs feature extraction and fusion through multi-resolution convolution to enhance the expression ability for complex wellbore waves and reflected waves. Since the time axis is regarded as the long axis and the space axis is the short axis, the overall data presents a narrow strip shape. In this case, using a rectangular convolution kernel is more adaptable to the data size. During the training process, the fidelity and similarity of the waveform are the key difficulties and focuses. To enhance the attention to key features, global context information is extracted through global average pooling at each level of feature extraction, and the channel weights are adjusted to enhance the attention to key features. Downsampling reduces the resolution (width and height) of the feature map through operations to compress information, reduce the computational amount, and retain the main features. Finally, the feature map for the decoding stage is output. The specific implementation uses max pooling MaxPool2d, that is, taking the maximum value within the local window of the feature map as the output. The expression is as follows:

[0188]

[0189] Among them, y(i, j) represents the output feature element after pooling at position (i, j), window(i,j) represents the 2x2 window centered on (i, j), x(s, t) represents the pixel value within the window, and the final feature map size is 1 / 2 of the input, and the width and height changes are Height in represents the height of the input feature map, Weight in represents the width of the input feature map, Height out represents the height of the output feature map after pooling, Weight out represents the width of the output feature map after pooling.

[0190] (3) Skip connection;

[0191] As Figure 7 shown, skip connections (Skip Connections) can effectively retain low-level feature information by fusing the feature maps of the encoder part with the feature maps of the corresponding layers in the decoder part. The encoder is used to extract the wellbore wave features of the original noisy data. The decoder part gradually restores the spatial resolution of the image. At the same time, through skip connections, feature fusion is performed with the corresponding layers of the encoder. The output feature map of the encoder is concatenated or added to the upsampled feature map in the decoder, enhancing the detail expression ability of the network when restoring the spatial resolution. It helps to slow down the loss of information in the deep network and improve the training efficiency and performance of the model.

[0192]

[0193] Among them, represents the input of the skip connection, represents the output calculated through the intermediate layer, Indicates that the final output contains the output after the original input is kernel-transformed.

[0194] (4) GRU layer;

[0195] In GRU, Update Z t and Reset R t control the dynamic adjustment of the hidden state. They are obtained by the weighted sum of the current input x t and the hidden state H at the previous time step t-1 , passed through W z and W R , and calculated by passing through the Sigmoid activation function σ. The expression is as follows:

[0196] Z t = σ(W z [H t-1 , x t ) (11)

[0197] R t = σ(W R [H t-1 , x t ) (12)

[0198] Among them, W z represents the weight matrix for calculating the update gate Z t , and W R represents the weight matrix for calculating the update gate R t . Z t (Update Gate) determines how much of the current hidden state is inherited from the previous state; the closer the value is to 1, the more information from the previous state is retained, and the closer it is to 0, the more attention is paid to the current input. R t (Reset Gate) controls the influence of the previous hidden state on the candidate state.

[0199] Candidate hidden state is calculated by combining the current input x t and the previous hidden state H t-1 . S represents the weight matrix. The expression is as follows:

[0200]

[0201] The final hidden state H t is obtained by the weighted sum of the candidate state and the previous hidden state H t-1 . The expression is as follows:

[0202]

[0203] Such as Figure 8As shown, in this embodiment, in the bidirectional GRU (BiGRU), the forward GRU and the reverse GRU process sequence information simultaneously and generate forward hidden states and reverse hidden states respectively. The two are combined through concatenation or weighted combination to obtain the final hidden state h t , and the expression is as follows:

[0204]

[0205]

[0206]

[0207] For wellbore wave noise, there are wellbore reflection waves in some data, and their development time is slightly later than that of wellbore waves. In order to better enable the network to recognize and judge the relationship between wellbore waves and their reflection waves, and at the same time achieve accurate prediction of them, adding a GRU layer with temporal dependence can effectively extract features in the time dimension, thereby further enhancing the model's representation ability for the complex relationship between wellbore waves and their reflection waves.

[0208] (5) Decoder part (Decoder);

[0209] The decoder is used to gradually map the deep features extracted by the encoder back to the spatial resolution of the original input to generate the final output. This process includes upsampling, feature map concatenation, and feature fusion.

[0210] The decoder performs bilinear interpolation upsampling multiple times through Up sample. After each upsampling, the number of channels is reduced through the Conv2d and BatchNorm2d modules. Combining the corresponding feature maps of the encoder (skip connection), the fusion of features is achieved through torch.cat. On the feature fusion result after each upsampling, a multi-resolution convolution module is applied to further extract and reconstruct features. These modules ensure that the decoder can gradually restore the spatial resolution of the data through convolution operations, while fusing low-level and high-level features and reducing the interference of irrelevant information. The specific calculation expression is as follows:

[0211] x concat = Concat(x up , x encoder ) (18)

[0212] x out = Activation(Conv2D(x concat ) (19)

[0213] xSE (c) = x(c)·σ(W 2 *ReLU(W 1 *GAP(x))) (20)

[0214] Among them, x concat represents the feature map fused after upsampling and skip connection on a certain layer of the decoder, x up represents the feature map after upsampling, x encoder represents the feature map of the corresponding layer in the encoder, x out represents the output map of the corresponding layer in the encoder; x SE (c) represents the c-th channel of the feature map processed by the SE module, x(c) represents the c-th channel of the input feature map, σ represents the Sigmoid activation function, W 2 represents the weight matrix of the fully connected layer, used to map the feature map generated by W 1 and the GAP operation to the final channel weights, W 1 represents the weight matrix of another fully connected layer, usually used to convert the feature map obtained by the GAP operation (compressed global feature) into an intermediate representation, GAP(x) represents global average pooling.

[0215] (6) Apparent Velocity Physical Constraint Layer (FK Constraint Layer);

[0216] As Figure 9 shown, the FK Constraint Layer calculates the difference between the two and its corresponding apparent velocity error by converting the predicted noise and label noise to the frequency - wavenumber (F - K) domain, so as to improve the accuracy of the model output with physical constraints.

[0217] First, the input data is transformed through the F - K transform to obtain the spectral value spec in the frequency f and wavenumber k x domain, and then the apparent velocity v is calculated according to the formula. The expression is as follows:

[0218]

[0219] By converting the input predicted noise and label noise data from the time - space domain to the frequency - wavenumber (F - K) domain, and calculating the error of the spectrum and apparent velocity. Here, the loss consists of two parts: the spectral difference between the predicted noise and label noise in the F - K domain and the apparent velocity difference between the predicted value and the label value. The expression is as follows:

[0220]

[0221] Among them, A represents the number of samples, B represents the batch size, C represents the number of channels, H i , W i represent the frequency dimension and wavenumber dimension in the F-K domain, represents the spectral value of the predicted noise, represents the spectral value of the labeled noise, f h represents the frequency component at the h-th time.

[0222] The apparent velocity is the ratio of the frequency f and the wavenumber k x The apparent velocity loss is the mean square error between the predicted value and the labeled value, and the expression is as follows:

[0223]

[0224]

[0225] Among them, represents the apparent velocity of the predicted noise, represents the apparent velocity of the labeled noise. When k x , w = 0, the corresponding position can be ignored to avoid the error of dividing by zero.

[0226] The frequency-wavenumber domain loss Total FK-loss is calculated by weighting two parts, FK loss and Velocity Loss The apparent velocity loss constrains the global propagation characteristics (apparent velocity) of the noise and ensures the kinematic information of the overall borehole wave. The FK spectral loss emphasizes the fitting of local spectral characteristics, making the frequency and wavenumber of the noise closer to the label. The weighted total loss comprehensively considers the global and local characteristics in the optimization to improve the generalization ability of the model.

[0227] Total FK-loss = Velocity Loss + FK loss (25)

[0228] The apparent velocity physical constraint layer ensures that the frequency-wavenumber characteristics of the noise are consistent with the actual physical characteristics through the joint constraint of the spectrum and the apparent velocity, thereby optimizing the performance of the model.

[0229] (7) Output layer;

[0230] The output layer is used to map the last feature map of the network to the required output space. The output layer usually converts the multi-channel feature map into the noise of the target output through a 1x1 convolution. Then, for the input feature map F, the conversion is as follows:

[0231]

[0232] Among them, the size of the feature map F is B×C×H×W, where C represents the number of input channels (Channels), H represents the height of the data, and W represents the width of the data; W conv represents the convolutional kernel weight, and b conv represents the bias of the convolutional layer, represents the predicted output.

[0233] Finally, the predicted noise signal is converted into two-dimensional data, and the predicted noise signal is output and physically reasonable constraints are imposed in the frequency domain through the Frequency Constraint Layer (FK Constraint Layer) to ensure that the model output satisfies both data characteristics and physical laws. At this time, the data distribution range is in the range of [-1, 1] after normalization. To restore the predicted output data to its original amplitude range, the original amplitude of the data is restored by multiplying by the maximum amplitude T used during normalization. The expression is as follows:

[0234] D = D′×T (27)

[0235] Among them, D′ represents the predicted output after normalization, T represents the absolute maximum amplitude of the corresponding label data during the normalization process, and D represents the data after denormalization, that is, the target predicted output is restored to the original amplitude range.

[0236] (8) Design of the loss function;

[0237] Since the apparent velocity of the borehole wave is constant, using the energy distinguishability between the signal and the noise to extract the effective signal and suppress the noise is the main idea of signal processing. Utilizing this apparent velocity characteristic, in this embodiment, the apparent velocities of the borehole wave and the noise are used to participate in the loss calculation so that the predicted noise can better learn the kinematic characteristics of the apparent velocity of the label borehole wave. Through F-K domain constraints, the learning ability of the model for frequency and wavenumber information is improved to ensure that the noise predicted by the model is consistent with the real noise in the spatio-temporal domain and the apparent velocity.

[0238] At the same time, considering the specific complex spatio-temporal structure of VSP data, the spatio-temporal domain can capture the signal changes in the time series and spatial distribution, thereby more accurately evaluating the difference between the predicted result and the real signal. In the borehole wave noise denoising task, the spatio-temporal domain loss can effectively handle the non-linear changes of the noise in time and space and improve the denoising effect. The predicted noise signal is compared with the real signal, and MSE, MAE, and R are comprehensively considered 2A three-part loss function is used to optimize the model. By defining the weights of each loss (Weight_mse, Weight_mae, Weight_r2), their contributions to the total loss are controlled, enabling it to simultaneously recover effective signals and suppress noise on the spatio-temporal scale, thereby obtaining a more accurate denoising effect. The calculation expression is as follows:

[0239]

[0240]

[0241]

[0242] Among them, MAE and MSE are used to measure the absolute error and squared difference between the predicted value and the label value, dealing with data outliers. represents the predicted value, and y i represents the true value, n represents the number of samples. NR and NT respectively represent the number of geophones and the number of signal samplings, m represents the m-th sampling point in the time series, E[S(x j , t m ) represents the average value of all data in a single shot gather record, S(x j , t m ) represents the label signal, that is, the true value of the seismic signal at position x j and time t m , S pre represents the predicted signal obtained through the model, and the magnitude of this value is close to 0. S pre and S represent the predicted and target shot gather records. R 2 reflects the goodness of fit of the prediction, and the value range is (-∞, 1].

[0243] To meet the settings of the denoising network, by means of few-shot learning, it gets rid of traditional parameter dependencies, fully explores the physical information of regular noise to construct constraint conditions, considering the Total in the time-space domain time loss and the Total in the frequency-spectrum - wavenumber domain FK-loss The two parts of the loss are weighted and added according to the dynamically adjusted weight ratio FK weight to form the total loss Total loss , and the expression is as follows:

[0244] Total loss = Total time loss + FK weight * Total FK-loss (31)

[0245] In the early stage of training, the weight of the frequency-domain loss is large, and the model pays more attention to the frequency constraint; as the training progresses, the weight decreases, and the model pays more attention to the fitting of time-domain features. By integrating the time-domain and frequency-domain losses, the weighting mechanism helps the model balance the optimization of time-domain and frequency-domain features at different training stages.

[0246] In this embodiment, the step S6 is specifically as follows:

[0247] Use the newly acquired actual VSP data in real time to test the trained improved U-Net network model (target noise suppression and effective signal prediction model) obtained in step S5 to verify its generalization performance in practical applications. By suppressing the borehole wave in the actual VSP data with noise, evaluate whether the model can effectively adapt to the actual geological conditions and accurately suppress the borehole wave, so as to prove the reliability and practicability of the model, and obtain high-fidelity effective signals on the premise of minimizing the loss of effective signals.

[0248] In this embodiment, a set of untrained VSP data is selected for generalization verification, and the energy statistics are completed by calculating the sum of the squares of the amplitude values. The energy of the original noisy data is about 9126, and the energy of the labeled borehole wave is about 1037.28. The effective signal can be obtained by subtracting the borehole wave noise from the original noisy wave field, and its energy is about 8107.18. The schematic diagram and the energy statistics results are as Figure 10 shown.

[0249] Based on the original noisy data, for the purpose of comparing the experimental effects, three methods are used for comparison to perform noise suppression and effective signal extraction, as Figure 11 , Figure 12 and Figure 13 shown, which respectively show the wave fields and energy statistics of the noisy, borehole wave noise, and effective signals under the three methods.

[0250] Among them, the three methods are: the existing filtering method (F-K filtering), the U-Net method, and the method of the present invention.

[0251] Figure 11 The wave fields before and after noise suppression obtained by using the existing filtering method (F-K filtering) are shown. It can be clearly seen that although the F-K filtering has a certain effect on suppressing the borehole wave, it is obvious that some useful shear wave signals are filtered out at the same time. At this time, the waveform fidelity of the obtained noise and effective signals is not high, which is also the disadvantage of the existing filtering method for suppressing the borehole wave.

[0252] After using the U-Net network for noise suppression, the predicted noise and effective signals shown in Figure 12 can be obtained. And when using the method of the present invention for noise suppression, the corresponding predicted noise and effective signals shown in Figure 13 can be obtained. By comparing Figure 12 andFigure 13 For the grayscale image and energy statistical data of the residual results, it can be found that the residual of the noise predicted by the U-Net method is 82.083, while the residual of the noise predicted by the method of the present invention is only 2.93. The former residual value is about 28.01 times that of the method of the present invention.

[0253] Further comparing the predicted noise results of the three methods, it can be found that the noise energy predicted by the U-Net method is significantly lower than that of the method of the present invention. The existing F-K filtering method will damage some shear wave signals when suppressing the borehole wave and cannot guarantee the waveform fidelity. In contrast, the method of the present invention can better maintain the waveform characteristics of the borehole wave when suppressing the borehole wave noise and shows good results in the prediction of kinematic characteristics. The analysis of the denoised effective signal obtained by subtracting the predicted noise from the original noisy data shows that the effective signal fidelity of the U-Net method and the existing F-K filtering method is significantly lower than that of the method of the present invention.

[0254] In addition, there are obvious differences in the energy changes of the three methods. The method of the present invention is significantly superior to the U-Net method and the existing F-K filtering method in terms of predicted energy value, signal fidelity, and denoising effect. As shown in Table 1, by comparing the noise and effective signal data of the original accurate labels and combining evaluation indexes such as SSIM (structural similarity index), IOU (intersection over union), and R 2 and other evaluation indexes for analysis, the method of the present invention shows obvious advantages, especially in maintaining the high fidelity of the effective signal and the accuracy of noise suppression.

[0255] Table 1

[0256] Network F-K filtering U-Net The method of the present invention <![CDATA[Wellbore wave R 2 > 0.9228 0.9375 0.9549 <![CDATA[Valid signal R 2 > 0.9316 0.9466 0.9591 Effective signal SSIM 0.8351 0.8371 0.8590 Noise SSIM 0.8102 0.8493 0.8687 Effective signal Iou 0.8101 0.8395 0.8364 Noise Iou 0.8024 0.8392 0.8407

[0257] As can be seen from the table, the method of the present invention not only effectively improves the accuracy of borehole wave suppression, but also greatly reduces signal distortion in the process of extracting effective signals, which is superior to the current mainstream noise suppression methods. Based on this experimental verification, it can be concluded that the method of the present invention has significant performance advantages in noise suppression and effective signal recovery and has strong generalization ability.

[0258] In summary, the method of the present invention is based on an improved U-Net network architecture, which combines an encoder-decoder structure with skip connections to extract rich features while retaining signal details. By performing apparent velocity constraint filtering in the F-K domain, accurate judgment and suppression of the borehole wave and its reflected wave are achieved. By introducing a GRU time series module, the feature learning ability of the model for the signals before and after the borehole wave is significantly improved, enabling the determination of whether there is a borehole reflected wave in the original noisy data, and making the identification, prediction, and suppression of the borehole wave more efficient. Moreover, the method of the present invention uses a residual learning mechanism to compare the features of the borehole wave with the original noisy data, learn and extract the features of the borehole wave, and achieve efficient suppression of the borehole wave from the noisy data. Compared with existing methods, residual learning can quickly extract the target effective signal with less damage to the effective signal, providing an efficient and simple means for accurate suppression of the borehole wave. Compared with existing methods, the present method performs better in terms of the phase consistency and waveform similarity of the prediction, further verifying the high fidelity of the method of the present invention for predicting the borehole wave and its reflected wave. Analysis of the prediction results shows that the method of the present invention can better achieve accurate suppression of the borehole wave and restoration of the effective signal, providing a new method for borehole wave processing and VSP data analysis.

[0259] Those of ordinary skill in the art will realize that the embodiments described herein are for helping the reader understand the principles of the present invention, and it should be understood that the protection scope of the present invention is not limited to such specific statements and embodiments. For those skilled in the art, various modifications and changes can be made to the present invention. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included within the scope of the claims of the present invention.

Claims

1. A VSP wellbore wave suppression method based on physical constraint improved U-Net, the specific steps are as follows: S1. Establish a velocity model, solve the elastic wave equation by forward simulation using the finite difference method, and then use a three-component geophone to record three-component gathers to generate pure and noise-free VSP data. Then, by adding background noise of different intensities and synthetic wellbore wave data of different intensities, a synthetic data set is obtained. in, The particle vibration directions received by the three-component geophone are horizontal and vertical, and the vibration components of the wave field along the Z direction and the X direction are obtained, generating pure and noise-free VSP data; Pure and noise-free VSP data is added with background noise of different intensities, and the wellbore wave related parameters are adjusted to simulate complex geological conditions and noise interference characteristics, generating wellbore waves and wellbore reflection wave data with different starting positions and energy intensities to construct a synthetic data set; Among them, the borehole wave related parameters include: the amplitude, waveform and spectrum characteristics of the borehole wave; S2. Preprocess the actual VSP data, extract the borehole wave noise data and the data without borehole wave effective signal by using the FK filtering method, and then use the τ-P transform domain filtering to solve the dispersion phenomenon caused by the FK filtering and obtain the actual data set, that is, more accurate borehole wave noise data and data without borehole wave effective signal; The actual VSP data is the VSP data including wellbore wave noise acquired in the actual work area; the preprocessing includes: removing bad tracks and horizontal correction; S3, fusing the synthetic data set of step S1 with the actual data set of step S2 to construct a total data set, and slicing the total data set using a sliding slicing method, that is, enhancing the full wavefield data by using the sliding slicing technology; S4, classify and sort the total data set after slicing in step S3, complete the construction of sample and label data sets, and then normalize the sample data and label data respectively; S5, constructing an improved U-Net network model, expanding the data normalized in step S4 into a four-dimensional structure as the original input, and dividing the normalized data into a training set and a validation set through K-Fold cross validation, using the training set to train the model, using the validation set to evaluate the model, and repeating the training steps according to the evaluation results until a trained improved U-Net network model is obtained, i.e., a target noise suppression and effective signal prediction model; The improved U-Net network model includes: an input layer, an encoder part, a skip connection, a decoder part, an apparent velocity physical constraint layer, a GRU layer and an output layer; The improved U-Net network model takes the original normalized noisy data after slicing and the wellbore wave noise data as input, combines the structure of U-Net, time series processing, and physical constraints of apparent velocity, uses residual learning, and denoises the wellbore wave noise by learning the wellbore wave characteristics; S6. Based on step S5, the trained improved U-Net network model is obtained, and the generalization of the model is verified using new actual VSP data collected in real time.

2. A VSP wellbore wave suppression method based on physical constraint improved U-Net according to claim 1, characterized in that: The step S3 is specifically as follows: The synthetic data set of step S1 is fused with the actual data set of step S2 to construct the total data set. The VSP data in the total data set is assumed to be a time series matrix D∈R composed of multiple receivers. M×N ; Among them, M represents the number of time sampling points, N represents the number of seismic traces, and R represents the complex domain; The sliding slicing method adopts a rectangular M×N sliding slicing method, each slice includes M rows and N columns, and the sliding step length in the X direction is S x =16, the sliding step length in the Y direction is S y =128, then the extraction formula for each small block is as follows: P i，j =D[j:j+M,i:i+N] (1) Among them, P i , j It represents a small block extracted from the matrix D. The position is determined by indexes i and j. The local area in the matrix D has been extracted. D[j:j+M, i:i+N] means extracting a submatrix from row j to row j+M-1 and from column i to column i+N-1 from the matrix D. Then follow S x =16 and S y The total data set is sliced ​​with a sliding step size of =128, and the overlap between each small slice is 50%, so as to obtain wellbore wave slice data and wellbore wave original slice data with a size of 256×32.

3. A VSP wellbore wave suppression method based on physical constraint improved U-Net according to claim 1, characterized in that: The step S4 is specifically as follows: S41, classify and sort the total data set after slicing in step S3, and complete the construction of sample and label data sets; Among them, the original noisy data is the sample, and the borehole wave noise data is the label; S42, based on step S41, normalizing the maximum value of the data amplitude, that is, normalizing the sample data and the label data respectively; Based on step S41, the absolute maximum amplitude values ​​in the input sample and label gather data are calculated respectively, and the input sample and label gather data are normalized to the range of [-1, 1] with their respective maximum values; Input data, i.e. D∈R M×N , calculate the absolute maximum value of the sample / label amplitude according to formula (2), and normalize it using formula (3), the expression is as follows: T=max(|D i , l |) (2) Among them, |D i，l | represents the absolute value of the element in the i-th row and the l-th column in the matrix D. The absolute value of the amplitude in the wave field. T represents its absolute maximum value. If T=0, the original wave field data is returned. If T≠0, each data is normalized to obtain the normalized matrix D′.

4. A VSP wellbore wave suppression method based on physical constraint improved U-Net according to claim 1, characterized in that: The step S5 is specifically as follows: The improved U-Net network model expands the normalized original noisy data and borehole wave noise data after slicing into a four-dimensional structure as the original input, and divides the normalized data into a training set and a validation set by using K-Fold cross validation; During the model training process, different hyperparameter settings are selected to determine the best training effect. The hyperparameters include: learning rate, batch size, and number of hidden layer units; K-Fold cross-validation divides the data set into K non-overlapping subsets. In each cross-validation, K-1 subsets are used as training sets, and the remaining subset is used as a validation set for model evaluation. Assuming that the data set has a total of L samples, the data set Z is divided into K subsets, and the size of each subset is L / K. Among them, K=5, the ratio of training set to validation set is 4:1; In the Kth validation, use subset Z K represents the validation set, and the remaining Z\Z K represents the training set; then the data training set Data validation set The mathematical expression is as follows: Repeat the validation process of formula (4) K times, using a different validation set each time, and finally summarize all the validation set evaluation results to obtain the overall performance of the model; Let X represent the total dataset after slicing and normalization, X train represents the training set, X val Represents the validation set, split idx Indicates partition index, the expression is as follows: X train =X[split idx ] X val =X[split idx :] (6) Among them, Q represents the total number of samples in the data set; After slicing and normalization, the total data set is normalized to match sample data and label data one by one to form training pairs. The data set is loaded by DataLoader for training and verification. It is provided in batches of fixed size and the order of data loading is disrupted. The batch data expression is as follows: Among them, Batch q represents the data of the qth batch, x high,i and x low,i Represent the sample and label data pairs respectively, B represents the size of each batch, q represents the index of the batch, that is, the current one is the qth batch, Indicates the total number of batches and rounds up; The random order of data is achieved by permuting the index π, as expressed below: Shuffle((X)=Permute(X,π) π represents a random permutation index (8).

5. A VSP wellbore wave suppression method based on physical constraint improved U-Net according to claim 1, characterized in that: In step S5, the improved U-Net network model is specifically as follows: (1) Input layer; The input layer is used to receive the original noisy VSP data and wellbore wave data after slicing. The original noisy data after slicing is used as samples, and the wellbore wave noise data is used as labels as training inputs. During input, for each two-dimensional data sample, it is expanded into a four-dimensional structure, and the samples and labels are matched one by one, and the order of data loading is disrupted. The four-dimensional structure dimensions represent batch size, number of channels, data height and width respectively; (2) Encoder part; The encoder part of the network includes: feature extraction and downsampling of the input data; Among them, each layer performs feature extraction and fusion through multi-resolution convolution, and at each level of feature extraction, global context information is extracted through global average pooling, and channel weights are adjusted; downsampling is to reduce the resolution of the feature map through operations while retaining the main features, and finally output the feature map for the decoding stage, which is implemented using the maximum pooling MaxPool2d, that is, taking the maximum value in the local window of the feature map as the output, the expression is as follows: Among them, y(i, j) represents the output feature element after pooling at position (i, j), window(i, j) represents the 2x2 window centered at (i, j), x(s, t) represents the pixel value in the window, and the final feature map size is 1 / 2 of the input, and the width and height change is Height in Indicates the input feature map height, Weight in Indicates the width and height of the input feature map. out Indicates the height of the output feature map after pooling, Weight out Indicates the width of the output feature map after pooling; (3) Skip connections; The jump connection is used to fuse the feature map of the encoder part with the feature map of the corresponding layer of the decoder part, retaining the low-level feature information. The expression is as follows: in, represents the input of the skip connection, Indicates that the output is calculated through the intermediate layer, Indicates that the final output contains the original input kernel transformed output; (4) GRU layer; In GRU, Update Z t and Reset R t Control the dynamic adjustment of the hidden state through the current input x t and the hidden state H at the previous time step t-1 The weighted sum of z and W R Pass it to the Sigmoid activation function σ and calculate it. The expression is as follows: Z t =σ(W z [H t-1 ,x t ]) (11) R t =σ(W R [H t-1 ,x t ]) (12) Among them, W z Indicates that the update gate Z is calculated t The weight matrix, W R Indicates that the update gate R is calculated t The weight matrix of t Determines how much the current hidden state inherits from the previous state. The closer the value is to 1, the more previous information is retained, and the closer it is to 0, the more attention is paid to the current input; R t Control the influence of the previous hidden state on the candidate state; Candidate hidden states Combined with the current input x t and the previous hidden state H t-1 Calculation, S represents the weight matrix, the expression is as follows: The final hidden state H t The candidate state and the previous hidden state H t-1 The weighted summation is as follows: In the bidirectional GRU, the forward GRU and the reverse GRU process the sequence information simultaneously, generating the forward hidden state and the reverse hidden state respectively. and the reverse hidden state The two are concatenated or weighted to obtain the final hidden state h t , the expression is as follows: (5) decoder part; The decoder is used to gradually map the deep features extracted by the encoder back to the spatial resolution of the original input to generate the final output, including upsampling, feature map concatenation and feature fusion. The decoder performs multiple bilinear interpolation upsampling through Up sample. After each upsampling, the channel is reduced in dimension through Conv2d and BatchNorm2d modules. Combined with the corresponding feature map of the encoder, feature fusion is achieved through torch.cat. On the feature fusion result after each upsampling, a multi-resolution convolution module is applied to further extract and reconstruct features. The specific calculation expression is as follows: x concat =Concat(x up ,x encoder ) (18) x out =Activation(Conv2D(x concat ) (19) x SE (c)=x(c).σ(W2*ReLU(W1*GAP(x))) (20) Among them, x concat It represents the feature map fused after upsampling and skip connection at a certain layer of the decoder, x up represents the feature map after upsampling, x encoder represents the feature map of the corresponding layer in the encoder, x out represents the output graph of the corresponding layer in the encoder; x SE (c) represents the c-th channel of the feature map after SE module processing, x(c) represents the c-th channel of the input feature map, σ represents the Sigmoid activation function, W2 represents the weight matrix of the fully connected layer, which is used to map the features generated by W1 and GAP operation to the final channel weight, W1 represents the weight matrix of another fully connected layer, which is used to convert the feature map obtained by GAP operation into an intermediate representation, and GAP(x) represents global average pooling; (6) Apparent velocity physical constraint layer; The apparent velocity physical constraint layer converts the prediction noise and label noise into the frequency-wavenumber domain, calculates the difference between the two and its corresponding apparent velocity error, and improves the accuracy of the model output with physical constraints. First, the input data is transformed by FK to obtain the frequency f and wave number k x The spectrum value spec of the domain, and then calculate the apparent velocity v, the expression is as follows: By converting the input prediction noise and label noise data from the time-space domain to the frequency-wavenumber domain, and then calculating the spectrum and apparent velocity errors, the loss is composed of the spectrum difference between the prediction noise and label noise in the FK domain and the apparent velocity difference between the calculated prediction value and the label value. The expression is as follows: Among them, A represents the number of samples, B represents the batch size, C represents the number of channels, and H represents the i , W i represents the frequency dimension and wave number dimension in the FK domain, represents the spectrum value of the predicted noise, represents the spectrum value of label noise, f h represents the hth frequency component at time; Apparent velocity is the sum of frequency f and wave number k x The speed loss is regarded as the mean square error between the predicted value and the label value, and the expression is as follows: in, represents the apparent velocity of the prediction noise, represents the apparent velocity of label noise, when k x , w =0, the corresponding position can be ignored to avoid the error of dividing by 0; Frequency-wavenumber domain loss Total FK-loss By FK loss and Velocity Loss The two parts are weighted and the expression is as follows: Total FK-loss =Velocity Loss +FK loss (25) (7) Output layer; The output layer is used to map the last feature map of the network to the required output space; the output layer converts the multi-channel feature map into the noise of the target output through a 1x1 convolution. The input feature map F is converted as follows: The size of the feature map F is B×C×H×W, where H represents the height of the data and W represents the width of the data. conv represents the convolution kernel weight, b conv represents the bias of the convolutional layer, Represents the predicted output; Finally, the predicted noise signal is converted into two-dimensional data, and the predicted noise signal is output and physical rationality constraints are imposed in the frequency domain through the apparent velocity physical constraint layer. At this time, the data distribution range is the normalized [-1,1] range, and then the predicted output data is restored to its original amplitude range, that is, the original amplitude of the data is restored by multiplying it by the maximum amplitude T used during normalization. The expression is as follows: D=D′×T (27) Where D′ represents the normalized prediction output, T represents the absolute maximum amplitude of the corresponding label data during the normalization process, and D represents the data after denormalization, that is, the target prediction output is restored to the original amplitude range; (8) Loss function design; Compare the predicted noise signal with the real signal and comprehensively design MSE, MAE and R 2 The model is optimized by using three loss functions, and the weights of each loss are defined: Weight_mse, Weight_mae, Weight_r2 to control their contribution to the total loss. The expression is as follows: in, Represents the predicted value, y i represents the true value, n represents the number of samples; NR and NT represent the number of detectors and the number of signal samples, respectively; m represents the mth sampling point in the time series, E[S(x j , t m ) represents the average value of all data in a single shot collection record, S(x j , t m ) represents the label signal, that is, at position x j and time t m The true value of the seismic signal on pre Represents the prediction signal obtained by the model. The value is close to 0. pre and S represent the prediction and target shot gather records; R 2 Reflects the goodness of fit of the prediction, with a value range of (-∞, 1]; Set the time and space domain Total timeloss Total spectrum - wave number domain FK-loss Two parts of loss, according to the dynamically adjusted weight ratio FK weight Weighted addition to form the total loss Total loss , the expression is as follows: Total loss =Total timeloss +FK weight *Total FK-loss (31)。

Images