余 騰 胡伍生 吳 杰 孫小榮 趙升峰

(1.宿遷學(xué)院建筑工程學(xué)院，江蘇宿遷 223800；2.東南大學(xué)交通學(xué)院，江蘇南京 210096； 3.南京市測繪勘察研究院有限公司，江蘇南京 210019)

近年來，伴隨GNSS(Global Navigation Satellite System)技術(shù)的快速發(fā)展和廣泛應(yīng)用，在鐵道沿線建立GNSS帶狀控制網(wǎng)已成常態(tài)。借用GNSS控制點(diǎn)與道路沿線部分水準(zhǔn)點(diǎn)重合，采用GNSS技術(shù)獲取道路沿線一定數(shù)量控制點(diǎn)的高程數(shù)據(jù)，并利用其幾何水準(zhǔn)測量高程資料，選擇科學(xué)的數(shù)學(xué)方法進(jìn)行擬合計(jì)算并對精度進(jìn)行分析[1-2]。

GNSS測量得到的高程是以WGS-84參考橢球?yàn)榛鶞?zhǔn)的大地高，而實(shí)際工程中使用的是以似大地水準(zhǔn)面為基準(zhǔn)的正常高，同一點(diǎn)位兩個(gè)數(shù)值的差值即為高程異常值ξ；如何選擇合適的函數(shù)模型對數(shù)據(jù)進(jìn)行擬合，求得每個(gè)點(diǎn)的高程異常值，進(jìn)而準(zhǔn)確地將大地高轉(zhuǎn)化成正常高是研究者關(guān)注的重點(diǎn)，也是實(shí)現(xiàn)精確的GNSS高程擬合的關(guān)鍵[3-6]。

1 高程系統(tǒng)基礎(chǔ)理論

1.1 水準(zhǔn)面基本概念

在大地測量與工程測量中，通常會涉及水準(zhǔn)面、大地水準(zhǔn)面、似大地水準(zhǔn)面幾個(gè)概念[7]。

水準(zhǔn)面是基于地球自轉(zhuǎn)的慣性離心力和地球引力交叉產(chǎn)生的重力作用而形成的一個(gè)處處與重力方向垂直的連續(xù)曲面，也是一個(gè)重力等位面。

大地水準(zhǔn)面是一個(gè)與平均海水面重合并延伸到大陸內(nèi)部形成的水準(zhǔn)面。因地表起伏和地球內(nèi)部質(zhì)量分布不勻，故大地水準(zhǔn)面是一個(gè)略有起伏的不規(guī)則曲面。該面包圍的形體近似于一個(gè)旋轉(zhuǎn)橢球，稱為“大地體”，大地水準(zhǔn)面是最接近地球整體形狀的重力位水準(zhǔn)面，也是正高系統(tǒng)的高程基準(zhǔn)面。

似大地水準(zhǔn)面是從地面點(diǎn)沿正常重力線量取正常高所得端點(diǎn)構(gòu)成的封閉幾何曲面。由于正高與大地水準(zhǔn)面的確定涉及到地球內(nèi)部密度的假定，在理論上存在著不嚴(yán)密性，為了便于計(jì)算，原蘇聯(lián)專家莫洛金斯基提出似大地水準(zhǔn)面，可以應(yīng)用地面測量數(shù)據(jù)直接確定地球表面形狀而不需要對地球密度作任何假設(shè)。似大地水準(zhǔn)面接近于水準(zhǔn)面，沒有物理意義，只是一個(gè)輔助面。

1.2 主要高程系統(tǒng)

(1)大地高是地面點(diǎn)沿參考橢球面法線到參考橢球面的距離，它克服了難以獲得地球內(nèi)部準(zhǔn)確質(zhì)量數(shù)據(jù)的困難，各個(gè)區(qū)域建立自己的參考橢球也有助于獲得本區(qū)域較高精度的高程值。地球上任意位置處的大地高指沿經(jīng)由待測點(diǎn)的旋轉(zhuǎn)球體的法線到旋轉(zhuǎn)體之間的長度。大地高沒有實(shí)際意義，只是數(shù)學(xué)層面上的定義。GNSS測量得到的高程是以WGS-84參考橢球?yàn)榛鶞?zhǔn)的大地高[8]。在WGS-84世界大地坐標(biāo)系下，地表P點(diǎn)的大地高用Hp表示(如圖1)。

圖1 WGS-84 世界大地坐標(biāo)系與大地高

(2)正高系統(tǒng)是以大地水準(zhǔn)面為基準(zhǔn)面的高程系統(tǒng)。可表示為

(1)

其中g(shù)m是點(diǎn)沿垂線到相應(yīng)基準(zhǔn)面間的平均重力值，dh是該方向上的高差，因?yàn)間m和該未知點(diǎn)所處地方的密度和深度有關(guān)，其密度深度數(shù)據(jù)難以精確測量[9]，因此，正高不能準(zhǔn)確得到，把大地水準(zhǔn)面與橢球面間的距離用N表示，有關(guān)系式

N=H-H正

(2)

(3)正高雖有完整的物理意義，鑒于每個(gè)區(qū)域的引力常數(shù)無法準(zhǔn)確測得，正高并不能準(zhǔn)確得到，為了便于計(jì)算研究，人們提出正常高系統(tǒng)。公式記為

(3)

式中，γm為平均正常重力值，則地表某點(diǎn)到大地水準(zhǔn)面和似大地水準(zhǔn)面間距離差為

(4)

式中，gm-rm為重力異常，似大地水準(zhǔn)面和參考橢球兩個(gè)基準(zhǔn)之間距離叫做高程異常，記為ξ

ξ=H-H常

(5)

總體來說，幾種主要高程系統(tǒng)關(guān)系可用圖2示意。

圖2 高程系統(tǒng)關(guān)系示意

2 GNSS高程擬合原理

直接獲取測區(qū)內(nèi)GNSS點(diǎn)的大地高后，選擇若干數(shù)量和位置均能滿足高程擬合要求的GNSS點(diǎn)，通過水準(zhǔn)測量方法對GNSS點(diǎn)進(jìn)行水準(zhǔn)聯(lián)測或在水準(zhǔn)點(diǎn)上布設(shè)GNSS點(diǎn)，獲取其正常高數(shù)據(jù)，計(jì)算點(diǎn)位的高程異常，選擇模型算法進(jìn)行高程擬合計(jì)算，求得測區(qū)內(nèi)其他GNSS點(diǎn)的正常高[10]。

先對GNSS觀測基線向量進(jìn)行解算，把坐標(biāo)作為未知參數(shù)進(jìn)行自由網(wǎng)平差，求出其三維地心坐標(biāo)，再在網(wǎng)中選取不少于3個(gè)聯(lián)測控制點(diǎn)(控制點(diǎn)經(jīng)度、緯度、大地高數(shù)據(jù)已知)。將這些點(diǎn)轉(zhuǎn)到相應(yīng)橢球的三維坐標(biāo)中，有

(6)

(7)

式中，Δx，Δy，Δz為平移參數(shù)，m是尺度比，εx，εy，εz為旋轉(zhuǎn)參數(shù)。求得GNSS控制點(diǎn)相應(yīng)坐標(biāo)系下的直角坐標(biāo)，再經(jīng)過如下變換

L=arctan(y/x)

(8)

(9)

可得到與已知控制點(diǎn)同一參考橢球下的經(jīng)緯度和大地高數(shù)據(jù)。在進(jìn)行坐標(biāo)系轉(zhuǎn)換時(shí)多采用7參數(shù)或4參數(shù)法，GNSS網(wǎng)點(diǎn)經(jīng)過轉(zhuǎn)換后和地面點(diǎn)間仍存在間隙，可在編制軟件過程中采取了一些措施來消除部分不兼容情況[12]。

3 GNSS高程擬合方法與精度評定

3.1 適用于道路工程的GNSS高程擬合方法選取

根據(jù)數(shù)據(jù)獲取和處理方法不同，GNSS高程可分為：GNSS三角高程、GNSS重力高程、GNSS水準(zhǔn)高程[13]。GNSS三角高程是利用已知的GNSS點(diǎn)高程數(shù)據(jù)，測出各點(diǎn)之間的豎直角和基線邊長，再通過三角高程計(jì)算公式求得它們之間的高差；GNSS重力高程測量利用測點(diǎn)區(qū)域重力資料，通過地球重力場模型求出點(diǎn)位的高程異常，用GNSS所測得的大地高減去高程異常，得到待測點(diǎn)的正常高。在實(shí)際工程中，重力數(shù)據(jù)相對匱乏，其適用性受較大限制。GNSS水準(zhǔn)高程是在獲得GNSS大地高和水準(zhǔn)測量正常高的基礎(chǔ)上，利用它們之間的轉(zhuǎn)換關(guān)系獲取其余點(diǎn)位正常高的方法，此方法方便高效，精度也有保障，適用范圍較廣。

常用GNSS高程擬合方法有：等值線圖內(nèi)插法、曲線擬合、曲面擬合(主要有多項(xiàng)式曲面擬合、多面函數(shù)擬合，適用于GNSS點(diǎn)呈面狀分布)、地球重力場模型法、神經(jīng)網(wǎng)絡(luò)法等。這類利用數(shù)學(xué)模型來計(jì)算區(qū)域內(nèi)高程異常值方法的主要工作就是構(gòu)建一個(gè)與測區(qū)似大地水準(zhǔn)面接近程度相對較高的數(shù)學(xué)模型。

在選用何種擬合方法時(shí)，應(yīng)充分考慮GNSS點(diǎn)的數(shù)量、密度和分布狀況[14]；對于道路工程呈帶狀狹長區(qū)域的特點(diǎn)，GNSS點(diǎn)呈線狀分布，沿途的似大地水準(zhǔn)面屬于連續(xù)的曲線，采用曲線擬合方法更為合適。曲線擬合的原理是：由GNSS控制點(diǎn)二維數(shù)據(jù)和高程異常來構(gòu)造一個(gè)插值函數(shù)，求出函數(shù)的必要參數(shù)，由GNSS所得二維數(shù)據(jù)來得到高程異常，最后將GNSS所得大地高減去高程異常算得正常高。

曲線擬合方法主要有：多項(xiàng)式曲線擬合、三次樣條曲線擬合、Akima曲線擬合。Akima曲線是在兩個(gè)實(shí)測點(diǎn)之間進(jìn)行內(nèi)插，它需要知道這兩個(gè)實(shí)測點(diǎn)觀測值和與這兩個(gè)實(shí)測點(diǎn)相近鄰的四個(gè)實(shí)測點(diǎn)上的觀測值，此法在已知點(diǎn)數(shù)量不多時(shí)并不適用，其擬合雖求解效率高，但忽略了一些污染數(shù)據(jù)的影響，有片面性，其擬合后線型平滑，但求導(dǎo)誤差較大[15]。多項(xiàng)式曲線擬合是線狀分布擬合的主要方法，三次樣條曲線在路線長和控制點(diǎn)多時(shí)可構(gòu)造高次的擬合多項(xiàng)式而避免過于震蕩，故以下著重研究多項(xiàng)式曲線擬合、三次樣條曲線擬合兩種方法。

3.2 GNSS高程擬合精度評定

常用的GNSS高程擬合精度評價(jià)主要有內(nèi)符合精度、外符合精度等[16]。

(1)內(nèi)符合精度

(10)

(2)外符合精度

(11)

4 工程算例

4.1 工程概述

工程算例數(shù)據(jù)來源于圖3中的藍(lán)色某段區(qū)間(吳集至龍廟)，路線總長約30.95 km。經(jīng)搜集已知水準(zhǔn)高程資料和現(xiàn)場踏勘確認(rèn)，選定沿線13個(gè)控制點(diǎn)，已知水準(zhǔn)高程由四等水準(zhǔn)觀測獲得，經(jīng)過GNSS接收機(jī)布網(wǎng)進(jìn)行靜態(tài)觀測，得到橢球高、正常高，出于保密需要，高程數(shù)據(jù)經(jīng)過系統(tǒng)處理，見表1。

圖3 鐵路整體走向示意

點(diǎn)號橢球高/m正常高/m點(diǎn)間距/km1216.998217.3422209.521 209.8773148.947149.3054154.963155.3175137.465137.8336127.350127.7207160.069160.4718150.886151.2989112.937113.35110118.143118.54811155.251155.67512163.171163.59413133.587134.0170.472.620.584.630.875.260.485.670.744.560.714.36

將數(shù)據(jù)繪制成如圖4所示的形式。

圖4 道路控制點(diǎn)高程折線示意

4.2 多項(xiàng)式曲線擬合

多項(xiàng)式曲線擬合把測區(qū)看成一條連續(xù)曲線，其本質(zhì)是一個(gè)m次方的代數(shù)多項(xiàng)式?？稍O(shè)高程異常值ξ與任意點(diǎn)A(x，y)有以下關(guān)系

ξ=f(x)-v

(12)

式中，f(x)為高程異常的數(shù)學(xué)模型，v為誤差。

設(shè)

f(x)=a0+a1x+a2x2+a3x3+……

(13)

當(dāng)有n個(gè)已知點(diǎn)時(shí)，上式寫成矩陣形式

V=XB-ξ

(14)

式中

(15)

在Σν2=min限制條件下，依B=(XTX)-1XTξ，可求出向量B和待定系數(shù)a0、a1、a2…an；即可計(jì)算出未知點(diǎn)高程異常值[17]。數(shù)據(jù)分別采用二次和三次多項(xiàng)式擬合，過程和結(jié)果大致如下。

①二次多項(xiàng)式擬合

依ξ=a0+a1x+a2x2，選取1，5，10，13作為公共點(diǎn)，其余作為檢核點(diǎn)，計(jì)算得常數(shù)值：ξ=-0.343 9-0.003 004x+0.000 007 014x2

帶入公式后可得各點(diǎn)高程異常和殘差。

②三次多項(xiàng)式擬合

依ξ=a0+a1x+a2x2+a3x3，選取同樣公共點(diǎn)和檢核點(diǎn)，算得常數(shù)值：ξ=-0.3463-0.001 867x+0.000 157 4x2+0.000 004 39x3

同理，帶入公式后可得各個(gè)點(diǎn)的高程異常值和殘差值。

4.3 三次樣條曲線擬合

當(dāng)高程異常值波動大、測量線路較長且有一定數(shù)量公共點(diǎn)的情況下，求待定系數(shù)的誤差會變大?？砂褱y線看成若干個(gè)分線段，將每段設(shè)為三次多項(xiàng)式函數(shù)，然后將多項(xiàng)式函數(shù)組成三次樣條函數(shù)。處理后不僅函數(shù)自身連續(xù)，其一階、二階導(dǎo)數(shù)也連續(xù)[18]。

ξ(x)=ξ(xi)+(x-xi)ξ(xi，xi+1)+

(x-xi)(x-xi+1)ξ(x，xi，xi+1)

(16)

其中x為待定點(diǎn)坐標(biāo)，ξ(xi，xi+1)是一階差商，ξ(xi，xi+1)=(ξi+1-ξi)/(xi+1-xi)；ξ(x，xi，xi+1)為二階差商，ξ(x，xi，xi+1)=1/6[ξ″(xi)+ξ″(x)+ξ″(xi+1]，ξ″(X)=(i=1，2，…，n-1)，其線性方程組具有對稱三角陣的系數(shù)矩陣

(xi-xi+1)ξ″(xi-1)+2(xi+1，xi-1)ξ″(xi)+

(xi+1，xi)ξ″(xi+1)=6[ξ″(x，xi+1)-ξ(x，xi)]

(17)

ξ(x0)=ξ″(xn)=0

(18)

用追趕法對其求解，可得ξ″(xi)和ξ(xi，xi+1)

ξ″(x)=ξ″(xi)+(x-xi)ξ″(xi，xi+1)

(19)

ξi，i+1=(ξi+1-ξi)/(xi+1-xi) (i=1，2,…，n-1)

(20)

Ki，i+1=(Ki+1-Ki)/(xi+1-xi) (i=1，2，…，n-1)

求出系數(shù)Ki后，就可用下式求得測區(qū)內(nèi)任一點(diǎn)的高程異常，并計(jì)算各點(diǎn)正常高

ξ=ξi+ξi，i+1(x-xi)+(x-xi)·

(x-xi+1){2Ki+Ki+1+Ki，i+1(x-xi)}/6

以采用7個(gè)控制點(diǎn)分為5段為例，設(shè)置每段的三次樣條函數(shù)(如圖5～圖9)，橫坐標(biāo)為5段中某分段的點(diǎn)位，縱坐標(biāo)為Matlab軟件自動計(jì)算的序列點(diǎn)特征量，擬合過程大致如下。

圖5 第一段樣條曲線

同理，得到第二段曲線函數(shù)

ξ=-0.001 5t3+0.013 5t2-0.031 2t-0.344 8

圖6 第二段樣條曲線

擬合得到第一條曲線函數(shù)

ξ=-0.001 6t3+0.009 7t2-0.008 8t-0.353 8

第三段曲線函數(shù)

ξ=0.000 4t3-0.007 6t2+0.046 3t-0.439 6

圖7 第三段樣條曲線

第四段曲線函數(shù)

ξ=-0.000 3t3+0.009 8t2-0.098 6t-0.038 7

圖8 第四段樣條曲線

第五段曲線函數(shù)

ξ=-0.000 3t3+0.009 4t2-0.0948t-0.0505

圖9 第五段樣條曲線

4.4 兩種方法擬合殘差比較

對道路控制點(diǎn)高程用多項(xiàng)式、三次多項(xiàng)式、三次樣條曲線方法進(jìn)行擬合分析，擬合值殘差如表2所示。

表2 多項(xiàng)式擬合殘差 m

可繪制成如圖10所示的形式。

圖10 高程擬合殘差

5 結(jié)論

本段線路長度約為30.95 km，相較于傳統(tǒng)水準(zhǔn)測量，選用準(zhǔn)確高效的GNSS高程擬合方法能極大地節(jié)省人力和物力。

選取多項(xiàng)式擬合和三次樣條曲線擬合方法進(jìn)行鐵路GNSS高程擬合適用且有效；對13個(gè)控制點(diǎn)的測試表明，其擬合誤差最大點(diǎn)的殘差值為-0.071 m，若經(jīng)簡易整體平差改正數(shù)按比例反向賦配，改正后的水準(zhǔn)高程誤差均在0.04 m以內(nèi)，對于Ⅱ級及以下等級鐵路基本可以滿足道路工程測量中地形測量、斷面測量和施工測量對于控制點(diǎn)高程的精度要求。

經(jīng)數(shù)據(jù)分析比較，二次多項(xiàng)式、三次多項(xiàng)式、三次樣條曲線三種方法的最大擬合殘差依次為：-0.071 m，-0.059 m，-0.042 m；三種方法的殘差均值依次為：-0.019 m，-0.015 m，-0.011 m；依據(jù)上述精度評定方法，三種方法的內(nèi)符合精度依次為0.066 m，0.055 m，0.045 m；外符合精度依次為0.070 m，0.059 m，0.044 m。由此可見，兩種GNSS高程擬合方法的點(diǎn)位擬合殘差平均值都在0.02 m以內(nèi)，擬合精度較好，且三次樣條曲線擬合效果更優(yōu)。究其原因，(1)它保留了多項(xiàng)式表達(dá)簡便的特點(diǎn)，又解決了單個(gè)多項(xiàng)式機(jī)械、不靈活的缺點(diǎn)；(2)當(dāng)線路較長時(shí)，即便是高次多項(xiàng)式，其擬合程度也不能明顯提高，而選取更多的插值節(jié)點(diǎn)，采用分段低次插值反而可行[19]。

對于較短或者更長距離的高程擬合還有待進(jìn)一步探討，其他插值方法和組合方法對于鐵道工程GNSS高程擬合的研究還需進(jìn)一步深入[11，14，15，19]。