Method for 3d inversion of airborne electromagnetic data based on compressed sensing and preconditioned stochastic gradient

Abstract

The present application belongs to the technical field of geophysical electromagnetic inversion, and specifically relates to an airborne electromagnetic three-dimensional inversion method based on compressed sensing and preconditioned random gradient, and comprises the following steps: step one, screening of available data in airborne electromagnetic measured data for inversion; step two, random sampling of measuring points by using a Poisson disc sampling method, forward calculation of the sampling points by using a finite volume method to obtain predicted data of the random measuring points, and reconstruction of the predicted data of all measuring points by using a compressed sensing method; and step three, data fitting of the measured data of step one and the predicted data of step two, and construction of an objective function of airborne electromagnetic three-dimensional regularization inversion in combination with model roughness, wherein the convergence speed can be comparable to that of traditional full-batch data three-dimensional inversion, and the inversion calculation efficiency is greatly improved while ensuring inversion accuracy.

Description

Technical Field

[0001] This invention relates to the field of geophysical electromagnetic inversion technology, specifically to an airborne electromagnetic three-dimensional inversion method based on compressed sensing and preconditional stochastic gradient. Background Technology

[0002] Airborne electromagnetic exploration (AEIA) is a geophysical exploration technique based on platforms such as fixed-wing aircraft, helicopters, or drones. It utilizes airborne or suspended ungrounded wires as magnetic emission sources and employs magnetic sensors to receive underground electromagnetic signals for exploration purposes. It features rapid sampling, high exploration efficiency, and eliminates the need for ground personnel. These advantages of AEIA perfectly solve the problem of exploration being impossible in uninhabited areas; therefore, its development is an inevitable trend.

[0003] The purpose of airborne electromagnetic inversion is to infer the true geoelectric model of underground space from the electromagnetic response received on air, thereby assisting geological interpreters in obtaining more accurate geological structure judgments. Extensive research has been conducted by scholars both domestically and internationally on airborne electromagnetic methods. Due to the massive amount of airborne electromagnetic data, data inversion and interpretation methods typically employ imaging or one-dimensional inversion. However, imaging accuracy is low, and one-dimensional inversion has significant limitations, currently unable to meet the requirements for characterizing fine underground structures. Therefore, various optimization algorithms have been applied to the three-dimensional airborne electromagnetic inversion problem. Currently popular algorithms include: the Gauss-Newton method (Mackie and Madden, 1993; Newman and Alumbaugh, 2000; Sasaki et al., 2013), the quasi-Newton method (Bae et al., 2012; Liu et al., 2013), and the nonlinear conjugate gradient method (Liu et al., 2013; Kamm and Pedersen, 2014), etc.

[0004] Traditional 3D inversion algorithms require calculations on all measurement points in the survey or target area, along with extensive adjoint forward modeling to obtain the sensitivity matrix, resulting in high time consumption and low efficiency. With the rise of machine learning theory, stochastic gradient descent has gained attention. This method randomly selects one or a batch of samples for training, replacing the traditional full-batch calculation, thus significantly improving training efficiency while achieving the training objective. However, due to the limited data volume, the gradient direction is not optimal, leading to slow inversion convergence and the introduction of gradient noise, potentially causing repeated deviations from the optimal solution. These problems can be effectively mitigated through preconditioning techniques. Therefore, a time-domain airborne electromagnetic 3D inversion method based on compressed sensing and preconditioned stochastic gradients can achieve convergence speeds comparable to traditional 3D inversion methods while significantly improving computational efficiency while maintaining inversion accuracy. Summary of the Invention

[0005] The purpose of this section is to outline some aspects of the embodiments of the present invention and to briefly describe some preferred embodiments. Simplifications or omissions may be made in this section, as well as in the abstract and title of this application, to avoid obscuring the purpose of these documents; however, such simplifications or omissions should not be construed as limiting the scope of the invention.

[0006] In view of the problems existing in the above and / or prior art, the present invention is proposed.

[0007] Therefore, the purpose of this invention is to provide an airborne electromagnetic three-dimensional inversion method based on compressed sensing and preconditional stochastic gradients, which has a convergence speed comparable to traditional full-batch data three-dimensional inversion, while greatly improving the inversion calculation efficiency while ensuring inversion accuracy.

[0008] To address the aforementioned technical problems, according to one aspect of the present invention, the present invention provides the following technical solution:

[0009] The airborne electromagnetic 3D inversion method based on compressed sensing and preconditional stochastic gradients includes the following steps:

[0010] Step 1: Select usable data from airborne electromagnetic field measurements for inversion;

[0011] Step 2: Randomly sample the measurement points using the Poisson disk sampling method, and perform forward modeling on the sampled points using the finite volume method to obtain the predicted data of the random measurement points. Then, reconstruct the predicted data of all measurement points using the compressed sensing method.

[0012] Step 3: Fit the measured data from Step 1 and the predicted data from Step 2, and construct the objective function for the three-dimensional regularized inversion of airborne electromagnetic systems by combining the model roughness.

[0013] Step 4: Differentiate both sides of the forward equation and use the adjoint forward model to calculate the sensitivity matrix;

[0014] Step 5: The gradient, i.e. the first derivative of the objective function, can be calculated by multiplying the transpose of the sensitivity matrix with the data fitting difference and adding the model conductivity parameter. The Hessian matrix, i.e. the second derivative of the objective function, can be obtained by combining the sensitivity matrix with the model covariance and the data covariance.

[0015] Step 6: Establish preconditional operators, form preconditional stochastic gradient-Gauss-Newton method inversion equations, calculate model update amount, and obtain new iterative model;

[0016] Step 7: Repeat steps 2 to 6 until the maximum number of iterations is reached or the termination condition is met, to obtain the final inversion result;

[0017] In step six, the inversion equation of the Gauss-Newton method with preconditioned stochastic gradients is:

[0018]

[0019] in This represents the stochastic gradient calculated from randomly sampled measurement points. As an intermediate variable, Let m represent the stochastic gradient calculated from random points of the Hessian matrix, where m = (m1, m2, ... m). M ) T , is the conductivity parameter vector of the M-dimensional model, and P is the preconditioner operator, a diagonal and positive definite matrix, which can be expressed as:

[0020] P ij =α / β+u ij (18)

[0021] Where β is the sampling rate, α is the weighting coefficient, and u ij For noise;

[0022] Finally, the model update amount is obtained by solving the inversion equation (17).

[0023] As a preferred embodiment of the airborne electromagnetic three-dimensional inversion method based on compressed sensing and preconditional stochastic gradient described in this invention, in step one, after removing atmospheric noise, background field, and motion noise from the airborne electromagnetic data, usable airborne electromagnetic data is selected for inversion.

[0024] As a preferred embodiment of the airborne electromagnetic three-dimensional inversion method based on compressed sensing and preconditional stochastic gradient described in this invention, in step two, the Poisson disk sampling divides the region of all measurement points into several non-overlapping circular regions of the same radius. Only one point is taken within each circle, completing random sampling of all measurement points. For the randomly undersampled points, three-dimensional forward modeling is performed using the finite volume method. The double curl equation of the electric field can be expressed as:

[0025]

[0026] Among them, E b and E s Let represent the electric field intensity of the background field and the secondary field, respectively; t represent time; and σ, ε, and μ represent the conductivity, permittivity, and permeability, respectively. This indicates the Hamiltonian operator, where the second time derivative term is caused by the displacement current. This term is much smaller than the conduction current, so it is ignored.

[0027] Discretize equation (1) using the regular grid finite volume method, and then integrate over each control volume to obtain the integral form of equation (1):

[0028]

[0029] in, Let V represent the surface area of ​​any control volume V, and n represent the normals of each surface of the control volume.

[0030] Equation (2) is discretized in time using the unconditionally stable back-electron Euler method, considering all emission sources, all components, and all time channels obtained by random sampling, and then rearranged into a unified form:

[0031]

[0032] in

[0033]

[0034] Where Nc is the number of measurement points randomly sampled in airborne electromagnetic space, Tn is the number of time channels used in the program's internal calculations, and C and D are parameters related to the grid size. and This represents the coefficients, electric field values, and right-hand side terms obtained from randomly sampled measurement points;

[0035] After the forward response calculation for random measurement points is completed, compressed sensing is used to reconstruct the data from all measurement points. Reconstruction is an optimization process, i.e., solving for:

[0036]

[0037] Where ε represents data noise, S is the measurement result obtained from the sampling matrix, and ψ is the sparse transformation matrix; that is, by minimizing the L1 norm, we continuously search for a more sparse matrix. Then, the forward modeling results of all measurement points in the survey area are obtained by using sparse inverse transform.

[0038] As a preferred embodiment of the airborne electromagnetic three-dimensional inversion method based on compressed sensing and preconditional stochastic gradients described in this invention, wherein: in step three, the objective function for regularized inversion is constructed based on airborne electromagnetic observation data and predicted electromagnetic data obtained through forward modeling:

[0039]

[0040] in and These are the data fitting term and the model constraint term, respectively, d = (d1, d2, ..., d...). N ) T It is an N-dimensional data vector, and N = nst × n, where nst and n represent the number of airborne electromagnetic measurement points and the number of time channels, respectively, and γ is a regularization factor;

[0041] For data fitting terms, the specific form is:

[0042]

[0043] Where f(m) represents the forward response corresponding to model parameter m, C d The covariance matrix of airborne electromagnetic data is represented by T, where T is the transpose operator.

[0044] These are model constraint terms, specifically in the form of:

[0045]

[0046] Where m0 represents the initial parameters of the model, C m It is the model covariance used to constrain the changes in model parameters during the inversion iteration process.

[0047] As a preferred embodiment of the airborne electromagnetic three-dimensional inversion method based on compressed sensing and preconditional stochastic gradient described in this invention, wherein: in step four, the calculated sensitivity matrix is ​​missing, and the process of calculating the sensitivity matrix is ​​as follows:

[0048] First, take the differential of both sides of the forward equation (3) with respect to the model parameter m, and after a simple transformation, obtain:

[0049]

[0050] The partial derivatives of the model response with respect to the model parameters are used as the sensitivity matrix J:

[0051]

[0052] Among them, F b and F s These represent the background field response and the anomalous field response generated by the anomalous body, respectively.

[0053] When the background field is a one-dimensional half-space or a layered medium response, it does not change with the conductivity at any point on the ground, therefore:

[0054]

[0055] in, This represents the forward coefficient matrix corresponding to the randomly sampled data. L is an intermediate parameter consisting of a coefficient matrix and a secondary electric field. It includes the interpolation process that directly obtains the measured time channel response at the receiver from the electric field at all times and locations obtained by forward modeling. Compared to the sensitivity matrix J calculated by traditional methods, columns are missing.

[0056] Denote intermediate variables Therefore:

[0057]

[0058] Therefore, V can be obtained by solving the system of adjoint forward equations (13), and the transpose of the sensitivity matrix can be obtained by substituting V into equation (12).

[0059] As a preferred embodiment of the airborne electromagnetic three-dimensional inversion method based on compressed sensing and preconditional stochastic gradient described in this invention, wherein: in step five, the gradient is obtained by taking the first derivative of the airborne electromagnetic regularized inversion objective function with respect to the model parameters:

[0060]

[0061] in, It is an intermediate variable, and

[0062] The Hessian matrix can be obtained by taking the second derivative of the airborne electromagnetic regularization inversion objective function with respect to the model parameters:

[0063]

[0064] The first term of the Hessian matrix corresponds to a full matrix and is not positive definite, making its solution complex. Therefore, in the Gauss-Newton inversion method, this term can be ignored, i.e., the expression for the Hessian matrix is:

[0065]

[0066] As a preferred embodiment of the airborne electromagnetic three-dimensional inversion method based on compressed sensing and preconditional stochastic gradient described in this invention, in step seven, if the iteration termination condition is not met during the iteration process, it is determined whether the reduction in the root mean square error of the data fitting between two adjacent iterations is less than the reduction threshold. If it is less than the reduction threshold, the regularization factor is updated and the iteration calculation continues.

[0067] As a preferred embodiment of the airborne electromagnetic three-dimensional inversion method based on compressed sensing and preconditional stochastic gradient described in this invention, in step seven, the iteration termination conditions include: (1) the root mean square error of the data fitting is less than the error threshold; (2) the airborne electromagnetic regularized inversion objective function converges; (3) the preset maximum number of iterations is reached; the iteration terminates as long as any one of the iteration termination conditions is met.

[0068] Compared with existing technologies, the present invention provides a fast three-dimensional inversion method for airborne electromagnetic fields based on compressed sensing and preconditional stochastic gradients in the time domain. The numerical simulation algorithm used for inversion is a fast forward modeling algorithm based on the three-dimensional time domain finite volume method using compressed sensing. This algorithm can simulate electromagnetic field signals under arbitrary conductivity structures and has the advantages of fast speed, high accuracy in simulating complex geological and topographical structures, and memory saving. At the same time, the inversion method is based on the Gauss-Newton method of preconditional stochastic gradients and proposes effective preconditional operator estimation methods, regularization constraint methods, and sensitivity solution methods. These methods can greatly improve the computational efficiency of three-dimensional inversion while ensuring accuracy and convergence speed. Attached Figure Description

[0069] To more clearly illustrate the technical solutions of the embodiments of the present invention, the present invention will be described in detail below with reference to the accompanying drawings and detailed embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort. Wherein:

[0070] Figure 1 This is a flowchart of a fast three-dimensional inversion method for time-domain airborne electromagnetic systems based on compressed sensing and preconditional stochastic gradients provided in an embodiment of the present invention.

[0071] Figure 2 This is a uniform half-space model provided in the embodiments of the present invention, wherein (a) is a schematic diagram of the uniform half-space model, the black dot in the figure is the location of the airborne electromagnetic receiving point, which is 40m away from the ground, and the horizontal line below represents the boundary line of the ground surface; (b) is the shape of the trapezoidal wave current emitted by the forward model VTEM.

[0072] Figure 3 This is a schematic diagram of the accuracy verification of the time-domain airborne electromagnetic three-dimensional finite volume method and the one-dimensional semi-analytical solution of the uniform half-space model provided in the embodiment of the present invention. (a) is a comparison result of the time-domain airborne electromagnetic three-dimensional finite volume method and the one-dimensional semi-analytical solution of the uniform half-space model; (b) is a single-point relative error analysis diagram of different time channels.

[0073] Figure 4 This is a schematic diagram of a single anomaly provided in an embodiment of the present invention, with color markings representing a logarithmic distribution;

[0074] Figure 5 This is an early time-channel reconstruction structure diagram provided by an embodiment of the present invention. The color scales of (a), (b), and (c) are logarithmic distributions, and the color scale of (d) is a linear distribution. Among them, (a) is the forward response of the early time channel; (b) is the forward response diagram obtained by random sampling at a sampling rate of 40%; (c) is the response obtained by compressed sensing reconstruction; and (d) is the single-point relative error between the response obtained by compressed sensing reconstruction and the original response.

[0075] Figure 6 These are schematic diagrams of multiple anomaly models and inversion results of various inversion methods provided in embodiments of the present invention. The color scales represent a logarithmic distribution. Among them, (a) is a schematic diagram of multiple anomaly models; (b) is the Gauss-Newton inversion result of the full batch data; (c) is the Gauss-Newton inversion result based on stochastic gradients; and (d) is the Gauss-Newton inversion result based on preconditional stochastic gradients.

[0076] Figure 7 These are iterative parameter diagrams for three inversion methods provided in this embodiment of the invention, wherein (a) is a diagram showing the RMS change during the inversion process of the three inversion methods; and (b) is a diagram showing the change of the objective function value during the inversion process of the three inversion methods.

[0077] Figure 8 Three-dimensional inversion results based on compressed sensing and preconditional stochastic gradient when using different weighting factors with a 40% sampling rate;

[0078] Figure 9 These are slices of inversion results at different depths based on compressed sensing and preconditional stochastic gradient inversion algorithms with weighting factors of 1.4 and 1.8, respectively;

[0079] Figure 10 This is the mean square error decay curve of the data fitting for the traditional full-batch data inversion 40% sampling preconditional stochastic gradient inversion method in this example. Detailed Implementation

[0080] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

[0081] Many specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways different from those described herein, and those skilled in the art can make similar extensions without departing from the spirit of the invention. Therefore, the invention is not limited to the specific embodiments disclosed below.

[0082] Secondly, the present invention is described in detail with reference to the schematic diagrams. When detailing the embodiments of the present invention, for ease of explanation, the cross-sectional views illustrating the device structure may be partially enlarged, not according to the usual scale. Furthermore, the schematic diagrams are merely examples and should not limit the scope of protection of the present invention. In addition, actual fabrication should include three-dimensional spatial dimensions of length, width, and depth.

[0083] To make the objectives, technical solutions, and advantages of the present invention clearer, the embodiments of the present invention will be described in further detail below with reference to the accompanying drawings.

[0084] This invention provides an airborne electromagnetic 3D inversion method based on compressed sensing and preconditional stochastic gradients. Its convergence speed is comparable to traditional full-batch data 3D inversion, while significantly improving computational efficiency while maintaining inversion accuracy. Please refer to [link to relevant documentation]. Figures 1-7 This includes the following steps:

[0085] S1, Filter available data from airborne electromagnetic field measurement data for inversion;

[0086] After preprocessing the airborne electromagnetic data, such as removing atmospheric noise, background field, and motion noise, usable airborne electromagnetic data are selected for inversion.

[0087] During the acquisition process, airborne electromagnetic data is subject to interference from various factors, resulting in signal distortion and warping, which poses significant difficulties for further inversion and interpretation. Therefore, it is necessary to preprocess the airborne electromagnetic data to obtain usable airborne electromagnetic data.

[0088] S2, the Poisson disk sampling method is used to randomly sample the measurement points, and the sampling points are subjected to three-dimensional forward modeling using the finite volume method to obtain the prediction data of the random measurement points. The compressed sensing method is then used to reconstruct the prediction data of all measurement points.

[0089] Poisson disk sampling divides the region of all measurement points into several non-overlapping circular regions of the same radius. Only one point is selected within each circle, thus completing random sampling of all measurement points. When performing finite volume forward modeling for these random sampling points, the double curl equation of the electric field can be expressed as:

[0090]

[0091] Among them, E b and E s Let represent the electric field intensity of the background field and the secondary field, respectively; t represent time; and σ, ε, and μ represent the conductivity, permittivity, and permeability, respectively. Let represent the Hamiltonian operator. The second time derivative term is the displacement current term, which is much smaller than the conduction current term, so this term is ignored.

[0092] Figure 2 This is a schematic diagram of the uniform half-space model provided in the embodiment. In diagram (a), the black dot represents the location of the airborne electromagnetic receiving point, 40m above the ground. The horizontal line below indicates the ground boundary. The air conductivity is set to 1×10⁻⁶. 8 (a) shows the trapezoidal current shape emitted by the forward VTEM, with an off-time of 7.516 ms. The half-space conductivity is 100 Ω·m.

[0093] Discretize equation (1) using the regular grid finite volume method, and then integrate over each control volume to obtain the integral form of equation (1):

[0094]

[0095] in, Let V represent the surface area of ​​any control volume V, and n represent the normals of each surface of the control volume.

[0096] In this embodiment, a trapezoidal wave is used as the transmitted waveform. Due to the significant electrical difference between air and underground media, if explicit time discretization is used, a sufficiently short time step Δt is required to meet the stability condition, which is very unfavorable for late-stage response calculation. Therefore, this embodiment uses the implicit back-pull Euler method for time discretization, which is unconditionally stable regardless of the time step Δt. That is, the first-order back-pull Euler method is applied to equation (2), considering all emission sources, all components, and all time channels obtained by random sampling, and they are organized into a unified form:

[0097]

[0098] in

[0099]

[0100] Where Nc is the number of measurement points randomly sampled in airborne electromagnetic space, Tn is the number of time channels used in the program's internal calculations, and C and D are parameters related to the grid size. and This represents the coefficients, electric field values, and right-hand side terms obtained from randomly sampled measurement points.

[0101] After the forward response calculation for random measurement points is completed, compressed sensing is used to reconstruct the data from all measurement points. Reconstruction is an optimization process, i.e., solving for:

[0102]

[0103] Where ε represents data noise, S is the measurement result obtained from the sampling matrix, and ψ is the sparse transformation matrix. That is, by minimizing the L1 norm, we continuously search for a more sparse matrix. Finally, the forward response of all measuring points in the survey area is obtained through sparse inverse transformation, which is used to calculate the inverse objective function.

[0104] To verify the accuracy of the forward modeling process and compressed sensing reconstruction, a method was first designed as follows: Figure 2 (a) shows the half-space model. The finite volume method used in this experimental example was employed for calculation, and the results were compared with the one-dimensional semi-analytical solution. The model parameters are as follows: air resistivity is taken as 1 × 10⁻⁶. 8Ω·m, half-space resistivity is 100Ω·m, transmitting coil diameter is 26m, transmitting as Figure 2 (b) shows a trapezoidal current with a flight altitude of 40m. The black dot is the location of the airborne electromagnetic receiving point.

[0105] The verification results of the uniform half-space are as follows Figure 3 As shown, (a) is a comparison of the time-domain airborne electromagnetic three-dimensional finite volume method and the one-dimensional semi-analytical solution for the uniform half-space model; (b) is a single-point relative error analysis diagram for different time channels; from Figure 3 It can be seen that the single-point relative error of both is less than 5% except for the last time trace, indicating that the forward modeling algorithm used in this experimental example has high accuracy.

[0106] To verify the effectiveness of the compression and reconstruction used in this experimental example, we designed the following... Figure 4 The single anomaly model shown, after random sampling, selects an early trace to calculate the forward response of 40% of the randomly sampled measurement points. Then, it reconstructs the entire area data through compression reconstruction. The model parameters are as follows: anomaly resistivity is 1 Ω·m, top surface burial depth is 32m, total thickness is 52m. Half-space resistivity is 100 Ω·m, transmitting coil diameter is 26m, and the transmitting... Figure 2 (b) shows a trapezoidal current, with a flight altitude of 40m and an air resistivity of 1×10⁻⁶. 8 Ω·m.

[0107] The verification results are as follows Figure 5 As shown, (a) is the forward modeling response of the early time trace; (b) is the forward modeling response obtained by random sampling at a 40% sampling rate; (c) is the response obtained by compressed sensing reconstruction; and (d) is the single-point relative error between the response obtained by compressed sensing reconstruction and the original response. Figure 5 It is evident that a random sampling rate of 40% can effectively calculate the three-dimensional forward modeling signals of all measurement points in the survey area, with reconstruction errors all within 4%. This indicates that the compression reconstruction algorithm used in this experimental example has high accuracy when the sampling rate is 40%, and can be used for inversion research.

[0108] S3, fit the measured data selected in S1 and the predicted data in S2, and construct the objective function for the three-dimensional regularized inversion of airborne electromagnetic systems by combining the model roughness.

[0109] Based on airborne electromagnetic observation data and predicted electromagnetic data obtained through forward modeling, the objective function for regularized inversion is constructed as follows:

[0110]

[0111] in and These are the data fitting term and the model constraint term, respectively, m = (m1, m2, ... mM ) T It is an M-dimensional model conductivity parameter vector, d=(d1,d2,…,d N ) T It is an N-dimensional data vector, and N = nst × n, where nst and n represent the number of airborne electromagnetic measurement points and the number of time channels, respectively, and γ is a regularization factor.

[0112] For data fitting terms, the specific form is:

[0113]

[0114] Where f(m) represents the forward response corresponding to model parameter m, C d The covariance matrix of airborne electromagnetic data is represented by T, where T is the transpose operator.

[0115] These are model constraint terms, specifically in the form of:

[0116]

[0117] Where m0 represents the initial parameters of the model, C m It is the model covariance used to constrain the changes in model parameters during the inversion iteration process.

[0118] S4, differentiate both sides of the forward equation, where the partial derivatives of the predicted data with respect to the model parameters are the sensitivity matrix (Jacobi matrix);

[0119] The sensitivity matrix calculated in this step is missing. The process for calculating the sensitivity matrix is ​​as follows:

[0120] First, take the differential of both sides of the forward equation (3) with respect to the model parameter m, and after a simple transformation, obtain:

[0121]

[0122] The partial derivatives of the model response with respect to the model parameters are used as the sensitivity matrix J:

[0123]

[0124] Among them, F b and F s These represent the background field response and the anomalous field response generated by the anomalous body, respectively.

[0125] When the background field is a one-dimensional half-space or a layered medium response, it does not change with the conductivity at any point on the ground, therefore:

[0126]

[0127] in, This represents the forward coefficient matrix corresponding to the randomly sampled data. It consists of intermediate parameters composed of a coefficient matrix and a secondary electric field. L includes the interpolation process that directly obtains the measured time channel response at the receiver from the electric field at all times and locations obtained by forward modeling. In contrast to the sensitivity matrix J calculated by traditional methods, which is missing column by column and padded with zeros, this method provides a more comprehensive approach.

[0128] Denote intermediate variables Therefore:

[0129]

[0130] Therefore, V can be obtained by solving the system of adjoint forward equations (13), and then V can be substituted into equation (12).

[0131] The transpose of the sensitivity matrix can then be obtained.

[0132] S5. The gradient of the objective function can be calculated by multiplying the transpose of the sensitivity matrix with the product of the fitting difference between the predicted data and the observed data, and adding the model conductivity parameter. The Hessian matrix can be obtained by combining the sensitivity matrix with the model covariance and numerical covariance.

[0133] The gradient of the airborne electromagnetic regularization inversion objective function is obtained by taking the first derivative with respect to the model parameters:

[0134]

[0135] in, It is an intermediate variable, and

[0136] The Hessian matrix can be obtained by taking the second derivative of the airborne electromagnetic regularization inversion objective function with respect to the model parameters:

[0137]

[0138] The first term of the Hessian matrix corresponds to a full matrix and is not positive definite, making its solution complex. Therefore, in the Gauss-Newton inversion method, this term can be ignored, i.e., the expression for the Hessian matrix is:

[0139]

[0140] S6. The Gauss-Newton method with preconditional stochastic gradients is used to calculate the model update amount and obtain a new iterative model.

[0141] The inversion equation for the Gaussian-Newton method with preconditioned stochastic gradients is:

[0142]

[0143] in This represents the stochastic gradient calculated from randomly sampled measurement points. It is an intermediate variable, and Let m represent the stochastic gradient calculated from random points of the Hessian matrix, where m = (m1, m2, ... m). M ) T , is the conductivity parameter vector of the M-dimensional model, and P is the preconditioner operator, a diagonal and positive definite matrix, which can be expressed as:

[0144] P ij =α / β+u ij ,i=j (18)

[0145] Where β is the sampling rate, α is the weighting coefficient, and u ij It is noise.

[0146] Finally, the model update amount is obtained by solving the inversion equation (17).

[0147] S7. Repeat S2 to S6 until the maximum number of iterations is reached or the termination condition is met, to obtain the final inversion reference model.

[0148] During the iteration process, if the iteration termination condition is not met, it is determined whether the decrease in the root mean square error of the data fitting between two adjacent iterations is less than the decrease threshold. If it is less than the decrease threshold, the regularization factor is updated and the iteration calculation continues. In this embodiment, the iteration termination condition includes: (1) the root mean square error of the data fitting is less than the error threshold, i.e., RMS < r0; (2) the airborne electromagnetic regularization inversion objective function converges, i.e., ||g k ||≤ε;(3) The iteration terminates as soon as any one of the preset maximum iteration numbers N is met, and the inversion solution that conforms to the real underground geoelectric structure is obtained. Among them, the formula for calculating RMS is:

[0149]

[0150] Where, ε i The value represents the data error, and n represents the number of measurement points.

[0151] During the iteration process, if the iteration termination condition is not met, it is determined whether the decrease in the root mean square error (RMS) of the data fitting between two adjacent iterations is less than the decrease threshold χ, where χ > 0 and is generally taken as 2%. If it is less than the decrease threshold, the regularization factor is updated and the iteration calculation continues.

[0152] This embodiment provides an inversion example to verify the correctness and efficiency of the inversion algorithm. Figure 6These are schematic diagrams of multiple anomaly models and inversion results provided in embodiments of the present invention. (a) shows schematic diagrams of multiple anomaly models; (b) shows the Gauss-Newton inversion results for the full batch of data; (c) shows the Gauss-Newton inversion results based on stochastic gradients; and (d) shows the Gauss-Newton inversion results based on preconditioned stochastic gradients. Figure 6 As shown in (a), multiple anomalies are distributed in different locations, each with a different shape. The resistivity of each anomaly is 1 Ω·m, the top surface is buried at a depth of 80 m, and the total thickness is 120 m. The resistivity of the half-space is 100 Ω·m, the diameter of the transmitting coil is 26 m, and the emission is as follows... Figure 2 (b) shows a trapezoidal current, with a flight altitude of 40m and an air resistivity of 1×10⁻⁶. 8 Ω·m. From Figure 6 As can be seen from (b), (c), and (d), all three inversion methods provide a good reflection of the morphology, location, and resistivity of the target anomaly.

[0153] The RMS and objective function values ​​of the three inversion methods change with the number of iterations as follows: Figure 7 As shown, (a) is a graph showing the RMS change during the inversion process of the three inversion methods; (b) is a graph showing the change of the objective function value during the inversion process of the three inversion methods. Meanwhile, a comparison of the parameters of the three inversion methods is shown in Table 1. Table 1 is a table of inversion iteration parameters for different inversion methods provided in this embodiment of the invention, including the RMS at termination of the three inversion methods, the number of iterations, the time taken, and the efficiency improvement compared to the Gauss-Newton method for full-batch data.

[0154] Inversion method RMS at termination Number of iterations Time / h efficiency / % All data inversion 1.78 4 70 -- SG inversion - 40% sampling 1.82 6 36 48.6 PSG inversion 1.92 4 20 71.4

[0155] Table 1

[0156] Depend on Figure 7 It can be seen that the convergence speed of the Gauss-Newton method based on preconditioned stochastic gradients is comparable to that of the Gauss-Newton method based on full-batch data. Table 1 shows that in the Gauss-Newton inversion of full-batch data, the RMS decreased to 1.78 after 4 iterations, with a total time of 70 hours. The Gauss-Newton method based on stochastic gradients requires 6 iterations, with the RMS decreasing to 1.82, but takes 36 hours, showing a certain improvement in efficiency compared to the traditional method. After adding preconditions, only 4 iterations are needed, with the RMS decreasing to 1.92, and the time taken is 20 hours. This indicates a significant improvement in both the number of iterations and the convergence time. Furthermore, the number of iterations is comparable to that of the Gauss-Newton inversion of full-batch data, with a time of approximately 20 hours. Compared to the 70 hours of the Gauss-Newton inversion of full-batch data, the efficiency is improved by approximately 71.4%, demonstrating a significant improvement in 3D inversion efficiency.

[0157] Example 2 is an experiment on the inversion of data from a single well-conducting anomaly model underground, based on the above implementation process. Figure 8The results are 3D inversion based on compressed sensing and preconditional stochastic gradient descent when using different weighting factors with a 40% sampling rate. The black boxes represent the locations of the true models. Figure 8 It can be seen that the preconditioning inversion algorithm can effectively recover the location and conductivity of the subsurface model when the weighting factors are 0.6, 1.0, and 1.6, verifying the effectiveness of the method. Table 2 compares the time consumption of compressed sensing and preconditioning stochastic gradient inversion with traditional full-data 3D inversion. It can be seen that the full-batch Gauss-Newton inversion takes 55.6 hours, while the weighting factor of 1.6 takes only 12.3 hours, improving efficiency by 78%.

[0158]

[0159] Table 2

[0160] Example 3 shows the results of airborne electromagnetic three-dimensional inversion of complex geoelectric structures based on the above implementation process. Figure 9 These are slices of inversion results at different depths, based on compressed sensing and preconditioned stochastic gradient inversion algorithms with weighting factors of 1.4 and 1.8 at a 40% sampling rate. The black boxes in the figures represent the actual target body locations. It can be seen that at depths of 50m and 80m, the inversion results can well recover shallow anomalies. At a depth of 150m, due to the small anomaly scale and large burial depth resulting in weak anomaly signals, the inversion results can still provide a good characterization of some large-scale anomalies. The overall inversion can reflect the main distribution, location, and conductivity of the anomalies. Figure 10 This example shows the mean square error decay curve of the data fitting for the traditional full-batch data inversion method with 40% sampling preconditional stochastic gradient inversion. It can be seen that 40% sampling achieves a convergence speed comparable to traditional full-data inversion, but the computational cost per iteration is only about 40% of the traditional method, significantly reducing the computational load of 3D inversion and improving the inversion speed. In this example, full-data inversion took 252 hours, and preconditional stochastic gradient inversion took 76 hours, resulting in an overall efficiency improvement of approximately 70%.

[0161] Although the present invention has been described above with reference to embodiments, various modifications can be made and components can be replaced with equivalents without departing from the scope of the invention. In particular, as long as there is no structural conflict, the features in the disclosed embodiments can be combined with each other in any manner. The lack of an exhaustive description of these combinations in this specification is merely for the sake of brevity and resource conservation. Therefore, the present invention is not limited to the specific embodiments disclosed herein, but includes all technical solutions falling within the scope of the claims.

Claims

1. A method for airborne electromagnetic three-dimensional inversion based on compressed sensing and preconditional stochastic gradients, characterized in that, Includes the following steps: Step 1: Select usable data from airborne electromagnetic field measurements for inversion; Step 2: Randomly sample the measurement points using the Poisson disk sampling method, and perform forward modeling on the sampled points using the finite volume method to obtain the predicted data of the random measurement points. Then, reconstruct the predicted data of all measurement points using the compressed sensing method. Step 3: Fit the measured data from Step 1 and the predicted data from Step 2, and construct the objective function for the three-dimensional regularized inversion of airborne electromagnetic systems by combining the model roughness. Step 4: Differentiate both sides of the forward equation and use the adjoint forward model to calculate the sensitivity matrix; Step 5: The gradient, i.e. the first derivative of the objective function, can be calculated by multiplying the transpose of the sensitivity matrix with the data fitting difference and adding the model conductivity parameter. The Hessian matrix, i.e. the second derivative of the objective function, can be obtained by combining the sensitivity matrix with the model covariance and the data covariance. Step 6: Establish preconditional operators, form preconditional stochastic gradient-Gauss-Newton method inversion equations, calculate model update amount, and obtain new iterative model; Step 7: Repeat steps 2 to 6 until the maximum number of iterations is reached or the termination condition is met, to obtain the final inversion result; In step six, the inversion equation of the Gauss-Newton method with preconditioned stochastic gradients is: in This represents the stochastic gradient calculated from randomly sampled measurement points. As an intermediate variable, Let m represent the stochastic gradient calculated from random points of the Hessian matrix, where m = (m1, m2, ... m). M ) T , is the conductivity parameter vector of the M-dimensional model, and P is the preconditioner operator, a diagonal and positive definite matrix, which can be expressed as: P ij =a / b+u ij (18) Where β is the sampling rate, α is the weighting coefficient, and u ij For noise; Finally, the model update amount is obtained by solving the inversion equation (17).

2. The airborne electromagnetic three-dimensional inversion method based on compressed sensing and preconditional stochastic gradient as described in claim 1, characterized in that, In step one, after removing atmospheric noise, background field, and motion noise from the airborne electromagnetic data, usable airborne electromagnetic data is selected for inversion.

3. The airborne electromagnetic three-dimensional inversion method based on compressed sensing and preconditional stochastic gradient as described in claim 2, characterized in that, In step two, the Poisson disk sampling divides the region of all measuring points into several non-overlapping circular regions of the same radius. Only one point is taken in each circle, completing random sampling of all measuring points. For the randomly undersampled points, three-dimensional forward modeling is performed using the finite volume method. The double curl equation of the electric field can be expressed as: Among them, E b and E s Let represent the electric field intensity of the background field and the secondary field, respectively; t represent time; and σ, ε, and μ represent the conductivity, permittivity, and permeability, respectively. This indicates the Hamiltonian operator, where the second time derivative term is caused by the displacement current. This term is much smaller than the conduction current, so it is ignored. Discretize equation (1) using the regular grid finite volume method, and then integrate over each control volume to obtain the integral form of equation (1): in, Let V represent the surface area of ​​any control volume V, and n represent the normals of each surface of the control volume. Equation (2) is discretized in time using the unconditionally stable back-electron Euler method, considering all emission sources, all components, and all time channels obtained by random sampling, and then rearranged into a unified form: in Where Nc is the number of measurement points randomly sampled in airborne electromagnetic space, Tn is the number of time channels used in the program's internal calculations, and C and D are parameters related to the grid size. and This represents the coefficients, electric field values, and right-hand side terms obtained from randomly sampled measurement points; After the forward response calculation for random measurement points is completed, compressed sensing is used to reconstruct the data from all measurement points. Reconstruction is an optimization process, i.e., solving for: Where ε represents data noise, S is the measurement result obtained from the sampling matrix, and ψ is the sparse transformation matrix; that is, by minimizing the L1 norm, we continuously search for a more sparse matrix. Then, the forward modeling results of all measurement points in the survey area are obtained by using sparse inverse transform.

4. The airborne electromagnetic three-dimensional inversion method based on compressed sensing and preconditional stochastic gradient as described in claim 3, characterized in that, In step three, the objective function for regularized inversion is constructed based on airborne electromagnetic observation data and predicted electromagnetic data obtained through forward modeling: in and These are the data fitting term and the model constraint term, respectively, d = (d1, d2, ..., d...). N ) T It is an N-dimensional data vector, and N = nst × n, where nst and n represent the number of airborne electromagnetic measurement points and the number of time channels, respectively, and γ is a regularization factor; For data fitting terms, the specific form is: Where f(m) represents the forward response corresponding to model parameter m, C d The covariance matrix of airborne electromagnetic data is represented by T, where T is the transpose operator. These are model constraint terms, specifically in the form of: Where m0 represents the initial parameters of the model, C m It is the model covariance used to constrain the changes in model parameters during the inversion iteration process.

5. The airborne electromagnetic three-dimensional inversion method based on compressed sensing and preconditional stochastic gradient as described in claim 4, characterized in that, In step four, the calculated sensitivity matrix is ​​missing. The process for calculating the sensitivity matrix is ​​as follows: First, take the differential of both sides of the forward equation (3) with respect to the model parameter m, and after a simple transformation, obtain: The partial derivatives of the model response with respect to the model parameters are used as the sensitivity matrix J: Among them, F b and F s These represent the background field response and the anomalous field response generated by the anomalous body, respectively. When the background field is a one-dimensional half-space or a layered medium response, it does not change with the conductivity at any point on the ground, therefore: in, This represents the forward coefficient matrix corresponding to the randomly sampled data. L is an intermediate parameter consisting of a coefficient matrix and a secondary electric field. It includes the interpolation process that directly obtains the measured time channel response at the receiver from the electric field at all times and locations obtained by forward modeling. Compared to the sensitivity matrix J calculated by traditional methods, columns are missing. Denote intermediate variables Therefore: Therefore, V can be obtained by solving the system of adjoint forward equations (13), and the transpose of the sensitivity matrix can be obtained by substituting V into equation (12).

6. The airborne electromagnetic three-dimensional inversion method based on compressed sensing and preconditional stochastic gradient as described in claim 5, characterized in that, In step five, the gradient of the airborne electromagnetic regularization inversion objective function is obtained by taking the first derivative with respect to the model parameters: in, It is an intermediate variable, and The Hessian matrix can be obtained by taking the second derivative of the airborne electromagnetic regularization inversion objective function with respect to the model parameters: The first term of the Hessian matrix corresponds to a full matrix and is not positive definite, making its solution complex. Therefore, in the Gauss-Newton inversion method, this term can be ignored, i.e., the expression for the Hessian matrix is:

7. The airborne electromagnetic three-dimensional inversion method based on compressed sensing and preconditional stochastic gradient as described in claim 1, characterized in that, In step seven, if the iteration termination condition is not met during the iteration process, it is determined whether the reduction in the root mean square error of the data fitting between two adjacent iterations is less than the reduction threshold. If it is less than the reduction threshold, the regularization factor is updated and the iteration calculation continues.

8. The airborne electromagnetic three-dimensional inversion method based on compressed sensing and preconditional stochastic gradient as described in claim 1, characterized in that, In step seven, the iteration termination conditions include: (1) the root mean square error of the data fitting is less than the error threshold; (2) the objective function of the aerospace electromagnetic regularization inversion converges; (3) the preset maximum number of iterations is reached; the iteration terminates as long as any one of the iteration termination conditions is met.

Images