无人机载荷图像地理信息拼接及验证算法*

梁中岩,戚红雨,王伟良,胡 杰

(1.中国电子科技集团公司第二十八研究所，江苏 南京 210007；2.中国人民解放军96902部队，北京 100015；3.中国人民解放军91604部队，山东 龙口 265700)

1 引言

由于灵活、机动的特点，无人机越来越多地应用于生产、生活的各个方面，如农业监控、电力监控和交通监控等[1]。侦察人员在得到无人机航拍图像的同时，往往需要知道关注目标的地理信息数据，如经度、纬度和高度等。现有的图像处理技术仅考虑图像中色彩、纹理等图像信息的处理，而未对地理信息进行处理，处理后的图像会丢失地理信息数据，最终使得处理后的图像仅有色彩、纹理信息，而不具有地理信息数据。为了满足侦察人员的需要，本文提出了一种图像与地理数据信息一体化的表示方法，使得地理数据信息处理与图像处理过程同步，实现图像与地理数据的一体化融合处理。同时，结合图像配准、融合和拼接技术，将地理信息融合到图像处理过程中，实现图像地理数据的同步处理，以便于侦察情报分析人员能够迅速从图像中直接获取感兴趣目标的地理位置。

2 相关工作

本文以图像拼接为例介绍带地理信息的无人机图像表示和处理算法。图像拼接的主要步骤包括图像配准和图像融合[2]。

图像配准主要通过检测图像中的特征点进行配准，常用的特征点检测算法主要是局部图像特征检测算法。近年来，有大量基于尺度不变特征变换SIFT(Scale Invariant Feature Transform)[3]、加速鲁棒特性SURF(Speeded Up Robust Features)[4]、有方向的加速分割测试特征和有旋转的二进制鲁棒独立基本特征ORB(Oriented FAST and Rotated BRIEF)[5]和关键点或角点检测Harris[6]的拼接算法被提出来，但基本都是仅仅处理图像本身。而在无人机遥感领域，无人机所拍摄的图像不仅包含图像本身，还带有各种地理信息和飞行参数。在无人机正视俯拍的情况下，可以利用飞机的载荷参数计算出图像中像素点的地理信息[7]，并将其与图像一同传回地面控制中心。在图像处理的同时，处理人员希望能够不丢失地理信息，例如，在拼接后的图像上可以更加直观地观察图像中的目标，对感兴趣的目标可以直接进行地理定位。

随着小型无人机的发展，无人机航拍图像的专业处理软件也蓬勃发展起来。其中，PhotoScan和Pix4D[8,9]可用于离线生成无人机航拍正射影像，然而在其处理过程中，多使用基于点云的稠密重建来生成全景影像，从而导致在图像插值的区域可能发生图像纹理扭曲的现象。虽然Arcgis[8]中也集成有无人机航拍图像的离线全景正射影像生成、无人机航迹点的地图显示等功能，然而该功能附属于Arcgis平台，主要支持民用的小型无人机数据与图像的传输、控制与处理，平台封闭性较强。

为了解决大量拼接计算时间长的问题，本文提出了基于分组控制的带地理信息的图像拼接算法，充分利用处理器多线程核心处理能力，对拼接图像的不同部分同时进行拼接，并在最终将各个处理器的处理结果拼接为最终拼接结果。

本文采用多维矩阵表示地理信息矩阵，参与图像拼接过程，与图像进行相同的仿射变换、平移操作，从而将地理信息拼接起来，并与原始图像一一对应。

3 地理信息数据表示方法

在计算机视觉领域，图像可由红、绿、蓝3个通道表示。由于地理信息与图像像素可一一对应，本文提出了一种使用图像存储结构存储并处理地理信息数据的方法。数字图像可由多通道无符号字符型(unsigned char，0～255)整数矩阵表示。本文采用双精度浮点多通道矩阵表示地理信息数据的主要信息，例如经度(单位为度(°)，-180°(不含)～180°)、纬度(单位为度(°)，-90°～90°)和高度(单位为m)。将地理信息存储为双精度浮点多通道矩阵，可利用矩阵处理算法对地理信息数据进行处理。

地理信息数据表示方法如图 1所示。

Figure 1 Geographic information data representation method图1 地理信息数据表示方法

地理信息数据多通道双精度浮点矩阵存储结构如图 2所示。

Figure 2 Multi-channel geographic information matrix storage structure图2 多通道地理信息矩阵存储结构

4 地理信息数据变换处理过程

地理信息数据变换处理过程如下所示：

(1) 使用多通道双精度浮点矩阵表示地理信息数据，记为矩阵Gorg；与之对应的图像数据矩阵记为Iorg。

(2) 通过归一化方法，将矩阵Gorg元素的取值范围转换到[0,1]，记为G，其中经度通道矩阵记为lon，纬度通道矩阵记为lat。

①经度通道归一化计算方法为：对于lon矩阵第y行x列的元素lonx,y=(lon0(x,y)+180)/360，其中lon0(x,y)为(-180,180]值域范围内的经度矩阵在y行x列的值。

②纬度通道归一化计算方法为：对于lat矩阵第y行x列的元素latx,y=(lat0(x,y)+90)/180，其中lat0(x,y)为[-90,90]值域范围内的纬度矩阵在y行x列的值。

③ 高度采用经纬度所在位置的数字高程模型DEM(Digital Elevation Model)高程数据，不参与计算。

(3) 通过归一化方法，将图像Iorg的取值范围转换到[0,1]，记为I。图像I归一化计算方法为：

(1)

(4) 从图像中提取特征，并在[0,1]进行图像处理变换。例如，对于仿射变换Γ，变换结果图像记为I′=ΓI。

(5) 矩阵G遵循相同的变换操作，包括平移、缩放、翻转、旋转、仿射等。例如，对于仿射变换Γ，变换结果矩阵记为G′=ΓG。

(6) 将结果图像I′和变换后的地理信息数据矩阵G′转换回原始值域空间。

5 基于ORB带地理信息的多线程无人机图像拼接算法

5.1 基于分组控制的带地理信息的图像拼接算法

为了解决大量拼接计算时间长的问题，本文提出基于ORB带地理信息的多线程无人机图像拼接算法，通过多线程分组控制大量图像的拼接过程，减少多幅航拍图像拼接过程中的误差累积，保证了无人机航拍图像大场景拼接的准确性。

无人机图像拼接任务往往面临大量的图像需要进行拼接处理，顺序拼接将耗时耗力，由于存在累积误差，拼接效果并不好。为了提高计算速度，充分利用计算资源，本文算法自动对拼接的图像进行分批，每个处理器核心运动一个线程，处理一批图像。多线程拼接过程如图 3所示。

每个线程处理的图像数计算方法为：

(2)

其中，l为每个线程处理的图像数，N为图像总数(N≥2)，b为不同批次之间重叠图像数(一般取b=1，取值范围为0≤b≤N-1，为了提高多线程拼接后对多个线程拼接结果进行合并的成功率，不同线程处理的图像之间需要有一定的重叠，重叠数并不是越大越好，典型值为1)，n为线程数(一般等于计算机CPU线程数)。图3中t为线程执行时间，tT(T=1,2,3,…,e)为线程执行时刻。

5.2 图像特征提取与变换矩阵计算

带地理信息的无人机图像处理算法首先需要通过图像处理算法对图像特征进行提取，而后在进行地理信息变换时共用图像处理过程计算得到的变换矩阵，如单应性矩阵H。本文以图像拼接为例介绍图像特征提取和变换矩阵计算过程。

Figure 3 Image mosaicking process with geographic information based on grouping control图3 基于分组控制的带地理信息的图像拼接过程

5.2.1 ORB特征提取

ORB算法[5]是基于关键点检测oFAST(orientation FAST)[10]和旋转敏感的二进制鲁棒独立基本特征BRIEF rBRIEF(rotation-aware Binary Robust Independent Elementary Features)[11]特征检测的一种新型特征提取算法。与SIFT算法相比，该算法具有计算量小、计算准确度满足使用需求的特点。

对图像进行ORB特征提取的步骤如下所示：

(1) 对图像构建金字塔。

(2) 用加速分割测试特征FAST(Features from Accelerated Segment Test)算法检测关键点的位置，这里的关键点指角点，角点为邻域内具有2个主方向的特征点。

(4) 对于选出的M个角点，利用强度中心(Intensity Centroid)算法[13]计算角点的方向，得到oFAST特征。

(5) 由于BRIEF算法是无向的，将(4)中计算出来的角点方向作为BRIEF的方向进行旋转，就得到了有向的BRIEF，并用贪婪学习算法筛选出具有高方差和高不相关的有向BRIEF，称之为rBRIEF。

(6) ORB特征就是oFAST和rBRIEF的组合。

5.2.2 特征匹配

ORB特征匹配步骤如下所示：

(1) 使用k近邻算法[14]计算任意2幅图像之间的特征点距离，记特征点最近距离为m0，次近距离为m1,如果r=m0/m1<δ，那么特征点为匹配点，匹配的特征点记为〈p0,p1〉；对任意图像进行特征点的匹配，并删除置信度低(即认为不是同一个全景图中的)的图像。

(2) 通过最小采样一致性算法RANSAC(RANdom SAmple Consensus)[15]进行优选。

(4) 将2幅可以拼接的图像放到拼接集合Φ={I0,I1}中，并不断扩展该集合形成最大的拼接集Φ={I0,I1,…,Im}。

5.2.3 光束法平差

由于多个单应性矩阵合成全景拼接图像时会造成累积误差，因此需要为每幅图像加上光束法平差值，以初始化图像为相同的旋转和焦距长度。使用光束法平差(Bundle Adjustment)[16]算法对所有图像进行相机参数K校正，以提高估计精度。光束法平差具有鲁棒性，其目标函数是一个映射误差的平方和函数,即每一个特征点都要映射到其他的图像中，使计算出的相机参数的误差的平方和最小。

5.3 带地理信息数据的图像变换及多线程拼接结果融合

根据相机参数和旋转矩阵，对图像进行单应性矩阵变换，再对图像进行初始拼接。地理信息矩阵也使用对应图像的旋转参数进行单应性矩阵变换，以对地理信息进行拼接。地理信息跟随图像进行相同的单应性矩阵变换，保证了变换后的图像有相同的地理信息。同时，由于经度为-180°(不含)～180°，纬度为-90°～90°，在插值变换时需要注意地理信息不能超出有效范围。

对拼接后的图像进行亮度的增量补偿以及基于图像金字塔的多波段融合[17]。对拼接后的地理信息进行有条件的插值融合。如果在拼接过程中，其中一幅图像全部或部分区域地理信息未知，而另一幅图像相同区域的地理信息已知，则直接将已知的地理信息作为拼接后图像的地理信息；如果拼接的2幅图像在融合点处的经纬度相差不大(小于0.000 2°)，则取两者的平均值作为拼接图像的经纬度。

带地理信息数据的图像变换步骤为：

(1) 对重叠区域进行分块估计平均光强。

(2) 用最大流算法检测缝隙。

(3) 进行多波段融合，将原图像分解成多个不同空间分辨率、不同尺度的子图像构成金字塔，然后由各层金字塔分别进行融合，最后组合得到拼接图像。

(4) 对拼接后的地理信息进行有条件的插值融合。

最后，对每一个处理器核心处理的拼接结果进行综合拼接融合，形成拼接结果图像。

6 基于SIFT带地理信息的多线程无人机图像拼接验证算法

对于图像拼接结果的验证，本文提出一种基于SIFT特征点对拼接后的地理信息进行验证的方法。该方法使用拼接前的图像与拼接后的图像进行SIFT特征匹配，用随机采样一致性RANSAC算法选出匹配度高的特征点，并比较拼接前和拼接后特征点地理信息的差异，如果其差异小于或等于0.000 1，则认为拼接的地理信息正确。具体步骤为：

(1) 提取拼接前图像的SIFT特征。

(2) 提取拼接后图像的SIFT特征。

(3) 使用RANSAC算法[15]选择高分特征点。

(4)比较匹配特征点之间的差值d，如果d≤0.0001(单位为°，在纬度尺度上大约11.1 m)，那么认为拼接结果是正确的。

7 实验结果与分析

为了验证本文算法的有效性，本文收集了不同场景下大量无人机航拍影像，分别对带地理信息的无人机图像进行性能测试、拼接结果测试和拼接结果验证。在实验之前，本文已通过无测距定位算法[7]对无人机图像中携带的地理信息进行了逐像素点定位，将每一像素对应的经纬度信息存储在多通道双精度浮点矩阵中。

7.1 性能测试

测试平台采用Intel i7-3770 3.4 GHz 4核8线程CPU，8 GB内存。本文分别使用单线程连续拼接算法、不使用多线程分组控制的本文算法(单线程)和使用多线程分组控制的本文算法做比较，测试结果如表 1所示,图像分辨率为720×576。

Table 1 Multi-thread and single-thread image mosaicking performance comparison表1 多线程与单线程拼接性能对比

从表 1可以看出，基于分组控制的带地理信息的图像拼接算法大幅提高了图像拼接速度和准确度，由于累积误差的存在，单线程拼接算法在拼接50幅以上图像时会因为畸变过大而失败。

单线程连续拼接是指多幅图像按拍摄顺序一幅一幅拼接，每次新拼接的图像只与之前的拼接结果做比较计算单应性矩阵H，而该矩阵变换模型的估计值与真实值之间始终存在误差。

分析单应性矩阵，第q幅图像与第r幅图像(r

(3)

其中,Hi,i+1为第i幅图像至第i+1幅图像的单应性矩阵估计值。由于Hi,i+1为矩阵变量，矩阵乘法是矩阵元素的加法和乘法的混合运算，随机误差在加法和乘法中传递。

Figure 4 Image mosaicking results of 100 UAV images图4 100幅无人机图像拼接结果

(1)随机误差在加法中的传递[19]。

(4)

即随机变量的和的方差等于各个方差之和。

(2)随机误差在乘法中的传递[19]。

对于X=A×B×C×…有:

(5)

即随机变量积的相对偏差的平方等于各个相对偏差的平方之和。

这种误差的累积效应会导致拼接误差越来越大[18,20]，最终得到的配准点无法用于计算单应性矩阵，从而导致拼接失败。拼接失败是指拼接程序不能完成指定图像数量的拼接，例如，表1中50幅单线程连续拼接失败的原因：未完成50幅图像的拼接，只拼接了前45幅。

7.2 图像拼接结果

使用本文基于分组控制的带地理信息的图像拼接算法拼接的100幅图像结果如图 4所示。

7.3 地理信息验证实验结果

基于SIFT特征点对拼接后的地理信息进行验证。对拼接前的图像与拼接后的图像进行SIFT特征匹配，用RANSAC算法选出匹配度高的特征点，并比较拼接前和拼接后特征点地理信息的差异，如果其差异小于或等于0.000 1°，则认为拼接的地理信息正确。验证结果如图 5所示，左图为右图拼接前的图像，右图为100幅图像拼接后的图像，共220个匹配点，纵向距离约4.463 km，部分匹配点差值如表 2所示(真实经纬度的整数部分分别替换为YY和XX，距离差值通过Haversine公式[21,22]计算)，均在有效误差范围内。

Figure 5 SIFT feature points matching verification图5 SIFT特征点匹配验证

从图 5和表 2可以看出，本文使用SIFT算法自动提取拼接前和拼接后图像的关键点，关键点的地理信息数据精度达到了算法要求，表示本文提出的算法在处理地理信息数据方面有效。

表2中距离差值通过Haversine公式[21,22]计算：

(6)

Table 2 Key points latitude and longitude matching results before and after mosaicking (partial matched points)表2 拼接前后图像匹配点部分关键点经纬度匹配结果(部分匹配点)

D=Δσ×RE

(7)

其中，假设地球是一个完美的球体，RE≈6371 km为地球的平均半径，Δφ为拼接前后图像中对应的2点纬度弧度差值的绝对值，Δγ为拼接前后图像中对应的2点经度弧度差值的绝对值，Δσ为拼接前后图像中对应的地球球面上2点分别与球心的连线所构成的圆心角，φA和φB分别为拼接前和拼接后对应2点的纬度，D为2点间的球面距离。

与其他算法的准确度比较如表3所示，其中本文算法数据是在220个匹配点上计算得出的。

从表3可以看出，文献[23]所述算法实验结果经纬度差值最大小于0.000 2°；文献[24]所述算法实验结果经纬度差值最大小于0.002°；文献[20]所述算法实验结果拼接平均距离差值在56.9 cm(0.000 569 km)，最大距离差值在文献[20]中未明确指出，表中数据为根据图表估算；本文算法实验结果经纬度差值最大小于或等于0.000 1°(在纬度尺度上大约11.1 m)优于文献[23,24]算法的计算结果，距离商用软件PhotoScan和Pix4D结果[25]还有一定的差距。

Table 3 Accuracy comparison of the proposed algorithm with other algorithms表3 本文算法与其他算法的准确度比较

8 结束语

本文首先提出了一种基于ORB带地理信息的无人机图像拼接算法。该算法利用地理信息数据和图像像素点之间一一对应的关系，并将地理信息的经度、纬度和高度使用多通道双精度浮点矩阵存储，使得图像和矩阵处理算法可以用于处理地理信息数据。同时，应用基于分组控制的多线程处理算法，提高了100幅图像拼接时的执行速度和准确性。最后，本文提出了一种基于SIFT算法的关键点拼接结果验证方法，利用该方法可以验证拼接结果的有效性。