鞏丹超，韓軼龍，黃 旭

1. 西安測繪研究所,陜西 西安 710054； 2. 地理信息工程國家重點實驗室, 陜西 西安 710054； 3. 武漢大學(xué)遙感信息工程學(xué)院，湖北 武漢 430079； 4. 武漢市工程科學(xué)技術(shù)研究院, 湖北 武漢 430019

立體影像密集匹配是在立體影像中，逐像素或者接近逐像素地尋找兩張影像之間的同名點，并根據(jù)同名光線對對相交的原理，生成三維點云的過程[1]。立體影像密集匹配具有點云密度大、匹配結(jié)果穩(wěn)定、成本低等突出優(yōu)勢，因此可以廣泛應(yīng)用于一些智能3D場景，如：測繪制圖、變化監(jiān)測、國防軍事、智慧城市、無人駕駛、虛擬現(xiàn)實、游戲動畫、搶險救災(zāi)等[2-6]。

立體影像密集匹配常采用固定匹配窗口來計算和比較同名點之間的特征相似性，其基本假設(shè)是局部匹配窗口內(nèi)，所有像素的視差均一致。該假設(shè)在大部分場景下適用，但是在建筑物邊緣區(qū)域，由于存在明顯的視差/高程階躍變化，該假設(shè)無法成立。在建筑物邊緣區(qū)域，固定匹配窗口會帶來較大的匹配歧義性，從而降低建筑物邊緣的匹配精度，具體表現(xiàn)為：建筑物邊緣不規(guī)則，建筑物邊緣外擴等問題，如圖1所示。圖1(a)、(c)表示原始影像。圖1(b)、(d)表示對應(yīng)的立體影像密集匹配結(jié)果。圖1(a)、(c)中紅色的線，以及圖1(b)、(d)中白色的線是對應(yīng)的建筑物直線邊緣特征。其中，密集匹配算法采用著名的半全局密集匹配算法(semi-global matching,SGM)[7-8]，匹配特征算子采用Census[9]，匹配窗口大小為9×9像素。從圖1可以看出，SGM密集匹配算法的匹配結(jié)果，在建筑物邊緣區(qū)域存在明顯的外擴問題，從而影響后續(xù)建筑物的三維重建精度。除了SGM算法以外，其他基于固定匹配窗口的全局/半全局密集匹配算法，也存在同樣的問題[10-13]。除了固定匹配窗口因素以外，匹配常用的視差平滑約束，同樣也是邊緣外擴的因素之一。

圖1 密集匹配算法在建筑物邊緣區(qū)域的匹配結(jié)果Fig.1 The dense matching result in the edge area of the buildings

為了獲取高精度的建筑物邊緣，目前主要有兩種解決方案。第1種方法是在密集匹配過程中加入線特征約束或者面特征約束，削弱建筑物邊緣處的匹配歧義性，從而提高建筑物邊緣處的匹配精度[14-15]。這類方法往往可以取得比傳統(tǒng)SGM更高的匹配精度，但是這些方法的線、面特征約束是根據(jù)自身算法的實際情況來定制設(shè)計的，不一定適用于其他密集匹配方法，比如，文獻[14]采用動態(tài)規(guī)劃思想來優(yōu)化密集匹配結(jié)果，并采用影像線特征來自適應(yīng)調(diào)整動態(tài)規(guī)劃方程中的懲罰項系數(shù)，以提高邊緣處的密集匹配精度；但是對于不含懲罰項系數(shù)的密集匹配方法(如：基于最小生成樹的密集匹配方法[16]，基于深度學(xué)習(xí)的密集匹配方法[17])，是無法使用這類線特征約束。此外，目前測繪單位和地理信息平臺存儲大量數(shù)字表面模型產(chǎn)品(digital surface model,DSM)，第1種方法只能用于生產(chǎn)新的DSM，無法對現(xiàn)有DSM產(chǎn)品進行邊緣優(yōu)化。第2種方法是通過后處理的方式，直接對密集匹配視差圖或者DSM中的建筑物邊緣進行優(yōu)化。這類后處理的方法，要么考慮灰度特征分布和實際視差/高程分布之間的相關(guān)性，修正所有像素的視差/高程，使其與對應(yīng)的影像灰度分布近似，從而達到邊緣優(yōu)化的目的；要么從影像中提取的線特征或者面特征，并修正視差/高程邊緣，使其與影像線特征或面特征盡量接近。因此，通過后處理來優(yōu)化邊緣的方法可以分為兩類：①基于灰度特征分布的邊緣銳化算法；②基于線面特征的優(yōu)化算法。

大多數(shù)基于灰度特征分布的邊緣銳化算法通過局部窗口內(nèi)的灰度特征，來改善視差/高程，從而也稱為局部邊緣銳化算法[18-24]。該類算法假設(shè)在局部窗口內(nèi)，灰度相似的像素，其視差/高程也是一致的。由于建筑物邊緣往往也是灰度階躍的邊緣，因此可以采用這類局部邊緣銳化算子進行濾波，通過灰度分布來重新調(diào)整建筑物邊緣，從而獲得更加精細的建筑物邊緣。但是受限于窗口大小，該類局部銳化算法很難改正建筑物邊緣處一些較大的誤匹配。為了解決邊緣處較大誤匹配的問題，需要從更加全局的角度來修正建筑物邊緣，也稱為全局邊緣銳化算法[25-26]。該類算法假設(shè)整張影像中，相鄰且灰度相似的像素，其視差/高程應(yīng)該盡量接近，從而將視差/高程優(yōu)化問題轉(zhuǎn)化為全局能量函數(shù)的最優(yōu)解計算問題。全局邊緣銳化算法能夠有效優(yōu)化建筑物邊緣，但是由于影像紋理(非邊緣區(qū)域)也具有灰度階躍特性，因此在影像紋理區(qū)域容易出現(xiàn)“銳化”假象。

基于線面特征的優(yōu)化算法首先從影像中提取線、面等幾何特征，然后改正建筑物邊緣的視差/高度，使得邊緣特征與影像的幾何特征相符。文獻[7]通過影像分割算法，提取影像中的面特征，并將每個面特征作為一個視差平面進行優(yōu)化，從而達到改善邊緣的目的。但是，過大的面特征不一定滿足視差平面假設(shè)，可能導(dǎo)致匹配精度下降。文獻[27]進一步采用直線特征提取算子(line segment detector,LSD)來提取建筑物的邊緣，并假設(shè)邊緣特征附近的局部區(qū)域滿足平面約束，從而實現(xiàn)邊緣優(yōu)化。該方法計算效率高，但是由于只考慮直線附近區(qū)域的視差/高程優(yōu)化，容易導(dǎo)致與直線區(qū)域以外的像素之間的視差/高程斷裂，并且容易將直線附近區(qū)域的地形地貌強行抹平。文獻[28]通過平移邊緣附近的像素，來改善建筑物邊緣，其本質(zhì)是舍棄直線邊緣附近的所有視差/高程，然后以直線邊緣為中心，分別利用邊緣以外的建筑物視差/高程信息和地面視差/高程信息，向中心邊緣進行內(nèi)插。但是，該方法同樣會將建筑物邊緣附近的地物抹平。文獻[12]計算邊緣像素的視差和匹配代價，通過取最小代價的視差，來優(yōu)化邊緣。文獻[29]通過超像素分割來確定物體邊緣，并通過加權(quán)中值濾波來優(yōu)化物體邊緣。

為了解決局部邊緣銳化算法難以解決較大誤匹配的問題，并解決基于線面特征優(yōu)化算法強行抹平附近地物的問題，本文提出了一種基于直線特征約束的建筑物邊緣全局優(yōu)化方法，能夠在解決較大的建筑物邊緣誤差的同時，保留建筑物附近的地物。該方法的貢獻在于，將建筑物邊緣優(yōu)化問題，轉(zhuǎn)化為一個新的全局能量函數(shù)的最優(yōu)解計算問題，并采用圖割算法進行全局最優(yōu)計算。全局能量函數(shù)由代價項和正則項兩部分組成?？紤]到保留建筑物邊緣附近的實際地物，代價項的設(shè)計主要根據(jù)估計視差/高程與原始視差/高程之間的距離，以及該像素的視差/高程置信度來共同決定的。正則項的設(shè)計，則是依據(jù)灰度相近的相鄰像素，其視差/高程也是相似的原則來設(shè)計。本文所提出的建筑物邊緣優(yōu)化方法，能夠有效改正建筑物邊緣的誤匹配問題，解決建筑物邊緣外擴的問題，最終重建高精度的建筑物三維邊緣。

1 基于線特征的建筑物密集匹配邊緣的全局優(yōu)化方法

本文重點解決密集匹配中，建筑物邊緣精度較低的問題。具體流程如圖2所示。算法的輸入是原始核線影像和視差圖(如圖2(a)所示)，或者正射影像和對應(yīng)的DSM，其中，圖2(a)中視差圖的視差范圍是34像素至96像素。首先，采用直線特征提取算子LSD[30]，來檢測并提取建筑物的邊緣，如圖2(b)所示，紅色、白色的線分別表示采用LSD算法檢測和提取的建筑物邊緣。兩條邊緣線均位于屋頂邊緣。然后，分別優(yōu)化每一條直線特征，將直線特征的優(yōu)化問題轉(zhuǎn)化為全局能量函數(shù)最優(yōu)解的計算問題。全局能量函數(shù)共有代價項和正則項兩部分構(gòu)成，如圖2(c)所示。其中，代價項的計算是根據(jù)估計視差/高程與原始視差/高程之間的距離Δd，當前像素與建筑物邊緣之間的距離s等因素，來共同決定的。Δd越小、s越小，則代價越小，反之則越大；正則項的計算，則是依據(jù)灰度相近的相鄰像素，其視差/高程也是相似的原則來設(shè)計。全局能量函數(shù)的具體形式，如1.2節(jié)式(1)所示。圖2(c)中的3個藍色的圓形表示邊緣附近的3個像素，黑色的連接線表示全局優(yōu)化中的視差平滑約束，線越粗表示約束越強。灰度相近的像素之間的約束較強，從而在后續(xù)全局優(yōu)化中鼓勵他們的視差/高程一致，反之，則灰度差異較大像素之間的約束較小。最后采用圖割算法求解整型最優(yōu)解，并采用基于圖像引導(dǎo)濾波的方法獲得連續(xù)最優(yōu)解。

圖2 本文算法流程Fig.2 The pipeline of the algorithm proposed in this paper

1.1 建筑物邊緣檢測

本文采用LSD算子來檢測和提取影像的直線特征。但是，受地物紋理、陰影等因素的影響，采用LSD算子提取的直線特征，大部分不屬于建筑物的邊緣?？紤]到建筑物邊緣的視差/高程階躍特性，本文采用線緩沖區(qū)方法[27]來從LSD直線特征中，進一步提取建筑物的邊緣。具體流程為：①首先，采用LSD算子檢測影像的直線特征，如圖3(b)所示；②以每一條直線為中心，建立緩沖區(qū)，統(tǒng)計直線兩邊的視差/高程情況；③如果直線位于視差/高程的階躍區(qū)域，則作為建筑物的邊緣，提取出來，如圖3(c)所示。本文算法實質(zhì)以直線邊緣為基本單元，優(yōu)化直線邊緣附近的高程/視差。因此，除了文獻[27]提出的基于緩沖區(qū)的直線邊緣提取方法之外，其他能夠提取直線邊緣的方法，同樣適用于本文算法，如手動提取直線。

圖3 建筑物邊緣檢測結(jié)果Fig.3 Building edge detection results

1.2 全局能量函數(shù)設(shè)計

提取出建筑物邊緣的直線特征之后，本文分別對每一條邊緣直線特征進行優(yōu)化。本文算法的核心思想是改變邊緣直線附近像素的視差/高程，從而使得改正后的視差邊緣/高程邊緣與影像的邊緣直線盡量接近。因此，首先需要尋找邊緣直線附近的像素。本文采用“內(nèi)/外緩沖區(qū)”的思路，來統(tǒng)計邊緣直線附近的像素，如圖4所示。白色直線表示建筑物的邊緣直線。紅色方框表示內(nèi)緩沖區(qū)。本文算法通過改變內(nèi)緩沖區(qū)內(nèi)所有像素的視差/高程，使得視差/高程邊緣與影像邊緣盡量接近。藍色方框表示外緩沖區(qū)。外緩沖區(qū)內(nèi)像素的視差/高程是不做任何變化的。外緩沖區(qū)的作用有兩個：①本文采用固定窗口來計算能量函數(shù)的代價項，因此，在計算代價項的過程中，會有一部分代價窗口超出內(nèi)緩沖區(qū)的范圍。超出的部分可以采用外緩沖區(qū)來計算。具體代價計算詳見1.2.1小節(jié)。②傳統(tǒng)的圖割算法只能提供整型最優(yōu)解，因此，在內(nèi)緩沖區(qū)和外緩沖區(qū)之間，存在明顯的視差/高程斷裂。在后續(xù)連續(xù)最優(yōu)解計算過程中，將考慮外緩沖區(qū)的像素，使得內(nèi)外緩沖區(qū)之間的邊界平滑過渡。外緩沖區(qū)的半徑rout設(shè)置為代價窗口的一半?？紤]到立體影像密集匹配窗口通常采用9×9窗口[31-32]，本文將內(nèi)緩沖區(qū)半徑rin設(shè)置為10。

圖4 邊緣直線的內(nèi)/外緩沖區(qū)Fig.4 Inner/outer buffer zone with straight edges

在統(tǒng)計邊緣直線附近的像素之后，本文通過改變內(nèi)緩沖區(qū)中像素的視差/高程，來達到邊緣優(yōu)化的目的?？傮w來說，本文的邊緣優(yōu)化準則包含兩個方面：①優(yōu)化后的視差/高程，與原始視差/高程盡量接近；②相鄰且灰度相近的像素，其視差也是相似的。根據(jù)上述兩個準則，可以將視差邊緣優(yōu)化問題，轉(zhuǎn)化為全局能量函數(shù)的最優(yōu)解計算問題，如式(1)所示

minE(D)=∑p∈BinSc(p)·Cost(p,d)+

∑p,q∈NP·wc(p,q)·we(p,q)·

min(|dp-dq|,σ)

(1)

式中，D表示內(nèi)緩沖區(qū)中所有像素的估計視差/高程的集合；E(D)表示全局能量函數(shù)的目標函數(shù)值；Bin表示內(nèi)緩沖區(qū)所有像素的集合；p表示內(nèi)緩沖區(qū)中的任意一個像素；Cost(p,d)表示像素p對應(yīng)于視差/高程d的代價，是根據(jù)估計視差/高程與原始視差/高程之間的距離來計算，詳見1.2.1小節(jié)介紹；Sc(p)表示像素p的置信度，表示像素p的代價計算的可靠性，用于降低誤匹配點的代價權(quán)重，具體詳見1.2.1小節(jié)介紹；N表示所有鄰域像素的集合；q表示p的鄰域像素；dp、dq表示像素p、q的估計視差/高程；σ表示截斷值，用于定義視差階躍條件的閾值；P表示懲罰項系數(shù)；wc(p,q)表示根據(jù)像素p、q的灰度差，計算出來的權(quán)重，一般灰度差越小，權(quán)越大；we(p,q)表示根據(jù)像素p、q的位置，定義的權(quán)重。如果p、q位于直線邊緣的同一側(cè)，則權(quán)we為1，否則，權(quán)we為0。we用于禁止直線兩邊的優(yōu)化結(jié)果相互影響，從而能夠獲得更加銳利的邊緣。本文將能量函數(shù)式(1)中的第1項稱為代價項，將第2項稱為正則項。代價項和正則項的具體介紹，詳見1.2.1小節(jié)和1.2.2小節(jié)。

1.2.1 代價項計算

大多數(shù)算法[7,27-28]采用平面模型來改變邊緣直線附近的視差/高程，因此，會強行刪掉邊緣直線附近的實際地物。為了解決這個問題，本文將估計視差/高程與原始視差/高程之間的距離，作為代價計算的一個依據(jù)。當距離越小，則代價越??；反之，距離越大，則代價越大，從而盡可能地保留直線附近的實際地物，如式(2)所示

(2)

但是，原始視差/高程同樣存在一定的誤匹配問題，單純采用上述距離測度計算的代價，有可能在優(yōu)化中將誤匹配點也保留下來了。因此，本文定義了視差/高程置信度準則，用于衡量內(nèi)緩沖區(qū)中每個像素是正確匹配點的概率。置信度越高則表示正確匹配點的概率越大，則在后續(xù)全局優(yōu)化中的權(quán)重越大；反之，則權(quán)重越低。本文定義了兩種置信度準則：第1種是像素距離邊緣直線的距離，考慮到匹配邊緣的外擴問題，距離越近則是正確匹配點的概率越小，距離越遠則是正確匹配點的概率越大；第2種是像素與周圍灰度相近像素之間的視差/高程差異，差異越小則是正確匹配點的概率越大，差異越大則是正確匹配點的概率越小。因此，綜上兩個準則，置信度的計算如式(3)所示

(3)

1.2.2 正則項計算

能量函數(shù)的正則項用于對緩沖區(qū)內(nèi)的所有像素引入視差/高程平滑約束。相鄰像素的灰度越接近，則視差/高程平滑約束越強，從而在能夠改正邊緣直線附近的誤匹配點的同時，盡可能多地保留緩沖區(qū)內(nèi)的細節(jié)信息?？偟膩碚f，本文定義的正則項與3個要素有關(guān)：①相鄰像素的灰度差；②相鄰像素是否位于直線的同一側(cè)；③相鄰像素的視差/高程差異。

相鄰像素灰度越接近，則正則項的視差/高程平滑約束越強。因此，本文采用高斯核函數(shù)來定義相鄰像素之間的灰度接近程度，如式(4)所示

(4)

式中，Ip、Iq分別表示像素p、q的灰度；σc表示與灰度相關(guān)的高斯核函數(shù)平滑因子，本文取10。當像素p、q之間的灰度越接近，灰度權(quán)wc越接近1。

考慮到有些情況下，邊緣直線兩側(cè)的灰度差異較小。在這種情況下，灰度權(quán)wc會較大，從而導(dǎo)致邊緣直線兩側(cè)的視差/高程平滑過渡，降低邊緣直線的銳利程度。因此，為了保證優(yōu)化后的視差/高程邊緣，與影像邊緣直線特征盡量接近，本文禁止直線兩側(cè)像素之間的正則平滑約束，如圖5所示。其中，白色直線表示建筑物直線邊緣；藍色圓圈代表緩沖區(qū)內(nèi)的像素；黃色線表示相鄰像素之間的正則項，線的粗細表示正則項平滑約束的大小；紅色線框表示內(nèi)緩沖區(qū)。在本文算法中，只有位于同一側(cè)的相鄰像素才有正則項的平滑約束。

如果相鄰像素屬于同一側(cè)緩沖區(qū)，則需要引入正則項平滑約束，將正則項中的同側(cè)權(quán)值we(p,q)定義為1；反之，則定義為0，如式(5)所示

(5)

1.3 連續(xù)最優(yōu)解計算

本文采用圖割算法[34]，獲取全局能量函數(shù)(式(1))的最優(yōu)解。但是，傳統(tǒng)的圖割算法只能獲得整型最優(yōu)解，從而導(dǎo)致緩沖區(qū)內(nèi)的像素和緩沖區(qū)外的像素出現(xiàn)輕微的視差/高程斷裂，如圖6(a)所示。其中，白色的線表示影像上的邊緣直線特征；紅色的線是剖分圖軌跡，對應(yīng)于右側(cè)的折線剖分圖。為了獲得連續(xù)解，需要綜合考慮外緩沖區(qū)像素的視差/高程，對內(nèi)緩沖區(qū)像素進行平滑濾波。為了保證濾波后的視差/高程邊緣的幾何精度不受影響，本文采用保邊濾波算子(如基于圖像引導(dǎo)的濾波[22])，來對緩沖區(qū)內(nèi)的像素進行平滑濾波?？紤]到一些房子的屋頂與地面的灰度差異不是特別明顯，即使采用保邊濾波算子，也無法在濾波過程中保持建筑物邊緣的尖銳性。因此，本文以邊緣直線特征為界限，分別對邊緣直線兩側(cè)進行濾波。當有濾波窗口中的像素超出邊緣直線時，超出部分的像素不參與濾波計算，如圖6(a)所示。圖6(a)中的綠色方格表示當前需要濾波的像素；藍色方格表示周圍的8鄰域像素。其中，最右側(cè)的3個方格超出邊緣直線部分，則這3個方格不參與濾波計算。

圖6 連續(xù)最優(yōu)解計算示意Fig.6 Continuous optimal solution calculation

基于圖像引導(dǎo)的濾波[22]是一種高效的保邊平滑算子，其核心思想是將局部窗口內(nèi)的像素視差/高程進行加權(quán)平均。窗口內(nèi)每個像素的權(quán)值大小取決于局部視差和灰度之間的線性模型。由于視差/高程邊緣往往也是灰度邊緣，因此基于圖像引導(dǎo)的濾波具有很好的保邊性質(zhì)。本文采用基于圖像引導(dǎo)的濾波算子，來進一步改善內(nèi)緩沖區(qū)中的整型視差/高程。在濾波過程中，外緩沖區(qū)中的像素同樣參與計算，從而能夠為內(nèi)緩沖區(qū)視差/高程的優(yōu)化提供邊界約束。基于圖像引導(dǎo)的濾波算子的定義如式(6)所示

(6)

2 試驗與結(jié)果分析

本文分別采用航拍立體像對數(shù)據(jù)集和衛(wèi)星立體像對數(shù)據(jù)集，從定量和定性兩個角度來分析本文算法的正確性和有效性。其中，航拍立體像對數(shù)據(jù)集包括International Society for Photogram-metry and Remote Sensing(ISPRS)提供的Vaihingen航拍數(shù)據(jù)集和Toronto航拍數(shù)據(jù)集；衛(wèi)星立體像對數(shù)據(jù)集包括約翰斯·霍普金斯大學(xué)和IAPRA提供的杰克遜維爾地區(qū)的WorldView-3衛(wèi)星影像數(shù)據(jù)集和對應(yīng)的激光點云(LiDAR)DSM。

2.1 航拍立體像對數(shù)據(jù)集

該部分試驗主要測試視差圖中建筑物邊緣的優(yōu)化效果。其中，Toronto數(shù)據(jù)集由Microsoft Vexcel’s UltraCam-D(UCD)相機拍攝，影像大小為2000×2000像素，地面分辨率為15 cm；Vaihingen數(shù)據(jù)集由Intergraph/ZI DMC相機拍攝，影像大小為2000×2000像素，地面分辨率為8 cm。兩個數(shù)據(jù)集共由9個立體像對組成。每個立體像對均重采樣成核線影像，并采用SGM半全局密集匹配算法[8]計算對應(yīng)的視差圖，并采用左右一致性檢測剔除誤匹配點，然后再進行視差內(nèi)插。具體測試數(shù)據(jù)集和對應(yīng)的原始視差圖如圖7所示。圖7的上半部分表示原始核線影像，下半部分表示SGM密集匹配結(jié)果。所有立體像對均包含較密集的房屋，而傳統(tǒng)的密集匹配算法在這些建筑物的邊緣區(qū)域的匹配精度較差。

圖7 航拍數(shù)據(jù)集與SGM密集匹配視差Fig.7 Aerial data and dense matching disparity using SGM

采用本文算法，最常用的雙邊濾波算法和加權(quán)中值濾波算法，分別對SGM密集匹配結(jié)果進行邊緣優(yōu)化，并對各個算法的邊緣優(yōu)化結(jié)果進行性能評估、對比和分析。其中，雙邊濾波算法[33]和加權(quán)中值濾波算法[24]均假設(shè)距離相近且灰度相近的像素，其視差應(yīng)該是一致的。由于建筑物邊緣往往也是灰度階躍的影像邊緣，因此，兩種方法均能對建筑物邊緣的匹配結(jié)果進行優(yōu)化，并且已經(jīng)廣泛應(yīng)用于一些立體影像密集匹配算法中[13,35]。為了與本文算法的參數(shù)盡量接近(緩沖區(qū)半徑為10)，雙邊濾波和加權(quán)中值濾波兩種方法均采用11×11的濾波窗口。由于本文算法僅僅優(yōu)化建筑物的邊緣，因此，本文僅在提取的建筑物邊緣區(qū)域(1.1節(jié)介紹的方法)分別評估3種方法的邊緣優(yōu)化結(jié)果。由于建筑物邊緣的誤匹配，絕大多數(shù)是由于建筑物外擴導(dǎo)致的，因此，本文根據(jù)直線緩沖區(qū)對應(yīng)地面那一側(cè)內(nèi)的誤匹配點數(shù)目，來評估各種算法的優(yōu)化性能，如圖8所示。其中，白色直線表示建筑物邊緣；黃色像素表示誤匹配點；藍色像素表示正確匹配點；紅色方框表示直線的緩沖區(qū)。本文根據(jù)圖8中的黃色像素和藍色像素數(shù)目，來綜合評價邊緣優(yōu)化算法的性能。

perbad=|Serror|/|Serror+Scorrect|

(7)

式中，|·|表示集合中所有元素的數(shù)目；perbad表示誤匹配像素百分比。perbad越小，表示邊緣優(yōu)化效果越好

σlow=

(8)

采用本文優(yōu)化算法，雙邊濾波和加權(quán)中值濾波，分別對圖7中9個立體像對的SGM密集匹配結(jié)果進行邊緣優(yōu)化，并采用式(7)和式(8)兩個測度，分別對原始匹配結(jié)果、本文優(yōu)化算法結(jié)果、雙邊濾波結(jié)果和加權(quán)中值濾波結(jié)果進行性能評定。評定結(jié)果如圖9(a)和(b)所示。

圖9 航拍數(shù)據(jù)集不同方法的性能對比結(jié)果Fig.9 Performance comparison results using different edge-refine methods on aerial dataset

由圖9可以看出，在兩種性能評定測度中，本文算法、雙邊濾波方法和加權(quán)中值方法均能夠提高建筑物邊緣的匹配效果。在誤差像素百分比測度perbad中，原始視差圖的誤匹配點百分比為25.53%，加權(quán)中值方法平均能夠?qū)⒔ㄖ镞吘壍恼`匹配像素百分比減少到20.04%，雙邊濾波方法平均能夠?qū)⒄`匹配像素百分比減少到14.21%，而本文算法平均能將誤匹配像素百分比減少到9.12%。在方差測度σlow中，原始視差圖建筑物邊緣區(qū)域的方差平均為5.86像素，加權(quán)中值方法平均能夠?qū)⒎讲顪p少到5.48像素，雙邊濾波方法平均能夠?qū)⒎讲顪p少到4.34像素，而本文算法平均能將方差減少到3.63像素。總體來說，本文算法在兩種測度中的整體表現(xiàn)均是最好的。具體到每個立體像對的邊緣優(yōu)化結(jié)果評定中，本文算法在9個像對中的表現(xiàn)，同樣也都是最好的。這是由于雙邊濾波方法和加權(quán)中值方法，均是采用局部窗口，用于優(yōu)化的信息較少；而本文方法用到了邊緣直線緩沖區(qū)內(nèi)的所有信息，并進行全局優(yōu)化，從而優(yōu)化效果更好。為了更全面地評估本文算法的性能，本文以第1對核線立體影像密集匹配的邊緣優(yōu)化結(jié)果為例，列出一些優(yōu)化前后的建筑物邊緣匹配結(jié)果的對比，以及對應(yīng)的建筑物剖面圖，如圖10所示。

圖10 本文算法的邊緣優(yōu)化效果Fig.10 Building edge refined results using our proposed algorithm

在圖(a)展示了原始影像、原始視差圖和本文算法優(yōu)化后視差圖的縮略圖結(jié)果。由于建筑物邊緣的誤匹配像素往往只有幾個像素的誤差，因此，在縮略圖中看不出太明顯的變化。為了進一步展示本文算法的效果，本文挑選了3個紅框所示的區(qū)域，并放大顯示，如圖10中(b)圖—(d)圖所示。在圖10(b)圖—(d)圖的所有原始影像中，紅色直線表示提取的影像直線邊緣；所有原始視差圖中，白色直線表示提取的影像直線邊緣，紅色直線表示建筑物剖面圖軌跡，對應(yīng)最右側(cè)的剖面圖，用于檢查優(yōu)化前后的建筑物邊緣變化。剖面圖中，紅色的線表示原始視差圖邊緣；綠色的線表示經(jīng)過本文算法修正后的邊緣；藍色的線表示影像直線邊緣。影像直線邊緣與地面之間的視差差異，主要是根據(jù)影像邊緣直線的位置，及其緩沖區(qū)兩側(cè)平均視差所共同決定的。原始SGM密集匹配結(jié)果，其建筑物邊緣距離影像直線邊緣會外擴若干個像素，從而降低建筑物邊緣的三維重建精度。本文算法能夠高精度地將建筑物邊緣修正到影像直線邊緣附近0～1個像素。在圖10(a)中，本文算法修正后的邊緣，與影像直線邊緣完全重合，而在圖10(b)、(c)中，本文算法修正后的邊緣，僅僅與影像直線邊緣差1個像素，從而充分說明本文算法的正確性和有效性。

2.2 衛(wèi)星立體像對數(shù)據(jù)集

該部分試驗主要測試DSM中建筑物邊緣的優(yōu)化效果。其中，WorldView-3衛(wèi)星數(shù)據(jù)集總共包含4個立體像對，空間分辨率為0.3 m，像對重疊區(qū)域約為0.5 km2。激光點云DSM的空間分辨率為0.5 m，范圍約為0.1 km2。首先，將每個立體像對重采樣成核線影像，采用SGM半全局密集匹配算法[8]計算對應(yīng)的視差圖，并采用左右一致性檢測剔除誤匹配點；然后根據(jù)視差圖生產(chǎn)對應(yīng)的DSM，并進行DSM內(nèi)插。為了便于利用激光點云DSM來評估精度，每個衛(wèi)星立體像對DSM和正射影像均切割成與激光點云DSM同樣的范圍，并將兩者的空間分辨率重采樣成0.5 m。具體測試數(shù)據(jù)集、對應(yīng)的原始DSM和激光點云DSM如圖11所示。圖11的上半部分表示衛(wèi)星正射影像，下半部分表示立體像對DSM和激光點云DSM。

采用本文算法和基于緩沖區(qū)平面模型的建筑物邊緣修正方法(plane based boundary refinement,PBR)[27]，分別對圖11中的4組原始DSM進行邊緣優(yōu)化，并以激光點云DSM作為高程真值，對各個算法的邊緣優(yōu)化結(jié)果進行性能評估、對比和分析。其中，基于緩沖區(qū)平面模型的建筑物邊緣修正方法(PBR)[27]假設(shè)緩沖區(qū)內(nèi)的高程/視差滿足平面約束，通過計算緩沖區(qū)內(nèi)的平面方程，來優(yōu)化建筑物的邊緣。由于兩種算法均僅僅優(yōu)化建筑物的邊緣，因此，本文僅僅在提取的建筑物邊緣區(qū)域(1.1節(jié)介紹的方法)評估兩種方法的邊緣優(yōu)化精度。為了客觀評價兩種算法的優(yōu)化效果，以激光點云DSM作為真值，定量分析優(yōu)化后DSM與原始DSM相比，精度的提升情況。總的來說，本文通過如下指標來評價各個算法的性能：①平均誤差變化，建筑物邊緣區(qū)域的平均誤差在優(yōu)化前后的變化(即優(yōu)化后的平均誤差減去原始平均誤差)；②小于2 m誤差百分比變化(bad-2)，統(tǒng)計誤差小于2 m的像素所占百分比在優(yōu)化前后的變化(即優(yōu)化后的百分比減去原始的百分比)；③小于3 m誤差百分比變化(bad-3)，統(tǒng)計誤差小于3 m的像素所占百分比在優(yōu)化前后的變化。采用這3種精度評價準則，分別對PBR和本文算法進行性能評定，結(jié)果如圖12所示。在圖12中，平均誤差變化越小表示優(yōu)化后的精度提升越大；而bad-2，bad-3越大，表示優(yōu)化后正確點的百分比提升越大。

圖11 衛(wèi)星數(shù)據(jù)集與對應(yīng)的數(shù)字表面模型Fig.11 Satellite dataset and the corresponding DSM

由圖12可以看出，本文算法的優(yōu)化效果明顯優(yōu)于PBR算法(除了(c)圖中的像對2)。在平均誤差變化測度中，PBR算法在部分像對的平均誤差變化大于0，說明PBR算法反而降低了原始DSM的精度，這是因為建筑物邊緣附近經(jīng)常存在其他地物(如樹木)，不符合平面約束條件，因此基于平面約束的PBR算法在這種情況下會降低DSM精度。而本文算法采用了更加靈活的約束，依靠相鄰且灰度相似的像素來優(yōu)化建筑物邊緣，從而能保留建筑物附近的不相似地物，最終提高DSM的精度。在大部分像對中，PBR算法的bad-2和bad-3兩個測度小于0，說明PBR算法反而減少了建筑物邊緣附近的正確點數(shù)目。而本文算法在大部分像對中的bad-2測度和bad-3測度大于0，說明經(jīng)過本文算法優(yōu)化，建筑物邊緣附近的正確點數(shù)目得到提升。因此，綜合高程誤差和正確點百分比兩個方面，均證明本文算法能夠提高衛(wèi)星DSM中的建筑物邊緣精度。為了更全面地評估本文算法的性能，本文列出一些優(yōu)化前后的建筑物邊緣對比，如圖13所示。

圖12 衛(wèi)星數(shù)據(jù)集的性能對比結(jié)果Fig.12 Performance comparison results on satellite dataset

在圖13中，正射影像中的紅色直線表示提取的邊緣直線，與DSM中的白色直線是一一對應(yīng)的關(guān)系。原始DSM中的紅色像素表示因為誤匹配而導(dǎo)致外擴的邊緣，而被紅色曲線所包圍的像素(如圖(c)所示)，表示因為遮擋區(qū)域高程內(nèi)插而導(dǎo)致的邊緣外擴問題。PBR算法通過強行將邊緣附近的區(qū)域擬合為平面，獲得與正射影像直線特征完全吻合的建筑物邊緣，但是這類過強的約束，同時也會破壞建筑物附近的地形地貌，反而降低了DSM的精度，如圖13(b)和(c)所示。本文算法能夠在優(yōu)化建筑物邊緣的同時，有效保留建筑物附近的不相似地物，如圖13(b)所示，并對遮擋區(qū)域的建筑物邊緣，有著明顯的改善效果(如圖13(c)中的紅色曲線區(qū)域所示)。但是，本文算法的效果依賴于圖像的質(zhì)量，當邊緣誤匹配像素的灰度信息，與周圍正確匹配像素的灰度信息存在較大差異時，本文算法無法對這類邊緣進行優(yōu)化，如圖13(a)中下方的建筑物邊緣所示。此外，盡管采用圖割算法進行全局優(yōu)化，本文算法的時間復(fù)雜度較低。以WorldView-3衛(wèi)星數(shù)據(jù)集為例，在單核CPU @2.60 GHz的配置下，采用本文算法處理一張512×512像素的影像，平均處理時間為0.699 s，平均每條邊緣的處理時間為0.012 s。這是因為本文算法僅僅處理建筑物邊緣附近的像素，而不是處理整張影像，因此，時間復(fù)雜度較低。

圖13中的邊緣優(yōu)化對比是以直線段為基本單元。為了更加全面地評估本文算法的邊緣優(yōu)化效果，本文以圖11中的立體像對4作為測試對象，測試整棟建筑物的直線邊緣優(yōu)化效果。測試中的立體像對DSM和正射影像采用切割前的產(chǎn)品，范圍約為0.5 km2(2221×2223像素)，采用文獻[27]的方法對建筑物進行直線邊緣檢測，部分沒檢測出來的直線，采用人工方法提取，從而形成完整的建筑物輪廓，總計提取出46棟建筑物的邊緣。然后，采用本文算法進行建筑物邊緣的全局優(yōu)化，優(yōu)化結(jié)果如圖14所示。圖14中，圖(a)表示優(yōu)化前后的全局DSM；圖14(b)—(d)表示局部若干棟建筑物的邊緣在優(yōu)化前/后的對比。圖14(b)—(d)的具體地理位置，如圖(a)中的紅框所示。由于本文算法只優(yōu)化建筑物邊緣，因此從整體上看，優(yōu)化前后的兩組全局DSM區(qū)別不大，如圖(a)所示。但是，在局部范圍內(nèi)，建筑物邊緣的優(yōu)化效果較為明顯，如圖14(b)—(d)所示。其中，優(yōu)化前建筑物邊緣有一定程度的外擴。在建筑物和地面之間存在平滑過渡區(qū)域，從而導(dǎo)致邊緣模糊化。經(jīng)過本文算法優(yōu)化后，建筑物在DSM中的邊緣具有較明顯的直線特征，從而提高了建筑物邊緣重建的精度。在單核CPU @2.60 GHz的配置下，采用本文算法處理該區(qū)域的時間為6.23 s，平均每棟建筑物的邊緣優(yōu)化時間為0.14 s。本文算法以直線邊緣為基本單元進行優(yōu)化，即一個個的直線進行優(yōu)化，而非DSM/視差圖的整體處理。因此，本文優(yōu)化方法的時間效率取決于建筑物直線邊緣的提取數(shù)目。考慮到建筑物邊緣及其附近的像素數(shù)目，在整張DSM中所占百分比較小，因此盡管采用了圖割解法，本文優(yōu)化算法的時間效率仍然較高，每棟建筑物的優(yōu)化效率可以達到亞秒級?？紤]到每條直線邊緣單獨處理，因此本文優(yōu)化算法非常容易用于并行框架，從而能夠進一步提高處理效率。

圖13 兩種算法的邊緣優(yōu)化效果對比Fig.13 The boundary refinement comparisons of the two algorithms

圖14 整棟建筑物優(yōu)化前后的邊緣對比效果Fig.14 The comparisons of building boundaries before/after the optimization

3 結(jié) 論

本文提出了一種基于直線特征的建筑物邊緣精化方法，用于提高建筑物邊緣的匹配精度或者高程精度。采用航拍影像數(shù)據(jù)集Vaihingen和Toronto，以及WorldView-3衛(wèi)星數(shù)據(jù)集來驗證本文算法的正確性和有效性。試驗結(jié)果表明：本文算法明顯優(yōu)于常用的邊緣銳化算子(雙邊濾波和加權(quán)中值)和一種新的基于緩沖區(qū)平面的邊緣優(yōu)化方法，能夠有效改善建筑物邊緣的精度，減少建筑物邊緣的誤匹配像素。因此，本文算法在一些建筑物三維重建場景(如：虛擬現(xiàn)實、智慧城市、城市管理等)，均能得到應(yīng)用。在未來，計劃結(jié)合深度學(xué)習(xí)技術(shù)，研究更精確的影像邊緣直線提取算法，從而進一步提升建筑物邊緣精化的效果。

致謝：感謝ISPRS提供的Vaihingen和Toronto航拍數(shù)據(jù)集；感謝約翰霍普金斯大學(xué)應(yīng)用物理實驗室，IARPA以及IEEE GRSS圖像分析和數(shù)據(jù)融合技術(shù)委員會提供的融合競賽數(shù)據(jù)。