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

（1.重慶市地理信息和遙感應用中心，重慶 401147；2.重慶郵電大學 通信與信息工程學院，重慶 400065）

重慶是“一帶一路”和長江經(jīng)濟帶的連接點，也是西部大開發(fā)的戰(zhàn)略支點，隨著“兩地”、“兩高”建設目標的推進，對高精度數(shù)字正射影像的需求與日俱增。重慶地貌復雜多樣，東、南、北三面中高山環(huán)繞，中西部丘陵廣布，河流、山脈沿線的地勢起伏較大，最大落差接近3 000 m[1-2]。另外，重慶的云霧天氣頻繁，導致可見光波段的穿透性較差[3]，嚴重影響了衛(wèi)星遙感、大飛機航拍的成像質(zhì)量。載人直升機航空攝影雖然能達到要求，但成本較高且具有一定的風險，因此無人機攝影測量已經(jīng)逐漸成為攝影測量數(shù)據(jù)獲取的重要手段。

受技術手段以及區(qū)域內(nèi)地形等因素影響，不同生產(chǎn)單位、不同批次的影像在分辨率方面的質(zhì)量差異較大。王海風等[4]通過分析傾斜相機的結構特點，給出了面向傾斜影像的列向水平分辨率、列向垂直分辨率、行向水平分辨率的計算公式。曹洪濤等[5]基于傾斜攝影的成像原理，推算了多視角相機的分辨率數(shù)學模型。但這些方法沒有顧及無人機的瞬時位姿影響，且缺少對地形起伏的考慮，所以對特殊地形條件下的無人機攝影測量無法完全適用。因此，目前重慶的無人機遙感影像檢查主要通過人工判讀的方法實現(xiàn)，需要消耗大量人力物力，且只能做到抽樣檢查，無法保證系統(tǒng)性、全面性。

本文以數(shù)字攝影測量基本原理為基礎，綜合考慮相機畸變、拍攝時的位置與姿態(tài)、地形起伏等，建立像素級分辨率評估計算的嚴密數(shù)學模型，并設計開發(fā)了一款評估軟件。通過實驗證明，本文方法具有自動化程度高、計算速度快、評估精度準等優(yōu)勢。

1 算法理論

本文提出的無人機影像空間分辨率評估算法，以帶畸變差改正的共線方程為基礎，在數(shù)字高程模型、無人機影像、傳感器之間建立數(shù)學模型，可以準確地分析像素級的分辨率情況。

經(jīng)典的畸變差改正模型[6]：

式中，x、y為像點坐標；Δx、Δy為像點坐標的改正值；x0、y0為像主點坐標；為像點到主點的距離；k1、k2為徑向畸變系數(shù)；p1、p2為切向畸變系數(shù)；α為CCD像元非正方形比例因子；β為CCD陣列非正交畸變系數(shù)。

引入畸變差改正的共線方程：

式中，f為影像的內(nèi)方位元素中的主距；Xs、Ys、Zs為攝站坐標；XA、YA、ZA為物方坐標；ai、bi、ci（i=1,2,3）為旋轉矩陣的9個組成元素，可由3個外方位角元素計算得到。

經(jīng)典的數(shù)字攝影測量有3種轉角系統(tǒng)[7]，分別如下：

根據(jù)r值的正負大小可以分為高度正相關、中度正相關、低度正相關、高度負相關、中度負相關、低度負相關（見表5）。

1）φ，ω，κ轉角系統(tǒng)，將偏角φ的旋轉軸（Y軸）作為主軸，傾角ω的旋轉軸作為第二軸，旋角κ的旋轉軸作為第三軸。則：R=RφRωRκ。

2）φ′，ω′，κ′轉角系統(tǒng)，將傾角ω′的旋轉軸作為主軸，偏角φ′、旋角κ′的旋轉軸分別作為第二軸、第三軸。則：R=Rφ′Rω′Rκ′。

3）A，ν，κ轉角系統(tǒng)，主要用于單像攝影測量。將Z軸作為主軸，在旋轉中具備空間方向不變的性質(zhì)。ν為攝影光束與鉛垂線的夾角，κ為像片主縱線與y軸之間的夾角。則：R=RARνRκ。

通過對主流軟件進行分析（如Pix4dMapper等），其采用了第二種轉角系統(tǒng)，即φ′，ω′，κ′轉角系統(tǒng)。旋轉矩陣計算方式如下：

在本文算法中，數(shù)字高程模型（DEM）提供（XA、YA、Za），外方位元素提供線元素（Xs、Ys、Zs）與角元素（φ′、ω′、κ′），內(nèi)方位元素提供（x0、y0、f），畸變校正提供（k1、k2、p1、p2、α、β）。一般而言，α與β可以缺省，也能保證足夠的精度。通過共線方程，可以建立數(shù)字高程模型坐標與圖像像素之間的對應關系。對于DEM上一點，對應像素的（x、y）坐標可由共線方程計算得到，此處地面分辨率ρgsd的計算公式：

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

2 軟件設計與實現(xiàn)

本文實現(xiàn)的無人機影像分辨率評估軟件如圖1所示，適用操作系統(tǒng)為Windows 7或Windows 10，編程語言為C#，開發(fā)工具為Visual Studio 2010，運行環(huán)境為.Net Framework 4.0。針對DEM數(shù)據(jù)的操作，基于ArcGIS Engine 10.0進行二次開發(fā)[8]。針對數(shù)字圖像的操作，調(diào)用OpenCV功能實現(xiàn)[9]。

圖1 無人機影像分辨率評估軟件

軟件實現(xiàn)步驟：

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函數(shù)得到（XA、YA），再利用GetPixelValue函數(shù)得到ZA。

4）根據(jù)公式（2）的共線方程，計算得到對應無人機影像的像素（x、y），并轉換為行列號（i、j）。利用公式（5）計算（i、j）處的分辨率ρgsd，并將ρgsd作為像素值。

5）重復步驟3，直到遍歷完成。生成無人機影像分辨率分布圖。

6）針對無人機影像分辨率分布圖中反投影法導致的一些孔隙，使用形態(tài)學處理中的膨脹算法[10]，可以較好地消除。

3 實驗與分析

本文使用Pix4dMapper的測試數(shù)據(jù)集，如圖所示。其中，圖2a為原始的無人機影像，圖2b為空三計算成果，圖2c為生成的正射影像圖，圖2d為生成的數(shù)字表面模型，可以取代DEM參與計算，能更好表現(xiàn)建筑、植被等凸起地物對分辨率地影響。

圖2 Pix4dMapper數(shù)據(jù)

內(nèi)方位元素以及畸變改正參數(shù)記錄在文件internal_camera_parameters中，外方位元素記錄在文件external_camera_parameters中，提取后保存為軟件可用的標準格式。

分辨率分布圖如圖3所示，像素灰度值差異代表了分辨率的大小，黑色區(qū)域表示無對應的地形數(shù)據(jù)。圖中，地形起伏已經(jīng)得到了反映，河流處的高度較小導致分辨率較低，所以灰度值較暗。堆體處的高度較大，導致分辨率較高，所以灰度值較亮。其次，影像邊緣距離拍攝中心較遠，所以分辨率會逐漸降低，灰度值會相應地逐漸變暗。

圖3 分辨率分布圖

軟件還支持輸出最小分辨率、最大分辨率、平均分辨率等統(tǒng)計信息，并可以按照質(zhì)檢需求進行拓展。

4 結 語

為了解決現(xiàn)有無人機影像空間分辨率評估算法在特殊地形條件下應用時的不足，本文綜合考慮地形起伏、外方位元素、內(nèi)方位元素、畸變改正數(shù)等影響因素，提出了一種高精度的像素級分辨率評估方法，并設計開發(fā)了評估軟件。

本文中使用反投影的方法，遍歷DEM尋找對應的影像像素位置，具有速度快的優(yōu)點，但生成的分辨率分布圖存在孔隙，雖然使用形態(tài)學膨脹算法可以有效改善，但進一步研究遍歷影像尋找對應DEM像素的正投影算法，仍具有重要意義。