基于能量守恒方法的重力場反演快速異構并行算法

譚勖立 王慶賓 范 雕 馮進凱 黃 炎 黃子炎

1 信息工程大學地理空間信息學院，鄭州市科學大道62號，450001

隨著衛星重力測量技術的發展，利用重力衛星觀測數據反演地球重力場已成為普遍的重力場研究手段[1]。針對重力衛星觀測特點，眾多學者提出不同的重力場反演方法，總體上分為空域法和時域法[2]，時域法又可分為能量守恒方法[3-5]、短弧法[1,6-8]、動力學法[9]等。時域法反演地球重力場時，在處理海量觀測數據、獲取觀測方程以及構建并解算法方程的過程中存在計算量巨大、計算耗時長的問題[10]。針對以上問題，許多學者基于并行計算平臺，分析反演方法的并行潛力，設計了并行算法，以提高計算效率。文獻[10]總結并行技術提高反演效率的關鍵任務，提出了OpenMP并行算法；文獻[11]實現了OpenMP并行快速解算重力場模型，反演得到60階模型，解算耗時減少75%；文獻[12-13]在曙光集群高性能計算系統上實現了MPI并行算法，解算120、240階模型時，法矩陣求逆耗時僅為229 s、7 359 s，解算總耗時分別為40 min、7 h；文獻[14]綜合OpenMP和MKL(Math Kernel Library)庫二者優勢，有效提高了反演效率；文獻[15]分析基于超級計算機平臺的并行技術應用于衛星重力測量的相關問題。

現今，反演重力場的并行算法研究主要針對OpenMP、MPI等CPU同構并行展開，存在并行算法硬件需求高、未充分挖掘反演方法并行潛力的問題。隨著GPU(graphic processing unit)的發展，GPU單個核心運算能力和芯片核心集成數量都得到大幅提升，CPU+GPU組成的異構并行結構成為常見的異構并行運算平臺，單臺個人電腦即可構成其主體部分。相較于CPU，GPU具有多核心、高集成的特點，能夠同時開辟上千條線程進行運算，從而實現密集運算的高度并行化。

基于以上分析，本文針對衛星跟蹤衛星觀測模式，在能量守恒方法的基礎上，提出一種重力場反演快速異構并行算法(后文簡稱快速反演算法)，利用CUDA在GPU端實現了并行計算觀測方程設計矩陣與自由項向量，結合MKL庫與分區平差法[16]、預處理共軛梯度法[17]在CPU端完成了低內存消耗下的法方程的快速構建和解算。用該算法處理GRACE-FO衛星2020-01-01～06-30期間共13.77 GB觀測數據，反演獲得120階重力場模型GM-GraceFO2020h(GM為gravity model縮寫，后綴GraceFO表示解算采用GRACE-FO數據，2020代表所用數據年份，h為heterogeneous parallelism的首字母)，并與現有GRACE模型和算法相對比。結果表明，快速反演算法所得模型與現有模型精度相當，且相較于傳統串行算法，反演耗時減少98.479%，內存消耗減小1個量級。算法在保證反演精度的同時，能有效提高計算效率，縮短反演時間，降低對計算平臺的硬件需求，可以為大規模處理衛星重力數據提供經濟、高效的實施方案。

1 能量守恒方法原理

能量守恒方法反演地球重力場，在慣性系下可得單星能量守恒積分式[5]：

(1)

Vconv=Vsolar-tide+Vlunar-tide+Vsolid-tide+

Vocean-tide+Vatmosphere-tide+Vpole-tide

(2)

雙星積分式由單星積分式作差獲得：

(3)

式中，下標1、2為單星的編號，下標12表示1、2星對應項相減。將擾動位T球諧展開：

(4)

綜合式(1)～(4)，顧及球坐標與直角坐標的轉換關系，基于單星和雙星能量守恒方程可構建統一的觀測方程：

(5)

2 算法設計

并行算法的設計主要圍繞兩個方面，分別為觀測方程的并行計算和法方程的并行計算，下面圍繞這兩點進行討論。

2.1 觀測方程的并行計算

觀測方程的并行計算由CUDA實現。CUDA將計算任務劃分為如圖1所示格網，每個GPU線程完成對應格網的任務，線程通過索引值index獲取對應歷元的觀測數據，索引值index的計算方式如下:

圖1 任務格網Fig.1 Task-grid

index=blockIdx.x*blockDim.x+threadIdx.x

(6)

式中，blockDim為線程塊格網的維度，blockIdx為線程塊在格網中的編號，threadIdx為線程在線程塊中的編號。

GPU可開辟與其核心數同樣多的線程。以GTX1080為例，其擁有2 560個核心，可開辟2 560個線程同時進行運算，故GPU可充分利用能量守恒方法設計矩陣各行向量的獨立性，實現觀測方程構建過程的高度并行化，在單位時間內完成更多的計算。同時，為減少計算過程中的循環次數，參考文獻[19]，將中間變量向量化：

(7)

進而設計矩陣行向量Ai及自由項元素li，皆可由向量運算獲得(其中符號“·”表示向量內積，符號“°”表示同維向量對應元素相乘)。

(8)

2.2 法方程的并行計算

MKL(math kernel library)庫是由Intel提供的函數庫，經過了高度優化，并實現了CPU上的并行計算，能提供包括線性代數在內的各類數學運算，具有較高的性能和穩定性。利用MKL庫中的cblas_dgemm函數可實現大型矩陣的乘法快速并行運算。

在構建法方程時，傳統算法需要將設計矩陣整體記錄在內存，當面對海量觀測數據時，設計矩陣將占用巨大的內存空間。以衛星182 d數據、60 s采樣間隔反演120階模型為例，此時將獲得262 080條數據記錄，設計矩陣需28.12 G內存空間進行存儲，個人電腦難以提供足夠大小的內存。針對大型設計矩陣在存儲時存在內存需求大的問題，可采用分區平差的方法，對觀測數據按時序進行分區處理。將觀測數據分區后，可得各分區的誤差方程：

(9)

式中，Ai、Bi和Li分別為i分區的設計矩陣和自由項向量。各分區分別構建法方程，消去區域性未知數Yi，得到僅含聯系未知數的約化法方程：

(10)

將各聯系未知數的約化法方程系數和自由項相加，獲得總體約化法方程的法矩陣與自由項：

(11)

解算總體約化法方程便能求得聯系未知數X，若需進一步求取區域性未知數，則將X回代至各分區法方程中即可。利用分區平差法，每次計算僅需原有內存需求的ni/n(n為總觀測數，ni為各分區觀測數)。與MKL庫結合后，在減小對硬件資源消耗的同時能夠快速完成法矩陣及自由項計算。

在解算法方程時，采用預處理共軛梯度法[17]，并引入MKL庫并行加速相關線性代數運算，可以兼顧解算精度和計算效率。

3 實驗結果與分析

3.1 實驗平臺與數據準備

實驗計算平臺為戴爾Precision 7530移動工作站，CPU為Intel Xeon(頻率為2.9 GHz，6核12線程)，內存容量32 GB，GPU為NVIDIA Quadro P3200(6 G顯存，1 792個CUDA核心)。程序編譯環境為Visual Studio 2015-Visual C++，CUDA10.0，OpenMP2.0并行庫與MKL19.0科學計算庫。實驗所涉及算法都在Debug模式下進行測試。

本文數據來源如下：1)由德國GFZ(German Research Centre for Geosciences) ISDC(Information System and Data Center for Geoscientific Data)數據中心提供的GRACE-FO衛星1 Hz精密軌道數據GNV1B、加速度計數據ACT1B、衛星姿態數據SCA1B以及0.5 Hz星間激光干涉測距數據LRI1B、0.2 Hz Ka/K波段測距數據KBR1B(測距數據主要用于雙星動能差的計算[6,17])，各類數據都以每天一個文件的形式保存；2)海洋潮汐數據來自于EOT11a模型；3)日月星歷數據由DE430模型計算獲得。

對于衛星數據，首先利用觀測數據標準差依據3σ準則進行粗差剔除；然后用5階牛頓插值填補跨度小于30 s的缺失歷元數據，對于無法內插的歷元則用0補足；同時舍棄缺失歷元數過半的數據，最終獲得166 d的可用數據，每天有86 400個觀測歷元，總數據量為13.77 GB。

3.2 反演精度分析

利用本文算法處理GRACE-FO衛星2020-01-01～06-30共166 d的可用數據，反演獲得120階重力場模型GM-GraceFO2020h。為檢驗反演精度，將ITU_GRACE16、HUST-Grace2016s、Tongji-Grace02s、GGM05S等GRACE重力場模型作為對比組，與GM-GraceFO2020h進行對比分析。同時選取2 190階模型EIGEN-6C4(該模型綜合了LAGEOS-1/2、GRACE、GOCE以及地形數據，具有較高的精度[20])作為基準，對比檢驗模型內符合精度。

用GM-GraceFO2020h與4個對比組模型(截斷至120階)分別計算全球1°×1°高程異常與重力異常格網，并求取各自相較于EIGEN-6C4的差值，統計結果如表1所示。

表1 各重力場模型與EIGEN-6C4模型間高程異常、重力異常差值

由表1可知，GM-GraceFO2020h與EIGEN-6C4的模型高程異常、重力異常差值相較于對比組，在數值上相近，差值最小值的差異分別在7 cm和1.2 mGal以內，最大值的差異分別在15 cm和11 mGal以內，平均值和標準差的差異分別在0.1 cm、0.001 mGal和0.2 cm、0.003 mGal以內。圖2為GM-GraceFO2020h與對比組模型相較于EIGEN-6C4的階方差。由圖2可見，GM-GraceFO2020h與GRACE重力場模型在各階位系數上的差異較小，總體上精度相近。圖3為GM- GraceFO2020h階方差曲線，其與Kaula曲線基本吻合，說明該模型與Kaula準則符合度較高，具有較高的內符合精度。圖4為該模型的重力異常、緯向垂線偏差和高程異常。

圖2 各模型與EIGEN-6C4位系數差的階方差Fig.2 Difference degree with EIGEN-6C4

圖3 GM-GraceFO2020h模型階方差Fig.3 Degree variance of GM-GraceFO2020h

圖4 GM-GraceFO2020h模型重力異常、 緯向垂線偏差、高程異常Fig.4 Gravity anomaly, deflection of the vertical and height anomaly derived from GM-GraceFO2020h

3.3 計算效率分析

從觀測方程構建、法矩陣計算、法方程解算3個方面對快速反演算法的計算效率進行分析。同時，引入并行加速比Sn[21]作為量化加速效果的參數：

(12)

式中，n表示參與運算的線程數，tn表示有n個線程參與時完成運算所需的時間。

3.3.1 觀測方程構建

在對數據進行60 s間隔采樣條件下，分別利用串行算法與快速反演算法完成如下計算任務：1)處理衛星23 d數據(共29 011條記錄)，構建60階模型3 716個待求參數的觀測方程；2)處理75 d數據(共100 960條記錄)，構建96階模型9 404個待求參數的觀測方程；3)處理166 d數據(共225 341條記錄)，構建120階模型14 636個待求參數的觀測方程。表2為各任務耗時情況。

表2 不同數據規模的計算任務耗時對比

對比結果顯示，傳統串行算法耗時隨著計算量的增大呈線性遞增趨勢；而本文提出的算法在構建觀測方法時，相較于串行算法，能有效提高計算效率，其并行加速比達到20以上，但隨著計算量的增加，加速比會逐漸減小到10左右。

3.3.2 法矩陣計算

分別用串行算法、OpenMP算法[11]、快速反演算法計算956階(設計矩陣維度為8 862×956)、3 716階(設計矩陣維度為29 011×3 716)、9 404階(設計矩陣維度為100 960×9 404)、14 636階(設計矩陣維度為225 341×14 636)法矩陣，時間消耗情況如表3，內存空間消耗情況如圖5。

表3 不同階法矩陣計算耗時情況

圖5 內存空間消耗情況Fig.5 Memory consumption

結合表3和圖5可見，快速反演算法完成各計算任務所需時間分別是串行算法的7.660%、2.753%、1.763%和1.537%，并行加速比在13～66之間，是OpenMP算法的52.078%、17.756%、11.869%和10.186%，其計算效率優于二者。本文算法的內存消耗是串行算法的1/20，較之小一個量級，也是OpenMP算法的1/2；同時，串行算法的內存消耗會隨著法矩陣和設計矩陣維度的增加而快速增大，而本文算法的增速較小。通過以上分析可得，快速反演算法綜合MKL庫與分區平差法二者優勢，相較于現有算法，在快速計算法矩陣的同時有效限制了對內存資源的消耗，大幅降低了模型反演對計算平臺硬件水平的要求，使得在內存資源有限的普通PC上處理大規模衛星重力數據成為可能。

3.3.3 法方程解算

分別利用串行算法(Gauss-Jordan法)、MKL算法[14](LU分解求逆法)、快速反演算法，求解不同維數的法方程，計算耗時情況如表4。

表4 不同維數法方程解算耗時對比

由表4可見，快速反演算法在法矩陣階數為116時的解算耗時略多于MKL算法0.004 s，其主要原因在于，計算量較小時內存讀寫耗時具有較大占比；而在其余情況下其耗時皆低于前兩種算法，說明該算法將MKL庫與預處理共軛梯度法相結合，能有效加快法方程的解算速度。但需要注意，共軛梯度法并不對法矩陣求逆，故需估計所求位系數的精度時，仍推薦使用MKL算法。

3.3.4 反演計算總體情況

分別利用串行算法和快速反演算法處理衛星166 d數據，反演120階重力場模型，計算耗時總體情況如表5所示。

表5 反演計算總體耗時

由表5可見，串行算法反演耗時超過19 h，而快速反演算法僅需約17 min，反演耗時縮減98.479%，并行加速比超過65。算法耗時最多的計算過程為法方程的構建，可采用僅計算法矩陣下三角部分、引入斯特拉森(Strassen)算法等方式來進一步優化算法，提高計算效率。

4 結 語

本文基于能量守恒方法，提出重力場反演快速異構并行算法。該算法實現了GPU并行計算設計矩陣與自由項向量，通過結合MKL庫與分區平差、預處理共軛梯度法來快速構建并解算法方程。利用本文算法處理GRACE-FO衛星2020-01-01～06-30期間數據，反演獲得120階重力場模型GM-GraceFO2020h。通過對比實驗，分析驗證了算法的反演精度與計算效率。

在反演精度方面，快速反演算法所得模型與現有GRACE重力場模型在前120階精度相當，具有較高的內符合精度。在計算效率方面，該算法在重力場反演各主要計算過程的計算效率都優于現有算法，相較于傳統串行算法，反演耗時縮減98.479%，內存消耗減小1個量級。綜上所述，本文提出的重力場反演快速異構并行算法，能在保證反演精度的同時，有效提高計算效率，減少反演耗時，降低硬件需求，可以為大規模處理衛星重力數據提供經濟高效的計算實施方案。

致謝：感謝德國GFZ ISDC數據中心提供GRACE-FO衛星軌道、姿態、激光干涉測距及加速度計數據。