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

鞏丹超，韓軼龍，黃 旭

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

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

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

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

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

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

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

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

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

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

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

1.1 建筑物邊緣檢測(cè)

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

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

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

提取出建筑物邊緣的直線特征之后，本文分別對(duì)每一條邊緣直線特征進(jìn)行優(yōu)化。本文算法的核心思想是改變邊緣直線附近像素的視差/高程，從而使得改正后的視差邊緣/高程邊緣與影像的邊緣直線盡量接近。因此，首先需要尋找邊緣直線附近的像素。本文采用“內(nèi)/外緩沖區(qū)”的思路，來統(tǒng)計(jì)邊緣直線附近的像素，如圖4所示。白色直線表示建筑物的邊緣直線。紅色方框表示內(nèi)緩沖區(qū)。本文算法通過改變內(nèi)緩沖區(qū)內(nèi)所有像素的視差/高程，使得視差/高程邊緣與影像邊緣盡量接近。藍(lán)色方框表示外緩沖區(qū)。外緩沖區(qū)內(nèi)像素的視差/高程是不做任何變化的。外緩沖區(qū)的作用有兩個(gè)：①本文采用固定窗口來計(jì)算能量函數(shù)的代價(jià)項(xiàng)，因此，在計(jì)算代價(jià)項(xiàng)的過程中，會(huì)有一部分代價(jià)窗口超出內(nèi)緩沖區(qū)的范圍。超出的部分可以采用外緩沖區(qū)來計(jì)算。具體代價(jià)計(jì)算詳見1.2.1小節(jié)。②傳統(tǒng)的圖割算法只能提供整型最優(yōu)解，因此，在內(nèi)緩沖區(qū)和外緩沖區(qū)之間，存在明顯的視差/高程斷裂。在后續(xù)連續(xù)最優(yōu)解計(jì)算過程中，將考慮外緩沖區(qū)的像素，使得內(nèi)外緩沖區(qū)之間的邊界平滑過渡。外緩沖區(qū)的半徑rout設(shè)置為代價(jià)窗口的一半。考慮到立體影像密集匹配窗口通常采用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)計(jì)邊緣直線附近的像素之后，本文通過改變內(nèi)緩沖區(qū)中像素的視差/高程，來達(dá)到邊緣優(yōu)化的目的?？傮w來說，本文的邊緣優(yōu)化準(zhǔn)則包含兩個(gè)方面：①優(yōu)化后的視差/高程，與原始視差/高程盡量接近；②相鄰且灰度相近的像素，其視差也是相似的。根據(jù)上述兩個(gè)準(zhǔn)則，可以將視差邊緣優(yōu)化問題，轉(zhuǎn)化為全局能量函數(shù)的最優(yōu)解計(jì)算問題，如式(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ū)中所有像素的估計(jì)視差/高程的集合；E(D)表示全局能量函數(shù)的目標(biāo)函數(shù)值；Bin表示內(nèi)緩沖區(qū)所有像素的集合；p表示內(nèi)緩沖區(qū)中的任意一個(gè)像素；Cost(p,d)表示像素p對(duì)應(yīng)于視差/高程d的代價(jià)，是根據(jù)估計(jì)視差/高程與原始視差/高程之間的距離來計(jì)算，詳見1.2.1小節(jié)介紹；Sc(p)表示像素p的置信度，表示像素p的代價(jià)計(jì)算的可靠性，用于降低誤匹配點(diǎn)的代價(jià)權(quán)重，具體詳見1.2.1小節(jié)介紹；N表示所有鄰域像素的集合；q表示p的鄰域像素；dp、dq表示像素p、q的估計(jì)視差/高程；σ表示截?cái)嘀担糜诙x視差階躍條件的閾值；P表示懲罰項(xiàng)系數(shù)；wc(p,q)表示根據(jù)像素p、q的灰度差，計(jì)算出來的權(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項(xiàng)稱為代價(jià)項(xiàng)，將第2項(xiàng)稱為正則項(xiàng)。代價(jià)項(xiàng)和正則項(xiàng)的具體介紹，詳見1.2.1小節(jié)和1.2.2小節(jié)。

1.2.1 代價(jià)項(xiàng)計(jì)算

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

(2)

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

(3)

1.2.2 正則項(xiàng)計(jì)算

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

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

(4)

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

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

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

(5)

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

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

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

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

(6)

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

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

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

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

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

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

perbad=|Serror|/|Serror+Scorrect|

(7)

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

σlow=

(8)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié) 論

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

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