大規(guī)模方程組參數(shù)估計與精度評定的改進算法*

陳秋杰 張興福 沈云中

1)同濟大學測繪與地理信息學院，上海 200092

2)同濟大學空間信息與可持續(xù)發(fā)展應用中心，上海 200092

3)廣東工業(yè)大學測繪工程系，廣州510006

1 引言

測量數(shù)據(jù)處理常采用最小二乘法進行參數(shù)估計與精度評定。有關(guān)最小二乘估計的研究非常多，但絕大部分僅關(guān)心如何實現(xiàn)參數(shù)估計及其精度評定，很少關(guān)心其計算效率。大型科學計算問題的待估參數(shù)數(shù)量巨大，相應的法方程矩陣十分龐大，為了提高法方程的解算效率，本文主要研究了:1)利用矩陣分解技術(shù)［1-3］提高法方程解算效率;2)利用多核CPU 并行計算。針對大規(guī)模法方程參數(shù)估計與精度評定，若按照傳統(tǒng)算法(算法1)，需要對法矩陣求逆，工作繁重。基于間接平差模型，使用Cholesky 分解技術(shù)，導出參數(shù)估計與精度評定串行算法(算法2);基于此，引入OpenMP［4］并行運算庫，對算法2并行化(算法3)。并以不同階次的地球重力場反演為例，對3 種算法的效率進行比較。

2 理論與方法

2.1 間接平差模型

對于間接平差模型，有

其中，l 為觀測向量，B 為設計矩陣，x 為待估參數(shù)，D為觀測值先驗信息。法方程及待估參數(shù)精度為:

其中，Nbb=BTPB，w=BTPl。

2.2 傳統(tǒng)的LS 參數(shù)估計與精度評定方法

對于間接平差模型，需要求解式(1)法方程并給出精度評定，傳統(tǒng)的方法需要采用相關(guān)算法對法矩陣求逆:

1)求解Nbb逆矩陣

顯然，算法1(傳統(tǒng)LS 估計)無論是進行參數(shù)求解還是精度評定，都需要求得法矩陣所對應的逆矩陣。但是，矩陣求逆其時間復雜度為O(n^3)，在高維線性方程組中時間消耗很大，然而使用Cholesky分解算法，不僅可避開求逆運算，同時能夠完成參數(shù)估計和精度評定。因此，本文將Cholesky 分解算法引入到法方程參數(shù)求解與精度評定中。

2.3 參數(shù)求解與精度評定

1)參數(shù)求解

式(1)中法方程系數(shù)矩陣Nbb為對稱正定矩陣，由Cholesky 分解原理，存在唯一的下三角矩陣G，使得:

式(2)即為經(jīng)典的Cholesky 分解算法，其實現(xiàn)過程為［5］:

對i=1，2，…，n 作

使用Cholesky 分解法，可得下三角矩陣G，則

(1)中法方程可分解為:

即:

由式(4)、(5)可得參數(shù)求解算法，

2)精度評定

對間接平差模型，為了對待估參數(shù)進行精度評定，若按傳統(tǒng)的精度評定，勢必涉及到逆矩陣求解問題。由于矩陣求逆效率低，并且精度評定一般只需要逆矩陣對角元，因此可以避開對法矩陣整體求逆，而只求解逆矩陣對角元，具體推導過程如下，

由式(1)、(3)可得:

由式(7)可知，為了求得Qxx，只需求得下三角矩陣G 的逆矩陣G－1，設G－1矩陣為:

式中fi=(ffliff2i… ffni)T(i=1，2，…，n)為G－1中第i 個列向量，且存在

GG－1=G(f1f2… fn)=I=(I1I2… In)

即:

從而有:

式中fi利用式(7)求得，同時由于一般精度評定只需要協(xié)因數(shù)Qxx對角元，將其融入式(7)可得Qxx對角元Qi為:

參數(shù)精度評定算法可寫為:

式(8)反映協(xié)因數(shù)陣第i 個對角元恰好為G－1矩陣第i 個列向量所有元素的平方和。由于僅求解協(xié)因數(shù)陣對角元，因此可節(jié)省工作量。

2.4 參數(shù)求解與精度評定算法并行化

使用Cholesky 分解進行參數(shù)求解與精度評定，其主要工作量在于大量的循環(huán)運算。因此算法并行化也主要是針對循環(huán)體工作量的分解與CPU 任務分配。本文以兩個N 維列向量點乘(C=A'B)為例簡要描述循環(huán)體的并行化方法，為了方便理解，算法采用偽代碼描述。對于不作并行化處理的循環(huán)體，其形式非常簡單(圖1)。

圖1 串行算法Fig.1 Serial algorithm

對C=A'B 的并行化處理，其本質(zhì)思想在于將循環(huán)體的工作量分成m(m≥1)份，然后將每份工作量分配到CPU 的各個核心。根據(jù)經(jīng)驗，m 最好與CPU 核心個數(shù)相同，若m 值大于CPU 核心數(shù)，CPU核心間將產(chǎn)生累贅的線程切換，帶來多余的時間消耗，影響并行效率。若CPU 核心數(shù)為4 個，CPU 核心id 分別為1、2、3、4，各個CPU 核心計數(shù)累加器為sum(i)，i=1、2、3、4，則C=A'B 其算法如圖2。

圖2 并行算法Fig.2 Parallel algorithm

當每個CPU 核心完成所分配的工作后，需要由CPU 主核匯總各個CPU 核心的計算結(jié)果:

C=sum(1)+sum(2)+sum(3)+sum(4)

3 模擬計算與分析

為了對改進的LS 算法進行效率分析，以動力學法反演地球重力場為例，分別采用算法1、算法2 和算法3 進行參數(shù)估計與精度評定。程序開發(fā)平臺為VC + + 6.0，采 用Intel C/C + + 9.1 編 譯 器 和OpenMP3.0 并行庫，實驗所用的處理器為Intel Core i5-2400 @ 3.10GHz 4 核CPU，內(nèi)存為8GB(DDR2 FB—DIMM)，操作系統(tǒng)為64 位Windows 7。

利用動力學法反演地球重力場過程中，在每個積分弧段結(jié)束后一般都將局部參數(shù)消去，而只保留與位系數(shù)有關(guān)的法方程，并將每一弧段的法方程進行疊加得到求解位系數(shù)的總法方程［6-8］，記為:

由式(9)可求得地球重力場位系數(shù)，內(nèi)符合精度信息為:

式中，Qi為Quu中第i 個對角元。隨著地球重力場反演階次的增加，待估參數(shù)的數(shù)量近似指數(shù)增長。

第一步:為了驗證三種算法的正確性，以模擬的30天衛(wèi)星軌道數(shù)據(jù)和星間距離變率數(shù)據(jù)(相關(guān)參數(shù)為:采樣率10s，積分弧長2 h，衛(wèi)星位置誤差3 cm，衛(wèi)星速度誤差1.0 ×10－4m/s，衛(wèi)星加速度計誤差3.0 ×10－10m/s2，衛(wèi)星間距離變率誤差3.0 ×10－7m/s)，反演120 階次地球重力場，其結(jié)果如圖3 所示。圖3 反映，應用三種不同算法所解算的大地水準面累積誤差均是一致的，且第120 階的大地水準面累積誤差約為10 cm，表明三種算法均是正確的。

圖3 應用三種不同算法反演120 階次地球重力場Fig.3 Earth’s gravity field model recovery up to 120 degree and order with three different kinds of algorithms

第二步:利用第一步中相關(guān)參數(shù)，模擬30天的星間距離變率及衛(wèi)星軌道數(shù)據(jù)，分別使用三種算法完成50、60、70、80、90、100、110 和120 階次的地球重力場模型的參數(shù)估計與精度評定，統(tǒng)計三種算法的時間消耗(表1，算法1 所采用的法矩陣求逆方法為高斯消去法)。表1 反映，對于50、60、70、80、90、100、110 和120 階次的地球重力場模型參數(shù)求解與精度評定，算法2 和算法3 均比算法1 效率高。算法3 對于50、60、70 階次的地球重力場模型參數(shù)解算與精度評定，其效率不如算法2，這是由于算法3采用并行算法，對于50、60、70 階次的地球重力場模型解算，其法矩陣維數(shù)依次為2 597、3 717 和5 037，均比較低，并行算法中線程的啟動存在一部分時間消耗，而對于小規(guī)模計算問題，該時間消耗相對顯得更加突出。因此，對于低維法方程，算法3 的效率反而不如算法2。表1 中，隨著重力場反演階次由80提高到120，法矩陣維數(shù)則由6 557 增長到14 637，算法3 相對于算法1 的效率由508.7% 提高到750.3%，然而算法2 相對于算法1 的效率卻穩(wěn)定在5 倍左右，可見算法3 在大規(guī)模法方程參數(shù)解算和精度評定中，效率十分顯著。

表1 三種不同算法參數(shù)估計與精度評定的時間消耗Tab.1 The time cost of parametric estimation and accuracy assessment of earth’s gravity field model up to 120 degree and order with three different kinds of algorithms

4 結(jié)束語

1)改進后的LS 算法與傳統(tǒng)LS 算法的計算結(jié)果一致，但其計算效率得到顯著提高。

2)引入OpenMP 并行庫的改進算法更適合大規(guī)模法方程參數(shù)求解與精度評定。

3)僅以地球重力場反演中的參數(shù)估計與精度評定為例進行算法效率分析，但該改進的LS 算法同樣適合其他大規(guī)模法方程參數(shù)估計與精度評定。

1 Davis T A and Hager W W.Dynamic supernodes in sparse Cholesky update/downdate and triangular solves， ACM Trans［J］.Math Software.，2009，35(4):1－16.

2 Chen Y，et al.Rajamanickam，ACM Trans.Algorithm 887:CHOLMOD，supernodal sparse Cholesky factorization and update/downdate［J］.Math Software.，2009，35(3):22－36.

3 Davis T A and Hager W W.Row modifications of a sparse Cholesky factorization，SIAM［J］.Journal on Matrix Analysis and Applications，2005，26(3):621－639.

4 鄒賢才，等.OpenMP 并行計算在衛(wèi)星重力數(shù)據(jù)處理中的應用［J］.測繪學報，2010，39(6):636－641.(Zou Xiancai，et al.Application of parallel computing with OpenMP in data processing for satellite gravity［J］.Acta Geodaetica et Cartographica Sinica，2010，39(6):636－641)

5 同濟大學計算數(shù)學教研室.現(xiàn)代數(shù)值計算［M］.上海:人民郵電出版社，2009.(Computational Mathematics Department of Tongji University.Advanced numerical computation［M］.Shanghai:Post ＆ Telecom Press，2009)

6 張興福，沈云中.應用低軌衛(wèi)星跟蹤數(shù)據(jù)反演地球重力場模型［D］.同濟大學，2007.(Zhang Xingfu and Shen Yunzhong.The Earth’s field model recovery on the basis of satellite-to-satellite tracking missions［D］.Tongji University，2007)

7 沈云中.應用CHAMP 衛(wèi)星星歷恢復地球重力場的研究［D］.中國科學院測量與地球物理研究所，2000.(Shen Yunzhong.Study of recovering gravitational potential model from the ephemerides of CHAMP［D］.The Institute of Geodesy and Geophysics，Chinese Academy of Sciences ，2000)

8 王正濤.衛(wèi)星跟蹤衛(wèi)星測量確定地球重力場的理論與方法［D］.武漢大學，2005.(Wang Zhengtao.Theory and methodology of Earth gravity field recovery by satellite to satellite tracking data［D］.Wuhan University，2005)