柯生學(xué)

（中國電建集團(tuán)西北勘測設(shè)計研究院有限公司，陜西 西安 710100）

0.引言

傳統(tǒng)水準(zhǔn)測量精度高，可靠性強，但往往耗費大量人力物力，測量速度慢，且受天氣氣候約束[1-4]。GNSS測量能在布網(wǎng)控制點上獲得正高，但我國使用的是正常高系統(tǒng)（基于似大地水準(zhǔn)面），為使測量結(jié)果能直接使用，必須計算高程異常值[5,6]。高程異常擬合需要獲得測區(qū)一定數(shù)量的控制點坐標(biāo)、大地高、正常高，若使用重力場模型還需要獲得對應(yīng)重力場模型計算的長波項，根據(jù)區(qū)域高程、高程異常值等條件選擇適合本區(qū)域的擬合模型[7-10]。傳統(tǒng)高程異常擬合模型受限于數(shù)學(xué)模型的特點，難以達(dá)到較高精度[11-13]。EGM2008 重力場模型作為美國發(fā)布的新一代重力場模型，可以精確計算待定點重力項高程異常，將之與移去-恢復(fù)法及常規(guī)數(shù)學(xué)模型結(jié)合，可以極大地提高高程異常擬合精度[14,15]。本文以西安某區(qū)域GNSS 水準(zhǔn)數(shù)據(jù)為例，通過多項式擬合法、最小二乘配置擬合法、半?yún)?shù)模型結(jié)合移去-恢復(fù)討論小區(qū)域高程異常擬合的最優(yōu)模型[16,17]。

1.高程異常擬合的理論與方法

高程異常擬合有多種理論，本文主要采用以下三類：多項式擬合法、最小二乘配置法、半?yún)?shù)模型。

1.1 多項式擬合法

多項式擬合法是通過選取一定的擬合系數(shù)數(shù)量來控制模型起伏的一種數(shù)學(xué)模型，最少選擇三個擬合系數(shù)，這時模型就是平面擬合模型。隨著選擇的擬合參數(shù)增加，模型起伏會有所變化，所需要的最少起算點也會隨之增加。一般選擇六參數(shù)多項式。擬合參數(shù)并非多多益善，因為擬合參數(shù)的增加會引起模型曲面的復(fù)雜起伏變化，有時反而會降低擬合精度。

多項式擬合模型可用式（1）表示：為高程異常；ε為殘差；ai為待擬合參數(shù)；B，L為大地經(jīng)緯度。誤差方程式（2）為：

擬合系數(shù)矩陣

根據(jù)以上式可知多項式擬合至少需要六組已知點的高程異常及大地經(jīng)緯度。將公共點數(shù)值帶入上式，可得待估擬合參數(shù)X 計算公式（3）為：

1.2 最小二乘配置法

最小二乘配置法是一種優(yōu)化模型，需要基于某個基礎(chǔ)模型進(jìn)行精度，一般優(yōu)化六參數(shù)多項式。六參數(shù)多項式僅考慮到模型趨勢項影響，未考慮非趨勢項即隨機誤差的影響。最小二乘配置法將非趨勢項作為優(yōu)化目標(biāo)，通過計算提高模型擬合精度。

最小二乘配置法可用公式（4）表示：

式中，L為已知點高程異常值，X為函數(shù)模型部分待擬合參數(shù)，B為函數(shù)模型部分系數(shù)矩陣；S為觀測信號；Δ 為觀測信號噪聲。以S'表示未知點高程異常觀測信號，式（5）可表示為：

式中，DLL=DΔΔ+CDZZCT；未測點的最或然值見式（9）：

1.3 半?yún)?shù)模型

半?yún)?shù)模型可用式（10）表示：

式中，L為已知點高程異常；B，X為參數(shù)模型函數(shù)系數(shù)矩陣與待估參數(shù)；S為半?yún)?shù)分量；Δ 為誤差分量。誤差方程見式（11）：

式中，P為權(quán)陣（正定矩陣）；α 為在極化過程中平衡V與S的純量因子，稱平滑因子；R為給定的正則化矩陣。由式（10）構(gòu)造極值函數(shù)，可解參數(shù)分量與非參數(shù)分量

式中，N=BTPB；M=P+αR-PBN-1BTP。

依據(jù)最小距離法選取平滑因子α，αD為曲線上的點到原點距離最小時對應(yīng)的α，當(dāng)Sn2(αD) +Vn2(αD)=min時，α取得最優(yōu)值，非參分量最優(yōu)值經(jīng)平差計算獲得。利用式計算傾向參數(shù)Y，結(jié)合S計算出X，內(nèi)插獲得待算點的高程異常值見式（15）：

2.高程異常擬合優(yōu)化基于—EGM2008 的移去恢復(fù)

高程異?？煞譃橐韵聨讉€部分：由重力場模型確定的中長波部分ξG、由Stokes理論確定的重力異常殘余分量ξS、由地形改正得到的短波部分ξD?！耙迫?恢復(fù)”法的基本思想是：移去高程異常中的長波部分或者短波部分，然后將剩余部分結(jié)合擬合模型進(jìn)行擬合計算，最后在待定點上恢復(fù)被移去的部分。由于缺少高分辨率的DTM數(shù)據(jù)，本文將短波部分與殘余分量合并作為殘差高程異常ξΔ=ξS+ξD。計算過程如下：

（1）根據(jù)GNSS 水準(zhǔn)數(shù)據(jù)計算每點高程異常ξ'；

（2）將ξ'與以EGM2008 模型計算的中長波部分ξG相減，得殘差高程異常ξΔ；

（3）選擇擬合模型對ξΔ進(jìn)行擬合，確定最優(yōu)待定參數(shù)，建立殘差高程異常擬合模型；

（5）將外推所得ξΔ'恢復(fù)中長波部分,得該點高程異常擬合值

（6）對模型進(jìn)行內(nèi)外符合精度評定以及殘差計算。

3.算例

3.1 工程概況

本文所用數(shù)據(jù)來自西安某區(qū)域，南北長約55 km，東西寬約30 km，地面標(biāo)高為+35.38m～+36.29 m。礦區(qū)共有3 個C 級GNSS 點、19 個D 級GNSS 點和24 個E級GNSS 點，總計46 個GNSS 點都進(jìn)行了相應(yīng)等級的水準(zhǔn)聯(lián)測，礦區(qū)GNSS 點分布（如圖1 所示）。

由圖1 可知，已知點分布不均勻，因此以3 個C 級點做起算點，在礦區(qū)已有基礎(chǔ)上，新增10 個GNSS 點進(jìn)行D 級GNSS 測量，大地高與正常高觀測結(jié)果（如表1 所示），由表可知所有點大地高變化不大，點位穩(wěn)定可靠。

3.2 方案設(shè)計

根據(jù)測區(qū)概況，共設(shè)計六種計算方案，分別為：（1）六參數(shù)多項式；（2）六參數(shù)多項式+移去恢復(fù)；（3）最小二乘配置；（4）最小二乘配置+移去恢復(fù)；（5）半?yún)?shù)；（6）半?yún)?shù)+移去恢復(fù)。分別計算六種方案的擬合殘差與內(nèi)外符合精度進(jìn)行對比。

3.3 實例計算與精度分析

圖1 礦區(qū)已知點分布

圖2 多項式殘差對比

表1 GNSS 點可靠性分析

續(xù)：表1

圖3 最小二乘配置殘差對比

圖4 半?yún)?shù)殘差對比

表2 計算結(jié)果

六種方案計算結(jié)果見表1，其中MAX 表示殘差最大值，MIN 表示殘差最小值，RMSE 表示內(nèi)/外符合精度結(jié)果。 殘差對比（如圖3、圖4、圖5 所示）：由表2 可知，半?yún)?shù)+移去恢復(fù)方案內(nèi)外符合精度均較小，分別為0.0038m 和0.0091m，六參數(shù)多項式模型內(nèi)外符合精度較大，分別為0.0246m 和0.0311m；六參數(shù)多項式出現(xiàn)了最大內(nèi)符合殘差，為0.0695m，最小二乘配置法出現(xiàn)了最大外符合殘差，為0.0650m，最小內(nèi)外符合殘差均出現(xiàn)在半?yún)?shù)模型，為-0.0618m 和-0.0688m。若不加入移去-恢復(fù)技術(shù)，半?yún)?shù)模型相對比其他兩種模型并無明顯擬合優(yōu)勢；最小二乘配置法相對六參數(shù)多項式的提升也比較有限。加入移去-恢復(fù)技術(shù)后半?yún)?shù)模型擬合精度提升比較明顯，內(nèi)外符合精度均在毫米級，可以很好地用于實際測量工作中。

4.結(jié)束語

通過三種方案計算表明，移去恢復(fù)技術(shù)可以有效地提高高程異常擬合精度，但在多項式與最小二乘配置中提升有限，在半?yún)?shù)模型中提升較為明顯；半?yún)?shù)結(jié)合重力場模型的方案效果較好，最大內(nèi)符合殘差為0.0.471m，最大外符合殘差為0.0153m，最小內(nèi)符合殘差為-0.0133m，最小外符合殘差為-0.0078m，最適用于區(qū)域高程異常擬合。在實際應(yīng)用中，若對高程異常求解要求精度不高，可直接采用多項式模型，其計算簡單、不需要使用重力數(shù)據(jù)可以快速計算待定點概略高程異常；若需要高精度高程異常則需使用半?yún)?shù)+移去恢復(fù)的方法，以獲得較好的擬合效果。