高景一号(SuperView-1)卫星2016-12-28发射，轨道高度530 km，幅宽12 km，卫星全色分辨率高达0.5 m，多光谱分辨率为2 m，能够彰显细腻的地物细节，适用于高精度地图制作、变化监测和影像深度分析^[1]。受制于星上存储能力和数据传输能力，影像空间分辨率的提高必然导致影像幅宽减小，因而覆盖同一区域所需要的影像数也越来越多^[2]。在利用多幅数字正射影像图(digital orthophoto map，DOM)生产标准分幅产品时，相邻影像需要采用镶嵌线进行镶嵌，研究镶嵌线的自动生成方法对于提高影像产品的生产效率具有重要的理论意义和实际应用价值。

目前镶嵌线的自动生成主要有4大类方法，分别是基于重叠区影像差异的方法、基于同名点的方法、基于形态学的方法和基于辅助数据的方法^[3-7]。基于重叠区影像差异的方法是对重叠区域中的每个像元计算差异值，从而形成一个二维差异矩阵，然后采用一定的搜索策略在该矩阵上选择一条两侧影像差异最小的最优路径作为镶嵌线，然而简单的搜索算法无法得到全局最优镶嵌线，极易穿越建筑物等影像差异较大区域；基于同名点的方法主要从几何特征考虑，所生成的镶嵌线是由重叠区域中匹配得到的一系列同名点连接而成，该方法采用尺度不变特征变换(scale-invariant feature transform，SIFT)、加速稳健特征(speeded up robust features，SURF)等匹配算法寻找同名点时，有可能匹配到非同名点，需要进行错误点剔除，会对结果造成影响，另一方面影像匹配的点一般是一些角点，而这些角点很多是建筑物的拐点，很有可能导致生成的镶嵌线穿越地物；基于形态学算法能很好地避开房屋等明显突出地表的实体，得到不穿越地物的最优路径，但算法复杂，受限于影像质量，不适用于大范围数据处理的开展^[6-7]；基于辅助数据的方法是利用辅助数据对镶嵌线的走向加以约束，辅助数据往往具有一些关于地物的先验知识，便于更好地对镶嵌线进行优化，但该类方法仍然需要与其他方法结合^[8]。

本文提出一种利用机载激光雷达(light detecting and ranging，LiDAR)点云辅助的镶嵌线自动选择方法，该方法使用原始点云和滤波后的点云进行高程插值得到数字表面模型(digital surface model，DSM)和数字高程模型(digital elevation model，DEM)，继而得到地表物体顶面高差数据，并在构建代价矩阵的过程中考虑地物高差的影响和镶嵌线穿过区域的重叠影像相似度，最终在代价矩阵中通过动态规划得到从起点到终点的最小代价路径，即为最优镶嵌线。

1 LiDAR数据辅助的自动镶嵌方法

本文方法的流程图如图 1所示。对机载LiDAR点云进行规则格网插值，得到反映地表高程的DSM。将机载LiDAR点云进行滤波，得到地面点和非地面点，仅对地面点进行规则格网插值，可生成反映地形高程的DEM。将DSM与DEM相减，可以得到地物高差图，它反映了房屋、植被等地物顶部相对于地面的高差，能够更好地引导镶嵌线绕开房屋、植被等障碍地物。而镶嵌线的自动生成，则通过构造代价图，并在代价图中通过动态规划求解起点到终点的最小代价路径实现^[9-10]。本文的代价图既考虑了地物的高度信息，也考虑了重叠区域影像的灰度和梯度信息。

1.1 地物高差图生成

DEM或DSM是通过内插规则格网点高程得到。首先，对机载点云数据$ \mathbb{G}$或滤波后的地面点数据${\mathbb{{{\tilde G}}}} $，构造不规则三角网，通过三角形确定的斜平面来表达地形的起伏情况。为了保证内插结果的稳定性，避免三角形相交处的高程突变导致内插结果的跳变，本文在内插高程时，既考虑格网点落入的三角形的高程，也考虑落入三角形相邻三角形的高程为：

$ {H_P} = \omega '{L_P} + \left( {1 - \omega '} \right){N_P} $

(1)

式中，P是某个格网点；ω′是协调因子，决定了相邻三角形格网对高程内插的影响，应设为一个小于1的正数，本文实验中设为0.5；L_P值是点P落入以点P₁、点P₂、点P₃为顶点的△L内插的结果，其内插高程值为：

$ {L_P} = \frac{{{\omega _1} + {H_{{P_1}}} + {\omega _2} + {H_{{P_2}}} + {\omega _3} + {H_{{P_3}}}}}{{{\omega _1} + {\omega _2} + {\omega _3}}} $

(2)

式中，H_P1、H_P2和H_P3是三角形的角点高程；ω₁、ω₂、ω₃为点P₁、点P₂、点P₃插值格网点P所对应的权重，取值为3个角点到插值点P的平面距离倒数。式(1)中，N_P是△L的3个相邻△L⁽¹⁾、△L⁽²⁾和△L⁽³⁾对应的内插结果，其计算方式为：

$ {H_P} = \frac{{{\omega _L}^{\left( 1 \right)} + {L_P}^{\left( 1 \right)} + {\omega _L}^{\left( 2 \right)} + {L_P}^{\left( 2 \right)} + {\omega _L}^{\left( 3 \right)} + {L_P}^{\left( 2 \right)}}}{{{\omega _L}^{\left( 1 \right)} + {\omega _L}^{\left( 2 \right)} + {\omega _L}^{\left( 3 \right)}}} $

(3)

式中，格网点在领域三角形上内插高程的方式同式(2)；ω_L⁽¹⁾、ω_L⁽²⁾、ω_L⁽³⁾为△L⁽¹⁾、△L⁽²⁾和△L⁽³⁾在格网点P高程内插中所对应的权重，取值为对应的3个三角形中心到格网点P的平面距离倒数^[11]。

为了得到地物高差信息，需要先对点云$ \mathbb{G} $进行滤波，得到地面点的点云$\mathbb{{{\tilde G}}} $。本文在试验中采用了布料模拟滤波(cloth simulation filter, CSF)方法^[12]进行滤波，该方法使用物理过程模拟地形表面，为假设的“布料”添加重力和结点间作用力，自下而上覆盖在地形上，最终形变结果作为地形表面。通过式(1)对点云$ \mathbb{G} $和滤波后的地面点的点云${\mathbb{{{\tilde G}}}} $分别进行高程内插^[13]，得到DSM和DEM。设置地物高差格网D，定义格网位置上DSM和DEM的高程值分别为H_DSM和H_DEM，将两者进行相减，即可得到地物高差格网点上的高程值：

$ \mathit{\boldsymbol{D}} = {\mathit{\boldsymbol{H}}_{{\rm{DSM}}}} - {\mathit{\boldsymbol{H}}_{{\rm{DEM}}}} $

(4)

1.2 镶嵌代价矩阵构造

在设计代价矩阵时，综合考虑高差矩阵D中的地物高度信息、以及重叠区域影像的灰度和梯度信息，其构造如下：

$ {\rm{Cost}}\left( p \right) = \left( {1 + {\gamma _1}{D^*}\left( p \right)} \right) \times \left( {{C_{{\rm{chon}}}}\left( p \right) + {\gamma _2}G\left( p \right)} \right) $

(5)

式中，p对应两景正射影像重叠区中的某个像素位置；γ₁为高差代价系数，本文设置为10；γ₂为灰度梯度代价系数，设为1；D^*(p)为归一化高差，其值为：

$ {D^*}\left( p \right) = \frac{{D\left( p \right) - {D_{\min }}}}{{{D_{\max }} - {D_{\min }}}} $

(6)

式中，D_max和D_min分别为高差格网D中的最大高差值和最小高差值；D(p)为点p在高差格网D中内插得到的高差值；C_chon(p)为Chon代价函数，其值为：

$ {C_{{\rm{chon}}}}\left( p \right) = \left( {1.0 - NCC\left( p \right)} \right)/2.0 $

(7)

式中，NCC(p)为以点p为中心的5像素×5像素窗口中，两景重叠影像I₁和I₂灰度值的归一化相关系数；式(5)中G(p)为图像梯度差异代价函数，其值为重叠影像上的梯度差异的模(此处采用2范数)：

$ G\left( p \right) - {\left\| {{\rm{Grad}}\left( {{I_1}\left( p \right)} \right) - {\rm{Grad}}\left( {{I_2}\left( p \right)} \right)} \right\|_2} $

(8)

在代价矩阵的基础上，已知镶嵌线的起点S和终点E，以代价最小为原则，采用Dijkstra算法进行动态规划路径搜索，进而得到最优镶嵌线。

2 试验与分析

本文采用4组高景一号影像对进行实验，其详细信息如表 1所示。高景一号数据的几何精校正和正射影像生产采用多源遥感影像处理系统(multi-source image processing，MIP)完成。在MIP系统处理流程中，卫星影像的几何精校正通过基于有理函数模型的区域网平差实现，从而保证相邻卫星影像在重叠区具有较高的配准精度。区域网平差时使用高精度全球导航卫星系统(global navigation satellite system, GNSS)外业控制点，使得卫星影像有较高的绝对定位精度，从而保证与机载LiDAR数据的配准。本次实验采用的机载LiDAR数据基本覆盖全广东省陆域，获取时间为2018-2019年，点密度为2~5点/m²，激光点的测距精度达到cm级，平面精度优于1 m。

由于本文仅关心双像镶嵌线的优化问题，因此试验中仅将4组卫星影像分别进行双像镶嵌测试，通过本文方法在机载LiDAR数据的辅助下获得的自动镶嵌线结果如图 2所示。图 2中，红线表示无LiDAR辅助的镶嵌线，蓝线表示LiDAR点云辅助的镶嵌线。为了对机载LiDAR数据辅助的效果进行评价，对比方法仅采用Chon代价构造代价矩阵，而不适应LiDAR数据，并通过Dijkstra算法进行最小代价路径搜索得到镶嵌线。两种方法得到的镶嵌线穿过房屋的次数统计如表 2所示。

由表 2可知，相比于不使用辅助数据的方法，使用本文提出的机载LiDAR数据辅助方法自动生成的镶嵌线穿过房屋的次数显著减少，不过仍然无法完全避免镶嵌线穿过房屋^[14]。部分镶嵌线细节对比截图如图 3所示，其中, 图 3(a)表示传统方法生成的镶嵌线；图 3(b)表示机载LiDAR数据辅助方法生成的镶嵌线。

由图 3可知，本文方法通过在目标函数中对穿过地物高差较大地区的镶嵌线设置较大的代价，使得镶嵌线基本位于裸地、道路等无遮挡的地表，从而有效避免自动生成的镶嵌线两侧出现明显地物差异。但由于LiDAR点云滤波可能漏掉一些建筑物，同时，本文方法未考虑侧视卫星影像中高大建筑物遮挡的因素，因此，尚无法完全避免镶嵌线穿过房屋。

3 结束语

本文提出了一种激光点云辅助的高景一号卫星影像镶嵌线自动生成方法，基于地物高度信息及重叠区域影像灰度和梯度信息差异构建代价矩阵，并对初始镶嵌线进行动态规划优化，使其有效避免穿过房屋等存在高程突变的地物。该方法已经成功应用到自然资源广东省卫星应用技术中心的影像处理系统，有效提高了卫星影像的镶嵌效率及质量。不过，本方法生成的镶嵌线对建筑物阴影及投影差导致的倾斜影像不能完全有效避开，这将是下一步要重点研究解决的问题。