基于地統計學的GPS高程異常擬合研究①

周長江，紀志鵬

（中國礦業大學 環境與測繪學院，江蘇 徐州221008）

0 引 言

GPS測量具有精度高、實時、快速等優點，廣泛用于工程測量、變形監測等領域。GPS在平面定位方面已經可以達到很高的精度，而GPS高程控制網精度比較低。GPS高程是以WGS－84為參考橢球的大地高，而實際需要的是以似大地水準面為基準的正常高，這兩者之間的差值叫做高程異常，若能求出高精度的GPS點的高程異常ζ就可以把大地高轉換到正常高。近年來，許多學者圍繞這一目標提出了很多方法，如多項式擬合、BP神經網絡算法［1]。但多項式擬合算法誤差較大，BP算法受初始權值的影響較大，不夠穩定。近年來，ArcGIS的應用范圍越來越廣，它具有強大的空間數據管理、分析、處理能力。ArcGIS中的Geostatistical Analysis Tools是基于地統計學的原理對數據進行處理和分析的工具［2]。用到的方法主要是空間內插方法，按其實現的數學原理分為：確定性插值、克里金插值。確定性插值方法有反距離加權值法、全局性插值法、局部性插值法、徑向基插值法。克里金插值又分為普通克里金、簡單克里金和泛克里金等。結合實例將地統計分析工具中的相關插值方法運用到GPS點的高程異常擬合中，并對它們擬合的精度進行分析比較。

1 確定性插值

1.1 反距離加權值法

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

式中：ζi為相鄰點的高程異常；pi為數據ζi的權值。

1.2 全局多項式插值法

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

1.3 局部多項式插值

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

1.4 徑向基函數插值

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

2 克里金插值法

克里金插值法又稱空間局部插值法，是以變異函數和結構分析為基礎，對有限區域內區域變化量進行無偏最優估計的一種方法。克里金方法又分為普通克里金、簡單克里金、泛克里金等。

2.1 普通克里金估計

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

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

2.2 泛克里金估計方法

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

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

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

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

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

3 空間探索性分析

3.1 數據分布檢驗

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

圖1 樣本數據的正態QQPlot分布圖

3.2 趨勢分析

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

圖2 趨勢分析

3.3 方差變異分析

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

圖3 半變異函數云圖

圖4 協方差函數云圖

3.4 粗差剔除

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

圖5 Voronoi Map

4 實例計算和精度比較

4.1 實例計算

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

利用ArcGIS中的地統計分析模塊中的插值工具，分別采用上述方法對樣本數據進行試驗。由于每種方法中參數選取的不同會對內插的精度產生影響，選取的原則就是使樣本數據通過轉換盡可能服從不同算法的前提假設，并使樣本數據的空間相關性達到最大，從而使擬合的結果均方根最小。為了便于比較，將每種插值方法中所有涉及到的參數進行組合分別算出結果，選取均方根最小的作為該種算法的最終結果進行比較。反復試驗得出當反距離加權值、全局插值方法中的power參數取為3次時均方根最小。局部插值法中power參數取為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時均方根最小。不同算法擬合結果如表1所示（高程異常單位為m，誤差單位為cm，限于篇幅所限，局部插值法數據沒有列出）。

圖6 已知數據點位分布和高程異常

表1 不同方法的結果比較

圖7 高程異常預測曲面

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

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.從實驗結果可以看出徑向基法比較穩定，且誤差比較小，精度比較高，克里金法精度略低于徑向基法。這主要是由于克里金方法對原始觀測數據進行了平滑處理，導致克里金估值的結果的空間結構在整體上擬合原始觀測數據不如徑向基法精度高。此外，本實驗采用樣本數據點比較多，且測區內高程異常曲面變化相對比較緩慢適合采用徑向基方法。

5 結 論

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

［1]楊明清，靳 蕃.用神經網絡方法轉換高程［J].測繪學報，1999，28（4）：301－307.

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

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

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