A VSP borehole wave suppression method based on improving U-Net with physical constraints

Abstract

The present invention discloses a VSP borehole wave suppression method based on improving U-Net with physical constraints. By introducing apparent velocity constraints in the F-K transform domain, physical constraint conditions are constructed. At the same time, combined with the GRU network, the temporal feature extraction and prediction of borehole waves are realized, and the mapping relationship between the noise image and the clean image is learned by the network. By training the network, the details and structure of the noise image are accurately restored, and the difference between the real denoised image and the network output is minimized to achieve the purpose of image denoising. The model retains the effective signal through the residual learning mechanism and simultaneously achieves precise suppression of borehole wave noise. The method of the present invention can accurately identify and effectively suppress borehole waves and their reflected waves in a complex VSP wavefield, significantly improve the signal-to-noise ratio, reduce the residual noise, improve the data processing efficiency while improving the data quality, provide technical support for high-precision seismic imaging under complex geological conditions, and provide a feasible solution for the popularization and application in the 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 uses sensors in deep underground wells to collect seismic signals, achieving 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 increasing attention. With the increasing complexity of exploration targets and the rising demand for wavefield processing accuracy and efficiency, how to effectively suppress 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, with 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 suppression of wellbore waves is often achieved 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 with relatively 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. The effective suppression of wellbore waves is still a difficult point in VSP signal fidelity processing and is a research direction 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 suppression of wellbore waves without manual parameter selection is achieved 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 wellbore 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 wellbore wave 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 dataset.

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

[0010] 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, obtaining the actual dataset, that is, more accurate wellbore wave noise data and effective signal data without wellbore waves.

[0011] Among them, the actual VSP data is the VSP wellbore wave noise data 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 perform slicing on the total dataset using the sliding slice method, that is, enhance the full-wavefield 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 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, 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 has been trained 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, temporal 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 acquired actual VSP data in real time to verify the model generalization.

[0018] Furthermore, step S3 is specifically as follows:

[0019] Fuse the synthetic dataset in step S1 and the actual dataset in step S2 to construct a total dataset, and set the VSP data in the total dataset to be a time series matrix D ∈ R 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 indices 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 dataset according to the sliding step sizes of S x = 16 and S y = 128. The overlap degree between each small slice is 50%, and borehole wave slice data with a size of 256×32 and original slice data containing borehole waves are obtained.

[0025] Furthermore, step S4 is specifically 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 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. It is assumed that the data set has L samples. The data set Z is divided 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, and split idx denotes the partitioning index. The expressions are 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. The DataLoader is used to load the dataset for training and validation, provided in batches of a fixed size Batch, and at the same time, the data loading order is shuffled. The batch data expressions are as follows:

[0047]

[0048] 0where 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, currently it is the q-th batch, represents the total number of batches and is rounded up.

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

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

[0051] Further, in the 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 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 global context information is extracted through global average pooling at each level of feature extraction, 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 changes in width and height are Height in represents the height of the input feature map, Weight in represents the width of the input feature map, Height out [[ID=3,2]]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 to 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, and 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, and 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, 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. The specific calculation expression is as follows:

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

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

[0082] x SE (c) = x(c).σ(W2 * ReLU(W1 * 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 being processed by the SE module, 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 the GAP operation to the final channel weights, W1 represents the weight matrix of another fully connected layer, which is used to transform the feature map obtained by the 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 converting 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 through the 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 converting 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 error of the apparent velocity, 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 calculation of 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 the 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 the label noise. When k x , wWhen it is equal to 0, the corresponding position can be ignored to avoid the error of division by zero.

[0095] Frequency-wavenumber domain loss Total FK-loss Through FK loss and Velocity Loss The two parts are weighted and calculated, and 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 required output space. The output layer 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:

[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, and 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 the physical rationality constraint is imposed in the frequency domain through the apparent velocity physical constraint layer. At this time, the data distribution range is in the range of [-1,1] after normalization. Then, the predicted output data is restored to its original amplitude range, that is, the original amplitude of the data is restored by multiplying 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] The predicted noise signal is compared with the true signal, and the MSE, MAE, and R 2 are comprehensively designed to optimize the model, and by defining the weights of each loss: Weight_mse, Weight_mae, Weight_r2, the contribution of each to the total loss is controlled. The expression is as follows:

[0106]

[0107]

[0108]

[0109] Among them, 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 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 parts of the loss in the time-space domain Total time loss and the frequency-spectrum 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, in combination with the GRU network, it realizes the extraction and prediction of the temporal characteristics of borehole waves, and utilizes 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, thereby achieving 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 VSP borehole wave suppression method based on a physically constrained improved U-Net according to 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 This 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 This is a schematic diagram of the predicted residual data, predicted wellbore wave data, predicted effective signal, and energy statistics of the U-Net method in the embodiment of the present invention.

[0125] Figure 13 This is a schematic diagram of the predicted residual data, predicted wellbore wave data, predicted effective signal, and energy statistics of the present method 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 through forward simulation by the finite difference method, then use a three-component geophone to record three-component gathers, generate pure noise-free VSP data, and then obtain a synthetic dataset by adding background noise of different intensities and synthetic wellbore wave data of 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 wellbore wave noise data and effective signal data without wellbore waves, and then use τ-P transform domain filtering to solve the dispersion phenomenon generated by F-K filtering to obtain an 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 and the validation set to evaluate the model. 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 f-k velocity physical 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, and establishes physical constraint conditions by introducing the f-k velocity constraint in the f-k transform domain, combining the differences in f-k velocity and the frequency spectrum in the f-k domain, and making it participate 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 acquired 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, with a horizontal coordinate of 3500 m and a depth range from 400 m to 3500 m, and the source is placed at the 5000 m position. A P-wave source with a dominant frequency of 30 Hz Ricker wavelet is used, and 161 geophones are arranged at 20 m intervals in the wellbore.

[0144] The elastic wave simulation is carried out in a two-dimensional underground medium to obtain the vibration components of the wave field along 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, and 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 is usually presented in the form of a spatio-temporal two-dimensional matrix. The time axis is called the long axis, and the space axis is the short axis, and there is a significant difference 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 indices 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 schematically as Figure 6 shown.

[0153] Among them, Figure 6 (a), (f) are respectively the sample and label schematic diagrams of slice data 1, Figure 6 (b), (g) are respectively the sample and label schematic diagrams of slice data 2, Figure 6 (c), (h) are respectively the sample and label schematic diagrams of slice data 3, Figure 6 (d), (i) are respectively the sample and label schematic diagrams of slice data 4, Figure 6 (e), (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 organize 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 to 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 the uneven gradient caused by the difference in dimension between different features, and then 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 formula (2), and perform normalization using formula (3). The expressions are as follows:

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

[0162]

[0163] where |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. Assuming 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 the evaluation results of all validation sets to obtain the overall performance of the model.

[0172] Assume X represents the total dataset after slicing and normalization, X train represents the training set, X val represents the validation set, split idx represents the division index, and 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 put in one-to-one correspondence 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 seismic trace number 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 putting the sample and the label in one-to-one correspondence, 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 borehole 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 at the same time. 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 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.

[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 borehole 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, it performs feature fusion with the corresponding layers of the encoder, and splices or adds the output feature map of the encoder with 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 through intermediate layer calculation, Indicates that the final output contains the original input and the output after nuclear transformation.

[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 weighted summation of the current input x t and the hidden state H at the previous time step t-1 and passed to the Sigmoid activation function σ for calculation. The expressions are as follows: z and W R The expressions are 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 expressions are as follows:

[0200]

[0201] The final hidden state H t is obtained by weighted summation of the candidate state and the previous hidden state H t-1 . The expressions are 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 respectively and reverse hidden states 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 ability to represent 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 [[ID=4I]]= Activation(Conv2D(x concat ) (19)

[0213] xSE (c) = x(c)·σ(W2 * ReLU(W1 * GAP(x))) (20)

[0214] 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 processed by the SE module, 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 operations to the final channel weights, W1 represents the weight matrix of another fully connected layer, which is usually used to transform the features obtained by the GAP operation (compressed global features) into intermediate representations, and 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 into 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 values 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 transforming 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 is composed 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 iDenote the frequency dimension and wavenumber dimension in the F-K domain, denote the spectral values of the predicted noise, denote the spectral values of the labeled noise, f h denote 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 (Velocity_Loss) is the mean square error between the predicted value and the labeled value, and the expression is as follows:

[0223]

[0224]

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

[0226] The loss in the frequency-wavenumber domain 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 conform to 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 target output noise 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, and the expression is as follows:

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

[0235] where 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) Loss function design;

[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. Using 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 labeled borehole wave. Through the F-K domain constraint, 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 the VSP data, the spatio-temporal domain can capture the signal changes in the time series and spatial distribution, so as to more accurately evaluate 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 the three parts of the loss function of MSE, MAE and R 2 are comprehensively considered to optimize the model. By defining the weights (Weight_mse, Weight_mae, Weight_r2) of each loss, their contributions to the total loss are controlled so that it can simultaneously restore the effective signal and suppress the noise on the spatio-temporal scale, thereby obtaining a more accurate denoising effect. Then 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, and handle data outliers. represents the predicted value, 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 by 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, get rid of traditional parameter dependencies through few-shot learning, fully explore the physical information of regular noise to construct constraint conditions, consider the Total in the time-space domain time loss and the Total in the frequency-spectrum wavenumber domain FK-loss two parts of losses, 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:

[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 frequency constraints; as the training progresses, the weight decreases, and the model pays more attention to the fitting of time domain features. Combining 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] The trained improved U-Net network model (target noise suppression and effective signal prediction model) obtained in step S5 is tested using newly acquired actual VSP data in real time to verify its generalization performance in practical applications. By suppressing the borehole waves in the actual noisy VSP data, it is evaluated whether the model can effectively adapt to the actual geological conditions and accurately suppress the borehole waves, thereby proving the reliability and practicability of the model. On the premise of minimizing the loss of effective signals, high-fidelity effective signals are obtained.

[0248] In this embodiment, a group 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 waves is about 1037.28. The effective signals can be obtained by subtracting the borehole wave noise from the original noisy wavefield, 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 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, respectively showing the wavefields 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 wavefields 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 waves, 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 waves.

[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 the grayscale images and energy statistics data of the residual results of Figure 12 and Figure 13 , it can be found that the residual of the predicted noise by the U-Net method is 82.083, while the residual of the predicted noise 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 comparison of the predicted noise results of the three methods reveals 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 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 kinematic feature prediction. 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 analyzing them in combination with evaluation indexes such as SSIM (structural similarity index), IOU (intersection over union), and R 2 and other evaluation indexes, 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 during the extraction of 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 uses an encoder-decoder structure combined with skip connections to extract rich features while retaining signal details. By performing apparent velocity constrained 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 judgment 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 compares the features of the borehole wave with the original noisy data through a residual learning mechanism, learns and extracts the features of the borehole wave, and realizes the efficient suppression of the borehole wave from the noisy data. Compared with the 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 the accurate suppression of the borehole wave. Compared with the existing methods, the method of the present invention 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 the prediction of the borehole wave and its reflected wave. The analysis of the prediction structure shows that the method of the present invention can better achieve the accurate suppression of the borehole wave and the restoration of the effective signal, providing a new method for the processing of the borehole wave and the 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 changes and modifications can be made to the present invention. Any modifications, equivalent replacements, improvements, 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 borehole wave suppression method based on improving U-Net with physical constraints, the specific steps are as follows: S1. Establish a velocity model, solve the elastic wave equation through forward simulation 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; 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; The pure and noise-free VSP data is added with background noise with 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; Among them, the relevant parameters of the borehole wave include: the amplitude, waveform and spectral characteristics of the borehole wave; 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; 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 traces and horizontal correction; 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-wavefield data through the sliding slice technique; S4. Classify and organize the total dataset after slicing in step S3, complete the construction of the sample and label datasets, and then perform normalization processing on the sample data and the label data respectively; 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 is trained, that is, the target noise suppression and effective signal prediction model; 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; The improved U-Net network model takes the original noisy data and borehole wave noise data after slicing and normalization as input, 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; S6. Based on the improved U-Net network model trained in step S5, use newly acquired actual VSP data in real time to verify the generalization of the model.

2. The VSP borehole wave suppression method based on improving U-Net with physical constraints according to claim 1, characterized in that The specific content of step S3 is as follows: Fuse the synthetic dataset in step S1 with the actual dataset in step S2 to construct a total dataset, and set the VSP data in the total dataset 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, where 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: P i，j = D[j:j + M, i:i + N] (1) where P i , j represents a small block extracted from matrix D, the position of which is determined by indices i and j, and a local area in matrix D has been extracted; D[j:j+M, i:i+N] represents a submatrix extracted from matrix D from row j to row j+M-1 and from column i to column i+N-1; Then, slice the total data set according to the sliding step sizes of S x = 16 and S y = 128, with an overlap of 50% between each small slice, to obtain wellbore wave slice data with a size of 256×32 and original slice data containing wellbore waves.

3. A VSP borehole wave suppression method based on improving U-Net by physical constraints according to claim 1, characterized in that The specific steps of step S4 are as follows: S41. Classify and sort the total data set after slicing in step S3 to complete the construction of the 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, normalize the maximum value of the data amplitude, that is, perform normalization processing on the sample data and the label data respectively; 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; The input data is \(D\in R\). M×N , calculate the absolute maximum value of the sample / tag amplitude according to Equation (2), and perform normalization using Equation (3). The expressions are as follows: T = max(|D i , l |) (2) where |D i，l | represents the absolute value of the amplitude in the absolute value wavefield of the element in the \(i\)-th row and \(l\)-th column of matrix \(D\), \(T\) represents its absolute maximum value. If \(T = 0\), the original wavefield data is returned. If \(T\neq0\), each data is normalized to obtain the normalized matrix \(D'\).

4. A VSP borehole wave suppression method based on improving U-Net with physical constraints according to claim 1, characterized in that The specific steps of step S5 are as follows: 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; 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; 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; it is assumed that the data set has L samples, and the data set Z is divided into K subsets, and the size of each subset is L / K; Among them, K = 5, and the division ratio of the training set to the validation set is 4:1; In the K-th verification, the subset Z is used K to represent the validation set, and the remaining Z\Z K represents the training set; then the data training set data validation set has the following mathematical expressions: Repeat the validation process in 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; Let \(X\) denote the total dataset after slicing and normalizing, \(X\) train denote the training set, \(X\) val denote the validation set, \(split\) idx denote the partitioning index, and 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; In the sliced and normalized total data set, the sample data and the label data are in one-to-one correspondence to form training pairs. The DataLoader is used to load the data set for training and validation, which is provided in batches of a fixed size Batch, and at the same time, the data loading order is shuffled; the expression of the batch data is as follows: 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, 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; The random order of the data is realized through the permutation index π, and the expression is as follows: Shuffle((X) = Permute(X, π) π represents the random permutation index (8).

5. A VSP borehole wave suppression method based on improving U-Net with physical constraints 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 sliced VSP original noisy data and borehole wave data. Using the sliced original noisy data as the sample and the borehole wave noise data as the label as the training input, when inputting, for each two-dimensional data sample, it is expanded into a four-dimensional structure, and the sample and the label are in one-to-one correspondence, and the data loading order is shuffled; Among them, the dimensions of the four-dimensional structure represent the batch size, the number of channels, the data height, and the 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 in each level of feature extraction, global context information is extracted through global average pooling, and the channel weights are adjusted; downsampling reduces the resolution of the feature map through operations while retaining the main features, and finally outputs 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: Among them, y(i, j) represents the output feature element after pooling at position (i, j), window(i, j) represents a 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, with the width and height changes 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; (3) Skip connection; 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 to retain the low-level feature information, and the expression is as follows: Among them, represents the input of the skip connection, represents the output calculated through the intermediate layer, represents that the final output contains the output after the original input is kernel-transformed; (4) GRU layer; In the 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 by weighted summation, via W z and W R passed to the Sigmoid activation function σ for calculation, and 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 represents the weight matrix for calculating the updated gate Z t , and W R represents the weight matrix for calculating the updated 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, and the closer it is to 0, the more attention is paid to the current input; R t controls the influence of the previous hidden state on the candidate state; Candidate hidden state Combined with the current input x t and the previous hidden state H t-1 Calculate, where S represents the weight matrix, and the expression is as follows: Final hidden state H t is obtained by weighted summation of the candidate state and the previous hidden state H t-1 as shown in the following expression: In a bidirectional GRU, the forward GRU and the reverse GRU process sequence information simultaneously, generating 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: (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, specifically including: upsampling, feature map concatenation and feature fusion; The decoder performs bilinear interpolation upsampling multiple times through Up sample. After each upsampling, the channels are reduced in dimension through the Conv2d and BatchNorm2d modules; combined with the corresponding feature map of the encoder, the feature fusion is realized through torch.cat; on the feature fusion result after each upsampling, a multi-resolution convolution module is applied to further extract and reconstruct the 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 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 being processed by the SE module, 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 for mapping the features generated by W1 and the GAP operation to the final channel weights, W1 represents the weight matrix of another fully connected layer for transforming the feature map obtained by the 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 improves the accuracy of the model output by physically constraining by converting the predicted noise and label noise into the frequency-wavenumber domain and calculating the difference between the two and its corresponding apparent velocity error; First, the input data is transformed by the F-K transform to obtain the spectral values spec in the frequency f and wavenumber k domain, and then the apparent velocity v is calculated. The expressions are as follows: x Domain spectral value spec, and then calculate the apparent velocity v. The expression is as follows: By converting 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 error of the apparent velocity. At this time, the loss is composed of the spectrum difference between the predicted noise and label noise in the F-K domain and the calculation of the apparent velocity difference between the predicted value and the label value, and the expression is as follows: 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 labeled noise, f h represents the frequency component at the h-th time; The apparent velocity is the ratio of the frequency f and the wave number k x The apparent velocity loss is the mean square error between the predicted value and the label value, and the expression is as follows: 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 division by zero; Total loss in frequency-wavenumber domain FK-loss Through FK loss and Velocity Loss The weighted calculation is performed in two parts, 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 target output noise through a 1x1 convolution. Then for the input feature map F, the conversion is as follows: Among them, 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; W conv represents the convolution kernel weight, and b conv represents the bias of the convolutional layer, represents the predicted output; 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], 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 by the maximum amplitude T used during normalization, and the expression is as follows: D = D' × T (27) 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; (8) Loss function design; Compare the predicted noise signal with the true signal, and comprehensively design the MSE, MAE and R 2 Three parts of the loss function to optimize the model, and by defining the weights of each loss: Weight_mse, Weight_mae, Weight_r2, to control their contribution to the total loss, the expression is as follows: Among them, represents the predicted value, 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 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 by the model, 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]; Set the time - space domain Total timeloss and the frequency - wavenumber domain Total FK-loss The two parts of the loss are weighted by the dynamically adjusted weight ratio FK weight and added together to form the total loss Total loss , and the expression is as follows: Total loss = Total timeloss + FK weight * Total FK-loss (31).

Images