技术领域
[0001] 本
发明涉及
图像处理技术领域，具体涉及到一种加强鲁棒性的无人机航拍图像拼接方法。
背景技术
[0002] 近年来，无人机航拍技术以其灵巧便捷、采集数据准确高效等优点，已成为各个领域获取信息的重要方式。由于无人机相机视
角和飞行高度的限制，单张航拍图像往往难以满足信息获取的需求，因此需要将两幅或多幅具有重叠部分的图像进行拼接融合，进而得到一幅高
精度、细节丰富的全景图像。
[0003]
图像配准是图像拼接中最重要的部分，其主要思想是，利用一组航拍图像中的重叠区域，寻找合适的数学模型和变换参数，建立图像之间的仿射关系，最终达到在同一
坐标系中进行特征点匹配的目的。根据所使用图像信息的区别，配准
算法可分为基于灰度的配准、基于特征的配准和基于图像理解与解释的配准三类，其中基于特征的图像配准算法更加稳定、快速，应用较为广泛，主要有Harris和SUSAN算子、SIFT、SURF、ORB和BRISK算法等等。尽管这些算法各有优点，但并不适用于实际的航拍图像配准，例如SIFT算法鲁棒性强但内存占用高、SURF算法综合性能强但容易丢失细节、ORB算法速度很快但鲁棒性较差。使用无人机获得的航拍图像往往
分辨率较高，而且
图像序列不在同一平面内，存在光照、尺度等的差异，这对航拍图像拼接算法的综合性能提出了很高的要求。
[0004] 本项目由国家自然科学基金项目(61773087)资助。
发明内容
[0005] 本发明的目的是针对当前无人机航拍图像拼接方法细节丢失严重、鲁棒性差的问题，提供了一种加强鲁棒性的无人机航拍图像拼接方法。
[0006] 为了达到上述目的，本发明采用的技术方案为：
[0007] 一种加强鲁棒性的无人机航拍图像拼接方法，包括如下步骤：
[0008] S1：输入两幅具有重叠区域的航拍图像，并进行预处理；
[0009] S2：使用边缘加强的SURF(Canny-SURF)特征检测算法处理航拍图像，得到两幅航拍图像的特征点；
[0010] S3：利用长距离FREAK(L-FREAK)算法对检测到的特征点进行特征描述，提取航拍图像特征；
[0011] S4：通过基于K近邻(k-Nearest Neighbor，KNN)算法对提取的航拍图像特征进行匹配；
[0012] S5：利用随机抽样一致(RANdom SAmple Consensus，RANSAC)算法对匹配后的特征点对进行提纯，提高匹配正确率，同时计算参考图像与待拼接图像之间的变换矩阵模型参数；
[0013] S6：使用坐标变换矩阵处理待拼接图像，并采用加权融合平均算法，使处理后的待拼接图像与参考图像在同一坐标系下进行加权平均融合，完成航拍图像拼接工作。
[0014] 进一步的，所述S1中，对图像进行预处理。采用Bilateral双边滤波对输入的航拍图像进行平滑处理，双边滤波表达式如下：
[0015]
[0016] 其中，k表示归一化因子，x和ξ分别为中心
像素点和邻域像素点，f(ξ)为邻域像素点的像素值，权重c(ξ,x)和s(f(ξ),f(x))与两像素点的几何和灰度距离负相关。
[0017] 进一步的，所述S2中，采用边缘加强的SURF(Canny-SURF)特征检测算法处理步骤S1得到的航拍图像，对航拍图像进行特征点检测。具体包括如下步骤：
[0018] S2.1：使用SURF算法检测图像特征点。
[0019] S2.1.1：构建Hessian矩阵。图像中的每一个特征点都可以求出一个特定的Hessian矩阵，其行列式为：
[0020]
[0021] 设l(x,y)为输入图像的像素点，则图像在尺度空间为σ的Hessian矩阵可以表示为：
[0022]
[0023] 其中函数L(x,σ)是高斯函数g(σ)的二阶导数 与像素点的卷积。
[0024] 为保证获取的特征点具备尺度不变特性，需要对输入图像进行滤波处理，降低图像灰度。SURF算法采用盒式滤波代替高斯滤波，通过加权系数来平衡精确值与近似值之间的误差，加权系数大小由尺度决定，最终得到Hessian矩阵行列式为：
[0025] det(H)＝LxxLyy-(ωLxy)2
[0026] 通过求解det(H)的特征值可以判断该点是否为特征点。
[0027] S2.1.2：构建尺度空间。尺度空间通常用高斯金字塔来表示，SURF算法在保证每组图像尺寸不变的前提下，逐步增大
滤波器的尺寸来构建高斯金字塔。
[0028] S2.1.3：特征点
定位。使用Hessian矩阵处理像素点之后，将得到的极值与其二维空间及邻域空间内的26个像素点进行比较，初步确定特征点
位置。然后使用差值运算进行特征点提纯，筛除定位不准确、不稳定的特性点，得到用于配准的特征点。
[0029] S2.1.4：特征点主方向分配。通过统计以特征点为中心的圆形区域内的Haar小波特征确定特征点主方向。具体来说，取圆心角为60°的扇形为单位，按照一定的间隔旋转一周，将每个单位扇形中的点在垂直和
水平方向的Haar小波响应矢量相加，取最大的作为特征点主方向。
[0030] S2.2：为了解决SURF算法边缘抑制较强，定位不准确等问题，使用Canny算法单独进行
边缘检测。
[0031] S2.2.1：梯度计算。像素梯度方向和幅值可以用x和y方向的一阶偏导的有限差分Gx和Gy确定：
[0032] 梯度方向：
[0033] 梯度幅值：
[0034] S2.2.2：特征点提纯。首先进行非极大值抑制，将目标像素在梯度方向θ上的幅值G与邻域内两个像素的幅值进行比较，若目标像素幅值最大则保留，否则舍弃。然后使用双
阈值法分割图像，去除由于噪声和灰度变化产生的虚假边缘。使用双阈值分割图像时，梯度值小于低阈值的点被抑制，剩余的像素点为强边缘和弱边缘，其中强边缘是真正的边缘，而弱边缘包含的一些噪声、灰度引起的虚假边缘需要去掉。
[0035] S2.2.3：滞后边缘
跟踪。通常强边缘和弱边缘中的真实边缘应该是相互连通的，而噪声引起的边缘一般是单独存在的。滞后边缘跟踪就是在弱边缘的3×3邻域内，若存在强边缘点就保留下来，否则抑制，通过不断递归跟踪，最终将边缘全部连接起来。
[0036] S2.3：使用Canny-SURF算法，即将所述S2.1和S2.2中得到的特征点集合并，并去除重复的特征点完成特征检测。
[0037] 进一步的，所述S3中利用长距离FREAK(L-FREAK)算法对所述S2中检测到的特征点进行描述，提取图像特征。具体步骤为：
[0038] S3.1：特征点
采样。L-FREAK算法是一种模仿人类视觉系统的二进制描述子，采样方式是以特征点为圆心，以指数递减的方式在特征点周围的圆形区域内寻找采样点。
[0039] S3.2：构建特征描述子。采样结束后，需要进行强度比对来获取二进制串生成FREAK描述子，描述子用F表示：
[0040]
[0041] 其中，α为
特征向量维度，N为期望的特征向量维数，Pα为采样点对，T(Pα)为采样点对的函数，定义如下：
[0042]
[0043] 其中， 和 分别Pα中两个采样点的灰度值。
[0044] 按照S3.1方式进行采样，每个特征点对应 个采样点对，会产生信息冗余，占用大量内存。所以在得到二进制描述子之后，需要对其进行筛选，提高描述子辨识度。通过构造矩阵，每列按照均值取出前521列作为特征描述子。
[0045] S3.3：基于长距离特征点的主方向确定。传统的FREAK算法是在在特征点周围选45个采样点提取方向信息，为了提高算法对旋转尺度较大的图片的鲁棒性，借鉴BRISK算法长距离点对的思想，L-FREAK算法仅使用远距离的30个采样点的方向信息进行主方向确定。这样特征点梯度O和角度信息θ的计算公式可以表示为：
[0046]
[0047] θ＝arctan(ox,oy)
[0048] 其中，M为选出的特征点对数，G为特征点对的集合，Po为采样点对的坐标信息，和 为采样点对Po中两个采样点的位置。
[0049] S3.4：特征粗匹配。对所述S3.3中生成的二进制描述子，通过参考图像与带拼接图像中的特征点之间的HANMING距离来衡量相似程度，得到粗略的匹配点对。
[0050] 进一步的，所述S4中的特征匹配，将经过S2和S3处理后得到的粗匹配点对使用KNN算法进行细匹配处理。具体为，设参考图像和待拼接图像的特征集合为A和B，对于A中任意描述子Ai，在B中找到与其距离最近和次近的描述子Bi和Bj，距离用d1和d2表示，若 则认为Ai对应的匹配点为Bi，阈值μ取值为2。
[0051] 进一步的，所述S5使用随机抽样一致算法(RANdom SAmple Consensus，RANSAC)对匹配后的特征点进行特征点提纯，并计算变换矩阵参数。具体步骤为：
[0052] S5.1：设t1＝(x1,y1,l)，t2＝(x2,y2,l)分别为同一景物在参考图像与待拼接图像中的投影，它们满足关系：
[0053] t2＝Ht1
[0054] 其中，矩阵H是从参考图像到待拼接图像的二维投影变换矩阵：
[0055] 这样从特征点(x1,y1)到(x2,y2)的透视变换公式为：
[0056]
[0057] 其中，m0,m1,m3,m4为图片的尺度和旋转量，m2为水平方向位移量，m5为垂直方向位移量，m6,m7为水平和垂直方向
变形量；共有8个待确定参数，需要4对匹配点进行求解，由于存在误匹配点，必须先经过提纯处理。
[0058] S5.2：首先任意选择4对匹配点(任意3对不共线)，求解变换矩阵H；然后利用求得的矩阵H计算其余匹配点在待拼接图像中的位置，并用d表示其于与实际对应点的垂直距离；设定一个阈值D，若满足d＜D则将其定义为内点；多次重复上述步骤，取其中包含内点最多的集合对应的H作为最优坐标变换矩阵。所述的阈值D一般根据所选的特定问题及数据集来确定。
[0059] 进一步的，所述S6中使用的加权平均融合方法具体为：
[0060] 假设f1(x,y)和f2(x,y)分别表示参考图像与经过S5投影变换之后的待拼接图像，则融合之后的图像f(x,y)可以表示为：
[0061]
[0062] 其中，a1和a2为像素点在两幅图像中分配到的权重，满足a1+a2＝1(a1,a2＞0)，a1可以由f1(x,y)到边界的距离与重叠区域的宽度的比值得到。a1由0变化到1，可以实现重叠区域光滑无缝融合，得到最终的全景航拍图像。
[0063] 本发明针对SURF算法在航拍图像拼接过程中出现的边缘定位不准确以及鲁棒性较差的问题，提出了一种加强鲁棒性的无人机航拍图像拼接算法。首先使用双边滤波来处理原始图像，利用Canny-SURF算法和L-FREAK二进制描述子进行特征点的检测和描述，加强边缘细节，将得到的特征点利用KNN算法进行粗匹配，然后采用RANSAC算法对匹配的特征点对提纯，并计算变换矩阵的参数，最后使用变换矩阵处理原始图像，并采用加权平均融合算法进行图像融合，完成航拍图像拼接。实验证明，本发明对拼接效果良好，对旋转变化较大的航拍图像具有更强的鲁棒性，同时相比其他经典算法综合性能更优，可以满足航拍图像拼接的实际要求。
附图说明
[0065] 图2是本发明
实施例所使用的航拍图像，其中(a)(b)(c)为参考图像，(d)(e)(f)为待拼接图像；另外(a)(d)、(b)(e)、(c)(f)分别为三组针对光照、尺度、旋转变化所挑选的航拍图像，用于验证算法的鲁棒性。
[0066] 图3是本发明实施例使用Canny-SURF算法和L-FREAK算法对图2中的三组航拍图像提取特征后得到的图像。
[0067] 图4是本发明实施例中对图3的航拍图像进行特征粗匹配后得到的图像。(a)为基于光照变化的航拍图像匹配结果，(b)为基于尺度变化的航拍图像匹配结果，(c)基于旋转变化的航拍图像匹配结果。
[0068] 图5是本发明实施例中使用RANSAC算法提纯之后的配准图像。(a)为基于光照变化的航拍图像配准结果，(b)为基于尺度变化的航拍图像配准结果，(c)基于旋转变化的航拍图像配准结果。
[0069] 图6是使用本发明所提算法拼接后最终得到的全景航拍图像。(a)为基于光照变化的航拍图像拼接结果，(b)为基于尺度变化的航拍图像拼接结果，(c)基于旋转变化的航拍图像拼接结果。
具体实施方式
[0070] 为了本技术领域的人员更好的理解本发明，下面结合附图和以下实施例对本发明作进一步详细描述。
[0071] 参阅图1，本实施例提供一种加强鲁棒性的无人机航拍图像拼接方法，包括如下步骤：
[0072] S1：如图2所示，输入两幅具有重叠区域的航拍图像，并进行预处理；
[0073] S2：使用边缘加强的SURF(Canny-SURF)特征检测算法处理航拍图像，得到两幅图像的特征点；
[0074] S3：利用长距离FREAK(L-FREAK)算法对检测到的特征点进行特征描述，提取图像特征，如图3所示；
[0075] S4：通过基于K近邻(k-Nearest Neighbor，KNN)算法对提取的航拍图像特征进行匹配，如图4所示；
[0076] S5：利用随机抽样一致(RANdom SAmple Consensus，RANSAC)算法对匹配后的特征点对进行提纯，提高匹配正确率，同时计算参考图像与待拼接图像之间的变换矩阵模型参数，得到的配准图像如图5所示；
[0077] S6：使用坐标变换矩阵处理待拼接图像，使之与参考图像在同一坐标系下进行加权平均融合，完成航拍图像拼接工作，如图6所示。
[0078] 具体的，所述S1中，航拍图像预处理。采用Bilateral双边滤波对图2所示的航拍图像进行平滑处理，双边滤波表达式如下：
[0079]
[0080] 其中x和ξ分别为中心像素点和邻域像素点，权重c(ξ,x)和s(f(ξ),f(x))与两像素点的几何和灰度距离负相关。
[0081] 所述S2中，采用边缘加强的SURF(Canny-SURF)算法对图像进行特征点检测。具体包括如下步骤：
[0082] S2.1：使用SURF算法检测图像特征点。
[0083] S2.1.1：构建Hessian矩阵。像素点l(x,y)在尺度空间为σ的Hessian矩阵可以表示为：
[0084]
[0085] 其中函数L(x,σ)是高斯函数g(σ)的二阶导数 与像素点的卷积。SURF算法采用了盒式滤波代替高斯滤波，最终得到Hessian矩阵行列式为：
[0086] det(H)＝LxxLyy-(ωLxy)2
[0087] 通过求解det(H)的特征值可以判断该点是否为特征点。
[0088] S2.1.2：构建尺度空间。尺度空间用高斯金字塔来表示，SURF算法在保证每组图像尺寸不变的前提下，逐步增大滤波器的尺寸来构建高斯金字塔。
[0089] S2.1.3：特征点定位。使用Hessian矩阵处理像素点之后，将得到的极值与其二维空间及邻域空间内的26个像素点进行比较，初步确定特征点位置。然后使用差值运算进行特征点提纯，筛除定位不准确、不稳定的特性点，得到用于配准的特征点。
[0090] S2.1.4：特征点主方向分配。取圆心角为60°的扇形为单位，按照一定的间隔旋转一周，将每个单位扇形中的点在垂直和水平方向的Haar小波响应矢量相加，取最大的作为特征点主方向。
[0091] S2.2：为了解决SURF算法边缘抑制较强，定位不准确等问题，使用Canny算法单独进行边缘检测。
[0092] S2.2.1：梯度计算。像素梯度方向和幅值可以用x和y方向的一阶偏导的有限差分Gx和Gy确定：
[0093] 梯度方向：
[0094] 梯度幅值：
[0095] S2.2.2：特征点提纯。首先进行非极大值抑制，将目标像素在梯度方向θ上的幅值G与邻域内两个像素的幅值进行比较，若幅值最大则保留，否则舍弃。然后使用双阈值法分割图像，去除由于噪声和灰度变化产生的虚假边缘。
[0096] S2.2.3：滞后边缘跟踪。滞后边缘跟踪就是在弱边缘的3×3邻域内，若存在强边缘点就保留下来，否则抑制，通过不断递归跟踪，最终将边缘全部连接起来。
[0097] S2.3：使用Canny-SURF算法，即将所述S2.1和S2.2中得到的特征点集合并，并去除重复的特征点完成特征检测。
[0098] 进一步的，所述S3中利用长距离FREAK(L-FREAK)算法对所述S2中检测到的特征点进行描述，提取图像特征。具体步骤为：
[0099] S3.1：特征点采样。L-FREAK算法采样方式是以特征点为圆心，以指数递减的方式在特征点周围的圆形区域内寻找采样点。
[0100] S3.2：构建特征描述子。采样结束后，进行强度比对来获取二进制串生成FREAK描述子，描述子用F表示：
[0101]
[0102] 其中，N为特征向量维数，T(Pα)为采样点对的函数，定义如下：
[0103]
[0104] 其中，I(Pα)为高斯模糊后采样点的灰度值。
[0105] 按照S3.1方式进行采样，每个特征点对应 个采样点对，会产生信息冗余，占用大量内存。通过构造矩阵，每列按照均值取出前521列作为特征描述子。
[0106] S3.3：基于长距离特征点的主方向确定。借鉴BRISK算法长距离点对的思想，L-FREAK算法仅使用远距离的30个采样点的方向信息进行主方向确定。特征点梯度O和角度信息θ的计算公式可以表示为：
[0107]
[0108] θ＝arctan(ox,oy)
[0109] 其中，M为选出的特征点对数，G为特征点对的集合，Po为采样点坐标信息，图3为提取到的特征点。
[0110] S3.4：特征粗匹配。对所述S3.3中生成的二进制描述子，通过参考图像与带拼接图像中的特征点之间的HANMING距离来衡量相似程度，得到粗略的匹配点对，如图4所示。
[0111] 进一步的，所述S4中的特征匹配，将图4中的粗匹配点对使用KNN算法进行细匹配处理。在待拼接图像中寻找位置最进和次近的对应点，若比值小于阈值则保留，本实施例中阈值μ取2。
[0112] 进一步的，所述S5使用RANSAC算法进行特征点提纯和变换矩阵参数计算。具体步骤为：
[0113] S5.1：矩阵H是从参考图像到待拼接图像的二维投影变换矩阵：
[0114]
[0115] 从特征点(x1,y1)到(x2,y2)的透视变换公式为：
[0116]
[0117] 共有8个待确定参数，需要4对匹配点进行求解，由于存在误匹配点，必须先经过提纯处理。
[0118] S5.2：首先任意选择4对匹配点(任意3对不共线)，求解变换矩阵H；然后利用求得的矩阵H计算其余匹配点在待拼接图像中的位置，并用d表示其于与实际对应点的垂直距离；设定一个阈值D，本实施例中取值为1，若满足d＜D则将其定义为内点；多次重复上述步骤，取其中包含内点最多的集合对应的H作为最优坐标变换矩阵，得到的配准图像如图5所示。
[0119] 进一步的，所述S6中使用加权平均融合方法具体为：
[0120] 使用融合算法之后的图像f(x,y)可以表示为：
[0121]
[0122] 其中，a1和a2为像素点在两幅图像中分配到的权重。a1由0变化到1，可以实现重叠区域光滑无缝融合，得到最终的全景航拍图像。
[0123] 仿真实验中使用的图片均为无人机实际拍摄图像，分辨率为4000*3000，如图2所示，每组两幅图像之间存在重叠场景，为了检验算法鲁棒性，三组图像分别是针对光照、旋转、尺度变化。图6为使用本发明所提算法与现有的SURF算法进行拼接得到的全景图片，可以看出本发明较好的完成了拼接任务，光滑过度、无明显拼接缝。
[0124] 本发明与基于SIFT算法、SURF算法的航拍图像拼接用时、鲁棒性对比表如表1、2、3所示，可以看出本发明的拼接速度和鲁棒性等综合性能要优于其他两种算法。
[0125] 表1特征点提取时间对比(s)
[0126]
[0127] 表2特征点匹配时间对比(s)
[0128]
[0129] 表3算法鲁棒性对比(％)
[0130]
[0131] 以上所述，仅为本发明的较佳实施例，并不用以限制本发明，本发明的
专利保护范围以
权利要求书为准，凡是运用本发明的
说明书及附图内容所作的等同结构变化，同理均应包含在本发明的保护范围内。
