A difference-constrained inversion method for three-dimensional airborne transient electromagnetic parameters

Abstract

The invention discloses a three-dimensional aviation transient electromagnetic induced polarization parameter difference constraint inversion method, comprising: constructing a data difference calculation model according to the correlation of aviation transient electromagnetic induced polarization parameters; constructing an objective function of aviation transient electromagnetic three-dimensional inversion based on the data difference calculation model, wherein the variables in the objective function are different induced polarization parameters, the calculation result is a fitting difference, and the constraint term in the objective function is constructed based on the data difference calculation model, and the model weighting term in the objective function is an improved Laplace operator; obtaining aviation transient electromagnetic observation data containing induced polarization effect, and iteratively performing three-dimensional inversion on the aviation transient electromagnetic observation data by Gauss-Newton method based on the objective function until the fitting difference meets a set threshold, thereby obtaining the inversion results of different induced polarization parameters of the underground medium.

Description

Technical Field

[0001] The invention belongs to the technical field of multi-parameter inversion in geophysical exploration, relates to a multi-parameter inversion method of induced polarization parameters in airborne transient electromagnetic exploration, and in particular to a difference-constrained inversion method of three-dimensional airborne transient electromagnetic induced polarization parameters. Background Art

[0002] In the inversion of IP parameters using airborne transient electromagnetic method, the inversion results of some IP parameters are often poor. It is difficult to accurately describe the position and morphological characteristics of the polarization target body of the underground medium from the inversion results, which will cause great difficulties for the later geological interpretation and the delineation of the scope of metal mines. Therefore, in order to improve the effect of multi-parameter inversion, it is necessary to improve the inversion process and reduce the negative impact caused by the multi-solution of the inverse problem as much as possible. At present, there is a lack of a systematic constrained inversion method to solve the problem of poor inversion effect of IP parameters using airborne transient electromagnetic method, that is, a method that can constrain the inversion by using the degree of difference between data to improve the inversion results. Summary of the invention

[0003] In order to solve the above-mentioned technical problems that the inversion effect of induced polarization multi-parameters in the prior art of airborne transient electromagnetic method is not good and there is a lack of effective constraints for induced polarization multi-parameter inversion, the present invention proposes a difference constrained inversion method for three-dimensional airborne transient electromagnetic induced polarization parameter inversion. This data difference constrained inversion method can utilize the similarity relationship between different physical parameters, narrow the scope of the solution space in the process of solving the inverse problem, obtain better induced polarization parameter inversion results, obtain more accurate dielectric induced polarization parameter information, and reduce the difficulty of subsequent geological interpretation.

[0004] To achieve the above object, the present invention provides a difference-constrained inversion method for three-dimensional aviation transient electromagnetic induced polarization parameters, comprising:

[0005] According to the correlation of aviation transient electromagnetic induced polarization parameters, a data difference calculation model is constructed;

[0006] Based on the data difference calculation model, the objective function of airborne transient electromagnetic three-dimensional inversion is constructed, wherein the variables in the objective function are different induced polarization parameters, the calculation results are fitting differences, and the constraint terms in the objective function are constructed based on the data difference calculation model, and the model weighting term in the objective function is an improved Laplace operator;

[0007] The airborne transient electromagnetic observation data containing induced polarization effect are obtained. Based on the objective function, the airborne transient electromagnetic observation data are iteratively inverted in three dimensions by the Gauss-Newton method until the fitting error meets the threshold, and the inversion results of different induced polarization parameters of the underground medium are obtained.

[0008] Optionally, the data difference calculation model includes:

[0009] When the correlation of IP parameters is unknown:

[0010]

[0011] When IP parameters are positively correlated:

[0012]

[0013] When IP parameters are negatively correlated:

[0014]

[0015] Where, s(m1, m2) represents the data difference between the first IP parameter and the second IP parameter, σ1 and σ2 represent the population standard deviations of the first IP parameter and the second IP parameter, respectively. They represent the mean values ​​corresponding to the first and second IP parameters respectively.

[0016] Optionally, the objective function is:

[0017]

[0018]

[0019] Where Φ(m) is the fitting error, They represent the observed data difference term, the model data difference term, and the electromagnetic parameter data difference term, respectively; λ1 and λ2 are different regularization factors; W d is the data weighting matrix; W m is the model weight matrix; s is the data difference matrix; d obs is the observed data; F(m) is the forward numerical simulation value; m is the model parameter vector; m0 is the reference model parameter vector; the superscript T indicates the transpose of the matrix;

[0020] Optionally, in the modified Laplacian, the diagonal elements are:

[0021]

[0022] The differential coefficients of the off-diagonal elements in three directions are:

[0023]

[0024] Where c represents the diagonal element, m(i, j, k) represents the electromagnetic parameter in the grid at the coordinate (i, j, k), i, j, k) represents the coordinate sequence of the grid in the x, y, and z directions respectively, and c x ,c yand c z are the differential coefficients in the x-, y- and z-directions respectively. Δx(i), Δy(j) and Δz(k) represent the sizes of the grids in the x-, y- and z-directions respectively. α is a constant. w is a weighting coefficient ranging from 0 to 1.

[0025] Optionally, the aviation transient electromagnetic induced polarization parameters include resistivity, polarizability, time parameter and frequency correlation coefficient.

[0026] Optionally, the process of performing a three-dimensional inversion of airborne transient electromagnetic observation data by the Gauss-Newton method includes:

[0027] Obtaining a data weighting matrix based on the observed data;

[0028] Based on the inversion area, the initial model of the underground space is constructed by the grid division method, and the forward numerical simulation value of the initial model is obtained, wherein the forward numerical simulation value is the transient electromagnetic response;

[0029] The initial model is updated to generate an updated model, wherein the process of updating the initial model includes: calculating the model weighting items according to the grid in the initial model; calculating the time domain partial derivative matrix of the observed data with respect to the model variables; calculating the data difference degree of the model variables according to the data difference degree calculation model, and calculating the partial derivative matrix of the data difference degree with respect to the model variables; obtaining the model update variables according to the time domain partial derivative matrix and the partial derivative matrix; updating the model variables in the initial model according to the model update variables to generate an updated model, wherein the model variables are different induced polarization parameters;

[0030] The updated model is forward modeled to generate forward numerical simulation values, and the fitting difference between the observed data and the forward modeling results is calculated through the objective function. If the fitting difference is greater than a threshold, the updated model is iteratively updated, forward modeled, and inverted for fitting difference calculation based on the process of updating the initial model until the fitting difference is less than or equal to the threshold, thereby obtaining inversion results of different induced polarization parameters, and the inversion results of the different induced polarization parameters are plotted to obtain the final inversion model of different induced polarization parameters.

[0031] Optionally, the process of acquiring the model update variable includes:

[0032]

[0033] Among them, Δm is the model update variable, J is the time domain partial derivative matrix of the observed data with respect to the model variable, and B s is the matrix of partial derivatives of the data variance with respect to the model variables.

[0034] Optionally, the process of obtaining the data weighting matrix includes:

[0035]

[0036] Compared with the prior art, the present invention has the following advantages and technical effects:

[0037] The data difference constrained inversion method for airborne transient electromagnetic induced polarization parameter inversion provided by the present invention can solve the problem of low resolution of underground medium polarization target body in the original traditional method to a certain extent by constraining the data similarity according to the correlation in different physical property parameters. The present invention corrects the inversion result through the data difference constraint to achieve the purpose of improving the polarization body inversion effect. BRIEF DESCRIPTION OF THE DRAWINGS

[0038] The drawings constituting a part of the present application are used to provide a further understanding of the present application. The illustrative embodiments and descriptions of the present application are used to explain the present application and do not constitute an improper limitation on the present application. In the drawings:

[0039] Figure 1 The implementation process of the data difference degree constrained inversion method of aviation transient electromagnetic induced polarization parameters according to the embodiment of the present invention;

[0040] Figure 2 An example calculation result diagram of data difference generated under different data structures in an embodiment of the present invention;

[0041] Figure 3 is a schematic diagram of a theoretical model of an embodiment of the present invention;

[0042] Figure 4 A comparison diagram of resistivity inversion results of an embodiment of the present invention;

[0043] Figure 5 A comparison diagram of polarizability inversion results of an embodiment of the present invention;

[0044] Figure 6 A comparison diagram of time constant inversion results of an embodiment of the present invention;

[0045] Figure 7 It is a comparison diagram of the frequency correlation coefficient inversion results of an embodiment of the present invention. DETAILED DESCRIPTION

[0046] It should be noted that, in the absence of conflict, the embodiments and features in the embodiments of the present application can be combined with each other. The present application will be described in detail below with reference to the accompanying drawings and in combination with the embodiments.

[0047] It should be noted that the steps shown in the flowcharts of the accompanying drawings can be executed in a computer system such as a set of computer executable instructions, and that, although a logical order is shown in the flowcharts, in some cases, the steps shown or described can be executed in an order different from that shown here.

[0048] The present invention provides a data difference constraint inversion method for three-dimensional aviation transient electromagnetic induced polarization parameters, which mainly includes the following contents:

[0049] 1. Establish the expression corresponding to the data difference calculation model between different physical parameters. The physical parameters m1 and m2 both contain n data. The expressions corresponding to the data difference calculation model are as follows:

[0050]

[0051] Where, s(m1, m2) represents the data difference between the first IP parameter and the second IP parameter, σ1 and σ2 represent the population standard deviations of the first IP parameter and the second IP parameter, respectively. They represent the mean values ​​of the first IP parameter and the second IP parameter respectively. represents the mean value of the physical property parameter, j is the data label of the induced polarization parameter, and n is the total number of induced polarization parameter data.

[0052] The above-mentioned physical property parameters refer to the three-dimensional aviation transient electromagnetic parameters in the present invention. Formula (1) is an expression for measuring the data difference when the correlation between physical property parameters cannot be clearly described, that is, when the physical property parameters are unknown. Formula (2) is an expression for the data difference when the two sets of physical property parameters are positively correlated. Formula (3) is an expression for the data difference when the two sets of physical property parameters are negatively correlated. Different expressions can be used to represent the degree of difference between data according to the specific data characteristics. The more similar the data structure is, the greater the data difference is. Figure 2 Formula (1) is used to show the calculation effect of data difference. The calculated data difference of two groups of data with similar data structures is small (e.g. Figure 2 (a), (b), (c), (d)) show the example diagrams of different data structure combinations and corresponding data difference curves. The data difference values ​​calculated for two sets of data with different data structures are large (e.g. Figure 2 (e) and (f) show examples of different data structure combinations and corresponding data discrepancy curves. Data discrepancy is the data discrepancy.

[0053] 2. Establish the objective function of airborne transient electromagnetic 3D inversion including data difference constraint:

[0054]

[0055] Where Φ(m) is the fitting error, They represent the observed data difference term, the model data difference term, and the electromagnetic parameter data difference term, respectively. λ is the regularization factor, and λ1 and λ2 are Corresponding different regularization factors; W d is the data weighting matrix; W m is the improved model weighting matrix; s is the data difference matrix; d obs is the observation data (airborne transient electromagnetic observation data); F(m) is the forward numerical simulation value (the transient electromagnetic response generated by the model); m is the model parameter (IP parameter) vector; m0 is the reference model parameter (variable of the initial space model, i.e., the IP parameter set initially) vector; the superscript T indicates the transpose of the matrix;

[0056] Adding the data difference constraint established in step 1 to the inversion can reduce the scope of the solution space and utilize the inherent correlation between the induced polarization parameters to obtain several groups of structurally similar induced polarization parameter models, so that the inversion results can more clearly and accurately describe the polariton information.

[0057] The improved model weight matrix is ​​an improved Laplace operator, and its diagonal elements are

[0058]

[0059] The off-diagonal elements are arranged in the form of central differences, and the difference coefficients in the three directions are c x ,c y and c z :

[0060]

[0061] Where c represents the diagonal element, m(i, j, k) represents the electromagnetic parameter in the grid at the coordinate (i, j, k), i, j, k represent the coordinate sequence of the grid in the x, y, z directions respectively, and c x ,c y and c z are the differential coefficients in the x-, y- and z-directions respectively. Δx(i), Δy(j) and Δz(k) represent the sizes of the i-th, j-th and k-th grids in the x-, y- and z-directions respectively. α is a relatively small constant. w is the weighting coefficient, which can be selected between 0 and 1.

[0062] In the inversion of induced polarization parameters, the absolute scale differences between various physical parameters are too large. The present invention takes into account the absolute differences of model variables in the improved Laplace operator. Adding it to the inversion can improve the morbidity of the inversion equation to a certain extent. Adding differential coefficients can also obtain a suitable inversion result as expected. If a smooth inversion result is expected, the weighting coefficient w can be close to 1. If the boundary information of the model is more emphasized, the weighting coefficient w can be set to a constant close to 0.

[0063] 3. Using the objective function in step 2, the Gauss-Newton method is used to perform three-dimensional inversion on the airborne transient electromagnetic observation data containing induced polarization effect. The airborne transient electromagnetic observation data (induced electromotive force attenuation curve) is input into the inversion program. When the iterative fitting error meets the set threshold, the inversion results of various induced polarization parameters of the underground medium are output.

[0064] like Figure 1 As shown, the present invention describes the implementation process of the above technical solution in detail in conjunction with the accompanying drawings.

[0065] Step 1: Establish data difference expressions between different physical parameters. For the Cole-Cole model used in the present invention to describe the IP characteristics of rocks and minerals, there are four IP parameters, namely, resistivity ρ0, polarizability m, time constant τ, and frequency correlation coefficient c. The relationship between the above IP parameters is not clear, so three sets of expressions are designed, namely:

[0066]

[0067] Among them, the subscript σ of different parameters represents the mean square error of the corresponding parameters, and the superscript “-” of the IP parameter represents the mean value of the corresponding IP parameter.

[0068] Step 2: Based on the data difference expression set in step 1, construct the objective function described by formula (6-9), and perform Gauss-Newton method inversion from step 3 to step 10.

[0069] Step 3: Input the airborne transient electromagnetic observation data (induced electromotive force attenuation curve) and the measurement point location into the inversion program, and determine the data weighting matrix based on the observation data, where the data weighting diagonal matrix is ​​expressed as:

[0070]

[0071] Step 4: Divide the inversion target area into hexahedral grids and set the initial model of the underground space. If there is a priori model, it can be set as the priori model. If there is no priori model, it can be set as a uniform half-space model. Perform forward calculation on the initial model to obtain the transient electromagnetic response of the initial model.

[0072] Step 5: Based on the hexahedral mesh divided in step 4, calculate the model weighting matrix described by formula (10-13).

[0073] Step 6: Calculate the time domain partial derivative matrix of the observed data with respect to the model variables (IP parameters).

[0074] Step 7: Calculate the data difference and the partial derivative matrix of the data difference with respect to the model variables according to the multi-property parameter model.

[0075] Step 8: According to the following iterative formula, the objective function equation group of the three-dimensional inversion is formed, and the objective function equation group of the three-dimensional inversion is solved to obtain the model update variable Δm:

[0076]

[0077] Among them, Δm is the model update variable, J is the time domain partial derivative matrix of the observed data with respect to the model variable, and B s is the matrix of partial derivatives of the data variance with respect to the model variables.

[0078] Step 9: Add the model update variable Δm to the original model and update the model:

[0079] m k+1 =m k +Δm (19)

[0080] Among them, m k is the corresponding model variable in the kth iteration process, and the model is updated by adjusting the model variables.

[0081] Step 10: Input the updated model into the forward modeling program and compare the data fitting difference between the forward modeling result data and the observed data. If the fitting difference is small enough and meets the set threshold, go to step 11; if it does not meet the inversion iteration termination condition, go to step 5.

[0082] Step 11: When the above threshold is met, the model variables corresponding to the current model, namely the IP parameters, are output as the inversion results of different IP parameters, and the inversion results of the IP parameters are plotted as the inversion model of the IP parameters.

[0083] Here, a set of theoretical models are calculated to compare the results of data difference constrained inversion and traditional inversion.

[0084] Two low-resistance prisms are buried in the uniform half space, one of which is a polarizability body and the other is non-polarizability. The surrounding rock resistivity is 200Ω·m, the polarizability is 0.01, the time constant is 0.0001s, and the frequency correlation coefficient is 0.01; the length and width of the target body are 80m and 60m respectively, the burial depth range is 66.2m-154.3m, the distance between the two target bodies is 120m, the target body resistivity is 50Ω·m, one of the polarizability bodies has a polarizability of 0.5, the time constant is 0.005s, the frequency correlation coefficient is 0.5, and the induced polarization parameters of the other polarizability body are the same as those of the surrounding rock. The radius of the transmitting coil is set to 15m, the transmitting current is 100A, the transmitting coil and the receiving coil are both located in the air at 30m, with a total of 121 measuring points, and the model is as follows: Figure 3 As shown, Figure 3 The figure shows a schematic diagram of the theoretical model, where the parts are labeled: (a) XOY plane view, (b) XOZ cross-section view.

[0085] The inversion calculation is performed using the implementation process of the present invention, and the results are as follows: Figure 4-Figure 7 As shown:

[0086] Figure 4 This is a comparison chart of resistivity inversion results, where the parts are labeled: (a) resistivity model, (b) conventional inversion result, (c) data difference constrained inversion result, and ρ represents resistivity. Figure 5 This is a comparison chart of polarizability inversion results, where the parts are labeled: (a) resistivity model, (b) conventional inversion result, (c) data difference constrained inversion result, m represents the corresponding polarizability, and the unit is dimensionless. Figure 6 This is a comparison chart of the time constant inversion results, where the parts are labeled: (a) resistivity model, (b) conventional inversion results, (c) data difference constrained inversion results, where T represents the time constant in seconds. Figure 7 The figure is a comparison chart of the frequency correlation coefficient inversion results, where the parts are labeled: (a) resistivity model, (b) conventional inversion result, (c) data difference constrained inversion result, where c is the frequency correlation coefficient, and the unit is dimensionless.

[0087] For zero-frequency resistivity, both conventional inversion and data difference constrained inversion can achieve good results. For the other three IP parameters, the conventional inversion has a small correction amount for the IP parameters. After adding the data difference constraint, the correction amount for the other three parameters in the inversion iteration process becomes larger. Although there is still a gap between the inversion results and the actual IP parameter values, judging by their relative sizes, the other three parameters can still individually characterize the position of the polarimetric body.

[0088] The calculation results of this model show that the data difference constrained inversion method can solve the problem of low resolution of underground medium polarization target volume in the original traditional method to a certain extent, and achieve the purpose of improving the polarization volume inversion effect.

[0089] The above are only preferred specific implementations of the present application, but the protection scope of the present application is not limited thereto. Any changes or substitutions that can be easily thought of by a person skilled in the art within the technical scope disclosed in the present application should be included in the protection scope of the present application. Therefore, the protection scope of the present application should be based on the protection scope of the claims.

Claims

1. A difference-constrained inversion method for three-dimensional airborne transient electromagnetic induced polarization parameters, characterized in that: include: According to the correlation of aviation transient electromagnetic induced polarization parameters, a data difference calculation model is constructed; Based on the data difference calculation model, the objective function of airborne transient electromagnetic three-dimensional inversion is constructed, wherein the variables in the objective function are different induced polarization parameters, the calculation results are fitting differences, and the constraint terms in the objective function are constructed based on the data difference calculation model, and the model weighting term in the objective function is an improved Laplace operator; Acquire airborne transient electromagnetic observation data containing induced polarization effect, and perform iterative three-dimensional inversion on the airborne transient electromagnetic observation data by Gauss-Newton method based on the objective function until the fitting error meets the threshold, and obtain the inversion results of different induced polarization parameters of the underground medium; The data difference calculation model includes: When the correlation of IP parameters is unknown: When IP parameters are positively correlated: When IP parameters are negatively correlated: Where, s(m1, m2) represents the data difference between the first IP parameter and the second IP parameter, σ1 and σ2 represent the population standard deviations of the first IP parameter and the second IP parameter, respectively. They represent the mean values ​​corresponding to the first and second IP parameters respectively.

2. The method according to claim 1, characterized in that The objective function is: Where Φ(m) is the fitting error, They represent the observed data difference term, the model data difference term, and the electromagnetic parameter data difference term, respectively; λ1 and λ2 are different regularization factors; W d is the data weighting matrix; W m is the model weight matrix; s is the data difference matrix; d obs is the observed data; F(m) is the forward numerical simulation value; m is the model parameter vector; m0 is the reference model parameter vector; the superscript T represents the transpose of the matrix.

3. The method according to claim 2, characterized in that In the modified Laplacian operator, the diagonal elements are: The differential coefficients of the off-diagonal elements in three directions are: Where c represents the diagonal element, m(i, j, k) represents the electromagnetic parameter in the grid at the coordinate (i, j, k), i, j, k represent the coordinate sequence of the grid in the x, y, z directions respectively, and c x ,c y and c z are the differential coefficients in the x-, y- and z-directions respectively. Δx(i), Δy(j) and Δz(k) represent the sizes of the grids in the x-, y- and z-directions respectively. α is a constant. w is a weighting coefficient ranging from 0 to 1.

4. The method according to claim 1, characterized in that: Aviation transient electromagnetic induced polarization parameters include resistivity, polarizability, time parameters and frequency correlation coefficient.

5. The method according to claim 3, characterized in that: The process of three-dimensional inversion of airborne transient electromagnetic observation data by Gauss-Newton method includes: Obtaining a data weighting matrix based on the observed data; Based on the inversion area, the initial model of the underground space is constructed by the grid division method, and the forward numerical simulation value of the initial model is obtained, wherein the forward numerical simulation value is the transient electromagnetic response; The initial model is updated to generate an updated model, wherein the process of updating the initial model includes: calculating the model weighting items according to the grid in the initial model; calculating the time domain partial derivative matrix of the observed data with respect to the model variables; calculating the data difference degree of the model variables according to the data difference degree calculation model, and calculating the partial derivative matrix of the data difference degree with respect to the model variables; obtaining the model update variables according to the time domain partial derivative matrix and the partial derivative matrix; updating the model variables in the initial model according to the model update variables to generate an updated model, wherein the model variables are different induced polarization parameters; The updated model is forward modeled to generate forward numerical simulation values, and the fitting difference between the observed data and the forward modeling results is calculated through the objective function. If the fitting difference is greater than a threshold, the updated model is iteratively updated, forward modeled, and inverted for fitting difference calculation based on the process of updating the initial model until the fitting difference is less than or equal to the threshold, thereby obtaining inversion results of different induced polarization parameters, and the inversion results of the different induced polarization parameters are plotted to obtain the final inversion model of different induced polarization parameters.

6. The method according to claim 5, characterized in that The process of obtaining the model update variables includes: Among them, Δm is the model update variable, J is the time domain partial derivative matrix of the observed data with respect to the model variable, and B s is the matrix of partial derivatives of the data variance with respect to the model variables.

7. The method according to claim 5, characterized in that The process of obtaining the data weighting matrix includes:

Images