CN111044041B - 基于重力场三维特征的重力辅助惯性导航适配区选取方法 - Google Patents

Abstract

本发明公开了一种基于重力场三维特征的重力辅助惯性导航适配区选取方法，能够提高适配区选取的有效性。该方法首先获得重力场采样数据，将重力场采样数据中的经纬度坐标转换为投影坐标x、y，将重力场采样数据中的重力异常值G进行因子变换，获得重力转化值z；x、y、z共同构成重力场三维数据；利用Delaunay三角剖分算法对重力场三维数据进行剖分，获得三角面构成的重力场三维模型；然后针对重力场三维模型的每个三角面，提取三角面矢量特征、三角面矢量变化特征和重力场复杂度特征，形成重力场多维特征；根据重力场多维特征选取重力辅助惯性导航适配区。

Description

基于重力场三维特征的重力辅助惯性导航适配区选取方法

技术领域

本发明属于惯性导航系统与重力导航系统交叉领域，具体涉及一种基于重力场三维特征的重力辅助惯性导航适配区选取方法。

背景技术

惯性导航系统是水下航行器广泛使用的导航技术，但其主要缺点是误差会随着时间推移而变大，不能满足水下运载体长航时的导航要求，必须通过其他导航方法来纠正，抑制位置误差的增长。重力场是一个稳定的物理场，具有良好的空间分布特征，近些年随着卫星测高和重力测量技术发展，提高了水下重力场的分辨率和精度。因此人们提出了水下重力辅助惯性导航方法，利用重力场匹配所得的导航定位信息，校正惯导误差，提升导航精度，进而实现水下运载体的长航时、高精度自主导航。

适配区的选取是水下重力匹配辅助惯性导航系统的关键技术之一，对航行区域适配能力分类及适配性的描述是影响定位精度的关键因素。水下重力匹配辅助惯性导航系统通常由海洋重力场背景图、海洋重力传感器、惯性导航系统、匹配算法四个模块组成，其中匹配算法是导航系统的核心技术，利用匹配算法的估计位置对惯性导航系统的导航信息进行反馈校正。在重力场变化不明显的区域，很难实现重力匹配，在重力场变化特征越是明显的区域，重力匹配效果越好，精度越高。传统的重力场适配性分析方法主要是是利用传统的数据统计特征参数或者统计特征参数的组合来表征局部区域的适配性，往往会导致统计结果不全面，不能全面、细致的描述重力场特征变化，容易产生适配区域漏选的情况，进而影响重力匹配效果。

传统适配区选取方法为多特征参数阈值法，仅对重力场进行简单的统计分析，特征参数包括标准偏差、经纬度相关性、重力异常变化系数、峰度系数和偏度系数等。并且，传统方法利用多重力场特征参数的阈值来选择适配区，往往需要经验判断，缺乏客观性。因此需要找到新的方法拟补目前方法的局限性，充分利用重力场的空间三维特征，为水下辅助惯性导航提供丰富的参考。

发明内容

有鉴于此，本发明提供了一种基于重力场三维特征的重力辅助惯性导航适配区选取方法，提高适配区选取的有效性。

为了解决上述技术问题，本发明是这样实现的：

一种基于重力场三维特征的重力辅助惯性导航适配区选取方法，其特征在于，包括：

步骤1、获得重力场采样数据，将重力场采样数据中的经纬度坐标转换为投影坐标x、y，将重力场采样数据中的重力异常值G进行因子变换，获得重力转化值z；x、y、z共同构成重力场三维数据；利用Delaunay三角剖分算法对重力场三维数据进行剖分，获得三角面构成的重力场三维模型；

步骤2、针对重力场三维模型的每个三角面，提取三角面矢量特征、三角面矢量变化特征和重力场复杂度特征，形成重力场多维特征；

步骤3、根据重力场多维特征选取重力辅助惯性导航适配区。

优选地，所述三角面矢量特征包括坡度和坡向；所述坡度为三角面的法线矢量与z轴的夹角；所述坡向为三角面的法线矢量在水平面xoy的投影与正北方向的夹角。

优选地，所述三角面矢量变化特征包括坡度变异系数和坡向变异系数；三角面A的所述坡度变异系数为以三角面A为中心的设定连续区域内所有三角面的坡度方差与坡度均值的比值；三角面A的坡向变异系数为以三角面A为中心的设定连续区域内所有三角面的坡向方差与坡向均值的比值。

优选地，以三角面A为中心的设定连续区域为：三角面A以及与三角面A相邻的三个三角面组成的区域。

优选地，所述重力场复杂度特征包括起伏度和粗糙度；所述起伏度为三角面内最大z值与最小z值之差；所述粗糙度为三角面面积S与其在水平面xoy上的投影面积S’之比。

优选地，所述步骤2形成重力场多维特征后，进一步利用主成分分析算法PCA对多维特征进行降维。

优选地，所述步骤3根据重力场多维特征选取重力辅助惯性导航适配区为：对重力场多维特征进行聚类，从而将重力场数据按照适配性分为K类，从K类中确定出惯性导航适配区。

聚类可以采用K-means算法或其他任何聚类方法。

优选地，所述步骤3进行K-means算法进行分类时，类别数K选取3，则分类后获得强适配区，弱适配区和非适配区。

优选地，所述从K类中确定出惯性导航适配区为：

根据地理位置和重力场多维特征所分的K个类别，绘制出K个分类区域；利用地理位置和重力异常值G绘制重力异常等值线；将重力异常等值线最密集处对应的分类区域确认为导航适配区；所述地理位置为经纬度或者为xy值。

优选地，所述步骤1对重力异常值G进行因子变换获得重力转化值z的公式为

有益效果：

与传统方法相比，本发明的主要贡献有以下几点：

1、创新性的将Delaunay三角剖分方法应用到了重力场三维模型的建立。利用重力场三维模型优势在于：一是对采样点分布没有具体明确的硬性要求，可以处理分析采样点分布不规则的重力场。二是三角形的形状和大小可以根据不同的情况灵活变化，所以基于不规则三角网的重力场三维模型可以高精度逼近复杂的重力场表面，精确的描述复杂的重力场起伏状态。

2、全面的适配性分析。该方法深度挖掘重力场的三维矢量特征，并用多维特征向量描述重力场，尤其是考虑了重力场的方向性，弥补了传统方法描述重力场特征单一性的不足。

3、该方法尝试利用了机器学习的方法，对重力场特征进行自动特征选取与聚类分析，减少了人工干预与人工经验判断，成功选取了一些新的优质适配区。

附图说明

图1为本发明适配区选取方法具体实施例流程图。

图2是重力场三维模型中基本单元——三角面结构图；

图3为本发明两个相邻三角面的示意图。

图4为一实例中传统特征参数阈值法适配区选取结果；

图5为一实例中本发明的选取结果。

具体实施方式

下面结合附图并举实施例，对本发明进行详细描述。

本发明提出了基于重力场三维特征的适配区选取方法，该方法基于Delaunay三角网理论构建重力场三维模型，全面挖掘重力场的三维矢量特征，构建重力场多维特征，根据重力场多维特征选取重力辅助惯性导航适配区。

本发明利用Delaunay三角剖分方法建立连续的重力场三维模型，三角形的形状和大小可以根据不同的情况灵活变化，所以基于不规则三角网的重力场三维模型可以高精度逼近复杂的重力场表面，精确的描述复杂的重力场起伏状态；从该重力场三维模型提取的矢量特征，能够更精确的表征重量场信息，尤其是考虑了重力场的方向性，弥补了传统方法描述重力场特征单一性的不足。

进一步地，本发明还尝试利用机器学习的方法，对重力场特征进行自动特征选取与聚类分析，减少了人工干预与人工经验判断，成功选取了一些新的优质适配区。

实验结果证明，与传统适配区选取方法相比较，基于重力场三维特征的适配区选取方法，能够挖掘出更多的重力场特征，能使得适配性分析更全面，能选取更多有效地适配区。

如图1所示，本发明基于重力场三维特征的重力辅助惯性导航适配区选取方法的具体流程如下：

步骤一、构建重力场三维模型。

重力场数据包括经纬度坐标(B,L)和重力异常值G。一般采用规则格网的形式获取采样点(B,L)的重力异常值G并存储，是一种典型的栅格数据结构，每个网格单元对应一个重力异常值。其缺点在于，一是不规则的重力场特性与规则的数据表示之间本身就不协调，二是规则的格网，是一个不连续的函数，虽然计算方便，但不利于复杂重力场的表示，不能准确表示重力场的结构特征和细部特征，不适用于重力场起伏复杂程度不同的地区。

为了解决上述问题，需要建立G值连续的重力场三维模型。因此考虑进行数字建模，将重力异常采样点按照某个规则连接成的面片的集合。在重力场三维模型建立过程中，需要解决两个问题：第一是统一单位，保证构建的三维模型是有意义的立方单位制；第二是构建的三维模型具有唯一性，保证构建的三维模型具有研究意义。下面的步骤12解决了统一单位的问题，步骤13解决了三维模型唯一性问题。

基于上述考虑，本步骤的重力场三维模型构建过程如下：

步骤11、获得重力场采样数据。

本发明不限制采集数据的路径，可以按照网格进行采集，或者按照任意路径进行采集。获得的重力场采样数据包括经纬度坐标(B,L)和重力异常值G。

步骤12、进行数据格式转换

本步骤需要将经纬度坐标(B,L)和重力异常值G统一到同一单位下。

针对经纬度坐标(B,L)：构建重力场三维模型首先需要将重力场采样点数据的经纬度坐标(B,L)转换为投影坐标系下的平面坐标(x,y)。本实施例采用墨卡托投影实现上述转换。

针对重力异常值G：重力场数据的投影坐标系(x,y)横纵坐标单位为米(m)，重力异常值G的单位为毫伽(mGal)，因为单位不同，所以需要一个转换因子Kz来构建重力异常值与高度之间的映射关系。

“重力异常空间校正”是用来描述重力采样点的位置变化对重力值的影响，由重力异常空间校正公式Δg＝kH＝0.3086H。可知，每1m的高程位移H产生约0.3086mGal的重力异常值变化Δg。可见K_z＝1/k，可以表达mGal到m的变换，可以应用到本发明将重力异常值G转换为以m为单位的重力转化值z。本发明将这种方法称为因子变换，本发明定义的因子为K_z＝1/k，k＝0.3068，变换公式如下：

z＝G·K_z

其中G为采样点重力异常值，z为因子变换后的重力转化值。

通过位置投影和重力异常值的因子变换，生成了oxyz坐标系下的重力场三维数据(x，y，z)，且该三维数据具有统一单位，这样生成的三维模型是有意义立方单位制。

步骤13、采用Delaunay三角剖分算法构建三维模型。

本发明利用Delaunay三角剖分算法对重力场三维数据(x，y，z)进行剖分，获得重力场三角网，即获得所述重力场三维模型。该重力场三维模型利用连续三角面来逼近重力场三维表面。在三角化算法中，Delaunay三角剖分算法的应用最为广泛，其具有很好的理论基础和数学特性。Delaunay三角剖分方法符合空外接圆准则和最大最小角准则，即将最近的三个节点组成三角形；任意四点不共圆；形成的三角形的最小角最大。如果将Delaunay三角剖分算法应用于重力场三角网的剖分，则得到的重力场三角网具有唯一性，从任何区域开始，生成的三角网最终都能达到一致，组成三角网的三角形互不相交、互不包含。基于Delaunay三角网的重力场三维模型GM可以表示为三角形T的集合。

采用现有的Lawson算法实现重力场数据的Delaunay三角剖分。图2所示为重力场三维模型的一个三角面，该三角面的基本元素有三个：节点、线和面。重力场采样点的空间位置和重力异常值决定了三角面的形状和大小，三角面是描述重力场三维表面的基本单元。图3示出了两个三角面形成的局部重力场变化情况，图中A₁B₁C₁和A₁D₁C₁为两个相邻三角面，其在xoy面上的投影为ABC和ADC，它们的法线矢量分别为n₁、n₂。基于这样的重力场空间矢量模型可以得到三角形的特性特征以及三角形之间的拓扑关系，进而可以对重力场进行三维空间分析，提取相应三维特征。

步骤2、针对重力场三维模型的每个三角面，提取三角面矢量特征、三角面矢量变化特征和重力场复杂度特征，形成重力场多维特征。

三角面矢量特征、三角面矢量变化特征和重力场复杂度特征分别从单个三角面、三角面变化以及重力场复杂度来描述重力场特征，这些特征能够很好地表征重力辅助惯性导航适配区信息。其中，三角面矢量特征从大小和方向两个维度对重力场微观尺度的特征进行了描述，即单个三角面的特征描述；三角面矢量变化特征则是对重力场中观尺度的特征进行了描述，即相邻三角面的变化特征描述；重力场复杂度则是对重力场宏观尺度的特征进行了描述，即重力场区域性的边和特征描述。而这些丰富的三维特征能够很好地表征重力辅助惯性导航适配区的适配性。

(1)本发明实施例中，三角面矢量特征包括坡度和坡向。

如图3所示，三维空间中的法向量在球面坐标系下需要两个互相独立的参量表征，一个是与水平面的倾斜角即俯仰角，另一个是水平面上的旋转角。法向量的偏移就体现在法向量在俯仰角及旋转角上的变化。三角形的坡度和坡向与三角形的法向量相对应。三角形平面方程可表示为：z＝a₀+a₁x+a₂y；a₀、a₁、a₂为曲线系数。

●坡度

坡度(slope)为三角面的法线矢量

与z轴即

的夹角，可以表示重力场表面某一位置的最陡下坡倾斜程度，是每个三角形中的最大重力异常变化率。坡度值范围为0°-90°。坡度值越小，重力场变化越平缓，坡度值越大，重力场变化越陡峭。

●坡向

坡向(aspect)为三角面的法线矢量

在水平面xoy的投影

与正北方向(x方向)的夹角，可以表示表面的水平方向，是坡度所面对的方向。坡向按照顺时针方向进行测量，角度范围介于0(正北)到360(仍是正北)之间，即完整的圆。不具有下坡方向的平坦区域将赋值为-1。坡向可以用于识别重力场表面某一位置处的最陡下坡方向。

aspect＝arctan(a₁/a₂)

(2)三角面矢量变化特征包括坡度变异系数和坡向变异系数。坡度变异系数与坡向变异系数反映了重力场坡度与坡向的变化趋势与走势。

三维重力场表面是三角剖面的集合面，每个三角面具有唯一坡度，对于每一个三角面来讲，在其坡面内，重力场没有发生变化，只有当从某一坡面转移到相邻坡面的时候，重力场才发生了变化。三角面的法向量反映了重力场曲面的倾斜程度。法向量区域一致区域，重力场曲面趋于平缓；法向量发生明显偏移的区域，重力场曲面呈现高低起伏。法向量趋于水平方向变化，重力场会呈陡峭的斜坡，法向量若趋于竖直方向变化，重力场趋于平坦，相邻三角形法向量偏移的趋势越明显，重力场的高低起伏变化越剧烈。因此三维重力场形态的变化部位主要集中在坡度与坡向发生变化的地方，通过计算坡度变异系数与坡向变异系数，可以得到坡度与坡向的变化情况。

●坡度变异系数

坡度变异系数(cv_slope)表征了当前三角面A与周围几个三角面之间的变化关系；又称坡度标准差率，反映重力场某一三角形面与其周围三角面坡度值之间的离散程度，变异系数越大，偏离程度越大，说明此处坡度变化越明显。

其中

slope_i为三角面i的坡度值。n为参与计算的三角面个数，参数计算的三角面区域为以当前三角面A为中心的一定连续区域内所有三角面的集合，包括三角面A。如果连续区域直选与三角面相邻的3个面，则n＝4，即连续区域是指三角面A以及与三角面A相邻的三个三角面组成的区域。

●坡向变异系数

坡向变异系数(cv_aspect)又称坡向标准差率，反映重力场某一三角形面与其相邻三角面坡向值之间的离散程度，变异系数越大，偏离程度越大，说明此处坡向变化越明显。

其中

aspact_i为三角形的坡向值。n为参与计算的三角面个数，参数计算的三角面区域为以当前三角面A为中心的一定连续区域内所有三角面的集合，包括三角面A。如果连续区域直选与三角面相邻的3个面，则n＝4，即连续区域是指三角面A以及与三角面A相邻的三个三角面组成的区域。

(3)重力场复杂度特征包括起伏度和粗糙度。重力场的复杂程度是必须考虑的因素，可以正确的反映重力场细部变化特征，重力场复杂度可以用起伏度，粗糙度来表示。

●起伏度

起伏度(relief)是指某个三角形内最大z值z_max与最小z值z_min之差，它是描述重力场特征的一个宏观性指标，起伏度反映的是一片区域内重力场起伏变化的特征，起伏度越大，则说明重力场特征信息越丰富。

relief＝z_max-z_min

●粗糙度

粗糙度为空间三角形面积S与其在水平面xoy上的投影面积S′之比，它是描述重力场特征的一个重要量化指标，粗糙度反映的是重力场表面的起伏变化和侵蚀程度。粗糙度越大，则说明重力场特征变化越丰富。

roughness＝S/S′＝1/cos(slope)

上述6个重力场三维特征因子构成重力场多维特征，对重力场的表面形态进行了更加合理、更加精确的描述。在实际应用中，还可以加入其它三角面特征参数，例如坡度方差、坡度标准差、坡向方差、坡度标准差、坡度变化率、破向变化率、斜边长度、边的变化情况等等，共同构成重力场多维特征。

步骤3、利用主成分分析算法(PCA)对多维特征进行降维和特征提取。

为了方便重力场特征分析与数据计算，首先确定多个特征因子之间的相关关系，本发明利用PCA方法对重力场多维特征向量T＝[x₁,x₂,…,x_m]进行降维处理和特征提取。本实施例中m＝6。将重力场m维特征映射到p维空间上(p≤m)，即投影到一个新的低维度的正交基空间中，重新构建出重力场的p维特征，用较少的新的P维特征向量来表示重力场，保留重力场特征向量中对方差贡献最大的特征，且新特征之间是彼此独立不相关的，达到对重力场特征的精简描述。

PCA算法的具体过程如下：

1)标准化处理

由于重力场的各个特征因子不在同一个数量级别，并且各个特征的单位不统一，因此首先要对综合特征向量T进行标准化处理。

2)主成分分析

首先建立T中变量的相关系数阵R，然后求R的特征值λ₁≥λ₂≥…≥λ_m>0及相应的单位特征向量，求出主成分，每一个主成分都是原始变量的线性组合，各主成分之间互不相关。

3)计算贡献率

计算主成分贡献率及累计贡献率，一般取累计贡献率达到95％的特征值λ₁,λ₂,…,λ_p，所对应得前p(p≤m)个主成分。

4)降维处理并特征提取

将重力场数据从原来的坐标系转换到新的坐标系下，第一个新坐标轴选择的是原始数据中方差最大的方向，第二个新坐标轴选取的是与第一个坐标轴正交且具有最大方差的方向，依次类推，可以取到p个坐标轴。将原来数据每个样本的特征值的个数降到了p个。最后获得具有最大特征值的特征向量构成的转换矩阵A_p，用重力场矩阵T乘以转换矩阵A_p，得到新的重力场特征向量T′。

步骤4、根据重力场多维特征选取重力辅助惯性导航适配区。

本实施例基于K-means算法对重力场多维特征进行聚类，从而将重力场数据按照适配性分为K类，从K类中确定出惯性导航适配区。

步骤3获得的重力场三维模型中包含M组特征，每组特征对应一个三角面，每组特征有p维元素。从M组特征随机选取K个特征向量作为初始的聚类中心，计算其余特征向量到这些聚类中心的相似性度量，并依据计算结果将其指派给相似度最高的那个聚类中心，更新聚类中心并重新分配特征向量，重复迭代更新聚类中心和特征向量分配环节，直到算法收敛，最终得到K个分组具有同一类别中特征向量之间相似度高、不同类别中的特征向量之间相似度小的特点。

通过对重力场特征进行聚类分析，最终将重力场数据按照适配性分为K类。需要结合一些先验知识确定K类中哪一类对应适配区。具体来说，存储一组数据(xy，经纬度B,L，重力异常值G，K分类值)，根据地理位置(B,L或x,y)和重力场多维特征所分的K个类别，绘制出K个分类区域，每一个分类区域用一个颜色或一种灰度表示。同时利用地理位置和重力异常值G绘制重力异常等值线。对比绘制出的K个分类区域和重力异常等值线，将重力异常等值线最密集处对应的分类区域确认为导航适配区。

本实施例的仿真实验中将重力场特征向量按适配性分为3类，包括强适配区、弱适配区、非适配性。那么，通过聚类获得3个分组，按照经纬度位置，将重力场背景图绘制成3种不同的颜色，用来表示区域的适配性。如图5所示，浅灰色区域的重力异常等值线密集，则该区域为适配区；大面积的深灰色区域重力异常等值线稀疏，为非适配区。而黑色线条区域为坡向变化明显的区域，虽然特征明显但区域连接性不强、区域面积不大，因此为弱适配区。与传统特征参数阈值法的结果(参见图4)相比，可以明显的看出，1)本方法选出的适配区覆盖了传统方法选出的适配区1-3，说明本方法的有效性。2)本方法选出的适配区范围更大，空间连续性更好，而传统方法都漏选了一部分重力场特征明显、适配性强的区域4-5，说明本方法的优越性。3)本文方法挖掘了重力场变化方向性特征，即坡向变化特征，对重力匹配导航具有重大意义，说明本方法具有独特性。

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种基于重力场三维特征的重力辅助惯性导航适配区选取方法，其特征在于，包括：

步骤1、获得重力场采样数据，将重力场采样数据中的经纬度坐标转换为投影坐标x、y，将重力场采样数据中的重力异常值G进行因子变换，获得重力转化值z；x、y、z共同构成重力场三维数据；利用Delaunay三角剖分算法对重力场三维数据进行剖分，获得三角面构成的重力场三维模型；

所述对重力异常值G进行因子变换获得重力转化值z的公式为

k＝0.3068；

步骤2、针对重力场三维模型的每个三角面，提取三角面矢量特征、三角面矢量变化特征和重力场复杂度特征，形成重力场多维特征；

所述三角面矢量特征包括坡度和坡向；所述坡度为三角面的法线矢量与z轴的夹角；所述坡向为三角面的法线矢量在水平面xoy的投影与正北方向的夹角；

所述三角面矢量变化特征包括坡度变异系数和坡向变异系数；三角面A的所述坡度变异系数为以三角面A为中心的设定连续区域内所有三角面的坡度方差与坡度均值的比值；三角面A的坡向变异系数为以三角面A为中心的设定连续区域内所有三角面的坡向方差与坡向均值的比值；

所述重力场复杂度特征包括起伏度和粗糙度；所述起伏度为三角面内最大z值与最小z值之差；所述粗糙度为三角面面积S与其在水平面xoy上的投影面积S’之比；

步骤3、根据重力场多维特征选取重力辅助惯性导航适配区。

2.如权利要求1所述的方法，其特征在于，以三角面A为中心的设定连续区域为：三角面A以及与三角面A相邻的三个三角面组成的区域。

3.如权利要求1所述的方法，其特征在于，所述步骤2形成重力场多维特征后，进一步利用主成分分析算法PCA对多维特征进行降维。

4.如权利要求1所述的方法，其特征在于，所述步骤3根据重力场多维特征选取重力辅助惯性导航适配区为：对重力场多维特征进行聚类，从而将重力场数据按照适配性分为K类，从K类中确定出惯性导航适配区。

5.如权利要求4所述的方法，其特征在于，所述步骤3利用K-means算法对重力场多维特征进行聚类，类别数K选取3，则分类后获得强适配区，弱适配区和非适配区。

6.如权利要求4或5所述的方法，其特征在于，所述从K类中确定出惯性导航适配区为：

根据地理位置和重力场多维特征所分的K个类别，绘制出K个分类区域；

利用地理位置和重力异常值G绘制重力异常等值线；

将重力异常等值线最密集处对应的分类区域确认为导航适配区；

所述地理位置为经纬度或者为xy值。

Images