基于地統(tǒng)計學(xué)的GPS高程異常擬合研究①

周長江，紀志鵬

（中國礦業(yè)大學(xué) 環(huán)境與測繪學(xué)院，江蘇 徐州221008）

0 引 言

GPS測量具有精度高、實時、快速等優(yōu)點，廣泛用于工程測量、變形監(jiān)測等領(lǐng)域。GPS在平面定位方面已經(jīng)可以達到很高的精度，而GPS高程控制網(wǎng)精度比較低。GPS高程是以WGS－84為參考橢球的大地高，而實際需要的是以似大地水準面為基準的正常高，這兩者之間的差值叫做高程異常，若能求出高精度的GPS點的高程異常ζ就可以把大地高轉(zhuǎn)換到正常高。近年來，許多學(xué)者圍繞這一目標提出了很多方法，如多項式擬合、BP神經(jīng)網(wǎng)絡(luò)算法［1]。但多項式擬合算法誤差較大，BP算法受初始權(quán)值的影響較大，不夠穩(wěn)定。近年來，ArcGIS的應(yīng)用范圍越來越廣，它具有強大的空間數(shù)據(jù)管理、分析、處理能力。ArcGIS中的Geostatistical Analysis Tools是基于地統(tǒng)計學(xué)的原理對數(shù)據(jù)進行處理和分析的工具［2]。用到的方法主要是空間內(nèi)插方法，按其實現(xiàn)的數(shù)學(xué)原理分為：確定性插值、克里金插值。確定性插值方法有反距離加權(quán)值法、全局性插值法、局部性插值法、徑向基插值法?？死锝鸩逯涤址譃槠胀死锝?、簡單克里金和泛克里金等。結(jié)合實例將地統(tǒng)計分析工具中的相關(guān)插值方法運用到GPS點的高程異常擬合中，并對它們擬合的精度進行分析比較。

1 確定性插值

1.1 反距離加權(quán)值法

反距離加權(quán)值法的基本原理就在于空間物體的距離越小，它們的性質(zhì)就越相似；距離越大則相似性就越小。它是以插值點與樣本點的距離為權(quán)重進行加權(quán)平均，離插值點越近的樣本點賦予的權(quán)重越大。用反距離加權(quán)值法進行高程異常擬合，其計算公式為

式中：ζi為相鄰點的高程異常；pi為數(shù)據(jù)ζi的權(quán)值。

1.2 全局多項式插值法

全局性插值法以整個測區(qū)的高程異常樣本數(shù)據(jù)集為基礎(chǔ)，用一個數(shù)學(xué)多項式來模擬計算預(yù)測值，可以視為用一個多項式曲面或是平面來進行全區(qū)域的擬合。但是，利用全局插值法生成的表面易受極高和極低樣點值的影響，擬合的曲面很少能夠與已知樣點完全重合，因此，它是一種非精確的方法，適用于地勢變化平緩的測區(qū)。

1.3 局部多項式插值

局部多項式采用多個多項式，每個多項式都處在特定的重疊鄰近區(qū)域內(nèi)，可以生成一個比較平滑的曲面，但它也是一個非精確的方法，適用于建立平滑表面和小范圍的變異。

1.4 徑向基函數(shù)插值

徑向基函數(shù)插值方法適用于樣本數(shù)據(jù)比較多，同時要求獲得平滑的曲面情況。對于變化緩慢的曲面，如求其高程異常，徑向基函數(shù)插值法能夠得到理想的結(jié)果，但是對于短距離內(nèi)，變化比較大的曲面，采用徑向基函數(shù)誤差就會很大。

2 克里金插值法

克里金插值法又稱空間局部插值法，是以變異函數(shù)和結(jié)構(gòu)分析為基礎(chǔ)，對有限區(qū)域內(nèi)區(qū)域變化量進行無偏最優(yōu)估計的一種方法。克里金方法又分為普通克里金、簡單克里金、泛克里金等。

2.1 普通克里金估計

普通克里金估計是一種內(nèi)蘊假設(shè)（或二階平穩(wěn)假設(shè)）下期望未知的區(qū)域變化量估值方法［3]，它是一種建立在樣本數(shù)據(jù)變化呈正態(tài)分布的前提下，運用普通克里金插值進行高程異常擬合，其估值公式為

其中：ζ＊（x0，y0）是在待估位置P（x0，y0）處的高程異常估值；ζ（xi，yi）是某個位置的已知高程異常值；λi為分配給ζ（xi，yi）的權(quán)重；n是樣本個數(shù)。普通克里金方法的最優(yōu)估計條件為估計誤差方差最小，其無偏估計條件為

2.2 泛克里金估計方法

泛克里金估計假設(shè)數(shù)據(jù)中存在主導(dǎo)趨勢，且該趨勢可以用一個多項式擬合或是確定的函數(shù)擬合，在進行克里金趨勢分析時，分析數(shù)據(jù)中存在的變化趨勢，獲得擬合模型，對殘差數(shù)據(jù)進行克里金分析，將趨勢面和殘差面分析的結(jié)果相加得到最終結(jié)果。泛克里金估計方法的估計公式為

式中：ζ＊（x0，y0）是在待估位置P（x0，y0）處的高程異常估值；ζ（xi，yi）是待估點附近某個位置的已知高程異常值；λi為分配給ζ（xi，yi）的權(quán)重；n是樣本個數(shù)。泛克里金方法的最優(yōu)估計條件為使估計誤差方差最小，其無偏估計條件

式中，fl（x，y）為P（x，y）的一次或二次函數(shù)。

2.3 基于克里金方法的GPS高程異常擬合步驟

采用克里金插值求GPS點的高程異常主要分為三個步驟：空間數(shù)據(jù)探索性分析、空間樣本點的結(jié)構(gòu)化量化分析、對未知點的高程異常進行預(yù)測。由于克里金方法建立在一定的假設(shè)基礎(chǔ)上，其在一定程度上要求其數(shù)據(jù)具有相同的變異性。簡單克里金要求數(shù)據(jù)服從正態(tài)分布，泛克里金插值假設(shè)數(shù)據(jù)存在主導(dǎo)趨勢。因此，首先要進行空間數(shù)據(jù)探索性分析來了解樣本數(shù)據(jù)的分布情況?？臻g樣本的結(jié)構(gòu)化分析指對樣本數(shù)據(jù)擬合一個空間獨立模型，進而可以用擬合的模型由未知點的平面位置對區(qū)域的GPS點進行高程異常預(yù)測。

3 空間探索性分析

3.1 數(shù)據(jù)分布檢驗

由于普通克里金建立在已知樣本數(shù)據(jù)服從正態(tài)分布的前提下，需要對數(shù)據(jù)分布進行分析。Arc－GIS地統(tǒng)計分析模塊的探索性空間數(shù)據(jù)分析（ESDA）工具允許用戶用多種方法對樣本數(shù)據(jù)進行分析。ESDA提供的方法有直方圖法（hisgogram）、正態(tài) QQPlot分布圖（voronoi map、normal QQPlot）、正交協(xié)方 差云（semivariogram／covariance cloud）等，可以用來查明數(shù)據(jù)分布、尋找局部和全局離群值、探查全局趨勢、檢測數(shù)據(jù)空間自相關(guān)以及數(shù)據(jù)間的協(xié)變。QQPlot分布圖是可以將樣本數(shù)據(jù)的分布與標準正態(tài)分布對比，利用樣本數(shù)據(jù)分布的分位數(shù)作出的圖形是否接近一條直線，從而來分析和評價現(xiàn)有數(shù)據(jù)是否近似服從正態(tài)分布。從圖1可以看出樣本分布分位數(shù)作出的圖形基本上在一條直線上，因此，樣本數(shù)據(jù)分布近似服從正態(tài)分布。

圖1 樣本數(shù)據(jù)的正態(tài)QQPlot分布圖

3.2 趨勢分析

泛克里金插值方法假設(shè)樣本數(shù)據(jù)中存在主導(dǎo)趨勢，地物的空間趨勢反映了物體在空間區(qū)域上變化的主體特征?？臻g趨勢面分析主要依靠樣本數(shù)據(jù)來擬合一個曲面，從而大致反映其空間分布變化情況。圖2反映了實驗樣本的空間變化情況。圖中豎軸Z代表高程異常值，X、Y代表點的平面位置。數(shù)據(jù)點被投影到兩個正交的面上，通過投影點可以做出一條最佳擬合曲線，用來模擬特定方向上的趨勢，如果擬合的線是直的說明沒有趨勢變化。

圖2 趨勢分析

3.3 方差變異分析

地理空間自相關(guān)是指時間序列相鄰數(shù)值間的相關(guān)關(guān)系，大部分的地理現(xiàn)象都具有空間相關(guān)特性，即距離越近的事物越相似。半變異／協(xié)方差函數(shù)云圖就是這種相似性的定量表示。其值越小，則就越相似?？臻g相關(guān)性僅與距離有關(guān)時稱為各向同性，但在實際應(yīng)用中各向異性更為普遍，即當考慮方向影響時，在某個方向距離更遠的事物確有更好的相關(guān)性。在實際操作中通過調(diào)節(jié)show search direction的角度值使半變異／協(xié)方差函數(shù)云的值最小 如圖3、圖4所示，在大約68°方向時半變異／協(xié)方差函數(shù)值趨近于零，此時，樣本數(shù)據(jù)間的高程異常之間的相關(guān)性最大。如果能夠探測出自相關(guān)中的方向效應(yīng)，就可以在擬合模型中考慮這個因素，從而提高擬合的準確性。

圖3 半變異函數(shù)云圖

圖4 協(xié)方差函數(shù)云圖

3.4 粗差剔除

已知樣本數(shù)據(jù)的好壞對后續(xù)的擬合精度有很大影響，因此，在擬合前要對數(shù)據(jù)做粗差檢驗，并剔除粗差點。在ArcGIS地統(tǒng)計分析工具中，粗差點的檢驗可以利用Voronoi圖、直方圖、半變異／協(xié)方差函數(shù)云三種工具進行檢驗。粗差點在直方圖中體現(xiàn)為一些孤立存在或被一些顯著不同的點值包圍的點。在Voronoi圖中則表現(xiàn)為眾多邊形中顏色與周圍顏色截然不同的點。圖5反映了采用Voronoi圖方法對樣本數(shù)據(jù)中粗差點的檢驗情況。從圖中可以看出一個綠色的多邊形與周圍多邊形的顏色均不一樣，說明其代表的樣本點有可能是含有粗差的點。當此樣本點在直方圖中表現(xiàn)為孤立存在的點，并且半變異／協(xié)方差函數(shù)云值不趨近于零偏離較遠時 ，就有理由相信此樣本點是粗差點，應(yīng)予剔除。

圖5 Voronoi Map

4 實例計算和精度比較

4.1 實例計算

實驗數(shù)據(jù)選自南方某地區(qū)聯(lián)測了水準的51個點，這些點的正常高、大地高和高程異常已知。把這些已知數(shù)據(jù)分為兩組，一組34個點作為已知樣本數(shù)據(jù)進行空間結(jié)構(gòu)化量化分析，另一組17個點作為檢驗點，用來檢驗擬合精度情況。GPS點點位分布情況和高程異常值如圖6所示。

利用ArcGIS中的地統(tǒng)計分析模塊中的插值工具，分別采用上述方法對樣本數(shù)據(jù)進行試驗。由于每種方法中參數(shù)選取的不同會對內(nèi)插的精度產(chǎn)生影響，選取的原則就是使樣本數(shù)據(jù)通過轉(zhuǎn)換盡可能服從不同算法的前提假設(shè)，并使樣本數(shù)據(jù)的空間相關(guān)性達到最大，從而使擬合的結(jié)果均方根最小。為了便于比較，將每種插值方法中所有涉及到的參數(shù)進行組合分別算出結(jié)果，選取均方根最小的作為該種算法的最終結(jié)果進行比較。反復(fù)試驗得出當反距離加權(quán)值、全局插值方法中的power參數(shù)取為3次時均方根最小。局部插值法中power參數(shù)取為2時均方差最小。徑向基方法中，當kernel function選為Thin plate spline時均方差最小。普通克里金中Transformation中選為Box－cox，Order of trend removal選為Const，Model選為Gaussan時均方根最小。泛克里金中Transformation選為 Box－cox，Order of trend removal選為Const，Model選為 Gaussan，Standard中的neighbors to include選為7時均方根最小。不同算法擬合結(jié)果如表1所示（高程異常單位為m，誤差單位為cm，限于篇幅所限，局部插值法數(shù)據(jù)沒有列出）。

圖6 已知數(shù)據(jù)點位分布和高程異常

表1 不同方法的結(jié)果比較

圖7 高程異常預(yù)測曲面

擬合高程異常求出后，生成高程異常的預(yù)測曲面，圖7是泛克里金法生成的高程異常預(yù)測曲面。

4.2 精度比較

把均方根（RMS）誤差作為衡量算法精度的指標［4]，不同算法的均方根誤差如表2所示。

表2 不同算法的均方根誤差

從表1、2可以看出反距離插值算法的精度最低，其誤差最大值為9.9cm，平均誤差為3.59cm.全局插值法誤差最大為4.9cm，最小為0，平均為1.97cm.徑向基最大誤差為4.4cm，平均誤差為1.4cm.普通克里格最大誤差為4.2cm，平均誤差為1.524cm.泛克里格法最大誤差為4.2cm，平均誤差為1.5cm.從實驗結(jié)果可以看出徑向基法比較穩(wěn)定，且誤差比較小，精度比較高，克里金法精度略低于徑向基法。這主要是由于克里金方法對原始觀測數(shù)據(jù)進行了平滑處理，導(dǎo)致克里金估值的結(jié)果的空間結(jié)構(gòu)在整體上擬合原始觀測數(shù)據(jù)不如徑向基法精度高。此外，本實驗采用樣本數(shù)據(jù)點比較多，且測區(qū)內(nèi)高程異常曲面變化相對比較緩慢適合采用徑向基方法。

5 結(jié) 論

ArcGIS中的地統(tǒng)計分析工具具有強大的空間數(shù)據(jù)分析處理能力，用地統(tǒng)計法進行GPS高程異常擬合優(yōu)點在于不用考慮GPS點的坐標的基準，可以由任何基準的平面坐標得到高程異常，而且還能夠達到較高的精度。在實際應(yīng)用中要注意不同算法參數(shù)選取對結(jié)果產(chǎn)生很大的影響。此外，地統(tǒng)計分析方法應(yīng)用于GPS高程異常擬合對于不同的樣本和樣本容量，這些算法的精度高低不同，因此，需具體問題具體分析，常選擇合適的算法。

［1]楊明清，靳 蕃.用神經(jīng)網(wǎng)絡(luò)方法轉(zhuǎn)換高程［J].測繪學(xué)報，1999，28（4）：301－307.

［2]池 建.精通ArcGIS地理信息系統(tǒng)［M].北京：清華大學(xué)出版社，2011.

［3]秦 昆.GIS空間分析理論與方法［M].武漢：武漢大學(xué)出版社，2010.

［4]魏旭東，曹先革，楊金玲.平坦地區(qū)GPS高程異常擬合研究［J].測繪與空間地理信息，2011，34（5）：59－62.