Gaussian mixture clustering constrained full waveform inversion method, system, equipment and medium

Abstract

This invention provides a Gaussian mixture clustering-constrained full waveform inversion method, system, equipment, and medium, belonging to the field of geophysical exploration technology. The method includes: constructing an initial velocity model; obtaining an observation data matrix based on the true velocity model; performing reverse time migration imaging using the initial velocity model and the observation data matrix to obtain a reflection amplitude image; generating a subject mask based on the Gaussian mixture model and the reflection amplitude image; obtaining a boundary mask based on the subject mask; processing the subject mask and the boundary mask in conjunction with time-domain full waveform inversion to obtain a weighted gradient; and iteratively updating the initial velocity model using a multi-scale frequency strategy and optimization algorithm based on the weighted gradient until convergence conditions are met, outputting the final velocity model. This invention improves the imaging accuracy of complex geological boundaries while reducing the dependence on the initial model.

Description

Technical Field

[0001] This invention relates to the field of geophysical exploration technology, specifically to a Gaussian mixture clustering constrained full waveform inversion method, system, equipment, and medium. Background Technology

[0002] Seismic exploration is the only geophysical method that simultaneously possesses both great depth and high resolution, making it the preferred exploration technique in the field of deep resource exploration. The main function of seismic exploration is to identify subsurface geological structures by processing seismic signals acquired by detectors deployed on the surface. Therefore, obtaining the distribution of subsurface medium physical parameters from seismic data is particularly important, and among these parameters, the propagation velocity of seismic waves in different media is a key factor determining the quality of migration imaging. However, deep, special spaces (such as salt caverns, abandoned mines, and aquifers) typically possess complex geological structures and physical properties, making it difficult for conventional seismic exploration methods to accurately detect and characterize the fine geological information of these spaces. With the increasing utilization of subsurface space, these problems have become increasingly prominent, urgently requiring innovative technological breakthroughs. Therefore, in-depth research into high-precision velocity modeling of deep, special spaces has become a key technical issue in the current field of seismic exploration.

[0003] Currently, the widely used methods for subsurface medium velocity modeling both domestically and internationally mainly include three types: velocity analysis, tomography, and full waveform inversion. Among these, velocity analysis and tomography only utilize travel time information from seismic wave propagation, neglecting many other important information. Facing the challenges of complex geological structures and physical property distributions in deep, special spaces, especially in unique subsurface spaces such as salt caverns, abandoned mines, and aquifers, they cannot effectively identify subtle subsurface structures. Therefore, their velocity modeling results cannot meet the requirements of high-precision seismic exploration. Full waveform inversion (FWI) technology utilizes all information from the pre-stack seismic wavefield (including kinematic and dynamic information such as amplitude, phase, and travel time), achieving high-precision velocity modeling of subsurface media. Furthermore, it can combine this with the physical property distribution characteristics of deep, special spaces to provide accurate velocity distribution maps. These maps can reveal the velocity distribution of the subsurface medium, thereby inferring changes in the structure of deep, special spaces.

[0004] Full waveform inversion faces three major technical challenges in high-precision velocity modeling: First, it is highly dependent on the initial model. If the initial velocity model deviates significantly from the actual underground velocity distribution, the inversion is prone to getting trapped in local extrema, causing the model to converge to a non-physical solution. Second, it has weak resistance to noise and ambiguity. Actual seismic data inevitably contains acquisition noise and interference waves, and the heterogeneity of the underground medium will exacerbate the ambiguity of the inversion, leading to false structures or velocity anomalies in the inversion results. Third, it is not adaptable to complex structures. In areas with strong heterogeneity such as faults, salt domes, and volcanic rocks, underground velocities exhibit drastic lateral changes. Traditional inversion methods struggle to accurately capture velocity abrupt change boundaries, easily causing fuzziness or distortion of structural details.

[0005] To address the aforementioned issues, existing technologies primarily employ three types of solutions, but all have significant limitations. First, they rely on conventional time-domain full-waveform inversion without introducing additional constraints, making the results entirely dependent on the quality of the initial model. This can easily lead to convergence bias in low signal-to-noise ratio data or complex structural regions. Second, while simple smoothing constraints or spatial regularization methods can suppress noise to some extent, they excessively smooth out the actual changes in underground velocities, resulting in blurred structural boundaries (such as reservoir top and bottom plates, faults). Third, some methods attempt to introduce geostatistical constraints (such as variograms), but these require pre-setting geological prior parameters, have poor adaptability to unknown structural regions, and struggle to handle the complex distribution characteristics of non-stationary velocity fields. Summary of the Invention

[0006] In view of this, it is necessary to provide a Gaussian mixture clustering constrained full waveform inversion method, system, device and medium to solve the technical problems in the existing technology, such as high dependence on the initial model, weak noise resistance, difficulty in accurately characterizing complex structural boundaries, and easy occurrence of local extrema, false anomalies or blurred boundaries in the inversion results.

[0007] To address the aforementioned technical problems, in a first aspect, the present invention provides a Gaussian mixture clustering constrained full waveform inversion method, comprising:

[0008] Construct an initial velocity model; The observation data matrix is ​​obtained based on the true velocity model; The initial velocity model and the observation data matrix are used to perform reverse time migration imaging to obtain a reflection amplitude image; A subject mask is generated based on the Gaussian mixture model and the reflection amplitude image; The boundary mask is obtained based on the main body mask; The weighted gradient is obtained by combining the main body mask and the boundary mask with time-domain full waveform inversion. A multi-scale frequency strategy and optimization algorithm are adopted to iteratively update the initial velocity model based on the weighted gradient until the convergence condition is met, and the final velocity model is output.

[0009] In one possible implementation, constructing the initial velocity model includes: Obtain the true velocity model of the target area; The initial velocity model is obtained by performing Gaussian smoothing filtering on the actual velocity model.

[0010] In one possible implementation, obtaining the observation data matrix based on the true velocity model includes: Based on the real velocity model, the acoustic wave equation is solved, and forward modeling is performed to obtain the simulated seismic wave field. A set proportion of Gaussian white noise is added to the simulated seismic wavefield to simulate the actual acquisition environment, forming the observation data matrix.

[0011] In one possible implementation, generating the subject mask based on the Gaussian mixture model and the reflection amplitude image includes: Based on the Gaussian mixture model, the pixels of the reflection amplitude image are probabilistically clustered. The probability density function of the Gaussian mixture model is: ; in, I This represents the amplitude value of a pixel extracted from the reflection amplitude image. This indicates the preset number of cluster categories. Indicates category index, Denotes the mixing weights of the k-th Gaussian component and satisfies , This represents the mean of the k-th Gaussian component. Let the covariance matrix of the k-th Gaussian component be denoted as . The probability density function is... The optimal values ​​of the parameters are obtained by solving the probability density function, and the posterior probability of each pixel belonging to each category is calculated based on the optimal values; the parameters include , and ; Compare the posterior probabilities of each pixel belonging to each category, and assign each pixel to the category with the highest posterior probability. Based on the partitioning results, the mean values ​​of the Gaussian components corresponding to each category are compared, and the category with the largest mean value is determined as the construction region. All pixels in the reflection amplitude image that are assigned to the constructed region are marked with a first value, and pixels that are assigned to other categories are marked with a second value, thereby generating the binarized subject mask.

[0012] In one possible implementation, obtaining the boundary mask based on the subject mask includes: The main body mask is optimized. The optimized main body mask is calculated together with the original main body mask to obtain the boundary mask; the first value of the boundary mask represents the constructed boundary region and the second value represents the non-boundary region.

[0013] In one possible implementation, the step of obtaining the weighted gradient by combining the main body mask and the boundary mask with time-domain full waveform inversion includes: Construct a spatial weighting function based on the main body mask and the boundary mask: ; in, The spatial weighting function is... For the main body mask, For the boundary mask, The gain coefficient and ; The weighted gradient is calculated by combining the residual gradient of the time-domain full waveform inversion data obtained by the adjoint state method with the spatial weighting function.

[0014] In one possible implementation, the use of a multi-scale frequency strategy and optimization algorithm to iteratively update the initial velocity model based on the weighted gradient includes: The observation data matrices of multiple frequency bands are selected in order from low frequency to high frequency, and iterative optimization of time domain full waveform inversion is performed sequentially. In each round of iterative optimization, the velocity model is updated using the weighted gradient of the current iteration using the iterative optimization algorithm; the update formula is:

[0015] in, The location of the underground space at the nth iteration The current speed at that location, The location of the underground space at the (n+1)th iteration The latest speed at the location, As a weighted gradient, the gain coefficient decreases exponentially with the number of iterations.

[0016] Secondly, the present invention also provides a Gaussian mixture clustering constrained full waveform inversion system, comprising: The initial model building module is used to build the initial velocity model; The observation data matrix synthesis module is used to obtain the observation data matrix based on the real velocity model. The acquisition module is used to perform reverse time migration imaging using the initial velocity model and the observation data matrix to obtain a reflection amplitude image; The main body mask generation module is used to generate a main body mask based on the Gaussian mixture model and the reflection amplitude image; A boundary mask generation module is used to obtain a boundary mask based on the main body mask; The gradient weighting processing module is used to obtain a weighted gradient by combining the main body mask and the boundary mask with time-domain full waveform inversion. The iterative inversion module is used to iteratively update the initial velocity model based on the weighted gradient using a multi-scale frequency strategy and optimization algorithm until the convergence condition is met, and output the final velocity model.

[0017] Thirdly, the present invention also provides an electronic device, including a memory and a processor, wherein, The memory is used to store programs; The processor, coupled to the memory, is used to execute the program stored in the memory to implement the steps in the Gaussian mixture clustering constrained full waveform inversion method described in any of the above implementations.

[0018] Fourthly, the present invention also provides a computer-readable storage medium for storing a computer-readable program or instruction, which, when executed by a processor, can implement the steps in the Gaussian mixture clustering constrained full waveform inversion method described in any of the above implementations.

[0019] The beneficial effects of this invention are as follows: The Gaussian mixture clustering constrained full-waveform inversion method provided by this invention firstly uses a main mask generated by a Gaussian mixture model to effectively distinguish between real reflection signals and noise interference, concentrating the inversion update energy on effective structural regions, thereby more accurately restoring the main structure and macroscopic trends of the subsurface velocity model. Furthermore, by utilizing the boundary mask derived from the main mask, the sensitivity of the wavefield to stratigraphic interfaces is particularly enhanced during the inversion process, helping to sharpen the boundaries of the velocity model and obtain clearer and more accurate structural contours and geological details such as faults and pinch-outs, significantly improving the spatial resolution of the imaging. The intelligently generated mask distinguishes between effective signals and noise, focusing the inversion energy on key structural regions, significantly improving the accuracy of restoring the main structure and macroscopic trends of the velocity model; simultaneously, the boundary mask sharpens the stratigraphic interfaces, making the imaging of details such as faults and pinch-outs clearer. Furthermore, by processing the weighted gradient based on the main body mask and boundary mask in conjunction with time-domain full waveform inversion, a priori spatial constraint is introduced into the gradient direction. The mask-weighted gradient introduces reliable prior spatial constraints, effectively suppressing erroneous updates caused by noise and blind zones, significantly alleviating cycle jump problems, and making the inversion process more stable and convergent. Furthermore, the mask and multi-scale strategy work together to avoid redundant computation in invalid regions, accelerate global convergence, and protect the stability of high-frequency inversion. Furthermore, although an initial velocity model is still required, the mask-guided gradient update enhances tolerance and correction capabilities for erroneous model regions, reduces the stringent dependence on the accuracy of the initial model, and strengthens the ability to correct erroneous models. Simultaneously, the mask originates from data-driven imaging results, reducing subjective bias from human interpretation and improving automation and robustness. Attached Figure Description

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

[0021] Figure 1 This is a schematic flowchart of an embodiment of the Gaussian mixture clustering constrained full waveform inversion method provided by the present invention; Figure 2 For the present invention Figure 1 A schematic diagram of an embodiment of S101; Figure 3 For the present invention Figure 1 A schematic diagram of an embodiment of S102; Figure 4 For the present invention Figure 1A schematic diagram of an embodiment of S104; Figure 5 For the present invention Figure 1 A schematic diagram of an embodiment of S105; Figure 6 For the present invention Figure 1 A schematic diagram of an embodiment of S106; Figure 7 For the present invention Figure 1 A schematic diagram of an embodiment of S107; Figure 8 This is a schematic diagram of the reverse-time migration imaging result obtained based on the initial velocity model and the observation data matrix in one embodiment of the present invention; Figure 9 This is a schematic diagram of the category labels obtained after performing Gaussian mixture model clustering on a reverse time-shifted image in one embodiment of the present invention; Figure 10 This is a schematic diagram of a binary subject mask generated based on clustering results and a boundary mask extracted after morphological processing in one embodiment of the present invention; Figure 11 This is a schematic diagram comparing the final velocity model obtained by GMM-constrained full waveform inversion with the actual model in one embodiment of the present invention; Figure 12 This is a schematic diagram comparing the traditional full waveform inversion result with the actual velocity curve on a typical profile in one embodiment of the present invention; Figure 13 This is a schematic diagram of the changes in the loss function value and gradient norm with the number of iterations during the GMM-constrained full waveform inversion process in one embodiment of the present invention; Figure 14 This is a schematic diagram comparing the velocity curve of the GMM-constrained full waveform inversion result with the actual velocity curve at the middle horizontal position of the model in one embodiment of the present invention. Figure 15 This is a schematic diagram comparing the velocity curve at the middle depth position of the model with the actual velocity curve in one embodiment of the present invention, showing the velocity variation with horizontal position in the GMM-constrained full waveform inversion result. Detailed Implementation

[0022] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, and not all of them. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative effort are within the scope of protection of the present invention.

[0023] In the description of the embodiments of the present invention, unless otherwise stated, "multiple" means two or more. "And / or" describes the relationship between related objects, indicating that there can be three relationships. For example, A and / or B can represent three situations: A exists alone, A and B exist simultaneously, and B exists alone.

[0024] The terms "first," "second," etc., used in the embodiments of this invention are for descriptive purposes only and should not be construed as indicating or implying their relative importance or implicitly specifying the number of technical features indicated. Therefore, a technical feature defined with "first" or "second" may explicitly or implicitly include at least one of that feature.

[0025] In this document, the term "embodiment" means that a particular feature, structure, or characteristic described in connection with an embodiment may be included in at least one embodiment of the invention. The appearance of this phrase in various places throughout the specification does not necessarily refer to the same embodiment, nor is it a separate or alternative embodiment mutually exclusive with other embodiments. It will be explicitly and implicitly understood by those skilled in the art that the embodiments described herein can be combined with other embodiments.

[0026] This invention provides a Gaussian mixture clustering constrained full waveform inversion method, system, device and medium, which are described below.

[0027] Figure 1 This is a schematic flowchart of an embodiment of the Gaussian mixture clustering constrained full waveform inversion method provided by the present invention, as shown below. Figure 1 As shown, the Gaussian mixture clustering constrained full waveform inversion method includes: S101. Construct the initial velocity model; S102. Obtain the observation data matrix based on the real velocity model; S103. Using the initial velocity model and the observation data matrix, perform reverse time migration imaging to obtain a reflection amplitude image; S104. Generate a subject mask based on the Gaussian mixture model and the reflection amplitude image; S105. Obtain the boundary mask based on the main body mask; S106. Based on the main body mask and the boundary mask, the weighted gradient is obtained by combining time-domain full waveform inversion. S107. Using a multi-scale frequency strategy and optimization algorithm, the initial velocity model is iteratively updated based on the weighted gradient until the convergence condition is met, and the final velocity model is output.

[0028] It should be noted that: Based on the geological understanding of the work area, well logging data, or existing seismic wave velocity information, an initial velocity model reflecting the approximate distribution of formation velocities and tectonic trends is established. This model will serve as the starting point for subsequent inversion. Forward simulation is performed using known or hypothesized real velocity models to generate a seismic wavefield observation data matrix recorded by multiple shots and receivers, which is then organized into a matrix form for comparison and constraint in subsequent inversion. Using the reverse time migration method, the observation data matrix is ​​backpropagated and cross-correlated with the forward propagation wavefield. Based on the initial velocity model, an amplitude image reflecting the subsurface reflection interface is calculated, which contains tectonic information but also noise and artifacts. Gaussian mixture model clustering analysis is performed on the amplitude value distribution of the reflection amplitude image to distinguish between major reflecting structural regions and noise background. Based on this, a binarized subject mask is generated to identify the main areas where effective reflection signals are located in the image. Morphological boundary detection or gradient calculation is performed on the subject mask to extract the edge information of the effective areas, generating a boundary mask. This mask will be used to strengthen the constraint on tectonic boundaries during inversion. In the time-domain full-waveform inversion process, a main body mask is used to constrain the model update region, and a boundary mask is used to enhance the gradient at the interface. By spatially weighting the inversion gradient, the contribution of effective tectonic regions is highlighted and the interference of irrelevant regions is suppressed. Data from different frequency bands are used stepwise from low frequency to high frequency to participate in the inversion. Combined with optimization algorithms (such as L-BFGS, conjugate gradient method, etc.), the velocity model is repeatedly updated according to the weighted gradient until the data residual is small enough or the model change tends to stabilize, and finally a high-precision velocity model that conforms to the actual geological conditions is obtained.

[0029] In summary, the Gaussian mixture clustering-constrained full-waveform inversion method provided in this invention firstly uses a main mask generated by a Gaussian mixture model to effectively distinguish between real reflection signals and noise interference, concentrating the inversion update energy on effective structural regions, thereby more accurately restoring the main structure and macroscopic trends of the subsurface velocity model. Furthermore, by utilizing a boundary mask derived from the main mask, the sensitivity of the wavefield to stratigraphic interfaces is particularly enhanced during the inversion process, helping to sharpen the boundaries of the velocity model and obtain clearer and more accurate structural contours and geological details such as faults and pinch-outs, significantly improving the spatial resolution of the imaging. The intelligently generated mask distinguishes between effective signals and noise, focusing the inversion energy on key structural regions, significantly improving the accuracy of restoring the main structure and macroscopic trends of the velocity model; simultaneously, the boundary mask sharpens stratigraphic interfaces, making the imaging of details such as faults and pinch-outs clearer. Furthermore, by processing the weighted gradient based on the main body mask and boundary mask in conjunction with time-domain full waveform inversion, a priori spatial constraint is introduced into the gradient direction. The mask-weighted gradient introduces reliable prior spatial constraints, effectively suppressing erroneous updates caused by noise and blind zones, significantly alleviating cycle jump problems, and making the inversion process more stable and convergent. Furthermore, the mask and multi-scale strategy work together to avoid redundant computation in invalid regions, accelerate global convergence, and protect the stability of high-frequency inversion. Furthermore, although an initial velocity model is still required, the mask-guided gradient update enhances tolerance and correction capabilities for erroneous model regions, reduces the stringent dependence on the accuracy of the initial model, and strengthens the ability to correct erroneous models. Simultaneously, the mask originates from data-driven imaging results, reducing subjective bias from human interpretation and improving automation and robustness.

[0030] In some embodiments of the present invention, such as Figure 2 As shown, step S101 includes: S201. Obtain the true velocity model of the target area; S202. Perform Gaussian smoothing filtering on the real velocity model to obtain the initial velocity model.

[0031] It should be noted that: initial velocity model This is the starting point of FWI, representing a priori estimate of subsurface velocity. Its acquisition methods typically include: conventional velocity analysis; tomographic imaging results; and artificially constructed smoothed models. This invention validates the model by extracting 600×250 sampling points from a two-dimensional Marmousi model (i.e., the true velocity model). The initial velocity model is generated after Gaussian smoothing of the two-dimensional Marmousi model with sigma=20. Simultaneously, it is input into both traditional full-waveform inversion and full-waveform inversion based on Gaussian mixture clustering. In the two-dimensional seismic velocity model, x represents the horizontal (lateral) spatial coordinate; z represents the vertical (depth) spatial coordinate. Therefore, the velocity model... Indicates the position as The velocity of the underground medium was determined. Specifically, the internationally recognized two-dimensional Marmousi velocity model, used for algorithm testing, was selected. From the complete Marmousi velocity model, a rectangular region with 600 horizontal sampling points and 250 vertical sampling points was extracted as the actual velocity model for this embodiment. For the real velocity model Apply a Gaussian smoothing filter (standard deviation σ=20).

[0032] In this embodiment, the real velocity model is processed by Gaussian smoothing filtering through S201 and S202 to construct an initial velocity model that deviates significantly from the actual underground conditions and has poor quality, so as to rigorously test the robustness and effectiveness of the subsequent inversion method under unfavorable initial conditions.

[0033] In some embodiments of the present invention, such as Figure 3 As shown, step 102 includes: S301. Solve the acoustic wave equation based on the real velocity model and perform forward modeling to obtain the simulated seismic wave field. S302. Add a set proportion of Gaussian white noise to the simulated seismic wavefield to simulate the actual acquisition environment and form the observation data matrix.

[0034] It should be noted that: based on the actual speed model For the parameters of the subsurface medium, the two-dimensional acoustic wave equation is numerically solved using the finite difference method. For example, 20 seismic sources are set up, and a Ricker wavelet with a dominant frequency of 25Hz is used as the excitation signal. 3% to 10% Gaussian white noise is added to the simulated clean seismic record to simulate uncontrollable factors such as environmental noise and acquisition interference. The noisy time-series data from each shot point and received by all geophones are organized according to three dimensions: shot point, geophone point, and time, forming an observation data matrix. .

[0035] In this embodiment, the present invention obtains a true velocity model using the standard Marmousi model, and accurately simulates the actual acquisition environment based on it to form an observation data matrix. It can accurately quantify and invert the output of the final velocity model and the true velocity model. The error between them (such as mean square error, MSE).

[0036] In some embodiments of the present invention, such as Figure 4 As shown, step S104 includes: S401. Based on the Gaussian mixture model, perform probability clustering on the pixels of the reflection amplitude image. The probability density function of the Gaussian mixture model is: (Equation 1); Wherein, I represents the amplitude value of the pixel extracted from the reflection amplitude image. This indicates the preset number of cluster categories. Indicates category index, Denotes the mixing weights of the k-th Gaussian component and satisfies , This represents the mean of the k-th Gaussian component. Let the covariance matrix of the k-th Gaussian component be denoted as . The probability density function is... S402. Obtain the optimal value of the parameters by solving the probability density function, and calculate the posterior probability of each pixel belonging to each category based on the optimal value; the parameters include , and ; S403. Compare the posterior probabilities of each pixel belonging to each category, and classify each pixel into the category with the highest posterior probability. S404. Based on the partitioning results, compare the mean values ​​of the Gaussian components corresponding to each category, and determine the category with the largest mean value as the construction region. S405. Mark all pixels in the reflection amplitude image that are assigned to the constructed region as a first value, and mark pixels that are assigned to other categories as a second value, thereby generating the binarized subject mask.

[0037] It should be noted that a Gaussian mixture model is used to analyze the reflection amplitude image. The amplitude value distribution of all pixels in the image is modeled. The Gaussian mixture model assumes that the pixel values ​​of the entire reflection amplitude image are composed of K different Gaussian distributions (i.e., components or categories) mixed in a certain proportion, and its joint probability density function is shown in Equation 1. The two-dimensional reflection amplitude image is then used as a model. All pixel values ​​are flattened into a one-dimensional vector, which serves as the input dataset for the GMM algorithm. Each data sample represents the amplitude value of a single pixel. The expectation-maximization algorithm is used to solve for the probability density function. That is, based on the current parameter estimate ( , , ), calculate the posterior probability of each pixel I_i belonging to the k-th Gaussian component, and use the calculated posterior probability as weights to re-estimate and update the parameters of the model. , , This process maximizes the likelihood probability of generating the observed data matrix under the current model. After iteration, the model parameters are obtained. , , The optimal estimate of ). For For each pixel (x, z) in the image, its amplitude value is substituted into the trained Gaussian Mixture Model (GMM), and the posterior probability P(k|I(x, z)) belonging to the k-th class is calculated using Bayes' theorem. For a pixel (x, z), its posterior probability P(k|I(x, z)) belonging to the K classes is compared, and the pixel is assigned to the class with the highest posterior probability value. The mean of the Gaussian components corresponding to the k clusters is examined, and the class with the larger mean μ_k is interpreted and determined as the construction region. The class with the smaller mean is interpreted as the background region. A new image with the original reflection amplitude is created. A two-dimensional matrix of identical size is used as the main mask, which is obtained by traversing all pixels (x, z). If a pixel is classified into the category determined as a "construction region" in S404, then in the main mask... The corresponding position is assigned the first value (usually 1). Otherwise, it is assigned the second value (usually 0). The final result is a binary image, i.e., the subject mask. In this context, the white area with a value of 1 represents the main structure automatically identified by the algorithm, while the black area with a value of 0 represents the background.

[0038] This invention currently primarily targets stratigraphic boundaries and background, thus dividing them into two categories: Category 1 mainly consists of boundaries or anomalies with high gradients and high amplitudes; Category 2 mainly consists of background information with low gradients and low amplitudes. However, this invention allows... Furthermore, the optimal number of clusters K is selected using BIC to address more complex geological conditions, which will not be discussed in detail here.

[0039] In this embodiment, by solving the probability density function to calculate the posterior probability of different categories, the background region and the constructed region can be effectively distinguished. This ensures that even under poor initial conditions, roughly correct spatial constraints can still be extracted, laying the foundation for successful subsequent inversion. Furthermore, based on the category classification results, all pixels in the reflection amplitude image that are classified into the constructed region are marked with a first value, and pixels classified into other categories are marked with a second value, thereby generating the binarized subject mask. Through inputting the RTM amplitude image → automatic GMM modeling and clustering → automatic determination of the constructed class based on statistical features → outputting the binarized subject mask, the gradient update intensity of different regions can be adjusted in a differentiated manner, effectively guiding subsequent inversion.

[0040] In some embodiments of the present invention, such as Figure 5 As shown, step S105 includes: S501. Optimize the main body mask; S502. Calculate the optimized main body mask and the original main body mask to obtain the boundary mask; the first value of the boundary mask represents the constructed boundary region and the second value represents the non-boundary region.

[0041] It should be noted that after clustering, the cluster labels of the pixels form a two-dimensional mask label image, which serves as the main body mask. In this context, the pixel values ​​of the subject mask are either 0 or 1 (1 indicates category 1 is the subject construction area, and 0 indicates category 2 is the background). To more accurately depict the construction boundaries, this invention modifies the subject mask... Gaussian filtering, morphological opening and closing operations are performed, followed by morphological dilation to uniformly expand all constructed regions with a value of 1 outward by one pixel, thus obtaining the dilated mask (i.e., the optimized main body mask). Then, the original main body mask is subtracted from the dilated mask. The boundary lines are extracted using this method, and the boundary mask can be represented as follows: (Equation 2). The resulting boundary mask. In the diagram, a pixel value of 1 represents a construction boundary, and a pixel value of 0 represents a non-boundary region.

[0042] In this embodiment, the original subject mask is purified. This eliminates small irregularities caused by clustering noise or imaging defects, resulting in a boundary mask with smooth and continuous regions. This ensures that the applied prior constraints are geologically reliable, thereby improving the rationality and stability of the inversion results.

[0043] In some embodiments of the present invention, such as Figure 6 As shown, step S106 includes: S601. Construct a spatial weighting function based on the main body mask and the boundary mask: (Equation 2); in, The spatial weighting function is... For the main body mask, For the boundary mask, The gain coefficient and .

[0044] S602. The weighted gradient is calculated by combining the residual gradient of the time-domain full waveform inversion data obtained by the adjoint state method with the spatial weighting function.

[0045] It should be noted that: based on the input body mask and boundary mask Construct the spatial weighting function as shown in Equation 2. In traditional time-domain full waveform inversion, without using any constraints, the residual gradient of time-domain full waveform inversion (TDFWI) can be obtained by the adjoint method. (Equation 3), the corresponding speed update amount is: ,in Indicates forward-modeling wave field, Indicates the accompanying wave field, Indicates step size factor, The unweighted gradient is the inverted gradient. This invention introduces a GMM mask constraint to spatially weight this gradient, obtaining a weighted gradient. Specifically, for the region where the main mask M=1: a gain is assigned. Promotes updates to major structural regions; for regions with boundary mask B=1: imparts stronger gain. This enhances the recovery of velocity jump positions; the original gradient update is maintained in the background region where M=0 and B=0. In each iteration, this invention applies spatial weights to the residual gradient obtained by the adjoint method, resulting in the final weighted gradient used for updating. .

[0046] In some embodiments of the present invention, such as Figure 7 As shown, step S107 includes: S701. Select the observation data matrix of multiple frequency bands in order from low frequency to high frequency, and perform iterative optimization of time domain full waveform inversion in sequence. S702. In each round of iterative optimization, the velocity model is updated using the weighted gradient of the current iteration using the iterative optimization algorithm; the update formula is:

[0047] in, The location of the underground space at the nth iteration The current speed at that location, The location of the underground space at the (n+1)th iteration The latest speed at the location, As a weighted gradient, the gain coefficient decreases exponentially with the number of iterations.

[0048] It should be noted that traditional FWI relies on data errors, while the uniform spatial application of gradients can lead to insufficient updates of the constructed entity, blurred boundaries, severe erroneous updates in the velocity background region, getting trapped in local extrema or exhibiting spurious velocity anomalies, and high sensitivity to the initial velocity model. Therefore, after completing the spatial mask construction and gradient weighting design based on GMM, this invention integrates the obtained spatial gradient weighting strategy into the overall solution framework of Time Domain Full Waveform Inversion (TDFWI) to construct a structure-aware inversion update method. The traditional FWI objective function (full waveform inversion objective function) without any prior constraints is: (Equation 4), where Indicates the observed earthquake records, For the current velocity model The simulation data obtained from forward modeling are used to select observation data matrices from multiple frequency bands in order from low frequency to high frequency. Iterative optimization of time-domain full waveform inversion is then performed sequentially. In each iteration, the velocity model is updated using the weighted gradient of the current iteration, and the final velocity model formula is: (Equation 7). Wherein... This is the learning rate. The location of the underground space at the nth iteration The current speed at that location, The location of the underground space at the (n+1)th iteration The latest speed at the location, This is a weighted gradient.

[0049] This invention uses the Deepwave framework in Python to implement TDFWI (Time Domain Full Waveform Inversion) and GMM-FWI (Gaussian Mixture Model Constrained Full Waveform Inversion), employing an optimization algorithm (L-BFGS optimizer) to avoid overfitting. Gradient weighting coefficients are used. , Gradually decays with iteration: (Equation 8), where, This is the decay constant, used to control the decay rate.

[0050] In this embodiment, the present invention initially obtains the amplitude image of the inversion region through reverse time migration, extracts the location information of the boundary and anomalous body regions in the image using Gaussian mixture clustering, and converts the location information into a mask as a priori spatial constraint to be embedded in the traditional two-dimensional time domain full waveform inversion. Through the design of regularization terms, the performance and recognition accuracy of the full waveform inversion are improved, and high-precision velocity modeling of deep underground space is realized, which can effectively solve the problems faced by traditional full waveform inversion.

[0051] This invention uses the Marmousi seismic velocity model as the research object, selecting a region of 600 (horizontal) × 250 (vertical) sampling points as the true velocity model. The spatial sampling interval is dx = dz = 5m, and the model velocity range is approximately 1500~4500m / s, including multiple layers of velocity discontinuous stratigraphic structures. The implementation process, based on a combination of reverse time migration imaging (RTM), Gaussian mixture clustering (GMM), and full waveform inversion (FWI), includes the following specific steps: An initial velocity model was constructed, using the selected Marmousi model region as the true velocity model. To simulate the situation of insufficient initial information in actual exploration, the true velocity model was Gaussian smoothed: first, the velocity was converted to a slower velocity (1 / velocity), a Gaussian filter with a standard deviation Sigma=20 was applied, and then it was converted back to the velocity domain. This processing filtered out high-frequency details and sharp boundaries, retaining only the macroscopic background velocity trend. The resulting smooth model served as the initial velocity model for subsequent inversion and reverse time migration.

[0052] The synthetic observation data matrix is ​​generated by performing forward modeling of the acoustic wave equations on a constructed realistic Marmousi model. Twenty seismic sources were set up with a spacing of 30m and a burial depth of 5m. The excitation signal used a Ricker wavelet with a dominant frequency of 25Hz. Six hundred receivers were horizontally deployed along the ground surface with a spacing of 5m and a burial depth of 5m. The time sampling interval was dt = 0.0004s, and the total recording time was 2 seconds. 3%–10% Gaussian white noise was added to the simulated data obtained from the forward modeling to simulate the actual acquisition environment, ultimately forming the observation data matrix. In the inversion, a sequence frequency f = {3, 5, 7, 10, 12, 15, 18, 20} Hz was used for multi-scale inversion, utilizing the stable convergence of low-frequency (<10Hz) data to recover large-scale structures, and high-frequency data to recover fine details.

[0053] Loading and processing the GMM boundary mask involves performing reverse time migration (RTM) imaging using the constructed initial velocity model and the generated observation data matrix to obtain a preliminary reflection amplitude image (which reflects the approximate outline of the structure, such as...). Figure 8 (As shown). The imaging results are saved as a 600×250 two-dimensional matrix. Subsequently, Gaussian Mixture Model (GMM) clustering is performed on this reverse-time-migrated amplitude image, with the number of classes set to 2 to distinguish between the main subsurface structures and the background area. The clustering process divides the pixels into different categories (clustering results are visualized as shown). Figure 9 As shown in the figure), the output is converted into a binary mask image. By performing morphological processing on this mask, the significantly marked boundary information can be extracted to form a boundary mask for constraint inversion (the mask extraction result is shown in the figure). Figure 10 (As shown).

[0054] The time-domain full waveform inversion based on hybrid constraints involves optimizing the velocity model through full waveform inversion after obtaining the boundary mask information. The specific process is as follows: 1) The initial velocity model is set to a smooth velocity field with values ​​ranging from 1400 to 1600 m / s; the inversion frequency is set to... Hz; 2) The L-BFGS optimizer, a gradient descent optimization method based on the Deepwave library of the deep learning platform, is adopted. The learning rate is dynamically adjusted between 1e-3 and 1e-4 to ensure stable convergence of the inversion process. 3) Use the example code TDFWI.py from the Deepwave library to iterate the initial velocity model, executing the inversion loop 200 times, and output the TDFWI results (for comparison with the real model, please refer to...). Figure 12 (As shown in the trend of the middle section), the mean square error (MSE) between the initial velocity model and the true model is 4219.329590, and the mean square error (MSE) between the inverted velocity model and the true model is 2826.648682, with an improvement rate of 33.01%. 4) The extracted boundary mask is introduced into the inversion as a spatial regularization term, and the weight factor of the boundary constraints is set as follows: In each iteration, the gradient is weighted using a mask, giving higher update weights to boundary regions and lower weights to background regions. This is the GMM-constrained full waveform inversion (GMM-FWI) described in this invention.

[0055] 5) Execute GMM-FWI, loop 200 times, and output the final velocity model (the final result is compared with the two-dimensional model of the real model, for example...). Figure 11 (As shown). The calculated MSE of the final model compared to the real model was 2076.961914, representing an improvement of 50.78% compared to the initial model.

[0056] 6) Output the loss value and gradient norm curve during the GMM-FWI inversion process (e.g.) Figure 13 As shown in the figure, both gradually stabilize with increasing iterations, indicating that the inversion process converges well. Furthermore, velocity curves at the mid-horizontal and depth positions of the model are extracted and compared (e.g., ...). Figure 14 and Figure 15 As shown in the figure, the velocity curve obtained by inversion can closely approximate the real velocity curve at different locations, indicating that the method of the present invention effectively restores the stratigraphic structure, anomalies and lateral variation characteristics of the strata.

[0057] To better implement the Gaussian mixture clustering constrained full waveform inversion method in the embodiments of the present invention, based on the Gaussian mixture clustering constrained full waveform inversion method, the embodiments of the present invention also provide a Gaussian mixture clustering constrained full waveform inversion system, which includes: The initial model building module is used to build the initial velocity model; The observation data matrix synthesis module is used to obtain the observation data matrix based on the real velocity model. The acquisition module is used to perform reverse time migration imaging using the initial velocity model and the observation data matrix to obtain a reflection amplitude image; The main body mask generation module is used to generate a main body mask based on the Gaussian mixture model and the reflection amplitude image; A boundary mask generation module is used to obtain a boundary mask based on the main body mask; The gradient weighting processing module is used to obtain a weighted gradient by combining the main body mask and the boundary mask with time-domain full waveform inversion. The iterative inversion module is used to iteratively update the initial velocity model based on the weighted gradient using a multi-scale frequency strategy and optimization algorithm until the convergence condition is met, and output the final velocity model.

[0058] The Gaussian mixture clustering constrained full waveform inversion system provided in the above embodiments can realize the technical solutions described in the above Gaussian mixture clustering constrained full waveform inversion method embodiments. The specific implementation principles of each module or unit can be found in the corresponding content of the above Gaussian mixture clustering constrained full waveform inversion method embodiments, which will not be repeated here.

[0059] The present invention also provides an electronic device. This electronic device includes a processor, a memory, and a display. Only some components of the electronic device have been shown above; however, it should be understood that it is not required to implement all the shown components, and more or fewer components may be implemented alternatively.

[0060] In some embodiments, the processor may be a central processing unit (CPU), a microprocessor, or other data processing chip, used to run program code stored in memory or process data, such as the Gaussian mixture clustering constrained full waveform inversion method in this invention.

[0061] In some embodiments, the processor may be a single server or a group of servers. The server group may be centralized or distributed. In some embodiments, the processor may be local or remote. In some embodiments, the processor may be implemented on a cloud platform. In one embodiment, the cloud platform may include a private cloud, a public cloud, a hybrid cloud, a community cloud, a distributed cloud, an intranet, a multi-cloud, or any combination thereof.

[0062] In some embodiments, the memory can be an internal storage unit of an electronic device, such as a hard drive or RAM. In other embodiments, the memory can also be an external storage device of the electronic device, such as a plug-in hard drive, a smart media card (SMC), a secure digital (SD) card, a flash card, etc.

[0063] Furthermore, memory can include both internal storage units and external storage devices within an electronic device. Memory is used to store application software and various types of data installed in the electronic device.

[0064] In some embodiments, the display may be an LED display, a liquid crystal display, a touch-sensitive liquid crystal display, or an OLED (Organic Light-Emitting Diode) touchscreen. The display is used to show information from an electronic device and to display a visual user interface. Components of the electronic device communicate with each other via a system bus.

[0065] In one embodiment, when the processor executes the Gaussian mixture clustering constrained full waveform inversion program in memory, the following steps can be performed: Construct an initial velocity model; The observation data matrix is ​​obtained based on the true velocity model; The initial velocity model and the observation data matrix are used to perform reverse time migration imaging to obtain a reflection amplitude image; A subject mask is generated based on the Gaussian mixture model and the reflection amplitude image; The boundary mask is obtained based on the main body mask; The weighted gradient is obtained by combining the main body mask and the boundary mask with time-domain full waveform inversion. A multi-scale frequency strategy and optimization algorithm are adopted to iteratively update the initial velocity model based on the weighted gradient until the convergence condition is met, and the final velocity model is output.

[0066] It should be understood that when the processor executes the Gaussian mixture clustering constrained full waveform inversion program in memory, in addition to the functions mentioned above, it can also perform other functions, as detailed in the description of the corresponding method embodiments above.

[0067] Furthermore, the embodiments of the present invention do not specifically limit the type of electronic device mentioned. The electronic device can be a mobile phone, tablet computer, personal digital assistant (PDA), wearable device, laptop computer, or other portable electronic device. Exemplary embodiments of portable electronic devices include, but are not limited to, portable electronic devices running iOS, Android, Microsoft, or other operating systems. The aforementioned portable electronic devices can also be other portable electronic devices, such as laptop computers with touch-sensitive surfaces (e.g., touch panels). It should also be understood that in some other embodiments of the present invention, the electronic device may not be a portable electronic device, but rather a desktop computer with a touch-sensitive surface (e.g., a touch panel).

[0068] Accordingly, this application also provides a computer-readable storage medium for storing computer-readable programs or instructions. When the programs or instructions are executed by a processor, they can implement the steps or functions of the Gaussian mixture clustering constrained full waveform inversion method provided in the above-described method embodiments.

[0069] Those skilled in the art will understand that all or part of the processes of the methods described in the above embodiments can be implemented by a computer program instructing related hardware (such as a processor, controller, etc.), and the computer program can be stored in a computer-readable storage medium. The computer-readable storage medium may be a disk, optical disk, read-only memory, or random access memory, etc.

[0070] The Gaussian mixture clustering constrained full waveform inversion method, system, device and medium provided by the present invention have been described in detail above. Specific examples have been used to illustrate the principle and implementation of the present invention. The description of the above embodiments is only for the purpose of helping to understand the method and core idea of ​​the present invention. At the same time, for those skilled in the art, there will be changes in the specific implementation and application scope based on the idea of ​​the present invention. Therefore, the content of this specification should not be construed as a limitation of the present invention.

Claims

1. A Gaussian mixture clustering constrained full waveform inversion method, characterized in that, include: Construct an initial velocity model; The observation data matrix is ​​obtained based on the true velocity model; The initial velocity model and the observation data matrix are used to perform reverse time migration imaging to obtain a reflection amplitude image; A subject mask is generated based on the Gaussian mixture model and the reflection amplitude image; The boundary mask is obtained based on the main body mask; The weighted gradient is obtained by combining the main body mask and the boundary mask with time-domain full waveform inversion. A multi-scale frequency strategy and optimization algorithm are adopted to iteratively update the initial velocity model based on the weighted gradient until the convergence condition is met, and the final velocity model is output.

2. The method according to claim 1, characterized in that, The construction of the initial velocity model includes: Obtain the true velocity model of the target area; The initial velocity model is obtained by performing Gaussian smoothing filtering on the actual velocity model.

3. The method according to claim 1, characterized in that, The observation data matrix obtained based on the true velocity model includes: Based on the real velocity model, the acoustic wave equation is solved, and forward modeling is performed to obtain the simulated seismic wave field. A set proportion of Gaussian white noise is added to the simulated seismic wavefield to simulate the actual acquisition environment, forming the observation data matrix.

4. The method according to claim 1, characterized in that, The step of generating the subject mask based on the Gaussian mixture model and the reflection amplitude image includes: Based on the Gaussian mixture model, the pixels of the reflection amplitude image are probabilistically clustered. The probability density function of the Gaussian mixture model is: ； in, I This represents the amplitude value of a pixel extracted from the reflection amplitude image. This indicates the preset number of cluster categories. Indicates category index, Denotes the mixing weights of the k-th Gaussian component and satisfies , This represents the mean of the k-th Gaussian component. Let the covariance matrix of the k-th Gaussian component be denoted as . The probability density function is... The optimal values ​​of the parameters are obtained by solving the probability density function, and the posterior probability of each pixel belonging to each category is calculated based on the optimal values; the parameters include , and ; Compare the posterior probabilities of each pixel belonging to each category, and assign each pixel to the category with the highest posterior probability. Based on the partitioning results, the mean values ​​of the Gaussian components corresponding to each category are compared, and the category with the largest mean value is determined as the construction region. All pixels in the reflection amplitude image that are assigned to the constructed region are marked with a first value, and pixels that are assigned to other categories are marked with a second value, thereby generating the binarized subject mask.

5. The method according to claim 1, characterized in that, The process of obtaining the boundary mask based on the main body mask includes: The main body mask is optimized. The optimized main body mask is calculated together with the original main body mask to obtain the boundary mask; the first value of the boundary mask represents the constructed boundary region and the second value represents the non-boundary region.

6. The method according to claim 1, characterized in that, The step of obtaining the weighted gradient by combining the main body mask and the boundary mask with time-domain full waveform inversion includes: Construct a spatial weighting function based on the main body mask and the boundary mask: ； in, The spatial weighting function is... For the main body mask, For the boundary mask, The gain coefficient and ; The weighted gradient is calculated by combining the residual gradient of the time-domain full waveform inversion data obtained by the adjoint state method with the spatial weighting function.

7. The method according to claim 6, characterized in that, The method employs a multi-scale frequency strategy and optimization algorithm to iteratively update the initial velocity model based on the weighted gradient, including: The observation data matrices of multiple frequency bands are selected in order from low frequency to high frequency, and iterative optimization of time domain full waveform inversion is performed sequentially. In each round of iterative optimization, the velocity model is updated using the weighted gradient of the current iteration using the iterative optimization algorithm; the update formula is: in, The location of the underground space at the nth iteration The current speed at that location, The location of the underground space at the (n+1)th iteration The latest speed at the location, As a weighted gradient, the gain coefficient decreases exponentially with the number of iterations.

8. A Gaussian mixture clustering constrained full waveform inversion system, characterized in that, include: The initial model building module is used to build the initial velocity model; The observation data matrix synthesis module is used to obtain the observation data matrix based on the real velocity model. The acquisition module is used to perform reverse time migration imaging using the initial velocity model and the observation data matrix to obtain a reflection amplitude image; The main body mask generation module is used to generate a main body mask based on the Gaussian mixture model and the reflection amplitude image; A boundary mask generation module is used to obtain a boundary mask based on the main body mask; The gradient weighting processing module is used to obtain a weighted gradient by combining the main body mask and the boundary mask with time-domain full waveform inversion. The iterative inversion module is used to iteratively update the initial velocity model based on the weighted gradient using a multi-scale frequency strategy and optimization algorithm until the convergence condition is met, and output the final velocity model.

9. An electronic device, characterized in that, Including memory and processor, among which, The memory is used to store programs; The processor, coupled to the memory, is used to execute the program stored in the memory to implement the steps in the Gaussian mixture clustering constrained full waveform inversion method according to any one of claims 1 to 7.

10. A computer-readable storage medium, characterized in that, Used to store computer-readable programs or instructions, which, when executed by a processor, can implement the steps in the Gaussian mixture clustering constrained full waveform inversion method according to any one of claims 1 to 7.

Images