二维地形匹配导航通过实时获取的二维地形数据与事先制备的带有地理信息的基准数据进行匹配以获取飞行器位置^[1-2]。为保证实时数据与机载计算机中的基准数据成功匹配并获得高精度的定位结果, 需要对候选地形区域的匹配性能进行分析和评价, 以甄选出飞行轨迹上匹配性能高的区域。该项技术被称为适配性分析^[3]。目前适配性分析研究的主流思想是对基准数据的特征进行提取和描述^[4-5], 通过经验判决、统计学习的方法建立特征量与适配性衡量指标的关系^[6-7]。现有方法在提取特征阶段皆进行了数据的遍历计算, 计算复杂度高; 构建的适配性决策准则均针对特定大小的实时数据与基准数据。而面向不同的飞行任务与应用场景时, 伴随着不同大小的实时图与基准图参数, 参数的多样性大大增加了适配性分析的计算量并难以制定统一的决策准则。因此, 面向二维地形匹配导航, 研究计算速度快、参数自适应的适配性分析方法具有重要意义。本文提出了一种基于显著性检测的地标点自动选取方法, 在秦岭山区、两广丘陵大范围的地形数据上进行了仿真实验分析。

1 基于显著性检测的地标点选取方法 1.1 基于地形分层统计的显著性检测

从统计学的角度考虑地形数据, 认为数据的稀缺性一定程度上代表了它的独特性, 其适配性能也相应好。占比低的数据以地势最低部分的山谷区域与地势最高部分的山峰区域为主, 并且呈现出一定的局部性^[8]。

设D(x, y)为二维DEM栅格数据, 设定θ_i为D在高程值域范围[min(D), max(D)]内的等间隔划分, 间隔为固定值θ_int。

(1)

对于每一个划分区间[θ_i, θ_i+1], 将数据D中θ_i < D(x, y)≤θ_i+1的像素赋值为1, 其余像素赋值为0, 生成对应的二值图像B_i, 如(2)式所示

(2)

基于每一幅二值图像B_i, 定义在该二值图像B_i上的稀缺性A(B_i)如(3)式所示

(3)

式中，‖B_i‖_p为二值图像B_i的p范数(p可取1或2, 文中实验p值为1)。B_i是二值图像, ‖B_i‖_p与图像中非零元素个数相关。公式(3)右端两项相加得到单幅二值图像B_i对应的稀缺值图像A(B_i)。最终基于分层统计稀缺性的显著图S_R(D)由N-1幅二值图像的A(B_i)累加而成。

(4)

1.2 基于抑制重复模式的显著性检测

频域检测常常用于数据去向关^[9], 此类方法不仅在检测效率方面具有很大的优势, 对重复模式的抑制也有良好效果。

设基准地形图为D, 对基准地形图像进行傅里叶变换转至频域

(5)

式中:F为图像傅里叶变换过程；Re, Im分别为频域的实部与虚部

(6)

(7)

A为幅度谱, P为相位谱。对幅度谱取对数, 随后进行高斯滤波, 得到新的幅度谱A_S, 如(8)式所示

(8)

式中:g为高斯滤波器; “*”为卷积操作

(9)

由于地形数据与自然图像都有如下特质:其幅度谱在低频部分的占比高, 高频部分的占比低。对幅度谱A_S的高斯滤波意味着平滑了幅度谱的低频部分, 同时增强幅度谱的高频。幅度谱A_S低频部分的数据在地形中分别代表了分布较为广泛的地形结构, 而高频部分代表了地形数据中较为稀缺的地形结构。故该操作有效抑制了地形中重复出现的地形结构。利用A_S与原有相位谱P, 进行傅里叶反变换, 得到显著图S

(10)

求取显著图S的幅度谱并进行平方, 如(11)式所示

(11)

随后对S_S进行归一化操作, 得到最终的抑制重复模式显著图S_N。

(12)

S_N中显著性目标的尺度与(9)式中高斯滤波函数g的尺度因子σ相关。当σ较小时, 在幅度谱上小范围的滤波, 抑制的是空域中较大尺度的重复模式, 同时凸显的也是空域中大尺度的显著目标(大范围结构、区域); 当σ较大时, 在幅度谱上较大范围的滤波, 抑制的是空域中较小尺度的重复模式, 同时凸显的也是空域中小尺度的显著目标(小范围结构、边界、个体像素)。针对地形匹配的显著性区域检测, 应当以地标点的大小为模板, 以实时图的大小为搜索范围, 检测出其中显著性好的区域作为地标点。由于地标点大小与g的滤波尺度因子σ成反相关, 而实时图大小应与g的滤波尺度因子σ成正相关。由此, 本文将实时图、地标点的大小引入(9)式中的高斯滤波函数, 若实时图D的大小为M×M, 地标点T的大小为m×m, (13)式中的σ为

(13)

式中，k为常值系数, 本文实验k值为1。由于实时图、地标点的大小是随着实际场景和应用条件变化的, 引入尺度因子σ的S_N也具有参数自适应性。

1.3 融合地形显著图

进行基于分层统计的显著性检测与基于去重复模式的显著性检测后, 得到相应的显著图S_R, S_N。其中, S_R为非归一化的显著图, 该显著图与DEM的高度范围有关, 范围越大分层越多, 累加得到的显著性值就越高, 可对不同地形的显著程度量化区别, 但无法检测地形中的重复结构; S_N为归一化的显著图, 可抑制地形数据中的重复模式, 并且对地标点和实时图大小的参数自适应, 但无法对不同地形的显著程度进行区分。公式(14)中将S_R, S_N相乘得到最终的融合显著图S_F。

(14)

融合显著图S_F中具有较高值的区域会被选为候选地标点, 候选地标点应在S_R, S_N中都具有较高响应, 因此结合了2种显著性检测的特性。通过非极大抑制方法选出所有的大于阈值T₁峰值位置, 将峰值的位置作为候选地标点的中心。

2 面向地形匹配导航的地标点自动选取

本文地标点选取方法以及仿真验证的整体流程如图 1所示。

1) 根据传感器参数得到实时二维地形数据的范围, 也即实时图的大小; 根据实时图的大小和搜索效率的需求, 确定所需要的地标点大小。

2) 在待测地形区域内遍历地截取实时图大小的地形数据为模板, 利用本文第一节提出的基于显著性检测的方法对待检测的地形进行显著图的计算, 并选取候选地标点。

3) 采用仿真验证的方法评估候选地标点的性能。采用相应的匹配算法和适配性评价指标对候选地标点在模拟实时图中进行匹配实验和性能评估。最终筛选出符合匹配性能要求的候选地标点作为最终确定的地标点。

二维地形匹配导航的匹配算法实质是模板匹配, 通常采用相关度量与距离度量方法, 归一化互相关算法NCC(normalized cross correlation)具有快速计算方案的优点成为最常用的匹配算法之一, 本文选取NCC相关匹配算法对地标点与模拟实时数据进行匹配。

定义E为地标点在模拟实时图上的单次匹配误差。

(15)

式中:(x, y)为通过NCC匹配算法的匹配位置；(x_T, y_T)为地标点实际位置, 定义匹配误差超过1个像素为错误匹配, 小于等于1个像素单位为正确匹配。

实验验证的评价指标采用匹配概率。匹配概率C_R为正确匹配次数N_C与所有匹配次数N_A的比率

(16)

3 仿真分析 3.1 数据及仿真平台

本文的DEM数据来自于地面分辨率约为30 m的ASTER-GDEM V2高程数据源。针对山地、丘陵两种典型的地形类型, 分别选取秦岭山地(33°N~34°N, 108°E~109°E)、两广丘陵(24°N~25°N, 107°E~108°E)区域内的DEM作为地形数据源。仿真平台采用i5-3230m 2.6 GHz CPU, 内存8 GB的笔记本电脑, 仿真环境为Matlab2017b。

3.2 仿真参数设置

设实时图大小为300×300, 地标点大小为50×50, 将2种大范围地形(3 600×3 600)从上至下, 从左至右依次裁剪成实时图大小(300×300)的144块地形区域, 采用本文基于显著性的地标点选取方法对2种地形区域中共288块地形区域进行计算, 并选取出相应的候选地标点。其中对分层统计显著图S_R的分层参数θ_int为30 m。本文中选取峰值阈值T₁为5×10^-3。实验验证部分对全部288块地形区域进行加噪处理, 在每块地形区域上叠加5种不同标准差的高斯噪声, 以模拟实际获取的实时地形数据。本文对2种地形分别叠加高斯噪声的标准差是平均地形标准差的10%~50%, 其中秦岭地区144块地形区域的平均地形标准差为244.11 m, 两广地区144块地形区域的平均地形标准差为133.37 m。统计分析地形区域内候选地标点匹配概率, 若在5种随机噪声情况下都能够正确匹配的候选地标点将成为最终“确定地标点”。本文引入“候选地标点正确率”作为地标点选取的评价指标, 该指标定义为经过验证的“确定地标点”与“候选地标点”的数量之比。

(17)

为了进一步对“确定地标点”与其他区域进行对比分析, 还选取显著图S_F中的显著图中较低位置的谷底区域作为“非地标点”, 对“非地标点”的匹配性能也做了同样的仿真验证, 由于“非地标点”的数量远大于“候选地标点”的数量, 为了保证统计数据上的公平, 进行统计计算时, 在显著图中位置较低的谷底区域采用随机选取的方法, 选取与该区域候选地标点个数相同的“非地标点”。

3.3 秦岭山区及两广丘陵的地标点自动选取

在秦岭地区144块地形区域中, 通过显著性方法选出的候选地标共673个, 经过仿真实验验证出其中497个地标点在5种高斯噪声情况下均可以正确匹配, 占候选地标点的73.85%;选取相同数量673个“非地标点”, 验证出其中有59个在5种高斯噪声情况下均可以正确匹配, 占非地标点的8.77%。在两广丘陵地区144块地形区域中, 通过显著性方法选出的候选地标共537个, 经过仿真实验验证出其中474个地标点在5种噪声情况下均可以正确匹配, 占候选地标点的88.27%, 选取相同数量537个“非地标点”, 验证出其中有41个在5种高斯噪声情况下均可以正确匹配, 占非地标点的7.64%, 相较秦岭地区的候选地标点与确定地标点数量都有所减少。

计算在秦岭山地的673个候选地标点与673个“非地标点”的平均匹配概率; 同样的, 求取两广地区的537个候选地标点与537个“非地标点”的平均匹配概率, 匹配性能指标如表 1所示。

在2种地形中, 通过本文方法选出的候选地标点的匹配概率都已到达80%以上, 而“非地标点”的匹配概率在40%以下。在2种地形上, 候选地标点的匹配性能指标均显著优于非地标点, 验证了基于显著性分析选出的候选地标点已具有较高的匹配性能。

采用文献[7]中的地形标准差σ_T、地形坡度均值E_S、地形坡度标准差σ_S 3种统计指标对“确定地标点”和“非地标点”的统计特性进行了计算。由表 2可见, “确定地标点”的3种统计在2种地形中都高于“非地标点”的统计指标值, 其中地形标准差σ_T、地形坡度均值E_S的差异较为显著。从表 2也可以看出本文选择出的地标点在地形起伏粗糙度的系列指标上具有可区分性。

在2种地形中生成的显著图像有如下规律:显著图S_R局部性较明显; S_N的显著区域较为分散独立; 根据融合显著图S_F选择出的地标点在秦岭山区都离散的分布于山谷和山顶区域, 而在两广丘陵区域趋向于分布在更为稀缺的山顶区域; 地标点的分布在秦岭地区的较两广丘陵的地标点更为密集。

对于尺度因子σ, 设计了另外2组不同地标点参数下所选取的候选地标点, 并进行了加噪验证实验(同4.3节)。3组地标点大小分别为50×50、100×100、150×150, 基准图大小均为300×300, 3组地标点适配性能的评价指标如表 3所示(指标R中列写了确定地标点/候选地标点的数量)。由表可见, 随着地标点尺寸的增大, 尺度因子σ减小, 所选地标点的数量也随之减少, 而从2种地形上看, 秦岭地区在相同尺度因子σ参数下的候选地标点与确定地标点都多于两广地区, 这与基于地形特征的地标点选取规律相吻合。候选地标点正确率均高于73%, 匹配概率均高于80%。

3.4 计算效率

本文对比了其他2种地标点的选择方法:第一种通过遍历匹配的方法, 对地形区域中的所有子区域进行匹配验证, 选择出在各种噪声情形下匹配性能好的子区域作为地标点, 该方法的实质就是将所有地形子区域进行实验验证; 第二种通过遍历分析所有地形子区域的地形统计特性, 并采用判定方法选出候选地标点, 随后通过实验验证确定出最终地标点。其中地形统计特征采用了文献[7]中的地形标准差σ、地形坡度均值E_S、地形坡度标准差σ_S。设(M, M)为实时图的尺寸, (N, N)为地标点的尺寸。在仿真验证部分的匹配算法采用NCC, 该方法复杂度为O(M²logM)。采用3种适配性分析及地标点选取方法的复杂度分析如表 4所示。

以秦岭地区一块区域300×300像素的区域作为测试对象, 设置地标点大小为50×50像素, 采用3种不同的方法选取地标点的用时如表 5所示, 其中假设特征判决方法所选的候选地标点数量与本文方法相同。遍历匹配的方法(I)最为耗时, 基于特征判定的方法(II)其次, 本文的方法用时最少且比其他2种方法提高了2至3个数量级的效率, 与算法复杂度分析结果相符合。在100 km²范围内的秦岭山地和两广丘陵, 本文自动选择地标点的方法均可在1 h内完成, 而遍历匹配的方法则需要数天甚至更长时间的计算。并且在实际应用中, 地标点的选取应紧密结合飞行器以及传感器参数, 才可确定出实时图的大小与所需地标点的大小, 而基于特征决策的地标点选取方法需要根据不同的应用场景重新制定不同参数的特征决策准则。本文的方法在执行效率和参数自适应方面都具有显著的优势。

4 结论

1) 在基于地形分层统计和频域去重复的显著图中可以有效选取出适配性高的地标点。显著图的生成对于实时图与地标点大小的参数自适应。

3) 相较传统地标点选取方法, 基于显著性检测的地标点选取方法可节省大量的预测与验证时间, 计算效率提升2至3个数量级。所选取的候选地标点在2种典型地形中不同参数下的正确率均高于73%, 平均匹配概率均高于80%。