周 偉

（1. 山東省臨沂市水利勘測(cè)設(shè)計(jì)院，山東 臨沂 276002）

網(wǎng)絡(luò)RTK技術(shù)的發(fā)展為采集數(shù)據(jù)提供了方便，通過(guò)流動(dòng)站GPS的觀測(cè)可以得到任意位置的平面坐標(biāo)（如2000國(guó)家大地坐標(biāo)）和大地高。似大地水準(zhǔn)面的高程可通過(guò)似大地水準(zhǔn)面精化的成果獲得，也可通過(guò)水準(zhǔn)測(cè)量擬合的方式來(lái)獲得[1-4]。由于似大地水準(zhǔn)面精化的成果一般都是保密的，很難獲得，因此利用插值和擬合方法求解高程異常就成為了必須的測(cè)量工作。在工程測(cè)量工作中，經(jīng)常采用工地校正或選取幾個(gè)點(diǎn)解算高程異常的方法，但在覆蓋區(qū)域較大的情況下（如臨沂市沂河及其支流水位的推算），該方法將給工作帶來(lái)不便?，F(xiàn)有的一些GPS自帶的處理軟件在求解轉(zhuǎn)換參數(shù)上會(huì)有離散點(diǎn)數(shù)量上的約束（如天寶GPS最多為19個(gè)點(diǎn)），當(dāng)點(diǎn)過(guò)多時(shí)則無(wú)法求解，需對(duì)項(xiàng)目進(jìn)行分塊求解，但利用不同的項(xiàng)目分塊將導(dǎo)致模型的不一致，在后期檢查、施工時(shí)也將加大工作難度。

為了獲取臨沂市整個(gè)區(qū)域的高程異常分布規(guī)律，本文介紹了多項(xiàng)式插值曲面、三次樣條函數(shù)插值曲面（Spline插值）、三次B樣條插值曲面以及多項(xiàng)式擬合曲面數(shù)學(xué)模型[5-8]，并利用臨沂市國(guó)土局提供的兩套共203個(gè)平高控制點(diǎn)（1980西安坐標(biāo)系/2000國(guó)家大地坐標(biāo)系，1985國(guó)家高程基準(zhǔn)的水準(zhǔn)高/大地高，平高同點(diǎn)）計(jì)算高程異常；再分別利用插值和擬合模型對(duì)高程異常進(jìn)行求解，分析這些模型與臨沂市精化GPS數(shù)據(jù)的一致性。

1 高程異常模型的計(jì)算方法

1.1 插值模型

高程異常插值曲面是對(duì)已知高程異常點(diǎn)的插值，得到未知區(qū)域的高程異常，本文主要討論多項(xiàng)式插值曲面和樣條函數(shù)插值曲面[5-8]模型的數(shù)學(xué)表示方法。

1.1.1 多項(xiàng)式插值曲面

1）線性插值曲面。該方法比較簡(jiǎn)單，計(jì)算公式為：

首先構(gòu)建三角網(wǎng)，并分別對(duì)兩個(gè)方向進(jìn)行求導(dǎo)，即可得到兩個(gè)方向的斜率，通過(guò)三角網(wǎng)的3個(gè)端點(diǎn)值就可求出方程，從而得到三角網(wǎng)內(nèi)任意點(diǎn)的值。該插值方法是通過(guò)已知端點(diǎn)得到的方程，運(yùn)算速度較快，得到的曲面不光滑，可用來(lái)查看高程異常曲面的插值效果。

2）高次多項(xiàng)式插值曲面。對(duì)于n＞1的情形，其計(jì)算公式為：

線性插值可表示為n=1的情形；當(dāng)n=3時(shí)，即三次多項(xiàng)式插值就需要求解10個(gè)未知數(shù)，次數(shù)越高，需要求解的未知數(shù)就越多。以n=3為例，計(jì)算未知數(shù)的過(guò)程依賴于插值數(shù)據(jù)的特性，常用方法為利用4個(gè)鄰近頂點(diǎn)的高度以及每個(gè)頂點(diǎn)的3個(gè)導(dǎo)數(shù)方程，一階導(dǎo)數(shù)分別表示兩個(gè)方向的表面斜率，二階導(dǎo)數(shù)表示同時(shí)在兩個(gè)方向的斜率，將頂點(diǎn)值代入多項(xiàng)式和導(dǎo)數(shù)方程，即可求解出這10個(gè)未知數(shù)。

1.1.2 三次樣條函數(shù)插值曲面

1）Spline插值。Spline插值是在多項(xiàng)式分段插值的基礎(chǔ)上發(fā)展而來(lái)的[5]，克服了曲線不光滑，只有一階導(dǎo)數(shù)連續(xù)的缺點(diǎn)。三次樣條函數(shù)的定義為：設(shè)則三次樣條函數(shù)且在每個(gè)區(qū)間為三次多項(xiàng)式，同時(shí)滿足三次樣條函數(shù)的精度依賴于節(jié)點(diǎn)的數(shù)量，若測(cè)量時(shí)觀測(cè)數(shù)據(jù)較少，則會(huì)引入誤差。該方法可表述為特定的方法求解各區(qū)間的三次多項(xiàng)式。

2）三次B樣條插值曲面。B樣條曲線是Bezier曲線的進(jìn)一步發(fā)展，也是樣條函數(shù)的進(jìn)一步發(fā)展。該方法克服了Bezier曲線階次過(guò)高、缺乏靈活性的缺點(diǎn)，可通過(guò)反求的方式得到經(jīng)過(guò)各端點(diǎn)并連續(xù)的方程[6-7]，目前已被廣泛應(yīng)用于工程制圖方面[8]。三次B樣條插值曲面模型如式（3）所示，其中0≤t≤1。

三次B樣條是B樣條為3次的特殊情況，通過(guò)相鄰4個(gè)頂點(diǎn)即可計(jì)算得到三次B樣條的曲面形狀；當(dāng)端點(diǎn)多于4個(gè)時(shí)則可通過(guò)反求的方式得到經(jīng)過(guò)各端點(diǎn)并連續(xù)的方程，從而得到區(qū)間內(nèi)的插值曲面模型。

1.2 多項(xiàng)式擬合曲面模型

多項(xiàng)式擬合曲面模型如式（2）所示，理想狀態(tài)下，已知點(diǎn)均通過(guò)擬合曲面；但實(shí)際上當(dāng)已知點(diǎn)的數(shù)量大于需要求解的未知參數(shù)時(shí)，總會(huì)有節(jié)點(diǎn)不通過(guò)曲面方程，即存在誤差。由于高程異常由兩個(gè)方向的值確定，即兩個(gè)自變量確定一個(gè)因變量的情況，為了得到區(qū)間內(nèi)最優(yōu)的曲面方程，就需要利用二元回歸的方法尋找一個(gè)最接近已知點(diǎn)的曲面，利用最小二乘方法進(jìn)行約束[5]，使高程異常誤差的均方根誤差最小，從而求出唯一最優(yōu)的曲面方程。

2 實(shí)驗(yàn)與結(jié)果分析

2.1 數(shù)據(jù)準(zhǔn)備

本文采用的數(shù)據(jù)包括臨沂市國(guó)土局提供的1980西安坐標(biāo)系的GPS C級(jí)點(diǎn)、II等三角點(diǎn)共203個(gè)，其中1985國(guó)家高程基準(zhǔn)的III等水準(zhǔn)點(diǎn)138個(gè)、GPS精化高程點(diǎn)65個(gè)（平高同點(diǎn)），以及2000國(guó)家大地坐標(biāo)系下這203個(gè)平高控制點(diǎn)。高程異常可表示為大地高減去正常高，利用已知數(shù)據(jù)求解得到203個(gè)點(diǎn)的高程異常值，作為精度統(tǒng)計(jì)的標(biāo)準(zhǔn)值。在現(xiàn)階段的測(cè)量工作中，推薦使用2000國(guó)家大地坐標(biāo)系，在實(shí)驗(yàn)中將2000國(guó)家大地坐標(biāo)作為自變量對(duì)高程異常曲面進(jìn)行插值和擬合。

2.2 數(shù)據(jù)實(shí)驗(yàn)

實(shí)驗(yàn)主要利用Matlab 2018b軟件進(jìn)行實(shí)現(xiàn)，首先分別利用線性插值、三次多項(xiàng)式插值和三次B樣條插值對(duì)138個(gè)III等水準(zhǔn)點(diǎn)進(jìn)行插值，插值曲面如圖1～3所示，并利用插值曲面函數(shù)對(duì)65個(gè)已知的GPS精化高程異常進(jìn)行計(jì)算，得到實(shí)際值與模型估算值的誤差；再分別利用1～5次多項(xiàng)式擬合方法對(duì)138個(gè)III等水準(zhǔn)點(diǎn)的高程異常曲面進(jìn)行求解，并利用求解的曲面模型對(duì)65個(gè)已知的GPS精化高程異常進(jìn)行計(jì)算，得到實(shí)際值與模型估算值的誤差。多項(xiàng)式擬合曲面如圖4～8所示。

圖1 線性插值曲面

圖2 三次多項(xiàng)式插值曲面

圖3 三次B樣條插值曲面

圖4 線性擬合曲面

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

圖6 三次多項(xiàng)式擬合曲面

圖7 四次多項(xiàng)式擬合曲面

圖8 五次多項(xiàng)式擬合曲面

2.3 數(shù)據(jù)分析

2.3.1 定性分析

1）通過(guò)圖2、3可直觀發(fā)現(xiàn)，多項(xiàng)式插值方法只對(duì)已知節(jié)點(diǎn)進(jìn)行插值，而樣條插值方法和擬合方法，在東坐標(biāo)和北坐標(biāo)的區(qū)間范圍內(nèi)均可估算高程異常曲面。

2）插值方法可以保證插值曲面經(jīng)過(guò)插值節(jié)點(diǎn)，但整體是不光滑的；擬合方法不能保證通過(guò)所有節(jié)點(diǎn)，但整體是比較光滑的。

3）多項(xiàng)式擬合曲面時(shí)，次數(shù)大于3時(shí)，在沒(méi)有點(diǎn)控制的邊緣區(qū)域?qū)a(chǎn)生較大的形變。

2.3.2 定量分析

1）擬合曲面的精度統(tǒng)計(jì)。由于插值方法通過(guò)插值節(jié)點(diǎn)，因此理論上插值曲面的節(jié)點(diǎn)誤差為零。擬合曲面不能保證所有節(jié)點(diǎn)都通過(guò)擬合曲面，存在誤差，多項(xiàng)式擬合曲面中誤差統(tǒng)計(jì)如表1所示。多項(xiàng)式擬合模型的高程精度至少可達(dá)0.1 m，利用三次多項(xiàng)式擬合時(shí)的高程精度為0.06 m。

表1 多項(xiàng)式擬合曲面中誤差統(tǒng)計(jì)/m

2）檢驗(yàn)點(diǎn)的精度統(tǒng)計(jì)。本文分別討論插值曲面和擬合曲面求解65個(gè)檢驗(yàn)點(diǎn)的精度，由于高次多項(xiàng)式擬合曲面在邊緣區(qū)域變形較大，且精度提升不明顯、多項(xiàng)式系數(shù)較多不方便計(jì)算，因此這里不再統(tǒng)計(jì)3～5次多項(xiàng)式擬合曲面的精度。插值和擬合方法精度統(tǒng)計(jì)如表2所示。

3）精度分析。由表1可知，多項(xiàng)式擬合方法可用于評(píng)價(jià)已知III等水準(zhǔn)數(shù)據(jù)的質(zhì)量，這里表明國(guó)土局提供的水準(zhǔn)數(shù)據(jù)質(zhì)量還是比較好的，精度優(yōu)于0.1 m。由表2可知，不同模型得到的精度是不同的，插值曲面的精度低于擬合曲面的精度；插值曲面并不是模型越復(fù)雜精度越高，從插值模型的精度統(tǒng)計(jì)來(lái)看，線性插值曲面精度優(yōu)于樣條插值曲面和三次多項(xiàng)式插值曲面；多項(xiàng)式次數(shù)不同擬合曲面的精度也不同，二次擬合曲面精度優(yōu)于線性擬合曲面精度；從樣條插值曲面和擬合曲面的最大殘差分布和中誤差統(tǒng)計(jì)情況來(lái)看，曲面的邊緣在沒(méi)有控制點(diǎn)的情況下，精度略低，為了達(dá)到更高的精度需適當(dāng)增加邊緣控制；二次多項(xiàng)式擬合曲面與臨沂市GPS精化高程的一致性最好，精度最高，檢查點(diǎn)中誤差為0.038 m，除了插值區(qū)域外的一個(gè)點(diǎn)，其余點(diǎn)的檢驗(yàn)殘差均優(yōu)于0.1 m，高程精度滿足1∶500～1∶2 000等大比例尺地形圖測(cè)量規(guī)范要求。

表2 插值和擬合曲面求解高程異常殘差統(tǒng)計(jì)

3 結(jié) 語(yǔ)

插值和擬合曲面模型是求解高程異常的一種方法，不同的區(qū)域可選擇不同的模型進(jìn)行高程異常解算。通過(guò)實(shí)驗(yàn)發(fā)現(xiàn)，臨沂地區(qū)利用多項(xiàng)式擬合曲面模型優(yōu)于插值曲面模型，插值和擬合精度是否能通過(guò)節(jié)點(diǎn)的增加而進(jìn)一步提高，將在下一步的工作中逐步完善。

擬合曲面模型減少了大面積測(cè)量時(shí)GPS數(shù)據(jù)處理軟件無(wú)法控制全區(qū)域而將項(xiàng)目分塊的麻煩，通過(guò)多項(xiàng)式擬合方法得到的模型是唯一最優(yōu)解，多項(xiàng)式系數(shù)是唯一確定的，這樣就可以將高程異常的求解轉(zhuǎn)化為平面坐標(biāo)和高程的函數(shù)關(guān)系。每次外業(yè)作業(yè)結(jié)束后可通過(guò)函數(shù)關(guān)系將GPS采集的 2000國(guó)家大地坐標(biāo)和大地高轉(zhuǎn)化為2000國(guó)家大地坐標(biāo)和正常高，同時(shí)保證了整個(gè)臨沂市高程異常數(shù)據(jù)的完整性，為臨沂市整個(gè)區(qū)域范圍內(nèi)的大面積測(cè)量提供了方便。