插值、平穩(wěn)假設(shè)、變異函數(shù)、基臺(tái)、塊金、克里格…地學(xué)計(jì)算概念及公式推導(dǎo)

1、image插值、平穩(wěn)假設(shè)、變異函數(shù)、基臺(tái)、塊金、克里格地學(xué)計(jì)算概念及公式推導(dǎo)1引言最近的幾篇博客，分別從多光譜與高光譜遙感的實(shí)際應(yīng)用出發(fā)，對(duì)影像前期處理與相關(guān)算法、反演操作等加以詳細(xì)介紹。而通過遙感手段獲取了豐富的各類地表信息數(shù)據(jù)后，如何對(duì)數(shù)據(jù)加以良好的數(shù)學(xué)處理與科學(xué)分析，同樣是我們需要重視的問題。因此，準(zhǔn)備由這一篇博客入手，新建一個(gè)專欄，逐篇地對(duì)地學(xué)計(jì)算方面的內(nèi)容加以初步總結(jié)。那么首先，我們就由地學(xué)計(jì)算的幾個(gè)基本概念入手，對(duì)相關(guān)理論方面的內(nèi)容加以一定了解。需要注意的是，以下內(nèi)容如果單獨(dú)來看或許有些不好理解，但一旦將其與實(shí)際應(yīng)用結(jié)合，便會(huì)豁然開朗。其中，具體的實(shí)際應(yīng)用部分我將會(huì)在后面的博客中涉

2、及。2空間插值空間數(shù)據(jù)的獲取是進(jìn)行空間分析的基礎(chǔ)與起源。為了提高研究結(jié)論的精度，我們亦總是希望能夠獲取研究區(qū)域內(nèi)更多、更全面的精確空間屬性數(shù)據(jù)信息。然而，在實(shí)際研究、工作中，由于人力、成本、資源等外部條件限制，我們不可能對(duì)全部未知區(qū)域加以采樣與測(cè)量，而往往只能得到研究區(qū)域內(nèi)有限數(shù)量的采樣點(diǎn)及其相關(guān)屬性數(shù)據(jù)。因此，往往可以考慮選取合適的空間采樣點(diǎn)，利用一定數(shù)學(xué)模型，依據(jù)已知采樣點(diǎn)各自對(duì)應(yīng)屬性數(shù)據(jù)對(duì)研究區(qū)域所有位置的未知屬性信息加以預(yù)測(cè)?？臻g插值(SpatialInterpolation)即可以實(shí)現(xiàn)這一需求。其是一種將離散采樣點(diǎn)測(cè)量數(shù)據(jù)轉(zhuǎn)換為連續(xù)數(shù)據(jù)曲面的常用方法，包括內(nèi)插(Interpolati

3、on)和外推(Extrapolation)兩種應(yīng)用形式。一般地，對(duì)樣本點(diǎn)范圍以內(nèi)(即所有采樣點(diǎn)最大外接矩形內(nèi)部)的空間進(jìn)行插值，才可稱作“內(nèi)插”(部分文獻(xiàn)亦直接用“插值”代替“內(nèi)插”)；反之則稱“外推”或“預(yù)測(cè)”，往往認(rèn)為外推的結(jié)果誤差較大?？臻g插值理論及其方法基于著名的“地理學(xué)第一定律(ToblersFirstLawofGeography)”，即一般地，距離越近的地物具有越高的相關(guān)性。這一至今已產(chǎn)生深遠(yuǎn)影響的地學(xué)定律最早由美籍瑞士地理學(xué)家WaldoR.Tobler教授于1970年提出。在各方法所對(duì)應(yīng)的數(shù)學(xué)計(jì)算原理層面，空間插值一般可以分為確定性插值方法(DeterministicInterp

4、olation)與地統(tǒng)計(jì)插值方法(Geostatistics，亦稱非確定性插值方法)兩種。其中，確定性插值方法基于研究區(qū)域內(nèi)各信息點(diǎn)之間相似程度或整個(gè)曲面的平滑程度，從而創(chuàng)建連續(xù)的擬合曲面；其依據(jù)插值計(jì)算時(shí)納入考慮的采樣點(diǎn)分布范圍，又可進(jìn)一步分為整體插值法與局部插值法。地統(tǒng)計(jì)插值方法則是基于研究區(qū)域內(nèi)各信息點(diǎn)的綜合統(tǒng)計(jì)學(xué)規(guī)律，以變異函數(shù)(Variogram)理論與結(jié)構(gòu)分析為基礎(chǔ)，實(shí)現(xiàn)其屬性的空間自相關(guān)性定量化，從而創(chuàng)建得出連續(xù)插值曲面。在所創(chuàng)建連續(xù)插值表面通過全部采樣點(diǎn)的與否層面，空間插值一般又可以分為精確性插值與非精確性插值兩種。其中，前者預(yù)測(cè)樣點(diǎn)的屬性數(shù)值與其各自實(shí)測(cè)值相等，即其采樣點(diǎn)屬性

5、數(shù)據(jù)全部落入預(yù)測(cè)結(jié)果曲面；后者預(yù)測(cè)樣點(diǎn)的屬性數(shù)值則往往不與其各自實(shí)測(cè)值相等，即其采樣點(diǎn)屬性數(shù)據(jù)一般不會(huì)落入預(yù)測(cè)結(jié)果曲面。因此，使用非精確性插值方法往往可以避免在預(yù)測(cè)表面中出現(xiàn)明顯的波峰或波谷，整體呈現(xiàn)出平緩態(tài)勢(shì)。3幾個(gè)重要假設(shè)地學(xué)計(jì)算中，幾個(gè)重要假設(shè)可以說是關(guān)鍵中的關(guān)鍵；而其往往也比較難以理解，不怕，我們慢慢往下看。3.1平穩(wěn)假設(shè)平穩(wěn)假設(shè)(StationaryAssumption)是指，一組觀測(cè)值的均值是始終固定的，其與觀測(cè)值所在位置無關(guān)；將既定的某個(gè)點(diǎn)集由某一研究區(qū)域內(nèi)某處移動(dòng)至另一處時(shí)，隨機(jī)函數(shù)的性質(zhì)保持不變。即隨機(jī)函數(shù)的分布規(guī)律不因位置的改變而改變，具有嚴(yán)格的平穩(wěn)性。平穩(wěn)性假設(shè)的公式表達(dá)

6、為：image其中，F(xiàn)_(x_1，,x_n)(z_1,z_n)表示位置在(x_1，,x_n)上的點(diǎn)集(z_1,z_n)對(duì)應(yīng)的隨機(jī)函數(shù)。3.2二階平穩(wěn)性假設(shè)二階平穩(wěn)性假設(shè)(SecondStationaryAssumption)又稱弱平穩(wěn)假設(shè)，其與協(xié)方差函數(shù)(CovarianceFunction)有關(guān)。這一假設(shè)認(rèn)為，隨機(jī)函數(shù)的均值為一常數(shù)，且任意兩個(gè)隨機(jī)變量之間的協(xié)方差僅僅依賴于其二者之間的距離與方向，而與其具體位置無關(guān)。將上述兩個(gè)條件用公式分別表達(dá)為：EZ(z)二肌QCovZ(x)rZ(x+h)=EZ(x)Z(x+h)-m2二image其中，EZ(x)為區(qū)域化變量Z(x)的數(shù)學(xué)期望，CovZ(x

7、),Z(x+h)為區(qū)域化變量Z(x)與Z(x+h)所對(duì)應(yīng)的協(xié)方差函數(shù)，C(h)為僅與h有關(guān)的協(xié)方差取值，m為一常數(shù)，h為滯后距。上述二階平穩(wěn)性假設(shè)是針對(duì)整個(gè)研究區(qū)域范圍而言。若區(qū)域化變量?jī)H僅在整個(gè)研究區(qū)域內(nèi)的某個(gè)有限區(qū)域中滿足上述條件，即條件僅在局部區(qū)域成立，則稱此區(qū)域化變量滿足準(zhǔn)二階平穩(wěn)性假設(shè)(QuasiSecondStationaryAssumption)。準(zhǔn)二階平穩(wěn)性假設(shè)可以視作一種折中方案，既考慮到平穩(wěn)范圍的大小，又顧及到有效數(shù)據(jù)的數(shù)量。3.3本征假設(shè)本征假設(shè)(IntrinsicHypothesis)又稱內(nèi)蘊(yùn)假設(shè)，其與變異函數(shù)有關(guān)。這一假設(shè)認(rèn)為，區(qū)域化變量的增量滿足以下兩個(gè)條件：在整個(gè)

8、研究區(qū)域內(nèi)，區(qū)域化變量增量的數(shù)學(xué)期望為0;且其方差函數(shù)存在，并只依賴于滯后距，而與所處位置無關(guān)。將上述兩個(gè)條件用公式分別表達(dá)為：Z(z)-Z(x+a)=00VarZ(x)-Z(x+K)=2y(h)image當(dāng)E(x)存在時(shí)，上述第一個(gè)公式可以寫作：EZ(x)=EZ(x+h)二image其中，VarZ(x)-Z(x+h)為區(qū)域化變量Z(x)與Z(x+h)所對(duì)應(yīng)的方差函數(shù)，Y(h)為區(qū)域化變量在滯后距為h時(shí)的變異函數(shù)，m為一常數(shù)，h為滯后距，其它符號(hào)同前述意義。同樣的，上述本征假設(shè)亦是針對(duì)整個(gè)研究區(qū)域范圍而言。若區(qū)域化變量?jī)H僅在整個(gè)研究區(qū)域內(nèi)的某個(gè)有限區(qū)域中滿足上述條件，則稱此區(qū)域化變量滿足準(zhǔn)本征

9、假設(shè)(QuasiIntrinsicHypothesis)。與準(zhǔn)二階平穩(wěn)性假設(shè)類似，準(zhǔn)本征假設(shè)亦可視作一種折中方案，其同樣既考慮到了本征image假設(shè)對(duì)應(yīng)范圍的大小，又顧及到了有效數(shù)據(jù)的數(shù)量。此外，本征假設(shè)是地統(tǒng)計(jì)學(xué)中對(duì)隨機(jī)函數(shù)的基本假設(shè)。3.4不同假設(shè)對(duì)比結(jié)合上述二階平穩(wěn)性假設(shè)與本征假設(shè)相關(guān)原理，可以看到兩種假設(shè)的討論對(duì)象具有一定區(qū)別：二階平穩(wěn)性假設(shè)更多是討論某一特定研究區(qū)域內(nèi)區(qū)域化變量自身【即Z(x)】的特征，而本征假設(shè)則是研究區(qū)域化變量所對(duì)應(yīng)增量【即Z(x)-Z(x+h)】的特征。般認(rèn)為，二階平穩(wěn)性假設(shè)對(duì)區(qū)域化變量的要求較之本征假設(shè)更為嚴(yán)格，即若某個(gè)研究區(qū)域內(nèi)的某一區(qū)域化變量滿足二階平穩(wěn)性

10、假設(shè)，則其一定滿足本征假設(shè)；反之則反，若僅知區(qū)域化變量滿足本征假設(shè)，則其不一定滿足二階平穩(wěn)性假設(shè)。再結(jié)合平穩(wěn)假設(shè)，上述三種假設(shè)的嚴(yán)格程度由大至小依次排列為：平穩(wěn)假設(shè)、二階平穩(wěn)性假設(shè)與本征假設(shè)。此外，結(jié)合二階平穩(wěn)性假設(shè)的兩個(gè)條件，還可以推出協(xié)方差函數(shù)與變異函數(shù)之間的關(guān)系：Y(Ji)=C(0)-image其中，廠(h)為區(qū)域化變量所對(duì)應(yīng)變異函數(shù)，C(O)為距離為0時(shí)此區(qū)域化變量所對(duì)應(yīng)的協(xié)方差取值【即基臺(tái)值y(8)】，C(h)為距離為h時(shí)此區(qū)域化變量所對(duì)應(yīng)的協(xié)方差取值。由這一關(guān)系可知，用以衡量某一區(qū)域化變量在相距為h的兩空間位置點(diǎn)分別取值的自相關(guān)性的指標(biāo)一一協(xié)方差函數(shù)與變異函數(shù)之間具有相互關(guān)系。因此

11、，在滿足二階平穩(wěn)性假設(shè)的條件下，若協(xié)方差函數(shù)平穩(wěn)，則可知變異函數(shù)平穩(wěn)，即其取值只與滯后距h有關(guān)。4變異函數(shù)克里格插值法需要借助空間數(shù)據(jù)的試驗(yàn)變異函數(shù)及其散點(diǎn)圖特點(diǎn)，因此變異函數(shù)的計(jì)算在克里格插值過程中發(fā)揮著重要作用；變異函數(shù)及其模型擬合對(duì)克里格插值結(jié)果精度具有較大影響。變異函數(shù)(Variogram)，又稱為半變異函數(shù)、半方差函數(shù)(Semi-variogram)等，其用以描述區(qū)域化變量的空間變化特征與強(qiáng)度，被定義為區(qū)域化變量增量平方的數(shù)學(xué)期望。在一維條件下，直接將區(qū)域化變量Z在位置(x)與(x+h)處的取值Z(x)與Z(x+h)之差的方差定義為變異函數(shù)，其因變量為距離h;而在二維或三維條件下，可

12、以將上述一維中具有單一方向的距離h進(jìn)一步引申為在任意方向a上的距離|h|。具體公式表達(dá)為如下。2*y(x,k)=VarZ(x)一Z(x+h)。image其中，廠(x,h)即為變異函數(shù)。由于公式中在其前具有一個(gè)系數(shù)“2”，因此其亦被稱作半變異函數(shù)。結(jié)合(準(zhǔn))二階平穩(wěn)性假設(shè)、(準(zhǔn))本征假設(shè)等地學(xué)基本假設(shè)，變異函數(shù)取值與區(qū)域化變量樣點(diǎn)所處位置x無關(guān)，僅僅與樣點(diǎn)之間的距離h有關(guān)，則變異函數(shù)可以寫作：呦#=Sz(Xi)z(Xi+評(píng)。i=limage其中，廠(h)#為區(qū)域化變量Z(x)的變異函數(shù)，N(h)為區(qū)域化變量樣點(diǎn)集中距離為h的點(diǎn)對(duì)數(shù)量，x_i為距離為h時(shí)所對(duì)應(yīng)的第i個(gè)點(diǎn)。在這里，距離h亦被稱作滯后

13、距(LagDistance)。般地，區(qū)域化變量變異函數(shù)圖像往往呈現(xiàn)出“先快速上升，再增速減緩，后趨于平穩(wěn)”的曲線特征。其具有三個(gè)十分重要的相關(guān)概念，分別為塊金常數(shù)(Nugget)、基臺(tái)值(Sill)與變程(Range)。塊金常數(shù)代表區(qū)域化變量的隨機(jī)性大小。由理論角度，在間距為0(即滯后距為零)時(shí)，區(qū)域化變量采樣點(diǎn)數(shù)值應(yīng)當(dāng)相等；而在間距無限趨近于0時(shí)，對(duì)應(yīng)變異函數(shù)數(shù)值應(yīng)當(dāng)亦向0趨近。但是，在實(shí)際研究中，試驗(yàn)變異函數(shù)在滯后距為0時(shí)，其取值并不為0，而是一個(gè)大于0的數(shù)值。這一數(shù)值便稱為塊金常數(shù)。一般地，上述塊金效應(yīng)的產(chǎn)生可以歸因于測(cè)量誤差，或小于采樣間隔距離處的空間變化。基臺(tái)值用以衡量區(qū)域化變量變化

14、幅度的大小。當(dāng)滯后距無限增大并到達(dá)某一程度后，試驗(yàn)變異函數(shù)若趨于平穩(wěn)，則這一平穩(wěn)水平所對(duì)應(yīng)的數(shù)值即為基臺(tái)值。然而，并不是所有的區(qū)域化變量均具有基臺(tái)值一一如無基臺(tái)值模型對(duì)應(yīng)的變異函數(shù)。變程用以衡量區(qū)域化變量自相關(guān)范圍的大小。當(dāng)滯后距無限增大并到達(dá)某一程度后，試驗(yàn)變異函數(shù)若趨于平穩(wěn)，則此時(shí)對(duì)應(yīng)的滯后距即為變程。其中，小于變程的距離所對(duì)應(yīng)的樣本位置與空間自相關(guān)，而大于變程的距離所對(duì)應(yīng)樣本位置不存在空間自相關(guān)。此外，變異函數(shù)還有其它相關(guān)指標(biāo)，如基臺(tái)值與塊金常數(shù)的差值一一偏基臺(tái)值(PartialSill)，用以衡量空間變異性程度的塊金常數(shù)與基臺(tái)值的比值塊金系數(shù)等?；诓煌瑓^(qū)域化變量對(duì)應(yīng)的變異函數(shù)特征，可

15、以將其分為不同類別。依據(jù)變異函數(shù)基臺(tái)值的有無，可以將模型分為有基臺(tái)值模型、無基臺(tái)值模型與孔穴效應(yīng)模型。其中，有基臺(tái)值模型可以依據(jù)變異函數(shù)的特征進(jìn)一步分為純塊金效應(yīng)模型(PureNuggetEffectModel)、球狀模型(SphericalModel)、指數(shù)模型(ExponentialModel)、高斯模型(CubicModel或GaussianModel)與線性有基臺(tái)模型(LinearwithSillModel)等。上述幾種模型中，較為常用的模型包括球狀模型、指數(shù)模型與高斯模型。此外，同樣依據(jù)變異函數(shù)的特征，無基臺(tái)值模型還可進(jìn)一步分為線性無基臺(tái)值模型、冪指數(shù)模型與對(duì)數(shù)模型等。同樣的，孔穴效

16、應(yīng)模型可分為基臺(tái)值模型和無基臺(tái)值模型。同時(shí)，針對(duì)某種區(qū)域化變量而言，其在不同方向、不同滯后距情況下可能受到不同因素影響；套合結(jié)構(gòu)可以很好解決這一問題。套合結(jié)構(gòu)可以表示為多個(gè)變異函數(shù)之和，每一個(gè)變異函數(shù)均代表著某種方向或某一尺度中的變異性，從而對(duì)區(qū)域化變量的特征加以更好概括。5克里格插值克里格插值法(KrigingMethod)又稱為空間局部插值法，是以上述變異函數(shù)理論及其結(jié)構(gòu)分析為基礎(chǔ)，在有限區(qū)域內(nèi)對(duì)區(qū)域化變量進(jìn)行線性無偏最優(yōu)估計(jì)(BestLinearUnbiasedPrediction，BLUP)的一種方法，在地統(tǒng)計(jì)學(xué)中也被稱為空間最優(yōu)無偏估計(jì)器(SpatialBLUP)。其中，上述“線性”

17、是指克里格插值法對(duì)未知點(diǎn)屬性數(shù)值的估計(jì)采用線性估計(jì)，其公式如下：1=1image其中，(z_0)是區(qū)域化變量在點(diǎn)(x_0,y_0)位置處的預(yù)測(cè)值，入為第i個(gè)已知點(diǎn)的權(quán)重系數(shù)，z_i為第i個(gè)已知點(diǎn)的實(shí)測(cè)值。上述權(quán)重幾可以使得各點(diǎn)處的預(yù)測(cè)值與實(shí)測(cè)值之間的方差最小，而求取這一權(quán)重則為克里格插值的主要內(nèi)容。上述“無偏”是指區(qū)域化變量在各點(diǎn)上估計(jì)量的數(shù)學(xué)期望等于其在同一位置上的真實(shí)值，公式如下:結(jié)合本征假設(shè)，無偏性可以表示為:E（岔一石）二0儀image71A1=1image上述“最優(yōu)”是指區(qū)域化變量在各點(diǎn)上估計(jì)量與其在同一位置上的真實(shí)值的方差最小，公式如下:min一z0）image其中，上述方差被稱作

18、估計(jì)方差或估值方差，是對(duì)估值準(zhǔn)確程度的一種定量表示；而在克里格插值方法中，又可以稱為克里金方差。后續(xù)將克里金方差記作o_kA2o經(jīng)過統(tǒng)計(jì)學(xué)相關(guān)推導(dǎo)，可以將克里金方差寫作:n71時(shí)=2*久：*y（xf,觀）-心*/.j*Y（xb丐）-y（x0,心嚴(yán)image由此轉(zhuǎn)換為在無偏條件約束下的最小值求解問題。引入拉格朗日乘子0,構(gòu)造拉格朗日函數(shù):nL（A1/12/A3j兀爭(zhēng)）二於-2*0歐*1=1image對(duì)權(quán)值與拉格朗日乘子求一階偏微分:偏微分求解結(jié)果為:image+(p=yfe-x0)J=lf2f-,n-*image將上述求和部分展開：h*心磯）+h水血勺）+九*仇觀）+卩二y（呵心W右*血羽）+血

19、*rfe-倉(cāng)）+九*yfeiO+卩二f（切丸W右*rfe,曲）+血事f（勿忑）+九歩譏呦屁）+卩二y（呦心W:;qI:白ahL水v（x.x+彳“水訶+/L*譏疋王J+m二v（x.xAimage可將上式進(jìn)一步寫為矩陣相乘形式，即化為:恥冊(cè)二恥image其中，A代表在原有變異函數(shù)矩陣中額外添加全1行與全1列（交界處1換為0）后的矩陣，幾代表各權(quán)重組成的列向量，0代表前述分析引入的拉格朗日乘子，B為各位置與待求解位置對(duì)應(yīng)距離的變異函數(shù)值組成的列向量，且在列尾增加一個(gè)1。由此，即將上述函數(shù)轉(zhuǎn)化為（n+1）個(gè)未知數(shù)、（n+1）個(gè)表達(dá)式組成的方程組；通過矩陣求逆，求解方程組即可得到待求解位置與其它已知點(diǎn)的權(quán)重。對(duì)每一個(gè)待插點(diǎn)進(jìn)行同樣操作，完成克里格插值。6回歸克里格正如上述分析，普通克里格（OrdinaryKriging）方法更多依靠所選采樣點(diǎn)數(shù)據(jù)對(duì)空間加以插值，其插值效果較依賴于采樣點(diǎn)的個(gè)數(shù)、密度，以及其數(shù)據(jù)的準(zhǔn)確性；另一方面，許多空間屬性往往受到其它環(huán)境變量的影響。例如，土壤有機(jī)碳含量與氣溫、海拔、降水量、坡度等多種環(huán)境因子在0.05水平顯著相關(guān)3,4;近地面氣溫與海拔、海陸距離、NDVI等環(huán)境因子在0.01水平顯著相關(guān)5。由此觀之，若簡(jiǎn)單地忽略環(huán)境要素對(duì)待插值空間屬性的影響，可能會(huì)降低最終插值結(jié)果精度。基于這種考慮，可以