Surface wave dispersion 3d tomography method based on adaptive dictionary constraint

Abstract

The application belongs to the technical field of seismology and geophysical tomography, and relates to a surface wave dispersion three-dimensional tomography method based on adaptive dictionary constraint, which comprises the following steps: performing three-dimensional grid dissection on a research area, constructing a model parameter vector for representing a three-dimensional shear wave velocity field, constructing a linearized sensitivity matrix equation between a travel time residual and a shear wave velocity disturbance, obtaining a sensitivity matrix and a column norm distribution of the sensitivity matrix, equivalently converting maximum joint posterior probability density function into minimization according to Bayes theorem, obtaining a joint inversion objective function composed of a travel time data fitting term, a global-assisted coupled constraint term and an image block sparse representation error term, and solving the joint inversion objective function by using an alternating minimization strategy, and outputting a final three-dimensional shear wave velocity structure model after iterative convergence. The optimal balance between data constraint and structure prior is achieved, and high-fidelity imaging of the three-dimensional shear wave velocity structure is realized.

Description

Technical Field

[0001] This application belongs to the field of seismology and geophysical tomography technology, specifically relating to a surface wave dispersion three-dimensional tomography method based on adaptive dictionary constraints. Background Technology

[0002] Surface wave dispersion three-dimensional direct tomography (SMT) technology directly constructs a three-dimensional velocity model of the subsurface medium by jointly inverting multi-period surface wave dispersion data. This effectively avoids the accumulation of errors in traditional two-step imaging methods, demonstrating significant advantages in acquiring high-resolution crustal velocity structures. However, the uneven distribution of stations in actual observations leads to significant differences in ray coverage density, resulting in highly uncertain imaging conditions. While traditional global smoothing regularization methods can address this, they cannot completely eliminate the inherent limitations of traditional methods. While norm-constrained stable inversion can overcompress velocity anomaly amplitudes and blur tectonic boundaries in ray-sparse regions, it severely limits the ability to identify fine structures such as fault zones and lithological interfaces. In recent years, sparse representation and dictionary learning methods have been gradually introduced into seismic tomography to overcome these limitations. However, existing works often use fixed pre-trained dictionaries or simple wavelet bases as priors for sparse transformation. These priors lack adaptability to the actual geological structure of the region to be inverted and cannot be dynamically updated with the velocity model during the inversion iteration process. This leads to a mismatch between sparse constraints and local structural features, making it difficult to achieve a smooth velocity background and discontinuous tectonic boundaries in strongly inhomogeneous media. Summary of the Invention

[0003] To address the problems in existing surface wave dispersion three-dimensional direct tomography methods, such as insufficient ray coverage leading to excessive suppression of abnormal velocity amplitudes and blurred boundary construction due to global smoothing regularization, this application provides a surface wave dispersion three-dimensional tomography method based on adaptive dictionary constraints.

[0004] This application provides a surface wave dispersion three-dimensional tomography method based on adaptive dictionary constraints, including: The study area is divided into three-dimensional meshes to construct a model parameter vector to characterize the three-dimensional shear wave velocity field. The model parameter vector consists of multiple shear wave velocity perturbations. The continuous slow field is discretized by ray tracing and bilinear interpolation to obtain the theoretical travel time. The observed travel time is calculated based on the observation data between station pairs. The travel time residual is constructed according to the theoretical travel time and the observed travel time. A linearized sensitivity matrix equation between the travel time residual and the shear wave velocity perturbation is constructed to obtain the sensitivity matrix and the norm distribution of the sensitivity matrix. Based on the sensitivity matrix, the likelihood function of the observed data is modeled as a Gaussian distribution. An auxiliary sparse model is introduced as a latent variable to establish a two-level Gaussian prior distribution, including a first-level Gaussian prior distribution and a second-level conditional prior distribution. The product of the likelihood function and the two-level Gaussian prior distribution constitutes the joint posterior probability density function. According to Bayes' theorem, maximizing the joint posterior probability density function is equivalently transformed into minimizing it, resulting in a joint inversion objective function composed of the weighted sum of three terms: the travel time data fitting term, the global-auxiliary coupling constraint term, and the image patch sparse representation error term. In the early stage of inversion, a warm-up iteration is performed, and the initial velocity background is constructed using Laplace smoothing regularization. After the warm-up, the dictionary sparse regularization is switched to, and the joint inversion objective function is solved using an alternating minimization strategy. After iterative convergence, the final three-dimensional shear wave velocity structure model is output.

[0005] Furthermore, when the likelihood function of the observed data is modeled as a Gaussian distribution, the mean of the Gaussian distribution is the sensitivity matrix multiplied by the difference between the global velocity model and the current velocity model. The difference represents the predicted travel time change value generated by the velocity disturbance of the current velocity model relative to the linearized expansion point after being mapped by the sensitivity matrix. The variance of the Gaussian distribution is the travel time observation error variance, which is used to characterize the noise level in the observed data. The weight of the likelihood function is controlled by the reciprocal of the travel time observation error variance.

[0006] Furthermore, the first-level Gaussian prior distribution in the two-level Gaussian prior distribution is a Gaussian distribution centered on the auxiliary sparse model for the global velocity model. The mean of the first-level Gaussian prior distribution is taken from the auxiliary sparse model, and the variance is taken from the prior variance parameter. The prior variance parameter quantitatively characterizes the range within which the global velocity model is allowed to fluctuate around the auxiliary sparse model. Based on the value of the prior variance parameter, it is used to select whether the construction boundary and velocity discontinuity features carried in the auxiliary sparse model are completely transmitted to the global velocity model, or whether the global velocity model deviates from the auxiliary sparse model.

[0007] Furthermore, the second-level conditional prior distribution is the conditional prior distribution of the auxiliary sparse model based on the sparse representation of the image patch dictionary, and it is constructed as follows: The lateral velocity distribution of the auxiliary sparse model at each depth layer is divided into several overlapping spatial image patches, each containing a fixed number of velocity grid points. The velocity value of each spatial image patch is extracted from the auxiliary sparse model using an extraction operator; Each spatial image patch is approximated by a sparse linear superposition of atoms in an overcomplete dictionary. Each spatial image patch corresponds to a sparse representation coefficient. The number of non-zero elements in the sparse representation coefficient must not exceed the preset sparsity upper limit. The sparse representation coefficients of each spatial image patch are statistically independent of each other. Using the dictionary sparse approximation error of each spatial image patch as the mean and the variance of the velocity approximation error of the spatial image patch as the variance, the second-level conditional prior distribution of the auxiliary sparse model is constructed as the product of the independent Gaussian distributions of each spatial image patch.

[0008] Furthermore, the travel time data fitting term is used to calculate the sum of squares of the difference between the observed travel time and the sensitivity matrix multiplied by the velocity disturbance, and the travel time data fitting term is weighted by the reciprocal of the variance of the travel time observation error; The global-auxiliary coupling constraint term is used to calculate the sum of squares of the velocity differences between the global velocity model and the auxiliary sparse model at each grid node, and is weighted by the reciprocal of the prior variance parameter. The image patch sparse representation error term is used to calculate the sum of squares of the differences between the velocity values ​​of each spatial image patch and the approximations of the sparse linear combination of dictionary atoms, and is accumulated across all spatial image patches. It is also weighted by the reciprocal of the variance of the spatial image patch velocity approximation error.

[0009] Furthermore, the joint inversion objective function is multiplied by the travel time observation error variance and normalized to normalize the coefficient of the travel time data fitting term to 1. A first tradeoff coefficient is defined as the coefficient of the global-auxiliary coupling constraint term to control the weight between the travel time data fitting term and the global-auxiliary coupling constraint term. A second tradeoff coefficient is defined as the coefficient of the image patch sparse representation error term to control the weight between the travel time data fitting term and the image patch sparse representation error term.

[0010] Furthermore, an alternating minimization strategy is employed to solve the joint inversion objective function, including: The joint optimization problem is broken down into three subproblems, which are solved iteratively, optimizing only one subproblem at a time while keeping the other two variables fixed: The first subproblem is the global velocity field update. With the auxiliary sparse model and sparse representation coefficients fixed, the LSMR algorithm is used to solve the augmented least squares equations with the sensitivity matrix as the core coefficient matrix, thus completing the update of the global velocity model. The second subproblem is the estimation of sparse representation coefficients. With the global velocity model and the auxiliary sparse model fixed, after the mean is removed from each spatial image patch, the sparse representation coefficients of each spatial image patch are solved one by one using the orthogonal matching pursuit algorithm. The third subproblem is solving the auxiliary sparse model. By fixing the global velocity model and the sparse representation coefficient matrix, the closed-form solution of the auxiliary sparse model is obtained using cyclic boundary conditions.

[0011] Furthermore, in the sparse representation coefficient estimation, the column norm distribution of the sensitivity matrix is ​​used to perform ray coverage checks on each spatial image block, and the proportion of pixels with zero ray density within the spatial image block is counted. When the proportion of pixels exceeds a preset threshold, the corresponding spatial image block is considered to be in a region without effective ray constraints and is removed from the sparse representation coefficient estimation.

[0012] Furthermore, during the solution process of the auxiliary sparse model, information sharing is achieved between spatial image blocks through overlapping regions. The velocity value of the overlapping part is simultaneously constrained by the dictionary reconstruction values ​​of multiple spatial image blocks. After projecting the reconstruction values ​​of each spatial image block back to the complete grid using cyclic boundary conditions, the final velocity value of each pixel is equal to the weighted average of the reconstruction values ​​of all spatial image blocks covering that pixel at that pixel position.

[0013] The closed-form solution of the auxiliary sparse model is solved independently for each depth layer without introducing depth direction regularization, and the solution of each depth layer does not affect each other; when depth direction regularization is introduced, the solution of each depth layer is coupled through a tridiagonal linear equation system.

[0014] Compared with the prior art, the advantages of this application are as follows: This application systematically overcomes the shortcomings of traditional smoothing methods in over-suppressing velocity discontinuities and blunting fine-scale structures in regions with insufficient ray coverage by using adaptive dictionary sparse regularization under the Bayesian framework. It achieves the optimal balance between data constraints and structural priors, and realizes high-fidelity imaging of three-dimensional shear wave velocity structures. Attached Figure Description

[0015] Figure 1 A flowchart of a surface wave dispersion three-dimensional tomography method based on adaptive dictionary constraints provided in this application embodiment; Figure 2 The distribution map of research area stations provided for embodiments of this application; Figure 3 The image shows the results of the checkerboard test at a depth of 5 km. (a) shows the phase velocity measurement and ray path distribution of the 4–6 s period Rayleigh wave; (b) and (f) are the theoretical input models; (c) and (g) are the comparison method and the inversion results of this application, respectively; (d) and (h) are magnified views of the areas shown in the pink dashed boxes in (c) and (g), respectively; and (e) shows the ray path density distribution. Figure 4The image shows the 40km depth checkerboard pattern reconstruction results provided in this application embodiment, where (a) is the Rayleigh wave phase velocity measurement and ray path distribution with a period of 39–41s; (b) and (f) are both theoretical input models; (c) and (g) are the comparison method and the inversion results of this application, respectively; (d) and (h) are enlarged views of the areas shown in the pink dashed boxes in (c) and (g), respectively; and (e) is the ray path density distribution. Figure 5 The image shows the peak test recovery results at a depth of 35 km provided for the embodiments of this application, where (a) is the Rayleigh wave phase velocity measurement and ray path distribution with a period of 34–36 s; (b) and (f) are theoretical input models; (c) and (g) are the inversion results of the comparison method and the method of this application, respectively; (d) and (h) are magnified views of the areas shown in the pink dashed boxes in (c) and (g), respectively; and (e) is the ray path density distribution. Figure 6 The following are the results of the one-dimensional velocity recovery curves extracted along the AA′, BB′ and CC′ measurement lines in the checkerboard test provided in the embodiments of this application: (a) is the one-dimensional velocity recovery curve of the AA′ measurement line, (b) is the root mean square error of the AA′ measurement line, (c) is the one-dimensional velocity recovery curve of the BB′ measurement line, (d) is the root mean square error of the BB′ measurement line, (e) is the one-dimensional velocity recovery curve of the CC′ measurement line, and (f) is the root mean square error of the CC′ measurement line. Figure 7 Figure 1 shows the effect of the sparsity upper limit T on the root mean square error in the chessboard recovery test at different depth layers provided in the embodiments of this application. Figure 8 The graph shows the effect of dictionary size Q on root mean square error in chessboard recovery tests at different depths, as provided in the embodiments of this application. Figure 9 First trade-off coefficient provided for embodiments of this application Second tradeoff coefficient The sensitivity analysis results are shown in the figure. (a) is the distribution of the root mean square error as a function of parameter combinations; (b) is the distribution of the root mean square error with a fixed second tradeoff coefficient. The time mean square error varies with λ1; (c) is the fixed first tradeoff coefficient. Root mean square error varies with the second weighting factor The changes.

[0016] Figure 10 Depth coupling coefficient in checkerboard recovery tests at different depths provided in embodiments of this application The effect of the root mean square error is shown in the figure. Figure 11The morphological evolution of 5km depth dictionary atoms provided in this application embodiment before and after learning by the ITKM algorithm, where (a) is the morphology before learning, (b) is the morphology after learning, and (c) is the morphological change of atoms arranged neatly in rows and columns before and after learning. Detailed Implementation

[0017] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the scope of this application.

[0018] See Figure 1 As shown, a surface wave dispersion three-dimensional tomography method based on adaptive dictionary constraints includes: S101. The study area is divided into three-dimensional meshes to construct a model parameter vector to characterize the three-dimensional shear wave velocity field. The model parameter vector consists of multiple shear wave velocity perturbations. The continuous slow field is discretized by ray tracing and bilinear interpolation to obtain the theoretical travel time. The observed travel time is calculated based on the observation data between station pairs. The travel time residual is constructed according to the theoretical travel time and the observed travel time. The linearized sensitivity matrix equation between the travel time residual and the shear wave velocity perturbation is constructed to obtain the sensitivity matrix and the norm distribution of the sensitivity matrix. S102, based on the sensitivity matrix, models the likelihood function of the observed data as a Gaussian distribution; introduces an auxiliary sparse model as a latent variable, and establishes a two-level Gaussian prior distribution including a first-level Gaussian prior distribution and a second-level conditional prior distribution; the product of the likelihood function and the two-level Gaussian prior distribution constitutes a joint posterior probability density function, and according to Bayes' theorem, maximizing the joint posterior probability density function is equivalently transformed into minimizing it, resulting in a joint inversion objective function composed of the weighted sum of squares of three terms: travel time data fitting term, global-auxiliary coupling constraint term, and image patch sparse representation error term; S103, in the early stage of inversion, a warm-up iteration is performed, and the initial velocity background is constructed by Laplace smoothing regularization. After the warm-up, the dictionary sparse regularization is switched, and the joint inversion objective function is solved by an alternating minimization strategy. S104, after iterative convergence, outputs the final three-dimensional shear wave velocity structure model.

[0019] In step S101, the study area is divided into three-dimensional meshes, and a total of [number] meshes are set. One horizontal grid node, This represents the number of nodes in the meridional grid. This represents the number of nodes in the latitudinal grid. The total number of horizontal grid nodes, and the number of each grid node. Included in the vertical direction depth sampling points Thus, an independent one-dimensional vertical profile is constructed. The model parameter vector is expressed as: , in, For model parameter vectors, Indicates the total number of model parameters. This represents the total number of depth sampling points. For the first Each horizontal grid point at depth The shear wave velocity at a given depth is the perturbation of the shear wave velocity relative to a reference velocity model. The reference velocity model is a pre-defined velocity structure established before inversion. The model parameter vector is composed of the perturbations of the shear wave velocity at each horizontal grid node relative to the reference velocity model at each depth sampling point. Based on observation data between stations, the station pair... Hetai Station Frequency between The travel time of the phase velocity of a Rayleigh wave can be expressed as the integral of its slowness along the path: , in On the representative path The phase velocity slowness value at point [point] For the station Hetai Station The propagation path between them For the station and Between, angular frequency The travel time of the Rayleigh wave phase velocity at a given location. After discretizing the continuous slow-motion field using bilinear interpolation, the th... Ray path The corresponding time-of-flight expression, i.e., the first... The ray path at angular frequency The theoretical timeout is: The theoretical travel time is obtained by discretization through ray tracing and bilinear interpolation under the current velocity model: for each ray path, the phase velocity slowness of each horizontal grid node is weighted and summed according to the bilinear interpolation weights along its propagation path.

[0020] in The phase velocity slowness of the k-th grid node is calculated using the Haskell-Thomson matrix method based on a one-dimensional hierarchical velocity model. Represents the forward calculation function, interpolation weights Defined as: , The summation range in the formula covers the four corner nodes involved in bilinear interpolation. , The coefficients are the bilinear interpolation shape functions. The corresponding path length weighting coefficients are used. Ray tracing is performed in a three-dimensional velocity field using a fast-advance method and is dynamically updated in each outer iteration. The outer iteration refers to each complete loop process of re-performing ray tracing using the updated three-dimensional velocity model, reconstructing the sensitivity matrix, and then solving the inversion problem again.

[0021] The aforementioned Continuous Slowness Field refers to the continuous spatial distribution function formed by the reciprocal of the propagation velocity of seismic waves in the subsurface medium at a specific frequency at each spatial node within the study area.

[0022] slow phase velocity disturbance right Linearization, i.e., taking the partial derivative with respect to the phase velocity slowness, directly yields the travel time residual: The travel time residual is the difference between the observed travel time and the theoretical travel time, and can be written as: , in, For the first The ray path at angular frequency The observation time at the location.

[0023] In the formula, ,in, For phase velocity disturbance, The phase velocity is the relationship between the phase velocity and the phase velocity slow speed, and the phase velocity perturbation is obtained. The dependence of each elasticity parameter can be expanded using integral form as follows: , in, , These represent the perturbations of P-wave velocity and density, respectively. Each partial derivative is obtained by applying a small perturbation to each depth node. The phase velocity is then approximated using a finite difference quotient. The phase velocity is determined by a one-dimensional vertical velocity profile below the grid nodes, which in turn is determined by the P-wave velocity. Shear wave velocity and density It is composed of three elastic parameters.

[0024] After the above steps, based on Brocher's empirical relationships (a set of empirical formulas connecting different elastic parameters of crustal rocks), and Represented as The function; after discretizing the depth direction, the travel time residual with respect to the shear wave velocity perturbation. The linearized relationship can be simplified into a linearized sensitivity matrix equation: , in Represents the time residual. Let be the sensitivity matrix. The linearized sensitivity matrix equation is an underdetermined problem and requires regularization constraints for stabilization. Each column of the sensitivity matrix corresponds to a model parameter (i.e., the shear wave velocity perturbation at a certain depth layer for a certain horizontal node), and each row corresponds to a ray path (a travel time observation data). The column norm is defined as the sum of the absolute values ​​of all elements in that column.

[0025] Traditional methods employ Laplace smoothing regularization, and the objective function of traditional methods... for: , Where L is a first-order or second-order spatial finite difference operator. This is the regularization tradeoff coefficient. However, this globally uniform smoothing strategy is prone to over-suppressing velocity discontinuities and passivating fine-scale structures in regions with insufficient ray coverage.

[0026] In step S102, the sensitivity matrix has been obtained. Based on the norm distribution, a Bayesian inversion framework is further constructed. Based on the sensitivity matrix, the likelihood function of the observed data is modeled as a Gaussian distribution. The mean of the Gaussian distribution is the sensitivity matrix multiplied by the difference between the global velocity model and the current velocity model. This difference represents the predicted travel time change of the current velocity model relative to the linearized expansion point after mapping through the sensitivity matrix. The variance of the Gaussian distribution is the travel time observation error variance, used to characterize the noise level in the observed data. The weights of the likelihood function are controlled by the reciprocal of the travel time observation error variance.

[0027] Specifically, this includes: defining the three-dimensional shear wave velocity field. For the global velocity model, the following content uses The global velocity model is shown in the table, the first In the next iteration, using the current velocity model As the linear expansion point, the likelihood function of the observed data is modeled as a Gaussian distribution:

[0028] in, Let be the likelihood function of the observed data. It follows a Gaussian normal distribution. It is the identity matrix. Let Variance be the variance of the travel time observation error. This represents the mean of a Gaussian distribution. For the global velocity model (used to represent the three-dimensional shear wave velocity field, all adopting...), The difference between the travel time model and the current velocity model is used to calculate the travel time data fitting term of the joint inversion objective function by removing the negative logarithm of the likelihood function of the observed data.

[0029] In one embodiment, the first-level Gaussian prior distribution in the two-level Gaussian prior distribution is a Gaussian distribution centered on the auxiliary sparse model for the global velocity model. The mean of the first-level Gaussian prior distribution is taken from the auxiliary sparse model, and the variance is taken from the prior variance parameter. The prior variance parameter quantitatively characterizes the range within which the global velocity model is allowed to fluctuate around the auxiliary sparse model. Based on the value of the prior variance parameter, it is used to select whether the construction boundary and velocity discontinuity features carried in the auxiliary sparse model are completely transmitted to the global velocity model, or whether the global velocity model deviates from the auxiliary sparse model.

[0030] To achieve coupling between the global velocity model and the auxiliary sparse model, the three-dimensional shear wave velocity field was analyzed. Introduced to assist sparse models The first-order Gaussian prior distribution centered at: , in The prior variance parameter characterizes the allowable fluctuation range of the global velocity model around the auxiliary sparse model.

[0031] Prior variance parameter When a smaller value is taken, the probability density of the first-order Gaussian prior distribution is highly concentrated near the auxiliary sparse model, and the global velocity model is strongly constrained around the auxiliary sparse model. In this case, the construction boundary and velocity discontinuity features carried in the auxiliary sparse model are completely transferred to the global velocity model, and the final inversion result mainly reflects the prior structural features.

[0032] Prior variance parameter When a larger value is taken, the probability density of the first-order Gaussian prior distribution tends to flatten, and the global velocity model has a greater degree of freedom deviating from the auxiliary sparse model. In this case, the inversion results mainly rely on the constraints of the travel-time observation data, and the prior structure exists only as a weak constraint.

[0033] The construction boundary and velocity discontinuity features carried in the auxiliary sparse model are passed to the final global velocity model through the following mechanism: First, the direct transmission effect of the mean value. Since the mean of the first-order Gaussian prior distribution is directly taken as the auxiliary sparse model, in the process of minimizing the joint inversion objective function, the second term... It is continuously minimized. When the velocity value of a grid node in the auxiliary sparse model changes abruptly due to the existence of the construction boundary, the velocity value of the corresponding node in the global velocity model will inevitably tend towards the same abrupt change value under the drive of the coupling constraints. This direct transfer effect ensures the precise spatial correspondence of the construction boundary. Second, the constraint strength control of the global-auxiliary coupling constraint term. In the normalized joint inversion objective function, the weight of the second term is controlled by the regularization coefficient. When When the velocity is large, the global velocity model is forced to closely follow all velocity variation features of the auxiliary sparse model, including local details and construction boundaries; when At smaller scales, the globally sparse model offers greater flexibility, referencing the structure of the auxiliary model only at larger scales. This can be achieved by adjusting... The value of can be used to achieve a balance between fully inheriting the auxiliary model structure and adequately fitting the observed data.

[0034] And in the bidirectional coupling feedback in subsequent alternating optimization.

[0035] The second-level conditional prior distribution is the conditional prior distribution of the auxiliary sparse model based on the sparse representation of the image patch dictionary. It is constructed as follows: The lateral velocity distribution of the auxiliary sparse model at each depth layer is divided into several overlapping spatial image patches, each containing a fixed number of velocity grid points. The velocity value of each spatial image patch is extracted from the auxiliary sparse model using an extraction operator; Each spatial image block is approximated by a sparse linear superposition of atoms in an overcomplete dictionary. Each spatial image block corresponds to a sparse representation coefficient. The number of non-zero elements in the sparse representation coefficient must not exceed the preset sparsity upper limit. The sparse representation coefficients of each spatial image block are statistically independent of each other. Using the dictionary sparse approximation error of each spatial image patch as the mean and the variance of the image patch velocity approximation error as the variance, the conditional prior distribution of the auxiliary sparse model is constructed as the product of the independent Gaussian distributions of each spatial image patch.

[0036] Specifically, this means: Setting up an auxiliary sparse model As latent variables. For auxiliary sparse velocity modulo... The lateral velocity distribution at each depth layer is divided into... There are three overlapping spatial image patches, each with a side length of 1. (Total) (Velocity grid points). Extraction operator is used. Extracting the first from the auxiliary sparse model In a spatial image patch Each velocity grid element, i.e. .

[0037] Assuming each spatial image patch can be derived from a dictionary (including) (atoms) It can be approximated by the sparse linear superposition of all atoms in the middle: , in For sparse representation coefficients, To measure sparsity, To count the number of non-zero elements, the coefficients of each spatial image patch are assumed to be statistically independent. Using the dictionary sparse approximation error of each spatial image patch as the mean and the variance of the velocity approximation error of each spatial image patch as the variance, the second-level conditional prior distribution of the auxiliary sparse model is constructed as the product of the independent Gaussian distributions of each spatial image patch. The second-level conditional prior distribution of the auxiliary sparse model is expressed as: , in The variance of the velocity approximation error of the spatial image patch. The sparse coefficient matrix of the image patch Indicates the first Sparse representation coefficients of spatial image patches.

[0038] Based on the above observation data, the likelihood function is modeled as a Gaussian distribution, with a first-order Gaussian prior distribution and a second-order conditional prior distribution. According to Bayes' theorem, the joint posterior probability is proportional to the product of the likelihood function and the prior probability. ; in, For the joint posterior probability density, Let be the likelihood function of the observed data. It is a first-order Gaussian prior distribution. This is a second-order conditional prior distribution. For a sparse prior distribution, the sparse prior distribution is the distribution of the sparse coefficient matrix. The applied prior probability distribution.

[0039] Using the conditional independence assumption, which states that under certain conditions, the probability distributions of different variables (or different data points) can be decomposed into the product of their independent components, taking the negative logarithm of each Gaussian distribution, and considering the sparse representation coefficients of the sparse prior distribution. The prior probability is obtained by taking the logarithm. Approximately This is approximated by the number of its non-zero elements. A sparse representation has more non-zero elements in its coefficients, resulting in a lower prior probability of the vector appearing. In the final joint inversion objective function, this will be further transformed into a hard constraint. , The sparsity is set as an upper limit to ensure that each spatial image patch is reconstructed from only a few dictionary atoms.

[0040] Maximizing the posterior probability is equivalent to minimizing the posterior probability, which refers to the probability obtained from the observations given the travel time. Under the condition that the global velocity model, the auxiliary sparse model, and the sparse coefficient matrix appear simultaneously, the joint posterior probability density yields the joint inversion objective function: , , in, For the time-fitting term, For global-auxiliary coupling constraint terms, The error term is used to represent the sparse representation of image patches.

[0041] make ,in, This is the first tradeoff coefficient, used to control the weights of the travel time data fitting term and the global-auxiliary coupling constraint term. This is the second tradeoff coefficient, used to control the weights of the travel time data fitting term and the global-auxiliary coupling constraint term. It will be used as... Normalizing the above equation, we get: , st , st in the above represents the upper limit constraint on sparsity, and st indicates that it is constrained by .

[0042] The following summary of the above process is provided to facilitate understanding of the solution in this embodiment: Bayesian statistical inference is used to handle the inversion problem of 3D seismic travel-time tomography. Bayes' theorem states that after obtaining observation data, the posterior probability of the model parameters is proportional to the product of the likelihood function of the observation data and the prior probability of the model parameters.

[0043] In this embodiment, the model parameters to be estimated include three variables: a global velocity model, an auxiliary sparse model, and a sparse coefficient matrix. The global velocity model is the final output 3D shear wave velocity model; the auxiliary sparse model is an introduced latent variable serving as a bridge between the global velocity model and the dictionary sparse representation; and the sparse coefficient matrix is ​​the set of sparse coefficients for all image patches. These three variables constitute the joint unknown parameter space for Bayesian inference.

[0044] According to Bayes' theorem, given travel time observation data Under these conditions, the joint posterior probability density function of the above three variables is proportional to the product of the likelihood function and the two-level Gaussian prior distribution: , It comprises four factors. The first factor is the likelihood function of the observed data, representing the probability of observing actual travel time data given a global velocity model. The second factor is the first-order Gaussian prior distribution, which couples the global velocity model with the auxiliary sparse model, forcing the global velocity model to fluctuate around the auxiliary sparse model. The third factor is the second-order conditional prior distribution, which imposes sparse representation constraints based on spatial image patches on the auxiliary sparse model. The fourth factor is the sparse prior distribution, which imposes sparsity constraints on the sparse coefficient vector itself.

[0045] The construction of the likelihood function for the observed data includes: firstly, constructing the likelihood function for the observed data based on the obtained sensitivity matrix and its column norm distribution. In the... In the next outer iteration, the current velocity model is used as the linearization expansion point to calculate the observed data. The residuals between the current velocity model predictions and the observed data are assumed to have zero mean noise and zero variance in travel time observation errors. If the data follows a Gaussian distribution and the errors between each observation point are independent, then the likelihood function of the observation data is also Gaussian.

[0046] By taking the negative logarithm of the likelihood function of the observed data and discarding the constant term, the fitting term for the travel time data is obtained.

[0047] Construction of the first-order Gaussian prior distribution: An auxiliary sparse model is introduced as a latent variable. This auxiliary sparse model has the exact same mesh dimension and total number of model parameters as the global velocity model. Centered on this auxiliary sparse model, a first-order Gaussian prior distribution around which the global velocity model fluctuates is constructed.

[0048] The mean of the first-order Gaussian prior distribution is taken as the auxiliary sparse model, and the variance is taken as the prior variance parameter. The first-order Gaussian prior distribution is: , The basic assumption of this first-order Gaussian prior distribution is that, under the condition of no travel-time observation data constraints, the global velocity model is most likely to take the velocity structure represented by the auxiliary sparse model, and the degree of deviation follows a Gaussian statistical law. Prior variance parameter The range of fluctuations around the auxiliary sparse model allowed by the global velocity model is quantitatively characterized.

[0049] Taking the negative logarithm of the first-order Gaussian prior distribution and discarding the constant term yields the global-auxiliary coupling constraint term.

[0050] Construction of the second-level conditional prior distribution: A second-level conditional prior constraint based on the sparse representation of the image patch dictionary is applied to the auxiliary sparse model. Specifically, the horizontal velocity distribution of the auxiliary sparse model at each depth layer is divided into several horizontally overlapping spatial image patches according to a preset spatial image patch side length. Each spatial image patch contains... The nth velocity grid point, and dictionary D is an overcomplete dictionary, are extracted from the auxiliary sparse model using an extraction operator. The velocity values ​​of a spatial image patch are given, and it is assumed that the spatial image patch can be approximated with high precision by a linear combination of a few atoms in the dictionary.

[0051] Assume that the sparse representation errors of each image patch are independent of each other, and that the errors follow a zero mean and a variance of . If the Gaussian distribution is given, then the second-order conditional prior distribution is the product of the independent Gaussian distributions of each spatial image patch.

[0052] Taking the negative logarithm of the second conditional prior distribution and discarding the constant term yields the image patch sparse representation error term.

[0053] In the joint inversion objective function, the travel-time data fitting term is used to calculate the sum of squares of the differences between the observed travel time and the sensitivity matrix multiplied by the velocity perturbation, and this term is weighted by the reciprocal of the travel-time observation error variance. The global-auxiliary coupling constraint term is used to calculate the sum of squares of the velocity differences between the global velocity model and the auxiliary sparse model at each grid node, and is weighted by the reciprocal of the prior variance parameter. The image patch sparse representation error term is used to calculate the sum of squares of the differences between the velocity values ​​of each spatial image patch and the dictionary atomic sparse linear combination approximation values, and accumulates them across all spatial image patches, weighted by the reciprocal of the image patch velocity approximation error variance. To eliminate dimensional differences and simplify hyperparameter tuning, the first term of the joint inversion objective function is normalized.

[0054] The joint inversion objective function is multiplied by the travel time observation error variance and normalized to make the coefficient of the travel time data fitting term normalized to 1, and a first tradeoff coefficient is defined. As a coefficient of the global-auxiliary coupling constraint term, it is used to control the weight between the travel time data fitting term and the global-auxiliary coupling constraint term, defining a second tradeoff coefficient. As a coefficient of the image patch sparse representation error term, it is used to control the weight between the travel time data fitting term and the image patch sparse representation error term.

[0055] In one embodiment, an alternating minimization strategy is used to solve the objective function, including: The joint optimization problem is broken down into the following three sub-problems; Solve iteratively, optimizing only one subproblem at a time while keeping the other two variables fixed: The first subproblem is the global velocity field update. With the auxiliary sparse model and sparse representation coefficients fixed, the LSMR algorithm is used to solve the augmented least squares equations with the sensitivity matrix as the core coefficient matrix, thus completing the update of the global velocity model. The second subproblem is the estimation of sparse representation coefficients. With the global velocity model and the auxiliary sparse model fixed, after the mean is removed from each spatial image patch, the sparse representation coefficients of each spatial image patch are solved one by one using the orthogonal matching pursuit algorithm. The third subproblem is solving the auxiliary sparse model. By fixing the global velocity model and the sparse representation coefficient matrix, the closed-form solution of the auxiliary sparse model is obtained using cyclic boundary conditions.

[0056] The joint inversion objective function is about the estimator The joint optimization is nonconvex. To estimate the global velocity model, To estimate the auxiliary sparse model, To estimate the sparse representation coefficients, we need a global velocity model, an auxiliary sparse model, and a sparse coefficient matrix. The coupling relationship between these three is as follows: the global velocity model appears in the first and second terms of the joint inversion objective function; the auxiliary sparse model appears in the second and third terms; and the sparse representation coefficients appear only in the third term. This coupling structure ensures that the joint optimization of the joint inversion objective function with respect to all variables is non-convex, but it is convex or sparsely solvable when optimizing each subproblem separately.

[0057] Therefore, in this embodiment, the joint optimization problem is broken down into three sub-problems, which are solved iteratively in sequence, optimizing only one sub-problem at a time while fixing the other two variables. Specifically, the first sub-problem is the global velocity field update, the second sub-problem is the sparse representation coefficient estimation, and the third sub-problem is the solution of the auxiliary sparse model. The three sub-problems are executed sequentially, constituting a complete internal iteration, which is repeated multiple times until convergence.

[0058] Therefore, an alternating minimization strategy is adopted to break it down into three sub-problems for step-by-step solution. Dictionary initialization uses an overcomplete DCT basis initialization. At that time, Laplace smoothing regularization is used to construct the initial velocity background; when When converting to dictionary-based sparse regularization, The number of warm-up iterations is preset based on the convergence results, ranging from 3 to 5.

[0059] The first subproblem is the global velocity field update, using a fixed auxiliary sparse model. With the sparse representation coefficients, the LSMR algorithm (least squared residual method) is used to solve the augmented least squares equations with the sensitivity matrix as the core coefficient matrix, thus updating the global velocity model. Fixed auxiliary sparse model. After the sparse representation coefficients, let Minimizing the first two terms of the joint inversion objective function for the global velocity model is equivalent to solving the following augmented least squares equations: , in, The update amount for the global velocity model is given by the LSMR algorithm in undamped mode, which completes the global velocity update as follows: .

[0060] When solving augmented least squares systems using the LSMR algorithm, the entire solution process proceeds in the following order: Construct the augmented system. Stack the sensitivity matrix and the weighted identity matrix vertically to form the augmented coefficient matrix; stack the travel-time residual vector and the weighted prior expectation update vertically to form the augmented observation vector. The size of this augmented system is the number of rows in the original sensitivity matrix plus the total number of model parameters, and the number of columns is the total number of model parameters.

[0061] Initialize the iteration variables. Set the initial solution to the zero vector, the initial residual to the augmented observation vector itself, and the initial search direction to the product of the identity matrix and the residual. Simultaneously, set the iteration counter to zero.

[0062] Enter the iterative loop. In each iteration, first calculate the product of the current augmented coefficient matrix and the current search direction, and then use this product to update the residual estimate and search direction. Specifically, through a series of orthogonal transformations and vector update operations, gradually construct the projection of the augmented coefficient matrix onto the Krylov subspace, and search for the solution that minimizes the sum of squared residuals in this subspace.

[0063] Update the current solution vector. Based on the information accumulated during the iteration process, update the estimated value of the model update amount. Due to the nature of the LSMR algorithm, each updated solution monotonically decreases the sum of squared residuals of the joint inversion objective function.

[0064] Check convergence criteria. Calculate the norm of the current residuals and determine if it has fallen below a preset tolerance threshold. Common convergence criteria include: the relative change in the residual norm is less than a preset threshold, the number of iterations reaches a preset upper limit, or the residual norm falls below the noise level of the observed data. If any convergence criterion is met, terminate the iteration; otherwise, increment the iteration counter and return to continue iteration.

[0065] The second subproblem is the estimation of sparse representation coefficients. With the global velocity model and auxiliary sparse model fixed, after removing the mean from each spatial image patch, the orthogonal matching pursuit algorithm is used to solve for the sparse representation coefficients of each spatial image patch independently. Before solving, mean centering is performed on each spatial image patch: , , in, Within a spatial image patch Average velocity of each pixel It is a vector of all 1s. This is the velocity vector of the image patch after the mean was removed.

[0066] Calculate the sparse representation coefficients for the mean-removed image patch: , , in, For the first Estimated sparse representation coefficients for spatial image patches.

[0067] The above equation is solved using the Orthogonal Matching Pursuit (OMP) algorithm, and the estimated sparse representation coefficient matrices of all spatial image patches are summarized and denoted as follows: , , The sparse coefficient matrix of the image patch Indicates the first The sparse representation coefficients of the _n spatial image patch. After restoring the mean, the _n The velocity of a spatial image block is approximately: .

[0068] The Orthogonal Matching Pursuit (OMP) algorithm iteratively selects the best-matching atoms from an overcomplete dictionary to sparsely represent spatial image patches. After each selection, the sparse representation coefficients and residuals are updated until a predetermined number of atoms are selected. In the sparse representation coefficient estimation, the norm distribution of the sensitivity matrix is ​​used to perform a ray coverage check on each spatial image patch. The proportion of pixels with zero ray density within the spatial image patch is statistically analyzed. When this proportion exceeds a preset threshold, the corresponding spatial image patch is considered to be in a region without effective ray constraints and is removed from the sparse representation coefficient estimation. During the patch-by-patch solution process of the OMP algorithm, the removed image patches are skipped; that is, they are not subjected to mean reduction, atom selection, coefficient updates, or residual calculations. This operation is linked to the alternating minimization solution process, and the removed image patches are not included in the calculation of the sparse coefficient estimation subproblem in subsequent outer iterations.

[0069] While estimating the sparse representation coefficients, the dictionary D is simultaneously optimized, allowing dictionary atoms to dynamically adapt to the local spatial characteristics of the current velocity model during the inversion process. Dictionary learning is based on the mean-reduced image patch velocity vectors. The process is performed using the ITKM algorithm. Before inversion initiation, the dictionary is initialized with an overcomplete DCT basis, and the atoms of the overcomplete DCT basis... The Each component Defined as: , Then, each atom was normalized. ).

[0070] set up For the first The support set of a spatial image patch (containing T atomic subscripts). For the support set in dictionary D The corresponding dictionary submatrix. The ITKM algorithm maximizes the atomic projections of all spatial image patches onto their respective support sets. Sum of norms: , in, To be the optimal dictionary To approximate the nonconvex optimization of the formula by the transpose of the dictionary submatrix corresponding to the support set, the following two alternating iterative steps are used, given the current dictionary. : Determine the support set: for each mean-reduced image patch velocity vector Find the way to project atoms Largest norm The set of atomic subscripts constitutes the support set: , in, This represents the number of iterations in the ITKM algorithm. For the current dictionary, supported by the set The transpose of the corresponding dictionary submatrix composed of atoms.

[0071] Atom update steps: For each atom In order to satisfy Summation is performed on spatial image patches in signed directions. atomic number Belongs to the A spatial image patch in the current dictionary The support set below And normalize it to a unit vector: , in, Indicates the first Atoms The value after iteration For the first Atoms The transpose of the iterated value.

[0072] During the inversion process, the dictionary is used at each interval The outermost layer is updated once using the ITKM algorithm, and each depth layer is updated once. Maintain a dictionary independently This allows for adaptation to the speed and structural characteristics of each layer. Practice shows that when the proportion of pixels with zero ray density within a spatial image patch exceeds 10%, that spatial image patch should be removed from dictionary training and sparse coefficient calculation.

[0073] The third subproblem is solving the auxiliary sparse model. With the global velocity model and sparse representation coefficient matrix fixed, a closed-form solution for the auxiliary sparse model is obtained using cyclic boundary conditions. During the solution process, information is shared between spatial image patches through overlapping regions. The velocity values ​​of the overlapping regions are simultaneously constrained by the dictionary reconstruction values ​​of multiple spatial image patches. After projecting the reconstruction values ​​of each spatial image patch back onto the complete grid using cyclic boundary conditions, the final velocity value of each pixel is equal to the weighted average of the reconstruction values ​​of all spatial image patches covering that pixel at that pixel location.

[0074] Specifically, this includes: the global velocity model Minimize the last two terms of the joint inversion objective function. Let Using wrap-around boundary conditions ensures that each pixel is exactly... The spatial image patch covers, that is, satisfies , To extract operators, To extract the projection operator of the operator, the closed-form solution of the auxiliary sparse model is obtained as follows: , when hour, Approaching the spatial image patch dictionary estimate, Indicates the estimation of the auxiliary sparse model; when This indicates that the contribution of sparse constraints is comparable to that of global constraints.

[0075] In the process of solving the sparse model, to eliminate boundary effects and ensure that each pixel in the entire study area receives uniform constraints, cyclic boundary conditions are used to segment the spatial image patches. Specifically, the study area is regarded as a ring-shaped topology in the horizontal direction. When a spatial image patch exceeds the right boundary, it continues to extend from the left; when it exceeds the lower boundary, it continues to extend from the top, so that pixels located at the edge of the study area have the same number of coverages as internal pixels. Through the above cyclic segmentation method, each horizontal position is covered by multiple spatial image patches, and adjacent spatial image patches form overlapping regions in the horizontal direction. The reconstructed values ​​of different spatial image patches in the overlapping region are derived from their respective selected dictionary atom combinations, and there are certain differences between these reconstructed values. The final pixel velocity value is determined by a weighted average: the final velocity value of each pixel is equal to the weighted average of the reconstructed values ​​of all spatial image patches covering that pixel at that pixel position, and the weights are determined by the overlap step size of the spatial image patches. In the central region of a spatial image patch, the contribution weight of that spatial image patch to the pixel velocity value is larger; in the region of the spatial image patch near the edge, the contribution weight gradually decreases. After the reconstructed values ​​of all spatial image patches are projected back to the complete grid, they are fused by weighted average to form an auxiliary velocity model with a final smooth transition.

[0076] The mathematical essence of this weighted average fusion process is as follows: by using cyclic boundary conditions, the sum of the self-products of the extraction operators is simplified to a constant multiplied by the identity matrix, thus representing the superposition of the reconstructed values ​​of each spatial image patch on the global grid as a uniform weighted average of the number of times each pixel is covered. This mechanism ensures that the final velocity value of each pixel is not determined solely by the dictionary reconstruction value of a single spatial image patch, but rather by a combination of independent estimates from all image patches covering that pixel. This overlapping region information-sharing mechanism eliminates stitching marks caused by image patch boundary division, giving the auxiliary velocity model spatial continuity. On the other hand, the common constraints of multiple spatial image patches suppress potential local errors or noise from a single image patch. Furthermore, since each spatial image patch independently solves for sparse coefficients, and information sharing in overlapping regions is achieved during the weighted averaging stage, this mechanism ensures model smoothness without affecting the computational efficiency of independent parallel solving of each image patch, making the entire auxiliary model solution process both accurate and efficient. When depth-direction continuity constraints are required, the information-sharing mechanism in the horizontal direction and the solution of the tridiagonal equations in the depth direction work together to achieve dual continuity and smoothness of the auxiliary velocity model in both the horizontal and vertical directions.

[0077] The aforementioned layer-by-layer independent solution does not involve continuity constraints in the depth direction. The closed-form solution of the auxiliary sparse model, without introducing depth direction regularization, is solved independently for each depth layer, with each layer's solution independent of the others. When depth direction regularization is introduced, the solutions for each depth layer are coupled through a tridiagonal linear equation system. When the depth coupling parameter... At this time, a first-order depth direction (Tikhonov) regularization term is introduced for each horizontal position in the joint inversion objective function: , in, For depth-direction regularization, These are the depth-direction regularization weight coefficients. This represents the total number of depth sampling points. Number the depth layers. For the first Auxiliary velocity values ​​for each depth layer For the first The auxiliary velocity value for each depth layer.

[0078] make , Horizontal position At the auxiliary velocity value, the layer-by-layer analytical solution of the estimated closed-form solution of the auxiliary sparse model is replaced with the solution at the horizontal position. Tridiagonal linear equations of order , It is a tridiagonal coefficient matrix. Horizontal position The right-hand term vector at position, where No. The components are: , in, The term vector on the right side is the first term. Each component.

[0079] Tridiagonal coefficient matrix for: , , in, The inner layer's main diagonal element. These are the main diagonal elements of the first and last layers.

[0080] The Thomas algorithm (also known as the catch-up algorithm) is used to solve the above tridiagonal linear equation system. The Thomas algorithm is specifically designed for solving tridiagonal linear equation systems. When, the tridiagonal coefficient matrix It degenerates into a diagonal matrix, at which point the result of solving layer by layer is completely equivalent to estimating the closed-form solution of the auxiliary sparse model.

[0081] The Thomas algorithm consists of two sequential stages: a forward elimination stage and a back-substitution stage. The forward elimination stage proceeds sequentially from the shallowest to the deepest layer, gradually eliminating the coupling coefficients on the lower diagonal and transforming the original tridiagonal equations into a simplified system containing only the main and upper diagonals. The back-substitution stage proceeds in reverse, from the deepest to the shallowest layer, calculating the auxiliary velocity values ​​for each depth layer.

[0082] The forward elimination phase starts from the first layer of the depth layer and proceeds sequentially to the last layer.

[0083] For the first layer, the intermediate coefficients of this layer are first calculated by dividing the negative value of the coupling coefficient between the first and second layers by the main diagonal element of the first layer. Simultaneously, the corrected right-hand side term of this layer is calculated by dividing the original right-hand side term of the first layer by the main diagonal element of the first layer.

[0084] For each intermediate layer, from the second layer to the penultimate layer, perform the following calculations layer by layer: First, calculate the denominator of the current layer, which equals the sum of the main diagonal elements of the current layer and the product of the coupling coefficient and intermediate coefficient of the previous layer. Then, divide the negative value of the coupling coefficient between the current layer and the next layer by the denominator to obtain the intermediate coefficient of the current layer. Finally, add the original right-hand side of the current layer to the product of the coupling coefficient and the modified right-hand side of the previous layer, and divide by the denominator to obtain the modified right-hand side of the current layer.

[0085] For the last layer, only its corrected right-hand side needs to be calculated: add the original right-hand side of the last layer to the product of the coupling coefficient of the previous layer and the corrected right-hand side of the previous layer, and then divide by the sum of the products of the main diagonal elements of the last layer, the coupling coefficient of the previous layer, and the intermediate coefficient of the previous layer.

[0086] After forward elimination is completed, the back substitution solution stage begins. This back substitution solution stage starts from the deepest layer and works backward to the shallowest layer.

[0087] First, the auxiliary velocity value of the last layer is directly equal to the right-hand side term of the last layer calculated in the forward elimination stage.

[0088] Then, starting from the second-to-last layer, the calculation proceeds upwards layer by layer to the first layer. For each layer, the auxiliary velocity value is equal to the right-hand side of the correction term for that layer minus the product of the intermediate coefficient of that layer and the auxiliary velocity value of the next layer.

[0089] From a physics perspective, the forward elimination stage is equivalent to propagating and accumulating the mutual influence between adjacent depth layers from shallow to deep. The corrected right-hand side of each layer includes not only its original right-hand side but also the cumulative influence of all layers above it. Intermediate coefficients record the effective coupling relationship between the current layer and the next. The back-substitution stage propagates known information backward from deep to shallow, determining the actual velocity values ​​of each layer. These two stages, one forward and one backward, comprehensively solve for the velocity values ​​of all depth layers.

[0090] The tridiagonal equations need to be solved independently for each horizontal grid position within the study area. Since the solutions for each horizontal position are completely independent, the Thomas algorithm can be executed iteratively for each horizontal position during the computation, or a parallel computation method can be used to solve the equations for all horizontal positions simultaneously.

[0091] Since most elements in the tridiagonal matrix are zero, the Thomas algorithm only processes the non-zero elements on the main diagonal and two adjacent diagonals. Therefore, its computational complexity is proportional to the number of depth layers. The forward elimination stage requires several multiplication and division operations for each layer, and the back substitution stage requires several multiplication and subtraction operations for each layer. The total computational complexity increases linearly with the number of depth sampling points, far exceeding the cubic complexity of ordinary Gaussian elimination.

[0092] The method described in this application is validated using a Rayleigh wave phase velocity dispersion dataset published for a specific basin. The dataset originates from background noise records from 448 stations, with observation periods ranging from 3 to 60 seconds. The stations are spatially unevenly distributed, resulting in significant differences in ray path coverage density within the study area. Figure 2 As shown, the white dashed line represents the main structural boundary, encompassing data acquired by both fixed and temporary seismic stations. The horizontal grid spacing for the study area is [value missing]. The vertical resolution is 5km for areas above 40km and 10km for areas below 40km, with the inversion depth range set at 0–60km. It should be noted that the actual data application in the aforementioned basin is merely an exemplary verification of the method described in this application and is not intended to limit its applicability. The method described in this application does not depend on specific regional geological conditions or station network layout. For other study areas, a three-dimensional velocity structure model of the corresponding area can be constructed simply by following the same grid partitioning strategy and inversion process based on the actual station distribution and dispersion observation data of that area.

[0093] To quantify the model reconstruction accuracy at each depth layer, for each depth layer... Calculate the root mean square error (RMSE) for each layer: , in For grid nodes The shear wave velocity value obtained from the inversion, For the meridional (X-direction) grid node numbering, Numbering of grid nodes in the latitudinal (Y direction) direction. Number the depth layers. To correspond to the actual shear wave velocity value of the theoretical model, Representing the The total number of valid grid points participating in the statistics for each layer. The set of valid grid points is determined as follows: ,in For the sensitivity matrix, the first The sum of the absolute values ​​of all elements in the column. This represents the maximum value within the global velocity model range. This is the sensitivity effectiveness threshold.

[0094] To verify the lateral resolution of the inversion scheme, a checkerboard synthesis test was designed and implemented. The input model consisted of alternating velocity perturbation blocks, each 0.6° x 0.6° in size, with relative velocity anomaly amplitude set to ±3%. Gaussian random noise with a standard deviation of 0.5% was superimposed on the synthesized travel time.

[0095] Both schemes employing the method described in this application use the same initial velocity model and mesh partitioning settings. Specifically, the smoothing regularization weight and damping coefficient of the comparison method (Direct Surface Wave Tomography) and the method described in this application are set to 7 and 16, respectively. The specific parameter configuration of the method described in this application is as follows: ,in, For dictionary size, This represents the upper limit of sparsity.

[0096] Figure 3 This is a map showing the restoration results of a checkerboard pattern test at a depth of 5km. Figure 3 (a) shows the phase velocity measurement and ray path distribution of Rayleigh waves with a period of 4–6 s; Figure 3 (b) and Figure 3 In this context, (f) represents the theoretical input model; Figure 3 (c) and Figure 3 (g) represents the comparison method and the inversion result of this application, respectively; Figure 3 (d) and Figure 3 (h) in the middle are respectively Figure 3 (c) and Figure 3 Enlarged view of the area shown in the pink dashed box in (g); Figure 3(e) represents the ray path density distribution. The black dashed line represents the main structural boundary, the triangles represent station locations, and the green solid line represents the profile. Location. At this shallow depth, the ray coverage is relatively sufficient, and both methods can reconstruct the overall shape of the checkerboard anomaly quite well. In comparison, the root mean square error of the method in this application is about 40% lower than that of the comparative method, and it also shows better amplitude recovery in the edge regions with sparse ray coverage. Figure 4 The results of the chessboard reconstruction at a depth of 40km are presented, in which... Figure 4 (a) shows the phase velocity measurement and ray path distribution of Rayleigh waves with a period of 39–41 s; Figure 4 (b) and Figure 4 In this context, (f) represents the theoretical input model; Figure 4 (c) and Figure 4 (g) represents the comparison method and the inversion result of this application, respectively; Figure 4 (d) and Figure 4 (h) in the middle are respectively Figure 4 (c) and Figure 4 Enlarged view of the area shown in the pink dashed box in (g); Figure 4 (e) represents the ray path density distribution. The black dashed line represents the main structural boundary, the triangles represent station locations, and the green solid line represents the profile. Location. As the ray coverage at this depth significantly decreases, the reconstruction quality of the contrast method deteriorates significantly, while the method of this application can still maintain the morphology of the checkerboard anomaly well, with a root mean square error reduction of about 38%, indicating that the method of this application can effectively improve the recovery accuracy of the anomaly under insufficient ray coverage.

[0097] Figure 5 The image presents the recovery results of the peak test at a depth of 35km, in which... Figure 5 (a) shows the phase velocity measurement and ray path distribution of Rayleigh waves with a period of 34–36 s; Figure 5 (b) and Figure 5 In this context, (f) represents the theoretical input model; Figure 5 (c) and Figure 5 In the figures (g), the inversion results are obtained by the comparative method and the method of this application, respectively. Figure 5 (d) and Figure 5 (h) in the middle are respectively Figure 5 (c) and Figure 5 Enlarged view of the area shown in the pink dashed box in (g); Figure 5 The (e) ray path density distribution is shown in the figure. The black dashed lines represent the main structural boundaries, the triangles represent station locations, and the green solid lines represent cross-sections. Location. The ray density is higher in the central part of the study area and relatively lower in the surrounding areas. In the sparsely covered areas (within the pink dashed box), the velocity recovery of the contrast method shows significant amplitude attenuation and tailing along the dominant ray direction, with a root mean square error of 58.00 m / s. In contrast, the method of this application has a significant suppression effect on tailing artifacts in the same area and maintains strong amplitude recovery capability, with the root mean square error reduced to 37.57 m / s, a reduction of about 35%, verifying that the method of this application has a substantial improvement in the recovery accuracy of isolated anomalies in underconstrained areas.

[0098] Figure 6 The results of the one-dimensional velocity recovery curves extracted along the AA′, BB′, and CC′ survey lines in the checkerboard test are presented. Figure 6 (a) in the figure is the one-dimensional velocity recovery curve of the AA′ survey line. Figure 6 In the figure, (b) represents the root mean square error of the AA′ survey line. Figure 6 (c) in the figure represents the one-dimensional velocity recovery curve of the BB′ survey line. Figure 6 In the figure, (d) represents the root mean square error of the BB′ survey line. Figure 6 (e) in the figure represents the one-dimensional velocity recovery curve of the CC′ survey line. Figure 6 In this context, (f) represents the root mean square error of the CC′ survey line at a depth of 5 km. On the cross-section, the one-dimensional velocity recovery curve of the method in this application is close to that of the one-dimensional velocity recovery curve of the actual velocity model. The root mean square error of the method in this application decreases from 30.78 m / s in the comparison method to 14.49 m / s, and the correlation coefficient increases from 0.92 to 0.98. (40 km depth) On the cross-section, the one-dimensional velocity recovery curve of the method in this application is close to that of the one-dimensional velocity recovery curve of the real velocity model. Although the overall accuracy decreases due to the reduced ray coverage, the method in this application still has stronger reconstruction capabilities, with the root mean square error reduced by 62% and the correlation coefficient increasing from 0.79 to 0.98. (Depth: 35 km) On the cross-section, the one-dimensional velocity recovery curve of the proposed method is close to that of the actual velocity model. The ray coverage in the edge region of the survey line is relatively sparse, and the comparative method exhibits significant attenuation of velocity amplitude and tailing distortion along the dominant direction in this area. The proposed method effectively suppresses the smearing artifacts in this region and has a better ability to recover the velocity amplitude of the anomaly. The root mean square error decreased from 58.00 m / s in the comparative method to 37.57 m / s, a reduction of approximately 35%. The results indicate that the regularization strategy proposed in this method can significantly improve the reconstruction accuracy of velocity anomalies in regions with insufficient ray constraints.

[0099] The above data shows that the method in this application has more accurate constraints on the amplitude and spatial location of velocity anomalies.

[0100] In checkerboard testing, adjacent perturbation blocks are tightly arranged, and the striped tails and morphological distortions of individual anomalies are often masked by the responses of surrounding anomalies. Therefore, discrete spike testing was further conducted to evaluate the ability of the two schemes to constrain the amplitude of isolated anomalies and suppress directional tail artifacts under non-uniform ray coverage conditions.

[0101] The spike test further validated the ability of both methods to identify isolated local velocity anomalies in a uniform background field.

[0102] To test the sensitivity of the proposed method to key regularization hyperparameters, parameter sensitivity experiments were conducted using a section checkerboard synthetic dataset. To conserve computational resources, the tests employed... The grid and over 20,000 dispersion observation data points were used, with the forward modeling structure and noise settings consistent with the complete synthesis experiment. The effects of sparsity, dictionary size, regularization parameters, and deep coupling coefficients on the inversion results were examined sequentially. When analyzing a particular parameter, all other parameters were fixed to the following baseline values: .

[0103] Figure 7 This diagram illustrates the impact of the sparsity ceiling T on the root mean square error (RMSE) in checkerboard recovery tests at different depths, including RMSE curves for 3km, 5km, 10km, 15km, 20km, 25km, 30km, 35km, and the mean. At that time, the fluctuation range of each layer and the root mean square error was relatively small, and the lowest root mean square error appeared at... Place; when At that time, the root mean square error increases with the upper limit of sparsity. The sparsity increases monotonically and its deeper effects are more pronounced, indicating that a moderate sparsity can balance the expressive power of image patches with the strength of sparsity constraints. An excessively high sparsity ceiling... On the contrary, it impairs the quality of recovery of abnormal bodies.

[0104] Figure 8 The influence of dictionary size Q on the root mean square error (RMSE) in checkerboard recovery tests at different depths is presented, and the RMSE is given as a function of dictionary size. The changing trends include the root mean square error (RMSE) curves for distances of 3km, 5km, 10km, 15km, 20km, 25km, 30km, and 35km, as well as the mean value. The RMSE varies with dictionary size. The increase shows an increasing trend. The lowest value was obtained at that location (32.83 m / s). The speed increased to 40.93 m / s, indicating that a smaller dictionary can effectively characterize local velocity features.

[0105] Figure 9 First trade-off coefficient Second tradeoff coefficient The sensitivity analysis results are shown in the figure. Among them, Figure 9 (a) in the figure is the distribution of the root mean square error as a function of parameter combinations. Figure 9 (a) shows the location of the optimal value; Figure 9 (b) in the equation represents the fixed second tradeoff coefficient. The root mean square error varies with λ1, with the second weighting factor fixed. The values ​​are 1, 4, 16, 64, and 256; Figure 9 (c) in the equation represents the fixed first tradeoff coefficient. Root mean square error varies with the second weighting factor Changes, fixing the first tradeoff coefficient The values ​​are 0.5, 2, 10, and 50. Within the current testing range, the minimum root mean square error is 32.78 ms / km, corresponding to... and .exist .

[0106] and Within the interval, the root mean square error variation does not exceed 10%, indicating that the inversion results have low overall sensitivity to both, and there is a relatively wide effective parameter range within the test range.

[0107] Figure 10 Depth coupling coefficient in chessboard recovery test at different depth layers The impact on root mean square error (RMSE) is shown, including RMSE variation curves for distances of 3km, 5km, 10km, 15km, 20km, 25km, 30km, and 35km, as well as the mean, with the baseline starting from 0. The deep coupling coefficient is also illustrated. The change in root mean square error at each depth layer as the value is gradually increased from 0 to 5. At this point, the algorithm degenerates into a layer-by-layer independent two-dimensional dictionary regularization scheme, which can be used as a benchmark. Under the current test conditions, the average root mean square error is... The velocity reaches a minimum of 32.11 m / s, a 4.8% reduction compared to the baseline, indicating that moderate interlayer coupling is beneficial for enhancing the vertical continuity of the velocity model without significantly compromising lateral resolution. At this time, excessive interlayer constraints will suppress lateral velocity variations, causing the average root mean square error to rise to 37.65 m / s. Overall, within the experimental framework, It can achieve a good balance between vertical continuity and lateral resolution.

[0108] To verify the adaptability of the dictionary adaptive learning process in this application to the local structural features of the velocity model, a 5km deep layer is used as an example to demonstrate the morphological evolution of dictionary atoms before and after learning by the ITKM algorithm.

[0109] Figure 11 The morphological evolution of dictionary atoms at a depth of 5km before and after learning by the ITKM algorithm is shown. Figure 10 (a) in the text represents the form before learning. Figure 10 (b) represents the learned form. Figure 10 (c) shows the morphological changes of atoms before and after learning, arranged neatly in rows and columns; the numbers in the figure represent the magnitude of these changes. The initial dictionary uses an overcomplete DCT (Discrete Cosine Transform) basis, with atoms arranged from low to high according to the overcomplete DCT basis, covering various spatial modes from smooth low-frequency to high-frequency oscillations. It can be seen that only one atom remains in a dormant state, while the remaining 15 have transitioned to an adapted state, becoming activated atoms. This indicates that the learned dictionary atoms have successfully adapted to the local structure of the current velocity model. Compared to the initial DCT atoms, the learned atoms exhibit more smooth gradients and local boundary features, with a reduction in high-frequency alternating modes, suggesting that the adaptive dictionary can more effectively characterize the local continuous and discontinuous structure of the velocity field in the study area.

[0110] The above description is merely a preferred embodiment of this application and is not intended to limit this application. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of this application should be included within the protection scope of this application.

Claims

1. A surface wave dispersion three-dimensional tomography method based on adaptive dictionary constraints, characterized in that, include: The study area is divided into three-dimensional meshes to construct a model parameter vector to characterize the three-dimensional shear wave velocity field. The model parameter vector consists of multiple shear wave velocity perturbations. The continuous slow field is discretized by ray tracing and bilinear interpolation to obtain the theoretical travel time. The observed travel time is calculated based on the observation data between station pairs. The travel time residual is constructed according to the theoretical travel time and the observed travel time. A linearized sensitivity matrix equation between the travel time residual and the shear wave velocity perturbation is constructed to obtain the sensitivity matrix and the norm distribution of the sensitivity matrix. Based on the sensitivity matrix, the likelihood function of the observed data is modeled as a Gaussian distribution. An auxiliary sparse model is introduced as a latent variable to establish a two-level Gaussian prior distribution, including a first-level Gaussian prior distribution and a second-level conditional prior distribution. The product of the likelihood function and the two-level Gaussian prior distribution constitutes the joint posterior probability density function. According to Bayes' theorem, maximizing the joint posterior probability density function is equivalently transformed into minimizing it, resulting in a joint inversion objective function composed of the weighted sum of three terms: the travel time data fitting term, the global-auxiliary coupling constraint term, and the image patch sparse representation error term. In the early stage of inversion, a warm-up iteration is performed, and the initial velocity background is constructed using Laplace smoothing regularization. After the warm-up, the dictionary sparse regularization is switched to, and the joint inversion objective function is solved using an alternating minimization strategy. After iterative convergence, the final three-dimensional shear wave velocity structure model is output.

2. The surface wave dispersion three-dimensional tomography method based on adaptive dictionary constraints according to claim 1, characterized in that, When the likelihood function of the observed data is modeled as a Gaussian distribution, the mean of the Gaussian distribution is the sensitivity matrix multiplied by the difference between the global velocity model and the current velocity model. The difference represents the predicted travel time change value generated by the velocity disturbance of the current velocity model relative to the linearized expansion point after being mapped by the sensitivity matrix. The variance of the Gaussian distribution is the travel time observation error variance, which is used to characterize the noise level in the observed data. The weight of the likelihood function is controlled by the reciprocal of the travel time observation error variance.

3. The surface wave dispersion three-dimensional tomography method based on adaptive dictionary constraints according to claim 1, characterized in that, The first-level Gaussian prior distribution in the two-level Gaussian prior distribution is a Gaussian distribution centered on the auxiliary sparse model for the global velocity model. The mean of the first-level Gaussian prior distribution is taken from the auxiliary sparse model, and the variance is taken from the prior variance parameter. The prior variance parameter quantitatively characterizes the range within which the global velocity model is allowed to fluctuate around the auxiliary sparse model. Based on the value of the prior variance parameter, it is used to select whether the construction boundary and velocity discontinuity features carried in the auxiliary sparse model are completely transmitted to the global velocity model, or whether the global velocity model deviates from the auxiliary sparse model.

4. The surface wave dispersion three-dimensional tomography method based on adaptive dictionary constraints according to claim 1, characterized in that, The second-level conditional prior distribution is the conditional prior distribution of the auxiliary sparse model based on the sparse representation of the image patch dictionary, and it is constructed as follows: The lateral velocity distribution of the auxiliary sparse model at each depth layer is divided into several overlapping spatial image patches, each containing a fixed number of velocity grid points. The velocity value of each spatial image patch is extracted from the auxiliary sparse model using an extraction operator; Each spatial image block is approximated by a sparse linear superposition of atoms in an overcomplete dictionary. Each spatial image block corresponds to a sparse representation coefficient, and the number of non-zero elements in the sparse representation coefficient must not exceed the preset sparsity limit. The sparse representation coefficients of each spatial image patch are statistically independent of each other; Using the dictionary sparse approximation error of each spatial image patch as the mean and the variance of the velocity approximation error of the spatial image patch as the variance, the second-level conditional prior distribution of the auxiliary sparse model is constructed as the product of the independent Gaussian distributions of each spatial image patch.

5. The surface wave dispersion three-dimensional tomography method based on adaptive dictionary constraints according to claim 4, characterized in that, The travel time data fitting term is used to calculate the sum of squares of the difference between the observed travel time and the sensitivity matrix multiplied by the velocity disturbance. The travel time data fitting term is weighted by the reciprocal of the variance of the travel time observation error. The global-auxiliary coupling constraint term is used to calculate the sum of squares of the velocity differences between the global velocity model and the auxiliary sparse model at each grid node, and is weighted by the reciprocal of the prior variance parameter. The image patch sparse representation error term is used to calculate the sum of squares of the differences between the velocity values ​​of each spatial image patch and the approximations of the sparse linear combination of dictionary atoms, and is accumulated across all spatial image patches. It is also weighted by the reciprocal of the variance of the spatial image patch velocity approximation error.

6. The surface wave dispersion three-dimensional tomography method based on adaptive dictionary constraints according to claim 5, characterized in that, The joint inversion objective function is multiplied by the travel time observation error variance and normalized to normalize the coefficient of the travel time data fitting term to 1. A first tradeoff coefficient is defined as the coefficient of the global-auxiliary coupling constraint term to control the weight between the travel time data fitting term and the global-auxiliary coupling constraint term. A second tradeoff coefficient is defined as the coefficient of the image patch sparse representation error term to control the weight between the travel time data fitting term and the image patch sparse representation error term.

7. The surface wave dispersion three-dimensional tomography method based on adaptive dictionary constraints according to claim 4, characterized in that, An alternating minimization strategy is used to solve the joint inversion objective function, including: The joint optimization problem is broken down into three subproblems, which are solved iteratively, optimizing only one subproblem at a time while keeping the other two variables fixed: The first subproblem is the global velocity field update. With the auxiliary sparse model and sparse representation coefficients fixed, the LSMR algorithm is used to solve the augmented least squares equations with the sensitivity matrix as the core coefficient matrix, thus completing the update of the global velocity model. The second subproblem is the estimation of sparse representation coefficients. With the global velocity model and the auxiliary sparse model fixed, after the mean is removed from each spatial image patch, the sparse representation coefficients of each spatial image patch are solved one by one using the orthogonal matching pursuit algorithm. The third subproblem is solving the auxiliary sparse model. By fixing the global velocity model and the sparse representation coefficient matrix, the closed-form solution of the auxiliary sparse model is obtained using cyclic boundary conditions.

8. The surface wave dispersion three-dimensional tomography method based on adaptive dictionary constraints according to claim 7, characterized in that, In the sparse representation coefficient estimation, the column norm distribution of the sensitivity matrix is ​​used to perform ray coverage checks on each spatial image block, and the proportion of pixels with zero ray density within the spatial image block is counted. When the proportion of pixels exceeds a preset threshold, the corresponding spatial image block is considered to be in a region without effective ray constraints and is removed from the sparse representation coefficient estimation.

9. The surface wave dispersion three-dimensional tomography method based on adaptive dictionary constraints according to claim 7, characterized in that, During the solution process of the auxiliary sparse model, information sharing is achieved between spatial image blocks through overlapping regions. The velocity value of the overlapping part is simultaneously constrained by the dictionary reconstruction values ​​of multiple spatial image blocks. After projecting the reconstruction values ​​of each spatial image block back to the complete grid using cyclic boundary conditions, the final velocity value of each pixel is equal to the weighted average of the reconstruction values ​​of all spatial image blocks covering that pixel at that pixel position.

10. The surface wave dispersion three-dimensional tomography method based on adaptive dictionary constraints according to claim 7, characterized in that, The closed-form solution of the auxiliary sparse model is solved independently for each depth layer without introducing depth direction regularization, and the solution of each depth layer does not affect each other; when depth direction regularization is introduced, the solution of each depth layer is coupled through a tridiagonal linear equation system.

Images