Multi-geophysical field joint inversion method based on physical property subspace decoupling

Abstract

The invention discloses a multi-geophysical field joint inversion method based on physical subspace decoupling, and belongs to the field of geophysical exploration. According to the invention, the problem of poor inversion effect caused by difficult characterization of multivariate and nonlinear physical relations of different types of rocks under complex geological conditions in the existing structure coupling joint inversion method is solved; according to the method, K subspace clustering is utilized to dynamically divide underground model units into different low-rank subspaces according to physical property parameter characteristics in each inversion iteration, and then a low-rank constraint based on a Schatten-p norm is applied to a physical property parameter matrix in each subspace; the clustering-coupling process and the inversion iteration are alternately carried out, so that the division of the physical property subspace and the model updating in each subspace are mutually promoted, the self-adaptive optimization is realized, and the clustering-coupling process and the inversion iteration are carried out at the same time; and finally, high-precision and high-resolution reconstruction of multi-type rock physical property distribution in the complex underground structure is realized.

Description

Technical Field

[0001] This invention relates to the field of geophysical exploration technology, specifically to a multi-geophysical field joint inversion method based on physical property subspace decoupling. Background Technology

[0002] Geophysical inversion often suffers from the problem of non-uniqueness of solutions. To improve detection accuracy, joint inversion methods are frequently employed, utilizing different physical properties (such as density and magnetism) to jointly characterize subsurface structures. Joint inversion mainly relies on coupled constraints of rock properties or coupled constraints of model structure. The former is limited by empirical formulas for physical properties that are difficult to apply universally, while the latter reduces ambiguity by seeking structural similarities between different physical property models, and is currently the more widely used approach.

[0003] In structural coupling constraint methods, cross-gradient and Gramian constraints are two commonly used techniques. Cross-gradient requires that the gradient directions of different physical property models be parallel, while Gramian constraints aim to promote the linear correlation of multiple physical property vectors in the overall space. However, when the actual underground conditions are complex and the relationship between the multiple physical property parameters of different rock types is not a single linear trend, the effectiveness of these traditional structural coupling methods will be limited.

[0004] The Earth's interior is a complex system. Different types of rocks, under the influence of different diagenetic environments and tectonic movements, develop their own unique combinations of multiple physical parameters, belonging to different "physical property subspaces." If these different physical property subspaces can be identified and decoupled, and the correlation of physical properties within each subspace can be enhanced during inversion, it is hoped that complex geological structures can be characterized more precisely. In this context, low-rank approximation and subspace clustering techniques (such as K-plane clustering) provide effective tools for identifying and decoupling these different low-rank subspaces from high-dimensional physical property data; therefore, we propose a multi-geophysical field joint inversion method based on physical property subspace decoupling. Summary of the Invention

[0005] The purpose of this invention is to provide a joint inversion method for multiple geophysical fields based on property subspace decoupling. By introducing the K-subspace clustering algorithm to adaptively decouple the multiple property parameter spaces, and constructing clustered Schatten-p low-rank constraint terms to enhance the structural correlation within each property subspace, the problem mentioned in the background art is solved.

[0006] To achieve the above objectives, the present invention provides the following technical solution: a multi-geophysical field joint inversion method based on property subspace decoupling, comprising the following steps:

[0007] Organize gravity and magnetic measurement data, and divide the inversion region into... A block arranged according to a rule;

[0008] Data matrix, sensitivity matrix and data standard deviation matrix are constructed based on gravity observation data and magnetic measurement observation data, respectively.

[0009] Construct a depth weighting matrix and set inversion parameters including the number of subspaces, cluster low-rank constraint weighting coefficients, depth constraint weighting coefficients, upper and lower bounds of physical property parameters, and Schatten-p norm parameters;

[0010] Initialize the density model, magnetization model, auxiliary iteration variables, and Nesterov momentum parameters;

[0011] K-subspace clustering is performed on the standardized property matrix to obtain the subspace labels of each block;

[0012] Calculate the gradient of the density versus magnetization model and update the auxiliary iteration variables;

[0013] Based on the clustering results, singular value decomposition is performed on the physical property parameter matrix in each subspace, and clustered Schatten-p norm near-end mapping is applied for update.

[0014] The updated physical property parameters are subjected to box-constrained projection, and the updated density and magnetization parameters are projected to the preset upper and lower bound intervals respectively through box projection.

[0015] The Nesterov acceleration strategy is used to update the auxiliary iteration variables for the next step;

[0016] Determine if the convergence condition is met. If it is, output the inversion result; otherwise, return to the K-subspace clustering step to continue iterating.

[0017] Furthermore, based on gravity observation data and magnetic measurement observation data, a data matrix, a sensitivity matrix, and a data standard deviation matrix are constructed, including:

[0018] Data matrices were generated using gravity and magnetic field observation data respectively. and The underground grid is divided and gravity and magnetic susceptibility matrices are generated separately. and Data standard deviation matrices were generated using gravity and magnetic data respectively. and .

[0019] Furthermore, in the step of constructing the depth-weighted matrix, the formula is used. Constructing a depth-weighted matrix and ;

[0020] in, Indicates the first The depth of each block; Represents a constant, usually . ; This represents the weighted decay index.

[0021] Furthermore, the density and magnetization models, auxiliary iteration parameters, and momentum parameters are initialized, including:

[0022] Set the initial density model Initial model of magnetization Density model auxiliary iteration quantity Magnetization intensity model auxiliary iteration quantity Initial values ​​of Nesterov momentum parameters and the current iteration number .

[0023] Furthermore, K-subspace clustering is performed on the standardized property matrix, including:

[0024] Construct the first The normalized property matrix of the next iteration ;

[0025] In the formula, Indicates the first The standardized property matrix constructed in the next inversion iteration; This represents the total number of grid cells after the inverted region is subdivided. The property matrix is ​​a A real number matrix with 2 rows and 2 columns, where each row represents a grid cell and each column represents a standardized physical property parameter;

[0026] For any column vector in the property matrix Standardized computing includes ;

[0027] in, , indicates the calculation of column vectors The arithmetic mean;

[0028] , indicates the calculation of column vectors The sample standard deviation;

[0029] Represents the first element in the standardized vector. The value of each element;

[0030] by The K-subspace clustering algorithm is used to generate a clustering index set for the input data. ;

[0031] in, The The element represents the element. The iteration of the ... Subspace labels with multiple physical parameters.

[0032] Furthermore, the K-subspace clustering steps include:

[0033] The density and magnetization auxiliary iteration variables are standardized to construct the property matrix;

[0034] Set the number of subspaces, the dimension of the subspaces, and the convergence threshold;

[0035] Initialize the subspace labels of each block and calculate the initial basis vectors of each subspace;

[0036] Iteratively update the subspace labels of each block and recalculate the basis vectors of each subspace;

[0037] The iteration stops when the labels no longer change or the rate of change of the objective function is less than the threshold, and the subspace index set is output.

[0038] Furthermore, the clustered Schatten-p norm near-end mapping update steps include:

[0039] The physical property parameter matrix is ​​divided into several sub-matrices based on the clustering index;

[0040] Perform singular value decomposition on each submatrix;

[0041] Solve the Schatten-p norm near-end mapping problem for each element in the singular value vector;

[0042] Reconstruct the submatrix using the updated singular values ​​and update the physical property parameter values ​​at the corresponding indices.

[0043] The Schatten-p norm near-end mapping problem is a scalar optimization problem, as detailed below:

[0044]

[0045] In the formula, These are the original singular values; For regularization parameters;

[0046] Furthermore, in the Nesterov acceleration step, the auxiliary iteration variables for the next step are calculated based on the changes in momentum and physical property parameters of the current iteration step.

[0047] Furthermore, the stopping criterion is that the L2 norm of the gradient vector is less than a preset tolerance, that is: .

[0048] Furthermore, the inversion parameters include upper and lower bounds for density constraints and upper and lower bounds for magnetization constraints, which are used to limit the reasonable range of physical property parameters in the box projection step.

[0049] Compared with the prior art, the beneficial effects of the present invention are:

[0050] 1. This invention introduces a clustered Schatten-p low-rank constraint based on K-subspace clustering into the joint inversion objective function, thereby achieving adaptive decoupling and coupling of the physical property relationships of multiple rock types in complex underground spaces. This overcomes the limitations of traditional cross-gradient or Gramian constraints that emphasize overall linear correlation, and can effectively characterize the specific coupling relationships of different rock types in their respective physical property subspaces, significantly improving the inversion resolution and reliability under complex geological structures.

[0051] 2. This invention combines property subspace clustering, low-rank constraint update of clustering, and Nesterov accelerated optimization strategy to form an efficient adaptive joint inversion framework. This framework not only dynamically optimizes the structural coupling relationship through iterative update of the property subspace, but also uses accelerated optimization algorithm to ensure convergence speed. It solves the problems of poor adaptability and low computational efficiency of traditional joint inversion methods in complex multi-rock type scenarios, thus providing a powerful technical means for deep mineral resource exploration and crustal structure research. Attached Figure Description

[0052] Figure 1 This is a flowchart of the multi-geophysical field joint inversion method based on property subspace decoupling of the present invention;

[0053] Figure 2 This is a schematic diagram of the joint inversion framework based on the decoupling and coupling of the property subspace of the present invention;

[0054] Figure 3 This is a schematic diagram of the theoretical model experiment of the present invention;

[0055] Wherein, (a) is the theoretical model, (b) and (c) are the forward gravity and magnetic anomalies with noise, respectively, (d) and (e) are the density and magnetization results of independent inversion, respectively, (f) and (g) are the density and magnetization results of inversion by the method in this paper, respectively, and (h) and (i) are the physical property distributions of the results of independent inversion and inversion by the method in this paper, respectively.

[0056] Figure 4 This is a schematic diagram of the Gramian-constrained joint inversion results of the theoretical model of this invention;

[0057] Wherein, (a), (b), and (c) are the density model, magnetization model, and spatial distribution of physical properties when the Gramian weight is 0.001, respectively; (d), (f), and (e) are the spatial distributions of density, magnetization, and physical properties when the Gramian weight is 0.01; (g), (h), and (i) are the spatial distributions of density, magnetization, and physical properties when the Gramian weight is 0.1; and (j), (k), and (l) are the spatial distributions of density, magnetization, and physical properties when the Gramian weight is 1. Detailed Implementation

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

[0059] To address the technical problem of poor inversion results in existing structural coupling joint inversion methods under complex geological conditions due to the difficulty in characterizing the multivariate and nonlinear physical property relationships of different rock types, please refer to [link to relevant documentation]. Figures 1-4 This embodiment provides the following technical solution:

[0060] The multi-geophysical field joint inversion method based on property subspace decoupling includes the following steps:

[0061] Organize gravity and magnetic measurement data, and divide the inversion region into... A block arranged according to a rule;

[0062] Data matrix, sensitivity matrix and data standard deviation matrix are constructed based on gravity observation data and magnetic measurement observation data, respectively.

[0063] Construct a depth-weighted matrix and set the number of subspaces. Clustering low-rank constraint weighting coefficients 0. Inversion parameters including depth constraint weighting coefficients, upper and lower bounds of physical property parameters, and Schatten-p norm parameters. Among them, the depth constraint weighting coefficients are divided into gravity depth constraint weighting coefficients. 0. Magnetic depth constraint weighting coefficient 0;

[0064] Initialize the density model, magnetization model, auxiliary iteration variables, and Nesterov momentum parameters;

[0065] K-subspace clustering is performed on the standardized property matrix to obtain the subspace labels of each block;

[0066] Calculate the gradient of the density versus magnetization model and update the auxiliary iteration variables. The gradient calculation is shown below:

[0067]

[0068] ;

[0069] Based on the clustering results, singular value decomposition is performed on the physical property parameter matrices in each subspace, and a clustered Schatten-p norm near-end mapping update is applied. The parameter matrix update is shown below:

[0070]

[0071] ;

[0072] The updated physical property parameters are subjected to box-constrained projection. The updated density and magnetization parameters are projected onto preset upper and lower bound intervals using box projection. The box projection is... and ;

[0073] Nest The ERov acceleration strategy updates the auxiliary iteration variables for the next step;

[0074] Determine if the convergence condition is met; if so, output the inversion result, i.e., the output density and magnetization model updated in the current iteration. and Otherwise, return to the K-subspace clustering step and continue iterating, where, and This is the model from the previous iteration.

[0075] The technical effects of the above solution are as follows: By introducing K-subspace clustering and Schatten-p low-rank constraints, automatic decoupling and coupling of complex underground multi-physical parameter spaces are achieved, thereby effectively overcoming the limitations of traditional cross-gradient or Gramian constraint methods in nonlinear, multi-trend physical property relationship scenarios. Furthermore, this method can adaptively identify physical property subspaces corresponding to different lithologies and enhance the correlation of multi-physical parameters in each subspace through low-rank approximation, thereby significantly improving the resolution and geological interpretation rationality of the joint inversion model. At the same time, by combining Nesterov acceleration and box-constrained projection techniques, computational efficiency can be improved while ensuring stable convergence of the inversion process, ultimately obtaining a multi-physical parameter model that is more consistent with actual geological conditions.

[0076] Based on gravity and magnetic measurement data, a data matrix, sensitivity matrix, and data standard deviation matrix were constructed, including:

[0077] Data matrices were generated using gravity and magnetic field observation data respectively. and The underground grid is divided and gravity and magnetic susceptibility matrices are generated separately. and Data standard deviation matrices were generated using gravity and magnetic data respectively. and .

[0078] The technical effects of the above-mentioned technical solution are as follows: by constructing data matrices, sensitivity matrices, and data standard deviation matrices for gravity and magnetic measurement data respectively, the error distribution of different observation data can be effectively quantified, and the physical correlation between surface observations and subsurface physical parameters can be accurately established through the sensitivity matrix. This achieves an independent and unified mathematical representation of multi-source geophysical data, thus laying an accurate data foundation for subsequent joint inversion.

[0079] In the step of constructing the depth-weighted matrix, the formula is used. Constructing a depth-weighted matrix and ;

[0080] in, Indicates the first The depth of each block; Represents a constant, usually . ; This represents the weighted decay index.

[0081] The technical effect of the above solution is as follows: by using a weighting function with depth attenuation characteristics to construct a depth weighting matrix, the ability to detect deep anomalies can be significantly enhanced, making the inversion results more reasonably distributed in the vertical direction, thereby obtaining a more geologically significant underground physical parameter model.

[0082] Initialize the density and magnetization model, auxiliary iteration parameters, and momentum parameters, including:

[0083] Set the initial density model Initial model of magnetization Density model auxiliary iteration quantity Magnetization intensity model auxiliary iteration quantity Initial values ​​of Nesterov momentum parameters and the current iteration number .

[0084] The technical effect of the above solution is that by initializing the density and magnetization model and its auxiliary iterations to zero, and setting the Nesterov momentum parameter to a standard initial value, the problem of local optimal solutions caused by improper selection of the initial model can be effectively avoided.

[0085] K-subspace clustering is performed on the standardized property matrix, including:

[0086] Construct the first The normalized property matrix of the next iteration ;

[0087] In the formula, Indicates the first The standardized property matrix constructed in the next inversion iteration; This represents the total number of grid cells after the inverted region is subdivided. The property matrix is ​​a A real number matrix with 2 rows and 2 columns, where each row represents a grid cell and each column represents a standardized physical property parameter;

[0088] For any column vector in the property matrix Standardized computing includes ;

[0089] in, , indicates the calculation of column vectors The arithmetic mean;

[0090] , indicates the calculation of column vectors The sample standard deviation;

[0091] Represents the first element in the standardized vector. The value of each element;

[0092] by The K-subspace clustering algorithm is used to generate a clustering index set for the input data. ;

[0093] in, The The element represents the element. The iteration of the ... Subspace labels with multiple physical parameters.

[0094] The technical advantages of the above solution are as follows: By standardizing the auxiliary iterative variables of density and magnetization and constructing a property matrix, and then using the K-subspace clustering algorithm to dynamically identify and decouple the low-rank property subspaces corresponding to different lithologies, the nonlinear coupling relationship of multiple property parameters in complex underground environments can be adaptively characterized. This effectively overcomes the dependence of traditional joint inversion methods on global linear correlation. Specifically, the standardization step calculates the arithmetic mean of each property parameter. and sample standard deviation This can eliminate dimensional differences and highlight structural features, while the index set generated by clustering can divide each grid cell into a specific subspace label, thus laying a solid foundation for the subsequent implementation of clustered low-rank constraints, significantly improving the joint inversion model's ability to identify multi-type lithological structures and the rationality of geological interpretation.

[0095] The K-subspace clustering steps include:

[0096] The density and magnetization auxiliary iteration variables are standardized to construct the property matrix;

[0097] Set the number of subspaces The subspace dimension and convergence threshold, where the subspace rank is... and 0;

[0098] Initialize the subspace labels of each block and calculate the initial basis vectors of each subspace;

[0099] Iteratively update the subspace labels of each block and recalculate the basis vectors of each subspace;

[0100] The iteration stops when the labels no longer change or the rate of change of the objective function is less than the threshold, and the subspace index set is output.

[0101] In this embodiment, the initial label For each , This represents the number of samples in cluster k during initialization; the samples of this cluster are stacked row-wise to obtain the row-stacked data matrix during initialization. Use the first r eigenvectors to set the basis of the subspace: and ,in This represents the expression that returns the first r orthogonal eigenvectors of its symmetric independent variable. Matrix operators; setting ;

[0102] In this embodiment, label assignment: for each Distribute using the following formula ;

[0103] In this embodiment, for each , build ,in It is the number of samples in the k-th cluster in this round; then let and ;

[0104] In this embodiment, stop iterative discrimination: Let ,like or If yes, then stop; otherwise, let Then return to the K subspace clustering step to continue iterating;

[0105] In this embodiment, the output subspace index set is... ,in .

[0106] The technical effects of the above solution are as follows: by systematically executing the complete process of standardization, parameter setting, initialization, iterative optimization and convergence judgment, the automatic low-dimensional subspace partitioning of the high-dimensional physical property parameter space is realized, thereby enabling adaptive identification and decoupling of physical property structure patterns corresponding to different lithologies. This provides an accurate clustering basis for subsequent implementation of clustering low-rank constraints, effectively overcomes the dependence of traditional joint inversion methods on global linear correlation, and improves the identification accuracy of multi-type lithological structures and the reliability of inversion results under complex geological conditions.

[0107] The steps for updating the clustered Schatten-p norm near-end mapping include:

[0108] The physical property parameter matrix is ​​divided into several sub-matrices based on the clustering index;

[0109] Perform singular value decomposition on each submatrix;

[0110] Solve the Schatten-p norm near-end mapping problem for each element in the singular value vector;

[0111] Reconstruct the submatrix using the updated singular values ​​and update the physical property parameter values ​​at the corresponding indices.

[0112] Among them, let For each type calculate , and have ;

[0113] Here and They respectively represent to put and Limited to the index set The sub-vector obtained above;

[0114] Each singular value is updated by solving the following scalar proximal problem. ;

[0115] The Schatten-p norm near-end mapping problem is a scalar optimization problem, as detailed below:

[0116]

[0117] In the formula, These are the original singular values; For regularization parameters; ;

[0118] And reconstruct using the following formula , and .

[0119] The technical effects of the above solution are as follows: by dividing the physical property parameter matrix into several sub-matrices according to the clustering index, and performing singular value decomposition on each sub-matrix, the Schatten-p norm near-end mapping problem is solved for each element in the singular value vector. The scalar optimization problem is used to achieve adaptive compression of the singular values. Finally, the updated singular values ​​are used to reconstruct the sub-matrices and update the corresponding physical property parameter values. This effectively suppresses noise interference while preserving the main structural features of each subspace, significantly enhancing the stability of the joint inversion results and the ability to identify different lithological physical property structures.

[0120] In the Nesterov acceleration step, based on the changes in momentum and physical property parameters in the current iteration step, auxiliary iteration variables for the next step are calculated, including the calculation of the dynamic update formula for the momentum parameter. ,

[0121] Extrapolation update formula for auxiliary iterative variables and .

[0122] The technical effects of the above solution are as follows: By using the momentum update mechanism designed in the Nesterov acceleration step, the dynamic update formula of the momentum parameter and the extrapolation update formula of the auxiliary iteration variable are used to effectively accelerate the convergence process of the joint inversion algorithm, and significantly improve the inversion calculation efficiency while maintaining the stability of the algorithm.

[0123] The stopping condition for iteration is that the L2 norm of the gradient vector is less than a preset tolerance, i.e.: Otherwise 1.

[0124] The technical effect of the above solution is as follows: by setting the condition that the gradient vector L2 norm is less than the preset tolerance as the criterion for stopping the iteration, a clear and reliable convergence criterion can be provided for the joint inversion process, thereby enabling adaptive judgment of the optimization degree of the inversion solution and effectively avoiding under-convergence or over-convergence problems.

[0125] The inversion parameters include upper and lower bounds for density constraints and magnetization constraints, used to limit the reasonable range of physical property parameters during the box projection step. The upper bound for density constraints is... The lower bound of the density constraint is The upper bound of the magnetization constraint is The lower bound of the magnetization constraint is .

[0126] The technical effects of the above solution are as follows: by setting upper and lower bounds for density constraints and upper and lower bounds for magnetization constraints as inversion parameters, and by forcing the physical property parameters to fall within a preset reasonable range in the box projection step, the inversion process is effectively constrained by prior geological knowledge. This not only ensures the geophysical rationality of the inversion results, but also improves the numerical stability and convergence efficiency of the algorithm.

[0127] Specifically, this embodiment also proposes a framework for the inversion method (such as...). Figure 2 As shown), specifically:

[0128] Divide the underground space into A multi-geophysical model is established by inverting a set of regularly arranged blocks to obtain multiple geophysical parameters for each block. For geophysical inversion of gravity and magnetic data, the density and magnetization of the i-th block can be expressed as... and Then a property matrix can be constructed.

[0129] (1)

[0130] To enhance the correlation between different rock types in their respective physical property subspaces, we propose a joint inversion optimization model as follows:

[0131] (2)

[0132] in,

[0133] (3)

[0134] and

[0135] (4)

[0136] These are the data fitting term and the model constraint term, respectively. , , , and These are the gravity data vector, sensitivity matrix, data weighting matrix, depth weighting term, and weight parameters. Correspondingly, , , , and These are the magnetic data vector, sensitivity matrix, data weighting matrix, and weight parameters. Sensitivity matrix and It is related to the spatial location of the sampling point and the rectangular block. and The j-th diagonal element is the reciprocal of the standard deviation of the j-th data element. and The j-th diagonal element is ,in and This represents the depth and depth-weighted parameters of the corresponding rectangular block, and and Recommended for and The third term in formula (2) is a coupling constraint term for multiple rock types, where... Low-rank structural coupling constraint terms representing a single property space

[0137] (5)

[0138] in Indicates belonging to the k-th property subspace Sub-data matrices in Describing the Schatten-p norm ( ), whose value is equal to the Lp norm of the singular values ​​of the matrix. Furthermore, in formula (2), K is the weighting parameter, and K represents the number of property subspaces.

[0139] After determining the objective function (2), it is assumed that the data lies on multiple hyperplanes with different trends in space, each representing a low-rank subspace. Subspace clustering methods can be used to identify and decouple data located in different low-rank subspaces. K-plane clustering is one of the commonly used low-rank approximation methods for subspace clustering. This method achieves clustering of high-dimensional data by iteratively fitting the low-dimensional linear subspace of each cluster with PCA and redistributing the samples to the subspace with the smallest projection error.

[0140] K-plane clustering minimizes points with multiple physical property parameters. With the material property subspace The total projection error realizes the material property subspace structure, and its objective function can be expressed as:

[0141] (6)

[0142] in express To the material property subspace The distance. The calculation process of formula (6) can be achieved through alternating iterations: First, each sample point is assigned to the nearest subspace; then, the point set belonging to the same subspace is updated, and the mean position and principal direction of the subspace are re-estimated, thus obtaining a new subspace representation. The above assignment and update process is carried out alternately until the sample partitioning and subspace structure converge. In the above process, the point label assignment follows the following formula:

[0143] (7)

[0144] in It is the center point of the k-th property subspace. It is the basis of the k-th property subspace, obtained from principal component analysis of the k-th property subspace data:

[0145] (8)

[0146] in yes The first d columns. Therefore, the k-th property subspace can be regarded as the first d columns. Zhang Cheng's d-dimensional space. It should be noted that the dimension of the property subspace is smaller than the dimension of the common property space, but its dimension can be set in advance. For example, when the value of d is set to 1, classification is performed according to a one-dimensional linear trend; when the value of d is set to 2, the property subspace can be divided according to a two-dimensional planar trend. Clearly, different dimensions of property subspaces can be set according to the actual situation.

[0147] Formulas (2) and (6) represent the decoupling and coupling processes of the property subspace, respectively. By using formula (6) to adaptively decouple the property subspace in each iteration of the inversion process, and using formula (2) to enhance the correlation of multiple physical parameters of the inversion model in their respective property subspaces, the structural coupling of multiple types of rocks in multiple property subspaces is finally realized.

[0148] The preferred embodiments of the present invention will now be described in detail. However, the scope of protection of the present invention is not limited thereto:

[0149] Example

[0150] Constructing models with different shapes and physical properties, such as Figure 3 As shown in (a), the density and magnetic parameters of models 1-4 are shown in Table 1. The magnetic inclination and declination of the geomagnetic field are set to 60° and -10°, respectively. The forward gravity and magnetic anomaly data with Gaussian white noise are shown in Table 1. Figure 3 (b) and Figure 3 As shown in (c).

[0151] Table 1 Figure 2 The physical properties of the theoretical model

[0152] <![CDATA[Density(g / cm 3 )]]> Magnetization (A / m) Model-1 1.5 3 Model-2 -1.8 2.2 Model-3 -1.8 2.2 Model-4 1.8 1.5

[0153] The underground space was divided into 40×40×20 blocks, and density and magnetization models were obtained using an independent inversion method. Figure 3 (d) and Figure 3 As shown in (e), the theoretical data were further calculated using the method presented in this paper. The number of property subspaces was set to 3, and the property values ​​of the initial model were set to 0. The calculated inversion density and magnetization models are shown below. Figure 3 (f) and Figure 3 As shown in (g). Figure 3 (h) and Figure 3 (i) The property distributions calculated by independent inversion and the inversion method presented in this paper are shown respectively. The property distributions of independent inversion are divergent, and the property distributions of Model 1 and Model 4 are mixed together, making them difficult to distinguish. The density slices of independent inversion show the presence of a weak positive density model volume around the negative density model, and the non-uniform distribution of weak density models in the non-model region. The property distributions calculated by the method presented in this paper have a clear multilinear trend, and can still be distinguished well even for Model 1 and Model 4, which have similar linear property parameter trends. The density model inverted by the method presented in this paper is clear and has a good spatial distribution correspondence with the magnetization model. The results of the Gramian-constrained joint inversion are shown in [reference needed]. Figure 4 Compared to Gramian-constrained joint inversion, the method presented in this paper achieves joint inversion of multiple rock types. Experiments in this paper verify the effectiveness of the proposed method.

[0154] The improved experimental results are primarily due to the fact that this invention, for the first time, simultaneously introduces two novel methods in joint inversion: low-rank subspace partitioning and low-rank structural constraints, used for subspace decoupling and recoupling, respectively. Subspace decoupling technology is proposed for the first time in this invention, and is not addressed in traditional joint inversion methods. The new method can partition multiple physical property parameters located in a common physical property space into their respective low-dimensional physical property subspaces. Furthermore, it utilizes clustered low-rank structural coupling constraints to enhance the correlation of multiple physical property parameters in each independent physical property subspace. Through alternating decoupling and coupling processes during iteration, it ultimately achieves adaptive coupling of multiple physical parameters in the multiple physical property subspaces. The new method can adaptively distinguish the physical property subspaces to which different types of rocks belong, and obtain multi-type rock models with strong correlations among multiple physical property parameters, without relying on the initial model or prior information on multiple physical property parameters. This new method provides a possibility for the joint inversion of complex geological structures.

[0155] Working principle: By utilizing standardized physical property matrices and dynamic K-subspace clustering, the physical property structure patterns corresponding to different lithologies can be automatically identified and decoupled. Then, through clustered singular value decomposition and Schatten-p near-end mapping, low-rank reconstruction of the physical property parameter matrix is ​​achieved in their respective subspaces. Simultaneously, depth weighting matrix is ​​combined to overcome sensitivity decay, box projection applies physical property boundary constraints, and Nesterov momentum accelerates the convergence process, thus forming a closed-loop iterative optimization process. Finally, the inversion results are output under the gradient 2 norm convergence criterion. The above design breaks through the dependence of traditional joint inversion on global linear correlation and realizes adaptive decoupling and coupling of complex multi-physical property parameter spaces, thereby effectively improving the identification accuracy of multi-type lithological structures, the numerical stability of inversion results, and the rationality of geological interpretation.

[0156] It should be noted that, in this document, relational terms such as "first" and "second" are used only to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such process, method, article, or apparatus.

[0157] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention.

Claims

1. A method for joint inversion of multiple geophysical fields based on physical property subspace decoupling, characterized in that, The method comprises the following steps: The gravity and magnetic observation data are arranged, and the inversion region is divided into regularly arranged blocks; Constructing data matrix, sensitivity matrix and data standard deviation matrix based on gravity observation data and magnetic observation data respectively; Constructing deep weighting matrix and setting inversion parameters including number of subspaces, clustering low-rank constraint weighting coefficient, deep constraint weighting coefficient, upper and lower bounds of physical property parameters and Schatten-p norm parameter; Initializing density model, magnetization model, auxiliary iteration variable and Nesterov momentum parameter; Performing K subspace clustering on the normalized physical property matrix to obtain subspace labels of each block; Calculating the gradient of the density and magnetization model and updating the auxiliary iteration variable; Performing singular value decomposition on the physical property parameter matrix in each subspace based on the clustering result and applying clustered Schatten-p norm proximal mapping update; Performing box constraint projection on the updated physical property parameters and projecting the updated density and magnetization parameters into the preset upper and lower bound interval through box projection; Updating the auxiliary iteration variable in the next step by using Nesterov acceleration strategy; If the convergence condition is met, output the inversion result, otherwise return to the K subspace clustering step for iteration.

2. The method of claim 1, wherein, Constructing data matrix, sensitivity matrix and data standard deviation matrix based on gravity observation data and magnetic observation data respectively, comprising: Generating data matrices from gravity and magnetic field observations and Subdividing the subsurface grid and generating gravity and magnetic sensitivity matrices respectively and Generating data standard deviation matrices from gravity and magnetic data respectively and .

3. The method of claim 1, wherein, In the step of constructing the depth weighting matrix, the formula constructing the depth weighting matrix and ; wherein, represents the depth of the th block; represents a constant, typically ; represents a weighted decay exponent.

4. The method of claim 1, wherein, Initializing density and magnetization model, auxiliary iteration variable and momentum parameter, comprising: density initial model magnetization initial model density model auxiliary iteration quantity magnetization model auxiliary iteration quantity nesterov momentum parameter initial value and current iteration number ​ 5. The method of claim 1, wherein, Performing K subspace clustering on the normalized physical property matrix, comprising: Constructing the first standardized property matrix of the second iteration ; In the formula, represents the standardized physical property matrix constructed in the nth inversion iteration; represents the total number of grid cells after the inversion region is divided; represents that the physical property matrix is a 2-row and 1-column real number matrix, each row represents a grid cell, and each column represents a standardized physical property parameter; represents that the physical property matrix is a 2-row and 1-column real number matrix, each row represents a grid cell, and each column represents a standardized physical property parameter;​ wherein, for any column vector in the property matrix , the standardization calculation comprises ; wherein denotes the arithmetic mean of the column vector ; , denotes the sample standard deviation of the column vector ; Represents the first element in the standardized vector. The value of each element; With Generating a cluster index set for data input using K-subspace clustering algorithm ; wherein, the first element of the first element of the the first element of the the first element of the 6. The method of claim 1, wherein, The K subspace clustering step comprises: Standardizing the density and magnetization auxiliary iteration variable to construct the physical property matrix; Setting the number of subspaces, subspace dimension and convergence threshold; Initializing the subspace labels of each block and calculating the initial basis vectors of each subspace; Iteratively updating the subspace labels of each block and recalculating the basis vectors of each subspace; Stopping iteration when the labels no longer change or the target function change rate is less than the threshold, and outputting the subspace index set.

7. The method of claim 1, wherein, The clustered Schatten-p norm proximal mapping update step comprises: Dividing the physical property parameter matrix into several submatrices according to the clustering index; Performing singular value decomposition on each submatrix; Solving the Schatten-p norm proximal mapping problem for each element in the singular value vector; Reconstructing the submatrix using the updated singular value and updating the physical property parameter value of the corresponding index; The Schatten-p norm proximal mapping problem is a scalar optimization problem, which is shown as follows: ； wherein is the original singular value; is a regularization parameter; .

8. The method of claim 1, wherein, In the Nesterov acceleration step, the auxiliary iteration variable in the next step is calculated according to the momentum parameter of the current iteration step and the physical property parameter change.

9. The method of claim 1, wherein, The stop iteration criterion is that the two-norm of the gradient vector is less than a preset tolerance, i.e. .

10. The method of claim 1, wherein, The inversion parameters include the upper and lower bounds of the density constraint and the upper and lower bounds of the magnetization constraint, which are used to limit the reasonable range of the physical property parameters in the box projection step.

Images