航空重力测 交叉点的非遍历逼近方法
屈进红1,2
, 姜作喜1,2
, 周锡华1,2,罗锋1,2,李芳1,2
1. 自然资源部航空地球物理与遥感地质重点实验室, 北京 100083;
2. 中国自然资源航空物探遥感中心, 北京 100083
基金项目：国家重点研发计划（2017YFC0601706；2017YFC0601705）
关键词：交叉点 快速逼近 航空重力 测 非遍历
引文格式：屈进红, 姜作喜, 周锡华, 等. 航空重力测 交叉点的非遍历逼近方法[J]. 测绘学 ，2022，51(1)：71-79. DOI:
10.11947/j.AGCS.2022.20210037
QU Jinhong, JIANG Zuoxi, ZHOU Xihua, et al. Non-ergodic approximation method for intersections of airborne gravity survey network[J]. Acta Geodaetica et Cartographica Sinica, 2022, 51(1): 71-79. DOI: 10.11947/j.AGCS.2022.20210037
阅读全文：
http://xb.sinomaps.com/article/2022/1001-1595/2022-1-71.htm
引 言
航空重力测量以飞机为载体，融合了重力传感器、惯性系统(INS)和差分全球定位系统(DGNSS)获取近地空中重力加速度的一种动态重力测量方法[1-6]。在中、高频重力场信息方面，弥补了地面重力、卫星测高在获取重力场信息方面的不足[7-8]，具有一定的优越性。由于其动态性特点，同一观测量无法进行重复观测，因此测线上往往布设切割线形成纵横交错的 状，根据交叉点重力值的符合程度来评估某个架次的测量质量[9]，直至评价整个测区的测量精度[10-15]。不仅如此，测线间存在的重力场水平差异通过与切割线上的交叉点重力差值进行水平调整，是数据处理中的关键性步骤[16]，对数据处理和地质解释具有重要的意义[17]。交叉点不符值的搜索除用于航空重力外，还可用于船载和其他机载的地球物理场测量中，甚至卫星测高中的平叉处理也涉及轨道交叉点的获取[18-19]，研究该方法具有十分重要的意义。
航空重力测 设计密集且数据采样率高，整体数据量庞大，如果采用逐一遍历或局部遍历搜索法，求取交叉点耗时长、效率低。目前常用的交叉点搜索方法主要是跳跃搜索法[20]、主副测线斜率法[21]、坐标平均法[22]、多项式拟合法[23]、滑动窗口法[24-25]及四叉树分解法[26]等。其中，跳跃搜索法、主副测线斜率法和坐标平均法，都是根据测 计算近似交叉点，再设置一定范围内进行遍历搜索出实际交叉点。多项式拟合法通过拟合测 飞行轨迹后计算近似交叉点，再设置搜索范围后遍历搜索，准确率比前几种方法有所提高，但搜索时间变长了。滑动窗口法对上述方法进行了改进，以一定搜索范围在测线上进行滑动遍历搜索，相交到切割线后确定搜索中心再以较小的搜索半径进行遍历搜索，准确率获得了保障，但没有提高搜索效率。四叉树分解法，通过矩形窗口在测线和切割线上滑动遍历锁定相交矩形范围后，再以四叉树结构分解重合区域，自适应加密剖分出遍历范围，是一次提升遍历搜索效率的尝试，但对于大型测 交叉点的搜索效率，仍无法突破瓶颈。
1 交叉点搜索方法
1.1 测 交叉点分析
航空重力测量布置测线和切割线相交叉的 状测 ，搜索交叉点就是获取它们相交于水平面位置的坐标、场值差或其他物理量差值。在航空重力测量作业中，飞机基本保持直线飞行，受风向、人为控制不当等原因对飞行产生扰动而发生偏离航迹称为偏航。因此，测 交叉点搜索方法设计时尽量考虑到极端偏航情况。
如图 1所示，假设航迹发生严重的偏航现象，A为偏航距。图 1模拟给出了测线line 1—line 5与切割线tie 1可能存在的5种位置情况：①line 1与tie 1正常相交；②line 2与tie 1真实相交，但没有被切割线两个端点相交到；③line 3与tie 1不相交，但端点延伸后被相交；④line 4与tie 1不相交；⑤line 5被tie 1的端点相交到，但实际没有相交。示意图中，实际相交的只有测线line 1、line 2，测线line 3、line 5也被误相交到，后续需剔除。搜索测 交叉点的核心就是判断测线与切割线是否真正相交，如果相交则获取二者相交的准确坐标。图 1中虚线外扩端点，就是为了解决因偏航可能被遗漏的实际交叉点，把所有测线和切割线端点延伸，延长量为规范中偏航阈值80m[27]的2~3倍，确保交叉点无遗漏。
图 1 测 中存在多种相交位置 Fig. 1 Multiple intersections at the boundary of the survey network
图选项
1.2 交叉点搜索思想
1.3 交叉点相交判断方法
图 2 两线段的向量叉乘 Fig. 2 Vector cross product of two line segments
图选项
图 2中两线段4个端点构成4个三角形：L1L2T1、L1L2T2、T1T2L1和T1T2L2，它们的叉乘面积按式(1)-式(4)计算
(1)
(2)
(3)
(4)
满足
符合交叉点存在条件后，计算交叉点位于切割线两端点之间的位置系数α为
(5)
通过式(5)可获得交叉点坐标P(x,y)，经度x、纬度y分别为
(6)
(7)
2 交叉点计算的程序实现
2.1 程序设计
按照上述搜索方法，设计了交叉点搜索流程，如图 3所示。装载测 后，对测线和切割线延伸的4个远端点组成线段后判断是否相交，存在相交后第1次获取近似交叉点坐标，再计算测线和切割线起始测点到近似交叉点的距离，根据各自的平均点距第1次快速逼近找出4个端点；若不相交接着把这4个端点再次进行线段相交，第2次获取近似交叉点坐标，重复上述过程第2次快速逼近找出4个端点；若不相交，程序进入微调逼近循环直至相交。
图 3 交叉点搜索流程 Fig. 3 Flow chart of intersections searching
图选项
图 3流程中，以微调逼近循环左侧为例，先以tie线段为基准，判断测线的相邻端点L1、L2是否与之相交，若不相交测线两个端点逐点向第2次近似交叉点坐标方向滑动，直至与tie线段相交；接着切割线的端点T1、T2，若不与line线段相交，两端点逐点向最新的近似交叉点坐标方向滑动，直至与line线段相交；因切割线端点的滑动，为防止发生测线端点与tie线段又不相交，再次重复上述过程，直至真正相交。而图 3中的微调逼近循环右侧，则用于在快速逼近后测线的两端点已经相交于tie线段，但切测线的两端点与line线段不相交的情况，实现过程与左侧类同。在交叉点相交后，再搜索下一条测线或切割线，直至整个测 的交叉点搜索完成。程序通过快速逼近和微调逼近的组合进行迭代式逼近，搜索方法思路清晰，程序比较容易实现。
2.2 实例分析
2.2.1 美国大地测量局EN01测
2013年3月美国大地测量局(National Geodetic Survey, NGS)对外发布航空重力测量EN01数据块Beta2版，此数据块国内学者研究较多[25-26, 28]，包含测点号、采样时间、经度、纬度和场值，共100830个数据点。表 1给出了EN01测 的设计参数，图 4为测 交叉点的分布图。
表 1 EN01数据块的设计参数 Tab. 1 Parameters of the data block EN01
| 参数名 | 参数值 |
| 测线间距/km | 10 |
| 切割线间距/km | 80 |
| 飞行高度/m | 6096 |
| 飞行速度/(km/h) | 407.44 |
| 采样率/Hz | 1 |
| 测点间距/m | 113 |
| 测线/条 | 26 |
| 切割线/条 | 5 |
| 交叉点数量 | 61 |
表选项
图 4 EN01数据块测 及交叉点 Fig. 4 Survey network and cross point diagram of data block EN01
图选项
表 2 EN01用不同搜索方法的统计信息 Tab. 2 Statistical information with different methods for EN01
| 搜索方法 | 搜索时间/s | 正确率/(%) |
| 逐一遍历法 | 41.9 | 100 |
| 优化逐一遍历法 | 2.1 | 100 |
| 坐标平均法 | 8.2 | 89.4 |
| 主副测线斜率法 | 7.8 | 89.5 |
| 滑动窗口法 | 7.8 | 100 |
| 四叉树分解法 | 0.28 | 100 |
| 0.00013 | 100 |
表选项
2.2.2 大型测 算例
表 3 与商业地球物理软件处理大型测 对比 Tab. 3 Compared with commercial geophysical software for large survey networks
| 搜索方法 | AS08 | TAZB | ||
| 搜索时间 /s |
正确率 /(%) |
搜索时间 /s |
正确率 /(%) |
|
| 优化逐一遍历法 | 5217 | 100 | 331.1 | 95.8 |
| Oasis Montaj8.5 | 8.1 | 100 | 83.1 | 100 |
| 0.031 | 100 | 0.061 | 100 | |
表选项
3 交叉点搜索效率分析
3.1 快速逼近设计
根据多个测 验证，常规测 可以设计成快速逼近循环直接逼近交叉点，但交叉点处在缺失数据的不连续测 中会使程序陷入死循环，另一方面，过多使用快速逼近也会影响到搜索效率，应视测 情况合理设置快速逼近的使用频率。
取两个国内航空重力测 为例，鞍山本溪直升机起伏飞行测 ANB的平均点距和标准差为：16.17±0.97m，对测 分别进行2~10倍加密和2~10倍抽稀，形成平均点距1.6~161.7m的共19个不同密度的测 。程序中分别设置快速逼近2、3和4次结合微调逼近后，统计这19个测 完成交叉点搜索时所用的逼近次数，以快速逼近2次的结果进行归一化，观察快速逼近3、4次的执行情况，如图 5所示。图中对平均点距在8.1m及以上测 ，设置3、4次的逼近次数归一化大于1，不如2次时的执行效率高。在平均点距8.1m以下的加密测 ，设置3、4次执行效率要好于2次，但也没有因为设置越高，执行效率被提高。
图 5 测 ANB快速逼近次数归一化 Fig. 5 Normalized graph of fast approximation times of ANB survey network
图选项
海南岛周边固定翼飞机测 HAIN的平均点距和标准差为：31.03±0.32m，同样进行2~10倍加密和2~10倍抽稀形成19个测 。然后进行快速逼近2、3、4次的逼近试验，以设置2次的结果进行归一化，发现设置2次时在各个测 密度中的执行效率都是最佳的，并且随着测 变疏它们之间的归一化结果趋于收敛，如图 6所示。HAIN测 测点距的标准差偏小，即测点间分布比较均匀，搜索结论与ANB测 出现截然相反的现象，测 加密下没有因为增加快速逼近的次数，而提高搜索效率反而出现下降。
图 6 测 HAIN快速逼近次数归一化 Fig. 6 Normalized graph of fast approximation times of HAIN survey network
图选项
对比不同情况的测 ，大量的试验都得出类似结果。程序中设置成快速逼近2次和微调逼近组合，在常规测 中就可以发挥出较好的搜索效率；如果测点距标准差偏大，同时测点距小于8m以下的测 ，宜采用快速逼近3次甚至4次，再结合微调逼近进行搜索。
3.2 EN01测 搜索分析
式中，M为所有测线的测点数和；nL为测线条数；N为所有切割线的测点数和；nT为切割线条数。
优化逐一遍历法：对切割线相交测线的近似交叉点坐标进行排序，切割线逐点遍历时只需遍历第1条测线，找到交叉点后把它排除在外，剩余测线重新排序，以此类推。避开了大量的重复遍历，缺点是对拼接测 、重复线混叠的时候不能保证100%搜索准确率。
表 4 EN01测 交叉点的搜索效率统计 Tab. 4 Search efficiency statistics for EN01 intersections
| 搜索方法 | 逼近次数 /次 |
平均逼近/ (次/个) |
搜索时间 /s |
| 逐一遍历 | 8.24×108 | 1.35×107 | 41.9 |
| 优化逐一遍历 | 4.65×107 | 7.62×105 | 2.1 |
| 400 | 6.56 | 0.00013 |
表选项
表 5 多个测 交叉点的搜索效率统计 Tab. 5 Search efficiency statistics for more network intersections
| 测 | 测 名 | 交叉点 /个 |
测线 /条 |
切割线 /条 |
测点数 /个 |
平均点距 /m |
逼近次数 /次 |
平均逼近 /(次/个) |
搜索时间 /s |
| 美国大地 测量局 航空重力 测 |
EN01 | 61 | 26 | 5 | 1.01×105 | 118.32±275.0 | 400 | 6.56 | 0.00013 |
| EN01′ | 56 | 24 | 5 | 9.51×104 | 117.43±16.64 | 266 | 4.75 | 0.0001 | |
| PN02 | 245 | 72 | 7 | 2.78×105 | 133.74±43.04 | 1167 | 4.76 | 0.00042 | |
| PN02′ | 238 | 70 | 7 | 2.71×105 | 133.73±7.72 | 1063 | 4.47 | 0.00039 | |
| AS08 | 310 | 48 | 13 | 6.14×106 | 6.14±43.5 | 199470 | 643.45 | 0.031 | |
| AS08′ | 231 | 45 | 9 | 5.43×106 | 6.06±0.59 | 1533 | 6.64 | 0.00055 | |
| 国内航空 重力测 |
ANB | 635 | 90 | 13 | 2.04×105 | 16.17±0.97 | 2971 | 4.68 | 0.0011 |
| HAIN | 2752 | 250 | 64 | 9.74×105 | 31.03±0.32 | 10732 | 3.90 | 0.0046 | |
| TAZB | 29517 | 809 | 170 | 5.47×106 | 30.60±0.93 | 136994 | 4.64 | 0.061 | |
| 船载测 | YYQC | 939 | 44 | 22 | 4.31×105 | 2.61±1.01 | 20492 | 21.82 | 0.0047 |
表选项
3.3 多个测 交叉点搜索分析
表 5中，美国大地测量局使用二轴阻尼稳定平台TAGS航空重力仪，因抗飞行颠簸性能差，人为切除了EN01、PN02和AS08测 中部分测线的扰动异常，从而造成测线中测点的不连续性。其中AS08测 切除了7条线的扰动异常，6个交叉点又位于剪切处，采样率为20Hz；剔除这7条线后为AS08′测 ，平均点距及标准差为6.06±0.59m。EN01和PN02测 的采样率均为1Hz，EN01′和PN02′为各剔除2条不连续测线后的测 。表 5、表 6中，除测 AS08、AS08′和YYQC使用快速逼近3次外，其余测 均使用2次，再采用微调逼近。
表 6 多个测 的逼近次数与搜索完成比例统计 Tab. 6 Statistics of the number of approximations and completion ratios formore survey networks
| 测 | 测 名 | 平均点距/m | 2次/(%) | 4次/(%) | 5次/(%) | 6次/(%) | 7次/(%) | 8次/(%) |
| 美国大地测 量局-航空 重力测 |
EN01 | 118.3±275.0 | 0.0 | 50.8 | 68.9 | 86.9 | 92.1 | 92.1 |
| EN01′ | 117.4±16.6 | 0.0 | 55.4 | 75.0 | 94.6 | 100.0 | – | |
| PN02 | 133.7±43.0 | 0.0 | 58.8 | 93.9 | 97.6 | 99.2 | 99.6 | |
| PN02′ | 133.7±7.7 | 0.0 | 59.7 | 95.0 | 98.7 | 100.0 | – | |
| AS08 | 6.14±43.5 | 0.0 | 1.1 | 不统计 | 54.0 | 79.9 | 87.8 | |
| AS08′ | 6.06±0.59 | 0.0 | 1.3 | 不统计 | 60.6 | 88.3 | 95.7 | |
| 国内航空 重力测 |
ANB | 16.17±0.97 | 0.0 | 51.2 | 85.7 | 96.5 | 98.7 | 100.0 |
| HAIN | 31.03±0.32 | 6.5 | 97.0 | 100.0 | – | – | – | |
| TAZB | 30.60±0.93 | 0.4 | 62.5 | 86.6 | 94.4 | 97.1 | 98.2 | |
| 船载测 | YYQC | 2.60±1.01 | 0.0 | 10.1 | 不统计 | 50.9 | 68.1 | 73.4 |
表选项
测 EN01、PN02，平均点距的离散程度偏大，每个交叉点的平均逼近次数搜索效率略低；AS08测 因有6个交叉点位于不连续处，搜索效率更低。而较规则测 EN01′、PN02′、ANB、HAIN、TAZB的平均逼近次数在3.9~4.75之间；平均点距小于8m的AS08′测 ，快速逼近调高至3次后的平均逼近次数为6.64，调高至4次时没有获得提升。
某海域船载重力测 YYQC的平均点距及标准差为2.61±1.01m，从标准差的偏离程度看出测点距比较离散或者测线之间的航行速度控制不稳定，测 的船测航迹保持也欠规整，如图 7所示。使用快速逼近3次后的每个交叉点平均逼近次数为21.82，调高至4次后，平均逼近次数下降至13.13，获得明显提升。
图 7 远洋船载重力测 Fig. 7 Ocean-going ship gravity survey network
图选项
表 6中，统计多个测 的交叉点逼近次数2~8次，分别对应交叉点搜索的完成情况。平均点距大于8m，较规则测 中只有大型拼接测 TAZB在逼近8次后完成98.2%，测 EN01′、PN02′、ANB和HAIN都在8次逼近内完成。表 5、表 6中得出认识：测 测点连续性、航迹规整、平均点距大小和它的离散程度适宜，对提高交叉点的搜索效率至关重要。
4 结论
(1) 针对常规航空重力测 ，交叉点搜索流程中快速逼近只需设置2次，就能实现快速搜索。
(2) 测 中测点不连续或者平均点距小于8m，且标准差偏大，建议快速逼近调升至3次，发挥出较好的搜索效率。
(3) 对于航行速度慢又不稳定，以及航迹不规整的船载重力测 ，为提高搜索速度，建议快速逼近设置成4次，弥补船载测 的不足。
初审：张艳玲
复审：宋启凡
终审：金 君
往期推荐
资讯
○ 直播预告 | 刘光明:《CGCS2000坐标转换》开讲
○ 山东科技大学测绘学科诚聘青年博士人才
○ 地理信息强国建设院士专家谈 | 秦昆，张永军，张祖勋 ：大力培养创新能力强的地理信息人才
○ 《北斗卫星导航标准体系（2.0版）》发布
会议
○ Geoinformatics 2022 暨CPGIS成立30周年年会（二号通知）
○ 会议通知 | 第三届IEEE地质、测绘与遥感国际学术会议
○ 关于召开大地测量与导航2022年综合学术年会的通知（第一轮）
○ 诚邀莅临丨“新测绘 新融合 新跨越”南方测绘2022用户大会
《测绘学 》
○ 测绘学 | 董杰：超导重力数据检测到的2011年日本东北大地震(Mw 9.0)震前重力异常及同震重力变化
○ 测绘学 | 徐焕：利用垂直重力梯度异常反演海底地形的解析方法
○ 测绘学 | 黄谟涛：重力异常向上延拓严密改化模型及向下延拓应用
○ 《测绘学 》2022年第2期目录
《测绘通 》
○《测绘通 》2021年第12期目录
○ 《测绘通 》2022年第1期目录
○ 地市级实景三维城市建设及应用
○ 图残差神经 络支持下的建筑物群组模式分类
《北京测绘》
○《北京测绘》2022年第1期摘要推荐
○《北京测绘》2021年第12期摘要推荐
○《北京测绘》2021年第11期摘要推荐
○《北京测绘》2021年第10期摘要推荐
《测绘科学技术学 》
○ 摘要 |《测绘科学技术学 》2021年第2期摘要推荐
○ 摘要 |《测绘科学技术学 》2021年第4期摘要推荐
○ 摘要 |《测绘科学技术学 》2021年第5期摘要推荐
○ 摘要 |《测绘科学技术学 》2021年第6期摘要推荐
《地球信息科学学 》
○ 《地球信息科学学 》2022年第2期佳文推介
○ 龚健雅院士：全球位置信息叠加协议与位置服务 技术研究进展与展望
○ 佳文推介 | 大数据城市通勤交通模型的构建与模拟应用
○ 专刊征稿：社会感知与地理大数据挖掘（征稿中）
《测绘工程》
○ 摘要 |《测绘工程》2021年第5期摘要推荐
○ 摘要 |《测绘工程》2021年第6期摘要推荐
○ 摘要 |《测绘工程》2022年第1期摘要推荐
○ 佳文推介 | 单目视觉技术在室内定位中的应用研究
《中国空间科学技术》
○ 火卫二地形地貌探测综述
○《中国空间科学技术》2022年第1期摘要
○《中国空间科学技术》2021年第6期摘要
○《中国空间科学技术》2021年第5期摘要推荐
《卫星应用》
○《卫星应用》2021年第12期摘要
○《卫星应用》2021年第11期摘要
○《卫星应用》2021年第10期摘要
○ 摘要 |《卫星应用》2021年第9期摘要推荐
《Journal of Geodesy and Geoinformation Science》
○《测绘学 （英文版）》专刊征稿 | 地图学与地球空间信息教育：理论与实践
○ 《测绘学 （英文版）》专刊征稿 | 用于三维地理信息的摄影测量和计算机视觉
○ Kexian WANG et al. | 《测绘学 （英文版）》（JGGS）精选论文
○ Li LIU et al. | 《测绘学 （英文版）》（JGGS）精选论文
《Satellite Navigation》
○ 徐元博士：面向有色测量噪声下UWB/INS组合行人导航的分布式卡尔曼滤波| SANA佳文速递
○ 杨飞博士：GNSS天顶对流层精化模型的构建与分析| SANA佳文速递
○ 姚铮教授：北斗PPP服务信号调制复用方案设计与质量评估| SANA佳文速递
○ 2022征文| SatNav“普适定位、室内导航与基于位置服务”专题
《自然资源遥感》
○ 《自然资源遥感》入驻“智绘科服”融媒体平台！
○ 《自然资源遥感》征稿：“海岸带空间资源及生态健康遥感监测”专栏
○ 摘要 |《自然资源遥感》2021年第3期摘要推荐
○ 摘要 |《自然资源遥感》2021年第4期摘要推荐
《Journal of Geovisualization and Spatial Analysis》
○《Journal of Geovisualization and Spatial Analysis》入驻“智绘科服”融媒体平台！
○ JGSA国际期刊2021年第5卷第2期论文摘要
○ 高被引论文推荐 | Journal of Geovisualization and Spatial Analysis
声明：本站部分文章内容及图片转载于互联 、内容不代表本站观点，如有内容涉及侵权，请您立即联系本站处理，非常感谢！