近景區(qū)域網(wǎng)平差的預處理共軛梯度稀疏解法

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

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

近景區(qū)域網(wǎng)平差的預處理共軛梯度稀疏解法

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

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

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

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

0 引言

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

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

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

1 PCG算法

1.1 共軛梯度算法原理

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

1.2 算法收斂性

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

(1)

式中Cm為次數(shù)為m的Chebyshev多項式。

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

2 快速穩(wěn)定解算策略

2.1 預處理

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

Nx=W,

(2)

選擇對稱正定矩陣P，轉(zhuǎn)化為等價方程

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)較小。在計算數(shù)學中，廣泛采用矩陣的不完全Cholesky分解、不完全LU分解等作為預處理矩陣[9]；但如果矩陣規(guī)模較大，這些方法反而會增加運算負擔。本文根據(jù)區(qū)域網(wǎng)法方程系數(shù)矩陣的特點，選擇主對角線元素占優(yōu)、與N對應的對角平方根矩陣作為預處理矩陣，即

P=(diagN)1/2。

(5)

2.2 稀疏存儲

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

2.3 PCG算法實現(xiàn)

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

初始解及預處理矩陣:x0,P diagN初始殘差:r0 (Nx0-W)解算:Py0=r0,得到y(tǒng)0初始下降方向:P0=-y0,迭代次數(shù)k 0循環(huán)開始 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循環(huán)結束

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

Fig.1 PCG iterative algorithm

3 實驗與分析

3.1 稀疏解算效率對比

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

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

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

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

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

3.2 PCG算法實驗

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

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

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

4 結論

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

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

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

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，F(xiàn)uxin123009,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

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

TP 751.1; P 234.1

A

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

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