A frequency domain viscous acoustic full waveform inversion method based on dynamic Gaussian stable Q compensation

Abstract

The present application relates to the field of geophysical exploration and seismic wave imaging technology, specifically to a frequency domain viscous acoustic full waveform inversion method based on dynamic Gaussian stable Q compensation. It aims to solve the problems of numerical overflow, Gibbs pseudo-ghost and frequency band adaptive control failure caused by the mismatch between the high-frequency exponential divergence of the compensation operator and the dynamic evolution demand of the multi-scale frequency band in the existing time domain inversion. The method divides the progressive inversion frequency band, constructs a frequency domain decoupling double-branch complex slowness model containing forward and accompanying branches; a second-order Gaussian penalty stabilizing factor is introduced into the accompanying branch and multiplied by the compensation amplification term, and the safety frequency threshold is dynamically bound to the lower limit of the current frequency band to balance the high-frequency divergence; the forward branch is used to calculate the forward wave field and calculate the residual, which is used as a virtual source to be substituted into the accompanying branch to obtain the accompanying compensation wave field, and the imaging gradient is constructed accordingly to update the velocity model, and the iteration is performed frequency by frequency until convergence.

Description

Technical Field

[0001] This invention relates to the field of geophysical exploration and seismic wave imaging technology, and more specifically, to a frequency domain viscous acoustic wave full waveform inversion method based on dynamic Gaussian stable Q compensation. Background Technology

[0002] Currently, mainstream viscous acoustic wave full-waveform inversion attenuation compensation techniques mainly focus on the time domain. The conventional strategy is to introduce attenuation operators with opposite signs for the time-domain adjoint wave field extrapolation equations such as fractional-order Laplace's equations, in order to recover the high-frequency components lost when seismic waves pass through strongly attenuating rock media such as shallow gas clouds and the reflection amplitude from deep strata during the time-stepping inverse process.

[0003] However, existing technologies are limited by the time-domain computational framework, and wavefield attenuation compensation faces highly ill-posed mathematical and engineering bottlenecks, specifically manifested in:

[0004] Error accumulation and high-frequency numerical overflow in time-domain step-by-step solutions: The time-domain compensation operator is essentially an exponentially amplifying term that monotonically increases with extrapolation time. In finite-difference extrapolation of discrete time-domain grids, high-frequency dispersion noise accumulates and diverges exponentially with each time step, easily leading to severe high-frequency numerical overflow, thus disrupting the correct velocity update gradient.

[0005] The inherent drawbacks of time-domain stabilization strategies: To suppress time-domain divergence, existing schemes often employ low-pass hard truncation or first-order decay smoothing strategies. However, hard truncation disrupts the continuity of the wavenumber spectrum and easily induces severe numerical artifacts such as the Gibbs phenomenon during spatiotemporal extrapolation; while the decay rate of the first-order smoothing operator in the time domain is the same as the growth rate of the compensation term, making it difficult to achieve sufficient and effective numerical balance in the extreme high-frequency region.

[0006] The core contradiction of the time-domain framework lies in its lack of multi-scale band adaptive control capability: In wideband multi-scale full-waveform inversion iterations, the inversion process progresses gradually from low to high frequency bands. Because time-domain simulations must include wideband wavefield components in a single calculation, it is difficult to decouple and independently control a single frequency or specific band as in the frequency domain. This results in fixed time-domain stable filter parameters failing to adapt to bandwidth variations across different frequency bands, easily falling into the technical contradiction of excessive smoothing in narrow bands leading to resolution loss or insufficient suppression in wide bands causing iterative divergence; furthermore, time-domain wideband inversion is accompanied by extremely large computational resource consumption. Summary of the Invention

[0007] The purpose of this invention is to overcome the shortcomings of the prior art and provide a frequency domain full waveform inversion method and system based on dynamic Gaussian stability factor Q compensation. This solves the technical problems of numerical overflow, Gibbs artifacts and frequency band adaptive control failure caused by the exponential divergence characteristics of the accompanying compensation operator and the mismatch between fixed filter parameters and the evolution requirements of multi-scale progressive frequency bands in the existing time domain viscous acoustic full waveform inversion.

[0008] To achieve the above objectives, the present invention employs the following technical means:

[0009] This invention provides a frequency domain full waveform inversion method based on dynamic Gaussian stability factor Q compensation, comprising the following steps:

[0010] Step S1: Obtain the seismic wavefield data, initial velocity model, and quality factor model of the area to be inverted, and divide it into multiple inversion frequency bands from low to high according to the multi-scale inversion strategy;

[0011] Step S2: Within the current inversion frequency band, construct a frequency domain decoupled dual-branch complex slowness model, which includes a forward branch for simulating natural physical attenuation and an adjoint branch for performing attenuation compensation.

[0012] A second-order Gaussian penalty stabilization factor is introduced into the adjoint branch, and the second-order Gaussian penalty stabilization factor is multiplied by the compensation amplification term of the adjoint branch to counteract the high-frequency numerical divergence in the extrapolation of the adjoint wave field.

[0013] The safe frequency threshold of the second-order Gaussian penalty stability factor is dynamically bound to the lower limit of the current inversion frequency band;

[0014] Step S3: Solve the wave equation using the forward branch to obtain the propagating wave field and calculate the data residual. Substitute the data residual as a virtual source into the adjoint branch for back extrapolation to obtain the adjoint compensation wave field. Construct an imaging gradient based on the propagating wave field and the adjoint compensation wave field, and update the subsurface velocity model according to the imaging gradient. Iterate band by band until the inversion converges.

[0015] In the above scheme, step S1 includes:

[0016] Step S1.1: Receive the seismic wavefield data to be inverted, apply a spatial smoothing operator to the background velocity model to smooth it in order to construct a velocity background model that preserves the macroscopic structure, and convert the velocity background model into an initial quality factor model based on the rock physics empirical formula of quality factor and velocity.

[0017] Step S1.2: Based on the multi-scale inversion strategy, the seismic wavefield data is divided into multiple progressive inversion frequency bands from low to high with partial overlap between adjacent frequency bands. After the current inversion frequency band converges iteratively, the updated velocity model is used as the initial velocity model for the next inversion frequency band until the highest frequency band converges to obtain the final inversion result.

[0018] In the above scheme, the complex slowness model of the forward branch described in step S2 satisfies the following wavefield spatial evolution relationship:

[0019]

[0020] in, For a forward-modeling wave field, The initial amplitude, For transmission distance, For phase propagation terms, This is the amplitude attenuation term. ω is the angular frequency.

[0021] In the above scheme, the complex slowness model of the adjoint branch described in step S2 satisfies the following wavefield spatial evolution relationship:

[0022]

[0023] in, Accompanying compensation wave field; For amplitude compensation; the second-order Gaussian penalty stability factor It is introduced into this evolutionary relationship.

[0024] In the above scheme, the second-order Gaussian penalty stability factor Represented as:

[0025]

[0026] in, For frequency, The safe frequency is equal to the lower limit of the current inversion frequency band. is the Gaussian feature smoothing parameter.

[0027] In the above scheme, the construction of the imaging gradient in step S3 includes: performing conjugate cross-correlation operation on the forward propagation wavefield and the accompanying compensation wavefield to obtain the initial imaging gradient, and using the diagonal pseudo Hessian matrix to perform source and receiver amplitude modulation correction on the initial imaging gradient to complete the preconditioning process.

[0028] In the above scheme, the depth weighting process includes: performing depth weighting on the gradient after the preconditioning process using a depth weighting operator that grows nonlinearly with the subsurface depth, so as to forcibly compensate for the reflected illumination of the deep medium after penetrating the low-quality factor attenuation anomaly.

[0029] In the above scheme, the adaptive edge-preserving smoothing process includes: constructing the upper limit of the current multi-scale inversion frequency band. Equivalent spatial smoothing radius of dynamic binding :

[0030]

[0031] in, A minimum smoothing radius is set to prevent over-sharpening. To regulate the characteristic constant, The average wave velocity of the current stratum is given by dh, where dh is the spatial step size of the discrete grid; using the... A Gaussian spatial smoothing operator, which is the standard deviation of the filtering characteristics, is used to perform two-dimensional convolution filtering on the depth-weighted gradient so that the filtering effect adaptively converges as the inversion frequency band moves towards higher frequencies, thus preserving the boundaries of deep high-frequency reflective strata.

[0032] In the above scheme, the updated underground velocity model in step S3 adopts a quasi-Newton optimization algorithm, and iterates until convergence within the current inversion frequency band before entering the next frequency band.

[0033] This invention also provides a frequency domain full waveform inversion system based on dynamic Gaussian stability factor Q compensation, comprising:

[0034] The data preprocessing and multi-scale partitioning module is used to acquire seismic wavefield data, initial velocity model, and quality factor model of the area to be inverted, and to partition it into multiple inversion frequency bands from low to high.

[0035] The dual-branch wavefield simulation module is used to construct a frequency-domain decoupled dual-branch complex slowness model within the current inversion frequency band, and multiply the second-order Gaussian penalty stabilization factor with the compensation amplification term in the adjoint branch. The safe frequency threshold of the second-order Gaussian penalty stabilization factor is dynamically bound to the lower limit of the current inversion frequency band.

[0036] The gradient preconditioning and update module is used to construct the imaging gradient based on the forward propagation wavefield and the associated compensation wavefield, and update the subsurface velocity model according to the imaging gradient, iterating band by band until the inversion converges.

[0037] Because the present invention employs the above-mentioned technical means, it has the following beneficial effects:

[0038] 1. This invention derives a second-order Gaussian stable compensation operator for multi-scale inversion frequency band adaptation and substitutes it into the adjoint wave field inverse time extrapolation compensation mechanism (step S2). This overcomes the numerical instability problem caused by the exponential amplification of the backpropagated wave field energy in the full waveform inversion of viscous media, realizes the stable acquisition of gradient in the full waveform inversion of viscous media, and improves the accuracy of the full waveform inversion.

[0039] In the full waveform inversion of viscoelastic media, to recover the deep attenuation energy, it is necessary to use the adjoint state wave equation to perform inverse-time propagation compensation for absorption attenuation. However, the traditional compensation operator grows exponentially monotonically with propagation distance and frequency. If this is not constrained, the high wavenumber noise in the discrete calculation will be amplified infinitely, causing the back-extended wavefield energy to diverge and triggering numerical overflow. If conventional low-pass hard-truncation filtering is used, it will cause discontinuous abrupt changes in the wavenumber spectrum, inducing non-physical numerical oscillations (Gibbs effect) during the extrapolation process, thereby introducing false artifacts in the imaging gradient. This invention introduces a second-order Gaussian stabilizing operator into the adjoint state wave equation, so that amplitude compensation can pass without loss within the effective working frequency band, accurately recovering the dynamic characteristics of the formation. In the numerical dispersion band exceeding the grid resolution limit, the Gaussian stabilizing factor provides rapid attenuation control, making the extrapolated energy smoothly and gradually return to zero. This mechanism, which embeds spatial constraints into the wave equation, not only cuts off the chain of numerical divergence accumulation at the source, but also eliminates Gibbs numerical oscillations caused by hard truncation of filtering, providing a solid foundation for high-fidelity backpropagation wavefield reconstruction and stable gradient acquisition.

[0040] 2. This invention dynamically binds the safe frequency threshold of the viscoelastic Gaussian stability compensation operator to the lower limit of the current inversion working frequency band, and combines this with a multi-scale seismic frequency band update method that progresses from low frequency to high frequency and partially overlaps (steps S1 and S2). This method overcomes the defect of the conventional fixed Q-value compensation parameter and solves the technical problem that the absorption and dissipation law cannot adaptively perform multi-scale inversion due to changes in seismic wavelength; thus, it realizes the dynamic coordination and adaptive update of the frequency domain viscous acoustic wave compensation characteristics and the underground multi-scale velocity modeling process.

[0041] Multi-scale full-waveform inversion requires a frequency band-by-frequency progression from low to high frequencies to construct a subsurface velocity model. Using fixed, stable filtering parameters can lead to resolution loss due to over-smoothing in low-frequency band inversion, while insufficient suppression can cause iterative divergence in high-frequency band inversion. This invention dynamically binds a safety frequency threshold to the lower limit of the current inversion frequency band, making the passband boundary of the stability factor no longer static but automatically shifting and evolving as the inversion frequency band progresses. In low-frequency band inversion, a lower safety threshold allows for a wider safety band to construct a macroscopic geological framework; in high-frequency band inversion, the safety threshold automatically increases, precisely protecting the effective high-frequency signal at the current scale. This dynamic binding mechanism allows the compensation stability to evolve automatically with the inversion process, perfectly matching the dynamic evolution requirements of multi-scale frequency bands.

[0042] 3. In the frequency domain viscous acoustic wave control equation, this invention physically decouples the forward propagation physical branch of natural physical attenuation from the inverse propagation accompanying branch of high-fidelity energy recovery, and integrates the dynamic Gaussian penalty stabilization factor with the absorption compensation term in the wave field evolution (steps S2, S3 and corresponding system modules). This overcomes the bottleneck of mutual constraint between numerical divergence of accompanying wave field backpropagation and high-resolution reconstruction of fine reflection boundaries in the full waveform inversion of viscous media, and realizes high-definition characterization of complex geological structures below strong attenuation shielding areas under the premise of computational stability.

[0043] If the frequency domain is used alone to decouple the two-branch complex slowness model, although the forward attenuation and the associated compensation can be decoupled in terms of physical mechanism, the absorption compensation operator in the associated state wave equation is amplified exponentially with the propagation depth as a reverse time energy gain term. This causes the high-frequency discrete artifacts in the finite difference grid to expand sharply, eventually leading to serious numerical divergence and computational collapse. If only the second-order Gaussian stabilizing operator is used independently to perform conventional filtering post-processing on the imaging gradient or wave field data, the filtering is not deeply integrated into the wave equation of wave field evolution. It cannot block the spatial accumulation chain of computational overflow during backpropagation extrapolation, and its static cutoff frequency is completely unable to meet the requirements of multi-scale inversion mid-frequency band progressive advancement and gradual reconstruction of the ground velocity model skeleton.

[0044] By organically combining the two, the Gaussian stabilizing operator with an adaptive cutoff frequency is no longer a simple spatial filter, but rather a constraint on the evolution of physical properties, deeply integrated into the adjoint state-wave equation. This design introduces a dynamic physical game between "absorption compensation recovery" and "high-frequency spatial constraints" in each spatial propagation step of the adjoint wave field backpropagation: within the current inversion working frequency band, the Gaussian stabilizing factor maintains lossless passage, allowing the Q-value-compensated formation reflection amplitude characteristics to be reproduced with high fidelity; in the numerical dispersion band exceeding the grid resolution limit, the Gaussian stabilizing operator plays a rapid attenuation control role, smoothly and progressively reducing the energy of high-frequency numerical artifacts to zero, thus blocking the accumulation process of backpropagation calculation divergence at its source. More importantly, because the cutoff frequency of the stabilizing operator is dynamically correlated with the lower limit of the current inversion working frequency band, the adaptive control boundary of this "physical game" can dynamically shift from deep to shallow and from coarse to fine as the multi-scale inversion frequency band is gradually deduced from low to high. In low-frequency inversion, a wider safety frequency passband is allowed to reconstruct the macroscopic velocity skeleton; in high-frequency detail imaging, the passband cutoff surface automatically shifts to physically preserve high-resolution fine geological boundaries. This combination produces a synergistic effect of single-band adaptive convergence and stable multi-band connection, not only overcoming the inherent technical contradiction in imaging of strongly attenuated media where "high-resolution absorption compensation is inevitably accompanied by numerical divergence, while excessive smoothing results in the loss of fine geological boundaries," but also providing highly feasible engineering technical support for accurate imaging of deep thin reservoirs in complex high-attenuation exploration areas. Attached Figure Description

[0045] Figure 1 This is a diagram of the overall system architecture.

[0046] Figure 2 This is a comparison image of the inversion results;

[0047] The above figure shows the velocity reconstruction results after the highest frequency band iteration of the full waveform inversion. Figure c reveals the systematic suppression of deep illumination by physical attenuation: at depths shallower than 0.8 km, the strata achieve acceptable reconstruction quality; when the seismic wave passes through the central low-Q gas-bearing region, the high-frequency effective energy is significantly absorbed, the velocity gradient of the anticline limb is blurred, the velocity jump near the fault plane degenerates into a gradual transition, and dense non-physical water ripple-like oscillations appear in the shallow layer. Figure d shows the dual repair mechanism of the adaptive Gaussian stability factor Λ(f): within the effective frequency band, Λ(f)≈1, and the compensation gain is preserved; outside the upper limit of the frequency band, Gaussian quadratic attenuation suppresses the high-frequency components. Therefore, the velocity profile also shows a relatively clear anticline structure and fault boundary below the gas cloud, while the background area is smooth and natural. Detailed Implementation

[0048] The embodiments of the present invention will be described in detail below. Although the present invention will be described and illustrated in conjunction with some specific embodiments, it should be noted that the present invention is not limited to these embodiments. On the contrary, any modifications or equivalent substitutions made to the present invention should be covered within the scope of the claims of the present invention.

[0049] Furthermore, to better illustrate the present invention, numerous specific details are set forth in the following detailed embodiments. Those skilled in the art will understand that the present invention can be practiced without these specific details.

[0050] To address the technical challenges in existing time-domain viscous acoustic full-waveform inversion methods, such as the difficulty in simultaneously considering adaptive frequency band control, the tendency of attenuation compensation to cause high-frequency exponential numerical divergence and numerical oscillations (Gibbs effect), and the extreme imbalance of weak illumination gradients in deep strata, this invention breaks through the limitations of the time-domain framework and provides a frequency-domain full-waveform inversion method and system based on dynamic Gaussian stability factor Q compensation.

[0051] I. Operating Environment and Core Physical Parameter Configuration

[0052] 1. Hardware and software operating environment: Development and deployment are based on scientific numerical computing platforms such as MATLAB. This embodiment adopts a full waveform inversion framework, and the hardware supports multi-core CPU / GPU line computing to accelerate wavefield extrapolation.

[0053] 2. Model Construction Verification: A two-dimensional Marmousi gas cloud anomaly model was used for simulation verification.

[0054] Spatial grid settings: After thinning, the number of grid points is... (depth Distance), spatial grid step size m. An additional 40 layers of PM absorption boundary are added outside the effective boundary.

[0055] Initial and background parameters: A Gaussian smoothing operator is used to process the real velocity model to obtain a smooth initial velocity background; quality factor. The model generates a background skeleton based on empirical formulas from well logging petrophysics, constrains the water layer as a non-attenuation layer, and implants low-density sediments above the deep structures. A strongly attenuated gas cloud anomaly.

[0056] 3. Observation system and calculation parameters:

[0057] Seismic source setup: The excitation source uses the dominant frequency. Hz Ricker wavelet, reference frequency Hz.

[0058] Time-domain sampling: time sampling rate s, total number of sampling time steps

[0059] Algorithm optimizer: The L-BFGS algorithm is called during the model update phase, with a maximum number of iterations per frequency band set to 30 and a convergence tolerance of [value missing]. .

[0060] II. Multi-scale frequency band allocation strategy

[0061] Band 1: [1, 3] Hz

[0062] Band 2: [2, 6] Hz

[0063] Band 3: [5, 10] Hz

[0064] Band 4: [9, 15] Hz

[0065] Band 5: [13, 20] Hz

[0066] Band 6: [16, 24] Hz

[0067] The full waveform inversion iteration is performed independently within each frequency band, and the velocity model after convergence of the current frequency band is used as the initial model for the next frequency band.

[0068] III. Core Scheduling Pseudocode for Full Waveform Inversion

[0069] The core steps of the algorithm execution logic and dynamic compensation scheduling mechanism of the system's main control program are shown in the following pseudocode:

[0070] / / Start multi-scale inversion loop

[0071] FOR EACH band [f_low, f_high] IN Freq_Bands:

[0072] WHILE NOT Converged:

[0073] / / Step 1: Forward Evolution and Residual Calculation

[0074] Fwd_Slowness = Compute_Fwd_Branch(Current_Vp, Q_macro)

[0075] SimData, Pfwd = Solve_Helmholtz(Fwd_Slowness, Source)

[0076] Residual = SimData - ObsData

[0077] / / Step 2: Calculate the dynamic Gaussian stability factor

[0078] alpha = Calculate_Adaptive_Decay_Rate(f_low, f_high)

[0079] Lambda_f = Gaussian_Penalty(f, f_safe=f_low, alpha)

[0080] / / Step 3: Accompanied by stable compensation reverse extrapolation

[0081] Adj_Slowness = Compute_Adj_Branch(Current_Vp, Q_macro, Lambda_f)

[0082] Padj = Solve_Helmholtz(Adj_Slowness, Residual)

[0083] / / Step 4: Gradient Synthesis and Edge Preservation Update

[0084] Raw_Gradient = CrossCorrelation(Pfwd, Padj)

[0085] Gradient = Apply_Depth_Weight(Raw_Gradient, Z)

[0086] Gradient = Adaptive_Edge_Preserving_Filter(Gradient, f_high)

[0087] Current_Vp = L_BFGS_Update(Current_Vp, Gradient)

[0088] END WHILE

[0089] END FOR

[0090] Output High_Resolution_Vp_Model

[0091] The above embodiments systematically combine deep weighting and extreme high-frequency edge-preserving control, and under the action of the dual-branch complex slowness model, can accurately reconstruct shallow large-scale models and stably characterize complex high-precision reflection interfaces under deep gas clouds.

[0092] Example 1

[0093] The method flow of this embodiment includes the following steps:

[0094] Step S1: Data preparation and frequency band allocation

[0095] Seismic wavefield data, initial velocity model, and macroscopic quality factor model of the area to be inverted are obtained, and multiple inversion frequency bands are divided from low to high according to the above multi-scale inversion strategy.

[0096] Step S2: Construct a frequency-domain decoupled bi-branch complex slowness model

[0097] Within the current inversion frequency band, a frequency-domain decoupled bi-branch complex slowness model is constructed to address the different physical requirements of the forward evolution of the wave field and its accompanying reverse propagation.

[0098] Forward modeling branch construction: Construct a complex slowness model that includes a natural physical decay term. Separate real wavenumbers are obtained by performing a first-order decoupled expansion of the complex wavenumbers. With attenuation coefficient ,Right now:

[0099] .

[0100] Based on this, the spatial evolution law of the forward-modeling wave field can be derived:

[0101]

[0102] in, For a forward-modeling wave field, The initial amplitude, For transmission distance, For phase propagation terms, This is the amplitude attenuation term. ω is the angular frequency. This branch is used to accurately simulate the physical energy dissipation and dispersion drift of seismic waves propagating forward in the subsurface medium.

[0103] Accompanying branch construction: Invert the sign of the decaying imaginary part of the forward complex wavenumber, i.e. The evolution of the accompanying compensation wave field at this time is as follows:

[0104]

[0105] in, To accompany the compensation wave field, This is the amplitude compensation term. The formula shows that the amplitude of the accompanying wave field is exponentially amplified and recovered during extrapolation. To perform stable and reliable energy compensation, a second-order Gaussian penalty stability factor adaptively bound to the current inversion frequency band is introduced. .

[0106] To achieve safe extrapolation of the aforementioned adjoint branch, this paper presents a detailed and targeted design for the second-order Gaussian penalty stabilization factor. Theoretical analysis shows that the attenuation coefficient... Substituting into the adjoint evolution equation, the overall amplitude amplification factor of the adjoint wave field during extrapolation is... It will be controlled by a high-frequency growth term that is extremely dependent on frequency:

[0107]

[0108] This formula reveals the direct root cause of numerical instability: in actual discrete grid computing, when the wave field contains high-frequency components or noise, this linear exponential growth term will significantly amplify the proportion of high-frequency artifacts, causing exponential divergence.

[0109] Based on this, the present invention creatively introduces a second-order Gaussian decay balancing mechanism, which introduces a second-order Gaussian term:

[0110]

[0111] At this point, the total energy amplification factor of the high-frequency wave field It becomes:

[0112]

[0113] Due to the quadratic decay term in the exponent The divergence rate under limiting conditions is much faster than that of the linear compensating growth term. The speed allows the overall amplification factor to adaptively converge to zero at high frequencies.

[0114] To accommodate the progressive frequency bands of multi-scale full-waveform inversion, this scheme will include a safe frequency. Dynamic binding is used as the lower limit of the current inversion frequency band, establishing the following piecewise second-order Gaussian stability factor:

[0115]

[0116] exist area, This ensures that the effective Q-compensation information in the current inversion band can fully participate in gradient construction. In this region, the quadratic term provides an extremely smooth transition curve, protecting the effective signal bandwidth while forcefully cutting off numerical divergence beyond the physical limits of the grid. Since the full-waveform inversion gradient is typically constructed from the cross-correlation of the forward wavefield and the accompanying wavefield, high-frequency instabilities in the accompanying propagation are directly mapped to the gradient calculation process. Introducing the aforementioned stable Q-compensation factor effectively reduces the proportion of high-frequency artifacts in the gradient, thereby improving the stability of the inversion iteration and the imaging quality.

[0117] Step S3: Solve the wave equation using the forward branch to obtain the propagating wavefield and calculate the data residuals. Then, substitute the residuals as a virtual source into the adjoint branch for back extrapolation to obtain a stable adjoint compensated wavefield. Finally, perform conjugate cross-correlation on the propagating and inverse wavefields to obtain the imaging gradient. After passing through preconditioning based on the pseudo-Hessian matrix, depth weighting, and adaptive edge-preserving smoothing, update the subsurface velocity model using a quasi-Newton optimization algorithm, and iterate band-by-band until the inversion converges.

[0118] Key Algorithm / Data Processing Logic

[0119] (1) Core creative thinking analysis: Based on the need to combat high-frequency numerical overflow and the seismic wave dynamics mechanism of multi-scale bandwidth evolution, in order to solve the technical problem of attenuation compensation divergence, this invention creatively combines the dual-branch complex slowness decoupling model with the second-order dynamic Gaussian penalty factor, resulting in a synergistic technical effect of lossless amplitude preservation at low frequencies and effective suppression at high frequencies.

[0120] Core Step 1: Independently decouple the complex slowness square of the frequency domain Helmholtz equation and construct a forward branch containing natural physical attenuation and an adjoint branch containing a stability compensation factor.

[0121] Core Step 2: Process the adjoint branch compensation operator by multiplying the compensation exponential growth term with the second-order Gaussian decay term. Since the quadratic term diverges much faster than the first-order compensation term under high-frequency limiting conditions, a higher-order convergence balance is achieved mathematically to counteract the linear divergence term, effectively suppressing numerical explosion from a mathematical perspective.

[0122] (2) Core formulas and algorithm logic:

[0123] 1. Construction of the Two-Branch Complex Slowness Model: The forward branch complex slowness square is constructed as a decoupled product of the phase velocity term and the attenuation term, where the imaginary part of the attenuation term is negative. The adjoint branch complex slowness square is constructed as a decoupled product of the phase velocity term and the compensation term, where the imaginary part of the compensation term is positive, and a stability factor is introduced. The reference frequency is preferably set near the dominant frequency of the current source wavelet to minimize dispersive phase error.

[0124] 2. Dynamic Second-Order Gaussian Stabilization Factor and Adaptive Smoothing Parameter Algorithm: The mathematical form of the stability factor is as described in step S2. Within the safe frequency range, it is defined as a constant of 1; above the safe frequency range, it is defined as an exponential function exhibiting second-order Gaussian decay. To adapt to different bandwidths, the Gaussian characteristic smoothing parameter... An adaptive calculation algorithm is adopted. The calculation logic is determined by the extreme safety threshold and the width of the high-frequency buffer band. This logic ensures that when the maximum effective frequency exceeds a specific buffer threshold, the compensation energy is forcibly suppressed to a reasonable safety range, and has a minimum amplitude limiting mechanism to prevent high-frequency oscillation.

[0125] 3. Gradient depth weighting and scale-adaptive edge-preserving logic: The calculated initial imaging gradient matrix undergoes wavefield synthesis edge-preserving processing.

[0126] 3.1 Depth-weighted enhancement: Construct a weighted exponential matrix that grows non-linearly with depth, and forcibly compensate for the weak reflected illumination gradient in the deep layers through dot product operations.

[0127] 3.2. High-Frequency Edge Preservation Control: When performing Gaussian smoothing to eliminate numerical grid footprints, the smoothing radius is dynamically adjusted according to the upper limit of the current inversion frequency band. When entering the ultra-high frequency inversion stage and approaching the theoretical resolution limit, the algorithm modulus significantly reduces the smoothing operator to better preserve the high-frequency reflection boundaries of deep strata.

Claims

1. A frequency-domain full-waveform inversion method based on dynamic Gaussian stability factor Q compensation, characterized in that, Includes the following steps: Step S1: Obtain the seismic wavefield data, initial velocity model, and quality factor model of the area to be inverted, and divide it into multiple inversion frequency bands from low to high according to the multi-scale inversion strategy; Step S2: Within the current inversion frequency band, construct a frequency domain decoupled dual-branch complex slowness model, which includes a forward branch for simulating natural physical attenuation and an adjoint branch for performing attenuation compensation. A second-order Gaussian penalty stabilization factor is introduced into the adjoint branch, and the second-order Gaussian penalty stabilization factor is multiplied by the compensation amplification term of the adjoint branch to counteract the high-frequency numerical divergence in the extrapolation of the adjoint wave field. The safe frequency threshold of the second-order Gaussian penalty stability factor is dynamically bound to the lower limit of the current inversion frequency band; Step S3: Solve the wave equation using the forward branch to obtain the propagating wave field and calculate the data residual. Substitute the data residual as a virtual source into the adjoint branch for back extrapolation to obtain the adjoint compensation wave field. Construct an imaging gradient based on the propagating wave field and the adjoint compensation wave field, and update the subsurface velocity model according to the imaging gradient. Iterate band by band until the inversion converges.

2. The method according to claim 1, characterized in that, Step S1 includes: Step S1.1: Receive the seismic wavefield data to be inverted, apply a spatial smoothing operator to the background velocity model to smooth it in order to construct a velocity background model that preserves the macroscopic structure, and convert the velocity background model into an initial quality factor model based on the rock physics empirical formula of quality factor and velocity. Step S1.2: Based on the multi-scale inversion strategy, the seismic wavefield data is divided into multiple progressive inversion frequency bands from low to high with partial overlap between adjacent frequency bands. After the current inversion frequency band converges iteratively, the updated velocity model is used as the initial velocity model for the next inversion frequency band until the highest frequency band converges to obtain the final inversion result.

3. The method according to claim 1, characterized in that, The complex slowness model of the forward branch described in step S2 satisfies the following wavefield spatial evolution relationship: in, For a forward-modeling wave field, The initial amplitude, For transmission distance, For phase propagation terms, This is the amplitude attenuation term. ω is the angular frequency.

4. The method according to claim 1, characterized in that, The complex slowness model of the adjoint branch described in step S2 satisfies the following wavefield spatial evolution relationship: in, Accompanying compensation wave field; For amplitude compensation; the second-order Gaussian penalty stability factor It is introduced into this evolutionary relationship.

5. The method according to claim 4, characterized in that, The second-order Gaussian penalty stability factor Represented as: in, For frequency, The safe frequency is equal to the lower limit of the current inversion frequency band. is the Gaussian feature smoothing parameter.

6. The method according to claim 1, characterized in that, The construction of the imaging gradient in step S3 includes: performing conjugate cross-correlation operation on the forward propagation wavefield and the accompanying compensation wavefield to obtain the initial imaging gradient, and using the diagonal pseudo-Hessian matrix to perform source and receiver amplitude modulation correction on the initial imaging gradient to complete the preconditioning process.

7. The method according to claim 6, characterized in that, It also includes depth-weighted processing, which involves applying a depth-weighted operator that grows non-linearly with subsurface depth to the gradient after the preconditioning process, in order to forcibly compensate for the reflected illumination of the deep medium after penetrating the low-quality factor attenuation anomaly.

8. The method according to claim 6, characterized in that, The adaptive edge-preserving smoothing process includes: constructing the upper limit of the current multi-scale inversion frequency band. Equivalent spatial smoothing radius of dynamic binding : in, A minimum smoothing radius is set to prevent over-sharpening. To regulate the characteristic constant, The average wave velocity of the current stratum is given by dh, where dh is the spatial step size of the discrete grid; using the... A Gaussian spatial smoothing operator, which is the standard deviation of the filtering feature, is used to perform two-dimensional convolution filtering on the depth-weighted gradient so that the filtering effect adaptively converges as the inversion frequency band moves towards higher frequencies, thus preserving the boundary of the deep high-frequency reflection strata.

9. The method according to claim 1, characterized in that, The updated underground velocity model described in step S3 uses a quasi-Newton optimization algorithm and iterates within the current inversion frequency band until convergence before moving to the next frequency band.

10. A frequency-domain full-waveform inversion system based on dynamic Gaussian stability factor Q compensation, characterized in that, include: The data preprocessing and multi-scale partitioning module is used to acquire seismic wavefield data, initial velocity model, and quality factor model of the area to be inverted, and to partition it into multiple inversion frequency bands from low to high. The dual-branch wavefield simulation module is used to construct a frequency-domain decoupled dual-branch complex slowness model within the current inversion frequency band, and multiply the second-order Gaussian penalty stabilization factor with the compensation amplification term in the adjoint branch. The safe frequency threshold of the second-order Gaussian penalty stabilization factor is dynamically bound to the lower limit of the current inversion frequency band. The gradient preconditioning and update module is used to construct the imaging gradient based on the forward propagation wavefield and the associated compensation wavefield, and update the subsurface velocity model according to the imaging gradient, iterating band by band until the inversion converges.

Images