0 引言

航空电磁法(airborne electromagnetic method)因其不受地形影响、测量面积大等优点而广泛应用于复杂地区的矿产勘查、地质填图、地下水勘查及环境检测等领域。1984年Stanmac和McPhar公司在加拿大成功实现了固定翼飞机AEM系统的首次飞行, 标志着AEM的诞生, 之后航空电磁法测量平台和系统配置不断创新和发展^[1]。我国于20世纪60年代初开始研制长导线式航空电磁系统, 之后又陆续研制了双频和三频硬架式航电仪^[2]。随着航空电磁技术的发展, 航空电磁正演计算方法一直是理论研究的重要内容。Newman等^[3]应用交错有限差分(FD)分别进行了二维和三维带地形的三维航空电磁响应。朱凯光等^[4]开展了一维层状模型航空电磁法正演模拟, 同时研究了航空电磁响应的影响因素和一定条件下有效探测深度。高亮等^[5]运用Impulse系统(频率域吊舱式直升机航空电磁测量系统)对频率域航空电磁系统进行了磁性条件下频率域航空电磁正演研究。李小康^[6]基于MPI(message passing interface, 用于并行计算的通信协议)运用有限元法并行计算了2.5维频率域航空响应。曲昕馨^[7]基于有限差分法实现了三维频率域航空电磁法正演模拟, 同时分析了线圈姿态对结果的影响并给出了校正方案。王卫平等^[8]在实现二维和三维航空电磁正演模拟的基础上研究了地形影响, 并使用地形校正法对典型地形地电模型进行校正。殷长春等^[9]运用非结构化网格的有限元法进行了2.5维航空电磁正演模拟, 分析总结了典型地形对航空响应的影响及其特征。黄威等^[10]利用三维矢量有限元正演模拟研究了覆盖层和垂直接触带等典型地电构造对航空电磁响应的影响特征。Li等^[11]运用等参有限元法进行了三维源/二维地电模型(2.5维)的频率域航空电磁法正反演算法研究。

由于航空电磁方法探测环境复杂, 地下介质的各向异性特征越来越受到重视。Yin等^[12-13]开发了用于计算层状大地任意各向异性中直流电场和磁场的算法; 处理了任意各向异性介质中正演问题。Avdeev等^[14]研究了各向异性对航空电磁测量的影响, 但其研究结果只适用于三维简单各向异性地质模型。Yin等^[15-16]研究各向异性介质对于一维介质模型的影响; 考虑了各向异性对机载电磁系统响应的影响, 并分析了实测航空电磁数据的各向异性主轴方向。Liu等^[17]运用有限差分法进行了任意三维各向异性模型的航空电磁响应计算, 取得了较好的应用效果。Yin等^[18]进行了起伏地形任意各向异性模型下的三维时域航空电磁响应正演计算。Huang等^[19]采用谱元法对航空电磁各向异性响应进行三维正演模拟, 得出了不同各向异性异常体的电磁响应的收敛条件及其识别方式。上述研究中, 各向异性介质电磁模拟研究主要以有限差分为主, 对于地下复杂结构模拟难以达到满意的模拟精度。本文采用矢量有限元法开展三维任意各向异性介质频率域航空电磁响应模拟算法, 利用矢量有限元法求解散射电场耦合方程, 通过将总场分解为一次场(背景场)和二次场(散射场)来进行各向异性介质中的频率域三维航空电磁正演模拟。针对三维模型计算量大的问题, 本文采用共享内存直接求解器PARDISO对大规模稀疏矩阵并行计算, 大大提高了计算效率。

1 正演理论

假设时谐因子为e^iωt, 其中, , ω为角频率, t为时间。为了消除源的奇异性, 对总场进行分解:

式中:E为电场强度; H为磁场强度; 下标p和s分别代表一次场和二次场。为了在保证精度的同时能提高计算效率, 本文以空气为介质的均匀全空间为背景对一次场采用解析式计算, 采用的源为垂直磁偶极源(对应水平共面装置, 简称HCP), 引入谢坤诺夫矢量势F, 通过解偏微分方程可得矢量势F表达式:

式中:m为磁距; r为计算坐标点到源点的距离, r=; k、μ₀分别为真空中的波数和磁导率; u_z为z轴方向上的单位向量。

由此得到一次场电场和磁场的表达式:

式中:x、y、z为所计算点的坐标值; u_x、u_y分别为x、y方向的单位向量; 为阻抗率, 为导纳率, 其中, μ为磁导率, ε为介电常数, σ为电导率。

对于二次场的计算, 下式为正演模拟中二次场形式下的麦克斯韦方程组:

(6)式和(7)式中:J_s=[(σ-σ_p)+iω(ε-ε_p)]E_p, 为等效场源项, σ_p为背景电导率。在各向异性介质中, σ及σ_p用一个3×3的张量矩阵表示^[20]。为实现任意各向异性电导率介质模型, 利用欧拉旋转从主轴坐标系转换到模拟坐标系, 主轴电导率张量用如下形式表示:

其下标x、y和z为笛卡尔坐标系的方向。欧拉旋转计算过程如下:

式中:D为主旋转矩阵; D_x、D_y、D_z分别是绕着x轴、y轴和z轴旋的对应矩阵(图 1), φ、ψ、χ分别为绕着x轴、y轴、z轴旋转的角度, 即

最后得到如下形式的任意各向异性介质中的电导率张量, 该张量具有对称和正定性质^[21]:

假定磁导率为真空中的磁导率μ₀, 不考虑存在位移电流的情况, 对式(6)两边同时取旋度, 将(7)式代入得到二次电场的双旋度方程:

采用矢量有限元法, 用六面网格将模型剖分, 并将模拟区域离散成N_e个单元, 对式(12)采用伽辽金加权余量法可得:

式中:R_k为单元的残数加权积分之和; N_k^e为第e个单元k棱边的插值函数; Ω_e为单元体积。

由于模拟区域, 同时利用矢量恒等式和高斯定理, 式(13)可以表示成:

式中, n_s为面外法线单位向量。

使用狄利克雷边界条件, n为边界上的法向量, 即假设在远离异常体的边界处二次异常场衰减为0, 棱边元赋值为0, 该条件使得电磁场呈指数衰减, 从而边界处不产生反射并远离异常区^[22]:

根据物性变化界面处电磁场切向连续性及所选用边界条件, 式(14)的第二项可以忽略, 改写成:

根据矢量有限元法, 即将自由度赋于棱边而不是单元节点, 可得到任意单元内场的展开式:

式中, N_i^e和E_i^e分别为第e个六面体单元第i条棱边的矢量基函数和切向电场值。通过xyz坐标系和ξηζ坐标系转换(图 2), 各棱边矢量基函数如下:

式中, l_i^e为该六面体单元第i条棱的边长。

E_s=[E^e], 将式(17)代入式(16)得:

K^e和M^e即为单元刚度矩阵, 具体如下:

式中:k=1, …, 12;l=1, …, 12。之后, 引入雅克行列式将函数转换为以ξ、η和ζ表示的函数, 则式(22)和(23)转化为

将单元刚度矩阵M^e和K^e按式(22)进行全局刚度矩阵合成, 可以得到大型线性方程组:

式中:A为稀疏复对称系数矩阵; x为带求解的场值; b为源端项。引入PARDISO求解器, 运用其内嵌套分割并行算法对方程共享内存并行直接求解。求解出电场值后, 磁场值可根据法拉第定律求得:

2 正演模拟结果 2.1 精度验证

本文借鉴Liu等^[17]的三维各向异性模型进行精度验证。三层介质模型如图 3所示, 其中各向异性低阻层厚度为100 m, 取主轴电导率元素为σ_xx=1 S/m、σ_yy=10 S/m、σ_zz=1 S/m, 并令主电导率张量绕y轴旋转45°, 上方为空气层, 下方为σ=0.01 S/m的各向同性介质高阻均匀半空间。采用水平磁偶极子源, 接收线圈和发射线圈均离地表 20 m, 发射频率为900 Hz。三维有限元计算结果与一维拟解析解结果对比如图 4所示。由图 4可见, 一维解析解和本文三维矢量有限元模拟结果能很好地吻合, 说明了本文算法的可靠性。

2.2 各向异性介质对航空电磁响应影响特征

为了进一步分析不同电阻率各向异性特征对航空电磁响应的影响, 建立如图 5所示的三维地电模型。该模型为一个高导异常体埋藏在均匀半空间中, 模型大小为1 540 m×1 540 m×1 540 m, 异常体大小为50 m×50 m×50 m, 异常体中心位于坐标系原点, 埋深为30 m, 三维模型剖分为42×42×46个单元, 包含扩边区域和目标区域。其中, 扩边区单元为呈2倍扩展的渐变粗网格, 目标区域为均匀细网格。同时, 针对空中场的(空间)衰减较地中有耗衰减得慢的问题, 对于空气层的扩边多, 地下区域扩边相对少, 尽可能减少截断边界对模拟区域的影响, 从而在保证了计算精度的同时提高计算速度。长源发射频率为900 Hz, 记录磁场H_z分量, 发射线圈和接收线圈距离为10 m, 离地面20 m, 点距为10 m, 每条测线共计21个测点。

首先, 当各向异性异常体的主电导率张量绕不同主轴旋转后, 得到倾斜各向异性电导率张量, 设各向异性异常体主电导率参量为

然后分别讨论以下不同各向异性情况对航空电磁响应影响规律:1)将(28)式的主电导率参量分别绕着x、z轴旋转0°、45°和90°, 分别研究磁场垂直分量的实部和虚部电磁响应特征; 2)各向异性围岩电阻率变化响应特征。

2.2.1 绕z轴旋转各向异性磁场响应特征

当异常体主电导率绕z主轴旋转0°、45°、90°时, 相当于目标体从各向异性介质变为倾斜各向异性介质的过程。图 6为上述电导率张量绕z轴旋转不同角度后得到的磁场H_z分量在xy面的响应特征。由于垂直磁偶极子所产生的电流在地层内主要沿xy面流动, 因此磁场响应H_z分量沿着电导率的旋转角度旋转, 场的分布形态和主轴的方位对应良好, 可用场值的形态判断出主轴方位。根据磁场特征, 可以确定异常体的各向异性主方向以及主电导率张量旋转角度。

图 7所示为绕z轴旋转HCP装置H_z分量观测记录结果。对比图中实虚分量各向同性(黑实线)与不同旋转角度各向异性(黑虚线)在信号响应振幅上有较大的差异, 可见各向异性参数对航空电磁观测结果影响很大。由各向同性变为各向异性条件下, 实分量响应信号由单峰值信号逐渐变为双峰值信号, 并且信号振幅有较大幅度的衰减(图 7a)。H_z虚分量相对影响较小, 信号相位没有明显变化(图 7b)。

2.2.2 绕x轴旋转各向异性磁场响应特征

当各向异性电性参数沿着x轴分别旋转0°、45°和90°。HCP装置磁场H_z分量响应如图 8所示。随着旋转角度的变化, HCP装置H_z分量的实部随着旋转角度达到90°时, 磁场分量响应从双峰异常变成了关于x轴对称的单峰异常, 而虚部变化异常较小。同理, 当沿y轴旋转时, HCP装置H_z分量的实部随着旋转角度的变化从单峰异常变成关于y轴对称的单峰异常。

图 9所示为绕x轴旋转HCP装置H_z分量观测记录结果。与上述绕z轴旋转结果类似, 实虚分量各向同性(黑实线)与不同旋转角度各向异性(黑虚线)在信号响应振幅上有较大的差异, 实分量在不同旋转角度情况下, 响应信号振幅影响更大, 在目标位置响应信号振幅随着旋转角度的变化, 振幅峰值由-400×10^-9 A/m(旋转90°)转变为-100×10^-9 A/m(旋转0°)。相比而言, 虚分量幅值变化率相对较小。

2.2.3 各向异性围岩与各向异性异常体的电磁磁场分布特征

当地下介质呈交替分布时, 围岩电性参数呈各向异性特征。在上述模型(图 5)的基础上改变围岩电性参数各向异性特征, 使围岩沿y方向电导率增大, 围岩主电导率张量为

同时, 各向异性异常体的主电导率张量分别绕z、x轴旋转45°, 计算得到异常体在各向异性围岩内的电磁响应(图 10)。对比图 6和图 8围岩为各向同性旋转相同角度的磁场H_z分量响应图可见, 围岩各向异性对磁场响应影响较大, 特别是在目标信号能量分布上产生较大变化, 由对称的双峰或单峰结构变为不规则的能量团分布。当围岩沿y轴电导率增大则在该方向呈良导性质, 而相应的x轴方向呈高阻性质, 由于电磁波在良导和高阻的情况下有不同的衰减速度, 因此可以看出磁场H_z分量平面分布沿着y方向延展。

图 11所示为HCP装置记录所得H_z响应曲线。图中黑实线为常规各向同性记录结果, 黑虚线为异常体各向异性旋转45°结果, 黑点划线为背景介质同时为各向异性记录结果。由图 11实虚分量绕x、z轴的旋转HCP装置H_z分量记录结果可见, 当背景介质存在各向异性特征时, 记录信号在形状上与背景各向同性-目标各向异性类似, 但在信号振幅值上有较大的差异, 与完全各向同性相比, 目标响应信号较弱, 说明各向异性造成目标体的识别能力降低。

3 结论与讨论

本文基于矢量有限元法并通过分离总场提出了任意各向异性情况下三维频率域航空电磁正演模拟算法, 通过典型的理论模型与一维解析解对比, 验证了该算法的精度和有效性。针对典型的各向异性电阻率模型, 分别计算得到了不同参数条件下的磁场响应特征并总结得出以下规律:

1) 磁场分布和主电导率张量的旋转方向之间有很强的相关性, 为识别各向异性异常体的主轴方向和主轴电导率提供依据, 根据这些性质可以进一步识别异常体的层理性质, 即其走向和倾角。

2) 当围岩呈各向异性时, 由于各向异性介质中的通道效应, 即电流会沿呈导电性质的方向聚焦, 磁场响应会随之改变, 记录信号相比于各向同性介质在响应振幅和相位上都有较大的变化, 且异常体响应信号变弱。

3) 航空电磁响应受到各向异性介质影响较大, 但其也有明显的特征, 根据其分布可以判断各向异性异常体的电性情况与层理方向, 同时也能判断围岩的电性情况。本文的计算结果和各向异性影响特征为实际解释航空电磁响应结果提供了依据。同时各向异性正演算法也为进一步开展航空电磁各向异性反演提供了基础。