基于城市遙感衛星影像對的立體匹配

趙 杰，陳小梅*，侯瑋旻，韓嘉威

（1.北京理工大學 光電學院，北京100081；2.北京理工大學 光電成像技術與系統教育部重點實驗室，北京100081）

1 引 言

城市三維重建對于城市規劃、智慧城市構建具有非常重要的意義，三維重建的核心技術是立體匹配［1］。研究人員提出了大量的立體匹配算法，其4個步驟主要包括：（1）匹配代價計算；（2）代價聚合；（3）視差計算；（4）視差細化［2］。

在匹配代價計算中，學者們提出了絕對值差之和（Sum of Absolute differences，SAD）、平方差之和（Sum of Squared differences，SSD）、歸一化互相關系數（Normalized Cross Correlation，NCC）、基于互信息的匹配代價計算方法以及基于Census變換的匹配代價計算方法等。其中，基于Census變換的匹配代價計算方法由于它對于光照度變化不敏感以及計算效率高得到了廣泛的使用［3-4］。但是Census算法對于噪聲非常敏感且在弱紋理區域匹配精度不足，城市衛星遙感影像上存在大量的噪聲以及弱紋理區域，因此該算法對于城市區域衛星遙感影像無法取得非常好的效果［4］。

根據代價聚合方法的不同，立體匹配算法可分為局部立體匹配算法與全局立體匹配算法［2］。局部立體匹配算法通過設定一鄰域窗口的形式，利用像素局部的信息來表征當前像素的匹配特征，增加了可靠性與穩定性。全局立體匹配算法在進行像素點匹配時利用了圖像的所有信息，因此在處理具有噪聲、重復紋理以及遮擋現象的圖像時表現更為出色［5］。不過，全局立體匹配算法的計算效率低于局部立體匹配算法［6］。2008年，Hirschmuller提出了半全局匹配算法（Semi-Global Matching，SGM）［7］，使用動態規劃算法作為代價聚合的方法，通過多個在一維方向上的優化來代替二維圖像整體代價優化策略，該算法兼顧了局部匹配算法的效率和全局匹配算法的效果［8］。但在以動態規劃思想作為代價聚合的過程中，為保護視差的連續性，對視差不連續的區域進行了懲罰，這導致視差突變的物體邊緣區域的視差值準確度下降［9］。對于視差細化部分，子像素擬合、一致性檢查、唯一性約束以及剔除小連通域等方法是常用的視差優化方法，但是城市遙感影像具有大量的建筑陰影以及弱紋理的建筑屋頂，上述方法對這些較低信噪比的區域難以產生值得信賴的視差值［7］。

綜上所述，本文設計了一個基于改進的Census變換的匹配代價函數，它保持了Census變換對光照變化不敏感的特性，同時具有較強的抗噪能力，對匹配點相似度的刻畫更加細膩。代價聚合時，在動態規劃優化方法的基礎上，增強了物體邊緣保護能力；在常規視差細化方法的基礎上，針對城市區域衛星遙感影像中存在的大量建筑陰影、弱紋理建筑屋頂等區域導致的匹配視差不準確現象提出了針對性的細化策略。

2 基于改進Census變換的匹配代價

基于Census變換的匹配代價計算通過利用鄰域像素相對于中心像素的亮度關系作為該中心像素的局部特征信息，利用漢明距離作為衡量左右兩幅圖像中對應像素的相似度度量［10］。Census變換的步驟為取中心像素p的一鄰域M，鄰域內的像素值若是大于中心像素值便記為1，否則記為0。表達式如下：

其中：C表示變換結果；p表示中心像素；qi表示鄰域像素；I(p)和I(qi)分別表示對應的灰度值。

將中心像素p與鄰域每一像素點進行比較得到反映p點局部特征的二進制串。通過比較左右兩幅圖像對應像素點的二進制串的亦或關系得到這兩點的相似度，即匹配代價。

該方法簡單迅速，利用相對亮度作為特征，因此對于兩幅匹配圖像之間的光照變化并不敏感。但是該方法只將中心點鄰域像素劃分為兩種情況，并且劃分過程過于依賴中心點像素，若中心像素灰度值因為噪聲等影響發生較大變化，變換結果會產生較大的偏差［11］。例如，某一中心像素灰度值略大于所有鄰域像素的灰度值，因此變換后的二進制串結果應是一串0，若噪聲使得中心像素灰度值小于鄰域像素灰度值，則會得到一串1的二進制串。此外，該方法能夠區分的視差層級完全局限于所取鄰域的大小，若取5×5則只能分辨25個層級，不能得到非常細膩的視差值。

本文提出一種n階加權的Census變換算法，以適應城市衛星遙感影像大量噪聲及弱紋理區域的問題。n階加權Census算法的思想依舊沿用Census變換中利用相對亮度作為特征以應對光照變化，但不再以中心像素作為標準，而是以整個鄰域內亮度均值作為比較標準，利用鄰域內像素灰度值的最大值、最小值以及設置的n個閾值將整個鄰域劃分為n個狀態。計算方法如下：

其中：Imin與Imax分別代表p點鄰域像素內的最小灰度值與最大灰度值；n代表n個狀態；Ith1，Ith2，Ith(n-1)分別代表第1個閾值、第2個閾值以及第n-1個閾值。

由于整個鄰域像素被劃分為n種狀態，所獲得的匹配代價比兩匹配點的相似性度量更加細膩。

基于越靠近中心像素，越能代表中心像素局部信息特征的事實，對n階Census變換得到的數串進行相似度度量時，加入距離權重信息。兩匹配像素在分別經過n階Census變換之后得到兩數串，兩數串內對應數值之差的和越小，兩像素越相似，因此定義以n階加權Census變換為基礎的匹配代價計算為：

其中：Xl為左圖像素p點經n階Census變換之后的數串，Xr為右圖匹配點經n階Census變換之后的數串，xl，i與xr，i分別為Xl與Xr內一數值，c為對應鄰域像素與p點之間的距離。

3 邊緣引導優化的匹配代價聚合

匹配代價計算通常使用鄰域像素的灰度信息來計算兩像素之間的匹配程度，這很容易受到圖像噪聲、弱紋理或重復紋理等因素的影響［12］。代價聚合則是通過建立像素之間的聯系，以一定的準則，例如相鄰像素應該具有連續的視差值，來優化第一步得到的匹配代價矩陣，使匹配代價值能夠準確地反映像素之間的相關性。SGM算法［7］采用在多個一維方向的優化代替在二維圖像整體進行優化的思想，兼顧了局部匹配算法的效率和全局匹配算法的效果，因此得到了廣泛的應用［13-14］。

SGM算法的匹配代價聚合方法通過動態規劃優化的能量函數為：

其中：第一項為視差為D的所有像素的匹配代價之和；第二項表示對于像素p鄰域N p內的所有像素q，若其視差相對于p點視差變化為1，則加一懲罰值P1，第三項表示對于像素p鄰域N p內的所有像素q，若其視差相對于p點視差變化大于1，則加一懲罰值P2。

該能量函數基于這樣一個假設，即某像素鄰域內的視差值應與該像素視差值相等或相差很小，因此對于視差變化較大的區域加一較大的懲罰項P2。當一個匹配使得所有匹配點對應的匹配代價和最小時，得到最優匹配。但是城市區域遙感影像包含大量的建筑分布，因此存在視差變化較大的區域，若通過式（4）優化則建筑邊緣處的匹配效果較差。

若能夠提前知道圖像中視差變化較大的位置，則可以對這些位置設置一個較小的P2懲罰值，從而解決視差變化區域難以得到較小匹配代價問題。視差變化較大的區域一般發生在建筑屋頂的邊緣位置。圖像的梯度信息一定程度上體現了圖像的邊緣信息，梯度較大的位置較有可能位于物體邊緣位置，因此可以通過求解圖像梯度的方法獲得建筑的邊緣，進而根據梯度幅值設置合理的P2懲罰值。P2的取值原則如下：

其中：P2c為大視差懲罰值，g為該像素位置處的梯度幅值，P1為式（4）中的懲罰值。

4 基于場景的視差優化

視差細化時，通常使用子像素擬合的方法使當前獲得的整像素級視差精度提高到子像素精度，使用一致性檢測和唯一性約束來提高匹配結果的置信度，通過剔除小連通域來減小錯誤視差區域［7］。對于大多數場景，這些方法能夠達到較好的視差細化效果，但是由于城市區域衛星遙感影像存在大量的建筑陰影和建筑屋頂的弱紋理，上述方法無法在陰影區域取得準確的視差，并獲得建筑屋頂的完整視差。因此，針對城市區域需要進行陰影去除和提高屋頂視差的完整度。

4.1 陰影檢測及去除算法

建筑陰影區域的信噪比非常低，幾乎無法獲得有助于匹配的信息。陰影與非陰影的邊緣區域因為像素亮度變化大，容易被認為是視差不連續區域，從而導致誤匹配。因此，通過檢測圖像中的陰影區域，進而對視差圖中陰影區域進行處理可減小陰影帶來的影響。

本文采用雙峰陰影檢測法。城市衛星遙感影像的直方圖具有明顯的雙峰現象，如圖1所示，其中小于峰谷的像素值可認為是場景中的陰影區域［15-16］。雙峰法陰影檢測結果如圖2所示。

圖1 城市衛星遙感影像直方圖Fig.1 Histogram of urban satellite remote sensing image

圖2 雙峰法陰影檢測結果Fig.2 Shadow detection result of double peak method

陰影一般是由太陽光照射到建筑表面投影到地面或者較低的建筑屋頂造成的，因此陰影區域的視差值應與其鄰域內高程值較小的視差值相同。本文使用陰影周圍其他像素點的視差信息對陰影區域視差值進行合理推斷，并進行填充修正。陰影區域視差值修正填充流程如下：

（1）獲取陰影區域周圍鄰域像素。首先，基于形態學的方法對檢測到的陰影連通域進行膨脹操作；然后，使用膨脹之后的陰影連通域減去原始陰影連通域，由此得到陰影區域周圍鄰域像素M。計算公式如下：

其中：A為陰影連通域，B為3×3矩陣模板，⊕為形態學膨脹操作，M為陰影區域周圍鄰域像素。

（2）分析統計陰影鄰域區域像素的視差數值。由式（1）獲得M陰影鄰域統計視差信息相同的相鄰像素的個數，其中具有相同視差的連續像素點個數占比超過一定閾值ε（本文設置為0.25），且高程信息較小的高程值即是合理的陰影區域填充視差值。圖3為陰影鄰域區域視差示意圖。圖4為去除陰影影響的圖像，由此可以看出該步驟對減小陰影影響具有非常好的效果。

圖3 陰影鄰域區域視差示意圖Fig.3 Disparity diagram of shaded neighborhood area

圖4 陰影區域視差修正結果Fig.4 Results of disparity correction in shadow area

4.2 基于建筑區域提取的建筑屋頂視差優化

由于城市內建筑物屋頂大多為弱紋理區域，該區域像素值變化不大，可提供給立體匹配算法的可用信息較少，造成視差圖中屋頂區域存在一些空洞，并且鄰域內視差分布不平整。這種情形下，可采用基于區域生長的分割算法。由于建筑屋頂大多為矩形，因此在分割算法中輔以直線段檢測算法，進一步分割邊緣位置，得到較為精確的建筑屋頂分割區域。城市中大多數的建筑屋頂都是水平平面，同一建筑屋頂具有相同的高程，也即擁有相同的視差。因此，以區域分割輔以直線段檢測的方法得到建筑屋頂區域；平均建筑屋頂區域視差得到平整的屋頂視差［17］。分割算法可以準確地檢測建筑屋頂區域，得到完整的建筑屋頂視差。

基于區域分割，輔以直線段檢測的建筑屋頂分割算法和建筑區域視差平整算法分別用偽代碼Algorithm 1，Algorithm 2表示，如圖5所示。

圖5 算法偽代碼Fig.5 Algorithm pseudocode

5 實驗及分析

本文提出的立體匹配算法以基于n階Census變換的匹配代價為代價計算方法，以保留物體邊緣信息動態規劃方法為代價聚合，并以WTA策略作為代價計算方法，同時針對城市陰影以及建筑屋頂進行視差優化。首先在Middle-Bury實驗數據集上對本算法的通用性及準確性進行了驗證，然后在真實的城市遙感影像對上進行了實驗。

5.1 算法代價計算及代價聚合改進效果的實驗

MiddleBury實驗數據數據集是2001年Schaestein等為統一評價已有的各種立體匹配算法的效果而創建的數據集。該數據集包含多個場景，每個場景包含左視圖和右視圖，以及對應的左視差圖和右視差圖［2］，可以方便對算法得到的視差圖與真實視差圖進行比較，因此在立體匹配算法的效果評價中得到了廣泛的應用。MiddleBury實驗數據集可在https：//vision.middlebury.edu/stereo/data/獲得。本文使用Middle-Bury實驗數據數據集中多個場景對本文算法的效果進行了驗證，并與網站中SGM算法的實驗結果進行比較。

MiddleBury實驗數據集中Piano場景的實驗結果如圖6所示。相較于SGM，本文提出的算法可得到更完整的物體輪廓和精準的物體邊緣視差圖。

圖6 MiddleBury數據集Piano場景圖像的實驗結果Fig.6 Experimental renderings of Piano scene images in MiddleBury dataset

在MiddleBury實驗數據集的網站中默認展示了各種立體匹配算法（包括SGM算法）誤差大于2的像素占總像素數的百分比。因此，這里將視差誤差小于2的像素認定為匹配正確，統計正確匹配像素占總像素數的比例作為準確率（acc），以反映計算所得視差和真實視差之間的差異性，計算方法如下：其中：cond Small(S，T，a，b)表示若S小于T為真，則取a，否則取b；d表示計算得到的視差；d G表示真實視差；M表示圖像像素個數。

表1展示了本文算法與SGM算法在Middle-Bury實驗數據多個場景下準確率的統計結果。由表1可知，在改進匹配代價函數以及代價聚合優化函數之后，本算法在多個場景下都取得了更高的準確率。在各場景下本算法的準確率平均值相較SGM算法提升了4.43%。

表1 Middlebury實驗數據集上算法的準確率Tab.1 Accuracy rate of algorithms on Middlebury experimental dataset （%）

5.2 真實遙感影像實驗

World View-2于2009年10月發射，是第一顆高分辨率8波段多光譜商業衛星，能夠提供0.5 m分辨率的立體影像產品。為了檢驗該算法在城市區域衛星遙感影像對上進行立體匹配的可行性，實驗使用的WorldView-2衛星遙感影像對如圖7所示。

圖7 World View-2國家圖書館區域遙感影像對Fig.7 WorldView-2 national library regional remote sensing image pair

為展示本文算法在應對場景亮度變化方面的效果，使用WorldView-2衛星拍攝的國家圖書館區域遙感影像對進行了實驗，該場景中包含明顯的亮度變化，嚴重影響立體匹配算法的準確性，實驗結果如圖8所示。

由圖8可以看出，以基于n階加權Census算法作為匹配代價函數以及保留物體邊緣的代價聚合策略的立體匹配算法，可以很好地應用于城市衛星遙感影像對的立體匹配中，同時本文算法比SGM算法具有更少的誤匹配點，圖像邊緣更細膩。

圖8 WorldView-2遙感影像對上不同算法的效果圖Fig.8 Algorithm renderings on WorldView-2 remote sensing image pair

為驗證本文所提的在建筑屋頂區域的視差優化方法，使用與圖7同一景影像中建筑密度以及復雜程度更高的場景進行了實驗，實驗原圖如圖9所示，該場景中還包含大量的陰影區域。建筑屋頂視差優化的實驗結果如圖10所示，可以看出，經過優化后的視差圖更加平滑，建筑屋頂完整度更高。

圖9 建筑屋頂視差修正實驗原圖Fig.9 Original images of building roof parallax in correction experiment

圖10 建筑屋頂區域視差優化效果圖Fig.10 Disparity optimization images of roof area of building

由于遙感影像對沒有標準的Ground Truth視差圖作為比較，對于城市區域遙感影像對的匹配結果只從主觀視覺進行了分析。為了定量分析其效果，本文利用基于有理函數的三維重建方法，由本文算法獲得的視差圖10（b）計算對應的數字 高 程 模 型（Digital Elevation Model，DEM）圖［18］，結果如圖11所示。

由先驗知識可知，城市建筑的屋頂大多為水平平面，其高程一致，因此可通過計算建筑屋頂平面的高程方差來評價算法的準確度。選取圖11中不同高程、不同建筑的a，b，c 3個區域進行分析，分別計算3個區域內的高程均值以及高程方差，結果如表2所示。

圖11 基于WorldView-2衛星遙感影像視差圖的DEM圖Fig.11 DEM map based on disparity map of WorldView-2 satellite remote sensing image

由表2可知，a區域的高程方差為1.739 0，b區域的高程方差為0.009 8，c區域方差為0.37，3個區域的方差平均值為0.71，這表示本算法可以得到比較平整的建筑屋頂重建效果，準確性較高。

表2 建筑屋頂高程數據統計信息Tab.2 Statistical information of building roof elevation data

6 結 論

本文為解決城市遙感影像多建筑陰影以及建筑屋頂無紋理導致立體匹配效果不理想的問題，提出了適用于城市遙感影像對的立體匹配方法。首先，提出了基于改進Census變換的代價計算方法。然后，利用建筑邊緣信息強化了算法在屋頂邊緣處匹配的穩定性。最后，給出了建筑陰影區域以及建筑屋頂區域的視差修正方法。實驗結果表明：在MiddleBury數據集上，本文算法的準確率比經典SGM算法高4.43%；在城市區域WorldView-2立體影像對上，建筑屋頂的高程方差為0.71。該方法滿足基于城市衛星遙感影像對獲取高精度視差圖的要求，能夠為城市三維重建提供良好的條件。