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

(1.北京大學(xué)空間信息集成與3S工程應(yīng)用北京市重點(diǎn)實(shí)驗(yàn)室，北京 100871；2.遼寧工程技術(shù)大學(xué)審計(jì)處，阜新 123009；3.武漢大學(xué)測(cè)繪學(xué)院，武漢 430079)

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

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

(1.北京大學(xué)空間信息集成與3S工程應(yīng)用北京市重點(diǎn)實(shí)驗(yàn)室，北京 100871；2.遼寧工程技術(shù)大學(xué)審計(jì)處，阜新 123009；3.武漢大學(xué)測(cè)繪學(xué)院，武漢 430079)

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

預(yù)處理；稀疏矩陣；預(yù)處理共軛梯度(PCG)；空中三角測(cè)量；光束法平差；從運(yùn)動(dòng)到結(jié)構(gòu)

0 引言

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

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

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

1 PCG算法

1.1 共軛梯度算法原理

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

1.2 算法收斂性

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

(1)

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

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

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

2.1 預(yù)處理

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

Nx=W,

(2)

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

P-1Nx=P-1W;

(3)

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

N0y=W0

(4)

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

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

P=(diagN)1/2。

(5)

2.2 稀疏存儲(chǔ)

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

2.3 PCG算法實(shí)現(xiàn)

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

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

圖1 預(yù)處理共軛梯度迭代算法

Fig.1 PCG iterative algorithm

3 實(shí)驗(yàn)與分析

3.1 稀疏解算效率對(duì)比

本文選用的法方程大小分別為45階、1 437階及15 945階3種規(guī)模，在系數(shù)矩陣分別為滿矩陣及稀疏矩陣的2種情況下，采用CG解法進(jìn)行測(cè)試。測(cè)試環(huán)境為Win7操作系統(tǒng)，4核Intel(R)Core(TM)i3 CPU 處理器，內(nèi)存為2G。處理程序?yàn)樽孕芯幹频囊惶谆贛atlab平臺(tái)的數(shù)字近景攝影測(cè)量光束法平差(DCBA3)代碼。解算效率對(duì)比情況見(jiàn)表1。

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

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

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

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

3.2 PCG算法實(shí)驗(yàn)

仍采用DCBA3程序，分別對(duì)上述數(shù)據(jù)應(yīng)用CG與PCG方法進(jìn)行測(cè)試，測(cè)試結(jié)果見(jiàn)表2。

表2 預(yù)處理共軛梯度稀疏解算實(shí)驗(yàn)Tab.2 PCG sparse solving experiment

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

4 結(jié)論

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

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

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

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

[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] 馮其強(qiáng),李廣云,李宗春.基于點(diǎn)松弛法的自檢校光束法平差快速計(jì)算[J].測(cè)繪科學(xué)技術(shù)學(xué)報(bào),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] 朱肇光.攝影測(cè)量學(xué)[M].2版.北京:測(cè)繪出版社,1995. Zhu Z G.Photogrammetry[M].2nd ed.Beijing:Surveying and Mapping Press,1995.

[8] 徐振亮.軸角描述的車載序列街景影像空中三角測(cè)量與三維重建方法研究[D].武漢：武漢大學(xué)，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] 吳建平,王正華,李曉梅.稀疏線性方程組的高效求解與并行計(jì)算[M].長(zhǎng)沙:湖南科學(xué)技術(shù)出版社,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.

(責(zé)任編輯： 劉心季)

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

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

TP 751.1; P 234.1

A

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

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