聯(lián)合EGM2008模型重力異常和GOCE觀測數(shù)據(jù)構(gòu)建超高階地球重力場模型SGG-UGM-1

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

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

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

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

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

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

1 基本原理

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

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

(1)

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

(2)

(3)

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

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

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

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

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

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

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

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

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

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

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

(4)

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

(5)

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

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

(6)

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

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

(7)

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

(8)

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

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

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

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

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

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

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

3 SGG-UGM-1超高階重力場模型的求解

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

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

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

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

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

(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)星觀測法方程，并求解得到220階的衛(wèi)星重力場模型，其數(shù)據(jù)處理方法的描述見下節(jié)。

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

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

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

本文使用的GOCE衛(wèi)星重力觀測數(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)星重力梯度觀測數(shù)據(jù)、梯度儀坐標(biāo)系相對于慣性系之間的姿態(tài)數(shù)據(jù)(四元素)、精密衛(wèi)星軌道數(shù)據(jù)、地固系與慣性系之間的轉(zhuǎn)換矩陣、加速度觀測數(shù)據(jù)、運動學(xué)軌道的方差-協(xié)方差數(shù)據(jù)等[18]。

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

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

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

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

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

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

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

4 SGG-UGM-1超高階重力場模型的精度分析

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

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

GOSG-EGM模型是利用本文計算得到的220階衛(wèi)星重力場模型與EGM2008重力場模型系數(shù)直接拼接得到的，其2～220階系數(shù)來自衛(wèi)星重力場模型，221～2159階系數(shù)來自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ù)的計算方法和使用數(shù)據(jù)類似，都是聯(lián)合LAGEOS、GRACE和GOCE 等數(shù)據(jù)求解的[9,22]，EGM2008的2～70階次系數(shù)主要是來自于GRACE觀測數(shù)據(jù)的貢獻；而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 重力場模型的系數(shù)階誤差RMSFig.5 Degree error RMS of the gravity field models

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 總 結(jié)

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

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

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

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

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

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

參考文獻：

[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] 李新星, 吳曉平, 李姍姍, 等. 塊對角最小二乘方法在確定全球重力場模型中的應(yīng)用[J]. 測繪學(xué)報, 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]. 北京: 測繪出版社, 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] 李新星. 超高階地球重力場模型的構(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. 物理大地測量學(xué)[M]. 盧?？? 胡國理, 譯. 北京: 測繪出版社, 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)密求解重力場的模擬研究[J]. 武漢大學(xué)學(xué)報(信息科學(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] 陳國良. 并行計算: 結(jié)構(gòu)·算法·編程[M]. 3版. 北京: 高等教育出版社, 2011.

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

[21] 鄒賢才, 李建成, 汪海洪, 等. OpenMP并行計算在衛(wèi)星重力數(shù)據(jù)處理中的應(yīng)用[J]. 測繪學(xué)報, 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.