顧及似大地水準(zhǔn)面趨勢(shì)變化的高程異常擬合方法

趙保成 李國(guó)忠 徐 堅(jiān) 徐 健 肖 瀟

(1. 長(zhǎng)江科學(xué)院 空間信息技術(shù)應(yīng)用研究所, 湖北 武漢 430010; 2. 武漢市智慧流域工程技術(shù)研究中心, 湖北 武漢 430010)

0 引言

2020年北斗3號(hào)全球組網(wǎng)成功,國(guó)際上基本形成了以美國(guó)全球定位系統(tǒng)(global positioning system,GPS)、俄羅斯格洛納斯、歐盟伽利略、中國(guó)北斗為主的四大全球衛(wèi)星導(dǎo)航定位系統(tǒng)。伴隨著各星座系統(tǒng)導(dǎo)航衛(wèi)星數(shù)目的增加以及迭代更新,全球?qū)Ш蕉ㄎ幌到y(tǒng)(global navigation satellite system,GNSS)在各行各業(yè)的應(yīng)用也將越來(lái)越廣泛,尤其對(duì)于測(cè)繪地理信息行業(yè)。GNSS測(cè)量具備全天時(shí)、全天候、累積誤差小、測(cè)站不受通視條件影響的優(yōu)點(diǎn),它的出現(xiàn)極大提升了測(cè)繪作業(yè)效率[1-3]。但是GNSS接收機(jī)直接獲取的測(cè)點(diǎn)高程信息是基于參考橢球面的大地高,而在實(shí)際工程應(yīng)用中,我國(guó)使用的卻是基于似大地水準(zhǔn)面的正常高,在不考慮垂線偏差影響的情況下,二者之間存在一個(gè)非固定差值即高程異常,因此,如何求取測(cè)點(diǎn)高程異常值成為GNSS高程數(shù)據(jù)成果直接應(yīng)用于實(shí)際工程中需首要解決的關(guān)鍵問(wèn)題。

針對(duì)上述關(guān)鍵問(wèn)題,眾多行業(yè)學(xué)者做出了諸多的嘗試。文獻(xiàn)[4]利用基于穩(wěn)健估計(jì)的總體最小二乘算法較好地解決了高程異常曲面擬合過(guò)程中已知點(diǎn)平面坐標(biāo)與高程異常值中均含有誤差的問(wèn)題;文獻(xiàn)[5]提出了一種基于多面函數(shù)和移動(dòng)曲面法的綜合高程異常擬合模型,在一定程度上提高了高程異常擬合精度;文獻(xiàn)[6]將多種常見(jiàn)高程異常曲面擬合模型進(jìn)行加權(quán)組合,對(duì)不同地形的適應(yīng)性進(jìn)行了分析;文獻(xiàn)[7]利用全球超高階地球重力場(chǎng)模型(Earth gravitationational model 2008,EGM2008)去除高程異常中的長(zhǎng)波項(xiàng),進(jìn)而使用加權(quán)組合模型對(duì)高程異常的剩余項(xiàng)擬合,擬合效率得到了一定的提升;文獻(xiàn)[8]提出了蟻群—遺傳算法結(jié)合多面函數(shù)擬合的改進(jìn)曲面擬合方法,在很大程度上提高了擬合模型的精度;文獻(xiàn)[9]基于二次曲面和EGM2008全球重力場(chǎng)模型,構(gòu)建了多種移去—恢復(fù)型函數(shù)擬合模型,探究了各種模型的適用性;文獻(xiàn)[10]利用礦區(qū)水準(zhǔn)觀測(cè)數(shù)據(jù)及重力場(chǎng)模型,綜合物理和幾何方法建立了礦區(qū)似大地水準(zhǔn)面模型;文獻(xiàn)[11]構(gòu)建了二次曲面—徑向基函數(shù)(radial basis function,RBF)神經(jīng)網(wǎng)絡(luò)組合的GPS高程轉(zhuǎn)換模型,提高了單一模型的擬合精度。

綜上所述,現(xiàn)有的高程異常擬合研究多是從擬合函數(shù)模型角度自身出發(fā),或是模型的改進(jìn),或是多種模型的加權(quán)組合,缺少了對(duì)參與高程異常曲面擬合的特征點(diǎn)的選取方法的研究,特征點(diǎn)若選取不當(dāng),將導(dǎo)致擬合曲面關(guān)鍵信息丟失,影響最終的擬合效果。

本文在總結(jié)了多位行業(yè)學(xué)者研究成果的基礎(chǔ)上,提出了一種顧及似大地水準(zhǔn)面趨勢(shì)變化的高程異常擬合方法,首先利用SGG-UGM-2全球重力場(chǎng)模型預(yù)先獲取測(cè)區(qū)似大地水準(zhǔn)面趨勢(shì)面特征點(diǎn),后對(duì)特征點(diǎn)進(jìn)行GNSS水準(zhǔn)聯(lián)測(cè),最后結(jié)合移去恢復(fù)法及二次多項(xiàng)式曲面進(jìn)行高程異常擬合。經(jīng)過(guò)實(shí)例驗(yàn)證,本文提出的高程異常擬合方法,參與高程異常曲面擬合的特征點(diǎn)選取更加合理,關(guān)鍵信息損失少,擬合精度能夠達(dá)到厘米級(jí),完全滿足絕大多數(shù)大比例尺地形圖測(cè)量要求。

1 高程系統(tǒng)

我國(guó)常用的高程系統(tǒng)有3種,大地高、正高和正常高。大地高H為測(cè)點(diǎn)沿法線至參考橢球面的距離,可由GNSS接收機(jī)測(cè)得;正高h(yuǎn)為測(cè)點(diǎn)沿鉛垂線至大地水準(zhǔn)面的距離,此值為理論值,一般不可直接測(cè)得;正常高h(yuǎn)r為測(cè)點(diǎn)沿鉛垂線至似大地水準(zhǔn)面的距離,一般通過(guò)水準(zhǔn)測(cè)量方法確定[12]。三者的幾何關(guān)系如圖1所示,數(shù)學(xué)關(guān)系見(jiàn)式(1)與式(2)。我國(guó)法定使用的1985國(guó)家高程基準(zhǔn)使用的高程系統(tǒng)即是正常高系統(tǒng)。式(1)中ζ表示高程異常,式(2)中N表示大地水準(zhǔn)面差距。

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

2 原理與方法

2.1 利用SGG-UGM-2全球重力場(chǎng)模型求解高程異常中長(zhǎng)波項(xiàng)

SGG-UGM-2全球重力場(chǎng)模型是由我國(guó)梁偉等學(xué)者于2020年采用球諧分析方法構(gòu)建并發(fā)布的一個(gè)最新的2 190階2 159次全球重力場(chǎng)模型,模型分辨率為5′,構(gòu)建此模型使用的數(shù)據(jù)包括全球衛(wèi)星重力觀測(cè)數(shù)據(jù)、衛(wèi)星測(cè)高數(shù)據(jù)和EGM2008全球重力場(chǎng)模型的重力異常數(shù)據(jù)[13],在我國(guó),此模型精度優(yōu)于其他全球重力場(chǎng)模型(如SGG-UGM-1、EGM2008等)。以SGG-UGM-2全球重力場(chǎng)模型為基礎(chǔ),利用Bruns公式[14]求解地球上任意點(diǎn)的高程異常中長(zhǎng)波項(xiàng)ζGM,公式如下:

(3)

2.2 二次多項(xiàng)式曲面擬合

曲面擬合常用的方法有二次多項(xiàng)式曲面擬合、多面函數(shù)法、曲面樣條函數(shù)法、移動(dòng)曲面法等。在實(shí)際工程應(yīng)用中,為方便計(jì)算,大多采用二次多項(xiàng)式曲面函數(shù)擬合高程異常短波項(xiàng),其數(shù)學(xué)表達(dá)如式(4)所示。

(4)

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

根據(jù)列立方程組可知,共有6個(gè)未知數(shù),因此為求解模型系數(shù),至少需要列立6個(gè)方程。當(dāng)已知點(diǎn)個(gè)數(shù)大于6時(shí),其表達(dá)如式(5)所示。

(5)

將式(5)變換為誤差方程形式，如式(6)所示,其中，ε為模型的擬合殘差。

(6)

根據(jù)最小二乘法原理,令模型擬合殘差平方和最小(i≥6),在此條件下,可以求得二次多項(xiàng)式曲面模型待定系數(shù)X,見(jiàn)式(7),其中式(8)中的P為權(quán)矩陣。

2.3 方法總體流程

根據(jù)物理大地測(cè)量學(xué)理論,任意測(cè)點(diǎn)高程異常可由為中長(zhǎng)波項(xiàng)與短波項(xiàng)構(gòu)成,中長(zhǎng)波分量變化較小,趨勢(shì)相對(duì)平緩,短波分量主要由地形起伏引起,變化相對(duì)較大[15-16],其數(shù)學(xué)關(guān)系如式(9)所示。

(9)

移去恢復(fù)法是高程異常擬合的經(jīng)典方法,該方法能夠有效地去除無(wú)關(guān)變量對(duì)擬合最終結(jié)果的影響,移去恢復(fù)法思想為：首先將高程異常中變化較小的中長(zhǎng)波分量移除,然后對(duì)短波項(xiàng)進(jìn)行曲面擬合并通過(guò)插值計(jì)算出待求點(diǎn)上的短波項(xiàng),進(jìn)而在待定點(diǎn)上恢復(fù)中長(zhǎng)波分量,最終獲取待求點(diǎn)上的高程異常。本文提出的顧及似大地水準(zhǔn)面趨勢(shì)變化的高程異常擬合方法對(duì)經(jīng)典的移去恢復(fù)法做了相應(yīng)的改進(jìn),其具體操作流程如圖2所示。

圖2 移去恢復(fù)流程圖

第一步：確定測(cè)區(qū)范圍并將其格網(wǎng)化,采用SGG-UGM-2全球重力場(chǎng)模型計(jì)算測(cè)區(qū)格網(wǎng)點(diǎn)的高程異常值(真實(shí)高程異常中的中長(zhǎng)波項(xiàng)),構(gòu)建基于高程異常中長(zhǎng)波項(xiàng)的似大地水準(zhǔn)面模型并將其圖形化顯示。

第二步：利用ArcMap軟件中的空間分析模塊對(duì)上一步得到的似大地水準(zhǔn)面曲面模型進(jìn)行地形表面分析,計(jì)算其高程異常極大值、極小值、坡向、坡度、曲率等參數(shù),合理選擇似大地水準(zhǔn)面趨勢(shì)變化較大的特征點(diǎn)作為GNSS水準(zhǔn)聯(lián)測(cè)點(diǎn),通過(guò)GNSS水準(zhǔn)方法求出特征點(diǎn)的真實(shí)高程異常值。

第三步：特征點(diǎn)的真實(shí)高程異常值減去第一步中計(jì)算得到的中長(zhǎng)波分量,剩余高程異常短波項(xiàng)。以特征點(diǎn)的平面坐標(biāo)作為輸入值,高程異常短波項(xiàng)作為輸出值,構(gòu)建二次多項(xiàng)式曲面擬合模型,根據(jù)模型曲面方程求解特征點(diǎn)覆蓋范圍內(nèi)的任意待求點(diǎn)的高程異常短波項(xiàng)。

第四步：利用SGG-UGM-2全球重力場(chǎng)模型計(jì)算出待求點(diǎn)的中長(zhǎng)波分量,再加上擬合步驟中該點(diǎn)高程異常短波項(xiàng),最終得到待求點(diǎn)的高程異常值。

3 實(shí)驗(yàn)驗(yàn)證

3.1 實(shí)驗(yàn)區(qū)概況

實(shí)驗(yàn)區(qū)位于長(zhǎng)江中游城市—武漢市,該區(qū)域?qū)儆诮瓭h平原,相對(duì)高差較小。實(shí)驗(yàn)區(qū)沿武漢市某河流岸線布置,長(zhǎng)度約70 km,寬度約0.3 km,實(shí)驗(yàn)區(qū)具有狹長(zhǎng)條形特點(diǎn)。實(shí)驗(yàn)區(qū)內(nèi)均勻地設(shè)置有若干等精度GNSS水準(zhǔn)聯(lián)測(cè)點(diǎn),點(diǎn)位平面按照《全球定位系統(tǒng)(GPS)測(cè)量規(guī)范》中D級(jí)點(diǎn)要求施測(cè),點(diǎn)位高程按照《國(guó)家三、四等水準(zhǔn)測(cè)量規(guī)范》中的三等點(diǎn)要求聯(lián)測(cè)。

3.2 實(shí)驗(yàn)區(qū)似大地水準(zhǔn)面趨勢(shì)分析

由于實(shí)驗(yàn)區(qū)呈線狀分布,因此選擇將實(shí)驗(yàn)區(qū)內(nèi)GNSS水準(zhǔn)聯(lián)測(cè)點(diǎn)自河流上游至下游串聯(lián)起來(lái),利用式(1)計(jì)算GNSS水準(zhǔn)聯(lián)測(cè)點(diǎn)的高程異常值真值,同時(shí)利用SGG-UGM-2全球重力場(chǎng)模型計(jì)算GNSS水準(zhǔn)聯(lián)測(cè)點(diǎn)高程異常中長(zhǎng)波分量,兩者對(duì)比結(jié)果如圖3所示。 GNSS水準(zhǔn)聯(lián)測(cè)點(diǎn)的高程異常值中的中長(zhǎng)波分量與高程異常的真值變化趨勢(shì)基本一致,基本反映了實(shí)驗(yàn)區(qū)內(nèi)似大地水準(zhǔn)面趨勢(shì)變化情況;實(shí)驗(yàn)區(qū)內(nèi)高程異常值均為負(fù)值,說(shuō)明實(shí)驗(yàn)區(qū)似大地水準(zhǔn)面整體位于橢球面(2000國(guó)家大地坐標(biāo)系橢球)之下,高程異常的最小值為-15.821 m,最大值為-14.254 m,相差1.567 m,沿河流上游至下游,自西向東,高程異常值的絕對(duì)值逐漸變小;實(shí)驗(yàn)區(qū)河流中上游段高程異常值變化近似呈線性規(guī)律,下游段高程異常值變化劇烈,整體呈現(xiàn)非線性變化特征。

圖3 實(shí)驗(yàn)區(qū)似大地水準(zhǔn)面趨勢(shì)

對(duì)于本實(shí)驗(yàn)區(qū)的高程異常擬合特征點(diǎn)選取工作：①需保證GNSS水準(zhǔn)聯(lián)測(cè)點(diǎn)覆蓋整個(gè)測(cè)區(qū)；②設(shè)置GNSS水準(zhǔn)聯(lián)測(cè)點(diǎn)的點(diǎn)數(shù)需確保高程異常擬合曲面方程有解;③根據(jù)似大地水準(zhǔn)面的地形分析結(jié)果,GNSS水準(zhǔn)聯(lián)測(cè)點(diǎn)需在似大地水準(zhǔn)面變化大的區(qū)域(實(shí)驗(yàn)區(qū)下游區(qū))盡量多選點(diǎn)。

3.3 實(shí)驗(yàn)方案設(shè)計(jì)

為驗(yàn)證本文提出的高程異常擬合方法的優(yōu)越性,采用二次多項(xiàng)式曲面擬合法、移去恢復(fù)型二次多項(xiàng)式曲面擬合法、顧及似大地水準(zhǔn)面趨勢(shì)變化的二次多項(xiàng)式曲面擬合法、顧及似大地水準(zhǔn)面趨勢(shì)變化的移去恢復(fù)型二次多項(xiàng)式曲面擬合法4種擬合方案進(jìn)行高程異常擬合效果對(duì)比實(shí)驗(yàn),將4種實(shí)驗(yàn)方案分別編號(hào)為方案1#、方案2#、方案3#、方案4#,如表1所示。

表1 擬合實(shí)驗(yàn)方案表

為減少無(wú)關(guān)變量對(duì)實(shí)驗(yàn)對(duì)比結(jié)果的影響,4種方案選擇GNSS水準(zhǔn)聯(lián)測(cè)點(diǎn)數(shù)量均為8個(gè)(最少點(diǎn)數(shù)為6),檢核點(diǎn)均為6個(gè),特征點(diǎn)采用黑色圓形表示,檢核點(diǎn)采用三角形表示,各方案的點(diǎn)位選擇空間分布如圖4所示。

(a)方案1#選點(diǎn)方案

(b)方案2#選點(diǎn)方案

(c)方案3#選點(diǎn)方案

(d)方案4#選點(diǎn)方案圖4 各實(shí)驗(yàn)選點(diǎn)方案

3.4 精度評(píng)定

中誤差對(duì)一組測(cè)量值中的異常值非常敏感,能較好地反映測(cè)量結(jié)果波動(dòng)大小,因此本次實(shí)驗(yàn)使用中誤差評(píng)定高程異常擬合精度,計(jì)算中誤差的表達(dá)見(jiàn)式(10)。

(10)

其中,mζ為高程異常擬合中誤差;Δ為在檢核點(diǎn)上曲面擬合模型插值計(jì)算的高程異常值與GNSS水準(zhǔn)計(jì)算得出的高程異常值(真值)的差值;n為檢核點(diǎn)的個(gè)數(shù)。

3.5 實(shí)驗(yàn)對(duì)比分析

對(duì)比4種方案的擬合效果,由圖5和表2的結(jié)果可得：

圖5 擬合殘差圖

表2 精度評(píng)定表 單位：m

(1)不考慮似大地水準(zhǔn)面趨勢(shì)變化選擇擬合特征點(diǎn)且不使用移動(dòng)恢復(fù)法的方案1#的高程異常擬合效果最差,最大殘差值為0.033 m,最小殘差值為0.012 m,中誤差為±0.025 m;選擇似大地水準(zhǔn)面趨勢(shì)面特征點(diǎn)參與曲面擬合且使用移去恢復(fù)法的方案4#擬合效果最好,最大殘差值為0.027 m,最小殘差值為0.000 m,中誤差為±0.014 m。

(2)對(duì)比方案1#與方案3#,二者均未采用移去恢復(fù)法,擬合中誤差均為±0.025 m,總體擬合效果相當(dāng);方案3#考慮了似大地水準(zhǔn)面趨勢(shì)變化,相比方案1#的擬合殘差的極大值和極小值都小,在未采用移去恢復(fù)法的情況下,考慮似大地水準(zhǔn)面趨勢(shì)變化有助于提高高程異常擬合精度,但是總體擬合效果并不明顯。

(3)對(duì)比方案2#與方案4#,二者均采用了移去恢復(fù)法,方案4#考慮了似大地水準(zhǔn)面趨勢(shì)變化,其擬合中誤差為±0.014 m,方案2#擬合中誤差為±0.017 m,方案4#比方案2#總體擬合精度提升了17.65%。因此,在采用移去恢復(fù)法情況下,考慮似大地水準(zhǔn)面趨勢(shì)變化的選點(diǎn)方案能夠顯著地提升高程異常擬合精度。

(4)由方案1#與方案2#、方案3#與方案4#實(shí)驗(yàn)結(jié)果對(duì)比表明：移去恢復(fù)法能夠有效地提升高程異常曲面擬合精度,其中方案2#比方案1#擬合精度提升32.00%,方案4#比方案3#擬合精度提升44.00%。

4 結(jié)束語(yǔ)

基于SGG-UGM-2全球重力場(chǎng)模型,結(jié)合移去恢復(fù)擬合法,本文提出了一種顧及似大地水準(zhǔn)面趨勢(shì)變化的高程異常擬合方法。在對(duì)陌生測(cè)區(qū)無(wú)任何高程先驗(yàn)知識(shí)的情況下,此方法能夠相對(duì)準(zhǔn)確地判斷出測(cè)區(qū)高程異常曲面擬合中關(guān)鍵已知點(diǎn)信息,從而減小曲面擬合過(guò)程中的精度損失。通過(guò)實(shí)例驗(yàn)證,本文提出的高程異常曲面擬合方法最大殘差值為0.027 m,最小殘差值為0.000 m,中誤差為±0.014 m,其完全能夠使實(shí)驗(yàn)區(qū)的GNSS高程測(cè)量替代四等水準(zhǔn)方法。