韓昊，李參海，丘曉楓

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

0 引言

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

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

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

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

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

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

1.1 影像預處理

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

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

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

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

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

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

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

網(wǎng)絡最終輸出的匹配代價由M來表示，代價M的計算方法如式(1)所示。

(1)

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

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

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

(2)

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

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

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

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

(3)

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

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

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

(4)

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

(5)

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

(6)

式中：h為高程值；C為點集；e為常數(shù);σ為標準差。

2 實驗分析

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

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

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

2.2 訓練數(shù)據(jù)與平臺

在采用預訓練模型進行遷移學習的基礎(chǔ)上，添加80%WHU MVS(wuhan university multi view stereo)航空影像立體數(shù)據(jù)集[21]與20%Middlebury立體數(shù)據(jù)集進行模型精細化訓練，共選取325組影像對。采用混合數(shù)據(jù)集是為了增加模型的泛化能力。本文從325訓練圖像對中提取了12 800個示例，一半標為正類，一半標為負類，并將圖像轉(zhuǎn)化為灰度圖像。通過減去平均值并除以像素灰度值的標準偏差來預處理每張圖像。訓練時網(wǎng)絡每批的輸入是128對圖像塊，在預測時，網(wǎng)絡兩端輸入一對核線立體像對輸出相似性分數(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張圖像，并有相對應的 1 776 個深度圖作為地面實況，且提供了總共1 760張沿飛行方向的視差圖。該數(shù)據(jù)集中涵蓋了六個代表性的子區(qū)域，這些區(qū)域涵蓋了不同的立體匹配場景類型，分別為工廠與郊區(qū)、樹木與道路區(qū)、居民區(qū)、復雜屋頂區(qū)、城鎮(zhèn)中心區(qū)、農(nóng)業(yè)用地和山區(qū)，可直接作為深度學習方法的訓練和測試集。

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

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

2.3 數(shù)字地表模型提取實驗

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

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

2.4 實驗結(jié)果分析

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

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

表2 DSM精度對比 m

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

3 結(jié)束語

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