基于位场曲面延拓的航磁成图数据精度评价方法研究

张义蜜， 熊盛青， 于长春， 王万银

1 长安大学重磁方法技术研究所， 长安大学地质工程与测绘学院，长安大学西部矿产资源与地质工程教育部重点实验室， 西安 710054 2 海洋油气勘探国家工程研究中心， 北京 100028 3 中国自然资源航空物探遥感中心， 自然资源部航空地球物理与遥感地质重点实验室， 北京 100083 4 加拿大纽芬兰纪念大学， 地球科学系， 圣约翰斯 AlB3X5

0 引言

航磁测量数据精度评价是非常必要的，但影响数据精度评价结果的因素有很多，Reeves(1993)指出利用切割线进行调平可能导致测线上的磁异常失真，若以调平后的数据进行精度评价可能会带入误差；现行的航空磁测技术规范中指出航空磁测总精度是由航空磁力仪系统测量和定位误差及各项校正不充分或不准确等误差组成，经过校正后进行评价得到航空磁测总精度，再经过精细调平等处理后评价得到成图总精度(中华人民共和国国土资源部，2010；薛典军，2001).分析得出，随着磁力仪灵敏度的提高，测量方法的改进，飞行高度是影响测量总精度的最主要原因，特别是测线和切割线交叉点处的高度偏差，即航磁测线和切割线测点的位置偏差及调平等处理会对评价结果造成一定的影响，在地形起伏大的山区影响更大(熊盛青等，2009；Domzalski，1957).前人关于航磁成图数据精度评价方法可以归纳为航磁数据经过精细调平等处理后，计算切割线(重复线)与测线磁场差值的均方根偏差，以此作为成图总精度，具体方法有《航空磁测技术规范》中提出的方法(下简称航磁规范法)(中华人民共和国国土资源部，2010；陈斌等，2010；李军峰等，2014；李志鹏等，2018；刘玲等，2018；崔志强等，2019；乔中坤等，2020)和重复线法(徐东礼等，2014，2016；梁建等，2018).

航磁数据调平是航空物探数据预处理的关键性步骤，对数据处理和解释具有重要意义.目前调平方法可分为切割线调平(Foster et al., 1970; Green, 1983; Mittal, 1984; 胥值礼等, 2010)、微调平(Minty, 1991; Ferraccioli et al., 1998; 骆遥等, 2012a,b)、伪切割线调平(朱月娥, 1994; 毛玉仙, 1995; Huang, 2008; White and Beamish, 2015)以及其他方法(Nelson, 1994; Pilkington and Roest, 1998; Fedi and Florio, 2003; Mauring and Kihle, 2006).调平方法众多，调平后计算得到的磁场成图数据精度也不尽相同，这将影响航磁成图数据精度评价的客观性和准确性.

综合以上，前人关于航磁成图数据评价方法存在以下两点需改进的地方：(1)舍弃磁场差值大于3倍预设精度或高差大于100 m的交叉点后，交叉点处高程不一致问题依然存在，精度评价过程中需解决交叉点处高程不一致问题；(2)因调平方法不同和调平次数不同导致最终的航磁成图精度不同，从而不能客观评价成图数据的质量情况.针对以上问题，本文基于位场理论，利用位场数据非规则网延拓方法，提出一种基于曲面延拓的航磁数据精度评价方法，改善了已有航磁成图数据精度评价方法存在的不足，从而可以更加客观、有效地评价航磁成图数据的精度.

1 原理与方法

1.1 航磁规范成图数据精度评价方法分析

航空磁测总精度是由航空磁力仪系统的测量和定位误差及各项校正不充分或不准确等误差的总和.在经过各项校正后，航空磁测总精度采用计算切割线与测线交点上磁场差值的均方根偏差σ来衡量，磁场水平精细调整前计算出的σ作为主要由定位引起的磁场测量误差.其计算公式为

(1)

其中，δi为第i(i=1,2,3,…,n)个交叉点处测线与切割线磁力异常差值；n为交叉点个数.σ和δi的单位均为nT.

在完成正常场、磁日变校正和整架次或每条测线粗略水平调整后，按式(1)计算调平前的总误差σ值.计算σ时，允许舍掉磁场梯度偏差大(600 nT·km-1，交叉点处磁场变化>3σ者)、磁场偏差>3σ或高程差大于100 m的交叉点，其余的交叉点均应参与σ的计算.在完成精细调平后应按上述舍点要求再计算σ，作为全区测量的成图总精度值.

以上为航磁规范成图数据精度评价方法，但若舍弃磁场偏差>3σ后，测线和切割线交叉点处磁场值为ai(i=1,2,…,n)，切割线交叉点处磁场值为bi(i=1,2,…,n)，满足|ai-bi|≤3σ(i=1,2,…,n)，则

(2)

即仅舍掉磁场差值>3σ的交叉点后，磁场总精度已经小于2.1σ.

舍弃磁场差值>3σ的交叉点后，若ai-bi≤σ占比λ，σ

(3)

若不“舍点”，测线和切割线的点位高差引起的磁异常偏差则会引入到精度评价结果中，从而影响最终结果的准确性.

1.2 基于位场曲面延拓精度评价方法

基于位场的等效性，通过设立等效源代替真实场源，拟合观测数据建立等效源模型，进而利用此模型进行位场数据的处理和转换.本文利用此思想，由经典场论可知，在无源空间，位场延拓可以表示为Laplace方程边值问题(Courant and Hilbert，1962)

(4)

由(4)式可得到半空间Dirichlet的解为(Blakely，1996)

(5)

其中，Uc(x,y,z)和Up(ξ,η,z0)为场值，A为积分核函数.如果Uc(x,y,z)和Up(ξ,η,z0)均为离散分布，并假设离散数据分别为M、N.则式(5)可以表示为

(6)

其中

(7)

式(7)常使用矩形积分公式进行离散化(刘福香等，2021).

式(6)采用矩阵形式可表示为

Uc=G·Up,

(8)

其中，Uc为观测数据序列；Up为等效平面数据序列；核函数矩阵G表示连接观测数据和等效平面数据的正演算子.通过共轭梯度法求解方程(8)，即可得到等效平面数据序列Up.

刘芬等(2019)研究等效平面高程应大于等于0.4倍点线距dx(若无点线距，则为观测数据点距的平均数)，若已知观测面高程zi(i=1,2,…,M)和待延拓面(点)高程zc,i(i=1,2,…,K)，即可得到等效平面高程z0：

z0=min[min(z),min(zc)]-0.4dx,

(9)

式中，函数min表示求最小值运算.

关于该方法计算效率问题，本文采用滑动窗口和OpenMP并行编程技术共同提高计算效率.在保证计算精度的条件下，刘芬等(2019)研究得到窗口半径一般选择20倍的延拓高度以上.而OpenMP是一个跨平台的多线程实现，主线程(顺序的执行指令)生成一系列的子线程，并将任务划分给这些子线程进行执行.这些子线程并行的运行，由运行时环境将线程分配给不同的处理器，因此可大大节省计算耗时.

通过以上步骤即可求得等效平面位场值，从而能够实现磁场数据向上延拓、向下延拓计算.通过将测线磁异常延拓至切割线测点(交叉点)处得到延拓磁异常，将延拓磁异常和切割线测点(交叉点)磁异常利用公式(1)计算航磁数据精度，后续将这种精度评价方法简称为曲面延拓法.

2 理论模型测试

为了对比航磁规范成图数据精度评价方法和本文提出的曲面延拓法的评价结果，设计理论模型，对比各个方法评价结果的差异性.在切割线测点处，实测值和延拓值之间的差值由航磁数据测量误差和延拓误差两部分构成.王万银等(1991)指出曲面延拓算法自身误差远小于测量误差；蒋涛(2018)在航空重力向上延拓外部估计中指出解析延拓误差可以忽略不计，本文首先测试非规则网延拓算法(曲面延拓)引起的误差量级，确保其不影响最终的数据精度评价结果.

2.1 非规则网延拓误差评价

本文设计了理论模型测试非规则网延拓引起的误差量级.组合模型体由五个直立六面体组成，空间位置和物性参数如表1所示.组合模型体在起伏观测面(图1a)引起的磁力异常如图1b所示.

图1 起伏观测面及磁力异常图(a) 观测面高程等值线图；(b)模型体在观测面高程处引起的磁力异常等值线图(注：图中白框实线为模型体水平展布位置；两条东西向和南北向剖线穿过异常幅值变化剧烈处).Fig.1 Contour map of observation surface fluctuation and magnetic anomaly(a) Observation surface elevation contour map; (b) Magnetic anomaly contour map caused by the model body at the observation surface elevation (Note: the solid line in the white box in the figure is the horizontal position of the model body. The two east-west and north-south cross-sections pass through the place where the abnormal amplitude changes sharply).

表1 组合模型体空间位置坐标及物性参数Table 1 The space position coordinates and physical property parameters of combined model

切割线和测线的分布关系有以下3种情况：切割线在测线下方、切割线在测线下方及上方、切割线在测线上方，因此需要分别测试这3种情况的延拓精度.将模型体在观测面(图1a)引起的磁力异常(图1b)分别延拓至1000 m平面(向下延拓)、1300 m平面(部分向上延拓、部分向下延拓)以及1600 m平面(向上延拓)得到延拓磁力异常，以此模拟切割线分布在测线的下方及上方等情况，并分别对比图1的两条剖线处的理论磁异常和延拓磁异常(图2、图3、图4).

图2 理论磁异常与延拓磁异常对比(切割线位于测线下方)(a) 剖线1测线和切割线高程曲线； (b) 剖线2测线和切割线高程曲线； (c) 剖线1切割线理论磁异常与延拓磁异常曲线；(d) 剖线2切割线理论磁异常与延拓磁异常曲线.Fig.2 Comparison of theoretical magnetic anomaly and extended magnetic anomaly (the cutting line is below the measuring line)(a) Profile 1 measuring line and cutting line elevation curve; (b) Profile 2 measuring line and cutting line elevation curve; (c) Profile 1 theoretical magnetic anomalies and extended magnetic anomaly for the cutting line; (d) Profile 2 theoretical magnetic anomalies and extended magnetic anomaly for the cutting line.

图3 理论磁异常与延拓磁异常对比(切割线位于测线上方及下方)(a) 剖线1测线和切割线高程曲线； (b) 剖线2测线和切割线高程曲线； (c) 剖线1切割线理论磁异常与延拓磁异常曲线；(d) 剖线2切割线理论磁异常与延拓磁异常曲线.Fig.3 Comparison of theoretical magnetic anomaly and extended magnetic anomaly (the cutting line is below and above the measuring line)(a) Profile 1 measuring line and cutting line elevation curve; (b) Profile 2 measuring line and cutting line elevation curve; (c) Profile 1 theoretical magnetic anomalies and extended magnetic anomaly for the cutting line; (d) Profile 2 theoretical magnetic anomalies and extended magnetic anomaly for the cutting line.

图4 理论磁异常与延拓磁异常对比(切割线位于测线上方)(a) 剖线1测线和切割线高程曲线； (b) 剖线2测线和切割线高程曲线； (c) 剖线1切割线理论磁异常与延拓磁异常曲线；(d) 剖线2切割线理论磁异常与延拓磁异常曲线.Fig.4 Comparison of theoretical magnetic anomaly and extended magnetic anomaly (the cutting line is above the measuring line)(a) Profile 1 measuring line and cutting line elevation curve; (b) Profile 2 measuring line and cutting line elevation curve; (c) Profile 1 theoretical magnetic anomalies and extended magnetic anomaly for the cutting line; (d) Profile 2 theoretical magnetic anomalies and extended magnetic anomaly for the cutting line.

2.2 航磁成图数据精度评价结果对比

计算切割线分布在测线下方、上方等情况下的延拓磁力异常与理论磁力异常的偏差，不同类型的数据偏差结果如表2所示.

从表2可以得到，延拓后的磁力异常与理论值偏差小于±0.36 nT(延拓高差范围：0～473 m；异常幅值范围：1544.30 nT)，相对均方根偏差不大于±0.04%.

表2 不同情况下延拓值与理论值偏差统计Table 2 The deviation statistics of the continuation value and the theoretical value of different situations

延拓绝对误差和相对均方根误差均较小，可以忽略不计，因此可以满足后续利用此延拓方法进行精度评价的需求.

设计组合模型对比航磁规范法和曲面延拓法的成图精度评价结果.组合模型空间位置坐标及物性参数见表1.向磁力异常(图1b)加入均值为0，标准差为±2nT的高斯白噪声模拟实测数据(图5)，测线和切割线分布如图6所示.

图5 加噪磁力异常等值线图Fig.5 Contour map of magnetic anomaly with noise

图6 测线和切割线分布图Fig.6 Distribution of measuring lines and cutting lines

通过交叉点自动搜索技术得到切割线与测线的交叉点坐标及磁力异常值，将测线磁力异常通过非规则网延拓计算得到切割线交叉点处延拓磁力异常.统计交叉点处高程差、测线与切割线交叉点磁力异常差值以及切割线处延拓磁力异常与实测磁力异常差值如图7所示.从图7中可以看出，交叉点处测线与切割线磁力异常差值波动明显大于切割线交叉点处延拓磁力异常与实测磁力异常差值的波动，这是因为交叉点处测线与切割线磁力异常差值其中包含了数据精度误差以及测点坐标偏差所引起的误差值.如果不顾测点坐标偏差引起的磁力异常而直接计算数据精度，那将是不准确的.相反，延拓至切割线处的磁力异常与原测点在同一位置，不存在测点坐标(水平坐标以及垂直坐标)引起的偏差.

图7 测线与切割线交叉点处高程差值曲线以及两种方法的磁异常差值曲线(a) 高程差值曲线； (b) 两种方法在交叉点处磁异常差值曲线.Fig.7 Elevation difference curves at the intersection of measuring line and cutting line and magnetic anomaly difference curves of the two methods(a) Elevation difference curve; (b) Difference curves of magnetic anomalies at the intersection of the two methods.

为了详细对比航磁规范法和曲面延拓法的评价结果，设计以下两种方案：

(1)舍弃不同高程差的交叉点，分别舍弃高程差为10，20，30，40，50，60以及70m，再进行精度评价结果对比(图8a)；

(2)舍弃超过不同倍数预设精度的交叉点，分别模拟预设精度为0.5 nT和1.0 nT，舍弃倍数分别为1，2，3，4，5以及6倍，再进行精度评价结果对比(图8b)；

图8 航磁规范法和曲面延拓法精度评价结果对比(a) 舍弃不同高程差的交叉点； (b) 舍弃超过不同倍数预设精度的交叉点.Fig.8 Comparison of the accuracy evaluation results of aeromagnetic norm method and surface extension method(a) Discard the intersection of different elevations differences; (b) Discard intersections that exceed the preset accuracy of different multiples.

并统计在两种方案下精度评价结果的差异性(表3).

表3 两种方案精度评价结果对比Table 3 Comparison of the accuracy evaluation results of the two schemes

从图8和表3中可以得到：舍弃不同高程差的交叉点以及舍弃超过不同倍数预设精度的交叉点后再进行精度评价，航磁规范法的结果自身相差较大，但曲面延拓法的结果相差不大，基本在0.5 nT附近趋于稳定.这说明曲面延拓法适应性强，能较客观、准确地评价数据的真实质量情况.

根据航磁规范说明，利用测线磁异常和切割线磁异常进行精度评价，需要舍弃高差大于100 m以及均方根偏差大于3倍预设精度(本次设计预设精度为0.5 nT和1.0 nT)的交叉点，之后对比航磁规范法和曲面延拓法的评价结果如表4所示.

表4 不同方法精度评价结果对比Table 4 Comparison of accuracy evaluation results of different methods

从表4中可以看出，舍弃高差超100 m并舍弃异常差值超过不同预设精度的交叉点后，航磁规范法评价结果相差较大，曲面延拓法评价结果相差不大.但此数据为同一套数据，精度评价结果不应差别太大，因此航磁规范法存在一些不合理；相反，曲面延拓法精度评价结果趋于一致，这和实际情况是比较吻合的.

3 实际资料处理

选取内蒙古某地实测航磁数据进行对比航磁规范法和曲面延拓法的评价结果.研究区地形较为简单，中部为独立山包，东部和南部有连续的山脉，西部和北部较为平坦.在此研究区，图9为地形高程，研究区中部隆起，四周变化较缓.飞机延地形缓起伏飞行，测点高程如图10所示.图11为测点离地高度统计图，基本呈现正态分布，平均离地高度为69.68 m.测线沿东西向分布，测区中两条切割线沿南北向分布.测点点距约为2 m，线距约为100 m.实测磁异常如图12所示.

图9 地形高程等值线图Fig.9 Topographic elevation contour map

图11 观测点离地高度直方图Fig.11 Histogram of observation point height above ground

图12 磁异常等值线图Fig.12 Contour map of magnetic anomaly

将测线与切割线交叉点处高差、磁异常偏差及延拓磁异常与切割线实测磁异常偏差绘制图13.图13中两个红圈粉色区域，高差较大，测线和切割线磁异常偏差也较大，这部分偏差主要是由高差引起的，如果直接用此偏差计算数据精度，那么结果显然是不客观的.通过延拓的方式，即可消除因高差变化引起的磁异常偏差，使得精度评价时所使用的数据点位一致(水平位置和垂直位置)，最终的评价结果会更加客观、有效.

根据航磁规范法说明，对此区域航磁异常成图数据进行精度评价，具体评价结果对比如表5所示.

图13 测线与切割线交叉点处高差、磁异常偏差及延拓磁异常与实测磁异常偏差图Fig.13 Line chart of elevation deviation, magnetic anomaly deviation and continuation of the intersection of measuring line and cutting line

表5 航磁数据精度评价结果对比Table 5 Comparison of accuracy evaluation results of aeromagnetic data

从表5中可以得出，在本研究区，无论是舍点还是不舍点处理，曲面延拓法的精度评价结果都相差不大，数值趋于稳定并和航磁规范法(舍点)的评价结果基本一致.但这并不能保证其他地区实测航磁成图数据评价也是类似的结果.因此，本次进行了如下试验：若在切割线磁异常中加入部分高斯噪声后(图14)，继续对比两种方法的评价结果如表6所示.

图14 切割线测点磁异常与加噪磁异常对比图Fig.14 Comparison of magnetic anomalies of cutting line measurement points and noise-added magnetic anomalies

表6 航磁成图数据精度评价结果对比Table 6 Comparison of accuracy evaluation results of aeromagnetic mapping data

由于对切割线磁异常进行了加噪处理，原评价方法会继续将偏差大的交叉点进行舍弃，从而不参与评价，那么最终的评价结果和未加噪的数据评价结果就会一致，未能反映出数据中加入了噪声.而本文方法对于未加入噪声评价结果为2.01 nT，加入高斯噪声的数据评价结果为2.88 nT，前后两次评价结果明显不同，从而能较好地反映出原始数据中加入了噪声.

综上所述，曲面延拓法在提高航空磁测成图数据评价的客观性和真实性方面具有一定的贡献.

由于该方法使用等效源反演进行位场数据延拓，因此在实际应用中计算效率值得关注.该研究区航磁测线数据共有162953个，切割线数据共有5185个，经搜索得到129个交叉点.以交叉点坐标为中心，10倍平均点线距为半径布设等效源进行位场数据延拓，再配合OpenMP并行编程计算航磁成图数据精度，在计算机(参数：Intel(R) Core(TM) i7-8565U CPU，内存16G)上总耗时为745 s，内存需要约为963 MB.以上列出本文方法在该研究区的计算时间成本，对以后的实际工作中使用该方法起到一定的参考作用.

4 结论及建议

本文研究了基于位场曲面延拓的航磁成图数据精度评价方法，采用将航磁数据经各项改正后的测线磁测数据通过非规则网延拓算法延拓至切割线测点(或交叉点)处，再计算切割线测点(或交叉点)处延拓值与测量值的均方根偏差，以此作为全区的成图数据精度.曲面延拓法不需要舍弃任何交叉点而且克服了测点坐标不一致(主要是高程坐标)的问题，可以更加客观地评价航磁成图数据的精度，也作为评价航磁数据精度的一种新方法.通过理论模型测试和实际资料试验显示，该方法的评价结果更加稳定，对提高航磁成图数据精度评价的客观性和真实性方面有一定的贡献.

通过本文研究，曲面延拓成图数据精度评价方法也可推广至航空重力测量、海洋重磁测量或航空重磁分量或张量的成图数据精度评价中去.此外，也可采用本文提出的方法对航磁数据进行延拓处理得到某一光滑曲面或平面上的磁异常数据，以此作为成图数据，克服由于测点高程差别较大造成的成图问题.

致谢在本文撰写过程中，评审专家和编辑部等提出了很多宝贵的建议，在此谨表谢意！限于水平和篇幅，本文疏漏和不对之处，敬请批评指正.