無人機影像空間分辨率評估軟件的設計與實現

曾安明，蒲德祥，唐 輝，楊 烽

（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像素的正投影算法，仍具有重要意義。