海量衛(wèi)星重力梯度觀(guān)測(cè)數(shù)據(jù)確定地球重力位模型的數(shù)值方法*

朱廣彬 常曉濤 鄒賢才 徐新禹 王建強(qiáng)

(1)國(guó)家測(cè)繪局衛(wèi)星測(cè)繪應(yīng)用中心，北京 100830 2)武漢大學(xué)測(cè)繪學(xué)院，武漢 430079 3)東華理工大學(xué)測(cè)繪工程學(xué)院，撫州344000)

海量衛(wèi)星重力梯度觀(guān)測(cè)數(shù)據(jù)確定地球重力位模型的數(shù)值方法*

朱廣彬1)常曉濤1)鄒賢才2)徐新禹2)王建強(qiáng)3)

(1)國(guó)家測(cè)繪局衛(wèi)星測(cè)繪應(yīng)用中心，北京 100830 2)武漢大學(xué)測(cè)繪學(xué)院，武漢 430079 3)東華理工大學(xué)測(cè)繪工程學(xué)院，撫州344000)

基于空域最小二乘法，對(duì)衛(wèi)星重力梯度數(shù)據(jù)確定地球重力場(chǎng)中的Cholesky分解法、預(yù)條件共軛梯度法以及OpenMP并行算法3種數(shù)值方法進(jìn)行比較與分析。研究表明，在計(jì)算機(jī)硬件資源有限的情況下，傳統(tǒng)的Cholesky分解法已經(jīng)無(wú)法滿(mǎn)足求解要求;預(yù)條件共軛梯度法的求解效率較之Cholesky分解法有改進(jìn)，但其以損失小量精度為代價(jià);OpenMP并行算法在不損失求解精度的條件下，可提高求解的效率。

衛(wèi)星重力梯度;Cholesky分解;預(yù)條件共軛梯度;OpenMP并行算法;數(shù)據(jù)處理

1 引言

地球重力位模型一般采用球諧系數(shù)進(jìn)行表達(dá)。基于空域最小二乘法或時(shí)域最小二乘法［1-5］，利用衛(wèi)星重力梯度數(shù)據(jù)確定地球重力場(chǎng)時(shí)，隨著求解階數(shù)的增加，未知數(shù)個(gè)數(shù)呈平方倍增長(zhǎng)。因此，如何對(duì)海量數(shù)據(jù)進(jìn)行快速有效的數(shù)值處理是需要重點(diǎn)考慮的內(nèi)容。一般來(lái)講，有以下3種處理方式：方程組的直接解法，例如Cholesky分解法、矩陣QR分解法等。此種方法經(jīng)過(guò)有限次運(yùn)算能夠獲得方程的精確解，但是總的運(yùn)算量和運(yùn)算時(shí)間并沒(méi)有得到優(yōu)減，故而該方法在求解大型方程組時(shí)實(shí)用性較差。第二種是方程組的迭代解法，即通過(guò)構(gòu)造一個(gè)向量的無(wú)窮序列，其極限對(duì)應(yīng)了方程組的精確解，從某一初始解開(kāi)始，通過(guò)逐次迭代，逐漸收斂至真解。該類(lèi)方法在有限迭代次數(shù)內(nèi)，無(wú)法得到精確解和精度信息，但運(yùn)算時(shí)間大大縮減。共軛梯度法、最速下降法、預(yù)條件共軛梯度法(PCCG)等均屬于該類(lèi)方法的范疇。第三種方法是基于并行計(jì)算思想，充分發(fā)揮計(jì)算機(jī)硬件結(jié)構(gòu)特點(diǎn)，對(duì)直接解法的計(jì)算結(jié)構(gòu)進(jìn)行優(yōu)化組合，從而提高程序的執(zhí)行效率。本文將對(duì)Cholesky分解法、預(yù)條件共軛梯度法、OpenMP(Open Multi-Processing)并行解法進(jìn)行數(shù)值分析和比較。

2 方法介紹

2.1 Cholesky分解法

基于空域最小二乘法或者時(shí)域最小二乘法，利用衛(wèi)星重力梯度數(shù)據(jù)恢復(fù)地球重力場(chǎng)時(shí)，均可以化為［3-6］：

其中，法方程N(yùn)=ATPA滿(mǎn)足正定對(duì)稱(chēng)特性，W= ATPl，A為設(shè)計(jì)矩陣，P為權(quán)陣，l為觀(guān)測(cè)向量，因此可以利用Cholesky分解法進(jìn)行直接求解。

2.2 PCCG方法

共軛梯度法的基本思想是利用迭代求得的增量改正上一次迭代的初始向量，所得結(jié)果作為下次迭代的初值。增量方向取逼近其解的最速化方向，且利用上一次的初始向量和增量按此梯度方向更新增量。為了進(jìn)一步提高迭代的收斂速度，可以選擇一個(gè)預(yù)條件陣Nbd應(yīng)用到共軛梯度方法中，該矩陣的逆要求便于計(jì)算，且與法方程矩陣的乘積的條件數(shù)要小于法方程矩陣自身的條件數(shù)，即：

PCCG方法的具體計(jì)算步驟如下［7］：

1)選擇參數(shù)初值x0==0，計(jì)算殘差向量與方向向量的初值為預(yù)條件矩陣逆陣，k=0;

3)計(jì)算新的參數(shù)向量xk+1=xk+αkpk;

4)計(jì)算新的殘差向量rk+1=rk-αkak;

8)k=k+1，回到迭代步驟2)。

前后兩次迭代解向量的差利用

2.3 OpenMP并行編程模型

OpenMP采用Fork-Join并行執(zhí)行模型［8］。當(dāng)程序開(kāi)始執(zhí)行時(shí)只有主線(xiàn)程存在，主線(xiàn)程會(huì)一直串行的執(zhí)行;當(dāng)遇到需要并行計(jì)算的并行域時(shí)，主線(xiàn)程會(huì)創(chuàng)建一隊(duì)并行的線(xiàn)程開(kāi)始并行執(zhí)行;當(dāng)派生線(xiàn)程在并行域中執(zhí)行完畢后，它們退出或者掛起，此時(shí)只有主線(xiàn)程在運(yùn)行(圖1)。OpenMP應(yīng)用編程接口可以根據(jù)不同并行域的需要?jiǎng)討B(tài)地改變線(xiàn)程數(shù)，且該并行結(jié)構(gòu)可以嵌入到別的并行結(jié)構(gòu)中去。從本質(zhì)上講，OpenMP的并行化是通過(guò)使用嵌入到源代碼中的編譯制導(dǎo)語(yǔ)句來(lái)實(shí)現(xiàn)的，其是一個(gè)外部編程模型，而不是自動(dòng)的編程模型。

圖1 OpenMP并行執(zhí)行模式Fig.1 OpenMP parallel execution mode

并行算法的加速比與效率是并行算法的兩個(gè)重要指標(biāo)，其表征了使用并行算法求解實(shí)際問(wèn)題所能獲得的好處［9］。并行算法的絕對(duì)加速比定義為：

式中，t1(n)為求解一規(guī)模為n的問(wèn)題的最優(yōu)串行算法在單處理器上的運(yùn)行時(shí)間，tp(n)為求解該問(wèn)題的并行算法在p個(gè)處理器上的運(yùn)行時(shí)間。最優(yōu)的串行算法很難找到，甚至不存在，一般利用實(shí)際所用的串行算法來(lái)代替。

并行算法的效率定義為：

式中，p為并行算法運(yùn)行所需的處理器數(shù)。

如果并行算法的加速比與處理器個(gè)數(shù)成正比，則稱(chēng)該并行算法具有線(xiàn)性加速比;如果Ep(n)＞1，則稱(chēng)之為超線(xiàn)性加速比。

3 數(shù)值分析

3.1 計(jì)算效率評(píng)測(cè)

由于利用重力梯度數(shù)據(jù)求解地球重力場(chǎng)是一個(gè)超大型的最小二乘問(wèn)題。這對(duì)計(jì)算機(jī)硬件的要求較高。表1列出了30天，5 s采樣的重力梯度觀(guān)測(cè)數(shù)據(jù)在不同階數(shù)下所對(duì)應(yīng)的計(jì)算機(jī)資源需求。

從表1可以看出，對(duì)于普通的小型服務(wù)器(本文計(jì)算均基于小型服務(wù)器P575-2進(jìn)行，具體參數(shù)見(jiàn)表2)而言，設(shè)計(jì)矩陣的存儲(chǔ)基本上是不可能的，只能通過(guò)分塊策略組成法方程。相對(duì)而言，預(yù)條件矩陣為一稀疏陣，采用CSR格式存儲(chǔ)可以大大節(jié)省內(nèi)存空間。下面對(duì) Cholesky分解法、PCCG方法和OpenMP并行解法的法方程求解時(shí)間進(jìn)行比較(不包括法方程的組成時(shí)間)，如表3所示。并行解法所采用的CPU數(shù)為8個(gè)。對(duì)于Cholesky分解法，由于計(jì)算時(shí)間過(guò)于冗長(zhǎng)，僅計(jì)算至120階。

表3顯示，對(duì)于大數(shù)據(jù)量的計(jì)算，3種方法中PCCG方法的計(jì)算時(shí)間最短，恢復(fù)200階的地球重力場(chǎng)僅需要24次迭代，耗時(shí)半個(gè)小時(shí)左右即可完成。比較兩種嚴(yán)密求解方法，并行算法大大提高了求解的效率，這對(duì)于求解超大型的線(xiàn)性估計(jì)問(wèn)題是非常好的一個(gè)選擇。圖2給出了并行解法加速比與CPU個(gè)數(shù)和模型階數(shù)的關(guān)系。其中，左圖的求解規(guī)模為120階地球重力位模型，右圖的CPU數(shù)固定為8個(gè)。圖3則顯示了相應(yīng)的并行算法效率值。

從圖3可以發(fā)現(xiàn)，設(shè)計(jì)的OpenMP并行解法具有超線(xiàn)性加速比。隨著CPU個(gè)數(shù)的增加，并行算法的效率值逐漸降低，但這種現(xiàn)象是建立在降低總解算時(shí)間的基礎(chǔ)之上。此外，隨著求解規(guī)模的增加，并行算法的效率值逐漸增加，充分說(shuō)明了OpenMP并行解法對(duì)于求解超大型線(xiàn)性問(wèn)題的高效性。

表1 不同階數(shù)對(duì)應(yīng)的計(jì)算機(jī)硬件資源需求Tab.1 Computer hardware resource dependence of different degree

表2 IBM P575-2服務(wù)器主要參數(shù)Tab.2 Main parameters of IBM P575-2 server

表3 不同數(shù)值解法的法方程求解時(shí)間比較Tab.3 Comparison among the computation time for normal equation solution with different numerical methods

3.2 PCCG方法與并行解法的精度比較

PCCG方法經(jīng)過(guò)24次迭代計(jì)算即可達(dá)到收斂，其求解效率較之Cholesky分解法和OpenMP并行解法都要高，對(duì)于超大型線(xiàn)性問(wèn)題尤為顯著。梯度模擬數(shù)據(jù)利用EIGEN-GL04C地球重力位模型計(jì)算獲得，階數(shù)計(jì)算至200階，然后基于空域最小二乘法［3-5］，分別利用預(yù)條件共軛梯度方法和OpenMP并行解法恢復(fù)地球重力場(chǎng)，階數(shù)同樣取至200階。結(jié)果見(jiàn)圖4～7。

從圖4可以看出，除低次項(xiàng)以外，PCCG方法能夠以10-16以上的精度恢復(fù)地球重力場(chǎng)，而并行算法則較之更高。相對(duì)而言，低次項(xiàng)的精度較低，這主要是由極空白問(wèn)題引起的。對(duì)比兩種方法的模型階誤差(圖5)可以看出，并行算法精度較之PCCG方法要高一個(gè)數(shù)量級(jí)左右，達(dá)到10-14以上，但PCCG方法收斂較為迅速，且精度能夠滿(mǎn)足求解的需要。對(duì)比兩極大地水準(zhǔn)面以及大地水準(zhǔn)面誤差緯向變化圖(圖6，7)可以發(fā)現(xiàn)，PCCG方法獲取的大地水準(zhǔn)面在極區(qū)引入了一定的誤差，最大在1 mm左右，隨著緯度的升高，其精度愈來(lái)愈低，但是其引入的誤差對(duì)全球重力場(chǎng)的解算可以忽略不計(jì);并行算法則不同，除兩極地區(qū)外，其他區(qū)域的精度分布較為均勻，除了計(jì)算機(jī)舍入誤差外，OpenMP并行解法計(jì)算嚴(yán)密，未引入其他誤差。

圖2 OpenMP并行解法加速比與CPU數(shù)和階數(shù)的關(guān)系Fig.2 Relation between OpenMP parallel algorithm speedup and CPU number as well as number of degree

圖3 OpenMP并行解法效率值與CPU數(shù)和階數(shù)的關(guān)系Fig.3 Relation between OpenMP parallel algorithm efficiency and CPU number as well as number of degree

圖4 PCCG方法與OpenMP并行解法的球諧系數(shù)誤差譜Fig.4 Error spectrum of spherical harmonic coefficients of PCCG method and OpenMP parallel algorithm

圖5 PCCG方法與并行算法的模型階誤差比較Fig.5 Comparison between model’s degree error RMS of PCCG method and OpenMP parallel algorithm

圖6 PCCG方法與OpenMP并行解法的大地水準(zhǔn)面誤差極區(qū)分布Fig.6 Distribution of geoid error RMS in the polar area of PCCG method and OpenMP parallel algorithm

圖7 PCCG方法與OpenMP并行解法的大地水準(zhǔn)面誤差緯向分布Fig.7 Latitudinal dependence of geoid error RMS of PCCG method and OpenMP parallel algorithm

3 結(jié)語(yǔ)

衛(wèi)星重力梯度數(shù)據(jù)確定地球重力位模型最終可化為一大型最小二乘求解問(wèn)題。針對(duì)GOCE衛(wèi)星重力梯度觀(guān)測(cè)數(shù)據(jù)的海量特征，對(duì)地球重力位模型的數(shù)值解法，包括Cholesky分解法、預(yù)條件共軛梯度法以及OpenMP并行解法3種方法，進(jìn)行了系統(tǒng)研究和詳細(xì)分析。研究表明，預(yù)條件共軛梯度方法能夠在滿(mǎn)足求解精度的要求下，有效地提高大型矩陣求逆的效率，但這也帶來(lái)了一定精度的損失;OpenMP并行算法具有簡(jiǎn)單通用、移植性和可擴(kuò)展性好、開(kāi)發(fā)快速的特點(diǎn)，能夠在不損失求解精度的條件下，有效提高計(jì)算效率。特別是在計(jì)算機(jī)硬件條件有限的情況下，OpenMP并行解法無(wú)疑是一個(gè)非常好的選擇。

1 Rummel R，et al.Spherical harmonic analysis of satellite gradiometry［R］.Delft：Netherland Geodetic Commission， 1993.

2 Koop R.Global gravity field modeling using satellite gravity gradiometry［R］.Delft：Netherland Geodetic Commission，1993.

3 徐新禹.衛(wèi)星重力梯度及衛(wèi)星跟蹤衛(wèi)星數(shù)據(jù)確定地球重力場(chǎng)的研究［D］.武漢大學(xué)，2008.(Xu Xinyu.Study of determining the Earth’s gravity field model from satellite gravity gradient and satellite-to-satellite tracking data［D］.Wuhan University，2008)

4 鐘波.基于GOCE衛(wèi)星重力測(cè)量技術(shù)確定地球重力場(chǎng)的研究［D］.武漢大學(xué)，2010.(Zhong Bo.Study on the determination of the Earth’s gravity field from satellite gravimetry mission GOCE［D］.Wuhan University，2010)

5 朱廣彬.衛(wèi)星重力梯度測(cè)量技術(shù)確定地球重力場(chǎng)的理論方法研究［D］.武漢大學(xué)，2010.(Zhu Guangbin.Research on the theory and methodology for the Earth’s gravity field determination using satellite gravity gradiometry measurement［D］.Wuhan University，2010)

6 徐士良.FORTRAN常用算法程序集［M］.北京：清華大學(xué)出版社，1992.(Xu Shiliang.Fortran algorithms commonly used procedures set［M］.Beijing：Tsinghua University Press，1992)

7 Ditmar P and Klees K.A method to compute the Earth’s gravity field from SGG/SST data to be acquired by the GOCE satellite［M］.Delft University Press，2002.

8 鄭鋒，李名世，蔡佳佳.基于OpenMP的并行遺傳算法探討［J］.心智與計(jì)算，2007，1(4)：396-402.(Zheng Feng，Li Mingshi and Cai Jiajia.Parallel genetic algorithms based on OpenMP［J］.Mind and Computation，2007，1 (4)：396-402)

9 李曉梅，吳建平.數(shù)值并行算法與軟件［M］.北京：科學(xué)出版社，2007.(Li Xiaomei and Wu Jianping.Numerical parallel algorithms and software［M］.Beijing：Science Press，2007)

ON NUMERICAL METHODS FOR DETERMINATION OF EARTH GRAVITY FIELD MODEL USING MASS SATELLITE GRAVITY GRADIOMETRY DATA

Zhu Guangbin1)，Chang Xiaotao1)，Zou Xiancai2)，Xu Xinyu2)and Wang Jianqiang3)

(1)Satellite Surveying and Mapping Application Center，SBSM，Bejing 100830
2)School of Geodesy and Geomatics，Wuhan University，Wuhan 430079
3)Institute of Surveying and Mapping，East China Institute of Technology，F(xiàn)uzhou344000)

On the basis of Space-Wise Least Square method，three numerical methods including Cholesky decomposition，Pre-conditioned conjugate gradient and Open Multi-Processing parallel algorithm are applied into the determination of gravity field with satellite gravity gradiometry data.The results show that，Cholesky decomposition method has been unable to meet the requirements of computation efficiency when the computer hardware is limited.Pre-conditioned conjugate gradient method can improve the computation efficiency of huge matrix inversion，but it also brings a certain loss of accuracy.The application of Open Multi-Processing parallel algorithm could achieve a good compromise between accuracy and computation efficiency.

satellite gravity gradiometry;Cholesky decomposition;pre-conditioned conjugate gradient;open multiprocessing parallel algorithm;data processing

1671-5942(2011)06-0140-05

2011-03-19

國(guó)家自然科學(xué)基金(40874012，40904003，40974016，41004007)

朱廣彬，男，1981年生，博士，主要從事衛(wèi)星大地測(cè)量學(xué)的研究.E-mail：whu_gbzhu@hotmail.com

P223

A