本文内容來源于《測繪學報》2022年第11期（審圖号GS京（2022）1229号）
無人機影像DSM自動生成随機傳播COLVLL算法
張春森1
, 葛英偉1, 郭丙軒2
, 張月瑩1
1. 西安科技大學測繪科學與技術學院，陝西 西安 710054;
2. 武漢大學測繪遙感資訊工程國家重點實驗室，湖北 武漢 430079
基金項目：國家自然科學基金(92038301)
摘要：針對現有密集比對方法在弱紋理及高差較大區域表現不佳, 以及密集比對融合生成DSM時導緻資訊丢失等問題, 本文提出基于随機傳播COLVLL的DSM生成方法。在對空三加密後的影像進行有效像對篩選的基礎上, 利用随機傳播機制對DSM像素區域進行掃描疊代, 結合VLL算法對随機生成的高程值進行疊代更新得到DSM。以弱紋理、大高差區域無人機影像為試驗資料與現有生成DSM商業軟體進行試驗對比, 并以ISPRS WGII/4提供的Vaihingen資料集為參考對本文方法生成DSM及真正射影像資料進行試驗分析, 結果均證明了本文方法的有效性和适用性。
關鍵詞：數字表面模型 無人機影像 密集比對 鉛垂線軌迹法 物方面元
引文格式：張春森, 葛英偉, 郭丙軒, 等. 無人機影像DSM自動生成随機傳播COLVLL算法[J]. 測繪學報，2022，51(11)：2346-2354. DOI: 10.11947/j.AGCS.2022.20210325
ZHANG Chunsen, GE Yingwei, GUO Bingxuan, et al. Automatic generation DSM of UAV image based on random propagation COLVLL algorithm[J]. Acta Geodaetica et Cartographica Sinica, 2022, 51(11): 2346-2354. DOI: 10.11947/j.AGCS.2022.20210325
閱讀全文：http://xb.chinasmp.com/article/2022/1001-1595/20221111.htm
引 言
數字表面模型(digital surface model，DSM)是指包含了地表建築物、樹木和橋梁等高度的地面高程模型。它在數字高程模型(digital elevation model，DEM)的基礎上，進一步涵蓋了除地面以外的其他高程的資訊。DSM資料廣泛應用于各行各業，如在森林地區可以用于檢測森林的生長情況，在城區可以用于檢查城市的發展情況，在軍事領域，可實作精确制導和自動化監控系統指揮作戰等[1]。通常通過影像密集比對算法求得視差或者深度圖，再利用視差或深度圖得到DSM。
無人機影像由于其影像幾何變形大、待比對資料量多等原因，在地物高差大、弱紋理區域影像比對會産生一系列問題[2]。如在地物遮擋處或高差大的地物上比對到的同名點過少，導緻解算影像視內插補點或深度值不穩定[3]。加之，現有密集比對算法在處理諸如弱紋理、大高差等對地觀測影像方面表現不佳。雖然對于密集比對國内外學者進行了大量相關研究與改進，文獻[4]提出半全局比對(semi-global matching, SGM)算法，利用雙重懲罰系數實作能量方程的最優求解，該算法通過計算8個方向或者16個方向的代價聚合實作全局能量最小化。文獻[5]在SGM算法中加入中值濾波對其進行改進，然而以上兩種方法在弱紋理區域均表現不佳。SGM采用代價聚合的能量函數來确定比對代價生成視差圖後融合得到DSM，但存在記憶體占用率高，視差圖中間結果流程複雜，計算量大，弱紋理、高差大時比對效果差。鉛垂線軌迹法(vertical line locus，VLL)是一種基于物方可直接得到地面點高程的影像比對方法，該方法能夠有效地避免輻射畸變對影像比對準确性的影響[6]，然而VLL法在地物深度不連續或高差較大的區域存在誤比對[7-9]。Colmap密集比對是基于随機物方patch面元掃描優化深度圖的方法，它基于相鄰像元的深度值法向量和目前像元的深度值法向量建構多個物方patch組合優化深度圖，并使用光學一緻性、幾何一緻性進行深度圖優化，通過過濾、融合得到最終深度圖[10-11]。基于物方patch面元的比對算法由文獻[12]提出，它通過引入随機采樣及随機傳播的機制，在影像間快速傳播并搜尋到相似影像塊進行比對。文獻[13]提出了PatchMatch雙目密集比對算法。近幾年來，随着深度學習的興起，卷積神經網絡被廣泛地應用于密集比對、同步定位與三維重建[14-17]。文獻[15]提出一種端到端的監督學習方式，通過建構代價體來表示左右影像的對應關系，利用3D卷積神經網絡對代價體進行學習得到視差圖。然而由于深度學習算法運作所需的内(顯)存較大，加之需耗費大量訓練和人工标注時間，其适用性和泛化性較差。
針對傳統密集比對方法在弱紋理及高差較大區域表現不佳，以及密集比對結果融合成DSM時導緻資訊丢失等問題，本文提出基于随機傳播COLVLL的無人機影像DSM快速生成方法，通過分析patch面元影像比對及随機傳播機制的思想，給出結合VLL算法對随機生成的高程值進行疊代更新的無人機影像DSM快速生成算法。分别選取城市高差大、弱紋理區域與一般區域無人機影像及ISPRS WGII/4提供的Vaihingen資料集進行試驗驗證，結果證明了本文算法的有效性。
1 本文算法
本文算法包括：①利用空三解算後的影像确定DSM測區生成範圍；②針對生成範圍進行掃描疊代，确定測區内關聯影像，利用候選影像篩選政策對候選影像進行篩選；③對候選影像進行高程、法向量等随機初始化；④對多個物方面元組合進行計算，利用VLL政策更新目前像素高程及法向量值；⑤在每次疊代結束後，利用相關系數比對計算得出最優組合，并進行疊代更新；⑥通過多次疊代更新最終得到DSM等步驟。算法流程如圖 1所示。
| |
| 圖 1 本文算法流程Fig. 1 Flowchart of the proposed algorithm |
| 圖選項 |
1.1 最佳候選影像篩選本文最佳候選影像篩選包括資料擷取、影像組像對打分、影像組像對篩選等步驟。(1) 資料擷取：準确的空三可提高比對初值的精度，并且保證比對的可靠性及降低疊代計算次數。針對空三加密後的資料确定DSM測區範圍，以及測區内最大、最小高程值。在最大最小高程範圍内，随機确定所有影像對應的高程值圖。(2) 影像組像對打分：以計算參考影像與每張比對影像的交會角得分、公共可視點個數作為其打分标準。設兩條光線在物方空間中的表示分别為U(a1, b1, c1)、V(a2, b2, c2)，則交會角θ為
(1)(3) 影像組像對篩選：基于交會角門檻值、公共可視點個數、比對影像門檻值，通過式(2)對比對影像組進行評價打分，确定最終篩選得到的最佳比對影像組
(2)式中，f為某一物方點；FR為基準影像可視物方點的集合；FI為其比對影像可視物方點的集合；Wa(f)為在物方點f處，比對影像I的交會角得分；Ws(f)為在物方點f處，比對影像I的分辨率相似性得分；Pc(f)為在物方點f處，比對影像I的公共可視點個數。1.2 随機PatchMatch掃描疊代PatchMatch算法最早被用于提高兩個圖像區域之間的近似最近鄰域的計算效率，通過大量随機采樣找到良好的patch比對。如圖 2所示，算法流程包括核線校正、深度資訊随機初始化、深度資訊全局傳播和深度資訊随機搜尋。
| 圖 2 PatchMatch算法流程Fig. 2 Flowchart of PatchMatch algorithm |
| 圖選項 |
PatchMatch算法使用比對代價估計的深度圖作為其融合密集點雲的中間結果，其深度圖的生成、優化、融合過程，流程複雜，計算量大，無法直接比對得到DSM成果。1.3 基于VLL物方高程自适應更新本文算法改進随機傳播深度值的計算方式，通過引入VLL算法，将深度值傳播改進為高程值傳播。利用随機PatchMatch對影像進行掃描疊代，同時結合VLL算法對随機生成的高程值進行疊代更新得到DSM。算法的核心是掃描疊代更新高程變化，自适應計算物方patch面元。如圖 3所示，物方地面三維點P(Xc, Yc, Zc)投影到參考影像I0上的像點為p0(x0, y0)，即以像點p0(x0, y0)為中心的影像視窗w0為物方patch面元的投影。過物點P(Xc, Yc, Zc)的物方patch面元由過物點P的法向量N所決定，而法向量N由法線方向角(α, β)确定。
| |
| 圖 3 像方與物方patch面元投影關系Fig. 3 The projection relationship between the image side and the object side patch |
| 圖選項 |
法向量N與法線方向角(α, β)之間的關系為
(3)若物方patch面元法向量為N(a, b, c)，則patch面元方程為
(4)過投影中心S0(XS0, YS0, ZS0)，以及參考像點p0(x0, y0)的投影光線方程為
(5)式中，(u0, v0, w0)為參考像點p0(x0, y0)在像空輔助坐标系下的坐标；λ為投影系數；(X, Y, Z)為投影光線與物方patch面元的交點坐标。像空間輔助坐标(u0, v0, w0)與參考像點p0(x0, y0)之間的關系為
(6)式中，R為由影像外方位角元素φ、ω、κ所構成的旋轉矩陣。物方patch面元坐标系與像方坐标系間的單應矩陣H為
(7)式中，K為影像内參數矩陣；(x0, y0)為物方patch面元中心坐标(自定義)；Rpatch是将patch方納入像方間坐标系下的旋轉矩陣；N為法向量；S為主距；C為影像外方位線元素；而(XspYspZsp)為平移向量。基于單應矩陣根據比對點位的變化自适應計算物方patch面元的高程及随機法向量的初值。通過疊代計算最優解，确定比對的高程值。歸納本文物方patch面元高程自适應計算步驟如下：(1) 通過物方patch面元坐标系與像方坐标系間的單應矩陣H實作像點p0與物方patch面元(法向量方向角為α、β)的投影。(2) 根據物方patch面元法向量随機值，求方向角α、β。(3) 在各候選影像上以像點p′i為中心，大小為w的視窗進行權重相關系數比對計算。(4) 根據得到的比對最優組合資訊更新目前點高程資訊。(5) 根據篩選出的最優組合确定目前像素的法向量、高程等資訊。(6) 繼續傳播疊代，最後得到目前比對點高程值。
2 試驗分析
2.1 試驗資料為驗證基于無人機影像資料采用本文算法在高差較大及弱紋理區域DSM生成效果，試驗資料(圖 4)分别選取安徽蚌埠某地18幅無人機影像，測區内高樓林立高差較大；四川綿陽某機場18幅無人機影像，測區内建築物影像呈弱紋理；江蘇揚中某地77幅無人機影像，測區為城市一般區域。3組影像資料相關資訊見表 1。算法實作基于Win7 64位作業系統，開發環境為VS 2017，采用C++語言程式實作。試驗硬體配置CPU為Intel Core 3.3 GHz，記憶體大小為16 GB。
| |
| 圖 4 試驗用無人機影像資料(部分)Fig. 4 UAV image data for test (part) |
| 圖選項 |
表 1 試驗資料相關資訊Tab. 1 Information about test data
表選項
2.2 大高差區域試驗針對第1組無人機影像資料，分别采用基于半全局SGM密集比對和本文算法生成DSM，以空三結束後的特征點作為控制點，統計兩組生成DSM後的平均誤差及均方根誤差(表 2)。圖 5為兩種算法得到的DSM結果對比，圖 6為兩種算法得到的陰影結果對比。其中，均方根誤差計算公式為
(8)
表 2 比對精度及效率比較Tab. 2 Comparison of matching accuracy and efficiency
表選項
| |
| 圖 5 大高差建築物區域DSM生成對比Fig. 5 Comparison of DSM generation in building area with large elevation difference |
| 圖選項 |
| |
| 圖 6 大高差建築物陰影對比Fig. 6 Shadow contrast of buildings with large elevation difference |
| 圖選項 |
改變面元視窗的尺寸大小，比較相關視窗尺寸對影像比對結果的影響(表 3)。
表 3 本文算法視窗尺寸對比對結果的影響
Tab. 3 Influence of window size on matching results of the proposed method
表選項
2.3 弱紋理區域試驗分别以SGM密集比對、視差圖融合得到的深度圖為基礎生成DSM，以及本文基于随機高程值傳播直接生成DSM方法對第2組弱紋理區域資料進行試驗(圖 7)。圖 8為兩種算法得到的DSM結果對比。圖 9為兩種算法得到的陰影結果對比。弱紋理區域比對精度及效率對比見表 4。
| |
| 圖 7 弱紋理影像(局部)Fig. 7 Weak texture image (local) |
| 圖選項 |
| |
| 圖 8 弱紋理區域DSM生成對比Fig. 8 Weak texture area DSM generation comparison |
| 圖選項 |
| |
| 圖 9 弱紋理區域陰影對比Fig. 9 Shadow contrast in weak texture areas |
| 圖選項 |
表 4 弱紋理區域比對精度及效率對比
Tab. 4 Comparison of matching accuracy and efficiency in weak texture areas
表選項
2.4 不同軟體DSM生成對比試驗基于第3組資料，采用本文算法與目前國内外主流商業軟體(SURE，PhotoScan，Pix4D)生成DSM進行試驗對比，如圖 10和圖 11所示。
| |
| 圖 10 不同算法DSM對比(整體效果)Fig. 10 Comparison of DSM effects of different methods(overall effect) |
| 圖選項 |
| |
| 圖 11 不同算法DSM對比(局部)Fig. 11 Comparison of DSM effects of different methods(local) |
| 圖選項 |
采用真正射影像(圖 12)在測區範圍選取一定數量特征點(10個)作為DSM定量分析檢查點(表 5)，其精度統計見表 6。
| |
| 圖 12 測區真正射影像Fig. 12 True orthophoto (TOP) of test area |
| 圖選項 |
表 5 定量分析檢查點統計(部分)Tab. 5 Quantitative analysis checkpoint statistics(part) m
表選項
表 6 不同軟體DSM生成精度統計
Tab. 6 Precision generated with different software DSM
m
表選項
2.5 基于ISPRS資料集的試驗為進一步驗證本文算法的有效性，使用ISPRS WGII/4提供的ISPRS Test Project on Urban Classification and 3D Building Reconstruc-tion Vaihingen資料集，對本文算法進行試驗驗證。該資料集包括由Leica ALS50系統擷取的機載雷射掃描器資料(平均點密度約為4個/m2)、數字表面模型(DSM)資料及真正射影像TOP資料(TOP和DSM的地面采樣距離均為9 cm)。研究區域包含典型的歐洲建築類型，其結構形态各異，包括山牆、坡屋頂及其混合結建構築。試驗以Vaihingen資料集中雷射掃描資料為真值，分别就本文算法生成的DSM與資料集中DSM，以及各自生成的真正射影像進行可視化對比與精度分析。(1) DSM生成對比。圖 13(a)、(b)、(c)是資料集中選取的3個研究區域1、2、3的DSM，圖 14(a)、(b)、(c)是使用本文算法生成對應區域DSM，圖 15(a)、(b)、(c)對應雷射掃描資料生成DSM。對比圖 13、圖 14、圖 15可以看出，由本文算法生成的DSM邊界清晰，更接近于雷射掃描資料生成DSM資料。
| |
| 圖 13 資料集DSMFig. 13 DSM of dataset |
| 圖選項 |
| |
| 圖 14 本文方法生成DSMFig. 14 DSM generated by the proposed algorithm |
| 圖選項 |
| |
| 圖 15 LAS資料生成DSMFig. 15 DSM generated by LAS |
| 圖選項 |
(2) 生成真正射影像對比。圖 16(a)、(b)是資料集中選取的兩個研究區域4和5比對結果生成真正射影像，圖 17(a)、(b)是采用本文算法生成對應區域的真正射影像，對比發現本文算法生成真正射影像效果較優。
| |
| 圖 16 資料集真正射影像Fig. 16 True orthophoto of dataset |
| 圖選項 |
| |
| 圖 17 本文算法生成真正射影像Fig. 17 True orthophoto generated by the paper algorithm |
| 圖選項 |
以Vaihingen資料集中雷射掃描資料為真值，分别與本文算法生成的DSM及資料集中比對生成(自帶)DSM進行對比，在DSM上随機擷取采樣資料點，計算不同地物處采樣點的均方根誤差與最大值，精度統計見表 7。由表 7可以看出，本文算法生成DSM的精度在山牆、坡屋頂及混合建築物等處均優于資料集比對生成DSM的精度。
表 7 與ISPRS資料對比精度分析
Tab. 7 Accuracy analysis compared with ISPRS data
m
表選項
2.6 試驗分析(1) 本文随機傳播COLVLL在patch方采用ZNCC雙邊濾波權重相關系數法比對得到較為準确的同名像點。由表 3可以看出，視窗大小在一定範圍内時(小于31×31像素)，視窗越大，比對到的正确值的機率也就越高。随着物方patch視窗的變大，計算的平均相關系數随之變大，即比對到正确值的精度在逐漸提高。試驗發現，當視窗尺寸取21×21像素時，可以取得較好的比對的精度與比對效率，說明patch面元視窗大小選取較為重要，如果選擇過小會導緻高程異常值多，過大會導緻比對效率過低。當patch視窗過大時，視窗内的紋理資訊增加，空間資訊複雜，導緻視窗相關系數求得的數值變小，比對像點的準确度降低，故視窗尺寸的選擇對比對結果的影響較大。(2) 由圖 5可以看出，本文算法能正确比對出高差較大的建築物，并且建築物邊緣完整。從建築物陰影圖上看(圖 6)，傳統方法比對的視差出現大量無效值，在進行深度圖融合及高程值融合時，前者導緻樓頂高程出現大量無效值。是以相較于其他方法，本文算法可以較好地處理高差較大區域。(3) 由于本文算法能夠較好地比對得到正确高程值，由圖 8弱紋理(原弱紋理影像如圖 7所示)影像DSM生成對比可以看出生成的DSM邊緣平滑，紋理清楚。同時由表 4可以看出，本文算法生成DSM精度及比對效率均高于傳統算法，說明本文算法可以較好地處理弱紋理區域影像。(4) 由表 6可以看出，本文算法DSM生成精度高于商業軟體，由圖 11可以看出，本文算法在建築物頂面處較其他3款商業軟體生成DSM更平滑，邊緣更清晰。(5) 由圖 13、圖 14、圖 15可以看出，本文算法生成的DSM邊界清晰，更接近于雷射掃描資料生成DSM資料。由表 7可以看出，本文算法在山牆、坡屋頂及混合建築物等處生成DSM精度均優于Vaihingen資料集比對生成DSM精度。由圖 16、圖 17可以看出，采用本文算法生成真正射影像優于Vaihingen資料集中比對結果生成真正射影像。
3 結論
借鑒Colmap基于随機物方patch面元掃描優化深度圖進而生成DSM的算法，本文提出了基于随機傳播COLVLL無人機影像DSM生成算法，即通過引入VLL算法，将深度值傳播改進為高程值傳播。利用随機PatchMatch對影像進行掃描疊代，同時結合VLL算法對随機生成的高程值進行多次疊代更新得到DSM。由于省略了深度圖等中間過程，是以克服了在弱紋理區域、高差範圍過大區域可能産生的誤比對，以及流程複雜，計算量大，記憶體占用率高等問題，在最終生成的DSM精度及效率方面較目前國内外主流商業軟體(SURE，PhotoScan，Pix4D)有了一定的改善。對于大場景超大資料量無人機影像DSM自動生成效率與精度的提升，将是筆者後續的研究工作。
作者簡介第一作者簡介：張春森(1963—)，男，博士，教授，研究方向為數字攝影測量、計算機視覺與遙感應用。E-mail：[email protected]通信作者：郭丙軒, E-mail：[email protected]
初審：張豔玲複審：宋啟凡
終審：金 君
資訊