无人机影像空间分辨率评估软件的设计与实现

曾安明，蒲德祥，唐 辉，杨 烽

（1.重庆市地理信息和遥感应用中心，重庆 401147；2.重庆邮电大学 通信与信息工程学院，重庆 400065）

重庆是“一带一路”和长江经济带的连接点，也是西部大开发的战略支点，随着“两地”、“两高”建设目标的推进，对高精度数字正射影像的需求与日俱增。重庆地貌复杂多样，东、南、北三面中高山环绕，中西部丘陵广布，河流、山脉沿线的地势起伏较大，最大落差接近3 000 m[1-2]。另外，重庆的云雾天气频繁，导致可见光波段的穿透性较差[3]，严重影响了卫星遥感、大飞机航拍的成像质量。载人直升机航空摄影虽然能达到要求，但成本较高且具有一定的风险，因此无人机摄影测量已经逐渐成为摄影测量数据获取的重要手段。

受技术手段以及区域内地形等因素影响，不同生产单位、不同批次的影像在分辨率方面的质量差异较大。王海风等[4]通过分析倾斜相机的结构特点，给出了面向倾斜影像的列向水平分辨率、列向垂直分辨率、行向水平分辨率的计算公式。曹洪涛等[5]基于倾斜摄影的成像原理，推算了多视角相机的分辨率数学模型。但这些方法没有顾及无人机的瞬时位姿影响，且缺少对地形起伏的考虑，所以对特殊地形条件下的无人机摄影测量无法完全适用。因此，目前重庆的无人机遥感影像检查主要通过人工判读的方法实现，需要消耗大量人力物力，且只能做到抽样检查，无法保证系统性、全面性。

本文以数字摄影测量基本原理为基础，综合考虑相机畸变、拍摄时的位置与姿态、地形起伏等，建立像素级分辨率评估计算的严密数学模型，并设计开发了一款评估软件。通过实验证明，本文方法具有自动化程度高、计算速度快、评估精度准等优势。

1 算法理论

本文提出的无人机影像空间分辨率评估算法，以带畸变差改正的共线方程为基础，在数字高程模型、无人机影像、传感器之间建立数学模型，可以准确地分析像素级的分辨率情况。

经典的畸变差改正模型[6]：

式中，x、y为像点坐标；Δx、Δy为像点坐标的改正值；x0、y0为像主点坐标；为像点到主点的距离；k1、k2为径向畸变系数；p1、p2为切向畸变系数；α为CCD像元非正方形比例因子；β为CCD阵列非正交畸变系数。

引入畸变差改正的共线方程：

式中，f为影像的内方位元素中的主距；Xs、Ys、Zs为摄站坐标；XA、YA、ZA为物方坐标；ai、bi、ci（i=1,2,3）为旋转矩阵的9个组成元素，可由3个外方位角元素计算得到。

经典的数字摄影测量有3种转角系统[7]，分别如下：

根据r值的正负大小可以分为高度正相关、中度正相关、低度正相关、高度负相关、中度负相关、低度负相关（见表5）。

1）φ，ω，κ转角系统，将偏角φ的旋转轴（Y轴）作为主轴，倾角ω的旋转轴作为第二轴，旋角κ的旋转轴作为第三轴。则：R=RφRωRκ。

2）φ′，ω′，κ′转角系统，将倾角ω′的旋转轴作为主轴，偏角φ′、旋角κ′的旋转轴分别作为第二轴、第三轴。则：R=Rφ′Rω′Rκ′。

3）A，ν，κ转角系统，主要用于单像摄影测量。将Z轴作为主轴，在旋转中具备空间方向不变的性质。ν为摄影光束与铅垂线的夹角，κ为像片主纵线与y轴之间的夹角。则：R=RARνRκ。

通过对主流软件进行分析（如Pix4dMapper等），其采用了第二种转角系统，即φ′，ω′，κ′转角系统。旋转矩阵计算方式如下：

在本文算法中，数字高程模型（DEM）提供（XA、YA、Za），外方位元素提供线元素（Xs、Ys、Zs）与角元素（φ′、ω′、κ′），内方位元素提供（x0、y0、f），畸变校正提供（k1、k2、p1、p2、α、β）。一般而言，α与β可以缺省，也能保证足够的精度。通过共线方程，可以建立数字高程模型坐标与图像像素之间的对应关系。对于DEM上一点，对应像素的（x、y）坐标可由共线方程计算得到，此处地面分辨率ρgsd的计算公式：

式中，ρccd为成像感光器件CCD的像元尺寸。

2 软件设计与实现

本文实现的无人机影像分辨率评估软件如图1所示，适用操作系统为Windows 7或Windows 10，编程语言为C#，开发工具为Visual Studio 2010，运行环境为.Net Framework 4.0。针对DEM数据的操作，基于ArcGIS Engine 10.0进行二次开发[8]。针对数字图像的操作，调用OpenCV功能实现[9]。

图1 无人机影像分辨率评估软件

软件实现步骤：

2）估算影像对应DEM范围，目的是减少程序遍历时的计算量。影像的长和宽分别为m，n个像素，则影像四个角的像平面坐标分为：左下角（-ρccd·m/2,-ρccd·n/2），左上角（-ρccd·m/2,+ρccd·n/2），右下角（+ρccd·m/2,-ρccd·n/2），右上角（+ρccd·m/2,+ρccd·n/2）。当ZA取DEM中高程最小值，利用公式（2）计算四个角对应的物方坐标，并求出最小外包盒对应的DEM像素行列号，分别是起始行号rb，终止行号re，起始列号cb，终止列号ce。

3）遍历DEM中第r（rb≤r≤re）行第c（cb≤c≤ce）列的像素，利用ArcGIS Engine的Pixel2Map函数得到（XA、YA），再利用GetPixelValue函数得到ZA。

4）根据公式（2）的共线方程，计算得到对应无人机影像的像素（x、y），并转换为行列号（i、j）。利用公式（5）计算（i、j）处的分辨率ρgsd，并将ρgsd作为像素值。

5）重复步骤3，直到遍历完成。生成无人机影像分辨率分布图。

6）针对无人机影像分辨率分布图中反投影法导致的一些孔隙，使用形态学处理中的膨胀算法[10]，可以较好地消除。

3 实验与分析

本文使用Pix4dMapper的测试数据集，如图所示。其中，图2a为原始的无人机影像，图2b为空三计算成果，图2c为生成的正射影像图，图2d为生成的数字表面模型，可以取代DEM参与计算，能更好表现建筑、植被等凸起地物对分辨率地影响。

图2 Pix4dMapper数据

内方位元素以及畸变改正参数记录在文件internal_camera_parameters中，外方位元素记录在文件external_camera_parameters中，提取后保存为软件可用的标准格式。

分辨率分布图如图3所示，像素灰度值差异代表了分辨率的大小，黑色区域表示无对应的地形数据。图中，地形起伏已经得到了反映，河流处的高度较小导致分辨率较低，所以灰度值较暗。堆体处的高度较大，导致分辨率较高，所以灰度值较亮。其次，影像边缘距离拍摄中心较远，所以分辨率会逐渐降低，灰度值会相应地逐渐变暗。

图3 分辨率分布图

软件还支持输出最小分辨率、最大分辨率、平均分辨率等统计信息，并可以按照质检需求进行拓展。

4 结 语

为了解决现有无人机影像空间分辨率评估算法在特殊地形条件下应用时的不足，本文综合考虑地形起伏、外方位元素、内方位元素、畸变改正数等影响因素，提出了一种高精度的像素级分辨率评估方法，并设计开发了评估软件。

本文中使用反投影的方法，遍历DEM寻找对应的影像像素位置，具有速度快的优点，但生成的分辨率分布图存在孔隙，虽然使用形态学膨胀算法可以有效改善，但进一步研究遍历影像寻找对应DEM像素的正投影算法，仍具有重要意义。