海量衛星重力梯度觀測數據確定地球重力位模型的數值方法*

朱廣彬 常曉濤 鄒賢才 徐新禹 王建強

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

海量衛星重力梯度觀測數據確定地球重力位模型的數值方法*

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

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

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

衛星重力梯度;Cholesky分解;預條件共軛梯度;OpenMP并行算法;數據處理

1 引言

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

2 方法介紹

2.1 Cholesky分解法

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

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

2.2 PCCG方法

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

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

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

3)計算新的參數向量xk+1=xk+αkpk;

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

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

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

2.3 OpenMP并行編程模型

OpenMP采用Fork-Join并行執行模型［8］。當程序開始執行時只有主線程存在，主線程會一直串行的執行;當遇到需要并行計算的并行域時，主線程會創建一隊并行的線程開始并行執行;當派生線程在并行域中執行完畢后，它們退出或者掛起，此時只有主線程在運行(圖1)。OpenMP應用編程接口可以根據不同并行域的需要動態地改變線程數，且該并行結構可以嵌入到別的并行結構中去。從本質上講，OpenMP的并行化是通過使用嵌入到源代碼中的編譯制導語句來實現的，其是一個外部編程模型，而不是自動的編程模型。

圖1 OpenMP并行執行模式Fig.1 OpenMP parallel execution mode

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

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

并行算法的效率定義為：

式中，p為并行算法運行所需的處理器數。

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

3 數值分析

3.1 計算效率評測

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

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

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

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

表1 不同階數對應的計算機硬件資源需求Tab.1 Computer hardware resource dependence of different degree

表2 IBM P575-2服務器主要參數Tab.2 Main parameters of IBM P575-2 server

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

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

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

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

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

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

圖4 PCCG方法與OpenMP并行解法的球諧系數誤差譜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并行解法的大地水準面誤差極區分布Fig.6 Distribution of geoid error RMS in the polar area of PCCG method and OpenMP parallel algorithm

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

3 結語

衛星重力梯度數據確定地球重力位模型最終可化為一大型最小二乘求解問題。針對GOCE衛星重力梯度觀測數據的海量特征，對地球重力位模型的數值解法，包括Cholesky分解法、預條件共軛梯度法以及OpenMP并行解法3種方法，進行了系統研究和詳細分析。研究表明，預條件共軛梯度方法能夠在滿足求解精度的要求下，有效地提高大型矩陣求逆的效率，但這也帶來了一定精度的損失;OpenMP并行算法具有簡單通用、移植性和可擴展性好、開發快速的特點，能夠在不損失求解精度的條件下，有效提高計算效率。特別是在計算機硬件條件有限的情況下，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 徐新禹.衛星重力梯度及衛星跟蹤衛星數據確定地球重力場的研究［D］.武漢大學，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衛星重力測量技術確定地球重力場的研究［D］.武漢大學，2010.(Zhong Bo.Study on the determination of the Earth’s gravity field from satellite gravimetry mission GOCE［D］.Wuhan University，2010)

5 朱廣彬.衛星重力梯度測量技術確定地球重力場的理論方法研究［D］.武漢大學，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］.北京：清華大學出版社，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］.心智與計算，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 李曉梅，吳建平.數值并行算法與軟件［M］.北京：科學出版社，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，Fuzhou344000)

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

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

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

P223

A