插值、平穩(wěn)假設、變異函數(shù)、基臺、塊金、克里格…地學計算概念及公式推導

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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