Three-dimensional magnetotelluric inversion method based on L1 norm and adaptive moment estimation algorithm

Abstract

The invention relates to the technical field of earth electromagnetism, in particular to a three-dimensional magnetotelluric inversion method based on an L1 norm and an adaptive moment estimation algorithm. Collecting magnetotelluric data of the target area as observation data; dividing the non-structural tetrahedral mesh of the target area to establish an initial resistivity model; performing forward modeling calculation according to the current resistivity model to obtain predicted electromagnetic response data of the target area; constructing an inversion objective function according to the data fitting term based on the L1 norm and the model constraint term based on the L2 norm; calculating a data fitting difference between the observation data and the predicted electromagnetic response data; if the data fitting difference does not meet the set condition, updating the current resistivity model by adopting an ADAM algorithm; and iterating in sequence until the data fitting difference meets a set condition, and outputting the resistivity model at the moment as an inversion result. According to the method, the influence of a big data fitting error in inversion can be effectively suppressed, and a higher inversion resolution is obtained.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of geoelectricity, and particularly to a three-dimensional magnetotelluric inversion method based on L1 norm and adaptive matrix estimation algorithm. BACKGROUND

[0002] Three-dimensional magnetotelluric (MT) inversion technology has developed into an important means for detecting deep crust and upper mantle structure. The existing mainstream methods include: nonlinear conjugate gradient method (NLCG): with good memory control ability, suitable for large-scale three-dimensional inversion, such as ModEM software; quasi-Newton method (QN) and limited memory quasi-Newton method (L-BFGS): with fast convergence, suitable for large-scale parallel inversion; Gauss-Newton method (GN): fast convergence speed and high accuracy, suitable for complex geological structure, but large memory consumption; Occam inversion method: based on model smoothing principle, suitable for noise-sensitive data inversion, widely used in one-dimensional and two-dimensional inversion in early stage, and then extended to three-dimensional; regularization techniques including Tikhonov regularization, minimum support, L1 norm, curvelet sparse transform, gradient filtering regularization, etc., effectively alleviate the instability of inverse problem, and enhance the stability and resolution of model.

[0003] Although the existing three-dimensional magnetotelluric inversion technology is relatively mature, there are still deficiencies in the following key aspects: first, the development of existing inversion technology in computational efficiency and memory has bottlenecks, for example, high-precision methods such as Gauss-Newton method require explicit calculation of Jacobian matrix, which has high memory occupation and poor scalability, and is not suitable for large data sets; NLCG is memory friendly, but has slow convergence speed and is easy to fall into local optimum, although L-BFGS improves this problem, but the gradient calculation of non-structured grid still depends on additional optimization; second, the regularization and structure expression ability of the model is limited, and traditional Tikhonov or least squares smoothing method is easy to cause abnormal boundary blur, although minimum support and sparse regularization are improved, but it is still difficult to balance smoothness and abnormal focusing, and the model analysis ability is still insufficient in complex structure (such as fault, volcanic channel) area; the data fitting error measurement method is single, and the existing inversion mainly uses L2 norm to measure data fitting error, which is sensitive to abnormal points, and cannot reflect the data distribution structure, the lack of structure-sensitive measurement method makes the inversion model easy to overfit noise or lose distribution characteristics. SUMMARY

[0004] In order to solve the problems in the prior art, the application provides a three-dimensional magnetotelluric inversion method based on L1 norm and adaptive moment estimation algorithm.

[0005] The application adopts the following technical scheme, the three-dimensional magnetotelluric inversion method based on L1 norm and adaptive moment estimation algorithm comprises: Collect magnetotelluric data of a target region as observation data; Divide the non-structured tetrahedral grid of the target region, and establish an initial resistivity model in the non-structured tetrahedral grid according to the observation data; Forward calculation is performed according to the current resistivity model to obtain predicted electromagnetic response data of the target region; A data fitting term based on L1 norm is established according to the observation data and the predicted electromagnetic response data; A model constraint term based on L2 norm is established according to the initial resistivity model and the current resistivity model; An inversion objective function is constructed according to the data fitting term based on L1 norm and the model constraint term based on L2 norm; The data fitting difference between the observation data and the predicted electromagnetic response data under the current resistivity model is calculated; If the data fitting difference does not meet the set condition, the ADAM algorithm is used to update the current resistivity model; The data fitting difference between the observation data and the predicted electromagnetic response data corresponding to the updated resistivity model is calculated, and iteration is performed in sequence until the data fitting difference meets the set condition, and the resistivity model at this time is output as the inversion result.

[0006] Further, the predicted electromagnetic response data of the target region is obtained by forward calculation according to the current resistivity model, specifically: The frequency domain electric field double spinor equation of the target region is established according to the current resistivity model using Maxwell's equation set; The frequency domain electric field of each tetrahedral element in the non-structured tetrahedral grid is interpolated using a vector basis function; The mass matrix and the stiffness matrix of each tetrahedral element are constructed based on the Galerkin method; The linear equations of the frequency domain electric field are established according to the mass matrix and the stiffness matrix of all the tetrahedral elements; The linear equations are solved by using a direct solver to obtain the frequency domain electric field value of the target region; The magnetic field of the target region is calculated according to the frequency domain electric field value of the target region through the Faraday electromagnetic induction law; The predicted electromagnetic response data of the target region are obtained according to the frequency domain electric field value and the magnetic field of the target region.

[0007] Further, the inversion objective function is constructed according to the data fitting term based on the L1 norm and the model constraint term based on the L2 norm, and is represented as: ; In the formula, the inversion objective function is represented as, is a covariance matrix of the observation data, is a covariance matrix of the predicted electromagnetic response data, is the observation data, is the predicted electromagnetic response data, is an initial resistivity model, is a current resistivity model, is a trade-off parameter, is the L1 norm, is the L2 norm.

[0008] Further, the ADAM algorithm is used to update the current resistivity model, including: The gradient of the objective function under the current resistivity model is calculated; The first-order momentum and the second-order momentum are calculated according to the gradient of the objective function; The current resistivity model is updated according to the first-order momentum and the second-order momentum to obtain an updated resistivity model.

[0009] Further, when the first-order momentum and the second-order momentum are calculated according to the gradient of the objective function, it further includes: The first-order momentum and the second-order momentum are corrected, and are respectively represented as: ; Wherein, the corrected first-order momentum is represented as, the first-order momentum of the i-th step is represented as, is a decay coefficient of the first-order momentum, represents the number of updates of the current resistivity model; ; Wherein, denotes the corrected second order momentum, denotes the second order momentum of the i-th step, is the decay coefficient of the second order momentum.

[0010] Further, the current resistivity model is updated according to the first order momentum and the second order momentum, denoted as: ; wherein, denotes the current resistivity model, denotes the updated resistivity model, denotes the corrected first order momentum, denotes the corrected second order momentum, is the global learning rate, is a constant to prevent the denominator from being zero.

[0011] Further, the data fitting difference between the observed data and the predicted electromagnetic response data under the current resistivity model is calculated, denoted as: ; wherein, is the data fitting difference, is the covariance matrix of the observed data, is the observed data, is the predicted electromagnetic response data, is the L1 norm, and N is the total data amount.

[0012] Further, the setting condition is that when the data fitting difference is less than a set threshold or does not change after multiple iterations, the iteration is stopped.

[0013] The present application has the beneficial effects that: the present application adopts non-structural grid vector finite element forward, which can significantly reduce the memory consumption, solve the storage bottleneck of Jacobian matrix of Gauss-Newton method, further combine L1 norm and L2 norm to construct the inversion objective function, establish the data fitting term by adopting L1 norm, which can effectively suppress abnormal value interference, reduce the high resistance body analysis error, at the same time, the model constraint term of L2 norm can well couple the geological prior information, the resolution of complex structure such as fault and intrusive body is significantly enhanced, further adopt ADAM algorithm to combine first order and second order momentum to update the inversion model, compared with traditional quasi-Newton algorithm and other methods, the model convergence speed is greatly improved, so that under the same computing resources, high resolution inversion of complex geological structure is realized. BRIEF DESCRIPTION OF DRAWINGS

[0014] In order to more clearly illustrate the technical solutions in the embodiments of the present application or the prior art, the following will briefly introduce the drawings needed to be used in the embodiments or prior art description. Obviously, the drawings in the following description only constitute some embodiments of the present application, and for those skilled in the art, other drawings can also be obtained from these drawings without creative labor.

[0015] Figure 1 A flowchart of a three-dimensional magnetotelluric inversion method based on L1 norm and adaptive matrix estimation algorithm of an embodiment of the present application is shown in Figure 2 A convergence path diagram of a gradient descent method of an embodiment of the present application is shown in Figure 3 A chessboard model and grid diagram of an embodiment of the present application is shown in Figure 4 A relationship diagram between inversion parameters and iteration numbers of an embodiment of the present application is shown in Figure 5 A horizontal slice diagram of a chessboard model inversion result of an embodiment of the present application is shown in Figure 6 A vertical slice diagram of a chessboard model inversion result of an embodiment of the present application is shown in Figure 7 A data fitting difference distribution diagram of different data fitting metrics and optimization methods in a chessboard model inversion test of an embodiment of the present application is shown in DETAILED DESCRIPTION

[0016] The technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments only constitute some embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor are within the protection scope of the present application.

[0017] A flowchart of a three-dimensional magnetotelluric inversion method based on L1 norm and adaptive matrix estimation algorithm of an embodiment of the present application is shown in Figure 1 as shown, comprising: Collecting magnetotelluric data of a target area as observation data; In the embodiment of the present application, firstly, a plurality of measuring points are arranged in a grid form according to the geological structure of the target area, each measuring point is equipped with a broadband magnetometer, a non-polarized electrode and a GPS clock and other data electromagnetic acquisition equipment, and the data acquisition time is set according to the actual situation, and the geoelectric data of the target area is obtained by continuously acquiring data; the geoelectric data acquired includes amplitude and phase data, and the standard deviation of the geoelectric data is further calculated after the geoelectric data is obtained, so as to facilitate data calculation in the subsequent inversion process.

[0018] The non-structural tetrahedral grid of the target area is divided, and an initial resistivity model is established in the non-structural tetrahedral grid according to the observation data; In the embodiment of the present application, the non-structural tetrahedral grid of the target area is divided based on the terrain data, geological data and the like of the target area, such as DEM elevation model, fault, lithological interface and the like, and a division algorithm is given in the embodiment of the present application, which can be a Delaunay triangulation algorithm, the target area is discretized into non-structural tetrahedral grid units by the algorithm, and a local encryption strategy is adopted for the near-surface complex terrain area and the abnormal body boundary, such terrain area is represented as a region with a grid unit edge length ≤1km in the tetrahedral grid, and the grid unit edge length gradient increases in the deep area, that is, the maximum edge length in the tetrahedral grid is ≤50km, so as to balance the calculation accuracy and efficiency.

[0019] For the initial resistivity model, the observation data acquired by the present embodiment is assigned to each unit of the tetrahedral grid by the Kriging difference method, and the geological prior information of the target area is taken as a hard constraint to limit the upper limit of the resistivity of each unit in the tetrahedral grid, and finally the resistivity data of each unit in the tetrahedral grid is obtained as the initial resistivity model.

[0020] The predicted electromagnetic response data of the target area is obtained by forward calculation according to the current resistivity model; In the embodiment of the present application, the process of obtaining the predicted electromagnetic response data of the target area according to the current resistivity model is as follows: The frequency domain electric field double spinor equation of the target area is established by using Maxwell equation set according to the current resistivity model; the frequency domain electric field of each tetrahedral unit in the non-structural tetrahedral grid is interpolated by using vector basis function; the mass matrix and stiffness matrix of each tetrahedral unit are constructed based on the Galerkin method; the linear equation set of the frequency domain electric field is established according to the mass matrix and stiffness matrix of all tetrahedral units; the linear equation set is solved by using a direct solver to obtain the frequency domain electric field value of the target area; the magnetic field of the target area is calculated by Faraday's law of electromagnetic induction according to the frequency domain electric field value of the target area; and the predicted electromagnetic response data of the target area is obtained according to the frequency domain electric field value and the magnetic field of the target area.

[0021] In one specific embodiment of the present application, the differential form of Maxwell's equations is expressed as: ; wherein, is the electric field intensity, with the unit of V / m, ; is the magnetic induction intensity, with the unit of T, ; is the magnetic field intensity, with the unit of A / m, ; is the current density, with the unit of A / m2, ; is the electric displacement vector, with the unit of C / m2, ; is the free charge density, with the unit of C / m3, ; is the vector differential operator, denotes the curl; assuming the time-harmonic factor of e-jwt, wherein, is the angular frequency, and ignoring the displacement current, the electric field double-curl equation in the frequency domain can be obtained by transformation, expressed as: ; wherein, is the magnetic permeability in vacuum, is the imaginary unit, is the electric conductivity, is the angular frequency.

[0022] In order to use the vector finite element method for the geomagnetic electromagnetic forward, it is necessary to interpolate the electric field of each element in the tetrahedral grid by using the vector basis function, expressed as: ; wherein, denotes the electric field of the element in the tetrahedral grid, denotes the electric field along the i-th boundary of the tetrahedral grid, the superscript e denotes the element in the tetrahedral grid, is the vector basis function, is the three-dimensional space coordinate, by using the vector basis function to interpolate the electric field of each element in the tetrahedral grid, the false interpretation in simulation can be effectively avoided, and the vector basis function can be defined as: ; wherein, denotes the scalar basis function, , denote the starting node and the ending node in the tetrahedral grid, respectively, represents the length of the i-th edge in the tetrahedral mesh, and the superscript e represents a unit in the tetrahedral mesh.

[0023] In the embodiment of the present application, based on the Galerkin method, a vector basis function is used as a weight function, and a weighted residual method is applied to a double spinor equation of an electric field in a frequency domain in a solution space Ω to obtain: In the formula, is a vector differential operator, represents a spinor, represents a solution space, is an electric field intensity, is a magnetic permeability in vacuum, is an imaginary unit, is an electric conductivity, is an angular frequency, is a vector basis function.

[0024] Further according to a vector identity, and substituting an expression of the electric field in each unit of the tetrahedral mesh through the vector basis function interpolation, a finite element equation of the unit e in the tetrahedral mesh is obtained as: In the formula, is a mass matrix, which is defined as , is a stiffness matrix, which is defined as In the formula, represents a volume of the tetrahedral structure, and the indices i, j represent boundary numbers.

[0025] Based on a relationship between local and global mesh units in the tetrahedral mesh, matrices of all units can be assembled, and a Dirichlet boundary condition is applied on an outer boundary. When a boundary distance from a three-dimensional anomaly is far enough (usually 3-5 times a skin depth), an electric field can be approximated by an induced field in a uniform half-space, thereby obtaining a linear equation group of the electric field in the frequency domain, which is expressed as: In the formula, , represents a coefficient matrix, wherein, is a stiffness matrix, is a mass matrix, is a magnetic permeability in vacuum, is an imaginary unit, is an angular frequency, is an electric field intensity of an edge to be solved, is a right end term related to a boundary condition.

[0026] ​​​Since the coefficient matrix of the whole tetrahedral mesh is symmetric and sparse, the compressed storage method is used to reduce the memory consumption in the embodiment of the application, and a direct solver is used to solve the equation, so that the electric field intensity E at any edge in the tetrahedral mesh is obtained, and then the electric field at any position is calculated through the frequency domain electric field double vorticity equation, and then the corresponding magnetic field intensity is calculated through the Faraday electromagnetic induction law At this time, the impedance of the magnetotelluric can be calculated through the following expression: ; In the formula, Z represents the impedance of the magnetotelluric, E represents the electric field intensity, H represents the magnetic field intensity, and the superscript -1 represents the inverse of the matrix, and the superscripts 1 and 2 of the electric field and the magnetic field represent different polarization modes, and represent the direction; the impedance of the magnetotelluric obtained by calculation can obtain the predicted electromagnetic response data.

[0027] An inversion objective function is constructed according to a data fitting term based on an L1 norm and a model constraint term based on an L2 norm; Since the three-dimensional inversion problem of the magnetotelluric is a typical underdetermined problem, the regularization method is used to construct the objective function in the embodiment of the application, in the embodiment of the application, the data fitting term based on the L1 norm is established according to the observation data and the predicted electromagnetic response data; the model constraint term based on the L2 norm is established according to the initial resistivity model and the current resistivity model, and the inversion objective function is constructed through the data fitting term based on the L1 norm and the model constraint term based on the L2 norm, and is represented as: ; In the formula, Z represents the inversion objective function, C represents the covariance matrix of the observation data, C represents the covariance matrix of the predicted electromagnetic response data, d represents the observation data, d represents the predicted electromagnetic response data, that is, the data fitting term, r represents the initial resistivity model, r represents the current resistivity model, that is, the model constraint term, is a weighting parameter, is the L1 norm, is the L2 norm.

[0028] calculating a data fitting difference between the observation data and the predicted electromagnetic response data under the current resistivity model; if the data fitting difference does not meet the set condition, updating the current resistivity model by using the ADAM algorithm; calculating a data fitting difference between the observation data and the predicted electromagnetic response data corresponding to the updated resistivity model, and iteratively calculating in sequence until the data fitting difference meets the set condition, and outputting the resistivity model at this time as the inversion result.

[0029] The data fitting term will generally converge from a value much larger than 1.0 to a set threshold value, and if the L2 norm is used to measure the data fitting term, the objective function will be highly sensitive to outliers, which are usually caused by serious observation errors. Ideally, outliers in the data can be assumed, and when the objective function is constructed, it is selected to remove or reduce the weight of the outliers. However, in practice, it is difficult to assign a small weight to each data with a large fitting difference, because many of them are of high quality. At this time, a large amount of numerical tests are required to identify the real outliers. Sometimes, a large data fitting difference may be caused by a large geological anomaly, in which case the weight of the data fitting term should not be changed to avoid inversion bias; another method is to use the L1 norm to measure the data fitting term, which is not overly sensitive to large outliers, thereby obtaining a more robust inversion process. Since there is no effective tool to optimize the non-smooth objective function, the L1 norm is mainly used to form the lasso problem in geophysics, and its application as a data fitting term is less. The emergence of the ADAM method has changed this situation, and therefore, the resistivity model for inversion in the objective function is optimized by using this method in subsequent research.

[0030] In the embodiment of the application, the data fitting difference is defined as: ; Wherein, is the data fitting difference, is the covariance matrix of the observation data, is the observation data, is the predicted electromagnetic response data, is the L1 norm, is the total data amount.

[0031] In the embodiment of the application, it is set that when the data fitting difference is less than the set threshold value or does not change after multiple iterations, the iteration is stopped and the final model is output. The embodiment of the application gives a set threshold value which can be .

[0032] In a specific embodiment of the present invention, the ADAM algorithm is an adaptive matrix estimation algorithm. Its core idea is to use the first and second momentum of the gradient to estimate the gradient, thereby more accurately adjusting the update step size of each parameter. The first momentum represents the exponentially weighted average of the gradient and is used to estimate the mean of the gradient, smoothing the gradient update process and reducing oscillations. Figure 2 As shown, the blue semi-circular dashed line represents the quadratic function approximation path of the finite-memory quasi-Newton algorithm, while the green convergence path exhibits less oscillation. This is attributed to the effect of the first moment estimated through gradient moving average. Furthermore, when the algorithm approaches the minimum, the second moment automatically adjusts the step size, thereby reducing the probability of skipping the minimum. The gradient of the objective function under the current resistivity model is calculated, and the process of updating the inverted resistivity model in the objective function using this algorithm is as follows: Calculate the first-order and second-order momentum based on the gradient of the objective function; update the current resistivity model based on the first-order and second-order momentum to obtain the updated resistivity model.

[0033] In this embodiment of the invention, the recursive formula for calculating the first-order momentum based on the gradient of the objective function is as follows: ; In the formula, Let represent the first-order momentum at step i. Indicates the first The first-order momentum of the step, The gradient of the objective function. The decay coefficient of the first-order momentum is usually taken as 0.9. To reduce the problem of the initial momentum being too small, the first-order momentum calculated in the ADAM algorithm is further corrected for deviation in this embodiment of the invention, as follows: ; in, This represents the corrected first-order momentum. Let represent the first-order momentum at step i. The decay coefficient of first-order momentum, This indicates the number of update steps for the current resistivity model.

[0034] For second-order momentum, the ADAM algorithm introduces an exponentially weighted average of the squared gradient as the second-order momentum to estimate the variance of the gradient, and its expression is: ; In the formula, Let represent the second momentum at step i. No. The second momentum of the step, For Hadama accumulation, is the decay coefficient of the second order momentum, usually taking 0.999, used to smooth the variance of the gradient, so as to dynamically adjust the learning rate, prevent the learning rate from being too large or too small due to the sharp change of the gradient; similarly, the second order momentum also carries out deviation correction, and the correction formula is: ; wherein, denotes the corrected second order momentum, denotes the second order momentum of the i-th step, is the decay coefficient of the second order momentum.

[0035] According to the above, the current resistivity model is updated according to the first order momentum and the second order momentum, and is denoted as: ; wherein, denotes the current resistivity model, denotes the updated resistivity model, denotes the corrected first order momentum, denotes the corrected second order momentum, is the global learning rate, is a constant to prevent the denominator from being zero.

[0036] In an experimental embodiment of the present application: As shown in Figure 3 , the present embodiment designs a three-dimensional magnetotelluric chessboard model for inversion test, which contains nine abnormal bodies, the top of which is buried at a depth of 25 kilometers, and the thickness is 100 kilometers, and the specific model can be seen from Figure 3 part (a), the size of all abnormal bodies is 100 kilometers x 100 kilometers x 50 kilometers, the interval is 20 kilometers, the background resistivity is 100 W x m, the resistivity of the abnormal bodies is 1000 W x m (blue) and 10 W x m (red) respectively, in x and y directions, there are a total of 289 measuring points uniformly distributed on the ground surface, the interval is 25 kilometers, the center area of the model is 500 kilometers x 500 kilometers x 150 kilometers, the size of the entire model including the extension is 10000 kilometers x 10000 kilometers x 10000 kilometers, in order to obtain accurate predicted electromagnetic response, the measuring points of the forward model also need to be locally encrypted, and finally 697,543 unstructured tetrahedral meshes are obtained, that is, Figure 3 part (b); and for the inversion model, the grid of the abnormal body is removed, and the measuring points are still locally encrypted, and finally 594077 unstructured tetrahedral meshes are generated, as shown in Figure 3Six frequencies were selected to generate full impedance data in the logarithmic equidistant range from 0.0001 to 0.1 Hz, and 3% Gaussian random noise was added as the input data of inversion, the inversion test was calculated on an Intel® Xeon® Gold 6258R CPU@2.7GHZ, 512GB memory workstation, the calculation time of each frequency was 30 seconds, and each iteration took about 6 minutes and 40 seconds, in this test, the initial resistivity model was a uniform half-space model with 100 W x m, lambda was set to 10, and the cooling factor was 0.5.

[0037] According to the above settings, the present embodiment first analyzes the convergence of the BFGS method and the ADAM method for the objective function in the form of L2-L2 and L1-L2, as shown in Figure 4 Regardless of which method and which form of objective function is selected, the inversion can converge to the predetermined threshold of data fitting, but the convergence speed of the ADAM method is significantly faster than that of the BFGS method, as shown in Figure 4 As shown in part (b) of FIG. 1, for the data fitting term of the L1 norm, the convergence curve of the ADAM method is a rapidly descending curve on a logarithmic scale in the early iterations, in contrast, the BFGS method has a slower descending speed in the initial iterations, and this phenomenon occurs because the objective function is mainly affected by the data fitting term in the early iterations, and the data fitting term of the norm is a non-smooth function without a second-order derivative, which makes it challenging for the BFGS method to accurately approximate the Hessian matrix, Figure 4 As shown in parts (a) and (b) of FIG. 1, the ADAM method has obvious oscillation in the convergence curve when it approaches convergence, which is caused by inaccurate step size estimation, although this phenomenon requires more iterations for the ADAM method, it sometimes helps to escape local minima and obtain good inversion results. Figure 4 Part (c) of FIG. 1 gives the variation curve of the model constraint term with the number of iterations for different results, and the final model constraint term variation tends to be consistent; in addition, Figure 4 Part (d) of FIG. 1 gives the variation curve of the regularization factor l with the number of iterations, it can be seen that when the L1 norm is used, l decreases lower, especially when the ADAM algorithm is used, the regularization factor l is constantly corrected to update the direct balance relationship between the model constraint term and the data fitting term.

[0038] Figure 5 and Figure 6 Further show the depth slices and vertical slices of the inversion results obtained by using different methods and different data fitting measures, respectively, wherein, Figure 5 Part (a) of FIG. 2 is the ADAM-L2 inversion result, Figure 5Part (b) in the table shows the BFGS-L2 inversion results. Figure 5 Part (c) in the table represents the ADAM-L1 inversion results. Figure 5 Part (d) in the table represents the BFGS-L1 inversion results. Figure 5 Part (e) in the model is the real model; Figure 6 In the image, the first row shows vertical slices of the real model, the second row shows vertical slices of the ADAM-L2 inversion result, the third row shows vertical slices of the ADAM-L1 inversion result, the fourth row shows vertical slices of the BFGS-L2 inversion result, and the fifth row shows vertical slices of the BFGS-L1 inversion result. It can be observed that for both data fit metrics, the BFGS method shows little difference in its ability to characterize the lower boundary of the outlier model, and both have relatively low resolution. The main reason for this is that the L2 norm tends to minimize the data fit terms of all measurement points to a similarity level, combined with the influence of L2 norm model constraint terms that may lead to smooth boundaries in the inversion result. When using the L1 norm to measure the data fit terms... Figure 4 The results show that the BFGS method struggles to fit data in the initial stages. Only after dozens of inversion iterations can the BFGS method achieve stable convergence, and the overall form of the objective function may approximate as a smooth function with a second derivative. During this process, the inversion objective function gradually transforms into a quadratic function, primarily reducing the poor fit on large datasets. Figure 7 As shown, the inversion results are similar to those obtained using the L2-L2 form objective function. In contrast, the ADAM method converges normally in both L1 and L2 norm cases and better reveals the differences in inversion results between the two data fit difference measures.

[0039] In addition, it can be found Figure 5 and Figure 6 The ADAM-L1 inversion results show the highest resolution, particularly in the recovery of the morphology and resistivity of high-resistivity anomalies. When using the L2 norm as a measure of data fitting, data with large fitting errors significantly impact the inversion process, which primarily focuses on reducing these errors. Although high-resistivity anomalies also contribute to data fitting errors, they become less significant in the inversion compared to low-resistivity anomalies, which generate larger errors. Modifying low-resistivity anomalies suffices for data fitting. In contrast, the L1 norm-based data fitting better balances the fitting errors across all measurement points and frequencies, resulting in more accurate recovery of high-resistivity anomalies in the inversion results.

[0040] To better analyze the differences between L1 norm and L2 norm in the data fitting term Figure 7The data fitting differences of each measuring point of the data fitting term of different methods and different norms are given. The errors of the data fitting are similar no matter which data fitting measurement method is selected, which indicates that the data fitting converges to a similar level. The data fitting term of L1 norm has more abnormal values than the data fitting term of L2 norm. This is because the L1 norm does not overfit the points with larger data fitting differences, can better maintain the original inversion weight of the data, and will not excessively amplify the influence of the points with larger fitting errors. If the data fitting term of L2 norm is used to modify the inversion weight based on the original data fitting error, the influence of the points with smaller fitting errors will be weakened, or the weak anomaly will be covered by the strong anomaly.

[0041] The above merely describes preferred embodiments of the present application and is not intended to limit the present application. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present application shall be included in the protection scope of the present application.

Claims

1. A three-dimensional magnetotelluric inversion method based on L1 norm and adaptive matrix estimation algorithm, characterized in that, The method comprises the following steps: Collecting magnetotelluric data of a target area as observation data; Dividing a non-structural tetrahedral grid of the target area, and establishing an initial resistivity model in the non-structural tetrahedral grid according to the observation data; Performing forward calculation according to the current resistivity model to obtain predicted electromagnetic response data of the target area; Establishing a data fitting term based on L1 norm according to the observation data and the predicted electromagnetic response data; Establishing a model constraint term based on L2 norm according to the initial resistivity model and the current resistivity model; Constructing an inversion objective function according to the data fitting term based on L1 norm and the model constraint term based on L2 norm; Calculating a data fitting difference between the observation data and the predicted electromagnetic response data under the current resistivity model; If the data fitting difference does not meet the set condition, updating the current resistivity model by using an ADAM algorithm; Calculating a data fitting difference between the observation data and predicted electromagnetic response data corresponding to the updated resistivity model, and iteratively calculating in sequence until the data fitting difference meets the set condition, and outputting the resistivity model at this time as an inversion result.

2. The L1 norm and adaptive matrix estimation algorithm based 3D magnetotelluric inversion method of claim 1, wherein: The method further comprises the following steps of performing forward calculation according to the current resistivity model to obtain predicted electromagnetic response data of the target area: Establishing a frequency domain electric field double spinor equation of the target area by using Maxwell equations according to the current resistivity model; Interpolating a frequency domain electric field of each tetrahedral element in the non-structural tetrahedral grid by using a vector basis function; Constructing a mass matrix and a stiffness matrix of each tetrahedral element based on a Galerkin method; Establishing a linear equation group of the frequency domain electric field according to the mass matrix and the stiffness matrix of all tetrahedral elements; Solving the linear equation group by using a direct solver to obtain frequency domain electric field values of the target area; Calculating a magnetic field of the target area by Faraday's law of electromagnetic induction according to the frequency domain electric field values of the target area; Obtaining predicted electromagnetic response data of the target area according to the frequency domain electric field values and the magnetic field of the target area.

3. The L1 norm and adaptive matrix estimation algorithm based 3D MT inversion method of claim 1, wherein: The method further comprises the following steps of constructing an inversion objective function according to the data fitting term based on L1 norm and the model constraint term based on L2 norm, and represented as: ； wherein, denotes an inversion objective function, is a covariance matrix of the observed data, is a covariance matrix of the predicted electromagnetic response data, is the observed data, is the predicted electromagnetic response data, is an initial resistivity model, is a current resistivity model, is a trade-off parameter, is an L1 norm, is an L2 norm.

4. The L1 norm and adaptive matrix estimation algorithm based 3D MT inversion method of claim 1, wherein: The method further comprises the following steps of updating the current resistivity model by using an ADAM algorithm, comprising: Calculating a gradient of the target function under the current resistivity model; Calculating a first momentum and a second momentum according to the gradient of the target function; Updating the current resistivity model according to the first momentum and the second momentum to obtain an updated resistivity model.

5. The L1 norm and adaptive matrix estimation algorithm based 3D MT inversion method of claim 4, wherein: The method further comprises the following steps of calculating a first momentum and a second momentum according to the gradient of the target function, and represented as: The method further comprises the following steps of updating the current resistivity model according to the first momentum and the second momentum, and represented as: ； wherein, denotes the first order momentum after correction, denotes the first order momentum of the i-th step, is the decay factor of the first order momentum, denotes the number of updates of the current resistivity model; ； wherein, denotes the modified second order momentum, denotes the second order momentum of the i-th step, is the decay factor of the second order momentum.

6. The L1 norm and adaptive matrix estimation algorithm based 3D MT inversion method of claim 4, wherein: The method further comprises the following steps of calculating a data fitting difference between the observation data and the predicted electromagnetic response data under the current resistivity model, and represented as: ； wherein, denotes the current resistivity model, denotes the updated resistivity model, denotes the corrected first order momentum, denotes the corrected second order momentum, is a global learning rate, is a constant to prevent division by zero.

7. The L1 norm and adaptive moment estimation algorithm based 3D MT inversion method of claim 1, wherein: The set condition is that the iteration is stopped when the data fitting difference is less than a set threshold or does not change after multiple iterations. ； wherein, is the difference for data fitting, is the covariance matrix of the observed data, is the observed data, is the predicted electromagnetic response data, , is the total amount of data.

8. The L1 norm and adaptive moment estimation algorithm based 3D MT inversion method of claim 1, wherein: ​

Images