长方体全张量磁梯度的两种正演公式对比

钟炀,管彦武,石甲强,肖锋

吉林大学地球探测科学与技术学院，长春130026

0 引言

磁法勘探探测技术的发展过程大致分为总磁场标量测量、磁场三分量测量以及全张量磁梯度测量3个阶段[1-2]。当前航空磁测主要利用组合式磁梯度计进行大工区的全张量磁梯度测量，全张量磁梯度数据较磁标量数据和三分量矢量数据具有抗干扰能力更强、分辨率更高、利于分辨场源磁化方向及小异常体边界特征等优点[2-7]。长方体是正演计算中常用的三度体模型[8-9]，因此推导长方体全张量磁梯度理论表达式具有重要的理论意义。

针对长方体模型磁场正演问题，近年来许多专家学者开展了大量研究。郭志宏等人首次指出此前未曾被人发现或重视的长方体ΔT场表达式中场值无法计算的解析奇点问题，并从直立长方体模型引力位出发，在求解积分过程中多次运用凑因式法及等价变换法消除导致奇点的多项式，推导出东南下坐标系的长方体ΔT场及其梯度场无解析奇点表达式[9]，但其推导过程在消除奇点的同时也增大了计算量；管志宁在郭志宏等人的研究基础上给出了长方体磁场三分量的计算公式[8]；为简化长方体模型正演计算公式，骆遥等人引入欧拉方程对长方体磁场理论表达式进行重新推导，得出了形式更加统一、简洁的长方体ΔT场及磁场三分量在上半无源空间的无解析奇点理论表达式[10]；针对前人理论公式仅能计算长方体模型上半无源空间，而无法适应起伏地形条件下的正演问题，匡星涛等人利用变量替换法，并单独考察让磁场表达式无意义的奇点处的磁场值，推导出适用于整个无源空间的长方体ΔT场无解析奇点表达式[11]；随着航空全张量磁梯度测量仪器不断取得重大突破，全张量磁梯度正反演理论也得到了发展。其中规则形体的全张量磁梯度正演算法受到广泛的重视，徐熠通过求解直立长方体模型引力位二阶偏导数，将其带入泊松公式，求解出长方体全张量磁梯度表达式[12]；干博、吕文杰分别根据骆遥等人推导出的长方体磁场三分量表达式，给出长方体模型全张量磁梯度表达式[13-14]；修春晓在骆遥等人推导出的长方体磁场三分量表达式基础上，给出基于体剖分模型的全张量磁梯度公式及计算效率较高的基于网格节点物性的长方体全张量磁梯度理论表达式[15]。

考虑到地球物理勘探常在场源外的上半空间进行，故本文在管志宁给出的东南下坐标系和骆遥等人给出的北东下坐标系的磁场三分量表达式的基础上，进一步推导了全张量磁梯度的计算公式，并进行了对比分析。沿用郭志宏及骆遥等人的长方体模型参数[10-16]，给出这两种计算公式的全张量磁梯度正演结果，它们完全相同。在三维物性反演中，反演的地质体往往比较复杂，通常需要剖分地下模型空间，利用正演公式进行迭代计算[17-19]，进而拟合观测异常，达到反演的目的。例如将地下地质体剖分为100×100×10的长方体组合模型，那么在正演公式中每增加一步计算，会增加10万次的计算[20]。因而有必要选择更为简洁的计算公式，为模型正演计算节约时间，提高反演效率。

1 东南下坐标系长方体全张量磁梯度计算公式

全张量磁梯度是磁场强度3个分量(Bx，By，Bz)在空间直角坐标系X、Y、Z的3个坐标轴方向的变化率。即：

由于磁场强度的旋度为零，所以全张量磁梯度矩阵为对称矩阵，即式(1)中，Bxy=Byx，Bxz=Bzx，Byz=Bzy；在无源空间中，磁标量位满足拉普拉斯方程，即U=0，故有Bxx+Byy+Bzz=0，所以在全张量磁梯度的9个分量中，只有5个独立分量[3]。在实际应用中，为了便于表达全张量磁梯度或验证磁标量位满足拉普拉斯方程，通常需要计算出上三角矩阵中的6个磁梯度分量。

首先，建立如图1所示东南下空间直角坐标系(X′,Y′,Z′)[8]，其中，X′正轴指向地理东方向，Y′正轴指向地理南方向，Z′轴铅直向下。

图1 东南下坐标系长方体模型Fig.1 Cuboid model in east-south-down coordinate system

在该坐标系中，直立长方体上半无源空间磁场三分量理论表达式为[16]:

(2)

(3)

(4)

对式(2)～(4)分别沿x′，y′，z′方向求偏导数，求得在东南下坐标系中的长方体全张量磁梯度表达式，即公式(5)～(10)：

(5)

(6)

(7)

(8)

(9)

(10)

2 北东下坐标系长方体全张量磁梯度计算公式

北东下坐标系，即NED(North East Down)坐标系，在航空航天、测绘和勘探等领域中经常使用该坐标系[24]。建立NED坐标系如图2所示。

图2 北东下坐标系长方体模型Fig.2 Cuboid model in north-east-down coordinate system

骆遥等人给出在该坐标系下直立长方体上半无源空间磁场三分量理论表达式[10]：

(11)

(12)

(13)

骆遥等人根据长方体重力场及其梯度满足构造指数为1的欧拉方程，并将Okabe给出的长方体重力场垂向一阶导数[25]带入其中，得到不含分子分母同时为零的解析奇点的引力位二阶导数，并将其带入泊松方程中，得到长方体无解析奇点的磁场三分量理论表达式。但推导出的引力位二阶导数中反正切项中仍存在分母为零的情况，再利用单变量阶梯函数在对称区间的三重积分恒为零的性质，对分母为零的反正切函数进行换元，求得形式更加整齐、简洁的长方体无解析奇点的磁场三分量表达式(11)～(13)。以上三式将模型角点到观测点的距离移至积分限中，利于简化计算，从而在编程时减少冗余的计算步骤。

对以上三式分别沿x，y，z方向求偏导数，即保持积分上下限不变，对式中的ξ，η，ζ求偏导数。求得在北东下坐标系中的长方体磁场全张量理论表达式，见式(14)～(19)。对比吕文杰及干博所推导的长方体全张量磁梯度公式[13-14]，他们没有将模型角点到观测点的距离移至积分限中，故所求长方体全张量磁梯度公式显得冗长。徐熠及修春晓所推导的长方体全张量磁梯度公式[12,15]中，均存在未进行因式分解至最简结果的项。

(14)

(15)

(16)

(17)

(18)

(19)

为方便描述，简称公式(14)～(19)为算法二。

这两种算法都是先计算出长方体引力位二阶导数，再根据泊松方程计算磁场三分量，进一步求方向导数，获得全张量磁梯度表达式。但是在求解引力位二阶导数的过程中，两者采用了不同的方法。算法一通过等价变换及凑因式法避免了前人[21-22]在求解积分过程中引入的奇点问题；算法二则是引用了Okabe给出的长方体重力场公式[25]，并带入到构造指数为1的欧拉方程中求解直立长方体引力位二阶导数。因而两种算法得到了形式不尽相同的长方体磁场三分量表达式，进而全张量磁梯度的表达式也简繁不同。

算法二即为北东下坐标系下长方体全张量磁梯度理论表达式，对比算法一。两种算法需要的参数个数相同，但算法二的公式形式相对简洁；算法二采用的坐标系更符合常规，且磁化偏角与地磁偏角的定义一致，便于使用；算法二设定的积分限在简化计算方面的优势也继续得到体现。为进一步对比两种算法计算量的差异，以同一长方体模型，计算观测面上同一点为例，统计了两种算法中每个磁梯度分量的计算量(表1)。其中，算法二全张量磁梯度公式加减法计算总量不到算法一的四分之一；乘除法计算总量约为算法一的一半。

表1 两种算法全张量磁梯度计算次数

Table 1 Calculation times of full tensor magnetic gradient with two algorithms

FTMG分量加减法计算次数乘除法计算次数算法一算法二算法一算法二Bxx1926013284Bxy2143014544Bxz1933012644Byy1906013074Byz1893012444Bzz116609474总计1 094270751370

3 长方体全张量磁梯度模型计算

为检验在两种坐标系下导出的长方体全张量磁梯度理论表达式的正确性，对算法一、算法二取完全相同的长方体模型进行正演计算，模型参数沿用郭志宏、骆遥等人的参数[10、16]。图3a～3f为利用算法一计算得到的全张量磁梯度等值线图，为方便对比，已转换到常用的北东下坐标系中显示。坐标转换方式为：

模型参数为：磁化方向I=50°，A=30°;磁化强度M=1 A/m;地磁场方向I0=60°，A0=10°;长方体中心埋深x0=10 000 m，y0=10 000 m，z0=3 000 m;长方体长a=8 000 m，宽b=4 000 m，高c=1 000 m;网格间距为100 m×100 m;计算高度z=0 m;等值线单位为nT/m。

图4为模型在北东下坐标系中利用算法二计算得到的全张量磁梯度等值线图。对比两图，可见同一长方体模型，两种算法的全张量磁梯度等值线图形态完全一致。为进一步对比两坐标系下全张量磁梯度差异，用图4减去图3所对应的张量值，并绘制成等值线图(图5)。

图5表明，两种计算公式得到的全张量磁梯度差值的量级均在10-17～10-15nT/m之间，相比于当前精度较高的航空全张量磁梯度测量的实际观测精度为10-2nT/m[26]，这两种算法的计算结果可视为无差别。此外，为对比计算效率，笔者建立了一个地下剖分为100×100×10的长方体组合模型，测区范围20 000 m×20 000 m，点线距均为100 m，X和Y方向网格点数均为201个。所用计算机CPU为Intel(R)Core(TM)i7-8750H，内存(RAM)8.0GB，Matlab版本2018a，算法一耗时1 896.5 s；而算法二耗时1 703.7 s。比较程序计算用时，算法二较算法一节省约10.2%的时间，可见两者计算量相差较大。主要原因是两种计算方法各张量表达式差异，算法二表达式相对简洁。

图3 由算法一计算的长方体全张量磁梯度等值线图Fig.3 Contour map of full tensor magnetic gradient of cuboid with algorithm 1

图4 由算法二计算的长方体全张量磁梯度等值线图(模型参数与图3相同)Fig.4 Contour map of full tensor magnetic gradient of cuboid with algorithm 2

图5 图4与图3对应张量之差等值线图Fig.5 Contour map of the difference between Fig.4 and Fig.3

4 结论

(1)分别推导了东南下和北东下两种坐标系下长方体模型的全张量磁梯度理论表达式，经验证两种算法的计算结果完全相同。

(2)从推导过程及公式形式上对比了两种算法，它们采用了不同的推导过程，故公式形式不尽相同，但获得了相同的计算结果；算法二中坐标系和磁偏角的定义更符合常规，便于使用；算法二计算效率高于算法一，在模型试算中，算法二节省10.2%的计算时间。