CN110598367A - 一种足迹引导的高效航空电磁法数值模拟方法 - Google Patents

Abstract

本发明提供一种足迹引导的高效航空电磁法数值模拟方法，该方法包括：使用均匀网格剖分航空电磁法探测区域；逐步增加均匀单元数量，确定航空电磁法观测装置的足迹大小；对第一个测站使用截断边界矢量有限元法计算其足迹内存在的散射电流，并将此步骤中计算得到的格林函数进行存储；计算上一步得到的足迹内散射电流在第一个测点的产生的电磁场响应，并将并将此步骤中计算得到的格林函数进行存储；对后续测站重复前两个步骤。本发明选择将探测区域按单个测站足迹逐次计算，极大地缩小了计算区域；使用均匀网格剖分，在后续测站使用截断边界矢量有限元法时能重复使用格林函数，很大程度上提高了后续测站的计算效率。

Description

一种足迹引导的高效航空电磁法数值模拟方法

技术领域

本发明涉及航空电磁法正演技术领域，具体的涉及一种足迹引导的高效航空电磁法数值模拟方法。

背景技术

航空电磁法观测装置使用机载平台，具有高效率的特点，探测区域往往覆盖几十平方公里到上千平方公里。传统数值模拟方法要求将全部探测区域作为计算区域，这种大区域网格计算将耗费大量的计算内存与时间。

截断边界矢量有限元法只要求将探测区域内异常体作为计算区域，具有效率高的特点。但是在面对大范围异常体存在时，需要进行大量的格林函数计算，使得该方法效率备受质疑。

因此，开发一种计算区域适当且计算效率高的模拟方法具有重要意义。

发明内容

为了解决航空电磁法探测区域范围大带来的计算困难，本发明提供了一种足迹引导的高效航空电磁法数值模拟方法，该方法结合航空电磁法观测装置具有有限足迹区域(即航空电磁法足迹特征)与截断边界矢量有限元法，将整个航空电磁法探测区域按测站细分为多个子区域逐一计算测站测点电磁场，具有计算区域适当且计算效率高的特征，具体技术方案如下：

本发明提供一种足迹引导的高效航空电磁法数值模拟方法，包括以下步骤：

步骤S100：使用均匀网格剖分航空电磁法探测区域；

步骤S200：计算所使用航空电磁法观测装置的足迹区域在x、y、z三个方向的延伸尺度；所述足迹区域在x与y方向网格数量相等，且z方向网格数量为x或y方向网格数量的1/4至1/2之间；网格具体数量由均匀半空间模型测点电磁场解析解与数值解比率决定；

步骤S300：使用截断边界矢量有限元法将第一个测站对应的足迹作为异常体，计算足迹内存在的散射电流，并将计算过程中产生的格林函数进行存储备用；

所述格林函数在使用截断边界矢量有限元时产生，是用来计算足迹内离散单元散射电流在截断边界矢量有限元法计算边界产生的电磁场；

步骤S400：使用第一测站足迹内离散单元散射电流计算测点电磁场响应，并对计算过程中产生的格林函数进行存储备用；

所述格林函数用来计算足迹内离散单元散射电流在测站观测点产生的电磁场；

步骤S500：对后续测站重复使用第一测站在S300与S400步骤所存储的格林函数，先按照S300计算测站足迹内散射电流，再按照S400计算站足迹内散射电流在测点电磁场。

以上技术方案中优选的，所述步骤S100中网格为规则的六面体网格。

以上技术方案中优选的，步骤S200中具体要求足迹区域以发射源在地表投影为中心，按4:4:1比例逐步增加x、y、z方向的网格数量，直到足迹内散射电流在测点产生的电磁场与测点解析解误差小于5％。

以上技术方案中优选的，所述步骤S200具体是：

步骤S210：计算航空电磁发射源在接收点产生的二次电磁场解析解；

计算航空电磁装置发射源正下方均匀半空间地下A区域的激发电流，并计算区域内激发电流在航空电磁测量装置产生的二次场数值解；此处A区域为400×400×100m³的区域；

步骤S220：计算得到的二次磁场数值解与得到的二次磁场解析解的相对误差；

步骤S230：相对误差大小判断，若磁场实部与虚部相对误差均不超过5％，则步骤S210中的A区域体积定义为该航空电磁装置足迹；

若磁场实部与虚部相对误差均超过5％，则按照4:4:1比例扩展重新定义S210中A区域并计算区域内的激发电流和区域内激发电流在航空电磁测量装置产生的二次场数值解，返回步骤S220。

以上技术方案中优选的，步骤S300中不论测站下方异常电导率规模，只需要将足迹区域作为计算目标区域，使用截断边界矢量有限元法计算足迹内存在的散射电流。

以上技术方案中优选的，所述步骤S300中截断边界矢量有限元法的使用包括：

步骤S310：将计算区域定义为足迹及其一个单元厚度包裹层，电场定义在计算区域剖分单元棱边中心点，使用矢量有限元理论建立电场所满足的待求解方程组；

步骤S320：计算格林函数联系足迹内单元中心电流与计算区域边界电场，将格林函数排列为矩阵形式，行数为计算区域边界棱边数，列数为足迹区域内部单元棱边数；

步骤S330：将步骤S320关系方程组系数矩阵存储备用；

步骤S340：形成计算区域边界棱边中心电场由足迹内部单元棱边中心电场关系式方程组；

步骤S350：将步骤S310中计算区域边界电场使用步骤S320表达式替换，形成截断边界矢量有限元法控制方程组；

步骤S360：求解方程组，得到足迹内单元棱边中心电场值，使用线性插值得到单元中心散射电流。

以上技术方案中优选的，所述步骤S400具体包括以下步骤：

步骤S410：计算联系第一测站测点与第一测站足迹内离散单元散射电流格林函数，并排列为行矩阵，矩阵列数为散射电流个数，即足迹离散单元个数；

步骤S420：将步骤S410形成的行矩阵存储备用；

步骤S430：将格林函数行矩阵乘以散射电流形成的列矩阵，得到测点磁场值。

以上技术方案中优选的，所述步骤S500包括以下步骤：

步骤S510：重复步骤S310，调用步骤S330存储的格林函数，重复步骤S340-S360，得到后续测站足迹内单元散射电流；

步骤S520：调用步骤S420存储的格林函数，重复步骤S430，得到后续测站测点磁场值。

应用本发明的技术方案，效果是：

1、本发明提供的足迹引导的高效航空电磁法数值模拟方法，摒弃了传统数值模拟方法要求将所有探测区域作为计算区域，选择将探测区域按单个测站足迹逐次计算，极大的缩小了计算区域。

2、本发明提供的足迹引导的高效航空电磁法数值模拟方法，使用均匀网格剖分，所有测站足迹由相同网格组成。在后续测站使用截断边界矢量有限元法时，可以重复使用格林函数第一测站存储的格林函数，很大程度上提高了后续测站计算效率。

具体请参考根据本发明的一种足迹引导的高效航空电磁法数值模拟方法提出的各种实施例的如下描述，将使得本发明的上述和其他方面显而易见。

附图说明

图1是本实施例中足迹引导的高效航空电磁法数值模拟方法的示意图；

图2是测站装置形式及下方均匀地下介质均匀剖分示意图(使用均匀网格剖分航空电磁法探测区域)；

图3a是同轴装置形式地下均匀介质内电流实部分布图；

图3b是同轴装置形式地下均匀介质内电流虚部分布图；

图4a是共面装置形式地下均匀介质内电流实部分布图；

图4b是共面装置形式地下均匀介质内电流虚部分布图；

图5是地下足迹区域尺度评估过程图(实线框是第一次评估使用400×400×100m³计算区域，虚线框是后续按照4:4:1扩展后计算区域)；

图6是不同方法计算大型水平延伸薄板异常体所需计算区域示意图(一号框是传统微分法所需计算区域，二号框是传统截断边界矢量有限元法所需计算区域，三号框是足迹引导的截断边界矢量有限元法所需计算区域)；

图7是后续测站足迹区域示意图(第一测站计算区域如实线黑框所示，第二测站计算区域如虚线黑框所示)；

图8是锯齿复杂三维模型示意图(一号框是本实施例第一测站计算区域，二号框是本实施例第45测站计算区域，三号框是传统截断边界矢量有限元法在任意测站都需要的计算区域)；

图9a是使用同轴装置时本实施例方法与传统截断边界矢量有限元法在锯齿复杂三维模型的磁场数值解虚部对比图(黑点是本实施例方法计算结果虚部，黑线是传统截断边界矢量有限元法计算结果虚部)；

图9b是使用同轴装置时本实施例方法与传统截断边界矢量有限元法在锯齿复杂三维模型的磁场数值解实部对比图(圆圈是本实施例方法计算结果实部，虚线是传统截断边界矢量有限元法计算结果实部)；

图10a是使用共面装置时本实施例方法与传统截断边界矢量有限元法在锯齿复杂三维模型的磁场数值解虚部对比图(黑点是本实施例方法计算结果虚部，黑线是传统截断边界矢量有限元法计算结果虚部)；

图10b是使用共面装置时本实施例方法与传统截断边界矢量有限元法在锯齿复杂三维模型的磁场数值解实部对比图(圆圈是本实施例方法计算结果实部，虚线是传统截断边界矢量有限元法计算结果实部)。

具体实施方式

构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。

实施例：

见图1，本发明提供的一种足迹引导的高效航空电磁法数值模拟方法，包括以下步骤：

步骤S100：使用均匀网格剖分航空电磁法探测区域，优选的，此处使用的均匀网格为规则的六面体网格，网格x、y和z方向尺度相等，大小等于在5、10和25三个数之间选取，尺度越小精度要高。

步骤S200：计算所使用航空电磁法观测装置的足迹区域在x、y、z三个方向的延伸尺度，具体包括如下步骤：

第一步：计算单测站发射源在测点产生的均匀半空间解析解与全空间解析解，使用半空间解析解减去全空间解析解得到均匀地下介质内电流在测点产生的二次磁场解析解；

第二步：计算航空电磁装置发射源在其正下方均匀地下A区域(400×400×100m³区域)均匀网格单元中心位置产生的激发电流；

第三步：计算所定义体积区域(400×400×100m³区域)均匀网格单元内激发电流在测量装置处产生的二次场进行叠加，得到测点二次磁场数值解；

第四步：计算第三步得到的二次磁场数值解与第一步得到的二次磁场解析解的相对误差；

第五步：相对误差大小判断，若磁场实部与虚部相对误差均不超过5％，则步骤第二步中的A区域体积定义为该航空电磁装置足迹；

若磁场实部与虚部相对误差均超过5％，则按照4:4:1比例扩展(按照x、y与z方向)扩展均匀单元数量，重新定义第二步中A区域并计算区域内的激发电流和区域内激发电流在航空电磁测量装置产生的二次场数值解，返回步骤第三步。

步骤S300：使用截断边界矢量有限元法将第一个测站对应的足迹作为目标区域，计算足迹内存在的散射电流，并将计算过程中产生的格林函数进行存储备用；

本实施例中截断边界矢量有限元法的使用包括以下步骤：

步骤S310：将计算区域定义为足迹及其一个单元厚度包裹层，电场定义在计算区域剖分单元棱边中心点，使用矢量有限元理论建立电场所满足的待求解方程组；

步骤S320：计算格林函数联系足迹内单元中心电流与计算区域边界电场，将格林函数排列为矩阵形式，行数为计算区域边界棱边数，列数为足迹区域内部单元棱边数；

步骤S330：将步骤S320关系方程组系数矩阵存储备用；

步骤S340：形成计算区域边界棱边中心电场由足迹内部单元棱边中心电场关系式方程组；

步骤S350：将步骤S310中计算区域边界电场使用步骤S320表达式替换，形成截断边界矢量有限元法控制方程组；

步骤S360：求解方程组，得到足迹内单元棱边中心电场值，使用线性插值得到单元中心散射电流。

步骤S400：使用测站足迹内离散单元散射电流计算测点电磁场响应，步骤如下：

步骤S410：计算联系第一测站测点与第一测站足迹内离散单元散射电流格林函数，并排列为行矩阵，矩阵列数为散射电流个数，即足迹离散单元个数；

步骤S420：将步骤S410形成的行矩阵存储备用；

步骤S430：将格林函数行矩阵乘以散射电流形成的列矩阵，得到测点磁场值。

步骤S500：对后续测站重复使用第一测站在S300与S400步骤所存储的格林函数，先按照S300计算测站足迹内散射电流，再按照S400计算站足迹内散射电流在测点电磁场。

优选的，具体S500包括以下步骤：

步骤S510：重复步骤S310，调用步骤S330存储的格林函数，重复步骤S340-S360，得到后续测站足迹内单元散射电流；

步骤S520：调用步骤S420存储的格林函数，重复步骤S430，得到后续测站测点磁场值。

应用本实施例的方案，详情如下：

1、使用均匀网格剖分探测区域

如图2所示，航空电磁法使用机载平台搭载观测装置，观测装置由发射线圈与接收线圈组成。当发射线圈与接收线圈为水平线圈时称为共面装置形式，当发射线圈与接收线圈为垂直线圈时称为同轴装置形式。航空电磁法距离地表高度在40m之间，航线从右至左，发射线圈源位于接收线圈前方7.5m，发射线圈发射频率为10000hz。为了计算该航空电磁法装置在探测区域各个测站的磁场数据，本实施例使用均匀网格地下介质。

2、计算航空电磁观测装置足迹区域

航空电磁法发射装置所使用的垂直线圈时(同轴装置)可以当做水平磁偶极子源处理，发射装置所使用的水平线圈时(共面装置)可以当做垂直磁偶极子源处理。在使用时谐因子e^iwt时，设置直角坐标，地表为xy面，z坐标向下。

2.1、均匀地下初始电场求取，如下：

空气中发射源为单位水平磁偶极子源(同轴装置)在接收点(地下网格单元中心)产生的电场表达式满足表达式1)-3)：

E_z＝0 3)；

空气中单位垂直磁偶极子源(共面装置)在地下产生的电场表达式满足表达式4)-6)：

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

2.2、均匀地下发射源产生电流分布，如下：

图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³区域内，在足迹评估过程中，首先选定该区域计算测点二次磁场数值解。

2.3、测点二次磁场数值解与解析解求取，如下：

本实施例选择均匀半空间模型，通过比较两种不同方法计算地下介质在测点产生的磁场来定义航空电磁法装置足迹区域大小，具体是：首先，计算得到发射源在地下选定区域引起的电流分布之后；然后，计算选定区域网格内电流在航空电磁观测点引起的二次磁场叠加得到测点二次磁场数值解；最后，得到的二次磁场数值解与解析解相比较，判断所选电流区域是否已达到所要求规模。其中判断准则为以二次磁场解析解为标准，二次磁场数值解误差小于5％，解析解由测点磁场均匀半空间解析解减去全空间解析解得到。

在求取所选区域内网格二次磁场数值解时，网格单元内部电流被视为均匀分布，在计算其在测点产生的磁场时，网格单元被当做电偶极子源。在航空电磁法同轴装置时，单位电偶极子源在全空间产生的磁场为表达式7)-9)：

使用表达式1)-3)计算得到地下网格单元中心电场分布后，使用表达式7)-9)计算每个单元中心电流在测点产生的磁场，通过累加所选区域内单元在测点磁场贡献得到测点二次磁场数值解。

在航空电磁法共面装置时，单位电偶极子源在全空间产生的磁场为表达式10)-12)：

使用表达式4)-6)计算得到地下网格单元中心电场分布后，使用表达式10)-12)计算每个单元中心电流在测点产生的磁场，通过累加所选区域内单元在测点磁场贡献得到测点二次磁场数值解。

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

在求取所选区域内网格二次磁场解析解时，在同轴装置下即发射装置使用x方向磁偶极子源，接收装置测量x方向磁场，单位磁偶极发射源在全空间条件下在测点产生的磁场表达式为表达式13)：

表达式13)中各参数与表达式7)-12)中各参数一致。

在均匀半空间条件下单位磁偶极发射源在测点产生的磁场表达式为

表达式14)中各参数与表达式1)-6)中各参数一致。使用表达式14)中减去表达式15)中可以得到航空电磁法同轴装置发射源在测量点产生的二次磁场解析解。

在共面装置下即发射装置使用z方向磁偶极子源，接收装置测量z方向磁场，单位磁偶极发射源在全空间条件下在测点产生的磁场为表达式15)：

表达式15)中各参数与表达式7)-12)中各参数一致。

在均匀半空间条件下单位磁偶极发射源在测点产生的磁场为表达式16)：

表达式16)中各参数与表达式1)-6)中各参数一致。使用表达式16)中减去表达式15)中可以得到航空电磁法共面装置发射源在测量点产生的二次磁场解析解。

2.4、自适应扩展网格区域求取足迹范围，如下：

图5展示了逐步扩大网格区域，使得测点二次磁场数值解逐渐接近解析解，求得航空电磁装置足迹区域的过程。第一次数值解计算使用400×400×100m³区域内网格，当区域网格内电流在测点二次磁场数值解与解析解差异大于5％时，按4:4:1比例扩展网格区域。

3、使用截断边界矢量有限元法计算第一测站足迹内散射电流，如下：

本实施例基于二次电场矢量有限元法，有限元离散网格单元棱边中心电场满足表达式17)：

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

根据Liu,R.,R.W.Guo,J.X.Liu,C.Y.Ma,and Z.W.Guo,“A hybrid solver basedon the integral equation method and vector finite-element method for 3Dcontrolled-source electromagnetic method modeling,”Geophysics,vol.83,no.5,pp.E319-E333,Sep.2018.矢量有限元离散网格中二次电场E_s满足表达式18)：

KE_s＝S+S^H 18)；

其中：K是整体刚度矩阵，S是与发射源在地下产生的初始电场相关项，S^H是电磁场在计算区域边界产生的坡印廷矢量项；当计算区域内部单元用I标示，计算区域边界单元用B标示，表达式18)可以分解为表达式19)：

其中：是关于内部单元的电磁场积分项，根据电磁场在电导率边界的切向连续性，内部会相互抵消。使用表达式19)的上半部分，计算区域内部二次电场满足表达式20)：

K_IIE_Is+K_IBE_Bs＝S_I 20)；

其中：K_II对称且稀疏；K_IB是不对称的。图6一号线框表示了矢量有限元法的计算区域，该计算区域边界远离异常体(水平黑色块状体)，使得E_BS＝0。

如图6二号框所示，截断边界矢量有限元法将矢量有限元法计算区域截断为异常体及其单元厚度包裹层。表达式20)中的E_Bs与E_Is通过格林函数建立联系，截断边界矢量有限元法计算区域边界二次电场E_Bs可以由内部单元散射电流表示为表达式21)：

E_Bs＝g^eeJ_I 21)；

其中：g^ee是联系内部单元中心处散射电流J_I与计算区域边界电场E_Bs。通过将内部单元棱边电场E_Is使用线性插值算子到单元中心，中心处散射电流J_I可表示为：

J_I＝V(σ-σ_*)N^e(E_IS+E_IP) 22)；

其中：E_IP是内部单元棱边初始电场，V是内部单元体积，N^e是线性插值算子，σ_*是均匀地下介质电导率；σ是地下介质电导率。然后，计算区域边界电场E_Bs可以由内部单元电场矩阵形式表达为表达式23)：

E_Bs＝G^ee(E_Is+E_IP) 23)；

其中：G^ee为表达式24)：

G^ee＝g^eeV(σ-σ_*)N^e 24)；

将表达式23)代入表达式20)中可得，计算区域内部单元二次电场控制方程组为表达式25)：

(K_II+K_IBG^ee)E_Is＝S_I-K_IBG^eeE_IP 25)；

通过求解方程组表达式25)，在得到E_Is后，使用表达式22)得到足迹内离散单元中心处电流强度J_I。

4、计算足迹区域内散射电流在第一测站测点的磁场，如下：

足迹内异常体离散单元中心电流J_I在测点产生的磁场可以视为，均匀地下电偶极源在空气中产生的磁场。使用共轴装置时，地下单位电偶极源在测点产生的二次磁场表达为为表达式26)-28)：

使用同面装置时，地下单位电偶极源在测点产生的二次磁场表达为表达式29)-31)：

其中：H下标表示电偶极子方向，上标表示产生的磁场方向。表达式26)-32)中各参数含义与表达式1)-6)一致。最终，测点二次磁场表达式可以写为矩阵形式表达式32)：

H_s＝g^emJ_I 32)；

其中g^em是联系均匀半空间地下电流与航空装置测点磁场的格林函数。

5、后续测站测点磁场计算，如下：

表达式21)与表达式32)中的矩阵g^ee与g^em计算包含第一类贝塞尔积分，此类积分计算需要耗费大量的时间。幸运的是，本发明足迹引导的截断边界矢量有限元法在对第一测站进行矩阵g^ee与g^em计算后，后续测站可以重复使用。如图7所示，第一测站计算区域为实线框，第二测站计算区域为虚线框，第二测站计算区域内部离散单元与计算区域边界之间相对位置关系，均可以从第一测站计算区域内部离散单元与计算区域边界之间相对位置关系中相应找到，鉴于此第二测站可以使用第一测站的g^ee。同样的，第二测站计算区域内离散单元与第二测站测点之间位置关系可以在第一测站计算区域内离散单元与第一测站测点之间位置关系中相应找到。鉴于此第二测站可以使用第一测站的g^em。以此类推，后续测站可以重复使用第一测站计算过程中产生的矩阵g^ee与g^em，由此很大程度上加快后续测站二次磁场响应计算。

以下结合具体实例对本发明的方法在不同发射装置及模型的应用进行详细说明。

本发明计算了一维层状模型，已验证本发明计算精度。航空电磁法使用图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方向磁场计算结果。

表1足迹引导的截断边界矢量有限元法数值解与解析解对比

同轴装置	实部(Hsx)	虚部(Hsx)	实部(Hsz)	虚部(Hsz)
					截断边界矢量有限元法	-3.899e-9	2.238e-9	-8.032e-10	2.198e-10
解析解	-3.940e-9	2.346e-9	-8.112e-10	2.313e-10
					相对误差(％)	-1.0	-4.6	-0.9	-4.9
共面装置	实部(Hsx)	虚部(Hsx)	实部(Hsz)	虚部(Hsz)
					截断边界矢量有限元法	8.117e-10	-2.235e-10	-7.977e-9	4.500e-9
解析解	8.119e-10	-2.313e-10	-7.928e-9	4.694e-9
					相对误差(％)	0.0	-3.3	0.6	-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二号框所示。

在使用传统截断边界矢量有限元法计算该模型时，该模型的计算区域包括两个导电异常体和包裹层，计算区域如图8三号框所示，计算区域尺度为834×21×20m³，表2对比了本发明方法与传统截断边界矢量有限元法在计算该模型单测站数据的计算效率。表2数据表明足迹引导的截断边界矢量有限元法形成的方程组自由度、需要的格林函数、计算时间明显小于传统截断边界矢量有限元法。两种方法在计算一次联系区域内部单元电场与边界节点电场之间的格林函数后，可以将其以矩阵形式存储。然而，由于足迹引导的截断边界矢量有限元法每次计算区域与接收装置相对位置不变，只需要计算一次联系计算区域内部单元电场与接收点磁场的格林函数，大大减少了计算时间。

图9a、图9b、图10a和图10b显示了基于相对复杂三维模型的截断边界矢量有限元法和足迹引导截断边界矢量有限元法(本发明)之间的数值比较。计算得到的测点二次磁场响应结果以发射源全空间磁场的百万分之一(ppm)绘制，并对这两种方法进行比较。对于同轴装置和共面装置，足迹引导的截断边界矢量有限元法计算得到的水平方向磁场(Hsx)与垂直方向磁场(Hsz)的实部与虚部结果与传统截断边界矢量有限元法的结果一致。该结果进一步验证了本发明计算结果正确，具有普遍适用性。

表2本发明足迹引导的截断边界矢量有限元法与传统截断边界矢量有限元法计算效率对比

	格林函数数量	方程组自由度	计算时间
				截断边界矢量有限元法	55823760	8280	147s
足迹引导截断边界矢量有限元法	13206900	3990	27s

本发明算法计算层状模型电磁响应具体实施如下：

步骤一：以边长为4m的相同大小立方体作为单元体，以发射源地下投影为中心分别在x与y轴方向将剖分网格50个，往z轴方向剖分网格25个，形成100×100×25规模网格，计算区域大小为400×400×100m³。

步骤二：计算400×400×100m³内，离散剖分网格中心发射源激发电流大小；计算各离散剖分网格中心发射源激发电流在接收点产生的二次磁场，并将各网格内电流在测点产生的磁场叠加，形成测点磁场数值解。

步骤三：计算均匀半空间地下介质在测点产生的磁场解析解，以解析解为标准，计算数值解误差。得到误差在5％以内，确定航空电磁法装置足迹区域大小为400×400×100m³。

步骤四：将无限延伸低阻板状模型，截断为400×400×4m³的有限体积板状体。足迹引导的截断边界矢量有限元法将计算区域定为400×400×12m³，使用4×4×4m³网格剖分，

步骤五：使用截断边界矢量有限元法计算400×400×12m³区域内离散单元中心散射电流。

步骤六：计算400×400×12m³区域内离散单元中心散射电流在测点产生的磁场，将每个单元产生的磁场叠加得到层状模型足迹引导的截断边界矢量有限元法数值解。

本领域技术人员将清楚本发明的范围不限制于以上讨论的示例，有可能对其进行若干改变和修改，而不脱离所附权利要求书限定的本发明的范围。尽管己经在附图和说明书中详细图示和描述了本发明，但这样的说明和描述仅是说明或示意性的，而非限制性的。本发明并不限于所公开的实施例。

通过对附图，说明书和权利要求书的研究，在实施本发明时本领域技术人员可以理解和实现所公开的实施例的变形。在权利要求书中，术语“包括”不排除其他步骤或元素，而不定冠词“一个”或“一种”不排除多个。在彼此不同的从属权利要求中引用的某些措施的事实不意味着这些措施的组合不能被有利地使用。权利要求书中的任何参考标记不构成对本发明的范围的限制。

Claims (8)

1.一种足迹引导的高效航空电磁法数值模拟方法，其特征在于，包括以下步骤：

步骤S100：使用均匀网格剖分航空电磁法探测区域；

步骤S200：计算所使用航空电磁法观测装置的足迹区域在x、y、z三个方向的延伸尺度；

所述足迹区域在x与y方向网格数量相等，z方向网格数量在x或y网格数量的1/4至1/2之间，网格具体数量由均匀半空间模型测点电磁场解析解与数值解比率决定；

步骤S300：使用截断边界矢量有限元法将第一个测站对应的足迹作为异常体，计算足迹内存在的散射电流，并将计算过程中产生的格林函数进行存储备用；

所述格林函数在使用截断边界矢量有限元时产生，是用来计算足迹内离散单元散射电流在截断边界矢量有限元法计算边界产生的电磁场；

步骤S400：使用第一测站足迹内离散单元散射电流计算测点电磁场响应，并对计算过程中产生的格林函数进行存储备用；

所述格林函数用来计算足迹内离散单元散射电流在测站观测点产生的电磁场；

步骤S500：对后续测站重复使用第一测站在S300与S400步骤所存储的格林函数，先按照S300计算测站足迹内存在的散射电流，再按照S400计算站足迹内离散单元散射电流在测点电磁场。

2.根据权利要求1所述的一种足迹引导的高效航空电磁法数值模拟方法，其特征在于，所述步骤S100中网格为规则六面体网格。

3.根据权利要求1所述的一种足迹引导的高效航空电磁法数值模拟方法，其特征在于，步骤S200中：要求足迹区域以发射源在地表投影为中心，按4:4:1比例逐步增加x、y、z方向的网格数量，直到足迹内散射电流在测点产生的电磁场与测点解析解误差不超过5％。

4.根据权利要求3所述的一种足迹引导的高效航空电磁法数值模拟方法，其特征在于，所述步骤S200具体是：

步骤S210：计算航空电磁发射源在接收点产生的二次电磁场解析解；计算航空电磁装置发射源正下方均匀半空间地下A区域的激发电流，并计算区域内激发电流在航空电磁测量装置产生的二次场数值解；此处A区域为400×400×100m³的区域；

步骤S220：计算得到的二次磁场数值解与得到的二次磁场解析解的相对误差；

步骤S230：相对误差大小判断，若磁场实部与虚部相对误差均不超过5％，则步骤S210中的A区域体积定义为该航空电磁装置足迹；

若磁场实部与虚部相对误差均超过5％，则按照4:4:1比例扩展重新定义S210中A区域并计算区域内的激发电流和区域内激发电流在航空电磁测量装置产生的二次场数值解，返回步骤S220。

5.根据权利要求4所述的一种足迹引导的高效航空电磁法数值模拟方法，其特征在于，步骤S300中不论测站下方异常电导率规模，只需要将足迹区域作为计算目标区域，使用截断边界矢量有限元法计算足迹内存在的散射电流。

6.根据权利要求4所述的一种足迹引导的高效航空电磁法数值模拟方法，其特征在于，所述步骤S300中截断边界矢量有限元法的使用包括以下步骤：

步骤S310：将计算区域定义为足迹及其一个单元厚度包裹层，电场定义在计算区域剖分单元棱边中心点，使用矢量有限元理论建立电场所满足的待求解方程组；

步骤S320：计算格林函数联系足迹内单元中心电流与计算区域边界电场，将格林函数排列为矩阵形式，行数为计算区域边界棱边数，列数为足迹区域内部单元棱边数；

步骤S330：将步骤S320关系方程组系数矩阵存储备用；

步骤S340：形成计算区域边界棱边中心电场由足迹内部单元棱边中心电场关系式方程组；

步骤S350：将步骤S310中计算区域边界电场使用步骤S320表达式替换，形成截断边界矢量有限元法控制方程组；

步骤S360：求解方程组，得到足迹内单元棱边中心电场值，使用线性插值得到单元中心散射电流。

7.根据权利要求6所述的一种足迹引导的高效航空电磁法数值模拟方法，其特征在于，

步骤S410：计算联系第一测站测点与第一测站足迹内离散单元散射电流格林函数，并排列为行矩阵，矩阵列数为散射电流个数，即足迹离散单元个数；

步骤S420：将步骤S410形成的行矩阵存储备用；

步骤S430：将格林函数行矩阵乘以散射电流形成的列矩阵，得到测点磁场值。

8.根据权利要求7所述的一种足迹引导的高效航空电磁法数值模拟方法，其特征在于，所述步骤S500包括以下步骤：

步骤S510：重复步骤S310，调用步骤S330存储的格林函数，重复步骤S340-S360，得到后续测站足迹内单元散射电流；

步骤S520：调用步骤S420存储的格林函数，重复步骤S430，得到后续测站测点磁场值。