楊幸彬，呂京國(guó),2，江 珊，張丹璐

1. 北京建筑大學(xué)測(cè)繪與城市空間信息學(xué)院，北京 100044； 2. 城市空間信息工程北京市重點(diǎn)實(shí)驗(yàn)室，北京 100044

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

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

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

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

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

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

1.1 連接點(diǎn)幾何約束下的影像精確定向

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

(1)

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

(2)

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

1.2 核線影像生成

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

1.3 影像密集匹配

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

(3)

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

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

圖2 視差空間上代價(jià)積聚過(guò)程Fig.2 Aggregation of costs in disparity space

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

Lrp,d=Cp,d+

(4)

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

(5)

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

D=mindSp,d

(6)

1.4 改進(jìn)的半全局匹配

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

文獻(xiàn)[7]提出tSGM算法，通過(guò)在金字塔影像上逐層進(jìn)行半全局匹配，利用本層匹配結(jié)果限制下一層匹配時(shí)的視差搜索邊界，從而動(dòng)態(tài)限定每個(gè)像素的視差搜索范圍。tSGM算法確定像素的視差搜索范圍時(shí)，首先在視差像素D(x)周?chē)翱谥兴阉饕暡钭畲笾礵max、最小值dmin，然后判斷最大、最小視差差值是否超出最大預(yù)設(shè)搜索范圍R。如果dmax-dmin

(7)

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

這種搜索窗口內(nèi)視差最大、最小值直接確定視差邊界的方法存在一定的不足。首先，在視差變化較大的區(qū)域(地物邊緣處)，由于像素周?chē)翱谥凶畲蟆⒆钚〉囊暡钪荡嬖谳^大的變化，因此會(huì)造成視差邊緣處的像素誤匹配概率的增大；其次，視差邊緣處的搜索范圍較大時(shí)，會(huì)影響代價(jià)矩陣整體的內(nèi)存開(kāi)銷(xiāo)，導(dǎo)致計(jì)算時(shí)間增加。

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

(8)

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

(9)

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

(10)

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

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

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

(11)

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

圖3為tSGM和顧及影像紋理信息的視差圖膨脹腐蝕算法獲取的視差上下邊界?？梢钥闯觯谝暡钭兓^大的像素處，顧及紋理信息的膨脹、腐蝕算法得到的視差上下搜索邊界范圍更小，最大范圍小于20個(gè)像素(見(jiàn)圖3中虛線框1)，而tSGM算法的最大視差搜索范圍基本達(dá)到了60個(gè)像素(見(jiàn)圖3中虛線框2)，說(shuō)明顧及影像紋理信息的視差圖膨脹腐蝕算法能夠更好地確定像素的視差搜索范圍。

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

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

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

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

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

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

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)

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

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

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

(13)

圖5顯示了改進(jìn)的半全局匹配算法和tSGM算法獲取的視差圖?？梢钥闯?，相比于tSGM算法，本文方法更好地保護(hù)了視差圖的邊緣和尖角(如橢圓1、2所示)，且有效像素?cái)?shù)要多于tSGM算法(如橢圓3所示)。

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

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

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

表1 計(jì)算時(shí)間對(duì)比

2 試驗(yàn)與分析

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

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

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

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

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

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

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

表2 WorldView 3立體像對(duì)重投影誤差

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

2.2 核線影像生成

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

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

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

2.3 DSM生成及精度分析

2.3.1 WorldView 3生成的DSM精度對(duì)比分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

表4 不同試驗(yàn)區(qū)域下的DSM精度評(píng)定結(jié)果

2.4 計(jì)算效率和內(nèi)存開(kāi)銷(xiāo)對(duì)比分析

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

表5 3種算法的計(jì)算性能對(duì)比結(jié)果

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

3 結(jié) 論

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