Nothing Special   »   [go: up one dir, main page]

CN114297905B - A Fast Numerical Simulation Method for Two-Dimensional Electromagnetic Field - Google Patents

A Fast Numerical Simulation Method for Two-Dimensional Electromagnetic Field Download PDF

Info

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
primary
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
Application number
CN202210227897.6A
Other languages
Chinese (zh)
Other versions
CN114297905A (en
Inventor
吉杭
郭荣文
柳建新
李伟
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Central South University
Original Assignee
Central South University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Central South University filed Critical Central South University
Priority to CN202210227897.6A priority Critical patent/CN114297905B/en
Publication of CN114297905A publication Critical patent/CN114297905A/en
Application granted granted Critical
Publication of CN114297905B publication Critical patent/CN114297905B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/30—Assessment of water resources

Landscapes

  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明提供了一种二维大地电磁场的快速数值模拟方法，该方法将有限单元法与傅里叶变换相结合，通过沿水平方向做傅里叶变换，将二维偏微分问题转换为不同波数间相互独立的一维常微分问题，并行度高，同时有效减少了计算量，提高了计算效率；为大规模条件下的大地电磁法的精细化数值模拟提供了条件；经过验证，采用本发明的方法计算的结果在Zhdanov et al.1997文献中收录的国际公开数据的误差棒范围内，且较为靠近均值，满足精度要求。

Figure 202210227897

The invention provides a fast numerical simulation method of two-dimensional magnetotelluric field. The method combines finite element method and Fourier transform, and transforms two-dimensional partial differential problem into different wave numbers by performing Fourier transform along the horizontal direction. The one-dimensional ordinary differential problem that is independent of each other has a high degree of parallelism, and at the same time effectively reduces the amount of calculation and improves the calculation efficiency; it provides conditions for the refined numerical simulation of the magnetotelluric method under large-scale conditions; after verification, the present invention is adopted The results calculated by the method are within the error bar range of the international public data included in the literature of Zhdanov et al.1997, and are relatively close to the mean value, which meets the accuracy requirements.

Figure 202210227897

Description

一种二维大地电磁场的快速数值模拟方法A Fast Numerical Simulation Method for Two-dimensional Electromagnetic Field

技术领域technical field

本发明涉及地球物理方法技术领域，具体涉及一种二维大地电磁场的快速数值模拟方法。The invention relates to the technical field of geophysical methods, in particular to a fast numerical simulation method of a two-dimensional magnetotelluric field.

背景技术Background technique

大地电磁法是一种频率域电磁勘探方法，其利用高空中垂直入射的天然电磁场（频段为10-4-104Hz）进而实现对地下地质构造、异常体进行探测的方法，具有操作简单、工作效率高、成本低廉、探测深度大等优势。目前，该方法已经被广泛应用于地质勘探、金属矿与油气资源探测以及地球动力学研究等领域，并取得了越来越明显的经济和社会效益。然而，此方法的成功开展依赖于对地下地质模型的正反演解释。The magnetotelluric method is a frequency-domain electromagnetic exploration method, which utilizes the natural electromagnetic field (frequency range is 10 -4 -10 4 Hz) of vertical incidence at high altitude to realize the detection of underground geological structures and abnormal bodies. It has the advantages of high work efficiency, low cost, and large detection depth. At present, this method has been widely used in geological exploration, metal ore and oil and gas resource exploration, and geodynamic research, and has achieved more and more obvious economic and social benefits. However, the successful development of this method relies on the forward and inverse interpretation of the subsurface geological model.

由于大地电磁正反演过程处理的数据量非常庞大，在实际的应用过程中往往需要消耗很长的计算时间和运算内存，且对计算机的要求较高，这在很大程度上限制了大地电磁方法在实际中的应用。在大地电磁的数值模拟过程中，需要求解关于电磁场的双旋度方程，但此方程不存在解析解，因此只能通过数值方法进行求解。Due to the huge amount of data processed in the magnetotelluric forward and inversion process, it often takes a long time to calculate and memory in the actual application process, and the requirements for the computer are relatively high, which limits the magnetotelluric to a large extent. application of the method in practice. In the process of numerical simulation of magnetotelluric, it is necessary to solve the double curl equation about the electromagnetic field, but there is no analytical solution to this equation, so it can only be solved by numerical method.

现有文献公开了很多种数值模拟方法在大地电磁正演中的应用，主要包括有常规的积分方程法(IE)、有限差分法(FD)、有限单元法(FE)等，这些数值模拟方法在计算精度和计算效率上各有优劣。积分方程法在计算大型异常体模型时，面临巨大的挑战；有限单元法网格剖分、系数矩阵的形成相对比较耗时，限制了其在大尺度电磁勘探中的应用；有限差分法使用的是结构化网格，很难处理起伏的地形和复杂异常体。而且这些数值模拟方法随着计算规模的增大，计算时间和运算内存呈非线性增加，极大的限制了大规模大地电磁法勘探的正反演解释。The existing literature discloses the application of many numerical simulation methods in MT forward modeling, mainly including the conventional integral equation method (IE), finite difference method (FD), finite element method (FE), etc. These numerical simulation methods There are advantages and disadvantages in calculation accuracy and calculation efficiency. The integral equation method faces huge challenges when calculating large-scale anomalous models; the finite element method is relatively time-consuming to mesh and form coefficient matrices, which limits its application in large-scale electromagnetic exploration; the finite difference method uses Structured mesh, difficult to deal with undulating terrain and complex anomalies. Moreover, these numerical simulation methods increase nonlinearly with the increase of computing scale, which greatly limits the forward and inverse interpretation of large-scale magnetotelluric exploration.

综上，有必要在保障一定精度的条件下，解决大地电磁正演数值模拟中计算效率和计算内存的问题，从而进一步实现大地电磁的快速正反演成像，使得大规模大地电磁的数据解释成为可能。To sum up, it is necessary to solve the problems of computational efficiency and computational memory in the forward numerical simulation of magnetotelluric under the condition of ensuring a certain accuracy, so as to further realize the fast forward and inversion imaging of magnetotelluric, and make large-scale magnetotelluric data interpretation become a possible.

发明内容SUMMARY OF THE INVENTION

本发明目的在于提供一种二维大地电磁场的快速数值模拟方法，旨在保障一定的精度的前提下，能够进一步的提高计算效率并减少一定的计算内存，具体技术方案如下：The purpose of the present invention is to provide a fast numerical simulation method of a two-dimensional magnetotelluric field, which can further improve the calculation efficiency and reduce a certain calculation memory under the premise of ensuring a certain accuracy. The specific technical scheme is as follows:

一种二维大地电磁场的快速数值模拟方法，包括以下步骤：A fast numerical simulation method for a two-dimensional magnetotelluric field, comprising the following steps:

步骤A1：由研究区域的形状、大小和电导率分布设计二维地电模型，并进行测点布置；将二维地电模型进行网格剖分，根据地下介质的电性分布对每个节点的电导率进行赋值，得到二维介质的电导率分布模型；Step A1: Design a two-dimensional geoelectric model based on the shape, size and conductivity distribution of the study area, and arrange the measuring points; divide the two-dimensional geoelectric model into grids, and analyze each node according to the electrical distribution of the underground medium. The conductivity is assigned to obtain the conductivity distribution model of the two-dimensional medium;

步骤A2：计算研究区域的背景电导率模型对应的一次场；Step A2: Calculate the primary field corresponding to the background conductivity model of the study area;

步骤A3：采用有限单元方法构造不同波数的线性方程组，具体是：Step A3: Use the finite element method to construct linear equations with different wave numbers, specifically:

采用一次场中一次电场替换空间域二次场控制方程中的总电场，然后进行水平方向的傅里叶变换得到空间波数域的二次场控制方程；The primary electric field in the primary field is used to replace the total electric field in the quadratic field governing equation in the space domain, and then the Fourier transform in the horizontal direction is performed to obtain the quadratic field governing equation in the spatial wavenumber domain;

对于空间波数域的二次场控制方程，利用伽辽金方法对每个单元进行分析，总体合成，并应用已知的边界条件得到不同波数的线性方程组；For the quadratic field governing equations in the spatial wavenumber domain, Galerkin's method is used to analyze each element, synthesize it as a whole, and apply the known boundary conditions to obtain linear equations with different wavenumbers;

步骤A4：采用追赶法求解不同波数的线性方程组，并对其解进行反傅里叶变换得到二次场，将二次场叠加一次场后获得总电场或总磁场，若为总磁场则进一步计算获得总电场，然后对总电场进行迭代修正；Step A4: Use the chasing method to solve linear equations with different wave numbers, and perform inverse Fourier transform on the solution to obtain a quadratic field. After superimposing the quadratic field once, the total electric field or total magnetic field is obtained. If it is a total magnetic field, further Calculate the total electric field, and then iteratively correct the total electric field;

步骤A5：根据前后两次迭代结果的相对残差判断是否达到收敛条件，若未达到收敛条件则重复步骤A3和步骤A4；Step A5: Judging whether the convergence condition is reached according to the relative residuals of the results of the two previous iterations, and if the convergence condition is not met, repeat Step A3 and Step A4;

若达到收敛条件，则用满足收敛条件的总电场代入空间域二次场控制方程中求得二次场，将其叠加至一次场获得最终总电场或最终总磁场；If the convergence conditions are met, substitute the total electric field satisfying the convergence conditions into the quadratic field control equation in the space domain to obtain the secondary field, and superimpose it into the primary field to obtain the final total electric field or the final total magnetic field;

步骤A6：获得最终总电场或最终总磁场后，计算测点上的视电阻率、阻抗和相位。Step A6: After obtaining the final total electric field or the final total magnetic field, calculate the apparent resistivity, impedance and phase on the measuring point.

以上技术方案中优选的，所述步骤A2中，若为TE极化模式则计算获得一次场中的一次电场，若为TM极化模式则计算获得一次场中的一次电场和一次磁场。Preferably in the above technical solutions, in the step A2, if the TE polarization mode is used, the primary electric field in the primary field is calculated and obtained, and if the TM polarization mode is used, the primary electric field and the primary magnetic field in the primary field are calculated and obtained.

以上技术方案中优选的，所述步骤A3中：在TE极化模式下，采用一次电场替代空间域二次电场控制方程中的总电场，然后进行水平方向的傅里叶变换得到空间波数域的二次电场控制方程；在TM极化模式下，采用一次电场替代空间域二次磁场控制方程中的总电场，然后进行水平方向的傅里叶变换得到空间波数域的二次磁场控制方程。Preferably in the above technical solutions, in the step A3: in the TE polarization mode, the primary electric field is used to replace the total electric field in the control equation of the secondary electric field in the space domain, and then the Fourier transform in the horizontal direction is performed to obtain the spatial wavenumber domain. The quadratic electric field governing equation; in the TM polarization mode, the primary electric field is used to replace the total electric field in the quadratic magnetic field governing equation in the space domain, and then the Fourier transform in the horizontal direction is performed to obtain the quadratic magnetic field governing equation in the spatial wavenumber domain.

以上技术方案中优选的，所述步骤A4中：在TE极化模式下反傅里叶变换后得到空间域的二次电场，叠加一次电场获得总电场；在TM极化模式下反傅里叶变换后得到空间域的二次磁场，叠加一次磁场得到总磁场，并进一步计算求得总电场。Preferably in the above technical solutions, in the step A4: in the TE polarization mode, the secondary electric field in the space domain is obtained after inverse Fourier transformation, and the total electric field is obtained by superimposing the primary electric field; in the TM polarization mode, the inverse Fourier transform After the transformation, the secondary magnetic field in the space domain is obtained, the primary magnetic field is superimposed to obtain the total magnetic field, and the total electric field is obtained by further calculation.

以上技术方案中优选的，所述步骤A5中：若达到收敛条件，TE极化模式下则用满足收敛条件的总电场代入空间域二次电场控制方程中求得二次电场，将其叠加至一次电场获得最终总电场；TM极化模式下则用满足收敛条件的总电场代入空间域二次磁场控制方程中求得二次磁场，将其叠加至一次磁场获得最终总磁场。Preferably in the above technical solutions, in the step A5: if the convergence condition is reached, in the TE polarization mode, the total electric field that satisfies the convergence condition is substituted into the control equation of the secondary electric field in the space domain to obtain the secondary electric field, and superimposed to The primary electric field is used to obtain the final total electric field; in the TM polarization mode, the total electric field satisfying the convergence condition is substituted into the quadratic magnetic field control equation in the space domain to obtain the secondary magnetic field, and it is superimposed to the primary magnetic field to obtain the final total magnetic field.

以上技术方案中优选的，所述步骤A3具体为：Preferably in the above technical solutions, the step A3 is specifically:

TE极化模式下，空间域二次电场的控制方程为：In the TE polarization mode, the governing equation of the quadratic electric field in the space domain is:

Figure 478070DEST_PATH_IMAGE001
（15），
Figure 478070DEST_PATH_IMAGE001
(15),

设

Figure 958730DEST_PATH_IMAGE002
，对公式（15）进行水平方向的傅里叶变换得空间波数域的二次电场控制方程：Assume
Figure 958730DEST_PATH_IMAGE002
, perform the horizontal Fourier transform of formula (15) to obtain the quadratic electric field control equation in the space wavenumber domain:

Figure 251171DEST_PATH_IMAGE003
（16），
Figure 251171DEST_PATH_IMAGE003
(16),

其中，

Figure 170586DEST_PATH_IMAGE004
为波数，
Figure 864872DEST_PATH_IMAGE005
为空间波数域的二次电场；
Figure 516433DEST_PATH_IMAGE006
为
Figure 561750DEST_PATH_IMAGE007
方向的总电场，
Figure 19276DEST_PATH_IMAGE008
为空间域的二次电场，
Figure 568069DEST_PATH_IMAGE009
为背景电导率，
Figure 718428DEST_PATH_IMAGE010
为异常电导率，
Figure 985461DEST_PATH_IMAGE011
为虚数单位，
Figure 541951DEST_PATH_IMAGE012
为磁导率，
Figure 273146DEST_PATH_IMAGE013
为角频率；in,
Figure 170586DEST_PATH_IMAGE004
is the wave number,
Figure 864872DEST_PATH_IMAGE005
is the secondary electric field in the space wavenumber domain;
Figure 516433DEST_PATH_IMAGE006
for
Figure 561750DEST_PATH_IMAGE007
The total electric field in the direction,
Figure 19276DEST_PATH_IMAGE008
is the secondary electric field in the space domain,
Figure 568069DEST_PATH_IMAGE009
is the background conductivity,
Figure 718428DEST_PATH_IMAGE010
is the abnormal conductivity,
Figure 985461DEST_PATH_IMAGE011
is an imaginary unit,
Figure 541951DEST_PATH_IMAGE012
is the magnetic permeability,
Figure 273146DEST_PATH_IMAGE013
is the angular frequency;

TM极化模式下，空间域二次磁场控制方程为：In the TM polarization mode, the governing equation of the quadratic magnetic field in the space domain is:

Figure 532089DEST_PATH_IMAGE014
（20），
Figure 532089DEST_PATH_IMAGE014
(20),

令

Figure 20840DEST_PATH_IMAGE015
，
Figure 85747DEST_PATH_IMAGE016
，并对公式（20）进行水平方向的傅里叶变换得空间波数域的二次磁场控制方程：make
Figure 20840DEST_PATH_IMAGE015
,
Figure 85747DEST_PATH_IMAGE016
, and perform the Fourier transform in the horizontal direction on formula (20) to obtain the quadratic magnetic field control equation in the space wavenumber domain:

Figure 609133DEST_PATH_IMAGE017
（21），
Figure 609133DEST_PATH_IMAGE017
(twenty one),

其中，

Figure 337180DEST_PATH_IMAGE018
，
Figure 578805DEST_PATH_IMAGE019
为背景电阻率，
Figure 181825DEST_PATH_IMAGE020
为异常电阻率；
Figure 559716DEST_PATH_IMAGE021
为空间波数域的二次磁场，
Figure 160462DEST_PATH_IMAGE022
为空间域的二次磁场，
Figure 686121DEST_PATH_IMAGE023
分别为沿y、z方向的总电场；in,
Figure 337180DEST_PATH_IMAGE018
,
Figure 578805DEST_PATH_IMAGE019
is the background resistivity,
Figure 181825DEST_PATH_IMAGE020
is abnormal resistivity;
Figure 559716DEST_PATH_IMAGE021
is the secondary magnetic field in the space wavenumber domain,
Figure 160462DEST_PATH_IMAGE022
is the secondary magnetic field in the space domain,
Figure 686121DEST_PATH_IMAGE023
are the total electric field along the y and z directions, respectively;

TE极化模式下使用步骤A2获得的沿x方向的一次电场

Figure 30515DEST_PATH_IMAGE024
替代公式（15）中的
Figure 528492DEST_PATH_IMAGE025
，并转换为公式（16）进行求解；Primary electric field along the x -direction obtained using step A2 in TE polarization mode
Figure 30515DEST_PATH_IMAGE024
Substitute Equation (15) for
Figure 528492DEST_PATH_IMAGE025
, and converted to formula (16) for solving;

TM极化模式下使用步骤A2获得的沿y、z方向的一次电场

Figure 96877DEST_PATH_IMAGE026
替代公式（20）中的
Figure 313095DEST_PATH_IMAGE027
，并转换为公式（21）进行求解；Primary electric field along the y and z directions obtained using step A2 in TM polarization mode
Figure 96877DEST_PATH_IMAGE026
Substitute the formula (20) for
Figure 313095DEST_PATH_IMAGE027
, and converted to formula (21) to solve;

对于替代后得到的公式（16）或公式（21），利用伽辽金方法对每个单元进行分析，并形成以节点上空间波数域的电磁场为未知量的代数方程组，并强加第一类边界条件，得到带宽为5且对角占优的线性方程组。For the formula (16) or formula (21) obtained after substitution, each unit is analyzed using the Galerkin method, and a system of algebraic equations is formed with the electromagnetic field in the spatial wavenumber domain on the node as the unknown, and the first type of algebraic equations are imposed. Boundary conditions, a system of linear equations with a bandwidth of 5 and diagonal dominance is obtained.

以上技术方案中优选的，所述步骤A4中根据公式（22）进行迭代修正：Preferably in the above technical solutions, in the step A4, iterative correction is performed according to formula (22):

Figure 195600DEST_PATH_IMAGE028
（22），
Figure 195600DEST_PATH_IMAGE028
(twenty two),

其中，

Figure 108936DEST_PATH_IMAGE029
为第
Figure 785905DEST_PATH_IMAGE030
次迭代进行修正后得到的电场；
Figure 223840DEST_PATH_IMAGE031
为第
Figure 972353DEST_PATH_IMAGE032
次迭代得到的未进行修正的电场，TE极化模式下
Figure 179343DEST_PATH_IMAGE031
为步骤A4获得的总电场
Figure 292793DEST_PATH_IMAGE033
，TM极化模式下
Figure 280340DEST_PATH_IMAGE034
为步骤A4获得的总电场
Figure 504648DEST_PATH_IMAGE035
和
Figure 831724DEST_PATH_IMAGE036
，且
Figure 178392DEST_PATH_IMAGE035
和
Figure 325340DEST_PATH_IMAGE037
分别采用公式（22）进行迭代更新；其中
Figure 41754DEST_PATH_IMAGE038
，
Figure 285654DEST_PATH_IMAGE039
。in,
Figure 108936DEST_PATH_IMAGE029
for the first
Figure 785905DEST_PATH_IMAGE030
The electric field obtained after the second iteration is corrected;
Figure 223840DEST_PATH_IMAGE031
for the first
Figure 972353DEST_PATH_IMAGE032
The uncorrected electric field from the second iteration, in TE polarization mode
Figure 179343DEST_PATH_IMAGE031
The total electric field obtained for step A4
Figure 292793DEST_PATH_IMAGE033
, in TM polarization mode
Figure 280340DEST_PATH_IMAGE034
The total electric field obtained for step A4
Figure 504648DEST_PATH_IMAGE035
and
Figure 831724DEST_PATH_IMAGE036
,and
Figure 178392DEST_PATH_IMAGE035
and
Figure 325340DEST_PATH_IMAGE037
Iteratively updated by formula (22) respectively; where
Figure 41754DEST_PATH_IMAGE038
,
Figure 285654DEST_PATH_IMAGE039
.

以上技术方案中优选的，所述步骤A5中，收敛的判断条件是：当两次迭代的相对残差

Figure 364075DEST_PATH_IMAGE040
时，迭代停止，
Figure 263897DEST_PATH_IMAGE041
为误差限。Preferably in the above technical solutions, in the step A5, the judgment condition for convergence is: when the relative residuals of the two iterations are
Figure 364075DEST_PATH_IMAGE040
, the iteration stops,
Figure 263897DEST_PATH_IMAGE041
is the error limit.

以上技术方案中优选的，所述步骤A6具体是：步骤A5获得最终总电场

Figure 830008DEST_PATH_IMAGE033
或最终总磁场
Figure 866097DEST_PATH_IMAGE042
后，使用数值方法得到其沿深度方向的偏导数，在TE极化模式下为
Figure 554567DEST_PATH_IMAGE043
，在TM极化模式下为
Figure 676107DEST_PATH_IMAGE044
；Preferably in the above technical solutions, the step A6 is specifically: step A5 obtains the final total electric field
Figure 830008DEST_PATH_IMAGE033
or the final total magnetic field
Figure 866097DEST_PATH_IMAGE042
Then, use the numerical method to obtain its partial derivative along the depth direction, in the TE polarization mode, it is
Figure 554567DEST_PATH_IMAGE043
, in the TM polarization mode is
Figure 676107DEST_PATH_IMAGE044
;

TE极化模式下有：In TE polarization mode, there are:

Figure 108226DEST_PATH_IMAGE045
（23），
Figure 108226DEST_PATH_IMAGE045
(twenty three),

TM极化模式下有：In TM polarization mode:

Figure 264400DEST_PATH_IMAGE046
（24），
Figure 264400DEST_PATH_IMAGE046
(twenty four),

其中，

Figure 61455DEST_PATH_IMAGE047
分别为TE极化模式下对应的阻抗、视电阻率、相位；
Figure 234073DEST_PATH_IMAGE048
分别为TM极化模式下对应的阻抗、视电阻率、相位；
Figure 266620DEST_PATH_IMAGE049
分别为虚部、实部。in,
Figure 61455DEST_PATH_IMAGE047
are the corresponding impedance, apparent resistivity, and phase in the TE polarization mode, respectively;
Figure 234073DEST_PATH_IMAGE048
are the corresponding impedance, apparent resistivity, and phase in TM polarization mode, respectively;
Figure 266620DEST_PATH_IMAGE049
are the imaginary part and the real part, respectively.

以上技术方案中优选的，所述步骤A2中：Preferably in the above technical solutions, in the step A2:

TE极化模式为：The TE polarization mode is:

Figure 277301DEST_PATH_IMAGE050
（2），
Figure 277301DEST_PATH_IMAGE050
(2),

TM极化模式为：The TM polarization modes are:

Figure 245257DEST_PATH_IMAGE051
（3），
Figure 245257DEST_PATH_IMAGE051
(3),

其中，

Figure 403706DEST_PATH_IMAGE052
分别为
Figure 115310DEST_PATH_IMAGE053
三个方向的总电场，
Figure 246077DEST_PATH_IMAGE054
分别为
Figure 384934DEST_PATH_IMAGE053
三个方向的总磁场；
Figure 263635DEST_PATH_IMAGE055
为总电导率；in,
Figure 403706DEST_PATH_IMAGE052
respectively
Figure 115310DEST_PATH_IMAGE053
The total electric field in the three directions,
Figure 246077DEST_PATH_IMAGE054
respectively
Figure 384934DEST_PATH_IMAGE053
The total magnetic field in three directions;
Figure 263635DEST_PATH_IMAGE055
is the total conductivity;

所述步骤A4中，TM极化模式下获得总磁场

Figure 778930DEST_PATH_IMAGE056
后，根据公式（3）进一步求解得到总电场
Figure 764204DEST_PATH_IMAGE035
和
Figure 136279DEST_PATH_IMAGE037
。In the step A4, the total magnetic field is obtained in the TM polarization mode
Figure 778930DEST_PATH_IMAGE056
Then, according to formula (3), the total electric field is obtained by further solving
Figure 764204DEST_PATH_IMAGE035
and
Figure 136279DEST_PATH_IMAGE037
.

应用本发明的技术方案，具有以下有益效果：Applying the technical scheme of the present invention has the following beneficial effects:

本发明提出了一种二维大地电磁场的快速数值模拟方法，该方法将有限单元法与傅里叶变换相结合，通过沿水平方向（y轴）做傅里叶变换，将二维偏微分问题转换为不同波数间相互独立的一维常微分问题，并行度高；其中，采用傅里叶变换方法实现了二维大地电磁数值模拟的降维，把复杂的二维问题转化为多个小问题，有效提高了计算效率，并减少了计算内存；而采用追赶法求解有限单元法离散后形成的定带宽方程组，可以实现高效求解。The invention proposes a fast numerical simulation method of the two-dimensional magnetotelluric field. The method combines the finite element method and the Fourier transform, and performs the Fourier transform along the horizontal direction (y-axis) to solve the two-dimensional partial differential problem. It is converted into a one-dimensional ordinary differential problem that is independent of each other between different wavenumbers, and has a high degree of parallelism. Among them, the Fourier transform method is used to realize the dimensionality reduction of the two-dimensional magnetotelluric numerical simulation, and the complex two-dimensional problem is transformed into a number of small problems. , which effectively improves the computational efficiency and reduces the computational memory; and the use of the catch-up method to solve the fixed-bandwidth equations formed after the discrete finite element method can be efficiently solved.

本发明的方法充分的利用了傅里叶变换的高效性和有限单元方法的准确性，有效提高了大地电磁法数值模拟的计算效率，并减少了计算时间，为大规模条件下的大地电磁法的精细化数值模拟提供了条件；经过验证，采用本发明的方法计算的结果在Zhdanov etal.1997文献中收录的国际公开数据的误差棒范围内，且较为靠近均值，满足精度要求。The method of the invention fully utilizes the high efficiency of the Fourier transform and the accuracy of the finite element method, effectively improves the calculation efficiency of the magnetotelluric method numerical simulation, and reduces the calculation time, and is suitable for the magnetotelluric method under large-scale conditions. The refined numerical simulation provides conditions; after verification, the results calculated by the method of the present invention are within the error bar range of the international public data included in the document of Zhdanov et al.

除了上面所描述的目的、特征和优点之外，本发明还有其它的目的、特征和优点。下面将参照图，对本发明作进一步详细的说明。In addition to the objects, features and advantages described above, the present invention has other objects, features and advantages. The present invention will be described in further detail below with reference to the drawings.

附图说明Description of drawings

构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：The accompanying drawings constituting a part of the present application are used to provide further understanding of the present invention, and the exemplary embodiments of the present invention and their descriptions are used to explain the present invention and do not constitute an improper limitation of the present invention. In the attached image:

图1是本发明二维大地电磁场的快速数值模拟方法的流程图；Fig. 1 is the flow chart of the fast numerical simulation method of two-dimensional magnetotelluric field of the present invention;

图2是本发明提供的背景电导率模型示意图；2 is a schematic diagram of a background conductivity model provided by the present invention;

图3是实施例1的验证案例中COMMEMI-2D1国际标准模型示意图；Fig. 3 is the schematic diagram of COMMEMI-2D1 international standard model in the verification case of embodiment 1;

图4a是TE极化模式下采用本发明方法与Zhdanov et al.1997的结果对比图；Fig. 4a is the result comparison diagram of adopting the method of the present invention and Zhdanov et al.1997 under the TE polarization mode;

图4b是TM极化模式下采用本发明方法与Zhdanov et al.1997的结果对比图。Figure 4b is a comparison diagram of the results obtained by using the method of the present invention and Zhdanov et al.1997 in the TM polarization mode.

具体实施方式Detailed ways

为了便于理解本发明，下面将对本发明进行更全面的描述，并给出了本发明的较佳实施例。但是，本发明可以以许多不同的形式来实现，并不限于本文所描述的实施例。相反地，提供这些实施例的目的是使对本发明的公开内容的理解更加透彻全面。In order to facilitate understanding of the present invention, the present invention will be described more fully below, and preferred embodiments of the present invention will be given. However, the present invention may be embodied in many different forms and is not limited to the embodiments described herein. Rather, these embodiments are provided so that a thorough and complete understanding of the present disclosure is provided.

除非另有定义，本文所使用的所有的技术和科学术语与属于本发明的技术领域的技术人员通常理解的含义相同。本文中在本发明的说明书中所使用的术语只是为了描述具体的实施例的目的，不是旨在于限制本发明。Unless otherwise defined, 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 terms used herein in the description of the present invention are for the purpose of describing specific embodiments only, and are not intended to limit the present invention.

实施例1：Example 1:

参见图1，本实施例提供了一种二维大地电磁场的快速数值模拟方法，具体流程如下：Referring to FIG. 1 , this embodiment provides a fast numerical simulation method for a two-dimensional electromagnetic field. The specific process is as follows:

步骤S1：构建二维介质的电导率分布模型，具体是：Step S1: build a conductivity distribution model of a two-dimensional medium, specifically:

由研究区域的形状、大小和电导率分布设计二维地电模型，并确定计算频率以及测点布置。The two-dimensional geoelectric model is designed according to the shape, size and conductivity distribution of the research area, and the calculation frequency and the arrangement of measuring points are determined.

将所构建的二维地电模型沿水平方向（y轴）、深度方向（z轴）进行网格剖分，并记录沿y轴、z轴方向所剖分的网格节点数分别为

Figure 941424DEST_PATH_IMAGE057
，y轴、z轴方向相邻节点间的长度分别记为
Figure 994831DEST_PATH_IMAGE058
；Mesh the constructed two-dimensional geoelectric model along the horizontal direction ( y -axis) and depth direction ( z -axis), and record the number of mesh nodes along the y -axis and z -axis, respectively:
Figure 941424DEST_PATH_IMAGE057
, the lengths between adjacent nodes in the y -axis and z -axis directions are respectively recorded as
Figure 994831DEST_PATH_IMAGE058
;

根据地下介质的电性分布，分别对每个节点的电导率进行赋值，进而完成二维介质的电导率分布模型。According to the electrical distribution of the underground medium, the electrical conductivity of each node is assigned a value, and then the electrical conductivity distribution model of the two-dimensional medium is completed.

步骤S2：选择需要计算的极化模式，具体为选择TE极化模式或TM极化模式；Step S2: selecting the polarization mode to be calculated, specifically selecting the TE polarization mode or the TM polarization mode;

对于大地电磁法，地下介质中的传导电流远远大于位移电流，所以忽略位移电流，Maxwell方程组（麦克斯韦方程组）可化为如下：For the magnetotelluric method, the conduction current in the underground medium is much larger than the displacement current, so ignoring the displacement current, Maxwell's equations (Maxwell's equations) can be transformed into the following:

Figure 896928DEST_PATH_IMAGE059
（1），
Figure 896928DEST_PATH_IMAGE059
(1),

其中：

Figure 377588DEST_PATH_IMAGE060
、
Figure 997925DEST_PATH_IMAGE061
分别表示电场和磁场，
Figure 589443DEST_PATH_IMAGE062
为磁导率，
Figure 113091DEST_PATH_IMAGE063
为角频率，且角频率与计算频率
Figure 826969DEST_PATH_IMAGE064
的关系为
Figure 606706DEST_PATH_IMAGE065
，
Figure 710839DEST_PATH_IMAGE066
为总电导率，
Figure 525211DEST_PATH_IMAGE067
为虚数单位。in:
Figure 377588DEST_PATH_IMAGE060
,
Figure 997925DEST_PATH_IMAGE061
represent the electric and magnetic fields, respectively,
Figure 589443DEST_PATH_IMAGE062
is the magnetic permeability,
Figure 113091DEST_PATH_IMAGE063
is the angular frequency, and the angular frequency is the same as the calculated frequency
Figure 826969DEST_PATH_IMAGE064
relationship is
Figure 606706DEST_PATH_IMAGE065
,
Figure 710839DEST_PATH_IMAGE066
is the total conductivity,
Figure 525211DEST_PATH_IMAGE067
is an imaginary unit.

对于二维的大地电磁问题，公式（1）可以解耦成两个独立的模式，即为TE和TM极化模式。For the two-dimensional magnetotelluric problem, equation (1) can be decoupled into two independent modes, namely TE and TM polarization modes.

TE极化模式为：The TE polarization mode is:

Figure 347673DEST_PATH_IMAGE068
（2），
Figure 347673DEST_PATH_IMAGE068
(2),

TM极化模式为：The TM polarization modes are:

Figure 677023DEST_PATH_IMAGE069
（3），
Figure 677023DEST_PATH_IMAGE069
(3),

式中，

Figure 875924DEST_PATH_IMAGE070
分别为
Figure 607119DEST_PATH_IMAGE071
三个方向的总电场，
Figure 662800DEST_PATH_IMAGE072
分别为
Figure 417129DEST_PATH_IMAGE073
三个方向的总磁场；x、y、z轴相互垂直。In the formula,
Figure 875924DEST_PATH_IMAGE070
respectively
Figure 607119DEST_PATH_IMAGE071
The total electric field in the three directions,
Figure 662800DEST_PATH_IMAGE072
respectively
Figure 417129DEST_PATH_IMAGE073
Total magnetic field in three directions; x , y , z axes are perpendicular to each other.

步骤S3：计算研究区域的背景电导率模型对应的一次场，具体是：Step S3: Calculate the primary field corresponding to the background conductivity model of the research area, specifically:

背景电导率模型，即一维水平均匀层状介质模型，一般是将电导率分布模型的异常体部分的电导率赋值为周围介质的电导率，其一次场可由解析解公式求解得出。The background conductivity model, that is, the one-dimensional horizontal uniform layered medium model, generally assigns the conductivity of the abnormal body part of the conductivity distribution model as the conductivity of the surrounding medium, and its primary field can be obtained by the analytical solution formula.

参见图2，假设有一

Figure DEST_PATH_IMAGE075AA
层水平均匀层状介质模型，且第1层介质为空气层。在TE极化模式下，我们假设天然场源（位于空气层的顶面，即
Figure 793622DEST_PATH_IMAGE076
处）为
Figure 317007DEST_PATH_IMAGE077
方向，则可将一次电场和一次磁场的形式解设为：Referring to Figure 2, suppose there is a
Figure DEST_PATH_IMAGE075AA
Layer-level uniform layered media model, and the first layer of media is the air layer. In the TE polarization mode, we assume a natural field source (located on the top surface of the air layer, i.e.
Figure 793622DEST_PATH_IMAGE076
place) for
Figure 317007DEST_PATH_IMAGE077
direction, the formal solutions of the primary electric field and the primary magnetic field can be set as:

Figure DEST_PATH_IMAGE078
（4），
Figure DEST_PATH_IMAGE078
(4),

其中，

Figure 543589DEST_PATH_IMAGE079
；
Figure 785214DEST_PATH_IMAGE080
为第
Figure 325917DEST_PATH_IMAGE081
层介质
Figure 766126DEST_PATH_IMAGE082
方向对应的一次电场、
Figure 366871DEST_PATH_IMAGE083
为第
Figure 830214DEST_PATH_IMAGE084
层介质
Figure 174608DEST_PATH_IMAGE085
方向对应的一次磁场，
Figure 970787DEST_PATH_IMAGE086
为第
Figure 742434DEST_PATH_IMAGE087
层介质
Figure 693073DEST_PATH_IMAGE088
方向对应的一次磁场，
Figure 575578DEST_PATH_IMAGE089
为虚数单位，
Figure 990379DEST_PATH_IMAGE090
为第
Figure 932927DEST_PATH_IMAGE091
层介质的背景电导率，
Figure 370862DEST_PATH_IMAGE092
为测点的深度，
Figure 119375DEST_PATH_IMAGE093
为第
Figure 326365DEST_PATH_IMAGE094
层介质顶面的深度，
Figure 439815DEST_PATH_IMAGE095
，
Figure 365046DEST_PATH_IMAGE096
为第
Figure 651670DEST_PATH_IMAGE081
层介质对应的公式（4）中的系数，
Figure 978747DEST_PATH_IMAGE097
为自然常数，为数学中的一常数，是一个无限不循环小数。in,
Figure 543589DEST_PATH_IMAGE079
;
Figure 785214DEST_PATH_IMAGE080
for the first
Figure 325917DEST_PATH_IMAGE081
Layer medium
Figure 766126DEST_PATH_IMAGE082
The primary electric field corresponding to the direction,
Figure 366871DEST_PATH_IMAGE083
for the first
Figure 830214DEST_PATH_IMAGE084
Layer medium
Figure 174608DEST_PATH_IMAGE085
The primary magnetic field corresponding to the direction,
Figure 970787DEST_PATH_IMAGE086
for the first
Figure 742434DEST_PATH_IMAGE087
Layer medium
Figure 693073DEST_PATH_IMAGE088
The primary magnetic field corresponding to the direction,
Figure 575578DEST_PATH_IMAGE089
is an imaginary unit,
Figure 990379DEST_PATH_IMAGE090
for the first
Figure 932927DEST_PATH_IMAGE091
the background conductivity of the layer medium,
Figure 370862DEST_PATH_IMAGE092
is the depth of the measuring point,
Figure 119375DEST_PATH_IMAGE093
for the first
Figure 326365DEST_PATH_IMAGE094
the depth of the top surface of the layer medium,
Figure 439815DEST_PATH_IMAGE095
,
Figure 365046DEST_PATH_IMAGE096
for the first
Figure 651670DEST_PATH_IMAGE081
The coefficients in formula (4) corresponding to the layer medium,
Figure 978747DEST_PATH_IMAGE097
is a natural constant, a constant in mathematics, an infinite non-repeating decimal.

电性分界面上切向电场和磁场满足连续边界条件，在分界面

Figure 997518DEST_PATH_IMAGE098
上得：On the electrical interface, the tangential electric and magnetic fields satisfy the continuous boundary condition, and at the interface
Figure 997518DEST_PATH_IMAGE098
Got:

Figure 970897DEST_PATH_IMAGE099
（5），
Figure 970897DEST_PATH_IMAGE099
(5),

令

Figure 998896DEST_PATH_IMAGE100
和
Figure 180479DEST_PATH_IMAGE101
，由公式（5）得：make
Figure 998896DEST_PATH_IMAGE100
and
Figure 180479DEST_PATH_IMAGE101
, obtained from formula (5):

Figure 698048DEST_PATH_IMAGE102
（6），
Figure 698048DEST_PATH_IMAGE102
(6),

其中，

Figure 332291DEST_PATH_IMAGE103
为第
Figure 960719DEST_PATH_IMAGE091
层介质底面的电磁反射系数，
Figure 262387DEST_PATH_IMAGE104
为第
Figure 888540DEST_PATH_IMAGE105
层介质顶面与底面间的垂向距离，R i 为第
Figure 72397DEST_PATH_IMAGE106
层介质
Figure 442199DEST_PATH_IMAGE107
与
Figure 598373DEST_PATH_IMAGE108
间的比值系数。in,
Figure 332291DEST_PATH_IMAGE103
for the first
Figure 960719DEST_PATH_IMAGE091
The electromagnetic reflection coefficient of the bottom surface of the layer medium,
Figure 262387DEST_PATH_IMAGE104
for the first
Figure 888540DEST_PATH_IMAGE105
The vertical distance between the top surface and the bottom surface of the layer medium, R i is the first
Figure 72397DEST_PATH_IMAGE106
Layer medium
Figure 442199DEST_PATH_IMAGE107
and
Figure 598373DEST_PATH_IMAGE108
ratio coefficient between.

Figure 395428DEST_PATH_IMAGE109
时电场和磁场有限，满足辐射边界条件，所以第N层介质的
Figure 568046DEST_PATH_IMAGE110
，由递推公式（6）从下至上可以计算出所有的R i 。
Figure 395428DEST_PATH_IMAGE109
When the electric field and magnetic field are limited, the radiation boundary conditions are satisfied , so the
Figure 568046DEST_PATH_IMAGE110
, all R i can be calculated from bottom to top by recursive formula (6).

空气层中，在模型顶面上

Figure 475959DEST_PATH_IMAGE111
可得：In the air layer, on top of the model
Figure 475959DEST_PATH_IMAGE111
Available:

Figure 486640DEST_PATH_IMAGE112
（7），
Figure 486640DEST_PATH_IMAGE112
(7),

由公式（7）以及

Figure 454596DEST_PATH_IMAGE113
得：By formula (7) and
Figure 454596DEST_PATH_IMAGE113
have to:

Figure 613045DEST_PATH_IMAGE114
（8），
Figure 613045DEST_PATH_IMAGE114
(8),

由边界条件公式（5），可得：From the boundary condition formula (5), we can get:

Figure 324649DEST_PATH_IMAGE115
（9），
Figure 324649DEST_PATH_IMAGE115
(9),

将公式（8）中的

Figure 189837DEST_PATH_IMAGE116
和
Figure 391011DEST_PATH_IMAGE117
带入公式（9），可从上至下依次求出各层的系数
Figure 974439DEST_PATH_IMAGE107
和
Figure 489734DEST_PATH_IMAGE118
，将
Figure 209428DEST_PATH_IMAGE107
和
Figure 581504DEST_PATH_IMAGE119
以及第
Figure 652228DEST_PATH_IMAGE084
层介质中对应的深度
Figure 705635DEST_PATH_IMAGE120
代入公式（4），可计算模型中TE极化模式下第
Figure 545415DEST_PATH_IMAGE121
层介质任意深度
Figure 586927DEST_PATH_IMAGE122
处的一次电场和一次磁场。Put in formula (8)
Figure 189837DEST_PATH_IMAGE116
and
Figure 391011DEST_PATH_IMAGE117
Bringing in formula (9), the coefficients of each layer can be calculated from top to bottom
Figure 974439DEST_PATH_IMAGE107
and
Figure 489734DEST_PATH_IMAGE118
,Will
Figure 209428DEST_PATH_IMAGE107
and
Figure 581504DEST_PATH_IMAGE119
and the
Figure 652228DEST_PATH_IMAGE084
Corresponding depth in layer medium
Figure 705635DEST_PATH_IMAGE120
Substitute into formula (4), the first in the TE polarization mode in the model can be calculated
Figure 545415DEST_PATH_IMAGE121
Layer medium arbitrary depth
Figure 586927DEST_PATH_IMAGE122
primary electric field and primary magnetic field.

在TM模式下，我们假设天然场源为

Figure 144947DEST_PATH_IMAGE123
方向，一次电场和一次磁场的形式解为：In TM mode, we assume that the natural field source is
Figure 144947DEST_PATH_IMAGE123
direction, the form of the primary electric field and the primary magnetic field are:

Figure 736465DEST_PATH_IMAGE124
（10），
Figure 736465DEST_PATH_IMAGE124
(10),

式中，

Figure 430752DEST_PATH_IMAGE125
为第
Figure 410209DEST_PATH_IMAGE091
层介质
Figure 189946DEST_PATH_IMAGE126
方向对应的一次电场，
Figure 585156DEST_PATH_IMAGE127
为第
Figure 399528DEST_PATH_IMAGE091
层介质
Figure 487570DEST_PATH_IMAGE128
方向对应的一次磁场，
Figure 754603DEST_PATH_IMAGE129
为第
Figure 953503DEST_PATH_IMAGE091
层介质
Figure 622382DEST_PATH_IMAGE130
方向对应的一次电场，各层的系数
Figure 474800DEST_PATH_IMAGE131
和
Figure 229130DEST_PATH_IMAGE108
以及一次场的求解方式与TE极化模式下相同，故此处不再阐述。In the formula,
Figure 430752DEST_PATH_IMAGE125
for the first
Figure 410209DEST_PATH_IMAGE091
Layer medium
Figure 189946DEST_PATH_IMAGE126
The primary electric field corresponding to the direction,
Figure 585156DEST_PATH_IMAGE127
for the first
Figure 399528DEST_PATH_IMAGE091
Layer medium
Figure 487570DEST_PATH_IMAGE128
The primary magnetic field corresponding to the direction,
Figure 754603DEST_PATH_IMAGE129
for the first
Figure 953503DEST_PATH_IMAGE091
Layer medium
Figure 622382DEST_PATH_IMAGE130
The primary electric field corresponding to the direction, the coefficient of each layer
Figure 474800DEST_PATH_IMAGE131
and
Figure 229130DEST_PATH_IMAGE108
And the solution method of the primary field is the same as that in the TE polarization mode, so it will not be described here.

步骤S4：采用有限单元方法构造不同波数的线性方程组，具体是：Step S4: Use the finite element method to construct linear equations with different wave numbers, specifically:

总场由一次场和二次场构成（一次场包括一次电场和一次磁场，二次场包括二次电场和二次磁场），其中总场对应总电导率

Figure 966141DEST_PATH_IMAGE132
，一次场对应背景电导率
Figure 53308DEST_PATH_IMAGE133
，二次场对应异常电导率
Figure 483153DEST_PATH_IMAGE134
，其中
Figure 724778DEST_PATH_IMAGE135
、
Figure 999902DEST_PATH_IMAGE133
和
Figure 705690DEST_PATH_IMAGE136
三者间的关系如式（11）：The total field consists of the primary field and the secondary field (the primary field includes the primary electric field and the primary magnetic field, and the secondary field includes the secondary electric field and the secondary magnetic field), where the total field corresponds to the total conductivity
Figure 966141DEST_PATH_IMAGE132
, the primary field corresponds to the background conductivity
Figure 53308DEST_PATH_IMAGE133
, the secondary field corresponds to anomalous conductivity
Figure 483153DEST_PATH_IMAGE134
,in
Figure 724778DEST_PATH_IMAGE135
,
Figure 999902DEST_PATH_IMAGE133
and
Figure 705690DEST_PATH_IMAGE136
The relationship between the three is as formula (11):

Figure 40856DEST_PATH_IMAGE137
（11），
Figure 40856DEST_PATH_IMAGE137
(11),

电导率分布模型的电性参数随y轴和z轴发生变化，由公式（2）可得，TE极化模式下总电场

Figure 769778DEST_PATH_IMAGE138
满足的控制方程为：The electrical parameters of the conductivity distribution model change with the y -axis and z -axis, and can be obtained from formula (2), the total electric field in the TE polarization mode
Figure 769778DEST_PATH_IMAGE138
The governing equations that are satisfied are:

Figure 848592DEST_PATH_IMAGE139
（12），
Figure 848592DEST_PATH_IMAGE139
(12),

一次电场

Figure 408886DEST_PATH_IMAGE140
满足的控制方程为：primary electric field
Figure 408886DEST_PATH_IMAGE140
The governing equations that are satisfied are:

Figure 180533DEST_PATH_IMAGE141
（13），
Figure 180533DEST_PATH_IMAGE141
(13),

公式（12）减去公式（13）可以得到：Subtracting equation (13) from equation (12) yields:

Figure 131172DEST_PATH_IMAGE142
（14），
Figure 131172DEST_PATH_IMAGE142
(14),

式中，

Figure 13677DEST_PATH_IMAGE143
为二次电场。In the formula,
Figure 13677DEST_PATH_IMAGE143
is the secondary electric field.

此处，如果直接对公式（14）进行傅里叶变换，

Figure 631740DEST_PATH_IMAGE144
这一项会由空间域的乘积关系变为波数域的褶积关系，故对其进行进一步变换可以得到空间域二次电场的控制方程为：Here, if the Fourier transform is directly performed on formula (14),
Figure 631740DEST_PATH_IMAGE144
This term will change from the product relationship in the space domain to the convolution relationship in the wavenumber domain, so by further transforming it, the governing equation of the quadratic electric field in the space domain can be obtained as:

Figure 574288DEST_PATH_IMAGE145
（15），
Figure 574288DEST_PATH_IMAGE145
(15),

式中，

Figure 12223DEST_PATH_IMAGE146
。In the formula,
Figure 12223DEST_PATH_IMAGE146
.

设

Figure 432840DEST_PATH_IMAGE147
，对公式（15）进行水平方向的傅里叶变换可得空间波数域二次电场控制方程：Assume
Figure 432840DEST_PATH_IMAGE147
, the Fourier transform of formula (15) in the horizontal direction can be obtained to obtain the quadratic electric field control equation in the space wavenumber domain:

Figure 262999DEST_PATH_IMAGE148
（16），
Figure 262999DEST_PATH_IMAGE148
(16),

式中，

Figure 376449DEST_PATH_IMAGE149
为波数，
Figure 301680DEST_PATH_IMAGE150
为空间波数域的二次电场。In the formula,
Figure 376449DEST_PATH_IMAGE149
is the wave number,
Figure 301680DEST_PATH_IMAGE150
is the quadratic electric field in the spatial wavenumber domain.

由公式（3）可得，TM极化模式下总磁场

Figure 588305DEST_PATH_IMAGE151
满足的控制方程为：From formula (3), the total magnetic field in TM polarization mode can be obtained
Figure 588305DEST_PATH_IMAGE151
The governing equations that are satisfied are:

Figure 649801DEST_PATH_IMAGE152
（17），
Figure 649801DEST_PATH_IMAGE152
(17),

其中，

Figure 934152DEST_PATH_IMAGE153
，为总电阻率。in,
Figure 934152DEST_PATH_IMAGE153
, is the total resistivity.

一次磁场

Figure 346679DEST_PATH_IMAGE154
满足的控制方程为：primary magnetic field
Figure 346679DEST_PATH_IMAGE154
The governing equations that are satisfied are:

Figure 640257DEST_PATH_IMAGE155
（18），
Figure 640257DEST_PATH_IMAGE155
(18),

公式（17）减去公式（18）可得二次磁场

Figure 821840DEST_PATH_IMAGE156
的控制方程：The secondary magnetic field can be obtained by subtracting the formula (18) from the formula (17)
Figure 821840DEST_PATH_IMAGE156
The governing equation of :

Figure 11513DEST_PATH_IMAGE157
（19），
Figure 11513DEST_PATH_IMAGE157
(19),

进一步化简可得空间域二次磁场控制方程（其中，

Figure 911336DEST_PATH_IMAGE158
，
Figure 336501DEST_PATH_IMAGE159
）：Further simplification can obtain the quadratic magnetic field governing equation in the space domain (where,
Figure 911336DEST_PATH_IMAGE158
,
Figure 336501DEST_PATH_IMAGE159
):

Figure 638169DEST_PATH_IMAGE160
（20），
Figure 638169DEST_PATH_IMAGE160
(20),

令

Figure 998743DEST_PATH_IMAGE161
，
Figure 385862DEST_PATH_IMAGE162
，并对公式（20）进行水平方向的傅里叶变换可得空间波数域二次磁场控制方程：make
Figure 998743DEST_PATH_IMAGE161
,
Figure 385862DEST_PATH_IMAGE162
, and perform the Fourier transform in the horizontal direction on the formula (20) to obtain the quadratic magnetic field control equation in the space wavenumber domain:

Figure 522708DEST_PATH_IMAGE163
（21），
Figure 522708DEST_PATH_IMAGE163
(twenty one),

其中，

Figure 413303DEST_PATH_IMAGE164
为空间波数域的二次磁场，
Figure 210358DEST_PATH_IMAGE165
，
Figure 819194DEST_PATH_IMAGE166
为背景电阻率，
Figure 851741DEST_PATH_IMAGE167
为异常电阻率。in,
Figure 413303DEST_PATH_IMAGE164
is the secondary magnetic field in the space wavenumber domain,
Figure 210358DEST_PATH_IMAGE165
,
Figure 819194DEST_PATH_IMAGE166
is the background resistivity,
Figure 851741DEST_PATH_IMAGE167
is the abnormal resistivity.

此时，TE极化模式下公式（15）中的

Figure 862422DEST_PATH_IMAGE168
是未知量，本实施例中采用一次电场
Figure 830378DEST_PATH_IMAGE169
进行替代，并转换为公式（16）进行求解；TM极化模式下公式（20）中的
Figure 988827DEST_PATH_IMAGE170
是未知量，本实施例中采用一次电场
Figure 700431DEST_PATH_IMAGE171
来进行替代，并转换为公式（21）进行求解。At this time, in the TE polarization mode, the equation (15)
Figure 862422DEST_PATH_IMAGE168
is an unknown quantity, the primary electric field is used in this embodiment
Figure 830378DEST_PATH_IMAGE169
Substitute and convert to Equation (16) for solution; in Equation (20) in TM polarization mode
Figure 988827DEST_PATH_IMAGE170
is an unknown quantity, the primary electric field is used in this embodiment
Figure 700431DEST_PATH_IMAGE171
to substitute and convert to Equation (21) to solve.

对于替代后得到的公式（16）和公式（21），利用伽辽金方法对每个单元进行分析，并形成以节点上空间波数域的电磁场为未知量的代数方程组，并强加第一类边界条件，可以得到带宽为5、对角占优的线性方程组。For the formula (16) and formula (21) obtained after substitution, each unit is analyzed using the Galerkin method, and an algebraic equation system is formed with the electromagnetic field in the spatial wavenumber domain on the node as the unknown quantity, and the first type of algebraic equations are imposed. Boundary conditions, a system of linear equations with a bandwidth of 5 and diagonal dominance can be obtained.

步骤S5：采用追赶法求解不同波数的线性方程组，并对其解进行反傅里叶变换，同时，对得到的

Figure 565619DEST_PATH_IMAGE172
进行修正，具体是：Step S5: Use the chasing method to solve the linear equations with different wave numbers, and perform inverse Fourier transform on the solutions.
Figure 565619DEST_PATH_IMAGE172
Make corrections, specifically:

根据步骤S4得到的线性方程组的特点，采用追赶法进行高效求解，并对其解进行反傅里叶变换，TE极化模式下得到空间域的二次电场

Figure 704476DEST_PATH_IMAGE173
，叠加一次电场
Figure 553484DEST_PATH_IMAGE174
即可得到总电场
Figure 803199DEST_PATH_IMAGE175
；TM极化模式下得到空间域的二次磁场
Figure 788473DEST_PATH_IMAGE176
，叠加一次磁场
Figure 98231DEST_PATH_IMAGE177
得到总磁场
Figure 526545DEST_PATH_IMAGE178
。对于TM极化模式，得到总磁场
Figure 579952DEST_PATH_IMAGE179
后可通过公式（3）中的后两式求解得到
Figure 419732DEST_PATH_IMAGE180
。According to the characteristics of the linear equations obtained in step S4, the chasing method is used to solve the problem efficiently, and the inverse Fourier transform is performed on the solution, and the secondary electric field in the space domain is obtained in the TE polarization mode.
Figure 704476DEST_PATH_IMAGE173
, superimposing an electric field
Figure 553484DEST_PATH_IMAGE174
the total electric field
Figure 803199DEST_PATH_IMAGE175
; Obtain the secondary magnetic field in the space domain in the TM polarization mode
Figure 788473DEST_PATH_IMAGE176
, superimposing a magnetic field
Figure 98231DEST_PATH_IMAGE177
get the total magnetic field
Figure 526545DEST_PATH_IMAGE178
. For the TM polarization mode, the total magnetic field is obtained
Figure 579952DEST_PATH_IMAGE179
Afterwards, it can be obtained by solving the last two equations in formula (3)
Figure 419732DEST_PATH_IMAGE180
.

由于迭代开始采用了一次电场

Figure 900392DEST_PATH_IMAGE181
代替公式（15）中的
Figure 458412DEST_PATH_IMAGE175
，或者是，采用一次电场
Figure 315510DEST_PATH_IMAGE182
代替公式（20）中的
Figure 9796DEST_PATH_IMAGE180
，无法获得精确的总场，因此直接进行求解得到的解是Born近似解，直接进行迭代很难实现稳定的收敛，针对此问题，本实施例采用了一种新的迭代计算电磁场的方法，其具体的迭代格式如下公式（22）：Since the iteration starts with an electric field
Figure 900392DEST_PATH_IMAGE181
in place of Equation (15)
Figure 458412DEST_PATH_IMAGE175
, or, using a primary electric field
Figure 315510DEST_PATH_IMAGE182
in place of Equation (20)
Figure 9796DEST_PATH_IMAGE180
, the exact total field cannot be obtained, so the solution obtained by directly solving is the Born approximate solution, and it is difficult to achieve stable convergence by direct iteration. To solve this problem, this embodiment adopts a new method for iterative calculation of the electromagnetic field. The specific iterative format is as follows (22):

Figure 661358DEST_PATH_IMAGE183
（22），
Figure 661358DEST_PATH_IMAGE183
(twenty two),

式中，

Figure 768991DEST_PATH_IMAGE184
为第
Figure 164200DEST_PATH_IMAGE185
次迭代进行修正后得到的电场；
Figure 712993DEST_PATH_IMAGE186
为第
Figure 801035DEST_PATH_IMAGE187
次迭代得到的未进行修正的电场，TE极化模式下
Figure 927123DEST_PATH_IMAGE188
为总电场
Figure 126023DEST_PATH_IMAGE189
，TM极化模式下
Figure 529322DEST_PATH_IMAGE190
为总电场
Figure 555309DEST_PATH_IMAGE191
和
Figure 309639DEST_PATH_IMAGE192
，且
Figure 46651DEST_PATH_IMAGE193
和
Figure 570036DEST_PATH_IMAGE194
分别采用公式（22）进行迭代更新；进一步的，
Figure 858935DEST_PATH_IMAGE195
，
Figure 834981DEST_PATH_IMAGE196
，是与一次场电导率
Figure 375684DEST_PATH_IMAGE197
（即背景电导率）、二次场电导率
Figure 19155DEST_PATH_IMAGE198
（即异常电导率）有关的张量。In the formula,
Figure 768991DEST_PATH_IMAGE184
for the first
Figure 164200DEST_PATH_IMAGE185
The electric field obtained after the second iteration is corrected;
Figure 712993DEST_PATH_IMAGE186
for the first
Figure 801035DEST_PATH_IMAGE187
The uncorrected electric field from the second iteration, in TE polarization mode
Figure 927123DEST_PATH_IMAGE188
is the total electric field
Figure 126023DEST_PATH_IMAGE189
, in TM polarization mode
Figure 529322DEST_PATH_IMAGE190
is the total electric field
Figure 555309DEST_PATH_IMAGE191
and
Figure 309639DEST_PATH_IMAGE192
,and
Figure 46651DEST_PATH_IMAGE193
and
Figure 570036DEST_PATH_IMAGE194
Iteratively updated by formula (22) respectively; further,
Figure 858935DEST_PATH_IMAGE195
,
Figure 834981DEST_PATH_IMAGE196
, is related to the primary field conductivity
Figure 375684DEST_PATH_IMAGE197
(i.e. background conductivity), secondary field conductivity
Figure 19155DEST_PATH_IMAGE198
(i.e. anomalous conductivity) related tensor.

在TE极化模式下，

Figure 619900DEST_PATH_IMAGE199
，且
Figure 348822DEST_PATH_IMAGE200
；TM极化模式下，
Figure 427636DEST_PATH_IMAGE201
，且
Figure 925614DEST_PATH_IMAGE202
，
Figure 493998DEST_PATH_IMAGE203
分别为
Figure 710216DEST_PATH_IMAGE204
方向上的单位向量。In TE polarization mode,
Figure 619900DEST_PATH_IMAGE199
,and
Figure 348822DEST_PATH_IMAGE200
; In TM polarization mode,
Figure 427636DEST_PATH_IMAGE201
,and
Figure 925614DEST_PATH_IMAGE202
,
Figure 493998DEST_PATH_IMAGE203
respectively
Figure 710216DEST_PATH_IMAGE204
unit vector in the direction.

根据公式（22）对计算完成后的电场

Figure 592721DEST_PATH_IMAGE205
进行修正更新（即TE模式下对
Figure 945205DEST_PATH_IMAGE206
进行修正更新，TM模式下分别对
Figure 245343DEST_PATH_IMAGE191
和
Figure 683278DEST_PATH_IMAGE207
进行修正更新），进而得到新的电场值
Figure 103895DEST_PATH_IMAGE208
。According to formula (22), the electric field after the calculation is completed
Figure 592721DEST_PATH_IMAGE205
Make a correction update (that is, in TE mode,
Figure 945205DEST_PATH_IMAGE206
Make corrections and updates, respectively in TM mode
Figure 245343DEST_PATH_IMAGE191
and
Figure 683278DEST_PATH_IMAGE207
to correct and update), and then get the new electric field value
Figure 103895DEST_PATH_IMAGE208
.

步骤S6：根据前后两次迭代结果的相对残差，判断是否达到收敛条件，若未达到收敛条件则重复步骤S4、步骤S5；Step S6: according to the relative residuals of the two iteration results before and after, determine whether the convergence condition is reached, and if the convergence condition is not met, repeat steps S4 and S5;

若收敛，TE极化模式下则用满足收敛条件的总电场

Figure 576465DEST_PATH_IMAGE206
代入公式（15）中求得二次电场
Figure 689914DEST_PATH_IMAGE209
，然后将二次电场
Figure 615145DEST_PATH_IMAGE209
叠加步骤S3获得的一次电场
Figure 839453DEST_PATH_IMAGE210
，获得最终总电场
Figure 228846DEST_PATH_IMAGE211
；TM极化模式下则用满足收敛条件的
Figure 513197DEST_PATH_IMAGE212
和
Figure 925723DEST_PATH_IMAGE194
代入公式（20）中求得二次磁场
Figure 688143DEST_PATH_IMAGE213
，将二次磁场
Figure 932043DEST_PATH_IMAGE213
叠加步骤S3获得的一次磁场
Figure 387295DEST_PATH_IMAGE214
，获得最终总磁场
Figure 287118DEST_PATH_IMAGE215
。If it converges, in the TE polarization mode, use the total electric field that satisfies the convergence condition
Figure 576465DEST_PATH_IMAGE206
Substitute into formula (15) to obtain the secondary electric field
Figure 689914DEST_PATH_IMAGE209
, and then the secondary electric field
Figure 615145DEST_PATH_IMAGE209
The primary electric field obtained in step S3 is superimposed
Figure 839453DEST_PATH_IMAGE210
, to obtain the final total electric field
Figure 228846DEST_PATH_IMAGE211
; in TM polarization mode, use the one that satisfies the convergence condition
Figure 513197DEST_PATH_IMAGE212
and
Figure 925723DEST_PATH_IMAGE194
Substitute into formula (20) to obtain the secondary magnetic field
Figure 688143DEST_PATH_IMAGE213
, the secondary magnetic field
Figure 932043DEST_PATH_IMAGE213
Superimpose the primary magnetic field obtained in step S3
Figure 387295DEST_PATH_IMAGE214
, to obtain the final total magnetic field
Figure 287118DEST_PATH_IMAGE215
.

具体的，收敛的判断条件是：当两次迭代的相对残差

Figure 853228DEST_PATH_IMAGE216
时，迭代停止，
Figure 453099DEST_PATH_IMAGE217
为误差限，本实施例设置为
Figure 79252DEST_PATH_IMAGE218
。Specifically, the judgment condition for convergence is: when the relative residuals of the two iterations are
Figure 853228DEST_PATH_IMAGE216
, the iteration stops,
Figure 453099DEST_PATH_IMAGE217
For the error limit, this example is set as
Figure 79252DEST_PATH_IMAGE218
.

步骤S7：计算对应测点上的视电阻率、阻抗和相位，具体是：Step S7: Calculate the apparent resistivity, impedance and phase on the corresponding measuring point, specifically:

当计算得到最终总电场

Figure 466371DEST_PATH_IMAGE206
或最终总磁场
Figure 429648DEST_PATH_IMAGE219
后，可以使用数值方法得到其沿深度方向的偏导数，在TE极化模式下为
Figure 585823DEST_PATH_IMAGE220
，在TM极化模式下为
Figure 382878DEST_PATH_IMAGE221
。When calculating the final total electric field
Figure 466371DEST_PATH_IMAGE206
or the final total magnetic field
Figure 429648DEST_PATH_IMAGE219
Afterwards, its partial derivative along the depth direction can be obtained numerically, which in the TE polarization mode is
Figure 585823DEST_PATH_IMAGE220
, in the TM polarization mode is
Figure 382878DEST_PATH_IMAGE221
.

TE极化模式下有：In TE polarization mode, there are:

Figure 257293DEST_PATH_IMAGE222
（23），
Figure 257293DEST_PATH_IMAGE222
(twenty three),

TM极化模式下有：In TM polarization mode:

Figure 430785DEST_PATH_IMAGE223
（24），
Figure 430785DEST_PATH_IMAGE223
(twenty four),

式中，

Figure 441467DEST_PATH_IMAGE224
分别为TE极化模式下对应的阻抗、视电阻率、相位；
Figure 471739DEST_PATH_IMAGE225
分别为TM极化模式下对应的阻抗、视电阻率、相位；
Figure 567871DEST_PATH_IMAGE226
、
Figure 13896DEST_PATH_IMAGE227
分别为虚部、实部。In the formula,
Figure 441467DEST_PATH_IMAGE224
are the corresponding impedance, apparent resistivity, and phase in the TE polarization mode, respectively;
Figure 471739DEST_PATH_IMAGE225
are the corresponding impedance, apparent resistivity, and phase in TM polarization mode, respectively;
Figure 567871DEST_PATH_IMAGE226
,
Figure 13896DEST_PATH_IMAGE227
are the imaginary part and the real part, respectively.

本实施例还提供了二维大地电磁场的快速数值模拟方法的验证案例：This embodiment also provides a verification case of the fast numerical simulation method of the two-dimensional magnetotelluric field:

为了对本实施例的方法的正确性进行验证，设计了如图3的COMMEMI-2D1国际标准模型（其为一低阻异常体模型），其详情如下：在地下介质中有一沿x轴方向无限延伸的厚板状体，其截面积为

Figure 705515DEST_PATH_IMAGE228
，埋深为
Figure 906690DEST_PATH_IMAGE229
；大地中的背景电阻率为
Figure 552435DEST_PATH_IMAGE230
，异常体电阻率为
Figure 802150DEST_PATH_IMAGE231
，设置上方的空气层电阻率为
Figure 787424DEST_PATH_IMAGE232
。将异常体的中心在地面上的投影点作为坐标原点
Figure 159499DEST_PATH_IMAGE233
，
Figure 230223DEST_PATH_IMAGE234
方向上在
Figure 283630DEST_PATH_IMAGE235
范围内均匀设置600个测点。In order to verify the correctness of the method in this embodiment, the COMMEMI-2D1 international standard model as shown in Figure 3 (which is a low-resistance anomaly body model) is designed. The details are as follows: In the underground medium, there is an infinite extension along the x axis direction A thick plate-like body with a cross-sectional area of
Figure 705515DEST_PATH_IMAGE228
, buried deep
Figure 906690DEST_PATH_IMAGE229
; the background resistivity in the earth is
Figure 552435DEST_PATH_IMAGE230
, the abnormal volume resistivity is
Figure 802150DEST_PATH_IMAGE231
, set the resistivity of the upper air layer to be
Figure 787424DEST_PATH_IMAGE232
. Use the projected point of the center of the abnormal body on the ground as the origin of coordinates
Figure 159499DEST_PATH_IMAGE233
,
Figure 230223DEST_PATH_IMAGE234
in the direction
Figure 283630DEST_PATH_IMAGE235
Set 600 measuring points evenly within the range.

对于网格的剖分，作以下的说明：若是采用

Figure 123410DEST_PATH_IMAGE236
的网格对研究区域进行剖分，则模拟区域为：
Figure 167852DEST_PATH_IMAGE234
方向为
Figure 460293DEST_PATH_IMAGE237
，
Figure 317390DEST_PATH_IMAGE238
方向为
Figure 73994DEST_PATH_IMAGE239
；若是采用
Figure DEST_PATH_IMAGE241A
的网格对研究区域进行剖分，则模拟区域为：
Figure 725555DEST_PATH_IMAGE242
方向为
Figure 833188DEST_PATH_IMAGE235
，
Figure 228398DEST_PATH_IMAGE243
方向为
Figure 777191DEST_PATH_IMAGE244
；若是采用
Figure DEST_PATH_IMAGE246A
的网格对研究区域进行剖分，则模拟区域为：
Figure 927549DEST_PATH_IMAGE234
方向为
Figure 591372DEST_PATH_IMAGE237
，
Figure 524693DEST_PATH_IMAGE247
方向为
Figure 193572DEST_PATH_IMAGE248
。上述的三种网格剖分方案均能有效模拟图3所示的低阻异常体模型。For the mesh division, make the following instructions: if using
Figure 123410DEST_PATH_IMAGE236
The grid of , divides the study area, then the simulation area is:
Figure 167852DEST_PATH_IMAGE234
direction is
Figure 460293DEST_PATH_IMAGE237
,
Figure 317390DEST_PATH_IMAGE238
direction is
Figure 73994DEST_PATH_IMAGE239
; if using
Figure DEST_PATH_IMAGE241A
The grid of , divides the study area, then the simulation area is:
Figure 725555DEST_PATH_IMAGE242
direction is
Figure 833188DEST_PATH_IMAGE235
,
Figure 228398DEST_PATH_IMAGE243
direction is
Figure 777191DEST_PATH_IMAGE244
; if using
Figure DEST_PATH_IMAGE246A
The grid of , divides the study area, then the simulation area is:
Figure 927549DEST_PATH_IMAGE234
direction is
Figure 591372DEST_PATH_IMAGE237
,
Figure 524693DEST_PATH_IMAGE247
direction is
Figure 193572DEST_PATH_IMAGE248
. The above three meshing schemes can effectively simulate the low-resistance anomaly body model shown in Figure 3.

首先对本实施例方法的正确性进行验证，采用本实施例中的二维大地电磁场的快速数值模拟方法计算COMMEMI-2D1国际标准模型，图4a和图4b分别为在TE极化模式、TM极化模式下的视电阻率与Zhdanov et al.1997文献中收录的国际公开数据的对比图，其中，剖分的网格大小为

Figure DEST_PATH_IMAGE250A
，计算频率为0.1
Figure 514832DEST_PATH_IMAGE251
。由图4a和图4b可知，本实施例方法的模拟结果均在国际公开数据的误差棒范围内，且较为靠近均值，从而验证了本实施例的方法对二维模型进行数值模拟的正确性。First, the correctness of the method in this embodiment is verified, and the COMMEMI-2D1 international standard model is calculated by using the fast numerical simulation method of the two-dimensional magnetotelluric field in this embodiment. The comparison chart of apparent resistivity under the mode and the international public data included in Zhdanov et al.1997, where the mesh size is
Figure DEST_PATH_IMAGE250A
, the calculation frequency is 0.1
Figure 514832DEST_PATH_IMAGE251
. It can be seen from Fig. 4a and Fig. 4b that the simulation results of the method of this embodiment are all within the error bar range of the international public data, and are relatively close to the mean value, which verifies the correctness of the method of this embodiment for numerical simulation of the two-dimensional model.

接下来，将本实施例的快速数值模拟方法与传统的有限差分法进行对比，具体如下：Next, the rapid numerical simulation method of this embodiment is compared with the traditional finite difference method, as follows:

表1.传统的有限差分法与快速数值模拟方法对比表Table 1. Comparison table between traditional finite difference method and fast numerical simulation method

Figure DEST_PATH_IMAGE253A
Figure DEST_PATH_IMAGE253A

表1为本实施例中快速数值模拟方法与传统的有限差分法在不同的剖分网格下的计算时间和消耗内存统计表，其中模拟的频率为0.1

Figure 65899DEST_PATH_IMAGE254
，选择的极化模式为TE极化。Table 1 is the calculation time and memory consumption statistics table of the fast numerical simulation method and the traditional finite difference method in the present embodiment under different meshes, and the frequency of the simulation is 0.1
Figure 65899DEST_PATH_IMAGE254
, and the selected polarization mode is TE polarization.

由表1可以看出，在相同的剖分网格大小的情况下，快速数值模拟方法的计算速度比传统的有限差分法快几个数量级，耗费的内存更小。同时，快速数值模拟方法的计算时间和消耗内存随着剖分网格的增大，近似呈线性规律缓慢增加，而传统的有限差分法呈非线性增加，表明网格剖分规模越大，快速数值模拟方法的优势越明显。因此，本实施例的方法对于开展大规模的大地电磁法的精细化数值模拟具有重要研究价值，可以有效的提高计算效率，并减少消耗内存。It can be seen from Table 1 that under the same mesh size, the calculation speed of the fast numerical simulation method is several orders of magnitude faster than the traditional finite difference method, and the memory consumption is smaller. At the same time, the calculation time and memory consumption of the fast numerical simulation method increase slowly with the increase of the mesh, while the traditional finite difference method increases nonlinearly, indicating that the larger the mesh size, the faster The advantages of numerical simulation methods are more obvious. Therefore, the method of this embodiment has important research value for developing a large-scale refined numerical simulation of the magnetotelluric method, which can effectively improve the computing efficiency and reduce the consumption of memory.

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。The above descriptions are only preferred embodiments of the present invention, and are not intended to limit the present invention. For those skilled in the art, the present invention may have various modifications and changes. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included within the protection scope of the present invention.

Claims (10)

1.一种二维大地电磁场的快速数值模拟方法，其特征在于，包括以下步骤：1. a fast numerical simulation method of two-dimensional magnetotelluric field, is characterized in that, comprises the following steps: 步骤A1：由研究区域的形状、大小和电导率分布设计二维地电模型，并进行测点布置；将二维地电模型进行网格剖分，根据地下介质的电性分布对每个节点的电导率进行赋值，得到二维介质的电导率分布模型；Step A1: Design a two-dimensional geoelectric model based on the shape, size and conductivity distribution of the study area, and arrange the measuring points; divide the two-dimensional geoelectric model into grids, and analyze each node according to the electrical distribution of the underground medium. The conductivity is assigned to obtain the conductivity distribution model of the two-dimensional medium; 步骤A2：计算研究区域的背景电导率模型对应的一次场；Step A2: Calculate the primary field corresponding to the background conductivity model of the study area; 步骤A3：采用有限单元方法构造不同波数的线性方程组，具体是：Step A3: Use the finite element method to construct linear equations with different wave numbers, specifically: 采用一次场中一次电场替换空间域二次场控制方程中的总电场，然后进行水平方向的傅里叶变换得到空间波数域的二次场控制方程；The primary electric field in the primary field is used to replace the total electric field in the quadratic field governing equation in the space domain, and then the Fourier transform in the horizontal direction is performed to obtain the quadratic field governing equation in the spatial wavenumber domain; 对于空间波数域的二次场控制方程，利用伽辽金方法对每个单元进行分析，总体合成，并应用已知的边界条件得到不同波数的线性方程组；For the quadratic field governing equations in the spatial wavenumber domain, Galerkin's method is used to analyze each element, synthesize it as a whole, and apply the known boundary conditions to obtain linear equations with different wavenumbers; 步骤A4：采用追赶法求解不同波数的线性方程组，并对其解进行反傅里叶变换得到二次场，将二次场叠加一次场后获得总电场或总磁场，若为总磁场则进一步计算获得总电场，然后对总电场进行迭代修正；Step A4: Use the chasing method to solve linear equations with different wave numbers, and perform inverse Fourier transform on the solution to obtain a quadratic field. After superimposing the quadratic field once, the total electric field or total magnetic field is obtained. If it is a total magnetic field, further Calculate the total electric field, and then iteratively correct the total electric field; 步骤A5：根据前后两次迭代结果的相对残差判断是否达到收敛条件，若未达到收敛条件则重复步骤A3和步骤A4；Step A5: Judging whether the convergence condition is reached according to the relative residuals of the results of the two previous iterations, and if the convergence condition is not met, repeat Step A3 and Step A4; 若达到收敛条件，则用满足收敛条件的总电场代入空间域二次场控制方程中求得二次场，将其叠加至一次场获得最终总电场或最终总磁场；If the convergence conditions are met, substitute the total electric field satisfying the convergence conditions into the quadratic field control equation in the space domain to obtain the secondary field, and superimpose it into the primary field to obtain the final total electric field or the final total magnetic field; 步骤A6：获得最终总电场或最终总磁场后，计算测点上的视电阻率、阻抗和相位。Step A6: After obtaining the final total electric field or the final total magnetic field, calculate the apparent resistivity, impedance and phase on the measuring point. 2.根据权利要求1所述的二维大地电磁场的快速数值模拟方法，其特征在于，所述步骤A2中，若为TE极化模式则计算获得一次场中的一次电场，若为TM极化模式则计算获得一次场中的一次电场和一次磁场。2. the fast numerical simulation method of two-dimensional magnetotelluric field according to claim 1, is characterized in that, in described step A2, if it is TE polarization mode, then calculate and obtain the primary electric field in the primary field, if it is TM polarization The mode calculates the primary electric field and primary magnetic field in the primary field. 3.根据权利要求2所述的二维大地电磁场的快速数值模拟方法，其特征在于，所述步骤A3中：在TE极化模式下，采用一次电场替代空间域二次电场控制方程中的总电场，然后进行水平方向的傅里叶变换得到空间波数域的二次电场控制方程；在TM极化模式下，采用一次电场替代空间域二次磁场控制方程中的总电场，然后进行水平方向的傅里叶变换得到空间波数域的二次磁场控制方程。3. The fast numerical simulation method of the two-dimensional magnetotelluric field according to claim 2, is characterized in that, in described step A3: in TE polarization mode, adopt primary electric field to replace total in the control equation of space domain quadratic electric field. electric field, and then perform the Fourier transform in the horizontal direction to obtain the quadratic electric field governing equation in the space wavenumber domain; in the TM polarization mode, the primary electric field is used to replace the total electric field in the quadratic magnetic field governing equation in the space domain, and then the horizontal direction The Fourier transform obtains the quadratic magnetic field governing equation in the spatial wavenumber domain. 4.根据权利要求3所述的二维大地电磁场的快速数值模拟方法，其特征在于，所述步骤A4中：在TE极化模式下反傅里叶变换后得到空间域的二次电场，叠加一次电场获得总电场；在TM极化模式下反傅里叶变换后得到空间域的二次磁场，叠加一次磁场得到总磁场，并进一步计算求得总电场。4. the fast numerical simulation method of two-dimensional magnetotelluric field according to claim 3, is characterized in that, in described step A4: obtain the secondary electric field of space domain after inverse Fourier transform under TE polarization mode, superposition The total electric field is obtained by the primary electric field; the secondary magnetic field in the space domain is obtained after the inverse Fourier transform in the TM polarization mode, the total magnetic field is obtained by superimposing the primary magnetic field, and the total electric field is obtained by further calculation. 5.根据权利要求4所述的二维大地电磁场的快速数值模拟方法，其特征在于，所述步骤A5中：若达到收敛条件，TE极化模式下则用满足收敛条件的总电场代入空间域二次电场控制方程中求得二次电场，将其叠加至一次电场获得最终总电场；TM极化模式下则用满足收敛条件的总电场代入空间域二次磁场控制方程中求得二次磁场，将其叠加至一次磁场获得最终总磁场。5. The fast numerical simulation method of the two-dimensional magnetotelluric field according to claim 4, it is characterized in that, in described step A5: if reach convergence condition, under TE polarization mode, then substitute the total electric field that satisfies convergence condition into space domain The secondary electric field is obtained from the quadratic electric field control equation, and it is superimposed to the primary electric field to obtain the final total electric field; in the TM polarization mode, the total electric field satisfying the convergence condition is substituted into the space domain quadratic magnetic field control equation to obtain the secondary magnetic field , and superimpose it into the primary magnetic field to obtain the final total magnetic field. 6.根据权利要求3-5任意一项所述的二维大地电磁场的快速数值模拟方法，其特征在于，所述步骤A3具体为：6. according to the fast numerical simulation method of the two-dimensional magnetotelluric field described in any one of claim 3-5, it is characterized in that, described step A3 is specifically: TE极化模式下，空间域二次电场的控制方程为：In the TE polarization mode, the governing equation of the quadratic electric field in the space domain is:
Figure FDA0003614737180000021
Figure FDA0003614737180000021
 设J＝-jωμσaEx，对公式(15)进行水平方向的傅里叶变换得空间波数域的二次电场控制方程：Suppose J=-jωμσ a E x , perform the Fourier transform in the horizontal direction on formula (15) to obtain the quadratic electric field control equation in the spatial wavenumber domain:
Figure FDA0003614737180000022
Figure FDA0003614737180000022
 其中，k为波数，
Figure FDA0003614737180000023
为空间波数域的二次电场；Ex为x方向的总电场，
Figure FDA0003614737180000024
为空间域的二次电场，σb为背景电导率，σa为异常电导率，j为虚数单位，μ为磁导率，ω为角频率；where k is the wave number,
Figure FDA0003614737180000023
is the secondary electric field in the spatial wavenumber domain; E x is the total electric field in the x direction,
Figure FDA0003614737180000024
is the secondary electric field in the space domain, σ b is the background conductivity, σ a is the abnormal conductivity, j is the imaginary unit, μ is the magnetic permeability, and ω is the angular frequency; TM极化模式下，空间域二次磁场控制方程为：In the TM polarization mode, the governing equation of the quadratic magnetic field in the space domain is:
Figure FDA0003614737180000025
Figure FDA0003614737180000025
 令
Figure FDA0003614737180000026
并对公式(20)进行水平方向的傅里叶变换得空间波数域的二次磁场控制方程：make
Figure FDA0003614737180000026
And perform the Fourier transform in the horizontal direction to formula (20) to obtain the quadratic magnetic field control equation in the space wavenumber domain:
Figure FDA0003614737180000031
Figure FDA0003614737180000031
 其中，ρ＝ρa+ρb，ρb为背景电阻率，ρa为异常电阻率；
Figure FDA0003614737180000032
为空间波数域的二次磁场，
Figure FDA0003614737180000033
为空间域的二次磁场，Ey、Ez分别为沿y、z方向的总电场；Among them, ρ=ρ a +ρ b , ρ b is the background resistivity, and ρ a is the abnormal resistivity;
Figure FDA0003614737180000032
is the secondary magnetic field in the space wavenumber domain,
Figure FDA0003614737180000033
is the secondary magnetic field in the space domain, E y and E z are the total electric field along the y and z directions, respectively; TE极化模式下使用步骤A2获得的沿x方向的一次电场
Figure FDA0003614737180000034
替代公式(15)中的Ex，并转换为公式(16)进行求解；Primary electric field along the x-direction obtained using step A2 in TE polarization mode
Figure FDA0003614737180000034
Substitute E x in formula (15) and convert to formula (16) for solving; TM极化模式下使用步骤A2获得的沿y、z方向的一次电场
Figure FDA0003614737180000035
替代公式(20)中的Ey、Ez，并转换为公式(21)进行求解；Primary electric field along the y and z directions obtained using step A2 in TM polarization mode
Figure FDA0003614737180000035
Substitute E y and E z in formula (20), and convert to formula (21) to solve; 对于替代后得到的公式(16)或公式(21)，利用伽辽金方法对每个单元进行分析，并形成以节点上空间波数域的电磁场为未知量的代数方程组，并强加第一类边界条件，得到带宽为5且对角占优的线性方程组。For formula (16) or formula (21) obtained after substitution, Galerkin method is used to analyze each element, and an algebraic equation system is formed with the electromagnetic field in the spatial wavenumber domain on the node as the unknown, and the first type of algebraic equations are imposed. Boundary conditions, a system of linear equations with a bandwidth of 5 and diagonal dominance is obtained. 7.根据权利要求6所述的二维大地电磁场的快速数值模拟方法，其特征在于，所述步骤A4中根据公式(22)进行迭代修正：7. The fast numerical simulation method of two-dimensional magnetotelluric field according to claim 6, is characterized in that, in described step A4, carry out iterative correction according to formula (22): E(n)＝αE(n)′+βE(n-1) (22)，E (n) = αE (n)′ + βE (n-1) (22), 其中，E(n)为第n次迭代进行修正后得到的电场；E(n)′为第n次迭代得到的未进行修正的电场，TE极化模式下E(n)′为步骤A4获得的总电场Ex，TM极化模式下E(n)′为步骤A4获得的总电场Ey和Ez，且Ey和Ez分别采用公式(22)进行迭代更新；其中
Figure FDA0003614737180000036
Among them, E (n) is the electric field obtained after the nth iteration after correction; E (n)′ is the uncorrected electric field obtained in the nth iteration, and E (n)′ in the TE polarization mode is obtained in step A4 The total electric field Ex of , E ( n )′ in the TM polarization mode is the total electric field E y and E z obtained in step A4, and E y and E z are iteratively updated using formula (22) respectively; where
Figure FDA0003614737180000036
 8.根据权利要求7所述的二维大地电磁场的快速数值模拟方法，其特征在于，所述步骤A5中，收敛的判断条件是：当两次迭代的相对残差
Figure FDA0003614737180000037
时，迭代停止，λ为误差限。8. The fast numerical simulation method of the two-dimensional magnetotelluric field according to claim 7, wherein in the step A5, the judgment condition for convergence is: when the relative residuals of two iterations are
Figure FDA0003614737180000037
When , the iteration stops, and λ is the error limit. 9.根据权利要求8所述的二维大地电磁场的快速数值模拟方法，其特征在于，所述步骤A6具体是：步骤A5获得最终总电场Ex或最终总磁场Hx后，使用数值方法得到其沿深度方向的偏导数，在TE极化模式下为
Figure FDA0003614737180000041
在TM极化模式下为
Figure FDA0003614737180000042
9. the fast numerical simulation method of the two-dimensional magnetotelluric field according to claim 8, is characterized in that, described step A6 is specifically: after step A5 obtains final total electric field E x or final total magnetic field H x , uses numerical method to obtain Its partial derivative along the depth direction in the TE polarization mode is
Figure FDA0003614737180000041
In TM polarization mode, it is
Figure FDA0003614737180000042
 TE极化模式下有：In TE polarization mode, there are:
Figure FDA0003614737180000043
Figure FDA0003614737180000043
 TM极化模式下有：In TM polarization mode:
Figure FDA0003614737180000044
Figure FDA0003614737180000044
 其中，ZTE、ρTE、φTE分别为TE极化模式下对应的阻抗、视电阻率、相位；ZTM、ρTM、φTM分别为TM极化模式下对应的阻抗、视电阻率、相位；Im、Re分别为虚部、实部；σ为总电导率。Among them, Z TE , ρ TE , and φ TE are the corresponding impedance, apparent resistivity, and phase in the TE polarization mode, respectively; Z TM , ρ TM , and φ TM are the corresponding impedance, apparent resistivity, and phase in the TM polarization mode, respectively. phase; Im and Re are imaginary and real parts, respectively; σ is the total conductivity. 10.根据权利要求7所述的二维大地电磁场的快速数值模拟方法，其特征在于，所述步骤A2中：10. The fast numerical simulation method of two-dimensional magnetotelluric field according to claim 7, is characterized in that, in described step A2: TE极化模式为：The TE polarization mode is:
Figure FDA0003614737180000051
Figure FDA0003614737180000051
 TM极化模式为：The TM polarization modes are:
Figure FDA0003614737180000052
Figure FDA0003614737180000052
 其中，Ex，Ey，Ez分别为x，y，z三个方向的总电场，Hx，Hy，Hz分别为x，y，z三个方向的总磁场；σ为总电导率；Among them, E x , E y , E z are the total electric field in the three directions of x, y and z respectively, H x , Hy , H z are the total magnetic field in the three directions of x, y and z respectively; σ is the total conductance Rate; 所述步骤A4中，TM极化模式下获得总磁场Hx后，根据公式(3)进一步求解得到总电场Ey和Ez。In the step A4, after the total magnetic field H x is obtained in the TM polarization mode, the total electric field E y and E z are obtained by further solving according to formula (3).
CN202210227897.6A 2022-03-10 2022-03-10 A Fast Numerical Simulation Method for Two-Dimensional Electromagnetic Field Active CN114297905B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210227897.6A CN114297905B (en) 2022-03-10 2022-03-10 A Fast Numerical Simulation Method for Two-Dimensional Electromagnetic Field

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210227897.6A CN114297905B (en) 2022-03-10 2022-03-10 A Fast Numerical Simulation Method for Two-Dimensional 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 A Fast Numerical Simulation Method for Two-Dimensional Electromagnetic Field

Country Status (1)

Country Link
CN (1) CN114297905B (en)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115130341B (en) * 2022-06-23 2024-04-12 中国人民解放军国防科技大学 TM polarization fast cross-correlation contrast source electromagnetic inversion method under uniform background

Citations (8)

* Cited by examiner, † Cited by third party
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 中南大学 A three-dimensional magnetotelluric forward modeling numerical simulation method
CN113792445A (en) * 2021-11-15 2021-12-14 中南大学 A 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)

* Cited by examiner, † Cited by third party
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

Patent Citations (8)

* Cited by examiner, † Cited by third party
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 中南大学 A three-dimensional magnetotelluric forward modeling numerical simulation method
CN113792445A (en) * 2021-11-15 2021-12-14 中南大学 A 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)

* Cited by examiner, † Cited by third party
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
CN108710153B (en) A Wavenumber Domain Method for Inverting Subsurface 3D Magnetic Distribution with Magnetic Full Tensor Gradient
JP7142968B2 (en) FULL WAVEFORM INVERSION METHOD, APPARATUS AND ELECTRONICS
Soloveichik et al. Finite-element solution to multidimensional multisource electromagnetic problems in the frequency domain using non-conforming meshes
Lu et al. 3D finite-volume time-domain modeling of geophysical electromagnetic data on unstructured grids using potentials
CN114970290B (en) Ground transient electromagnetic method parallel inversion method and system
CN109977585B (en) A high-precision magnetotelluric forward modeling method
Zhang et al. 3D inversion of time-domain electromagnetic data using finite elements and a triple mesh formulation
Xue et al. 3D pseudo‐seismic imaging of transient electromagnetic data–a feasibility study
CN104360404A (en) Magnetotelluric regularization inversion method based on different constraint conditions
CN115755199B (en) Three-dimensional electromagnetic inversion smoothing regularization method for unstructured grid
Le Bouteiller et al. A discontinuous Galerkin fast-sweeping eikonal solver for fast and accurate traveltime computation in 3D tilted anisotropic media
CN105717547A (en) Anisotropy medium magnetotelluric meshless value simulation method
CN105354421A (en) Magnetotelluric meshless numerical simulation method for random conductive medium model
Faucher et al. Adjoint-state method for Hybridizable Discontinuous Galerkin discretization, application to the inverse acoustic wave problem
CN114297905B (en) A Fast Numerical Simulation Method for Two-Dimensional Electromagnetic Field
CN111079278A (en) Three-dimensional time-domain hybridization discontinuous Galerkin method plus electromagnetic source term processing method
CN113569447A (en) A Fast Forward Modeling Method for Transient Electromagnetic 3D Based on Schur Compensation
CN114970289B (en) Three-dimensional magnetotelluric anisotropy forward modeling numerical simulation method, equipment and medium
CN110826283A (en) Preprocessor and three-dimensional finite difference electromagnetic forward modeling calculation method based on preprocessor
CN118566993B (en) Aviation electromagnetic three-dimensional random inversion method based on multi-scale analysis
CN115017782B (en) Three-dimensional Natural Source Electromagnetic Field Calculation Method Considering Medium Anisotropy
CN113627027B (en) Method and system for simulating electromagnetic field value of non-trivial anisotropic medium
CN114065577A (en) A three-dimensional forward modeling method for DC resistivity wavelet Galerkin
Zhao et al. U-net-based pseudoseismic imaging for the short-offset transient electromagnetic method
Zhou et al. A regularized magnetotelluric inversion with a minimum support gradient constraint

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