A smoothing regularization method for 3D electromagnetic inversion on unstructured grids

Abstract

The present invention discloses a practical unstructured grid three-dimensional electromagnetic inversion smoothing regularization method, comprising: obtaining measured electromagnetic data; applying a smoothing regularization operator to construct a three-dimensional electromagnetic inversion regularization objective function, wherein the smoothing regularization operator is constructed based on a model parameter vector under a logarithmic constraint domain; forward calculating the electromagnetic response data of the current geological model; combining the electromagnetic response data and the measured electromagnetic data to calculate the first-order partial derivative of the three-dimensional electromagnetic inversion regularization objective function relative to the model parameter vector under the logarithmic constraint domain, and obtaining the gradient of the current geological model; updating the model parameter vector under the logarithmic constraint domain according to the gradient of the current geological model; converting the updated model parameter vector into the model parameter of the non-logarithmic constraint domain; iterative inversion until the iteration termination condition is met, and outputting the final updated model parameters. The method ensures the stability, convergence, and practicality of three-dimensional electromagnetic inversion using unstructured tetrahedral grids and improves the inversion efficiency.

Description

Technical Field

[0001] The invention belongs to the technical field of geophysical electromagnetic exploration inversion, and in particular relates to a non-structured grid three-dimensional electromagnetic inversion smoothing regularization method. Background Art

[0002] Electromagnetic exploration is an important branch of geophysical exploration. It uses the electrical property differences between most rocks and minerals in underground media as a physical basis. By observing and studying the spatial and temporal distribution laws and changing characteristics of artificial or natural electromagnetic fields caused by electrical property differences, the geoelectric structure distribution of underground media can be determined. This exploration method has the advantages of not being shielded by high-resistance layers and having a large exploration depth. It has been widely used in geoscience fields such as mineral exploration, engineering geophysical exploration, and deep earth structure research. Early processing methods for electromagnetic data were restricted by limited computing hardware and software conditions. The main methods include (1) numerical simulation, such as building a theoretical model for one-dimensional numerical simulation to study the response distribution characteristics and infer which geoelectric structure affects the measured data; (2) one-dimensional or two-dimensional inversion imaging. This method describes the electromagnetic response of the model by constructing a correct electromagnetic field equation, updates the model parameters through an iterative algorithm to make the model's electromagnetic response fit the measured data, and maps the model parameters to obtain the geometric information and electrical information of the inferred underground geoelectric structure. However, the real underground medium should be three-dimensional, and simplified assumptions will cause some results to deviate from the actual situation. Nowadays, with the rapid development of computing software and hardware, three-dimensional electromagnetic inversion research has gradually become a hot topic in this field. This inversion can obtain more accurate geoelectric structure information. However, the three-dimensional electromagnetic inversion problem is more ill-posed than the one-dimensional and two-dimensional inversion. According to the solution of ill-posed problems (Tikhonov AN and Arsenin VY, 1977), solving this problem requires designing a suitable regularization method so that the inversion can converge stably.

[0003] The regularization methods used in the inversion of traditional structured grids are mainly first-order or second-order difference operators, which have achieved good inversion results, such as: three-dimensional conjugate gradient magnetotelluric inversion (Newman GA and Alumbaugh DL, 2000), three-dimensional inversion processing of static effect magnetotelluric data (Sasaki Y., 2004), and computational methods for electromagnetic inversion problems (Egbert GD and Kelbert A., 2012). However, structured grids have limited drawing capabilities when dealing with actual undulating terrain and complex geological models. Unstructured tetrahedral grids have more flexible bevels, and using them for inversion can solve such problems well. Therefore, they are now being studied and applied by more and more scholars, such as: three-dimensional minimized structured magnetotelluric data inversion based on the finite element method of tetrahedral grids (Jahandari H. and Farquharson CG, 2017), L-BFGS three-dimensional magnetotelluric inversion based on unstructured finite elements (Cao XY, Yin CC and Zhang B. et al., 2018), and three-dimensional magnetotelluric inversion of triaxial anisotropic media based on data space (Xie JT, Cai, HZ and Hu, XY et al., 2022).

[0004] However, for unstructured tetrahedral grids, the interface between cells is not perpendicular to the axis, and the conventional model constraint difference operator will result in scattered distribution of model parameters, which will make the inversion unstable. Although the final result can be smoothed by mapping technology, the efficiency of three-dimensional electromagnetic inversion is greatly reduced. Therefore, it is necessary to design a more effective regularization method to make the changes of model parameters smooth in the small area around each cell. Existing unstructured grid regularization methods, such as the volume-weighted model constraint regularization method mentioned in the three-dimensional inversion of transient electromagnetic data based on unstructured tetrahedral grids (Liu YH, Yin C.C. and Qiu CK et al., 2019), refer to the controlled source electromagnetic and magnetotelluric data inversion code (Key K., 2016) for the discretization of triangular grids. The model constraint is used to perform volume-weighted discrete two-norm model constraints on the three-dimensional model discretized by tetrahedral grids. It is a good constraint for grids with relatively uniform volume, but the effect is not ideal for grids with obvious changes in field volume size. The reason is that all regularization methods depend on the selection of regularization factors, and the sensitivity between adjacent units cannot be guaranteed to be smooth in the later stage of inversion. Summary of the invention

[0005] In view of the above, the object of the present invention is to provide a smoothing regularization method for unstructured grid three-dimensional electromagnetic inversion, comprising the following steps:

[0006] Step 1, obtaining measured electromagnetic data;

[0007] Step 2, applying a smoothing regularization operator to construct a three-dimensional electromagnetic inversion regularization objective function, wherein the smoothing regularization operator is constructed based on a model parameter vector under a logarithmic constraint domain;

[0008] Step 3, forward calculation of electromagnetic response data of the current geological model;

[0009] Step 4, combining the electromagnetic response data and the measured electromagnetic data to calculate the first-order partial derivative of the three-dimensional electromagnetic inversion regularization objective function with respect to the model parameter vector under the logarithmic constraint domain, and obtain the gradient of the current geological model;

[0010] Step 5, updating the model parameter vector under the logarithmic constraint domain according to the gradient of the current geological model;

[0011] Step 6, converting the model parameter vector updated in step 5 into model parameters in the non-logarithmic constraint domain;

[0012] Step 7, iterate steps 2-6 until the iteration termination condition is met, and output the final updated model parameters.

[0013] Compared with the prior art, the present invention has the following beneficial effects:

[0014] A smoothing regularization operator is constructed based on the model parameter vector under the logarithmic constraint domain, and a three-dimensional electromagnetic inversion regularization objective function is constructed based on the smoothing regularization operator. During the calculation, the smoothing constraints between the model parameters are not directly displayed, and the model constraint terms are implicitly calculated in the form of a transposed matrix of the inverse of the smoothing matrix, thereby reducing the dependence of the three-dimensional electromagnetic inversion on the regularization factor and making the inversion result smoother, thereby ensuring the stability, convergence, and practicality of the three-dimensional electromagnetic inversion using unstructured tetrahedral grids and improving the inversion efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0015] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work.

[0016] Figure 1 The present invention provides a practical flow chart of a smoothing regularization method for three-dimensional electromagnetic inversion of unstructured grids.

[0017] Figure 2It is a schematic diagram of the geological model provided by an embodiment of the present invention. The black dot on the Z=0km plane is the position of the observation point on the ground, the Z axis is downward and positive, the red and blue squares are the buried positions of the anomalies, the red ones are low-resistance anomalies, and the blue ones are high-resistance anomalies.

[0018] Figure 3 It is an inversion parameter diagram provided by an embodiment of the present invention, wherein (a) is a diagram showing the change of RMS during the inversion process; (b) is a diagram showing the change of the linear search step during the existing volume weighted regularized inversion process; (c) is a diagram showing the change of the regularization factor λ during the inversion process; and (d) is a diagram showing the change of the linear search step during the smooth regularized inversion process proposed by the present invention.

[0019] Figure 4 It is a schematic diagram of the initial gradient distribution of inversion provided by an embodiment of the present invention, wherein (a) is the initial gradient distribution diagram in the existing volume weighted regularized inversion process, and (b) is the initial gradient distribution diagram in the smoothed regularized inversion process proposed by the present invention.

[0020] Figure 5 It is a three-dimensional inversion result slice diagram provided by an embodiment of the present invention. The one marked with original in the figure is the existing volume weighted regularized inversion result slice, and the one marked with p=1.0 is the smoothed regularized inversion result slice proposed by the present invention. DETAILED DESCRIPTION

[0021] To make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific implementation methods described herein are only used to explain the present invention and do not limit the scope of protection of the present invention.

[0022] In order to solve the problem that when using unstructured tetrahedral grid for three-dimensional electromagnetic inversion, the interface between two adjacent units is not perpendicular to the axis, resulting in scattered distribution of model parameters due to uneven grid, which affects the convergence and stability of the inversion and reduces the efficiency of three-dimensional electromagnetic inversion.

[0023] At the same time, in order to solve the problem that smooth regularization methods such as volume weighted model constraint methods cannot well overcome the scattered distribution of model parameters due to their reliance on regularization factors during the inversion process, the embodiment of the present invention proposes a practical unstructured grid three-dimensional electromagnetic inversion regularization method, which does not directly display the smooth constraints between model parameters, but achieves the purpose of smooth inversion by implicitly calculating the model constraint terms by smoothing the sensitivity of the unstructured unit itself. This method has less reliance on the regularization factor and can effectively improve the convergence, stability and computational efficiency of three-dimensional electromagnetic inversion based on unstructured grids.

[0024] Figure 11 is a flow chart of a practical unstructured grid three-dimensional electromagnetic inversion smoothing regularization method provided by an embodiment of the present invention. Figure 1 As shown, the three-dimensional electromagnetic inversion smoothing regularization method provided in the embodiment includes the following steps:

[0025] S1: Read in measured electromagnetic data.

[0026] Different electromagnetic methods correspond to different electromagnetic data. Magnetotelluric sounding data is the total impedance tensor, controlled source audio magnetotelluric data is the apparent resistivity and phase, marine controlled source electromagnetic method is the electric field component, transient electromagnetic method is the rate of change of magnetic induction intensity, etc. However, no matter what type of data is, it should correspond to a measured electromagnetic data vector d obs , whose size is the total amount of data n.

[0027] S2: Apply the smoothing regularization operator to construct the regularized objective function for three-dimensional electromagnetic inversion.

[0028] In the embodiment, the smoothing regularization operator is constructed based on the model parameter vector under the logarithmic constraint domain. The constructed three-dimensional electromagnetic inversion regularization objective function Φ is:

[0029]

[0030] Among them, d obs is the measured electromagnetic data, f(m) is the electromagnetic response data obtained by forward calculation of the current geological model, m is the electrical parameter vector of the model, referred to as the model parameter, generally the conductivity of the unit, C d is the error covariance matrix of the measured electromagnetic data and the electromagnetic response data, represents the square of the L2 norm, λ is the regularization factor, and y is the model parameter vector under the logarithmic constraint domain, which is expressed as:

[0031] y=C w (x inv -x pri )

[0032] Among them, the vector x = T (m), and T () is a conversion function used to control the upper limit b and lower limit a of the model parameter m. For the i-th unit parameter m i , T(m i )=ln(m i -a)-ln(bm i ), the superscript inv indicates that the model parameters are obtained by the current inversion, the superscript pri indicates that the model parameters are the initial model parameters, C w Represents the smoothing matrix. Since y is implicitly calculated, the initial y is equal to 0. During the iteration, the gradient is used to obtain the step size to update y, so the smoothing matrix C w The form is not important.

[0033] S3: Finite element forward modeling calculates the electromagnetic response of the current geological model.

[0034] The finite element method is an effective method for solving partial differential equations. The electromagnetic field equations satisfied by different electromagnetic field data are not the same, but this does not affect the implementation of the smoothing regularization operator proposed in this embodiment. The purpose of this step is to obtain the electromagnetic response data f(m) of the current geological model to compare with the measured electromagnetic data and evaluate whether the model obtained by the current inversion is close to the real underground medium. However, in order to fully illustrate the inversion method proposed in this embodiment, here we only take the finite element simulation of the isotropic magnetotelluric field as an example to briefly introduce the finite element forward simulation process:

[0035] Starting from Maxwell's equations, we can deduce the frequency domain electric field double curl equation that satisfies the isotropic medium magnetotelluric field:

[0036]

[0037] in ω, μ0 refer to the imaginary unit, angular frequency, and magnetic permeability in vacuum (μ0 = 4π×10 -7 T·m / A); E is the electric field strength, unit is V / m, H is the magnetic field strength, unit is A / m; the conductivity σ of the isotropic medium, According to the Galerkin method, the vector interpolation basis function N is used as the weighting function to apply the vector identity to the electric field double curl equation in the solving unit space Ω. The finite element governing equation is obtained:

[0038]

[0039] Following the unit matrix name used in mechanics by the finite element method, the integral term M related to the electrical parameters e is called the mass matrix, The integral term S related to the curl of the interpolation basis function e is called the stiffness matrix, Here, the superscript e indicates that the current physical quantity is a physical quantity inside any unit, the subscript s indicates the local number of the sth edge, and the subscript t indicates the local number of the tth edge corresponding to the sth edge. V is the volume of the tetrahedral unit.

[0040] In order to solve the equation, it is necessary to assemble all units into a general matrix according to the correspondence between the local unit number and the global unit number and apply boundary conditions. At the same time, for the full tensor impedance response Z of the magnetotelluric field, it is necessary to solve the linear equations under two types of orthogonal polarization field sources respectively: KE = b. Here is the overall finite element matrix, E is the electric field vector of the edge to be solved, and b is the right-hand side of the equation after the boundary conditions are applied. The electric field vector E on all unit edges under two types of polarization sources is obtained by solving the equations with the MUMPS direct solver, and then the interpolation operator is defined by Faraday's law of electromagnetic induction to obtain the magnetic field intensity vector H. Then, the full impedance tensor response Z of the magnetotelluric at the current model observation point is obtained by interpolation.

[0041]

[0042] The superscript -1 indicates the inversion operator, while the superscripts 1 and 2 of the electromagnetic fields E and H represent the results obtained by solving the equations under two different polarization sources.

[0043] S4: Calculate the gradient g of the current geological model.

[0044] In the embodiment, the gradient g refers to the first-order partial derivative of the three-dimensional electromagnetic inversion regularization objective function Φ with respect to the model parameter vector y under the logarithmic constraint domain. According to the chain derivative rule, the gradient g of the current geological model is expressed as:

[0045]

[0046] Where v refers to the normalized data fitting difference vector, J is the Jacobian matrix, which stores the first-order derivative of the data response function to the model parameter vector y. The superscript T is the transposition operator. p is a diagonal matrix. For the i-th unit, the diagonal matrix C p The elements of are:

[0047]

[0048] Among them, M represents the total number of units, and exp represents exponential operation. The role of this diagonal matrix is ​​to make the model parameters after the change within a certain range, so as to make the three-dimensional inversion converge more stably. The range (a, b) selected in the specific implementation can be given in combination with the relevant geological materials and geophysical data of the survey area.

[0049] The key to smooth regularization lies in the transposed matrix of the inverse of the smooth matrix How to give, in the embodiment, the inverse distance weighted operator is used to define the transposed matrix of the inverse of the smoothing matrix It is expressed as:

[0050]

[0051] in, is the centroid distance between the ith unit and its jth adjacent unit, N(i) is the number of adjacent units of the ith unit, M represents the total number of units, and the superscript p represents the weight factor of distance weighting, which is a real number greater than 0. Generally, the smaller p is, the greater the influence of the weighted point on the closer distance is.

[0052] When calculating, the inverse of the smoothing matrix is ​​transposed So that the vector J T The elements of each unit in v and its adjacent units tend to be smooth in a small area, while J T v can be obtained by using the adjoint forward modeling method.

[0053] S5: According to the gradient of the current geological model, the L-BFGS algorithm is used to obtain the model parameter vector y in the step-size update transformation domain. k .

[0054] In the embodiment, the calculation of the model parameter update amount is based on the limited memory quasi-Newton algorithm, also known as L-BFGS, with reference to the limited memory quasi-Newton method for large-scale optimization (Liu DC and Nocedal J., 1989).

[0055] The magnetotelluric inversion problem is nonlinear. After being discretized by Taylor expansion, it can be solved by a linear solution algorithm. The solution process is essentially the process of minimizing the objective function Φ. According to Newton's calculation method, the local minimum of a function can be obtained by calculating its first-order derivative and second-order derivative. For three-dimensional electromagnetic inversion problems, the second-order derivative matrix of the objective function with respect to the model parameters is also known as the Hessian matrix, and its computational complexity is astonishing. The L-BFGS algorithm iteratively obtains the approximate inverse matrix D of the Hessian matrix by giving an initial positive definite matrix related to the objective function gradient g, ensuring the descending direction of the inversion, and only storing a finite number of column vector pairs to construct the inverse matrix of the Hessian matrix, instead of storing the inverse matrix of the Hessian matrix, saving a lot of memory and avoiding complex calculations. It is an efficient algorithm for solving the solution to minimize the objective function. The specific process is as follows:

[0056] 1) According to the initial model gradient g0, the inverse matrix of the initial Hessian matrix is ​​set to D0 = I, where I is the unit matrix. At the same time, the iteration step limit δ is set, which is usually a positive integer between 5 and 10;

[0057] 2) Calculate p k =-D k g k , the subscript k indicates that the quantity is in the kth iteration. The iterative formula of D is:

[0058]

[0059] in Model parameter displacement ak =y k+1 -y k , gradient displacement b k =g k+1 -g k Whenever the total number of iterations exceeds δ, the initial vector pair a is deleted from memory. k and b k .

[0060] 3) In a given step interval [α min ,α max ]Inside the linear search, the step length α that satisfies the Wolfe condition is k :

[0061]

[0062] here represents the gradient operator, 1>β2>β1>0.

[0063] 4) New model parameter y k+1 =y k +α k p k .

[0064] S6: Convert the model parameter vector in the transformation domain updated in S5 into the model parameter m in the non-logarithmic constraint domain k .

[0065] In the embodiment, the updated model parameter m k , expressed as:

[0066]

[0067] in, for The transposed matrix of y is the one with rows and columns swapped. k Represents the model parameter vector at the kth iteration. Through this transformation, the inversion can indirectly update the electrical model parameter m required for the actual interpretation, which is generally the conductivity vector of the unit.

[0068] S7: Calculate the RMS.

[0069] For the updated model parameter m k Re-forward modeling is performed to obtain the electromagnetic field response f(m) at the observation point. In order to describe the new electromagnetic response data f(m k ) and measured electromagnetic data d obs The degree of fit is calculated based on the covariance operator matrix Define RMS, the root mean square error of the data, as follows:

[0070]

[0071] The subscript c represents the index of the data element in the data. The elements of the vector er are the estimated noise or error limits of the data.

[0072] S8: Determine whether the three-dimensional electromagnetic inversion reaches the termination condition and whether the regularization factor λ needs to be updated.

[0073] When the termination condition is reached, it is considered that a solution that conforms to the actual underground electrical structure is obtained. At this time, the current inversion model and response are output for geological interpretation, engineering geophysical exploration, and deep earth structure research.

[0074] If the termination condition is not met, first determine whether the regularization factor λ needs to be updated. The judgment condition is whether the reduction between RMS and the RMS of the previous iteration is less than TRMS and whether CL is greater than TCL. If so, CL=0 and the regularization factor λ is updated according to the proportional coefficient k, that is, the new regularization factor λ'=k·λ. Otherwise, continue to iterate and execute S3-S8 until the three-dimensional electromagnetic inversion reaches the termination condition, that is, RMS is less than the expected RMS.

[0075] Here, the expected RMS, the expected RMS reduction TRMS, and the "cooling" step size TCL are all constants given in advance before inversion. The "cooling" step size TCL means that the regularization factor λ is updated only when the RMS cannot be sufficiently reduced after a certain step size TCL is iterated. The expected RMS reduction is generally 0.02.

[0076] Taking the three-dimensional inversion of magnetotelluric data as an example, the three-dimensional electromagnetic inversion smoothing regularization method proposed in the present invention is used to carry out numerical example experiments, and compared with the existing model volume weighted regularization method to illustrate the advantages and effectiveness of the three-dimensional electromagnetic inversion smoothing regularization method proposed in the present invention. In the uniform half space of 250km×250km×150km, there are high and low resistance alternating cube anomalies, such as Figure 2 As shown, the red indicates a low-resistance anomaly with a resistivity of 10Ω·m; the blue indicates a high-resistance anomaly with a resistivity of 1000Ω·m. They are each buried 25km deep, with a size of 100km×100km×50km, and are evenly distributed at 20km intervals. The resistivity of the surrounding rock is 100Ω·m, and the resistivity of the overlying air is 10 8 Ω·m. 289 magnetotelluric measuring points were evenly arranged on the ground with a spacing of 25km. The entire model includes a cube with an extended area of ​​10000km. The forward model is discretized using 706060 unstructured tetrahedral units, in which the measuring point positions and the anomaly positions are locally encrypted. The magnetotelluric impedance tensor response data at the measuring points of the model are obtained by forward calculation using the Galerkin finite element algorithm, and then random Gaussian noise with an impedance amplitude of 3% is added to the data to simulate the measured electromagnetic data.

[0077] The initial geological model of the inversion is the forward model without the abnormal body. The inversion grid contains a total of 594,383 non-structural tetrahedral units, and the inversion area size is 250km×250km×150km. The initial regularization factor λ is 1000, the proportional coefficient k is 0.5, the expected CL is 10, the expected RMS is 1.0, and the expected RMS reduction is 0.02. The synthesized measured data are compared and inverted using the three-dimensional electromagnetic inversion regularization inversion method proposed by the present invention and the existing volume weighted regularization inversion method. The initial conditions and control conditions of the two types of inversion are consistent. When applying the three-dimensional electromagnetic inversion regularization inversion proposed by the present invention, the weight factor P is taken as 1.0.

[0078] The results of the two types of inversion are as follows Figure 3 , Figure 4 and Figure 5 As shown, the three-dimensional electromagnetic inversion regularization method proposed in the present invention is better than the volume weighted regularization method. Figure 3 is the iteration of the parameters in the inversion, Figure 3 (a) is the RMS iteration situation, original indicates the iteration situation of the volume weighted regularization method, and p=1.0 is the iteration situation of the three-dimensional electromagnetic inversion regularization proposed by the present invention. It can be seen that the method of the present invention reduces the RMS to the expected level with fewer iterations. Figure 3 (b) and Figure 3 (d) shows the iterative step sizes of the two types of inversions. By comparison, it can be seen that the inversion method of the present invention has only a few iterative step sizes that are unstable during the entire inversion process, and the step size is 1 in most iterations, indicating that the inversion stability of the inversion method of the present invention is better than that of the existing volume weighted regularization method. Figure 3 (c) shows the change of the inversion regularization factor λ. It can be seen that the inversion of the two methods is normal.

[0079] Figure 4 is the initial gradient distribution of the two types of inversion methods, from Figure 4 (a) and Figure 4 From the comparison in (b), it can be seen that the inversion gradient of the method of the present invention is smoother and more focused than the gradient of the volume weighted regularization method, which is more beneficial for electromagnetic three-dimensional inversion to achieve good convergence and stability.

[0080] Figure 5The horizontal and vertical slices of the three-dimensional inversion results are compared. By comparing the horizontal slices of the two types of inversion results, it can be seen that the existing volume weighted inversion results have a sudden change in model parameters between adjacent units in a very small area, while the resistivity values ​​obtained by the inversion results of the method of the present invention are within the correct range and are smooth and continuous, and accurately depict the horizontal position and morphology of the abnormal body. By comparing the y-direction slices of the two types of inversion results, it can also be seen that the inversion results obtained by the method of the present invention are more focused and have a stronger ability to depict the vertical position and morphology of the abnormal body.

[0081] It should be noted that, for the sake of completeness of the comparative experiment description of the embodiment, the objective function formula of the model volume weighted regularization method is given here:

[0082] Φ=Φ d +λΦ m

[0083]

[0084]

[0085]

[0086] The specific implementation methods described above provide a detailed description of the technical solutions and beneficial effects of the present invention. It should be understood that the above is only the most preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, supplements and equivalent substitutions made within the scope of the principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A smoothing regularization method for unstructured grid three-dimensional electromagnetic inversion, characterized in that: The following steps are involved: Step 1, obtaining measured electromagnetic data; Step 2, applying a smoothing regularization operator to construct a three-dimensional electromagnetic inversion regularization objective function, wherein the smoothing regularization operator is constructed based on a model parameter vector under a logarithmic constraint domain; Step 3, forward calculation of electromagnetic response data of the current geological model; Step 4, combining the electromagnetic response data and the measured electromagnetic data to calculate the first-order partial derivative of the three-dimensional electromagnetic inversion regularization objective function with respect to the model parameter vector under the logarithmic constraint domain, and obtain the gradient of the current geological model; Step 5, updating the model parameter vector under the logarithmic constraint domain according to the gradient of the current geological model; Step 6, converting the model parameter vector updated in step 5 into model parameters in the non-logarithmic constraint domain; Step 7, iterate steps 2-6 until the iteration termination condition is met, and output the final updated model parameters.

2. The unstructured grid three-dimensional electromagnetic inversion smoothing regularization method according to claim 1, characterized in that: The constructed three-dimensional electromagnetic inversion regularization objective function Φ is: Among them, d obs is the measured electromagnetic data, f(m) is the electromagnetic response data obtained by forward modeling of the current geological model, m is the model parameter, C d is the error covariance matrix of the measured electromagnetic data and the electromagnetic response data, represents the square of the L2 norm, λ is the regularization factor, and y is the model parameter vector under the logarithmic constraint domain, which is expressed as: y=C w (x inv -x pri ) Among them, the vector x = T (m), and T () is a conversion function used to control the upper limit b and lower limit a of the model parameter m. For the i-th unit parameter m i , T(m i )=ln(m i -a)-ln(bm i ), the superscript inv indicates that the model parameters are obtained by the current inversion, the superscript pri indicates that the model parameters are the initial model parameters, C w represents the smoothing matrix.

3. The unstructured grid three-dimensional electromagnetic inversion smoothing regularization method according to claim 1, characterized in that: The finite element forward modeling method is used to calculate the electromagnetic response data of the current geological model.

4. The unstructured grid three-dimensional electromagnetic inversion smoothing regularization method according to claim 2, characterized in that: By taking the first-order partial derivative of the regularized objective function of the three-dimensional electromagnetic inversion with respect to the model parameter vector under the logarithmic constraint domain, according to the chain derivative rule, the gradient g of the current geological model is expressed as: Where v refers to the normalized data fitting difference vector, J is the Jacobian matrix, which stores the first-order derivative of the data response function to the model parameter vector y. The superscript T is the transposition operator. p is a diagonal matrix. For the i-th unit, the diagonal matrix C p The elements of are: Where M represents the total number of units, and exp represents exponential operation; When calculating, the inverse of the smoothing matrix is ​​transposed So that the vector J T The elements of each unit in v and its adjacent units tend to be smooth in a small area, while J T v is obtained by using the adjoint forward modeling method.

5. The unstructured grid three-dimensional electromagnetic inversion smoothing regularization method according to claim 4, characterized in that: Define the inverse of the smoothing matrix using the inverse distance weighted operator The form is as follows: in, is the centroid distance between the ith unit and its jth adjacent unit, N(i) is the number of adjacent units of the ith unit, M represents the total number of units, and the superscript p represents the weight factor of distance weighting, which is a real number greater than 0. The smaller p is, the greater the influence of the weighted point on the closer distance is.

6. According to the unstructured grid three-dimensional electromagnetic inversion smoothing regularization method of claim 2, the model parameters updated in step 5 are converted into model parameters m in the non-logarithmic constraint domain. k , expressed as: in, for The transposed matrix of y is the one with rows and columns swapped. k Represents the model parameter vector at the kth iteration.

7. The unstructured grid three-dimensional electromagnetic inversion smoothing regularization method according to claim 1, characterized in that: The iteration termination condition includes the electromagnetic observation data d obs and the electromagnetic response data f(m k ) is expressed as: Among them, the subscript c represents the index of the data element in the data, the elements of the vector er are the data estimation noise, and n is the total number of data elements.

8. The unstructured grid three-dimensional electromagnetic inversion smoothing regularization method according to claim 1, characterized in that: According to the gradient of the current geological model, the L-BFGS algorithm is used to obtain the model parameter vector under the logarithmic constraint domain with a step size update.

Images