高分辨率遙感影像DSM的改進半全局匹配生成方法

楊幸彬，呂京國,2，江 珊，張丹璐

1. 北京建筑大學測繪與城市空間信息學院，北京 100044； 2. 城市空間信息工程北京市重點實驗室，北京 100044

從衛(wèi)星立體影像獲取精確數(shù)字表面模型的關鍵在于影像間的密集匹配，通常按照匹配優(yōu)化策略的不同主要分為局部和全局匹配算法兩種[1]。局部匹配算法通過計算像素的最小匹配代價獲取最佳匹配點，其特點是速度快，但在重復紋理或低紋理區(qū)域的匹配效果較差，因此不能滿足生成高精度數(shù)字表面模型的需求。全局匹配算法[2-3]匹配效果較好，但由于算法本身的計算復雜度較高、計算時間較長，因此很難應用到大像幅的遙感影像匹配。文獻[4—5]提出的半全局匹配算法在匹配效果和計算效率方面有較好的平衡，而且比較容易實現(xiàn)并行加速，因此被廣泛應用于近景[6]、航空[7]和衛(wèi)星[8-10]影像匹配中。

目前，將半全局匹配算法應用于衛(wèi)星影像匹配的思路主要有兩種：基于物方的OSGM匹配(object-based semi-global matching)和基于像方的SGM匹配。文獻[11—12]提出的OSGM匹配算法在物方空間上計算代價矩陣，然后進行代價積聚和最優(yōu)視差選擇得到最終的匹配結果，這種方法省去了核線影像生成的步驟，可以直接應用于多視影像的匹配中。基于像方的SGM匹配算法[13]在影像面上計算代價矩陣，屬于嚴格的逐像素匹配算法。文獻[10]利用衛(wèi)星影像數(shù)據(jù)對OSGM和SGM算法的匹配效果進行了對比分析，結果表明基于像方SGM算法獲取的DSM精度要高于OSGM，且SGM算法可以利用左右視差一致性檢查剔除誤匹配，因此生成的DSM粗差較少，此外，SGM算法的懲罰系數(shù)P2可以根據(jù)影像灰度信息自適應調(diào)整，也有效地提高了DSM精度。

由于衛(wèi)星影像的幅寬較大，直接應用SGM算法存在一定的困難。文獻[14—15]提出利用影像分塊約束視差搜索范圍，基于SGM算法逐塊獲取匹配結果；文獻[16]將已有的DEM高程信息反算至核線影像，限制SGM算法的視差空間大小。這兩種方法在一定程度上解決了大像幅影像匹配時的難點問題，然而算法本身的計算效率并不高，且每個像素都需要在相同的視差搜索范圍內(nèi)計算匹配代價，造成了誤匹配概率的增大。文獻[7,17]提出的tSGM算法在金字塔影像上逐層進行半全局匹配，利用上層匹配結果動態(tài)限定每個像素的視差搜索范圍，大大地減少了算法的內(nèi)存開銷和計算時間，同時提高了匹配準確度。文獻[18]將tSGM算法應用于衛(wèi)星立體影像匹配取得了較好的結果，但當影像重復紋理信息較多時，tSGM算法在核線范圍內(nèi)獲取的頂層匹配結果會受到一定的影響，且會逐層傳遞到下一層匹配結果上。另外，tSGM算法通過搜索窗口內(nèi)最大、最小視差來確定視差搜索范圍，視差圖邊緣像素的搜索范圍可能超出視差預設最大值，導致計算內(nèi)存增加，視差邊緣處誤匹配的概率也會增大。在視差后處理階段，tSGM算法采用傳統(tǒng)的中值濾波算法不能將噪聲同細節(jié)區(qū)別出來，因此視差圖的邊緣和尖角會受到破壞，導致圖像細節(jié)模糊。文獻[19]利用改進的導向中值濾波替代傳統(tǒng)的中值濾波算法有效保護了視差圖的邊緣信息，其不足之處在于濾波過程中閾值的設置不具備較好的普適性。

針對上述問題，本文提出一種基于改進半全局匹配的高分辨率遙感影像DSM生成方法。該方法利用影像分塊策略約束半全局匹配時的初始視差搜索范圍，對核線影像塊建立金字塔后逐層進行半全局匹配，并引入顧及影像紋理信息的視差圖膨脹、腐蝕算法和加權中值濾波算法保護視差圖邊緣信息，同時減少匹配過程中的內(nèi)存開銷和計算時間。利用WorldView3和資源三號立體影像匹配結果驗證了本文方法的有效性。

1 改進半全局匹配的高分辨率遙感影像DSM提取流程

本文基于改進半全局匹配的高分遙感影像DSM生成方法主要包括：影像精確定向、核線影像生成、影像密集匹配和前方交會生成DSM 4個部分，流程如圖1所示。

1.1 連接點幾何約束下的影像精確定向

高分辨率遙感衛(wèi)星影像通常采用有理函數(shù)模型來描述像點與地面點之間的對應關系，由于衛(wèi)星影像的RPC模型存在一定的系統(tǒng)誤差，需要對有理函數(shù)模型進行系統(tǒng)誤差補償。文獻[20]采用像方系統(tǒng)誤差補償方法獲得了較好的定向結果，其實質(zhì)是通過對原始的有理函數(shù)模型增加補償項，消除像點坐標上的系統(tǒng)誤差，從而提高基于有理函數(shù)模型的目標定位精度，其基本原理為

(1)

式中，(Δr,Δc)為像點坐標(r,c)的改正值；(Xn,Yn,Zn)為像點對應的地面點坐標；pi(i=1,2,3,4)為關于Xn、Yn、Zn的三階有理多項式。本文首先獲取影像間匹配點，得到立體影像間相對幾何約束關系后對定向模型進行系統(tǒng)誤差補償，幾何約束關系定義為影像面上的仿射變換

(2)

式中，(e0,e1,e2,f0,f1,f2)是仿射變換系統(tǒng)誤差補償參數(shù)。在式(1)和式(2)的基礎上，采用文獻[21]中的方法對每個連接點建立誤差方程，然后采用直接列改化法方程的策略消除一類連接點地面坐標未知數(shù)，再基于傳統(tǒng)光束法平差求解即可得到精確的定向模型。

1.2 核線影像生成

為了提高影像匹配的精度并減少匹配時間，在立體影像匹配前需要對影像進行核線校正，將二維匹配搜索降低為一維搜索。通常情況下，高分遙感影像的核線類似于雙曲線，但在局部范圍內(nèi)可以近似看作直線[22]。鑒于此，本文利用立體像對的連接點計算得到影像間的仿射變換模型，在模型的基礎上對原始影像進行分塊處理，對于左影像塊利用仿射變換關系可找到相應的右影像塊，分塊完成后利用文獻[23]中的投影軌跡法生成核線影像對，以保證相鄰核線影像塊生成的視差圖之間不產(chǎn)生空隙，塊與塊之間保證一定的重疊。通過對影像進行分塊，一方面可以減小核線影像的誤差，另一方面，通過將視差搜索限定在一個較小的區(qū)間，可以減小重復紋理區(qū)域因視差搜索范圍過大造成的誤匹配。

1.3 影像密集匹配

原始的SGM算法[4]通過建立全局能量函數(shù)，使能量函數(shù)達到最小時獲取潛在的最優(yōu)視差。能量函數(shù)為

(3)

全局能量函數(shù)ED由數(shù)據(jù)項和平滑項組成。式(3)右側第1項為數(shù)據(jù)項，表示所有像素點的匹配代價之和，Cp,Dp表示左影像像素p視差為Dp時與右影像點的匹配代價[24]。第2項為平滑項，其目的是對像點p臨近視差差值等于1的像素q增加懲罰數(shù)P1。這使得相鄰像素點的視差盡可能的一致，其中T[·]為一個判斷函數(shù)，當判斷條件為真時返回1，否則為0。第3項為平滑項，對像素點p臨近視差差值大于1的像素點q給予更大的懲罰數(shù)P2。通常P2隨影像梯度的變化而改變，P2不能過大，否則視差邊緣會被過度平滑。

SGM算法首先計算左影像上每個像素在視差搜索范圍內(nèi)的匹配代價，組成代價矩陣。代價矩陣維度為r×c×(dmax-dmin)，其中r為影像行數(shù)，c為影像列數(shù)，dmax、dmin表示視差搜索的上下界。計算完代價矩陣后，對代價矩陣上的每個點考慮從多個方向(8或者16方向)進行積聚，如圖2所示。

圖2 視差空間上代價積聚過程Fig.2 Aggregation of costs in disparity space

對像素點p，其在r方向上的積聚函數(shù)為

Lrp,d=Cp,d+

(4)

式中，右側第1項表示像素點p視差為d時的匹配代價；第2項表示像素p在r方向上的前一個像素代價積聚的最小值(包含懲罰系數(shù))；最后一項減去積聚方向上的最小代價值保證了Lr(p,d)

(5)

完成代價矩陣的積聚后，在視差搜索范圍內(nèi)利用勝者為王(winner take all，WTA)算法獲取積聚值最小位置即為當前像素的最優(yōu)視差，最優(yōu)視差選擇的公式見式(6)，D為最終的視差圖

D=mindSp,d

(6)

1.4 改進的半全局匹配

對于大像幅的高分遙感影像而言，將整條核線作為視差搜索范圍時應用半全局匹配算法存在較大的困難。一方面代價矩陣計算和代價積聚過程所需的內(nèi)存開銷較大，計算時間也相對較長；另一方面，由于每個像素都在整條核線上計算匹配代價，當影像中存在較多低紋理或重復紋理信息時，會造成誤匹配概率的增大。

文獻[7]提出tSGM算法，通過在金字塔影像上逐層進行半全局匹配，利用本層匹配結果限制下一層匹配時的視差搜索邊界，從而動態(tài)限定每個像素的視差搜索范圍。tSGM算法確定像素的視差搜索范圍時，首先在視差像素D(x)周圍窗口中搜索視差最大值dmax、最小值dmin，然后判斷最大、最小視差差值是否超出最大預設搜索范圍R。如果dmax-dmin

(7)

確定了本層影像的視差搜索邊界后，將視差上下邊界放大一倍作為下一層核線影像匹配時的視差搜索范圍，然后逐層進行半全局匹配得到最終的匹配結果。

這種搜索窗口內(nèi)視差最大、最小值直接確定視差邊界的方法存在一定的不足。首先，在視差變化較大的區(qū)域(地物邊緣處)，由于像素周圍窗口中最大、最小的視差值存在較大的變化，因此會造成視差邊緣處的像素誤匹配概率的增大；其次，視差邊緣處的搜索范圍較大時，會影響代價矩陣整體的內(nèi)存開銷，導致計算時間增加。

本文提出一種顧及影像紋理信息的視差圖膨脹、腐蝕算法約束視差搜索范圍。其基本原理是在獲取窗口中視差最大、最小值的同時，記錄最大、最小視差值對應核線影像上的像素灰度，依據(jù)灰度變化量與視差值的變化量存在正相關的關系，將視差上下邊界搜索范圍定義為

(8)

式中，dmax、dmin分別為視差像素D(x)膨脹、腐蝕后的視差值，其計算公式為

(9)

式中，ρ表示弱化因子(本文取0.2)，其作用是減弱灰度差對搜索邊界的影響；t_diff(x,xmax)、t_diff(x,xmin)分別為當前像素x與最大、最小視差值對應像素灰度差的絕對值，計算公式為

(10)

式中，D為視差圖像素值；s為結構元素；D?s為視差圖腐蝕操作；D⊕s為視差圖膨脹操作；dmax(i,j)、dmin(i,j)分別為膨脹、腐蝕后視差圖上索引為(i,j)的視差值；w為結構元素對應的窗口大小。

2）水域、草地和林地都有不同程度的減少，林地主要轉(zhuǎn)換為草地（約105.8 km2），水域和草地主要轉(zhuǎn)換為耕地（分別約為136.98 km2、112.19 km2）；

將紋理信息納入視差搜索范圍的判定中，可以進一步縮小視差搜索范圍，尤其是在視差變化較大的區(qū)域，能夠?qū)⒁暡钏阉鞣秶薅ㄔ谝粋€較為合理的區(qū)間內(nèi)，修改后的視差搜索范圍定義為

(11)

式中，Tmax(x)、Tmin(x)分別為視差上下邊界。需要說明的是，當搜索窗口內(nèi)存在多個最大或者最小的視差像素時，此時計算最大的視差tmax作為視差上邊界搜索范圍，最小的視差tmin作為視差下邊界搜索范圍。

圖3為tSGM和顧及影像紋理信息的視差圖膨脹腐蝕算法獲取的視差上下邊界。可以看出，在視差變化較大的像素處，顧及紋理信息的膨脹、腐蝕算法得到的視差上下搜索邊界范圍更小，最大范圍小于20個像素(見圖3中虛線框1)，而tSGM算法的最大視差搜索范圍基本達到了60個像素(見圖3中虛線框2)，說明顧及影像紋理信息的視差圖膨脹腐蝕算法能夠更好地確定像素的視差搜索范圍。

圖3 tSGM和本文方法的視差搜索范圍對比Fig.3 Contrast of tSGM and the proposed method of parallax search

結合顧及影像紋理信息的視差圖膨脹腐蝕算法，本文提出一種改進的半全局算法。首先利用影像分塊策略約束初始視差搜索范圍，減小影像重復紋理信息造成的誤匹配，然后對核線影像塊建立金字塔，在金字塔影像上逐層進行半全局匹配，匹配過程中引入顧及影像紋理信息的視差圖膨脹、腐蝕算法約束每一層影像的視差搜索范圍，逐層得到最終匹配結果，改進的半全局匹配算法流程如圖4所示。

在頂層金字塔核線影像上進行半全局匹配時，初始視差搜索范圍根據(jù)影像縮放比例進行設置(例如原始視差搜索范圍為128，對影像塊構建兩層金字塔后初始視差搜索范圍為32)，設置了初始視差搜索范圍后，采用9×7窗口的Census變換[25]作為代價函數(shù)計算代價矩陣，然后從16個方向?qū)Υ鷥r矩陣進行代價積聚，最后選取像素的最優(yōu)視差形成初始視差圖。

為了確定layer-1層影像的視差搜索范圍，首先利用顧及影像紋理信息的視差圖膨脹、腐蝕算法對layer層視差圖進行膨脹操作得到視差上界搜索范圍tmax，腐蝕后得到視差下界搜索范圍tmin，膨脹腐蝕的窗口大小為7×7，然后對視差搜索邊界放大一倍后作為layer-1層影像的視差搜索上下界。例如對layer層視差值為d的像素p，其在第layer-1層影像上的視差搜索范圍為[Tmin(p),Tmax(p)]，其中Tmin(p)=2(d-tmin(p))，Tmax(p)=2(d+tmax(p))。需要說明的是，對于layer層視差圖上的無效像素，將其膨脹后的視差值減去32作為視差搜索下邊界范圍。

圖4 改進的半全局匹配算法流程Fig.4 Flowchart of modified semi-global matching

利用上述方法確定layer-1層的視差搜索邊界后，在視差搜索范圍內(nèi)計算動態(tài)代價矩陣，然后從8個方向?qū)Υ鷥r矩陣進行積聚。積聚過程中由于當前像素和臨近像素的視差搜索區(qū)可能僅有部分重疊或者沒有重疊，可能導致Lr(p-r,d+k)不存在，因此將代價積聚公式修改為

if
d>Tmax(p-r) ordLr(p,d)=C(p,d)+P2
else
Lr(p,d)=C(p,d)+min(Lr(p-r,d),
Lr(p-r,d-1)+P1,Lr(p-r,d+1)+P1,
Lr(p-r,i)+P2)-minkLr(p-r,k)
end

(12)

上述匹配過程得到的每一層視差圖可能含有匹配粗差，因此每一層影像匹配結束時需要交換左右影像獲取視差圖，剔除左右視差差值大于1的像素，最后對視差圖進行濾波后得到最終的視差結果。對視差圖濾波處理時，本文對文獻[19]中的導向中值濾波算法進行改進，不采用固定閾值判斷鄰域范圍內(nèi)有效像素，而是在計算鄰域范圍內(nèi)的視差中值時，將像素的灰度差值信息視為權值(灰度差越大則對應的視差像素權值就越小)，采用加權中值濾波保護視差圖細節(jié)信息。加權中值濾波算法的原理為

dis(r,c)=weighted_meddis(i,j)=

medwij◇dis(i,j)，(i,j)∈W(r,c)

(13)

圖5顯示了改進的半全局匹配算法和tSGM算法獲取的視差圖?？梢钥闯觯啾扔趖SGM算法，本文方法更好地保護了視差圖的邊緣和尖角(如橢圓1、2所示)，且有效像素數(shù)要多于tSGM算法(如橢圓3所示)。

表1為相同硬件測試環(huán)境下3種算法的匹配時間統(tǒng)計結果(3種算法的代價矩陣計算、代價積聚和最優(yōu)視差選擇過程進行了并行加速[26])。從匹配效率上看，本文匹配算法效率要高于tSGM和SGM算法。需要說明的是，試驗時僅對影像建立了兩層金字塔，當金字塔層數(shù)增多時本文方法和tSGM算法計算時間相差不大。

采用改進的半全局匹配算法獲取視差圖后，利用核線影像與原始影像之間的對應關系，逐像素地將視差圖計算出的同名像點反算至原始影像，再基于前方交會原理[27]即可獲取物方密集點云。

圖5 tSGM和本文方法計算的視差圖Fig.5 Disparity maps computed by SGM and the proposed method

表1 計算時間對比

2 試驗與分析

試驗選取的第1組數(shù)據(jù)為阿曼首都馬斯喀特地區(qū)的WorldView 3下視與后視立體像對，影像分辨率為0.4 m，獲取時間為2014年10月，影像像幅大小為12 984×12 928，立體影像的基高比為0.725 8。由于缺少該區(qū)域精確的DSM作為基準參考，因此將本文算法和商業(yè)軟件ENVI5.4的SGM密集匹配模塊以及ERDAS2016的Tridicon SGM模塊進行相對精度和計算效率的對比分析。立體影像和本文方法生成的DSM如圖6所示。

圖6 WorldView 3全色立體像對和DSMFig.6 Panchromatic stereo images of WorldView-3 and DSM

第2組數(shù)據(jù)為法國圣馬克西姆地區(qū)的資源三號前視與后視立體影像，影像分辨率為3.5 m，獲取時間為2014年8月，像幅大小為16 306×16 384。該數(shù)據(jù)配套了同地區(qū)航空影像獲取的高精度DSM，分辨率為10 m，可作為絕對高程參考數(shù)據(jù)評價本文方法生成的DSM高程精度。資源三號立體影像和本文方法生成的DSM如圖7所示(本文方法生成的DSM在水域和影像無效值區(qū)域進行了掩膜處理)。

基于Win10 64位系統(tǒng)采用VC++2013開發(fā)了密集匹配算法，并使用OpenMP對算法的金字塔影像構建、顧及紋理信息的視差圖膨脹、腐蝕算法、代價矩陣計算、代價積聚、最優(yōu)視差選擇、視差一致性檢查和加權中值濾波過程進行了并行加速。試驗使用的硬件平臺為華碩PU403UF筆記本，處理器為Intel i7-6500U、主頻2.5 GBHz，內(nèi)存8 GB。

圖7 ZY-3全色立體像對和DSMFig.7 Panchromatic stereo images of ZY-3 and DSM

2.1 連接點幾何約束的RPC區(qū)域網(wǎng)平差

為驗證本文基于連接點幾何約束的區(qū)域網(wǎng)平差方法的有效性，在WorldView 3立體影像四角選取4個控制點，同時利用SIFT特征提取算子獲取21對同名像點，采用區(qū)域網(wǎng)平差方法求解系統(tǒng)誤差補償參數(shù)，然后利用補償?shù)挠欣砗瘮?shù)模型將同名點交會出的地面點反投影至左右影像，統(tǒng)計像點反投影誤差的最大值、平均值和標準差(表2)。

表2 WorldView 3立體像對重投影誤差

從表2的試驗結果可以看出，經(jīng)過連接點幾何約束的RFM區(qū)域網(wǎng)平差后，像點坐標的反投影誤差在行方向上標準差僅為0.1個像素，列方向上最大為0.18個像素，達到了亞像素級精度，可以滿足生成核線影像的要求。

2.2 核線影像生成

為驗證本文方法生成核線影像的精度，在WorldView 3影像上選取地形高差較大的城區(qū)和高差較小的平地區(qū)域的核線影像，利用SIFT算子提取核線影像上同名像點，并采用最小二乘匹配精化后得到精確的同名點坐標。通過比較同名像點的y視差差值作為誤差結果，對誤差進行統(tǒng)計分析，剔除不符合正態(tài)分布的粗差點后計算中誤差作為核線影像精度評價指標。表3為上下視差的精度評價結果。

表3 核線影像上下視差統(tǒng)計結果

從表3可以看出，本文方法在城區(qū)和平坦地區(qū)生成的核線影像上下視差標準差在0.3個像素左右，可以滿足核線影像的要求。

2.3 DSM生成及精度分析

2.3.1 WorldView 3生成的DSM精度對比分析

本文對核線影像塊構建兩層金字塔后進行密集匹配。為保證ENVI5.4的SGM密集匹配模塊、ERDAS2016的Tridicon SGM模塊和本文匹配方法對CPU具有相同的利用率，將本文密集匹配算法設置為4線程并行加速(Envi和ERDAS的密集匹配模塊默認為4線程并行)，Envi和ERDAS的匹配參數(shù)按默認的設置。最后，將3種方法生成的密集點云統(tǒng)一柵格化為1 m 分辨率的DSM后進行對比分析。圖8、9為平坦地區(qū)和城區(qū)的DSM對比結果。需要說明的是，本文方法生成的DSM在影像水域和無效值區(qū)域進行了掩膜處理。

從圖8平坦地區(qū)的對比結果可以看出，本文方法和ERDAS生成的DSM均能保持完整地表形態(tài)，而Envi生成的DSM在屋頂上存在缺失現(xiàn)象，造成房屋不完整(圖8中橢圓區(qū)域1)。從DSM的平滑性可以看出，ERDAS生成的DSM平滑性最好，本文方法結果介于ERDAS和Envi之間，Envi生成的DSM在屋頂和路面上存在較多的噪聲(圖8中橢圓區(qū)域1、2)。這主要是由于屋頂和路面屬于重復紋理區(qū)域，在這些區(qū)域的像素匹配代價較為相似，匹配錯誤率較高。比較圖9城區(qū)的DSM可以看出，3種方法生成的DSM在水域均存在較多噪聲，原因和屋頂與路面類似。從圖9中建筑物的DSM可以看出，由于采用顧及紋理信息的視差圖膨脹腐蝕算法和加權中值濾波算法，本文方法生成的DSM能夠更好地保持建筑物的邊緣信息(圖9中橢圓區(qū)域1、2、3)。在建筑物遮擋的陰影區(qū)域，本文方法存在匹配缺失現(xiàn)象(圖中橢圓區(qū)域1、2)，這種現(xiàn)象在ENVI和ERDAS生成的DSM上并不明顯。將ENVI和ERDAS生成的DSM與正射影像圖疊加后發(fā)現(xiàn)，DSM上建筑物區(qū)域的邊界范圍擴大，猜測可能是ENVI和ERDAS的密集匹配算法對建筑物陰影區(qū)域的DSM進行了填充處理。

圖8 平坦地區(qū)DSM對比分析Fig.8 Comparison of DSM in planer area

圖9 城區(qū)DSM對比分析Fig.9 Comparison of DSM in city area

為評價本文方法生成的DSM在垂直方向上的匹配效果，選取典型地物的高程斷面進行3種方法的對比分析。圖10為連續(xù)的房屋斷面，可通過4個連續(xù)的弧形屋頂和兩個斜坡屋頂輪廓反映DSM的垂直精度；圖11為平坦路面的斷面，從斷面曲線的趨勢和平滑性可以間接反映DSM的垂直精度；圖12為弧形棚房斷面，可以根據(jù)屋頂斷面曲線輪廓反映DSM垂直精度。

圖10 3種方法獲取的連續(xù)房屋斷面和局部細節(jié)圖Fig.10 The profiles of continuous house section as well as superposition of details derived from three methods

圖11 3種方法獲取的路面斷面和局部細節(jié)圖Fig.11 The profiles of pavement section as well as superposition of details derived from three methods

圖12 3種方法獲取的棚房斷面疊加和局部細節(jié)圖Fig.12 The profiles of shed building section as well as superposition of details derived from three methods

從圖10弧形屋頂和斜坡屋頂?shù)恼w斷面曲線可以看出，本文方法和ERDAS生成的DSM能夠較好地反映屋頂?shù)恼w輪廓，而Envi結果在斜坡屋頂上存在較多的誤匹配導致屋頂形態(tài)丟失，detail 1和detail 2為弧形屋頂和斜坡屋頂?shù)木植糠糯髤^(qū)域。從局部細節(jié)可以看出，ENVI生成的DSM平滑性不佳，且有一定的粗差產(chǎn)生，ERDAS在房屋遮擋的陰影區(qū)域和屋頂邊緣處(圖10中箭頭指向處)效果不好，本文方法在這些區(qū)域均表現(xiàn)出了較好的匹配結果。這主要是由于本文方法采用顧及紋理信息的視差圖膨脹、腐蝕算法更好地確定了這些區(qū)域的視差搜索范圍，使得誤匹配率降低，而且在視差圖后處理時利用加權中值濾波算法，較好地保持了屋頂?shù)倪吘壭畔?。從圖11路面斷面曲線的整體平滑度可以看出，本文方法獲取的DSM斷面曲線平滑度要優(yōu)于Envi，ERDAS的平滑效果最好，但從局部細節(jié)可以看出，ERDAS的斷面曲線整體趨勢存在一定的偏差，可能是由于ERDAS對DSM進行了過渡平滑處理的緣故。比較圖12棚房的斷面曲線可以看出，3種算法均能較好地反映出屋頂?shù)恼w輪廓，從detail 1和detail 2房屋邊緣的局部細節(jié)可以看出，本文方法較好地保護了房屋的邊緣信息，在DSM的整體效果和局部效果上都具有較好的結果。

2.3.2 ZY-3生成的DSM精度分析

為進一步定量評價本文方法生成DSM的絕對精度，選取資源三號影像上平原、丘陵和山地區(qū)域的DSM與參考DSM進行誤差對比分析，絕對值誤差分布見圖13—圖15，圖16為誤差直方圖統(tǒng)計結果。

從圖13—圖15絕對值誤差分布可以看出，本文方法獲取的DSM在大部分區(qū)域內(nèi)誤差較小，絕對值誤差均值最大的為山區(qū)，這主要是由于山體陰影處匹配誤差較大。對圖16中誤差進行統(tǒng)計分析，剔除不符合正態(tài)分布的誤差后，計算一定范圍內(nèi)誤差的均值、絕對值均值和均方根誤差來定量評定本文方法生成的DSM精度，結果如表4所示。從DSM的均方根誤差結果可以看出，本文方法生成的DSM在平原地區(qū)均方根誤差較小，精度較高，在山區(qū)接近于1.5倍GSD；從誤差的平均值可以看出，3個地區(qū)的誤差均值均表現(xiàn)為負值，推測其原因可能是影像定向結果中存在一定的系統(tǒng)誤差，從而傳遞到了最終生成的DSM中。

圖13 平原地區(qū)DSM絕對值誤差分布圖Fig.13 The absolute error map of DSM in plain area

圖15 山地區(qū)域DSM絕對值誤差分布圖Fig.15 The absolute error map of DSM in mountain area

圖16 不同試驗區(qū)域下的高程誤差直方圖Fig.16 Elevation error histogram under different experimental areas

表4 不同試驗區(qū)域下的DSM精度評定結果

2.4 計算效率和內(nèi)存開銷對比分析

為定量評價本文密集匹配算法的性能，在相同的硬件測試環(huán)境下將本文方法和Envi與ERDAS密集匹配方法進行計算效率和內(nèi)存開銷的對比分析，其中計算效率通過計算時間來反映，內(nèi)存開銷通過比較計算過程中內(nèi)存穩(wěn)定時的最大值，結果見表5。

表5 3種算法的計算性能對比結果

由于本文方法和ERDAS的輸入輸出數(shù)據(jù)相同，因此具有較好的對比性，由表5的對比結果可以看出，本文方法的計算效率要高于ERDAS，但計算時的內(nèi)存開銷也相對較多。與Envi對比后發(fā)現(xiàn)，本文方法計算效率要略低于Envi，但在內(nèi)存開銷上要優(yōu)于Envi，由于Envi的密集匹配模塊自帶了區(qū)域網(wǎng)平差功能，因此綜合計算效率要高于本文方法。從整體上看，本文方法在計算效率和內(nèi)存開銷方面有較好的平衡。

3 結 論

本文提出了一種基于改進半全局匹配的高分辨率遙感影像DSM生成方法，針對影像定向模型中存在的系統(tǒng)誤差，利用影像間連接點幾何約束關系對有理函數(shù)模型系統(tǒng)誤差進行了補償，獲得了亞像素級的定位精度；生成核線影像的過程中，通過對原始影像分塊后基于投影軌跡法獲取到了0.3個像素精度的核線影像；在密集匹配階段，為了減小代價空間大小和計算時間，對影像建立金字塔后逐層進行半全局匹配，匹配過程中引入顧及影像紋理信息的視差圖膨脹腐蝕算法約束視差搜索范圍，更好地保留地物的邊緣特性；最后，采用加權中值濾波算法對視差圖進行濾波，基于前方交會原理得到DSM。利用WorldView 3和資源三號立體影像進行了試驗，結果表明：本文算法獲取的DSM無論在高程精度或地物邊緣特性的保留上都有較好的結果，計算效率和內(nèi)存開銷方面也具有較好的平衡。