武龙飞 卢会龙 刘金号 李玉飞
(中冀建勘集团有限公司, 河北 石家庄 050000)
影像镶嵌是不同时间或不同视角获取的多幅影像通过自动计算镶嵌线,将影像镶嵌获取覆盖整个测区的影像技术[1]。影像镶嵌能够将不同视角、不同时间拍摄的具有重叠区域的多景影像镶嵌在一起,突破了影像拍摄时间、空间的局限性,广泛应用于医学影像处理、遥感影像处理、虚拟现实分析、计算机视觉等领域[2]。影像镶嵌技术需要将具有重叠区域的影像转换到相同坐标系下再进行镶嵌。正射影像纠正能够将影像坐标转换的指定坐标系下,能够满足镶嵌处理的要求。但由于地势起伏的影响,纠正后的正射影像中存在一定未消除投影误差,导致镶嵌结果存在几何错位现象。因此,为提高正射影像镶嵌线质量和效率同时保证自动获取一条避免建筑物被分开成为一个热点问题[3]。
影像镶嵌技术由Milgram[4]首次提出,之后自动选取镶嵌线一直是众多学者研究的一个重要方向。Davis[5]在影像重叠区域内采用迪杰斯特拉(Dijkstra)算法搜索差值影像的最优镶嵌线,但Dijkstra算法采用穷举策略,计算量很大。王作勇等人[6]提出A*搜索算法的镶嵌线自动提取方法,该方法利用建筑物辅助数据避免了镶嵌线穿过建筑物。程多祥等人[7]采用蚁群算法选取最优镶嵌线,该算法依赖于蚂蚁数目,局部搜索能力较弱。Kerschner[8]提出基于Voronoi图的镶嵌网络生成与优化的方法,利用数字表面模型(digital surface model,DSM)数据对镶嵌网络进行优化调整,改方法能有效地避免网络定点落在树木、建筑物等存在高差的区域。袁修孝等人[9]通过半全局立体匹配算法获取立体影像的视差图,结合形态学滤波算法分离出非地面点与地面点,采用贪吃蛇搜索算法获取镶嵌线,该方法计算视差图运算量比较大而且对影像辐射信息要求高。宫思伟等人[10]和陈松等人[11]分别利用DSM和激光雷达(Light Detection And Ranging,LiDAR)点云数据先进行非地面区域提取,再搜索镶嵌线,其适用性受到辅助数据的获取与精度的限制。
针对上述存在一定的问题,有必要研究一种满足快速生产数字正射影像(digital orthophoto map,DOM)产品的镶嵌线智能检测方法。本文提出一种基于A*算法的数学形态学镶嵌线检测方法,根据重叠区域的差值影像搜索一条避开房屋等非地面区域的最短路径。
计算镶嵌线时,应尽量避免穿越高大地物如房屋、树冠以及色调差异较大的区域,同时避免生成镶嵌影像上存在地物接边不一致以及色调存在差异等问题,确保生成影像的质量满足应用需求。本文利用影像的灰度差值表示重叠区域的影像差异,计算公式为
(1)
式中,g1(x,y)、g2(x,y)分别为两幅待镶嵌影像在重叠区域中像素(x,y)处灰度值。
对于彩色影像,本文取各波段差值的最大值表示影像差异,公式为
(2)
由于正射影像是根据数字高程模型(digital elevation model,DEM)纠正得到,高出地面的地物在正射影像上存在投影误差,灰度差值会较大。因此差值影像可以较好地反映高大地物和色调差异较大的区域,使搜索到的镶嵌线更好地满足无缝镶嵌的要求。
由于遥感影像地物复杂性,导致差值影像上存在地物边界不清或地物边界不连续现象。在实际应用中通常在差值影像上检测镶嵌线,会存在镶嵌线穿过建筑物造成建筑物接边存在误差或建筑物表达不完整现象,难以获取理想的检测结果。本文对重叠区域的差值影像采用数学形态学方法进行改进处理,确保优化后的差值影像即能正确的表示影像差异又能够保证生成地物的完整性。
数学形态学的膨胀和腐蚀两个基本运算,主要用于影像上特征尺寸的增大和减少。对重叠区域的差值影像采用正方形结构元素处理时,膨胀和腐蚀运算定义如下:
对于重叠区域的差值影像ΔD的膨胀运算为取局部窗口内的差值影像最大值,具体公式为
(3)
式中,(xp,yp)为差值影像上以点(x,y)中心5×5窗口内任一点。
腐蚀运算为取局部窗口内的差值影像最大值,具体公式为
(4)
通过膨胀运算对重叠区域差值影像处理能够将临近断开的地物进行连接,填充整个地物内部的细小孔洞,确保镶嵌影像上地物的完整性,避免镶嵌线穿过存在高差的建筑物导致接边处存在误差。
1.3.1 A*算法原理
A*算法是一种静态路网中求最短路径的启发式搜索算法[12],其效果相当于找到最佳路径的Dijkstra算法,在简单情况下,其速度相当于广度优先遍历(breath first search,BFS)算法,因此改算法具有Dijkstra算法的准确性,也具有BFS算法的效率,Dijkstra算法以初始点为目标节点进行贪婪搜索会产生大量的盲目搜索,而BFS算法是以目标点为结点进行贪婪搜索在遇到障碍物时不能及早躲避,而A*算法是将两者结合以初始点和目的地之和最短进行搜索,选择最优路径,减少搜索范围、提高搜索效率。
A*算法的估价函数f(n)表示为
(5)
式中,f(n)表示为起点过节点n到目标点的估价函数;g(n)表示起始点到节点n的实际代价;h(n)是从节点n到目标点的估计代价。
对A*算法中h(n)的可纳性准则满足搜索到最小点估值,即要求
(6)
1.3.2 镶嵌线检测
镶嵌线检测可以简化为在重叠区域从起始点到目标点搜索一条影像差异最小、最短路径优化问题。陈松[11]等选取镶嵌线时结合DSM数据避免障碍区域,采用欧式距离计算g(n)和h(n)值,公式为
(7)
(8)
式中,(xsta,ysta)为起点的像点坐标;(xgo,ygo)为目标点的像点坐标;(xn,yn)为过起始点与目标点路径上任意中间节点n的坐标;D为影像格网分辨率。
通过式(5)、式(7)、式(8)组成A*算法镶嵌线检测估价函数,仅能沿水平与垂直方向进行移动,算法运行效率较低且难以满足影像无缝镶嵌的需求。为提高搜索效率,对四邻域搜索方向进行改进,添加对角线方向搜索采用八邻域方法确定邻接顶点。原有的式(7)、式(8)难以精确表示估价函数值。针对八邻域搜索特点,本文对g(n)与h(n)的计算方法进行重新定义为
式中,hd(n)表示节点n沿对角线方向到目标点的距离;hs(n)表示除对角线行走后还需走直线距离;w(n-1,n)表示节点n与节点n-1两节点间距离权重为
(11)
式中,T为灰度阈值。
为了减少搜索区域,本文将影像差值小于阈值T的区域内搜索镶嵌线,减少搜索区域的存储空间,提高算法运算效率。A*算法在搜索过程中通过建立OPEN表记录所有被考虑寻找最优路径节点,CLOSE表记录被舍弃的路径。在优化重叠区域的差值影像上如图1,设A为镶嵌线的起点,B为目标点,绿色格子表示自由区域,灰色格子为建筑物,采用A*算法搜索从A到B最优镶嵌线的步骤如下:
(1)设定初始灰度阈值T,对于无人机影像,一般阈值T取30即可;
(2)将起始点A存入OPEN表,由于查找OPEN表中只有A一个节点,故将A从OPEN表中移出,存入CLOSE表。然后将节点A相邻节点加入OPEN表内
(3)检查与当前节点相邻且不在CLOSE表中的节点,计算该点与相邻节点g(n)、h(n)、f(n)值,g(n)为当前节点n到起始点A的值,h(n)表示为当前节点n通过可通区域到目标点的值,f(n)为二者的和。
(4)选取OPEN表中f(n)值最小的节点,移除OPEN表同时设置为父节点并加入CLOSE表内。若父节点的相邻节点为障碍物时,不操作。否则查看是否在OPEN表内,如果不在OPEN表内,将该子节点加入OPEN表内。若子节点在OPEN表内,则计算g(n)值并与原开启列表的g(n)值比较大小。若当前路径的g(n)值大于原开启列表g(n)值,则按照原始列表关系,h(n)值不变。最后,将OPEN表中f(n)值最小的点设置为父节点,直到检测到目标点B则完成镶嵌线的检测;
(5)由目标点B反向追溯父节点到A,生成最优的镶嵌线路径(橙色格子)。
由于遥感影像数据量较大,若直接在重叠区域的差值影像上进行A*算法搜索,计算量和所需内存依然较大,本文结合影像匹配过程中金字塔匹配策略思路,采用金字塔策略进行由粗到细的镶嵌线检索[13]。影像金字塔是遥感影像多分辨率分析处理中常用重要结构,通过将原始影像逐步降低分辨率获取多次度的遥感影像数据集,从上到下影像分辨率逐渐变高影像细节也逐渐清晰,如图1所示。
图1 金字塔结构
金字塔分层搜索策略是一种由粗到精的搜索策略,首先,根据需要在分辨率较低的影像上进行初始搜索;然后,将搜索到的镶嵌线同比例变换到上以分辨率影像上;最后,在镶嵌线相邻区域以当前节点为父节点在一定区域内进行搜索细化,以获取更加精准的镶嵌线。初始搜索图层的选择是根据原始影像(第零层影像)大小来确定的,对于航空影像,一般可采用1∶3的分层策略,即在第三层影像进行初始搜索。采用金字塔分层策略搜索镶嵌线,既可以避免影像较大情况下,中间数据过多造成内存不足的现象,又可以提高搜索效率,满足快速镶嵌的要求。
为了验证本文算法的有效性,在Windows10 64位操作系统、i5-8400QM CPU、8 GB内存的台式机下,采用VS 2010环境下采用C++语言完成文中提出镶嵌线检测算法实验。实验影像为某城市地区一对数字正射影像,影像大小为5 548×3 863,重叠区域大小为5 352×2 269。
为了验证差值影像优化的有效性,分别以优化前和优化后的差值影像计算代价矩阵检测镶嵌线,如图2所示。图3(a)为两影像直接作差生成的差值影像的二值化结果,阈值T=30,图3(b)为图3(a)进行数学形态学膨胀处理后的结果,膨胀模板为5×5的正方形结构元素。图4(a)、图4(b)分别是以图3(a)、图3(b)为差值影像进行镶嵌线检测后的镶嵌结果。图5(a)、图5(b)分别是图4(a)、图4(b)的局部放大图。
(a)上影像
(b)下影像图2 试验影像
(a)优化前
(b)优化后图3 差值影像比较
(a)优化前
(b)优化后图4 镶嵌结果比较
(a)优化前
(b)优化后图5 局部镶嵌结果比较
图3(a)中建筑物的边界不连续而且车辆形状表示不准确,导致检测到的镶嵌线穿过了2个建筑物和2辆汽车,产生了几何错位现象,如图5(a)所示。优化后的图3(b)中地物表示准确,搜索到的镶嵌线能够很好地避开建筑物区域和车辆,如图5(b)所示,保证了镶嵌影像的质量,说明本文对差值影像的优化是有效的。
为了分析本文提出算法的效率和质量,分别利用基于四邻域和八邻域的A*算法在数据一上检测镶嵌线,并统计消耗时间和镶嵌线节点数目,如表1所示。
表1 镶嵌线效率与节点数量统计
由表1可知,金字塔分层搜索策略使四邻域搜索和八邻域的搜索效率均提高了30倍左右。同时,本文提出的基于八邻域的A*搜索算法速度比四邻域搜索算法快,而且在不采用金字塔分层策略搜索时,即搜索影像越大时,速度优势越明显。由于本文方法可以沿对角线进行搜索,搜索到的镶嵌线更短,镶嵌节点更少,可以减小后续镶嵌缝消除处理的计算量。可以看出,本文提出的八邻域搜索算法的效率和质量均优于四邻域搜索算法。
为了分析本文方法的有效性,分别利用本文方法和商业软件OrthoVista对正射影像进行镶嵌,如图6所示。白色线为本文方法检测的镶嵌线,青色线为OrthoVista检测的镶嵌线。表2为两种方法搜索镶嵌线的效率与质量统计数据。
图6 对比实验
表2 镶嵌线效率与质量统计
由表2可知,本文方法和OrthoVista软件的效率相当,对数据进行处理时,OrthoVista检测的镶嵌线穿过了3个建筑物,而本文方法检测的镶嵌线绕过了建筑物,主要落在道路和草地上。相比之下,本文方法检测的镶嵌线具有更高的质量。
本文提出一种基于A*算法的数学形态学镶嵌线检测方法,通过对差值影像进行数学形态学膨胀处理,能使镶嵌线更好地避开房屋、树冠等高出地面的地物,同时利用八邻域和金字塔分层搜索策略,提高了镶嵌线自动选取的效率。该方法能够较快地检测出一条较高质量的最短镶嵌线,镶嵌线比较简洁,能够减少后续镶嵌缝消除处理的计算量。