樣條函數擬合GPS高程探討①

張海艷，張振濤，蘇貴波，李清生

（1.國家測繪局第一地形測量隊，陜西 西安710054；2.中國人民解放軍61491部隊，河南 許昌461131）

0 引 言

全球定位系統（GPS）是隨著現代科學技術的迅速發展，而建立起來的新一代精密定位系統［1－2]。由于該系統不受氣象條件的影響、自動化程度高，且定位精度良好，已經廣泛應用于各個領域，水平定位精度在短時間內即可達到±5mm＋1ppm［3].目前，GPS測高的能力得到了普遍的關注，研究GPS測定水準高程的方法將有助于更好的發揮GPS的作用。為了使GPS測量成果能夠在解決高程問題方面得到較好的應用，綜合利用水準測量和重力測量資料，確定大地水準面高程，綜合利用大地水準面資料確定點的高程。

1 數學模型

解決GPS測定高程問題的關鍵是如何確定高程異常差［4]。

1.1 三次樣條函數曲線擬合法

樣條數學含義為：曲線直到二階導數連續，沒有多余的拐點，曲率變化比較均勻［5]。當測線較長時，節點處的高程異常值變化會較大，當插值點過多時，若采用多項式曲線擬合法，會出現不光滑、不連續的現象，影響插值精度。但是，若采用三次樣條函數曲線擬合的方法就會彌補上述不足［6－7]，使插值精度得到保證。三次樣條函數是光滑連續的，可以更好的處理較長曲線的擬合問題［7]。

節點處的函數值和二階導數或函數值和一階導數來建立三次樣條插值函數的表達式及連續性方程和邊界條件給定及其算法。

記S（x）在節點xi處的函數值，一階導數和二階導數為

在每個小區間［xi－1，xi]上，S（x）都是不超過3次的多項式，則S′（x）在小區間上是不超過1次的多項式，那么S＂（x）在兩節點處的值應為常數，即y"＝0.此時，稱為自然插值三次樣條函數。假定節點的二階導數為零，進行擬合運算

1.1.1 邊界條件

由式（6）決定的M 連續性方程，都是n＋1個未知數的n－1個線性代數方程，顯然，僅有這些方程組是無解的，要想唯一確定方程組的解，必須再附加兩個方程，這里給出一種附加條件的方法，通常是在區間節點處給出，稱為邊界條件。

1.1.2 算法

為了求出三次樣條函數s（x），關鍵在于求解出式（4），在三對角系數矩陣中，（i＝1，2，…，n） ，且0≤λ0，μ0，λn，μn≤1主對角線元素等于2，對角嚴格占優，因為，這樣方程組的解存在并且唯一。為了節省計算時間和存儲量，一般都采用追趕法求解［8－10]。

1.2 多項式曲線擬合法［4，6，11]

設ξi和xi（或yi）存在的函數關系（i＝0，1，2，…，n）可以用以下m（m≤n）項多項式擬合

在節點處的離差Ri＝ξi（x）－ξi平方和最小的條件下：

求解出式中各系數，即可按ξm（x）＝a0＋a1x＋a2x2＋…＋amxn求出測線方向上（或方向線左右任一點的）ξ值。

2 實例驗算及分析比較

2.1 東北某地GPS測線的擬合概算

為了檢驗曲線擬合法在長距離長間距的擬合效果，我們選了東北某地由六個點組成的一條測線，如圖1所示。

測線總長474km，測線上各點都聯測了四等幾何水準，分別用三次樣條函數和多項式曲線擬合法進行了擬合，所得中誤差如表1所示。

比較樣條函數和多項式擬合的中誤差，可以看出，運用樣條函數的擬合效果比多項式擬合的精度高。

圖1 東北某地GPS測線

表1 東北某地樣條函數和多項式擬合中誤差對比

2.2 王灘GPS網的擬合概算

王灘試驗GPS網位于唐山市樂亭縣王灘鎮以南，瀕臨渤海灣，測區面積約170km2，全網共31個點，最長邊為23.5km，最短邊為956m，除30和31號點為已知坐標點遠離測區未作水準聯測外，其余點都實測了四等水準，如圖2所示。

圖2 王灘試驗GPS網

在擬合計算時，僅選了11－12－7這條測線上的點，以及其左右兩側的點，運用三次樣條函數及多項式曲線擬合法進行擬合，如表2、表3所示。

2.2.1 節點均勻方案

選取測線上分布較均勻的11、14、12、4、7等點進行擬合計算，得出各點的中誤差如表2所示。

表2 節點均勻中誤差對比

2.2.2 節點非均勻方案

選取測線上分布不均勻的11、15、12、3、7等點進行擬合計算，得出各點的中誤差如表3所示。

表3 節點不均勻中誤差對比

從表中可以看出，樣條函數擬合的精度高于多項式擬合的精度，而且均勻方案中點的擬合精度普遍高于不均勻方案點的擬合精度。

2.3 廈門GPS網擬合概算

廈門GPS網共113個點，為Ⅱ、Ⅲ等全面網，見圖3，地處東經117°53′～118°23′，北緯24°24′～24°54′，位于閩南金三角，面積約1 500km2.測區地勢為西北高、東南低，且瀕臨臺灣海峽。西北部最高點海拔564m，東南最低點海拔5m，平均海拔93m.

圖3 廈門GPS網

在驗證計算中，選用了測線67～80～92，并且以67、66、65、71、76、77、80、92為節點，以69、70、72、75、78、81、84、89、87、90、91為插值點，選用四種擬合方案進行擬合計算，具體結果如表4、5、6、7所示。

2.3.1 插值點相同而節點不同（均勻方案取70、81、87、91為插值點，67、71、80、92為節點）（如表4所示）

表4 插值點相同而節點不同中誤差對比

2.3.2 插值點相同而節點不同（非均勻方案取70、81、87、91為插值點，67、66、77、80為節點）（如表5所示）

表5插值點相同而節點不同中誤差對比

2.3.3 節點同而插值點不相同（均勻方案取70、81、87、91為插值點，67、80、92為節點）（如表6所示）

表6 節點同而插值點不相同中誤差對比

2.3.4 節點同而插值點不相同（非均勻方案取69、75、84、91為插值點，67、80、92為節點）（如表7所示）

表7 節點同而插值點不相同中誤差對比

當插值點相同而節點不同，其擬合效果將隨著節點選取的均勻程度提高而提高；而當節點相同，插值點不同時，插值點的擬合效果也隨其均勻程度的提高而提高。而且，當插值點距節點越近時其擬合效果越好，即插值點的擬合效果受節點的影響較大。因此，節點精度的高低對插值點的擬合精度有直接影響。

2.4 方案分析比較

2.4.1 已知點位分布規律

1）控制點（即節點）數相同，點位分布不同的比較

三次樣條函數曲線擬合時，非均勻分布的方案：在2.3.4節中，點69、75、84、91的中誤差分別為；m樣＝±30.5cm、±21.2cm、±8.7cm、±9.3cm，m多＝±31.5cm、±23.3cm、±18.5cm、±24.6cm.而均勻方案2.3.1和2.3.3節中，m樣＝±7.3cm、±8.3cm，m多＝±7.7cm、±8.4cm.在2.2節中，總體擬合效果都較好，最小的樣條擬合中誤差達±1.6cm，最小多項式擬合中誤差為 ±4.4cm；而最大的樣條擬合中誤差為±2.4cm，最大多項式擬合中誤差達到了±23.3cm，由此可以看出樣條函數擬合的效果要比多項式擬合的效果好得多。再由2.2節中第二種方案可看出，當節點與插值點都均勻時，其擬合的效果都較好；在2.1節中，由于測線較長，其 擬 合 效 果 不 太 好，m樣＝ ±41.4cm，m多＝±42.3cm.可見，已知控制點分布均勻，各種擬合法的精度都高于非均勻分布的精度。

2）擬合方法相同，已知點數不同的情況

擬合點的選擇對GPS水準高程精度的影響較大，一般應選擇一定的擬合點且均勻分布于測區。在2.3節中，可以看出，擬合點的增加對精度影響并不顯著，用適當的擬合點也能獲得滿意的結果。可見，增加已知點數來提高擬合精度是有限的。

2.4.2 擬合方法優劣的比較

由2.2節中二種方案，可以看出，在均勻分布情況下，兩種擬合法的精度基本相當，約為（0.2～0.3）±3cm.但一般情況下，三次樣條函數曲線擬合要優于多項式擬合。

3 結 論

1）在GPS高程擬合中，對于局部小區域內的線狀GPS網，在平原和丘陵地區，兩種方法擬合的精度都可以達到±5cm以內，而三次樣條函數曲線擬合法更適合于線狀GPS網的測區。

2）兩種擬合法，采用已知控制點的個數要適當，以略多于待定系數個數為宜，一般為網中總點數的20%～30%，點的布設以均勻分布最佳，且不宜外插擬合。

3）測線過長時，擬合精度不高，應加密GPS點，充分利用重力水準資源，選點要均勻。

4）獲取GPS數據宜采用雙頻接收機，消弱電離層影響，選擇最佳衛星幾何圖形，選擇衛星位近和觀測時段，以獲取精確的數據。

5）由于插值點距離節點越近，則節點對插值點高程影響越大，所以要盡量提高已知點的精度。

［1]劉基余，李征航.全球定位系統原理與應用［M].北京：測繪出版社，1993.

［2]周忠謨，易杰軍，周 琪.GPS衛星測量原理與應用［M].北京：測繪出版社，1997.

［3]吳俐民.GPS測定高程的方法研究［J].城市勘測，1997（4）：17－21.

［4]于來法.用數值擬合法確定GPS正常高的研究［J].解放軍測繪學報，1996，6（2）：90－94.

［5]蘇步青，劉鼎元.計算幾何［M].上海：上海科學技術出版社，1981.

［6]徐紹銓，李征航.擬合法求定GPS點的正常高［J].武測科技，1992（2）：15－21.

［7]陸付民.樣條函數在垂直變形分析中的應用［J].勘察科學技術，1990（5）：53－55.

［8]諸梅芳，耿美英，薛 毅，等.實習與計算方法［M].北京：北京工業大學出版社，1996.

［9]蔣爾雄，趙風光.數值逼近［M].上海：復旦大學出版社，1996.

［10]黃友謙.曲線曲面的數值表示和逼近［M].上海：上海科學技術出版社，1984.

［11]陶本藻.GPS水準似大地水準面擬合和正常高計算［J].測繪通報，1992（4）：14－18.