基于OPENMP的高性能GNSS法方程解算方法研究*

楊 凱

(重慶市國(guó)土資源和房屋勘測(cè)規(guī)劃院，重慶 400020)

1 引言

目前，隨著全球或區(qū)域范圍內(nèi)大規(guī)模GNSS 基準(zhǔn)站網(wǎng)的不斷建設(shè)，基準(zhǔn)站數(shù)目不斷增加，如何高效地處理?yè)碛袛?shù)百、甚至上千個(gè)GNSS 基準(zhǔn)站的大規(guī)模基準(zhǔn)網(wǎng)數(shù)據(jù)，是一個(gè)亟需解決的重要問(wèn)題。大規(guī)模GNSS 基準(zhǔn)站網(wǎng)數(shù)據(jù)處理策略目前主要分為PPP解、非差網(wǎng)解和雙差網(wǎng)解模式等。全球公認(rèn)的高精度數(shù)據(jù)處理軟件包括GIPSY、GAMIT、BERNESE、EPOS 等，其中GAMIT 只能解算不超過(guò)100 個(gè)站的基準(zhǔn)站網(wǎng)數(shù)據(jù)［1］，BERNESE 和EPOS 等雖然可以同時(shí)解算超過(guò)100 個(gè)基準(zhǔn)站的觀測(cè)數(shù)據(jù)，但解算時(shí)間過(guò)長(zhǎng)、占用的計(jì)算機(jī)硬件資源過(guò)多，限制了其在大規(guī)模GNSS 基準(zhǔn)站網(wǎng)數(shù)據(jù)處理中的應(yīng)用［2］。

針對(duì)上述技術(shù)難點(diǎn)，本文利用IGS 全球連續(xù)跟蹤站的實(shí)際數(shù)據(jù)，論述并給出了解決同時(shí)處理的測(cè)站數(shù)量不能超過(guò)100 個(gè)的方法，比較分析了濾波算法與最小二乘估計(jì)算法的優(yōu)劣，對(duì)比分析了現(xiàn)有法方程解算方法的效率，實(shí)現(xiàn)了基于OPENMP 的并行喬里斯基分解法方程求逆方法，并通過(guò)實(shí)際算例驗(yàn)證了其解算的高效率。

2 測(cè)站數(shù)量限制的解決方法

首先，在軟件編寫(僅考慮FORTRAN 語(yǔ)言)時(shí)需要考慮到內(nèi)存溢出問(wèn)題。當(dāng)數(shù)組規(guī)模太大時(shí)很容易引起內(nèi)存溢出從而導(dǎo)致編譯失敗。因此，在測(cè)站數(shù)量超過(guò)100 個(gè)時(shí)，應(yīng)按照FORTRAN 95 格式編寫代碼，并使用FORTRAN 95 標(biāo)準(zhǔn)中規(guī)定的動(dòng)態(tài)數(shù)組。

其次，在FORTRAN 語(yǔ)言中對(duì)于文件的訪問(wèn)(包括讀/寫)是通過(guò)所謂“文件號(hào)”進(jìn)行的。若將同一個(gè)文件號(hào)賦予多個(gè)文件，會(huì)導(dǎo)致文件訪問(wèn)出錯(cuò)。特別是在大規(guī)模GNSS 基準(zhǔn)站網(wǎng)軟件的編寫中，訪問(wèn)的文件數(shù)目將達(dá)到數(shù)百甚至數(shù)千個(gè)，一旦出錯(cuò)會(huì)導(dǎo)致整個(gè)數(shù)據(jù)處理失敗。對(duì)此需采用取如下措施加以避免:1)將文件按類型分類，例如第一類可以是IGS精密星歷文件、鐘差文件、系統(tǒng)參數(shù)控制文件等，該類文件數(shù)目有限，可分配10 ～100 的文件號(hào);第二類文件可以是測(cè)站RINEX 格式觀測(cè)文件，該類文件數(shù)目較多，可以分配200 ～2 000(上限值按基準(zhǔn)站數(shù)目而定)的文件號(hào);第三類文件可以是測(cè)站信息文件，例如包含了測(cè)站周跳、野值信息的LOG 文件，可分配3 000 ～5 000 的文件號(hào)。2)在每一類文件中，可自行編寫一個(gè)獲得空閑文件號(hào)的子程序?qū)ξ募?hào)加以自動(dòng)分配，以防止文件號(hào)分配混亂。

3 現(xiàn)有參數(shù)估計(jì)方法比較分析

目前參數(shù)估計(jì)方法主要包括濾波算法和最小二乘估計(jì)算法兩大類。濾波算法一般用以解決包含了狀態(tài)參數(shù)的參數(shù)估計(jì)問(wèn)題，具有實(shí)時(shí)輸出狀態(tài)參數(shù)及其方差，同時(shí)占用計(jì)算機(jī)內(nèi)存資源較少的優(yōu)點(diǎn)。濾波算法中的卡爾曼濾波［3］在眾多應(yīng)用領(lǐng)域中得到了廣泛的應(yīng)用并出現(xiàn)了多種變體。特別的，在GNSS 精密定軌定位領(lǐng)域，均方根信息濾波(Square-Root Information Filter，簡(jiǎn)稱SRIF)算法［4］，已在美國(guó)噴氣動(dòng)力實(shí)驗(yàn)室研發(fā)的高精度GNSS 精密數(shù)據(jù)處理軟件GIPSY 中得到應(yīng)用［5］，在處理GNSS 觀測(cè)數(shù)據(jù)時(shí)能夠有效避免濾波發(fā)散，具有較高的數(shù)值穩(wěn)健性和高效性［6，7］。

對(duì)于大規(guī)模GNSS 基準(zhǔn)站網(wǎng)的數(shù)據(jù)處理來(lái)說(shuō)，雖然在解決動(dòng)態(tài)誤差參數(shù)估計(jì)的問(wèn)題上濾波算法(特別是均方根信息濾波)具有非常優(yōu)秀的表現(xiàn)，但GNSS 事后精密定位所關(guān)心的是待估參數(shù)多為固定偏差(例如測(cè)站坐標(biāo))，一般不關(guān)心中間過(guò)程不確定因素影響的狀態(tài)參數(shù)(例如接收機(jī)鐘差、對(duì)流層延遲參數(shù)、太陽(yáng)光壓參數(shù)等)求解。同時(shí)，濾波算法是逐歷元求解一次法方程，SRIF 則需要逐歷元進(jìn)行信息矩陣的更新，當(dāng)測(cè)站超過(guò)100 個(gè)時(shí)，每歷元進(jìn)行一次法方程求解或信息矩陣更新，將耗費(fèi)大量的計(jì)算時(shí)間。在GIPSY 軟件中，為了提高效率，將觀測(cè)數(shù)據(jù)的采樣間隔壓縮為每300 秒一個(gè)歷元［8］。

本文實(shí)現(xiàn)了利用SRIF 算法進(jìn)行參數(shù)估計(jì)，并采用2009年5月2日的IGS 全球連續(xù)跟蹤站網(wǎng)數(shù)據(jù)(采樣率為300s)分別在50、100、150 和200 站的不同測(cè)站數(shù)目情況下利用SRIF 進(jìn)行解算，各種情況下的測(cè)站分布如圖1 所示，得到的每歷元SRIF 解算所花費(fèi)時(shí)間如圖2 所示。

由圖2 可知，隨著測(cè)站數(shù)目的增加，SRIF 算法單歷元解算時(shí)間增長(zhǎng)非常迅速。當(dāng)測(cè)站數(shù)目為150個(gè)時(shí)，單歷元解算時(shí)間超過(guò)5 分鐘;而測(cè)站數(shù)目為200 個(gè)時(shí)，單歷元解算時(shí)間超過(guò)15 分鐘，由此推斷，解算采樣率為300s 的全天觀測(cè)數(shù)據(jù)時(shí)間將超過(guò)72小時(shí)。這對(duì)于大規(guī)模GNSS 基準(zhǔn)站網(wǎng)數(shù)據(jù)解算來(lái)說(shuō)，是難以容忍的。

解算是在作者的個(gè)人筆記本電腦上運(yùn)行的，計(jì)算機(jī)軟、硬件平臺(tái)信息如表1 所示。

表1 計(jì)算機(jī)軟、硬件平臺(tái)Tab.1 Computer hardware and software platform

另一方面，可以看到，由于最小二乘估計(jì)算法不需要逐歷元求解法方程，只需在觀測(cè)時(shí)段的每個(gè)歷元進(jìn)行法方程疊加，直到最后一個(gè)歷元才顯式求解GNSS 精密定位所關(guān)心的確定性參數(shù)和狀態(tài)參數(shù)，因此在計(jì)算效率上具有獨(dú)特的優(yōu)點(diǎn)。

4 現(xiàn)有法方程解算方法比較

法方程解算(主要是法方程求逆過(guò)程)是影響數(shù)據(jù)處理效率的關(guān)鍵因素之一。一般來(lái)說(shuō)，對(duì)于GNSS 精密數(shù)據(jù)處理，法方程中最多的參數(shù)是整周模糊度參數(shù)，其次為ZTD 參數(shù)，二者占所有參數(shù)的比例可能超過(guò)90%。法方程求逆過(guò)程是一個(gè)標(biāo)準(zhǔn)的線性代數(shù)問(wèn)題，可以用任意標(biāo)準(zhǔn)方法進(jìn)行解算。對(duì)于大規(guī)模GNSS 基準(zhǔn)站網(wǎng)的高性能數(shù)據(jù)處理，需要找出最優(yōu)法方程求逆算法。

另一方面，對(duì)于法方程求逆這樣一個(gè)數(shù)值計(jì)算問(wèn)題，不僅運(yùn)用不同的方法求逆效率有差異，即便同一個(gè)方法、不同編程實(shí)現(xiàn)方式的解算效率同樣可能相差較大。現(xiàn)有的高性能線性代數(shù)函數(shù)庫(kù)，包括BLAS、LAPACK、MKL 等，其運(yùn)算效率已經(jīng)被很多學(xué)者加以驗(yàn)證［9］。因此，使用這些已有的函數(shù)庫(kù)，可以安全可靠地讓測(cè)繪科技人員將主要精力放在與GNSS 相關(guān)的算法研究上，而不必在與計(jì)算機(jī)底層硬件相關(guān)的計(jì)算效率問(wèn)題上過(guò)分糾纏。

為了研究分析適用于最小二乘估計(jì)方法的最優(yōu)法方程解算方法，本文實(shí)現(xiàn)了利用高斯消元法、LU分解法、Bunch-Kaufman 選主元法和喬里斯基分解法來(lái)解算法方程，并采用2009年5月2日的IGS 全球連續(xù)跟蹤站網(wǎng)數(shù)據(jù)(采樣率為300s)分別在50、100、150 和200 站的不同測(cè)站數(shù)目情況下進(jìn)行解算，所用數(shù)據(jù)如圖1 所示，解算時(shí)的計(jì)算機(jī)軟、硬件平臺(tái)如表1 所示，采用不同方法解算法方程的效率如表2 所示。

表2 不同方法法方程求解時(shí)間對(duì)比(單位:s)Tab.2 Comparison among processing time with different normal equation methods(unit:s)

由表2 可知，對(duì)于最小二乘估計(jì)法方程的解算，高斯消元法耗時(shí)最多，接近LU 分解法所用時(shí)間的20 倍。LU 分解法、Bunch-Kaufman 選主元法和喬里斯基分解法所用時(shí)間基本在同一量級(jí)，喬里斯基分解法所用時(shí)間最少，約為L(zhǎng)U 分解法所用時(shí)間的一半。在200 個(gè)測(cè)站情況下，喬里斯基分解法解算時(shí)間694s，僅為同等情況下高斯消元法的2.4%。

究其原因，喬里斯基分解法的高效率，來(lái)源于其與最小二乘估計(jì)法方程特征的匹配。最小二乘估計(jì)法方程的系數(shù)矩陣，是一個(gè)實(shí)正定對(duì)稱矩陣，喬里斯基分解法可以充分利用其正定對(duì)稱特性，而LU 分解法(適用于所有可逆實(shí)矩陣)和Bunch-Kaufman 選主元法(適用于所有可逆實(shí)對(duì)稱)則不能充分利用這一特征，因此喬里斯基分解法的解算效率相對(duì)最高。

5 基于OPENMP 的高性能并行計(jì)算

值得注意的是數(shù)據(jù)處理僅占用了計(jì)算機(jī)處理器的一個(gè)線程，即所謂的單線程解算。然而，目前的計(jì)算機(jī)技術(shù)已經(jīng)進(jìn)入多核處理器時(shí)代，Intel 和AMD公司面向個(gè)人消費(fèi)者推出了眾多款式的多核桌面/移動(dòng)處理器。因此，在進(jìn)行大規(guī)模GNSS 基準(zhǔn)站網(wǎng)的高性能數(shù)據(jù)處理時(shí)，應(yīng)考慮同時(shí)使用多個(gè)線程進(jìn)行并行計(jì)算，以提升解算效率。

一般來(lái)說(shuō)，并行計(jì)算是由并行計(jì)算機(jī)來(lái)實(shí)現(xiàn)的，目前主要的并行計(jì)算機(jī)包括多計(jì)算機(jī)系統(tǒng)和集中式多處理器系統(tǒng)。并行算法在并行計(jì)算機(jī)上需要通過(guò)并行編程來(lái)實(shí)現(xiàn)，當(dāng)前使用較多的并行編程環(huán)境主要有消息傳遞編程接口(Message Passing Interface，簡(jiǎn)稱MPI)以及共享存儲(chǔ)編程接口(Open Multi-Processing，簡(jiǎn)稱OPENMP)兩種［10］。

其中，MPI 更適合于多臺(tái)計(jì)算機(jī)節(jié)點(diǎn)組成的計(jì)算機(jī)集群系統(tǒng)，而OPENMP 更適用于多個(gè)計(jì)算機(jī)內(nèi)核/線程的單臺(tái)計(jì)算機(jī)(節(jié)點(diǎn))。本文暫不涉及基于MPI的并行計(jì)算。在筆者的個(gè)人筆記本計(jì)算機(jī)上實(shí)現(xiàn)了基于LAPACK/MKL 的OPENMP 并行化計(jì)算，所用的計(jì)算機(jī)軟、硬件平臺(tái)如表1 所示。其實(shí)現(xiàn)方式包括:1)在Intel Fortran 11.2 編譯器中打開(kāi)-openmp 選項(xiàng)以支持OPENMP 并行化計(jì)算;2)對(duì)上述四種最小二乘法方程解算函數(shù)源代碼做適當(dāng)修改以實(shí)現(xiàn)邏輯上的并行運(yùn)算，具體步驟限于篇幅暫不贅述。

在將上文中的四種最小二乘法方程求解方法進(jìn)行并行化改造之后，采用2009年5月2日的全球IGS 連續(xù)運(yùn)行跟蹤站數(shù)據(jù)進(jìn)行了算例分析，所用測(cè)站數(shù)目分別為50 站、100 站、150 站與200 站，采樣率為300s。測(cè)站分布如圖1 所示。喬里斯基分解法在并行化之后的解算效率如圖3 所示，解算時(shí)同時(shí)使用了四個(gè)線程進(jìn)行并行計(jì)算。

圖3 法方程解算時(shí)間對(duì)比圖Fig.3 Comparison chart of normal equation computing time(unit:s)

結(jié)合表2 和圖3 可知，在實(shí)現(xiàn)了并行化計(jì)算后，法方程解算效率提升超過(guò)400%，在200 個(gè)測(cè)站情況下，喬里斯基分解法求解法方程用時(shí)130s，僅為串行計(jì)算時(shí)694s 解算時(shí)間的18%，可見(jiàn)并行計(jì)算對(duì)解算效率的巨大影響。

6 結(jié)論

1)對(duì)于大規(guī)模GNSS 基準(zhǔn)站網(wǎng)全網(wǎng)整體處理時(shí)，不能同時(shí)解算100 個(gè)測(cè)站以上數(shù)據(jù)的問(wèn)題可以通過(guò)編程技巧解決;

2)大規(guī)模GNSS 基準(zhǔn)站網(wǎng)的事后靜態(tài)數(shù)據(jù)處理，采用最小二乘算法進(jìn)行參數(shù)估計(jì)在解算效率上具有獨(dú)特的優(yōu)勢(shì);

3)法方程的解算，建議采用BLAS/LAPACK/MKL 等高性能線性代數(shù)函數(shù)庫(kù)以提升解算效率;

4)基于OPENMP 的并行化喬里斯基分解方法，可充分利用多核/多線程處理器的優(yōu)勢(shì)進(jìn)行法方程求逆解算，從而可以大幅提高個(gè)人計(jì)算機(jī)平臺(tái)的解算效率。

1 Herring T A，King R W and McClusky S C.GAMIT reference manual［EB/OL］.2009［2011－02－10］.http://www－gpsg.mit.edu/ ～simon/gtgk/docs.htm.

2 Ge M，et al.A new data processing strategy for huge GNSS global networks［J］.Journal of Geodesy，2006，80(4):199－203.

3 Kalman R E.A new approach to linear filtering and prediction problems［J］.Journal of Basic Engineering，1960，82(1):35－45.

4 Bierman G J.The treatment of bias in the square-root information filter/smoother［J］.Journal of Optimization Theory and Applications，1975，16(1):165－178.

5 Webb F H and Zumberge J F.An Introduction to GIPsy/oasIs-Ⅱ［R］.JPL Publ D-11088，Jet Propul Lab.，Pasadena，California，1993.

6 趙齊樂(lè)，等.均方根信息濾波和平滑及其在低軌衛(wèi)星星載GPS 精密定軌中的應(yīng)用［J］.武漢大學(xué)學(xué)報(bào)(信息科學(xué)版)，2006，(1):12－15.(Zhao Qile，et al.Applications of square-root information filtering and smoothing on orbit determination of LEO satellites with on bord GPS data［J］.Geomatics and Information Science of Wuhan University，2006，(1):12－15.)

7 施闖，等.衛(wèi)星導(dǎo)航系統(tǒng)綜合分析處理軟件PANDA 及研究進(jìn) 展［J］.航天 器 工 程，2009，(4):64－70.(Shi Chuang，et al.PANDA:Comprehensive processing software for satellite navigation systems and its research progress［J］.Spacecraft Egineering，2009，(4):64－70.)

8 Lichten S M.Towards GPS orbit accuracy of tens of centimeters［J］.Geophysical Research Letters，1990，17(3):215－218.

9 K Gstr M B and Poromaa P.LAPACK-style algorithms and software for solving the generalized Sylvester equation and estimating the separation between regular matrix pairs［J］.ACM Transactions on Mathematical Software，1996，22(1):78－103.

10 Quinn M J.Parallel programming in C with MPI ＆ openMP［M］.New York:McGraw-Hill Press，2003.