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

賈建峰

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

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

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

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

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

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

1 基本原理

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

1.1變差函數(shù)

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

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

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

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

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

1.2普通克立格法

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

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

2 案例分析

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

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

以20個(gè)樣本點(diǎn)作為已知數(shù)據(jù)集的計(jì)算為例，首先啟動(dòng)Arcgis9.0軟件的地質(zhì)統(tǒng)計(jì)學(xué)模塊，點(diǎn)擊Geostatistical Analyst，選擇Geostatistical Wizard（地統(tǒng)計(jì)分析向?qū)В棾鋈缦聦?duì)話框：

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

圖2Geostatistical Method欄中選擇Ordinary Kriging方法下的Prediction Map，其他默認(rèn)后點(diǎn)擊Next，彈出圖3對(duì)話框：

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

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

2.2四種計(jì)算方案對(duì)比

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

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

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

3 結(jié)論

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

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

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

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