地質(zhì)統(tǒng)計學(xué)Kriging法在GPS高程擬合中的應(yīng)用探討

賈建峰

摘要：本文以某工業(yè)園區(qū)項目的GPS水準測量數(shù)據(jù)，利用地質(zhì)統(tǒng)計學(xué)中的Ordinary Kriging方法，進行了GPS高程異常擬合的不同方案試驗研究，結(jié)果發(fā)現(xiàn)普通克立格方法合理擬合方案下的估計精度可滿足測量規(guī)范要求，推薦應(yīng)用于GPS高程異常擬合工程。

關(guān)鍵詞：地質(zhì)統(tǒng)計學(xué) GPS高程擬合

GPS高程轉(zhuǎn)換實質(zhì)上主要是通過求取地面點的高程異常值，然后將外業(yè)GPS測量的大地高減去內(nèi)業(yè)獲取的高程異常值即可得到實用的正常高。而高程異常是地球重力場的重要參數(shù)之一，從理論上講，實現(xiàn)GPS大地高向正常高轉(zhuǎn)換的最好方法是綜合利用GPS測量數(shù)據(jù)、重力測量數(shù)據(jù)、地球重力場模型和其他數(shù)學(xué)方法，而對于一般單位而言，在無法獲取必要的重力資料的情況下，數(shù)學(xué)擬合方法仍然是目前進行GPS高程轉(zhuǎn)換的首選方案。

關(guān)于GPS高程的轉(zhuǎn)換問題，測繪界的許多專家學(xué)者提出了多種解決辦法。這些方法歸納起來主要有：重力法、GPS水準法（數(shù)學(xué)模型擬合法）、數(shù)學(xué)模型抗差估計法、數(shù)學(xué)模型優(yōu)化方法、GPS三角高程法、平差轉(zhuǎn)換法、整體平差法、神經(jīng)網(wǎng)絡(luò)法等。如常用的移動加權(quán)平均法、曲線擬合法、曲面擬合法均屬于數(shù)學(xué)模型擬合法。以上諸多方法在應(yīng)用過程中均取得了一定的成效，但各自也都存在一些缺點。主要由于受地球區(qū)域密度分布異常、地表地形的復(fù)雜程度等因素的影響，有的方法理論精度很高，但與實際情況的吻合程度不高；而有的方法計算過程比較復(fù)雜，且實際效果也不夠理想。因此，如何來進行高精度的GPS高程擬合與轉(zhuǎn)換仍然是測量工程實踐中需要不斷研究和探討的問題。

本文在介紹地質(zhì)統(tǒng)計學(xué)基本理論及普通克立格Ordinary Kriging方法的原理及應(yīng)用基礎(chǔ)上，以某測區(qū)的實測GPS高程異常數(shù)據(jù)為例，進行了基于普通克立格方法的GPS高程擬合的相關(guān)試驗分析與研究，按照不同已知高程異常點數(shù)據(jù)個數(shù)設(shè)計了四種方案，并對四種方案的擬合結(jié)果進行對比分析，證實了普通克立格方法用于GPS高程異常擬合的可行性，并探討了普通克立格方法高程異常擬

合的精度，以及已知高程異常點數(shù)量、分布與高程擬合精度的關(guān)系。

1 基本原理

地質(zhì)統(tǒng)計學(xué)是在經(jīng)典統(tǒng)計方法的基礎(chǔ)上，充分考慮地質(zhì)變量的空間變化特征-相關(guān)性和隨機性，并以反映地質(zhì)現(xiàn)象區(qū)域化的隨機函數(shù)-變差函數(shù)作為基本工具，來研究地質(zhì)和采礦等工作中的各種問題。在結(jié)構(gòu)分析的基礎(chǔ)上采用各種Kriging法，來估計實際問題的。克立格法（Kriging）是一種推求最優(yōu)、線性、無偏估計量的方法。即Kriging法是在考慮信息樣品的形狀、大小及其與待估塊段間的空間分布位置等幾何特征以及品位的空間結(jié)構(gòu)之后，為了達到線性、無偏和最小估計方差的估計，而對每一樣品值分別賦予一定的權(quán)系數(shù)，最后進行加權(quán)平均來估計塊段品位的方法。

1.1變差函數(shù)

變差函數(shù)既能描述區(qū)域化變量的結(jié)構(gòu)性變化，又能捕述其隨機性變化，且它的計算是許多其他地質(zhì)統(tǒng)計學(xué)計算的基礎(chǔ)。

變差函數(shù)的定義如下：變差函數(shù)為區(qū)域化變量Z（x）和Z（x+h）的增量平方的數(shù)學(xué)期望，即區(qū)域化變量增量的方差。當區(qū)域化變量滿足二階平穩(wěn)假設(shè)、本征假設(shè)或（準）二階平穩(wěn)假設(shè)、（準）本征假設(shè)時，變差函數(shù)只依賴于步長h，與x取值無關(guān)，就可以得到一維實驗變差函數(shù)的計算公式：

式（1）即為Matheron推薦的傳統(tǒng)計算公式，式中：N（h）是x軸上相隔h的點的對數(shù)，Z（x_i）和Z（x_i+h）是觀測值Z（x）和Z（x+h）的N（h）對現(xiàn)實。如果Z（x）是定義在二維、三維空間的區(qū)域化變量，則x是二維、三維空間中的點，h是二維、三維空間中的向量。

變差函數(shù)不僅是Kriging估值、條件模擬等計算的基礎(chǔ)，更重要的是它反映和刻畫了區(qū)域化變量的許多性質(zhì)，是分析區(qū)域化變量空間變異性的重要工具。一般情況下，代表高程異常的變差函數(shù)事先未知，可選用實用變差函數(shù)模型進行描述，然后根據(jù)實測資料進行擬合。設(shè)隨機場是二階平穩(wěn)且為統(tǒng)計各向同性的，則常用的理論變差函數(shù)與交叉協(xié)方差函數(shù)理論球狀模型如下：

式（2）中：C₀為塊金常數(shù)（Nugget Effect），塊金常數(shù)反映了變差函數(shù)在空間上的突變性，產(chǎn)生這種現(xiàn)象的原因可能是由于空間不相關(guān)的觀測誤差以及小于取樣間距尺度上的空間變異引起的；C₀+C為基臺值（Sill），C為拱高。變差函數(shù)基臺值的大小可反映變量在該方向上變化幅度或總的變異程度的大小。塊金常數(shù)大小可反映區(qū)域化變量隨機性大小。塊金值與基臺值之比C₀/（C₀+C）反映塊金方差占總空間異質(zhì)性變異的大小。

1.2普通克立格法

將所研究的區(qū)域化變量，如高程異常記為U（x），且U（x）滿足二階平穩(wěn)假設(shè)或本征假設(shè)，E[U（x）]=m為未知常數(shù)，對區(qū)域內(nèi)任一點x₀處的變量U（x₀）進行最優(yōu)、無偏、線性估計的方法稱為普通克立格法，且x₀處的估計量可表示為：

對高程異常的變異函數(shù)和結(jié)構(gòu)分析結(jié)果表明高程異常存在空間相關(guān)性，可視為區(qū)域化變量，則可利用普通克立格方法進行內(nèi)插或外推。該方法的實質(zhì)是根據(jù)高程異常這一區(qū)域化變量的原始數(shù)據(jù)和變異函數(shù)的結(jié)構(gòu)特點，利用已知樣點對未知樣點的高程異常進行線性無偏、最優(yōu)估計。估計時考慮了已知樣點與未知高程異常樣點的相互空間位置關(guān)系，以及變異函數(shù)提供的結(jié)構(gòu)信息。

2 案例分析

從某測區(qū)的靜態(tài)GPS控制網(wǎng)獲得100個GPS高程異常數(shù)據(jù)點，分別選用20、30、40、50個樣本點作為已知數(shù)據(jù)集，其余為未知點，用于外符合精度的檢驗。本文借助Arcgis9.0軟件成熟的地質(zhì)統(tǒng)計學(xué)模塊，計算并擬合變差函數(shù)，并進行高程異常的空間性相關(guān)分析；最后利用普通克立格方法擬合高程異常，并進行精度分析與方案對比。探討普通克立格方法在GPS高程異常擬合中的問題。

2.1Arcgis9.0軟件的地質(zhì)統(tǒng)計學(xué)模塊簡介

以20個樣本點作為已知數(shù)據(jù)集的計算為例，首先啟動Arcgis9.0軟件的地質(zhì)統(tǒng)計學(xué)模塊，點擊Geostatistical Analyst，選擇Geostatistical Wizard（地統(tǒng)計分析向?qū)В棾鋈缦聦υ捒颍?/p>

圖1Method欄中選擇Kriging方法，Dataset1中Inputdata選擇事先做好的已知高程異常Excel數(shù)據(jù)表，Attribute選擇“高程異常”，點擊Next，彈出圖2對話框：

圖2Geostatistical Method欄中選擇Ordinary Kriging方法下的Prediction Map，其他默認后點擊Next，彈出圖3對話框：

圖3選擇變異函數(shù)模型為球狀模型，不考慮帶狀各向異性，調(diào)節(jié)Number of Lags和Lag Size，獲得左上角的綠色變異函數(shù)擬合曲線，用于克立格估計。20個已知點數(shù)據(jù)變異函數(shù)塊金值0，基臺值0.13476，變程0.24143。獲得理想的變異函數(shù)模型后，點擊Next，進行交叉驗證，彈出圖4對話框：

由圖4可見20個高程異常點的預(yù)測值與已知值的分布曲線接近于綠色理論曲線，可由預(yù)測值與已知值的差值計算分析估計誤差，除有一點高程異常擬合誤差接近0.40m外，其余21個點擬合誤差均小于0.10m，與普通水準測量精度相當。點擊Finish后，完成克立格估計，生成測區(qū)高程異常差值圖5。

2.2四種計算方案對比

四種方案變差函數(shù)圖計較，如下圖：

變差函數(shù)的計算采用球狀模型擬合，四種方案的變差函數(shù)均具有較強的空間的結(jié)構(gòu)性。變差函數(shù)塊金值C₀均為零，表明高程異常測量誤差可以忽略不計，高程異常這一區(qū)域化變量的隨機性較弱。C₀+C為基臺值，變差函數(shù)基臺值的大小可反映變量在該方向上變化幅度或總的變異程度的大小，表1反映出高程異常具有很強的空間相關(guān)性，其空間相關(guān)半徑（變程值）約在0.277附近。

表2顯示四種方案的高程異常擬合外符合精度0.017m-0.047m，當已知點數(shù)量由20%增加至50%時，外符合精度逐漸提高，均為厘米水平，達到了四等水準測量的精度。而內(nèi)符合精度0.104m-0.044m，當已知點數(shù)量由20%增加至50%時，內(nèi)符合精度逐漸提高，除20%已知點方案外，其余均為厘米水平；值得注意的是當已知點為40%時，再增加已知點數(shù)量對精度的提高并無實際意義。

3 結(jié)論

本文建立了基于Kriging方法的高程異常變差函數(shù)模型，利用普通克立格公式進行插值計算，最后利用實測數(shù)據(jù)進行了交叉驗證與精度分析。借助ArcGIS地質(zhì)統(tǒng)計學(xué)模塊中克立格插值工具，采用普通克立格（Ordinary Kriging）方法對本例數(shù)據(jù)進行高程異常擬合的試驗計算，得到以下結(jié)論：

（1）高程異常數(shù)據(jù)的統(tǒng)計分布接近正態(tài)分布，高程異常的變差函數(shù)具有空間較強的空間結(jié)構(gòu)性，可將其視為區(qū)域化變量，利用最佳擬合的變差函數(shù)，對未知高程異常點進行普通克立格插值。

（2）從四種方案的高程異常擬合圖以及交叉誤差表，可直觀地分析擬合成果以及誤差分布情況，進一步統(tǒng)計分析后，發(fā)現(xiàn)普通克立格方法具有較高的內(nèi)、外符合精度，其外符合精度均達到了厘米級水平，與四等水準測量精度相當。

（3）隨著已知數(shù)據(jù)的增多，高程異常估計精度逐漸提高，但當已知點為總采樣點數(shù)的40%時，再增加已知點數(shù)量對高程異常擬合精度的提高已無實際意義。建議做高程擬合的已知點數(shù)量控制在總點數(shù)的40%以內(nèi)即可。