基于迭代傳播的小基高比立體匹配方法

門朝光，邊繼龍，李香

(哈爾濱工程大學 計算機科學與技術學院，黑龍江 哈爾濱，150001)

立體匹配是計算機視覺領域中的一個熱點問題，它通過匹配2幅或多幅圖像獲得對應點視差，再通過三角測量原理計算出景物的深度，在計算機視覺、立體測繪、航空攝影測量、遙感等領域有著廣泛的應用。目前，雖然基于局部和基于全局的立體匹配方法[1]取得了一定進展，但由于匹配過程是一個病態過程，并且存在著大量的不確定因素，使得立體匹配仍然是一個活躍的研究領域。由于重復紋理、非紋理區域和遮擋的存在使像素點之間的匹配具有很大程度的不確定性，特別當立體像對不是在同一時間獲取時，移動或消失的目標會使匹配變得更加困難。為解決上述因素的影響，在立體視覺領域中出現了另一個分支即小基高比立體匹配[2-8]。基于小基高比的立體匹配是指在攝影基線與攝像機高度比值較小，對目標物體進行拍照獲得立體像對，再利用立體匹配技術獲得對應點視差，最后，根據三角測量原理計算出物體的高程信息。基于小基高比的立體匹配算法中存在 1個相關基本等式[2-4](Central equation of correlation)，此等式把測量視差和真實視差聯系起來，可在匹配結束后對相關基本等式進行求解以恢復其真實視差。通過求解此等式可以解決小基高比立體匹配中存在的邊界膨脹問題即黏合現象。相關基本等式求解問題屬于數學上的反問題，具有不適定性[9]，需要對相關基本等式進行正則化處理，使問題變得適定可解。本文通過在目標泛函中引入全變差項進行正則化處理，把等式求解問題轉化為求解目標泛函的極值問題。目前，數字圖像處理領域中變分法受到了許多研究者的關注，它廣泛應用于圖像復原[10-11]、圖像融合[12]、圖像分割[13]和立體匹配[14]中。基于變分法的優勢在于它能提供一個良好的數學形式體系，而且可以生成亞像素級視差。變分法中的正則化參數選取關系到解的逼近性與數值穩定性，其選取方法可分為先驗和后驗2類策略。其中，先驗策略是在求出正則解之前就已將正則參數確定。這種方法往往具有理論分析價值，在實際中常常難以驗證其賴以施用的條件；后驗策略則是在計算正則解的過程中根據一定的原則確定與原始數據的誤差水平相匹配的正則參數，此方法比較實用。本文采用后驗策略根據Morozov偏差原理推導出正則參數選擇公式，并設計一個迭代過程進行求解。根據Morozov原理求解正則參數需要了解系統的誤差水平，然而，在立體匹配中我們無法事先獲知該誤差水平。若根據經驗設置此誤差水平，則在匹配不同的立體像對時需要選擇不同的經驗值，這會影響算法的適應性，為此本文根據啟發式信息估計該誤差水平。最后，根據延遲擴散定點迭代思想[15](lagged diffusivity fixed point iteration，FP)獲得目標泛函的迭代傳播等式，通過共軛梯度法進行求解。

1 相關性分析

假設u和u?代表立體像對中左右圖像，并滿足經典立體模型：

其中：函數 λ(x)在空間上緩慢變化；視差函數 ε(x)是一有界函數，描述立體像對u和之間的幾何畸變。在獲取立體像對的過程中，若2次成像之間的角度較大即大基高比，則立體像對會受光照變化(例如云層遮擋)、物體遮擋和非朗伯平面等因素的影響而導致它們之間存在較大的輻射差異，致使假設立體模型存在較大的偏差；若2次成像之間的角度較小即小基高比，則可以認為2次成像時的光照條件近似一樣而且小角度可以減少物體遮擋，減弱非朗伯平面對輻射差異的

此立體模型假設立體像對u和u?之間僅存在幾何畸變，不存在輻射差異，這種模型只有在基高比 b/h較小時才能成立，而且基高比越小模型越精確。

基于小基高比的立體匹配方法[2-8]屬于局部匹配方法，該方法的匹配精度受支撐窗口大小影響。當窗口大小增加時匹配精度隨之增加，但同時大窗口也會導致算法在視差不連續處產生誤匹配，在立體匹配中把這種誤匹配稱為黏合現象，這種現象會導致視差圖中的前景區膨脹失真。如圖1所示，左圖像中存在3個點a，b和c，它們分別與右圖像中的′，′和′ 3個點相對應，其視差分別為0，db和dc。在匹配時由于b點的支撐窗口包含了建筑物邊界致使算法在b點產生錯誤的視差db=dc，從而導致前景區的建筑物產生膨脹。影響，在這種條件下生成的立體像對可以近似滿足這一模型。在小基高比條件下，立體模型可以簡化為：

圖1 黏合現象Fig.1 Adhesion phenomenon

假設在小基高比模型(2)中真實視差 ε(x)足夠小，通過測量視差對這個模型進行一階近似，可以發現通過相關基本等式(4)能將最大化式(3)所獲得的測量視差與真實視差聯系起來。

式中：函數 φ ( x0- x )表示以x0為中心的窗口函數；ε(x)表示真實視差函數；0()d x為圖像u在x0點的相關密度函數[4]，其數學表達式為：

相關密度函數只依賴于參考圖像 u，它表達了窗口中心x0周圍的邊緣方向。通過求解相關基本等式可以解決小基高比立體匹配中存在的黏合現象。Delon等[4]通過 delta函數近似式(4)中的相關密度函數將等式(4)簡化為ε(x1)= m(x0)，再通過質心校正公式(6)計算出質心的灰度：

然后，在局部窗口內尋找一個與此灰度x1相似的像素位置作為質心位置并將測量視差賦予此位置。雖然通過簡單的質心校正取得了一定的成效，但根據該校正公式計算的灰度確定質心位置會導致誤差傳播。為此，本文提出通過變分原理求解該等式來解決小基高比立體匹配中存在的黏合現象，該方法可以有效減少匹配中存在的黏合現象，提高視差圖的準確率。

2 模型分析及求解

2.1 模型分析

小基高比立體匹配中存在的相關基本等式(4)可以表示成形如式(7)的第一類算子方程，其數學表達式為：

由于式(7)中的算子K是1個緊算子，因此，導致求解相關基本等式問題成為1個病態問題。為使等式變得適定可解，經常使用Tikhonov正則化方法對等式進行正則化處理。根據Tikhonov正則化方法，將求解等式(7)問題轉化為求解目標泛函極值問題，其目標泛函為：

式中：?ε為ε的梯度，其數學表達式為

圖2 懲罰函數的比較Fig.2 Comparison of different penalty functions

第1種懲罰函數式(9)(如圖2(a)所示)在視差函數ε不連續時經常導致虛假振蕩(spurious oscillations)；第2種懲罰函數式(10)要求視差函數ε處處平滑，然而在物體邊界(即視差不連續處)常常會違背這一假設，因此，這2種懲罰函數都不太適合此應用。為保持函數的不連續性，防止在邊緣處過度平滑，許多研究者引入了全變差[10-12]來解決此問題，其全變差的數學表達式為：

此懲罰函數的形狀如圖2(b)所示。由于此懲罰函數不要求視差函數ε處處平滑，因此該函數是一個很好的候選正則函數。然而該函數的歐式范數在0處不可微，為避免這一問題，正則函數修改為：

其修改后的懲罰函數形狀如圖 2(c)所示。在 β=0時，可以簡化為全變差函數

根據修改的全變差懲罰函數式(13)相應的目標泛函可以表達為：

然后，通過對目標泛函E(ε)求導獲得該泛函的歐拉-拉格朗日方程：

式中：K*代表伴隨算子，其數學表達式為

L(ε)為微分算子，它對函數ω的作用可以通過下式表示：

最后，根據Vogel提出的延遲擴散定點迭代思想最小化目標泛函式(14)。其FP迭代可以表達為：

每次迭代必須通過求解線性傳播式(18)來獲得下一步迭代εk+1，迭代過程的傳播率取決于前一迭代εk。

2.2 數值求解

經過離散化處理之后，式(18)轉化為超大稀疏線性方程組：

由于系數矩陣A為超大稀疏矩陣，直接通過高斯消除法求解線性方程組的方法不可行，因此，本文通過迭代的共軛梯度法來求解線性方程組，其共軛梯度法如下。

算法1(共軛梯度法[16])：

k=0

r0=b-Ax0

while rk≠0

k=k+1

if k=1

p1=r0

else

βk=(rk-1)Trk-1/(rk-2)Trk-2

pk= rk-1+βkpk

end

αk=(rk-1)Trk-1/(pk)TApk

xk=xk-1+αkpk

rk=rk-1-αkApk

end

x=xk

3 正則化處理

3.1 正則參數選擇

由于正則化參數α的選取關系到解的逼近性與數值穩定性，為此，根據Morozov原理[17]選擇此參數以獲得較為精確的視差。Morozov原理闡述了怎樣利用觀測數據選擇正則參數 α，它在線性反問題領域中受到了很大關注。由于正則化產生的誤差等于觀測數據所引起的誤差，因此，可以根據這一原理選擇正則參數α：

在計算正則參數的過程中需要利用以下幾個法則，本文在此給出這些法則和其證明過程。令F(α)是關于正則參數的函數，其數學表達式為：

法則1 對任意α＞0，式(21)存在惟一最小化解，它可以描述成下面系統的解：

或者用內積可表示為：

法則2 F(α)的一階導數和二階導數可以表示為：

證明：

根據式(23)可得：

因此，可得：

法則 3 F(α)是關于正則參數的函數，對任意α＞0，存在如下等式：

證明 根據法則1，有如下等式成立：

對上述等式兩邊關于α求導得：

根據法則2，上述等式可以轉化為：

對等式兩邊積分得：

根據 F(α)和 F′(α)可以將式(20)轉化為：

為了獲得合理的正則參數，根據Morozov原理推導出正則參數選擇方法。為推導模型函數，對式(25)的第3項進行如下近似：

式中：T1是待定常數。通過上述近似，式(25)可以簡化為：

對等式(28)兩邊積分可得：

為簡化正則參數α的計算，令m(0)=0，其模型函數m(α)的表達式可以簡化為：

對模型函數m(α)關于α求導得：

最后，通過迭代方法求解正則參數α，具體算法如下。

(1) 根據式(21)和式(24)計算 F(αk)和 F′(αk)，再根據式(32)計算Ck和Tk。

通過求解式(32)可得：

(2) 根據式(30)和式(31)計算 m(αk)與 m'(αk)。

(3) 根據 Morozov偏差方程(26)計算正則參數αk+1。

3.2 估計誤差水平

計算正則參數時需要使用觀測誤差水平δ，然而，立體匹配中往往無法事先獲知該誤差水平。為了能精確的估計正則參數，在計算正則參數之前首先利用啟發式方法估計該誤差水平。根據Morozov原理有如下不等式關系成立：

利用式(34)可推導出誤差觀測水平的存在區間：

在數值實現過程中，利用模型函數m(α)預測觀測誤差δ：

具體實現過程如下：

(1) 給定比例常數 σ ∈ [ 0 .5,1]和初始正則參數 α0＞0。

(2) 利用式(21)和式(24)計算 F(α0)和 F'(α0)。然后，利用式(37)計算式(30)中的參數C和T。

再利用式(38)計算正則參數α1：

式中：y0為曲線 y=m(α)上點(α0,m(α0))的切線與 y 軸的截距，其數學表達式為

(3) 根據步驟(2)中計算的正則參數α1，利用式(21)和式(24)計算 F(α1)和 F'(α1)，再利用式(25)計算參數C0，然后，利用式(29)和式(40)計算參數C1和T1。

利用式(29)計算m(0)和m(1)，誤差δ表示為：

3.3 算法總體流程

基于迭代傳播的小基高比立體匹配方法是在Delon等[4]提出的基于多分辨率的小基高比立體匹配方法的基礎上通過對匹配中存在的相關基本等式進行迭代求解來恢復其真實視差值。該方法可以減少小基高比立體匹配過程中因匹配窗口跨越邊界而導致的黏合現象，立體匹配中的這種黏合現象會導致視差圖中前景區產生“膨脹”。算法具體流程如下：

(1) 利用小基高比立體匹配算法[4]獲得初始測量視差ε0，并給定比例常數σ和初始正則參數α0。

(2) 利用3.2節中的算法獲得誤差水平δk。

(3) 利用3.1節中的算法獲得正則參數αk。

(4) 利用共軛梯度法求解線性方程組(19)獲得視差 εk+1，若或者 k ≥kmax則停止迭代；否則，k=k+1，轉入步驟(2)繼續迭代。

4 實驗分析

本實驗使用 CNES提供的航空攝影像對Toulouse(如圖3(a)，(b)所示)對本文算法的合理性進行驗證。這幅立體像對是1幅航空影像圖，基高比約為0.045，地面分辨率為 R=0.5，其獲取時間間隔為 20 min，這導致立體像對中存在明顯的運動和陰影移動，增加了視差估計的難度。實際上，理想的小基高比立體像對的2次獲取時間間隔較小，這樣可以有效的減少立體像對間的輻射差異和景物運動對匹配產生的不利影響。

圖3 Toulouse立體像對Fig.3 Toulouse stereo images

理論上講，小基高比條件下形成的立體像對非常相似，具有較小的視差搜索范圍。這幅立體像對的視差范圍為[-2，2]，可近似模擬小基高比條件下的立體像對，這種相似性使匹配過程變得更容易。圖 3(a)和(b)所示為這幅立體像對的左右圖像，看出它們非常相似，具有較小的幾何畸變與視差搜索范圍。其真實視差圖如3(c)所示。圖4(a)為Delon等[4]提出的基于多分辨率的小基高比立體匹配方法獲得的實驗結果，通過與真實視差圖對比可以看出該方法已經獲得了較好的效果，但是，在背景處還是存在一些明顯的誤匹配，而且在物體邊界處伴隨著較多黏合現象(如圖 4(b)所示)。圖4(c)所示為由本文算法生成的視差圖。從圖4可見：由本文算法生成的視差圖中場景各主要組成部分都清晰地檢測出來，特別是建筑物四周只發生微小變化的那些部分也都被清晰地檢測出來。通過對比圖4(b)和 4(d)可以看出：本文算法的黏合現象要明顯少于文獻[4]中的黏合現象，而且該算法生成的視差圖的準確率高達95%以上。

為驗證小基高比立體匹配的高程精度，本實驗使用由空間機電研究所提供的1幅帶有真實高程信息的小基高比立體像對(如圖 5(a)和(b)所示)。該立體像對的基高比約為0.05，地面分辨率為0.2。在該實驗中，根據本文算法生成的視差圖(如圖 5(c)所示)計算幾個主要建筑物的高程信息，這些建筑被標記在圖5(d)中，對應建筑物的視差、計算高程和真實高程如表1所示。高程計算公式如下：式中：δ為視差；B/H為基高比；R為地面分辨率。從表1可知：該算法具有較高的高程精度，其平均高程誤差為0.33 m，像元匹配差異精度為0.082 5，優于1/10個像元。

圖4 黏合現象比較Fig.4 Comparison of adhesion phenomenon

表1 高程信息Table 1 Information of DEM

圖5 高程精度驗證Fig.5 Verification of DEM accuracy

5 結論

(1) 根據啟發式信息估計視差誤差水平，解決了先驗參數估計問題，使匹配算法具有更強的適應性。

(2) 所提出的迭代正則參數選擇方法解決了正則解的穩定性與逼近性問題，使匹配算法獲得了更加精確的視差。

(3) 通過共軛梯度法求解目標泛函的迭代傳播等式加快了匹配算法的收斂速度。

(4) 該方法能有效地減少小基高比立體匹配中的黏合現象，提高視差圖的準確率，其視差圖的準確率可達95%以上，且亞像素精度可優于1/10個像元。

[1] 羅三定, 陳海波. 基于區域增長的自適應窗口立體匹配算法[J]. 中南大學學報: 自然科學版, 2005, 36(6): 1042-1047.LUO San-ding, CHEN Hai-bo. Stereo matching algorithm of adaptive window based on region growing[J]. Journal of Central South University: Science and Technology, 2005, 36(6):1042-1047.

[2] Delon J. Fine comparison of images and other problems[D].Paris: ENS Cachan. Mathematics Department, 2004: 89-108.

[3] Facciolo G. Variational adhesion correction with image based regularization for digital elevation models[D]. Montevicleo:University of the Republic. Engineering College, 2005: 37-49.

[4] Delon J, Rouge B. Small baseline stereovision[J]. Journal of Mathematical Imaging and Vision, 2007, 28(3): 209-223.

[5] Igual L, Preciozzi J, Garrido L, et al. Automatic low baseline stereo in urban areas[J]. Inverse Problems and Imaging, 2007,1(2): 319-348.

[6] Morgan G L K, LIU Jian-guo, YAN Hong-shi. Sub-pixel stereo matching for DEM generation from narrow baseline stereo imagery[C]//Proceedings of International Geoscience and Remote Sensing Symposium. Boston: IEEE, 2008: 1284-1287.

[7] Sabater N, Blanchet G, Moisan L,et al. Review of low-baseline stereo algorithms and benchmarks[C]//Proceedings of Image and Signal Processing for Remote Sensing XVI. Toulouse: IEEE,2010: 1-12.

[8] Morgan G L K, LIU Jian-guo, YAN Hong-shi. Precise subpixel disparity measurement from very narrow baseline stereo[J].IEEE Transactions on Geoscience and Remote Sensing, 2010,48(9): 3424-3433.

[9] 王彥飛. 反問題的計算方法及其應用[M]. 北京: 高等教育出版社, 2007: 3-15.WANG Yan-fei. Computational method for inverse problem and their application[M]. Beijing: High Education Press, 2007: 3-15.

[10] Aujol J F. Some first-order algorithms for total variation based image restoration[J]. Mathematical Imaging and Vision, 2009,34(3): 307-327.

[11] Beck A, Teboulle M. Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[J]. IEEE Transactions on Image Processing, 2009,18(11): 2419-2434.

[12] Mrityunjay K, Sarat D. A total variation-based algorithm for pixel-Level image fusion[J]. IEEE Transactions on Image Processing, 2009, 18(9): 2137-2143.

[13] Chung G, Vese L A. Image segmentation using a multilayer level-set approach[J]. Computing and Visualization in Science,2009, 12(6): 267-285.

[14] Benari R, Sochen N. Stereo matching with Mumford-Shah regularization and occlusion handling[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2010, 32(11):2071-2084.

[15] Vogel C R, Oman M E. Iterative methods for total variation denoising[J]. SIAM Journal on Scientific Computing, 1996,17(1): 227-238.

[16] Golub G H, Loan C F V. Matrix computations[M]. London:Johns Hopkins Univ, 1996: 520-532.

[17] Kunisch K, ZOU Jun. Iterative choices of regularization parameters in linear inverse problems[J]. Inverse Problems, 1998,14: 1247-1264.