CN114297905B - Quick numerical simulation method of two-dimensional earth electromagnetic field - Google Patents
Quick numerical simulation method of two-dimensional earth electromagnetic field Download PDFInfo
- Publication number
- CN114297905B CN114297905B CN202210227897.6A CN202210227897A CN114297905B CN 114297905 B CN114297905 B CN 114297905B CN 202210227897 A CN202210227897 A CN 202210227897A CN 114297905 B CN114297905 B CN 114297905B
- Authority
- CN
- China
- Prior art keywords
- electric field
- field
- polarization mode
- total
- magnetic field
- 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 90
- 238000004088 simulation Methods 0.000 title claims abstract description 42
- 230000005672 electromagnetic field Effects 0.000 title claims abstract description 17
- 238000004364 calculation method Methods 0.000 claims abstract description 21
- 230000005684 electric field Effects 0.000 claims description 105
- 230000010287 polarization Effects 0.000 claims description 83
- 238000012937 correction Methods 0.000 claims description 12
- 230000002159 abnormal effect Effects 0.000 claims description 10
- 238000011160 research Methods 0.000 claims description 7
- 238000006467 substitution reaction Methods 0.000 claims description 4
- 230000009466 transformation Effects 0.000 claims description 4
- 230000002547 anomalous effect Effects 0.000 claims description 3
- 230000002194 synthesizing effect Effects 0.000 claims description 2
- 238000012795 verification Methods 0.000 abstract description 4
- 230000008569 process Effects 0.000 description 6
- 230000008901 benefit Effects 0.000 description 4
- 230000005477 standard model Effects 0.000 description 3
- 238000001514 detection method Methods 0.000 description 2
- 238000006073 displacement reaction Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000035699 permeability Effects 0.000 description 2
- 125000002015 acyclic group Chemical group 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 239000002184 metal Substances 0.000 description 1
- 230000005855 radiation Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 101150097388 rgy gene Proteins 0.000 description 1
- 230000001131 transforming effect Effects 0.000 description 1
Images
Classifications
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/30—Assessment of water resources
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
The invention provides a rapid numerical simulation method of a two-dimensional geodetic electromagnetic field, which combines a finite element method with Fourier transform, converts a two-dimensional partial differential problem into a one-dimensional constant differential problem with mutually independent different wave numbers by performing Fourier transform along the horizontal direction, has high parallelism, effectively reduces the calculated amount and improves the calculation efficiency; conditions are provided for the refined numerical simulation of the geoelectromagnetic method under the large-scale condition; through verification, the result calculated by adopting the method is within the error bar range of international published data included in the Zhdannov et al 1997 document, is closer to the mean value, and meets the precision requirement.
Description
Technical Field
The invention relates to the technical field of geophysical methods, in particular to a quick numerical simulation method of a two-dimensional earth electromagnetic field.
Background
The magnetotelluric method is a frequency domain electromagnetic prospecting method which utilizes a natural electromagnetic field (with a frequency band of 10) vertically incident in the air-4-104Hz) and then realizes the detection of underground geological structures and abnormal bodies, has simple operation, high working efficiency, low cost and deep detectionHigh degree and the like. At present, the method is widely applied to the fields of geological exploration, metal ore and oil gas resource detection, geodynamics research and the like, and obtains more and more obvious economic and social benefits. However, successful development of this approach relies on a forward and backward interpretation of the subsurface geological model.
Because the data volume processed by the magnetotelluric forward inversion process is very huge, very long calculation time and calculation memory are often consumed in the actual application process, and the requirement on a computer is high, so that the application of the magnetotelluric method in the actual process is limited to a great extent. In the numerical simulation process of magnetotelluric, a double rotation equation about an electromagnetic field needs to be solved, but the equation has no analytic solution, so that the solution can be only carried out by a numerical method.
The prior literature discloses applications of a plurality of numerical simulation methods in magnetotelluric forward modeling, which mainly include a conventional integral equation method (IE), a finite difference method (FD), a finite element method (FE), and the like, and the numerical simulation methods are respectively good and bad in computational accuracy and computational efficiency. The integral equation method faces huge challenges when calculating a large abnormal body model; the grid subdivision of the finite element method and the formation of the coefficient matrix are relatively time-consuming, so that the application of the finite element method in large-scale electromagnetic exploration is limited; the finite difference method uses a structured grid, and is difficult to process undulating terrain and complex anomalies. And with the increase of the calculation scale, the calculation time and the operation memory are increased in a nonlinear way, so that the forward and backward interpretation of the large-scale geoelectromagnetic exploration is greatly limited.
In summary, it is necessary to solve the problems of calculation efficiency and calculation memory in magnetotelluric forward numerical simulation under the condition of ensuring a certain accuracy, so as to further realize rapid forward/reverse imaging of magnetotelluric, and make large-scale magnetotelluric data interpretation possible.
Disclosure of Invention
The invention aims to provide a quick numerical simulation method of a two-dimensional earth electromagnetic field, aiming at further improving the calculation efficiency and reducing a certain calculation memory on the premise of ensuring certain precision, and the specific technical scheme is as follows:
a method for rapid numerical simulation of a two-dimensional earth electromagnetic field, comprising the steps of:
step A1: designing a two-dimensional geoelectric model according to the shape, size and conductivity distribution of the research area, and arranging measuring points; mesh subdivision is carried out on the two-dimensional earth electric model, and the conductivity of each node is assigned according to the electric distribution of the underground medium, so that a conductivity distribution model of the two-dimensional medium is obtained;
step A2: calculating a primary field corresponding to a background conductivity model of the research area;
step A3: constructing a linear equation set with different wave numbers by adopting a finite element method, which specifically comprises the following steps:
replacing a total electric field in a space domain secondary field control equation by a primary electric field in a primary field, and then performing Fourier transform in the horizontal direction to obtain a secondary field control equation of a space wave number domain;
for a secondary field control equation of a space wave number field, analyzing and totally synthesizing each unit by using a Galerkin method, and obtaining linear equation sets of different wave numbers by applying known boundary conditions;
step A4: solving linear equation sets with different wave numbers by adopting a catch-up method, performing inverse Fourier transform on the solution to obtain a secondary field, superposing the secondary field on a primary field to obtain a total electric field or a total magnetic field, further calculating to obtain the total electric field if the total magnetic field is obtained, and then performing iterative correction on the total electric field;
step A5: judging whether a convergence condition is reached or not according to the relative residual errors of the iteration results of the previous iteration and the next iteration, and repeating the step A3 and the step A4 if the convergence condition is not reached;
if the convergence condition is met, substituting the total electric field meeting the convergence condition into a space domain secondary field control equation to obtain a secondary field, and superposing the secondary field to the primary field to obtain a final total electric field or a final total magnetic field;
step A6: and after the final total electric field or the final total magnetic field is obtained, calculating apparent resistivity, impedance and phase on the measuring point.
In the above technical solution, preferably, in the step a2, the primary electric field in the primary field is obtained by calculation in the TE polarization mode, and the primary electric field and the primary magnetic field in the primary field are obtained by calculation in the TM polarization mode.
Preferably, in the above technical solution, in the step a 3: in the TE polarization mode, a primary electric field is adopted to replace a total electric field in a space domain secondary electric field control equation, and then Fourier transformation in the horizontal direction is carried out to obtain a secondary electric field control equation of a space wave number domain; in a TM polarization mode, a primary electric field is adopted to replace a total electric field in a space domain secondary magnetic field control equation, and then Fourier transformation in the horizontal direction is carried out to obtain a secondary magnetic field control equation of a space wave number domain.
Preferably, in the above technical solution, in the step a 4: obtaining a secondary electric field of a space domain after inverse Fourier transform in a TE polarization mode, and superposing the primary electric field to obtain a total electric field; and obtaining a secondary magnetic field of a space domain after inverse Fourier transform in a TM polarization mode, superposing the primary magnetic field to obtain a total magnetic field, and further calculating to obtain the total electric field.
Preferably, in the above technical solution, in the step a 5: if the convergence condition is met, substituting the total electric field meeting the convergence condition into a space domain secondary electric field control equation to obtain a secondary electric field in the TE polarization mode, and superposing the secondary electric field to the primary electric field to obtain a final total electric field; and under the TM polarization mode, substituting the total electric field meeting the convergence condition into a space domain secondary magnetic field control equation to obtain a secondary magnetic field, and superposing the secondary magnetic field to the primary magnetic field to obtain a final total magnetic field.
Preferably in the above technical solution, the step a3 specifically is:
in the TE polarization mode, the control equation of the spatial domain secondary electric field is:
is provided withFourier transform in the horizontal direction into two in the spatial wavenumber domain for equation (15)The sub-electric field control equation:
wherein,in terms of the wave number, the number of waves,a secondary electric field in the spatial wavenumber domain;is composed ofThe total electric field in the direction of the field,is a secondary electric field in the spatial domain,in order to be a background conductivity,in order to be an abnormal electrical conductivity,is the unit of an imaginary number,in order to have a magnetic permeability,is the angular frequency;
in the TM polarization mode, the spatial domain secondary magnetic field control equation is:
order to,And performing Fourier transform in the horizontal direction on the formula (20) to obtain a secondary magnetic field control equation of a space wave number domain:
wherein,,in the case of the background resistivity,is the anomalous resistivity;is a secondary magnetic field in the spatial wavenumber domain,is a secondary magnetic field in the spatial domain,are respectively an edgey、zA total electric field of direction;
edge obtained using step A2 in TE polarization modexDirectional primary electric fieldInstead of in equation (15)And is converted into publicSolving the formula (16);
edge obtained using step A2 in TM polarization modey、zDirectional primary electric fieldInstead of in equation (20)And then converted into a formula (21) for solving;
and (3) analyzing each unit by using a Galerkin method for the formula (16) or the formula (21) obtained after substitution, forming an algebraic equation system taking the electromagnetic field of the space wave number domain on the node as an unknown quantity, and imposing a first class boundary condition to obtain a linear equation system with the bandwidth of 5 and the diagonal dominance.
Preferably, in the above technical solution, the iterative correction is performed in step a4 according to formula (22):
wherein,is as followsThe electric field obtained after correction is carried out in the secondary iteration;is as followsElectric field obtained by sub-iteration without correction, in TE polarization modeTotal electric field obtained for step A4In TM polarization modeTotal electric field obtained for step A4Andand is andandrespectively adopting a formula (22) to carry out iterative updating; wherein ,。
In the above technical solution, preferably, in the step a5, the convergence determination condition is: relative residual error of two iterationsAnd then, the iteration is stopped,is the error limit.
Preferably in the above technical solution, the step a6 specifically is: step A5 obtaining the final total electric fieldOr the final total magnetic fieldThen, the partial derivative along the depth direction is obtained by using a numerical method, and the partial derivative is obtained in the TE polarization modeIn the TM polarization mode is;
In the TE polarization mode:
in the TM polarization mode:
wherein,respectively corresponding impedance, apparent resistivity and phase under the TE polarization mode;respectively corresponding impedance, apparent resistivity and phase under a TM polarization mode;respectively, imaginary part and real part.
Preferably, in the above technical solution, in the step a 2:
the TE polarization mode is:
the TM polarization mode is:
wherein,are respectively asThe total electric field in the three directions is,are respectively asTotal magnetic field in three directions;is the total conductivity;
in the step A4, obtaining the total magnetic field in TM polarization modeThen, the total electric field is obtained by further solving according to the formula (3)And。
the technical scheme of the invention has the following beneficial effects:
the invention provides a quick numerical simulation method of a two-dimensional earth electromagnetic field, which combines a finite element method with Fourier transform, converts a two-dimensional partial differential problem into a one-dimensional ordinary differential problem which is independent among different wave numbers by performing Fourier transform along the horizontal direction (y axis), and has high parallelism; the Fourier transform method is adopted to realize the dimension reduction of the two-dimensional magnetotelluric numerical simulation, and the complex two-dimensional problem is converted into a plurality of small problems, so that the calculation efficiency is effectively improved, and the calculation memory is reduced; and the fixed bandwidth equation set formed after the finite element method is dispersed is solved by adopting a catch-up method, so that efficient solution can be realized.
The method of the invention fully utilizes the high efficiency of Fourier transform and the accuracy of the finite element method, effectively improves the calculation efficiency of the numerical simulation of the magnetotelluric method, reduces the calculation time and provides conditions for the refined numerical simulation of the magnetotelluric method under the large-scale condition; through verification, the result calculated by adopting the method is in the error bar range of international published data included in the Zhdannov et al 1997 document, is closer to the mean value, and meets the precision requirement.
In addition to the objects, features and advantages described above, other objects, features and advantages of the present invention are also provided. The present invention will be described in further detail below with reference to the drawings.
Drawings
The accompanying drawings, which are incorporated in and constitute a part of this application, illustrate embodiments of the invention and, together with the description, serve to explain the invention and not to limit the invention. In the drawings:
FIG. 1 is a flow chart of a method for rapid numerical simulation of a two-dimensional magnetotelluric field of the present invention;
FIG. 2 is a schematic diagram of a background conductivity model provided by the present invention;
FIG. 3 is a schematic view of the COMMEM-2D 1 international standard model in the verification case of example 1;
FIG. 4a is a graph comparing the results of using the method of the present invention with that of Zhdannov et al 1997 in TE polarization mode;
FIG. 4b is a graph comparing the results of using the method of the present invention with that of Zhdannov et al 1997 in TM polarization mode.
Detailed Description
In order that the invention may be more fully understood, a more particular description of the invention will now be rendered by reference to specific embodiments thereof that are illustrated in the appended drawings. This invention may, however, be embodied in many different forms and should not be construed as limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete.
Unless defined otherwise, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention belongs. The terminology used in the description of the invention herein is for the purpose of describing particular embodiments only and is not intended to be limiting of the invention.
Example 1:
referring to fig. 1, the embodiment provides a method for fast numerical simulation of a two-dimensional earth electromagnetic field, which includes the following specific processes:
step S1: constructing a conductivity distribution model of a two-dimensional medium, which specifically comprises the following steps:
a two-dimensional earth model is designed from the shape, size and conductivity distribution of the study area, and the calculation frequency and the measuring point arrangement are determined.
The constructed two-dimensional geoelectrical model is arranged along the horizontal direction (yAxis), depth direction: (zAxis) was meshed and the edge was recordedyA shaft,zThe number of the mesh nodes divided in the axial direction is respectively ,yA shaft,zThe lengths between adjacent nodes in the axial direction are respectively recorded as;
And respectively assigning the conductivity of each node according to the electrical distribution of the underground medium, thereby completing a conductivity distribution model of the two-dimensional medium.
Step S2: selecting a polarization mode to be calculated, specifically selecting a TE polarization mode or a TM polarization mode;
for geoelectromagnetic methods, the conduction current in the subsurface medium is much larger than the displacement current, so neglecting the displacement current, the Maxwell equation set (Maxwell equation set) can be:
wherein:、respectively representing an electric field and a magnetic field,in order to have a magnetic permeability,is the angular frequency, and the angular frequency and the calculated frequencyIn a relationship of,As a result of the total electrical conductivity,in units of imaginary numbers.
For the two-dimensional magnetotelluric problem, equation (1) can be decoupled into two independent modes, namely TE and TM polarization modes.
The TE polarization mode is:
the TM polarization mode is:
in the formula,are respectively asThe total electric field in the three directions is,are respectively asTotal magnetic field in three directions;x、y、zthe axes are mutually perpendicular.
Step S3: calculating a primary field corresponding to a background conductivity model of the research area, specifically:
the background conductivity model, i.e. the one-dimensional horizontal uniform layered medium model, generally assigns the conductivity of the abnormal body part of the conductivity distribution model to the conductivity of the surrounding medium, and the primary field of the model can be solved by an analytical solution formula.
Referring to FIG. 2, assume that there is oneThe layered medium model is uniform in layer level, and the layer 1 medium is an air layer. In the TE polarization mode, we assume the natural field source (on top of the air layer, i.e. on top of the air layer)Where) isThe direction, then, can solve the form of the primary electric field and the primary magnetic field as:
wherein,;is a firstLayer mediumPrimary electric field with corresponding direction,Is as followsLayer mediumThe primary magnetic field with the corresponding direction is provided,is as followsLaminated mediumThe primary magnetic field with the corresponding direction is provided,is a unit of an imaginary number, and is,is as followsThe background conductivity of the layer medium,in order to measure the depth of the spot,is as followsThe depth of the top surface of the layer medium,,is as followsThe layer medium corresponds to the coefficients in equation (4),is a natural constant, is a constant in mathematics, and is an infinite acyclic decimal.
The tangential electric field and magnetic field on the electric interface satisfy the continuous boundary conditionObtaining:
wherein,is as followsThe electromagnetic reflection coefficient of the bottom surface of the layer medium,is a firstThe vertical distance between the top surface and the bottom surface of the layer medium,R i is as followsLayer mediumAndthe coefficient of the ratio therebetween.
The time electric field and the magnetic field are limited, and the radiation boundary condition is satisfied, soNOf a layer mediumAll can be calculated from the recursion formula (6) from bottom to topR i 。
from the boundary condition equation (5), we can obtain:
will be given in formula (8)Andsubstituting into equation (9), the coefficients of each layer can be determined sequentially from top to bottomAndwill beAndand the firstCorresponding depth in the layer mediumSubstituting into equation (4), the first in TE polarization mode in the model can be calculatedAny depth of layer mediumA primary electric field and a primary magnetic field.
In TM mode, we assume that the natural field source isThe direction, primary electric field and primary magnetic field are solved in the form:
in the formula,is as followsLayer mediumThe primary electric field with the corresponding direction is provided,is as followsLayer mediumA primary magnetic field with a direction corresponding to the direction,is a firstLayer mediumPrimary electric field with corresponding direction, coefficient of each layerAndand the primary field is solved in the same manner as in TE polarization mode, and will not be described here.
Step S4: constructing a linear equation set with different wave numbers by adopting a finite element method, which specifically comprises the following steps:
the total field is composed of a primary field and a secondary field (the primary field comprises a primary electric field and a primary magnetic field, and the secondary field comprises a secondary electric field and a secondary magnetic field), wherein the total field corresponds to the total conductivityPrimary field versus background conductivitySecondary field corresponding to abnormal conductivityWherein、Andthe relationship between the three is shown as formula (11):
electrical parameter of conductivity distribution modelyShaft andzaxis change, as obtained from equation (2), total electric field in TE polarization modeThe control equation satisfied is:
equation (12) minus equation (13) yields:
Here, if the fourier transform is directly performed on the formula (14),the term is changed from the product relationship of the space domain to the convolution relationship of the wavenumber domain, so that the control equation of the secondary electric field of the space domain can be obtained by further transforming the term:
is provided withAnd carrying out Fourier transform in the horizontal direction on the formula (15) to obtain a space wave number domain secondary electric field control equation:
in the formula,is a wave number of the wave number,is a secondary electric field in the spatial wavenumber domain.
From equation (3), the total magnetic field in TM polarization modeThe control equation satisfied is:
the secondary magnetic field is obtained by subtracting the formula (18) from the formula (17)The control equation of (1):
further simplification can lead to a space-domain quadratic magnetic field control equation (where,,):
order to,And Fourier transform in the horizontal direction is carried out on the formula (20) to obtain a space wave number domain secondary magnetic field control equation:
wherein,is a secondary magnetic field in the spatial wavenumber domain,,in the case of the background resistivity,is the anomalous resistivity.
At this time, in the TE polarization mode, in the formula (15)Is unknown, the primary electric field is used in this embodimentReplacing, and converting into a formula (16) for solving; formula under TM polarization mode(20) In (1)Is unknown, the primary electric field is used in this embodimentTo be substituted and converted to equation (21) to be solved.
For the formula (16) and the formula (21) obtained after substitution, each unit is analyzed by using a Galerkin method, an algebraic equation system taking an electromagnetic field of a space wave number domain on a node as an unknown quantity is formed, a first class boundary condition is imposed, and a linear equation system with the bandwidth of 5 and the diagonal dominance can be obtained.
Step S5: solving linear equation sets with different wave numbers by adopting a catch-up method, performing inverse Fourier transform on the solution, and simultaneously performing inverse Fourier transform on the obtained solutionAnd correcting, specifically:
according to the characteristics of the linear equation set obtained in the step S4, efficient solution is carried out by adopting a catch-up method, inverse Fourier transform is carried out on the solution, and a secondary electric field of a space domain is obtained in a TE polarization modeSuperimposing a primary electric fieldThe total electric field can be obtained(ii) a Obtaining secondary magnetic field of space domain in TM polarization modeSuperimposing a primary magnetic fieldObtain the total magnetic field. For TM polarization mode, the total magnetic field is obtainedThe last two equations in the formula (3) can be solved to obtain。
Since the iteration begins with a primary electric fieldInstead of in equation (15)Or, using a primary electric fieldInstead of in equation (20)For this problem, the present embodiment adopts a new method for iteratively calculating an electromagnetic field, and a specific iteration format thereof is as follows formula (22):
in the formula,is as followsThe electric field obtained after correction is carried out in the secondary iteration;is as followsElectric field obtained by sub-iteration without correction, in TE polarization modeIs the total electric fieldIn TM polarization modeIs the total electric fieldAndand is andandrespectively adopting a formula (22) to carry out iterative updating; further, in the above-mentioned case, ,is related to the primary field conductivity(i.e., background conductivity), secondary field conductivity(i.e., abnormal conductivity) related tensors.
In the TE polarization mode, the polarization direction of the polarization direction,and is made of(ii) a In the TM polarization mode, the polarization of the light beam,and is and,are respectively asUnit vector in direction.
According to the formula (22) to the calculated electric fieldMake correction update (i.e. in TE modePerforming correction update, respectively in TM modeAndperforming a correction update) to obtain a new electric field value。
Step S6: judging whether a convergence condition is reached or not according to the relative residual errors of the iteration results of the previous iteration and the next iteration, and if the convergence condition is not reached, repeating the step S4 and the step S5;
if convergent, in TE polarization mode, with total electric field satisfying convergence conditionsSubstituting into equation (15) to obtain the secondary electric fieldThen applying a secondary electric fieldSuperimposing the primary electric field obtained in step S3To obtain the final total electric field(ii) a For TM polarization mode, satisfying convergence conditionsAndsubstituting into equation (20) to obtain a secondary magnetic fieldA secondary magnetic fieldSuperimposing the primary magnetic field obtained in step S3To obtain the final total magnetic field。
Specifically, the convergence determination condition is: relative residual error of two iterationsAnd then, the iteration is stopped,for error limitation, the present embodiment is configured as。
Step S7: calculating apparent resistivity, impedance and phase on a corresponding measuring point, specifically:
when the final total electric field is obtained by calculationOr the final total magnetic fieldThen, the partial derivative in the depth direction can be obtained by using a numerical method, which is in the TE polarization modeIn the TM polarization mode is。
In the TE polarization mode:
in the TM polarization mode:
in the formula,respectively corresponding impedance, apparent resistivity and phase under the TE polarization mode;respectively the corresponding impedance, apparent resistivity and electrical resistivity in TM polarization mode,A phase;、respectively, imaginary part and real part.
The embodiment also provides a verification case of the rapid numerical simulation method of the two-dimensional magnetotelluric field:
in order to verify the correctness of the method of this embodiment, a commmi-2D 1 international standard model (which is a low-resistance abnormal body model) as shown in fig. 3 is designed, and the details thereof are as follows: has a rim in the underground mediumxA thick plate body extending infinitely in the axial direction and having a cross-sectional area ofBuried deep into(ii) a Background resistivity in the earth isAbnormal volume resistivity ofThe resistivity of the air layer above is set to. Using the projection point of the center of the abnormal body on the ground as the coordinate origin,In the direction ofAnd 600 measuring points are uniformly arranged in the range.
ForThe mesh division is explained as follows: if it is adoptedThe grid(s) is (are) used to subdivide the study area, and then the simulation area is (are):in a direction of,In the direction of(ii) a If it is adoptedThe grid(s) is (are) used to subdivide the study area, and then the simulation area is (are):in the direction of,In a direction of(ii) a If it is adoptedThe grid(s) is (are) used to subdivide the study area, and then the simulation area is (are):in the direction of,In the direction of. The three mesh generation schemes can effectively simulate the low-resistance abnormal body model shown in the figure 3.
Firstly, the correctness of the method of the embodiment is verified, the COMMEM-2D 1 international standard model is calculated by adopting the rapid numerical simulation method of the two-dimensional earth electromagnetic field in the embodiment, and fig. 4a and 4b are respectively a comparison graph of apparent resistivity in a TE polarization mode and a TM polarization mode and international published data recorded in the ZHdanov et al 1997 document, wherein the mesh size is divided intoThe calculation frequency is 0.1. As can be seen from fig. 4a and 4b, the simulation results of the method of this embodiment are both within the range of the error bars of the international published data and are closer to the mean value, thereby verifying the correctness of the method of this embodiment in performing numerical simulation on the two-dimensional model.
Next, comparing the fast numerical simulation method of this embodiment with the conventional finite difference method, the following are specific:
TABLE 1 comparison table of traditional finite difference method and fast numerical simulation method
Table 1 is a statistical table of computation time and memory consumption of the fast numerical simulation method and the conventional finite difference method in different division grids in this embodiment, where the simulation frequency is 0.1The selected polarization mode is TE polarization.
As can be seen from table 1, under the same subdivision grid size, the computation speed of the fast numerical simulation method is several orders of magnitude faster than that of the conventional finite difference method, and the consumed memory is smaller. Meanwhile, the calculation time and the consumed memory of the rapid numerical simulation method are approximately linearly and slowly increased along with the increase of the mesh division, while the traditional finite difference method is nonlinearly increased, which shows that the larger the mesh division scale is, the more obvious the advantages of the rapid numerical simulation method are. Therefore, the method of the embodiment has an important research value for developing the fine numerical simulation of the large-scale magnetotelluric method, can effectively improve the calculation efficiency, and reduces the consumption of the memory.
The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention, and various modifications and changes may be made by those skilled in the art. Any modification, equivalent replacement, or improvement made within the spirit and principle of the present invention should be included in the protection scope of the present invention.
Claims (10)
1. A method for rapid numerical simulation of a two-dimensional earth electromagnetic field, comprising the steps of:
step A1: designing a two-dimensional earth model according to the shape, size and conductivity distribution of the research area, and arranging measuring points; mesh subdivision is carried out on the two-dimensional earth electric model, and the conductivity of each node is assigned according to the electric distribution of the underground medium, so that a conductivity distribution model of the two-dimensional medium is obtained;
step A2: calculating a primary field corresponding to a background conductivity model of the research area;
step A3: constructing a linear equation set with different wave numbers by adopting a finite element method, which specifically comprises the following steps:
replacing a total electric field in a space domain secondary field control equation by a primary electric field in a primary field, and then performing Fourier transform in the horizontal direction to obtain a secondary field control equation of a space wave number domain;
for a secondary field control equation of a space wave number field, analyzing and totally synthesizing each unit by using a Galerkin method, and applying known boundary conditions to obtain linear equation sets with different wave numbers;
step A4: solving linear equation sets with different wave numbers by adopting a catch-up method, performing inverse Fourier transform on the solution to obtain a secondary field, superposing the secondary field on a primary field to obtain a total electric field or a total magnetic field, further calculating to obtain the total electric field if the total magnetic field is obtained, and then performing iterative correction on the total electric field;
step A5: judging whether a convergence condition is reached or not according to the relative residual errors of the iteration results of the previous iteration and the next iteration, and repeating the step A3 and the step A4 if the convergence condition is not reached;
if the convergence condition is met, substituting the total electric field meeting the convergence condition into a space domain secondary field control equation to obtain a secondary field, and superposing the secondary field to the primary field to obtain a final total electric field or a final total magnetic field;
step A6: and after the final total electric field or the final total magnetic field is obtained, calculating apparent resistivity, impedance and phase on the measuring point.
2. The method for rapid numerical simulation of a two-dimensional magnetotelluric field according to claim 1, wherein in step A2, the primary electric field and the primary magnetic field in the primary field are obtained by calculation in TE polarization mode, and in TM polarization mode.
3. The method for rapid numerical simulation of a two-dimensional geoelectromagnetic field according to claim 2, wherein in said step a 3: in the TE polarization mode, a primary electric field is adopted to replace a total electric field in a space domain secondary electric field control equation, and then Fourier transformation in the horizontal direction is carried out to obtain a secondary electric field control equation of a space wave number domain; in a TM polarization mode, a primary electric field is adopted to replace a total electric field in a space domain secondary magnetic field control equation, and then Fourier transformation in the horizontal direction is carried out to obtain a secondary magnetic field control equation of a space wave number domain.
4. A method for the rapid numerical simulation of a two-dimensional geoelectromagnetic field according to claim 3, wherein in said step A4: obtaining a secondary electric field of a space domain after inverse Fourier transform in a TE polarization mode, and superposing the primary electric field to obtain a total electric field; and obtaining a secondary magnetic field of a space domain after inverse Fourier transform in a TM polarization mode, superposing the primary magnetic field to obtain a total magnetic field, and further calculating to obtain the total electric field.
5. The method for rapid numerical simulation of a two-dimensional magnetotelluric field according to claim 4, wherein in step A5: if the convergence condition is met, substituting the total electric field meeting the convergence condition into a space domain secondary electric field control equation to obtain a secondary electric field in the TE polarization mode, and superposing the secondary electric field to the primary electric field to obtain a final total electric field; and under the TM polarization mode, substituting the total electric field meeting the convergence condition into a space domain secondary magnetic field control equation to obtain a secondary magnetic field, and superposing the secondary magnetic field to the primary magnetic field to obtain a final total magnetic field.
6. The method for rapid numerical simulation of a two-dimensional magnetotelluric field according to any one of claims 3 to 5, wherein the step A3 is specifically as follows:
in the TE polarization mode, the control equation of the spatial domain secondary electric field is:
let J be-J ω μ σaExThe equation (15) is subjected to a secondary electric field control equation in which fourier transform in the horizontal direction is converted into a space wave number domain:
wherein k is a wave number,a secondary electric field in the spatial wavenumber domain; exIs the total electric field in the x-direction,is a secondary electric field of the spatial domain, σbAs background conductivity, σaIs the abnormal conductivity, j is the imaginary unit, mu is the magnetic conductivity, omega is the angular frequency;
in the TM polarization mode, the spatial domain secondary magnetic field control equation is:
order toAnd performing Fourier transform in the horizontal direction on the formula (20) to obtain a secondary magnetic field control equation in a space wave number domain:
where ρ is ρa+ρb,ρbAs background resistivity, paIs the anomalous resistivity;is a secondary magnetic field in the spatial wavenumber domain,a secondary magnetic field of a spatial domain, Ey、EzTotal electric fields along the y and z directions respectively;
primary electric field in x-direction obtained using step A2 in TE polarization modeInstead of E in the formula (15)xAnd then converted into a formula (16) for solving;
in TM polarization mode using the y, z directions obtained in step A2Primary electric field ofSubstitution of E in equation (20)y、EzAnd then converted into formula (21) to be solved;
and (3) analyzing each unit by using a Galerkin method for the formula (16) or the formula (21) obtained after substitution, forming an algebraic equation system taking the electromagnetic field of the space wave number domain on the node as an unknown quantity, and imposing a first class boundary condition to obtain a linear equation system with the bandwidth of 5 and the diagonal dominance.
7. The method for rapid numerical simulation of a two-dimensional magnetotelluric field according to claim 6, wherein step A4 is iteratively modified according to equation (22):
E(n)=αE(n)′+βE(n-1) (22),
wherein E is(n)An electric field obtained after correction is carried out on the nth iteration; e(n)′For the electric field obtained in the nth iteration without correction, E in TE polarization mode(n)′Total electric field E obtained for step A4xIn TM polarization mode E(n)′Total electric field E obtained for step A4yAnd EzAnd E isyAnd EzRespectively adopting a formula (22) to carry out iterative updating; wherein
9. The two-dimensional earth of claim 8The method for rapid numerical simulation of electromagnetic fields is characterized in that the step A6 specifically comprises the following steps: step A5 obtaining the final total electric field ExOr the final total magnetic field HxThen, the partial derivative along the depth direction is obtained by using a numerical method, and the partial derivative is obtained in the TE polarization modeIn the TM polarization mode is
In the TE polarization mode:
in the TM polarization mode:
wherein Z isTE、ρTE、φTERespectively corresponding impedance, apparent resistivity and phase under the TE polarization mode; zTM、ρTM、φTMRespectively corresponding impedance, apparent resistivity and phase under a TM polarization mode; im and Re are respectively an imaginary part and a real part; σ is the total conductivity.
10. The method for rapid numerical simulation of a two-dimensional geoelectromagnetic field according to claim 7, wherein in said step A2:
the TE polarization mode is:
the TM polarization mode is:
wherein E isx,Ey,EzTotal electric field in three directions of x, y and z, respectively, Hx,Hy,HzTotal magnetic fields in x, y and z directions respectively; σ is the total conductivity;
in the step A4, obtaining the total magnetic field H in the TM polarization modexThen, the total electric field E is obtained by further solving according to the formula (3)yAnd Ez。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210227897.6A CN114297905B (en) | 2022-03-10 | 2022-03-10 | Quick numerical simulation method of two-dimensional earth electromagnetic field |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210227897.6A CN114297905B (en) | 2022-03-10 | 2022-03-10 | Quick numerical simulation method of two-dimensional earth electromagnetic field |
Publications (2)
Publication Number | Publication Date |
---|---|
CN114297905A CN114297905A (en) | 2022-04-08 |
CN114297905B true CN114297905B (en) | 2022-06-03 |
Family
ID=80978631
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202210227897.6A Active CN114297905B (en) | 2022-03-10 | 2022-03-10 | Quick numerical simulation method of two-dimensional earth electromagnetic field |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN114297905B (en) |
Families Citing this family (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN115130341B (en) * | 2022-06-23 | 2024-04-12 | 中国人民解放军国防科技大学 | TM polarization rapid cross-correlation contrast source electromagnetic inversion method under uniform background |
Citations (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
EP2302360A2 (en) * | 2009-09-24 | 2011-03-30 | ASML Netherlands B.V. | Methods and apparatus for modeling electromagnetic scattering properties of microscopic structures and methods and apparatus for reconstruction of microscopic structures |
CN106199742A (en) * | 2016-06-29 | 2016-12-07 | 吉林大学 | A kind of Frequency-domain AEM 2.5 ties up band landform inversion method |
CN109977585A (en) * | 2019-04-04 | 2019-07-05 | 中南大学 | A kind of high-precision magnetotelluric the Forward Modeling |
CN113221393A (en) * | 2021-01-29 | 2021-08-06 | 吉林大学 | Three-dimensional magnetotelluric anisotropy inversion method based on non-structural finite element method |
CN113534270A (en) * | 2021-07-20 | 2021-10-22 | 中铁二院工程集团有限责任公司 | Semi-aviation transient electromagnetic conductivity-depth imaging method and equipment |
CN113553748A (en) * | 2021-09-22 | 2021-10-26 | 中南大学 | Three-dimensional magnetotelluric forward modeling numerical simulation method |
CN113792445A (en) * | 2021-11-15 | 2021-12-14 | 中南大学 | Three-dimensional magnetotelluric numerical simulation method based on integral equation method |
CN114065585A (en) * | 2021-11-22 | 2022-02-18 | 中南大学 | Three-dimensional electrical source numerical simulation method based on coulomb specification |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN114036745A (en) * | 2021-11-08 | 2022-02-11 | 中南大学 | Anisotropic medium three-dimensional magnetotelluric forward modeling method and system |
-
2022
- 2022-03-10 CN CN202210227897.6A patent/CN114297905B/en active Active
Patent Citations (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
EP2302360A2 (en) * | 2009-09-24 | 2011-03-30 | ASML Netherlands B.V. | Methods and apparatus for modeling electromagnetic scattering properties of microscopic structures and methods and apparatus for reconstruction of microscopic structures |
CN106199742A (en) * | 2016-06-29 | 2016-12-07 | 吉林大学 | A kind of Frequency-domain AEM 2.5 ties up band landform inversion method |
CN109977585A (en) * | 2019-04-04 | 2019-07-05 | 中南大学 | A kind of high-precision magnetotelluric the Forward Modeling |
CN113221393A (en) * | 2021-01-29 | 2021-08-06 | 吉林大学 | Three-dimensional magnetotelluric anisotropy inversion method based on non-structural finite element method |
CN113534270A (en) * | 2021-07-20 | 2021-10-22 | 中铁二院工程集团有限责任公司 | Semi-aviation transient electromagnetic conductivity-depth imaging method and equipment |
CN113553748A (en) * | 2021-09-22 | 2021-10-26 | 中南大学 | Three-dimensional magnetotelluric forward modeling numerical simulation method |
CN113792445A (en) * | 2021-11-15 | 2021-12-14 | 中南大学 | Three-dimensional magnetotelluric numerical simulation method based on integral equation method |
CN114065585A (en) * | 2021-11-22 | 2022-02-18 | 中南大学 | Three-dimensional electrical source numerical simulation method based on coulomb specification |
Non-Patent Citations (2)
Title |
---|
A hybrid solver based on the integral equation method and vector finite-element method for 3D controlled-source electromagnetic method modeling;Rongwen Guo et al.;《GEOPHYSICS》;20181231;第83卷(第5期);第1-41页 * |
基于各向异性介质的三维可控源电磁法快速正演研究;柳建新等;《中国地球科学联合学术年会 2020》;20201018;第2315-2317页 * |
Also Published As
Publication number | Publication date |
---|---|
CN114297905A (en) | 2022-04-08 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Schwarzbach et al. | Finite element based inversion for time-harmonic electromagnetic problems | |
Key et al. | A parallel goal-oriented adaptive finite element method for 2.5-D electromagnetic modelling | |
Rücker et al. | Three-dimensional modelling and inversion of dc resistivity data incorporating topography—I. Modelling | |
Fan et al. | Multidomain pseudospectral time-domain simulations of scattering by objects buried in lossy media | |
Guo et al. | An efficient multigrid solver based on a four-color cell-block Gauss-Seidel smoother for 3D magnetotelluric forward modeling | |
CN109977585B (en) | High-precision magnetotelluric forward modeling method | |
CN113553748B (en) | Three-dimensional magnetotelluric forward modeling numerical simulation method | |
CN113158527B (en) | Method for calculating frequency domain electromagnetic field based on implicit FVFD | |
Minisini et al. | Local time stepping with the discontinuous Galerkin method for wave propagation in 3D heterogeneous media | |
Pan et al. | 2.5-D and 3-D DC resistivity modelling using an extrapolation cascadic multigrid method | |
Zhang et al. | 3D inversion of time-domain electromagnetic data using finite elements and a triple mesh formulation | |
CN114297905B (en) | Quick numerical simulation method of two-dimensional earth electromagnetic field | |
Song et al. | High-frequency wavefield extrapolation using the Fourier neural operator | |
Long et al. | Three-dimensional forward modelling of gravity data using mesh-free methods with radial basis functions and unstructured nodes | |
Pardo et al. | Fourier series expansion in a non-orthogonal system of coordinates for the simulation of 3D-DC borehole resistivity measurements | |
CN105354421A (en) | Magnetotelluric meshless numerical simulation method for random conductive medium model | |
CN115292973B (en) | Arbitrarily sampled space wave number domain three-dimensional magnetic field numerical simulation method and system | |
CN111158059A (en) | Gravity inversion method based on cubic B spline function | |
Vatankhah et al. | Large-scale focusing joint inversion of gravity and magnetic data with Gramian constraint | |
Li et al. | An efficient algebraic multi-resolution sampling approach to 3-D magnetotelluric modelling | |
Han et al. | Efficient three-dimensional inversion of magnetotelluric data using approximate sensitivities | |
CN117972282A (en) | Electromagnetic propagation analysis method and system in anisotropic layered medium | |
Chen et al. | 3-D marine controlled-source electromagnetic modelling in an anisotropic medium using a Wavelet–Galerkin method with a secondary potential formulation | |
Bucha et al. | Cap integration in spectral gravity forward modelling up to the full gravity tensor | |
CN115017782A (en) | Three-dimensional natural source electromagnetic field calculation method considering medium anisotropy |
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 |