基于超高階地球重力場模型的GNSS 高程轉化方法

趙保成，徐健，徐堅

(1.長江科學院 空間信息技術應用研究所,武漢 430010；2.武漢市智慧流域工程技術研究中心,武漢 430010)

0 引言

近年來，隨著衛(wèi)星導航定位技術和移動網(wǎng)絡通訊技術的快速發(fā)展，基于連續(xù)運行參考基站(CORS)的網(wǎng)絡實時動態(tài)(RTK)技術發(fā)展勢頭迅猛，千尋知寸FindCM、中國移動CORS、北斗CORS、厘清/Locate-CM 以及各省的CORS 服務等極大地豐富了測繪技術人員的技術手段.借助于CORS 服務，在移動網(wǎng)絡覆蓋區(qū)域，單人使用單臺全球衛(wèi)星導航系統(tǒng)(GNSS)接收機便能夠瞬時獲取在指定參考坐標系下的厘米級三維(3D)定位結果，相比于傳統(tǒng)GNSS 應用的電臺“1+N”模式，該服務使測繪作業(yè)效率得到了極大地提升.但不可忽略的是，無論采用何種CORS 服務，衛(wèi)星接收機直接獲取的高程信息始終是基于參考橢球面的大地高，而在我國的高程體系中使用的卻是基于似大地水準面的正常高，二者高程系統(tǒng)不一致導致GNSS 在高程測量中的應用受到了一定的限制.為了使GNSS 測得的大地高直接應用于工程建設中，必須要進行高程數(shù)據(jù)轉換[1-3].

GNSS 高程轉換方法主要有兩種.一種是GNSS水準高程法，其原理為利用均勻分布于測區(qū)的同名點的大地高和水準高(正常高)成果通過數(shù)學函數(shù)模型確定局部似大地水準面，由于一般測區(qū)內(nèi)很少布設有高精度的GNSS 水準點，如果從外部引測，時間及人工成本高，因此該方法的適用范圍有限；另一種方法是GNSS 重力高程法，其原理是利用地面或者衛(wèi)星重力資料求解待測點的高程異常值，結合GNSS 接收機獲取高精度大地高，進而求得待測點正常高[4].

隨著全球衛(wèi)星重力測量技術的不斷發(fā)展，利用地面重力數(shù)據(jù)結合衛(wèi)星重力數(shù)據(jù)獲取的超高階地球重力場模型的精度和分辨率在不斷提高.基于此，眾多行業(yè)學者嘗試利用超高階地球重力場模型進行了諸多系統(tǒng)GNSS 高程轉換方法研究.文獻[5]利用EGM2008 地球重力場模型去除高程異常中的長波項，進而使用加權組合模型對高程異常的剩余項擬合，擬合效率得到了一定的提升；文獻[6]基于二次曲面和EGM2008 地球重力場模型，構建了多種移去-恢復型函數(shù)擬合模型，探究了各種模型的適用性；文獻[7]利用礦區(qū)水準觀測數(shù)據(jù)和EGM2008 地球重力場模型，綜合物理和幾何方法建立了礦區(qū)局部似大地水準面模型；文獻[8]利用EIGEN-6C4 地球重力場模型進行GNSS 高程擬合，擬合精度得到了進一步提升；文獻[9]詳細闡述了利用地球重力場模型位系數(shù)計算高程異常的數(shù)學模型，同時實現(xiàn)了全球高程異常的可視化；文獻[10]利用EIGEN-6C4、EIGEN-6C2、EGM2008 地球重力場模型結合移去-恢復法進行了高程轉換方法研究，得出了在我國顧及EIGEN-6C4地球重力場模型的GNSS 高程轉換精度更高的結論；文獻[11]通過GNSS 高程轉換精度分析，得出了XGM2019 重力場模型與我國似大地水準面符合程度更高的結論；文獻[12]討論了在不同地形條件下利用EGM2008 地球重力場模型轉換大地高為正常高的方法；文獻[13]利用常規(guī)水準限差評定了基于EGM2008 地球重力場模型的GNSS 擬合高程等級及適用性.

綜上，目前GNSS 高程轉換方法的研究存在兩方面問題，一是研究雖引入了高階地球重力場模型，但多是針對單種地球重力場模型并結合不同插值擬合方法，鮮有對利用目前已存的各種超高階地球重力場模型進行高程轉換方法效果的綜合評估；二是插值擬合方法多是從擬合函數(shù)的數(shù)學模型角度出發(fā)，而對于參與插值擬合的特征點的選取方法及原則并未做深入研究.本文在總結了多位行業(yè)學者研究成果的基礎上，針對截至2022年已存的5 種超高階地球重力場模型，結合改進的“移去-擬合-恢復”法的GNSS 高程轉換方法開展了研究，本研究針對在缺乏高程基礎控制資料及地面重力資料測區(qū)的高程系統(tǒng)的建立工作上提供解決方案及思路.

1 高程系統(tǒng)簡介

目前，常用的高程系統(tǒng)有三種，分別是以參考橢球面為基準面的大地高系統(tǒng)、以大地水準面為基準面的正高系統(tǒng)及以似大地水準面為基準面的正常高系統(tǒng).大地高H為地面點沿法線至參考橢球面的距離；正高h為地面點沿鉛垂線至大地水準面的距離；正常高hr為地面點沿鉛垂線至似大地水準面的距離[14].在不考慮垂線偏差影響情況下，三種高程系統(tǒng)之間的幾何關系如圖1 所示，數(shù)學關系由式(1)～(2)所示，式中ζ 表示高程異常，N表示大地水準面差距.我國法定使用的1985 國家高程基準的高程系統(tǒng)即是正常高系統(tǒng).

圖1 高程系統(tǒng)關系圖

2 原理與方法

2.1 超高階地球重力場模型

1)地球重力場模型(EGM2008)

EGM2008 地球重力場模型是由美國國家地理空間情報局通過先進建模技術與算法，以PGM2007B為參考，結合Topex 衛(wèi)星測高數(shù)據(jù)、GRACE 衛(wèi)星重力數(shù)據(jù)、地面重力數(shù)據(jù)、地形數(shù)據(jù)，得到拓展至2 190階球諧系數(shù)的全球超高階地球重力場模型.EGM2008地球重力場模型的地球覆蓋率高達83.80%，模型的基本分辨率為5'×5'，最高分辨率為1'×1'.

2)利用新技術改進的歐洲地球重力模型(EIGEN-6C4)

EIGEN-6C4 地球重力場模型于2014年發(fā)布，其球諧系數(shù)的階數(shù)和次數(shù)均為2 190.數(shù)據(jù)源主要為1985—2010年的LAGEO 數(shù)據(jù)、2003—2012年的GRACERL03GRGS數(shù)據(jù)、完整的GOCE-SGG 數(shù)據(jù)、DTU10 地面重力數(shù)據(jù).

3)SGG-UGM-2 地球重力場模型

SGG-UGM-2 是由我國梁偉等學者于2020年采用球諧分析方法構建并發(fā)布的一個新的2 190 階地球重力場模型，使用的數(shù)據(jù)包括衛(wèi)星重力觀測數(shù)據(jù)、衛(wèi)星測高數(shù)據(jù)和EGM2008 模型重力異常數(shù)據(jù).

4)XGM2019e_2159 地球重力場模型

該模型于2019年由慕尼黑技術大學天文和物理大地測量學研究所發(fā)布，其階次完全至2 159(球諧系數(shù)的階擴展至2 190)，空間分辨率5'×5'，數(shù)據(jù)來源包括GOCO 06 衛(wèi)星重力模型、陸地和海洋重力異常和地面上的地形導出重力.

5)GECO 地球重力場模型

該模型于2015年發(fā)布，其球諧系數(shù)的階為2 190，主要采用簡單譜組合原理將GO_CONS_GCF_2_TIM_R5 模型和EGM2008 模型結合組成.

2.2 改進的移去-擬合-恢復法

根據(jù)物理大地測量學理論[15]，地球任意點高程異常通常由2 部分構成，分別為中長波項與短波項，如式(3)所示，高程異常的中長波分量變化較小，趨勢相對平緩，短波分量主要由地形起伏引起，變化相對較大[16].

“移去-擬合-恢復法”是GNSS 高程轉換的常規(guī)方法，其基本思想為：首先將參與擬合點的高程異常中變化較平穩(wěn)的中長波分量移除，只對變化較大的短波項進行擬合，然后利用空間插值計算出待求點上的短波項，最后在待定點上恢復中長波分量.

然而，在擬合的步驟中，參與擬合點的選取往往具有一定的盲目性，極易造成擬合曲面關鍵信息的損失，最終導致高程轉換精度降低.為了解決上述問題，本文提出了一種改進的“移去-擬合-恢復”法，其具體操作流程如圖2 所示.

圖2 改進的移去-擬合-恢復法流程圖

1)確定測區(qū)范圍，將測區(qū)網(wǎng)格化處理，格網(wǎng)尺寸根據(jù)項目成果要求的比例尺精度需求確定，利用全球重力場模型計算測區(qū)格網(wǎng)點的高程異常值，并將其作為真實高程異常值的中長波分量，進而構建基于高程異常中長波分量的似大地水準面模型并將其柵格化顯示.

2)利用ArcMap 軟件中的空間分析模塊對上一步得到的似大地水準面曲面模型(基于高程異常中長波分量)進行地形表面及3D 分析，計算其高程異常曲面的極大值、極小值、坡向、坡度、曲率、粗糙度等參數(shù)，進而選擇特征點作為GNSS 水準聯(lián)測點參與后續(xù)的曲面擬合計算.對于特征點選取有以下原則：所選點位在似大地水準面上趨勢變化較大；點位的個數(shù)需保證曲面擬合參數(shù)方程有解；點位的位置需覆蓋整個測區(qū).

3)利用GNSS 水準聯(lián)測法獲取特征點的高精度大地高、正常高，求得真實高程異常，進而移去擬合點上的高程異常中長波分量，將參與擬合點的平面坐標作為輸入值，高程異常短波項作為輸出值，構建二次曲面擬合模型.

4)利用地球重力場模型計算出待求點的高程異常中長波分量，再加上擬合步驟中該點高程異常短波項擬合值，最終得到待求點的高程異常值.

2.3 利用地球重力場模型求解高程異常中長波項

可以利用2.1 節(jié)介紹的5 種超高階地球重力場模型獲取地表任意點重力勢的泛函值，例如重力異常、高程異常、垂線偏差等[15].利用Bruns 公式求解地球上任意點的高程異常中長波項ζGM，公式為

式中：ρ為計算點地心向徑(計算點到地心的距離)；Ψ為余緯(緯度90°)；λ為大地經(jīng)度；GM為包含大氣層在內(nèi)的地心引力常數(shù)；a為地球橢球長半徑；、nm為完全規(guī)格化位系數(shù)；(sin Ψ)為完全規(guī)格化締合函數(shù)；γ為計算點的正常重力位.

2.4 二次多項式曲面擬合

高程異常短波項的曲面擬合方法較多，常用的方法有二次曲面擬合、多面函數(shù)法、曲面樣條函數(shù)法、移動曲面法等.在實際工程應用中，為方便計算，大多采用二次曲面函數(shù)擬合高程異常短波項，其數(shù)學表達式為

式中：ζ短波為特征點的高程異常短波項；(xi,yi)為特征點在某參考系統(tǒng)內(nèi)的平面直角坐標；(a0,a1,a2,a3,a4,a5)為模型待定系數(shù).

根據(jù)曲面擬合函數(shù)模型可知，共有6 個未知參數(shù)，因此為求解模型系數(shù)，至少需要列立6 個方程.當已知點個數(shù)大于6 時，表達式為

將上式變換為誤差方程形式如式(7)所示，其中ε為模型的擬合殘差：

根據(jù)最小二乘法原理，令模型擬合殘差平方和最小(i≥6)，在此條件下，可以求得二次曲面模型待定系數(shù)X，如式(13)所示，其中，P為權矩陣，在實際應用中，常使用單位權矩陣代入：

3 算例驗證

3.1 實驗區(qū)概況

實驗區(qū)位于湖北省某縣城城區(qū)，此區(qū)域內(nèi)相對高差較小，區(qū)內(nèi)均勻地設置有20 個等精度GNSS 水準聯(lián)測點，點位平面按照《全球定位系統(tǒng)(GPS)測量規(guī)范》中D 級要求施測，高程值采用《國家三、四等水準測量規(guī)范》進行三等水準聯(lián)測.具體的點位分布如圖3 所示.

圖3 實驗點位分布圖

3.2 實驗方案設計及擬合精度評定

按照2.2 節(jié)中介紹的擬合點選點方法選取圖3中用紅色三角形表示的6 個點作為擬合點，選取黑色圓形表示的14 個點作為檢核點，開展高程異常擬合實驗.實驗方法采用改進的“移去-擬合-恢復”法，分別采用2.1 節(jié)介紹的5 種超高階地球重力場模型進行高程異常中長波項的求解，擬合方法均采用二次曲面擬合法.

擬合效果一般通過對比擬合值與真值的符合程度來評估，主要指標有中誤差、極差、偏度值、峰度值等.中誤差對一組測量值中的異常值非常敏感，能較好地反映測量結果波動大小，因此本文主要使用中誤差評定高程異常擬合精度，極差、偏度值、峰度值作為輔助參考值.中誤差計算式為

式中：mζ為高程異常擬合中誤差；Δ為在檢核點上曲面擬合模型插值計算的高程異常值與GNSS 水準計算得出的高程異常值(真值)的差值；n為檢核點的個數(shù).

3.3 實驗對比分析

如圖4 所示，本實驗區(qū)內(nèi)使用5 種地球重力場模型求解的高程異常均為負值，基本與真實高程異常值的變化趨勢一致，其中SGG-UGM-2 模型計算得到的高程異常與真實高程異常變化趨勢符合程度更高.

圖4 模型計算高程異常圖

由表1 所示，5 種地球重力場模型計算的高程異常值與真值殘差均較小，中誤差均不超過±0.012 0 m，殘差最大值為GECO 模型計算得出的結果，殘差值為0.094 3 m.其中，SGG-UGM-2 地球重力場模型計算得出的高程異常值與真值符合程度更高，中誤差為±0.009 3 m.

表1 模型高程異常精度評定表 m

由圖5 及表2 可知，5 種方案的擬合中誤差分別為±5.785 4 mm、±5.810 0 mm、±5.815 3 mm、±4.786 6 mm、±7.622 5 mm，擬合效果均比較優(yōu)秀，中誤差均不超過±8.000 0 mm，在此表明了基于超高階地球重力場模型結合改進的“移去-擬合-恢復”法進行GNSS 高程轉換的正確性與優(yōu)越性.其中，利用XGM2019e_2159重力場模型的方案4 在各方案中的中誤差、極差、偏度值、峰度值均為最小，分別為±4.786 6 mm、18.875 7 mm、-0.648 8、0.887 8，從而證明了利用該模型計算結果相對準確和可靠.而模型EGM2008、EIGEN-6C4和SGG-UGM-2 擬合精度大體一致，模型GECO 擬合精度相對最差.

圖5 模型計算高程異常擬合殘差圖

表2 5 種方案的擬合精度評定表

4 結束語

基于5 種不同的超高階地球重力場模型，本文在湖北省某縣城城區(qū)開展了GNSS 高程轉換方法實驗研究，通過對實驗結果對比分析，表明了基于超高階地球重力場模型結合改進的“移去-擬合-恢復”法的GNSS 高程轉換方法的可行性與可靠性.在不考慮我國高程基準與全球高程基準之間高差的影響下，在實驗測區(qū)內(nèi)，SGG-UGM-2 地球重力場模型計算得出的高程異常值與真值符合程度在5 種模型中最高.在改進的“移去-擬合-恢復”法實驗中，基于XGM2019e_2159地球重力場模型的高程轉換效果最好，高程異常擬合中誤差、極差、偏度值和峰度值最小.

本文提出了基于地球重力場模型改進的 “移去-擬合-恢復”方法，相比于僅依賴數(shù)學模型的擬合方法更具有實際的物理意義.除此之外，該方法還在曲面擬合特征點的選擇上更具有科學性和合理性，擬合殘差更小.

在絕大多數(shù)的工程建設中，高程控制測量往往采用水準法或三角高程法，雖然保證了精度，但是上述兩種方法的時間及人力成本較高.對于面積廣闊、處于偏遠地區(qū)、嚴重缺乏已有高程資料的測區(qū)的高程測量工作，本文的研究成果或將具有一定的參考指導意義.