Hybrid Lp regularization magnetotelluric two-dimensional inversion method based on adaptive weight

Abstract

The invention discloses a hybrid Lp regularization magnetotelluric two-dimensional inversion method based on adaptive weight, and the method comprises the steps: reading magnetotelluric observation data, carrying out the quadrilateral finite element grid discretization of a whole inversion region, building an inversion grid, and building an initial model and a prior model of the inversion conductivity of the inversion region; based on the initial model and the prior model of the inversion conductivity and the magnetotelluric observation data, constructing a mixed Lp regularization inversion objective function comprising an L1 regularization item, an L2 regularization item and a data fitting item; in an inversion iteration process, adaptively adjusting weight factors of the L2 regularization item and the L1 regularization item; and based on the inversion grid, performing magnetotelluric two-dimensional inversion according to the inversion objective function mixed with Lp regularization, the initial model and magnetotelluric observation data, and outputting optimal model parameters of Gaussian Newton inversion to obtain underground electrical structure information. According to the invention, the resolution and accuracy of the inversion result are improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of geophysical electromagnetic method inversion, and particularly relates to a mixed Lp regularization two-dimensional magnetotelluric inversion method. BACKGROUND

[0002] Magnetotelluric method reconstructs the underground electrical structure by observing the electromagnetic field excited by natural field source, has the advantages of low cost and large detection depth, and is an important means for resource exploration and deep structure research. With the gradual depletion of shallow mineral resources, China is shifting the focus of resource exploration to deep land and deep sea, and the exploration area is more complex, so it is urgent to develop high-resolution inversion methods. Two-dimensional magnetotelluric inversion plays an important role in actual exploration as an important means of data interpretation. However, the magnetotelluric inversion problem itself has strong nonlinearity, multiple solutions and ill-posedness, and traditional methods often face problems such as insufficient resolution, blurred boundaries and noise sensitivity under complex geological conditions.

[0003] The invention patent with application number 202310543810.0 discloses a magnetotelluric inversion method based on stratum adaptive encryption. The magnetotelluric inversion method based on stratum adaptive encryption adjusts the inversion grid adjustment indicator factor of the objective function data fitting difference term, and is based on different grid adjustment indicator factors and a variety of inversion grid adaptive encryption strategies. The dependence of the inversion result on the inversion grid is reduced, the convergence speed of the inversion is improved, the number of unknowns in the inversion process is reduced, and the data can be better fitted. However, the above-mentioned inversion does not consider the applicability of different regularizations to different geological models, and only considers the traditional L2 regularization, which realizes a stable solution by minimizing the second-order derivative of the model parameter. Although this method can suppress noise interference, excessive smoothing leads to the inability to effectively depict sharp electrical boundaries such as faults and dikes. In order to improve the resolution of the inversion result, some sparse inversions (such as L1 regularization) are applied to the magnetotelluric inversion problem. Sparse regularization can restore clear geoelectric interfaces, but it is more susceptible to data noise, is prone to produce false abnormal bodies, and has poor numerical calculation stability. In addition, there may be smooth structures and abrupt interfaces in the real geological structure, and a single regularization constraint cannot meet the reconstruction requirements of different structures. SUMMARY

[0004] In view of the low resolution of the existing magnetotelluric inversion method, which cannot realize accurate depiction of the electrical structure in the exploration of metal mineral resources and geothermal energy and other resources, the application provides a two-dimensional magnetotelluric inversion method based on adaptive weight mixed Lp regularization, which takes magnetotelluric data as observation data, constructs an objective function based on mixed Lp regularization, and minimizes the objective function by using a Gauss-Newton optimization method. The model update amount is obtained by using explicit solution of the Gauss-Newton normal equation. In the inversion process, the weights of L2 regularization and L1 regularization are updated by using adaptive weight. The step length is calculated by using Armijo line search, and then a new model is obtained. When the root mean square error (RMS) is less than a set threshold or the maximum iteration number is reached, the inversion calculation is ended, and the final inversion model is output, so as to obtain the underground electrical structure information. The application improves the resolution and accuracy of the inversion result, and can more effectively identify the electrical characteristics of the underground structure.

[0005] In order to achieve the above purpose, the technical scheme of the application is as follows:

[0006] A two-dimensional magnetotelluric inversion method based on adaptive weight mixed Lp regularization, comprising the following steps:

[0007] S1, reading magnetotelluric observation data, discretely constructing an inversion grid by using a quadrilateral finite element grid for the whole inversion region, and constructing an initial model and a prior model of inversion conductivity for the inversion region;

[0008] S2, based on the initial model and the prior model of inversion conductivity and the magnetotelluric observation data, constructing an inversion objective function of mixed Lp regularization including an L1 regularization term, an L2 regularization term and a data fitting term; in the inversion iteration process, the weight factors of the L2 regularization term and the L1 regularization term are adaptively adjusted;

[0009] S3, based on the inversion grid, performing two-dimensional magnetotelluric inversion according to the inversion objective function of mixed Lp regularization, the initial model and the magnetotelluric observation data, outputting the optimized model parameters of Gauss-Newton inversion, and obtaining underground electrical structure information.

[0010] Specifically, the expression of the inversion objective function of mixed Lp regularization is:

[0011]

[0012] wherein, is a data fitting term, and are L2 regularization and L1 regularization terms respectively; W d and W m are data weighting matrix and model weighting matrix respectively, d obsis the magnetotelluric observation data, F(m) is the forward calculation data, m is the inversion model parameter, and m pre is the prior model, λ1 is the weight factor of the L2 regularization term, and λ2 is the weight factor of the L1 regularization term.

[0013] Specifically, the method for outputting the optimized model parameters of Gauss-Newton inversion is:

[0014] Set the maximum number of inversion iterations k max The number and mean square error convergence threshold ε;

[0015] S31, calculate the forward response F(m k ) and the gradient vector g of the inversion objective function k , the Gauss-Newton method is used to minimize the inversion objective function, and the Gauss-Newton normal equation is solved explicitly using direct solution to obtain the model update δm k ;

[0016] S32, using Armijo line search to calculate the model update step size α k , determine the search direction, and update the step size α based on the model k and model update δm k , update the current model to obtain the inversion model;

[0017] S33, calculating the forward response corresponding to the inversion model, and then calculating the root mean square error between the forward calculation data of the inversion model and the observation data;

[0018] S34. If the value of the root mean square error is less than the set mean square error convergence threshold ε, or the maximum number of iterations k is reached max When , the inversion iteration ends and the optimized model parameters of Gauss-Newton inversion are output; otherwise, the process jumps to step S31 to perform the next Gauss-Newton iteration;

[0019] The current model refers to the parameterized model existing before each iteration starts. The current model of the 0th iteration is the input initial model. The current model participating in the nth iteration is obtained based on the inversion model of the n-1th iteration, where n is a natural number greater than 1.

[0020] Specifically, the method for adaptively adjusting the weight factors of the L2 regularization term and the L1 regularization term is:

[0021]

[0022] Among them, λ0 is the initial regularization factor, k is the number of iterations, L2 regularization The initial value of L2 regularization for the kth iteration a value of a weight factor of the L1 regularization, is the L1 regularization an initial value of a weight factor of the L1 regularization, a value of a weight factor of the L1 regularization of the kth iteration, a data fitting term of the k-1th iteration, an L2 regularization term of the k-1th iteration, an L1 regularization term of the k-1th iteration, the variable c is updated by the following formula:

[0023]

[0024] wherein an initial value of the variable c is c0∈[0, 1], and d is a decay factor, and a value of the decay factor d is greater than 1;

[0025] For the L1 regularization alone, the weight λ2 of the L1 regularization is equal to 0, and the weight factor λ1 of the L2 regularization term is updated by the following strategy:

[0026]

[0027] For the L1 regularization alone, the weight factor λ1 of the L2 regularization term is equal to 0, and the weight factor λ2 of the L1 regularization term is updated by the following strategy:

[0028]

[0029] Specifically, an expression of the Gauss-Newton normal equation is as follows:

[0030] H k δm k = -g k

[0031] wherein δm k is a model update amount, g k is a gradient vector of the inversion objective function, and H k is a Hessian matrix.

[0032] Specifically, the method for calculating the model update step α k by using the Armijo line search is as follows: a unit step is tested, and it is determined whether the following inequality is satisfied:

[0033]

[0034] wherein, is a gradient operator, Φ is the inversion objective function, c1 is a variable, and the model update step α kThe initial value of the step size is 1, and if the current search does not satisfy the formula, the step size is gradually reduced in a half-decay manner until the Armijo condition is satisfied.

[0035] Specifically, the calculation formula of the root mean square error is:

[0036]

[0037] where S is the total number of data, is the error of the kth data, is the kth component of the magnetotelluric observation data, (F(m k )) k is the kth component of the forward data.

[0038] Specifically, the data source of the magnetotelluric observation data is model-synthesized data or field-measured data.

[0039] When the data source of the magnetotelluric observation data is model-synthesized data, a finite element numerical method is used for calculation to obtain the magnetotelluric observation data.

[0040] When the data source of the magnetotelluric observation data is field-measured data, the processed observation data is directly read.

[0041] The magnetotelluric observation data includes apparent resistivity and phase.

[0042] Specifically, the method for calculating the magnetotelluric observation data by using the finite element numerical method is:

[0043] The partial differential equation to be solved for two-dimensional forward is:

[0044]

[0045] where Ω is the research area, is the boundary of the research area, is the gradient operator.

[0046] For TE mode:

[0047]

[0048] where E x is the electric field in the x direction, E0 is the electric field at the boundary, μ is the magnetic permeability in vacuum, ω is the angular frequency, i is a pure imaginary number, and σ is the conductivity in the medium.

[0049] For TM mode:

[0050] ​

[0051] Among them, H x is the magnetic field in the x direction, H0 is the magnetic field at the boundary The magnetic field at

[0052] The node basis functions of the quadrilateral are used as weight functions, and the weighted residual equation is obtained using the Galerkin weighted residual method:

[0053]

[0054] Among them, n e is the total number of cells in the grid, N i is the i-th node basis function of the unit quadrilateral;

[0055] The global coordinate system is transformed into the local coordinate system to simplify the subsequent calculation of the unit matrix. In the ξη coordinate system, the four basis functions of the quadrilateral unit are:

[0056]

[0057] After performing element analysis on each quadrilateral element, the following linear equations are obtained:

[0058] k e e e =b e

[0059] Among them, b e is the right end vector of the 4×1 quadrilateral element e, e e is the electric field value at the four nodes corresponding to the quadrilateral element e, k e is the unit matrix corresponding to the quadrilateral unit e;

[0060] For the unit matrix k e Perform global assembly and add Dirichlet boundary conditions at the domain boundaries to obtain the forward linear equations:

[0061] Ke=b

[0062] Among them, the dimension of the coefficient matrix K is NP, NP is the total number of nodes in the forward modeling grid, e is the electric field vector to be determined, and b is the right-hand side vector; the forward linear equations are solved directly; after completing the solution of the forward linear equations for the two polarization modes, the magnetotelluric observation data whose data source is the model synthesis data are obtained.

[0063] Specifically, the unit matrix k e The components are expressed as:

[0064]

[0065] Wherein, x, y are horizontal and vertical coordinates in the global coordinate system respectively, ξ, η are horizontal and vertical coordinates in the local coordinate system respectively, and a and b are the length and width of the quadrilateral element e respectively.

[0066] Compared with the prior art, the present application has the following beneficial effects:

[0067] The present application dynamically adjusts the weights of L2 and L1 regularization by introducing an adaptive weight mechanism, rather than relying on empirical selection, which has the following advantages:

[0068] It can reduce human bias and improve the objectivity and reliability of the inversion results.

[0069] By reasonably allocating the weights of L2 and L1 regularization, overfitting can be effectively prevented, and the generalization ability of the model can be improved.

[0070] The adaptive regularization strategy can be applied to complex geological models, making the changes and characteristics of the geological structure more clear.

[0071] In summary, the present application can improve the resolution of the two-dimensional magnetotelluric inversion results and obtain a more accurate geoelectric model, which is helpful for locating potential mineral resources, geothermal energy and the like, and has important application value. BRIEF DESCRIPTION OF DRAWINGS

[0072] 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 some embodiments of the present application, and for those skilled in the art, other drawings can also be obtained without creative labor on the basis of these drawings.

[0073] Figure 1 The present application is a method flowchart.

[0074] Figure 2 It is a theoretical model diagram in an embodiment of the present application, wherein the black triangle represents the distribution of measuring points.

[0075] Figure 3 It is an inversion resistivity distribution diagram of a theoretical model in an embodiment of the present application based on mixed Lp regularization of L2 regularization, L1 regularization and adaptive weight; wherein, Figure 3 (a) is the L2 regularization inversion result (λ2≡0, λ1 is selected as shown in the following formula (17), wherein λ0=50, d=2), Figure 3 (b) is the L1 regularization inversion result (λ1≡0, λ2 is selected as shown in the following formula (18), wherein λ0=50, d=2), Figure 3-(c) is the hybrid Lp regularization inversion result of the adaptive weight of the present invention (the selection of λ1 and λ2 is shown in the following formula (15), where λ0 = 50, ).

[0076] Figure 4 is the inversion result of magnetotelluric observation data based on mixed Lp regularization with different weights in one embodiment of the present invention, where the data source is model-synthesized data. Figure 4 -(a) is the hybrid Lp regularization inversion result of the adaptive weight of the present invention (the selection of λ1 and λ2 is shown in the following formula (15), where λ0=50, ), Figure 4 -(b) is the hybrid Lp regularization inversion result with adaptive weights (λ1 and λ2 are chosen as shown in formula (15) below, where λ0 = 50 and c is a fixed constant of 0.5), Figure 4 -(c) is the hybrid Lp regularized inversion result with adaptive weights (λ1 and λ2 are chosen as shown in Equation (15) below, where λ0 = 50 and c is a fixed constant of 0.3), Figure 4 -(d) is the mixed Lp regularization inversion result with fixed weights (λ1≡1,λ2≡100), Figure 4 -(e) is the mixed Lp regularization inversion result with fixed weights (λ1≡10,λ2≡10), Figure 4 -(f) is the mixed Lp regularization inversion result with fixed weights (λ1≡10, λ2≡1).

[0077] Figure 5 : is a graph showing the change of the root mean square error of data based on different regularizations with the number of iterations in one embodiment of the present invention. Figure 5 -(a) is Figure 3 The corresponding L1 regularization, L2 regularization and the hybrid Lp regularization of the present invention are the data root mean square error changes with the number of iterations, Figure 5 -(b) Figure 4 The corresponding curve of the root mean square error of the data for the 6 mixed Lp regularizations changes with the number of iterations.

[0078] Figure 6 is the inversion result of the magnetotelluric measured data based on L2 regularization in one embodiment of the present invention (λ2≡0, the selection of λ1 is shown in the following formula (17), where λ0=10 5 ,d=2).

[0079] Figure 7 is the inversion result of the magnetotelluric measured data based on L1 regularization in one embodiment of the present invention (λ1≡0, λ2 is selected as shown in the following formula (18), where λ0=10 5 ,d=2).

[0080] Figure 8 is the inversion result of the mixed Lp regularization based on adaptive weight of the magnetotelluric measured data in an embodiment of the present application (the selection of λ1 and λ2 is shown in the following formula (15), where λ0 = 1000, ).

[0081] Figure 9 is the root mean square error value of the single point data of the magnetotelluric measured data based on L2 regularization, L1 regularization and mixed Lp regularization in an embodiment of the present application. DETAILED DESCRIPTION

[0082] 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 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 the other embodiments obtained by those skilled in the art without creative labor fall within the protection scope of the present application.

[0083] As shown in Figure 1 , a two-dimensional inversion method of magnetotelluric based on mixed Lp regularization of adaptive weight includes the following steps:

[0084] S1, reading magnetotelluric observation data d obs , constructing an initial model m0 and a prior model m of inversion conductivity of the inversion region, and constructing an inversion grid by quadrilateral finite element grid discretization of the entire inversion region. pre .

[0085] When reading the magnetotelluric observation data, the data source is model synthesized data or field measured data. If the data is field observation data, only the processed observation data needs to be read;

[0086] If the data is model synthesized data, in the present embodiment, the model synthesized data is obtained by the theoretical model as shown in Figure 2 , and the finite element numerical method needs to be used for calculation, so as to obtain the magnetotelluric observation data d obs . The steps include:

[0087] Firstly, the partial differential equation to be solved for two-dimensional forward is given as follows:

[0088]

[0089] Wherein, Ω is the research region, is the boundary of the research region, is the gradient operator, and u, u0, τ, λ have different meanings for different polarization modes (TE mode and TM mode), and the specific meanings are as follows.

[0090] For TE mode (Transverse Electric Mode):

[0091]

[0092] where E x is the electric field in x direction, E0 is the electric field at the boundary , μ is the magnetic permeability in vacuum, ω is the angular frequency, i is the imaginary unit, and σ is the electrical conductivity in the medium.

[0093] For TM mode (Transverse Magnetic Mode):

[0094]

[0095] where H x is the magnetic field in x direction, H0 is the magnetic field at the boundary .

[0096] Using the Galerkin weighted residual method, the weighted residual equation is obtained as follows:

[0097]

[0098] where n e is the total number of elements in the mesh, N i is the i-th node function of the quadrilateral element. To simplify the calculation of the element matrix, the global coordinate system (x, y) is transformed to the local coordinate system (ξ, η). In the ξη coordinate system, the four basis functions of the quadrilateral element are:

[0099]

[0100] After element analysis of each quadrilateral element, the following linear equation system is obtained:

[0101] k e e e = b e (6)

[0102] where b e is the 4x1 right-hand vector corresponding to the quadrilateral element e, e e is the electric field value at the four nodes corresponding to the quadrilateral element e, and k e is the element matrix corresponding to the quadrilateral element e, whose dimension is 4x4. The components of the element matrix k e are represented as:

[0103]

[0104] where x, y are the horizontal and vertical coordinates in the global coordinate system, ξ, η are the horizontal and vertical coordinates in the local coordinate system, and a and b are the length and width of the quadrilateral element e, respectively.

[0105] The unit matrix k e is globally assembled and the Dirichlet boundary conditions are added at the region boundaries to obtain the forward linear equations:

[0106] Ke=b (8)

[0107] where the dimension of the coefficient matrix K is NP, NP is the total number of nodes in the forward grid, e is the electric field vector to be solved, and b is the right end vector.

[0108] The direct solution is used to solve the forward linear equations of formula (8). After solving the forward linear equations of the two polarization modes, the magnetotelluric observation data with model-synthesized data as the source, including apparent resistivity and phase data, are obtained.

[0109] In this embodiment, 2% Gaussian random noise is added to the magnetotelluric observation data (apparent resistivity and phase) with model-synthesized data as the source, and the purpose is to simulate the uncertainty of the actual field collected magnetotelluric observation data. The field collected magnetotelluric observation data will be affected by environmental noise, instrument error, etc. By adding noise to the model-synthesized magnetotelluric observation data, the robustness of the inversion algorithm can be tested.

[0110] In this embodiment, the initial model and the prior model of the inversion conductivity of the inversion region are constructed as a uniform half-space (such as a uniform half-space of 100 Ω.m, and the establishment of the specific prior model can also refer to specific geological information).

[0111] S2, based on the initial model m0 and the prior model m pre of the inversion conductivity and the magnetotelluric observation data, an inversion objective function of mixed Lp regularization including L1 regularization term, L2 regularization term and data fitting term is constructed.

[0112] The expression of the inversion objective function of the mixed Lp regularization is:

[0113]

[0114] where, is the data fitting term, and are the L2 regularization and L1 regularization terms, respectively. d and W m are the data weighting matrix and the model weighting matrix, respectively, d obsis the magnetotelluric observation data, F(m) is the forward calculation data, m is the inversion model parameter, m pre is the prior model, λ1 is the weight factor of the L2 regularization term, and λ2 is the weight factor of the L1 regularization term.

[0115] S3, based on the inversion grid, the mixed Lp regularization inversion objective function, the initial model and the magnetotelluric observation data are used to carry out magnetotelluric two-dimensional inversion, and the optimized model parameter of the Gauss Newton inversion is output, so that the underground electrical structure information is obtained;

[0116] The present application combines the traditional L2 regularization and L1 regularization to play the advantages of both. However, how to select the weight factors of different regularizations is crucial to the inversion result. At present, most of the researches on mixed regularization of magnetotelluric inversion are based on fixed weight factors, and the weight of different regularizations depends on experience selection, so that the advantages of mixed regularization cannot be fully played. Therefore, the present application proposes an adaptive weight technology, which adaptively adjusts the weight factors of the L2 regularization term and the L1 regularization term in the inversion iteration process, so as to take into account the advantages of the two kinds of regularizations, fully play the advantages of mixed regularization, and obtain high-resolution inversion results.

[0117] The maximum iteration times k max and the mean square error convergence threshold ε are set.

[0118] S31, the forward response F(m k ) corresponding to the current model and the gradient vector g k of the inversion objective function are calculated. k

[0119] The mixed Lp regularization inversion objective function is minimized by using the Gauss Newton method, and the forward function F(m) at the current model m k needs to be linearized:

[0120] F(m k +δm k )≈F(m k )+J k δm k (10)

[0121] Wherein, δm k is the model update amount, F(m k +δm k ) is the forward calculation data at m k +δm k , and J k is the sensitivity matrix of each iteration. The sensitivity matrix J k ​The calculation employs direct solution.

[0122] After linearization of the forward function F(m), the gradient vector (first derivative) g k and the Hessian matrix (second derivative) H k of the inversion objective function Φ(m) are calculated, which are expressed as follows:

[0123]

[0124] where W y and W z are the first-order difference matrices in the y and z directions, respectively, and R j (j = y, z) are diagonal matrices as follows:

[0125]

[0126] where x j1 , x j2 ,..., x jM are the elements of the vector M is the number of inversion parameters, and η is a very small positive number, such as η = 0.0001, used to control the smoothness of the L1 regularization inversion. A larger value of η corresponds to a smoother inversion result, while a smaller value of η corresponds to a more focused inversion result.

[0127] The Gauss-Newton normal equation is then obtained as follows:

[0128] H k δm k = -g k (14)

[0129] where δm k is the model update, g k is the gradient vector of the inversion objective function, and H k is the Hessian matrix.

[0130] The Gauss-Newton normal equation is solved by the direct method to obtain the model update δm k .

[0131] The method for adaptively adjusting the weight factors of the L2 regularization term and the L1 regularization term during the inversion iteration process is as follows:

[0132]

[0133] where λ0 is the initial regularization factor, k is the iteration number, is the initial value of the L2 regularization , is the L2 regularization The value of the weight factor, L1 regularization The initial value of L1 regularization for the kth iteration The value of the weight factor, is the data fitting term for the k-1th iteration, is the L2 regularization term for the k-1th iteration, is the L1 regularization term for the k-1th iteration, and the variable c is updated as follows:

[0134]

[0135] Among them, the initial value of the variable c is c0∈[0,1], d is the attenuation factor, and the value of the attenuation factor d is greater than 1.

[0136] For a single L2 regularization, which is actually equivalent to the L1 regularization weight λ2 of the mixed regularization, it is always equal to 0. Then the weight factor λ1 of the L2 regularization term is updated using the following strategy:

[0137]

[0138] For separate L1 regularization, the weight factor λ1 of the L2 regularization term is always equal to 0, and the weight factor λ2 of the L1 regularization term is updated using the following strategy:

[0139]

[0140] S32, using Armijo line search to calculate the model update step size α k , determine the search direction, and update the step size α based on the model k and model update δm k , update the current model and obtain the inversion model.

[0141] Armijo line search is used to calculate the model update step size α k The core idea is to test from the unit step size and determine whether the following inequality is satisfied:

[0142]

[0143] Among them, Φ is the inversion objective function, and the variable c1 is generally set to 10 -4 , model update step α k The initial value of is 1. If the current search does not satisfy the above formula, the step size needs to be gradually reduced (halved decay) until the Armijo condition is satisfied.

[0144] The updated model is expressed as:

[0145] mk+1 = m k + a k dm k (20)

[0146] S33, calculate the forward response corresponding to the inversion model, and further calculate the root mean square error between the forward calculation data of the inversion model and the observation data.

[0147] The calculation formula of the root mean square error is:

[0148]

[0149] Wherein, S is the total number of data, is the error of the kth data, is the kth component of the magnetotelluric observation data, (F(m k )) k is the kth component of the forward data.

[0150] S34, if the value of the root mean square error is less than the set mean square error convergence threshold ε, or the maximum iteration number k max is reached, the inversion iteration is ended, and the optimized model parameter of the Gauss-Newton inversion is output, otherwise, jump to step S31 to perform the next Gauss-Newton iteration.

[0151] The current model refers to the parameterized model before each iteration starts, and the current model of the 0th iteration is the input initial model, and the current model participating in the nth iteration is obtained according to the inversion model of the n-1th iteration, and n is a natural number greater than 1.

[0152] Figure 3 (c) and Figure 4 (a) is the mixed Lp regularization inversion optimized model parameter distribution of the adaptive weight of the application. Figure 3 -(a), 3-(b) and Figure 4 -(b)- Figure 4 -(f) distribution is the resistivity distribution diagram of the regularization inversion result of L2 regularization, L1 regularization and other mixed strategies. It can be seen that the adaptive strategy of the application has obvious advantages compared with single mixed regularization and other mixed regularization, and the optimization result of the application is more close to the real model electrical distribution, and the resolution is good.

[0153] Figure 5 -(a) is Figure 3 the root mean square error of the data corresponding to L1 regularization, L2 regularization and mixed Lp regularization of the application changes with the number of iterations, Figure 5 -(b) is Figure 4The RMS error of the data corresponding to the mixed Lp regularization in 6 varies with the iteration number. The experimental results show that basically all the regularization strategies can converge to the threshold 1.02 in a small number of iterations, which shows the superior convergence performance of the Gauss-Newton algorithm, and the adaptive mixed Lp regularization strategy only needs to iterate 3 times to converge.

[0154] For the measured data, the optimized model parameters (reference Figure 6 、 Figure 7 and Figure 8 ) of the Gauss-Newton inversion. Figure 6 is the L2 regularization inversion result based on the magnetotelluric measured data, Figure 7 is the L1 regularization inversion result based on the magnetotelluric measured data, Figure 8 is the adaptive mixed Lp regularization inversion result based on the magnetotelluric. As can be seen from the figure, the inversion result of L2 regularization presents a smooth and diffuse electrical structure, the L1 regularization inversion result obtains a relatively compact and focused inversion result, and the inversion result of the adaptive weight mixed Lp regularization of the present application is relatively smooth for some electrical structure, and is relatively focused for some electrical structure. It is shown that the present application can be applied to actual scenarios containing both smooth features and block features. Figure 9 is the single point data RMS error value based on L2 regularization, L1 regularization and mixed Lp regularization of the magnetotelluric measured data. In the inversion case of the measured data, the RMS of all regularizations corresponding to the single point almost decreases to below 3, and the RMS of some points corresponding to the mixed Lp regularization decreases due to L1 and L2 regularization, which shows that the convergence performance of the present application is good.

[0155] The above inversion test shows that the proposed inversion algorithm can reconstruct both smooth structure and block structure. For smooth structure, it maintains smoothness in the gradual change area and can enhance the boundary resolution in the sudden change area. Further, it can improve the resolution of the magnetotelluric inversion result.

[0156] The above only 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 two-dimensional magnetotelluric inversion method based on hybrid Lp regularization with adaptive weights, characterized in that: The following steps are involved: S1. Read the magnetotelluric observation data, discretize the entire inversion area into a quadrilateral finite element grid to construct an inversion grid, and build the initial model and prior model of the inversion conductivity of the inversion area; S2. Based on the initial model and prior model of the inverted conductivity and the magnetotelluric observation data, a mixed Lp regularized inversion objective function including the L1 regularization term, the L2 regularization term, and the data fitting term is constructed; during the inversion iteration process, the weight factors of the L2 regularization term and the L1 regularization term are adaptively adjusted; S3. Based on the inversion grid, perform two-dimensional magnetotelluric inversion according to the hybrid Lp regularized inversion objective function, the initial model and magnetotelluric observation data, output the optimized model parameters of Gauss-Newton inversion, and obtain the underground electrical structure information.

2. The magnetotelluric two-dimensional inversion method based on hybrid Lp regularization with adaptive weights according to claim 1, characterized in that: The expression of the inversion objective function of the mixed Lp regularization is: in, is the data fitting term, and are L2 regularization and L1 regularization terms respectively; W d and W m are the data weighting matrix and the model weighting matrix, d obs is the magnetotelluric observation data, F(m) is the forward calculation data, m is the inversion model parameter, and m pre is the prior model, λ1 is the weight factor of the L2 regularization term, and λ2 is the weight factor of the L1 regularization term.

3. The magnetotelluric two-dimensional inversion method based on hybrid Lp regularization with adaptive weights according to claim 2, characterized in that: The method for outputting the optimized model parameters of Gauss-Newton inversion is: Set the maximum number of inversion iterations k max The number and mean square error convergence threshold ε; S31, calculate the forward response F(m k ) and the gradient vector g of the inversion objective function k , the Gauss-Newton method is used to minimize the inversion objective function, and the Gauss-Newton normal equation is solved explicitly by direct solution to obtain the model update δm k ; S32, using Armijo line search to calculate the model update step size α k , determine the search direction, and update the step size α based on the model k and model update δm k , update the current model to obtain the inversion model; S33, calculating the forward response corresponding to the inversion model, and then calculating the root mean square error between the forward calculation data of the inversion model and the observation data; S34. If the value of the root mean square error is less than the set mean square error convergence threshold ε, or the maximum number of iterations k is reached max When , the inversion iteration ends and the optimized model parameters of Gauss-Newton inversion are output; otherwise, the process jumps to step S31 to perform the next Gauss-Newton iteration; The current model refers to the parameterized model existing before each iteration starts. The current model of the 0th iteration is the input initial model. The current model participating in the nth iteration is obtained based on the inversion model of the n-1th iteration, where n is a natural number greater than 1.

4. The magnetotelluric two-dimensional inversion method based on hybrid Lp regularization with adaptive weights according to claim 3, characterized in that: The method for adaptively adjusting the weight factors of the L2 regularization term and the L1 regularization term is: Among them, λ0 is the initial regularization factor, k is the number of iterations, L2 regularization The initial value of L2 regularization for the kth iteration The value of the weight factor, L1 regularization The initial value of L1 regularization for the kth iteration The value of the weight factor, is the data fitting term for the k-1th iteration, is the L2 regularization term for the k-1th iteration, is the L1 regularization term for the k-1th iteration, and the variable c is updated as follows: Among them, the initial value of variable c is c0∈[0,1], d is the attenuation factor, and the value of the attenuation factor d is greater than 1; For separate L2 regularization, the weight λ2 of L1 regularization is always equal to 0, then the weight factor λ1 of the L2 regularization term is updated using the following strategy: For separate L1 regularization, the weight factor λ1 of the L2 regularization term is always equal to 0, and the weight factor λ2 of the L1 regularization term is updated using the following strategy:

5. The magnetotelluric two-dimensional inversion method based on hybrid Lp regularization with adaptive weights according to claim 3 or 4, characterized in that: The expression of the Gauss-Newton normal equation is: H k δm k =-g k Among them, δm k is the model update amount, g k is the gradient vector of the inversion objective function, H k is the Hessian matrix.

6. The magnetotelluric two-dimensional inversion method based on hybrid Lp regularization with adaptive weights according to claim 5, characterized in that: The Armijo line search calculation model update step size α k The method is to test from the unit step size and determine whether the following inequality is satisfied: in, is the gradient operator, Φ is the inversion objective function, c1 is the variable, and the model update step size is α k The initial value of is 1. If the current search does not satisfy the above formula, the step size is gradually reduced by half decay until the Armijo condition is satisfied.

7. The magnetotelluric two-dimensional inversion method based on hybrid Lp regularization with adaptive weights according to claim 6, characterized in that: The calculation formula of the root mean square error is: Among them, S is the total number of data, is the error of the kth data, is the kth component of the magnetotelluric observation data, (F(m k )) k is the kth component of the forward modeling data.

8. The magnetotelluric two-dimensional inversion method based on hybrid Lp regularization with adaptive weights according to claim 1 or 7, characterized in that: The data source of the magnetotelluric observation data is model-synthesized data or field-measured data; When the data source of the magnetotelluric observation data is model-synthesized data, a finite element numerical method is used for calculation to obtain the magnetotelluric observation data; When the data source of the magnetotelluric observation data is field measured data, the processed observation data is directly read; The magnetotelluric observation data includes apparent resistivity and phase.

9. The magnetotelluric two-dimensional inversion method based on hybrid Lp regularization with adaptive weights according to claim 8, characterized in that: The method for calculating using the finite element numerical method to obtain magnetotelluric observation data is as follows: The partial differential equation required to be solved in the two-dimensional forward model is: Where Ω is the study area, is the boundary of the study area, is the gradient operator; For TE mode: u=E x ,u0=E0, λ=σ Among them, E x is the electric field in the x direction, E0 is the electric field at the boundary The electric field at , μ is the magnetic permeability in vacuum, ω is the angular frequency, i is a pure imaginary number, and σ is the conductivity of the medium. For TM mode: u=H x ,u0=H0, λ=iωμ Among them, H x is the magnetic field in the x direction, H0 is the magnetic field at the boundary The magnetic field at The node basis functions of the quadrilateral are used as weight functions, and the weighted residual equation is obtained using the Galerkin weighted residual method: Among them, n e is the total number of cells in the grid, N i is the i-th node basis function of the unit quadrilateral; The global coordinate system is transformed into the local coordinate system to simplify the subsequent calculation of the unit matrix. In the ξη coordinate system, the four basis functions of the quadrilateral unit are: After performing element analysis on each quadrilateral element, the following linear equations are obtained: k e e e =b e Among them, b e is the right end vector of the 4×1 quadrilateral element e, e e is the electric field value at the four nodes corresponding to the quadrilateral element e, k e is the unit matrix corresponding to the quadrilateral unit e; For the unit matrix k e Perform global assembly and add Dirichlet boundary conditions at the domain boundaries to obtain the forward linear equations: Ke=b Among them, the dimension of the coefficient matrix K is NP, NP is the total number of nodes in the forward modeling grid, e is the electric field vector to be determined, and b is the right-hand side vector; the forward linear equations are solved directly; after completing the solution of the forward linear equations for the two polarization modes, the magnetotelluric observation data whose data source is the model synthesis data are obtained.

10. The magnetotelluric two-dimensional inversion method based on hybrid Lp regularization with adaptive weights according to claim 9, characterized in that: The unit matrix k e The components are expressed as: Among them, x and y are the horizontal and vertical coordinates in the global coordinate system, ξ and η are the horizontal and vertical coordinates in the local coordinate system, and a and b are the length and width of the quadrilateral element e, respectively.

Images