基于数学形态学改进的无人机影像配准算法

陆 可, 郑伯桢, 卢春盛, 王中元

(1.中国矿业大学 江苏省资源环境信息工程重点实验室,江苏 徐州 221116; 2.中国矿业大学 环境与测绘学院,江苏 徐州 221116;3.韶关市国土资源技术中心,广东 韶关 512026)

近年来,无人机飞行平台凭借其机动灵活、成本低廉、操作便捷等优势,在变化检测、大比例尺地形图测绘、环境监测等诸多领域得到广泛应用[1]。影像的配准是上述各类应用的基础与前提,配准的精度将直接影响各类应用的结果。然而无人机航拍时容易受到旋转、光照、飞行姿态、相机畸变等因素影响,且无人机影像空间分辨率高,利用基于特征的方法进行配准时,特征点对提取数量庞大,给影像配准带来了困难。因此,研究如何提高无人机影像配准的精度与效率具有重要意义。

图像配准方法按处理思路可以分为基于灰度的配准方法和基于特征的配准方法2类。基于灰度的配准方法可以视作模版匹配的过程,计算简便、容易实现,但是对不同影像上的灰度差异非常敏感,在实际应用中往往难以取得满意的配准结果;而基于特征的配准算法对光照和几何变化具有较强的鲁棒性,且配准效率较高,是图像配准领域的研究热点。

基于特征的配准可以分为基于线特征、面特征及点特征的方法,无人机影像配准中最常用的配准算法大都基于点特征。基于点特征的配准方法中最具代表性的是文献[2]提出的尺度不变特征变换 (scale invariant feature transform,SIFT) 算法,该方法对旋转、仿射变化、尺度缩放、光照变化具有一定程度的稳定性。为了提高SIFT的性能,研究者们提出了一系列基于SIFT的改进算法,如加速稳健性特征 (speeded-up robust features,SURF) 算法[3]、主成分分析尺度不变特征变换 (principal component analysis SIFT,PCA-SIFT) 算法[4]、基于仿射变换的尺度不变特征变换 (affine-SIFT,ASIFT) 算法[5]等,这类算法都是通过高斯模糊构造尺度空间,容易使得对象的边界信息被过度平滑,导致细节特征丢失。文献[6]通过非线性滤波构建尺度空间,提出了KAZE算法,克服了SIFT算法的这一缺点;文献[7]提出加速KAZE(accelerated KAZE,AKAZE)算法,用快速显示扩散(fast explicit diffusion,FED)求解偏微分方程,提高了运算效率,同时构造了改进的局部差异二进制 (modified local difference binary,MLDB) 特征描述符,增加了对旋转与尺度变化的鲁棒性。此外,有研究者通过多特征联合提取来满足特定配准场景。文献[8]利用SURF算法提取特征点,然后用二进制稳健基元独立特征 (binary robust independent elementary features,BRIEF) 算法进行特征描述,满足了无人机影像实时拼接的要求;文献[9]融合AKAZE算法与稳定的RootSIFT描述符进行图像配准,并融合了最近邻 (k-nearest neighbors,KNN) 算法、双向匹配及余弦相似度约束进行特征匹配,该算法更好地适应了无人机影像的匹配,但是增加了算法复杂度,难以满足实时性要求;文献[10]在SIFT基础上加入归一化互信息进行精配准,克服了影像尺度与光谱差异的问题。

虽然上述算法在提高匹配精度方面取得了一定的效果,但是在匹配效率上仍存在不足。另外,文献[10]发现树冠区域内特征点的误差对配准会产生较大影响。为解决上述问题,本文提出一种基于数学形态学改进的无人机影像配准算法。首先根据可见光波段差异植被指数(visible-band difference vegetation index,VDVI)以及数学形态学原理,对树冠区域内特征点进行剔除;然后将双向匹配算法与KNN算法相结合,对特征点进行粗匹配;最后利用渐进一致采样(progressive sample consensus,PROSAC)算法[11]进行精匹配,进一步提高配准精度。

1 基于点特征的影像配准原理

基于点特征的影像配准是最常用的配准方法,通过提取影像中鲁棒性强的特征点,建立参考图与待配准图之间的映射关系,然后计算几何变换模型,完成影像的配准。其主要步骤可以概括为特征点提取、特征描述及特征匹配3步。

1.1 特征点提取

特征点是图像灰度在x、y轴方向都有明显变化的像素点[12]。近年来,利用尺度空间性质检测特征点是一种主流的特征点提取方法。尺度空间是图像领域模拟人眼在不同距离辨别物体的方法,尺度大则观察范围大,物体模糊;尺度小则观察范围小,可以看到细节特征。因为计算机在分析场景时,预先并不知道图像中物体的尺度,所以需要将对象进行多尺度分解。高斯核函数是唯一的尺度不变核函数, SIFT算法在金字塔框架下构造高斯尺度空间[13],SURF算法采用盒子滤波近似高斯导数;但高斯模糊会导致边界信息丢失,并且同一尺度下的几组照片会平滑掉相同的细节和噪声,影响特征点的定位精度。AKAZE算法采用FED来构建一个各向异性扩散的非线性尺度空间,这样既保留了物体的边界信息,又提高了构建速度。尺度空间构建效果如图1所示。

在不同尺度的尺度空间效果图中,图像上的某个结构会在某个尺度上达到最大响应。通过Hessian矩阵[14]求取响应值是最常用的方法,其中最大响应值的点就是极值点,而极值点就是潜在的特征点。因此,通过计算滤波后图像的局部响应值并求极值,可以获得特征点的值,求取方式如图2所示。

图1 尺度空间效果图 图2 极值点选取

若该点在相邻的3个尺度中3×3×3邻域内是最大或最小值,则认为该点是极值点。通过确定极值点的亚像素精度,就可以估计特征点位置。响应值计算公式为:

(1)

其中:i为尺度等级;Li为图像;σi,norm为步长。

极值点亚像素精度计算公式为:

(2)

1.2 特征描述

在检测到图像上的特征点后,为了将这些特征点区分开,通常使用特征点邻域像素信息来对该点进行描述。描述符可以分为传统描述符与二进制描述符。传统描述符的代表是SIFT,是一种浮点型描述符,浮点型描述符鲁棒性强,但运算速度慢且占用过多内存;而二进制描述符在尽量不丢失必要信息的前提下,减少了内存占用率,并大大提高了后续特征匹配的速度。

其中,MLDB描述符是一种改进的局部差分二进制(local difference binary,LDB)描述符。LDB通过在特征点附近划分多种不同大小的窗口,分别求取网格灰度均值信息以及水平、垂直方向上的梯度信息,将灰度信息和梯度信息值小的区域赋值为0,大的区域赋值为1,形成的子集连接成LDB的描述符。MLDB在LDB基础上将LDB窗口旋转到特征点的主方向,使描述符具有旋转不变性,并对网格像素进行降采样,用采样后的平均值代替细分网格的均值,这样对尺度变化也有很好的鲁棒性。

1.3 特征匹配

为了将参考图与待配准图像的特征点进行对应,需要在特征描述符基础上通过一定的距离准则寻找相似度最高的特征点对。特征匹配分为粗匹配与精匹配2个步骤。

1.3.1 特征粗匹配

KNN的原理为:在一个特征空间中,如果样本中的k个最相邻的样本大多属于同一类别,那么该样本也属于该类别,通常将对象间的欧氏距离作为相似性度量。选取参考影像中的特征点A,在待配准影像中根据欧氏距离确定最邻近点B1和次邻近点B2,设A与B1之间距离为d1,A与B2的距离为d2,阈值为T。若d1/d2>T,则认为特征点A周围有大量相似点,不具有代表性,将其剔除;若d1/d2

传统的特征点粗匹配都是单向匹配,但是欧氏距离在各维向量中分布不一,容易影响匹配精度。因此,本文在完成单向匹配基础上,以待配准影像为基准,再求取B1在参考影像中对应的最邻近特征点A1,当A1=A时,保留该点对。双向匹配有效提高了特征点粗匹配精度。

1.3.2 特征精匹配

树冠区域SIFT特征点提取如图3所示。

图3 树冠区域SIFT特征点提取

在完成特征点粗匹配之后采用PROSAC算法对特征点进行特征点提纯,PROSAC是针对随机抽样一致性 (RANdom SAmple Consensus,RANSAC) 算法的改进。RANSAC算法从整个数据集中随机采样进行计算,得到模型参数,而PROSAC先将特征点对按照相似度进行降序排序,取相似度较强的点对进行迭代计算,得到模型参数,因此该算法提高了正确样本的抽取几率,同时也节省了运算时间。

在利用已有配准算法进行无人机影像配准时发现,无人机影像具有丰富的结构和纹理信息,提取的特征点数量非常庞大,并且位于树冠区域的特征点十分密集。对图3中2幅影像提取SIFT特征点,左右两幅影像分别提取12 452、12 571个特征点。但是经过RANSAC精匹配后只有22个点对匹配成功,且大多位于轮廓上。

树冠区域特征点匹配结果如图4所示,其中特征点的连线交错相连,可以直观地看出正确匹配点对很少。因此,如果能够预先剔除树冠区域内的特征点再进行特征描述与匹配,那么将大大提高匹配的速度与精度。

图4 树冠区域特征点匹配结果

2 改进的点特征配准算法

本文改进算法流程如图5所示。

图5 本文改进算法流程

卫星影像由于成像高度很高,可以近似看作正射影像；而无人机影像成像高度较低,当地形起伏较大时,特征点高差会呈现显著的变化,导致影像的配准结果不理想。树冠区域特征点与地面点高差大,且纹理信息丰富导致特征点密集,容易受外界因素影响,大都是无效的特征点。如果能提取出树冠范围,然后剔除该区域内特征点,那么将大幅提高配准的精度与速率。采用植被指数方法可以将树冠区域大致分离出来,但在无人机影像上树冠区域是按像素分布,区域破碎且不连贯。为了能提取出完整、准确的树冠区域,本文采用数学形态学(简称“形态学”)的先验知识来改进特征点剔除算法。在形态学中,图像被视为一个集合,结构元素是基础研究对象,通过交集和并集运算达到图像分割的目的。其中,腐蚀与膨胀操作是形态学运算的基础,将腐蚀与膨胀操作结合可以实现不同形式的运算。

2.1 数学形态学原理

在形态学中,通常使用结构元素逐像素扫描目标图像,并根据结构元素与图像关系确定对应结构元素中心位置的像素值,结构元素通常是较小的矩阵。腐蚀操作中,当结构元素全都在目标区域中,结构元素的原点对应像素点属于该区域;反之认为结构元素的原点对应像素点不属于该区域。与此类似,膨胀的基本思想是:当结构元素中任一点处于目标区域,认为结构元素原点对应的像素点属于该区域;反之认为不属于该区域。腐蚀和膨胀算法的运算公式分别为:

C1=A1ΘB1={(x,y)|(B1)(x,y)⊆A1}

(3)

C2=A2⨁B2={(x,y)|(B2)(x,y)∩A2≠∅}

(4)

其中:C1、C2分别为腐蚀、膨胀后的结果;A1、A2为目标对象;B1、B2为结构元;(x,y)为结构元的原点。腐蚀与膨胀效果如图6所示。图6中,红色是腐蚀掉的像素,蓝色是膨胀后的像素。腐蚀操作可以将图像的边界点消除,使图像沿着边界向内收缩,借此可以进行噪声去除、元素分割等;膨胀操作则会将图像边界向外扩张,连通2个对象。

图6 腐蚀与膨胀效果图

2.2 基于数学形态学改进的特征点剔除算法

腐蚀和膨胀是形态学的基础,将两者组合就可以实现形态学开、闭运算等不同形式的操作。开运算是先对图像腐蚀后膨胀,在不改变物体轮廓的情况下,平滑物体边缘,消除小噪声点。闭运算是先对图像进行膨胀再腐蚀,可以消除区域内的“空洞”,或者是去除物体上的小黑点。

受此启发,在特征点提取后,以VDVI为基础,通过形态学开、闭运算可以从无人机影像中获取完整、简明的树冠范围,然后将该区域内特征点剔除。算法具体步骤如下:

(1) 将无人机影像R、G、B 3个通道信息进行拆分。

(2) 计算图像的植被指数IVDV。因为无人机影像没有近红外波段,所以通过VDVI[15]进行计算,计算公式为:

(5)

其中:IVDV取值范围为[-1,1];ρred、ρgreen、ρblue分别为红、绿、蓝3个波段的反射率或像元值。

(3) 利用大津阈值法进行分割,将无人机影像分为树冠区域(前景)与非树冠区域(背景)。

(4) 确定合适的结构元素,将小于结构元的对象视为噪声点并移除。

(5) 利用形态学算法去除植被区域边缘的“毛刺”,填补区域中的“空洞”。

(6) 获取完整的树冠范围,剔除树冠内的特征点。

特征点剔除的具体流程如图7所示。

图7 特征点剔除流程

图7a、图7e所示分别为特征点剔除前、后的结果。从图7e可以看出,树冠区域内的特征点已全部被剔除掉。

在无效特征点剔除后,对剩余特征点进行描述与特征匹配。其中,特征粗匹配本文采用双向匹配方法,特征精匹配用PROSAC算法。最后,通过计算两幅影像单应矩阵将影像配准到同一区域,单应矩阵H计算公式为:

(6)

其中:(x1,y1)为参考影像上特征点坐标;(x,y)为待配准影像上特征点坐标;hmn为单应矩阵的自由度参数。

3 实验结果与分析

3.1 实验数据与设备配置

为验证提出算法的有效性,选取4组无人机影像进行实验,实验数据如图8所示。

图8 实验影像对

本文研究的无人机影像配准主要针对城市区域,图8所示为4组典型的城市场景影像对，其中左图为参考图,右图为待配准图,分辨率均为6 000×4 000像素。为便于后续变化检测研究,选取的影像是重叠度较高的同一区域2期影像,间隔时间2周。本文所有算法均在Python3.7环境下编译完成,部分开源代码调用OpenCV库实现并在此基础上进行改进。实验设备为64位Win10系统的笔记本电脑,其CPU为i5-8265U,内存为8 GiB。

图8a的影像代表典型的城市住宅区,建筑物密集,左侧为低层,右侧为高层建筑,小区中道路两侧有绿化植被;图8b的影像是城市高架公路,道路两侧以及环形匝道中有成片高大乔木,树冠面积较大;图8c的影像中有多种场景元素,包括建筑物、主干道路、河湖、小型公园,这可以验证算法在复杂场景下的表现;图8d的影像是在建区域,几乎没有树冠区域,通过此区域可以验证本文算法在没有植被覆盖区域的表现。

3.2 定性分析

为验证改进算法效果,首先对配准前、后影像进行定性分析。以SIFT算法为例,结合本文提出的形态学算法进行两期影像配准,其中特征点匹配情况如图9所示,配准结果如图10所示。

图9 特征点匹配情况

从图9可以看出:住宅区由于建筑物比较密集,且2周内几乎没有变化,特征点分布多,在剔除树冠区域内特征点后,特征点匹配效果很好;主干道路区域较空旷,特征点分布较少且位于树冠区域内的已被剔除,但是道路上有同款汽车导致特征点存在误匹配情况;混合区特征点数量与分布适中,特征点匹配效果好;在建区域特征点数量与分布适中,但变化较明显,会存在部分特征点误匹配情况。

图10 影像配准结果

图10中,为了更直观地查看配准效果,在ArcGIS中利用卷帘工具进行对比,其中配准后的影像置于图层顶层,参考影像置于底层并用灰度图像展示。通过对建筑物、道路、河流等的边界位置进行比对可以得出结论,配准后的影像与原影像基本重合,配准效果较好。

3.3 定量分析

为了对本文提出的算法性能进行定量分析，评价其性能提升的幅度，分别对4个区域数据集中的影像对采用SIFT、SURF、BRISK[16]、AKAZE 4种算法进行配准。实验分为2组，其中一组采用原算法匹配，作为对照组，另一组采用本文提出的形态学算法进行改进。具体的评价指标为特征点的正确匹配率(correct matching ratio，CMR)RCM、均方根误差(root mean square errors,RMSE)ERMS及运算时间t。RCM、ERMS计算公式分别为：

RCM=n/N

(7)

(8)

其中:n为正确匹配点对数量;N为所有匹配点对数量;ERMS取值范围为(0,1),值越小精度越高;(xi,yi)为参考图像中特征点的坐标;P(xi′,yi′)为待配准影像中特征点转换后的坐标。

4种算法改进前、后对4个区域数据集进行配准,评价结果见表1所列。由表1可知,从CMR看,改进后算法在4个区域数据集中均有不同程度的提升,SIFT算法特征点正确匹配率提高了约6.47%,SURF算法提高了约7.09%,BRISK算法提高了约4.56%,AKAZE算法提高了约0.69%。其中,AKAZE算法在提取特征点时,采用FED来构建一个各向异性扩散的非线性尺度空间,提取的特征点鲁棒性较强且数量适中,通过观察AKAZE特征点分布情况,树冠区域分布较少,其算法本身性能很高,因此本文提出的改进算法效果不明显。但是对于SIFT、SURF这种特征点分布密集且鲁棒性不是很高的特征匹配算法而言,采用本文提出的改进算法后,其特征点匹配正确率有较大提升。

表1 4种算法改进前、后3个评价指标结果对比

为验证改进算法对配准精度的提升,在参考影像上选取10个均匀分布的特征点,然后记录其在配准前、后两幅图像上的坐标信息,利用(8)式计算出RMSE。由表1可知:从RMSE 看,4种算法的配准精度也有一定幅度的提升,SIFT、SURF、BRISK算法配准精度分别提升了约0.17、0.20、0.14像素;其中对住宅区和混合区的影像配准结果基本控制在1.0像素内,道路区域以及在建工地区域的影像配准精度在1.5像素内,可以满足后续变化检测等处理的要求。

在运行速度方面,由于无人机影像分辨率较高,运算时间较长,为了方便计算,本文在算法运行前先对无人机影像进行降采样。对改进前、后算法对比后发现,对住宅区和混合区配准的4种算法改进后速度得到一定的提升,这是因为这2片区域的树冠面积较大,剔除这些树冠区域内的特征点再进行特征描述和匹配,可以大幅减少运算时间。而在道路附近和在建区中,特征点分布较少,在建区甚至没有树冠区域,而本文算法在特征粗匹配中采用双向匹配,使运算时间稍有增加,在后续的研究中,应该继续改进。

总体而言,本文提出的基于形态学改进算法对影像配准的精度有了一定的提升,其中对SIFT、SURF这类提取特征点较密集的传统配准算法提升效果明显;当影像上树冠区域密集时提升了算法的运算速度,对于没有树冠区域分布的影像,时间复杂度略有增加。

4 结 论

基于特征的影像配准算法中,特征描述子的提取与匹配是影像配准中的关键,决定了配准的精度与效率。本文从树冠区域内特征点匹配结果中发现,该区域内的特征点对大多是错误匹配,因此在特征点提取时通过形态学处理提取树冠区域,然后剔除该区域内的特征点再进行特征描述和匹配,提高了算法的特征匹配精度,也有效提高了配准结果的精度。该算法对SIFT、SURF这类特征点鲁棒性较差的算法,改进效果很好;而对于AKAZE这类特征点数量适中且鲁棒性强的算法,改进效果不是很明显。在特征匹配方面,本文结合相关研究,通过双向匹配算法提高了粗匹配精度,但是增加了运算时间,可以根据任务情况进行选用;通过PROSAC算法提高了精匹配的效率。与已有的算法相比,本文提出的改进算法在特征点匹配精度与配准精度上有了一定的提升,运算速度也基本满足要求。