插值和擬合方法在高程異常模型建立中的應(yīng)用

周 偉

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

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

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

1 高程異常模型的計算方法

1.1 插值模型

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

1.1.1 多項式插值曲面

1）線性插值曲面。該方法比較簡單，計算公式為：

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

2）高次多項式插值曲面。對于n＞1的情形，其計算公式為：

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

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

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

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

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

1.2 多項式擬合曲面模型

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

2 實驗與結(jié)果分析

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

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

2.2 數(shù)據(jù)實驗

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

圖1 線性插值曲面

圖2 三次多項式插值曲面

圖3 三次B樣條插值曲面

圖4 線性擬合曲面

圖5 二次多項式擬合曲面

圖6 三次多項式擬合曲面

圖7 四次多項式擬合曲面

圖8 五次多項式擬合曲面

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

2.3.1 定性分析

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

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

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

2.3.2 定量分析

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

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

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

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

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

3 結(jié) 語

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

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