WO2021068527A1

WO2021068527A1 - Numerical simulation method for footprint-guided high-efficiency airborne electromagnetic survey

Info

Publication number: WO2021068527A1
Application number: PCT/CN2020/093495
Authority: WO
Inventors: 刘嵘; 柳卓; 柳建新; 王建新; 郭荣文
Original assignee: 中南大学
Priority date: 2019-10-12
Filing date: 2020-05-29
Publication date: 2021-04-15
Also published as: CN110598367A; ZA202108251B

Abstract

Provided is a numerical simulation method for a footprint-guided high-efficiency airborne electromagnetic survey. The method comprises: partitioning an airborne electromagnetic survey region using uniform grids; gradually increasing the number of uniform units, and determining a footprint size of an airborne electromagnetic survey observation device; calculating a scattering current existing in a footprint of a first survey station by using a truncated boundary vector finite element method, and storing a Green's function obtained from the calculation in this step; calculating an electromagnetic field response generated at a first survey point by the scattering current within the footprint obtained in the previous step, and storing a Green's function obtained from the calculation in this step; and repeating the two steps above for a subsequent survey station. In the present invention, calculation is performed for a survey region successively for footprints of individual survey stations, thereby significantly reducing a calculation region. In addition, the invention uses uniform grids for partitioning to enable a Green's function to be repeatedly used for a subsequent survey station when using the truncated boundary vector finite element method, thereby greatly improving calculation efficiency for subsequent survey stations.

Description

一种足迹引导的高效航空电磁法数值模拟方法A footprint-guided high-efficiency aviation electromagnetic method numerical simulation method

技术领域Technical field

本发明涉及航空电磁法正演技术领域，具体的涉及一种足迹引导的高效航空电磁法数值模拟方法。The invention relates to the technical field of aviation electromagnetic forward modeling, in particular to a footprint-guided high-efficiency aviation electromagnetic method numerical simulation method.

背景技术Background technique

航空电磁法观测装置使用机载平台，具有高效率的特点，探测区域往往覆盖几十平方公里到上千平方公里。传统数值模拟方法要求将全部探测区域作为计算区域，这种大区域网格计算将耗费大量的计算内存与时间。The aerial electromagnetic observation device uses an airborne platform, which is characterized by high efficiency, and the detection area often covers tens of square kilometers to thousands of square kilometers. Traditional numerical simulation methods require that the entire detection area be used as the calculation area. This large-area grid calculation will consume a lot of calculation memory and time.

截断边界矢量有限元法只要求将探测区域内异常体作为计算区域，具有效率高的特点。但是在面对大范围异常体存在时，需要进行大量的格林函数计算，使得该方法效率备受质疑。The truncated boundary vector finite element method only requires the abnormal body in the detection area as the calculation area, which has the characteristics of high efficiency. However, in the face of the existence of a large range of abnormal bodies, a large number of Green's function calculations are required, which makes the efficiency of this method questionable.

因此，开发一种计算区域适当且计算效率高的模拟方法具有重要意义。Therefore, it is of great significance to develop a simulation method with an appropriate calculation area and high calculation efficiency.

发明内容Summary of the invention

为了解决航空电磁法探测区域范围大带来的计算困难，本发明提供了一种足迹引导的高效航空电磁法数值模拟方法，该方法结合航空电磁法观测装置具有有限足迹区域(即航空电磁法足迹特征)与截断边界矢量有限元法，将整个航空电磁法探测区域按测站细分为多个子区域逐一计算测站测点电磁场，具有计算区域适当且计算效率高的特征，具体技术方案如下：In order to solve the calculation difficulties caused by the large range of the aviation electromagnetic method detection area, the present invention provides a footprint-guided high-efficiency aviation electromagnetic method numerical simulation method, which combines the aviation electromagnetic method observation device with a limited footprint area (that is, the aviation electromagnetic method footprint Features) and the truncated boundary vector finite element method, the entire aerial electromagnetic detection area is subdivided into multiple sub-areas according to the station to calculate the electromagnetic field of the station measurement point one by one, which has the characteristics of proper calculation area and high calculation efficiency. The specific technical solutions are as follows:

本发明提供一种足迹引导的高效航空电磁法数值模拟方法，包括以下步骤：The present invention provides a footprint-guided high-efficiency aviation electromagnetic method numerical simulation method, which includes the following steps:

步骤S100：使用均匀网格剖分航空电磁法探测区域；Step S100: Use a uniform grid to divide the aerial electromagnetic method to detect the area;

步骤S200：计算所使用航空电磁法观测装置的足迹区域在x、y、z三个方向的延伸尺度；所述足迹区域在x与y方向网格数量相等，且z方向网格数量为x或y方向网格数量的1/4至1/2之间；网格具体数量由均匀半空间模型测点电磁场解析解与数值解比率决定；Step S200: Calculate the extension scales of the footprint area of the aviation electromagnetic observation device used in the x, y, and z directions; the number of grids in the x and y directions of the footprint area is equal, and the number of grids in the z direction is x or Between 1/4 and 1/2 of the number of grids in the y direction; the specific number of grids is determined by the ratio between the analytical solution and the numerical solution of the electromagnetic field at the measurement point of the uniform half-space model;

步骤S300：使用截断边界矢量有限元法将第一个测站对应的足迹作为异常体，计算足迹内存在的散射电流，并将计算过程中产生的格林函数进行存储备用；Step S300: Use the truncated boundary vector finite element method to take the footprint corresponding to the first station as an abnormal body, calculate the scattered current in the footprint, and store the Green function generated during the calculation process for future use;

所述格林函数在使用截断边界矢量有限元时产生，是用来计算足迹内离散单元散射电流在截断边界矢量有限元法计算边界产生的电磁场；The Green function is generated when the truncated boundary vector finite element is used, and is used to calculate the electromagnetic field generated by the boundary of the truncated boundary vector finite element method in the calculation of the scattered current of the discrete elements in the footprint;

步骤S400：使用第一测站足迹内离散单元散射电流计算测点电磁场响应，并对计算过程中产生的格林函数进行存储备用；Step S400: Calculate the electromagnetic field response of the measuring point by using the scattered current of the discrete unit in the footprint of the first measuring station, and store the Green function generated in the calculation process for backup;

所述格林函数用来计算足迹内离散单元散射电流在测站观测点产生的电磁场；The Green's function is used to calculate the electromagnetic field generated by the scattered current of the discrete unit in the footprint at the observation point of the station;

步骤S500：对后续测站重复使用第一测站在S300与S400步骤所存储的格林函数，先按照S300计算测站足迹内散射电流，再按照S400计算站足迹内散射电流在测点电磁场。Step S500: Repeat the Green's function stored in steps S300 and S400 at the first station for subsequent stations, first calculate the scattered current in the station footprint according to S300, and then calculate the electromagnetic field of the station footprint at the measurement point according to S400.

对于上述本发明步骤S300和S400记载的格林函数可知，实质上属于两个不同用途的格林函数，因此，为了更好的区分，对于步骤S300中记载的格林函数，命名为第一格林函数，对于步骤S400中记载的格林函数，命名为第二格林函数。Regarding the Green's function recorded in steps S300 and S400 of the present invention, it can be seen that they essentially belong to two Green's functions with different purposes. Therefore, in order to better distinguish, the Green's function recorded in step S300 is named the first Green's function. The Green function recorded in step S400 is named the second Green function.

以上技术方案中优选的，所述步骤S100中网格为规则的六面体网格。Preferably in the above technical solutions, the grid in the step S100 is a regular hexahedral grid.

以上技术方案中优选的，步骤S200中具体要求足迹区域以发射源在地表投影为中心，按4:4:1比例逐步增加x、y、z方向的网格数量，直到足迹内散射电流在测点产生的电磁场与测点解析解误差小于5％。Preferably in the above technical solutions, step S200 specifically requires that the footprint area is centered on the projection of the emission source on the ground, and the number of grids in the x, y, and z directions is gradually increased at a ratio of 4:4:1 until the scattered current in the footprint is measured The error between the electromagnetic field generated by the point and the analytical solution of the measurement point is less than 5%.

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

步骤S210：计算航空电磁发射源在接收点产生的二次电磁场解析解；Step S210: Calculate the analytical solution of the secondary electromagnetic field generated by the aviation electromagnetic emission source at the receiving point;

计算航空电磁装置发射源正下方均匀半空间地下A区域的激发电流，并计算区域内激发电流在航空电磁测量装置产生的二次场数值解；此处A区域为400×400×100m ³的区域； Calculate the excitation current in the uniform half-space underground area A directly below the emission source of the aviation electromagnetic device, and calculate the numerical solution of the secondary field generated by the excitation current in the area in the aviation electromagnetic measurement device; here, the area A is an area of 400×400×100m ³ ；

步骤S220：计算得到的二次磁场数值解与得到的二次磁场解析解的相对误差；Step S220: the relative error between the calculated numerical solution of the secondary magnetic field and the obtained analytical solution of the secondary magnetic field;

步骤S230：相对误差大小判断，若磁场实部与虚部相对误差均不超过5％，则步骤S210中的A区域体积定义为该航空电磁装置足迹；Step S230: Determine the magnitude of the relative error. If the relative error between the real and imaginary parts of the magnetic field does not exceed 5%, the volume of area A in Step S210 is defined as the footprint of the aviation electromagnetic device;

若磁场实部与虚部相对误差均超过5％，则按照4:4:1比例扩展重新定义S210中A区域并计算区域内的激发电流和区域内激发电流在航空电磁测量装置产生的二次场数值解，返回步骤S220。If the relative error between the real part and the imaginary part of the magnetic field exceeds 5%, the area A in S210 will be redefined according to the 4:4:1 ratio expansion and the excitation current in the area and the secondary excitation current generated by the aeronautical electromagnetic measuring device in the area will be calculated. Field value solution, return to step S220.

以上技术方案中优选的，步骤S300中不论测站下方异常电导率规模，只需要将足迹区域作为计算目标区域，使用截断边界矢量有限元法计算足迹内存在的散射电流。Preferably in the above technical solutions, regardless of the scale of the abnormal conductivity below the station in step S300, only the footprint area needs to be used as the calculation target area, and the truncated boundary vector finite element method is used to calculate the scattered current in the footprint.

以上技术方案中优选的，所述步骤S300中截断边界矢量有限元法的使用包括：Preferably in the above technical solutions, the use of the truncated boundary vector finite element method in the step S300 includes:

步骤S310：将计算区域定义为足迹及其一个单元厚度包裹层，电场定义在计算区域剖分单元棱边中心点，使用矢量有限元理论建立电场所满足的待求解方程组；Step S310: The calculation area is defined as the footprint and an element thickness wrapper, the electric field is defined at the center point of the edge of the division element of the calculation area, and the vector finite element theory is used to establish a set of equations to be solved that the electric field satisfies;

步骤S320：计算格林函数联系足迹内单元中心电流与计算区域边界电场，将格林函数排列为矩阵形式，行数为计算区域边界棱边数，列数为足迹区域内部单元棱边数；Step S320: Calculate the Green function to connect the central current of the unit in the footprint with the boundary electric field of the calculation area, arrange the Green function in a matrix form, the number of rows is the number of edges of the calculation area, and the number of columns is the number of edges of the elements in the footprint area;

步骤S330：将步骤S320关系方程组系数矩阵存储备用；Step S330: Store the coefficient matrix of the relational equation group in step S320 for use;

步骤S340：形成计算区域边界棱边中心电场由足迹内部单元棱边中心电场关系式方程组；Step S340: forming a relational equation group for the electric field at the edge center of the calculation area boundary and the electric field at the edge center of the internal element of the footprint;

步骤S350：将步骤S310中计算区域边界电场使用步骤S320表达式替换，形成截断边界矢量有限元法控制方程组；Step S350: Replace the calculated area boundary electric field in step S310 with the expression in step S320 to form a truncated boundary vector finite element method control equation set;

步骤S360：求解方程组，得到足迹内单元棱边中心电场值，使用线性插值得到单元中心散射电流。Step S360: Solve the equations to obtain the electric field value at the center of the edge of the element in the footprint, and use linear interpolation to obtain the scattered current at the center of the element.

以上技术方案中优选的，所述步骤S400具体包括以下步骤：Preferably in the above technical solutions, the step S400 specifically includes the following steps:

步骤S410：计算联系第一测站测点与第一测站足迹内离散单元散射电流格林函数，并排列为行矩阵，矩阵列数为散射电流个数，即足迹离散单元个数；Step S410: Calculate and connect the measuring point of the first measuring station with the scattered current Green's function of the discrete unit in the footprint of the first measuring station, and arrange it into a row matrix. The number of matrix columns is the number of scattered currents, that is, the number of discrete units of the footprint;

步骤S420：将步骤S410形成的行矩阵存储备用；Step S420: store the row matrix formed in step S410 for future use;

步骤S430：将格林函数行矩阵乘以散射电流形成的列矩阵，得到测点磁场值。Step S430: Multiply the row matrix of Green's function by the column matrix formed by the scattered current to obtain the magnetic field value of the measuring point.

以上技术方案中优选的，所述步骤S500包括以下步骤：Preferably in the above technical solutions, the step S500 includes the following steps:

步骤S510：重复步骤S310，调用步骤S330存储的格林函数，重复步骤S340-S360，得到后续测站足迹内单元散射电流；Step S510: Repeat step S310, call the Green function stored in step S330, repeat steps S340-S360, and obtain the scattered current of the unit in the footprint of the subsequent station;

步骤S520：调用步骤S420存储的格林函数，重复步骤S430，得到后续测站测点磁场值。Step S520: Call the Green function stored in step S420, repeat step S430, and obtain the magnetic field value of the subsequent measuring station.

应用本发明的技术方案，效果是：Applying the technical scheme of the present invention, the effect is:

1、本发明提供的足迹引导的高效航空电磁法数值模拟方法，摒弃了传统数值模拟方法要求将所有探测区域作为计算区域，选择将探测区域按单个测站足迹逐次计算，极大的缩小了计算区域。1. The footprint-guided high-efficiency aviation electromagnetic method numerical simulation method provided by the present invention abandons the traditional numerical simulation method that requires all detection areas as calculation areas, and chooses to calculate the detection areas one by one according to the footprint of a single station, which greatly reduces the calculation area.

2、本发明提供的足迹引导的高效航空电磁法数值模拟方法，使用均匀网格剖分，所有测站足迹由相同网格组成。在后续测站使用截断边界矢量有限元法时，可以重复使用格林函数第一测站存储的格林函数，很大程度上提高了后续测站计算效率。2. The footprint-guided high-efficiency aviation electromagnetic method numerical simulation method provided by the present invention uses uniform grid division, and all station footprints are composed of the same grid. When the truncated boundary vector finite element method is used in subsequent stations, the Green function stored in the first station can be reused, which greatly improves the calculation efficiency of subsequent stations.

具体请参考根据本发明的一种足迹引导的高效航空电磁法数值模拟方法提出的各种实施例的如下描述，将使得本发明的上述和其他方面显而易见。For details, please refer to the following descriptions of various embodiments proposed by a footprint-guided high-efficiency aviation electromagnetic method numerical simulation method of the present invention, which will make the above and other aspects of the present invention obvious.

附图说明Description of the drawings

图1是本实施例中足迹引导的高效航空电磁法数值模拟方法的示意图；Fig. 1 is a schematic diagram of the numerical simulation method of the high-efficiency aviation electromagnetic method guided by the footprint in this embodiment;

图2是测站装置形式及下方均匀地下介质均匀剖分示意图(使用均匀网格剖分航空电磁法探测区域)；Figure 2 is a schematic diagram of the form of the station device and the uniform dissection of the uniform underground medium below (using a uniform grid to divide the detection area of the aerial electromagnetic method);

图3a是同轴装置形式地下均匀介质内电流实部分布图；Figure 3a is the layout of the real part of the current in the underground homogeneous medium in the form of a coaxial device;

图3b是同轴装置形式地下均匀介质内电流虚部分布图；Figure 3b is the layout of the imaginary part of the current in the underground homogeneous medium in the form of a coaxial device;

图4a是共面装置形式地下均匀介质内电流实部分布图；Figure 4a is the layout of the real part of the current in the underground homogeneous medium in the form of a coplanar device;

图4b是共面装置形式地下均匀介质内电流虚部分布图；Figure 4b is the layout of the imaginary part of the current in the underground homogeneous medium in the form of a coplanar device;

图5是地下足迹区域尺度评估过程图(实线框是第一次评估使用400×400×100m ³计算区域，虚线框是后续按照4:4:1扩展后计算区域)； Figure 5 is a diagram of the evaluation process of the underground footprint area scale (the solid line box is the 400×400×100m ³ calculation area used in the first evaluation, and the dashed box is the subsequent calculation area after expansion according to 4:4:1);

图6是不同方法计算大型水平延伸薄板异常体所需计算区域示意图(一号框是传统微分法所需计算区域，二号框是传统截断边界矢量有限元法所需计算区域，三号框是足迹引导的截断边界矢量有限元法所需计算区域)；Figure 6 is a schematic diagram of the calculation area required by different methods to calculate the large horizontally extending thin plate abnormal body (Box 1 is the calculation area required by the traditional differential method, Box 2 is the calculation area required by the traditional truncated boundary vector finite element method, and Box 3 is the calculation area required by the traditional truncated boundary vector finite element method. Footprint guided truncated boundary vector finite element method required calculation area);

图7是后续测站足迹区域示意图(第一测站计算区域如实线黑框所示，第二测站计算区域如虚线黑框所示)；Figure 7 is a schematic diagram of the footprint area of the subsequent station (the calculation area of the first station is shown by the solid black box, and the calculation area of the second station is shown by the dashed black box);

图8是锯齿复杂三维模型示意图(一号框是本实施例第一测站计算区域，二号框是本实施例第45测站计算区域，三号框是传统截断边界矢量有限元法在任意测站都需要的计算区域)；Figure 8 is a schematic diagram of a zigzag complex three-dimensional model (Box No. 1 is the calculation area of the first station in this embodiment, Box No. 2 is the calculation area of the station No. 45 in this embodiment, and Box No. 3 is the traditional truncated boundary vector finite element method. The calculation area required by the station);

图9a是使用同轴装置时本实施例方法与传统截断边界矢量有限元法在锯齿复杂三维模型的磁场数值解虚部对比图(黑点是本实施例方法计算结果虚部，黑线是传统截断边界矢量有限元法计算结果虚部)；Figure 9a is a comparison diagram of the imaginary part of the magnetic field numerical solution of the zigzag complex three-dimensional model using the method of this embodiment and the traditional truncated boundary vector finite element method when using a coaxial device (the black dot is the imaginary part of the calculation result of the method of this embodiment, and the black line is the traditional The imaginary part of the calculation result of the truncated boundary vector finite element method);

图9b是使用同轴装置时本实施例方法与传统截断边界矢量有限元法在锯齿复杂三维模型的磁场数值解实部对比图(圆圈是本实施例方法计算结果实部，虚线是传统截断边界矢量有限元法计算结果实部)；Fig. 9b is a comparison diagram of the numerical solution real part of the magnetic field of the zigzag complex three-dimensional model between the method of this embodiment and the traditional truncated boundary vector finite element method when a coaxial device is used (the circle is the real part of the calculation result of the method of this embodiment, and the dashed line is the traditional truncated boundary Real part of vector finite element method calculation results);

图10a是使用共面装置时本实施例方法与传统截断边界矢量有限元法在锯齿复杂三维模型的磁场数值解虚部对比图(黑点是本实施例方法计算结果虚部，黑线是传统截断边界矢量有限元法计算结果虚部)；Figure 10a is a comparison diagram of the imaginary part of the magnetic field numerical solution of the complex three-dimensional model of the sawtooth between the method of this embodiment and the traditional truncated boundary vector finite element method when the coplanar device is used (the black dot is the imaginary part of the calculation result of the method of this embodiment, and the black line is the traditional The imaginary part of the calculation result of the truncated boundary vector finite element method);

图10b是使用共面装置时本实施例方法与传统截断边界矢量有限元法在锯齿复杂三维模型的磁场数值解实部对比图(圆圈是本实施例方法计算结果实部，虚线是传统截断边界矢量有限元法计算结果实部)。Figure 10b is a comparison diagram of the numerical solution real part of the magnetic field of the zigzag complex three-dimensional model between the method of this embodiment and the traditional truncated boundary vector finite element method when a coplanar device is used (the circle is the real part of the calculation result of the method of this embodiment, and the dashed line is the traditional truncated boundary The real part of the calculation result of the vector finite element method).

具体实施方式Detailed ways

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

实施例：Examples:

见图1，本发明提供的一种足迹引导的高效航空电磁法数值模拟方法，包括以下步骤：As shown in Figure 1, a footprint-guided high-efficiency aviation electromagnetic method numerical simulation method provided by the present invention includes the following steps:

步骤S100：使用均匀网格剖分航空电磁法探测区域，优选的，此处使用的均匀网格为规则的六面体网格，网格x、y和z方向尺度相等，大小等于在5、10和25三个数之间选取，尺度越小精度要高。Step S100: Use a uniform grid to subdivide the airborne electromagnetic method to detect the area. Preferably, the uniform grid used here is a regular hexahedral grid with the same scale in the x, y and z directions, and the size is equal to 5, 10 and Choose between 25 three numbers, the smaller the scale, the higher the accuracy.

步骤S200：计算所使用航空电磁法观测装置的足迹区域在x、y、z三个方向的延伸尺度，具体包括如下步骤：Step S200: Calculate the extension scale of the footprint area of the aviation electromagnetic observation device used in the x, y, and z directions, which specifically includes the following steps:

第一步：计算单测站发射源在测点产生的均匀半空间解析解与全空间解析解，使用半空间解析解减去全空间解析解得到均匀地下介质内电流在测点产生的二次磁场解析解；The first step: Calculate the uniform half-space analytical solution and the full-space analytical solution generated by the emission source of a single station at the measurement point, and use the half-space analytical solution to subtract the full-space analytical solution to obtain the secondary current generated by the uniform underground medium at the measurement point Analytical solution of magnetic field;

第二步：计算航空电磁装置发射源在其正下方均匀地下A区域(400×400×100m ³区域)均匀网格单元中心位置产生的激发电流； Step 2: Calculate the excitation current generated by the emission source of the aviation electromagnetic device at the center of the uniform grid unit in the uniform underground A area (400×400×100m ^{3 area) directly below it;}

第三步：计算所定义体积区域(400×400×100m ³区域)均匀网格单元内激发电流在测量装置处产生的二次场进行叠加，得到测点二次磁场数值解； The third step: Calculate the secondary field generated by the excitation current at the measuring device in the uniform grid unit of the defined volume area (400×400×100m ³ area) to superimpose, and obtain the numerical solution of the secondary magnetic field at the measuring point;

第四步：计算第三步得到的二次磁场数值解与第一步得到的二次磁场解析解的相对误差；Step 4: Calculate the relative error between the numerical solution of the secondary magnetic field obtained in the third step and the analytical solution of the secondary magnetic field obtained in the first step;

第五步：相对误差大小判断，若磁场实部与虚部相对误差均不超过5％，则步骤第二步中的A区域体积定义为该航空电磁装置足迹；Step 5: Determine the relative error size. If the relative error between the real and imaginary parts of the magnetic field does not exceed 5%, the volume of area A in the second step of the step is defined as the footprint of the aviation electromagnetic device;

若磁场实部与虚部相对误差均超过5％，则按照4:4:1比例扩展(按照x、y与z方向)扩展均匀单元数量，重新定义第二步中A区域并计算区域内的激发电流和区域内激发电流在航空电磁测量装置产生的二次场数值解，返回步骤第三步。If the relative error between the real and imaginary parts of the magnetic field exceeds 5%, expand the number of uniform units according to the 4:4:1 ratio (according to the x, y and z directions), redefine the area A in the second step and calculate the area The numerical solution of the secondary field generated by the excitation current and the excitation current in the region in the aviation electromagnetic measuring device, return to the third step.

步骤S300：使用截断边界矢量有限元法将第一个测站对应的足迹作为目标区域，计算足迹内存在的散射电流，并将计算过程中产生的格林函数进行存储备用；Step S300: Use the truncated boundary vector finite element method to take the footprint corresponding to the first station as the target area, calculate the scattered current in the footprint, and store the Green function generated during the calculation process for use;

本实施例中截断边界矢量有限元法的使用包括以下步骤：The use of the truncated boundary vector finite element method in this embodiment includes the following steps:

步骤S400：使用测站足迹内离散单元散射电流计算测点电磁场响应，步骤如下：Step S400: Calculate the electromagnetic field response of the measuring point by using the scattered current of the discrete elements in the footprint of the measuring station. The steps are as follows:

优选的，具体S500包括以下步骤：Preferably, the specific S500 includes the following steps:

应用本实施例的方案，详情如下：The details of applying the scheme of this embodiment are as follows:

1、使用均匀网格剖分探测区域1. Use a uniform grid to divide the detection area

如图2所示，航空电磁法使用机载平台搭载观测装置，观测装置由发射线圈与接收线圈组成。当发射线圈与接收线圈为水平线圈时称为共面装置形式，当发射线圈与接收线圈为垂直线圈时称为同轴装置形式。航空电磁法距离地表高度在40m之间，航线从右至左，发射线圈源位于接收线圈前方7.5m，发射线圈发射频率为10000hz。为了计算该航空电磁法装置在探测区域各个测站的磁场数据，本实施例使用均匀网格地下介质。As shown in Figure 2, the aviation electromagnetic method uses an airborne platform to carry an observation device, which is composed of a transmitting coil and a receiving coil. When the transmitting coil and the receiving coil are horizontal coils, it is called a coplanar device form, and when the transmitting coil and the receiving coil are vertical coils, it is called a coaxial device form. The altitude of the airborne electromagnetic method is between 40m and the route from right to left. The source of the transmitting coil is located 7.5m in front of the receiving coil, and the transmitting frequency of the transmitting coil is 10000hz. In order to calculate the magnetic field data of each station in the detection area of the aerial electromagnetic method device, this embodiment uses a uniform grid underground medium.

2、计算航空电磁观测装置足迹区域2. Calculate the footprint area of the aerial electromagnetic observation device

航空电磁法发射装置所使用的垂直线圈时(同轴装置)可以当做水平磁偶极子源处理，发射装置所使用的水平线圈时(共面装置)可以当做垂直磁偶极子源处理。在使用时谐因子e ^iwt时，设置直角坐标，地表为xy面，z坐标向下。 The vertical coil (coaxial device) used in the aerial electromagnetic transmitter can be treated as a horizontal magnetic dipole source, and the horizontal coil used in the transmitter (coplanar device) can be treated as a vertical magnetic dipole source. When using the time-harmonic factor e ^iwt , set the rectangular coordinates, the ground surface is the xy plane, and the z coordinate is downward.

2.1、均匀地下初始电场求取，如下：2.1. The calculation of the uniform underground initial electric field is as follows:

空气中发射源为单位水平磁偶极子源(同轴装置)在接收点(地下网格单元中心)产生的电场表达式满足表达式1)-3)：The emission source in the air is a unit horizontal magnetic dipole source (coaxial device). The electric field expression generated at the receiving point (the center of the underground grid unit) satisfies the expression 1)-3):

E _z＝0 3)； E _z =0 3);

空气中单位垂直磁偶极子源(共面装置)在地下产生的电场表达式满足表达式4)-6)：The expression of the electric field generated by the unit vertical magnetic dipole source (coplanar device) in the air satisfies expression 4)-6):

E _z＝0 6)； E _z =0 6);

表达式1)-6)中：i是单位虚数，w是与发射源频率对应的角速度，Δx与Δy分别是接收点(地下网格中心点)与航空电磁发射源的x与y坐标差，r与ρ分别是接收点(地下网格中心点)与偶极子源之间的距离与水平投影距离，λ是空间波数，J ₀与J ₁分别是第一类0阶与1阶贝塞尔函数，z与z'分别为接收点(地下网格中心点)与发射源点的纵坐标，μ是空气磁导率，

σ是地下均匀介质电导率。 In expressions 1)-6): i is the unit imaginary number, w is the angular velocity corresponding to the frequency of the emission source, Δx and Δy are the difference between the x and y coordinates of the receiving point (the center point of the underground grid) and the aviation electromagnetic emission source, respectively, r and ρ are the distance between the receiving point (the center point of the underground grid) and the dipole source and the horizontal projection distance, λ is the spatial wave number, and J ₀ and J ₁ are the first type 0-order and 1st-order Bessels, respectively Er function, z and z'are the ordinates of the receiving point (the center point of the underground grid) and the emitting source point, respectively, μ is the air permeability,

σ is the conductivity of the underground homogeneous medium.

2.2、均匀地下发射源产生电流分布，如下：2.2. The current distribution generated by the uniform underground emission source is as follows:

图3a、图3b、图4a和图4b分别展示了水平磁偶极子源(同轴装置)与垂直磁偶极子源(共面装置)电流幅值分布的水平切片图与竖直切片图。图3a与图4a第一列展示了地下z＝20m剖面，地下水平面电流分布，电流大于10 ^-9.5A/m主要位于400×400m ²水平区域内。图3b与图4b第二列展示了源水平偶极子方向竖直切面电流分布，电流大于10 ^-9.5A/m主要在地下100m深度内。图3a、图3b、图4a和图4b表明了航空电磁法发射源引起的地下电流主要集中在发射源正下方400×400×100m ³区域内，在足迹评估过程中，首先选定该区域计算测点二次磁场数值解。 Figure 3a, Figure 3b, Figure 4a and Figure 4b respectively show the horizontal slice diagram and the vertical slice diagram of the current amplitude distribution of the horizontal magnetic dipole source (coaxial device) and the vertical magnetic dipole source (coplanar device) . The first column of Figure 3a and Figure 4a shows the underground z=20m section, and the current distribution on the ground level. Currents greater than 10 ^-9.5 A/m are mainly located in the horizontal area of ^{400×400m 2.} The second column of Fig. 3b and Fig. 4b shows the current distribution of the vertical section in the direction of the source horizontal dipole. The current is greater than 10 ^-9.5 A/m mainly in the depth of 100m underground. Figure 3a, Figure 3b, Figure 4a and Figure 4b show that the underground current caused by the aviation electromagnetic emission source is mainly concentrated in the 400×400×100m ³ area directly below the emission source. In the process of footprint assessment, this area is first selected for calculation Numerical solution of the secondary magnetic field at the measuring point.

2.3、测点二次磁场数值解与解析解求取，如下：2.3. The numerical solution and analytical solution of the secondary magnetic field at the measuring point are as follows:

本实施例选择均匀半空间模型，通过比较两种不同方法计算地下介质在测点产生的磁场来定义航空电磁法装置足迹区域大小，具体是：首先，计算得到发射源在地下选定区域引起的电流分布之后；然后，计算选定区域网格内电流在航空电磁观测点引起的二次磁场叠加得到测点二次磁场数值解；最后，得到的二次磁场数值解与解析解相比较，判断所选电流区域是否已达到所要求规模。其中判断准则为以二次磁场解析解为标准，二次磁场数值解误差小于5％，解析解由测点磁场均匀半空间解析解减去全空间解析解得到。In this embodiment, the uniform half-space model is selected, and the magnetic field generated by the underground medium at the measuring point is compared by comparing two different methods to define the size of the footprint area of the aviation electromagnetic method device. The specifics are: First, the calculation of the emission source caused by the selected underground area After the current distribution; then, calculate the secondary magnetic field superimposition caused by the current in the selected area grid at the aviation electromagnetic observation point to obtain the numerical solution of the secondary magnetic field at the measurement point; finally, the obtained secondary magnetic field numerical solution is compared with the analytical solution to determine Whether the selected current area has reached the required scale. The judgment criterion is that the secondary magnetic field analytical solution is the standard, and the numerical solution error of the secondary magnetic field is less than 5%. The analytical solution is obtained by subtracting the full-space analytical solution from the uniform half-space analytical solution of the magnetic field at the measuring point.

在求取所选区域内网格二次磁场数值解时，网格单元内部电流被视为均匀分布，在计算其在测点产生的磁场时，网格单元被当做电偶极子源。在航空电磁法同轴装置时，单位电偶极子源在全空间产生的磁场为表达式7)-9)：When calculating the numerical solution of the grid secondary magnetic field in the selected area, the internal current of the grid cell is regarded as a uniform distribution, and when calculating the magnetic field generated at the measuring point, the grid cell is regarded as an electric dipole source. In the airborne electromagnetic coaxial device, the magnetic field generated by the unit electric dipole source in the whole space is the expression 7)-9):

使用表达式1)-3)计算得到地下网格单元中心电场分布后，使用表达式7)-9)计算每个单元中心电流在测点产生的磁场，通过累加所选区域内单元在测点磁场贡献得到测点二次磁场数值解。After calculating the electric field distribution at the center of the underground grid cell using expressions 1)-3), use expressions 7)-9) to calculate the magnetic field generated by the center current of each cell at the measuring point, by accumulating the cell at the measuring point in the selected area The magnetic field contribution obtains the numerical solution of the secondary magnetic field at the measuring point.

在航空电磁法共面装置时，单位电偶极子源在全空间产生的磁场为表达式10)-12)：In the airborne electromagnetic coplanar device, the magnetic field generated by the unit electric dipole source in the whole space is the expression 10)-12):

使用表达式4)-6)计算得到地下网格单元中心电场分布后，使用表达式10)-12)计算每个单元中心电流在测点产生的磁场，通过累加所选区域内单元在测点磁场贡献得到测点二次磁场数值解。After calculating the electric field distribution in the center of the underground grid cell using expressions 4)-6), use expressions 10)-12) to calculate the magnetic field generated by the center current of each cell at the measuring point, by accumulating the cell at the measuring point in the selected area The magnetic field contribution obtains the numerical solution of the secondary magnetic field at the measuring point.

表达式7)-12)中：

i是单位虚数；σ ₀是空气电导率；w是与发射源频率对应的角速度；μ是空气磁导率；H下标表示电偶极子方向，上标表示产生的磁场方向；Δx、Δy、Δz分别是地下单元网格中心点与航空装置测量点之间的x、y、z坐标差；r是地下单元网格中心点与航空装置接收点之间的距离。 In expressions 7)-12):

i is the unit imaginary number; σ ₀ is the air conductivity; w is the angular velocity corresponding to the frequency of the emission source; μ is the air permeability; the H subscript indicates the direction of the electric dipole, and the superscript indicates the direction of the generated magnetic field; Δx, Δy , Δz are the x, y, z coordinate differences between the grid center point of the underground unit and the measurement point of the aeronautical device; r is the distance between the grid center point of the underground unit and the receiving point of the aeronautical device.

在求取所选区域内网格二次磁场解析解时，在同轴装置下即发射装置使用x方向磁偶极子源，接收装置测量x方向磁场，单位磁偶极发射源在全空间条件下在测点产生的磁场表达式为表达式13)：When obtaining the analytical solution of the grid secondary magnetic field in the selected area, under the coaxial device, that is, the transmitting device uses the x-direction magnetic dipole source, and the receiving device measures the x-direction magnetic field. The unit magnetic dipole transmitting source is in full space conditions. The following expression for the magnetic field generated at the measuring point is expression 13):

表达式13)中各参数与表达式7)-12)中各参数一致。Each parameter in expression 13) is consistent with each parameter in expression 7)-12).

在均匀半空间条件下单位磁偶极发射源在测点产生的磁场表达式为The expression of the magnetic field generated by the unit magnetic dipole emission source at the measuring point under the uniform half-space condition is

表达式14)中各参数与表达式1)-6)中各参数一致。使用表达式14)中

减去表达式15)中

可以得到航空电磁法同轴装置发射源在测量点产生的二次磁场解析解。 Each parameter in expression 14) is consistent with each parameter in expressions 1)-6). Use expression 14)

Minus expression 15)

The analytical solution of the secondary magnetic field generated by the emission source of the aerial electromagnetic method coaxial device at the measuring point can be obtained.

在共面装置下即发射装置使用z方向磁偶极子源，接收装置测量z方向磁场，单位磁偶极发射源在全空间条件下在测点产生的磁场为表达式15)：Under the coplanar device, that is, the transmitting device uses the z-direction magnetic dipole source, and the receiving device measures the z-direction magnetic field. The magnetic field generated by the unit magnetic dipole transmitting source at the measuring point under the full space condition is expression 15):

表达式15)中各参数与表达式7)-12)中各参数一致。Each parameter in expression 15) is consistent with each parameter in expression 7)-12).

在均匀半空间条件下单位磁偶极发射源在测点产生的磁场为表达式16)：Under the uniform half-space condition, the magnetic field generated by the unit magnetic dipole emission source at the measuring point is Expression 16):

表达式16)中各参数与表达式1)-6)中各参数一致。使用表达式16)中

减去表达式15)中

可以得到航空电磁法共面装置发射源在测量点产生的二次磁场解析解。 Each parameter in expression 16) is consistent with each parameter in expression 1)-6). Use expression 16)

Minus expression 15)

The analytical solution of the secondary magnetic field generated by the emission source of the airborne electromagnetic coplanar device at the measuring point can be obtained.

2.4、自适应扩展网格区域求取足迹范围，如下：2.4. Adaptively expand the grid area to obtain the footprint range, as follows:

图5展示了逐步扩大网格区域，使得测点二次磁场数值解逐渐接近解析解，求得航空电磁装置足迹区域的过程。第一次数值解计算使用400×400×100m ³区域内网格，当区域网格内电流在测点二次磁场数值解与解析解差异大于5％时，按4:4:1比例扩展网格区域。 Figure 5 shows the process of gradually expanding the grid area so that the numerical solution of the secondary magnetic field at the measuring point gradually approaches the analytical solution, and obtaining the footprint area of the aviation electromagnetic device. The first numerical solution calculation uses a 400×400×100m ³ area grid. When the difference between the numerical solution and the analytical solution of the secondary magnetic field of the current in the area grid at the measuring point is greater than 5%, the network is expanded according to the ratio of 4:4:1 Grid area.

3、使用截断边界矢量有限元法计算第一测站足迹内散射电流，如下：3. Use the truncated boundary vector finite element method to calculate the scattered current in the footprint of the first station, as follows:

本实施例基于二次电场矢量有限元法，有限元离散网格单元棱边中心电场满足表达式17)：This embodiment is based on the quadratic electric field vector finite element method, and the finite element discrete grid element edge center electric field satisfies expression 17):

其中：

是哈密顿算符，E _p是发射源在均匀地下产生的初始电场，可以由表达式1)-6)计算；E _s是地下介质电导率异常体在离散网格上产生的二次电场，σ _*是均匀地下介质电导率；σ是地下介质电导率，其他参数同表达式1)-6)中参数的含义。 among them:

Is Hamiltonian, E _p is the initial field emission sources uniformly generated in the ground, can be calculated by Expression 1) -6); E _s is a secondary electric field generated subsurface medium conductivity anomaly in the discrete grid, σ _* is the conductivity of the uniform underground medium; σ is the conductivity of the underground medium, and other parameters have the same meaning as those in expressions 1)-6).

根据Liu,R.,R.W.Guo,J.X.Liu,C.Y.Ma,and Z.W.Guo,“A hybrid solver based on the integral equation method and vector finite-element method for 3D controlled-source electromagnetic method modeling,”Geophysics,vol.83,no.5,pp.E319-E333,Sep.2018.矢量有限元离散网格中二次电场E _s满足表达式18)： According to Liu, R., RWGuo, JXLiu, CYMa, and ZWGuo, “A hybrid solver based on the integral equation method and vector finite-element method for 3D controlled-source electromagnetic method modeling,” Geophysics, vol. 83, no. 5 ,pp.E319-E333,Sep.2018. The secondary electric field E _s in the vector finite element discrete grid satisfies expression 18):

KE _s＝S+S ^H 18)； KE _s = S+S ^H 18);

其中：K是整体刚度矩阵，S是与发射源在地下产生的初始电场相关项，S ^H是电磁场在计算区域边界产生的坡印廷矢量项；当计算区域内部单元用I标示，计算区域边界单元用B标示，表达式18)可以分解为表达式19)： Among them: K is the overall stiffness matrix, S is the item related to the initial electric field generated by the emission source underground, and ^SH is the Poynting vector item generated by the electromagnetic field at the boundary of the calculation area; when the internal unit of the calculation area is marked with I, the calculation area boundary The unit is marked with B, and expression 18) can be decomposed into expression 19):

其中：

是关于内部单元的电磁场积分项，根据电磁场在电导率边界的切向连续性，

内部会相互抵消。使用表达式19)的上半部分，计算区域内部二次电场满足表达式20)： among them:

Is the integral term of the electromagnetic field of the internal unit, according to the tangential continuity of the electromagnetic field at the conductivity boundary,

Internally will cancel each other out. Using the upper half of expression 19), the internal secondary electric field in the calculation area satisfies expression 20):

K _IIE _Is+K _IBE _Bs＝S _I 20)； K _II E _Is +K _IB E _Bs =S _I 20);

其中：K _II对称且稀疏；K _IB是不对称的。图6一号线框表示了矢量有限元法的计算区域，该计算区域边界远离异常体(水平黑色块状体)，使得E _Bs＝0。 Among them: K _{II is} symmetric and sparse; K _IB is asymmetric. The line frame No. 1 in Fig. 6 shows the calculation area of the vector finite element method. The boundary of the calculation area is far away from the abnormal body (horizontal black block), so that E _Bs =0.

如图6二号框所示，截断边界矢量有限元法将矢量有限元法计算区域截断为异常体及其单元厚度包裹层。表达式20)中的E _Bs与E _Is通过格林函数建立联系，截断边界矢量有限元法计算区域边界二次电场E _Bs可以由内部单元散射电流表示为表达式21)： As shown in the second box of Figure 6, the truncated boundary vector finite element method cuts the vector finite element method calculation area into anomalous body and its element thickness wrap. _{The relationship between E Bs} and E _Is in expression 20) is established through Green's function. The truncated boundary vector finite element method calculates the secondary electric field E _{Bs at the} boundary of the region, which can be expressed by the internal element scattering current as expression 21):

E _Bs＝g ^eeJ _I 21)； E _Bs = g ^ee J _I 21);

其中：g ^ee是联系内部单元中心处散射电流J _I与计算区域边界电场E _Bs。通过将内部单元棱边电场E _Is使用线性插值算子到单元中心，中心处散射电流J _I可表示为： Among them: g ^ee is the connection between the scattered current J _I at the center of the internal unit and the electric field E _{Bs at the} boundary of the calculation area. By applying a linear interpolation operator to the edge electric field E _{Is of the} internal element to the center of the element, the scattered current J _I at the center can be expressed as:

J _I＝V(σ-σ _*)N ^e(E _Is+E _IP) 22)； J _I =V(σ-σ _* ) ^Ne (E _Is +E _IP ) 22);

其中：E _IP是内部单元棱边初始电场，V是内部单元体积，N ^e是线性插值算子，σ _*是均匀地下介质电导率；σ是地下介质电导率。然后，计算区域边界电场E _Bs可以由内部单元电场矩阵形式表达为表达式23)： Among them: E _IP is the initial electric field at the edge of the internal unit, V is the internal unit volume, ^Ne is the linear interpolation operator, σ _* is the conductivity of the uniform underground medium; σ is the conductivity of the underground medium. Then, the calculation area boundary electric field E _Bs can be expressed in the form of the internal unit electric field matrix as Expression 23):

E _Bs＝G ^ee(E _Is+E _IP) 23)； E _Bs =G ^ee (E _Is +E _IP ) 23);

其中：G ^ee为表达式24)： Among them: G ^ee is the expression 24):

G ^ee＝g ^eeV(σ-σ _*)n ^e 24)； G ^ee =g ^ee V(σ-σ _* )n ^e 24);

将表达式23)代入表达式20)中可得，计算区域内部单元二次电场控制方程组为表达式25)：Substituting expression 23) into expression 20), it can be obtained that the control equations for the secondary electric field of the internal unit of the calculation area are expression 25):

(K _II+K _IBG ^ee)E _Is＝S _I-K _IBG ^eeE _IP 25)； (K _II +K _IB G ^ee )E _Is =S _I -K _IB G ^ee E _IP 25);

通过求解方程组表达式25)，在得到E _Is后，使用表达式22)得到足迹内离散单元中心处电流强度J _I。 By solving the equation system expression 25), after obtaining E _Is , use expression 22) to obtain the current intensity J _I at the center of the discrete element in the footprint.

4、计算足迹区域内散射电流在第一测站测点的磁场，如下：4. Calculate the magnetic field of the scattered current at the measuring point of the first station in the footprint area, as follows:

足迹内异常体离散单元中心电流J _I在测点产生的磁场可以视为，均匀地下电偶极源在空气中产生的磁场。使用共轴装置时，地下单位电偶极源在测点产生的二次磁场表达为为表达式26)-28)： _{The magnetic field generated by the center current J I of the} discrete element of the abnormal body in the footprint at the measuring point can be regarded as the magnetic field generated in the air by a uniform underground electric dipole source. When using a coaxial device, the secondary magnetic field generated by the underground unit electric dipole source at the measuring point is expressed as expression 26)-28):

使用同面装置时，地下单位电偶极源在测点产生的二次磁场表达为表达式29)-31)：When using the same plane device, the secondary magnetic field generated by the underground unit electric dipole source at the measuring point is expressed as expression 29)-31):

其中：H下标表示电偶极子方向，上标表示产生的磁场方向。表达式26)-32)中各参数含义与表达式1)-6)一致。最终，测点二次磁场表达式可以写为矩阵形式表达式32)：Among them: the H subscript indicates the direction of the electric dipole, and the superscript indicates the direction of the generated magnetic field. The meaning of each parameter in expression 26)-32) is consistent with expression 1)-6). Finally, the expression of the secondary magnetic field at the measuring point can be written as a matrix expression 32):

H _s＝g ^emJ _I 32)； H _s =g ^em J _I 32);

其中g ^em是联系均匀半空间地下电流与航空装置测点磁场的格林函数。 Where g ^em is the Green's function linking the underground current in a uniform half-space with the magnetic field at the measuring point of the aeronautical device.

5、后续测站测点磁场计算，如下：5. The magnetic field calculation of the subsequent measuring stations is as follows:

表达式21)与表达式32)中的矩阵g ^ee与g ^em计算包含第一类贝塞尔积分，此类积分计算需要耗费大量的时间。幸运的是，本发明足迹引导的截断边界矢量有限元法在对第一测站进行矩阵g ^ee与g ^em计算后，后续测站可以重复使用。如图7所示，第一测站计算区域为实线框，第二测站计算区域为虚线框，第二测站计算区域内部离散单元与计算区域边界之间相对位置关系，均可以从第一测站计算区域内部离散单元与计算区域边界之间相对位置关系中相应找到，鉴于此第二测站可以使用第一测站的g ^ee。同样的，第二测站计算区域内离散单元与第二测站测点之间位置关系可以在第一测站计算区域内离散单元与第一测站测点之间位置关系中相应找到。鉴于此第二测站可以使用第一测站的g ^em。以此类推，后续测站可以重复使用第一测站计算过程中产生的矩阵g ^ee与g ^em，由此很大程度上加快后续测站二次磁场响应计算。 ^{The calculation of matrices g ee} and g ^{em in} Expression 21) and Expression 32) includes Bessel integrals of the first kind, and such integral calculation requires a lot of time. Fortunately, the truncated boundary vector finite element method guided by the footprint of the present invention can be reused in subsequent stations after the ^{matrix g ee} and g ^{em are calculated for the first station.} As shown in Figure 7, the calculation area of the first station is a solid line frame, and the calculation area of the second station is a dashed frame. The relative positional relationship between the discrete units within the calculation area of the second station and the boundary of the calculation area can be determined from the The relative positional relationship between the discrete unit within the calculation area of a station and the boundary of the calculation area can be found accordingly. In view of this, the second station can use the g ^{ee of the} first station. Similarly, the positional relationship between the discrete unit in the calculation area of the second station and the measurement point of the second station can be found in the position relationship between the discrete unit and the measurement point of the first station in the calculation area of the first station. In view of this, the second station can use the g ^{em of the} first station. By analogy, subsequent stations can reuse the matrices g ^ee and g ^em generated in the calculation process of the first station, which greatly speeds up the calculation of the secondary magnetic field response of the subsequent station.

以下结合具体实例对本发明的方法在不同发射装置及模型的应用进行详细说明。The application of the method of the present invention to different launching devices and models will be described in detail below with reference to specific examples.

本发明计算了一维层状模型，已验证本发明计算精度。航空电磁法使用图2所示装置，发射线圈与接收线圈相距7.5m，发射源接收装置距离地面40m。首先使用图5所示，首先使用均匀立方体网格剖分400×400×100m ³大小的计算区域，评估该区域内离散单元电流在测点产生的二次磁场与解析解误差，从而确定该区域是否满足足迹区域要求。计算结果表明在使用共轴装置时，400×400×100m ³体积中的激发电流在测点x方向二次磁场响应数值解实部误差0％，虚部误差2％；在使用同面装置时，400×400×100m ³体积中的激发电流在测点z方向二次磁场响应数值解实部误差3％，虚部误差2％。由此判断400×400×100m ³体积为该航空电磁法装置的足迹区域。图6虚线框为该航空电磁法装置足迹区域，其内部三号线框为足迹引导的截断边界矢量有限元法对层状介质模型的计算区域。表1给出了层状模型足迹引导的截断边界矢量有限元法(本发明方法)数值解与解析解在不同装置形式下的精度对比。由表1数据可知不同装置形式下，数值实部与虚部误差都保持在5％以内，符合航空电磁法精度要求。其中Hsx是x方向磁场计算结果，Hsz是z方向磁场计算结果。 The invention calculates a one-dimensional layered model, which has verified the calculation accuracy of the invention. The aviation electromagnetic method uses the device shown in Figure 2, the transmitting coil and the receiving coil are separated by 7.5m, and the transmitting source receiving device is 40m away from the ground. First use Figure 5, first use a uniform cubic grid to divide a 400×400×100m ³ calculation area, evaluate the secondary magnetic field generated by the discrete unit current in the area at the measuring point and the analytical solution error, so as to determine the area Whether to meet the footprint area requirements. The calculation results show that when using a coaxial device, ^{the excitation current in a volume of 400×400×100m 3} in the x direction of the measuring point has a secondary magnetic field response with a numerical solution error of 0% for the real part and 2% for the imaginary part; when using the same-plane device , The excitation current in the volume of 400×400×100m ³ in the z direction of the measuring point, the real part error of the numerical solution of the secondary magnetic field response is 3%, and the imaginary part error is 2%. Therefore, it is judged that the ^{volume of 400×400×100m 3} is the footprint area of the aerial electromagnetic method device. The dashed frame in Figure 6 is the footprint area of the aerial electromagnetic method device, and the internal line frame No. 3 is the calculation area of the layered medium model by the truncated boundary vector finite element method guided by the footprint. Table 1 shows the accuracy comparison between the numerical solution and the analytical solution of the truncated boundary vector finite element method (the method of the present invention) guided by the layered model footprint under different device forms. From the data in Table 1, it can be seen that under different device forms, the errors of the real and imaginary parts of the values are kept within 5%, which meets the accuracy requirements of the aviation electromagnetic method. Where Hsx is the calculation result of the magnetic field in the x direction, and Hsz is the calculation result of the magnetic field in the z direction.

表1 足迹引导的截断边界矢量有限元法数值解与解析解对比Table 1 Comparison of numerical solution and analytical solution of the truncated boundary vector finite element method guided by the footprint

同轴装置Coaxial device	实部(Hsx)Real part (Hsx)	虚部(Hsx)Imaginary part (Hsx)	实部(Hsz)Real part (Hsz)	虚部(Hsz)Imaginary part (Hsz)
截断边界矢量有限元法Truncated boundary vector finite element method	-3.899e-9-3.899e-9	2.238e-92.238e-9	-8.032e-10-8.032e-10	2.198e-102.198e-10
解析解Analytical solution	-3.940e-9-3.940e-9	2.346e-92.346e-9	-8.112e-10-8.112e-10	2.313e-102.313e-10
相对误差(％)Relative error(%)	-1.0-1.0	-4.6-4.6	-0.9-0.9	-4.9-4.9
共面装置Coplanar device	实部(Hsx)Real part (Hsx)	虚部(Hsx)Imaginary part (Hsx)	实部(Hsz)Real part (Hsz)	虚部(Hsz)Imaginary part (Hsz)
截断边界矢量有限元法Truncated boundary vector finite element method	8.117e-108.117e-10	-2.235e-10-2.235e-10	-7.977e-9-7.977e-9	4.500e-94.500e-9

解析解Analytical solution	8.119e-108.119e-10	-2.313e-10-2.313e-10	-7.928e-9-7.928e-9	4.694e-94.694e-9
相对误差(％)Relative error(%)	0.00.0	-3.3-3.3	0.60.6	-4.1-4.1

本发明计算了相对复杂三维模型，以验证本发明对复杂地质模型计算的高效性与准确性。复杂模型如图8所示，地下存在两个电导率不同的锯齿状低阻异常体。航空电磁法图2所示装置，第一测站位于x＝739.5m最后一个测站位于x＝211.5m，共45个测站，相邻测站相距12m。由于本次航空电磁法装置与层状模型一致，故装置的足迹区域同样为400×400×100m ³。本次航空电磁法第一测站使用足迹引导的截断边界矢量有限元法时，计算区域如图8一号框所示，区域大小为405×21×20m ³，选择3×3×4m ³网格剖分该区域。第二测站位于727.5m，使用足迹引导的截断边界矢量有限元法时，将第一测站剖分网格向左移动4个单元，得到第二测站计算区域。同理，后续测站陆续移动网格，在45号测站时，其计算区域如图8二号框所示。 The present invention calculates a relatively complex three-dimensional model to verify the high efficiency and accuracy of the present invention for the calculation of complex geological models. The complex model is shown in Figure 8. There are two jagged low-resistance anomalies with different electrical conductivity underground. The device shown in Figure 2 of the airborne electromagnetic method, the first station is at x=739.5m and the last station is at x=211.5m. There are 45 stations in total, and the distance between adjacent stations is 12m. Since the aerial electromagnetic method installation is consistent with the layered model, the footprint area of the installation is also 400×400×100m ³ . When using the truncated boundary vector finite element method guided by the footprint at the first station of the airborne electromagnetic method, the calculation area is shown in box No. 1 in Fig. 8. The area size is 405×21×20m ³ , and the 3×3×4m ³ network is selected. The grid is divided into the area. The second station is located at 727.5m. When using the truncated boundary vector finite element method guided by the footprint, the mesh of the first station is moved 4 units to the left to obtain the calculation area of the second station. In the same way, subsequent stations move the grid one after another. At station 45, the calculation area is shown in the second box in Figure 8.

在使用传统截断边界矢量有限元法计算该模型时，该模型的计算区域包括两个导电异常体和包裹层，计算区域如图8三号框所示，计算区域尺度为834×21×20m ³，表2对比了本发明方法与传统截断边界矢量有限元法在计算该模型单测站数据的计算效率。表2数据表明足迹引导的截断边界矢量有限元法形成的方程组自由度、需要的格林函数、计算时间明显小于传统截断边界矢量有限元法。两种方法在计算一次联系区域内部单元电场与边界节点电场之间的格林函数后，可以将其以矩阵形式存储。然而，由于足迹引导的截断边界矢量有限元法每次计算区域与接收装置相对位置不变，只需要计算一次联系计算区域内部单元电场与接收点磁场的格林函数，大大减少了计算时间。 When the traditional truncated boundary vector finite element method is used to calculate the model, the calculation area of the model includes two anomalous conductive bodies and the wrapping layer. The calculation area is shown in the third box in Figure 8. The calculation area size is 834×21×20m ³ Table 2 compares the calculation efficiency of the method of the present invention and the traditional truncated boundary vector finite element method in calculating the data of the single station of the model. The data in Table 2 show that the degree of freedom of the equations formed by the truncated boundary vector finite element method guided by the footprint, the required Green's function, and the calculation time are significantly less than the traditional truncated boundary vector finite element method. After calculating the Green's function between the electric field of the element inside the contact area and the electric field of the boundary node in the two methods, it can be stored in a matrix form. However, because the truncated boundary vector finite element method guided by the footprints calculates the relative position of the area and the receiving device every time, it only needs to calculate the Green's function linking the electric field of the unit inside the calculation area with the magnetic field of the receiving point, which greatly reduces the calculation time.

图9a、图9b、图10a和图10b显示了基于相对复杂三维模型的截断边界矢量有限元法和足迹引导截断边界矢量有限元法(本发明)之间的数值比较。计算得到的测点二次磁场响应结果以发射源全空间磁场的百万分之一(ppm)绘制，并对这两种方法进行比较。对于同轴装置和共面装置，足迹引导的截断边界矢量有限元法计算得到的水平方向磁场(Hsx)与垂直方向磁场(Hsz)的实部与虚部结果与传统截断边界矢量有限元法的结果一致。该结果进一步验证了本发明计算结果正确，具有普遍适用性。Figures 9a, 9b, 10a and 10b show the numerical comparison between the truncated boundary vector finite element method based on a relatively complex three-dimensional model and the footprint guided truncated boundary vector finite element method (the present invention). The calculated results of the secondary magnetic field response of the measuring point are plotted in parts per million (ppm) of the full-space magnetic field of the emission source, and the two methods are compared. For coaxial devices and coplanar devices, the real and imaginary parts of the horizontal magnetic field (Hsx) and vertical magnetic field (Hsz) calculated by the footprint-guided truncated boundary vector finite element method are the same as those of the traditional truncated boundary vector finite element method. The results are consistent. This result further verifies that the calculation result of the present invention is correct and has universal applicability.

表2 本发明足迹引导的截断边界矢量有限元法与传统截断边界矢量有限元法计算效率对比Table 2 Comparison of calculation efficiency between the truncated boundary vector finite element method guided by the present invention and the traditional truncated boundary vector finite element method

To	格林函数数量Number of Green's functions	方程组自由度Degrees of Freedom of the System of Equations	计算时间calculating time
截断边界矢量有限元法Truncated boundary vector finite element method	5582376055823760	82808280	147s147s
足迹引导截断边界矢量有限元法Footprint guided truncated boundary vector finite element method	1320690013206900	39903990	27s27s

本发明算法计算层状模型电磁响应具体实施如下：The specific implementation of the algorithm of the present invention to calculate the electromagnetic response of the layered model is as follows:

步骤一：以边长为4m的相同大小立方体作为单元体，以发射源地下投影为中心分别在x与y轴方向将剖分网格50个，往z轴方向剖分网格25个，形成100×100×25规模网格，计算区域大小为400×400×100m ³。 Step 1: Take a cube of the same size with a side length of 4m as the unit body, and divide the grid into 50 grids in the x and y axis directions and 25 grids in the z axis direction with the underground projection of the emission source as the center. 100×100×25 scale grid, the calculation area size is 400×400×100m ³ .

步骤二：计算400×400×100m ³内，离散剖分网格中心发射源激发电流大小；计算各离散剖分网格中心发射源激发电流在接收点产生的二次磁场，并将各网格内电流在测点产生的磁场叠加，形成测点磁场数值解。 Step 2: Calculate the ^{excitation current of the emission source in the center of the discrete grid within 400×400×100m 3} ; calculate the secondary magnetic field generated by the excitation current of the emission source in the center of each discrete grid at the receiving point, and combine the grids The magnetic field generated by the internal current at the measuring point is superimposed to form a numerical solution of the magnetic field at the measuring point.

步骤三：计算均匀半空间地下介质在测点产生的磁场解析解，以解析解为标准，计算数值解误差。得到误差在5％以内，确定航空电磁法装置足迹区域大小为400×400×100m ³。 Step 3: Calculate the analytical solution of the magnetic field generated by the uniform half-space underground medium at the measuring point, and use the analytical solution as the standard to calculate the numerical solution error. The error is within 5%, and the footprint area of the aviation electromagnetic device is determined to be 400×400×100m ³ .

步骤四：将无限延伸低阻板状模型，截断为400×400×4m ³的有限体积板状体。足迹引导的截断边界矢量有限元法将计算区域定为400×400×12m ³，使用4×4×4m ³网格剖分， Step 4: Cut the infinitely extending low-resistance plate-shaped model into a finite-volume plate-shaped body ^{of 400×400×4m 3.} The footprint-guided truncated boundary vector finite element method defines the calculation area as 400×400×12m ³ , and uses 4×4×4m ³ grid division.

步骤五：使用截断边界矢量有限元法计算400×400×12m ³区域内离散单元中心散射电流。 Step 5: Use the truncated boundary vector finite element method to calculate the scattered current at the center of the discrete element in the area of ^{400×400×12m 3.}

步骤六：计算400×400×12m ³区域内离散单元中心散射电流在测点产生的磁场，将每个单元产生的磁场叠加得到层状模型足迹引导的截断边界矢量有限元法数值解。 Step 6: Calculate the magnetic field generated by the scattered current from the center of the discrete unit in the 400×400×12m ³ area at the measuring point, and superimpose the magnetic field generated by each unit to obtain the numerical solution of the truncated boundary vector finite element method guided by the layered model footprint.

本领域技术人员将清楚本发明的范围不限制于以上讨论的示例，有可能对其进行若干改变和修改，而不脱离所附权利要求书限定的本发明的范围。尽管己经在附图和说明书中详细图示和描述了本发明，但这样的说明和描述仅是说明或示意性的，而非限制性的。本发明并不限于所公开的实施例。It will be clear to those skilled in the art that the scope of the present invention is not limited to the examples discussed above, and it is possible to make several changes and modifications to it without departing from the scope of the present invention defined by the appended claims. Although the present invention has been illustrated and described in detail in the drawings and specification, such description and description are only illustrative or illustrative, and not restrictive. The invention is not limited to the disclosed embodiments.

通过对附图，说明书和权利要求书的研究，在实施本发明时本领域技术人员可以理解和实现所公开的实施例的变形。在权利要求书中，术语“包括”不排除其他步骤或元素，而不定冠词“一个”或“一种”不排除多个。在彼此不同的从属权利要求中引用的某些措施的事实不意味着这些措施的组合不能被有利地使用。权利要求书中的任何参考标记不构成对本发明的范围的限制。By studying the drawings, the description and the claims, those skilled in the art can understand and realize the modifications of the disclosed embodiments when implementing the present invention. In the claims, the term "comprising" does not exclude other steps or elements, and the indefinite article "a" or "an" does not exclude a plurality. The fact that certain measures are cited in mutually different dependent claims does not mean that a combination of these measures cannot be used to advantage. Any reference signs in the claims do not constitute a limitation on the scope of the present invention.

Claims

一种足迹引导的高效航空电磁法数值模拟方法，其特征在于，包括以下步骤：A footprint-guided high-efficiency aviation electromagnetic method numerical simulation method is characterized in that it comprises the following steps:

步骤S100：使用均匀网格剖分航空电磁法探测区域；Step S100: Use a uniform grid to divide the aerial electromagnetic method to detect the area;

步骤S200：计算所使用航空电磁法观测装置的足迹区域在x、y、z三个方向的延伸尺度；Step S200: Calculate the extension scale of the footprint area of the aviation electromagnetic observation device used in the three directions of x, y, and z;

所述足迹区域在x与y方向网格数量相等，z方向网格数量在x或y网格数量的1/4至1/2之间，网格具体数量由均匀半空间模型测点电磁场解析解与数值解比率决定；The footprint area has the same number of grids in the x and y directions, and the number of grids in the z direction is between 1/4 and 1/2 of the number of x or y grids. The specific number of grids is analyzed by the electromagnetic field of the uniform half-space model. The ratio of solution to numerical solution is determined;

步骤S300：使用截断边界矢量有限元法将第一个测站对应的足迹作为异常体，计算足迹内存在的散射电流，并将计算过程中产生的第一格林函数进行存储备用；Step S300: Use the truncated boundary vector finite element method to take the footprint corresponding to the first station as an abnormal body, calculate the scattered current in the footprint, and store the first Green's function generated during the calculation process for future use;

所述第一格林函数在使用截断边界矢量有限元时产生，是用来计算足迹内离散单元散射电流在截断边界矢量有限元法计算边界产生的电磁场；The first Green's function is generated when a truncated boundary vector finite element is used, and is used to calculate the electromagnetic field generated by the boundary of the truncated boundary vector finite element method in the calculation of the scattered current of the discrete elements in the footprint;

步骤S400：使用第一测站足迹内离散单元散射电流计算测点电磁场响应，并对计算过程中产生的第二格林函数进行存储备用；Step S400: Calculate the electromagnetic field response of the measuring point by using the scattered current of the discrete elements in the footprint of the first measuring station, and store the second Green's function generated in the calculation process for use;

所述第二格林函数用来计算足迹内离散单元散射电流在测站观测点产生的电磁场；The second Green's function is used to calculate the electromagnetic field generated by the scattered current of the discrete unit in the footprint at the observation point of the station;

步骤S500：对后续测站重复使用第一测站在S300与S400步骤所存储的格林函数，先按照S300计算测站足迹内存在的散射电流，再按照S400计算站足迹内离散单元散射电流在测点电磁场。Step S500: Repeat the Green's function stored in steps S300 and S400 at the first station for subsequent stations, first calculate the scattered current in the station footprint according to S300, and then calculate the scattered current in the discrete unit in the station footprint according to S400. Point the electromagnetic field.
根据权利要求1所述的一种足迹引导的高效航空电磁法数值模拟方法，其特征在于，所述步骤S100中网格为规则六面体网格。The footprint-guided high-efficiency aviation electromagnetic method numerical simulation method according to claim 1, wherein the grid in the step S100 is a regular hexahedral grid.
根据权利要求1所述的一种足迹引导的高效航空电磁法数值模拟方法，其特征在于，步骤S200中：要求足迹区域以发射源在地表投影为中心，按4:4:1比例逐步增加x、y、z方向的网格数量，直到足迹内散射电流在测点产生的电磁场与测点解析解误差不超过5％。The footprint-guided high-efficiency aviation electromagnetic method numerical simulation method according to claim 1, characterized in that, in step S200: the footprint area is required to be centered on the projection of the emission source on the ground surface, and x is gradually increased in a ratio of 4:4:1 The number of grids in the, y, and z directions until the electromagnetic field generated by the scattered current in the footprint at the measuring point and the analytical solution error of the measuring point are not more than 5%.
根据权利要求3所述的一种足迹引导的高效航空电磁法数值模拟方法，其特征在于，所述步骤S200具体是：The footprint-guided high-efficiency aviation electromagnetic method numerical simulation method according to claim 3, wherein the step S200 is specifically:

步骤S210：计算航空电磁装置发射源在接收点产生的二次电磁场解析解；计算航空电磁装置发射源正下方均匀地下A区域的激发电流，并计算A区域内激发电流在航空电磁测量装置产生的二次电磁场数值解；此处A区域为400×400×100m ³的区域； Step S210: Calculate the analytical solution of the secondary electromagnetic field generated by the emission source of the aviation electromagnetic device at the receiving point; calculate the excitation current in the uniform underground area A directly below the emission source of the aviation electromagnetic device, and calculate the excitation current in the area A generated by the aviation electromagnetic measuring device Numerical solution of the secondary electromagnetic field; here, the area A is an area of 400×400×100m ³ ;

步骤S220：计算得到的二次电磁场数值解与得到的二次电磁场解析解的相对误差；Step S220: the relative error between the calculated numerical solution of the secondary electromagnetic field and the obtained analytical solution of the secondary electromagnetic field;

步骤S230：相对误差大小判断，若磁场实部与虚部相对误差均不超过5％，则步骤S210中的A区域体积定义为该航空电磁装置足迹；Step S230: Determine the magnitude of the relative error. If the relative error between the real and imaginary parts of the magnetic field does not exceed 5%, the volume of area A in Step S210 is defined as the footprint of the aviation electromagnetic device;

若磁场实部与虚部相对误差均超过5％，则按照4:4:1比例扩展重新定义S210中A区域并重新计算A区域内的激发电流和A区域内激发电流在航空电磁测量装置产生的二次电磁场数值解，返回步骤S220。If the relative error between the real and imaginary parts of the magnetic field exceeds 5%, redefine the area A in S210 according to the 4:4:1 ratio expansion and recalculate the excitation current in the A area and the excitation current in the A area generated by the aviation electromagnetic measuring device For the numerical solution of the secondary electromagnetic field of, return to step S220.
根据权利要求4所述的一种足迹引导的高效航空电磁法数值模拟方法，其特征在于，步骤S300中不论测站下方异常电导率规模，只需要将足迹区域作为计算目标区域，使用截断边界矢量有限元法计算足迹内存在的散射电流。A footprint-guided high-efficiency airborne electromagnetic method numerical simulation method according to claim 4, wherein in step S300, regardless of the scale of the abnormal conductivity below the station, only the footprint area is used as the calculation target area, and the truncated boundary vector is used. The finite element method calculates the scattered current present in the footprint.
根据权利要求4所述的一种足迹引导的高效航空电磁法数值模拟方法，其特征在于，所述步骤S300中截断边界矢量有限元法的使用包括以下步骤：The footprint-guided high-efficiency aviation electromagnetic method numerical simulation method according to claim 4, wherein the use of the truncated boundary vector finite element method in the step S300 includes the following steps:

步骤S310：将计算区域定义为足迹及其一个单元厚度包裹层，电场定义在计算区域剖分单元棱边中心点，使用矢量有限元理论建立电场所满足的待求解方程组；Step S310: The calculation area is defined as the footprint and an element thickness wrapper, the electric field is defined at the center point of the edge of the division element of the calculation area, and the vector finite element theory is used to establish a set of equations to be solved that the electric field satisfies;

步骤S320：计算第一格林函数，联系足迹内单元中心电流与计算区域边界电场，并将以第一格林函数排列为矩阵形式，行数为计算区域边界棱边数，列数为足迹区域内部单元棱边数；Step S320: Calculate the first Green's function, connect the center current of the unit in the footprint and the boundary electric field of the calculation area, and arrange the first Green's function as a matrix, the number of rows is the number of edges of the calculation area, and the number of columns is the internal unit of the footprint area Edge number

步骤S330：将步骤S320得到的第一格林函数矩阵存储备用；Step S330: store the first Green's function matrix obtained in step S320 for use;

步骤S340：形成计算区域边界棱边中心电场与足迹内部单元棱边中心电场关系式方程组；Step S340: Form a system of equations for the relationship between the center electric field at the edge of the calculation area and the center electric field at the edge of the footprint;

步骤S350：将步骤S310中计算区域边界电场使用步骤S340得到的关系方程组表达式替换，形成截断边界矢量有限元法控制方程组；Step S350: Replace the calculated regional boundary electric field in step S310 with the relational equation set expression obtained in step S340 to form a truncated boundary vector finite element method control equation set;

步骤S360：求解方程组，得到足迹内单元棱边中心电场值，使用线性插值得到单元中心散射电流。Step S360: Solve the equations to obtain the electric field value at the center of the edge of the element in the footprint, and use linear interpolation to obtain the scattered current at the center of the element.
根据权利要求6所述的一种足迹引导的高效航空电磁法数值模拟方法，其特征在于，The footprint-guided high-efficiency aviation electromagnetic method numerical simulation method according to claim 6, characterized in that:

步骤S410：计算联系第一测站测点与第一测站足迹内离散单元散射电流第二格林函数，并排列为行矩阵，矩阵列数为散射电流个数，即足迹离散单元个数；Step S410: Calculate and connect the second Green's function of the scattered currents of the discrete units in the footprint of the first station and the first station, and arrange them into a row matrix. The number of columns in the matrix is the number of scattered currents, that is, the number of discrete units of the footprint;

步骤S420：将步骤S410形成的第二格林函数行矩阵存储备用；Step S420: store the second Green's function row matrix formed in step S410 for use;

步骤S430：将第二格林函数行矩阵乘以散射电流形成的列矩阵，得到测点磁场值。Step S430: Multiply the second Green's function row matrix by the column matrix formed by the scattered current to obtain the magnetic field value of the measuring point.
根据权利要求7所述的一种足迹引导的高效航空电磁法数值模拟方法，其特征在于，所述步骤S500包括以下步骤：The footprint-guided high-efficiency aviation electromagnetic method numerical simulation method according to claim 7, wherein the step S500 comprises the following steps:

步骤S510：重复步骤S310，调用步骤S330存储的第一格林函数矩阵，重复步骤S340-S360，得到后续测站足迹内单元散射电流；Step S510: Repeat step S310, call the first Green's function matrix stored in step S330, and repeat steps S340-S360 to obtain the cell scattered current in the footprint of the subsequent station;

步骤S520：调用步骤S420存储的第二格林函数行矩阵，重复步骤S430，得到后续测站测点磁场值。Step S520: Call the second Green's function row matrix stored in step S420, repeat step S430 to obtain the magnetic field value of the subsequent measuring station.