于艷超，許捍衛(wèi)，張明希

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

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

于艷超1，許捍衛(wèi)1，張明希1

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

考慮GPS高程異常具有一定的空間相關(guān)性，可將ArcGIS地統(tǒng)計中的確定性插值方法應(yīng)用到GPS高程異常擬合中。結(jié)合Python腳本實現(xiàn)基于加權(quán)總體最小二乘的GPS高程異常擬合算法。結(jié)果表明，該算法的點位擬合精度和檢核精度均優(yōu)于ArcGIS地統(tǒng)計確定性插值。

GPS高程異常；地統(tǒng)計分析；總體最小二乘；Python

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

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

1 加權(quán)總體最小二乘趨勢面擬合算法

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

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

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

式中，

為擬合點數(shù)。

根據(jù)最小二乘估計準則，得最小二乘解為：

單位權(quán)中誤差為：

1.2 加權(quán)總體最小二乘算法

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

隨機模型為：

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

2 算例分析

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

圖1 技術(shù)路線

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

圖2 QQPlot分布

圖3 泰森多邊形

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

圖4 空間趨勢面分析

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

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

地統(tǒng)計分析插值與加權(quán)總體最小二乘擬合算法精度對比，如表1所示，其中l(wèi)b2、gb2、IDW和WTLS分別代表局部多項式、全局多項式、反距離權(quán)重和加權(quán)總體最小二乘。

表1 地統(tǒng)計分析插值與加權(quán)總體最小二乘擬合算法精度對比

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

3 結(jié) 語

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

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

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

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

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

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

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

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

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

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

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

P216

B

1672-4623（2014）04-0049-03

10.11709/j.issn.1672-4623.2014.04.017

于艷超，碩士，主要研究方向為地理信息系統(tǒng)開發(fā)與應(yīng)用。

2013-10-25。

項目來源：水利部公益性行業(yè)科研專項經(jīng)費資助項目（201201025）；國家自然科學(xué)基金資助項目（41101374、41101308）。