張春森，葛英偉，郭丙軒，張月瑩

1. 西安科技大學測繪科學與技術學院，陜西 西安 710054； 2. 武漢大學測繪遙感信息工程國家重點實驗室，湖北 武漢 430079

數字表面模型(digital surface model，DSM)是指包含了地表建筑物、樹木和橋梁等高度的地面高程模型。它在數字高程模型(digital elevation model，DEM)的基礎上，進一步涵蓋了除地面以外的其他高程的信息。DSM數據廣泛應用于各行各業(yè)，如在森林地區(qū)可以用于檢測森林的生長情況，在城區(qū)可以用于檢查城市的發(fā)展情況，在軍事領域，可實現精確制導和自動化監(jiān)控系統(tǒng)指揮作戰(zhàn)等[1]。通常通過影像密集匹配算法求得視差或者深度圖，再利用視差或深度圖得到DSM。

無人機影像由于其影像幾何變形大、待匹配數據量多等原因，在地物高差大、弱紋理區(qū)域影像匹配會產生一系列問題[2]。如在地物遮擋處或高差大的地物上匹配到的同名點過少，導致解算影像視差值或深度值不穩(wěn)定[3]。加之，現有密集匹配算法在處理諸如弱紋理、大高差等對地觀測影像方面表現不佳。雖然對于密集匹配國內外學者進行了大量相關研究與改進，文獻[4]提出半全局匹配(semi-global matching,SGM)算法，利用雙重懲罰系數實現能量方程的最優(yōu)求解，該算法通過計算8個方向或者16個方向的代價聚合實現全局能量最小化。文獻[5]在SGM算法中加入中值濾波對其進行改進，然而以上兩種方法在弱紋理區(qū)域均表現不佳。SGM采用代價聚合的能量函數來確定匹配代價生成視差圖后融合得到DSM，但存在內存占用率高，視差圖中間結果流程復雜，計算量大，弱紋理、高差大時匹配效果差。鉛垂線軌跡法(vertical line locus，VLL)是一種基于物方可直接得到地面點高程的影像匹配方法，該方法能夠有效地避免輻射畸變對影像匹配準確性的影響[6]，然而VLL法在地物深度不連續(xù)或高差較大的區(qū)域存在誤匹配[7-9]。Colmap密集匹配是基于隨機物方patch面元掃描優(yōu)化深度圖的方法，它基于相鄰像元的深度值法向量和當前像元的深度值法向量構建多個物方patch組合優(yōu)化深度圖，并使用光學一致性、幾何一致性進行深度圖優(yōu)化，通過過濾、融合得到最終深度圖[10-11]?；谖锓絧atch面元的匹配算法由文獻[12]提出，它通過引入隨機采樣及隨機傳播的機制，在影像間快速傳播并搜索到相似影像塊進行匹配。文獻[13]提出了PatchMatch雙目密集匹配算法。近幾年來，隨著深度學習的興起，卷積神經網絡被廣泛地應用于密集匹配、同步定位與三維重建[14-17]。文獻[15]提出一種端到端的監(jiān)督學習方式，通過構建代價體來表示左右影像的對應關系，利用3D卷積神經網絡對代價體進行學習得到視差圖。然而由于深度學習算法運行所需的內(顯)存較大，加之需耗費大量訓練和人工標注時間，其適用性和泛化性較差。

針對傳統(tǒng)密集匹配方法在弱紋理及高差較大區(qū)域表現不佳，以及密集匹配結果融合成DSM時導致信息丟失等問題，本文提出基于隨機傳播COLVLL的無人機影像DSM快速生成方法，通過分析patch面元影像匹配及隨機傳播機制的思想，給出結合VLL算法對隨機生成的高程值進行迭代更新的無人機影像DSM快速生成算法。分別選取城市高差大、弱紋理區(qū)域與一般區(qū)域無人機影像及ISPRS WGII/4提供的Vaihingen數據集進行試驗驗證，結果證明了本文算法的有效性。

1 本文算法

本文算法包括：①利用空三解算后的影像確定DSM測區(qū)生成范圍；②針對生成范圍進行掃描迭代，確定測區(qū)內關聯影像，利用候選影像篩選策略對候選影像進行篩選；③對候選影像進行高程、法向量等隨機初始化；④對多個物方面元組合進行計算，利用VLL策略更新當前像素高程及法向量值；⑤在每次迭代結束后，利用相關系數匹配計算得出最優(yōu)組合，并進行迭代更新；⑥通過多次迭代更新最終得到DSM等步驟。算法流程如圖1所示。

圖1 本文算法流程Fig.1 Flowchart of the proposed algorithm

1.1 最佳候選影像篩選

本文最佳候選影像篩選包括數據獲取、影像組像對打分、影像組像對篩選等步驟。

(1) 數據獲?。簻蚀_的空三可提高匹配初值的精度，并且保證匹配的可靠性及降低迭代計算次數。針對空三加密后的數據確定DSM測區(qū)范圍，以及測區(qū)內最大、最小高程值。在最大最小高程范圍內，隨機確定所有影像對應的高程值圖。

(2) 影像組像對打分：以計算參考影像與每張匹配影像的交會角得分、公共可視點個數作為其打分標準。

設兩條光線在物方空間中的表示分別為U(a1,b1,c1)、V(a2,b2,c2)，則交會角θ為

(1)

(3) 影像組像對篩選：基于交會角閾值、公共可視點個數、匹配影像閾值，通過式(2)對匹配影像組進行評價打分，確定最終篩選得到的最佳匹配影像組

gR(I)=∑f∈FR∩FIWa(f)Ws(f)pc(f)

(2)

式中，f為某一物方點；FR為基準影像可視物方點的集合；FI為其匹配影像可視物方點的集合；Wa(f)為在物方點f處，匹配影像I的交會角得分；Ws(f)為在物方點f處，匹配影像I的分辨率相似性得分；Pc(f)為在物方點f處，匹配影像I的公共可視點個數。

1.2 隨機PatchMatch掃描迭代

PatchMatch算法最早被用于提高兩個圖像區(qū)域之間的近似最近鄰域的計算效率，通過大量隨機采樣找到良好的patch匹配。如圖2所示，算法流程包括核線校正、深度信息隨機初始化、深度信息全局傳播和深度信息隨機搜索。

圖2 PatchMatch算法流程Fig.2 Flowchart of PatchMatch algorithm

PatchMatch算法使用匹配代價估計的深度圖作為其融合密集點云的中間結果，其深度圖的生成、優(yōu)化、融合過程，流程復雜，計算量大，無法直接匹配得到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由法線方向角(α,β)確定。

法向量N與法線方向角(α,β)之間的關系為

(3)

圖3 像方與物方patch面元投影關系Fig.3 The projection relationship between the image side and the object side patch

若物方patch面元法向量為N(a,b,c)，則patch面元方程為

a(X-Xc)+b(Y-Yc)+c(Z-Zc)=0

(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面元的高程及隨機法向量的初值。通過迭代計算最優(yōu)解，確定匹配的高程值。

歸納本文物方patch面元高程自適應計算步驟如下：

(1) 通過物方patch面元坐標系與像方坐標系間的單應矩陣H實現像點p0與物方patch面元(法向量方向角為α、β)的投影。

(2) 根據物方patch面元法向量隨機值，求方向角α、β。

(4) 根據得到的匹配最優(yōu)組合信息更新當前點高程信息。

(5) 根據篩選出的最優(yōu)組合確定當前像素的法向量、高程等信息。

(6) 繼續(xù)傳播迭代，最后得到當前匹配點高程值。

2 試驗分析

2.1 試驗數據

為驗證基于無人機影像數據采用本文算法在高差較大及弱紋理區(qū)域DSM生成效果，試驗數據(圖4)分別選取安徽蚌埠某地18幅無人機影像，測區(qū)內高樓林立高差較大；四川綿陽某機場18幅無人機影像，測區(qū)內建筑物影像呈弱紋理；江蘇揚中某地77幅無人機影像，測區(qū)為城市一般區(qū)域。3組影像數據相關信息見表1。算法實現基于Win7 64位操作系統(tǒng)，開發(fā)環(huán)境為VS 2017，采用C++語言程序實現。試驗硬件配置CPU為Intel Core 3.3 GHz，內存大小為16 GB。

表1 試驗數據相關信息

圖4 試驗用無人機影像數據(部分)Fig.4 UAV image data for test (part)

2.2 大高差區(qū)域試驗

針對第1組無人機影像數據，分別采用基于半全局SGM密集匹配和本文算法生成DSM，以空三結束后的特征點作為控制點，統(tǒng)計兩組生成DSM后的平均誤差及均方根誤差(表2)。圖5為兩種算法得到的DSM結果對比，圖6為兩種算法得到的陰影結果對比。其中，均方根誤差計算公式為

(8)

表2 匹配精度及效率比較

改變面元窗口的尺寸大小，比較相關窗口尺寸對影像匹配結果的影響(表3)。

表3 本文算法窗口尺寸對匹配結果的影響

圖5 大高差建筑物區(qū)域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

2.3 弱紋理區(qū)域試驗

分別以SGM密集匹配、視差圖融合得到的深度圖為基礎生成DSM，以及本文基于隨機高程值傳播直接生成DSM方法對第2組弱紋理區(qū)域數據進行試驗(圖7)。圖8為兩種算法得到的DSM結果對比。圖9為兩種算法得到的陰影結果對比。弱紋理區(qū)域匹配精度及效率對比見表4。

圖7 弱紋理影像(局部)Fig.7 Weak texture image (local)

圖8 弱紋理區(qū)域DSM生成對比Fig.8 Weak texture area DSM generation comparison

圖9 弱紋理區(qū)域陰影對比Fig.9 Shadow contrast in weak texture areas

表4 弱紋理區(qū)域匹配精度及效率對比

2.4 不同軟件DSM生成對比試驗

基于第3組數據，采用本文算法與當前國內外主流商業(yè)軟件(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)在測區(qū)范圍選取一定數量特征點(10個)作為DSM定量分析檢查點(表5)，其精度統(tǒng)計見表6。

圖12 測區(qū)真正射影像Fig.12 True orthophoto (TOP) of test area

表5 定量分析檢查點統(tǒng)計(部分)

表6 不同軟件DSM生成精度統(tǒng)計

2.5 基于ISPRS數據集的試驗

為進一步驗證本文算法的有效性，使用ISPRS WGII/4提供的ISPRS Test Project on Urban Classification and 3D Building Reconstruc-tion Vaihingen數據集，對本文算法進行試驗驗證。該數據集包括由Leica ALS50系統(tǒng)獲取的機載激光掃描儀數據(平均點密度約為4個/m2)、數字表面模型(DSM)數據及真正射影像TOP數據(TOP和DSM的地面采樣距離均為9 cm)。研究區(qū)域包含典型的歐洲建筑類型，其結構形態(tài)各異，包括山墻、坡屋頂及其混合結構建筑。試驗以Vaihingen數據集中激光掃描數據為真值，分別就本文算法生成的DSM與數據集中DSM，以及各自生成的真正射影像進行可視化對比與精度分析。

(1) DSM生成對比。圖13(a)、(b)、(c)是數據集中選取的3個研究區(qū)域1、2、3的DSM，圖14(a)、(b)、(c)是使用本文算法生成對應區(qū)域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)是數據集中選取的兩個研究區(qū)域4和5匹配結果生成真正射影像，圖17(a)、(b)是采用本文算法生成對應區(qū)域的真正射影像，對比發(fā)現本文算法生成真正射影像效果較優(yōu)。

圖16 數據集真正射影像Fig.16 True orthophoto of dataset

圖17 本文算法生成真正射影像 Fig.17 True orthophoto generated by the paper algorithm

以Vaihingen數據集中激光掃描數據為真值，分別與本文算法生成的DSM及數據集中匹配生成(自帶)DSM進行對比，在DSM上隨機獲取采樣數據點，計算不同地物處采樣點的均方根誤差與最大值，精度統(tǒng)計見表7。由表7可以看出，本文算法生成DSM的精度在山墻、坡屋頂及混合建筑物等處均優(yōu)于數據集匹配生成DSM的精度。

表7 與ISPRS數據對比精度分析

2.6 試驗分析

(1) 本文隨機傳播COLVLL在patch方采用ZNCC雙邊濾波加權相關系數法匹配得到較為準確的同名像點。由表3可以看出，窗口大小在一定范圍內時(小于31×31像素)，窗口越大，匹配到的正確值的概率也就越高。隨著物方patch窗口的變大，計算的平均相關系數隨之變大，即匹配到正確值的精度在逐漸提高。試驗發(fā)現，當窗口尺寸取21×21像素時，可以取得較好的匹配的精度與匹配效率，說明patch面元窗口大小選取較為重要，如果選擇過小會導致高程異常值多，過大會導致匹配效率過低。當patch窗口過大時，窗口內的紋理信息增加，空間信息復雜，導致窗口相關系數求得的數值變小，匹配像點的準確度降低，故窗口尺寸的選擇對匹配結果的影響較大。

(2) 由圖5可以看出，本文算法能正確匹配出高差較大的建筑物，并且建筑物邊緣完整。從建筑物陰影圖上看(圖6)，傳統(tǒng)方法匹配的視差出現大量無效值，在進行深度圖融合及高程值融合時，前者導致樓頂高程出現大量無效值。因此相較于其他方法，本文算法可以較好地處理高差較大區(qū)域。

(3) 由于本文算法能夠較好地匹配得到正確高程值，由圖8弱紋理(原弱紋理影像如圖7所示)影像DSM生成對比可以看出生成的DSM邊緣平滑，紋理清楚。同時由表4可以看出，本文算法生成DSM精度及匹配效率均高于傳統(tǒng)算法，說明本文算法可以較好地處理弱紋理區(qū)域影像。

(4) 由表6可以看出，本文算法DSM生成精度高于商業(yè)軟件，由圖11可以看出，本文算法在建筑物頂面處較其他3款商業(yè)軟件生成DSM更平滑，邊緣更清晰。

(5) 由圖13、圖14、圖15可以看出，本文算法生成的DSM邊界清晰，更接近于激光掃描數據生成DSM數據。由表7可以看出，本文算法在山墻、坡屋頂及混合建筑物等處生成DSM精度均優(yōu)于Vaihingen數據集匹配生成DSM精度。由圖16、圖17可以看出，采用本文算法生成真正射影像優(yōu)于Vaihingen數據集中匹配結果生成真正射影像。

3 結 論

借鑒Colmap基于隨機物方patch面元掃描優(yōu)化深度圖進而生成DSM的算法，本文提出了基于隨機傳播COLVLL無人機影像DSM生成算法，即通過引入VLL算法，將深度值傳播改進為高程值傳播。利用隨機PatchMatch對影像進行掃描迭代，同時結合VLL算法對隨機生成的高程值進行多次迭代更新得到DSM。由于省略了深度圖等中間過程，因此克服了在弱紋理區(qū)域、高差范圍過大區(qū)域可能產生的誤匹配，以及流程復雜，計算量大，內存占用率高等問題，在最終生成的DSM精度及效率方面較當前國內外主流商業(yè)軟件(SURE，PhotoScan，Pix4D)有了一定的改善。對于大場景超大數據量無人機影像DSM自動生成效率與精度的提升，將是筆者后續(xù)的研究工作。