楊 凱

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

1 引言

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

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

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

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

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

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

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

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

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

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

解算是在作者的個人筆記本電腦上運行的，計算機軟、硬件平臺信息如表1 所示。

表1 計算機軟、硬件平臺Tab.1 Computer hardware and software platform

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

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

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

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

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

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

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

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

5 基于OPENMP 的高性能并行計算

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

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

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

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

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

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

6 結(jié)論

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

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

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

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

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 趙齊樂，等.均方根信息濾波和平滑及其在低軌衛(wèi)星星載GPS 精密定軌中的應用［J］.武漢大學學報(信息科學版)，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)星導航系統(tǒng)綜合分析處理軟件PANDA 及研究進 展［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.