基于神經(jīng)網(wǎng)絡(luò)密集匹配的資源三號(hào)DSM提取

韓昊，李參海，丘曉楓

(1.遼寧工程技術(shù)大學(xué) 測(cè)繪與地理科學(xué)學(xué)院，遼寧 阜新 123000；2.自然資源部國(guó)土衛(wèi)星遙感應(yīng)用中心，北京 100048)

0 引言

立體匹配是一種從不同視角影像中尋找像點(diǎn)對(duì)應(yīng)關(guān)系的技術(shù)，通過(guò)立體匹配對(duì)多視航空航天遙感圖像進(jìn)行三維場(chǎng)景重建一直是攝影測(cè)量與遙感中的核心問(wèn)題[1]。隨著我國(guó)資源三號(hào)系列立體衛(wèi)星的陸續(xù)成功發(fā)射組網(wǎng)[2]，衛(wèi)星的立體匹配及數(shù)字地表模型生成成為了研究熱點(diǎn)。幾十年來(lái)，立體匹配作為一個(gè)經(jīng)典的研究主題，傳統(tǒng)上被描述為一個(gè)多階段優(yōu)化問(wèn)題，包括匹配代價(jià)計(jì)算、代價(jià)聚合、視差優(yōu)化和后處理。匹配成本計(jì)算是立體匹配的第一步，它為左圖像塊和可能的相應(yīng)右圖像塊提供了初始相似性度量。傳統(tǒng)的立體匹配方法利用像素周圍圖像塊的低級(jí)特征來(lái)測(cè)量差異。通常使用一些常見(jiàn)的局部描述符，例如絕對(duì)差(absolute difference，AD)、CENSUS、BRIEF、歸一化互相關(guān)(normalized cross-correlation，NCC)或它們的組合(例如AD-CENSUS)[3]。

已有的立體匹配算法可以分為局部匹配算法、全局匹配算法、半全局匹配算法以及基于深度學(xué)習(xí)的立體匹配[4]。局部方法運(yùn)行速度非常快，但是質(zhì)量通常較差。全局匹配算法匹配效果較局部方法效果更好，但由于全局匹配算法本身的計(jì)算復(fù)雜度高、時(shí)間長(zhǎng)，因此很難應(yīng)用到大像幅的遙感影像匹配。Hirschmüller[5]提出的半全局匹配算法在匹配效果和計(jì)算效率方面有較好的平衡，而且比較容易實(shí)現(xiàn)并行加速，因此被廣泛應(yīng)用于近景、航空影像匹配中[6]，取得了很好的效果。立體衛(wèi)星影像具有較大幾何和輻射差異，在使用半全局立體匹配(sim-global matching，SGM)算法匹配時(shí)通常會(huì)出現(xiàn)誤匹配和視差空洞。

近些年以神經(jīng)網(wǎng)絡(luò)為首的機(jī)器學(xué)習(xí)算法在多個(gè)領(lǐng)域取得了突破性進(jìn)展，計(jì)算機(jī)視覺(jué)與深度學(xué)習(xí)的結(jié)合給遙感帶來(lái)了新的研究熱點(diǎn)。在近景影像立體匹配方面，基于深度學(xué)習(xí)的立體匹配算法已經(jīng)勝過(guò)Kitti和Middlebury立體數(shù)據(jù)集上的每個(gè)SGM派生算法。通常將基于深度學(xué)習(xí)的立體匹配方法分為三類[7]：非端到端學(xué)習(xí)算法、端到端學(xué)習(xí)算法以及無(wú)監(jiān)督學(xué)習(xí)算法。對(duì)于非端到端立體聲方法，引入了卷積神經(jīng)網(wǎng)絡(luò)(convolutional neural network，CNN)來(lái)替換傳統(tǒng)立體匹配算法中的一個(gè)或多個(gè)流程。Zbontar等[8]首先成功地將匹配代價(jià)替換為深度學(xué)習(xí)的方法，并且在準(zhǔn)確性和速度方面都比傳統(tǒng)方法獲得了更可觀的效果。Shaked等[9]設(shè)計(jì)了一種新的高速網(wǎng)絡(luò)架構(gòu)，用于計(jì)算每種可能差異下的匹配成本。Seki等[10]提出了一種SGM-Net來(lái)計(jì)算SGM的代價(jià)。通過(guò)精心設(shè)計(jì)和監(jiān)督神經(jīng)網(wǎng)絡(luò)，還可以通過(guò)端到端的深度學(xué)習(xí)方法獲得精細(xì)的視差，無(wú)需進(jìn)行后處理。隨著Mayer等[11]的成功，端到端立體匹配網(wǎng)絡(luò)在立體匹配算法中變得越來(lái)越流行，端到端的網(wǎng)絡(luò)通常網(wǎng)絡(luò)層數(shù)多、訓(xùn)練時(shí)間長(zhǎng)、計(jì)算時(shí)間較慢，不適合數(shù)據(jù)量較大的衛(wèi)星遙感影像。在過(guò)去的幾年中，基于空間變換和視圖合成，已經(jīng)提出了幾種用于立體匹配的無(wú)監(jiān)督學(xué)習(xí)方法。Garg等[12]使用從泰勒展開(kāi)中得到的不完全可區(qū)分的圖像重建損失訓(xùn)練神經(jīng)網(wǎng)絡(luò)進(jìn)行深度估計(jì)。但迄今為止，無(wú)監(jiān)督的深度解決方案雖然產(chǎn)生了令人鼓舞的初步結(jié)果，但是目前在實(shí)際使用中仍不能獲得可靠的信息。

目前深度學(xué)習(xí)的立體匹配方法已經(jīng)被標(biāo)準(zhǔn)測(cè)試集成功驗(yàn)證，并廣泛應(yīng)用于車載立體匹配[13]以及無(wú)人機(jī)避障[14]，在遙感影像的立體影像處理尚不成熟。劉瑾等[15]研究了深度學(xué)習(xí)的方法在航空遙感影像密集匹配上的性能。本文將基于孿生網(wǎng)絡(luò)的立體匹配方案用于國(guó)產(chǎn)線陣衛(wèi)星影像，并與基于半全局匹配(sim-global block matching，SGBM)生成數(shù)字地表模型(digital surface model，DSM)的算法效果進(jìn)行對(duì)比，證明了本文算法提取效果較好，提取結(jié)果能夠達(dá)到與商業(yè)軟件相近。

1 基于孿生網(wǎng)絡(luò)密集匹配的DSM提取流程

提取流程共分為三大部分：預(yù)處理部分、密集匹配部分以及數(shù)字地表模型提取部分。

1.1 影像預(yù)處理

預(yù)處理分別應(yīng)用于每個(gè)經(jīng)正射糾正后的立體對(duì)，共分為四個(gè)步驟，如圖1所示。其目的是產(chǎn)生兩個(gè)主要輸出：①立體校正網(wǎng)格，該網(wǎng)格可將一對(duì)圖像重新采樣為完美對(duì)齊的幾何形狀；②視差范圍的估計(jì)值。

步驟1：此步驟包括根據(jù)傳感器模型RPC(rational polynomial coefficients)和初始低分辨率數(shù)字地表模型(例如SRTM)估算立體校正網(wǎng)格。這些網(wǎng)格用于對(duì)立體對(duì)中的圖塊逐塊重采樣，本文使用一種基于核線約束的幾何形狀的近似迭代[16]，這種方法遞歸地估計(jì)兩個(gè)糾正網(wǎng)格。這些網(wǎng)格可與連接點(diǎn)一起將輸入圖像重新排列為核線影像。

1.2 孿生網(wǎng)絡(luò)密集匹配生成視差圖

在進(jìn)行影像預(yù)處理后，計(jì)算左右影像視差圖。通過(guò)立體匹配計(jì)算視差圖，然后使用傳感器模型對(duì)視差進(jìn)行空三計(jì)算以獲得3D點(diǎn)。將匹配所得的目標(biāo)地形區(qū)域(默認(rèn)為成對(duì)覆蓋的最大區(qū)域面積)依據(jù)視差圖劃分為地形圖塊，并形成遮擋掩碼，重疊地形圖塊僅計(jì)算一次。

密集匹配的流程可分為代價(jià)計(jì)算、代價(jià)聚合、視差計(jì)算、視差精化4個(gè)步驟。本文算法將孿生神經(jīng)網(wǎng)絡(luò)MC-CNN(matching-cost CNN)計(jì)算匹配代價(jià)引入資源三號(hào)密集匹配生成數(shù)字地表模型流程中。MC-CNN包括兩種結(jié)構(gòu)的網(wǎng)絡(luò)：Fast和Slow結(jié)構(gòu)，前者的網(wǎng)絡(luò)結(jié)構(gòu)處理速度更快，但所生成視差值精度稍遜于后者。由于Slow網(wǎng)絡(luò)對(duì)訓(xùn)練數(shù)據(jù)量和計(jì)算機(jī)內(nèi)存均有較高要求，本文采用Fast網(wǎng)絡(luò)作為實(shí)驗(yàn)網(wǎng)絡(luò)。網(wǎng)絡(luò)框架如圖2所示。

1)代價(jià)計(jì)算。這里引入孿生網(wǎng)絡(luò)模型用來(lái)學(xué)習(xí)相似性度量(similarity measure)，即進(jìn)行代價(jià)計(jì)算。輸入為左右一組圖塊，輸出是它們之間的相似性。兩側(cè)卷積層用于分別提取左右圖塊高維特征，兩側(cè)網(wǎng)絡(luò)模型共享權(quán)值。模型結(jié)構(gòu)末尾設(shè)置了一個(gè)特征提取器及歸一化層，目的是將圖像塊表示成特征向量的形式，然后對(duì)這兩個(gè)特征向量進(jìn)行比較。決策網(wǎng)絡(luò)以點(diǎn)積作為相似性度量，對(duì)提取的共有特征計(jì)算點(diǎn)積并輸出結(jié)果。根據(jù)點(diǎn)積的大小反映左右兩樣本之間的相似性。

網(wǎng)絡(luò)的訓(xùn)練通過(guò)輸入相同圖像位置為中心的樣本對(duì)來(lái)計(jì)算損失，其中一個(gè)樣本屬于正類別，一個(gè)樣本屬于負(fù)類別。通過(guò)使鉸鏈損失函數(shù)最小化，進(jìn)行網(wǎng)絡(luò)訓(xùn)練。鉸鏈損失函數(shù)定義為：max(0，m+s--s+)。s-為負(fù)樣本的輸出；s+為正樣本的輸出；余量m為正實(shí)數(shù)。也就是說(shuō)，正樣本要比負(fù)樣本至少大m，loss的值才為0。通過(guò)最小化鉸鏈損失來(lái)訓(xùn)練網(wǎng)絡(luò)，實(shí)驗(yàn)中將m設(shè)置為0.2。

網(wǎng)絡(luò)最終輸出的匹配代價(jià)由M來(lái)表示，代價(jià)M的計(jì)算方法如式(1)所示。

(1)

2)代價(jià)聚合。在計(jì)算匹配代價(jià)后需進(jìn)行代價(jià)聚合，采用可變十字交叉代價(jià)聚合(cross-based cost aggregation，CBCA)[18]方式。CBCA 的目標(biāo)是找到像素p周圍顏色相似的像素，并按照一定的規(guī)則將它們的代價(jià)值聚合成像素p的最終代價(jià)值，聚合窗口如圖3所示。圖3中像素p的支持區(qū)域是合并其垂直臂上的所有像素的水平臂，pl為左側(cè)水平聚合臂，q是p的垂直臂上的某個(gè)像素，p的聚合區(qū)域就是所有q(包括p自身)的水平臂的并集。

3)視差計(jì)算及精化。視差圖D通過(guò)贏家通吃策略(winner takes all，WTA)進(jìn)行計(jì)算，即尋找使M(p，d)最小的視差d(式(2))。

D(p1→2(xe，ye))=argmindM(p，d)

(2)

式中：D(p1→2(xe，ye))為視差圖；p為像點(diǎn)；d為視差值。MC-CNN在視差優(yōu)化方法上采用了與SGM相同的策略，匹配后經(jīng)過(guò)亞像素增強(qiáng)、中值濾波以及雙邊濾波進(jìn)行處理。本文采用窗口為5×5的中值濾波器。

1.3 生成數(shù)字地表模型

1)視差圖轉(zhuǎn)換為3D點(diǎn)云。從視差圖p1→2，使用g1和g2獲得傳感器像空間中的同源點(diǎn)集H，計(jì)算如式(3)所示。

H(x，y)=(g1(xe，ye)，g2(xe+d1→2(xe，ye)，ye))

(3)

根據(jù)這些點(diǎn)，使用前方交會(huì)計(jì)算出對(duì)應(yīng)視線，并將這些線與3D點(diǎn)位置相交。出于數(shù)值精度的考慮，此計(jì)算先在ECEF坐標(biāo)中完成，后轉(zhuǎn)換為WGS84坐標(biāo)[19]。通過(guò)將此交點(diǎn)應(yīng)用于所有視差圖中已計(jì)算的視差，形成一個(gè)3D點(diǎn)云P=(x，y，h)k。

2)柵格化。最終的數(shù)字地表模型是通過(guò)點(diǎn)云的柵格化生成的，即將上文計(jì)算的3D點(diǎn)云轉(zhuǎn)換為地理定位的柵格圖像。對(duì)于地形網(wǎng)格的每個(gè)像元，根據(jù)所定義的分辨率形成像元中心，對(duì)于像元中心(cx，cy)，尋找一組與其歐氏距離小于k·res的點(diǎn)C，計(jì)算如式(4)、式(5)所示。

dc(x，y)=‖(x，y)-(cx，cy)‖

(4)

C=(x，y，h)∈P，dc(x，y)

(5)

式中：dc為視差；C為點(diǎn)集；res代表分辨率；k為用戶設(shè)定的參數(shù)；h為中心點(diǎn)的高程值，通過(guò)在一定范圍內(nèi)對(duì)3D點(diǎn)進(jìn)行高斯加權(quán)來(lái)獲得。高斯加權(quán)方法如式(6)所示，并最終生成數(shù)字地表模型柵格圖。

(6)

式中：h為高程值；C為點(diǎn)集；e為常數(shù);σ為標(biāo)準(zhǔn)差。

2 實(shí)驗(yàn)分析

2.1 孿生網(wǎng)絡(luò)模型及參數(shù)設(shè)置

本文模型訓(xùn)練采用了佐治亞理工學(xué)院[20]提供的5×5窗口孿生神經(jīng)網(wǎng)絡(luò)預(yù)訓(xùn)練模型，并在預(yù)訓(xùn)練模型基礎(chǔ)上進(jìn)行精細(xì)化二次訓(xùn)練，以節(jié)約訓(xùn)練時(shí)間。分支網(wǎng)絡(luò)超參數(shù)如表1所示。

表1 分支網(wǎng)絡(luò)超參數(shù)

2.2 訓(xùn)練數(shù)據(jù)與平臺(tái)

在采用預(yù)訓(xùn)練模型進(jìn)行遷移學(xué)習(xí)的基礎(chǔ)上，添加80%WHU MVS(wuhan university multi view stereo)航空影像立體數(shù)據(jù)集[21]與20%Middlebury立體數(shù)據(jù)集進(jìn)行模型精細(xì)化訓(xùn)練，共選取325組影像對(duì)。采用混合數(shù)據(jù)集是為了增加模型的泛化能力。本文從325訓(xùn)練圖像對(duì)中提取了12 800個(gè)示例，一半標(biāo)為正類，一半標(biāo)為負(fù)類，并將圖像轉(zhuǎn)化為灰度圖像。通過(guò)減去平均值并除以像素灰度值的標(biāo)準(zhǔn)偏差來(lái)預(yù)處理每張圖像。訓(xùn)練時(shí)網(wǎng)絡(luò)每批的輸入是128對(duì)圖像塊，在預(yù)測(cè)時(shí)，網(wǎng)絡(luò)兩端輸入一對(duì)核線立體像對(duì)輸出相似性分?jǐn)?shù)。

1) WHU MVS數(shù)據(jù)集。WHU MVS航拍影像數(shù)據(jù)集涵蓋了貴州省梅丹縣約6.7 km×2.2 km的地區(qū)，距離地面 550 m 處拍攝地面分辨率為0.1 m，大小為 5 376像素×5 376像素的1 776張圖像，并有相對(duì)應(yīng)的 1 776 個(gè)深度圖作為地面實(shí)況，且提供了總共1 760張沿飛行方向的視差圖。該數(shù)據(jù)集中涵蓋了六個(gè)代表性的子區(qū)域，這些區(qū)域涵蓋了不同的立體匹配場(chǎng)景類型，分別為工廠與郊區(qū)、樹(shù)木與道路區(qū)、居民區(qū)、復(fù)雜屋頂區(qū)、城鎮(zhèn)中心區(qū)、農(nóng)業(yè)用地和山區(qū)，可直接作為深度學(xué)習(xí)方法的訓(xùn)練和測(cè)試集。

2) Middlebury數(shù)據(jù)集。Middlebury 立體數(shù)據(jù)集中的影像對(duì)是由Middlebury學(xué)院拍攝的多組室內(nèi)場(chǎng)景。受控的照明條件使得其比 KITTI 數(shù)據(jù)集和Driving數(shù)據(jù)集具有更高的密度和視差精度。Middlebury 數(shù)據(jù)集中的圖像顯示了具有不同復(fù)雜度的靜態(tài)室內(nèi)場(chǎng)景，包括重復(fù)結(jié)構(gòu)、遮擋、線狀物體以及無(wú)紋理區(qū)域。Middlebury2003、2005、2006與2014 數(shù)據(jù)集共包含 65 個(gè)高分辨率立體影像對(duì)，并使用結(jié)構(gòu)光掃描儀獲得的相機(jī)校準(zhǔn)參數(shù)和真實(shí)視差圖。

3) 訓(xùn)練平臺(tái)。本文孿生網(wǎng)絡(luò)的訓(xùn)練和測(cè)試在Ubuntu系統(tǒng)下NVIDIA Tesla K80顯卡上進(jìn)行，采用Jupyter Notebook IDE、PyTorch1.5框架，精細(xì)化訓(xùn)練共迭代50次。由于本文在預(yù)訓(xùn)練模型基礎(chǔ)上精細(xì)化訓(xùn)練，學(xué)習(xí)速率設(shè)置為0.000 3。

2.3 數(shù)字地表模型提取實(shí)驗(yàn)

資源三號(hào)01、02、03星分別于2012、2016、2020年發(fā)射升空，填補(bǔ)了中國(guó)立體衛(wèi)星領(lǐng)域的空白。本文選取02星于2019年拍攝的河北邯鄲地區(qū)一組質(zhì)量較好的前后視影像進(jìn)行實(shí)驗(yàn)，影像分辨率2.5 m，實(shí)驗(yàn)區(qū)影像如圖4、圖5所示。實(shí)驗(yàn)區(qū)分為復(fù)雜山區(qū)和多地形影像兩組，多地形影像包含了山地、建筑物、河流、橋梁、平地等多種地貌。

數(shù)字地表模型生成及評(píng)價(jià)實(shí)驗(yàn)均在Ubuntu系統(tǒng)下，Intel Core i7-660U 2.60 GHz CPU上進(jìn)行，使用Python編程語(yǔ)言進(jìn)行程序編寫(xiě)。

2.4 實(shí)驗(yàn)結(jié)果分析

1) DSM定性結(jié)果分析。首先對(duì)所生成視差圖進(jìn)行定性分析，分別對(duì)影像中的典型復(fù)雜山地區(qū)域、丘陵區(qū)域、建筑物區(qū)域、平坦地區(qū)生成的未經(jīng)填充且未添加DEM輔助[22]的5 m分辨率數(shù)字地表模型效果進(jìn)行對(duì)比。從圖6、圖7、圖8中可以看出，本文方法通過(guò)卷積操作深度提取了匹配塊的特征，相較SGBM方法在實(shí)驗(yàn)區(qū)的山區(qū)能夠通過(guò)密集匹配生成連續(xù)視差，具有更好的DSM提取效果，山體部分匹配效果更加理想，SGBM算法在高分辨率數(shù)字地表模型提取時(shí)存在大量的誤匹配和視差空洞。但兩種方法均未在連續(xù)陰影及視差較大的區(qū)域取得良好效果。本文在低矮建筑物、橋梁匹配相較SGBM算法細(xì)節(jié)也更為突出，在平原地區(qū)的低矮建筑物DSM提取邊緣更加清晰，密集建筑物誤匹配現(xiàn)象更少。

2)DSM精度分析。在精度分析方面，針對(duì)每種地表類型選取適量可靠檢查點(diǎn)，將本文方法生成的DSM中無(wú)視差空洞部分與PCI Geomatica 2018軟件生成的5 m分辨率DSM進(jìn)行精度對(duì)比。PCI Geomatica是世界級(jí)的專業(yè)地理信息服務(wù)商PCI公司的旗艦產(chǎn)品，是目前最優(yōu)秀的衛(wèi)星數(shù)字地表模型生產(chǎn)軟件之一，其生成的DSM精度與效果通常優(yōu)于MicMac、ASP等同類軟件[23]。PCI生成的DSM如圖9所示。圖9中可明顯看出，在視差空洞部分，即圖中圓圈圈出部分，軟件對(duì)所生成的DSM進(jìn)行了填充或模糊處理，且PCI的DSM生成添加了低分辨率DEM作為輔助填充。精度對(duì)比由平均誤差、均方根誤差兩個(gè)指標(biāo)對(duì)比完成，結(jié)果如表2所示。

表2 DSM精度對(duì)比 m

從表2中數(shù)據(jù)可看出，本文方法生成的5 m分辨率DSM在各種地物類別中的誤差均值接近于0，均方根誤差較小，其中山地誤差較大，這主要是由于山地通常匹配難度高、視差大、遮擋較多。平均誤差數(shù)值在0±1.5 m內(nèi)，證明與PCI所生成DSM差別很小，驗(yàn)證了本文方法生成DSM的可靠性。

3 結(jié)束語(yǔ)

本文提出了一種基于孿生神經(jīng)網(wǎng)絡(luò)的資源三號(hào)衛(wèi)星影像立體匹配方法，通過(guò)孿生神經(jīng)網(wǎng)絡(luò)密集匹配實(shí)現(xiàn)了國(guó)產(chǎn)資源三號(hào)系列02星影像的立體匹配，并實(shí)現(xiàn)了高精度、高分辨率的DSM提取。運(yùn)用孿生網(wǎng)絡(luò)提取匹配塊的深度特征，并計(jì)算匹配測(cè)度，提高了資源三號(hào)影像的密集匹配效果。在具有復(fù)雜地貌的立體影像上進(jìn)行實(shí)驗(yàn)，與傳統(tǒng)方法做了比較，驗(yàn)證了本文方法在衛(wèi)星遙感影像立體匹配生成DSM方面效果與商業(yè)軟件效果相近，且匹配效果優(yōu)于基于半全局匹配的方法。本文的研究為開(kāi)展1∶2.5萬(wàn)以及更大比例尺DSM的生產(chǎn)提供了新的技術(shù)思路。