CN115755199A - Practical unstructured grid three-dimensional electromagnetic inversion smoothing and regularization method - Google Patents
Practical unstructured grid three-dimensional electromagnetic inversion smoothing and regularization method Download PDFInfo
- Publication number
- CN115755199A CN115755199A CN202211282575.8A CN202211282575A CN115755199A CN 115755199 A CN115755199 A CN 115755199A CN 202211282575 A CN202211282575 A CN 202211282575A CN 115755199 A CN115755199 A CN 115755199A
- Authority
- CN
- China
- Prior art keywords
- inversion
- electromagnetic
- model
- regularization
- data
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Images
Classifications
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/30—Assessment of water resources
Landscapes
- Geophysics And Detection Of Objects (AREA)
Abstract
The invention discloses a practical smoothing and regularizing method for three-dimensional electromagnetic inversion of a non-structural grid, which comprises the following steps: acquiring actually measured electromagnetic data; constructing a three-dimensional electromagnetic inversion regularization target function by using a smoothing regularization operator, wherein the smoothing regularization operator is constructed based on a model parameter vector under a logarithmic constraint domain; calculating the electromagnetic response data of the current geological model in a forward mode; calculating first-order partial derivatives of the three-dimensional electromagnetic inversion regularization target function relative to a model parameter vector under a logarithmic constraint domain by combining electromagnetic response data and actually measured electromagnetic data to obtain the gradient of the current geological model; updating a model parameter vector under a logarithmic constraint domain according to the gradient of the current geological model; converting the updated model parameter vector into a model parameter of a non-log-constrained domain; and (5) carrying out iterative inversion until an iteration termination condition is met, and outputting finally updated model parameters. The method ensures the stability, convergence and practicability of the three-dimensional electromagnetic inversion by utilizing the non-structural tetrahedral mesh and improves the inversion efficiency.
Description
Technical Field
The invention belongs to the technical field of geophysical electromagnetic prospecting inversion, and particularly relates to a practical smoothing and regularizing method for three-dimensional electromagnetic inversion of a non-structural grid.
Background
The electromagnetic exploration method is an important branch in the field of geophysical exploration, and determines the geoelectric structure distribution of an underground medium by observing and researching the space and time distribution rule and the change characteristics of an artificial or natural electromagnetic field caused by electrical property difference by taking the electrical property difference existing among a plurality of rock ores in the underground medium as a physical basis. The exploration method has the advantages of no shielding by a high-resistance layer, large exploration depth and the like, and is widely applied to the field of geoscience such as mineral exploration, engineering geophysical prospecting, earth deep structure research and the like. Early processing means of electromagnetic data are limited by the conditions of computing software and hardware, and the main methods include (1) numerical simulation, for example, by constructing a theoretical model to perform one-dimensional numerical simulation to study the response distribution characteristics and deduce which kind of electrical structure the actually measured data is affected by; (2) The method comprises the steps of constructing a correct electromagnetic field equation to describe the electromagnetic response of a model, updating model parameters through an iterative algorithm to enable the electromagnetic response of the model to be fitted with measured data, and mapping model parameters and deducing geometric information and electrical information of an underground electrical structure. But the real subsurface medium should be three-dimensional, and the simplifying assumption may cause some of the results to deviate from reality. Nowadays, with the rapid development of computing software and hardware, the research of three-dimensional electromagnetic inversion gradually becomes a hot point in the field, and the inversion can obtain more accurate geoelectrical structure information. However, the three-dimensional electromagnetic inversion problem is more ill-conditioned than the one-and two-dimensional inversion, and as can be seen from solutions of ill-defined problems (Tikhonov a.n. and Arsenin v.y., 1977), solving this problem requires designing an appropriate regularization method so that the inversion can converge stably.
Regularization methods adopted in inversion of the traditional structured grid are mainly first-order or second-order difference operators, and the regularization methods achieve good inversion effects, such as: three-dimensional conjugate gradient magnetotelluric inversion (Newman g.a. and Alumbaugh d.l., 2000), three-dimensional inversion processing of static effect magnetotelluric data (Sasaki y, 2004), computational methods of electromagnetic inversion problems (Egbert g.d. and Kelbert a, 2012). However, the structured grid has limited drawing capability when processing actual relief topography and complex geologic body models, and the unstructured tetrahedral grid has more flexible oblique edges, and the inversion performed by using the structured grid can well solve the problems, so that the structured grid is researched and applied by more and more researchers, such as: finite element method three-dimensional minimum structure magnetotelluric data inversion based on tetrahedral mesh (Jahandari h. And farquhardson c.g., 2017), L-BFGS three-dimensional magnetotelluric inversion based on non-structural finite elements (Cao x.y., yin c.c. and Zhang b. Etc., 2018), and triaxial anisotropic medium magnetotelluric three-dimensional inversion based on data space (Xie j.t., cai, h.z. And Hu, x.y. Etc., 2022).
However, for an unstructured tetrahedral mesh, the interface between units is not perpendicular to the axial direction, and the phenomenon of scattered distribution of model parameters can occur by adopting a conventional model constraint difference operator, so that the inversion tends to be unstable, and although the final result can be smoothly processed by the mapping technology, the efficiency of three-dimensional electromagnetic inversion is greatly reduced. Therefore, there is a need to design more efficient regularization methods such that the variation of the model parameters tends to be smooth in a small field around each cell. The existing unstructured grid regularization methods, such as: the regularization methods based on volume weighting for model constraint mentioned in transient electromagnetic data three-dimensional inversion (Liu y.h., yin c.c. and Qiu c.k. and so on, 2019) based on non-structural tetrahedral mesh refer to the discrete model inversion of triangular mesh in the inversion code (Key k., 2016) of controllable source electromagnetic and magnetotelluric data, the discrete two-norm model constraint of area weighting is adopted to perform the discrete two-norm model constraint of volume weighting on the discrete three-dimensional model of tetrahedral mesh, which is good for mesh with uniform volume, but the effect is not ideal for mesh with obviously changed field volume, because all regularization methods depend on the selection of regularization factors, and when the late inversion stage ensures that the sensitivity between adjacent units is smooth.
Disclosure of Invention
In view of the above, the present invention provides a practical method for smoothing and regularizing three-dimensional electromagnetic inversion of unstructured grids, comprising the following steps:
step 2, a smooth regularization operator is applied to construct a three-dimensional electromagnetic inversion regularization target function, wherein the smooth regularization operator is constructed based on a model parameter vector under a logarithmic constraint domain;
step 3, calculating the electromagnetic response data of the current geological model in a forward mode;
step 4, calculating a first-order partial derivative of the three-dimensional electromagnetic inversion regularization target function relative to a model parameter vector under a logarithmic constraint domain by combining the electromagnetic response data and the actually measured electromagnetic data to obtain the gradient of the current geological model;
step 6, converting the model parameter vector updated in the step 5 into a model parameter of a non-logarithmic constraint domain;
and 7, iterating the steps 2-6 until an iteration termination condition is met, and outputting a finally updated model parameter.
Compared with the prior art, the invention has the beneficial effects that at least:
the method comprises the steps of constructing a smooth regularization operator based on model parameter vectors in a logarithmic constraint domain, constructing a three-dimensional electromagnetic inversion regularization target function based on the smooth regularization operator, during calculation, not directly displaying smooth constraints among model parameters, implicitly calculating a model constraint term by adopting a form of calculating a transpose matrix of a smooth matrix inverse, reducing the dependence of the three-dimensional electromagnetic inversion on regularization factors, enabling an inversion result to be smoother, and accordingly ensuring the stability, convergence and practicability of the three-dimensional electromagnetic inversion by utilizing a non-structural tetrahedral mesh and improving the inversion efficiency.
Drawings
In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the embodiments or the prior art will be briefly described below, it is obvious that the drawings in the following description are only some embodiments of the present invention, and for those skilled in the art, other drawings can be obtained according to these drawings without creative efforts.
Fig. 1 is a flowchart of a practical three-dimensional electromagnetic inversion smoothing regularization method for an unstructured grid according to an embodiment of the present invention.
Fig. 2 is a schematic diagram of a geological model provided in an embodiment of the present invention, where a black point on a Z =0km plane is a position of an observation point on the ground, a Z-axis direction is positive, a red and blue square is a buried position of an anomalous body, red is a low-resistance anomalous body, and blue is a high-resistance anomalous body.
FIG. 3 is a graph of inversion parameters provided by an embodiment of the present invention, wherein (a) is a graph of RMS variation during inversion; (b) A change graph of the linear search step length in the existing volume weighting regularization inversion process is obtained; (c) is a variation graph of the regularization factor lambda in the inversion process; (d) The invention provides a change graph of linear search step length in the smoothing and regularization inversion process.
Fig. 4 is a schematic diagram of inversion initial gradient distribution provided in an embodiment of the present invention, where (a) is an initial gradient distribution map in an existing volume-weighted regularization inversion process, and (b) is an initial gradient distribution map in a smooth regularization inversion process provided in the present invention.
FIG. 5 is a slice diagram of the results of the three-dimensional inversion provided by the embodiment of the invention. The conventional volume weighted regularized inversion result slice is marked with original in the figure, and the smoothed regularized inversion result slice provided by the invention is marked with p = 1.0.
Detailed Description
In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention will be further described in detail with reference to the accompanying drawings and examples. It should be understood that the detailed description and specific examples, while indicating the scope of the invention, are intended for purposes of illustration only and are not intended to limit the scope of the invention.
The method aims to solve the problems that when three-dimensional electromagnetic inversion is carried out by using a non-structural tetrahedral mesh, the convergence and the stability of the inversion are influenced and the efficiency of the three-dimensional electromagnetic inversion is reduced due to the fact that the interface between two adjacent units is not vertical to the axial direction and model parameters are scattered due to non-uniform meshes.
Meanwhile, the method also aims to solve the problem that the dependence on regularization factors cannot well overcome the scattered distribution of model parameters in the inversion process of a smooth regularization method such as a volume weighted model constraint method. The embodiment of the invention provides a practical regularization method for three-dimensional electromagnetic inversion of a non-structural grid, which does not directly display the smoothness constraint among the parameters of a used model, and implicitly calculates the model constraint item through the sensitivity of a smooth non-structural unit to achieve the aim of smooth inversion.
Fig. 1 is a flowchart of a practical three-dimensional electromagnetic inversion smoothing regularization method for an unstructured grid according to an embodiment of the present invention. As shown in fig. 1, an embodiment provides a smoothing and regularizing method for three-dimensional electromagnetic inversion, which includes the following steps:
s1: and reading in measured electromagnetic data.
Different electromagnetic methods correspond to different electromagnetic data, magnetotelluric sounding data are full impedance tensors, controllable source audio magnetotelluric data are apparent resistivity and phase, an ocean controllable source electromagnetic method is electric field component, a transient electromagnetic method is change rate of magnetic induction intensity and the like, but whatever type of data corresponds to a measured electromagnetic data vector d obs The size is the total data amount n.
S2: and applying a smooth regularization operator to construct a three-dimensional electromagnetic inversion regularization target function.
In an embodiment, the smoothing regularization operator is constructed based on a model parameter vector under a log-constrained domain. The constructed three-dimensional electromagnetic inversion regularization target function phi is as follows:
wherein d is obs For actually measured electromagnetic data, f (m) is electromagnetic response data obtained by forward modeling calculation of the current geological model, m is an electric parameter vector of the model, which is called model parameter for short, and is generally the conductivity of the unit，C d An error covariance matrix of the measured electromagnetic data and the electromagnetic response data,
the square of the L2 norm is expressed, λ is the regularization factor, and y is the model parameter vector under the log-constrained domain, expressed as:
y＝C w (x inv -x pri )
where the vector x = T (m) and T () is a transfer function controlling the upper and lower limits b, a of the model parameter m for the ith unit parameter m i ，T(m i )＝ln(m i -a)-ln(b-m i ) The superscript inv represents that the model parameters are obtained by current inversion, the superscript pri represents that the model parameters are initial model parameters, C w A smoothing matrix is represented. Since y is implicitly calculated, the initial y is equal to 0, and y is updated by the gradient solving step size in the iterative process, the matrix C is smoothed w The form of (a) is not critical.
And S3, forward computing the electromagnetic response of the current geological model by using a finite element.
The finite element method is an effective method for solving partial differential equations, the electromagnetic field equations which are satisfied by different electromagnetic field data are different, but the implementation of the smoothing regularization operator provided by the embodiment is not influenced, and the purpose of the step is to obtain electromagnetic response data f (m) of the current geological model to compare with actually measured electromagnetic data and evaluate whether the model obtained by current inversion is close to a real underground medium. However, in order to fully explain the inversion method proposed in this embodiment, only the finite element simulation of the isotropic earth electromagnetic field is taken as an example, and a finite element forward simulation flow is briefly introduced:
the frequency domain electric field double rotation degree equation which meets the isotropic medium earth electromagnetic field can be derived from the Maxwell equation:
wherein
ω，μ 0 Respectively, the imaginary unit, the angular frequency and the magnetic permeability (mu) in vacuum 0 ＝4π×10 -7 T.m/A); e is the electric field intensity with the unit of V/m, H is the magnetic field intensity with the unit of A/m; the electrical conductivity a of the isotropic medium,
is a rotation operator. According to the Galerkin method, vector interpolation basis function N is used as a weighting function, and a vector identity equation is applied to an electric field double-rotation equation in the solution of unit space omega
Obtaining a finite element control equation:
calling the element matrix in mechanics by finite element, integral term M related to electric parameter e Referred to as the quality matrix, is,
s, t =1, …,6, and an integral term S related to the rotation of the interpolation basis function e Referred to as a matrix of stiffness values,
s, t =1, …,6. The superscript e here indicates that the current physical quantity is a physical quantity inside any one cell, the subscript s indicates the local number of the s-th edge, and the subscript t indicates the local number of the corresponding t-th edge of the s-th edge. V is the volume of the tetrahedral unit.
In order to solve the equation, all the units need to be assembled into an overall matrix according to the corresponding relation between the local number and the global number of the unit, and a boundary condition is applied, and simultaneously, for the full tensor impedance response Z of the earth electromagnetic field, linear equation sets under two types of orthogonal polarization field sources need to be solved respectively: KE = b. Here, the
And E is an edge electric field vector to be solved, and b is an equation right-end term after the boundary condition is applied. And solving an equation set through an MUMPS direct solver to obtain electric field vectors E on all unit edges under the two types of polarization sources, and defining an interpolation operator through a Faraday electromagnetic induction law to obtain a magnetic field intensity vector H. And then interpolating to obtain the full impedance tensor response Z of the magnetotelluric at the observation point of the current model.
The superscript-1 indicates the inverse, while the superscripts 1 and 2 for the electromagnetic fields E and H, respectively, represent the solution to the equation for two different sources of polarization.
S4: the gradient g of the current geological model is calculated.
In the embodiment, the gradient g refers to a first-order partial derivative of a three-dimensional electromagnetic inversion regularization target function Φ relative to a model parameter vector y under a logarithmic constraint domain, and the gradient g of the current geological model is obtained according to a chain derivative rule and is expressed as:
where v refers to the normalized data fitting difference vector,
j is Jacobian matrix, the first derivative of the data response function to the model parameter vector y is stored, superscript T is the transpose operator, C p As diagonal matrix, for the ith cell, diagonal matrix C p The elements of (A) are:
where M denotes the total number of cells and exp denotes the exponential operation. The function of the diagonal matrix is to make the model parameters within a certain range after being changed, so as to make the three-dimensional inversion more stably converged, and the range (a, b) in the specific implementation can be selected and given by combining geological materials and geophysical data related to a measured area.
The key to smooth regularization is the transpose of the inverse of the smooth matrix
How given, the inverse distance weighting operator is used in the embodiment to define the transpose of the inverse of this smooth matrix
Expressed as:
wherein, the first and the second end of the pipe are connected with each other,
the centroid distance between the ith unit and the jth adjacent unit is defined, N (i) is the number of adjacent units of the ith unit, M represents the total number of units, superscript p represents a weighting factor for distance weighting, a real number greater than 0 is taken, and generally, the smaller p is, the more the influence of the closer distance is exerted on a weighting point.
In the calculation, the inverted matrix is inverted by the smooth matrix
So that vector J T v the elements of each unit and its neighbors tend to be smooth in a small domain, and J T v can be obtained by a forward method.
S5: calculating a model parameter vector y under the step updating conversion domain by adopting an L-BFGS algorithm according to the gradient of the current geological model k 。
In the examples, the calculation of the model parameter update is based on a finite memory quasi-newton algorithm, also known as L-BFGS, referred to as the finite memory quasi-newton method for large scale optimization (Liu d.c. and Nocedal j., 1989).
The magnetotelluric inversion problem is nonlinear, the linear solving algorithm can be adopted for solving after Taylor expansion discretization, the solving process is substantially the process of minimizing an objective function phi, and the local minimum of a function can be obtained by solving a first derivative and a second derivative of the function according to a Newton calculation method. For three-dimensional electromagnetic inversion problems, the second derivative matrix of the objective function with respect to the model parameters, also known as hessian matrix, is prohibitively computationally intensive. The L-BFGS algorithm iteratively obtains an inverse matrix D of an approximate Hessian matrix by giving an initial positive definite matrix related to the gradient g of the objective function, the descending direction of inversion is ensured, the inverse matrix of the Hessian matrix is constructed by only storing column vector pairs with limited steps, the inverse matrix of the Hessian matrix is not stored, a large amount of memory is saved, complex calculation is avoided, and the L-BFGS algorithm is an efficient algorithm for solving the solution of the minimized objective function. The specific process is as follows:
1) Gradient g according to initial model 0 Setting an inverse D of the initial Hessian matrix 0 I, I is the identity matrix. Meanwhile, setting an iteration step number limit delta which is usually a positive integer between 5 and 10;
2) Calculating p k ＝-D k g k The subscript k indicates that the quantity is at the kth iteration. The iterative formula of D is:
wherein
Displacement a of model parameter k ＝y k+1 -y k Gradient displacement b k ＝g k+1 -g k . Deleting the original vector pair a in the memory whenever the total number of iteration steps exceeds delta k And b k 。
3) At a given step interval [ alpha ] min ,α max ]The step length alpha of the internal linear search satisfying the Wolfe condition k ：
Here, the
Representing gradient operators, 1 > beta 2 ＞β 1 ＞0。
4) New model parameters y k+1 ＝y k +α k p k 。
S6: converting the model parameter vector under the conversion domain updated by S5 into the model parameter m of the non-logarithmically constrained domain k 。
In an embodiment, the updated model parameter m k Expressed as:
wherein the content of the first and second substances,
is composed of
Transposed matrices, i.e. interchanging rows and columns, y k The model parameter vector at the k-th iteration is represented, and the inversion can indirectly update the electrical model parameter m required by the actual interpretation through the conversion, which is generally the conductivity vector of the unit.
S7: the RMS is calculated.
For the updated model parameter m k And performing forward modeling once again to obtain the electromagnetic field response f (m) at the observation point. To describe the new electromagnetic response data f (m) k ) With measured electromagnetic data d obs According to the covariance operator matrix
The RMS, i.e. the root mean square error of the data, is defined as follows:
where subscript c denotes 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.
S8: and judging whether the three-dimensional electromagnetic inversion reaches a termination condition and whether the regularization factor lambda needs to be updated.
And when the termination condition is reached, the solution which accords with the underground real geoelectric structure is considered to be obtained, and the current inversion model and the response are output to be used for the research of geological interpretation, engineering geophysical prospecting, earth deep structure and the like.
If the termination condition is not met, judging whether the regularization factor lambda needs to be updated or not, wherein the judgment condition is that whether the reduction amount between the RMS and the RMS of the last iteration is smaller than TRMS and whether CL is larger than TCL or not, if CL =0 and updating the regularization factor lambda according to a proportionality coefficient k, the new regularization factor lambda' = k & lambda is obtained, otherwise, continuously iterating and executing S3-S8 until the three-dimensional electromagnetic inversion reaches the termination condition, and the RMS is smaller than the expected RMS.
Here the expected RMS, the expected RMS reduction TRMS and the "cool down" step length TCL are all constants that are pre-given prior to inversion. The "cool down" step size TCL means that the regularization factor λ is selected to be updated only if the RMS cannot be sufficiently reduced after the iteration over a certain step size TCL. The expected RMS reduction is typically 0.02.
In the following, by taking magnetotelluric data three-dimensional inversion as an example, numerical example experiments are carried out by applying the three-dimensional electromagnetic inversion smoothing regularization method provided by the invention, and compared with the existing model volume weighting regularization method, the advantages and effectiveness of the three-dimensional electromagnetic inversion smoothing regularization method provided by the invention are clarified. A high-low resistance alternating cube anomaly is scattered in a uniform half space of 250km multiplied by 150km, as shown in fig. 2, wherein red indicates the low resistance anomaly and the resistivity is 10 Ω · m; blue is a high-resistance anomalous body, and the resistivity is 1000 Ω · m. The buried depth is 25km, the buried depth is 100km multiplied by 50km, the buried depth is 20km and the buried depth is equal to 20km, the surrounding rock resistivity is 100 omega.m, the overlying air resistivity is 10 8 Omega.m. 289 magnetotelluric measuring points are uniformly arranged on the ground at an interval of 25km. The whole model comprises a cube with an edge expanding area of 10000km, and the forward model is discretized by utilizing 706060 unstructured tetrahedron units, wherein the measuring point position and the abnormal body position are locally encrypted. And forward calculation is carried out by utilizing a Galerkin finite element algorithm to obtain magnetotelluric impedance tensor response data on the model measuring points, and then random Gaussian noise with 3% impedance amplitude is added to the data to simulate the actually measured electromagnetic data.
The inversion initial geological model removes abnormal bodies for the forward model, the inversion grid comprises 594383 non-structural tetrahedral units, and the size of the inversion area is 250km multiplied by 150km. The initial regularization factor λ is 1000, the scaling factor k is 0.5, CL is expected to be 10, RMS is expected to be 1.0, and RMS reduction is expected to be 0.02. The synthesized measured data is compared and inverted by the three-dimensional electromagnetic inversion regularization inversion method and the existing volume weighting regularization inversion method, and the initial conditions and the control conditions of the two types of inversion are consistent. When the regularized inversion of the three-dimensional electromagnetic inversion provided by the invention is applied, the weight factor P is 1.0.
The results of the two types of inversion are shown in fig. 3, fig. 4 and fig. 5, and the effect of the three-dimensional electromagnetic inversion regularization method provided by the invention is superior to that of the volume weighting regularization method. Fig. 3 is an iteration case of parameters in inversion, fig. 3 (a) is an RMS iteration case, original indicates an iteration case of the volume weighted regularization method, and p =1.0 indicates an iteration case of the regularization of the three-dimensional electromagnetic inversion proposed by the present invention, and it can be seen that the RMS is reduced to a desired value with a smaller number of iterations in the method of the present invention. In fig. 3, (b) and fig. 3, (d) are the two types of inversion iteration step size cases, and a comparison shows that the inversion stability of the method is better than that of the existing volume weighted regularization method, because the method of the present invention has unstable iteration step sizes for a few times in the whole inversion process, and the step sizes in most iterations are 1. Fig. 3 (c) shows the inversion regularization factor λ variation, and it can be seen that the inversion conditions of both methods are normal.
Fig. 4 is the initial gradient distribution of the two types of inversion methods, and it can be seen from the comparison between (a) in fig. 4 and (b) in fig. 4 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 to the electromagnetic three-dimensional inversion to obtain good convergence and stability.
FIG. 5 shows horizontal and vertical slices of three-dimensional inversion results, and comparing the horizontal slices of two types of inversion results shows that the existing volume weighted inversion result has sudden changes of model parameters between adjacent units in a small area, and the resistivity value obtained by the inversion result of the method is in a correct range, smooth and continuous, and simultaneously accurately depicts the horizontal position form of an abnormal body. The same can be seen by comparing the y-direction slices of the two types of inversion results, the inversion results obtained by the method are more focused, and the position and shape depicting capability on the abnormal body in the vertical direction is stronger.
It should be noted that, for completeness of description of the comparative experiment of the embodiment, an objective function formula of the applied model volume weighted regularization method is given here:
Φ＝Φ d +λΦ m
the above-mentioned embodiments are intended to illustrate the technical solutions and advantages of the present invention, and it should be understood that the above-mentioned embodiments are only the most preferred embodiments of the present invention, and are not intended to limit the present invention, and any modifications, additions, equivalents, etc. made within the scope of the principles of the present invention should be included in the scope of the present invention.
Claims (8)
1. A practical smoothing and regularizing method for three-dimensional electromagnetic inversion of unstructured grids is characterized by comprising the following steps:
step 1, acquiring actually measured electromagnetic data;
step 2, a smoothing regularization operator is applied to construct a three-dimensional electromagnetic inversion regularization target function, wherein the smoothing regularization operator is constructed on the basis of a model parameter vector under a logarithmic constraint domain;
step 3, calculating the electromagnetic response data of the current geological model in a forward mode;
step 4, calculating a first-order partial derivative of the three-dimensional electromagnetic inversion regularization target function relative to a model parameter vector under a logarithmic constraint domain by combining the electromagnetic response data and the actually measured electromagnetic data to obtain the gradient of the current geological model;
step 5, updating model parameter vectors under a logarithmic constraint domain according to the gradient of the current geological model;
step 6, converting the model parameter vector updated in the step 5 into a model parameter of a non-logarithmic constraint domain;
and 7, iterating the steps 2-6 until an iteration termination condition is met, and outputting the finally updated model parameters.
2. The practical unstructured grid three-dimensional electromagnetic inversion smoothing regularization method according to claim 1, characterized in that the constructed three-dimensional electromagnetic inversion regularization target function Φ is:
wherein d is obs For actually measured electromagnetic data, f (m) is electromagnetic response data obtained by forward modeling calculation of the current geological model, m is a model parameter, C d An error covariance matrix of the measured electromagnetic data and the electromagnetic response data,
expressing the square of L2 norm, lambda is a regularization factor, and y is in a logarithmically constrained domainIs expressed as:
y＝C w (x inv -x pri )
where the vector x = T (m) and T () is a transfer function controlling the upper and lower limits b, a of the model parameter m for the ith unit parameter m i ，T(m i )＝ln(m i -a)-ln(b-m i ) The superscript inv represents that the model parameters are obtained by current inversion, the superscript pri represents that the model parameters are initial model parameters, C w A smoothing matrix is represented.
3. The practical unstructured grid three-dimensional electromagnetic inversion smoothing and regularization method according to claim 1, characterized in that a finite element forward modeling method is used to calculate electromagnetic response data of the current geological model.
4. The practical unstructured grid three-dimensional electromagnetic inversion smoothing and regularization method according to claim 2, characterized in that the gradient g of the current geological model is obtained by first-order partial derivatives of the three-dimensional electromagnetic inversion regularization target function relative to the model parameter vector under the logarithmic constraint domain according to a chain derivation rule as follows:
where v refers to the normalized data fitting difference vector,
j is Jacobian matrix, the first derivative of the data response function to the model parameter vector y is stored, superscript T is the transpose operator, C p As diagonal matrix, for the ith cell, diagonal matrix C p The elements of (A) are:
where M denotes the total number of cells and exp denotes the exponential operation.
In the calculation, the inverted matrix is inverted by the smooth matrix
So that the vector J T v the elements of each unit and its neighbors tend to be smooth in a small domain, and J T v can be obtained by a forward method.
5. The method according to claim 4, wherein the inverse distance weighting operator is used to define a transpose of the smooth matrix inverse
The form is as follows:
wherein the content of the first and second substances,
the centroid distance between the ith unit and the jth adjacent unit is defined, N (i) is the number of adjacent units of the ith unit, M represents the total number of units, superscript p represents a weighting factor for distance weighting, a real number greater than 0 is taken, and generally, the smaller p is, the more the influence of the closer distance is exerted on a weighting point.
6. The practical smoothing and regularizing method for three-dimensional electromagnetic inversion of unstructured grids as set forth in claim 2, wherein the model parameters updated in step 5 are converted into model parameters m of a non-log constrained domain k Expressed as:
wherein, the first and the second end of the pipe are connected with each other,
is composed of
Transposed matrices, i.e. interchanging rows and columns, y k Representing the model parameter vector at the k-th iteration.
7. The method of claim 1, wherein the iteration termination condition comprises electromagnetic observation data d obs Electromagnetic response data f (m) with the k-th iteration k ) RMS error, expressed as:
the subscript c represents the index of the data elements in the data, the elements of the vector er are the data estimation noise, and n is the total amount of the data elements.
8. The practical unstructured grid three-dimensional electromagnetic inversion smoothing and regularization method according to claim 1, characterized in that a L-BFGS algorithm is adopted to solve model parameter vectors under a step updating logarithmic constraint domain according to the gradient of a current geological model.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211282575.8A CN115755199A (en) | 2022-10-19 | 2022-10-19 | Practical unstructured grid three-dimensional electromagnetic inversion smoothing and regularization method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211282575.8A CN115755199A (en) | 2022-10-19 | 2022-10-19 | Practical unstructured grid three-dimensional electromagnetic inversion smoothing and regularization method |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN115755199A true CN115755199A (en) | 2023-03-07 |
Family
ID=85352150
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202211282575.8A Pending CN115755199A (en) | 2022-10-19 | 2022-10-19 | Practical unstructured grid three-dimensional electromagnetic inversion smoothing and regularization method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115755199A (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116908928A (en) * | 2023-05-15 | 2023-10-20 | 重庆大学 | Stratum adaptive encryption-based magnetotelluric inversion method |
-
2022
- 2022-10-19 CN CN202211282575.8A patent/CN115755199A/en active Pending
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116908928A (en) * | 2023-05-15 | 2023-10-20 | 重庆大学 | Stratum adaptive encryption-based magnetotelluric inversion method |
| CN116908928B (en) * | 2023-05-15 | 2024-03-26 | 重庆大学 | Stratum adaptive encryption-based magnetotelluric inversion method |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN113221393B (en) | Three-dimensional magnetotelluric anisotropy inversion method based on non-structural finite element method | |
| Key et al. | Adaptive finite-element modeling using unstructured grids: The 2D magnetotelluric example | |
| CN106980736B (en) | A kind of ocean controllable source electromagnetic method finite element forward modeling method of anisotropic medium | |
| CN112949134B (en) | Earth-well transient electromagnetic inversion method based on non-structural finite element method | |
| Rezaie et al. | Fast 3D focusing inversion of gravity data using reweighted regularized Lanczos bidiagonalization method | |
| CN113255230B (en) | Gravity model forward modeling method and system based on MQ radial basis function | |
| Wang et al. | 2D joint inversion of CSAMT and magnetic data based on cross-gradient theory | |
| Cao et al. | 3D magnetotelluric inversions with unstructured finite-element and limited-memory quasi-Newton methods | |
| CN115755199A (en) | Practical unstructured grid three-dimensional electromagnetic inversion smoothing and regularization method | |
| CN115238550A (en) | Self-adaptive unstructured grid landslide rainfall geoelectric field numerical simulation calculation method | |
| Wang et al. | Three-dimensional scalar controlled-source audio-frequency magnetotelluric inversion using tipper data | |
| CN113569447B (en) | Transient electromagnetic three-dimensional fast forward modeling method based on Shu' er compensation method | |
| Persova et al. | Geometric 3-D inversion of airborne time-domain electromagnetic data with applications to kimberlite pipes prospecting in a complex medium | |
| CN116466402B (en) | Electromagnetic inversion method based on geological information and electromagnetic data combined driving | |
| Chen et al. | An application of Coons patch to generate grid-based digital elevation models | |
| CN114970289B (en) | Three-dimensional magnetotelluric anisotropy forward modeling numerical simulation method, equipment and medium | |
| Peng et al. | Efficient 3D inversion of potential field data using fast proximal objective function optimization algorithm | |
| Boonchaisuk et al. | Two-dimensional direct current (DC) resistivity inversion: Data space Occam's approach | |
| Yang et al. | Three‐dimensional inversion of semi‐airborne transient electromagnetic data based on finite element method | |
| Zhong et al. | Electrical resistivity tomography with smooth sparse regularization | |
| Lu et al. | 3D large-scale transient electromagnetic modeling using a Shift-and-Invert Krylov subspace method | |
| CN114200541A (en) | Three-dimensional gravity-magnetic joint inversion method based on cosine dot product gradient constraint | |
| Huo et al. | 3D nonlinear conjugate gradient inversion for frequency-domain airborne EM based on vector finite element method | |
| CN107710023B (en) | Efficient solution of inverse problems | |
| Gomes et al. | Efficient Calculation of the Jacobian Matrix for 3D inversion of MT Data |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |