半方差函數也稱為半變異函數,它是地統計學中研究土壤變異性的關鍵函數.
2.1.1半方差函數的定義和參數
如果隨機函數Z(x)具有二階平穩性,則半方差函數((h)可以用Z(x)的方差S2和空間協方差C(h)來定義:((h)= S2-C(h)
((h)反映了Z(x)中的空間相關部分,它等於所有以給定間距h相隔的樣點測值之差平方的數學期望:
(1)
實際可用:
(2)
式中N(h)是以h為間距的所有觀測點的成對數目.某個特定方向的半方差函數圖通常是由((h)對h作圖而得.在通常情況下,半方差函數值都隨著樣點間距的增加而增大,並在壹定的間距(稱為變程,arrange)升大到壹個基本穩定的常數(稱為基臺,sill).
土壤性質的半方差函數也可能持續增大,不表現出確定的基臺和變程,這時無法定義空間方差,說明存在有趨勢效應和非平穩性.另壹些半方差函數則可能完全缺乏空間結構,在所用的采樣尺度下,樣品間沒有可定量的空間相關性.
從理論上講,實驗半方差函數應該通過坐標原點,但是許多土壤性質的半方差函數在位置趨於零時並不為零.這時的非零值就稱為"塊金方差(Nugget variance)"或"塊金效應".它代表了無法解釋的或隨機的變異,通常由測定誤差或土壤性質的微變異所造成.
對於平穩性數據,基底方差與結構方差之和約等於基臺值.
2.1.2 方差函數的理論模型
土壤在空間上是連續變異的,所以土壤性質的半方差函數應該是連續函數.但是,樣品半方差圖卻是由壹批間斷點組成.可以用直線或曲線將這些點連接起來,用於擬合的曲線方程就稱為半方差函數的理論模型.在土壤研究中常用的模型有:
①線性有基臺模型:
式中C1/a是直線的斜率.這是壹維數據擬合的最簡單模型:
((h)=C0 +C1·h/a 0在極限情況下,C1/a可以為0,這時就有純塊金效應模型:
((h)=C0, h>0 (4)
((0)=0 h=0
②球狀模型
((h)= C0 +C1[1.5h/a-0.5(h/a)3] 0a (5)
((0)=0 h=0
③指數模型
((h)=C0+C1[1-exp-h/a ] h>0 (6)
((0)=0 h=0
④雙曲線模型
(7)
⑤高斯模型
((h)=C0+C1[1-exp(-h2/a2)] h>0 (8)
((0)=0 h=0
選定了半方差函數的擬合模型後,通常是以最小二乘法計算方程的參數,並應用Ross等的最大似然程序(MLP),得到效果最好的半方差方程.
2.1.3 模型的檢驗(cross-validation,又稱作jacknifing)
為了檢驗所選模型三個參數的合理性,必須作壹定的檢驗.但是到現在為止還沒有壹個有效的方法檢驗參數的置信區間;同時,由於我們不知道半方差模型的確切形式,所選定的模型只是半方差函數的近似式,故無法以確切的函數形式對模型參數進行統計檢驗.交叉驗證法的檢驗方法,壹種間接的結合普通克立格的方法,為檢驗所選模型的參數提供了壹個途徑.這個方法的優點是在檢驗過程中對所選定的模型參數不斷進行修改,直至達到壹定的精度要求.
交叉驗證法的基本思路是:依次假設每壹個實測數據點未被測定,由所選定的半方差模型,根據n-1個其它測定點數據用普通克立格估算這個點的值.設測定點的實測值為,估算值為,通過分析誤差,來檢驗模型的合理性.
2.1.4半方差函數的模型的選取原則和參數的確定
半方差函數的模型的選取原則是:首先根據公式計算出((h)的散點圖,然後分別用不同類型的模型來進行擬合,得到模型的參數值及離差平方和,首先考慮離差平方和較小的模型類型,其次,考慮塊金值和獨立間距,最後用交叉驗證法來修正模型的參數.
2.2 Kriging最優內插估值法
如果區域化變量滿足二階平穩或本征假設,對點或塊段的估計可直接采用點克立格法(Puctual Kriging )或者塊段克立格法(Block Kriging).這兩種方法是最基本的估計方法,也稱普通克立格法(Origing Kriging,簡稱OK).
半方差圖除用於分析土壤特性空間分布的方向性和相關距離外,還可用於對未測點的參數進行最優內插估值和成圖,該法原理如下:
Kriging最優內插法的原理
設x0為未觀測的需要估值的點,x1, x2,…, xN 為其周圍的觀測點,觀測值相應為y(x1 ),y(x2),…,y(xN).未測點的估值記為 (x0),它由相鄰觀測點的已知觀測值加權取和求得:
(9)
此處,(i為待定加權系數.
和以往各種內插法不同,Kriging內插法是根據無偏估計和方差最小兩項要求來確定上式中的加權系數(i的,故稱為最優內插法.
1. 無偏估計 設估值點的真值為y(x0).由於土壤特性空間變異性的存在,以及, y(x0)均可視為隨機變量.當為無偏估計時,
(10)
將式(9)代入(10)式,應有
(11)
2. 估值和真值y(x0)之差的方差最小.即
(12)
利用式(3-10),經推導方差為
(13)
式中,((xi,xj)表示以xi和xj兩點間的距離作為間距h時參數的半方差值,((xi, x0)則是以xi和x0兩點之間的距離作為間距h時參數的半方差值.觀測點和估值點的位置是已知的,相互間的距離業已知,只要有所求參數的半方差((h)圖,便可求得各個((xi,xj)和((xi,x0)值.
因此,確定式(9)中各加權系數的問題,就是在滿足式(11)的約束條件下,求目標函數以式(13)表示的方差為最小值的優化問題.求解時可采用拉格朗日法,為此構造壹函數,(為待定的拉格朗日算子.由此,可導出優化問題的解應滿足:
i=1,2,N (14)
由式(14)和式(11)組成n+1階線性方程組,求解此線性方程組便可得到n個加權系數(i和拉格朗日算子(.該線性方程組可用矩陣形式表示:
(15)
式中,( ij為((xi,xj)的簡寫.
求得各(i值和(值後,由式(9)便可得出x0點的最優估值y(x0).而且還可由式(13)求出相應該估值的方差之最小值(2min.將式(14)代入式(13),最小方差值還可由下式方便地求出:
(16)
上述最優化問題求解還可用其他方法,在應用Kriging內插法時還有其他方面的問題,在此都不壹壹列舉了.