于艷超，許捍衛(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)計(jì)中的確定性插值方法應(yīng)用到GPS高程異常擬合中。結(jié)合Python腳本實(shí)現(xiàn)基于加權(quán)總體最小二乘的GPS高程異常擬合算法。結(jié)果表明，該算法的點(diǎn)位擬合精度和檢核精度均優(yōu)于ArcGIS地統(tǒng)計(jì)確定性插值。

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

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

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

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

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

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

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

式中，

為擬合點(diǎn)數(shù)。

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

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

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

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

隨機(jī)模型為：

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

2 算例分析

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

圖1 技術(shù)路線

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

圖2 QQPlot分布

圖3 泰森多邊形

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

圖4 空間趨勢(shì)面分析

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

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

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

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

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

3 結(jié) 語

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

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

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

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

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

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

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

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

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

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

[8] 魯鐵定,周世?。傮w最小二乘的迭代解法[J]．武漢大學(xué)學(xué)報(bào)：信息科學(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

于艷超，碩士，主要研究方向?yàn)榈乩硇畔⑾到y(tǒng)開發(fā)與應(yīng)用。

2013-10-25。

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