...聯合EGM2008模型重力異常和GOCE觀測數據構建超高階地球重力場...

　　《測繪學報》

　　構建與學術的橋梁 拉近與權威的距離

　　聯合EGM2008模型重力異常和GOCE觀測數據構建超高階地球重力場模型SGG-UGM-1

　　梁偉1
, 徐新禹1,2

　　

　　 , 李建成1,2 , 朱廣彬3

　　1. 武漢大學測繪學院, 湖北 武漢 430079;
2. 武漢大學地球空間環境與大地測量教育部重點實驗室, 湖北 武漢 430079;
3. 國家測繪地理信息局衛星測繪應用中心, 北京 100048

　　收稿日期：2017-05-22；修回日期：2017-10-28

　　基金項目：國家自然科學基金（41774020；41210006；41404020）；國家高分專項高分遙感測繪應用示範系統項目（AH1701-4）

　　第一作者簡介：梁偉(1991—)，男，博士生，研究方向為地球重力場模型的構建。E-mail: [e-mail]

　　通信作者：徐新禹, E-mail: [e-mail]

　　摘要：本文研究了聯合衛星觀測數據和重力異常數據確定超高階重力場模型的理論方法，並使用EGM2008模型重力異常和GOCE（gravity field and ocean circulation explorer）觀測數據構建了重力場模型SGG-UGM-1。重點研究了由球面格網重力異常快速構建超高階重力場模型的塊對角最小二乘方法，將OpenMP技術引入到塊對角最小二乘中以提高計算效率，並基於模擬數據驗證了方法及算法和軟體模塊的正確性。採用本文制定的聯合解算策略，利用GOCE重力衛星觀測數據構建的220階次法方程和EGM2008模型重力異常構建的2159階次塊對角法方程，聯合求解了2159階次的重力場模型SGG-UGM-1。將SGG-UGM-1與EGM2008、EIGEN-6C2、EIGEN-6C4等超高階模型在頻譜域內進行了比較分析，結果表明SGG-UGM-1相對參考模型的係數誤差較小，且在220階次內的係數精度相比EGM2008模型有了提高。採用中國與美國的GPS/水準數據和毛烏素測區的航空重力觀測數據對這些模型進行了外符合精度的檢驗。檢核結果表明，在中國區域，SGG-UGM-1模型大地水準面的精度在EIGEN-6C2和EIGEN-6C4兩個模型之間，優於GOSG-EGM模型和EGM2008模型，與美國區域幾個模型的精度相當。利用毛烏素測區的航空重力數據對幾個模型進行了檢核，結果表明SGG-UGM-1模型計算的重力擾動精度與EGM2008、EIGEN-6C4模型相當，優於GOSG-EGM模型和EIGEN-6C2模型。

　　關鍵詞：SGG-UGM-1 超高階重力場模型 塊對角最小二乘方法 OpenMP並行計算

　　The Determination of an Ultra-high Gravity Field Model SGG-UGM-1 by Combining EGM2008 Gravity Anomaly and GOCE Observation Data

　　LIANG Wei1
, XU Xinyu1,2

　　

　　 , LI Jiancheng1,2 , ZHU Guangbin3

　　Abstract: The theory and methods of the determination of an ultra-high gravity field model by combination of satellite observation data and gravity anomaly data are studied.And an ultra-high gravity field model named SGG-UGM-1 is computed using EGM2008 derived gravity anomaly and GOCE observation data.The block-diagonal least squares (BDLS) method for quickly estimating an ultra-high gravity field model is researched and the corresponding software module is validated by numerical experiments.OpenMP technique is introduced into the BDLS program, which improves computing efficiency dramatically.An ultra-high gravity model SGG-UGM-1 complete to degree and order 2159 is derived using the proposed calculation strategies.The fully occupied normal equation system up to degree and order 220 formed by GOCE satellite data and the block-diagonal normal equation system up to degree and order 2159 formed by EGM2008 gravity anomaly data are used for the combination.Comparison among the models SGG-UGM-1 and EGM2008, EIGEN-6C2, EIGEN-6C4, GOSG-EGM in frequency domain has been done, which shows that SGG-UGM-1 is close to the reference models and that the coefficients lower than degree 220 of SGG-UGM-1 are more accurate than that of EGM2008.The models are also validated by GPS-leveling data in China and America and airborne gravity data in Maowusu surveying area.The results show that in China the accuracy level of SGG-UGM-1 derived geoid is between EIGEN-6C2 and EIGEN-6C4, and better than GOSG-EGM and EGM2008, while in America they are almost the same.In Maowusu area the accuracy level of SGG-UGM-1 derived gravity disturbance is almost at the same accuracy level with EGM2008 and EIGEN-6C4 and better than GOSG-EGM and EIGEN-6C2.

　　Key words: SGG-UGM-1 ultra-high gravity field model block diagonal least squares method OpenMP parallel computing technique

　　構建超高階地球重力場模型一直是物理大地測量學領域的熱點研究問題。採用重力梯度測量技術的GOCE重力衛星於2009年成功發射，僅利用其觀測數據恢復的全球重力場模型階次達到了280階[1]，模型中長波分量的精度更是有了2個數量級的提高，尤其在地面重力數據空白區，相比由GRACE和地面觀測數據聯合反演的EIGEN-5C等模型，在100 km空間解析度上，GOCE任務反演的模型精度有了明顯提升[2]。這一進展使得聯合GOCE衛星觀測數據和地面觀測數據獲得全球更高精度的超高階重力場模型成為可能。反演的超高階重力場模型可用於區域高精度大地水準面的確定[3]、全球高程基準的統一[4]、地球內部結構反演[5]等科學研究及應用。

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

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

　　本文將研究BDLS方法的基本原理，並將OpenMP並行計算技術引入到BDLS方法的求解中，研製相應的軟體模塊，並基於數值模擬試驗，分析BDLS方法求解2160階重力場模型的精度及軟體模塊的可靠性。然後聯合由GOCE任務獨立建立的衛星觀測法方程和EGM2008模型格網重力異常法方程，採用最小二乘方法估計一個2159階次的重力場模型，並利用獨立的GPS/水準數據對模型精度進行分析與檢驗。

　　1 基本原理1.1 由重力異常數據快速求解模型係數的BDLS方法

　　格網平均重力異常觀測值的觀測方程為[6]

　　
(1)

　　式中，Δσi=Δλ(cos θi-cos θi+1)；

　　

　　 ；(i, j)表示第i行j列個格網；(r, θ, λ)為球坐標[15]；GM為地球引力常數；R為參考球半徑；rij為r在格網上的平均值；Pnm(cos θ)為規格化的締合勒讓德函數；IPnmi為Pnm(cos θ)sin θ的積分值；(Cnm, Cnm)為擾動位球諧係數；(n, m)為勒讓德函數及模型係數的階次；Δgij為格網重力異常平均值。

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

　　
(2)

　　式中，L是觀測值Δgij；v是觀測值的改正數；A是基於式(1)建立的設計矩陣； 是待求的擾動位球諧係數(Cnm, Snm)。

　　在LS平差準則下，可得式(3)，其中P是觀測值Δgij的權陣，式(3)即為使用LS方法求解重力場模型的基本公式

　　
(3)

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

　　(1) 格網數據分布於旋轉曲面上(例如球面上)；

　　(2) 格網數據覆蓋整個曲面且經度方向的解析度一致；

　　(3) 格網數據的精度與經度無關，即數據的權重不依賴於經度的變化；

　　(4) 格網數據的精度關於赤道對稱，即緯度θ與-θ處的數據精度相同。

　　例如當係數先按次排列，然後同次的係數中Cnm係數在前，Snm係數在後，同次的Cnm和Snm均按階排列，則法方程將是嚴格的塊對角結構。例如：6階重力場模型係數的法方程矩陣結構如圖 1所示。

　　

圖 1 法方程矩陣的結構示意圖Fig. 1 The schematic diagram of the normal matrix灰色部分代表非零元素，其餘為零元素

　　圖選項

　　

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

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

　　1.2 多源重力觀測數據的最小二乘聯合求解

　　本文聯合重力異常數據和GOCE衛星重力數據的方法是最小二乘聯合平差方法[18-19]，它的基本原理如下：

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

　　
(4)

　　式中，v1、v2為觀測量的改正數；A1、A2為設計矩陣、 為待求參數；l1、l2為觀測量。

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

　　
(5)

　　式中，P1、P2分別為數據L1和L2的權陣。

　　同時也可聯合數據L1和L2，建立它們聯合觀測的誤差方程如下

　　
(6)

　　式中，vC=[v1 v2]T；AC=[A1 A2]T；lC=[l1 l2]T。

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

　　
(7)

　　式中，PC為數據的權陣。

　　當數據L1和L2不相關時，式(7)中權陣PC= ，其中w1和w2分別為數據L1和L2的加權因子，此時有[18-19]

　　
(8)

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

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

　　2 基於OpenMP技術的塊對角最小二乘軟體模塊設計及驗證

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

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

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

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

　　

圖 2 模型的係數階誤差RMSFig. 2 Degree error RMS of the model coefficients

　　圖選項

　　

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

　　EGM2008重力場模型在構建時使用了精度較高且覆蓋較全的地面重力數據集[8]，本文使用其計算球面上的全球重力異常格網平均值，並聯合該數據和GOCE衛星數據解算了2159階次重力場模型SGG-UGM-1。

　　3.1 超高階重力場模型解算的數據處理流程

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

　　

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

　　圖選項

　　

　　(1) 使用EGM2008模型計算參考球面(球半徑為6 378 136.3 m)上5′×5′格網解析度的重力異常數據ΔgEGM；由該數據使用BDLS方法構建並求解2159階塊對角法方程，得到2~2159階位係數，這部分係數命名為BDLS，使用BDLS方法構建法方程時認為格網數據的權相等。

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

　　(3) 使用GOCE衛星的衛星重力梯度(satellite gravity gradiometry，SGG)數據和衛星跟蹤衛星(satellite-to-satellite tracking，SST)數據構建220階衛星觀測法方程，並求解得到220階的衛星重力場模型，其數據處理方法的描述見下節。

　　(4) 步驟(2)建立了101~250階係數的法方程，步驟(3)建立了2~220階係數的法方程，使用1.2節描述的聯合平差方法，聯合這兩個法方程得到2~250階係數的聯合觀測法方程。3.3節中有關於法方程聯合過程更為詳細的描述。

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

　　3.2 GOCE衛星觀測法方程的構建

　　本文使用的GOCE衛星重力觀測數據為2009年11月1日—2012年5月31日期間的SGG數據和2009年11月1日—2010年7月5日期間的SST數據，數據的採樣間隔為1 s。其中包括：衛星重力梯度觀測數據、梯度儀坐標系相對於慣性系之間的姿態數據(四元素)、精密衛星軌道數據、地固系與慣性系之間的轉換矩陣、加速度觀測數據、運動學軌道的方差-協方差數據等[18]。

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

　　3.3 GOCE衛星觀測法方程和格網重力異常法方程的聯合解算

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

　　

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

　　圖選項

　　

　　使用EGM2008模型的係數誤差信息，可以計算得到模型在2190階的重力異常累計誤差約為5 mGal。本文使用5 mGal作為重力異常數據的單位權中誤差，並基於該單位中誤差以及上文中GOCE數據的單位權方差得到重力異常法方程與GOCE觀測法方程的加權因子。

　　聯合觀測法方程分為兩部分，其中251~2159階係數法方程與重力異常法方程中對應部分相同，而2~250階係數法方程是使用1.2節中的最小二乘聯合平差方法聯合101~250階重力異常數據塊對角法方程和GOCE全階次法方程(220階)得到的。構建重力異常數據101~250階塊對角法方程時需要在原始重力異常數據中移除2~100階及251~2159階係數的貢獻。

　　模型的220~250階係數是聯合GOCE數據和重力異常數據聯合求解的，這樣可以保證聯合模型的這部分係數及其驗後方差的頻譜有更好的連續性。重力異常數據求解低階係數的精度較差，重力異常低階次係數法方程與衛星法方程聯合可能會「汙染」衛星重力的低階係數，所以為減弱其影響，2~250階聯合法方程中只採用了重力異常的101~250階係數塊對角法方程。

　　4 SGG-UGM-1超高階重力場模型的精度分析4.1 SGG-UGM-1模型的係數誤差

　　EIGEN-6C4模型是當前精度較高的超高階重力場模型[11]。本文以此為參考模型，認為該模型的係數為「真值」計算SGG-UGM-1、GOSG-EGM、EGM2008和EIGEN-6C2[22]模型的「誤差」，並對它們的誤差進行對比分析，以在頻域內分析SGG-UGM-1的特性，評價模型的精度與合理性。

　　GOSG-EGM模型是利用本文計算得到的220階衛星重力場模型與EGM2008重力場模型係數直接拼接得到的，其2~220階係數來自衛星重力場模型，221~2159階係數來自EGM2008模型。圖 5給出了SGG-UGM-1、GOSG-EGM、EGM2008和EIGEN-6C2模型係數的階誤差RMS。

　　

圖 5 重力場模型的係數階誤差RMSFig. 5 Degree error RMS of the gravity field models

　　圖選項

　　

　　由圖 5可知，在2~70階的頻率範圍內，EIGEN-6C2的誤差最小，EGM2008次之，而SGG-UGM-1和GOSG-EGM較大。這是由於EIGEN-6C2與EIGEN-6C4的2~70階係數的計算方法和使用數據類似，都是聯合LAGEOS、GRACE和GOCE等數據求解的[9, 22]，EGM2008的2~70階次係數主要是來自於GRACE觀測數據的貢獻；而SGG-UGM-1和GOSG-EGM的2~70階次係數主要由GOCE數據求解，因此其與EIGEN-6C4模型的差異大於EGM2008、EIGEN-6C2。

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

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

　　

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

　　圖選項

　　

　　頻譜域的分析結果說明，本文求解的SGG-UGM-1模型精度可靠，且相比EGM2008模型在220階次內的係數精度有了明顯改善。

　　4.2 中國和美國GPS/水準數據的檢核結果

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

　　表 1 中國區域的GPS/水準數據檢核結果Tab. 1 Validation of the gravity field models using GPS-leveling data in China

　　

m模型最大值最小值平均值標準差EGM20081.729-1.5350.2390.240GOSG-EGM0.704-0.7930.2400.169SGG-UGM-10.744-0.6180.2460.162EIGEN-6C20.808-0.5920.2430.167EIGEN-6C40.729-0.6980.2430.157

　　表選項

　　表 2 美國區域的GPS/水準數據檢核結果Tab. 2 Validation of the gravity field models using GPS-leveling data in America

　　

m模型最大值最小值平均值標準差EGM20080.360-1.396-0.5110.284GOSG-EGM0.392-1.470-0.5110.291SGG-UGM-10.317-1.407-0.5110.280EIGEN-6C20.398-1.415-0.5120.283EIGEN-6C40.397-1.392-0.5120.282

　　表選項

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

　　4.3 毛烏素測區航空重力數據的檢核結果

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

　　表 3 航空重力數據的檢核結果Tab. 3 Validation of the gravity field models using airborne gravity data

　　

mGal模型最大值最小值平均值標準差EGM20084.694-25.565-10.2065.556GOSG-EGM3.089-28.550-9.8836.506SGG-UGM-15.046-26.715-10.0955.792EIGEN-6C27.779-23.414-7.5625.851EIGEN-6C44.179-25.227-10.1185.554

　　表選項

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

　　5 總結

　　本文使用EGM2008模型的5′×5′格網重力異常構建了2159階次的塊對角法方程，並使用最小二乘方法聯合GOCE重力衛星任務220階次的全階次法方程求解了2159階次的重力場模型SGG-UGM-1，使用EIGEN-6C4等重力場模型、中國與美國的GPS水準/數據和毛烏素測區的航空重力數據分析了SGG-UGM-1模型的內外符合精度。主要的結論如下：

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

　　(2) 採用最小二乘聯合平差方法不僅可以保證求解模型在聯合頻段內信號和誤差頻譜的連續性，同時求解的模型SGG-UGM-1精度優於直接拼接獲得的GOSG-EGM模型。

　　(3) 本文使用的GPS/水準數據對模型的檢核結果表明，SGG-UGM-1大地水準面在中國區域的精度介於EIGEN-6C2和EIGEN-6C4兩個模型之間，優於GOSG-EGM模型，遠優於EGM2008模型；在美國區域這些模型的精度相當。

　　(4) 本文使用的毛烏素測區航空重力數據對模型的檢核結果表明，SGG-UGM-1模型計算的重力擾動精度與EGM2008、EIGEN-6C4模型相當，優於GOSG-EGM模型和EIGEN-6C2模型。

　　模型的檢核結果說明本文構建超高階重力場模型策略的合理性，為進一步研究超高階重力場模型的構建奠定了基礎。但在聯合衛星觀測數據和重力異常數據時，本文參考了EIGEN-6C系列模型構建的經驗，後續應該進一步深入研究重力異常法方程與衛星觀測法方程的最優聯合問題。SGG-UGM-1模型的低階次係數的精度較EGM2008有了較大提升，但由於本文使用的重力異常數據是由EGM2008模型計算得到的，所以模型在高頻段沒有提升，且與EGM2008模型有很強的相關性。更為合理的做法是從全球高精度的地面重力、航空重力、衛星測高、地形等觀測數據出發，完全獨立自主處理各類數據得到全球高解析度的重力異常格網數據集，然後據此反演超高階重力場模型，這也是筆者後續的研究重點。

　　【引文格式】梁偉, 徐新禹, 李建成, 等. 聯合EGM2008模型重力異常和GOCE觀測數據構建超高階地球重力場模型SGG-UGM-1[J]. 測繪學報，2018，47(4)：425-434. DOI: 10.11947/j.AGCS.2018.20170269

　　

　　

　　往期 精彩回顧

　　「2017年度中國遙感領域十大事件」揭曉

　　同行評議的「醜陋面」！

　　漲知識！全球有多少CORS站？

　　2018年泰晤士世界大學200強榜單，中國入榜12所，快來看看都有誰？

　　楊元喜院士: 北鬥三號試驗系統性能分析【附PPT】

　　數學能力是中興與華為的唯一區別

　　關於在學術論文署名中常見問題或錯誤的誠信提醒

　　新時代之光 | 武大六院士當選「荊楚楷模」

　　足球大亨許家印1000億殺入信息產業，重點布局商業航天、智能交通、AI

　　寧津生院士：測繪界的十幾位院士正在研究這個問題

　　震撼！中國北鬥導航系統為何令歐美恨之入骨？

　　世界一流學科發布，哈佛再佔15個，中國高校9個學科世界第一

　　發表中文論文就很low？再談「中文期刊」與「文化自信」

　　李德仁院士：老師教我做人做學問

　　8個地球的科學冷知識顛覆你的世界觀！

　　關於稿件「時間」安排那些事兒~

　　世界上最有趣最冷門的地圖，刷新你的世界觀！

　　黃昕：當你動筆，成敗已定——來自IEEE評審專家的體會與思考

　　院士論壇| 高俊：圖到用時方恨少, 重繪河山待後生——《測繪學報》60年紀念與前瞻

　　

　　權威 | 專業 | 學術 | 前沿

　　微信投稿郵箱 | [e-mail]

　　微信公眾號中搜索「測繪學報」，關注我們，長按上圖二維碼，關注學術前沿動態。

　　歡迎加入《測繪學報》作者QQ群： 297834524

　　進群請備註：姓名+單位+稿件編號