基于移動二次曲面擬合的GNSS點高程異常計算

吳新輝

(漳州市測繪設(shè)計研究院，福建 漳州 363000)

包括北斗衛(wèi)星導航系統(tǒng)在內(nèi)的全球?qū)Ш叫l(wèi)星系統(tǒng)(GNSS)及連續(xù)運行參考站(CORS)已逐步成為測繪等相關(guān)領(lǐng)域所采用的主要測量手段[1]。目前我國使用的高程是正常高，GNSS獲得的高程為大地高，大地高數(shù)據(jù)對于高程控制測量只有幾何意義，只有獲得大地水準面精化模型成果，將測量所得的大地高轉(zhuǎn)換為正常高，才有實際應用意義[2]。

GNSS直接獲取的高程為大地高H，而我國采用正常高h，所以大地高數(shù)據(jù)對于高程控制測量僅具有幾何意義。因此，在實際生產(chǎn)應用過程中，應采用正常高才能起到高程控制作用，所以如何將獲取的大地高快速、準確地轉(zhuǎn)換為正常高已成為測繪工作中的一項重要內(nèi)容。

當前測繪生產(chǎn)中，一般采用靜態(tài)測量或者精密單點定位技術(shù)(PPP)獲得一些數(shù)量且分布合理的點的大地高，同時采用傳統(tǒng)水準測量方式進行高程聯(lián)測，然后擬合出該“空白區(qū)”的似大地水準面，即高程異常值，最后得到該區(qū)域其他待求點的正常高[3]。該方法被廣泛應用于地形平坦地區(qū)，且精度可達到四等水準或更高[4]，但效率較低，人力、物力投入較多，生產(chǎn)周期長，對于地形起伏較大、時間緊迫的任務不適用。隨著重力場模型的精度及分辨率越來越高，采用少量GNSS水準點進行高程異常擬合快速獲取正常高成為可能。

EGM2008重力場模型在中國大陸區(qū)域具有一定的適應性，且精度優(yōu)于其他模型[5]。本研究在EGM2008重力場模型的基礎(chǔ)上，采用漳州市區(qū)實測C級控制網(wǎng)數(shù)據(jù)，比較了移動曲面擬合模型中平面法、二次曲面法、加權(quán)平均法這3種方法在殘余高程異常計算中的應用，結(jié)果表明在EGM2008重力場模型的基礎(chǔ)上，對殘余高程異常進行擬合獲得的高程異常值，可以滿足一定范圍內(nèi)的測繪工作需求。

1 原理與方法

1.1 技術(shù)流程

正常高是以似大地水準面為基準面的，兩者之間的差值為高程異常ζ，三者關(guān)系如下：

h=H-ζ。

(1)

當前建立一個區(qū)域的似大地水準面模型，流行的方法是移去-恢復法，即對測區(qū)范圍內(nèi)離散的控制點進行GNSS水準聯(lián)測后，把控制點的高程異常ζ分解成中長波分量(ζGM)、短波分量(ζT)和殘余部分(ζres)，其關(guān)系如下：

ζ=ζGM+ζT+ζres。

(2)

其中，用重力場模型高程異常來表示中長波部分高程異常，剩余地形模型計算短波分量高程異常。

在計算時先把中長波和短波部分移去，再對殘余部分用數(shù)學或統(tǒng)計的方法進行擬合，然后將擬合結(jié)果與移去的中長波和短波部分相疊加，得到最終的高程異常，最后擬合出該區(qū)域的似大地水準面模型。

在似大地水準面精化研究中，劉瑛等[6]對多個重力場模型進行了對比，研究表明EGM2008重力場模型精度要優(yōu)于其他模型，并肯定了其在中國大陸區(qū)域的適用性。本研究亦選用EGM2008重力場模型。

在不考慮短波分量部分時，可以將短波分量部分歸入殘余高程異常，即

ζ=ζGM+ζres。

(3)

在對殘余高程異常部分進行擬合時，移動曲面擬合模型是日常測繪中較為常用的一種方法。

1.2 計算高程異常長波分量

EGM2008重力場模型是美國國家地理空間情報局(NGA)在充分利用GRACE衛(wèi)星重力、TOPEX衛(wèi)星測高、5′×5′地面重力等多種數(shù)據(jù)的基礎(chǔ)上，研制發(fā)布的2 190次球諧系數(shù)的全球超高階地球重力場模型。該模型的基本格網(wǎng)分辨率為5′，目前格網(wǎng)分辨率已經(jīng)達到2.5′×2.5′和更加精細的1′×1′，是當前計算重力高程異常的理想模型。目前可以通過NGA網(wǎng)站下載格網(wǎng)分辨率1′×1′的數(shù)據(jù)，利用Alltrans EGM2008 Calculator軟件可以直接計算高程異常ζGM。

1.3 殘余高程異常項處理

殘余高程異常項擬合方法較多，主要有二次曲面函數(shù)法、多面函數(shù)法、神經(jīng)網(wǎng)絡法、移動曲面法，每種方法都有其優(yōu)缺點。比如：二次曲面法僅能滿足地形起伏小且控制區(qū)域小的一般工程需求；多面函數(shù)由于其參數(shù)的選取而對擬合精度影響較大且工作量也大；神經(jīng)網(wǎng)絡法的算法復雜，在工程應用中使用神經(jīng)網(wǎng)絡法尚有諸多問題有待進一步研究[7]。考慮到模型復雜性、精度和可靠性，本研究的殘余高程異常擬合采用距離定權(quán)的移動曲面模型。該模型是在每一個待求點的鄰域半徑R范圍內(nèi)擬合出一個曲面，同時引入距離定權(quán)的權(quán)函數(shù)，在該區(qū)域內(nèi)，已知點與待求點距離越近權(quán)越大，反之越小[8-11]。移動曲面模型又包含平面法、二次曲面法和加權(quán)平均法。

1.3.1平面法

平面法模型如下：

ζresi=a0+a1xi+a2yi，i=0，1，2，…，n。

(4)

平面法要求已知點大于3個，利用最小二乘法的參數(shù)平差原理求解未知數(shù)a0、a1、a2，則式(4)的矩陣形式如下：

V=AX+L，

(5)

根據(jù)最小二乘原理VTPV可得

ATPAX+ATPL=0。

(6)

設(shè)ATPA=N、ATPL=U，則式(6)可改寫成如式(7)的法方程：

NX+U=0，

(7)

式中：N為系數(shù)陣；U為常數(shù)項。當N為非奇異陣時，可以求出未知數(shù)X：

X=-N-1U。

(8)

將已知數(shù)據(jù)代入函數(shù)，解出未知數(shù)a0、a1、a2，這樣就確定了一個擬合平面。把待求點的平面坐標代入式(4)，可求出待求點的殘余高程異常ζresi，其與重力場模型高程異常ζGMi之和即為待求點的高程異常ζi。最后，將ζi代入式(1)可得到正常高hi。

1.3.2二次曲面法

二次曲面法模型如下：

(9)

與平面法模型類似，由大于6個點的已知點可以求解出未知數(shù)a0、a1、a2、a3、a4、a5，擬合得到二次曲面，進而求解出殘余高程異常ζresi，再恢復重力場模型高程異常ζGMi，最后得到正常高hi。

1.3.3加權(quán)平均法

加權(quán)平均法模型待求點的殘余高程異常ζres可由周圍已知點計算，公式如下：

(10)

式中：n為設(shè)定半徑區(qū)域內(nèi)的已知點個數(shù)；Pi為各待求點的權(quán)。求出ζres后，疊加待求點的重力高程異常項，即可求出待求點的高程異常值。

權(quán)P是與待求點、已知點的距離有關(guān)的值，具體關(guān)系如下：

(11)

式中：S為擬合半徑，是待求點與已知點的水平距離；e為定權(quán)常數(shù)，是任意常數(shù)，根據(jù)移動曲面模型的特點，本研究的定權(quán)常數(shù)e取0.01；n為距離定權(quán)次方，一般取大于零的數(shù)。

另外，為防止S的距離過小導致P值無限大，規(guī)定當S

1.4 精度評定

采用最小二乘法中的單位權(quán)中誤差作為高程擬合精度分析的標準，如式(12)所示：

(12)

式中：v表示已知高程異常(H-h)與擬合所得高程異常(ζGM+ζres)的差值；n為檢查點個數(shù)。

2 實例與精度評定

2.1 實驗區(qū)情況

漳州市位于福建省最南部，東臨臺灣海峽，東北部與廈門市和泉州市接壤，西北部與龍巖市相連，西南部與廣東省梅州市和潮州市毗鄰。漳州市主要由山地、丘陵、平原組成，地勢從西北向東南傾斜。

圖1 已知點和檢查點分布圖 Fig.1 Distribution map of given points and check points

圖2 實驗區(qū)高程異常暈渲圖

利用漳州市C級GNSS控制網(wǎng)303個控制點作為實驗數(shù)據(jù)，所有控制點均包含CGCS2000坐標系下平面坐標、大地高和1985國家高程基準下的正常高成果，其中274個控制點作為已知點、29個控制點作為檢查點，已知點和檢查點的分布如圖1所示。根據(jù)303個控制點數(shù)據(jù)求得的高程異常，繪制漳州市高程異常暈渲圖，如圖2所示。

2.2 高程異常長波項計算

在高程異常長波項計算過程中，將303個已知點用Alltrans EGM2008 Calculator軟件計算長波項高程異常ζGM，表1列舉了部分已知GNSS水準點的重力場模型ζGM、大地高H、正常高h等信息。殘余高程異常暈渲圖見圖3，殘余高程異常值最小為0.092 m，最大為1.08 m。

表1 部分已知GNSS水準點信息Tab.1 Partially given GNSS bench marking point information m

圖3 殘余高程異常暈渲圖Fig.3 Dizziness map of residual elevation abnormal

2.3 殘余高程異常擬合

依次對274個已知點用移動曲面模型法中的平面法、二次曲面法和加權(quán)平均法進行殘余高程異常值擬合，其中待求點視為已知點的距離d取100 m，定權(quán)常數(shù)e取0.01，距離定權(quán)次方n依次取1、2、3，移動曲面擬合搜索半徑依次取25 km、30 km、35 km。

2.4 檢查點精度評價

為檢驗上述3種方法及不同的擬合參數(shù)求出的待定點高程異常值的精度，需要對待定點擬合高程異常值進行精度評定。本研究對前面選取的29個檢查點進行檢驗計算，具體實測高程異常值與擬合高程異常值較差的最小值(取較差絕對值)、最大值(取較差絕對值)和中誤差如表2至表4所示。

表2 檢查點移動曲面之平面法精度統(tǒng)計Tab.2 Accuracy statistics of plane method of check point moving surface

表3 檢查點移動曲面之二次曲面法精度統(tǒng)計Tab.3 Accuracy statistics of quadric surface method for checkpoint moving surface

表4 檢查點移動曲面之加權(quán)平均法精度統(tǒng)計Tab.4 Accuracy statistics of weighted average method of checkpoint moving surface

從表2至表4可以看出，29個檢查點的實測高程異常值與擬合高程異常值較差，采用二次曲面法，距離定權(quán)次方n=2和n=3時效果較好，最大值均小于10 cm，中誤差均小于4 cm，完全滿足CJJ/T 8—2011《城市測量規(guī)范》[13]中圖根點的精度要求。

3 結(jié)語

采用移去-恢復法的高程擬合思想，用EGM2008模型計算高程異常長波分量，利用在地形起伏較大且控制面積較大的區(qū)域效果較好的移動曲面模型進行殘余高程異常值的擬合。實踐表明，該模型是一種有效的GNSS高程擬合方法，特別是利用移動曲面之二次曲面法，根據(jù)實驗區(qū)內(nèi)的檢查點，實測高程異常值與擬合高程異常值較差最大值小于10 cm、中誤差小于4 cm。