CN115113286A - 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
- CN115113286A CN115113286A CN202210799540.5A CN202210799540A CN115113286A CN 115113286 A CN115113286 A CN 115113286A CN 202210799540 A CN202210799540 A CN 202210799540A CN 115113286 A CN115113286 A CN 115113286A
- Authority
- CN
- China
- Prior art keywords
- inversion
- data
- frequency domain
- model
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 94
- 230000004927 fusion Effects 0.000 title claims abstract description 29
- 239000011159 matrix material Substances 0.000 claims abstract description 36
- 230000006870 function Effects 0.000 claims description 30
- 238000004364 calculation method Methods 0.000 claims description 24
- 230000005684 electric field Effects 0.000 claims description 15
- 230000008569 process Effects 0.000 claims description 14
- 230000004044 response Effects 0.000 claims description 10
- 238000005457 optimization Methods 0.000 claims description 7
- 238000005259 measurement Methods 0.000 claims description 6
- 230000010354 integration Effects 0.000 claims description 3
- 230000035699 permeability Effects 0.000 claims description 3
- 238000009499 grossing Methods 0.000 claims description 2
- 230000002194 synthesizing effect Effects 0.000 claims description 2
- 238000010276 construction Methods 0.000 abstract description 2
- 238000004422 calculation algorithm Methods 0.000 description 8
- 230000002159 abnormal effect Effects 0.000 description 5
- 230000008859 change Effects 0.000 description 5
- 230000002547 anomalous effect Effects 0.000 description 3
- 238000009795 derivation Methods 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 238000012545 processing Methods 0.000 description 3
- 238000004088 simulation Methods 0.000 description 3
- 238000002939 conjugate gradient method Methods 0.000 description 2
- 238000001514 detection method Methods 0.000 description 2
- 238000003384 imaging method Methods 0.000 description 2
- 238000012804 iterative process Methods 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 238000007792 addition Methods 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000021615 conjugation Effects 0.000 description 1
- 238000001816 cooling Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000002474 experimental method 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
- 230000001788 irregular Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000001737 promoting effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
Images
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 three-dimensional inversion method based on multi-component frequency domain aviation electromagnetic data fusion, which comprises the following steps of firstly, based on the flexible modeling characteristic of a tetrahedral mesh, developing three-dimensional forward modeling of frequency domain aviation electromagnetic finite elements under any complex structure and undulating terrain conditions, and providing effective precondition guarantee for refined three-dimensional inversion; secondly, approximation of a hessian matrix is achieved based on a finite memory quasi-Newton method, and the aviation electromagnetic inversion speed is effectively accelerated; the inversion of the method is more accurate and more underground structure information can be obtained by comparing the single data inversion result with the multiple data fusion inversion result, and a new thought is provided for the fine construction explanation based on the frequency domain aviation electromagnetic data.
Description
Technical Field
The invention belongs to the technical field of geophysical signal processing, and particularly relates to a three-dimensional inversion method based on multi-component frequency domain aviation electromagnetic data fusion.
Background
Resource exploration needs to carry out secondary fine detection on marching and near-surface areas in areas with complex geological conditions. The aeroelectromagnetic method takes an airplane as a carrier, has the advantages of high exploration speed, low consumption cost and the like, can quickly realize regional electromagnetic data scanning acquisition work, and becomes one of the preferred technical methods for exploration work in regions with complicated geological conditions. 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 necessary guarantee is provided for the near-surface fine detection. However, as the structure of the exploration target becomes complex and refined, this presents a significant challenge to the processing of massive airborne electromagnetic data under high sampling rate conditions.
The interpretation of the airborne electromagnetic actual measurement data mainly ensures that the rapid processing imaging technology and the one-dimensional inversion of mass data are mainly performed, and the one-dimensional imaging and the inversion are based on the assumption of a layered or geoelectric model and are not enough to meet the fine inversion of a complex geoelectric model. Therefore, the three-dimensional aeroelectromagnetic inversion capable of recovering the fine and complex structure of the geoelectrical model is the key technical means for solving the interpretation of the aeroelectromagnetic mass data and is also the key technical problem for promoting complex geological conditions and near-surface aeroelectromagnetic exploration.
In addition, most frequency domain electromagnetic data acquisition systems take a helicopter nacelle as a platform and simultaneously carry two sets of coil pairs which are vertical, coaxial and horizontal and coplanar to acquire two sets of orthogonal components, and in data interpretation at present, data inversion is generally only based on each independent coil pair, which hardly meets the requirement of fine structure interpretation, so that the development of fusion inversion of data of the two sets of coil pairs is imperative.
The forward modeling is part of the inversion core technique. Among the forward algorithms, the finite difference algorithm is relatively simple and easy to program, and is most used in the current three-dimensional forward modeling, and then is an integral equation method and a finite volume method, and finally is a finite element method. The finite element method can adapt to various irregular grids such as tetrahedrons, deformed hexahedrons and the like, but the finite element method has high calculation memory consumption and limits the application of the finite element method in early three-dimensional forward modeling. As computer hardware has evolved, finite element methods based on unstructured grids have become the focus of current research. Currently, three-dimensional inversion algorithms are mature, and mainly include a nonlinear conjugate gradient method, a gaussian-newton method, a quasi-newton method, and a finite-memory quasi-newton method. The nonlinear conjugate gradient method avoids the Hessian matrix solution, greatly reduces the calculated amount, but has slow convergence speed. The Gaussian-Newton method ignores the second order term of the Hessian matrix and improves the convergence rate. The quasi-Newton method adopts an iteration method to obtain an approximate Hessian matrix in inversion, and compared with the traditional quasi-Newton algorithm, the finite memory quasi-Newton method effectively reduces the calculation memory and is suitable for the inversion of mass electromagnetic data. The current three-dimensional inversion method mainly adopts a finite difference method based on a structural grid, however, for an area with strong topographic relief, an air and ground interface needs to be split finely enough, the number of grids is large, memory consumption is huge, calculation time is long, and the application of the method in a complex terrain and a structural area is limited.
Moreover, the existing finite difference method and finite volume algorithm adopting regular grids cannot fit the terrain, the calculation error is large, and meanwhile, the existing inversion technology is independent inversion and strong in non-uniqueness. Both data cannot be utilized simultaneously.
Disclosure of Invention
In view of the above, the present 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 matches with a non-structural tetrahedral mesh, can solve the problem of efficiently fitting terrain, and has high calculation accuracy; meanwhile, the fusion inversion of two measured data, namely HCP and VCX, is adopted, so that the non-uniqueness of the three-dimensional inversion can be effectively reduced, and more information of the underground structure can be obtained.
In order to achieve the above object, an embodiment of the present invention provides a three-dimensional inversion method based on multi-component frequency domain airborne electromagnetic data fusion, including the following steps:
and 5, performing optimization calculation on the inversion objective function by using a quasi-Newton method with limited memory, adding partial derivatives of data fitting items into a finite element forward equation solution in the optimization process, and updating the model through iterative updating solution to realize the final inversion model.
In step 2 of an embodiment, when a finite element forward modeling based on an unstructured tetrahedral mesh is performed on an initial model in a frequency domain aeroelectromagnetic system, a frequency domain secondary 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 Conductivity for background field;
to ensure uniqueness of the quadratic field solution, at the boundary
Using dirichlet first class boundary conditions:
the method comprises the following steps of multiplying a formula (2) by a vector basis function N point, performing step-by-step integration by utilizing a Green first formula to obtain a coefficient matrix equation of each tetrahedral mesh unit, synthesizing the coefficient matrixes of all the tetrahedral mesh units into an overall matrix K, and obtaining a finite element forward equation of the frequency domain electromagnetic method:
KΕ s ＝S (3)
wherein S is a source item in a tetrahedral mesh unit, and a secondary magnetic field H of the frequency aviation electromagnetic response is obtained by solving a finite element forward equation s 。
wherein phi d (m) represents a data fit term, N is the number of data,
is expressed as measured multi-component frequency domain aeroelectromagnetic data, d n Secondary magnetic field H calculated for forward modeling s A magnetic field component of e n And m is the geological attribute of the model, namely the variable obtained by inversion.
Ф m (m)＝α r ‖Rm‖ 2 +α s ‖W s (m-m pri )‖ 2 (5)
wherein alpha is r And alpha s Respectively weight, R is a model smoothing matrix, m is a geological property of the model, m pri For a priori geological properties as model prior information, W s A prior constraint matrix of the model;
inverse objective function Φ, expressed as:
Ф＝Ф d (m)+λФ m (m)
＝Ф d (m)+λ(α r ‖Rm‖ 2 +α s ‖W s (m-m pri )‖ 2 ) (6)
Φ d (m) fitting dataAn item.
In step 5 of an embodiment, when the finite-memory quasi-newton method is used to solve and calculate the inversion objective function, in order to calculate the gradient of the data fitting term in the inversion objective function, the actually measured multi-component frequency domain airborne electromagnetic data is used
Secondary magnetic field H obtained by forward modeling s Magnetic field component d of n Separately arranged according to the real part and the imaginary part to obtain:
wherein N is the data index, N is the data number, Δ d n Is composed of
And d n The difference between the two;
and the gradient of the data fit term is expressed as:
for the kth tetrahedral mesh cell there are:
wherein, M is the total number of the tetrahedral mesh units,
an operator representing the real part of the complex number;
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 an embodiment, in the process of solving the inversion objective function by using the finite-memory quasi-newton method, after the partial derivative of the data fitting term is added to the finite-element forward equation solution, forward calculation is performed at each frequency, wherein one forward calculation is to calculate the secondary electric field E in the finite-element forward equation s Secondary magnetic field H of sum frequency aviation electromagnetic response s (ii) a Another forward calculation is accompanied by different frequency aviation electromagnetic response secondary magnetic field H s The adjoint forward calculation of the components to solve the overall matrix K.
In step 5 of an embodiment, in the process of solving the inversion objective function by using the quasi-newton method with limited memory, the hessian matrix is replaced by an approximate matrix constructed by the quasi-newton method, so that the approximation of the hessian matrix is realized, and the aviation electromagnetic inversion speed is accelerated.
In step 5 of an embodiment, during the process of solving the inversion objective function by using a finite-memory quasi-newton method, a model constraint term is directly calculated according to model prior information and an updated model obtained by the last solution.
Compared with the prior art, the invention has the beneficial effects that at least:
according to the technical scheme provided by the embodiment, a finite element method based on a non-structural tetrahedral mesh is combined with a finite memory quasi-Newton method to carry out fusion three-dimensional inversion of multi-component aviation electromagnetic data, firstly, based on the flexible modeling characteristic of the tetrahedral mesh, frequency domain aviation electromagnetic finite element three-dimensional forward simulation under the conditions of any complex structure and undulating terrain can be carried out, and effective precondition guarantee is provided for refined three-dimensional inversion; secondly, approximation of a hessian matrix is achieved based on a finite memory quasi-Newton method, and the aviation electromagnetic inversion speed is effectively accelerated; the comparison of a single data inversion result and a plurality of data fusion inversion results shows that the provided aviation electromagnetic data fusion three-dimensional inversion method based on the multi-component frequency domain is more accurate in inversion, can obtain more underground structure information, and provides a new idea for the fine construction and explanation of aviation electromagnetic data based on the frequency domain.
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 three-dimensional inversion method based on multi-component frequency domain airborne electromagnetic data fusion 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 an HCP device provided by an embodiment;
FIG. 4 is an inversion result of measured data collected by a VCX device according to an embodiment;
FIG. 5 is an inversion result of measured data collected by an aviation electromagnetic HCP device and measured data collected by an aviation electromagnetic VCX device according to the embodiment;
FIG. 6 is an inversion convergence graph of measured data collected for an HCP device provided by an embodiment;
FIG. 7 is an inversion convergence diagram of measured data collected by a VCX device according to an embodiment;
fig. 8 is an inversion convergence diagram of measured data acquired by an airborne electromagnetic HCP device and measured data acquired by an airborne electromagnetic VCX device according to an embodiment.
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.
FIG. 1 is a flowchart of a three-dimensional inversion method based on multi-component frequency domain airborne electromagnetic data fusion, provided by the embodiment. As shown in fig. 1. The embodiment provides a three-dimensional inversion method based on multi-component frequency domain aviation electromagnetic data fusion, which comprises the following steps:
In the embodiment, an initial model is constructed according to the inversion test geological structure, and the geological properties of the initial model are configured, wherein the geological properties comprise magnetic field information such as conductivity, resistivity and the like. In view of the flexible modeling characteristic based on the tetrahedral mesh and in view of the structural complexity and the topographic relief of the geological structure, the initial model is divided by adopting the non-structural tetrahedral mesh so as to facilitate the three-dimensional forward modeling of the aeronautical electromagnetic finite element in the frequency domain.
And 2, carrying out finite element forward modeling on the initial model based on the non-structural tetrahedral mesh.
In the embodiment, when the finite element forward modeling is carried out on the initial model in the frequency domain aeroelectromagnetic system, the secondary field and the total field are measured. And the literature indicates that the quadratic field has a smaller influence range (football print) than the total field, so that the quadratic field equation is solved in the finite element forward process. Assuming a positive time-harmonic e iωt Then, the frequency domain secondary field electric field satisfies the following helmholtz equation:
wherein E 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 Background field conductivity. To ensure uniqueness of the quadratic field solution, here at the boundary
Using dirichlet first class boundary conditions:
(1) the formula (2) and the formula (1) together form a forward theoretical boundary value problem. (1) The formula and the vector basis function are multiplied by N points, and the Green first formula is used for step-by-step integration to obtain
dV is the rightmost term of the integral operator, representing any minute volume within the region omega. By decomposing the calculation region of equation (3) into a series of tetrahedral mesh cells and calculating the integral of equation (3) within each tetrahedral mesh cell, a small coefficient matrix equation can be obtained
Where i is the index of the tetrahedral network element, K e Is a coefficient matrix composed of the sum of a density matrix and a stiffness matrix, S e For tetrahedral network element inliers, E e s is the secondary electric field of the tetrahedral network unit, the coefficient matrixes of all the tetrahedral network units are synthesized into a total matrix K, and a finite element forward equation of the frequency domain electromagnetic method can be obtained
KΕ s ＝S (5)
Wherein S is a source term matrix, and the equation (5) is solved through a multi-core parallel direct solver. Then, the secondary magnetic field H of the frequency aviation electromagnetic response s Can be obtained by faraday's theorem.
In frequency domain airborne electromagnetic data simulation, data of a fixed airborne electromagnetic HCP device or VCA device are processed, wherein the airborne electromagnetic HCP device means that a transmitting coil and a receiving coil are both vertical magnetic dipoles, and the airborne electromagnetic VCA device means that the transmitting coil and the receiving coil are both horizontal magnetic dipoles. In the embodiment, the measured data collected by the aviation electromagnetic HCP device and the measured data collected by the aviation electromagnetic VCX device aiming at the same geological structure are simultaneously obtained to form measured multi-component frequency domain aviation electromagnetic data, and the data is 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 target function based on the data fitting item and the model constraint item.
In the embodiment, an inversion target function Φ of the nonlinear frequency domain three-dimensional aeroelectromagnetic method constructed based on the data fitting term and the model constraint term is expressed as:
Ф＝Ф d (m)+λФ m (m) (6)
where λ is the balance data fitting term Φ d (m) and the model constraint term phi m The real number of (m), called the regularization factor. The three-dimensional aeroelectromagnetic inversion problem is equivalent to minimizing the objective function problem described above. Reasonable regularization factor may guarantee phi m (m) and phi d (m) does not over fit in Φ. According to equation (6), the gradient of the objective function can be expressed as:
and based on a non-structural finite element method, the calculation of the three-dimensional aeroelectromagnetic target function and the gradient thereof is described in detail. Phi ( d (m) may be expressed as:
wherein N is the number of data,
is expressed as measured multi-component frequency domain aeroelectromagnetic data, d n Secondary magnetic field H found for forward modeling s Of the magnetic field component H, i.e. the magnetic field component H x 、H y And H z ，e n As data error, m is the geology of the modelAnd (4) sex, namely obtaining the variable obtained by inversion.
Model constraint term phi m (m) may be expressed as:
Φ m (m)＝α r Φ r (m)+α s Φ s (m) (9)
wherein phi m (m)＝||Rm|| 2 (10)
Φ s (m)＝||W s (m-m pri )|| 2 (11)
Wherein alpha is r And alpha s Are respectively weight, R is a model constraint matrix, m is a geological attribute of the model, m is pri For a priori geological properties as model prior information, W s A prior constraint matrix of the model;
the objective function can be rewritten as:
Φ(m)＝Φ d (m)+λ(α r Φ r (m)+α s Φ s (m)) (12)
for unstructured tetrahedral mesh, the difference matrix R in equation (10) does not need to be solved by display, and W can be solved by the difference between each tetrahedral mesh unit and the adjacent tetrahedral mesh unit s Which can be expressed as the inverse of the volume of the current model, then there are:
wherein, V i Cell volume, N, for the ith tetrahedral mesh cell (i) Is the number of adjacent cells of the ith tetrahedral mesh cell, M is the total number of tetrahedral mesh cells, d ij Is the distance of the centers of adjacent cells, w j Is the weight coefficient of the jth neighboring cell, Δ m ij Model parameter difference, w, of neighboring cells j And Δ m ij Are defined as follows:
let the ith grid cell center point be
Then:
in summary, the gradient of the objective function can be expressed as:
and 5, performing optimal calculation on the inversion objective function by using a finite memory quasi-Newton method to update the model and realize the final inversion model.
In the embodiment, in order to calculate the gradient of a data fitting item in an inversion target function, the actually measured multi-component frequency domain aeroelectromagnetic secondary magnetic field data is subjected to
Secondary magnetic field H obtained by forward modeling s Magnetic field component d of n Respectively arranging according to the real part and the imaginary part separately, wherein the first N data are the real part of the magnetic field component, and the last N data are the imaginary part of the magnetic field component, so as to obtain:
wherein N is the data index, N is the data number,
representing plural forms
d n ' denotes d in complex form n ，Δd n Is composed of
And d n The difference, and the gradient of the data fit term can be expressed as:
thus, for the ith tetrahedral unit there is:
wherein, is complex conjugation.
In an embodiment, these partial derivatives, which can be shown in equations (20) and (21), are applied to the forward equation solution E during the optimization process. Is provided with
For a given frequency and a magnetic field component H of the magnetic field at point j interpolated from the electric field vector E at the emission source x,j And then:
in the same way, can define
Similar to equation (21) there are:
forward modeling of aviation electromagnetic methodEquation (5) for m k The derivation can be:
from the above derivation, equation (24) can be rewritten as:
wherein,
dependent on different aeronautical electromagnetic response components, e.g. for H x Corresponding to
Is composed of
Defining a single frequency all source simulation is called a forward run, equation (25) requires two forward runs at each frequency. One of the forward runs is calculation E and response H; it is additionally necessary twice for the right-hand term to be g T Is performed to solve an overall matrix K in which the coefficients g are relatively stable T Comprises the following steps:
thus, let v T Satisfies the following conditions:
v T ＝g T K -1 (27)
thus, a companion forward can be defined:
Kv＝g (28)
the complete calculation of equation (25) can be performed by equation (28).
Fitting term phi to data d Model constraint term phi m Can be directly analyzed to obtain:
the model constraint term is directly calculated according to the model prior information and the updated model obtained by the last solution.
For an objective function phi with M unknowns, the number can reach tens of thousands to millions in the multi-component aviation electromagnetic method inversion, in order to accelerate convergence, an inversion process of a finite-memory-based quasi-Newton optimization algorithm suitable for three-dimensional aviation electromagnetic data is adopted, and the algorithm is developed from a quasi-Newton method. The quasi-Newton method adopts an approximate matrix B in the k +1 iteration process of the model k To replace Hessian matrix H k According to the nature of the Taylor series, the objective function Φ (m) is at m k+1 The quadratic approximation of (d) can be expressed as:
the derivation of the above formula is:
g(m)≈g k+1 +H k+1 (m-m k+1 )，·················· (32)
let m be m k ，s k ＝m k+1 -m k ，y k ＝g k+1 -g k Then, there are:
H k+1 s k ≈y k ，························ (33)
the quasi-newton method requires a structure:
B k+1 s k ＝y k .························ (34)
the calculation relationship 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 the memory and improve the efficiency on the basis of the Newton method, and the transformation of the equation (35) by using a Sherman-Morrison formula can obtain:
by F k+1 To represent
Order to
And I is an identity matrix. Then the above equation reduces to:
considering the approximate hessian inverse matrix in the previous several adjacent iterations, an initial F is given first 0 The latest m-times correction relation is kept in the finite memory quasi-Newton method, and the following can be obtained:
it can be seen that the limited memory quasi-Newton method, stores only the column vector pairs s k And { y } k F is not stored k+1 。
The specific process can be summarized as follows:
1) setting inversion termination condition, setting precision epsilon > 0 and target root mean square value r 0 Let k equal to 1;
2) given an initial positive definite diagonal matrix
Wherein the scale factor
3) Computing
And RMS, if
Or RMS ≦ r 0 If the iterative inversion is ended, m is output k ；
4) Let n be min { k, m-1}, pair
Updating for n +1 times to obtain F k+1 The expression of (1);
5) finding Δ m using double-loop recursion k ＝-F k g k ；
6) Setting beta 'to be more than 0 and less than 0.5 and beta' to be less than 1 according to the Wolfe condition, and searching the step length alpha in a linear mode k ：
Let m k+1 ＝m k +α k d k ；
7) Let k be k +1, return to step 2).
According to the process, the finite-memory Newton optimization algorithm is simulated, and therefore the three-dimensional inversion of the multi-component frequency domain aviation electromagnetic data is achieved.
The aviation electromagnetic data fusion three-dimensional inversion method based on the multi-component frequency domain provided by the embodiment sufficiently explores the advantages of data collected by two devices, namely a HCP device and a VCA device, develops two data fusion inversion technologies, calculates the modes of the two devices and solves the sensitivity of the two devices in the forward and concomitant forward processes, performs optimized calculation by using a finite-memory quasi-Newton method, and performs fusion inversion on the two data in an ordered arrangement manner.
Examples of the experiments
In order to test the effectiveness of the method for three-dimensional inversion based on 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 the background space, and 5% of Gaussian noise is added into the forward-modeled magnetic field data to simulate the measured data. FIG. 2 shows the abnormal body and 225 measuring point distribution of the model, wherein, (a) is a top view of the model, and (b) is a transverse section view of the model, and the distance between each measuring point and a lateral line is 20 m. The resistivity of the background half space is 100 Ω · m, the anomalous body is located 20m directly below the center of the measurement point, the size of the anomalous body is 100m × 100m × 25m, and the resistivity of the anomalous body is 10 Ω · m. The grid used by the forward theoretical model contains 431223 cells, calculating the HCP and VCX responses at 2 frequencies, 900Hz and 5000Hz respectively. The grid used for the model inversion contains 214325 cells. The invention uses the same cooling principle to change lambda. Initial lambda was chosen to be 0.01, an additional regularization factor alpha r And alpha s They were all 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 in combining HCP and VCX data, separate HCP and VCX inversions were set simultaneously for comparison.
The experimental examples were performed on an Intel (R) Xeon (R) Gold 6254CPU 3.1GHz workstation. Fig. 3 shows the HCP inversion alone, where (a) is the y-z section plot (x ═ 0 m); (b) x-z slice (y-0 m); (c) the actual anomaly location is represented in the boxed area for an x-y slice (z-30 m). Analyzing FIG. 3, it can be seen that HCP inversion can approximately restore the true horizontal position of the anomaly. There are, however, two disadvantages, the first: the inverted slice shows that the depth of the abnormal body is difficult to reach the true depth, and only the shallow abnormal body is effectively recovered; secondly, the method comprises the following steps: the contour of the horizontal plane of the abnormal body is deformed to approximate a circle smaller than the true range. Two disadvantages of HCP inversion directly affect the accuracy of data inversion.
Fig. 4 shows a single VCX inversion, where (a) is a y-z slice (x ═ 0 m); (b) x-z slice images (y ═ 0 m); (c) the actual anomaly location is represented in the boxed area for an x-y slice (z-30 m). Analysis of FIG. 4 reveals that the VCX inversion depth is deeper than the HCP and that the anomaly profile is closer to true. The two defects existing in single HCP inversion are overcome, but the VCX inversion slice has some small high-resistance false anomalies in a shallow part, and the anomalies bring errors to data inversion interpretation.
Fig. 5 is the result of the joint inversion of fused HCP and VCX, where (a) is a y-z section (x ═ 0 m); (b) x-z slice images (y ═ 0 m); (c) the actual anomaly location is represented in the boxed area for an x-y slice (z-30 m). By analyzing fig. 5, it can be found that the joint inversion is superior to the single HCP inversion in both inversion depth and horizontal profile, and overcomes the disadvantages of the single HCP inversion, and meanwhile, the joint inversion does not have a false high-resistance anomaly in a shallow part, which illustrates that the joint inversion can also overcome the disadvantages of the false anomaly of the single VCX inversion.
The effectiveness of the method for three-dimensional inversion based on multi-component frequency domain aviation electromagnetic data fusion provided by the embodiment is further illustrated. FIG. 6 shows the HCP inversion convergence plot, where (a) shows the data fitting difference Φ of the inversion over 112 iterations for the HCP example d The decrease from 100.1 to 2.6 shows the roughness Φ of the model r And phi s And finally, stable convergence; (b) showing that there are 12 HCP inversions that require changing the step size α of the inversion, (c) showing that there are 6 HCP inversions that change the regularization factor λ.
FIG. 7 shows the HCP inversion convergence plot, where (a) shows the data fit difference Φ for the inversion of the VCX example over 66 iterations d From 110.2 to 24.6, the roughness of the model Φ r And phi s And finally, stable convergence; (b) the HCP inversion is shown to have 14 times of step size alpha required to change the inversion; (c) the VCX inversion is shown with 7 changes in the regularization factor λ.
FIG. 8 shows a convergence plot for the HCP and VCX joint inversion, where (a) shows the data fit difference Φ for the inversion of the example over 48 iterations d From 110.2 to 26.1, the roughness of the model is reduced to phi r And phi s And finally, stable convergence; (b) shows that there are 5 times of combined HCP and VCX inversions, the step size alpha of inversion needs to be changed, (c) shows that there are 5 HCP and VCX inversionsThe regularization factor is changed sub-times. First, HCP inversion converges slowest in iteration number, 112 iterations, and VCX inversion 68 times, whereas joint inversion is fastest, requiring only 48 times. Secondly, VCX changes up to 14 times from the number of times the step size is changed in the iterative process of inversion, HCP12 times, and the combined HCP and VCX inversion changes 5 times, because if the step size is not changed, each iteration requires one forward, the iteration of changing the step size is required, the required forward number is more than 1, and the effectiveness of the combined inversion is also proved from the number of times the step size is changed. Finally, in terms of the change times of the regularization factors, VCX changes most, HCP changes are the least, and the combination inversion changes least, and the less change times represent the better inversion stability, so that the effectiveness of the aviation electromagnetic data fusion three-dimensional inversion method based on the multi-component frequency domain is also proved.
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 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 performing non-structural tetrahedral mesh division on the initial model;
step 2, carrying out finite element forward modeling based on the non-structural tetrahedral mesh on the initial model;
step 3, acquiring actually measured multi-component frequency domain aviation electromagnetic data, comprising: actual measurement data collected by the aviation electromagnetic HCP device and actual measurement data collected by the aviation electromagnetic VCX device;
step 4, constructing a data fitting item according to the multi-component frequency domain aeroelectromagnetic data, constructing a model constraint item according to model prior information, and constructing an inversion target 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 with limited memory, adding partial derivatives of data fitting items into a finite element forward equation solution in the optimization process, and updating the model through iterative updating solution to realize the final inversion model.
2. The method for three-dimensional inversion based on multi-component frequency domain airborne electromagnetic data fusion according to claim 1, wherein in step 2, when finite element forward modeling based on non-structural tetrahedral mesh is performed on an initial model in a frequency domain airborne electromagnetic system, a frequency domain secondary 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 Conductivity for background field;
to ensure uniqueness of the quadratic field solution, at the boundary
Using dirichlet first class boundary conditions:
the method comprises the following steps of multiplying a formula (2) by a vector basis function N point, performing step-by-step integration by utilizing a Green first formula to obtain a coefficient matrix equation of each tetrahedral mesh unit, synthesizing the coefficient matrixes of all the tetrahedral mesh units into an overall matrix K, and obtaining a finite element forward equation of the frequency domain electromagnetic method:
KΕ s ＝S (3)
wherein S is a source item in a tetrahedral mesh unit, and a secondary magnetic field H of the frequency aviation electromagnetic response is obtained by solving a finite element forward equation s 。
3. The method for three-dimensional inversion based on multi-component frequency domain aviation electromagnetic data fusion according to claim 1, wherein in step 4, a data fitting term constructed according to measured data is expressed as:
wherein phi d (m) represents a data fit term, N is the number of data,
the table is measured multi-component frequency domain aviation electromagnetic data, d n Secondary magnetic field H found for forward modeling s A magnetic field component of e n And m is the geological attribute of the model, namely the variable obtained by inversion.
4. The method for three-dimensional inversion based on multi-component frequency domain aviation electromagnetic data fusion according to claim 1 or 3, characterized in that 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 weight, R is a model smoothing matrix, m is a geological property of the model, m pri For a priori geological properties as model prior information, W s A prior constraint matrix of the model;
inverse objective function Φ, expressed as:
Ф＝Ф d (m)+λФ m (m)
＝Ф d (m)+λ(α r ‖Rm‖ 2 +α s ‖W s (m-m pri )‖ 2 ) (6)
Φ d (m) is a data fit term.
5. The method for three-dimensional inversion based on multi-component frequency domain airborne electromagnetic data fusion according to claim 1, wherein in step 5, when a finite memory quasi-Newton method is used for solving and calculating an inversion target function, in order to calculate the gradient of a data fitting item in the inversion target function, actually measured multi-component frequency domain airborne electromagnetic data is subjected to calculation
Secondary magnetic field H obtained by forward modeling s Magnetic field component d of n Separately arranging according to the real part and the imaginary part to obtain:
wherein N is the data index, N is the data number, Δ d n Is composed of
And d n The difference between the two;
and the gradient of the data fit term is expressed as:
for the kth tetrahedral mesh cell there are:
wherein, M is the total number of the tetrahedral mesh units,
an operator representing the real part of the complex number;
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.
6. The method for three-dimensional inversion based on multi-component frequency domain airborne electromagnetic data fusion according to claim 1, wherein in the step 5, in the process of solving the inversion objective function by using a finite-memory quasi-Newton method, forward calculation is performed at each frequency after partial derivatives of data fitting terms are added to finite-element forward equation solutions, wherein one forward calculation is used for calculating the secondary electric field E in the finite-element forward equation s Secondary magnetic field H of sum frequency aviation electromagnetic response s (ii) a Another forward calculation is accompanied by different frequency aviation electromagnetic response secondary magnetic field H s The adjoint forward calculation of the components to solve the overall matrix K.
7. The method for three-dimensional inversion based on multi-component frequency domain airborne electromagnetic data fusion of claim 1, wherein in step 5, in the process of solving the inversion objective function by using a finite-memory quasi-Newton method, the Hessian matrix is replaced by an approximate matrix constructed by the quasi-Newton method, so that the approximation of the Hessian matrix is realized, and the airborne electromagnetic inversion speed is accelerated.
8. The method for three-dimensional inversion based on multi-component frequency domain airborne electromagnetic data fusion according to claim 1, wherein in the step 5, a finite memory quasi-Newton method is used for solving the inversion objective function, and the model constraint term is directly calculated according to model prior information and the updated model obtained by the last solution.
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 true CN115113286A (en) | 2022-09-27 |
| CN115113286B 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) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN118244370A (en) * | 2024-05-28 | 2024-06-25 | 吉林大学 | Three-dimensional inversion method for transient electromagnetic tunnel advanced detection based on finite element algorithm |
Citations (7)
| 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 |
| CN106199742A (en) * | 2016-06-29 | 2016-12-07 | 吉林大学 | A kind of Frequency-domain AEM 2.5 ties up band landform inversion method |
| US20170003409A1 (en) * | 2015-07-02 | 2017-01-05 | Volkan AKCELIK | Krylov-Space-Based Quasi-Newton Preconditioner for Full-Wavefield Inversion |
| 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 |
-
2022
- 2022-07-06 CN CN202210799540.5A patent/CN115113286B/en active Active
Patent Citations (7)
| 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 |
| US20170003409A1 (en) * | 2015-07-02 | 2017-01-05 | Volkan AKCELIK | Krylov-Space-Based Quasi-Newton Preconditioner for Full-Wavefield Inversion |
| 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 |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN118244370A (en) * | 2024-05-28 | 2024-06-25 | 吉林大学 | Three-dimensional inversion method for transient electromagnetic tunnel advanced detection based on finite element algorithm |
| CN118244370B (en) * | 2024-05-28 | 2024-07-23 | 吉林大学 | Three-dimensional inversion method for transient electromagnetic tunnel advanced detection based on finite element algorithm |
Also Published As
| Publication number | Publication date |
|---|---|
| CN115113286B (en) | 2023-09-15 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106199742B (en) | A kind of dimension of Frequency-domain AEM 2.5 band landform inversion method | |
| Dines et al. | Analysis of electrical conductivity imaging | |
| CN112949134B (en) | Earth-well transient electromagnetic inversion method based on non-structural finite element method | |
| Persova et al. | Methods and algorithms for reconstructing three-dimensional distributions of electric conductivity and polarization in the medium by finite-element 3D modeling using the data of electromagnetic sounding | |
| CN110852025B (en) | Three-dimensional electromagnetic slow diffusion numerical simulation method based on super-convergence interpolation approximation | |
| Donepudi et al. | A novel implementation of multilevel fast multipole algorithm for higher order Galerkin's method | |
| CN114970290B (en) | Ground transient electromagnetic method parallel inversion method and system | |
| CN110244351A (en) | A kind of Uniform Construction inversion method of different constraint Geophysical Inverse Problems | |
| CN105717547A (en) | Anisotropy medium magnetotelluric meshless value simulation method | |
| CN112733364B (en) | Foil cloud scattering rapid calculation method based on impedance matrix partitioning | |
| Cao et al. | 3D magnetotelluric inversions with unstructured finite-element and limited-memory quasi-Newton methods | |
| Zakaria | The finite-element contrast source inversion method for microwave imaging applications | |
| CN115113286B (en) | Aviation electromagnetic data fusion three-dimensional inversion method based on multi-component frequency domain | |
| CN116842813B (en) | Three-dimensional geoelectromagnetic forward modeling method | |
| CN116090283A (en) | Aviation electromagnetic three-dimensional inversion method based on compressed sensing and preconditioned random gradient | |
| Lee et al. | Application of whale optimization algorithm to inverse scattering of an imperfect conductor with corners | |
| LU505850B1 (en) | Fast seismic wave travel time calculation method for tunnel detection | |
| CN112526621B (en) | Ground-air electromagnetic data slow diffusion multi-parameter extraction method based on neural network | |
| Lu et al. | Three-dimensional magnetotelluric inversion using L-BFGS | |
| Lin et al. | Three-dimensional conjugate gradient inversion of magnetotelluric sounding data | |
| CN115563791B (en) | Magnetotelluric data inversion method based on compressed sensing reconstruction | |
| CN113627027B (en) | Method and system for simulating electromagnetic field value of non-trivial anisotropic medium | |
| Zhong et al. | Electrical resistivity tomography with smooth sparse regularization | |
| Yin et al. | A fast 3D gravity forward algorithm based on circular convolution | |
| Cheng et al. | 3D Step-by-step inversion strategy for audio magnetotellurics data based on unstructured mesh |
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 |