趙鍇坤 朱炬波 谷德峰 韋春博

1 中山大學(xué)物理與天文學(xué)院天琴中心,廣東省珠海市大學(xué)路2號(hào),519082

高精度、高分辨率的重力場(chǎng)模型在氣候變化、資源勘探、災(zāi)害預(yù)防、國(guó)防安全等領(lǐng)域有著非常重要的作用。重力場(chǎng)解算的難度主要來(lái)源于龐大的計(jì)算量與巨大的內(nèi)存占用。法方程的構(gòu)建是重力場(chǎng)反演過(guò)程中的關(guān)鍵步驟,而法矩陣的規(guī)模與重力場(chǎng)階數(shù)相關(guān)[1],是直接影響反演速度與內(nèi)存占用的重要因素。當(dāng)重力場(chǎng)階數(shù)過(guò)大時(shí),法方程的求解速度會(huì)大幅下降,故高階重力場(chǎng)的解算尤為困難[2]。隨著計(jì)算機(jī)技術(shù)的發(fā)展,已有一些新的技術(shù)與高性能算法可用于重力場(chǎng)的快速反演。鄒賢才等[3]使用OpenMP對(duì)重力場(chǎng)解算過(guò)程進(jìn)行并行化處理,計(jì)算速度提升近6倍;陳秋杰等[4]使用MKL(math kernel library)科學(xué)計(jì)算庫(kù)與OpenMP并行計(jì)算庫(kù),顯著提高了重力場(chǎng)解算的效率;周浩等[5-6]基于MPI并行技術(shù)對(duì)法方程構(gòu)建、大規(guī)模矩陣相乘、求逆等關(guān)鍵節(jié)點(diǎn)引入并行計(jì)算,解算120階重力場(chǎng)模型僅耗時(shí)40 min。

除CPU平臺(tái)的并行計(jì)算方法與高性能算法外,近年來(lái)計(jì)算機(jī)圖形學(xué)也得到快速發(fā)展,圖形處理器(GPU)的計(jì)算能力甚至已經(jīng)超過(guò)CPU。GPU編程工具的進(jìn)步也使得各種算法成功移植到GPU上,混合計(jì)算平臺(tái)已經(jīng)實(shí)現(xiàn)了針對(duì)不同學(xué)科應(yīng)用的大幅提速[7]。Hou等[8]提出一種結(jié)合CUDA和OpenMP的混合并行方法,能處理大規(guī)模的數(shù)據(jù);劉世芳等[9]實(shí)現(xiàn)了基于MPI和CUDA的稠密對(duì)稱矩陣三角化算法,取得了較好的加速效果。

上述研究均顯示出并行計(jì)算在重力場(chǎng)相關(guān)研究中能起到顯著作用。為提高使用最小二乘法反演重力場(chǎng)模型的解算速度,本文針對(duì)復(fù)雜的重力場(chǎng)反演流程,聯(lián)合MPI與CUDA設(shè)計(jì)合理的并行方案,同時(shí)實(shí)現(xiàn)內(nèi)存耗用與計(jì)算速度等多方面的優(yōu)化。

1 基于動(dòng)力積分法的重力場(chǎng)模型解算方法

動(dòng)力積分法是較早應(yīng)用于重力場(chǎng)反演的基礎(chǔ)方法之一,具有很高的求解精度,其計(jì)算流程可以概括為:

1)使用先驗(yàn)力學(xué)模型計(jì)算各種攝動(dòng)加速度對(duì)待估參數(shù)的偏導(dǎo)數(shù),并通過(guò)積分器獲得位置對(duì)各待估參數(shù)的偏導(dǎo),構(gòu)建設(shè)計(jì)矩陣;

2)通過(guò)積分器獲得參考軌道序列,并與實(shí)測(cè)數(shù)據(jù)作差,獲得O-C殘差序列b;

3)使用最小二乘法求解方程Ax=b,獲得衛(wèi)星初始狀態(tài)改進(jìn)量ΔX和重力場(chǎng)位系數(shù)的改進(jìn)量Δp;

4)對(duì)初始狀態(tài)和先驗(yàn)重力場(chǎng)模型進(jìn)行改進(jìn),重復(fù)上述流程,直到殘差滿足收斂條件。

下面將詳細(xì)闡述法方程的構(gòu)建與求解過(guò)程,以便對(duì)線性的重力場(chǎng)反演任務(wù)進(jìn)行分配與規(guī)劃。

首先,將設(shè)計(jì)矩陣A按照參數(shù)類型進(jìn)行分塊,表示為:

A=[LG]

(1)

式中,L為低軌衛(wèi)星位置對(duì)局部參數(shù)的偏導(dǎo)數(shù),局部參數(shù)是指如衛(wèi)星初始狀態(tài)向量、加速度計(jì)測(cè)定的非保守力等與時(shí)間有關(guān)的參數(shù);G為低軌衛(wèi)星位置對(duì)全局參數(shù)的偏導(dǎo)數(shù),全局參數(shù)是指如重力場(chǎng)位系數(shù)等與時(shí)間無(wú)關(guān)的參數(shù)。通過(guò)該方式對(duì)A進(jìn)行分塊,有利于多弧段法方程的合并。

將分塊后的單歷元觀測(cè)方程構(gòu)建為法方程:

ATAx=ATy

(2)

(3)

式中,ΔX為局部參數(shù)的改正量,Δp為全局參數(shù)的改正量。上式簡(jiǎn)記為:

(4)

由于單個(gè)弧段的所有歷元的待估參數(shù)完全一致,故單弧段的法方程形式與式(4)一致。但多弧段間的局部參數(shù)存在區(qū)別,需采用式(5)實(shí)現(xiàn)多弧段法方程合并:

(5)

使用最小二乘法將式(5)整理成法方程形式并化簡(jiǎn)為:

(6)

消去局部參數(shù)即可解得重力場(chǎng)位系數(shù)改進(jìn)量:

(7)

將式(7)代入式(6)即可求得衛(wèi)星初始狀態(tài)向量的改進(jìn)量:

(8)

動(dòng)力法反演重力場(chǎng)模型的詳細(xì)流程見(jiàn)圖1。

max(ΔX)表示衛(wèi)星狀態(tài)修正量在各個(gè)方向上的最大值,下同

經(jīng)測(cè)試,解算過(guò)程中耗時(shí)最多的環(huán)節(jié)為數(shù)值積分器計(jì)算參考軌道、大規(guī)模矩陣相乘以及大規(guī)模矩陣求逆。除計(jì)算速度外,由于設(shè)計(jì)矩陣大小近似為3倍歷元個(gè)數(shù)乘以n位系數(shù),而n位系數(shù)又與重力場(chǎng)模型階數(shù)n相關(guān)(式(9)),故當(dāng)重力場(chǎng)模型階數(shù)n過(guò)高時(shí),設(shè)計(jì)矩陣的大小會(huì)劇增,需要兼顧計(jì)算過(guò)程中的內(nèi)存及顯存耗用,故有必要根據(jù)不同需求分別制定不同的并行加速策略:

n位系數(shù)=(n+3)(n-1)

(9)

2 針對(duì)復(fù)雜計(jì)算邏輯的MPI并行策略

MPI提供的各接口可以使不同進(jìn)程間的信息相互傳遞[5],進(jìn)而實(shí)現(xiàn)復(fù)雜任務(wù)的并行化處理。使用MPI實(shí)現(xiàn)按弧段并行構(gòu)建設(shè)計(jì)矩陣及求解法方程,具體任務(wù)分配及運(yùn)行流程見(jiàn)圖2。

圖2 MPI任務(wù)分配及計(jì)算流程

在MPI中,各進(jìn)程編號(hào)通常為從0開(kāi)始的整數(shù)。為便于表述,稱0號(hào)進(jìn)程為主進(jìn)程,其他進(jìn)程為分進(jìn)程。由于軌道數(shù)據(jù)、加速度計(jì)數(shù)據(jù)、姿態(tài)數(shù)據(jù)等均按日期存儲(chǔ),故需要使用主進(jìn)程讀取各軌道文件,按弧段整理完畢后使用MPI_Send命令將數(shù)據(jù)發(fā)送到各分進(jìn)程,再啟動(dòng)積分器執(zhí)行后續(xù)任務(wù)。針對(duì)內(nèi)存有限的單機(jī)計(jì)算,則可以按天數(shù)分配任務(wù),為配合部分對(duì)數(shù)據(jù)連貫性有要求的函數(shù)(如批量讀取加速度計(jì)數(shù)據(jù)的函數(shù)),直接由各分進(jìn)程讀取連續(xù)數(shù)日的數(shù)據(jù)文件,并按弧段進(jìn)行存儲(chǔ)。為使各分進(jìn)程的任務(wù)均勻分配,分進(jìn)程總數(shù)可定為任務(wù)總數(shù)的因數(shù)。

在計(jì)算法方程前,式(1)中低軌衛(wèi)星位置對(duì)重力場(chǎng)位系數(shù)的偏導(dǎo)數(shù)矩陣G的內(nèi)存耗用最大。以單弧段3 h共2 160個(gè)歷元為例,不同重力場(chǎng)模型階數(shù)n對(duì)應(yīng)的矩陣G的規(guī)模見(jiàn)表1。

表1 不同階數(shù)n對(duì)應(yīng)的矩陣G的大小

一般而言,集群的單個(gè)節(jié)點(diǎn)內(nèi)存在16 GB以上,故在解算法方程前可以保證內(nèi)存充足,無(wú)需通過(guò)分塊處理、寫(xiě)文件存儲(chǔ)等方式防止內(nèi)存溢出。

在使用最小二乘法求解時(shí),需要通過(guò)式(2)構(gòu)建法方程,單弧段法方程的構(gòu)建結(jié)果見(jiàn)式(4)。當(dāng)位系數(shù)個(gè)數(shù)充足時(shí),式(4)中的Ngg大小將超過(guò)矩陣G。不同重力場(chǎng)模型階數(shù)n對(duì)應(yīng)的Ngg大小見(jiàn)表2。

表2 不同階數(shù)n對(duì)應(yīng)的矩陣Ngg的大小

顯然,當(dāng)重力場(chǎng)模型階數(shù)達(dá)到一定程度時(shí),Ngg占用的內(nèi)存會(huì)過(guò)大,需要使用合適的方法來(lái)控制其大小。

3 矩陣相乘的并行加速方法

矩陣相乘是重力場(chǎng)解算流程中耗時(shí)最多、內(nèi)存占用最大的環(huán)節(jié)。針對(duì)大規(guī)模矩陣相乘,可利用CUDA將程序編譯為核函數(shù)(kernel)并發(fā)送至設(shè)備(device),進(jìn)而在GPU上實(shí)現(xiàn)矩陣相乘的并行運(yùn)算。CUDA是由NVIDIA推出的運(yùn)算平臺(tái)和編程模型,包含有CUDA指令集架構(gòu)(ISA)與GPU內(nèi)部的并行計(jì)算引擎,使開(kāi)發(fā)者編寫(xiě)的程序能在GPU上高速運(yùn)行。另外,當(dāng)重力場(chǎng)階數(shù)過(guò)大時(shí),反演流程中的各矩陣,尤其是Ngg會(huì)占用大量顯存,而顯存容量一般較小,會(huì)嚴(yán)重限制可解算重力場(chǎng)階數(shù)的上限。本節(jié)將利用CUDA和MPI實(shí)現(xiàn)大規(guī)模矩陣相乘的并行加速與內(nèi)存壓縮,同時(shí)解決耗時(shí)過(guò)長(zhǎng)、內(nèi)存與顯存占用過(guò)大的問(wèn)題。

法方程構(gòu)建過(guò)程中涉及多個(gè)矩陣,其中以Ngg規(guī)模最大、占用內(nèi)存最多,因此本文主要針對(duì)Ngg矩陣進(jìn)行研究。

3.1 使用CUDA實(shí)現(xiàn)并行加速

矩陣G的大小為3倍歷元個(gè)數(shù)乘以n位系數(shù),是設(shè)計(jì)矩陣A的主要成分。由式(9)可以推算出串行計(jì)算Ngg=GTG的乘法次數(shù)約為3×歷元個(gè)數(shù)×n4。而使用CUDA并行計(jì)算時(shí),能使用大量線程(thread)同時(shí)計(jì)算,每個(gè)線程負(fù)責(zé)矩陣中單個(gè)元素的計(jì)算,使乘法次數(shù)下降到3×歷元個(gè)數(shù),大幅減少了計(jì)算次數(shù)。

整個(gè)過(guò)程在設(shè)備上進(jìn)行,計(jì)算完畢后再將結(jié)果轉(zhuǎn)移到主機(jī)。在GPU上計(jì)算的函數(shù)稱為核函數(shù),其邏輯見(jiàn)圖3。設(shè)線程的id為idx,則idx號(hào)線程負(fù)責(zé)的元素位置為(x/列數(shù),y%列數(shù)),%為取模運(yùn)算符;再將GT的第x行與G的第y列對(duì)應(yīng)元素相乘后求和,即可求得Ngg的(x,y)號(hào)元素的值。

圖3 CUDA計(jì)算矩陣相乘的簡(jiǎn)單邏輯

針對(duì)線性代數(shù)運(yùn)算,CUDA提供了專門(mén)的子程序庫(kù)cuBLAS(CUDA basic linear algebra subroutine library),其level-3函數(shù)cublasDgemm適用于double型常規(guī)矩陣相乘運(yùn)算,比未經(jīng)處理的簡(jiǎn)單矩陣并行相乘的核函數(shù)計(jì)算更快,后文的實(shí)驗(yàn)也將使用cuBLAS中的函數(shù)進(jìn)行計(jì)算。

針對(duì)前述內(nèi)存占用過(guò)大的問(wèn)題,可根據(jù)Ngg的特殊性進(jìn)行優(yōu)化,主要有以下2點(diǎn):

1)Ngg具有對(duì)稱性,對(duì)于C=AAT型的矩陣相乘,cuBLAS提供了cublasDsyrk函數(shù)進(jìn)行運(yùn)算,并將結(jié)果矩陣以上三角或下三角的形式存儲(chǔ)。但由于cublasDsyrk本身沒(méi)有將結(jié)果再對(duì)稱地運(yùn)算,故可在CPU上完成對(duì)稱操作。經(jīng)過(guò)該處理后,顯存耗用將減少一個(gè)矩陣G的大小,且計(jì)算量將減半。

2)當(dāng)Ngg占用內(nèi)存接近或超過(guò)顯存上限時(shí),需要對(duì)矩陣G進(jìn)行合適的分塊處理。分塊數(shù)量以及分塊方式的差異都會(huì)對(duì)矩陣相乘的計(jì)算速度產(chǎn)生影響,可根據(jù)實(shí)際情況靈活調(diào)整。

3.2 使用MPI減少內(nèi)存占用

在矩陣相乘環(huán)節(jié),采用矩陣分塊計(jì)算的方法僅能控制顯存,而無(wú)法壓縮矩陣占用的內(nèi)存。為進(jìn)一步提高相同硬件配置下MPI并行進(jìn)程的總數(shù)或能解算的重力場(chǎng)階數(shù)上限,在上節(jié)的基礎(chǔ)上,使用MPI進(jìn)一步對(duì)重力場(chǎng)反演任務(wù)進(jìn)行劃分,通過(guò)去除直接計(jì)算Ngg的步驟實(shí)現(xiàn)分進(jìn)程的內(nèi)存壓縮。主要流程(圖4)如下:

圖4 使用MPI減少分進(jìn)程內(nèi)存峰值的流程

4)代入式(7)計(jì)算最終結(jié)果。

該計(jì)算流程省去了分進(jìn)程中Ngg占用的大部分內(nèi)存,當(dāng)重力場(chǎng)階數(shù)較高時(shí)能將單進(jìn)程內(nèi)存峰值減少一半以上。同時(shí),由于該過(guò)程不涉及文件讀寫(xiě),計(jì)算速度不會(huì)有明顯的下降。

4 算例分析

本節(jié)主要測(cè)試算法的加速效果。使用30 d仿真軌道數(shù)據(jù)模擬重力場(chǎng)月解的計(jì)算場(chǎng)景,解算不同階數(shù)的時(shí)變重力場(chǎng)模型,“真實(shí)”重力場(chǎng)模型選用GGM05S,并不加入任何誤差到仿真軌道數(shù)據(jù)中,先驗(yàn)重力場(chǎng)模型選用EGM96,弧長(zhǎng)統(tǒng)一為3 h,共240個(gè)弧段。

使用的計(jì)算機(jī)配置如下:處理器為Intel(R) Xeon(R) W-2235 CPU @ 3.80 GHz,內(nèi)存為32 GB,六核心,顯卡為quadro P1000,顯存為4 GB,CUDA核心共640個(gè)。

4.1 MPI并行加速效果

本文程序基于Windows 10系統(tǒng)開(kāi)發(fā),使用微軟的Microsoft MPI進(jìn)行并行計(jì)算,為適配Visual Studio 2010,選用的MPI版本為8.1。

以30階重力場(chǎng)解算為例,使用已經(jīng)搭建完畢的重力場(chǎng)解算程序計(jì)算,僅修改并行的進(jìn)程數(shù)量進(jìn)行效率比對(duì),統(tǒng)計(jì)結(jié)果見(jiàn)表3與圖5。

表3 啟動(dòng)不同數(shù)量分進(jìn)程時(shí)的重力場(chǎng)解算耗時(shí)比較

圖5 啟動(dòng)不同數(shù)量分進(jìn)程時(shí)的重力場(chǎng)解算耗時(shí)比較

結(jié)果顯示,隨著并行進(jìn)程數(shù)量的增多,計(jì)算速度逐漸加快,但由于MPI并行本質(zhì)上是邏輯的并行而非硬件層面的并行,計(jì)算速度仍然受到CPU性能的約束,所以進(jìn)程數(shù)量越多,單進(jìn)程的計(jì)算速度越慢。當(dāng)進(jìn)程數(shù)量超過(guò)CPU核心數(shù)目時(shí),單進(jìn)程的效率將大幅降低。在本實(shí)驗(yàn)中,分進(jìn)程數(shù)量從5增加至10時(shí),單進(jìn)程計(jì)算耗時(shí)明顯增加,而整體計(jì)算耗時(shí)并沒(méi)有顯著減少,內(nèi)存峰值卻加倍,造成了不必要的資源占用。故使用MPI進(jìn)行并行計(jì)算時(shí),應(yīng)綜合考慮實(shí)際任務(wù)需求與硬件性能進(jìn)行適當(dāng)?shù)倪x擇,建議總進(jìn)程數(shù)量不超過(guò)CPU核數(shù)。

4.2 CUDA對(duì)矩陣相乘的加速效果

本文使用的CUDA版本為8.0,主要用于大規(guī)模矩陣相乘環(huán)節(jié)。測(cè)試所使用矩陣G的大小均為2 160×N,N表示重力場(chǎng)位系數(shù)個(gè)數(shù)。分別計(jì)算30階、60階和120階的重力場(chǎng)模型對(duì)應(yīng)的Ngg,使用的方法包括串行計(jì)算和CUDA并行計(jì)算,后者又包括簡(jiǎn)單核函數(shù)并行以及cuBLAS庫(kù)函數(shù)并行,計(jì)算結(jié)果見(jiàn)表4。

表4 不同方法計(jì)算矩陣相乘耗時(shí)比較

結(jié)果顯示,使用CUDA并行計(jì)算矩陣相乘的速度遠(yuǎn)快于串行計(jì)算,且矩陣越大,加速效果越明顯。但顯存容量通常要比內(nèi)存更小,在使用CUDA時(shí)應(yīng)綜合考慮顯存限制與計(jì)算性能,對(duì)大矩陣進(jìn)行合適的分塊處理。

4.3 30階和120階重力場(chǎng)解算實(shí)例

使用上述方法搭建動(dòng)力學(xué)法反演重力場(chǎng)模型的并行計(jì)算平臺(tái)。使用時(shí)長(zhǎng)為30 d的仿真軌道數(shù)據(jù)進(jìn)行測(cè)試,完成SST-HL模式下30階重力場(chǎng)模型的解算,共迭代3次。使用GGM05S重力場(chǎng)模型進(jìn)行外符合精度評(píng)估,計(jì)算出的位系數(shù)階誤差RMS見(jiàn)圖6。

圖6 30階重力場(chǎng)解算結(jié)果階誤差RMS

結(jié)果顯示,30階重力場(chǎng)模型的解算精度達(dá)到10-14量級(jí),且從第2次迭代開(kāi)始收斂?？烧J(rèn)為本文搭建的重力場(chǎng)并行解算平臺(tái)有足夠的低階重力場(chǎng)解算能力。本次共啟用10個(gè)MPI分進(jìn)程,單個(gè)弧段耗時(shí)3～6 s,迭代3次耗時(shí)約430 s。串行模式下單弧段耗時(shí)105～110 s,迭代3次耗時(shí)約81 600 s。對(duì)比可知,本文的并行計(jì)算程序在計(jì)算30階重力場(chǎng)模型時(shí)可加速約180倍。

將重力場(chǎng)模型階數(shù)提升到120階,此時(shí)串行模式已經(jīng)無(wú)法完成解算任務(wù)。使用與上一節(jié)仿真實(shí)驗(yàn)相同的條件,計(jì)算出的位系數(shù)階誤差RMS見(jiàn)圖7。

圖7 120階重力場(chǎng)解算結(jié)果階誤差RMS

結(jié)果顯示,經(jīng)過(guò)3次迭代得到的重力場(chǎng)模型的階誤差RMS可達(dá)10-13量級(jí),證明了計(jì)算方法的正確性。該算例啟用5個(gè)分進(jìn)程進(jìn)行計(jì)算,每個(gè)弧段耗時(shí)100～200 s不等,全過(guò)程耗時(shí)約6.3 h。同時(shí),該結(jié)果僅由單機(jī)計(jì)算,對(duì)于高性能的集群同樣適用。

5 結(jié) 語(yǔ)

本文針對(duì)基于最小二乘法的重力場(chǎng)模型解算過(guò)程,聯(lián)合MPI與CUDA實(shí)現(xiàn)了并行加速。主要得出以下幾點(diǎn)結(jié)論:

1)使用MPI和CUDA可分別在全局和局部實(shí)現(xiàn)重力場(chǎng)解算過(guò)程的并行加速。前者負(fù)責(zé)重力場(chǎng)反演任務(wù)拆分,加速比與進(jìn)程數(shù)正相關(guān),但建議進(jìn)程數(shù)量不超過(guò)CPU核數(shù);后者負(fù)責(zé)法矩陣相關(guān)計(jì)算,加速比可達(dá)380以上。

2)使用單機(jī)解算30階重力場(chǎng)時(shí)加速比可達(dá)180;解算120階重力場(chǎng)并迭代3次僅耗時(shí)約6.3 h,且該方法在高性能集群上同樣適用。