水下导航一直以来都是导航领域的重要研究方向。随着丰富的海洋资源越来越受到各国的关注，海上冲突不断发生，我国临海海域的主权争议问题也越来越突出，这使得水下运载体导航技术的发展显得极为迫切[1]。目前，我国水下运载体导航定位最常采用的导航系统是惯性导航系统[2,3]。惯性导航系统是1种自主式无源导航定位系统，能同时给出位置、姿态等多种信息，且短时导航精度较高[4,5,6]。但其存在定位误差随时间积累的问题，当其应用于长航时，必须通过其他导航方法对其发散的误差进行修正[7,8,9]。重力匹配导航是根据地球重力场的变化来进行定位的，不需要发射和接收信号，不易受外界干扰，具有隐蔽性、实时性好、精度高等优点[10,11,12]，是最适合水下辅助惯性导航的导航方式之一。

水下重力导航系统通过水下运载体自身搭载的重力传感器，测量所处位置的重力特征信息，将实时测量值与预先储存的海洋重力场基准图进行对比，通过一定的匹配准则，对运载体的位置进行估计[13,14]，以修正惯导系统累积导航的误差，提高水下运载体导航定位信息精度。重力匹配导航系统的定位效果与水下运载体所处匹配区域的重力场特征密切相关。当所匹配区域重力场特征丰富，即重力特征值差异大时，各位置点易识别匹配，定位效果好；反之，当匹配区域重力特征不明显时，容易导致误匹配，定位效果差。因此，为提高系统的匹配效率，在进行重力匹配前，必须要结合水下运载体航行区域的重力场分布特性，扬长避短，选取出有利于重力匹配的区域[15,16,17]。本文对重力匹配导航适配区选取的关键技术及研究进展进行分析。

1、重力场基准图插值重构方法

重力适配区的选取，就是对已测量得到的重力场基准图进行适配区与非适配区的划分。目前，我国积累了较丰富的海洋重力信息数据，可以通过船载重力测量、航空重力测量及卫星测高反演这3种方式获取重力场基准图[18]，但实际测量中，由于时间及空间的约束，这3种方式获得的重力场基准图，均无法同时满足高分辨率和大面积的要求。而在重力匹配时，重力场基准图分辨率越低，重力传感器数据采集处位于重力场基准图网格点处的概率就越低。若此时单纯以格网平均重力值进行重力图的搜索，其引起的误匹配概率会随着分辨率的降低而提高[19]。因此，为提高匹配精度，在进行适配区选取前，应适当对重力场基准图进行插值重构，提高其分辨率。

哈尔滨工程大学、海军工程大学、解放军信息工程大学、海军海洋测绘研究所和中科院测地所等单位已开展了重力场基准图插值重构方法的研究。文献[20,21,22]提出将地质统计学中的克里金算法用于重力网格数据加密，并验证Kriging算法明显优于距离平方反比法。文献[18]研究了高精度重力异常数据库的生成技术，比较了克里格（Krige）方法、改进的2次曲面谢帕德方法、径向基函数方法3种不同格网插值方法，提出了利用光滑因子和边缘数据延拓，分别对MQS和RBF方法进行改进的算法。文献[23]提出了1种采用加权非线性拟合细分技术，提高数字重力场图分辨率的方法。文献[24]将计算机构图中的孔斯（Coons）曲面建模，引入到导航用海洋重力异常图的加密重构中，建立不同边界曲线的双1次Coons曲面和双3次C1Coons曲面重力异常模型。文献[25,26]分别从小波分析及重力场的模型/空间/时间尺度的角度，研究重力场的多尺度特征及其对匹配定位的影响。文献[27]基于Shepard插值模型的基本原理，从权函数的构造和采样点的选取2个方面，对Shepard函数逼近模型和算法实现进行了改进。文献[28]针对重力场不满足线性或一定程度多项式等信号行为的假设，利用平面傅里叶级数和径向基函数进行重力建模，克服了传统插值方法存在的问题。文献[29]提出了基于斐波那契数列寻优的2维高斯样条函数逼近局部重力异常场的方法。文献[30]提出1种基于双调和算子格林函数计算的格林样条插值法，其为用中心点位于各观测数据点的多个格林函数进行加权叠加而解析计算出插值曲面（曲线）的全局插值方法。文献[31]为解决重力场高频区插值误差，提出由小波分析根据已知基准图确定重点测量区域，由随机采样策略在重点区域增加测点以降低插值误差的方法。文献[32]针对反距离加权插值算法，仅利用海域重力观测值的空间信息，没有考虑重力观测值物理特性的缺陷，提出1种借助EGM2008重力场模型构建海域重力异常模型的反距离插值算法。目前研究较多的重力场插值的方法有Kriging空间预测算法、改进的Shepard算法、径向基函数插值算法以及格林样条函数插值算法。由各插值算法的插值实验结果可知，Kriging插值算法的球状变异函数、指数变异函数更适合于重力场插值[13]。径向基函数法的逆多重2次曲面、多重2次曲面以及薄板样条基函数的重力场插值精度较高，但因重力场的相关性，基函数的平滑因子对插值精度影响较大。格林样条函数插值算法亦可得到较好的重力场内插与外推效果。若在选择最优变异函数及平滑因子的情况下，径向基函数插值算法精度最优，格林样条函数插值算法略低，最次为改进的Shepard算法及克里格插值算法。

下面简要介绍几种主要的插值方法：假设待插值点位置为X0，计算区域内已知点位置为Xi，已知点个数为L,g为重力异常值。

1）距离倒数乘方法。距离倒数乘方法是最基础的重力基准图插值方法，其通过距离远近来对待插值点进行估计[33]，其计算公式为

式中：为未知点X0处进行插值得到的重力异常值；为第i个已知点Xi处的重力异常值；hi0为距离倒数；β为权重乘方；δ为光滑参数。在计算权重时，距离倒数加权插值法只考虑已知点与待插值点的距离远近。

2)Kriging方法。Kriging法是从地质统计学借鉴而来的1种最优内插算法。其用数据的相关性代替数据的欧氏距离，根据相关程度对样本点进行加权求和以估计未知点处的数据值，实质为利用构建的变异函数模型进行Kriging计算。待插值点的估计值为

式中iλ为已知点上的权系数。Kriging方法有2个假设条件

式中：γ为变异函数；ψ为与方差最小有关的拉格朗日常数。

3）改进的2次曲面Shepard方法。改进的2次曲面Shepard方法插值重构效果较好，其与距离倒数乘方法相似，限定了邻域范围并通过2次曲面拟合函数对待插值进行估计，其计算方法为

式中iG(x)的计算方法为

式中：ai1、ia2、ai3、ia4和ia5通过各已知点的最佳距离加权最小二乘逼近。

4）曲线和曲面拟合法。常用的曲线和曲面拟合方法有：双线性内插法、2次曲线法、3次曲线法、加权2次曲线法[34,35]和基于Coons曲面的重力异常模型插值重构的方法等[24]。

2、重力特征统计参数

在获取高精度高分辨率的重力场基准图后，需要对基准图的重力信息分布特征进行描述与评估。常用的传统评价重力区域的统计参量包括重力场标准差、重力场变异系数、重力场经纬度粗糙度、重力场经纬度相关系数、重力场经纬度坡度、重力场坡度、重力场信息熵、重力场偏态系数和重力场峰态系数等[10,36,37,38]。它们从不同的角度反映了区域内重力异常值的离散程度、光滑程度、相关程度、集中趋势等重力场固有属性，对研究重力场特征和匹配精度的相互关系具有重要意义。

1）重力场标准差。标准差是常用的统计分布测量参数，σ可以衡量数据序列的离散程度越大，计算公式为

式中：为均值；m、n为统计区域大小；g(i,j)为统计区域内数据；(i,j)为数据标号。

2）重力场经纬度粗糙度。粗糙度是衡量曲面表面光滑程度的度量值，粗糙度越大，重力场重力异常变化越大[39]。经度方向的绝对粗糙度rλ和纬度方向的绝对粗糙度的表达式分别为：

3）重力异常经纬度相关系数。重力场重力异常相关系数反映了相邻点的重力异常序列的线性相关程度[39]。当相关系数值较大时，相邻点相关程度较大，在识别定位时容易发生混淆；反之，相邻点相关程度较小，重力异常变化较大，在识别定位时比较准确，适合作为适配区。经度方向相关系数Rλ和纬度方向相关系数的表达式分别，即：

4）重力场坡度。重力场坡度S为重力场曲面上某1点的法线方向和垂直方向的夹角[39]为

式中：经度方向坡度为

纬度方向坡度为

5）偏态系数。偏态系数Cs是说明统计区域内数据分布不对称程度的统计参数，Cs的绝对值越大，表明数据序列的偏斜程度越大，sC>0表示数据序列分布右偏，sC<0表示数据序列分布左偏。计算公式为

除以上常见的重力统计特征参数外，文献[40]提出了改进标准差特征参数，以及引入了图像匹配中边缘提取方法的边缘密度作为特征衡量指标。文献[41]提出了特征丰度特征参数指标用以衡量重力场变化特征丰富程度，定义为粗糙度与标准差的比值。文献[42]引入了图像匹配中的相关度指标，以衡量实时值与重力场基准图之间的相互作用关系，定义为基准图上某重力p点的相关曲面CFδ(p)的相关峰的尖锐性；还提出了跟踪度指标，用以度量重力基准图中的特征显著性，反映了格点所在局部计算区域的纹理特征，以及在匹配跟踪时，局部区域的抗噪性和稳定性。文献[43]引入了编码失真量[44]以及莫里兹提出的，用以表征重力异常协方差函数局部特征的相关长度[45]。文献[46]引入重力梯度能量，反映重力梯度值总量的大小。文献[47]对全张量重力梯度5个独立分量信息的多个特征参数进行提取和统计，分析了特征参数对匹配性能的影响。文献[48]引入重力异常可导航值，作为描述重力异常背景场的特征，计算了西太平洋海域1'×1'重力异常可导航值。文献[49]引入平均重力差参数，来表示某一区域的重力变化水平。

但以上统计特征参数方法，大部分为借鉴地形导航性分析思想，没有关注重力适配区的匹配定位效果，存在方向性差异的问题，即不能衡量重力适配区域的方向适配性。大量实验表明，通过传统特征参数选出的适配区，并不都在各个方向上适合匹配，大多数适配区为在一定方向范围内适合匹配。如有的区域，标准差较大的但只在经度方向上或纬度方向上粗糙度较大，另1个方向则比较小，即其重力场特征仅仅在某个方向上较为明显[50]。在实际航行时，有很大概率以不适合匹配的方向通过此区域，匹配效果差，甚至不如部分非适配区，如此容易产生误导。而反之，有些适配区其可能为在某一方向范围极为适合匹配，但传统统计参数却无法评估与挖掘出此特征，会将其忽略，视为非适配区，浪费导航区域。为解决传统重力统计特征参数没有考虑到区域在方向上是否适合匹配的问题，文献[51]发现适配区域的中心线与图形骨架[52]的概念相吻合，引入快速欧氏距离场算法和简化方法提取得到适配区简化骨架。文献[53]基于分形理论，通过分析重力异常序列3维曲面图的频谱特征，提出了等方性系数，用以衡量重力适配区是否在各个方向都适配，即等方向适配性。文献[54]对坡度智能衡量经度与纬度方向做出改进，提出局部计算范围内的所有格网点坡度的平方和开方，为新的全方向坡度值。文献[55]提出另1种基于虚拟航向的特征计算方法，给出了传统的重力特征参数在0°～360°航向角下的计算公式。文献[56]对方向适配性进行一步研究，受方向估计方法启发并以降低误匹配率为目的，提出了1个衡量导航点适配性的特征矢量参数，即重力异常差合向量，其模值可表征某重力网格点与局部区域的重力异常值差异程度，其方向为该重力网格点局部重力异常值差异程度最大方向。这些考虑方向适配性的指标较传统方法有先进性，但暂时仍不完善，有待进一步深入研究。

3、重力适配区选取方法

在对重力导航区域进行适配性分析后，需要提取重力适配区以供重力匹配。由于当前衡量重力场适配性的指标多种多样，且不同指标反映了重力场不同方面的特性[57]；因此，如何综合运用重力场不同方面的特征信息，成为适配区选择中的关键问题。国内能查阅到的关于重力场适配区选择方法方面的研究资料，相比国外要丰富得多。目前，相关学者已经提出了各种适配区选取方法，极少部分为基于单一的特征参数的选取方法，大部分为基于多种参数进行融合的选取方法。通过利用多指标综合评价方法对多种重力特征参数进行融合，以得到可综合评价重力导航区域的特征参数。目前在各领域应用较为广泛的多指标综合评价方法有主成分分析法、层次分析法、综合指数法、熵值法、因子分析法、模糊综合评价法、专家评价法、功效系数法、经济分析法、多目标效用综合法、全概率评分法等多种方法[58]。重力适配区技术相关研究人员通过将各种融合方法引入，以选取合理有效的重力适配区。

最传统的多指标融合方法为阈值法，即通过求取各种重力特征参数与其对应阈值比较的交集，作为最终选取准则。如文献[39]经过等值线图相关分析后，选取重力场标准差和经纬度方向相关系数作为重力特征参数，并通过经验确定阈值，应用阈值法得到最终选取标准。文献[59]基于数理统计选取重力梯度标准差、能量、绝对粗超度作为指标，并通过阈值法进行适配区选取。但阈值法中，参数的阈值需要反复实验确定，区域不同则参数的阈值也不同，参数阈值对适配区的选择非常关键。当参数阈值条件选择宽松时，会选择出一些无效的区域；当参数阈值条件选择严格时，又会少选一些区域，使适配区面积过小。这样使得阈值法的选取结果很难恰好满足要求，往往存在或多或少的不足。

其后研究趋势便转移为更高级的各种主客观多指标融合方法。层次分析法作为1种定性与定量分析相结合，能注重不同特征参数所反映的重力场本质特征，以及各特征参数之间的逻辑关系的主观多准则决策方法，被重力适配区研究人员广泛应用。文献[41]利用层次分析法，将重力异常熵、重力场标准差、经度或纬度方向粗糙度、坡度、相关系数以及重力场特征丰度进行融合，得到综合特征参数进行适配区选取。文献[60]则通过层次分析法，将重力场标准差、经度或纬度方向粗糙度、经度方向坡度、相关系数进行组合。文献[58]基于反演重力图的多项统计特征及匹配仿真结果，使用层次分析法得到综合评价指标。文献[61]基于主成分分析方法，可对相关性的数据进行降维，以获得1组独立的综合指标，并可在某种程度上挖掘出数据中隐藏信息的优势，提出了主成分分析与层次分析法相组合的重力梯度匹配区选择准则，首先利用主成分分析法确定对适配性影响较大的关键特征参数指标，再利用层次分析法对这些指标进行融合。与此相似，文献[62]提出了1种新型主成分加权平均归一化法，融合了重力场标准差、坡度标准差、粗糙度、重力异常差异熵和分形维数等特征。除此之外，文献[63]提出基于基本的量纲和谐π定理建立量纲公式，并利用重力特征参数进行基于SPSS(统计产品与服务解决方案）软件的回归分析，确定量纲式指数，得到综合特征参数，并分析重力场特征曲线拐点与重力匹配区域的关系，根据拐点判断得到滤波估计误差下降的区域[64]。但其存在需要效果好的适配区作为前提，且存在算法公式适用性不足的问题。文献[65]提出利用模糊决策理论方法，求解出指标的隶属度矩阵和权阵，以得到评价区域的综合参数。其后又提出基于信息熵的求解经纬度粗糙度的熵和，再乘以标准差得到最终的重力场特征参数信息熵，作为适配区选取参数[66]。文献[67]基于熵值法可以突出重力场的局部差异性，并且根据各个区域的实际重力异常数据求得最优权重，反映重力场评价指标信息熵值的效用价值，避免主观认为的影响因素的特点，利用熵值法对各特征参数赋权值，得到综合特征参数。但这些适配区选取方法，由于是对传统的特征参数进行融合，因此必然存在忽略导航区域方向适配性的问题。

除以上单纯的对多种重力特征参数进行融合作为适配区选取的方法外，相关学者还借鉴了图像、机器学习等领域的适配区分割方法。文献[68,69]利用等值线方法，对重力场特征数据进行区域分割，提取栅格特征数据区域的矢量边缘,并进行聚类处理。最后，利用凸壳算法提取这些区域的凸壳范围线,获得重力适配区。文献[42,70]引入水平集理论，采用1种基于芒福德-沙阿模型的全局化水平集的图像分割方法，对重力场可匹配区域进行分割。文献[51]提出计算适配区骨架，并由简化骨架点及其距离值生成局部适配区域，具有一定参考价值，但所选区域面积较小，较分散，连续性不好。其后又提出将局部重力图转换为8位灰度图像，提取灰度直方图复杂性、索贝尔（Sobel）算子、灰度共生矩阵法等图像纹理特征，并通过使用投影寻踪模型，得到适配性综合评价结果[71]。文献[47]对全张量重力梯度5个独立分量信息的多个特征参数进行提取和统计，基于支持向量机的方法[72]，运用K-fold交叉验证法和格网寻优法，得到最优支持向量机模型参数。

4、结束语

重力场特征作为无源且导航精度高的自主导航方式，国外相关研究已经较为成熟且证明了其应用价值。因此，重力匹配导航的研究具有极大意义。而重力场适配区的选取，是重力匹配导航的基础，只有在重力特征丰富的区域进行重力匹配，重力导航才能发挥作用。本文详细介绍了重力适配区选取技术中的重力场基准图的插值重构、重力场特征参数的提取、重力场综合特征参数构建以及其他适配区选取方法等关键技术。

现重力场适配区的研究，大都借鉴于相似领域如高程导航、地形导航、地磁导航等，也有部分借鉴于图像匹配领域等。在将来的研究中，针对适配区重力特征参数，可以研究的方向有：(1)如何深入且准确评价重力场方向适配性，以为航迹规划与重力匹配提供基础，并须研究如何与航迹规划或重力匹配相结合进行应用；(2)现大部分适配区特征都为针对于重力异常值，而针对重力梯度值的研究较少；(3)当前重力图为2维网格图，适配区特征都为1维特征，将来可以进行更高维的特征的提取。

针对适配区选取方法，将来可以研究的方向有：(1)借鉴更多多指标融合方法，不拘泥于传统常见的几种方法；(2)借鉴除高程、地形、地磁等相似领域外的其他领域的特征提取与特征区域获取方法；(3)尝试更多的适配区智能划分方法；(4)更注重适配区选取后如何与航迹规划和重力匹配相结合，而不只是单纯地进行特征分析。

参考文献：

[1]杨峻巍.水下航行器导航及数据融合技术研究[D].哈尔滨:哈尔滨工程大学,2012.

[10]李姗姗.水下重力辅助惯性导航的理论与方法研究[D].郑州:信息工程大学,2010.

[13]姚剑奇.水下重力辅助导航定位方法研究[D].哈尔滨:哈尔滨工程大学,2015.

[18]吴太旗,黄谟涛,欧阳永忠,等.高精度海洋重力异常格网插值技术研究[J].测绘科学,2008,33(5):70-72.

[19]付梦印,刘飞,袁书明,等.水下惯性/重力匹配自主导航综述[J].水下无人系统学报,2017,25(2):35-47.

[20]张京娟,刘繁明,郝燕玲.重力场的克立格估值研究[J].哈尔滨工程大学学报,2003,24(3):249-252.

[21]成怡,刘繁明,郝燕玲.海洋重力图加密仿真实验研究[J].中国惯性技术学报,2007,15(2):206-208.

[22]李森,安智明.Kriging模型在区域重力数据插值中的应用[J].测绘通报,2013(10):66-69.

[23]施闻明,杨晓东.重力匹配导航中的重力图细分重构方法研究[J].导航,2006,42(2):83-88.

[24]李姗姗,吴晓平,赵东明.导航用海洋重力异常图的孔斯曲面重构方法[J].测绘学报,2010,39(5):508-515.

[25]闫利,吴华玲.重力场基准图的多尺度分析研究[J].测绘通报,2009(2):10-13.

[26]朱庄生,周朋.重力辅助惯性导航中的重力场多尺度特性研究[J].地球物理学进展,2011,26(5):1868-1873.

[27]徐遵义,姜玉祥,赵亮,等.改进的Shepard算法及其在重力异常插值中的应用[J].武汉大学学报(信息科学版),2010,35(4):102-105.

[29]童余德,边少锋,蒋东方.基于高斯样条函数的局部重力异常场解析重构[J].测绘学报,2012,41(5):128-134.

[30]许丰,贺日政,张贵宾.格林样条插值算法及其应用[J].地球物理学进展,2013,28(4):91-98.

[31]刘繁明,姚剑奇,荆心.由小波与随机采样提升重力场插值精度的策略[J].大地测量与地球动力学,2015,35(6):948-952.

[32]曾敏,汪振华,邓凯亮,等.借助EGM2008的海域重力异常反距离加权插值算法[J].海洋测绘,2018,38(2):5-7.

[33]韩雨蓉.水下导航重力匹配算法研究[D].北京:北京理工大学,2017.

[34]刘毅,张彩明.约束5点决定二次曲线的研究及其在参数插值中的应用[J].计算机研究与发展,2005,42(12):2161-2168.

[35]刘晓艳,邓重阳.非均匀三次B样条曲线插值的Jacobi-PIA算法[J].计算机辅助设计与图形学学报,2015,27(3):485-491.

[36]李姗姗,吴晓平,张传定,等.我国重力场新的统计特征参数的计算分析[J].地球物理学报,2010,53(5):1099-1108.

[38]朱宇炜.基于重力信息的航迹规划方法研究[D].北京:北京理工大学,2016.

[39]程力,张雅杰,蔡体菁.重力辅助导航匹配区域选择准则[J].中国惯性技术学报,2007,15(5):559-563.

[40]马旭东,闫利,乔永杰.重力异常辅助惯性导航中匹配区提取方法[J].测绘科学技术学报,2009,26(6):399-402.

[41]李艳,刘繁明.重力图的可导航性分析[J].应用科技,2010,37(12):60-63.

[42]吴华玲,徐效波,刘波.重力场可匹配区域选择研究[J].测绘科学,2014,39(5):84-88.

[43]王伟.重力异常滤波匹配导航算法及适配性问题研究[D].郑州:信息工程大学,2015.

[44]胡小平,吴美平.水下地磁导航技术[M].北京:国防工业出版社,2013:94-109.

[45]陆仲患.地球重力场理论与方法[M].北京:解放军出版社,1996:307-309.

[46]李开寒.重力梯度辅助导航匹配区选择准则研究[D].武汉:武汉科技大学,2015.

[48]王虎彪,王勇,许大欣,等.重力垂直梯度导航值特征分布与可导航性分析[J].测绘学报,2010,39(4):364-369.

[50]马越原,欧阳永忠,黄谟涛,等.基于重力场特征参数信息熵的适配区选择方法[J].中国惯性技术学报,2016,24(6):763-768.

[51]刘繁明,姚剑奇,李艳.一种骨架提取的海洋重力适配区域局部分析方法[J].武汉大学学报(信息科学版),2014,39(4):428-434.

[57]董翠军.水下重力场辅助惯性导航匹配算法与适配区选择研究[D].武汉:武汉大学,2017.

[58]张辰,马林,张桂敏,等.一种新的重力辅助导航海区适配性分析方法[EB/OL].[2019-10-28].

[60]蔡体菁,陈鑫巍.基于层次分析法的重力匹配区域选择准则[J].中国惯性技术学报,2013,21(1):93-96.

[62]李钊伟,郑伟,房静,等.基于新型主成分加权平均归一化法优选水下重力匹配导航适配区[J].地球物理学报,2019,62(9):3269-3278.

[63]夏冰,蔡体菁.基于SPSS的重力匹配区域选择算法[J].中国惯性技术学报,2010,18(1):81-84.

[64]夏冰,蔡体菁.重力辅助导航的拐点判据算法明[J].舰船电子工程,2009,29(6):80-82.

[65]马越原,欧阳永忠,曲政豪,等.基于模糊决策理论的重力辅助导航区域选择方法[J].海洋测绘,2016,36(6):7-10.

[66]马越原,欧阳永忠,黄谟涛,等.基于重力场特征参数信息熵的适配区选择方法[J].中国惯性技术学报,2016,24(6):763-768.

[68]马旭东,闫利,乔永杰.重力异常辅助惯性导航中匹配区提取方法[J].测绘科学技术学报,2009,26(6):399-402.

[72]张凯,赵建虎,王锲.基于支持向量机的水下地形匹配导航中适配区划分方法研究[J].大地测量与地球动力学,2013,33(6):73-78.

王博,周明龙.水下重力辅助导航适配区选取的研究进展[J].导航定位学报,2020,8(03):32-39.