近景區域網平差的預處理共軛梯度稀疏解法

徐振亮， 李艷煥， 閆利， 晏磊

(1.北京大學空間信息集成與3S工程應用北京市重點實驗室，北京 100871；2.遼寧工程技術大學審計處，阜新 123009；3.武漢大學測繪學院，武漢 430079)

近景區域網平差的預處理共軛梯度稀疏解法

徐振亮1， 李艷煥2， 閆利3， 晏磊1

(1.北京大學空間信息集成與3S工程應用北京市重點實驗室，北京 100871；2.遼寧工程技術大學審計處，阜新 123009；3.武漢大學測繪學院，武漢 430079)

針對大規模、近病態法的近景區域網平差法方程快速解算問題，提出基于預處理共軛梯度(preconditioned conjugate gradient，PCG)法的稀疏解算方法。首先，通過選擇與法方程系數矩陣對應的對角平方根矩陣作為預處理矩陣，以改變待估參數向量的坐標基，進而改善法方程系數矩陣性態，達到利用PCG提高收斂速度和解算精度的目的；然后，通過應用稀疏矩陣提高平差法方程系數矩陣的儲存與求解效率。實驗結果證明，該方法不影響攝影測量中區域網平差中多類、多尺度參數同時解算的收斂域，不但具有很高的解算精度，而且速度較快。

預處理；稀疏矩陣；預處理共軛梯度(PCG)；空中三角測量；光束法平差；從運動到結構

0 引言

當前，數字攝影測量技術正在穩步邁進大數據處理時代，比如在空三處理中的光束法平差環節，影像外方位元素眾多，并且多數情況下定向參數之間存在很強的相關性，可能會導致法方程系數矩陣出現病態，為方程的快速、高精度解算帶來很大困難。

目前，對光束法平差技術的研究主要集中在2個方面： ①如何提高大規模平差速度。例如， Lourakis等[1]在稀疏光束法平差(sparse bundle adjustment，SBA)中使用稀疏存儲與矩陣分解技術來求解法方程，降低了計算機對矩陣的存儲與求解負擔。Cornou等[2]提出了減少優化的未知數個數。Bartoli[3]提出將非線性轉化為準線性求解以簡化求解模型，該方法基于計算機視覺，準確地說不屬于嚴格解算模型，解算精度不高。馮其強等[4]則提出了對三維點逐點解算、對攝像機逐個解算的點松弛法快速計算方法，避免了組建整體誤差方程和大型矩陣運算，明顯提高了計算速度；但這種分塊求解對物方點噪聲、攝像機外方位元素噪聲和像方點噪聲的魯棒性不強。②如何提高解算質量。如在針對系數矩陣近病態問題處理方面，對測量領域嶺估計(計算機視覺領域一般采用L-M)方法進行改進，但這種方法的解算結果是有偏估計，破壞了解的統計特性。Olsson等[5]和Kahl等[6]通過邊界約束非凸二次規劃方法，使光束法平差能夠全局收斂；但該方法沒有一般性，在求解過程中依然存在解算效率低的問題。在實際處理中，廣泛應用的仍是經典的光束法平差，即逐次分塊約化法求解[7]。由于在近景攝影測量中攝影方式靈活多變，導致構成的法方程系數矩陣不具有規則帶狀結構，因此難以沿用經典平差方法求解。方程組超線性迭代解算技術仍然作為一種有效方法被廣泛研究應用。

針對上述情況，本文采用預處理共軛梯度(preconditioned conjugate gradient，PCG)稀疏解法解決大規模區域網平差的快速、高質量解算問題。先根據法方程系數矩陣的特點進行預處理，再利用共軛梯度法整體稀疏解算。實驗結果證明，該算法不僅可以改善因法方程系數矩陣病態導致的解算不穩定問題，而且可以保證解算的精度及效率。

1 PCG算法

1.1 共軛梯度算法原理

共軛梯度(conjugate gradient，CG)算法是在每一迭代步利用當前處的最速下降方向生成凸二次函數f的Hesse矩陣G(相當于函數f的二階導數)的共軛方向、并求解f在Rn上的極小點的方法。由于該方法與最速下降法和牛頓法相比避免了二次導數的計算，因而降低了計算量和存儲量，運算效率大為提升。

1.2 算法收斂性

‖x-x(m)‖A≤‖x-x(0)‖A[Cm(1+2η)]-1,

(1)

式中Cm為次數為m的Chebyshev多項式。

由該定義可知，CG法的收斂速度與系數矩陣特征值的分辨率有關。對于n階矩陣A，在條件最壞的情況下，理論上CG法也可以在n步之內得到精確解。如果A的特征值分辨率差異不大，那么CG法可以穩定收斂，并且可以提高解的質量。預處理的目的就在于此。

2 快速穩定解算策略

2.1 預處理

近景影像區域網光束法平差法方程系數矩陣N為一大規模稀疏矩陣，并且條件數k(N)較大，即矩陣特征值尺度差異較大，若直接利用CG法會使解算效率低下；因此需要構造預處理矩陣P以降低原N的條件數，改善特征值分辨率，從而減少迭代次數，達到提高解算穩定的目的[8]。對于光束法平差的法方程

Nx=W,

(2)

選擇對稱正定矩陣P，轉化為等價方程

P-1Nx=P-1W;

(3)

若令N0=P1/2NP-1/2,W0=P1/2W,y=P1/2x,則方程

N0y=W0

(4)

與方程(2)等價。式(2)-(4)中的變量是正則化后的向量。

構造的預處理矩陣P應滿足4項要求： ①對稱正定；②與N的稀疏性相近；③易于求解；④P-1N的特征值分布比較集中，使k(P-1N)較小。在計算數學中，廣泛采用矩陣的不完全Cholesky分解、不完全LU分解等作為預處理矩陣[9]；但如果矩陣規模較大，這些方法反而會增加運算負擔。本文根據區域網法方程系數矩陣的特點，選擇主對角線元素占優、與N對應的對角平方根矩陣作為預處理矩陣，即

P=(diagN)1/2。

(5)

2.2 稀疏存儲

對于稀疏狀矩陣的處理，目前主要采用近似最小度排列技術(approximate minimum degree，AMD)，詳見參考文獻[10-11]。

2.3 PCG算法實現

經過預處理后的法方程即可按照共軛梯度算法解算，圖1給出整個算法實現過程的偽代碼。

初始解及預處理矩陣:x0,P diagN初始殘差:r0 (Nx0-W)解算:Py0=r0,得到y0初始下降方向:P0=-y0,迭代次數k 0循環開始 rk≠0αk rTkykPTkNpkxk+1 (xk+αkPk)rk+1 (rk+αkNPk)解算 Nyk+1 rk+1βk+1 rTk+1yk+1rTkykβk+1 (-yk+1+βk+1Pk)k k+1循環結束

圖1 預處理共軛梯度迭代算法

Fig.1 PCG iterative algorithm

3 實驗與分析

3.1 稀疏解算效率對比

本文選用的法方程大小分別為45階、1 437階及15 945階3種規模，在系數矩陣分別為滿矩陣及稀疏矩陣的2種情況下，采用CG解法進行測試。測試環境為Win7操作系統，4核Intel(R)Core(TM)i3 CPU 處理器，內存為2G。處理程序為自行編制的一套基于Matlab平臺的數字近景攝影測量光束法平差(DCBA3)代碼。解算效率對比情況見表1。

表1 稀疏解算效率比較Tab.1 Comparison of sparse solving efficiency

對于樣例中的54景近景影像數據，最低3°重疊，最高23°重疊，總共5 207個物方點，待估參數為54×6+5 207×3=15 945個。因此，法方程階數高且為稀疏狀(圖2)。

圖2 法方程系數矩陣非零元素分布Fig.2 Nonzero values of normal equation coefficient matrix

在15 945個待估參數的法方程系數矩陣(15 945×15 945)中，絕大多數為0元素，非零元素為932 715個，占整個矩陣空間的0.37%。測試結果表明，采用稀疏矩陣解算確實比滿矩陣效果明顯，在矩陣規模較小情況下，解算效率可提高一倍多。對于大規模法方程來說，采用滿矩陣方法會受限于計算機配置而無法解算，而對于稀疏解算來說是可以完成的。因此，在處理區域網平差問題時，采用稀疏矩陣可明顯提高解算效率。

3.2 PCG算法實驗

仍采用DCBA3程序，分別對上述數據應用CG與PCG方法進行測試，測試結果見表2。

表2 預處理共軛梯度稀疏解算實驗Tab.2 PCG sparse solving experiment

從表2可以看出，2種方法平差后殘差都有所減少，說明CG法和PCG法平差結果是有效的；另外，PCG法解算精度明顯優于CG法，說明構造的預處理矩陣達到了改善法方程系數矩陣特征值分布的目的，提高了解算的穩定性與精確性。從處理的時間上來看， CG法略短于PCG法，但差別不明顯。

4 結論

1)對于大規模近景攝影測量空三平差解算，必須采用稀疏解算技術。

2)與傳統迭代方法比較，由于CG法及PCG法省略了二次導數的計算，因而降低了計算量和存儲量。光束法法方程的共軛梯度解算優勢明顯。

3)通過選擇法方程系數矩陣對應的對角矩陣平方根來構造預處理矩陣，能夠發揮降低矩陣條件數、改善特征值分布的優勢，提高處理速度和精度；并且，該處理方法不改變多類、多尺度參數收斂域，能夠保證方程快速、穩定收斂。

4)PCG方法用于近景攝影測量空三平差處理是可行的。

[1] Lourakis M,Argyros A.The Design and Implementation of a Generic Sparse Bundle Adjustment Software Package Based on the Levenberg Marquardt Algorithm[R].ICS/FORTH Technical Report,No340,2004.

[2] Cornou S,Dhome M,Sayd P,et al.Bundle Adjustment:A Fast Method with Weak Initialisation[G].Cardiff:BMVC,2002:223-232.

[3] Bartoli A.A unified framework for quasi-linear bundle adjustment[C]//Proceedings of the 16th International Conference on Pattern Recognition.Quebec City,Quebec,Canada:IEEE,2002,2:560-563.

[4] 馮其強,李廣云,李宗春.基于點松弛法的自檢校光束法平差快速計算[J].測繪科學技術學報,2008,25(4):300-302. Feng Q Q,Li G Y,Li Z C,et al.Speedy calculation of self calibration bundle adjustment in digital industrial photogrammetry[J].Journal of Geomatics Science and Technology,2008,25(4):300-302.

[5] Olsson C,Kahl F,Oskarsson M.Optimal estimation of perspective camera pose[J].International Conference on Pattern Recongnition,2006,2:5-8.

[6] Kahl F,Henrion D.Globally optimal estimates for geometric reconstruction problems[J].International Journal of Computer Vision,2007,74(1):3-15.

[7] 朱肇光.攝影測量學[M].2版.北京:測繪出版社,1995. Zhu Z G.Photogrammetry[M].2nd ed.Beijing:Surveying and Mapping Press,1995.

[8] 徐振亮.軸角描述的車載序列街景影像空中三角測量與三維重建方法研究[D].武漢：武漢大學，2014. Xu Z L.Research on Aerial Triangulation Angle/Axis Representation and 3D Reconstruction for Vehicle-borne Street-level Image Sequence[D].Wuhan：Wuhan University,2014.

[9] 吳建平,王正華,李曉梅.稀疏線性方程組的高效求解與并行計算[M].長沙:湖南科學技術出版社,2004. Wu J P,Wang Z H,Li X M.Efficient Solving Sparse Linear Equations with Parallel Computing[M].Changsha:Hunan Science and Technology Press,2004.

[10]Davis T.Direct Methods for Sparse Linear Systems[M].Philadelphia:SIAM,2006.

[11]Dellaert F,Kaess M.Square root SAM:Simultaneous localization and mapping via square root information smoothing[J].International Journal of Robotics Research,2006,25(12):1181-1204.

(責任編輯： 劉心季)

PCG sparse algorithm for close-range block bundle adjustment

XU Zhenliang1， LI Yanhuan2， YAN Li3， YAN Lei1

(1.SpatialInformationIntegrationandApplicationsBeijingKeyLaboratory,PekingUniversity,Beijing100871,China；2.AuditOffice,LiaoningTechnicalUniversity，Fuxin123009,China；3.SchoolofGeodesyandGeomatics,WuhanUniversity,Wuhan430079,China)

Aimed at tackling the fast solver problem for the large-scale and nearly pathological close-range block sparse bundle adjustment normal equation，the authors propose a solution method based on the preconditioned conjugate gradient(PCG)sparse algorithm. Firstly，the normal equation coefficient matrix corresponding to diagonal matrix square root is selected as the preconditioning matrix by changing the coordinate base of the parameter vector to be estimated, which can improve the behavior of the normal equation coefficient matrix so as to achieve the purpose of improving the convergence rate of the conjugate gradient method. Then, through the application of the sparse matrix, the efficiency of storage can be improved and the adjustment of normal equation coefficient matrix can be achieved. Experiment results show that the method proposed in this paper has the advantage that any scalar change in variables has no effect on the range of convergence of the iterative technique，and hence it has not only high accuracy of calculation but also faster speed.

preconditioning；sparse matrix；preconditioned conjugate gradient(PCG)；aerial triangulation；bundle adjustment；structure from motion

2013-11-21；

2014-02-08

國家自然科學基金項目“高分測繪衛星動態成像質量與非平穩像移的多參數耦合機理研究”( 編號： 41271456)、“內視場拼接相機的數字基高比模型與精度評價機理”(編號： 11174017)和北京市共建項目“極坐標體系下的局部優化實時定位技術和三維重建”(編號： SYS1000010402)共同資助。

TP 751.1; P 234.1

A

1001-070X(2015)01-0044-04

徐振亮(1982-),男，博士后,講師,現主要從事近景攝影測量與計算機視覺方面的科研工作。Email：xuzhenliang@pku.edu.cn。