基于神經網絡密集匹配的資源三號DSM提取

韓昊，李參海，丘曉楓

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

0 引言

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

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

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

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

1 基于孿生網絡密集匹配的DSM提取流程

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

1.1 影像預處理

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

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

1.2 孿生網絡密集匹配生成視差圖

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

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

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

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

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

(1)

2)代價聚合。在計算匹配代價后需進行代價聚合，采用可變十字交叉代價聚合(cross-based cost aggregation，CBCA)[18]方式。CBCA 的目標是找到像素p周圍顏色相似的像素，并按照一定的規則將它們的代價值聚合成像素p的最終代價值，聚合窗口如圖3所示。圖3中像素p的支持區域是合并其垂直臂上的所有像素的水平臂，pl為左側水平聚合臂，q是p的垂直臂上的某個像素，p的聚合區域就是所有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在視差優化方法上采用了與SGM相同的策略，匹配后經過亞像素增強、中值濾波以及雙邊濾波進行處理。本文采用窗口為5×5的中值濾波器。

1.3 生成數字地表模型

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

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

(3)

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

2)柵格化。最終的數字地表模型是通過點云的柵格化生成的，即將上文計算的3D點云轉換為地理定位的柵格圖像。對于地形網格的每個像元，根據所定義的分辨率形成像元中心，對于像元中心(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為用戶設定的參數；h為中心點的高程值，通過在一定范圍內對3D點進行高斯加權來獲得。高斯加權方法如式(6)所示，并最終生成數字地表模型柵格圖。

(6)

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

2 實驗分析

2.1 孿生網絡模型及參數設置

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

表1 分支網絡超參數

2.2 訓練數據與平臺

在采用預訓練模型進行遷移學習的基礎上，添加80%WHU MVS(wuhan university multi view stereo)航空影像立體數據集[21]與20%Middlebury立體數據集進行模型精細化訓練，共選取325組影像對。采用混合數據集是為了增加模型的泛化能力。本文從325訓練圖像對中提取了12 800個示例，一半標為正類，一半標為負類，并將圖像轉化為灰度圖像。通過減去平均值并除以像素灰度值的標準偏差來預處理每張圖像。訓練時網絡每批的輸入是128對圖像塊，在預測時，網絡兩端輸入一對核線立體像對輸出相似性分數。

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

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

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

2.3 數字地表模型提取實驗

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

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

2.4 實驗結果分析

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

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

表2 DSM精度對比 m

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

3 結束語

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