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

周 偉

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

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

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

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

1.1 插值模型

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

1.1.1 多項式插值曲面

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

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

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

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

1.1.2 三次樣條函數插值曲面

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

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

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

1.2 多項式擬合曲面模型

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

2 實驗與結果分析

2.1 數據準備

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

2.2 數據實驗

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

圖1 線性插值曲面

圖2 三次多項式插值曲面

圖3 三次B樣條插值曲面

圖4 線性擬合曲面

圖5 二次多項式擬合曲面

圖6 三次多項式擬合曲面

圖7 四次多項式擬合曲面

圖8 五次多項式擬合曲面

2.3 數據分析

2.3.1 定性分析

1）通過圖2、3可直觀發現，多項式插值方法只對已知節點進行插值，而樣條插值方法和擬合方法，在東坐標和北坐標的區間范圍內均可估算高程異常曲面。

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

3）多項式擬合曲面時，次數大于3時，在沒有點控制的邊緣區域將產生較大的形變。

2.3.2 定量分析

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

表1 多項式擬合曲面中誤差統計/m

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

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

表2 插值和擬合曲面求解高程異常殘差統計

3 結 語

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

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