Multi-geophysical field joint inversion method based on low-rank structure constraint

Abstract

The application discloses a kind of based on low rank structure constraint multi-geophysical field joint inversion method, it is related to geophysical field joint inversion field, constructs low rank structure coupling joint inversion mathematical optimization model, utilizes singular value and rank to represent physical property parameter space structure, establishes including data fitting term, model constraint term and low rank structure constraint target function, adopts Schatten-p norm to carry out model multi-property low rank constraint, obtains optimization model optimal solution using nonlinear optimization method, establishes the multi-property attribute model with low rank physical property structure, not only can effectively represent and evaluate the global correlation of multi-physical property model in physical property space, but also can reduce the rank of multi-physical property in physical property space by reduction, gradually enhance the model multi-property correlation in iterative calculation process, reduce geophysical inversion multi-solution and interpretation uncertainty, finally establish the multi-property attribute model with structure similarity.

Description

Technical Field

[0001] The present invention belongs to the field of joint inversion of geophysical fields, and specifically to a multi-geophysical field joint inversion method based on low-rank structural constraints. Background Art

[0002] Joint inversion uses multiple physical property parameters to describe the underground spatial structure, obtains a multi-physical attribute model with structural similarity, realizes the complementary advantages of multiple geophysical methods, and reduces the multi-solution and interpretation uncertainty of inversion. It has been an important development direction and frontier hot issue in the field of geophysics in recent years. The optimal joint inversion constraint method includes physical property coupling constraint and structural coupling constraint. The physical property coupling constraint requires the establishment of a multi-physical property relationship of rock samples. Different research areas, different types and different genesis rocks have different physical property structures. When the multi-physical property relationship cannot be accurately established, it will have a greater impact on the joint inversion results of the rock physical property coupling constraint. The existing structural coupling constraint function is relatively simple and cannot evaluate the correlation of multiple physical properties from a global perspective. Summary of the Invention

[0003] The purpose of the present invention is to overcome the defects in the prior art and provide a multi-geophysical field joint inversion method based on low-rank structural constraints. Singular values ​​and ranks are used to characterize the spatial structure of physical property parameters, and a target function is established that includes data fitting terms, model constraint terms and low-rank structural constraint terms. The Schatten-p norm is used to perform low-rank constraints on the model's multi-properties. The nonlinear optimization method is used to obtain the optimal solution of the optimization model, and a multi-property attribute model with a low-rank physical property structure is established.

[0004] To solve the above technical problems, the present invention adopts the following technical solution: a multi-geophysical field joint inversion method based on low-rank structure constraints, comprising the following steps:

[0005] Step 1: Divide the underground space into regular blocks or irregular polyhedrons, and establish the relationship between rock properties and near-surface geophysical field responses. This relationship can be expressed using a linear equation as follows:

[0006]

[0007] where (d1, d2, ..., d I )、(G1,G2,…,G I ) and (m1,m2,…,m I ) represent the response data, sensitivity matrix and parameter vector of I physical property model respectively, through the data d1, d2, ..., d I Get models m1,m2,…,m with structural similarity I ;

[0008] Step 2: Establish a mathematical model for low-rank structural coupling joint inversion:

[0009]

[0010] Where μ is the regularization parameter, and They are data fitting term, model constraint term and low-rank structure coupling term, and their specific function expressions are as follows:

[0011]

[0012] Among them, D I and W I Represent the data standard deviation and depth weighting matrix of the first physical attribute model, ‖·‖ Sp represents the Schatten-p norm of the matrix, where p is a constant ranging from 1 to 2, λ1,λ2,…,λ I represents the regularization parameter.

[0013] Preferably, a nonlinear optimization method is used to minimize the low-rank structure coupling joint inversion mathematical optimization model, which is characterized by comprising the following steps:

[0014] S1: Data preprocessing, setting the initial value parameter k=1;

[0015] S2: Calculate the standard deviation matrix D1, D2, ..., D respectively I and depth weighted matrices W1, W2, ..., W I ;

[0016] S3: Calculate the sensitivity matrix G1, G2, ..., G I ;

[0017] S4: Calculate the Jacobian matrix

[0018] S5: Calculate the objective function gradient matrix

[0019] S6: Perform model update m k+1 =m k +δd, if the model does not converge, re-execute S3-S6 with parameter k=k+1;

[0020] S7: If model m k+1 =m k +δd converges, output m k+1 .

[0021] Preferably, step S1 data preprocessing includes: inputting geophysical data (d1, d2, ..., d I ), set the regularization parameters λ1,λ2,…,λ Iand μ, divide the underground space and initialize the model parameters Create a matrix

[0022] Preferably, the Jacobian matrix calculation process in step S4 includes:

[0023] S41: Calculate the gradient of the data fitting term:

[0024]

[0025] S42: Calculate the gradient of the model constraint:

[0026]

[0027] S43: Perform singular value decomposition and calculate Schatten-p norm;

[0028] S44: Calculation Jacobian matrix

[0029] S45: Calculate the Jacobian matrix of the objective function

[0030] Preferably, the singular value decomposition in step S43 includes: establishing a structural coupling function using the Schatten-p norm, wherein the Schatten-p norm is obtained using a singular value decomposition algorithm, and forming a physical property matrix M = [m1, ..., m I ], the singular value decomposition of the physical property matrix is:

[0031] M=UΣV T ;

[0032] Where U and V represent the singular value vectors, Σ=diag(σ1,…,σ I ) is a diagonal matrix, the diagonal elements σ1,…,σ I are the singular values ​​of the physical property matrix M, the singular values ​​σ1,…,σ I Can be calculated by M T The eigenvalues ​​of M are obtained, and the Schatten-p norm of the matrix M is obtained:

[0033] Preferably, establishing the linear variation trend of multi-physical data is an important step in studying the statistical characteristics of rock physical properties. Minimizing the data projection variance is the primary condition for determining the direction of the main basis vectors. The position and direction of the main basis vectors of the physical property subspace data are obtained, and the optimization model is established:

[0034]

[0035] Where u and v represent unit vectors. The closed-form solution of the optimization model is to calculate the singular value decomposition of the matrix M, select the left and right singular vectors u1 and v1 corresponding to the maximum singular value σ1, establish the direction of the main basis vectors, and use the singular value decomposition to obtain the linear change direction of the multi-physical data of the inversion model.

[0036] Preferably, when high-dimensional data is nested in a low-dimensional subspace, the subspace dimension is the data dimension, that is, the data subspace dimension dim(M) is equal to the rank of the data matrix rank(M), which is equal to the number of non-zero singular values, that is, the norm of the singular value matrix ||Σ||0. Singular value decomposition is used to quantitatively characterize the linear correlation of multi-physical data. The lower the data dimension, the lower the rank, and fewer linear equations can be used to represent the global data, which means that the data has a stronger correlation. For example, when the data is distributed as a straight line, a linear equation can be used to represent the global data, so the data dimension and rank are equal to 1. When the data is distributed on a plane, a set of orthogonal linear equations are used to represent the data on the plane, and the data dimension and rank are equal to 2. Therefore, singular value decomposition can be used to quantitatively characterize the linear correlation of multi-physical data. Singular value decomposition is used to quantitatively characterize the linear correlation of multi-physical data.

[0037] Beneficial effect: The multi-geophysical field joint inversion method based on low-rank structural constraints disclosed in the present invention is a new joint inversion structural coupling idea compared with the existing known technologies. It uses singular values ​​to characterize the spatial distribution characteristics of multiple physical attributes of the model, constructs an objective function containing low-rank structural constraints, realizes the joint inversion of multiple geophysical fields, and improves the correlation of multiple physical attribute models. BRIEF DESCRIPTION OF THE DRAWINGS

[0038] The accompanying drawings are used to provide further understanding of the present invention and constitute a part of the specification. They are used to explain the present invention together with the embodiments of the present invention and do not constitute a limitation of the present invention.

[0039] In the attached figure:

[0040] Figure 1 This is the overall flow chart of the low-rank structure-constrained multi-geophysical field joint inversion method of the present invention;

[0041] Figure 2 It is a schematic diagram representing the spatial distribution characteristics of multiple physical properties of the model of the present invention;

[0042] Figure 3 This is a schematic diagram representing the linear change trend of multiple physical property data of the present invention;

[0043] Figure 4 It is a schematic diagram of the two-dimensional theoretical model and forward gravity and magnetic anomaly of the present invention;

[0044] Figure 5Schematic diagram of independent inversion and joint inversion results of the present invention;

[0045] Figure 6 It is a physical property relationship diagram of the independent inversion and joint inversion models of the present invention. DETAILED DESCRIPTION

[0046] The technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. The following text is only used to describe the implementation method of the multi-geophysical field joint inversion method based on low-rank structure constraints of the present invention, and does not strictly limit the scope of protection specifically requested by the present invention.

[0047] In addition, the technical solutions between the various embodiments can be combined with each other, but they must be based on the fact that ordinary technicians in this field can implement them. When the combination of technical solutions is mutually contradictory or cannot be implemented, it should be deemed that such a combination of technical solutions does not exist and is not within the scope of protection required by the present invention.

[0048] In the first embodiment, Figure 1 、 Figure 2 and Figure 3 As shown, the present invention discloses a multi-geophysical field joint inversion method based on low-rank structure constraints, comprising the following steps:

[0049] Step 1: Divide the underground space into regular blocks or irregular polyhedrons, and establish the relationship between rock properties and near-surface geophysical field responses. This relationship can be expressed using a linear equation as follows:

[0050]

[0051] where (d1, d2, ..., d I )、(G1,G2,…,G I ) and (m1,m2,…,m I ) represent the response data, sensitivity matrix and parameter vector of I physical property model respectively, through the data d1, d2, ..., d I Get models m1,m2,…,m with structural similarity I ;

[0052] Step 2: Establish a mathematical model for low-rank structural coupling joint inversion:

[0053]

[0054] Where μ is the regularization parameter, and They are data fitting term, model constraint term and low-rank structure coupling term, and their specific function expressions are as follows:

[0055]

[0056]

[0057] Among them, D I and W I denote the data standard deviation and depth weighting matrix of the first physical attribute model, respectively, ||·|| Sp represents the Schatten-p norm of the matrix, where p is a constant ranging from 1 to 2, λ1,λ2,…,λ I represents the regularization parameter.

[0058] In the first embodiment, a nonlinear optimization method is used to minimize a low-rank structure coupling joint inversion mathematical optimization model, which is characterized by comprising the following steps:

[0059] S1: Data preprocessing, setting the initial value parameter k=1;

[0060] S2: Calculate the standard deviation matrix D1, D2, ..., D respectively I and depth weighted matrices W1, W2, ..., W I ;

[0061] S3: Calculate the sensitivity matrix G1, G2, ..., G I ;

[0062] S4: Calculate the Jacobian matrix

[0063] S5: Calculate the objective function gradient matrix

[0064] S6: Perform model update m k+1 =m k +δd, if the model does not converge, re-execute S3-S6 with parameter k=k+1;

[0065] S7: If model m k+1 =m k +δd converges, output m k+1 .

[0066] In the first embodiment, step S1 data preprocessing includes: inputting geophysical data (d1, d2, ..., d I ), set the regularization parameters λ1,λ2,…,λ I and μ, divide the underground space and initialize the model parameters Create a matrix

[0067] In the first embodiment, the Jacobian matrix calculation process in step S4 includes:

[0068] S41: Calculate the gradient of the data fitting term:

[0069]

[0070] S42: Calculate the gradient of the model constraint:

[0071]

[0072] S43: Perform singular value decomposition and calculate Schatten-p norm;

[0073] S44: Calculation Jacobian matrix

[0074] S45: Calculate the Jacobian matrix of the objective function

[0075] In the first embodiment, the singular value decomposition in step S43 includes: establishing a structural coupling function using the Schatten-p norm, wherein the Schatten-p norm is obtained using the singular value decomposition algorithm, and forming a physical property matrix M = [m1, ..., m I ], the singular value decomposition of the physical property matrix is:

[0076] M=UΣV T ;

[0077] Where U and V represent the singular value vectors, Σ=diag(σ1,…,σ I ) is a diagonal matrix, the diagonal elements σ1,…,σ I are the singular values ​​of the physical property matrix M, the singular values ​​σ1,…,σ I Can be calculated by M T The eigenvalues ​​of M are obtained, and the Schatten-p norm of the matrix M is obtained:

[0078] In the first embodiment, Figure 3 As shown in the figure, establishing the linear variation trend of multi-physical data is an important step in studying the statistical characteristics of rock physical properties. Taking the minimum data projection variance as the primary condition for determining the direction of the main basis vectors, the position and direction of the main basis vectors of the physical property subspace data are obtained, and the optimization model is established:

[0079]

[0080] Where u and v represent unit vectors. The closed-form solution of the optimization model is to calculate the singular value decomposition of the matrix M, select the left and right singular vectors u1 and v1 corresponding to the maximum singular value σ1, establish the direction of the main basis vectors, and use the singular value decomposition to obtain the linear change direction of the multi-physical data of the inversion model.

[0081] In embodiment one, when high-dimensional data is nested in a low-dimensional subspace, the subspace dimension is the data dimension, that is, the data subspace dimension dim(M) is equal to the rank of the data matrix rank(M), which is equal to the number of non-zero singular values, that is, the norm of the singular value matrix ||Σ||0. Singular value decomposition is used to quantitatively characterize the linear correlation of multi-physical data. The lower the data dimension, the lower the rank, and fewer linear equations can be used to represent the global data, which means that the data has a stronger correlation. For example, when the data is distributed as a straight line, a linear equation can be used to represent the global data, so the data dimension and rank are equal to 1. When the data is distributed on a plane, a set of orthogonal linear equations are used to represent the data on the plane, and the data dimension and rank are equal to 2. Therefore, singular value decomposition can be used to quantitatively characterize the linear correlation of multi-physical data. Singular value decomposition is used to quantitatively characterize the linear correlation of multi-physical data.

[0082] In the second embodiment, Figure 4 、 Figure 5 and Figure 6 As shown in the figure, the effectiveness of the joint inversion method of the present invention is verified by using two-dimensional gravity and magnetic data. Figure 4 As shown in (c), a double-tilted plate model is established, where the density and magnetic properties of model 1 are 1 g / cm 3 and 2A / m, the density and magnetic properties of model 2 are 3g / cm 3 and 2.5A / m, Figure 4 (a) and Figure 4 (b) The forward modeling of gravity anomalies and magnetic anomalies by the theoretical model, respectively. For comparison, the independent gravity and magnetic anomaly inversion method and the multi-geophysical field joint inversion method based on low-rank structure constraints of the present invention are used for comparison.

[0083] In the second embodiment, the underground space is divided into 50×32 rectangular grids, and gravity data and magnetic data are further inverted separately. Figure 5 (a) and Figure 5 (b) are the density and magnetization distributions shown respectively. The density model obtained by inversion alone has a good correspondence with the theoretical model, while the magnetization model diffuses to the deep, which is different from the theoretical model. The method of the present invention is used to jointly invert gravity data and magnetic data, and the inverted density and magnetization distributions are shown respectively. Figure 5 (c) and Figure 5 As shown in (d), the inversion results show that the density and magnetization intensity models obtained by the method of the present invention have similar spatial distribution relationships.

[0084] In Example 2, the comparison of the physical property relationship of the independent inversion model and the physical property relationship of the inversion model of the present invention shows that, Figure 6 As shown in (a), the physical property distribution calculated by the independent inversion method is relatively divergent, such as Figure 6 As shown in (b), the physical property relationship of the inversion model of the method of the present invention has a strong linear correlation, which verifies the effectiveness of the multi-geophysical field joint inversion method based on low-rank structural constraints of the present invention.

[0085] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. For those skilled in the art, the present invention may have various changes and variations. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A multi-geophysical field joint inversion method based on low-rank structural constraints, characterized by: The following steps are involved: Step 1: Divide the underground space into regular blocks or irregular polyhedrons, and establish the relationship between rock properties and near-surface geophysical field responses. This relationship can be expressed using a linear equation as follows: where (d1, d2, ..., d I )、(G1,G2,...,G I ) and (m1,m2,…,m I ) represent the response data, sensitivity matrix and parameter vector of I physical property model respectively, through the data d1, d2, ..., d I Get models m1,m2,…,m with structural similarity I ; Step 2: Establish a mathematical model for low-rank structural coupling joint inversion: Where μ is the regularization parameter, and They are data fitting term, model constraint term and low-rank structure coupling term, and their specific function expressions are as follows: Among them, D I and W I denote the data standard deviation and depth weighting matrix of the first physical attribute model, respectively, ||·|| Sp represents the Schatten-p norm of the matrix, where p is a constant ranging from 1 to 2, λ1,λ2,…,λ I represents the regularization parameter.

2. The multi-geophysical field joint inversion method based on low-rank structural constraints according to claim 1 is characterized in that: A nonlinear optimization method is used to minimize a low-rank structure coupling joint inversion mathematical optimization model, which is characterized by comprising the following steps: S1: Data preprocessing, setting the initial value parameter k=1; S2: Calculate the standard deviation matrix D1, D2, ..., D respectively I and depth weighted matrices W1, W2, ..., W I ; S3: Calculate the sensitivity matrix G1, G2, ..., G I ; S4: Calculate the Jacobian matrix S5: Calculate the objective function gradient matrix S6: Perform model update m k+1 =m k +δd, if the model does not converge, re-execute S3-S6 with parameter k=k+1; S7: If model m k+1 =m k +δd converges, output m k+1 .

3. The multi-geophysical field joint inversion method based on low-rank structural constraints according to claim 2 is characterized in that: Step S1 data preprocessing includes: inputting geophysical data (d1, d2, ..., d I ), set the regularization parameters λ1,λ2,…,λ I and μ, divide the underground space and initialize the model parameters Create a matrix 4. The multi-geophysical field joint inversion method based on low-rank structural constraints according to claim 2 is characterized in that: The Jacobian matrix calculation process in step S4 includes: S41: Calculate the gradient of the data fitting term: S42: Calculate the gradient of the model constraint: S43: Perform singular value decomposition and calculate Schatten-p norm; S44: Calculation Jacobian matrix S45: Calculate the Jacobian matrix of the objective function 5. The multi-geophysical field joint inversion method based on low-rank structural constraints according to claim 4 is characterized in that: The singular value decomposition in step S43 includes: using the Schatten-p norm to establish a structural coupling function, the Schatten-p norm is obtained by the singular value decomposition algorithm, and the physical property matrix M = [m1, ..., m I ], the singular value decomposition of the physical property matrix is: M=UΣV T ; Where U and V represent the singular value vectors, Σ=diag(σ1,…,σ I ) is a diagonal matrix, the diagonal elements σ1,…,σ I are the singular values ​​of the physical property matrix M, the singular values ​​σ1,…,σ I Can be calculated by M T The eigenvalues ​​of M are obtained, and the Schatten-p norm of the matrix M is obtained:

6. The multi-geophysical field joint inversion method based on low-rank structural constraints according to claim 5 is characterized in that: Taking the minimum data projection variance as the primary condition for determining the direction of the main basis vectors, the main basis vector positions and directions of the physical property subspace data are obtained, and the optimization model is established: Where u and v represent unit vectors. The closed-form solution of the optimization model is to calculate the singular value decomposition of the matrix M, select the left and right singular vectors u1 and v1 corresponding to the maximum singular value σ1, establish the direction of the main basis vectors, and use the singular value decomposition to obtain the linear change direction of the multi-physical data of the inversion model.

7. The multi-geophysical field joint inversion method based on low-rank structural constraints according to claim 6 is characterized in that: When high-dimensional data is nested in a low-dimensional subspace, the subspace dimension is the data dimension, that is, the data subspace dimension dim(M) is equal to the rank of the data matrix rank(M) is equal to the number of non-zero singular values, that is, the norm of the singular value matrix ||Σ||0. Singular value decomposition is used to quantitatively characterize the linear correlation of multi-physical data.

Images