A physical information circulation network W-VSP velocity inversion method based on imaging domain depth feature constraint

Abstract

The application discloses a physical information cycle network W-VSP velocity inversion method based on an imaging domain depth feature constraint, is applied to the field of seismic data processing and geophysical exploration, and aims at the non-convexity of a target function and extreme sensitivity of data quality of traditional data domain full waveform inversion; first, an RTM imaging network of PIRNN is established, a physical process of sound wave propagation is embedded into a recurrent neural network, forward propagation is realized to obtain simulation records by wave equation forward modeling, and a receiver-depth imaging domain gather is obtained by cross-correlation imaging conditions, and space-time wave field gradient information is recorded; then, a U-Net network is used for trace gather feature identification in the RTM prestack imaging domain, and a curved event skeleton is automatically picked up; then, a hybrid loss function fusing data domain residual and imaging domain flatness is constructed to realize representation of the trace gather feature; finally, space-time field gradient information of wave field propagation recorded by the PIRNN is used for backward propagation, and an adaptive gradient fusion algorithm is introduced to balance macroscopic constraints in the imaging domain and data domain detail updating, so that velocity inversion is realized.

Description

Technical Field

[0001] This invention belongs to the field of seismic data processing and geophysical exploration, and specifically relates to a high-precision noise-resistant velocity model inversion technology. Background Technology

[0002] High-precision subsurface velocity models are the cornerstone of successful seismic imaging, reservoir characterization, and oil and gas exploration and development. Among various seismic acquisition systems, the Walkaway-VSP (W-VSP) has a unique geometric configuration with the source at the surface and the geophone downhole. This configuration effectively avoids complex near-surface interference fields and theoretically has a higher signal-to-noise ratio (SNR) and higher resolution tectonic illumination than surface seismic data. To fully utilize the potential of W-VSP data, Full Waveform Inversion (FWI), a data fitting technique that utilizes full wavefield information (amplitude, phase, and travel time), is widely used to construct high-fidelity velocity models. However, traditional FWI faces significant challenges when processing real-world W-VSP data with strong noise.

[0003] The core challenge of traditional data-domain FWI lies in the non-convexity of the objective function and its extreme sensitivity to data quality. Traditional FWI is typically constructed in the data domain, using the L2 norm as the objective function. This method is highly susceptible to getting trapped in local minima, leading to the "cycle-skipping" phenomenon. More seriously, in actual W-VSP acquisition, data is often distorted by strong random noise, pipe waves, and environmental interference. Under low signal-to-noise ratio conditions, L2 norm-based inversion will forcibly fit these non-physical noise components, resulting in severe distortion of gradient calculations, ultimately causing inversion divergence or producing high wavenumber artifacts.

[0004] To overcome the limitations of data-domain FWI, the academic community has proposed imaging-domain methods (such as Wave-equation migration velocity analysis, WEMVA). These methods aim to evaluate velocity errors by assessing the focusing quality of the migration imaging or the flatness of the reflection phase axis. However, existing imaging-domain methods heavily rely on the accuracy of feature extraction. In low signal-to-noise ratio environments, traditional methods based on similarity or curvature calculation operators often fail to identify continuous valid phase axes from a noisy background, leading to the failure of the constructed objective function.

[0005] Meanwhile, while deep learning technology has shown potential in physics inversion, purely data-driven methods lack physical constraints and have poor generalization ability. There is an urgent need for a W-VSP velocity modeling method that can maintain inversion stability in noisy environments and effectively extract physical features to constrain the inversion process. Summary of the Invention

[0006] To address the aforementioned technical problems, this invention proposes a physical information cyclic network (W-VSP) velocity inversion method based on imaging domain depth feature constraints. This method is applicable to the RZ imaging domain feature inversion strategy of W-VSP observation systems and can effectively avoid the defect of traditional data domain FWI being overly sensitive to noise.

[0007] The technical solution adopted in this invention is: a W-VSP velocity inversion method based on imaging domain depth feature constraints, comprising:

[0008] S1. Based on the current velocity model, define the underground grid space and absorbing boundary conditions, then perform time-stepping simulation, inject the source wavelet, and use the acoustic wave equation to perform forward modeling to obtain synthetic data records.

[0009] S2. Obtain the waveform residual by calculating the difference between the synthetic data and the observed data;

[0010] S3. Generating the source wavefield based on the forward modeling of the source wavelet. The backpropagating wave field, obtained by combining waveform residuals or observation data from the detector location back underground, is then used to obtain the backpropagating wave field. , position Image value at Represented as source wave field The detector wave field propagating in the opposite direction The cross-correlation; by superimposing the imaging results of all detectors to form a 3D imaging data volume, and at the horizontal X coordinate position of the VSP well shaft. Slices were extracted from the area to obtain a two-dimensional depth-domain imaging gather that reflects the detector-depth relationship. ;

[0011] S4. Single-channel depth domain imaging gather The probability map based on the output of the pre-trained U-Net network. The algorithm performs a column-by-column search for the probability maxima and applies three-point parabolic interpolation for sub-pixel thinning to obtain continuous floating-point depth coordinates. ;

[0012] S5, Based on continuous floating-point depth coordinates Obtain the depth trajectory of the in-phase axis Based on in-phase axis depth trajectory Construct an objective function for imaging domain flatness that combines first-order slope and second-order curvature;

[0013] S6. Based on the objective function of imaging domain flatness, the fused total gradient of velocity model iterative update is obtained, and finally a high-precision velocity model is obtained.

[0014] The beneficial effects of this invention are as follows: The intelligent feature picking mechanism based on U-Net designed in this invention achieves robust estimation of the curvature of the in-phase axis in low signal-to-noise ratio RTM images, solving the problem of failure of the traditional WEMVA method under strong noise; the adaptive gradient fusion strategy effectively balances macroscopic velocity updates and detailed characterization. Compared with the traditional FWI, this invention significantly improves the robustness and convergence stability of the inversion, and can successfully reconstruct high-precision underground velocity structures under strong interference environments. Attached Figure Description

[0015] Figure 1 A schematic diagram illustrating the principle of constructing common imaging point gathers in the depth domain under different velocity errors;

[0016] Among them, (a) is the straight shape of the depth domain common imaging point gather under no error, (b) is the upward curve shape of the depth domain common imaging point gather under 0.6 times error, and (c) is the downward curve shape of the depth domain common imaging point gather under 1.2 times error.

[0017] Figure 2 This is a diagram illustrating the hybrid inversion network architecture and core flowchart of the present invention.

[0018] Figure 3 The diagram shows the feature extraction and inference results of the U-Net network under various typical noise and speed error conditions.

[0019] Among them, (a) is the feature extraction result of the common imaging point gather in the depth domain under the condition of no velocity error, (b) is the feature extraction result of the common imaging point gather in the depth domain under the condition of 1.2 times velocity error, (c) is the feature extraction result of the common imaging point gather in the depth domain under the condition of 0.8 times velocity error, (d) is the feature extraction result of the common imaging point gather in the depth domain under the condition of strong noise and no velocity error, (e) is the feature extraction result of the common imaging point gather in the depth domain under the condition of strong noise and 1.2 times velocity error, and (f) is the feature extraction result of the common imaging point gather in the depth domain under the condition of strong noise and 0.8 times velocity error.

[0020] Figure 4 A schematic diagram of the horizontal layered velocity real model and initial model for noise-free feasibility testing;

[0021] Wherein, (a) is the VSP observation model, and (b) is the initial model;

[0022] Figure 5A comparison chart of data domain loss and mixed loss inversion results under noise-free environment;

[0023] Among them, (a) is the inversion result based on traditional data domain loss, and (b) is the inversion result based on imaging domain smoothness loss of the present invention;

[0024] Figure 6 Comparison of depth domain common imaging point gathers before and after noise-free inversion;

[0025] Wherein, (a) is the depth domain common imaging point gather before inversion, and (b) is the depth domain common imaging point gather after inversion;

[0026] Figure 7 This is a schematic diagram of a horizontal layered model under strong noise testing (SNR=5dB);

[0027] Wherein, (a) is the VSP observation model under strong noise test, and (b) is the initial model under strong noise test;

[0028] Figure 8 A comparison of the final inversion velocity models of different methods under strong noise conditions;

[0029] Among them, (a) is the inversion result under a noise-free environment (as a benchmark reference), (b) is the inversion result based on traditional data domain loss under a strong noise environment, and (c) is the inversion result based on the hybrid loss of the present invention under a strong noise environment.

[0030] Figure 9 Comparison of depth domain imaging gathers before and after inversion under strong noise conditions;

[0031] Wherein, (a) is the depth domain imaging gather before inversion under strong noise environment, and (b) is the depth domain imaging gather after inversion under strong noise environment;

[0032] Figure 10 A comparison of the final imaging results of each method under strong noise conditions;

[0033] Among them, (a) is the final imaging result under noise-free environment (as a reference), (b) is the final imaging result under strong noise environment based on traditional data domain loss, and (c) is the final imaging result under strong noise environment based on the hybrid loss of the present invention.

[0034] Figure 11 This is a schematic diagram showing the actual velocity and initial velocity of the inclined layered model.

[0035] Where (a) is the actual velocity of the inclined layered model, and (b) is the initial velocity of the inclined layered model;

[0036] Figure 12Comparison of velocity inversion results for the inclined layered model under strong noise environment;

[0037] Among them, (a) is the inversion result under a noise-free environment (as a benchmark reference), (b) is the inversion result based on traditional data domain loss under a strong noise environment, and (c) is the inversion result based on the hybrid loss of the present invention under a strong noise environment.

[0038] Figure 13 Comparison of depth domain imaging gathers for tilted layered models under high-noise environments;

[0039] Wherein, (a) is the initial depth domain common imaging point gather before inversion, and (b) is the depth domain common imaging point gather after inversion;

[0040] Figure 14 Comparison of the final imaging results of the tilted layered model under strong noise environment;

[0041] Among them, (a) is the final imaging result under noise-free environment (as a reference), (b) is the final imaging result under strong noise environment based on traditional data domain loss, and (c) is the final imaging result under strong noise environment based on the hybrid loss of the present invention. Detailed Implementation

[0042] To facilitate understanding of the technical content of this invention by those skilled in the art, the following description, in conjunction with the accompanying drawings, further illustrates the invention.

[0043] like Figure 2 As shown, the present invention provides a physical information recurrent network (W-VSP) velocity inversion method based on imaging domain depth feature constraints. First, the subsurface velocity model to be inverted is set as... The observed W-VSP seismic record is In the first The next iteration mainly includes the following steps:

[0044] Step 1: Forward Modeling and Data Domain Residual Calculation

[0045] Based on the current speed model A physical information circulation network framework is established. First, a perfectly matched layer (PML) is set up, defining the subsurface grid space and absorbing boundary conditions. The current velocity model is mapped into the grid space as medium parameters to eliminate artificial boundary reflections (e.g., at the edge of a finite computational grid, wavefield propagation will produce non-physical spurious echoes; absorbing boundary conditions act as an energy-absorbing layer at the boundary, allowing the wavefield to smoothly pass through without reflection). Then, time-stepping simulation is performed, injecting the source wavelet into the source space where the wave equation is located. By using the acoustic wave equation for forward modeling, the current velocity model can be obtained. Controlled synthetic data recording .

[0046] The forward modeling process of the acoustic wave equation satisfies:

[0047]

[0048] in, The horizontal position coordinates are For depth coordinates, For time; For sound wave field, For the speed of sound waves, This is the source term (i.e., the source wavelet); Let be the Laplace operator, and let represent the second derivative of the wave field in space. This is the source term (source wavelet) injected into the wave equation, and This refers to the source wave field generated by forward modeling driven by the velocity model for this source term.

[0049] By calculating synthetic data With observation data The difference is used to obtain the waveform residual. In the formula, Represents seismic wavefield data. The composite data obtained from forward simulation. This is actual observation data; Indicates the difference. This represents the waveform residual between the two;

[0050] Then, the L2 loss term in the data domain is calculated using the waveform residuals. :

[0051] , It represents the square of the L2 norm, used to measure the overall energy of the waveform residual.

[0052] Step 2: Construct depth domain imaging gathers

[0053] like Figure 2 The system flowchart shown uses reverse time migration (RTM) imaging conditions (where observation data is used as boundary conditions in RTM, extrapolating time backward from the detector location towards the subsurface to construct the detector wavefield) to establish the relationship between observation data and migration images. This specifically includes the following two steps:

[0054] (1) Cross-correlation: position Image value at It can be represented as the source wave field. The detector wave field propagating in the opposite direction The cross-correlation. Among them, the source wave field. It is generated by the forward modeling of the source wavelet, while the wave field propagates in the reverse direction. It is obtained by transmitting waveform residual data or observation data, which contains waveform information, back to the underground from the detector location.

[0055]

[0056] in, The coordinates or index number of the detector's location. This represents the total recording time.

[0057] (2) Superposition: The imaging results of all detectors are superimposed to form a 3D imaging data volume, and the data is placed in a horizontal position. A slice is extracted at the horizontal X-coordinate of the VSP wellbore to obtain a two-dimensional depth-domain imaging gather that reflects the detector-depth relationship. .

[0058] Step 3: In-phase axis feature picking based on U-Net

[0059] Single-channel depth domain imaging gather The input is into the pre-trained U-Net network. The U-Net used in this invention is based on a typical convolutional encoder-decoder architecture (reference: Ronneberger et al., 2015, "U-Net: Convolutional networks for biomedical image segmentation").

[0060] Pre-training setup: The network is trained using 2000 random layered velocity models generated by forward modeling of the wave equation, supplemented with a depth-domain imaging gather dataset constructed using random noise of varying intensities. The training objective is to learn the physical distribution of the phase axis. Training labels are provided by continuous pixel paths (i.e., the true depth coordinate sequence of the phase axis skeleton, serving as high-quality pseudo-labels to supervise network learning) automatically tracked by a dynamic programming algorithm. Binary cross-entropy (BCE) is used as the loss function until the network converges and meets the stopping condition.

[0061] Network inference: Probability graph of network output after Sigmoid activation This represents the probability that each pixel belongs to the in-phase axis skeleton. Subsequently, the maximum probability is searched column by column, and sub-pixel thinning is performed using three-point parabolic interpolation to obtain continuous floating-point depth coordinates. This means that the detector is in a horizontal position. At that location, the specific depth value corresponding to the phase axis. .

[0062] Step 4: Construct the hybrid objective function

[0063] Based on the in-phase axis depth trajectory extracted from the network (i.e., coordinates) Construct an objective function for imaging domain smoothness that combines first-order slope and second-order curvature:

[0064] First-order slope loss (calculated by the depth difference between adjacent traces, used to eliminate overall slope error and restore macroscopic background velocity):

[0065]

[0066] Second-order curvature loss (measures the rate of change of slope, used to eliminate higher-order residual time difference and curvature caused by high-frequency noise):

[0067]

[0068] in This represents the total number of downhole detectors. , , They represent the first time. , No. and the The floating-point depth picked up by the network on each detector channel.

[0069] The total loss of imaging domain smoothness is the weighted sum of the two:

[0070]

[0071] In the formula, and The pre-set empirical weighting coefficients control the proportions of slope and curvature in the flatness constraint, respectively.

[0072] Final hybrid objective function Defined as the union of the data domain and the imaging domain:

[0073]

[0074] in, This represents the underground velocity model to be inverted. The core of the optimization in the objective function is to find an optimal velocity model. To minimize the error. During the algorithm iteration process, the first... The model for the next iteration is denoted as . This represents the speed of sound propagation in the sound wave equation.

[0075] Step 5: Gradient fusion and velocity model update

[0076] make Let the gradient of the objective function in the data domain with respect to the velocity model be... This represents the gradient of the objective function in the imaging domain with respect to the velocity model. To prevent any one gradient from dominating the update direction, the weights are dynamically adjusted using the ratio of the norms of the two gradients. :

[0077]

[0078] in, This is a preset target scaling factor, which is empirically chosen based on the magnitude of the data gradient in practical applications. For example, 0.7 is used to align the initial magnitudes of the two gradients. To prevent division by zero of extremely small constants, for example .

[0079] The total gradient of the fusion was calculated. Then, the speed model is iteratively updated using the Adam optimization algorithm:

[0080]

[0081] in, and These represent the speed models before and after the update, respectively. For learning rate, The Adam optimization algorithm is based on the fused total gradient. The calculated update step size matrix.

[0082] The process continues until the convergence requirement or the maximum number of iterations is reached, at which point the final high-precision speed model is output.

[0083] Figure 1 A physical benchmark for the inversion task was established. By demonstrating the morphological changes (straight, upturned, downturned) of the common imaging point gathers in the depth domain under different velocity errors (no error, 0.6 times error, 1.2 times error), the mapping relationship between gather flatness and velocity accuracy was intuitively demonstrated, providing a physical basis for the subsequent establishment of the flatness objective function.

[0084] Figure 3 The reliability of this invention under extreme environments was verified. By demonstrating the inference results of U-Net under various typical noise levels (such as SNR=5dB) and different velocity deviations, it was proved that this invention can overcome the limitations of traditional mathematical operators, robustly extract the in-phase axis skeleton under strong interference, and ensure the accuracy of the source of the inversion gradient.

[0085] Figures 4 to 6This constitutes an ideal experiment in a noise-free environment. Figure 4 A benchmark model for testing is given. Figure 4 In the middle (a), the VSP observation model is shown. Figure 4 (b) is the initial model obtained by adding Gaussian blur to the observation model; Figure 5 and Figure 6 The comparison results prove that "imaging domain flatness" as an independent objective function still possesses the physical integrity to drive the velocity model to converge accurately without relying on waveform matching, thus verifying the effectiveness of the core constraints of the present invention.

[0086] Figures 7 to 10 A comparative experiment was conducted to address the common problem of low signal-to-noise ratio in actual production. Figure 7 A strong noise (SNR=5dB) experimental environment was set up; Figures 8 to 10 The comparison results significantly demonstrate the superiority of the proposed solution over the traditional L2-FWI: when the traditional method causes the model to diverge due to overfitting noise, the proposed solution can still reconstruct a clear layered structure and accurate background velocity, proving the robustness of the proposed solution under strong interference environment.

[0087] Figures 11 to 14 Extend the application scenarios to complex models with tilted interfaces. Figure 11 Increased construction complexity; Figures 12 to 14 The inversion results show that even under conditions of complex wavefield propagation paths and strong noise, the present invention can still accurately capture the kinematic characteristics of tilted phase axes through U-Net, and achieve high-precision reconstruction of complex underground structures, verifying the generalization application value of the present invention in actual complex geological tasks.

[0088] Figure 4 , Figure 7 , Figure 11 In this context, Receivers refers to the locations where receivers are deployed, and Sources refers to the locations where seismic sources are deployed.

[0089] Those skilled in the art will recognize that the embodiments described herein are for the purpose of helping to understand the principles of the invention, and should be understood that the scope of protection of the invention is not limited to such specific statements and embodiments. Various modifications and variations can be made to the invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the invention should be included within the scope of the claims of the invention.

Claims

1. A velocity inversion method for a physical information recurrent network (W-VSP) based on imaging domain depth feature constraints, characterized in that, include: S1. Define the underground grid space and absorbing boundary conditions, map the current velocity model to the grid space as medium parameters, then perform time-stepping simulation, inject the source wavelet at the source space location where the wave equation is located, and use the acoustic wave equation to perform forward modeling to obtain synthetic data. S2. Obtain the waveform residual by calculating the difference between the synthetic data and the observed data; S3. Generating the source wavefield based on the forward modeling of the source wavelet. The backpropagating wave field, obtained by combining waveform residuals or observation data from the detector location back underground, is then used to obtain the backpropagating wave field. , position Image value at Represented as source wave field The detector wave field propagating in the opposite direction The cross-correlation; by superimposing the imaging results of all detectors to form a 3D imaging data volume, and at the horizontal X coordinate position of the VSP well shaft. Slices were extracted from the area to obtain a two-dimensional depth-domain imaging gather that reflects the detector-depth relationship. ; S4. Single-channel depth domain imaging gather The probability map based on the output of the pre-trained U-Net network. The algorithm performs a column-by-column search for the probability maxima and applies three-point parabolic interpolation for sub-pixel thinning to obtain continuous floating-point depth coordinates. ; S5, Based on continuous floating-point depth coordinates Obtain the depth trajectory of the in-phase axis Based on in-phase axis depth trajectory Construct an objective function for imaging domain flatness that combines first-order slope and second-order curvature; S6. Based on the objective function of imaging domain flatness, the fused total gradient of velocity model iterative update is obtained, and finally a high-precision velocity model is obtained.

2. The W-VSP velocity inversion method based on imaging domain depth feature constraints for physical information recurrence networks according to claim 1, characterized in that, The U-Net in step S4 employs a typical convolutional encoder-decoder architecture.

3. The W-VSP velocity inversion method based on imaging domain depth feature constraints for physical information recurrence networks according to claim 2, characterized in that, The training dataset used in step S4 for pre-training is a dataset constructed from several random layered velocity models generated based on the forward modeling of the wave equation and deep domain imaging gathers constructed by adding random noise of different intensities.

4. The W-VSP velocity inversion method based on imaging domain depth feature constraints for physical information recurrence networks according to claim 3, characterized in that, The training objective during the pre-training process is to learn the physical distribution law of the in-phase axes. The training labels are provided by the continuous pixel paths automatically tracked by the dynamic programming algorithm, and the training loss function is the binary cross-entropy.

5. The W-VSP velocity inversion method based on imaging domain depth feature constraints for physical information recurrence networks according to claim 4, characterized in that, The objective function for imaging domain flatness in step S5 is expressed as: ； in, The objective function for imaging domain smoothness is represented by... This represents the L2 loss term in the data domain. This represents the total loss of smoothness in the imaging domain. , and These are empirical weighting coefficients. This represents the first-order slope loss. This represents the second-order curvature loss; Indicates the weight.

6. The W-VSP velocity inversion method based on imaging domain depth feature constraints for physical information recurrence networks according to claim 5, characterized in that, It is obtained by calculating using waveform residuals.

7. The W-VSP velocity inversion method based on imaging domain depth feature constraints for physical information recurrence networks according to claim 6, characterized in that, Also includes the Dynamic adjustments will be made, and the adjustment process is as follows: Calculate separately , The gradient of the velocity model; The weights are dynamically adjusted using the ratio of their gradient norms. .

Images