OpenMP并行算法在衛(wèi)星重力場(chǎng)模型反演中的應(yīng)用*

周 浩 鐘 波 羅志才 張 坤

(1)武漢大學(xué)測(cè)繪學(xué)院，武漢 430079

2)地球空間環(huán)境與大地測(cè)量教育部重點(diǎn)實(shí)驗(yàn)室，武漢430079)

OpenMP并行算法在衛(wèi)星重力場(chǎng)模型反演中的應(yīng)用*

周 浩1)鐘 波1，2)羅志才1，2)張 坤1)

(1)武漢大學(xué)測(cè)繪學(xué)院，武漢 430079

2)地球空間環(huán)境與大地測(cè)量教育部重點(diǎn)實(shí)驗(yàn)室，武漢430079)

利用衛(wèi)星重力數(shù)據(jù)反演地球重力場(chǎng)需要解決重力場(chǎng)模型的高效計(jì)算問題。分析了最小二乘直接法求解重力場(chǎng)模型涉及的密集型計(jì)算任務(wù)，基于OpenMP實(shí)現(xiàn)了衛(wèi)星重力場(chǎng)模型直接求解的并行算法。利用30天、5秒采樣間隔的沿軌擾動(dòng)位T和徑向擾動(dòng)重力梯度Trr數(shù)據(jù)，分別反演了60階次的衛(wèi)星重力場(chǎng)模型，計(jì)算結(jié)果表明，OpenMP并行算法能夠有效提高直接法求解衛(wèi)星重力場(chǎng)模型的計(jì)算效率，并具有很好的穩(wěn)定性。

并行算法;重力場(chǎng)模型;衛(wèi)星重力;重力梯度;OpenMP

1 引言

衛(wèi)星重力測(cè)量技術(shù)為解決全球高覆蓋率、高空間分辨率和高時(shí)間重復(fù)率的重力測(cè)量開辟了新途徑，特別是新一代衛(wèi)星重力計(jì)劃的成功實(shí)施，已使得地球重力場(chǎng)模型的中長(zhǎng)波部分精度提高了約2個(gè)量級(jí)［1-3］。

新一代衛(wèi)星重力測(cè)量數(shù)據(jù)具有高采樣率、全球觀測(cè)的特點(diǎn)。利用衛(wèi)星重力觀測(cè)值反演地球重力場(chǎng)將面臨海量數(shù)據(jù)的處理問題，特別是對(duì)于高階重力場(chǎng)模型的反演，其計(jì)算量巨大，計(jì)算過程極為耗時(shí)。為了利用海量衛(wèi)星重力觀測(cè)值快速求解高精度的地球重力場(chǎng)模型，許多學(xué)者一方面是對(duì)現(xiàn)有重力場(chǎng)恢復(fù)方法進(jìn)行改進(jìn)，提出了一些快速解算法［4-6］;另一方面就是采用并行算法，即基于并行計(jì)算機(jī)，編寫合適的并行算法來提高計(jì)算速度［7-9］。本文對(duì)最小二乘直接求解過程中的密集型計(jì)算任務(wù)進(jìn)行了分析，并基于OpenMP實(shí)現(xiàn)其并行求解，最后通過數(shù)值計(jì)算對(duì)其有效性進(jìn)行了分析與驗(yàn)證。

2 重力場(chǎng)求解的并行算法［9－14］

重力場(chǎng)最小二乘直接求解的并行算法的計(jì)算步驟如圖1所示。從圖1可以看出，構(gòu)建設(shè)計(jì)矩陣A、構(gòu)建法方程系數(shù)陣N和b，以及解算法方程3個(gè)部分為計(jì)算的“熱點(diǎn)”，應(yīng)該引入并行算法進(jìn)行優(yōu)化。下面將基于OpenMP實(shí)現(xiàn)各個(gè)計(jì)算“熱點(diǎn)”的并行。

圖1 直接解法計(jì)算流程圖Fig.1 Flow chart of the computation with direct method

2.1 構(gòu)建設(shè)計(jì)矩陣A的并行

由于設(shè)計(jì)矩陣A的行與行之間無相關(guān)性，且每行的構(gòu)建具有相似性，為此，程序中使用BUILD_A_ THREAD子函數(shù)來構(gòu)建矩陣A中的一行，然后根據(jù)OpenMP語法需求，加入組合并行共享任務(wù)制導(dǎo)語句!$OMP PARALLEL DO/!$OMP END PARALLEL DO實(shí)現(xiàn)構(gòu)建設(shè)計(jì)矩陣A的并行。具體形式為：

!$OMP PARALLEL DO

DO I=1，NOBS

CALL BUILD_A_THREAD(……)

END DO

!$OMP END PARALLEL DO

其中，NOBS為觀測(cè)值個(gè)數(shù)，即設(shè)計(jì)矩陣A的行數(shù)。

2.2 矩陣相乘的并行

由圖1可知，根據(jù)最小二乘構(gòu)建法方程包含兩個(gè)矩陣相乘運(yùn)算，均是計(jì)算的“熱點(diǎn)”，有必要實(shí)現(xiàn)矩陣相乘的并行，以加速構(gòu)建法方程。構(gòu)建法方程可用直接法、分塊法、斯特拉森法等矩陣相乘的并行算法。

1)直接法

設(shè)矩陣A=(aij)為m×n矩陣，矩陣B=(bij)為n×l矩陣，其乘積C=AB=(cij)為m×l矩陣，其中：

式(1)各次計(jì)算中都要用到i、j、k，所以在并行過程中，需要各個(gè)節(jié)點(diǎn)都有一份這些參量的拷貝。根據(jù)OpenMP語法需求，需要將i、j、k設(shè)置為PRIVATE屬性，同時(shí)使用并行共享任務(wù)制導(dǎo)語句! $OMP PARALLEL DO/!$OMP END PARALLEL DO實(shí)現(xiàn)并行。具體形式為：

!$OMP PARALLEL

!$OMP DO PRIVATE(i，j，k)

DO i=1，M

DO j=1，L

C(i，j)=0.0

DO k=1，N

C(i，j)=C(i，j)+A(i，k)*B(k，j)

END DO

END DO

END DO

!$OMP END DO

!$OMP END PARALLEL

2)分塊法

行列分塊法的具體計(jì)算過程是首先將矩陣A、B、C分解為：

然后由A、B陣中的分塊相乘得到C陣中的對(duì)應(yīng)分塊，例如C11=A11B11+A12B21。至于分塊矩陣乘法的并行，其處理方式與直接法相同。

3)斯特拉森(Strassen)算法

利用此種方法主要是為了減少矩陣乘法中的運(yùn)算量，降低時(shí)間復(fù)雜度。與分塊法一樣，首先按照式(2)將矩陣A、B、C進(jìn)行行列分解，然后引入中間變量：

利用中間變量相加得到最后結(jié)果：

斯特拉森算法的實(shí)質(zhì)是反復(fù)使用二階矩陣乘法，使計(jì)算復(fù)雜度從O(n3)降為O(n2.81)。為了保證在M陣計(jì)算完成以后才計(jì)算C陣，可以引入隱含并行同步語句!$OMP BARRIER，具體形式為：

!$OMP PARALLEL

!$OMP DO PRIVATE(i，j，k)

實(shí)現(xiàn)公式(3)

!$OMP END DO

!$OMP BARRIER

!$OMP DO PRIVATE(i，j，k)

實(shí)現(xiàn)公式(4)

!$OMP END DO

!$OMP END PARALLEL

2.3 求解法方程的并行算法

根據(jù)計(jì)算流程，得到法方程中各個(gè)矩陣后，則可對(duì)法方程進(jìn)行求解。并用直接法、平方根法和改進(jìn)的平方根法實(shí)現(xiàn)相應(yīng)的OpenMP并行算法。

1)直接法

直接法分為矩陣求逆和矩陣相乘。其中，采用高斯-約當(dāng)方法對(duì)矩陣求逆。由于高斯-約當(dāng)法求逆是按照行序i進(jìn)行的，故按照下面方式實(shí)現(xiàn)并行：

DO i=1，U

處理對(duì)角線元素

!$OMP PARALLEL SHARED(N，i)

!$OMP DO PRIVATE(j)

!$OMP FLUSH(N)

處理主行元素

!$OMP END DO

!$OMP BARRIER

!$OMP DO PRIVATE(j，k)

!$OMP FLUSH(N)

處理非主行非主列元素

!$OMP END DO

!$OMP BARRIER

!$OMP DO PRIVATE(k)

!$OMP FLUSH(N)

處理主列元素

!$OMP END DO

!$OMP BARRIER

!$OMP END PARALLEL

END DO

其中，U為法方程系數(shù)陣N的維數(shù)。加入!$OMP FLUSH和!$OMP BARRIER語句，是為了保證數(shù)據(jù)的實(shí)時(shí)刷新與嚴(yán)格同步［13］。

2)平方根法［15］

由于法方程系數(shù)陣N=ATA為對(duì)稱正定陣，故可以利用平方根法求解法方程。法方程系數(shù)陣N可以作喬列斯基(Cholesky)分解，如果設(shè)定一次乘法和加法運(yùn)算時(shí)間或一次除法運(yùn)算時(shí)間為一個(gè)單位時(shí)間，則采用平方根方法的計(jì)算時(shí)間為(n3+3n2+ 2n)/6，計(jì)算整個(gè)方程的時(shí)間為(n3+9n2+2n)/6。

3)改進(jìn)的平方根法

!$OMP PARALLEL DO SHARED(N，j)num_ threads(U-j-1)

DO j=1，U

!$OMP PARALLEL DO SHARED(N，p，j) num_threads(U-j-1)

DO p=j+1，U

DO q=j+1，p

N(p，q)=N(p，q)-N(p，j)*N(q，j)/N(j，j)

END DO

END DO

!$OMP PARALLEL DO SHARED(N，j)

DO i=j+1，U

N(i，j)=N(i，j)/N(j，j)

END DO

END DO

利用上述代碼實(shí)現(xiàn)N陣的分解后，利用簡(jiǎn)單的矩陣乘法即可獲得法方程的解。

3 算例與分析

計(jì)算采用Dell R900共享內(nèi)存并行計(jì)算機(jī)，處理器為4顆4核的Intel E7430，主頻2.13GHz，內(nèi)存64GB，編譯環(huán)境為Intel Visual Fortran 10.1?；贕OCE衛(wèi)星模擬軌道數(shù)據(jù)，采用60階次的EIGENGL04C模型分別模擬了沿軌30天、5秒采樣間隔(總共518 400個(gè)觀測(cè)數(shù)據(jù))的擾動(dòng)位T和徑向擾動(dòng)重力梯度Tzz數(shù)據(jù)，并進(jìn)行重力場(chǎng)反演和精度比較分析。

3.1 設(shè)計(jì)矩陣計(jì)算與分析

首先設(shè)定不同的線程數(shù)來構(gòu)建設(shè)計(jì)矩陣A(維數(shù)為518 400×3 721)。以Tzz數(shù)據(jù)計(jì)算為例，獲得計(jì)算所需時(shí)間如表1所示。

表1 構(gòu)建設(shè)計(jì)矩陣A的時(shí)間分析(單位：s)Tab.1 Runtime of building design matrix A(unit：s)

表1中，“加速比”是指在相同任務(wù)量下的并行執(zhí)行時(shí)間與串行執(zhí)行時(shí)間之比。若加速比與線程數(shù)的比例能夠接近于1，則表明該算法能夠達(dá)到線性加速比，并行效率高［9，14］。由表1中第1～3列可知，對(duì)構(gòu)建設(shè)計(jì)矩陣A的并行基本上可達(dá)到線性加速比，由此可知構(gòu)建設(shè)計(jì)矩陣A的并行計(jì)算效率較高;而第4列和第5列之所以未能夠達(dá)到線性加速比，是因?yàn)楸疚牟捎玫臄?shù)據(jù)量較小，導(dǎo)致有些線程還沒有參與計(jì)算就已結(jié)束了計(jì)算任務(wù)。

3.2 矩陣乘法的計(jì)算與分析

設(shè)定不同的線程數(shù)，采用直接法、分塊法和Strassen法分別實(shí)現(xiàn)矩陣乘法的并行計(jì)算，數(shù)據(jù)計(jì)算結(jié)果如圖2所示。

圖2 不同矩陣相乘的并行計(jì)算效率比較Fig.2 Comparison among the efficiency of parallel methods for matrix multiplication

在評(píng)價(jià)并行算法時(shí)，我們通常假定任務(wù)量W一定，若隨著線程數(shù)p的增加，計(jì)算所需時(shí)間t縮短，且近似滿足t=W/p，則該并行算法是成功的。由圖2可知，3種方法的并行效率都趨近于反比例函數(shù)，這說明3種矩陣乘法的并行算法均是成功的。

從圖2還可知：在相同線程數(shù)、相同任務(wù)量的前提下，分塊法計(jì)算時(shí)間略小于直接法，而Strassen法的計(jì)算時(shí)間則最小。

3.3 法方程解算與分析

采用直接法、平方根法和改進(jìn)的平方根法，使用不同階數(shù)的法方程陣，分別實(shí)現(xiàn)了法方程求解的并行算法(計(jì)算中設(shè)定12個(gè)線程數(shù))，計(jì)算結(jié)果如圖3所示。

圖3 不同求解法方程方法的并行效率比較Fig.3 Comparison among the parallel efficiency of different methods for solving the normal equation

從圖3可以看出：在相同線程數(shù)的前提下，隨著任務(wù)量的增加，計(jì)算所需時(shí)間迅速增加;在解算所需時(shí)間上，改進(jìn)的平方根法相比于直接法和平方根法具有極大的優(yōu)勢(shì)，這說明改進(jìn)的平方根法的時(shí)間復(fù)雜度最小。但為了獲得更高精度的解算結(jié)果和精度評(píng)定的參數(shù)協(xié)因數(shù)陣，本文選用直接法對(duì)法方程進(jìn)行求解。

3.4 重力場(chǎng)反演計(jì)算與分析

分別采用518 400個(gè)觀測(cè)歷元的擾動(dòng)位T和二階徑向擾動(dòng)重力梯度Tzz數(shù)據(jù)來反演60階次的地球重力場(chǎng)模型，其中3個(gè)計(jì)算“熱點(diǎn)”分別使用構(gòu)建設(shè)計(jì)矩陣的并行算法、Strassen矩陣相乘的并行算法、直接法求解法方程的并行算法，計(jì)算時(shí)間如表2所示。

表2 反演衛(wèi)星重力場(chǎng)模型計(jì)算時(shí)間表(單位：s)Tab.2 Computational runtime of gravity model inversion (unit：s)

從表2可知，在本文計(jì)算條件下，使用串行算法利用一個(gè)月的觀測(cè)數(shù)據(jù)反演60階的重力場(chǎng)模型，約需要12小時(shí)，當(dāng)采用12個(gè)線程的OpenMP計(jì)算時(shí)，計(jì)算時(shí)間可節(jié)約8個(gè)小時(shí)。若使用更長(zhǎng)時(shí)間跨度的觀測(cè)數(shù)據(jù)，理論上節(jié)省的計(jì)算時(shí)間將更長(zhǎng)。

采用模型階誤差RMS分別對(duì)擾動(dòng)位T和二階徑向擾動(dòng)重力梯度Tzz數(shù)據(jù)反演的重力場(chǎng)模型精度進(jìn)行評(píng)價(jià)：

圖4 利用T和Tzz求解模型的階誤差RMSFig.4 Degree error RMS of different models for solving the T and Tzz

4 結(jié)束語

基于OpenMP實(shí)現(xiàn)了由最小二乘直接法求解衛(wèi)星重力場(chǎng)模型的并行算法。通過數(shù)值計(jì)算分析表明：采用按行存儲(chǔ)方式實(shí)現(xiàn)的構(gòu)建設(shè)計(jì)矩陣的并行算法，基本可達(dá)到線性加速比;采用Strassen法實(shí)現(xiàn)的矩陣乘法并行計(jì)算能夠減小計(jì)算的復(fù)雜度，極大提高計(jì)算效率;采用直接解法求解法方程的并行算法能提高計(jì)算效率，并可保證求解結(jié)果的準(zhǔn)確性和可靠性。因此，基于OpenMP的并行算法是反演地球重力場(chǎng)模型的有效工具。由于OpenMP是基于共享存儲(chǔ)系統(tǒng)上的并行編程標(biāo)準(zhǔn)，系統(tǒng)開銷較大，將在后續(xù)工作中利用OpenMP與MPI的聯(lián)合編程，實(shí)現(xiàn)更長(zhǎng)時(shí)間序列數(shù)據(jù)反演超高階次地球重力場(chǎng)的并行運(yùn)算。

1 寧津生.衛(wèi)星重力探測(cè)技術(shù)與地球重力場(chǎng)研究［J］.大地測(cè)量與地球動(dòng)力學(xué)，2002，(1)：1-5.(Ning Jinsheng.The satellite gravity surveying technology and research of Earth’s gravity field［J］.Journal of Geodesy and Geodynamics，2002，(1)：1-5)

2 許厚澤.衛(wèi)星重力研究：21世紀(jì)大地測(cè)量研究的新熱點(diǎn)［J］.測(cè)繪科學(xué)，2001，26(3)：1-3.(Xu Houze.Satellite gravity mission-new hot point in Geodesy［J］.Science of Surveying and Mapping，2001，26(3)：1-3)

3 孫文科.低軌道人造衛(wèi)星(CHAMP、GRACE、GOCE)與高精度地球重力場(chǎng)——衛(wèi)星重力大地測(cè)量的最新發(fā)展及其對(duì)地球科學(xué)的重大影響［J］.大地測(cè)量與地球動(dòng)力學(xué)，2002，(1)：92-100.(Sun Wenke.Satellite in low orbit (CHAMP，GRACE，GOCE)and high precision Earth gravity field：the latest progress of satellite gravity geodesy and its great influence on geosciences［J］.Journal of Geodesy and Geodynamics，2002，(1)：92-100)

4 Sanso F and Tscherning C C.Fast spherical collocation：theory and examples［J］.Journal of Geodesy，2003，77：101 -112.

5 Klees R，et al.Efficient gravity field recovery from GOCE gravity gradient observations［J］.Journal of Geodesy，2000，74：561-571.

6 Sneeuw N.A semi-analytical approach to gravity field analysisfrom satellite observations［D］.University of München，2000.

7 Klees R，et al.GOCE gravity field recovery using massive parallel computing［A］.Gravity，Geoid and Geodynamics 2000［C］.F Sansò(eds.)：IAG symposium，2002，123：109-116.

8 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.

9 Wittwer T.An introduction to parallel programming［M/ OL］.http：//www.vssd.nl/hlf/a019.htm.

10 羅志才.利用衛(wèi)星重力梯度數(shù)據(jù)確定地球重力場(chǎng)的理論和方法［D］.武漢測(cè)繪科技大學(xué)，1996.(Luo Zhicai.The theory and methodology for determining the earth’s gravity field using satellite gravity gradiometry data［D］.Wuhan Technical University of Surveing and Mapping，1996)

11 王正濤.衛(wèi)星跟蹤衛(wèi)星測(cè)量確定地球重力場(chǎng)的理論與方法［D］.武漢大學(xué)，2005.(Wang Zhengtao.Theory and methodology of earth gravity field recovery by satellite-tosatellite tracking data［D］.Wuhan University，2005)

12 鐘波.基于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)

13 Hermanns M.Parallel programming in Fortran 95 using OpenMP［M/OL］.http：//www.openmp.org/presentations/miguel/F95_OpenMPv1_v2.pdf.

14 陳國(guó)良，等.并行算法實(shí)踐［M］.北京：高等教育出版社，2004.(Chen Guoliang，et al.Parallel algorithm practice［M］.Beijing：Higher Education Press，2004)

15 朱方生，李大美，李素貞.計(jì)算方法［M］.武漢：武漢大學(xué)出版社，2003.(Zhu Fangsheng，Li Damei and Li Shuzhen.Numerical method［M］.Wuhan：Wuhan University Press，2003)

APPLICATION OF PARALLEL ALGORITHMS BASED ON OpenMP TO SATELLITE GRAVITY FIELD RECOVERY

Zhou Hao1)，Zhong Bo1，2)，Luo Zhicai1，2)and Zhang Kun1)

(1)School of Geodesy and Geomatics，Wuhan University，Wuhan 430079 2)Key Lab.of Geospace Environment and Geodesy，Ministry of Education，Wuhan University，Wuhan 430079)

The massive satellite data presents a significant computational challenge to estimate the gravity field.Analyzing the computing-intensive tasks in recovering the earth gravity field model based on the direct leastsquares algorithms，on applying parallel computing techniques capable of rigorously inverting geopotential coefficients using direct method on the basis of OpenMP is focused.Eventually，the earth’s gravity fields complete to degree and order 60 are recovered respectively on the basis of monthly along track disturbing potential observations and gravity gradient observations(with 5 s sampling interval).The simulation shows that introducing the parallel algorithms based on OpenMP into determination of the earth’s gravity field model can be more effective，and is of great reliability.

parallel algorithms;gravity filed model;satellite gravimetry;gravity gradient;OpenMP

1671-5942(2011)05-0123-05

2011-06-17

國(guó)家自然科學(xué)基金(40874002，41104014);國(guó)家863計(jì)劃項(xiàng)目(2008AA12Z105);測(cè)繪遙感信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室專項(xiàng)科研經(jīng)費(fèi)資助項(xiàng)目(2009);地球空間壞境與大地測(cè)量教育部重點(diǎn)實(shí)驗(yàn)室開放基金(10-02-09)

周浩，男，1987年生，碩士研究生，主要從事衛(wèi)星重力學(xué)研究.E-mail：lengyeshanren@126.com

P223，P228

A