基于位场曲面延拓的航磁成图数据
张义蜜1，2，4，熊盛青3，于长春3，王万银1，2
长安大学重磁方法技术研究所，长安大学地质工程与测绘学院，长安大学西部矿产资源与地质工程教育部重点实验室 海洋油气勘探国家工程研究中心 中国自然资源航空物探遥感中心，自然资源部航空地球物理与遥感地质重点实验室 加拿大纽芬兰纪念大学
导读:
为了评价航磁成图数据精度，通常在完成精细调平等处理后依据航空磁测总精度评价方法进行评价。但经过精细调平或反复调平可以使成图数据精度达到很高，不能反应实际情况。针对这一问题，本文提出了一种基于位场曲面延拓的航磁成图数据精度评价方法，该方法是将航磁数据经各项改正后、精细调平前的测线数据通过非规则网延拓算法延拓至切割线测点（或交叉点）处，再计算二者的均方根偏差来进行精度评价。该精度评价方法不需要经过精细调平也可避免多次调平，不需要舍弃交叉点，并且克服了交叉点处点位坐标不一致的问题，可以客观地评价航磁成图数据的精度，也作为评价航磁数据精度的一种新方法，同时也可用该曲面延拓方法进行曲面数据处理成图。通过模型测试和实际资料处理，验证了该方法的正确性和实用性。
0 引言
航磁测量数据精度评价是非常必要的，但影响数据精度评价结果的因素有很多，Reeves指出利用切割线进行调平可能导致测线上的磁异常失真，若以调平后的数据进行精度评价可能会带入误差；现行的航空磁测技术规范中指出航空磁测总精度是由航空磁力仪系统测量和定位误差及各项校正不充分或不准确等误差组成，经过校正后进行评价得到航空磁测总精度，再经过精细调平等处理后评价得到成图总精度。分析得出，随着磁力仪灵敏度的提高，测量方法的改进，飞行高度是影响测量总精度的最主要原因，特别是测线和切割线交叉点处的高度偏差，即航磁测线和切割线测点的位置偏差及调平等处理会对评价结果造成一定的影响，在地形起伏大的山区影响更大。前人关于航磁成图数据精度评价方法可以归纳为航磁数据经过精细调平等处理后，计算切割线（重复线）与测线磁场差值的均方根偏差，以此作为成图总精度，具体方法有《航空磁测技术规范》中提出的方法（下简称航磁规范法）和重复线法。
航磁数据调平是航空物探数据预处理的关键性步骤，对数据处理和解释具有重要意义。目前调平方法可分为切割线调平、微调平、伪切割线调平以及其他方法。调平方法众多，调平后计算得到的磁场成图数据精度也不尽相同，这将影响航磁成图数据精度评价的客观性和准确性。
综合以上，前人关于航磁成图数据评价方法存在以下两点需改进的地方：（1）舍弃磁场差值大于3倍预设精度或高差大于100m的交叉点后，交叉点处高程不一致问题依然存在，精度评价过程中需解决交叉点处高程不一致问题；（2）因调平方法不同和调平次数不同导致最终的航磁成图精度不同，从而不能客观评价成图数据的质量情况。针对以上问题，本文基于位场理论，利用位场数据非规则网延拓方法，提出一种基于曲面延拓的航磁数据精度评价方法，改善了已有航磁成图数据精度评价方法存在的不足，从而可以更加客观、有效地评价航磁成图数据的精度。
1.1 航磁规范成图数据精度评价方法分析
航空磁测总精度是由航空磁力仪系统的测量和定位误差及各项校正不充分或不准确等误差的总和。在经过各项校正后，航空磁测总精度采用计算切割线与测线交点上磁场差值的均方根偏差σ来衡量，磁场水平精细调整前计算出的σ作为主要由定位引起的磁场测量误差。其计算公式为
在完成正常场、磁日变校正和整架次或每条测线粗略水平调整后，按式（1）计算调平前的总误差σ值。计算σ时，允许舍掉磁场梯度偏差大（600nT·km-1，交叉点处磁场变化>3σ者）、磁场偏差>3σ或高程差大于100m的交叉点，其余的交叉点均应参与σ的计算。在完成精细调平后应按上述舍点要求再计算σ，作为全区测量的成图总精度值。
以上为航磁规范成图数据精度评价方法，但若舍弃磁场偏差>3σ后，测线和切割线交叉点处磁场值为ai（i=1，2，…，n），切割线交叉点处磁场值为bi（i=1，2，…，n），满足|ai-bi|≤3σ（i=1，2，…，n），则
即仅舍掉磁场差值>3σ的交叉点后，磁场总精度已经小于2.1σ。
舍弃磁场差值>3σ的交叉点后，若ai-bi≤σ占比λ，σ<ai-bi<3σ占比1-λ，则有
当时，λ=7/8=87.5%。即，若ai-bi≤σ占比大于87.5%时，磁场总精度就会小于预设精度σ。从上面的分析得到，“舍点”后进行精度评价的做法足以影响最终结果的客观性。
若不“舍点”，测线和切割线的点位高差引起的磁异常偏差则会引入到精度评价结果中，从而影响最终结果的准确性。
基于位场的等效性，通过设立等效源代替真实场源，拟合观测数据建立等效源模型，进而利用此模型进行位场数据的处理和转换。本文利用此思想，由经典场论可知，在无源空间，位场延拓可以表示为Laplace方程边值问题：
由（4）式可得到半空间Dirichlet的解为
式（7）常使用矩形积分公式进行离散化。
式（6）采用矩阵形式可表示为：
刘芬等研究等效平面高程应大于等于0.4倍点线距dx（若无点线距，则为观测数据点距的平均数），若已知观测面高程zi（i=1，2，…，M）和待延拓面（点）高程zc，i（i=1，2，…，K），即可得到等效平面高程z00：
关于该方法计算效率问题，本文采用滑动窗口和OpenMP并行编程技术共同提高计算效率。在保证计算精度的条件下，刘芬等研究得到窗口半径一般选择20倍的延拓高度以上。而OpenMP是一个跨平台的多线程实现，主线程（顺序的执行指令）生成一系列的子线程，并将任务划分给这些子线程进行执行。这些子线程并行的运行，由运行时环境将线程分配给不同的处理器，因此可大大节省计算耗时。
2.1 非规则网延拓误差评价
本文设计了理论模型测试非规则网延拓引起的误差量级。组合模型体由五个直立六面体组成，空间位置和物性参数如表1所示。组合模型体在起伏观测面（图1a）引起的磁力异常如图1b所示。
表1 组合模型体空间位置坐标及物性参数
图1 起伏观测面及磁力异常图
（a）观测面高程等值线图；（b）模型体在观测面高程处引起的磁力异常等值线图（注：图中白框实线为模型体水平展布位置；两条东西向和南北向剖线穿过异常幅值变化剧烈处）
切割线和测线的分布关系有以下3种情况：切割线在测线下方、切割线在测线下方及上方、切割线在测线上方，因此需要分别测试这3种情况的延拓精度。将模型体在观测面（图1a）引起的磁力异常（图1b）分别延拓至1000m平面（向下延拓）、1300m平面（部分向上延拓、部分向下延拓）以及1600m平面（向上延拓）得到延拓磁力异常，以此模拟切割线分布在测线的下方及上方等情况，并分别对比图1的两条剖线处的理论磁异常和延拓磁异常（图2、图3、图4）。
图2 理论磁异常与延拓磁异常对比（切割线位于测线下方）
（a）剖线1测线和切割线高程曲线；（b）剖线2测线和切割线高程曲线；（c）剖线1切割线理论磁异常与延拓磁异常曲线；（d）剖线2切割线理论磁异常与延拓磁异常曲线。
图3 理论磁异常与延拓磁异常对比（切割线位于测线上方及下方）
（a）剖线1测线和切割线高程曲线；（b）剖线2测线和切割线高程曲线；（c）剖线1切割线理论磁异常与延拓磁异常曲线；（d）剖线2切割线理论磁异常与延拓磁异常曲线。
图4 理论磁异常与延拓磁异常对比（切割线位于测线上方）
（a）剖线1测线和切割线高程曲线；（b）剖线2测线和切割线高程曲线；（c）剖线1切割线理论磁异常与延拓磁异常曲线；（d）剖线2切割线理论磁异常与延拓磁异常曲线。
计算切割线分布在测线下方、上方等情况下的延拓磁力异常与理论磁力异常的偏差，不同类型的数据偏差结果如表2所示。
从表2可以得到，延拓后的磁力异常与理论值偏差小于±0.36nT（延拓高差范围：0~473m；异常幅值范围：1544.30nT），相对均方根偏差不大于±0.04%。
2.2 航磁成图数据精度评价结果对比
设计组合模型对比航磁规范法和曲面延拓法的成图精度评价结果。组合模型空间位置坐标及物性参数见表1。向磁力异常（图1b）加入均值为0，标准差为±2nT的高斯白噪声模拟实测数据（图5），测线和切割线分布如图6所示。
图5 加噪磁力异常等值线图
图6 测线和切割线分布图
通过交叉点自动搜索技术得到切割线与测线的交叉点坐标及磁力异常值，将测线磁力异常通过非规则网延拓计算得到切割线交叉点处延拓磁力异常。统计交叉点处高程差、测线与切割线交叉点磁力异常差值以及切割线处延拓磁力异常与实测磁力异常差值如图7所示。从图7中可以看出，交叉点处测线与切割线磁力异常差值波动明显大于切割线交叉点处延拓磁力异常与实测磁力异常差值的波动，这是因为交叉点处测线与切割线磁力异常差值其中包含了数据精度误差以及测点坐标偏差所引起的误差值。如果不顾测点坐标偏差引起的磁力异常而直接计算数据精度，那将是不准确的。相反，延拓至切割线处的磁力异常与原测点在同一位置，不存在测点坐标（水平坐标以及垂直坐标）引起的偏差。
图7 测线与切割线交叉点处高程差值曲线以及两种方法的磁异常差值曲线
（a）高程差值曲线；（b）两种方法在交叉点处磁异常差值曲线。
为了详细对比航磁规范法和曲面延拓法的评价结果，设计以下两种方案：
（1）舍弃不同高程差的交叉点，分别舍弃高程差为10，20，30，40，50，60以及70m，再进行精度评价结果对比（图8a）；
（2）舍弃超过不同倍数预设精度的交叉点，分别模拟预设精度为0.5nT和1.0nT，舍弃倍数分别为1，2，3，4，5以及6倍，再进行精度评价结果对比（图8b）；
并统计在两种方案下精度评价结果的差异性（表3）。
图8 航磁规范法和曲面延拓法精度评价结果对比
（a）舍弃不同高程差的交叉点；（b）舍弃超过不同倍数预设精度的交叉点。
表3 两种方案精度评价结果对比
从图8和表3中可以得到：舍弃不同高程差的交叉点以及舍弃超过不同倍数预设精度的交叉点后再进行精度评价，航磁规范法的结果自身相差较大，但曲面延拓法的结果相差不大，基本在0.5nT附近趋于稳定。这说明曲面延拓法适应性强，能较客观、准确地评价数据的真实质量情况。
根据航磁规范说明，利用测线磁异常和切割线磁异常进行精度评价，需要舍弃高差大于100m以及均方根偏差大于3倍预设精度（本次设计预设精度为0.5nT和1.0nT）的交叉点，之后对比航磁规范法和曲面延拓法的评价结果如表4所示。
表4 不同方法精度评价结果对比
选取内蒙古某地实测航磁数据进行对比航磁规范法和曲面延拓法的评价结果。研究区地形较为简单，中部为独立山包，东部和南部有连续的山脉，西部和北部较为平坦。在此研究区，图9为地形高程，研究区中部隆起，四周变化较缓。飞机延地形缓起伏飞行，测点高程如图10所示。图11为测点离地高度统计图，基本呈现正态分布，平均离地高度为69.68m。测线沿东西向分布，测区中两条切割线沿南北向分布。测点点距约为2m，线距约为100m。实测磁异常如图12所示。
图9 地形高程等值线图
图10 测点高程等值线图
图11 观测点离地高度直方图
图12 磁异常等值线图
将测线与切割线交叉点处高差、磁异常偏差及延拓磁异常与切割线实测磁异常偏差绘制图13。图13中两个红圈粉色区域，高差较大，测线和切割线磁异常偏差也较大，这部分偏差主要是由高差引起的，如果直接用此偏差计算数据精度，那么结果显然是不客观的。通过延拓的方式，即可消除因高差变化引起的磁异常偏差，使得精度评价时所使用的数据点位一致（水平位置和垂直位置），最终的评价结果会更加客观、有效。
图13 测线与切割线交叉点处高差、磁异常偏差及延拓磁异常与实测磁异常偏差图
根据航磁规范法说明，对此区域航磁异常成图数据进行精度评价，具体评价结果对比如表5所示。
表5 航磁数据精度评价结果对比
从表5中可以得出，在本研究区，无论是舍点还是不舍点处理，曲面延拓法的精度评价结果都相差不大，数值趋于稳定并和航磁规范法（舍点）的评价结果基本一致。但这并不能保证其他地区实测航磁成图数据评价也是类似的结果。因此，本次进行了如下试验：若在切割线磁异常中加入部分高斯噪声后（图14），继续对比两种方法的评价结果如表6所示。
图14 切割线测点磁异常与加噪磁异常对比图
表6 航磁成图数据精度评价结果对比
由于对切割线磁异常进行了加噪处理，原评价方法会继续将偏差大的交叉点进行舍弃，从而不参与评价，那么最终的评价结果和未加噪的数据评价结果就会一致，未能反映出数据中加入了噪声。而本文方法对于未加入噪声评价结果为2.01nT，加入高斯噪声的数据评价结果为2.88nT，前后两次评价结果明显不同，从而能较好地反映出原始数据中加入了噪声。
综上所述，曲面延拓法在提高航空磁测成图数据评价的客观性和真实性方面具有一定的贡献。
本文研究了基于位场曲面延拓的航磁成图数据精度评价方法，采用将航磁数据经各项改正后的测线磁测数据通过非规则网延拓算法延拓至切割线测点（或交叉点）处，再计算切割线测点（或交叉点）处延拓值与测量值的均方根偏差，以此作为全区的成图数据精度。曲面延拓法不需要舍弃任何交叉点而且克服了测点坐标不一致（主要是高程坐标）的问题，可以更加客观地评价航磁成图数据的精度，也作为评价航磁数据精度的一种新方法。通过理论模型测试和实际资料试验显示，该方法的评价结果更加稳定，对提高航磁成图数据精度评价的客观性和真实性方面有一定的贡献。
