CN115113286B - Aviation electromagnetic data fusion three-dimensional inversion method based on multi-component frequency domain - Google Patents
Aviation electromagnetic data fusion three-dimensional inversion method based on multi-component frequency domain Download PDFInfo
- Publication number
- CN115113286B CN115113286B CN202210799540.5A CN202210799540A CN115113286B CN 115113286 B CN115113286 B CN 115113286B CN 202210799540 A CN202210799540 A CN 202210799540A CN 115113286 B CN115113286 B CN 115113286B
- Authority
- CN
- China
- Prior art keywords
- inversion
- data
- frequency domain
- model
- aviation electromagnetic
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 103
- 230000004927 fusion Effects 0.000 title claims abstract description 29
- 239000011159 matrix material Substances 0.000 claims abstract description 37
- 230000006870 function Effects 0.000 claims description 30
- 238000004364 calculation method Methods 0.000 claims description 23
- 230000008569 process Effects 0.000 claims description 16
- 230000005684 electric field Effects 0.000 claims description 15
- 238000005457 optimization Methods 0.000 claims description 9
- 230000004044 response Effects 0.000 claims description 9
- 238000005259 measurement Methods 0.000 claims description 7
- 239000013598 vector Substances 0.000 claims description 7
- 230000010354 integration Effects 0.000 claims description 3
- 230000035699 permeability Effects 0.000 claims description 3
- 238000005303 weighing Methods 0.000 claims description 3
- 238000009499 grossing Methods 0.000 claims description 2
- 230000002194 synthesizing effect Effects 0.000 claims description 2
- 230000005389 magnetism Effects 0.000 claims 1
- 238000012876 topography Methods 0.000 abstract description 5
- 238000004422 calculation algorithm Methods 0.000 description 9
- 230000002159 abnormal effect Effects 0.000 description 5
- 238000005516 engineering process Methods 0.000 description 4
- 238000002939 conjugate gradient method Methods 0.000 description 2
- 238000010276 construction Methods 0.000 description 2
- 238000009795 derivation Methods 0.000 description 2
- 238000001514 detection method Methods 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 238000003384 imaging method Methods 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 238000004088 simulation Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 230000005856 abnormality Effects 0.000 description 1
- 238000007792 addition Methods 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000001816 cooling Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000009826 distribution Methods 0.000 description 1
- PCHJSUWPFVWCPO-UHFFFAOYSA-N gold Chemical compound [Au] PCHJSUWPFVWCPO-UHFFFAOYSA-N 0.000 description 1
- 239000010931 gold Substances 0.000 description 1
- 229910052737 gold Inorganic materials 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
- 238000003860 storage Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V3/00—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
- G01V3/38—Processing data, e.g. for analysis, for interpretation, for correction
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V3/00—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
- G01V3/15—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation specially adapted for use during transport, e.g. by a person, vehicle or boat
- G01V3/16—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation specially adapted for use during transport, e.g. by a person, vehicle or boat specially adapted for use from aircraft
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V3/00—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
- G01V3/15—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation specially adapted for use during transport, e.g. by a person, vehicle or boat
- G01V3/165—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation specially adapted for use during transport, e.g. by a person, vehicle or boat operating with magnetic or electric fields produced or modified by the object or by the detecting device
-
- 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
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- Remote Sensing (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Aviation & Aerospace Engineering (AREA)
- Electromagnetism (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
The invention discloses a multi-component frequency domain based aviation electromagnetic data fusion three-dimensional inversion method, which comprises the steps of firstly, developing frequency domain aviation electromagnetic finite element three-dimensional forward modeling under any complex structure and fluctuating topography conditions based on flexible modeling characteristics of tetrahedron grids, and providing effective precondition guarantee for refined three-dimensional inversion; secondly, based on a quasi-Newton method of a limited memory, the approximation of a Heisen matrix is realized, and the aviation electromagnetic inversion speed is effectively accelerated; the comparison of the single data inversion result and the multiple data fusion inversion results shows that the inversion method is more accurate, more underground structure information can be obtained, and a new thought is provided for the fine structural interpretation of aviation electromagnetic data based on a frequency domain.
Description
Technical Field
The invention belongs to the technical field of geophysical signal processing, and particularly relates to a multi-component frequency domain based aviation electromagnetic data fusion three-dimensional inversion method.
Background
The resource exploration needs to carry out secondary fine detection on army and near-surface areas in areas with complex geological conditions. The aviation electromagnetic method takes the aircraft as a carrier, has the advantages of high exploration speed, low consumption cost and the like, can rapidly realize regional electromagnetic data scanning acquisition work, and becomes one of the preferred technical methods for geological condition complex regional exploration work. Meanwhile, with the continuous upgrading of an aviation electromagnetic hardware system, the transverse resolution and the sampling rate of the aviation electromagnetic hardware system are greatly improved, and the near-surface fine detection is provided with necessary guarantee. However, as the structure of the exploration target becomes complex and refined, this presents a serious challenge for massive aviation electromagnetic data processing at high sampling rates.
The interpretation of aviation electromagnetic actual measurement data mainly ensures that the imaging technology and one-dimensional inversion of mass data are mainly processed, and the one-dimensional imaging and inversion are based on the assumption of a layered or ground model, so that the precise inversion of a complex ground model is not enough. Therefore, the three-dimensional aviation electromagnetic inversion capable of recovering the fine and complex structure of the ground model is a key technical means for solving the explanation of aviation electromagnetic mass data, and is also a key technical problem for advancing complex geological conditions and near-surface aviation electromagnetic investigation.
In addition, most frequency domain electromagnetic data acquisition systems take a helicopter nacelle as a platform and simultaneously carry two groups of vertical coaxial and horizontal coplanar coil pairs to acquire two groups of orthogonal components, and in data interpretation, the data inversion is generally only based on each group of independent coil pairs, so that the requirement of fine structure interpretation is hardly met, and therefore, the development of fusion inversion of the data of the two groups of coil pairs is imperative.
Forward is an inversion core technology part. Among the forward algorithms, the finite difference algorithm is simpler and is convenient to program, and is most used in the current three-dimensional forward algorithm, and then the integral equation method and the finite volume method and finally the finite element method are adopted. The finite element method can adapt to various irregular grids such as tetrahedrons, deformation hexahedrons and the like, but the finite element method has larger memory consumption in calculation and limits the application of the finite element method in early three-dimensional forward modeling. With the development of computer hardware, finite element methods based on unstructured grids have become a hotspot of current research. Currently, three-dimensional inversion algorithms are relatively mature and mainly include nonlinear conjugate gradient methods, gaussian-Newton methods, quasi-Newton methods and finite memory quasi-Newton methods. The nonlinear conjugate gradient method avoids solving the hessian matrix, greatly reduces the calculated amount, and has slower convergence speed. The Gaussian-Newton method ignores the second order term of the Hessen matrix, and improves the convergence rate. The quasi-Newton method adopts an iteration method to obtain an approximate Hessen matrix in inversion, and compared with the traditional quasi-Newton method, the finite memory quasi-Newton method reduces the calculation memory more effectively and is suitable for massive electromagnetic data inversion. The current three-dimensional inversion method mainly adopts a limited difference method based on structural grids, however, for the region with strong topography fluctuation, the air and ground interface needs to be split sufficiently finely, the number of grids is large, the memory consumption is huge, the calculation time is long, and the application of the three-dimensional inversion method in complex topography and construction regions is limited.
Moreover, the existing algorithm such as the finite difference method and the finite volume of the regular grid cannot fit the terrain, the calculation error is large, and meanwhile, the existing inversion technology is single inversion and is strong in non-uniqueness. Two kinds of data cannot be utilized at the same time.
Disclosure of Invention
In view of the above, the invention aims to provide a three-dimensional inversion method based on multi-component frequency domain aviation electromagnetic data fusion, which adopts a finite element method and is matched with a non-structural tetrahedral grid, so that the problem of efficiently fitting terrains can be solved, and the calculation accuracy is high; meanwhile, the HCP and VCX measured data are adopted for fusion inversion, so that the non-uniqueness of three-dimensional inversion can be effectively reduced, and more underground structure information can be obtained.
In order to achieve the above object, an embodiment of the present invention provides a multi-component frequency domain based aviation electromagnetic data fusion three-dimensional inversion method, which includes the following steps:
step 1, constructing an initial model to be inverted, and carrying out unstructured tetrahedral mesh division on the initial model;
step 2, finite element forward modeling based on unstructured tetrahedral meshes is carried out on the initial model;
step 3, acquiring actually-measured multi-component frequency domain aviation electromagnetic data, including: the device comprises actual measurement data collected by an aviation electromagnetic HCP device and actual measurement data collected by an aviation electromagnetic VCX device;
step 4, constructing a data fitting item according to the multi-component frequency domain aviation electromagnetic data, constructing a model constraint item according to the model prior information, and constructing an inversion objective function based on the data fitting item and the model constraint item;
and 5, performing optimization calculation on the inversion objective function by using a quasi-Newton method of a limited memory, adding partial derivatives of the data fitting term into a finite element forward equation solution in the optimization process, and updating the model through iterative updating solution to realize a final inversion model.
In step 2 of one embodiment, when performing finite element forward modeling based on unstructured tetrahedral mesh on an initial model in a frequency domain aviation electromagnetic system, a frequency domain second-order field electric field satisfies the following helmholtz equation:
wherein E is s Is a secondary electric field, E p Is the primary electric field, ω is the angular velocity, μ is the magnetic permeability, σ is the electrical conductivity, σ p Is background field conductivity;
to ensure uniqueness of the quadratic field solution, at the boundaryThe dirichlet first class boundary conditions are used:
the finite element forward equation of the frequency domain electromagnetic method can be obtained by multiplying the formula (2) with the vector basis function N point and performing step-by-step integration by utilizing the first green formula to obtain a coefficient matrix equation of each tetrahedron grid unit, and synthesizing the coefficient matrices of all the tetrahedron grid units into an overall matrix K:
KE s ＝S (3)
s is an endogenous term of a tetrahedral grid unit, and a secondary magnetic field H of frequency aviation electromagnetic response is obtained by solving a finite element forward equation s 。
Step 4 of one embodiment, a data fitting term constructed from measured data, expressed as:
wherein phi is d (m) represents a data fitting term, N is the number of data,the table isActually measured multi-component frequency domain aviation electromagnetic data, d n A secondary magnetic field H determined for forward modeling s E, e n And (3) data errors, wherein m is a geological attribute of the model, namely a variable obtained by inversion.
Step 4 of one embodiment, constructing a model constraint term Φ according to the model prior information m (m) expressed as:
Ф m (m)＝α r ‖Rm‖ 2 +α s ‖W s (m-m pri )‖ 2 (5)
wherein alpha is r And alpha s Respectively weighing, R is a model smoothing matrix, m is a geological attribute of a model, and m pri For a priori geologic attribute as model a priori information, W s A constraint matrix is used for model priori;
inversion objective function Φ, expressed as:
Ф＝Ф d (m)+λФ m (m)
＝Ф d (m)+λ(α r ‖Rm‖ 2 +α s ‖W s (m-m pri )‖ 2 ) (6)
Φ d (m) a data fitting term.
In step 5 of one embodiment, when the solution calculation is performed on the inversion objective function by using the quasi-newton method with limited memory, in order to calculate the gradient of the data fitting term in the inversion objective function, the actually measured multi-component frequency domain aviation electromagnetic data is obtainedAnd a secondary magnetic field H obtained by forward modeling s Of (d) is a magnetic field component d n And respectively separating and arranging according to a real part and an imaginary part to obtain:
where N is the data index, N is the data number, Δd n Is thatAnd d n A difference between;
whereas the gradient of the data fitting term is expressed as:
for the kth tetrahedral grid cell there is:
wherein, is complex conjugate, M is the total number of tetrahedral grid units,an operator representing the complex number as a real part;
these partial derivatives of equation (9) are applied to the finite element forward equation solution E, which is the electric field vector, and the initial model is updated by iterative update solution.
In step 5 of one embodiment, in the process of solving the inversion objective function by using the quasi-newton method of the finite memory, after the partial derivatives of the data fitting term are added to the finite element forward equation solution, forward calculation is performed at each frequency, where one forward calculation is to calculate the quadratic electric field E in the finite element forward equation s And a secondary magnetic field H of frequency aviation electromagnetic response s The method comprises the steps of carrying out a first treatment on the surface of the Another forward calculation is to accompany the aviation electromagnetic response secondary magnetic field H with different frequencies s The accompanying forward computation of the components to solve the overall matrix K.
In step 5 of one embodiment, in the process of solving the inversion objective function by using the quasi-newton method with limited memory, the approximated matrix constructed by the quasi-newton method is used to replace the hessian matrix, so as to achieve approximation of the hessian matrix and accelerate the aviation electromagnetic inversion speed.
In step 5 of one embodiment, in the process of solving the inversion objective function by using the quasi-newton method with limited memory, the model constraint term is directly calculated according to the model prior information and the updated model obtained by the last solution.
Compared with the prior art, the invention has the beneficial effects that at least the following steps are included:
in the technical scheme provided by the embodiment, the fusion three-dimensional inversion of the multi-component aviation electromagnetic data is carried out based on a finite element method of an unstructured tetrahedron grid and a quasi-Newton method of a finite memory, firstly, the frequency domain aviation electromagnetic finite element three-dimensional forward modeling under the conditions of any complex structure and relief topography can be carried out based on the flexible modeling characteristics of the tetrahedron grid, and an effective precondition guarantee is provided for the refined three-dimensional inversion; secondly, based on a quasi-Newton method of a limited memory, the approximation of a Heisen matrix is realized, and the aviation electromagnetic inversion speed is effectively accelerated; the inversion method based on the multi-component frequency domain aviation electromagnetic data fusion is more accurate in inversion and can obtain more underground structure information through comparison of a single data inversion result and a plurality of data fusion inversion results, and a new thought is provided for fine structural interpretation based on the frequency domain aviation electromagnetic data.
Drawings
In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings that are required in the embodiments or the description of the prior art will be briefly described, and it is obvious that the drawings in the following description are only some embodiments of the present invention, and other drawings may be obtained according to these drawings without inventive effort for a person skilled in the art.
FIG. 1 is a flow chart of a multi-component frequency domain based aviation electromagnetic data fusion three-dimensional inversion method provided by an embodiment;
FIG. 2 is a model diagram of an experimental example provided in the examples;
FIG. 3 is an inversion of measured data collected for a HCP device, provided in the examples;
FIG. 4 is an inversion of measured data acquired for a VCX device provided by an embodiment;
FIG. 5 is an inversion of two data provided by the example for the measured data collected by the aero-electromagnetic HCP device and the measured data collected by the aero-electromagnetic VCX device;
FIG. 6 is an inversion convergence graph of measured data collected for a HCP device, provided in an embodiment;
FIG. 7 is an inversion convergence graph of measured data acquired for a VCX device provided in an embodiment;
fig. 8 is an inversion convergence graph of two data, namely, measured data collected for an aero-electromagnetic HCP device and measured data collected for an aero-electromagnetic VCX device, provided in an embodiment.
Detailed Description
The present invention will be described in further detail with reference to the drawings and examples, in order to make the objects, technical solutions and advantages of the present invention more apparent. It should be understood that the detailed description is presented by way of example only and is not intended to limit the scope of the invention.
Fig. 1 is a flowchart of an embodiment of a multi-component frequency domain based aviation electromagnetic data fusion three-dimensional inversion method. As shown in fig. 1. The method for fusion three-dimensional inversion of aviation electromagnetic data based on the multi-component frequency domain provided by the embodiment comprises the following steps:
step 1, constructing an initial model to be inverted, and carrying out unstructured tetrahedral mesh division on the initial model.
In an embodiment, an initial model is constructed according to a geologic structure of an inversion test, and geologic properties of the initial model including magnetic field information such as conductivity, resistivity and the like are configured. In view of flexible modeling characteristics based on tetrahedral grids, and in view of structural complexity and topography fluctuation of a geological structure, an initial model is divided by adopting unstructured tetrahedral grids, so that frequency domain aviation electromagnetic finite element three-dimensional forward modeling is facilitated.
And 2, performing finite element forward modeling based on the unstructured tetrahedral mesh on the initial model.
In an embodiment, in a frequency domain avionics system, the initial model is subjected to finite element forward modelingThe total field over the secondary field ratio is measured. And it is described that the quadratic field has a smaller influence range (footprint) than the total field, and therefore the quadratic field equation is solved in a finite element forward process. Assume a positive time harmonic e iωt Then the frequency domain second field electric field satisfies the following helmholtz equation:
wherein E is s Is a secondary electric field, E p Is the primary electric field, ω is the angular velocity, μ is the magnetic permeability, σ is the electrical conductivity, σ p Is background field conductivity. To ensure uniqueness of the quadratic field solution, here at the boundaryThe dirichlet first class boundary conditions are used:
(1) The formula and the formula (2) together form a forward theory boundary value problem. (1) The sum vector basis function N is multiplied by the point, and the first green formula is utilized to perform step integration, thus obtaining
dV is the rightmost term of the integral operator and represents any tiny volume within the regional omega. Decomposing the calculation region of the formula (3) into a series of tetrahedral grid units, and calculating the integral of the formula (3) in each tetrahedral grid unit to obtain a small coefficient matrix equation
Where i is the index of the tetrahedral network element, K e To be composed of density matrix and rigidityCoefficient matrix formed by sum of degree matrix S e Is an endogenous term of the tetrahedral network element,the coefficient matrix of all tetrahedral network units is synthesized into a total matrix K for the secondary electric field of the tetrahedral network units, and the finite element forward equation of the frequency domain electromagnetic method can be obtained
KΕ s ＝S (5)
S is a source term matrix, and equation (5) is solved by a multi-core parallel direct solver. Thereafter, a secondary magnetic field H of the frequency aeroelectromagnetic response s Can be obtained by faraday's theorem.
Step 3, acquiring actually-measured multi-component frequency domain aviation electromagnetic data, including: the device comprises actual measurement data collected by an aviation electromagnetic HCP device and actual measurement data collected by an aviation electromagnetic VCX device.
In frequency domain aero-electromagnetic data simulation, data of a fixed aero-electromagnetic HCP device or a VCX device are processed, wherein the aero-electromagnetic HCP device refers to a transmitting coil and a receiving coil which are both vertical magnetic dipoles, and the aero-electromagnetic VCX device refers to a transmitting coil and a receiving coil which are both horizontal magnetic dipoles. In an embodiment, the measured data collected by the aviation electromagnetic HCP device and the measured data collected by the aviation electromagnetic VCX device for the same geological structure are obtained simultaneously to form the measured multi-component frequency domain aviation electromagnetic data, and the measured multi-component frequency domain aviation electromagnetic data are used for inversion calculation of the geological structure.
And 4, constructing a data fitting item according to the multi-component frequency domain aviation electromagnetic data, constructing a model constraint item according to the model prior information, and constructing an inversion objective function based on the data fitting item and the model constraint item.
In an embodiment, an inversion objective function Φ of a nonlinear frequency domain three-dimensional aviation electromagnetic method constructed based on a data fitting term and a model constraint term is expressed as:
Ф＝Ф d (m)+λФ m (m) (6)
wherein lambda is the equilibrium data fitting term phi d (m) and modelConstraint term phi m And (m) a real number, called regularization factor. The three-dimensional aeroelectromagnetic inversion problem is equivalent to minimizing the objective function problem described above. Reasonable regularization factors can ensure phi m (m) and phi d (m) no excessive fitting in Φ. According to equation (6), the gradient of the objective function can be expressed as:
the calculation of the three-dimensional aviation electromagnetic method objective function and the gradient thereof is described in detail based on the unstructured finite element method. Phi (phi) d (m) can be expressed as:
wherein N is the number of data,the table is actually measured multi-component frequency domain aviation electromagnetic data, d n A secondary magnetic field H determined for forward modeling s Of (a), i.e. magnetic field component H x 、H y And H z ，e n And (3) data errors, wherein m is a geological attribute of the model, namely a variable obtained by inversion.
Model constraint term phi m (m) can be expressed as:
Φ m (m)＝α r Φ r (m)+α s Φ s (m) (9)
wherein phi is m (m)＝||Rm|| 2 (10)
Φ s (m)＝||W s (m-m pri )|| 2 (11)
Wherein alpha is r And alpha s Respectively weighing, R is a model constraint matrix, m is a geological attribute of a model, and m pri For a priori geologic attribute as model a priori information, W s A constraint matrix is used for model priori;
the objective function may be rewritten as:
Φ(m)＝Φ d (m)+λ(α r Φ r (m)+α s Φ s (m)) (12)
for unstructured tetrahedral meshes, the differential matrix R in equation (10) does not need to be obtained by display, and W can be obtained by solving the differences between tetrahedral mesh cells and adjacent cells of the tetrahedral mesh one by one s Can be expressed by the inverse of the volume of the current model, then there are:
wherein V is i For the cell volume corresponding to the ith tetrahedral mesh cell, N (i) The number of adjacent cells of the ith tetrahedral mesh unit is M, the total number of the tetrahedral mesh units is d ij Is the distance between the centers of adjacent units, w j For the weight coefficient of the j-th adjacent unit, deltam ij Model parameter differences, w, of adjacent cells j And Δm ij Is defined as:
let the ith grid cell center point beThen:
in summary, the gradient of the objective function can be expressed as:
and 5, performing optimization calculation on the inversion objective function by using a Newton method with limited memory so as to update the model and realize a final inversion model.
In an embodiment, in order to calculate the gradient of the data fitting term in the inversion objective function, the measured multi-component frequency domain aviation electromagnetic secondary magnetic field dataAnd a secondary magnetic field H obtained by forward modeling s Of (d) is a magnetic field component d n The method comprises the steps of respectively separating and arranging according to a real part and an imaginary part, wherein the first N data are the real parts of magnetic field components, and the last N data are the imaginary parts of the magnetic field components, so that the method comprises the following steps of:
where N is a data index, N is a data number,representing plural forms->d n ' represents d in complex form n ，Δd n Is->And d n The difference, and the gradient of the data fitting term, can be expressed as:
thus, for the ith tetrahedral unit there is:
wherein, is complex conjugate.
In an embodiment, these partial derivatives, as shown in equation (20) and equation (21), may be applied to the forward equation solution E during the optimization process. Is provided withFor a given frequency and the magnetic field component H of the j-site magnetic field interpolated from the electric field vector E at the source x,j Then:
similarly, it can be defined thatSimilar to formula (21), there are:
forward equation (5) of aeroelectromagnetic method about m k The derivation can be obtained:
from the above derivation, equation (24) can be rewritten as:
wherein,,dependent on different aeroelectromagnetic response components, e.g. for H x Corresponding->Is->Defining the simulation of all sources at a single frequency is called forward, equation (25) requires two forward calculations at each frequency. Wherein one forward is calculated E and response H; in addition, the right end item g is required to be twice T To solve the overall matrix K, wherein the relatively stable coefficients g T The method comprises the following steps:
thus, let v be T The method meets the following conditions:
v T ＝g T K -1 (27)
thus, the companion forward can be defined:
Kv＝g (28)
the complete calculation of equation (25) can be completed by equation (28).
Comparison data fitting term Φ d Model constraint term Φ m Can be directly resolved to obtain:
and directly calculating a model constraint term according to the model prior information and an updated model obtained by the last solving.
For the objective function phi with M unknowns, in the inversion of a multi-component aeroelectromagnetic method, the quantity can reach tens of thousands to millions, and in order to accelerate convergence, the method adoptsThe inversion process is suitable for the inversion process of the quasi-Newton optimization algorithm based on the limited memory of the three-dimensional aviation electromagnetic data, and the algorithm is developed from the quasi-Newton method. The quasi-Newton method adopts an approximate matrix B in the k+1th iteration process of the model k Instead of the hessian matrix H k According to the nature of the Taylor series, the objective function Φ (m) is at m k+1 The quadratic approximation at this point can be expressed as:
the above method is to find the derivative:
g(m)≈g k+1 +H k+1 (m-m k+1 )，·····················(32)
let m=m k ，s k ＝m k+1 -m k ，y k ＝g k+1 -g k The following steps are:
H k+1 s k ≈y k ，····················(33)
the quasi-newton method requires construction:
B k+1 s k ＝y k .··························(34)
the calculated relation in the iterative process of the quasi-Newton method can be obtained through further arrangement, namely
The quasi-Newton method with limited memory can further save memory and improve efficiency based on Newton method, and the transformation of (35) by using Sherman-Morrison formula can be obtained:
by F k+1 Representation ofLet->I is an identity matrix. The above simplification is:
taking into account the approximate hessian inverse matrix in the previous several adjacent iterations, an initial F is given first 0 The most recent m correction relations are preserved in the limited memory quasi-newton method, and the following can be obtained:
it can be seen that the Newton method is a finite memory method, and only the pair of column vectors { s }, is stored k Sum { y } k No storage of F k+1 。
The specific flow can be summarized as follows:
1) Setting inversion termination conditions, setting the precision epsilon to be more than 0 and setting a target root mean square value r 0 Let k=1;
2) Given an initial positive diagonal matrixWherein the proportion factor->γ k ＝1；
3) Calculation ofAnd RMS, if->Or RMS.ltoreq.r 0 The iteration inversion is terminated, and m is output k ；
4) Let n=min { k, m-1}, pairUpdating n+1 times to obtain F k+1 Is an expression of (2);
5) Determination of Δm using double-loop recursion k ＝-F k g k ；
6) According to Wolfe condition, setting 0 < beta '< 0.5, beta' < beta < 1, searching step alpha in linear mode k ：
Let m k+1 ＝m k +α k d k ；
7) Let k=k+1, return to step 2).
According to the process, the three-dimensional inversion of the multi-component frequency domain aviation electromagnetic data is realized by combining a limited memory quasi-Newton optimization algorithm.
According to the multi-component frequency domain-based aviation electromagnetic data fusion three-dimensional inversion method provided by the embodiment, the advantages of two devices, namely a HCP device and a VCA device, are fully explored, two data fusion inversion technologies are developed, in the forward and accompanying forward processes, two device modes are calculated respectively, sensitivity is solved, optimization calculation is carried out by using a Newton method with limited memory, and two data are orderly arranged for fusion inversion.
Experimental example
In order to test the effectiveness of the three-dimensional inversion method based on the multi-component frequency domain aviation electromagnetic data fusion provided by the embodiment, the experimental example simultaneously carries out forward modeling on the measured data of the HCP and VCX devices, the model is set to have a low-resistance abnormal body in a background space, and 5% Gaussian noise is added into the forward modeling magnetic field data to simulate the measured data. Fig. 2 shows the anomaly and 225 point distributions of the present model, wherein (a) is a top view of the model and (b) is a transverse slice of the model, and the distances between the points and the lateral lines are 20m. The resistivity of the background half space is 100 omega m, the abnormal body is positioned at 20m under the center of the measuring point, and the abnormal body is positioned at the position of the measuring pointThe size of the body was 100m×100m×25m, and the abnormal resistivity was 10Ω·m. The grid used by forward modeling of the theoretical model contained 431223 units, and the HCP and VCX responses were calculated at 2 frequencies, 900Hz and 5000Hz, respectively. The mesh used for model inversion contains 214325 cells. The invention uses the same cooling principle to change lambda. The initial lambda is chosen to be 0.01, the additional regularization factor alpha r And alpha s All were set to 0.1. The prior model and the reference model are set to a half space of 100 Ω·m. To illustrate the superiority of the present invention combining HCP and VCX data, separate HCP and VCX inversion was set up as a comparison.
Experimental examples were performed on an Intel (R) Xeon (R) Gold 6254 CPU 3.1ghz workstation. FIG. 3 shows an HCP inversion alone, where (a) is a y-z section plot (x=0m); (b) is an x-z slice (y=0m); (c) For the x-y slice (z= -30 m), the box area represents the actual outlier location. Analyzing FIG. 3, it can be seen that HCP inversion can generally restore true horizontal positions of anomalies. However, there are two disadvantages, first: the inverted slice shows that the depth of the abnormal body is difficult to reach the real depth, and only the abnormality of the shallow part is effectively recovered; second,: the profile of the horizontal plane of the anomaly is deformed, approximating a circle less than the true range your. Two disadvantages of HCP inversion will directly affect the accuracy of the data inversion.
FIG. 4 shows VCX inversion alone, where (a) is a y-z section (x=0m); (b) is an x-z slice (y=0m); (c) For the x-y slice (z= -30 m), the box area represents the actual outlier location. Analysis of fig. 4 shows that the VCX inversion depth is deeper than the HCP and that the profile of the anomaly is more nearly true. Two disadvantages of single HCP inversion are overcome, but VCX inversion slices have some small high-resistance false anomalies in the shallow portion, which can introduce errors into the data inversion interpretation.
FIG. 5 is a joint inversion result of a fusion HCP and VCX, wherein (a) is a y-z section plot (x=0m); (b) is an x-z slice (y=0m); (c) For the x-y slice (z= -30 m), the box area represents the actual outlier location. Analyzing fig. 5, it can be found that the joint inversion is superior to the single HCP inversion in both the inversion depth and the horizontal profile, and overcomes the disadvantages of the single HCP inversion, and meanwhile, the joint inversion also has no false high-resistance anomalies in the shallow part, which illustrates that the joint inversion can also overcome the disadvantages of the false anomalies of the single VCX inversion.
In order to further illustrate the effectiveness of the multi-component frequency domain based aviation electromagnetic data fusion three-dimensional inversion method provided by the embodiment. FIG. 6 shows an HCP inversion convergence plot wherein (a) shows the data fitting difference Φ of the inversion of the HCP calculation over 112 iterations d From 100.1 down to 2.6, the roughness Φ of the model is shown r And phi is s Finally, stable convergence is achieved; (b) The HCP inversion is shown with 12 steps a that require changing the inversion and (c) the HCP inversion is shown with 6 changes to the regularization factor λ.
FIG. 7 shows a HCP inversion convergence plot in which (a) shows the data fitting difference Φ for an inversion of the VCX algorithm over 66 iterations d From 110.2 down to 24.6, the roughness Φ of the model r And phi is s Finally, stable convergence is achieved; (b) The HCP inversion is shown with a step size α of 14 times that the inversion needs to be changed; (c) shows that the VCX inversion has 7 changes in regularization factor λ.
FIG. 8 shows a HCP and VCX joint inversion convergence plot in which (a) shows the data fitting difference Φ for an example over 48 iterations d From 110.2 down to 26.1, the roughness Φ of the model r And phi is s Finally, stable convergence is achieved; (b) The HCP and VCX joint inversion is shown with 5 steps α that require changing the inversion, and (c) the HCP and VCX inversion is shown with 5 changes in regularization factor. First, HCP inversion converges most slowly in terms of iteration number, 112 iterations, 68 VCX inversions, and joint inversion is fastest, requiring only 48 times. Secondly, the number of step changes in the inversion iteration process is up to 14, 12, and 5 for the combined HCP and VCX inversion, because if the step length is not changed, each iteration needs one forward, the number of required forward is greater than 1, and the effectiveness of the combined inversion is also demonstrated from the number of step changes. Finally, the VCX changes the most, the HCP the least, the joint inversion changes in terms of the number of changes in the regularization factorThe smaller number of changes represents better inversion stability, and from this point, the effectiveness of the multi-component frequency domain-based aviation electromagnetic data fusion three-dimensional inversion method is also demonstrated.
The foregoing detailed description of the preferred embodiments and advantages of the invention will be appreciated that the foregoing description is merely illustrative of the presently preferred embodiments of the invention, and that no changes, additions, substitutions and equivalents of those embodiments are intended to be included within the scope of the invention.
Claims (7)
1. A three-dimensional inversion method based on multi-component frequency domain aviation electromagnetic data fusion is characterized in that,
step 1, constructing an initial model to be inverted, and carrying out unstructured tetrahedral mesh division on the initial model;
step 2, finite element forward modeling based on the unstructured tetrahedral mesh is performed on the initial model, specifically, when finite element forward modeling based on the unstructured tetrahedral mesh is performed on the initial model in a frequency domain aviation electromagnetic system, the following Helmholtz equation is satisfied by a frequency domain secondary field electric field:
wherein E is s Is a secondary electric field, E p Is the primary electric field, ω is the angular velocity, μ is the magnetic permeability, σ is the electrical conductivity, σ p Is background field conductivity;
to ensure uniqueness of the quadratic field solution, at the boundaryThe dirichlet first class boundary conditions are used:
the finite element forward equation of the frequency domain electromagnetic method can be obtained by multiplying the formula (2) with the vector basis function N point and performing step-by-step integration by utilizing the first green formula to obtain a coefficient matrix equation of each tetrahedron grid unit, and synthesizing the coefficient matrices of all the tetrahedron grid units into an overall matrix K:
KE s ＝S (3)
s is an endogenous term of a tetrahedral grid unit, and a secondary magnetic field H of frequency aviation electromagnetic response is obtained by solving a finite element forward equation s ；
Step 3, acquiring actually-measured multi-component frequency domain aviation electromagnetic data, including: the device comprises actual measurement data collected by an aviation electromagnetic HCP device and actual measurement data collected by an aviation electromagnetic VCX device;
step 4, constructing a data fitting item according to the multi-component frequency domain aviation electromagnetic data, constructing a model constraint item according to the model prior information, and constructing an inversion objective function based on the data fitting item and the model constraint item;
and 5, performing optimization calculation on the inversion objective function by using a quasi-Newton method of a limited memory, adding partial derivatives of the data fitting term into a finite element forward equation solution in the optimization process, and updating the model through iterative updating solution to realize a final inversion model.
2. The method for three-dimensional inversion based on multicomponent frequency domain aviation electromagnetic data fusion according to claim 1, wherein step 4, a data fitting term constructed according to measured data is expressed as:
wherein phi is d (m) represents a data fitting term, N is the number of data,the table is actually measured multi-component frequency domain aviation electromagnetic data, d n Secondary magnetism obtained for forward modelingField H s E, e n And (3) data errors, wherein m is a geological attribute of the model, namely a variable obtained by inversion.
3. The method for fusion three-dimensional inversion based on multi-component frequency domain aviation electromagnetic data according to claim 1 or 2, wherein step 4, a model constraint term Φ is constructed according to model prior information m (m) expressed as:
Ф m (m)＝α r ‖Rm‖ 2 +α s ‖W s (m-m pri )‖ 2 (5)
wherein alpha is r And alpha s Respectively weighing, R is a model smoothing matrix, m is a geological attribute of a model, and m pri For a priori geologic attribute as model a priori information, W s A constraint matrix is used for model priori;
inversion objective function Φ, expressed as:
Ф＝Ф d (m)+λФ m (m)
＝Ф d (m)+λ(α r ‖Rm‖ 2 +α s ‖W s (m-m pri )‖ 2 ) (6)
Φ d (m) a data fitting term.
4. The three-dimensional inversion method based on multi-component frequency domain aviation electromagnetic data fusion according to claim 1, wherein in step 5, when the inversion objective function is solved and calculated by utilizing the quasi-newton method of finite memory, in order to calculate the gradient of the data fitting term in the inversion objective function, the actually measured multi-component frequency domain aviation electromagnetic data is obtainedAnd a secondary magnetic field H obtained by forward modeling s Of (d) is a magnetic field component d n And respectively separating and arranging according to a real part and an imaginary part to obtain:
where N is the data index, N is the data number, Δd n Is thatAnd d n A difference between;
whereas the gradient of the data fitting term is expressed as:
for the kth tetrahedral grid cell there is:
wherein, is complex conjugate, M is the total number of tetrahedral grid units,an operator representing the complex number as a real part;
these partial derivatives of equation (9) are applied to the finite element forward equation solution E, which is the electric field vector, and the initial model is updated by iterative update solution.
5. The three-dimensional inversion method based on multi-component frequency domain aviation electromagnetic data fusion according to claim 1, wherein in step 5, in the process of solving the inversion objective function by using the quasi-newton method of finite memory, after the partial derivatives of the data fitting term are added to the finite element forward equation solution, forward calculation is performed at each frequency, wherein one forward calculation is to calculate the quadratic electric field E in the finite element forward equation s Secondary magnetic field for sum frequency aero-electromagnetic responseH s The method comprises the steps of carrying out a first treatment on the surface of the Another forward calculation is to accompany the aviation electromagnetic response secondary magnetic field H with different frequencies s The accompanying forward computation of the components to solve the overall matrix K.
6. The three-dimensional inversion method based on multi-component frequency domain aviation electromagnetic data fusion according to claim 1 is characterized in that in step 5, in the process of solving an inversion objective function by using a quasi-Newton method with limited memory, an approximation matrix constructed by the quasi-Newton method is adopted to replace a Hessen matrix, so that approximation of the Hessen matrix is realized, and aviation electromagnetic inversion speed is accelerated.
7. The three-dimensional inversion method based on multi-component frequency domain aviation electromagnetic data fusion according to claim 1, wherein in step 5, in the process of solving an inversion objective function by using a quasi-Newton method with limited memory, model constraint terms are directly calculated according to model prior information and an updated model obtained by the last solving.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210799540.5A CN115113286B (en) | 2022-07-06 | 2022-07-06 | Aviation electromagnetic data fusion three-dimensional inversion method based on multi-component frequency domain |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210799540.5A CN115113286B (en) | 2022-07-06 | 2022-07-06 | Aviation electromagnetic data fusion three-dimensional inversion method based on multi-component frequency domain |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN115113286A CN115113286A (en) | 2022-09-27 |
| CN115113286B true CN115113286B (en) | 2023-09-15 |
Family
ID=83332678
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210799540.5A Active CN115113286B (en) | 2022-07-06 | 2022-07-06 | Aviation electromagnetic data fusion three-dimensional inversion method based on multi-component frequency domain |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115113286B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN118244370B (en) * | 2024-05-28 | 2024-07-23 | 吉林大学 | Three-dimensional inversion method for transient electromagnetic tunnel advanced detection based on finite element algorithm |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106199742A (en) * | 2016-06-29 | 2016-12-07 | 吉林大学 | A kind of Frequency-domain AEM 2.5 ties up band landform inversion method |
| CN110058317A (en) * | 2019-05-10 | 2019-07-26 | 成都理工大学 | Aviation transient electromagnetic data and aviation magnetotelluric data joint inversion method |
| CN111474594A (en) * | 2020-05-27 | 2020-07-31 | 长安大学 | Three-dimensional time domain aviation electromagnetic fast inversion method |
| CN112949134A (en) * | 2021-03-09 | 2021-06-11 | 吉林大学 | Earth-well transient electromagnetic inversion method based on non-structural finite element method |
| CN114460654A (en) * | 2022-02-22 | 2022-05-10 | 成都理工大学 | Semi-aviation transient electromagnetic data inversion method and device based on L1L2 mixed norm |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20090083006A1 (en) * | 2007-09-20 | 2009-03-26 | Randall Mackie | Methods and apparatus for three-dimensional inversion of electromagnetic data |
| US10838093B2 (en) * | 2015-07-02 | 2020-11-17 | Exxonmobil Upstream Research Company | Krylov-space-based quasi-newton preconditioner for full-wavefield inversion |
-
2022
- 2022-07-06 CN CN202210799540.5A patent/CN115113286B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106199742A (en) * | 2016-06-29 | 2016-12-07 | 吉林大学 | A kind of Frequency-domain AEM 2.5 ties up band landform inversion method |
| CN110058317A (en) * | 2019-05-10 | 2019-07-26 | 成都理工大学 | Aviation transient electromagnetic data and aviation magnetotelluric data joint inversion method |
| CN111474594A (en) * | 2020-05-27 | 2020-07-31 | 长安大学 | Three-dimensional time domain aviation electromagnetic fast inversion method |
| CN112949134A (en) * | 2021-03-09 | 2021-06-11 | 吉林大学 | Earth-well transient electromagnetic inversion method based on non-structural finite element method |
| CN114460654A (en) * | 2022-02-22 | 2022-05-10 | 成都理工大学 | Semi-aviation transient electromagnetic data inversion method and device based on L1L2 mixed norm |
Also Published As
| Publication number | Publication date |
|---|---|
| CN115113286A (en) | 2022-09-27 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Ford et al. | A nonhydrostatic finite-element model for three-dimensional stratified oceanic flows. Part I: Model formulation | |
| CN106199742B (en) | A kind of dimension of Frequency-domain AEM 2.5 band landform inversion method | |
| CN112949134B (en) | Earth-well transient electromagnetic inversion method based on non-structural finite element method | |
| CN106646645B (en) | A kind of gravity forward modeling accelerated method | |
| CN109946518B (en) | Power harmonic signal analysis method and analysis equipment based on Bayes method | |
| CN110346835B (en) | Magnetotelluric forward modeling method, forward modeling system, storage medium, and electronic device | |
| CN110852025B (en) | Three-dimensional electromagnetic slow diffusion numerical simulation method based on super-convergence interpolation approximation | |
| CN114970290B (en) | Ground transient electromagnetic method parallel inversion method and system | |
| Ren et al. | Recursive analytical formulae of gravitational fields and gradient tensors for polyhedral bodies with polynomial density contrasts of arbitrary non-negative integer orders | |
| CN115113286B (en) | Aviation electromagnetic data fusion three-dimensional inversion method based on multi-component frequency domain | |
| Ying et al. | The phase flow method | |
| CN115755199A (en) | Practical unstructured grid three-dimensional electromagnetic inversion smoothing and regularization method | |
| Abubakar et al. | Three-dimensional visco-acoustic modeling using a renormalized integral equation iterative solver | |
| Yao et al. | 3D finite-element modeling of Earth induced electromagnetic field and its potential applications for geomagnetic satellites | |
| CN116090283A (en) | Aviation electromagnetic three-dimensional inversion method based on compressed sensing and preconditioned random gradient | |
| CN115292973B (en) | Arbitrarily sampled space wave number domain three-dimensional magnetic field numerical simulation method and system | |
| CN112526621B (en) | Ground-air electromagnetic data slow diffusion multi-parameter extraction method based on neural network | |
| Pan et al. | Three-dimensional forward modelling of gravity field vector and its gradient tensor using the compact difference schemes | |
| Zaynidinov et al. | Digital processing of two-dimensional signals in the basis of Haar wavelets | |
| Wu et al. | Solving parametric elliptic interface problems via interfaced operator network | |
| Vujević et al. | Potential distribution for a harmonic current point source in horizontally stratified multilayer medium | |
| Luo et al. | Reconstructing airborne gravity gradient data based on Fourier transform | |
| Baxansky et al. | Calculating geometric properties of three-dimensional objects from the spherical harmonic representation | |
| Yin et al. | A fast 3D gravity forward algorithm based on circular convolution | |
| CN114818416A (en) | Magnetotelluric simulation method and device based on non-fitting grid and storage medium |
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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |