趙保成，徐健，徐堅

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

0 引言

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

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

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

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

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

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

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

2 原理與方法

2.1 超高階地球重力場模型

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

EGM2008 地球重力場模型是由美國國家地理空間情報局通過先進建模技術(shù)與算法，以PGM2007B為參考，結(jié)合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)利用新技術(shù)改進的歐洲地球重力模型(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年采用球諧分析方法構(gòu)建并發(fā)布的一個新的2 190 階地球重力場模型，使用的數(shù)據(jù)包括衛(wèi)星重力觀測數(shù)據(jù)、衛(wèi)星測高數(shù)據(jù)和EGM2008 模型重力異常數(shù)據(jù).

4)XGM2019e_2159 地球重力場模型

該模型于2019年由慕尼黑技術(shù)大學天文和物理大地測量學研究所發(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 模型結(jié)合組成.

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

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

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

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

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

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

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

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

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為權(quán)矩陣，在實際應用中，常使用單位權(quán)矩陣代入：

3 算例驗證

3.1 實驗區(qū)概況

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

圖3 實驗點位分布圖

3.2 實驗方案設(shè)計及擬合精度評定

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

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

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

3.3 實驗對比分析

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

圖4 模型計算高程異常圖

由表1 所示，5 種地球重力場模型計算的高程異常值與真值殘差均較小，中誤差均不超過±0.012 0 m，殘差最大值為GECO 模型計算得出的結(jié)果，殘差值為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，在此表明了基于超高階地球重力場模型結(jié)合改進的“移去-擬合-恢復”法進行GNSS 高程轉(zhuǎn)換的正確性與優(yōu)越性.其中，利用XGM2019e_2159重力場模型的方案4 在各方案中的中誤差、極差、偏度值、峰度值均為最小，分別為±4.786 6 mm、18.875 7 mm、-0.648 8、0.887 8，從而證明了利用該模型計算結(jié)果相對準確和可靠.而模型EGM2008、EIGEN-6C4和SGG-UGM-2 擬合精度大體一致，模型GECO 擬合精度相對最差.

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

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

4 結(jié)束語

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

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

在絕大多數(shù)的工程建設(shè)中，高程控制測量往往采用水準法或三角高程法，雖然保證了精度，但是上述兩種方法的時間及人力成本較高.對于面積廣闊、處于偏遠地區(qū)、嚴重缺乏已有高程資料的測區(qū)的高程測量工作，本文的研究成果或?qū)⒕哂幸欢ǖ膮⒖贾笇б饬x.