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

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

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

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

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

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

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

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

1 引言

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

2 方法介紹

2.1 Cholesky分解法

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

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

2.2 PCCG方法

共軛梯度法的基本思想是利用迭代求得的增量改正上一次迭代的初始向量，所得結(jié)果作為下次迭代的初值。增量方向取逼近其解的最速化方向，且利用上一次的初始向量和增量按此梯度方向更新增量。為了進(jìn)一步提高迭代的收斂速度，可以選擇一個預(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)程序開始執(zhí)行時只有主線程存在，主線程會一直串行的執(zhí)行;當(dāng)遇到需要并行計(jì)算的并行域時，主線程會創(chuàng)建一隊(duì)并行的線程開始并行執(zhí)行;當(dāng)派生線程在并行域中執(zhí)行完畢后，它們退出或者掛起，此時只有主線程在運(yùn)行(圖1)。OpenMP應(yīng)用編程接口可以根據(jù)不同并行域的需要動態(tài)地改變線程數(shù)，且該并行結(jié)構(gòu)可以嵌入到別的并行結(jié)構(gòu)中去。從本質(zhì)上講，OpenMP的并行化是通過使用嵌入到源代碼中的編譯制導(dǎo)語句來實(shí)現(xiàn)的，其是一個外部編程模型，而不是自動的編程模型。

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

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

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

并行算法的效率定義為：

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

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

3 數(shù)值分析

3.1 計(jì)算效率評測

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

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

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

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

表1 不同階數(shù)對應(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ù)值解法的法方程求解時間比較Tab.3 Comparison among the computation time for normal equation solution with different numerical methods

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

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

從圖4可以看出，除低次項(xiàng)以外，PCCG方法能夠以10-16以上的精度恢復(fù)地球重力場，而并行算法則較之更高。相對而言，低次項(xiàng)的精度較低，這主要是由極空白問題引起的。對比兩種方法的模型階誤差(圖5)可以看出，并行算法精度較之PCCG方法要高一個數(shù)量級左右，達(dá)到10-14以上，但PCCG方法收斂較為迅速，且精度能夠滿足求解的需要。對比兩極大地水準(zhǔn)面以及大地水準(zhǔn)面誤差緯向變化圖(圖6，7)可以發(fā)現(xiàn)，PCCG方法獲取的大地水準(zhǔn)面在極區(qū)引入了一定的誤差，最大在1 mm左右，隨著緯度的升高，其精度愈來愈低，但是其引入的誤差對全球重力場的解算可以忽略不計(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é)語

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

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ù)確定地球重力場的研究［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)星重力測量技術(shù)確定地球重力場的研究［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)星重力梯度測量技術(shù)確定地球重力場的理論方法研究［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

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

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

P223

A