Three-dimensional magnetotelluric inversion method based on l1 norm and adaptive matrix estimation algorithm

Abstract

The present application relates to the technical field of geoelectromagnetism, in particular to a three-dimensional magnetotelluric inversion method based on L1 norm and adaptive moment estimation algorithm. Collect magnetotelluric data of a target region as observation data; divide the non-structured tetrahedral grid of the target region to establish an initial resistivity model; perform forward calculation according to the current resistivity model to obtain predicted electromagnetic response data of the target region; construct an inversion objective function according to a data fitting term based on L1 norm and a model constraint term based on L2 norm; calculate the data fitting difference between the observation data and the predicted electromagnetic response data; if the data fitting difference does not meet the set condition, update the current resistivity model by using ADAM algorithm; iteratively, until the data fitting difference meets the set condition, and output the resistivity model at this time as the inversion result. The present application can effectively suppress the influence of large data fitting error in inversion, and obtain higher inversion resolution.

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): 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): 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 structural expression ability of the model is limited, and traditional Tikhonov or least squares smoothing method is easy to cause abnormal boundary blur, although there are minimum support and sparse regularization improvements, but it is still difficult to balance the 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. The method collects magnetotelluric data of a target region as observation data; a non-structured tetrahedral grid of the target region is divided to establish an initial resistivity model; forward calculation is performed according to the current resistivity model to obtain predicted electromagnetic response data of the target region; an inversion objective function is constructed according to a data fitting term based on L1 norm and a model constraint term based on L2 norm; a data fitting difference between the observation data and the predicted electromagnetic response data 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 iteration is sequentially performed until the data fitting difference meets the set condition, and the resistivity model at this time is output as the inversion result. The application can effectively suppress the influence of large data fitting error in inversion and obtain higher inversion resolution.

[0005] The application adopts the following technical scheme, a three-dimensional magnetotelluric inversion method based on L1 norm and adaptive moment estimation algorithm, comprising:

[0006] Collect magnetotelluric data of a target region as observation data;

[0007] A non-structured tetrahedral grid of the target region is divided, and an initial resistivity model is established in the non-structured tetrahedral grid according to the observation data;

[0008] Forward calculation is performed according to the current resistivity model to obtain predicted electromagnetic response data of the target region;

[0009] A data fitting term based on L1 norm is established according to the observation data and the predicted electromagnetic response data;

[0010] A model constraint term based on L2 norm is established according to the initial resistivity model and the current resistivity model;

[0011] 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;

[0012] A data fitting difference between the observation data and predicted electromagnetic response data under the current resistivity model is calculated;

[0013] If the data fitting difference does not meet the set condition, the ADAM algorithm is used to update the current resistivity model;

[0014] A data fitting difference between the observation data and predicted electromagnetic response data corresponding to the updated resistivity model is calculated, and the iteration is sequentially performed until the data fitting difference meets the set condition, and the resistivity model at this time is output as the inversion result.

[0015] Further, forward calculation is performed according to the current resistivity model to obtain predicted electromagnetic response data of the target region, specifically:

[0016] The frequency-domain electric field double-curl equations of the target region are established according to the current resistivity model and Maxwell equations;

[0017] The frequency-domain electric field of each tetrahedral element in the unstructured tetrahedral mesh is interpolated by using a vector basis function;

[0018] The mass matrix and the stiffness matrix of each tetrahedral element are constructed based on the Galerkin method;

[0019] The linear equations of the frequency-domain electric field are established according to the mass matrix and the stiffness matrix of all tetrahedral elements;

[0020] The linear equations are solved by using a direct solver to obtain the frequency-domain electric field value of the target region;

[0021] 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;

[0022] 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.

[0023] 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:

[0024] ;

[0025] In the formula, the inversion objective function, is the covariance matrix of the observation data, is the covariance matrix of the predicted electromagnetic response data, is the observation data, is the predicted electromagnetic response data, is the initial resistivity model, is the current resistivity model, is the trade-off parameter, is the L1 norm, is the L2 norm.

[0026] Further, the ADAM algorithm is used to update the current resistivity model, including:

[0027] The gradient of the objective function under the current resistivity model is calculated;

[0028] The first-order momentum and the second-order momentum are calculated according to the gradient of the objective function;

[0029] The current resistivity model is updated according to the first-order momentum and the second-order momentum to obtain an updated resistivity model.

[0030] Further, when calculating the first-order momentum and the second-order momentum according to the gradient of the objective function, the following is further included:

[0031] The first-order momentum and the second-order momentum are bias-corrected, and are respectively represented as:

[0032] ;

[0033] wherein, represents the corrected first-order momentum, represents the first-order momentum at the i-th step, is an attenuation coefficient of the first-order momentum, represents the number of updates of the current resistivity model;

[0034] ;

[0035] wherein, represents the corrected second-order momentum, represents the second-order momentum at the i-th step, is an attenuation coefficient of the second-order momentum.

[0036] Further, the current resistivity model is updated according to the first-order momentum and the second-order momentum, and is represented as:

[0037] ;

[0038] wherein, represents the current resistivity model, represents the updated resistivity model, represents the corrected first-order momentum, represents the corrected second-order momentum, is a global learning rate, is a constant to prevent the denominator from being zero.

[0039] Further, a data fitting difference between the observed data and the predicted electromagnetic response data under the current resistivity model is calculated, and is represented as:

[0040] ;

[0041] wherein, is the data fitting difference, is a covariance matrix of the observed data, is the observed data, is the predicted electromagnetic response data, is an L1 norm, and N is the total amount of data.

[0042] 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.

[0043] The beneficial effects of the present application are: the present application adopts a non-structural grid vector finite element forward, which can significantly reduce the memory consumption, solve the storage bottleneck of Jacobian matrix such as Gauss-Newton method, further combine L1 norm and L2 norm to construct the inversion objective function, establish the data fitting term by using L1 norm, which can effectively suppress the 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 using ADAM algorithm combined with first order and second order momentum to update the inversion model, compared with the 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

[0044] In order to more clearly illustrate the technical solutions of the embodiments of the present application or the prior art, the drawings needed to be used in the embodiments or prior art description will be briefly introduced below. Obviously, the drawings in the following description are only some embodiments of the present application, and other drawings can be obtained by those skilled in the art without creative labor.

[0045] Figure 1 The flow chart of the three-dimensional magnetotelluric inversion method based on L1 norm and adaptive matrix estimation algorithm of the embodiment of the present application is shown.

[0046] Figure 2 The convergence path diagram of the gradient descent method of the embodiment of the present application is shown.

[0047] Figure 3 The chessboard model and grid diagram of the embodiment of the present application is shown.

[0048] Figure 4 The relationship diagram between the inversion parameters and the number of iterations of the embodiment of the present application is shown.

[0049] Figure 5 The horizontal slice diagram of the chessboard model inversion result of the embodiment of the present application is shown.

[0050] Figure 6 The vertical slice diagram of the chessboard model inversion result of the embodiment of the present application is shown.

[0051] Figure 7 The data fitting difference distribution diagram of different data fitting metrics and optimization methods in the chessboard model inversion test of the embodiment of the present application is shown. DETAILED DESCRIPTION

[0052] With reference to the accompanying drawings, the technical solutions in the embodiments of the present application will be clearly and completely described below, obviously, the described embodiments are only a part of the 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 fall within the scope of the present application.

[0053] The three-dimensional magnetotelluric inversion method based on the L1 norm and adaptive matrix estimation algorithm in the embodiments of the present application is shown in the flowchart as Figure 1 shown, comprising:

[0054] Collecting magnetotelluric data of a target area as observation data;

[0055] In the embodiments of the present application, first, a plurality of measuring points are arranged in a grid form according to the geological structure of the target area, and each measuring point is simultaneously provided with data electromagnetic acquisition equipment such as a broadband magnetometer, a non-polarized electrode and a GPS clock, and the data collection time is set according to the actual situation, and the magnetotelluric data of the target area is obtained by continuously collecting data; the collected magnetotelluric data includes amplitude and phase data, and the standard deviation of the magnetotelluric data is further calculated after the magnetotelluric data is obtained, so as to facilitate data calculation in the subsequent inversion process.

[0056] 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;

[0057] In the embodiments 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, lithology interface and the like, and a division algorithm given in the embodiments of the present application can be a Delaunay triangulation algorithm, which discretizes the target area into non-structural tetrahedral grid units, and adopts a local encryption strategy 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.

[0058] For the initial resistivity model, the observation data collected by the present application 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.

[0059] According to the current resistivity model, the predicted electromagnetic response data of the target area is obtained by forward calculation;

[0060] In the embodiment of the present application, the process of obtaining the predicted electromagnetic response data of the target region according to the current resistivity model is specifically as follows:

[0061] The frequency domain electric field double curl equation of the target region is established according to the current resistivity model by using Maxwell equations; the frequency domain electric field of each tetrahedral element in the unstructured tetrahedral mesh is interpolated by 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 equation set of the frequency domain electric field is established according to the mass matrix and the stiffness matrix of all tetrahedral elements; the linear equation set is 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 by Faraday's law of electromagnetic induction; and 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.

[0062] In one specific embodiment of the present application, the differential form of Maxwell equations is represented as:

[0063] ;

[0064] In the formula, is the electric field intensity, and the unit is V / m ; is the magnetic induction intensity, and the unit is T ; is the magnetic field intensity, and the unit is A / m ; is the current density, and the unit is A / m2 ; is the electric displacement vector, and the unit is C / m2 ; is the free charge density, and the unit is C / m3 ; is the vector differential operator, represents the curl; assuming the time-harmonic factor is , wherein, is the angular frequency, and the displacement current is ignored, and the electric field double curl equation in the frequency domain is obtained by transformation combined with the constitutive relation, and is represented as:

[0065] ;

[0066] In the formula, is the magnetic permeability in vacuum, is the imaginary unit, is the electrical conductivity, is the angular frequency.

[0067] In order to use the vector finite element method for the magnetotelluric forward modeling, the electric field of each element in the tetrahedral mesh is interpolated by using a vector basis function, and is represented as:

[0068] ;

[0069] wherein, represents the electric field of a unit e in a tetrahedral mesh, represents the electric field along the i-th boundary in a tetrahedral mesh, and the superscript e represents a unit in the tetrahedral mesh, is a vector base function, is a three-dimensional space coordinate, by using the vector base function to interpolate the electric field of each unit in the tetrahedral mesh, false interpretation in simulation can be effectively avoided, and the vector base function may be defined as:

[0070] ;

[0071] wherein, represents a scalar base function, , respectively represent a starting node and an ending node in the tetrahedral mesh, represents the length of the i-th edge in the tetrahedral mesh, and the superscript e represents a unit in the tetrahedral mesh.

[0072] In the embodiment of the application, based on the Galerkin method, the vector base function is used as a weight function, the weighted residual method is applied to the electric field double vorticity equation in the frequency domain in a solution space Omega, and the following formula is obtained:

[0073] ;

[0074] wherein, is a vector differential operator, represents a vorticity, 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 base function.

[0075] Further according to the vector identity, and by substituting the expression of the electric field of each unit in the tetrahedral mesh interpolated by the vector base function, the finite element equation of the unit e in the tetrahedral mesh is obtained as follows:

[0076] ;

[0077] wherein, is a mass matrix, which is defined as , is a stiffness matrix, which is defined as ​wherein, represents the volume of the tetrahedron structure, and indexes i,j represent the boundary number.

[0078] Based on the relationship between the local and global grid cells in the tetrahedral mesh, the matrix of all cells can be assembled, and the Dirichlet boundary condition is applied on the outer boundary. When the boundary is far enough from the 3D anomaly body (usually 3-5 times the skin depth), the electric field can be approximated by the induced field in a homogeneous half-space, thereby obtaining the linear equation set of the frequency domain electric field, expressed as:

[0079] ;

[0080] In the formula, , represents the coefficient matrix, wherein, is the stiffness matrix, is the mass matrix, is the magnetic permeability in vacuum, is the imaginary unit, is the angular frequency, is the electric field intensity of the edge to be solved, is the right end term related to the boundary condition.

[0081] Since the coefficient matrix of the tetrahedral mesh as a whole is symmetric and sparse, the compressed storage method is adopted in the embodiment of the present application to reduce memory consumption, and the direct solver is used to solve the equation. After obtaining the electric field intensity E at any edge in the tetrahedral mesh, the electric field at any position can be calculated through the frequency domain electric field double spinor equation, and then the corresponding magnetic field intensity can be calculated through the Faraday's law of electromagnetic induction At this time, the impedance of the magnetotelluric can be calculated by the following expression:

[0082] ;

[0083] In the formula, represents the impedance of the magnetotelluric, is the electric field intensity, is the magnetic field intensity, and the superscript 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; by calculating the impedance of the magnetotelluric, the predicted electromagnetic response data can be obtained.

[0084] According to the data fitting term based on the L1 norm and the model constraint term based on the L2 norm, an inversion objective function is constructed;

[0085] Since the magnetotelluric three-dimensional inversion problem is a typical underdetermined problem, the regularization method is used to construct the objective function in the embodiment of the present application, in the embodiment of the present application, the data fitting term based on L1 norm is established according to the observed data and the predicted electromagnetic response data; the model constraint term based on L2 norm is established according to the initial resistivity model and the current resistivity model, and the inversion objective function is constructed by the data fitting term based on L1 norm and the model constraint term based on L2 norm, and is expressed as:

[0086] ;

[0087] In the formula, indicates the inversion objective function, is the covariance matrix of the observed data, is the covariance matrix of the predicted electromagnetic response data, is the observed data, is the predicted electromagnetic response data, that is, the data fitting term, is the initial resistivity model, is the current resistivity model, that is, the model constraint term, is the weighting parameter, is the L1 norm, is the L2 norm.

[0088] The data fitting difference between the observed data and the predicted electromagnetic response data under the current resistivity model is calculated; if the data fitting difference does not satisfy the set condition, the ADAM algorithm is used to update the current resistivity model; the data fitting difference between the observed data and the predicted electromagnetic response data corresponding to the updated resistivity model is calculated, and the iteration is performed in turn until the data fitting difference satisfies the set condition, and the resistivity model at this time is output as the inversion result.

[0089] The data fitting term usually converges from a value much larger than 1.0 to a set threshold, 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 selected to be removed or reduced in weight when constructing the objective function, however, it is actually impossible 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 needed to identify the real outliers, sometimes, a large data fitting difference may be caused by a large geological anomaly, in this 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, thus 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, and the emergence of the ADAM method changes this situation, therefore, the embodiment of the present application uses this method to optimize the resistivity model in the inversion of the objective function in subsequent research.

[0090] In the embodiment of the present application, the definition of the data fitting difference is:

[0091] ;

[0092] 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, is the total amount of data.

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

[0094] In one specific embodiment of the present application, the ADAM algorithm is an adaptive matrix estimation algorithm, the core idea of which is to estimate the gradient by using the first order momentum and the second order momentum of the gradient, so as to more accurately adjust the update step of each parameter, the first order momentum represents the exponentially weighted average of the gradient, which is used to estimate the mean value of the gradient to smooth the gradient update process and reduce the oscillation phenomenon, such as Figure 2As 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:

[0095] The first-order and second-order momentum are calculated based on the gradient of the objective function; the current resistivity model is updated based on the first-order and second-order momentum to obtain the updated resistivity model.

[0096] 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:

[0097] ;

[0098] 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:

[0099] ;

[0100] 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.

[0101] 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, expressed as:

[0102] ;

[0103] 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:

[0104] ;

[0105] 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.

[0106] 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:

[0107] ;

[0108] 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.

[0109] In one experimental embodiment of the present application:

[0110] As Figure 3 shown, 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 the part (a) of Figure 3 , 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 the x and y directions, a total of 289 measuring points are 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, the part (b) of Figure 3 ; 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 Wx m, lambda was set to 10, and the cooling factor was 0.5.

[0111] According to the above settings, the present embodiment first analyzes the convergence of the BFGS method and the ADAM method for the objective functions in the forms of L2-L2 and L1-L2, as shown in Figure 4 No matter 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 obviously 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 curve that rapidly decreases on a logarithmic scale in early iterations, and in comparison, the BFGS method has a slower decrease in the initial iterations, and this phenomenon occurs because 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 may require more iterations for the ADAM method, it can sometimes help 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, and 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.

[0112] 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.

[0113] 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.

[0114] 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.

[0115] The above merely describes the preferred embodiments of the present application and is not used 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; Specifically: According to the current resistivity model, Maxwell's equations are used to establish a frequency domain electric field double spinor equation of the target area; A vector basis function is used to interpolate the frequency domain electric field of each tetrahedral element in the non-structural tetrahedral grid; Based on the Galerkin method, a mass matrix and a stiffness matrix of each tetrahedral element are constructed; According to the mass matrix and the stiffness matrix of all tetrahedral elements, a linear equation set of the frequency domain electric field is established; The linear equation set is solved by using a direct solver to obtain the frequency domain electric field value of the target area; According to the frequency domain electric field value of the target area, the magnetic field of the target area is calculated through the Faraday's law of electromagnetic induction; According to the frequency domain electric field value and the magnetic field of the target area, the predicted electromagnetic response data of the target area are obtained; 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; A 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 the 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.

2. The L1 norm and adaptive matrix estimation algorithm based 3D magnetotelluric inversion method of claim 1, wherein: The inversion objective function is constructed according to the data fitting term based on L1 norm and the model constraint term based on L2 norm, and is expressed 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.

3. The L1 norm and adaptive matrix estimation algorithm based 3D MT inversion method of claim 1, wherein: The ADAM algorithm is used to update the current resistivity model, which comprises: 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 the updated resistivity model.

4. The L1 norm and adaptive matrix estimation algorithm based 3D MT inversion method of claim 3, wherein: When the first-order momentum and the second-order momentum are calculated according to the gradient of the objective function, it further comprises: The first-order momentum and the second-order momentum are bias corrected, and are respectively expressed 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.

5. The L1 norm and adaptive matrix estimation algorithm based 3D MT inversion method of claim 3, wherein: The current resistivity model is updated according to the first-order momentum and the second-order momentum, and is expressed 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 the denominator from being zero.

6. The L1 norm and adaptive moment estimation algorithm based 3D MT inversion method of claim 1, wherein: The data fitting difference between the observation data and the predicted electromagnetic response data under the current resistivity model is calculated, and is expressed as: ； wherein, is the data misfit, is the covariance matrix of the observed data, is the observed data, is the predicted electromagnetic response data, is the L1 norm, is the total amount of data.

7. The L1 norm and adaptive moment estimation algorithm based 3D MT inversion method of claim 1, wherein: The set 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.

Images