0 引言

固定翼航空电磁探测系统以固定翼飞机为载体，采用偶极-偶极收发装置，发射线圈安置于飞机周围，接收线圈由电缆线拖拽悬挂于飞机下方，可测量三分量电磁响应，具有勘探效率高、探测深度大、勘探面积广等优势，广泛应用于地质勘查中^{[1, 2, 3]}。固定翼航空电磁数据量庞大，反演解释耗时费力，一般采用电导率深度成像 ( conductivity-depth imaging，CDI) 技术快速处理海量航空电磁数据，得到地下电导率分布图，从而确定地下异常分布^{[4, 5]}。

目前国内外航空电磁探测的电导率深度成像方法很多，如Liu等^[6]提出了基于薄板模型的成像方法，其视电导率及视深度由最小二乘迭代法导出；Wolfgram等^[7]采用解卷积算法将任意波形的电磁响应转换为阶跃响应，通过τ域分解，得到地下电导率分布；Huang等^[8]提出了基于假层-半空间查表的电导率深度成像方法，引入电导率为0的假想层代表人文圈；朱凯光等^[9]利用神经网络获得更为准确的视电导率；陈小红等^[10]采用改进后的一阶、二阶电导率深度成像算法，提高了对地层电导率的灵敏性；毛立峰^[11]建立了σ-β数据表查询视电导率，消除了电磁响应的二值性。

在上述电导率深度成像算法中，直升机时间域电磁探测系统一般采用单分量测量、单分量成像；而传统固定翼航空电磁探测系统由于无法准确获取多分量电磁响应，因此，一般采用总场数据(即dB_total/dt=[(dB_x/dt)²+(dB_y/dt)²+(dB_z/dt)²])进行成像解释。虽然总场数据综合利用了多个分量的响应，但由于进行了多对一的转换，损失了多分量信息。近年来，国内外学者对电磁多分量响应特性进行了大量研究^{[12, 13, 14]}，发现水平分量响应(即x分量)有助于判断异常体的中心位置，垂直分量响应(即z分量)与水平分量相比具有较高的晚期道信噪比以及较好的垂向分辨率等。因此，若对航空电磁数据的多分量响应进行联合解释，不仅可以降低单一分量响应解释时描述的自由度，还能够实现对异常目标体几何形态特性的联合解释，提高成像精度。目前随着激光、惯导等技术的发展，航空电磁系统能够准确记录接收线圈姿态，从而通过姿态角度校正算法^[15]，获取多分量电磁响应，为多分量解释提供基础。

笔者根据磁场对导电异常体探测与单调性的优势^[16]提出了一种基于B场响应的双分量联合电导率深度成像算法。该算法基于固定翼时间域航空电磁正演理论，采用层状大地模型，推导出固定翼航空电磁系统x与z分量的B场响应，建立水平分量磁场(即B_x)-垂直分量磁场(即B_z)-电导率-飞行高度数据表，给出双分量B_x-B_z联合查表算法，最后通过理论仿真验证双分量查表法联合电导率深度成像算法的有效性。

建立固定翼航空电磁响应正演计算坐标系，如图 1所示，以发射线圈中心在地面的投影为坐标系原点O。假定固定翼航空电磁系统保持平稳飞行状态，发射线圈(Tx)中心距地面高度为h，发射磁偶极矩为m=I·A_Tx(I为发射电流强度，A_Tx为发射线圈的有效面积)，接收线圈(Rx)位于坐标系(x,y,z)处，发射线圈与接收线圈的水平距离r=,垂直距离为Δ,悬线长度L=。若大地模型为n层层状大地，其第j层的电导率、磁导率及厚度分别为σ_j，μ_j，d_j。

图 1 固定翼时间域航空电磁系统坐标系Fig.1 Coordinate of the fixed-wing time-domain electromagnetic system

图选项

根据电磁计算理论^[17]，接收线圈处的二次场频率域水平分量磁场响应B_x(w)和垂直分量磁场响应B_z(w)的表达式分别为

点击浏览公式

其中：J₀和J₁分别为零阶和一阶贝塞尔函数；w为角频率；λ为积分变量；μ₀为空气中的磁导率，μ₀=4×10^-7π H/m；反射系数r_TE为

其中:μ_j=μ₀。分别将式(1)、式(2)变换到s(复频率)域中，并除以s得到x分量和z分量的磁场瞬态响应表达式：

点击浏览公式

分别对式(3)、(4)做拉氏逆变换(L^-1)得到系统双分量B场响应，作为电导率深度成像的数据表来源：

点击浏览公式

式中，零阶和一阶贝塞尔函数的无限积分(汉克尔变换)利用数值滤波算法(Guptasarma120点和140点)^{[18, 19]}计算。

设定固定翼航空电磁探测系统保持平稳飞行状态，飞行高度 在80~180 m等间隔选取101个，半空间模型的电导率在0.001~20.000 s/m按等对数间隔选取55个，形成101×55组半空间探测模型。利用上述正演算法，计算14个时间道(0.01~10.00 ms)各模型对应的双分量B场响应。按时间道划分，建立14个B_x-B_z-电导率-飞行高度数据表，表中任一点为由对应飞行高度及半空间模型确定的x和z分量B场响应，记为(B_x,B_z,σ,h)。图 2给出了第k(1≤k≤14)道的部分B_x-B_z-电导率-飞行高度数据表示意图，其中电导率σ₁ < σ₂ < σ₃，飞行高度h₁>h₂>h₃>h₄>h₅>h₆>h₇>h₈>h₉。

和为实际飞行测量的第k道x和z分量的B场响应。图 2 第k道Bx-Bz-电导率-飞行高度数据表示意图Fig.2 Diagram of Bx-Bz-conductivity-altitude for channel k

图选项

传统的总场查表法与单分量查表法电导率深度成像^{[20, 21, 22]}采用总场响应或单分量响应进行匹配搜索查询视电导率，在搜索过程中需要提供实际的飞行高度以选择所需的数据表，而实际飞行探测过程中很难获取准确的飞行高度，给电导率的查询带来一定困难。

B_x-B_z分量联合查表法建立了B_x-B_z-电导率-飞行高度数据表，将飞行高度引入数据表中，在无法准确提供高度信息时不但可以快速获得视电导率，还可以得到对应的高度信息。同时采用x和z两个分量的电磁响应进行匹配搜索，不但增加了搜索时的约束条件，而且综合利用了x和z两个分量电磁响应的特性，从而提高成像的精确度。总场查表虽然也综合利用了多分量的响应，但对于两组不同的B_x、B_y和B_z分量，只要满足

就会得到相同的电导率查询结果。因此，总场查表虽然利用了多个分量，但由于进行了二对一的转换，损失了多分量测量信息。

在第k道对应的表中(图 2)查表搜索和，确定其在数据表中的具体位置，将其所在四边形的顶点分别记为(B_x1，B_z1,δ₁，h₁)，(B_x2，B_z2,δ₂，h₂),(B_x3，B_z3,δ₃，h₃),(B_x4，B_z4,δ₄，h₄),放大如图 3所示。其中，h₁=h₄，h₂=h₃，δ₁=δ₂，δ₃=δ₄。由于数据表的间隔较小，故具有准线性的特征，根据图 3中的几何关系，利用式(7)、式(8)线性插值获得第k道对应的视电导率以及飞行高度。反复进行上述步骤，即可得到电磁响应对应各道的视电导率及高度信息。

图 3 插值算法示意图Fig.3 Diagram of interpolation algorithm

图选项

点击浏览公式

其中：q₁——q₄为线段长度。

式中，t_k为第k个时间道对应的时间。在相邻2个时间道深度变化的基础上定义成像深度，即第k个时间道的成像深度是由上一时间道的成像深度加上利用趋肤深度公式计算得到的相邻两个时间道深度变化值的绝对值，从而获得第k个时间道对应的成像深度D_k^[11]：

点击浏览公式

综上，B_x-B_z分量查表法联合电导率深度成像的步骤如下：

1)利用正演算法，建立B_x-B_z-电导率-飞行高度数据表表格；

2)通过磁场双分量联合查表与插值算法，确定视电导率；

3)利用扩散深度公式得到视深度，并计算成像深度，最终导出双分量联合电导率深度成像结果。

设定固定翼航空电磁系统的z分量发射线圈高度为h=120 m，三分量接收线圈位于(-70,0,70) m处，系统参数同前。针对不同层状大地模型对B_x-B_z分量联合查表法、总场查表法与单分量查表法CDI进行对比研究。

以三层大地模型为例，选取了K型、H型典型地电模型作为参考模型对算法进行对比分析。表 1为K、H型大地模型参数。

图 4给出了K、H型两种三层大地模型采用B_x-B_z分量联合查表法、总场查表法与单分量查表法的CDI结果。可以看到，总场查表法成像结果与B_x和B_z单分量的查表结果一致，大致给出了三层大地模型的变化趋势，而B_x-B_z分量联合查表算法对K型大地模型(图 4a)中间高导层的电导率和深度的成像结果与理论模型接近，即使H型大地的高阻层(图 4b)，B_x-B_z分量联合查表算法也获得了优于总场与单分量查表法的成像结果。而由于查表算法是利用半空间大地模型的电磁响应作为数据表，而观测响应是含有异常体的三维大地模型;因此，利用观测值在半空间模型表中查询时，可能会出现部分点不在数据表中，使得成像结果不太光滑。

图 4 三层典型大地模型CDI结果Fig.4 Results of four CDI methods on three-layer models

图选项

综上，总场查表法虽然利用了B_x和B_z两个分量进行电导率深度成像，但是由于进行了二对一的转换，仅获得了与单分量查表法近似的成像结果，损失了多分量测量获取地下信息的能力。而B_x-B_z分量联合查表算法充分利用B_x和B_z两个分量的信息，进行联合查表，获得了与真实大地模型更为接近的地下信息，提高了成像精度。

为评价B_x-B_z分量查表法联合电导率深度成像算法在实际中的应用，仿真构建了准二维大地模型，由500个5层大地模型相互连接构成。如图 5所示，在测点100~250处的地下埋深2 m有一个厚度为5 m、电导率为0.2 s/m的异常，同时在埋深57 m处有一个厚度40 m、电导率0.2 s/m的异常，周围基岩的电导率为0.02 s/m。对上述大地模型进行正演计算获得其电磁响应剖面数据，并加入均方根相对误差为5%的白噪声，仿真野外实测数据。

图 5 理论大地模型Fig.5 Pseudo-2D theoretical model

图选项

为进一步对比成像效果，分别采用B_x分量查表法、B_z分量查表法、总场响应查表法及B_x-B_z分量联合查表法对数据进行电导率深度成像。图 6分别展示了上述四种方法对于含噪准二维数据的CDI结果。

a.Bx分量查表法;b. Bz分量查表法;c. 总场响应查表法;d. Bx-Bz联合查表法。图 6 四种方法的电导率深度成像结果对比图Fig.6 Results of four CDI methods on Pseudo-2D model

图选项

可以看到，B_x分量查表法(图 6a)、B_z分量查表法(图 6b)与总场查表法(图 6c)的准二维含噪数据的CDI结果都较好地反映了埋深2 m的导电异常体，屏蔽了其下埋深57 m处的导电异常体，而图 6d给出的B_x-B_z分量联合CDI结果却较为清晰地揭示了导电薄层下的低阻异常体，与理论模型一致。笔者计算了测点175处电导率深度成像结果的相对误差平均值得知，B_x分量查表、B_z分量查表、总场查表法以及B_x-B_z分量联合查表法的相对误差分别为22.07%、21.65%、21.84%、14.72%，即本文研究算法的成像精度较单分量电导率深度成像提高了7%左右。

1)本文研究的基于磁场双分量的联合电导率深度成像算法，充分利用了多分量电磁探测数据，取得了优于总场与单分量查表法的成像结果。

2)一维大地模型的仿真结果表明该算法提高了对导电薄层的识别能力；准二维大地模型加噪数据的成像结果进一步证明在低阻覆盖层存在的情况下，联合电导率深度成像算法能够增强对其下导电异常体的识别能力，综合成像精度较单分量电导率深度成像提高了7%左右，为固定翼航空瞬变电磁资料解释提供了新思路。

吉林大学地球信息探测仪器教育部重点实验室航空电磁组的全体成员对此文进行了帮助和指导，在此表示衷心的感谢！