在进行航磁测量时，飞机不可能始终处在同一高度的平面.观测面起伏变化导致测得的重磁场形态复杂化，对解释工作造成干扰.而对于磁场进行处理转换往往都是基于平面数据，因此对于磁场测量数据由起伏观测面转化到水平面的曲化平处理也就显得尤为重要.

在位场数据曲化平处理方面，等效源法(杜维本，1982)、有限元法(程振炎，1981)和边界元法(Xu，2001)一般需要解线性方程组，对于大数据量的实际资料由于需要计算机有很大的内存，并且耗费时间，这些方法并不实用.Taylor级数法(刘金兰等，2007)由于计算高阶导数放大异常高频成分，该方法也不宜处理起伏变化较大的观测面上的位场数据.迭代法曲化平(徐世浙，2006)计算过程简单，在条件较好时取得了很好的效果.但由于迭代曲化平过程中间最重要一步是将曲面上的数据下延到地形最低点所在的平面上，所以，迭代法曲化平要求必须是场源较深的弱地形条件(姚长利等，2012)，因此受地形的起伏变化的程度、下延到地形最低点的高差等影响较大.根据位场的性质则可以将以上方法应用于磁梯度数据的曲化平处理(李海侠等，2013；谢汝宽等，2015).

以上方法主要应用在单参量(磁异常分量或其梯度)的曲化平中，现在航磁梯度测量逐渐展开，磁法勘探也由过去的单参量测量发展为目前的多参量测量.在进行磁场分量测量时，也可以同时测量其梯度.梯度测量的意义在于它能弥补磁场分量测量空间分辨能力较低的缺点，对于浅表小异常有很好的反映，有利于对磁场进行三维描述分析.在能够实现多参量测量的年代，研究如何综合利用实际测量的多参量数据进行有效的曲化平也就变得很有必要.

姚长利曾利用测量磁场及垂向梯度进行曲化平处理.提出由观测曲面上的位场及其垂向一阶导数，利用样条函数计算该曲面上的位场垂向二阶导数及垂向三阶导数(姚长利等，1997)，再由位场的三阶Taylor展开式进行曲化平换算，其计算速度快，计算结果比较稳定.充分发挥了垂向梯度反映地质体磁场垂向变化特点的作用.为综合利用位场及其梯度进行曲化平处理提供了很好的思路.上述方法中的垂向梯度，需要事先获得.

曲化平的直接求解法(管志宁，2005)是在已知观测面上的磁场及其法向导数的情况下，利用位函数的基本积分公式，直接计算出观测面以上空间中任意点的磁场.如何基于这个法向导数，直接进行曲化平，正是本文的研究工作.

位函数的基本积分公式为

(1)

式(1) 表示由边界位场值求解调和域中的位场值.已知观测面上的磁场及其法向导数，可以直接利用式(1) 计算出观测面以上空间的磁场值.直接求解法中曲面上的法向导数可以通过实测的三轴向梯度值转换得到.

边界元法曲化平是通过利用积分节点对积分区域的立体张角ω_j建立线性方程组，需要解线性方程组求出磁场在曲面上的法向导数，对于大数据量的实际资料需要计算机有很大的内存，并且耗费时间.但边界元曲化平方法提供了利用剖分三角元和高斯数值积分的方法，给出了计算积分方程的系数的方法.如果有实测的梯度数据，那么法向导数可以由实测的梯度转换得到.

在磁场及其梯度已知的情况下，将得到一种简便的综合利用磁场分量(ΔT)及其梯度进行曲化平的方法.只要获得了高精度的曲化平磁异常结果，即可对此平面上的磁场进一步在频率域内计算出平面上的梯度，由此可同时实现磁异常和磁梯度的曲化平处理.可以看出该方法关键在于利用梯度数据对ΔT数据进行曲化平处理.

边界单元法的基本方法为：

(1) 用三角单元对地表边界Γ_s进行剖分，使每个单元Γ_e的地形近似平面，假定单元内的u和是线性变化的，用线性等参单元进行一次插值，通过将形函数与单元上各节点的值相乘得到单元内各高斯积分点的值.公式为

(2)

其中N_j、N_k、N_m是形函数，公式为

(3)

(2) 单元积分，公式为

(4)

其中：

(5)

将这些积分使用高斯数值积分进行，用7点高斯积分公式：

(6)

其中r_iq是i点到单元Γ_e上高斯积分点q的距离，cos(r_iq, n_Γe)是r_iq与Γ_e的外法线方向n_Γe的夹角的余弦，w_q是加权系数，Δ是三角形面积，7点的高斯积分系数N_j、N_k、N_m和加权系数w_q可以见表 1：

表 1 (Table 1)

表 1 三角元上7点高斯积分系数和加权系数 Table 1 The Gauss integral coefficient and the weighted coefficient of 7-points on the triangular element

表 1 三角元上7点高斯积分系数和加权系数 Table 1 The Gauss integral coefficient and the weighted coefficient of 7-points on the triangular element

(3) 最后将各单元积分相加，公式为

(7)

其中F_ij、D_ij是节点j周围诸单元f_ij、d_ij之和，共n个节点.

本文通过修改Matlab一个内置函数的系数，假定一个曲面函数做为航磁测量中飞机高度组成的曲面，函数为

(8)

在x和y取值范围为(-10⁴，10⁴)时，取网格间距为200，其曲面形态如图 1，由图 1可以看出该曲面模型z有三个极大值的坐标为：(0，4400，2870)、(5600，0，6460) 和(-1800，-1400，3000)，z最小值坐标为(5000，800，-5200).

图 1

Figure 1

图 1 假定的曲面示意图 Figure 1 Sketch map of the assumed curved surface

组合地质体模型为一个直立长方体和三个球体，其平面投影图如图 2.三个球体分别位于图 7所示的高度曲面极大值的下方，坐标为：(0，4400，2400)、(5600，0，5800) 和(-1800，-1400，2600)，三个球体半径均为100，保证了每个球体在任一方向都不“冒出”曲面.直立长方体长×宽×高为：5000×5000×1000 m，上顶面设置为-5500 m，低于图 1中曲面最小值.曲面最小值基本处于长方体南边缘上，距离较近.

图 2

Figure 2

图 2 模型平面投影图 Figure 2 Map for plane projection of model

图 3

Figure 3

图 3 正演计算ΔT和梯度的三分量等值线图 (a)ΔT；(b)北方向梯度；(c)东方向梯度；(d)垂向梯度. Figure 3 The contour map for ΔT and three component of gradient in forward calculation (a)ΔT; (b)Northward gradient ; (c)Eastward gradient ; (d)Vertical gradient .

图 4

Figure 4

图 4 用曲面函数求曲面的法向导数图(a)和高度为6600 m的理论ΔT等值线图(b) Figure 4 The normal derivative determined from surface function (a) and the ideal magnetic anomaly ΔT at the height of 6600m (b)

图 5

Figure 5

图 5 三角元向“下”的法线方向 Figure 5 Normal direction of triangle element directed downward

图 6

Figure 6

图 6 求节点上的法方向 Figure 6 Method of determining the normal direction on the node

图 7

Figure 7

图 7 利用梯度计算的法向导数(a)以及曲化平后的ΔT等值线图(b) Figure 7 The normal derivative determined from gradient (a) and the reduction magnetic anomaly ΔT (b)

设置地磁倾角I=45°，地磁偏角D=0°；所有模型磁化倾角i=45°，磁化偏角d=0°；模型磁化率均设定为200×10^-3(SI)，可以得到ΔT及其三个方向梯度、和，如图 3.

从图 3可以看出由于长方体南边缘距离飞行高度曲面最小值处很近，在该位置处所引起的异常达到了2000 nT，在画等值线图时，1号球体和3号球体引起的异常在图中只有微弱反映，而2号球体的磁异常在等值线图中则近乎湮没.从三个梯度的等值线图上看，三个球体的异常均得到了凸显，尤其是水平东向梯度；对于长方体在曲面最小值处引起的异常则只有较弱的反映.

在曲面函数z为已知的情况下，计算z在x、y方向导数分别为和，则曲面的法向量为.对其求单位向量n_j后，利用梯度三个轴向分量值可以求得曲面法方向的磁场梯度—即磁场的在曲面上的法向导数，可以表示为

(9)

式中：为单元上节点j上的外法向导数；n_j为节点上的单位法向量；、和为单元节点j上实测的x、y、z三个轴向的梯度，T表示求向量转置.

图 4a为利用梯度分量以及曲面函数的导数计算得到的曲面上的法线方向的梯度，将其作为法向导数的理论值.由于函数曲面极大值为6478.4，本文在此计算了高度为6600时ΔT值作为后期曲化平处理的理论值，见图 4b.

在实际测量中，要求出高度所构成的曲面函数是比较困难的.因此，整个求解问题就简化为利用测量的梯度计算曲面上磁场法向导数的问题，即确定曲面节点处法方向的问题.

节点上法方向的确定有多种方法，本文提供一种简单有效的方法.因为在积分过程中需要用到三角元的法方向单位向量及面积，因此可以利用每个三角元的法方向及面积求节点上的法方向.

每个三角元的法方向的确定很简单，使用三角单元三个顶点坐标(x_j、y_j、z_j)、(x_k、y_k、z_k)和(x_m、y_m、z_m)确定的平面单元，求其向“下”的单位法向量(向“下”的法方向即磁场调和域的外法线方向).

在图 5中所示的坐标系中三角元的单位法向量n_Γe可以表示为

(10)

(11)

式(10) 和(11) 中Δ为三角元的面积，det表示求行列式的值，A、B、C分别为

(12)

那么可以用以下方式求得节点处单位法向量，公式为

(13)

(14)

式(13) 表示求向量a的单位向量.由于三角元的法向量(detA detB detC)的模量为其面积的2倍，那么向量a就是节点j周围6个三角元的法向量的向量和，即周围每个三角元的单位法向量乘以其面积后求和，面积为每个三角元单位法向量的加权系数，再求向量a的单位向量.式(13) 和(14) 实际上是利用节点周围每个三角元的单位法向量及每个三角元的面积为权重计算出节点处的单位法向量(图 6).

在进行曲化平处理时，网格距设定为200 m，扩充边界为-20000~+20000 m(本文并未研究用何种方法进行扩边处理会更好一些，只是将正演计算的面积进行了放大).

利用本文提出的方法计算得到法向导数值，可见图 7a，可以看出其形态与图 4a中的理论曲线形态是一致的.与法向导数(图 4a)理论值进行比较，计算法向导数的误差范围为-0.0052~+0.0175 nT/m，均方差为±5.7565×10^-4 nT/m，可以看出该方法的有效性.

利用计算得到的法向导数曲化平到高度为6600的ΔT可见图 7b.对比图 7b和图 4b可以发现，曲化平后的曲线形态与理论曲线形态相同.曲化平后数据与理论值的误差范围为1.0201~3.1930 nT，均方差为2.1413 nT.所以曲化平的结果相当好.计算最大值为40.5803 nT，理论最大值为38.0227 nT，误差为2.5576 nT；计算最小值为-15.2796 nT，理论最小值为-17.2302 nT，误差为1.9506 nT.曲化平后的曲线两个极值和理论数据的两个极值误差很小.

通常认为在网格距越小，则拟合曲面的程度越好，那么曲化平的效果也会越好.在试验了多个不同的网格距和不同的边界扩充后，本文认为和边界扩充的程度也有一定关系，只有在无限扩充边界时网格距越小曲化平将越接近理论值.排除以上因素，从理论上讲只要曲面上的法向导数计算正确，即曲面的法方向计算正确，那么整个曲化平处理是没有问题的.

利用单参量进行曲化平处理或多或少都存在一些问题.本文通过模型计算证明了所提出的利用边界元方法中剖分的三角元的法向量计算曲面上节点处磁场的法向导数的方法是有效的.并利用位函数基本公式—格林公式—直接进行曲化平处理，该方法充分的利用了磁场梯度所反映的磁场在各个方向上变化特征的作用，因此在非弱地形条件下的中—高山区进行梯度测量可以实现高精度的曲化平处理结果.该方法原理清晰，程序实现简单.

本文利用的是正演计算理论数据，但在实际的梯度测量过程中由于飞机姿态的影响，实测的垂向梯度并不一定对应曲面的法向导数，如果把实测的垂向梯度看作曲面的法向导数，将磁场值及其垂向导数代入积分公式实现向上延拓，计算也就会有很大误差.在利用梯度前需要将测量的梯度转换到既定的坐标系中，即要对测量梯度进行姿态改正(梁建，2016)后再计算法向导数.

从理论上讲网格距越小，剖分的三角元对曲面的拟合程度越好，但整个处理过程的计算时间也会相应的加长.另外对磁场和梯度数据以及高度数据的扩边处理也需要进一步研究，使用不同扩边方法可能会降低时耗.网格距与扩边之间的关系控制好，可以使得化平之后的磁场能达到最佳效果.