基于超級計算機平臺的并行解技術在衛星重力測量中的應用*

聶琳娟 申文斌 王正濤 金濤勇

(1)武漢大學測繪學院，武漢 430079 2)湖北水利水電職業技術學院，武漢430070)

基于超級計算機平臺的并行解技術在衛星重力測量中的應用*

聶琳娟1，2)申文斌1)王正濤1)金濤勇1)

(1)武漢大學測繪學院，武漢 430079 2)湖北水利水電職業技術學院，武漢430070)

分析超級計算機平臺的并行解技術應用于衛星重力測量中的相關問題，對涉及的矩陣運算并行化給出了數值計算和分析，并利用衛星重力擾動位觀測基于最小二乘直接解法，比較了OpenMP和MPI兩種并行化技術的計算效率。

衛星重力;并行計算;重力位模型;OpenMP;MPI

1 前言

隨著衛星重力探測計劃的成功實施，地球重力場的研究成為國際大地測量學的熱點問題。如何高效深入地利用當前(CHAMP、GRACE)和未來(GOCE)衛星任務獲取的數據成為真正具有挑戰性的任務。隨著衛星觀測數目的不斷增多，數據處理將成為地球重力場恢復中的一個關鍵問題。對于GRACE衛星的數據，如果30秒的數據采樣，一個月將有86 400個觀測值，如果計算一年平均重力場模型，觀測值將達到1 036 800個，眾所周知，對于觀測方程的系數陣的形成，勒讓德函數的計算是最耗時間的過程，在GS80服務器上，如果球諧展開到90階，一個觀測將花費1.5秒計算時間，如果展開到120階，將花費4秒，因此，對于GRACE衛星任務，預期恢復重力場模型到120階，一個月的數據形成設計矩陣將花費96小時，這還未包括法方程的形成與解算。如果將來的衛星重力梯度計劃GOCE順利實施，期望恢復重力場位系數到250階，對應的未知數個數為62 997個，如果所有數據均以雙精度保存，法方程矩陣的存儲空間要求31.75 G，雖然目前計算機的硬盤存儲不再是關鍵問題，但程序運行需要31.75 G的內存空間，對于目前任何一臺單機都無法滿足要求，因此無論是從計算時間考慮，還是從內存管理考慮，都必須尋找解決問題的新方法，并行計算技術提供了解決這個問題的唯一途徑。

2 衛星重力解算地球重力場的并行化設計分析

利用衛星重力數據解算地球重力場模型引入并行算法的主要目的:在計算規模一定的情況下，盡量縮短計算時間;在給定的硬件設備限制范圍內，擴大問題求解的規模。

算法的評價通常用執行該算法所需的時間和所需的計算機內存容量(空間)來評價，即用算法的時間復雜度和空間復雜度衡量其效率和計算難度。由于我們的問題總是針對一定的規模進行，因此對于時間復雜度給予更多的關注。眾所周知，同一算法的計算機執行時間與計算機型號、程序語言及程序員的技能等有關，因此，這個時間無法精確度量，我們需要基于算法本身來確定算法的時間復雜度，而不受機型以及這個算法如何編碼等因素的影響。時間復雜度簡單地定義為:從輸入數據到計算出結果所需的時間，記為O(n)。表1給出了3類程序時間復雜度。

表1 三類程序時間復雜度計算Tab.1 Time complexity of three types of program

并行算法設計還需要研究并行程序中數據的相關性。能否將順序執行的程序轉換成語義等價的、可并行執行的程序，主要取決于程序的結構形式，特別是其中的數據相關性。如果在機群上程序P1和P2的數據是彼此相關的，則程序不可能設計為并行結構。通常利用伯恩斯坦準則判斷兩個程序段能否并行執行。如果用Ii表示程序Pi操作的讀入數據，Oi表示程序操作的輸出數據，則下列3種情況可以保證程序的并行化設計:

1)I1∩O2=Φ:P1的輸入數據與P2的輸出數據集不相交;

2)I2∩O1=Φ:P2的輸入數據與P1的輸出數據集不相交;

3)O1∩O2=Φ:P1與P2的輸出數據集不相交。

如果保證了算法中存在不相關的進程，從理論上講問題的解決可以通過并行程序實現。并行程序是通過并行語言來表達的，通常產生并行語言有以下3種途徑:

1)設計全新的并行語言;

2)擴展原有串行語言的語法成分，使其具有并行特征;

3)不改變串行語言，而為串行語言提供可調用的并行庫。

最完美的并行語言完全脫離串行語言的束縛，從語言成分上直接支持并行計算，可以使并行程序的書寫更方便更自然，相應的并行程序也更容易在并行機上實現。但是由于并行計算至今還沒有像串行計算那樣統一的馮·諾伊曼模型可供遵循，因此并行機、并行模型、并行算法和并行語言的設計和開發千差萬別，沒有統一標準，導致設計全新的并行語言實現起來難度和工作量都很大。

對于第二種方法，將對串行語言的并行擴充作為原來串行語言的注釋得到新的并行程序，若用原來的串行編譯器編譯，標注的并行擴充部分將不起作用，若使用擴充后的并行編譯器編譯，則該并行編譯器就會根據標注的要求將原來串行執行的部分轉化為并行執行。相對于設計全新的并行語言，此種方法難度有所降低，但仍需重新開發編譯器，使其能夠支持擴充的并行部分，仍具有相當的復雜度。

對于已有的大量串行程序，最簡單的并行化方法是直接加入對已有的并行庫的調用，這是一種對原來的串行程序設計改動最小的并行化方法，這樣原來的串行編譯器也能夠使用，不需要任何修改，新的程序也將具有并行計算能力。當前最流行的MPI并行程序設計就屬于這種方式。對于這3種并行語言的實現方法，目前最常用的是第三種方法。

由以上分析，我們必須針對重力場恢復算法進行并行處理的可行性研究，尋找其中可以并行計算的階段，即此程序段間的數據之間存在互不相關性，同時保證所給出的并行程序結構具有最小的時間復雜度。衛星重力觀測計算重力場模型可采用動力學法、能量法或加速度法，其核心思想即利用不同的衛星觀測數據或其推導值可以得到各種重力觀測量，據此建立起與地球重力場模型位系數之間的函數關系，即觀測方程，從而可用最小二乘方法形成法方程，直接對法方程系數陣求逆可得位系數最佳估值。對于上述計算過程，經分析可以得出衛星重力場模型解算中通常存在的3種類型并行結構:

1)海量觀測數據的并行計算:將讀入的觀測數據等分地分配到各個CPU，各個CPU協同分別計算設計矩陣，最終通過迭加得到設計矩陣A，即利用并行設計將一個任務分成多個子任務，達到減少計算時間的目的;

2)矩陣與矩陣乘和矩陣與向量乘;

3)線性方程組的求解(包括法方程直接求逆或共軛梯度迭代求解)。

對于后兩種并行計算的實現，可以在原有代碼中加入調用已有的并行庫語句，經過并行編譯器編譯，就可以實現并行計算，而不需要對原有程序作大的改動。

3 矩陣運算的MPI并行化設計

由分析可以看出，提高重力場解算效率的關鍵問題之一就是矩陣運算，尤其對于衛星重力數據觀測導致的稠密矩陣，主要包括兩種運算:矩陣相乘、矩陣與向量乘。由于矩陣與向量乘可以轉換為矩陣與一維矩陣乘，因此，本節給出矩陣相乘在GS80并行機上通過調用MPI庫的并行化設計。

給出一般性的假定:A為M×N階矩陣，B為N ×L階矩陣，C為M×L階矩陣，即求解矩陣乘積C =A·B。

由于GS80并行機共有8顆CPU，為了保證最大限度地利用計算資源，令程序進程數NPROCS等于8。在通常情況下，矩陣的維數不為8的整倍數，為描述簡單，假定M和L均為8的整倍數，則矩陣A、B和C的子塊大小分別為MLOC×N，N×LLOC，MLOC×L，A和C按行等分成子塊存儲在8個進程中，B按列等分成子塊存儲在對應的8個進程中，其中MLOC=M/8，LLOC=L/8。圖1給出了矩陣乘法并行化分塊策略。

圖1 矩陣乘法并行化分塊策略Fig.1 Parallel partitioned strategy for matrix multiplication

圖1中分塊矩陣對應的具體存儲方式為:

各個矩陣元素分別對應存儲在進程k的A、B和C數組中。如果保持矩陣A和C的子塊不動，矩陣B的子塊在各個進程間循環移動，即可完成矩陣的并行乘法，程序設計算法流程見圖2。

為了解程序并行化對計算效率的提高程度，在GS80服務器上安裝了MPICH1.2.6版本，它支持所有的MPI-1.2特征和部分的MPI-2的特征，MPICH編譯安裝成功后，所有MPI并行程序的編譯與原先的FORTRAN編譯命令參數一致，只是編譯器由原先的f77/f90變為mpif77，對應程序的運行需要調用mpirun命令，并且需要指定并行的CPU個數，例如mpirun-np 8 program.bat命令指定使用CPU個數為8顆。表2給出了不同大小矩陣相乘在不同并行進程數目下所花費的時間與每秒百萬次浮點運算MFLOPS(Million Floating-point Instruction Per Second)值。

從表2可以看出，對于小運算量，并行算法并不能體現并行的優點，甚至會降低計算速度，這主要是因為數據在總線的傳輸時間大于CPU的計算時間，但對于大計算量的工作，并行算法將以幾何級數提高計算效率(TIME:20s-10s-5s)，同時計算能力也大大增強(MFLOPS:50-100-200)。

圖2 矩陣乘法并行化流程圖Fig.2 Flow chart of matrix multiplication in parallel

表2 GS80服務器的并行計算效率Tab.2 Parallel computing efficiency of GS80 server

4 最小二乘直接解法并行分析

最小二乘直接解法形式上最為簡單，并且能夠提供誤差估計信息，但由于其對于內存的苛刻要求，使得200階以上重力場模型位系數在普通微機上的解算陷入困境。除了內存限制外，構建法方程矩陣及其求逆是一項十分耗時的工作，這些不利因素都可以通過并行計算得到解決。

對于最小二乘直接解法的流程，有3個關鍵的密集型計算任務。第一個是以行形式計算設計矩陣A，第二個是矩陣乘和向量乘分別得到法方程矩陣N和b，最后一個是線性方程系統的求逆。圖3給出了最小二乘直接解法流程圖。

圖3 最小二乘直接法流程圖Fig.3 Fflow chart of least squares direct method

作為實驗算例，給出直接最小二乘算法的耗時，利用16 200個衛星重力擾動位觀測，基于擾動位的球諧展開求解50階重力位系數，程序中涉及的矩陣相乘等的數學運算均直接調用數學函數庫，包括Intel MKL、Goto BLAS、BLACS和ScaLAPACK。通過在程序中設置時間斷點，獲取每一步的計算時間，當僅選取一個CPU時，相當于串行計算，具體耗時為:構建A，耗時0.936 s;構建N耗時12.739 s;構建b耗時0.049 s;解方程耗時0.778 s，總時間耗時14.531 s。從統計結果可以看出，對于直接解最小二乘法，大部分的時間用在計算法方程矩陣N上，構建設計矩陣A和解線性方程系統也需要花費一定的時間。總時間則由觀測個數和未知數個數決定，但各部分的時間消耗比例基本保持不變，不隨問題的大小而改變。下面分別就OpenMP和MPI兩種并行化方法進行詳細的對比與分析。

4.1 OpenMP并行化分析

使用的OpenMP并行環境是Intel Fortran自帶的，通過在編譯命令里增加-openmp參數開啟并執行并行運算。從分析和OpenMP的特性可以看出，設計矩陣A，法方程矩陣N和線性方程系統解算3部分工作可以進行并行計算。

對于多核心的CPU，改變默認線程數n的命令是:

BASH環境:export OMP_NUM_THREADS=n

C-Shell環境:setenv OMP_NUM_THREADS=n

解算問題的規模保持不變，分別設置線程數n =1，2，3，4，計算結果如表3所示。為直觀表示，圖4給出了對比結果。

表3 基于OpenMP環境的運行時間(單位：s)Tab.3 Run-time based on OpenMP parallel environment (unit：s)

圖4 基于OpenMP環境的運行時間Fig.4 Run-time based on OpenMP parallel environment

從圖4可以看出，OpenMP的多線程并行計算有明顯的效果，尤其對于計算法方程矩陣N，計算時間成比例遞減，效果最好;對于計算設計矩陣A和解線性系統也有明顯改善。從圖4還可以看出，由于除N的計算外，其余幾項任務對線程數不敏感，因此當線程數提高到3以上，總的計算時間并不能節約太多。

由此可以看出，基于OpenMP的并行化實現非常容易，只要程序中存在可以被并行化的部分，編譯器將自動分析并給出并行化執行方案，用戶不需要進行交互操作，也不需要對進程執行進行管理，但是這種方法局限于共享內存構架，只適用于SMP和DSM機器。

4.2 MPI并行化分析

相比于OpenMP線程級并行化，MPI具有更好的擴展性，它屬于進程級的并行計算，適合于各種構架，尤其是分布式存儲的特點，除了SMP和DSM機器外，還被廣泛應用于集群系統。

我們使用的MPI并行環境是Open MPI 1.2.6。它整合了其他MPI的優點，其目標是成為最好的MPI運行庫。編譯MPI程序需要執行下列命令:

利用上述命令可以自動鏈接所需要的MPI庫函數，程序的執行必須要求下面的格式:

mpirun--hostfile hosts-np n./program

n表示CPU個數或節點數。

基于MPI并行化計算的最小二乘算法的流程為:

1)初始化進程節點;

2)初始化矩陣和對應的內存分配;

3)分布式計算設計矩陣A;

4)分布式計算矩陣相乘和向量乘積;

5)線性方程系統的并行解;

6)參數估值的匯集;

7)退出各計算節點。

其中矩陣相乘和線性方程系統的解由ScaLAPACK的相關子程序計算。解得的結果分別存儲在分布式向量b中，輸出時，需要將所有節點的系數歸集到一個節點上。

解算問題的規模保持不變，分別設置n=1，2，3，4，計算結果列于表4;圖5給出了對比結果。

表4 基于MPI并行環境的運行時間(單位：s)Tab.4 Run-time based on MPI parallel environment(unit：s)

圖5 基于MPI并行環境的運行時間Fig.5 Run-time based on MPI parallel environment

從圖5可以看出，MPI的多線程并行計算同樣有明顯的效果，可大大加速法方程矩陣N的生成，對于設計矩陣A，由于處理的問題太小，其計算時間沒有成比例的減少，反而隨著CPU數目的增多耗時越來越大。

除了上述兩種并行化設計外，分布式共享內存(OpenMP+MPI)的并行計算機結合了前兩者的特點，是當今新一代并行計算機的一種重要發展方向。

由表5和圖6可以看出，相比于OpenMP，MPI由于涉及分布式內存的信息傳遞，對于這樣的小規模運算，并未能表現出優勢，相反，其運算速度比OpenMP還慢。

表5 不同計算環境運行時間(50階)(單位：s)Tab.5 Run-time based on different computing environments(unit：s)

圖6 不同計算環境運行時間Fig.6 Run-time based on different computing environments

5 結語

對于衛星重力場恢復解算中涉及的3種類型并行結構可以通過不同并行算法實現，其中海量觀測數據的并行計算需要單獨設計并行程序，根據CPU的個數決定并行進程的數目，在不同的CPU節點分別計算設計矩陣，完成各個節點計算后再累加求和即可得到設計矩陣A;對于矩陣相關運算和線性方程組的求解，則可以通過MPI庫函數的調用實現。通過并行解技術以及超級計算機的使用，海量的衛星重力觀測數據處理和更高階地球重力場模型球諧系數的解算將變得更加容易。與此同時，對于不同的并行算法，需要考慮各自的最優適用范圍，利用OpenMP可以簡單的實現算法的并行化，對于小規模的計算，其速度要高于MPI并行程序，主要是MPI涉及消息傳遞，總線速度的快慢對其影響較大;對于循環算法的并行化，最優的選擇是采用純MPI環境，注意避免更多的消息傳遞，尤其對于較慢的總線界面，例如傳統的網絡;當需要更少的消息傳遞時，多線程的MPI程序的效率要遠遠高于單線程程序。

1 Gropp W，Lusk E and Skjellum A.Using MPI:Portable parallel programming with the message-passing interface［M］.Cambridge，MA:MIT Press，1994.

2 Pail R and Plank G.Assessment of three numerical solution strategies for gravity field recovery from GOCE satellite gravity gradiometry implemented on a parallel platform［J］.Journal of Geodesy，2002，76:462-474.

3 Pail R and Plank G.Comparison of numerical solution strategies for gravity field recovery from GOCE SGG observations implemented on a parallel platform［J］.Advances in Geosciences，2003，1:39-45.

4 Xiao H D and Lu Y.Parallel computation for spherical harmonic synthesis and analysis［J］.Computersamp;Geosciences，2007，33(3):311-317.

5 Xie J.Implementation of parallel least-square algorithm for gravity field estimation［R］.Dept.of Geodetic Science，The Ohio State University，Columbus，Ohio，2005.

6 陳國良.并行計算——結構、算法、編程［M］.北京:高等教育出版社，1999.(Chen Guoliang.Parallel computing— structures，algorithms and programming［M］.Beijing: Higher Education Press，1999)

7 李迎春.利用衛星重力梯度測量數據恢復地球重力場的理論與方法［D］.解放軍信息工程大學，2004.(Li Yingchun.Theories and methods of the earth’s gravity field recovery by using of satellite gravity gradiometry data［D］.PLA Information Engineering University，2004)

8 寧津生.衛星重力探測技術與地球重力場研究［J］.大地測量與地球動力學，2002，(1):1-5.(Ning Jinshedng.The satellite gravity surveying technology and research of earth’s gravity field［J］.Journal of Geodesy and Geodynamics，2002，(1):1-5)

9 黃謨濤，等.超高階地球位模型的計算與分析［J］.測繪學報，2001，30(3):208-213.(Huang Motao，et al.Analysis and computation of ultra high degree geopotential model［J］.Acta Geodaetica et Cartographica Sinica，2001，30 (3):208-213)

APPLICATIONS OF PARALLEL SOLUTIONS TECHNOLOGY IN SATELLITE GRAVITY MEASUREMENT

Nie Linjuan1，2)，Shen Wenbin1)，Wang Zhengtao1)and Jin Taoyong1)

(1)School of Geodesy and Geomatics，Wuhan University，Wuhan 430079 2)Hubei Water Resources Technical College，Wuhan 430070)

The parallel algorithm based on supercomputer for solving the Earth’s gravity field model was investigated.Shared and distributed parallel environment(OpenMPamp;MPI)applying to the least squares method are analysed and compared.Different threads，different node to determine the gravity field model at different levels of time consumption and efficiency are compared with each other.

satellite gravity;parallel computing;geopotenial model;OpenMP;MPI

1671-5942(2012)02-0064-06

2011-08-13

國家自然科學基金(41074014);湖北省高等學校青年教師深入企業行動計劃項目(XD2010461);教育部博士點新教師基金(200804861022)

聶琳娟，女，博士生，講師.主要從事大地測量專業理論方法的教學與科研工作.E-mail:ljnie@whu.edu.cn

P223

A