Electrical source airborne transient electromagnetic three-dimensional focused inversion method

Abstract

The application discloses a kind of electric source ground air transient electromagnetic three-dimensional focusing inversion methods, comprising the following steps: (1) selecting work area data file;(2) reading parameter file;(3) according to parameter setting initial model;(4) setting inversion termination condition;(5) calculating frequency domain three-dimensional forward response again with G-S transformation is converted into time domain response;(6) by Gauss Newton inversion method obtains L1 regular inversion objective function;(7) obtains model correction Δm;(8) update model vector;Calculate frequency domain three-dimensional forward response again with G-S transformation is converted into time domain response;(9) judge whether to meet stop condition, if not satisfied, return to step (5), otherwise execute step (10);(10) terminate iteration, output current model m.The application can restore model real electrical characteristics better, with the effect of boundary sharpening, make up for the defect of current ground air transient electromagnetic detection data inversion interpretation system.

Description

Technical Field

[0001] This invention belongs to the field of geophysical time-domain electromagnetic exploration, and relates to three-dimensional inversion methods and technologies for ground and semi-airborne transient electromagnetic data, particularly a three-dimensional focusing inversion method for electrical source ground-to-air transient electromagnetic data. Background Technology

[0002] Ground-to-air transient electromagnetic methods are an important branch of transient electromagnetic exploration methods based on electrical sources. For a long time, due to limitations in transient electromagnetic forward modeling theory, data volume, and computation time, practical interpretation techniques for ground-based and semi-airborne transient electromagnetic sources have remained primarily one-dimensional. Three-dimensional characteristics are more prevalent in natural media, making three-dimensional inversion interpretation techniques more reasonable. However, when interpreting data in three dimensions, computational efficiency and model interpretation accuracy are insufficient. Among these challenges, model structure reconstruction and electrical interface identification are key research areas. To address these issues, this invention, based on three-dimensional forward modeling theory, improves computational efficiency and adapts to the practical needs of large datasets, building upon existing conventional three-dimensional inversion research. It also improves the inversion model's expectation to enhance electrical interface resolution. These two aspects are crucial to the current development of transient electromagnetic three-dimensional inversion technology and are of great significance to the widespread practical application of inversion algorithms. Summary of the Invention

[0003] The purpose of this invention is to overcome the shortcomings of the prior art and provide a three-dimensional focusing inversion method for electrical source ground-to-space transient electromagnetic data, which can better restore the true electrical characteristics of the model, has the effect of boundary sharpening, and makes up for the defects of the current ground-to-space transient electromagnetic data inversion and interpretation system.

[0004] The objective of this invention is achieved through the following technical solution: a three-dimensional focusing inversion method for transient electromagnetic sources at ground and space, comprising the following steps:

[0005] (1) Select the work area data file, including the parameter file and the observation data file d;

[0006] (2) Read the parameter file, including: receiving coil area, response time, wire source length, offset, peak current, background resistivity, and response time;

[0007] (3) Set the initial model according to the parameters, including: 3D mesh generation and physical property filling;

[0008] (4) Set the inversion termination condition, including: maximum number of iterations I max Minimum fitting error Rms; Initialize iteration count i = 1; Initialize model m: The elements in model m are resistivity. If the initial model m does not include prior information, then arbitrarily set a resistivity value. If there is prior information, then set the resistivity value according to the prior information.

[0009] (5) Calculate the three-dimensional forward response in the frequency domain, then use the adjoint forward method to solve the frequency domain rate Jacobian matrix and update the data item Hessian matrix; and obtain the transient response of the magnetic field in the time domain by converting the frequency domain response to the time domain response, and calculate the root mean square error Rms1.

[0010] (6) The frequency domain Jacobian matrix is ​​converted into the time domain Jacobian matrix by GS transformation, and the L1 canonical inversion objective function is obtained by Gauss-Newton inversion method.

[0011] (7) The L1 regularization term in the L1 regularization inversion objective function is calculated by using the iterative reweighted least squares method, and then the model correction amount Δm is obtained by solving the linear equation system by using the conjugate gradient method.

[0012] (8) Update the model vector: m = m + Δm; recalculate the three-dimensional forward response in the frequency domain, obtain the transient response of the magnetic field in the time domain by converting the frequency domain response to the time domain response, and calculate the root mean square error Rms2;

[0013] (9) Calculate the root mean square error ΔRms: ΔRms=|Rms1-Rms2| / Rms1, and determine whether the root mean square error is less than the minimum fitting error Rms, or whether the number of iterations i is less than the maximum number of iterations I. max If none of them are satisfied, let i = i + 1, take the updated model in (8) as the current model, and return to step (5); otherwise, execute step (10);

[0014] (10) Terminate the iteration and output the current model m.

[0015] The specific implementation method of the three-dimensional forward response is as follows: the frequency domain electric field double curl equation is rewritten into a weighted residual equation using the Galerkin method, and the electric field is obtained by solving the equation system using the bistable conjugate gradient method with incomplete LU decomposition; finally, the auxiliary magnetic field is calculated using the difference method.

[0016] The cell weighted residual equation is:

[0017]

[0018] in, This represents the electric field in the frequency domain to be solved. This indicates the curl of the electric field in the frequency domain; Let f represent the background field (obtained directly from the analytical expression), ∂f / ∂ω represent the vector interpolation function, ∂f / ∂ω represent the curl, ω be the angular frequency, μ₀ be the free permeability, σ be the conductivity, dΩ be the volume integral, and Δσ be the residual conductivity. Indicates the background electric field;

[0019] The first term on the left-hand side of equation (1) can be expressed as:

[0020]

[0021] K le It is a coefficient matrix;

[0022] The second term on the left is represented as:

[0023]

[0024] This is a first-order interpolation function, where the superscript 'e' represents a cell in the 3D mesh, and the subscripts 'p' and 'q' represent the edge numbers. K... 2e It is a coefficient matrix;

[0025] The right-hand side is represented as:

[0026]

[0027] Sorting and assembling all the element meshes results in a large sparse linear system of equations. Solving this system yields the edge electric field components of the entire mesh, expressed as follows:

[0028] KE a =B (5)

[0029] Where K is the global stiffness matrix, composed of all element left-hand terms, obtained through formulas (2) and (3), E a The electric field coefficient to be solved is B, which is composed of the right-hand terms of all cells and is obtained by formula (4).

[0030] Finally, through E a The magnetic field H is obtained using a cubic spline interpolation algorithm. a .

[0031] The L1 regularized inversion objective function is expressed as follows:

[0032]

[0033] Equation (6) can be expanded as follows:

[0034]

[0035] The first term on the right-hand side of equation (7) is the data item, and the second term is the L1 regularization term, where N m x represents the number of underground grid cells. j W is the j-th element of the regularization vector; d The data is weighted matrix, with the superscript T indicating matrix transpose; λ is the regularization factor, d is the observation data file, and m... ref This is the reference model vector.

[0036] The iterative reweighted least squares method is as follows: In equation (7), the L1 regularization term is approximately represented by matrix expansion as follows:

[0037]

[0038] Weighted matrix V m element V m It changes with each iteration and is related to the result of the previous iteration; that is, it is a reweighted relationship after each iteration; it can be expressed in the following form:

[0039]

[0040] Where ε is a small value;

[0041] Taking the derivative of equation (7) with respect to the model increment, we get:

[0042]

[0043] The square brackets on the left side of formula (10) represent the Hessian matrix of the data item, J. jcb W is a Jacobian matrix; m Let m be the constraint matrix of model m;

[0044] The model correction amount Δm is obtained through formula (10).

[0045] The beneficial effects of this invention are as follows: The transient electromagnetic three-dimensional focusing inversion imaging method for electric sources provided by this invention is a ground-to-air time-domain electromagnetic detection method for electric source emission and secondary field signals received by ground or airborne UAVs. This invention is the first to introduce the iterative reweighted least squares method to solve the L1 canonical zero-point non-differentiability problem in the field of semi-airborne transient electromagnetic three-dimensional inversion. Simultaneously, the inversion can better restore the true electrical characteristics of the model and has a boundary sharpening effect, thus overcoming the shortcomings of current ground-to-air transient electromagnetic detection data inversion interpretation systems. Attached Figure Description

[0046] Figure 1 A flowchart for the three-dimensional focusing inversion of transient electromagnetic sources in the air-to-ground region;

[0047] Figure 2 This is a schematic diagram of the grid division of the observation area at the measuring points;

[0048] Figure 3 A stereoscopic image of L1 and L2 canonical inversion of semi-airborne transient electromagnetic L2;

[0049] Figure 4 Comparison of results from different work on a certain survey line;

[0050] Figure 5 Comparison of L1 and L2 canonical three-dimensional inversion results of transient electromagnetic data from the grounding wire source of the ink mine. Detailed Implementation

[0051] This invention presents a three-dimensional focusing inversion method for transient electromagnetic events between ground and air, including improvements to the regularization term of the objective function and the selection of regularization factors and reference models. In the selection of three-dimensional inversion techniques, this invention employs the second-order convergent Gauss-Newton method. Its overall objective function, the Hessian approximation matrix, can be considered as the sum of the approximation Hessian term and the exact regularization term Hessian term in the data fitting project, which is beneficial for achieving robust inversion or focusing inversion strategies with automatic variable weights. To enhance the inversion model's ability to distinguish electrical abrupt changes, this invention adopts the L1 norm for the regularization term of the objective function and introduces an iterative reweighted least squares strategy within the Gauss-Newton iterative framework for approximate solution. Its iterative stability is superior to optimization strategies that directly minimize the L1 norm approximation or use subgradient search. The Gauss-Newton-type Hessian approximation of the data fitting utilizes the Jacobian matrix inner product form, employing the adjoint forward method to first solve the frequency domain Jacobian matrix in parallel to transform it to the time domain. To minimize the initial model dependency of the inversion algorithm, a smooth model obtained from inversion with a large regularization factor is used as the initial model strategy for inversion with a small regularization factor. Based on the above research, inversion of the data can better approximate the true electrical structure of the model. The invention will now be further described in detail with reference to the accompanying drawings and specific embodiments.

[0052] like Figure 1 As shown, a three-dimensional focusing inversion method for ground-space transient electromagnetic sources includes the following steps:

[0053] (1) Select the work area data file, including the parameter file and the observation data file d;

[0054] (2) Read the parameter file, including: receiving coil area, response time, wire source length, offset, peak current, background resistivity, and response time;

[0055] (3) Set the initial model according to the parameters, including: 3D mesh generation and physical property filling;

[0056] (4) Set the inversion termination condition, including: maximum number of iterations I max Minimum fitting error Rms; Initial iteration count i = 1; Initialize model m: The elements in model m are electrical conductivity. If no prior information is added to the initial model m, an arbitrary electrical conductivity value is set. If there is prior information, the electrical conductivity value is set according to the prior information. Prior information refers to geological data, which is used to set the initial values ​​of the reference model.

[0057] (5) Calculate the three-dimensional forward response in the frequency domain, then use the adjoint forward method to solve the frequency domain rate Jacobian matrix and update the data item Hessian matrix; and obtain the transient response of the magnetic field in the time domain by converting the frequency domain response to the time domain response (GS transformation), and calculate the root mean square error Rms1.

[0058] The specific implementation method of the three-dimensional forward response is as follows: the frequency domain electric field double curl equation is rewritten into a weighted residual equation using the Galerkin method. The equation system is solved by the bistable conjugate gradient method with incomplete LU decomposition to obtain the main field (frequency domain electric field). Finally, the auxiliary field (frequency domain magnetic field) is calculated using the difference method. The total frequency domain magnetic field is obtained by GS transformation to obtain the time domain magnetic field transient response dBz / dt.

[0059] The cell weighted residual equation is:

[0060]

[0061] in, This represents the electric field in the frequency domain to be solved. This indicates the curl of the electric field in the frequency domain; Let f represent the background field (obtained directly from the analytical expression), and let f represent the vector interpolation function. denoted by curl; ω is the angular frequency; μ0 is the free permeability; σ is the conductivity; dΩ is the volume integral; Δσ is the residual conductivity. Indicates the background electric field;

[0062] The first term on the left-hand side of equation (1) can be expressed as:

[0063]

[0064] K le It is a coefficient matrix;

[0065] The second term on the left is represented as:

[0066]

[0067] This is a first-order interpolation function, where the superscript 'e' represents a cell in the 3D mesh, and the subscripts 'p' and 'q' represent the edge numbers. K... 2e It is a coefficient matrix;

[0068] The right-hand side is represented as:

[0069]

[0070] Sorting and assembling all the element meshes results in a large sparse linear system of equations. Solving this system yields the edge electric field components of the entire mesh, expressed as follows:

[0071] KE a =B (5)

[0072] Where K is the global stiffness matrix, composed of the left-hand side terms of all elements (Equations 2 and 3), E aLet B be the electric field coefficient to be solved, composed of the right-hand terms of all cells (Formula 4). Finally, by applying E... a The frequency domain magnetic field H is obtained using a cubic spline interpolation algorithm. a The obtained frequency domain magnetic field H a and background magnetic field H p The results are summed and then transformed by GS to obtain the transient response dBz / dt of the magnetic field in the time domain. The response dBz / dt is denoted as F(m).

[0073] (6) The frequency domain Jacobian matrix is ​​converted to the time domain Jacobian matrix by the GS transform, and the L1 canonical inversion objective function is obtained by the Gauss-Newton inversion method. The L1 canonical inversion objective function is expressed as:

[0074]

[0075] Equation (6) can be expanded as follows:

[0076]

[0077] The first term on the right-hand side of equation (7) is the data item, and the second term is the L1 regularization term, where N m x represents the number of underground grid cells. j W is the j-th element of the regularization vector; d The data is weighted matrix, with the superscript T indicating matrix transpose; λ is the regularization factor, d is the observation data file, and m... ref This is the reference model vector.

[0078] (7) The L1 regularization term in the L1 regularization inversion objective function is calculated by using the iterative reweighted least squares method, and then the model correction amount Δm is obtained by solving the linear equation system by using the conjugate gradient method.

[0079] The iterative reweighted least squares method is as follows: In equation (7), the L1 regularization term is approximately represented by matrix expansion as follows:

[0080]

[0081] Weighted matrix V m element V m It changes with each iteration and is related to the result of the previous iteration; that is, it is a reweighted relationship after each iteration; it can be expressed in the following form:

[0082]

[0083] Where ε is a small value, here taken as 1e-3.

[0084] Weighted matrix V mThe adjustment allows the regularization term to retain the properties of the L1 norm while resolving the issue of non-differentiability in its solution domain. The L1 regularization term in the objective function makes it easier for the model's gradient vector to obtain sparse solutions. The inversion process preserves the most basic features of the model, and the boundaries gradually sharpen, thus achieving a focusing effect.

[0085] Taking the derivative of equation (7) with respect to the model increment, we get:

[0086]

[0087] The square brackets on the left side of formula (10) represent the Hessian matrix of the data item, J. jcb The Jacobian matrix is ​​obtained from step (5); W d W is a weighted matrix for the data. m Let V be the constraint matrix of model m, where the superscript T denotes matrix transpose; λ is the regularization factor, and V... m The iterative reweighting matrix is ​​given by Δm, where Δm is the model correction, F(m) is the model response (in 3D forward modeling), d is the observation data file, and m is the model vector (with elements representing conductivity σ). ref This is the reference model vector.

[0088] The model correction amount Δm is obtained through formula (10). The inversion process is to add the solved Δm and m to obtain a new model vector after correction, and then use the inversion step (9) to determine whether the result meets the conditions. The inversion process ends when the result after correcting m meets the termination condition of step (9), and the final model vector is output.

[0089] (8) Update the model vector: m = m + Δm; recalculate the three-dimensional forward response in the frequency domain, obtain the transient response of the magnetic field in the time domain by converting the frequency domain response to the time domain response (GS transformation), and calculate the root mean square error Rms2;

[0090] (9) Calculate the root mean square error ΔRms: ΔRms=|Rms1-Rms2| / Rms1, and determine whether the root mean square error is less than the minimum fitting error Rms, or whether the number of iterations i is less than the maximum number of iterations I. max If none of them are satisfied, let i = i + 1, take the updated model in (8) as the current model, and return to step (5); otherwise, execute step (10);

[0091] (10) Terminate the iteration and output the current model m.

[0092] The corresponding software module for the electric source ground-space transient electromagnetic three-dimensional focusing inversion imaging, which is supported by the above-mentioned methods and technologies, includes: a system function module and a bottom support module. The system function module includes a data file management module, a preprocessing module, a forward modeling library module, and an electromagnetic response value search and mapping module. The bottom support module includes a data file I / O module, an embedded database module, and a general mathematical library module. The bottom support module provides general function functions to the system function module.

[0093] Example 1: Three-dimensional focusing inversion of rule model data

[0094] The theoretical model has a background resistivity of 100 Ω·m, into which a low-resistivity anomaly of 20 Ω·m is embedded. This anomaly measures 300 m × 300 m × 100 m, with its top surface buried at a depth of 200 m. A schematic diagram of the grid partitioning of the observation area (excluding extended regions) is shown below. Figure 2 As shown in the figure, (a) represents the measuring points on the XY cross-section, and (b) represents the measuring points on the XZ (YZ) cross-section; the light color represents 100 Ω·m, and the dark color represents 20 Ω·m. The total grid is divided into 23×23×21, with 81 measuring points in total, 9 grids each in the horizontal x and y directions of the target area. In the extended area, 7 grids are set in the positive and negative x and y directions respectively, extending to more than 5km. In the air, 6 grids are set to extend to more than 3km, and in the underground, 6 grids are set to extend to less than 6km. The designed line source length is 1000m, the current is 1A, the origin of the coordinate system is used as the center point of the line source, the minimum offset distance from the line source to the observation area is 300m, and the receiving coil heights are 0m on the ground (ground transient electromagnetic) and 50m in the air (semi-airborne transient electromagnetic).

[0095] Figure 3 The figures show three-dimensional inversion diagrams for semi-airborne transient electromagnetic L1 and L2 canonical models. (a) is an L1 canonical view with the surface removed from the top view; (b) is an L2 canonical view with the surface removed from the top view; (c) is an L1 canonical view with one corner removed from the side view; and (d) is an L2 canonical view with one corner removed from the side view. The figures demonstrate that the three-dimensional focused inversion can restore the true electrical characteristics of the model.

[0096] Example 2: Three-dimensional focusing inversion of measured mineral data in a certain area

[0097] The ground transient electromagnetic survey lines in this area are generally laid out in a regular pattern, with a line source length of 2km, a total of 5 survey lines, a line spacing of 200m, and a measurement point spacing of 80m. The target area of ​​the 3D model is divided into regular hexahedrons according to the survey line layout, with an outer extension of 10km.

[0098] Figure 4 The table shows the characteristics of the natural electric field potential curve, the known ore body, and the 3D focused inversion results. The top part shows the natural potential curve, the middle part shows the geological exploration results, and the bottom part shows the L1 regular 3D inversion results. Figure 5 The image shows a comparison of the L1 and L2 canonical 3D inversion results of transient electromagnetic data from the grounding wire source in the Mo mine. Left: L1 canonical, right: L2 canonical. The image demonstrates that the inversion algorithm proposed in this invention yields results in a horizontal position consistent with the known ore body location and traditional electrical resistivity tomography (EDT) methods.

[0099] In summary, the natural electric field data also show a good correspondence with the three-dimensional inversion results of transient electromagnetic data.

[0100] Those skilled in the art will recognize that the embodiments described herein are intended to help the reader understand the principles of the invention, and should be understood that the scope of protection of the invention is not limited to such specific statements and embodiments. Those skilled in the art can make various other specific modifications and combinations based on the technical teachings disclosed in this invention without departing from the spirit of the invention, and these modifications and combinations are still within the scope of protection of this invention.

Claims

1. A three-dimensional focusing inversion method for ground-space transient electromagnetic sources, characterized in that, Includes the following steps: (1) Select the work area data file, including the parameter file and the observation data file d; (2) Read the parameter file, including: receiving coil area, response time, wire source length, offset, peak current, background resistivity, and response time; (3) Set the initial model according to the parameters, including: 3D mesh generation and physical property filling; (4) Set the inversion termination condition, including: maximum number of iterations. Minimum fitting error Rms; Initialize iteration count i=1; Initialize model :Model The elements in the model are resistivity, if the initial model If no prior information is included, an arbitrary resistivity value is set; if prior information is included, the resistivity value is set according to the prior information. (5) Calculate the frequency domain three-dimensional forward response, and then use the adjoint forward method to solve the frequency domain rate Jacobian matrix and update the data item Hessian matrix; obtain the time domain magnetic field transient response by converting the frequency domain response to the time domain response, and calculate the root mean square error Rms1; the specific implementation method of the three-dimensional forward response is as follows: rewrite the frequency domain electric field double curl equation into a weighted residual equation by the Galerkin method, and solve the equation system by the bistable conjugate gradient method of incomplete LU decomposition to obtain the electric field; finally, use the difference method to calculate the auxiliary magnetic field; The weighted residual equation is: （1）； in, This represents the electric field in the frequency domain to be solved. This indicates the curl of the electric field in the frequency domain; Indicates the background magnetic field. Represents a vector interpolation function. Indicates curl; Angular frequency, The permeability of free space, Indicates electrical conductivity. Represents the volume integral. Indicates residual conductivity. Indicates the background electric field; (6) The frequency domain Jacobian matrix is ​​converted into the time domain Jacobian matrix by GS transformation, and the L1 canonical inversion objective function is obtained by Gauss-Newton inversion method; (7) The L1 regularization term in the L1 regularization inversion objective function is calculated by using the iterative reweighted least squares method, and then the model correction is obtained by solving the linear equation system using the conjugate gradient method. ; (8) Update model vectors: The frequency domain three-dimensional forward response is calculated again, and the time domain magnetic field transient response is obtained by converting the frequency domain response to the time domain response. The root mean square error Rms2 is then calculated. (9) Calculate the root mean square error ΔRms: ΔRms = |Rms1 - Rms2| / Rms1, and determine whether the root mean square error ΔRms is less than the minimum fitting error Rms, or whether the number of iterations i is less than the maximum number of iterations. If none of them are satisfied, then let If the updated model in (8) is used as the current model, return to step (5); otherwise, execute step (10). (10) Terminate the iteration and output the current model m.

2. The method for three-dimensional focusing inversion of electric source ground-space transient electromagnetic fields according to claim 1, characterized in that, The specific implementation method of the three-dimensional forward response is as follows: The first term on the left side of equation (1) is expressed as: （2）； It is a coefficient matrix; The second term on the left is represented as: （3）； This is a first-order interpolation function. The superscript 'e' represents a cell in the 3D mesh, and the subscripts 'p' and 'q' represent the edge numbers. It is a coefficient matrix; The right-hand side is represented as: （4）； Sorting and assembling all the element meshes results in a large sparse linear system of equations. Solving this system yields the edge electric field components of the entire mesh, expressed as follows: （5）； in The overall stiffness matrix is ​​composed of the left-hand terms of all elements and is obtained through formulas (2) and (3). Let be the electric field coefficient to be solved. It consists of the right-hand items of all cells, obtained through formula 4; Finally, through the... The frequency domain magnetic field is obtained using a cubic spline interpolation algorithm. , frequency domain magnetic field and background magnetic field The results are summed, and then the time-domain magnetic field transient response dBt / dt is obtained through the GS transform. The response dBt / dt is denoted as... .

3. The method according to claim 2, wherein, The L1 regularized inversion objective function is expressed as follows: （6）； Equation (6) can be expanded as follows: （7）； The first term on the right-hand side of equation (7) is the data item, and the second term is the L1 regularization term, where N m The number of underground grid cells. The j-th element of the regularization vector; This is a weighted matrix of data, where the superscript T denotes the matrix transpose. As a regularization factor, For observation data files, This is the reference model vector.

4. The method for three-dimensional focusing inversion of electric source ground-space transient electromagnetic fields according to claim 3, characterized in that, The iterative reweighted least squares method is as follows: In equation (7), the L1 regularization term is approximately represented by matrix expansion as follows: （8）； Weighted matrix elements in It changes with each iteration and is related to the result of the previous iteration; that is, it is a reweighted relationship after each iteration; it can be expressed in the following form: （9）； in It is a small value; Taking the derivative of equation (7) with respect to the model increment, we get: （10）； The square brackets on the left side of formula (10) represent the Hessian matrix of the data items. It is a Jacobian matrix; Let m be the constraint matrix of model m; The model correction amount is obtained through formula (10). .

Images