基于Python的GPS高程異常擬合研究

于艷超，許捍衛，張明希

（1. 河海大學 地球科學與工程學院，江蘇 南京 210098）

基于Python的GPS高程異常擬合研究

于艷超1，許捍衛1，張明希1

（1. 河海大學 地球科學與工程學院，江蘇 南京 210098）

考慮GPS高程異常具有一定的空間相關性，可將ArcGIS地統計中的確定性插值方法應用到GPS高程異常擬合中。結合Python腳本實現基于加權總體最小二乘的GPS高程異常擬合算法。結果表明，該算法的點位擬合精度和檢核精度均優于ArcGIS地統計確定性插值。

GPS高程異常；地統計分析；總體最小二乘；Python

擬合方法是工程領域進行GPS高程轉換的首選方案，其基本原理是在中小測區將似大地水準面看作平面或曲面，用多項式函數擬合法進行高程轉換[1-3]。考慮到GPS三維坐標中平面坐標與高程以及GPS測量所得大地高與水準測量所得正常高的不等精度性，通過對觀測向量以及系數矩陣定權，采用加權總體最小二乘（weighted total least-squares，WTLS）的GPS高程曲面擬合算法，無論數學的嚴密性，還是模型參數的精確性以及空間數據的適用性，均優于ArcGIS地統計中的確定性插值方法[4]。

Python作為一種通用腳本語言，可以像Matlab一樣隨寫隨運行[5]。從ArcGIS10．0開始，ESRI在所有產品中都集成了基于Python的站點包ArcPy。作為ArcGIS系統擴展和二次開發的語言，其靈活性和重用性可以很方便地實現空間數據的處理[6]。ArcPy提供了大量的類和函數，可以輕松執行ArcGIS工具，并創建原生對象，如幾何、柵格、空間參考等[7]。本文基于Python，實現了基于加權總體最小二乘的GPS高程異常擬合算法。

1 加權總體最小二乘趨勢面擬合算法

1.1 基于最小二乘的GPS高程異常擬合算法

GPS高程異常擬合的基本原理是在中小測區將似大地水準面看作一個連續的平面或曲面。本文采用二次多項式曲面擬合：

式中，a0～a5為待求參數。利用線性回歸模型進行參數估計的Guass-Markov 模型為:

式中，

為擬合點數。

根據最小二乘估計準則，得最小二乘解為：

單位權中誤差為：

1.2 加權總體最小二乘算法

從最小二乘算法模型與解算過程不難看出，最小二乘法只考慮了高程異常的誤差，而未考慮系數矩陣的誤差。然而，從組成系數矩陣A的元素可見，除了第一列，其他列均是與大地高坐標有關聯的，忽略系數矩陣的誤差不僅不符合實際，而且從理論上講也是不夠嚴密的。系數矩陣A中含有誤差的模型[8]為:

隨機模型為：

GPS測量所得的三維坐標中平面坐標的精度要高于高程，且高程異常值為大地高與正常高之間的差值，這就導致平面點位X和Y的觀測精度相當，而高程異常ζ的精度則不同。定義每個觀測點之間相互獨立且具有相同的平面測量精度，則σx=σy，系數矩陣A的列向量權陣表示系數矩陣中每列觀測值之前的權比關系，也可以表示為不同類觀測值之間的權比關系。所以，根據協因素傳播律以及系數矩陣的數值形式，可得：

2 算例分析

本文選取某城市四等水準聯測的58組城市E級GPS控制網高程測量數據為數據源，根據已知數據的點位坐標和大地高、正常高的測量值，計算高程異常。數據點分布在21．5 km×27．5 km的范圍內，點位密度滿足GPS高程異常曲面擬合的要求。高程異常值從0．656到1．769，變化范圍較小。在進行曲面擬合時，為了檢驗模型的可靠性，根據數據點的分布情況，從中選取分布較為分散的10個點作為檢核點，其余48個作為模型擬合的數據點。由于測點的平面坐標數值較大，為了避免在運算時出現病態矩陣的情況，需要對測點平面坐標數據進行中心化處理。依據這48個點建立誤差方程，求解模型參數并計算檢核點中誤差。技術路線如圖1。

圖1 技術路線

首先利用QQPlot分布圖和泰森多邊形對數據進行探索性分析。從QQPlot分布圖（圖2）可以看出，樣本分位數作出的圖形基本在一條直線上。因此，該樣本數據的分布近似服從正態分布。粗差點在泰森多邊形中表現為眾多多邊形中顏色與周圍截然不同的點。從圖3可以看出，沒有與周圍顏色不一樣的多邊形，說明泰森多邊形代表的樣本數據中不含粗差點。

圖2 QQPlot分布

圖3 泰森多邊形

空間趨勢面分析主要依靠樣本數據來擬合一個曲面，從而大致反映其空間分布變化情況。圖4表明，樣本數據在ZX方向和ZY方向呈現二次擬合曲線。

圖4 空間趨勢面分析

最后，采用地統計分析模塊中的插值工具對樣本數據進行實驗。本實驗中采用全局多項式插值、反距離權重和局部多項式插值。每種方法中，參數選取不同會對內插精度產生影響，因此選取的原則是使樣本數據通過轉換盡可能服從不同算法的前提假設，并使樣本數據的空間相關性達到最大，從而使擬合的結果均方根最小。本文開發的加權總體最小二乘趨勢面擬合算法工具界面如圖5所示。

圖5 加權總體最小二乘趨勢面擬合算法工具界面

地統計分析插值與加權總體最小二乘擬合算法精度對比，如表1所示，其中lb2、gb2、IDW和WTLS分別代表局部多項式、全局多項式、反距離權重和加權總體最小二乘。

表1 地統計分析插值與加權總體最小二乘擬合算法精度對比

從表1可以看出，WTLS方法的檢核中誤差較lb2、gb2和IDW有所改善，WTLS方法擬合的精度相比于lb2、gb2和IDW分別提高了26%、5．1%和63%。由此可見，加權總體最小二乘這種顧及系數陣誤差和觀測值誤差，并對不同觀測值賦予不同的權值參與平差計算的方法更合理、更精確。

3 結 語

1）當系數陣含有觀測隨機誤差時，應該使用總體最小二乘方法進行參數的求解，而當系數陣列向量之間具有相關性或觀測值之間不等權時，就要考慮采用加權總體最小二乘方法。在應用加權總體最小二乘方法時要注意定權方法，確保給定權值的合理性，這樣得出的結果才是最優的。

2）以GPS高程異常擬合為例，引入地統計分析工具對空間數據分析和處理，剔除粗差，進行趨勢面分析。將確定性插值法的計算結果與加權總體最小二乘算法處理效果比較，為ArcGIS地統計模塊的算法設計提供參考。結合開源的Python工具包（NumPy、SciPy），開發基于加權總體最小二乘的趨勢面擬合算法，以Python腳本的形式集成到ArcToolbox中，并通過ArcGIS Server發布，腳本簡單易用，減少了處理海量數據時人機交互的工作量。

3）曲面擬合是三維建模逆向工程的重要組成部分，主要研究利用空間散亂數據點來重構符合要求的曲面。當空間散亂點具有空間相關性時，就可以采用基于ArcGIS地統計工具的空間插值方法。但是當測量數據精度達不到要求時，擬合后的曲面誤差較大。用擬合效果更好的基于加權總體最小二乘的二次曲面建模方法，可以得到更優的曲面擬合精度和更好的未知空間點預測值。

[1] 魏立峰,何建國． GPS高程擬合似大地水準面的方法[J]．地理空間信息,2010,8(4): 72-73

[2] 曹先革,劉金鵬．基于MATLAB的GPS高程擬合程序設計[J]．地理空間信息,2009,7(2): 22-24

[3] 王小輝,王琪潔,丁元蘭，等． 基于二次曲面和BP神經網絡組合模型的GPS高程異常擬合[J]．大地測量與地球動力學,2012,32(6): 103-105

[4] 王樂洋,許才軍．總體最小二乘研究進展[J]．武漢大學學報：信息科學版,2013,38(7): 850-856

[5] 張若愚．Python科學計算[M]．北京：清華大學出版社，2012

[6] 彭海波,向洪普．基于Python的空間數據批量處理方法[J]．測繪與空間地理信息,2011,34(4):81-82

[7] Dobesova Z． Programming Language Python for Data Processing[C]．ICECE 2011, 2011

[8] 魯鐵定,周世健．總體最小二乘的迭代解法[J]．武漢大學學報：信息科學版，2010,35(11): 1 351-1 354

P216

B

1672-4623（2014）04-0049-03

10.11709/j.issn.1672-4623.2014.04.017

于艷超，碩士，主要研究方向為地理信息系統開發與應用。

2013-10-25。

項目來源：水利部公益性行業科研專項經費資助項目（201201025）；國家自然科學基金資助項目（41101374、41101308）。