梁 偉，徐新禹,2，李建成,2，朱廣彬

1. 武漢大學(xué)測(cè)繪學(xué)院，湖北 武漢 430079； 2. 武漢大學(xué)地球空間環(huán)境與大地測(cè)量教育部重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430079； 3. 國(guó)家測(cè)繪地理信息局衛(wèi)星測(cè)繪應(yīng)用中心，北京 100048

構(gòu)建超高階地球重力場(chǎng)模型一直是物理大地測(cè)量學(xué)領(lǐng)域的熱點(diǎn)研究問(wèn)題。采用重力梯度測(cè)量技術(shù)的GOCE重力衛(wèi)星于2009年成功發(fā)射，僅利用其觀測(cè)數(shù)據(jù)恢復(fù)的全球重力場(chǎng)模型階次達(dá)到了280階[1]，模型中長(zhǎng)波分量的精度更是有了2個(gè)數(shù)量級(jí)的提高，尤其在地面重力數(shù)據(jù)空白區(qū)，相比由GRACE和地面觀測(cè)數(shù)據(jù)聯(lián)合反演的EIGEN-5C等模型，在100 km空間分辨率上，GOCE任務(wù)反演的模型精度有了明顯提升[2]。這一進(jìn)展使得聯(lián)合GOCE衛(wèi)星觀測(cè)數(shù)據(jù)和地面觀測(cè)數(shù)據(jù)獲得全球更高精度的超高階重力場(chǎng)模型成為可能。反演的超高階重力場(chǎng)模型可用于區(qū)域高精度大地水準(zhǔn)面的確定[3]、全球高程基準(zhǔn)的統(tǒng)一[4]、地球內(nèi)部結(jié)構(gòu)反演[5]等科學(xué)研究及應(yīng)用。

基于球面重力異常數(shù)據(jù)反演超高階全球重力場(chǎng)模型常用的方法主要有數(shù)值積分(numerical quadrature，NQ)方法、最小二乘配置(least squares collocation，LSC)方法和最小二乘(least squares，LS)方法[6]。NQ方法是基于面球諧函數(shù)的正交性，采用積分方法獨(dú)立求解模型的每個(gè)系數(shù)，其計(jì)算簡(jiǎn)單，但無(wú)法直接給出模型系數(shù)的驗(yàn)后方差。相比NQ方法，LS方法不僅可評(píng)估模型系數(shù)的精度，而且可以聯(lián)合解算多類(lèi)觀測(cè)數(shù)據(jù)。同時(shí)，在LS方法中，如果位系數(shù)按照特定的順序排列，法方程矩陣會(huì)具有塊對(duì)角特點(diǎn)，因此可以分塊求解位系數(shù)，這種方法稱(chēng)為塊對(duì)角最小二乘(block diagonal least squares，BDLS)方法[6]。該方法可避免超高維矩陣的構(gòu)建及求逆，計(jì)算量大大減少，已成為目前求解超高階重力場(chǎng)模型的主要方法，國(guó)際上精度得到認(rèn)可的EGM96[7]、EGM2008[8]、EIGEN-6C4[9]等模型在構(gòu)建時(shí)均使用了該方法。文獻(xiàn)[10]分析了BDLS方法在重力場(chǎng)模型確定中的應(yīng)用，并通過(guò)模擬試驗(yàn)比較了BDLS方法與NQ方法求解360階重力場(chǎng)模型的精度。

考慮到衛(wèi)星重力觀測(cè)數(shù)據(jù)和重力異常數(shù)據(jù)在反演全球重力場(chǎng)模型時(shí)存在頻譜互補(bǔ)性，前者對(duì)重力場(chǎng)的中長(zhǎng)波信號(hào)更為敏感，后者對(duì)短波及甚短波信號(hào)更為敏感，因此聯(lián)合衛(wèi)星重力數(shù)據(jù)和重力異常數(shù)據(jù)可以反演得到高精度高分辨率的重力場(chǎng)模型，這也是當(dāng)前求解超高階重力場(chǎng)模型的主要方法。EGM2008模型構(gòu)建時(shí)聯(lián)合了GRACE衛(wèi)星觀測(cè)法方程和重力異常數(shù)據(jù)[8]，EIGEN-6C4模型不僅使用了這些數(shù)據(jù)，還使用了GOCE觀測(cè)數(shù)據(jù)，相比EGM2008模型在一些區(qū)域的精度有明顯提高[9]。文獻(xiàn)[11]聯(lián)合EGM2008模型和GOCE-TIM-R5模型解算了2160階次重力場(chǎng)模型GECO；文獻(xiàn)[12]采用全矩陣求逆的最小二乘方法解算了720階次的GOCO05c模型，模型精度在15′的空間分辨率上(720階次)與EIGEN-6C4相當(dāng)；文獻(xiàn)[13]基于衛(wèi)星、陸地、海洋等觀測(cè)數(shù)據(jù)，使用球諧分析方法確定了2159階次的UGM08模型；文獻(xiàn)[14]采用譜權(quán)組合方法聯(lián)合EIGEN-6S衛(wèi)星重力場(chǎng)模型和格網(wǎng)重力異常構(gòu)建了2160階次的重力場(chǎng)模型。

本文將研究BDLS方法的基本原理，并將OpenMP并行計(jì)算技術(shù)引入到BDLS方法的求解中，研制相應(yīng)的軟件模塊，并基于數(shù)值模擬試驗(yàn)，分析BDLS方法求解2160階重力場(chǎng)模型的精度及軟件模塊的可靠性。然后聯(lián)合由GOCE任務(wù)獨(dú)立建立的衛(wèi)星觀測(cè)法方程和EGM2008模型格網(wǎng)重力異常法方程，采用最小二乘方法估計(jì)一個(gè)2159階次的重力場(chǎng)模型，并利用獨(dú)立的GPS/水準(zhǔn)數(shù)據(jù)對(duì)模型精度進(jìn)行分析與檢驗(yàn)。

1 基本原理

1.1 由重力異常數(shù)據(jù)快速求解模型系數(shù)的BDLS方法

格網(wǎng)平均重力異常觀測(cè)值的觀測(cè)方程為[6]

(1)

由式(1)，可以建立如下誤差方程

(2)

(3)

使用LS方法構(gòu)建法方程時(shí)，當(dāng)重力異常數(shù)據(jù)及其精度的空間分布滿(mǎn)足下面4個(gè)條件時(shí)，法方程為稀疏矩陣，如按照一定順序排列待求解參數(shù)，法方程矩陣N為塊對(duì)角形式[16]：

(1) 格網(wǎng)數(shù)據(jù)分布于旋轉(zhuǎn)曲面上(例如球面上)；

(2) 格網(wǎng)數(shù)據(jù)覆蓋整個(gè)曲面且經(jīng)度方向的分辨率一致；

(3) 格網(wǎng)數(shù)據(jù)的精度與經(jīng)度無(wú)關(guān)，即數(shù)據(jù)的權(quán)重不依賴(lài)于經(jīng)度的變化；

(4) 格網(wǎng)數(shù)據(jù)的精度關(guān)于赤道對(duì)稱(chēng)，即緯度θ與-θ處的數(shù)據(jù)精度相同。

灰色部分代表非零元素，其余為零元素圖1 法方程矩陣的結(jié)構(gòu)示意圖Fig.1 The schematic diagram of the normal matrix

如果法方程矩陣為塊對(duì)角形式，求解參數(shù)時(shí)可以不構(gòu)建和求解整個(gè)法方程，而是根據(jù)法方程具體的分塊結(jié)構(gòu)，單獨(dú)構(gòu)建法方程矩陣N中的每個(gè)塊對(duì)角部分及與其對(duì)應(yīng)的U，形成子法方程，對(duì)其單獨(dú)求解，所有子法方程解的組合與原法方程的解完全一致。這種求解參數(shù)的方法稱(chēng)為塊對(duì)角最小二乘(BDLS)方法。求解超高階重力場(chǎng)模型時(shí)，參數(shù)個(gè)數(shù)較多，法方程龐大，現(xiàn)有的計(jì)算條件很難滿(mǎn)足嚴(yán)格LS方法的計(jì)算需求，而B(niǎo)DLS方法可以避免大型矩陣的構(gòu)建及求逆，計(jì)算量遠(yuǎn)小于LS方法。

值得注意的是，對(duì)于實(shí)際重力觀測(cè)數(shù)據(jù)，可以對(duì)其進(jìn)行處理使其滿(mǎn)足條件(1)和(2)，但是處理后的數(shù)據(jù)的精度特點(diǎn)很難滿(mǎn)足條件(3)和(4)，此時(shí)法方程不是嚴(yán)格塊對(duì)角結(jié)構(gòu)，而是塊對(duì)角占優(yōu)的結(jié)構(gòu)。此時(shí)，如果僅取法方程的塊對(duì)角部分采用BDLS方法會(huì)帶來(lái)計(jì)算誤差。文獻(xiàn)[17]通過(guò)模擬試驗(yàn)探討了重力異常數(shù)據(jù)精度的空間分布不滿(mǎn)足條件(3)和(4)時(shí)BDLS方法的計(jì)算誤差。結(jié)果表明，當(dāng)全球重力異常數(shù)據(jù)精度的空間分布存在差異，且其差異達(dá)到10 mGal(1 mGal=10-5m/s2)的情況下，采用BDLS方法仍能以較高的精度求解重力場(chǎng)模型，因此本文將不考慮重力異常數(shù)據(jù)精度的空間差異性問(wèn)題。

1.2 多源重力觀測(cè)數(shù)據(jù)的最小二乘聯(lián)合求解

本文聯(lián)合重力異常數(shù)據(jù)和GOCE衛(wèi)星重力數(shù)據(jù)的方法是最小二乘聯(lián)合平差方法[18-19]，它的基本原理如下：

可以分別根據(jù)觀測(cè)數(shù)據(jù)L1和L2的觀測(cè)模型建立它們的誤差方程如式(4)所示

(4)

在最小二乘準(zhǔn)則下可以建立數(shù)據(jù)L1和L2的法方程如式(5)所示，求解法方程可以得到待求參數(shù)

(5)

式中，P1、P2分別為數(shù)據(jù)L1和L2的權(quán)陣。

同時(shí)也可聯(lián)合數(shù)據(jù)L1和L2，建立它們聯(lián)合觀測(cè)的誤差方程如下

(6)

式中，vC=[v1v2]T；AC=[A1A2]T；lC=[l1l2]T。

在最小二乘準(zhǔn)則下可以建立它們的聯(lián)合觀測(cè)法方程，如式(7)所示

(7)

式中，PC為數(shù)據(jù)的權(quán)陣。

(8)

由式(8)可知，當(dāng)數(shù)據(jù)L1和L2不相關(guān)時(shí)，使用最小二乘聯(lián)合平差方法聯(lián)合數(shù)據(jù)L1和L2構(gòu)建的法方程和直接將各類(lèi)觀測(cè)數(shù)據(jù)構(gòu)建的法方程按照加權(quán)因子相加形成的法方程相同[18-19]。即當(dāng)觀測(cè)數(shù)據(jù)不相關(guān)時(shí)，直接把法方程加權(quán)相加就可實(shí)現(xiàn)多源數(shù)據(jù)的最小二乘聯(lián)合處理，這種方法簡(jiǎn)單有效。

在最小二乘聯(lián)合平差中的一種特殊情況是數(shù)據(jù)L1和L2求解的參數(shù)x1和x2可能不完全一樣，它們構(gòu)建的法方程的參數(shù)和維數(shù)不同，這樣不能直接將法方程相加。此時(shí)可以分別由數(shù)據(jù)L1和L2構(gòu)建參數(shù)x1和x2的參數(shù)集合x(chóng)C的誤差方程。此時(shí)建立的參數(shù)xC的誤差方程與原來(lái)參數(shù)x1和x2的誤差方程的區(qū)別在于設(shè)計(jì)矩陣中增加了0元素，這些0元素是與觀測(cè)量無(wú)關(guān)的參數(shù)的系數(shù)。進(jìn)而建立參數(shù)xC的法方程，此時(shí)將法方程相加可以得到它們對(duì)參數(shù)xC的聯(lián)合觀測(cè)法方程。

2 基于OpenMP技術(shù)的塊對(duì)角最小二乘軟件模塊設(shè)計(jì)及驗(yàn)證

構(gòu)建超高階重力場(chǎng)模型需要求解的參數(shù)較多，雖然BDLS方法的計(jì)算量遠(yuǎn)小于LS方法，但是一般計(jì)算機(jī)在單核單進(jìn)程下使用BDLS方法求解2160階重力場(chǎng)模型的時(shí)耗仍很長(zhǎng)。為提高BDLS方法的計(jì)算效率，本文將OpenMP并行計(jì)算引入BDLS方法的解算中，在高性能計(jì)算平臺(tái)上實(shí)現(xiàn)超高階重力場(chǎng)模型的快速解算。

并行計(jì)算是指同時(shí)使用多種計(jì)算資源完成計(jì)算任務(wù)，是縮短計(jì)算時(shí)間的有效手段。常用的并行計(jì)算有OpenMP多線(xiàn)程技術(shù)和MPI多進(jìn)程技術(shù)。OpenMP可為程序中的并行計(jì)算區(qū)域開(kāi)辟多個(gè)線(xiàn)程，每個(gè)線(xiàn)程獨(dú)立運(yùn)行，多個(gè)線(xiàn)程協(xié)同完成并行計(jì)算區(qū)域的計(jì)算任務(wù)[20-21]。

使用BDLS方法求解超高階重力場(chǎng)模型時(shí)，需要依次構(gòu)建法方程矩陣N中的塊對(duì)角部分和U中與其對(duì)應(yīng)的部分，并對(duì)其單獨(dú)求解。法方程的每個(gè)塊對(duì)角部分的構(gòu)建和求解可單獨(dú)進(jìn)行，在程序設(shè)計(jì)時(shí)一般用DO循環(huán)完成。OpenMP技術(shù)可以為相互獨(dú)立的DO循環(huán)開(kāi)辟多個(gè)線(xiàn)程，將多個(gè)循環(huán)分配到多個(gè)線(xiàn)程上運(yùn)行，充分利用高性能計(jì)算平臺(tái)的計(jì)算資源，以提高計(jì)算效率。

根據(jù)BDLS方法和OpenMP技術(shù)的特點(diǎn)，本文設(shè)計(jì)了求解超高階重力場(chǎng)模型的OpenMP并行計(jì)算軟件模塊。為驗(yàn)證BDLS方法和編寫(xiě)軟件模塊的正確性，設(shè)計(jì)了以下模擬試驗(yàn)：使用EGM2008模型(截?cái)嗟?159階)模擬計(jì)算了半徑為6 378 136.3 m球面上分辨率為5′×5′的格網(wǎng)重力異常平均值數(shù)據(jù)，基于該數(shù)據(jù)使用編寫(xiě)的BDLS并行計(jì)算軟件模塊求解了2159階的重力場(chǎng)模型，命名為T(mén)estModel，它的系數(shù)階誤差RMS如圖2所示(見(jiàn)后文)。試驗(yàn)使用的計(jì)算機(jī)為HP服務(wù)器，所用計(jì)算節(jié)點(diǎn)上有20個(gè)主頻為2.4 GHz的Intel Xeon CPU，在單節(jié)點(diǎn)上運(yùn)行，使用了80個(gè)線(xiàn)程，計(jì)算耗時(shí)7 h。從圖2可以看出，求解的模型系數(shù)的誤差小于10-11量級(jí)，試驗(yàn)結(jié)果驗(yàn)證了BDLS方法和相應(yīng)軟件模塊的正確性，可以高效精確地求解超高階重力場(chǎng)模型。

3 SGG-UGM-1超高階重力場(chǎng)模型的求解

EGM2008重力場(chǎng)模型在構(gòu)建時(shí)使用了精度較高且覆蓋較全的地面重力數(shù)據(jù)集[8]，本文使用其計(jì)算球面上的全球重力異常格網(wǎng)平均值，并聯(lián)合該數(shù)據(jù)和GOCE衛(wèi)星數(shù)據(jù)解算了2159階次重力場(chǎng)模型SGG-UGM-1。

3.1 超高階重力場(chǎng)模型解算的數(shù)據(jù)處理流程

本文構(gòu)建超高階重力場(chǎng)模型的解算流程如圖3 所示，具體描述如下：

(1) 使用EGM2008模型計(jì)算參考球面(球半徑為6 378 136.3 m)上5′×5′格網(wǎng)分辨率的重力異常數(shù)據(jù)ΔgEGM；由該數(shù)據(jù)使用BDLS方法構(gòu)建并求解2159階塊對(duì)角法方程，得到2～2159階位系數(shù)，這部分系數(shù)命名為BDLS，使用BDLS方法構(gòu)建法方程時(shí)認(rèn)為格網(wǎng)數(shù)據(jù)的權(quán)相等。

(2) 根據(jù)步驟(1)求得系數(shù)中2～100階以及251～2159階系數(shù)計(jì)算重力異常，然后在數(shù)據(jù)ΔgEGM中減去這部分重力異常得到殘余重力異常數(shù)據(jù)ΔgR。ΔgR數(shù)據(jù)中僅有101～250階系數(shù)的貢獻(xiàn)，再使用ΔgR數(shù)據(jù)構(gòu)建101～250階系數(shù)的塊對(duì)角法方程。

(3) 使用GOCE衛(wèi)星的衛(wèi)星重力梯度(satellite gravity gradiometry，SGG)數(shù)據(jù)和衛(wèi)星跟蹤衛(wèi)星(satellite-to-satellite tracking，SST)數(shù)據(jù)構(gòu)建220階衛(wèi)星觀測(cè)法方程，并求解得到220階的衛(wèi)星重力場(chǎng)模型，其數(shù)據(jù)處理方法的描述見(jiàn)下節(jié)。

(4) 步驟(2)建立了101～250階系數(shù)的法方程，步驟(3)建立了2～220階系數(shù)的法方程，使用1.2節(jié)描述的聯(lián)合平差方法，聯(lián)合這兩個(gè)法方程得到2～250階系數(shù)的聯(lián)合觀測(cè)法方程。3.3節(jié)中有關(guān)于法方程聯(lián)合過(guò)程更為詳細(xì)的描述。

(5) 求解步驟(4)建立的2～250階系數(shù)聯(lián)合觀測(cè)法方程，可以得到2～250階系數(shù)，這部分系數(shù)與步驟(1)中求解的系數(shù)中251～2159階部分一同組成2～2159階重力場(chǎng)模型系數(shù)。

3.2 GOCE 衛(wèi)星觀測(cè)法方程的構(gòu)建

本文使用的GOCE衛(wèi)星重力觀測(cè)數(shù)據(jù)為2009年11月1日—2012年5月31日期間的SGG數(shù)據(jù)和2009年11月1日—2010年7月5日期間的SST數(shù)據(jù)，數(shù)據(jù)的采樣間隔為1 s。其中包括：衛(wèi)星重力梯度觀測(cè)數(shù)據(jù)、梯度儀坐標(biāo)系相對(duì)于慣性系之間的姿態(tài)數(shù)據(jù)(四元素)、精密衛(wèi)星軌道數(shù)據(jù)、地固系與慣性系之間的轉(zhuǎn)換矩陣、加速度觀測(cè)數(shù)據(jù)、運(yùn)動(dòng)學(xué)軌道的方差-協(xié)方差數(shù)據(jù)等[18]。

基于上面的觀測(cè)數(shù)據(jù)，建立GOCE衛(wèi)星觀測(cè)法方程的簡(jiǎn)要描述如下：首先，對(duì)原始數(shù)據(jù)進(jìn)行粗差探測(cè)與剔除、時(shí)變重力場(chǎng)信號(hào)的分離、數(shù)據(jù)內(nèi)插分段等預(yù)處理；直接利用沿軌的引力梯度對(duì)角線(xiàn)數(shù)據(jù)(Vxx，Vyy，Vzz)建立觀測(cè)方程及相應(yīng)的法方程，再使用方差分量估計(jì)方法估計(jì)引力梯度張量對(duì)角線(xiàn)分量各自的單位權(quán)中誤差，3個(gè)分量的單位權(quán)中誤差分別為2.5、3.3和4.3 mE?？紤]到引力梯度觀測(cè)值的有色噪聲特點(diǎn)，在最小二乘過(guò)程中引入通帶為0.005～0.041 Hz的ARMA(auto regressive moving-average)濾波器。采用點(diǎn)域加速度方法由PKI(precise kinematic orbit)軌道觀測(cè)值建立相應(yīng)的SST法方程，求解SST法方程得到觀測(cè)值殘差，并基于殘差估計(jì)了單位權(quán)方差?；赟GG各分量的單位權(quán)方差以及SST的單位權(quán)方差，得到了它們之間的加權(quán)因子，最后基于加權(quán)因子將SGG數(shù)據(jù)各分量的法方程以及SST法方程相加得到GOCE衛(wèi)星觀測(cè)法方程，并求解得到相應(yīng)的衛(wèi)星重力場(chǎng)模型，求解時(shí)采用Tikhonov正則化方法解決極空白問(wèn)題引起的法方程不穩(wěn)定性問(wèn)題。文獻(xiàn)[18]中有更為詳細(xì)的數(shù)據(jù)處理方法和過(guò)程的描述。

圖3 SGG-UGM-1模型的解算流程Fig.3 Calculation process of SGG-UGM-1 model

3.3 GOCE 衛(wèi)星觀測(cè)法方程和格網(wǎng)重力異常法方程的聯(lián)合解算

聯(lián)合求解GOCE衛(wèi)星法方程與重力異常塊對(duì)角法方程可以得到最優(yōu)的模型解，但是由于兩個(gè)法方程在形式、階次、頻譜精度上有較大差異，聯(lián)合解算時(shí)需要采取有效的策略。根據(jù)文獻(xiàn)[9]，本文制定的法方程聯(lián)合解算策略如圖4所示。

使用EGM2008模型的系數(shù)誤差信息，可以計(jì)算得到模型在2190階的重力異常累計(jì)誤差約為5 mGal。本文使用5 mGal作為重力異常數(shù)據(jù)的單位權(quán)中誤差，并基于該單位中誤差以及上文中GOCE 數(shù)據(jù)的單位權(quán)方差得到重力異常法方程與GOCE觀測(cè)法方程的加權(quán)因子。

聯(lián)合觀測(cè)法方程分為兩部分，其中251～2159階系數(shù)法方程與重力異常法方程中對(duì)應(yīng)部分相同，而2～250階系數(shù)法方程是使用1.2節(jié)中的最小二乘聯(lián)合平差方法聯(lián)合101～250階重力異常數(shù)據(jù)塊對(duì)角法方程和GOCE全階次法方程(220階)得到的。構(gòu)建重力異常數(shù)據(jù)101～250階塊對(duì)角法方程時(shí)需要在原始重力異常數(shù)據(jù)中移除2～100階及251～2159階系數(shù)的貢獻(xiàn)。

模型的220～250階系數(shù)是聯(lián)合GOCE數(shù)據(jù)和重力異常數(shù)據(jù)聯(lián)合求解的，這樣可以保證聯(lián)合模型的這部分系數(shù)及其驗(yàn)后方差的頻譜有更好的連續(xù)性。重力異常數(shù)據(jù)求解低階系數(shù)的精度較差，重力異常低階次系數(shù)法方程與衛(wèi)星法方程聯(lián)合可能會(huì)“污染”衛(wèi)星重力的低階系數(shù)，所以為減弱其影響，2～250階聯(lián)合法方程中只采用了重力異常的101～250階系數(shù)塊對(duì)角法方程。

4 SGG-UGM-1超高階重力場(chǎng)模型的精度分析

4.1 SGG-UGM-1模型的系數(shù)誤差

EIGEN-6C4模型是當(dāng)前精度較高的超高階重力場(chǎng)模型[11]。本文以此為參考模型，認(rèn)為該模型的系數(shù)為“真值”計(jì)算SGG-UGM-1、GOSG-EGM、EGM2008和EIGEN-6C2[22]模型的“誤差”，并對(duì)它們的誤差進(jìn)行對(duì)比分析，以在頻域內(nèi)分析SGG-UGM-1的特性，評(píng)價(jià)模型的精度與合理性。

GOSG-EGM模型是利用本文計(jì)算得到的220階衛(wèi)星重力場(chǎng)模型與EGM2008重力場(chǎng)模型系數(shù)直接拼接得到的，其2～220階系數(shù)來(lái)自衛(wèi)星重力場(chǎng)模型，221～2159階系數(shù)來(lái)自EGM2008模型。圖5給出了SGG-UGM-1、GOSG-EGM、EGM2008和EIGEN-6C2模型系數(shù)的階誤差RMS。

由圖5可知，在2～70階的頻率范圍內(nèi)，EIGEN-6C2的誤差最小，EGM2008次之，而SGG-UGM-1和GOSG-EGM較大。這是由于EIGEN-6C2與EIGEN-6C4的2～70階系數(shù)的計(jì)算方法和使用數(shù)據(jù)類(lèi)似，都是聯(lián)合LAGEOS、GRACE和GOCE 等數(shù)據(jù)求解的[9,22]，EGM2008的2～70階次系數(shù)主要是來(lái)自于GRACE觀測(cè)數(shù)據(jù)的貢獻(xiàn)；而SGG-UGM-1和GOSG-EGM的2～70階次系數(shù)主要由GOCE 數(shù)據(jù)求解，因此其與EIGEN-6C4模型的差異大于EGM2008、EIGEN-6C2。

圖2 模型的系數(shù)階誤差RMSFig.2 Degree error RMS of the model coefficients

圖4 SGG-UGM-1的法方程聯(lián)合解算示意圖Fig.4 The schematic diagram of normal equation combination in SGG-UGM-1

圖5 重力場(chǎng)模型的系數(shù)階誤差RMSFig.5 Degree error RMS of the gravity field models

在70～220階，SGG-UGM-1和EIGEN-6C2與EIGEN-6C4更為接近，主要原因是兩個(gè)模型在此頻段內(nèi)均有GOCE和地面重力數(shù)據(jù)的貢獻(xiàn)，而GOSG-EGM僅有GOCE衛(wèi)星的貢獻(xiàn)，從170階開(kāi)始它的誤差迅速增大。同時(shí)，在170～220階部分，因?yàn)镾GG-UGM-1是聯(lián)合GOCE和地面重力數(shù)據(jù)解算的，因此比GOSG-EGM模型的系數(shù)誤差平滑，在220階沒(méi)有出現(xiàn)突變現(xiàn)象。

圖6為SGG-UGM-1等重力場(chǎng)模型的大地水準(zhǔn)面階誤差，反映了系數(shù)的誤差。圖6中模型的大地水準(zhǔn)面系數(shù)誤差與系數(shù)相對(duì)參考模型的系數(shù)誤差呈現(xiàn)出了一定的相似性。由圖6可知，SGG-UGM-1的系數(shù)誤差在220階附近達(dá)到最大，然后開(kāi)始下降，它的誤差譜連續(xù)，沒(méi)有出現(xiàn)突變現(xiàn)象。相比EGM2008、GOCE-EGM、BDLS模型，SGG-UGM-1模型低階次的系數(shù)誤差得到了明顯改善，這也說(shuō)明了本文模型聯(lián)合求解方法的合理性。

圖6 重力場(chǎng)模型的大地水準(zhǔn)面階誤差Fig.6 Degree geoid error of the gravity field models

頻譜域的分析結(jié)果說(shuō)明，本文求解的SGG-UGM-1模型精度可靠，且相比EGM2008模型在220階次內(nèi)的系數(shù)精度有了明顯改善。

4.2 中國(guó)和美國(guó)GPS/水準(zhǔn)數(shù)據(jù)的檢核結(jié)果

為了分析SGG-UGM-1模型的外符合精度，本文使用中國(guó)區(qū)域的649個(gè)GPS/水準(zhǔn)點(diǎn)數(shù)據(jù)[23]和美國(guó)區(qū)域的6169個(gè)GPS/水準(zhǔn)點(diǎn)數(shù)據(jù)[24]對(duì)模型進(jìn)行外部檢核，并與GOSG-EGM、EGM2008、EIGEN-6C2及EIGEN-6C4模型的檢核結(jié)果進(jìn)行對(duì)比分析，結(jié)果見(jiàn)表1、表2。

表1中國(guó)區(qū)域的GPS/水準(zhǔn)數(shù)據(jù)檢核結(jié)果

Tab.1 Validation of the gravity field models using GPS-leveling data in China m

表2美國(guó)區(qū)域的GPS/水準(zhǔn)數(shù)據(jù)檢核結(jié)果

Tab.2 Validation of the gravity field models using GPS-leveling data in America m

如表1所示，在中國(guó)區(qū)域SGG-UGM-1大地水準(zhǔn)面的精度遠(yuǎn)優(yōu)于EGM2008，略?xún)?yōu)于EIGEN-6C2和GOSG-EGM模型，略遜于EIGEN-6C4模型。如表2所示，在美國(guó)區(qū)域SGG-UGM-1模型的精度最高，但其他幾個(gè)模型的精度基本相當(dāng)。SGG-UGM-1模型在美國(guó)和中國(guó)區(qū)域的精度均優(yōu)于EGM2008模型，這主要是因?yàn)镾GG-UGM-1模型使用了GOCE數(shù)據(jù)，GOCE數(shù)據(jù)提高了模型低階次系數(shù)的精度，這與上文的模型系數(shù)分析的結(jié)果一致。尤其與EGM2008模型在地面重力數(shù)據(jù)空白或稀疏區(qū)對(duì)比，SGG-UGM-1模型在這些區(qū)域(如中國(guó)大陸)的精度提升更為明顯[25]。SGG-UGM-1模型的大地水準(zhǔn)面在中國(guó)區(qū)域和美國(guó)區(qū)域的精度檢核結(jié)果均優(yōu)于GOSG-EGM，充分體現(xiàn)了最小二乘聯(lián)合平差方法的優(yōu)勢(shì)。

4.3 毛烏素測(cè)區(qū)航空重力數(shù)據(jù)的檢核結(jié)果

為了進(jìn)一步分析SGG-UGM-1模型的外符合精度，本文使用毛烏素測(cè)區(qū)的航空重力擾動(dòng)觀測(cè)數(shù)據(jù)對(duì)模型進(jìn)行外部檢核，并與GOSG-EGM、EGM2008、EIGEN-6C2及EIGEN-6C4模型的檢核結(jié)果進(jìn)行對(duì)比分析，結(jié)果如表3示。毛烏素測(cè)區(qū)位于陜北地區(qū)，測(cè)區(qū)范圍為38°48′N(xiāo)—39°16′N(xiāo)、109°17′E—110°27′E，測(cè)區(qū)大小為50 km×100 km。

表3航空重力數(shù)據(jù)的檢核結(jié)果

如表3所示，在毛烏素測(cè)區(qū)，SGG-UGM-1模型計(jì)算重力擾動(dòng)的精度略遜于EGM2008和EIGEN-6C4模型，優(yōu)于GOSG-EGM和EIGEN-6C2模型，總體來(lái)說(shuō)幾個(gè)模型的精度相當(dāng)。相比GPS水準(zhǔn)檢核的結(jié)果，雖然SGG-UGM-1模型的低階次系數(shù)相較EGM2008模型有了較大提升，但由于其高階次系數(shù)主要是由EGM2008模型重力異常求解的，且重力擾動(dòng)的誤差主要來(lái)自高階次系數(shù)，所以SGG-UGM-1航空重力擾動(dòng)的精度相較EGM2008沒(méi)有提升。

5 總 結(jié)

本文使用EGM2008模型的5′×5′格網(wǎng)重力異常構(gòu)建了2159階次的塊對(duì)角法方程，并使用最小二乘方法聯(lián)合GOCE重力衛(wèi)星任務(wù)220階次的全階次法方程求解了2159階次的重力場(chǎng)模型SGG-UGM-1，使用EIGEN-6C4等重力場(chǎng)模型、中國(guó)與美國(guó)的GPS水準(zhǔn)/數(shù)據(jù)和毛烏素測(cè)區(qū)的航空重力數(shù)據(jù)分析了SGG-UGM-1模型的內(nèi)外符合精度。主要的結(jié)論如下：

(1) BDLS方法可以快速精確地確定超高階重力場(chǎng)模型，采用OpenMP技術(shù)可以大幅度提高BDLS軟件模塊的計(jì)算效率。

(2) 采用最小二乘聯(lián)合平差方法不僅可以保證求解模型在聯(lián)合頻段內(nèi)信號(hào)和誤差頻譜的連續(xù)性，同時(shí)求解的模型SGG-UGM-1精度優(yōu)于直接拼接獲得的GOSG-EGM模型。

(3) 本文使用的GPS/水準(zhǔn)數(shù)據(jù)對(duì)模型的檢核結(jié)果表明，SGG-UGM-1大地水準(zhǔn)面在中國(guó)區(qū)域的精度介于EIGEN-6C2和EIGEN-6C4兩個(gè)模型之間，優(yōu)于GOSG-EGM模型，遠(yuǎn)優(yōu)于EGM2008模型；在美國(guó)區(qū)域這些模型的精度相當(dāng)。

(4) 本文使用的毛烏素測(cè)區(qū)航空重力數(shù)據(jù)對(duì)模型的檢核結(jié)果表明，SGG-UGM-1模型計(jì)算的重力擾動(dòng)精度與EGM2008、EIGEN-6C4模型相當(dāng)，優(yōu)于GOSG-EGM模型和EIGEN-6C2模型。

模型的檢核結(jié)果說(shuō)明本文構(gòu)建超高階重力場(chǎng)模型策略的合理性，為進(jìn)一步研究超高階重力場(chǎng)模型的構(gòu)建奠定了基礎(chǔ)。但在聯(lián)合衛(wèi)星觀測(cè)數(shù)據(jù)和重力異常數(shù)據(jù)時(shí)，本文參考了EIGEN-6C系列模型構(gòu)建的經(jīng)驗(yàn)，后續(xù)應(yīng)該進(jìn)一步深入研究重力異常法方程與衛(wèi)星觀測(cè)法方程的最優(yōu)聯(lián)合問(wèn)題。SGG-UGM-1模型的低階次系數(shù)的精度較EGM2008有了較大提升，但由于本文使用的重力異常數(shù)據(jù)是由EGM2008模型計(jì)算得到的，所以模型在高頻段沒(méi)有提升，且與EGM2008模型有很強(qiáng)的相關(guān)性。更為合理的做法是從全球高精度的地面重力、航空重力、衛(wèi)星測(cè)高、地形等觀測(cè)數(shù)據(jù)出發(fā)，完全獨(dú)立自主處理各類(lèi)數(shù)據(jù)得到全球高分辨率的重力異常格網(wǎng)數(shù)據(jù)集，然后據(jù)此反演超高階重力場(chǎng)模型，這也是筆者后續(xù)的研究重點(diǎn)。

參考文獻(xiàn)：

[1] BROCKMANN J M, ZEHENTNER N, H?CK E, et al. EGM_TIM_RL05: An Independent Geoid with Centimeter Accuracy Purely Based on the GOCE Mission[J]. Geophysical Research Letters, 2014, 41(22): 8089-8099.

[2] PAIL R, GOIGINGER H, SCHUH W D, et al. Combined Satellite Gravity Field Model GOCO01S Derived from GOCE and GRACE[J]. Geophysical Research Letters, 2010, 37(20): L20314.

[3] FEATHERSTONE W E. GNSS-based Heighting in Australia: Current, Emerging and Future Issues[J]. Journal of Spatial Science, 2008, 53(2): 115-133.

[4] RUMMEL R. Height Unification Using GOCE[J]. Journal of Geodetic Science, 2012, 2(4): 355-362.

[5] MCKENZIE D, YI Weiyong, RUMMEL R. Estimates ofTefrom GOCE Data[J]. Earth and Planetary Science Letters, 2014, 399(2): 116-127.

[6] PAVLIS N K. Modeling and Estimation of a Low Degree Geopotential Model from Terrestrial Gravity Data[R]. Report No.386. Columbus: The Ohio State University, 1988.

[7] LEMONIE F G, KENYON S C, FACTOR J K, et al. The Development of the Joint NASA GSFC and the National Imagery and Mapping Agency (NIMA) Geopotential Model EGM96[R]. Greenbelt: NASA Goddard Space Flight Center, 1998.

[8] PAVLIS N K, HOLMES S A, KENYON S C, et al. The Development and Evaluation of the Earth Gravitational Model 2008 (EGM2008)[J]. Journal of Geophysical Research, 2012, 117(3): 205-213.

[9] F?RSTE C, BRUINSMA S L, ABRIKOSOV O, et al. EIGEN-6C4 the Latest Combined Global Gravity Field Model Including GOCE Data up to Degree and Order 2190 of GFZ Potsdam and GRGS Toulouse[EB/OL]. http:∥doi.org/10.5880/icgem.2015.1.

[10] 李新星, 吳曉平, 李?yuàn)檴? 等. 塊對(duì)角最小二乘方法在確定全球重力場(chǎng)模型中的應(yīng)用[J]. 測(cè)繪學(xué)報(bào), 2014, 43(8): 778-785. DOI:10.13485/j.cnki.11-2089.2014.0110.

LI Xinxing, WU Xiaoping, LI Shanshan, et al. The Application of Block-diagonal Least-squares Methods in Geopotential Model Determination[J]. Acta Geodaetica et Cartographica Sinica, 2014, 43(8): 778-785. DOI:10.13485/j.cnki.11-2089.2014.0110.

[11] GILARDONI M, REGUZZONI M, SAMPIETRO D. GECO: A Global Gravity Model by Locally Combining GOCE Data and EGM2008[J]. Studia Geophysica et Geodaetica, 2016, 60(2): 228-247.

[12] FECHER T, PAIL R, GRUBER T, et al. GOCO05c: A New Combined Gravity Field Model Based on Full Normal Equations and Regionally Varying Weighting[J]. Surveys in Geophysics, 2017, 38(3): 571-590.

[13] 王正濤, 黨亞民, 晁定波. 超高階地球重力位模型確定的理論與方法[M]. 北京: 測(cè)繪出版社, 2011.

WANG Zhengtao, DANG Yamin, CHAO Dingbo. Theory and Methodology of Ultra-high-degree Geopotential Model Determination[M]. Beijing: Surveying and Mapping Press, 2011.

[14] 李新星. 超高階地球重力場(chǎng)模型的構(gòu)建[D]. 鄭州: 信息工程大學(xué), 2013.

LI Xinxing. Building of an Ultra-high-degree Geopotential Model[D]. Zhengzhou: Information Engineering University,

2013.

[15] 海斯卡涅W A, 莫里茲H. 物理大地測(cè)量學(xué)[M]. 盧?？? 胡國(guó)理, 譯. 北京: 測(cè)繪出版社, 1979.

HEISKANEN W A, MORITZ H. Physical Geodesy[M]. LU Fukang, HU Guoli, trans. Beijing: Surveying and Mapping Press, 1979.

[16] COLOMBO O L. Numerical Methods for Harmonic Analysis on the Sphere[R]. Columbus: The Ohio State University,

1981.

[17] LIANG Wei, LI Jiancheng, XU Xinyu, et al. Analysis of the Impact on the Gravity Field Determination from the Data with the Ununiform Noise Distribution Using Block-diagonal Least Squares Method[J]. Geodesy and Geodynamics, 2016, 7(3): 194-201.

[18] XU Xinyu, ZHAO Yongqi, REUBELT T, et al. A GOCE Only Gravity Model GOSG01S and the Validation of GOCE Related Satellite Gravity Models[J]. Geodesy and Geodynamics, 2017, 8(4): 260-272.

[19] 鐘波, 寧津生, 羅志才, 等. 聯(lián)合GOCE衛(wèi)星軌道和重力梯度數(shù)據(jù)嚴(yán)密求解重力場(chǎng)的模擬研究[J]. 武漢大學(xué)學(xué)報(bào)(信息科學(xué)版), 2012, 37(10): 1215-1220.

ZHONG Bo, NING Jinsheng, LUO Zhicai, et al. Simulation Study of Rigorous Gravity Field Recovery by Combining GOCE Satellite Orbit and Gravity Gradient Data[J]. Geomatics and Information Science of Wuhan University, 2012, 37(10): 1215-1220.

[20] 陳國(guó)良. 并行計(jì)算: 結(jié)構(gòu)·算法·編程[M]. 3版. 北京: 高等教育出版社, 2011.

CHEN Guoliang. Parallel Computing: Architecture, Algorithm and Programming[M]. 3rd ed. Beijing: Higher Education Press, 2011.

[21] 鄒賢才, 李建成, 汪海洪, 等. OpenMP并行計(jì)算在衛(wèi)星重力數(shù)據(jù)處理中的應(yīng)用[J]. 測(cè)繪學(xué)報(bào), 2010, 39(6): 636-641.

ZOU Xiancai, LI Jiancheng, WANG Haihong, et al. Application of Parallel Computing with OpenMP in Data Processing for Satellite Gravity[J]. Acta Geodaetica et Cartographica Sinica, 2010, 39(6): 636-641.

[22] F?ERSTE C, BRUINSMA S L, FLECHTNER F, et al. A Preliminary Update of the Direct Approach GOCE Processing and a New Release of EIGEN-6C[C]∥Proceedings of American Geophysical Union 2012 Fall Meeting. San Francisco, USA: American Geophysical Union, 2012.

[23] LI Jiancheng, JIANG Weiping, ZOU Xiancai, et al. Evaluation of Recent GRACE and GOCE Satellite Gravity Models and Combined Models Using GPS/leveling and Gravity Data in China[M]∥MARTI U. Gravity, Geoid and Height Systems. Cham, Switzerland: Springer, 2014: 67-74.

[24] MILBERT D G. Documentation for the GPS Benchmark Data Set of 23-July-98[R]. Milan: IGeS International Geoid Service, 1998: 29-42.

[25] RUMMEL R, YI Weiyong, STUMMER C. GOCE Gravitational Gradiometry[J]. Journal of Geodesy, 2011, 85(11): 777-790.