陳秋杰 張興福 沈云中

1)同濟(jì)大學(xué)測(cè)繪與地理信息學(xué)院，上海 200092

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

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

1 引言

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

2 理論與方法

2.1 間接平差模型

對(duì)于間接平差模型，有

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

其中，Nbb=BTPB，w=BTPl。

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

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

1)求解Nbb逆矩陣

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

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

1)參數(shù)求解

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

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

對(duì)i=1，2，…，n 作

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

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

即:

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

2)精度評(píng)定

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

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

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

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

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

即:

從而有:

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

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

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

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

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

圖1 串行算法Fig.1 Serial algorithm

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

圖2 并行算法Fig.2 Parallel algorithm

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

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

3 模擬計(jì)算與分析

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

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

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

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

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

圖3 應(yīng)用三種不同算法反演120 階次地球重力場(chǎng)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 階次的地球重力場(chǎng)模型的參數(shù)估計(jì)與精度評(píng)定，統(tǒng)計(jì)三種算法的時(shí)間消耗(表1，算法1 所采用的法矩陣求逆方法為高斯消去法)。表1 反映，對(duì)于50、60、70、80、90、100、110 和120 階次的地球重力場(chǎng)模型參數(shù)求解與精度評(píng)定，算法2 和算法3 均比算法1 效率高。算法3 對(duì)于50、60、70 階次的地球重力場(chǎng)模型參數(shù)解算與精度評(píng)定，其效率不如算法2，這是由于算法3采用并行算法，對(duì)于50、60、70 階次的地球重力場(chǎng)模型解算，其法矩陣維數(shù)依次為2 597、3 717 和5 037，均比較低，并行算法中線程的啟動(dòng)存在一部分時(shí)間消耗，而對(duì)于小規(guī)模計(jì)算問題，該時(shí)間消耗相對(duì)顯得更加突出。因此，對(duì)于低維法方程，算法3 的效率反而不如算法2。表1 中，隨著重力場(chǎng)反演階次由80提高到120，法矩陣維數(shù)則由6 557 增長(zhǎng)到14 637，算法3 相對(duì)于算法1 的效率由508.7% 提高到750.3%，然而算法2 相對(duì)于算法1 的效率卻穩(wěn)定在5 倍左右，可見算法3 在大規(guī)模法方程參數(shù)解算和精度評(píng)定中，效率十分顯著。

表1 三種不同算法參數(shù)估計(jì)與精度評(píng)定的時(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)改進(jìn)后的LS 算法與傳統(tǒng)LS 算法的計(jì)算結(jié)果一致，但其計(jì)算效率得到顯著提高。

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

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

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 并行計(jì)算在衛(wèi)星重力數(shù)據(jù)處理中的應(yīng)用［J］.測(cè)繪學(xué)報(bào)，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 同濟(jì)大學(xué)計(jì)算數(shù)學(xué)教研室.現(xiàn)代數(shù)值計(jì)算［M］.上海:人民郵電出版社，2009.(Computational Mathematics Department of Tongji University.Advanced numerical computation［M］.Shanghai:Post ＆ Telecom Press，2009)

6 張興福，沈云中.應(yīng)用低軌衛(wèi)星跟蹤數(shù)據(jù)反演地球重力場(chǎng)模型［D］.同濟(jì)大學(xué)，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 沈云中.應(yīng)用CHAMP 衛(wèi)星星歷恢復(fù)地球重力場(chǎng)的研究［D］.中國科學(xué)院測(cè)量與地球物理研究所，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)星測(cè)量確定地球重力場(chǎng)的理論與方法［D］.武漢大學(xué)，2005.(Wang Zhengtao.Theory and methodology of Earth gravity field recovery by satellite to satellite tracking data［D］.Wuhan University，2005)