1、1 第二講第二講 克里金方法（克里金方法（Kriging）, 是以南非礦業是以南非礦業工程師工程師D.G.Krige (克里格克里格)名字命名的一項名字命名的一項實用空間估計技術，是實用空間估計技術，是地質統計學地質統計學 的重要的重要組成部分，也是地質統計學的核心。組成部分，也是地質統計學的核心。 2主要是為解決礦床儲量計算和誤差估計問題而主要是為解決礦床儲量計算和誤差估計問題而發展起來的發展起來的 由法國巴黎國立高等礦業學院由法國巴黎國立高等礦業學院G馬特隆教授于馬特隆教授于1962年所創立。年所創立。 3H. S. Sichel (1947)D.G. Krige (1951)Krigin

2、g法法（克里金法，克立格（克里金法，克立格法）法）:“根據樣品空間位置不同、樣根據樣品空間位置不同、樣品間相關程度的不同，對每個樣品品間相關程度的不同，對每個樣品品位賦予不同的權，進行滑動加權品位賦予不同的權，進行滑動加權平均，以估計中心塊段平均品位平均，以估計中心塊段平均品位”G. Materon(1962)提出了提出了“地質統計學地質統計學”概念概念 (法文法文Geostatistique)發表了專著發表了專著應用地質統計學論應用地質統計學論。闡明了一整套區域化變量的理論，闡明了一整套區域化變量的理論，為地質統計學奠定了理論基礎。為地質統計學奠定了理論基礎。區域化變量理論區域化變量理論克里

3、金估計克里金估計隨機模擬隨機模擬應用統計學方法研究金礦品位應用統計學方法研究金礦品位1977年我國開始引入年我國開始引入 4 克里金插值方法克里金插值方法niiixzxz10*井眼地震（普通克里金）（應用（應用隨機函數隨機函數理論）理論） 不僅考慮待估點位置與不僅考慮待估點位置與 已知數據位置的相互關已知數據位置的相互關 系，而且還考慮變量的系，而且還考慮變量的 空間相關性。空間相關性。 5 為一個實值變量，可根據概率分布取不同的值。為一個實值變量，可根據概率分布取不同的值。每次取值（觀測）結果每次取值（觀測）結果z為一個確定的數值，稱為為一個確定的數值，稱為隨機變量隨機變量Z的的一個實現。一

4、個實現。P1. 隨機變量隨機變量6連續變量：連續變量：累積分布函數（cdf） cumulative distribution function)(Pr);(zuZobzuF條件累積分布函數（ccdf）后驗 conditional cumulative distribution function)( |)(Pr)( |;(nzuZobnzuF離散變量（類型變量）：離散變量（類型變量）：)( |)(Pr)( |;(nkuZobnkuFZ (u)PP不同的取值方式：估計（estimation） 模擬（simulation）7連續型地質變量連續型地質變量構造深度構造深度砂體厚度砂體厚度有效厚度有效厚度

5、孔隙度孔隙度滲透率滲透率含油飽和度含油飽和度離散型地質變量離散型地質變量（范疇變量）（范疇變量）砂體砂體相相 流動單元流動單元隔夾層隔夾層斷層斷層類型變量類型變量8設設離散型隨機變量離散型隨機變量的所有可能取值為的所有可能取值為 x1，x2，其相應的概率為，其相應的概率為P (=xk)= pk, k=1,2,. 隨機變量的特征值：隨機變量的特征值：(1)數學期望數學期望 是隨機變量是隨機變量的整體代表性特征數。的整體代表性特征數。 則當級數 絕對收斂時，稱此級數的和為的數學期望，記為E()，或E。 E() = 1kkpxk1kkpxk9設連續型隨機變量的可能取值區間為(-，+)， p(x)為其

6、概率密度函數，若無窮積分 絕對收斂，則稱它為的數學期望，記為E()。 dxxxp)(E() = dxxxp)( 數學期望是隨機變量的最基本的數字特征， 相當于隨機變量以其取值概率為權的加權平均數。 從矩的角度說，數學期望是的一階原點矩。對于一組樣本：對于一組樣本：NzmNii)(110 為隨機變量的離散性特征數。若數學期望E-E()2存在，則稱它為的方差，記為D()，或Var()，或2。 = 222 )(E -)( )(E-E)(ED 從矩的角度說，方差是的二階中心矩。 (2)方差方差 其簡算公式為 D()=E(2) E()2D()= E-E()2方差的平方根為標準差，記為 11 研究范圍內的

7、一組隨機變量。研究范圍內的一組隨機變量。),(研究范圍uuZ)(uZ簡記為)( |)(,)(Pr)( |,;,(1111nzuZzuZobnzzuuFKKKK 隨機場：隨機場：當隨機函數依賴于多個當隨機函數依賴于多個自變量時，稱為隨機場。自變量時，稱為隨機場。如具有三個自變量如具有三個自變量(空間空間點的三個直角坐標點的三個直角坐標)的隨的隨機場機場2. 隨機函數隨機函數條件累積分布函數（ccdf）P12 二個隨機變量二個隨機變量，的協方差為二維隨機變量的協方差為二維隨機變量(，)的二階混合中心矩的二階混合中心矩11，記為，記為Cov(，)，或，或,。 協方差協方差(Variance)：Cov

8、(,) = , = E-E()-E()其簡算公式為其簡算公式為 Cov(,) = E ()-E() E() 隨機函數的特征值隨機函數的特征值 13P 任何統計推斷（cdf，數學期望等）均要求重復取樣。 但在儲層預測中，一個位置只能有一個樣品。 同一位置重復取樣，得到cdf，不現實14考慮鄰近點，推斷待估點 空間一點處的觀測值可解釋為一個隨機變量在該點 處的一個隨機實現。 空間各點處隨機變量的集合構成一個隨機函數。區域化變量： 能用其空間分布來表征一個自然現象的變量。（將空間位置作為隨機函數的自變量）（可以應用隨機函數理論解決插值和模擬問題）15考慮鄰近點，推斷待估點 -空間統計推斷要求平穩假設

9、),;,(),;,(1111KKKKzzhuhuFzzuuF 嚴格平穩嚴格平穩);();(zhuFzuF對于單變量而言：可從研究區內所有數據的累積直方圖推斷而得 （將鄰近點當成重復取樣點）太強的假設，不符合實際P16 當區域化變量Z(u)滿足下列二個條件時，則稱其為二階平穩或弱平穩： EZ(u) = EZ(u+h) = m(常數) xh 隨機函數在空間上的變化沒有明顯趨勢，隨機函數在空間上的變化沒有明顯趨勢，圍繞圍繞m值上下波動。值上下波動。 在整個研究區內有在整個研究區內有Z(u)的數學期望存在，的數學期望存在， 且等于常數，即：且等于常數，即：二階平穩二階平穩17 在整個研究區內，在整個研

10、究區內，Z(u)的協方差函數存在且平穩的協方差函數存在且平穩 (即只依賴于滯后即只依賴于滯后h，而與，而與u無關無關)， 即即 CovZ(u),Z(u+h) = EZ(u)Z(u+h)-EZ(u)EZ(u+h) = EZ(u)Z(u+h)- = C(h) 特殊地，當h=0時，上式變為VarZ(u)=C(0), 即方差存在且為常數。 協方差不依賴于空間絕對位置，而依賴于相對位置協方差不依賴于空間絕對位置，而依賴于相對位置 ， 即具有空間的平穩不變性。即具有空間的平穩不變性。uu+h18 在整個研究區內有在整個研究區內有 EZ(u)-Z(u+h) = 0 本征假設本征假設 當區域化變量Z(u)的增

11、量Z(u)-Z(u+h)滿足下列二條件時，稱其為滿足本征假設或內蘊假設。可出現EZ(u)不存在， 但EZ(u)-Z(u+h)存在并為零的情況存在并為零的情況 intrinsic hypotheseEZ(u)可以變化,但EZ(u)-Z(u+h)=0（比二階平穩更弱的平穩假設）19 增量增量Z(u)-Z(u+h)的方差函數的方差函數 (變差函數,Variogram) 存在且平穩存在且平穩 (即不依賴于即不依賴于u)，即：，即： VarZ(u)-Z(u+h) = EZ(u)-Z(u+h)2-EZ(u)-Z(u+h)2 = EZ(u)-Z(u+h)2 = 2(u,h) = 2(h), 相當于要求：相當

12、于要求：Z(u)的變差函數存在且平穩。的變差函數存在且平穩。20例：物理學上的著名的布朗運動是一種呈現出無限離散性的物理現象，其隨機函數的理論模型就是維納-勒維(Wiener-Levy)過程(或隨機游走過程)。布朗運動: 可出現協方差函數不存在，但變差函數存在的情況。 既不能確定驗前方差，也不能確定協方差函數。 但是其增量卻具有有限的方差： VarZ(x)-Z(x+h) = 2 = A|h| (其中，A是個常數)， 變差函數= |h|，且隨著|h|線性地增大。2A)(h21 若區域化變量若區域化變量Z(x)在整個區域內不滿足二階平在整個區域內不滿足二階平穩穩(或本征假設或本征假設) ，但在有限

13、大小的鄰域內是二階平，但在有限大小的鄰域內是二階平穩穩(或本征或本征)的，則稱的，則稱Z(x)是準二階平穩的是準二階平穩的(或準本征或準本征的的)。準二階平穩假設及準本征假設準二階平穩假設及準本征假設22 設 為區域上的一系列觀測點， 為相應的觀測值。區域化變量在 處的值 可采用一個線性組合來估計： nxx,1 nxzxz,10 x0*xzZ*(x0) niiixzxz10*min00*00*0 xZxZVarxZxZE無偏無偏最優最優無偏性和估計方差最小被作為 選取的標準 i-以普通克里金為例23從本征假設出發, 可知 為常數,有 xZE 0*11000mmxZxZExZxZEniiniii

14、可得到關系式: 11nii（1）無偏條件）無偏條件Z*(x0)（在搜尋鄰域內為常數，不同鄰域可以有差別）24njxZxZEnijj, 1, 021200*（2）估計方差最小）估計方差最小min200*200*00*xZxZExZxZExZxZE2k應用拉格朗日乘數法求條件極值Z*(x0)25niijniijinjxxCxxC1011, 1進一步推導，可得到n+1階的線性方程組, 即克里金方程組 當隨機函數不滿足二階平穩,而滿足內蘊（本征）假設時,可用變差函數來表示克里金方程組如下:niijniijinjxxxx1011, 1Z*(x0)26最小的估計方差,即克里金方差可用以下公式求解: nii

15、ikxxCxxC1000200102xxxxniiikZ*(x0)27 變差函數變差函數(或叫或叫變程方差函數變程方差函數，或，或變異函數變異函數)是是地質統計學所特有的基本工具。它既能描述區域化地質統計學所特有的基本工具。它既能描述區域化變量的空間結構性變化，又能描述其隨機性變化。變量的空間結構性變化，又能描述其隨機性變化。躍遷現象1. 變差函數的概念與參數變差函數的概念與參數 28),(hx 假設空間點假設空間點x只在一維的只在一維的x軸上變化，則將區域化軸上變化，則將區域化變量變量Z(x)在在x，x+h兩點處的兩點處的值之差值之差的方差之半定義的方差之半定義為為Z(x)在在x軸方向上的變

16、差函數，記為軸方向上的變差函數，記為一維情況下的定義：一維情況下的定義：VarZ(x)-Z(x+h) EZ(x)-Z(x+h)2-EZ(x)-Z(x+h)2 ),(hx21=21半變差函數(或半變異函數)29 在在二階平穩假設，或作本征假設二階平穩假設，或作本征假設，此時：，此時： 地質統計學中最常用的基本公式之一。EZ(x)-Z(x+h) = 0hVarZ(x)-Z(x+h) EZ(x)-Z(x+h)2-EZ(x)-Z(x+h)2 ),(hx21=21EZ(x)-Z(x+h)2 ),(hx21=則：30)()0()(hCCh（二階平穩假設條件下邊查函數與寫防查的關系）31變程變程(Range

17、) ：指區域化變量在空間上具有相關性的指區域化變量在空間上具有相關性的范圍。在變程范圍之內，數據具有相關性；而在變范圍。在變程范圍之內，數據具有相關性；而在變程之外，數據之間互不相關，即在變程以外的觀測程之外，數據之間互不相關，即在變程以外的觀測值不對估計結果產生影響。值不對估計結果產生影響。 32具不同變程具不同變程的克里金插的克里金插值圖象值圖象33塊金值塊金值(Nugget) ：變差函數如果在原點間斷，在地質統計學中稱：變差函數如果在原點間斷，在地質統計學中稱為為“塊金效應塊金效應”，表現為在很短的距離內有較大的空間變異性，表現為在很短的距離內有較大的空間變異性，無論無論h多小，兩個隨機

18、變量都不相關多小，兩個隨機變量都不相關 。它可以由測量誤差引起，。它可以由測量誤差引起，也可以來自礦化現象的微觀變異性。在數學上，塊金值也可以來自礦化現象的微觀變異性。在數學上，塊金值c0相當于相當于變量純隨機性的部分。變量純隨機性的部分。 34 如果品位完全是典型的隨機變量，則不論如果品位完全是典型的隨機變量，則不論觀測尺度大小，所得到的實驗變差函數曲線總觀測尺度大小，所得到的實驗變差函數曲線總是接近于純塊金效應模型。是接近于純塊金效應模型。 當采樣網格過大時，將掩蓋小尺度的結構，當采樣網格過大時，將掩蓋小尺度的結構，而將采樣尺度內的變化均視為塊金常數。這種而將采樣尺度內的變化均視為塊金常數

19、。這種現象即為塊金效應的尺度效應。現象即為塊金效應的尺度效應。塊金效應的尺度效應塊金效應的尺度效應12111333335基臺值基臺值(Sill)：代表變量在空間上的總變異性大小。即為變代表變量在空間上的總變異性大小。即為變差函數在差函數在h大于變程時的值，為大于變程時的值，為塊金值塊金值c0和和拱高拱高cc之和。之和。 拱高拱高為在取得有效數據的尺度上，可觀測得到的變異性幅為在取得有效數據的尺度上，可觀測得到的變異性幅度大小。當塊金值等于度大小。當塊金值等于0時，基臺值即為拱高。時，基臺值即為拱高。 = C(0) C(h)(h36幾何各向異性：幾何各向異性：變差函數變差函數在空間各個方向上的在

20、空間各個方向上的變程變程不同不同，但，但基臺值不變基臺值不變（即（即變化程度相等）。這種情變化程度相等）。這種情況能用一個簡單的幾何坐況能用一個簡單的幾何坐標變換將各向異性結構變標變換將各向異性結構變換為各向同性結構。換為各向同性結構。帶狀各向異性：帶狀各向異性：不同方向不同方向的變差函數具有的變差函數具有不同的基不同的基臺值臺值，其中，其中變程可以不同，變程可以不同，也可以相同也可以相同。這種情況不。這種情況不能通過坐標的線性變換轉能通過坐標的線性變換轉化為各向同性，因而結構化為各向同性，因而結構套合是比較復雜的。套合是比較復雜的。 地質變量相關性的各向異性地質變量相關性的各向異性12111

21、3333(2)372. 變差函數的理論模型變差函數的理論模型設Z(x)為滿足本征假設的區域化變量，則常見的理論變差函數有以下幾類：球狀模型球狀模型指數模型指數模型高斯模型高斯模型冪函數模型冪函數模型空洞效應模型空洞效應模型38 接近原點處，變差函接近原點處，變差函 數呈線性形狀，在變數呈線性形狀，在變 程處達到基臺值。程處達到基臺值。 原點處變差函數的切原點處變差函數的切 線在變程的線在變程的2/3處與處與 基臺值相交。基臺值相交。 ahcahahahchahSphch,2123003球狀模型：球狀模型： c為基臺值，為基臺值，a為變程，為變程，h為滯后距。為滯后距。39指數模型：指數模型：

22、ahcahExpch3exp1 變差函數漸近地逼近變差函數漸近地逼近 基臺值。基臺值。 在實際變程處，變差在實際變程處，變差 函數為函數為0.95c。 模型在原點處為直線。模型在原點處為直線。40高斯模型：高斯模型： 223exp1ahch 變差函數漸近地逼近變差函數漸近地逼近 基臺值。基臺值。 在實際變程處，變差函在實際變程處，變差函 數為數為0.95c。 模型在原點處為拋物線。模型在原點處為拋物線。 41冪函數模型：冪函數模型： hch. 冪函數模型為一種無基冪函數模型為一種無基臺值的變差函數模型。這臺值的變差函數模型。這是一種特殊的模型。是一種特殊的模型。 當當 =1時，變差函數為一時，

23、變差函數為一直線，即為線性模型，這直線，即為線性模型，這一模型即為著名的一模型即為著名的布朗運布朗運動（隨機行走過程）動（隨機行走過程）的變的變差函數模型；差函數模型； 當當 1時，變差函數為拋時，變差函數為拋物線形狀，為物線形狀，為分數布朗運分數布朗運動動(fBm)的變差函數模型。的變差函數模型。 布朗運動布朗運動分數布朗運動分數布朗運動分數布朗運動分數布朗運動 h2111h42空洞效應模型空洞效應模型(Hole Effect)： 2cosexp1.bhahch 變差函數并非單調增加，變差函數并非單調增加， 而顯示出一定周期性的而顯示出一定周期性的 波動。波動。 模型可以有基臺值，也模型可以

24、有基臺值，也 可以無基臺值；可以有可以無基臺值；可以有 塊金值，也可以無塊金塊金值，也可以無塊金 值。值。 空洞效應在地質上多沿空洞效應在地質上多沿 垂向上出現，如富礦層垂向上出現，如富礦層 與貧礦層互層、砂巖與與貧礦層互層、砂巖與 泥巖頻繁薄互層等等。泥巖頻繁薄互層等等。（b為富礦化帶重復距離）)(hh43 通過區域化變量的空間觀測值來通過區域化變量的空間觀測值來構建相應的變構建相應的變差函數模型差函數模型, 以表征該變量的主要結構特征。以表征該變量的主要結構特征。(求變求變差差) (1)數據準備數據準備 區域化變量的選取區域化變量的選取、 數據質量檢查及校正數據質量檢查及校正、 數據的變換

25、數據的變換（如對滲透率進行對數變換）、（如對滲透率進行對數變換）、 數據的統計數據的統計（如分相對儲層參數計算平均值、（如分相對儲層參數計算平均值、 方差，作直方圖、相關散點圖等）、方差，作直方圖、相關散點圖等）、 叢聚數據的解串叢聚數據的解串等。等。3. 區域化變量的區域化變量的44(2)(2)實驗變差函數的計算實驗變差函數的計算 實驗變差函數是指應用觀測值計算的變差函實驗變差函數是指應用觀測值計算的變差函數。對于不同的滯后距數。對于不同的滯后距h h，可算出相應的實驗變，可算出相應的實驗變差函數差函數。 )(*h=N(h)1i2iih)Z(x-)Z(xN(h)21一維實驗變差函數的計算公式

26、(i=1,N(h) Z(xi)-Z(xi+h)2的算術平均值一半即為一個h的變差函數值45對不同的滯后h，進行計算，得出各個h的變差函數值 )(*h=N(h)1i2iih)Z(x-)Z(xN(h)21h3h5hh46設Z(x)為一維區域化變量，滿足本征假設，又已知Z(1)=2，Z(2)=4，Z(3)=3，Z(4)=1，Z(5)=5，Z(6)=3，Z(7)=6，Z(8)=4， ， ) 1 (*)2(*)3(*例：例：試求：試求：)(*h=N(h)1i2iih)Z(x-)Z(xN(h)2147721) 1 (*)2(*)3(*=22+12+22+42+22+32+22 = 1442= 3.0062

27、112+32+22+22+12+12 = 1220= 1.6752112+12+02+52+12 = 1028= 2.80482D情況情況（1）分不同方向，進行1D變差函數計算3D情況情況： 增加垂向方向（2）確定主變程方向 次變程方向角度容限步長容限h3h5hh四方向試算（考慮主變程方向的 走向、傾向和傾角）49(3)(3)理論變差函數的最優擬合與結構套合理論變差函數的最優擬合與結構套合 選擇合適的理論變差函數模型，同時還需進選擇合適的理論變差函數模型，同時還需進行結構套合，從而得到一條反映不同層次（或行結構套合，從而得到一條反映不同層次（或不同空間規模）結構的、統一的、最優的變差不同空間規

28、模）結構的、統一的、最優的變差函數曲線。函數曲線。球狀模型球狀模型指數模型指數模型高斯模型高斯模型冪函數模型冪函數模型空洞效應模型空洞效應模型50 復雜的區域化變量往往包含各種尺度上的多層次、多方向的變化性，反映在變差函數上即為多層次結構。將不同結構組合為統一結構的過程稱為“結構套合”結構套合結構套合各層次套合各層次套合例如，對于200米寬的河道，在h50m的觀測尺度上可以將其與河道間的變化性區分出來，但卻無法區分層理和礦物成分的變化性 (即無法找出更細微的結構來)，它們在50m尺度得到的結構上只能作為“塊金效應”出現。若觀測尺度為500米，河道的變化也只能作為“塊金效應”。 12111333

29、3大尺度的變化性總是包含著小尺度的變化性，但卻不能從大尺度的變化性中區分出小尺度的變化性。51)()()()(210rrrr)(0r=。0, 0, 00rCr)(1r= 1 1131311 ar ,Car0 ),2123(ararC)(2r 2 2232322 ar ,Car0 ),r21-r23(aaC=代表微觀變化性的變程極小的球狀模型，可近似地看作純塊金效應型 球狀模型，沒有塊金常數，基臺值為C1，變程為a1 ，反映了小規模范圍的變化 球狀模型，沒有塊金常數，基臺值為C2，變程較大，為a2 ，反映了大規模范圍的變化 可以用反映各種不同尺度變化性的多個變差函數之和來表示一個套合結構。（各層

30、次理論模型可以不一樣）)(ri可以是不同模型的變差函數52其中 21aa 則套合結構的表達式為 )(r2210213232210133221122110a r ,CCCa ),2123(0 ,r )(21)(230, 03raararCCCaraCaCraCaCCr=。0, 0, 00rCr 1 1131311 ar ,Car0 ),2123(ararC 2 2232322 ar ,Car0 ),r21-r23(aaC)(0r)(1r)(2r=53對于幾何各向異性，先根據異向比壓縮 距離軸，使之成為各向同性的模型；對于帶狀各向異性，運用模型疊加的方法加以處理。先用壓縮距離軸的辦法，使其變程變為

31、相同，然后再把具有相同變程的兩個球狀模型疊加起來，構成一個新的球狀模型 各方向套合各方向套合（將各向異性套合為各向同性，以便于 在克里金估計時，不同方向均可用統一 的結構模型計算實際的變差函數值）54(4)(4)變差函數參數的最優性檢驗：變差函數參數的最優性檢驗： 變差函數是否符合實際，應該進行檢驗。變差函數是否符合實際，應該進行檢驗。一種實用的檢驗方法為一種實用的檢驗方法為“交叉驗證法交叉驗證法”（Cross-validationCross-validation），檢驗標準是在各實測），檢驗標準是在各實測點，根據周圍點計算的點，根據周圍點計算的克里金估計值與該實測克里金估計值與該實測值的誤差

32、平方值的誤差平方平均最小。平均最小。 估計誤差的平方估計誤差的平方與與克里金估計方差克里金估計方差之比越接之比越接近近1 1，則說明變差函數與實際的符合程度越高。，則說明變差函數與實際的符合程度越高。實際上，這種方法在檢驗變差函數的同時，也實際上，這種方法在檢驗變差函數的同時，也在檢驗所使用的克里金估計方法的適用性。在檢驗所使用的克里金估計方法的適用性。Z*(x0)55niijniijinjxxCxxC1011, 1（以普通克里金為例）i求取變差函數（或協方差）；求取變差函數（或協方差）；解克里金方程組解克里金方程組56 設有一個油藏，在平面上S1，S2，S3，S4處有四個井點，其孔隙度值分別

33、為Z1，Z2，Z3，Z4。據此估計S0點處的孔隙度值Z0 設孔隙度Z(x)是二階平穩的。其在平面上的二維變差函數是一個各向同性的球狀模型，其參數為：塊金值C02，變程a200，拱高C20，即： 實例實例200,h 22,200,h0 ),(200)h21-200h2320(20,h 0,(h)3357Z0的估計量為 41iii*0ZZ普通克里金方程組的矩陣形式為 K =M2 0 1 1 1 1 1 C C C C 1 C C C C1 C C C C1 C C C C K, 444342413433323124232221141312114321 ,1 CCCCM04030201221MKni

34、ijniijinjxxCxxC1011, 1(求解）58 0 1 1 1 1 1 C C C C 1 C C C C1 C C C C1 C C C C K, 444342413433323124232221141312114321)()0()(hChC求解求解: Cij ,1 CCCCM040302012C11C12C01試求200,h 22,200,h0 ),(200)h21-200h2320(20,h 0,(h)33?59 C11 = C22 = C33 = C44 = C(0) =2 = C0+C = 22,由于C(h)=C(0) - (h)=22 - (h)當ij時，Cij=C(|S

35、i-Sj|)=22 - (|Si-Sj|).于是，C12=C21=C04=22- )250(,84. 9)200250(21)20025023(202223,22. 1)50150(22223113CC,98. 4)50100(2222024114CCC,32. 2)100100(22223223CC,28. 0)100150(22224224CC0)50200(22224334CC,66.12)50(2201C,72. 1)150(2203C 60將以上數值代入普通克里金方程組解的矩陣形式中，得 19.841.724.9812.66 0 1 1 1 1 1 22 0 0.28 4.98 1

36、0 22 2.32 1.221 0.28 2.32 22 9.841 4.98 1.22 9.84 22 14321經計算得: =0.5182, =0.0220, =0.0886, =0.3712。 1234Z00.5182Z1+0.0220Z2+0.0886Z3+0.3712Z421MK61 搜索鄰域搜索鄰域注意注意1： 搜索鄰域中的數據點才參加估計節省CPU和內存局域平穩 搜索橢圓或橢球的選擇方法與選擇變差函數橢圓或橢球相同。62注意注意2： 參與計算的數據點不能太多，否則計算太慢 一般軟件中都內置或可選最大的 數據點數目（與待估點最近的數據點），如10。注意注意3： 防止數據叢聚帶來的數

37、據代表性不強井眼井眼垂向數據太密，若待估點與該井近，則可能忽視鄰井數據八分搜尋，保證各象限均有代表數據63 若搜尋范圍無數據，則應用邊際概率。若搜尋范圍無數據，則應用邊際概率。64x0 簡單克里金簡單克里金(SK) (SK) 普通克里金普通克里金(OK) (OK) 泛克里金泛克里金(UK) (UK) 協同克里金協同克里金(CK) (CK) 貝葉斯克里金（貝葉斯克里金（BKBK） 指示克里金指示克里金(IK) (IK) 65所有克里金估計都應用線性回歸算法，形式為：m為期望)()()(1*umuZumZnSK求取權系數的克里金方程組的非平穩形式nnuuCuuCu1), 2 , 1( ),(),(

38、)(求(n+1)個m(u), 求(n+1)(n+1) 個C(u,u)66二階平穩假設二階平穩假設EZ(u) = EZ(u+h) = m(常數)C(u,u+h) = C(h) nnSKumuZuuZ11*)(1)()()(簡單克里金估計的平穩形式：)()()(1*umuZumZnSKEZ(u) = EZ(u+h) = m(常數)67應用條件：應用條件： 隨機函數二階平穩隨機函數二階平穩 隨機函數的期望值 m為常數并已知已知不能用于具有局部趨勢的情況nnuuCuuCu1), 2 , 1( ),()()(簡單克里金方程組的平穩形式：nnuuCuuCu1), 2 , 1( ),(),()(C(u,u+

39、h) = C(h) （C與位置有關）（C與位置無關）68 nuZuuZ1*)(nnunuuCuuuCu111)(, 1)()(69應用要求：應用要求： 隨機函數二階平穩或符合內蘊假設隨機函數二階平穩或符合內蘊假設 隨機函數的期望值 m在搜尋鄰域內穩定但未知 協方差平穩 與簡單克里金相比，普通克里金相當于在每一與簡單克里金相比，普通克里金相當于在每一個位置個位置u，重新估計，重新估計 m。 由于普通克里金估計常使用滑動數據鄰域，由于普通克里金估計常使用滑動數據鄰域，相當于均值相當于均值m隨位置可變，即隨位置可變，即Z*(u)，此時，實際，此時，實際上是一種非平穩算法，對應于變化的均值和平穩上是一

40、種非平穩算法，對應于變化的均值和平穩的協方差。的協方差。70 u uu umZE非平穩隨機函數的漂移函數非平穩隨機函數的漂移函數(drift),簡稱為漂移或趨勢簡稱為漂移或趨勢 )()()(uuuRmZ隨機函數隨機函數 = 趨勢趨勢 + 殘差殘差區域化變量Z(X)是非平穩的，即EZ(x)=m(x) Kriging with a trend model (KT)具有趨勢的克里金具有趨勢的克里金71ma fkkkK( )( )uu0用光滑的確定性函數來模擬，或用擬合方法 趨勢函數趨勢函數一維的線性趨勢 maa x( )u 01二維的二次趨勢: maa xa ya xa ya xy( )u 0123

41、242572R( )u用均值為0、協方差函數為 的平穩隨機函數來模擬。CR( )hZZKTKTn*()( )( ) ()uuu1泛克里金估計值泛克里金估計值: 殘差殘差73()()( )()( )()(), ,( )()( ), ,KTRnkkkKRKTkknCfCnffkKuuuuuuuuuu1011201 為權值 ()KTk( )u是與(K+1)個權值的限制條件相對應的(K+1)個拉格朗日參數泛克里金方程組泛克里金方程組 CR( )h為殘差協方差函數 74)()()(10uuuYaamZEZZKTKTn*()( )( ) ()uuu1估計值估計值 當K = 1時，線性趨勢函數為ma fkk

42、kK( )( )uu0趨勢函數可理解為二級變量75 （1）外部變量必須在空間光滑）外部變量必須在空間光滑 地變化，否則可能導致地變化，否則可能導致KT 線性系統不穩定；線性系統不穩定； （2）在主變量的所有數據點）在主變量的所有數據點u 處和待估計的處和待估計的 位置位置u處，外部變量都必須是已知的。處，外部變量都必須是已知的。 nKTnKTRnRKTYYnCYC1)(1)(101)()()()(1)(, 1)()()()()()(uuuuuuuuuuuu克里金方程組： 可理解為地震 數據（如深度）（K=0時，？）76 利用幾個變量之間的空間相關性，對其中的一個或幾利用幾個變量之間的空間相關性

43、，對其中的一個或幾個變量進行空間估計，尤其適用于被估計變量的觀察數據個變量進行空間估計，尤其適用于被估計變量的觀察數據較少的情況較少的情況 。 mjjjniiiyxZ110*協同克里金估計值（初始變量和二級變量）協同克里金估計值（初始變量和二級變量） -隨機變量在位置0處的估計值；-初始變量的n個樣本數據；-二級變量的m個樣本數據；-需要確定的協同克里金加權系數。0*Znxx,1myy,1naa,1及m,17701, 2 , 1, 2 , 1,1102110111mjjniijmjjiinijiijmjjiinijiimjyxCyyCyxCnixxCxyCxxC協同克里金方程組協同克里金方程組

44、傳統普通協克里金傳統普通協克里金78標準化普通協克里金標準化普通協克里金)(110*YXmmyxZmjjjniii111mjjniimX = Ex(u)mY = Ey(u)79 為協同克里金的簡化形式，即如果二級變量密為協同克里金的簡化形式，即如果二級變量密集取樣時，只保留與估計點同位的二級變量。集取樣時，只保留與估計點同位的二級變量。)()()()()(1uYuuZuuZjniii對應的協同克里金方程組只要求知道Z - 協方差函數以及Z-Y 互協方差函數CZ（h）)()(hChCZZY)0(/ )0()0(ZYZYCCP（同位兩種數據的 相關系數）（方差函數）同位協同克里金同位協同克里金 C

45、ollocated Cokriging80H.Omre在（1987）把線性貝葉斯理論用于克里金估計技術，提出了貝葉斯克里金估計技術。他構想了一個模型，把用于空間估計的數據分為兩類： 觀察數據：是指那些精度比較高，但數量比較少的數據 猜測數據：是指那些精度比較低，但分布廣泛的數據 在觀測數據比較多的地方，估計結果主要受觀測數據的影響；在觀測數據比較少的地方，則主要受猜測數據的影響。 顯然，井數據和地震數據的關系符合貝葉斯估計中觀測數據和猜測數據的關系。 81 設Z(x)，xA，是觀察數據的區域化變量。 設M(x)，xA，是猜測數據的區域化變量。 Z*(x0) = EZ(x0) = a0+M(x0

46、) EM(x) M(x) ，xAM(x)是對Z(x)的一種猜測，誤差為a0 x0a0 ？82 設已得到設已得到Z(x)，xA的一組的一組(N個個)觀察值觀察值 Z(xi); i=1,2,N。 定義一個新的隨機函數：定義一個新的隨機函數： ZT(x)Z(x)-M(x)，xAZT(xi) = Z(xi) -M(x)， i=1,2,NZ(x0)的貝葉斯克里金估計量為的貝葉斯克里金估計量為NiMiiBKxxZxZxZT100*0*)()()()(x0對這個N個觀察值有（相當于誤差a0 ）（誤差的隨機函數）83基于無偏性和估計方差最小兩個條件：基于無偏性和估計方差最小兩個條件：min00*00*0 xZ

47、xZVarxZxZE利用拉格朗日乘數法可得到貝葉斯克里金方程組：利用拉格朗日乘數法可得到貝葉斯克里金方程組： 1,00|1|jjiMiMZjjiMjiMZjaxxxxxxxxaNj, 2 , 1 Z(x,x) = ZM(x-x)+ M(x,x)84將數據按照不同的門檻值編碼為將數據按照不同的門檻值編碼為1或或0的過程。的過程。 對于模擬目標區內的對于模擬目標區內的每一類相，當它出現于某每一類相，當它出現于某一位置時，指示變量為一位置時，指示變量為1，否則為否則為0。A (100) B (010) A (100) C(001) 類型變量的指示變換：類型變量的指示變換： 01ui 變量 u 屬于范

48、疇A 其它指示變換指示變換1982年由AGJournel(儒爾奈耳)教授提出 85(00111) (00001) (01111) (00011) 首先將連續變量截斷首先將連續變量截斷為類型變量，然后進為類型變量，然后進行指示變換。行指示變換。如： z = 10， 15， 20， 25, 30 zxzzxzzui01; 連續變量的指示變換連續變量的指示變換 86設沿空間某一方向，在間距為h的5對樣品點處觀測了Z(x)及Z(x+h)的值 (=1,2，5)。 1 2 3 4 5 Z(x) 0 2 3 6 9 Z(x+h) 1 6 7 8 8 設指示XI(x; z) = 設指示Y=I(x+h; z)

49、= z)Z(, 0z)Z(, 1當當zh)Z(x, 0zh)Z(x, 1當當87假定只選定了5個門限值：0，2，3，6，9 XI(x;z) Y=I(x+h;z) 當 z= 當 z= 0 2 3 6 9 0 2 3 6 9 1 1 1 1 1 1 0 1 1 1 1 2 0 1 1 1 1 0 0 0 1 1 3 0 0 1 1 1 0 0 0 0 1 4 0 0 0 1 1 0 0 0 0 1 5 0 0 0 0 1 0 0 0 0 1 1 2 3 4 5 Z(x) 0 2 3 6 9 Z(x+h) 1 6 7 8 8 z)Z(, 0z)Z(, 1當當zh)Z(x, 0zh)Z(x, 1當當8

50、8 指示指示(函數函數)的數學期望的數學期望 當x固定時，若z再給定，則I(x；z)就是個隨機變量，就有數學期望： EI(x；z)1PI(x；z)=1+0PI(x；z)=0 = PI(x；z)=1=P =F(x；z)， zZ(x) z+ 在x點處區域化變量Z(x)的先驗分布函數F(x；z)，z+就是x點處指示函數的數學期望 EI(x；z)i*(x；z) = ),();(11zFzxinAAAnn89(1) z=3時，設指示隨機變量X I (x;3)E(X)=EI(x;3) = 51)(6 . 053)3 ;(51xmxiVar(X)= Var I(x;3)=mx(1-mx)=0.60.4=0.

51、24= 2xXI(x;z) Y=I(x+h;z) 當 z= 當 z= 0 2 3 6 9 0 2 3 6 9 1 1 1 1 1 1 0 1 1 1 1 2 0 1 1 1 1 0 0 0 1 1 3 0 0 1 1 1 0 0 0 0 1 4 0 0 0 1 1 0 0 0 0 1 5 0 0 0 0 1 0 0 0 0 1 試求指示隨機變量的數學期望和方差：試求指示隨機變量的數學期望和方差：90(2) z=6時，設指示隨機變量YI (x+h;6)E(Y)= EI(x+h;6)= 51)(4 . 052)6 ;(51YhmxiVar(Y)= Var I(x+h;6)= mY(1-mY)=0.

52、40.6=0.24= 2YXI(x;z) Y=I(x+h;z) 當 z= 當 z= 0 2 3 6 9 0 2 3 6 9 1 1 1 1 1 1 0 1 1 1 1 2 0 1 1 1 1 0 0 0 1 1 3 0 0 1 1 1 0 0 0 0 1 4 0 0 0 1 1 0 0 0 0 1 5 0 0 0 0 1 0 0 0 0 1 91設設Z(X)是個一維區域化變量，在等間距的是個一維區域化變量，在等間距的10個點個點 處有處有10個個觀測值：觀測值：Z(1)=3，Z(2)=5， Z(3)=6，Z(4)=2，Z(5)=7，Z(6)=1，Z(7)=4，Z(8)=8，Z(9)=9，Z(1

53、0)=7，設門限值，設門限值Z分別等于分別等于2 , 3, 4, 5, 6, 7, 8時，求指示變差函數的估計值時，求指示變差函數的估計值 (h;z)，h=1，2，3，4，5。計算計算*I92)5 ;()5 ;(212hxIxIEI(h;5)= *I)(12)5 ;()5 ;()(21hnhxixihn(h;5) = 首先，計算i(x;5)：i(x1;5)1,i(x2;5)1,i(x3;5)0,i(x4;5)=1, i(x5;5)0, i(x6;5)=1, i(x7;5)1, i(x8;5)0,i(x9;5)0,i(x10;5)0,*I;2778. 0185)001011110(9212222

54、22222(1;5) = zxzzxzzxia01;93 1 2 3 4 5 6 7 8 9 10 i(x;2) i(x;3) i(x;4) i(x;5) i(x;6) i(x;7) i(x;8) h z 2 3 4 5 6 7 8 1 2 3 4 5 列表計算各指示值i(x;z) 根據i(x;z)值算出不同z值的 (h;z)值 *IZ(X)356271489794h z 2 3 4 5 6 7 8 1 0.2222 0.2778 0.2778 0.2778 0.1667 0.1111 0.1111 2 0.1250 0.1875 0.3125 0.2500 0.2500 0.1875 0.0

55、625 3 0.2857 0.2143 0.2143 0.2857 0.2143 0.1429 0.0714 4 0.2500 0.3333 0.4167 0.3333 0.2500 0.1667 0.0833 5 0.2000 0.1000 0.2000 0.1000 0.2000 0.2000 0.1000 根據i(x;z)值算出的不同z值的 (h;z)值 *I 計算的各指示值i(x;z) 1 2 3 4 5 6 7 8 9 10 i(x;2) 0 0 0 1 0 1 0 0 0 0 i(x;3) 1 0 0 1 0 1 0 0 0 0 i(x;4) 1 0 0 1 0 1 1 0 0 0

56、 i(x;5) 1 1 0 1 0 1 1 0 0 0 i(x;6) 1 1 1 1 0 1 1 0 0 0 i(x;7) 1 1 1 1 1 1 1 0 0 1 i(x;8) 1 1 1 1 1 1 1 1 0 1 95指示克里金指示克里金 線性估計量線性估計量 i*(x;z) = nzxizx1);();(其中，x(=1,2,，n)點處的Z(x)值已知， (x;z) (=1,2,，n)為IK權系數 (00111) (00001) (01111) (00011) （普通指示克里金）96nnIIzxnzxxCzzxxCzx111);(.,2 , 1),;()();();(nnIIzxnzxxzzxxzx111);(.,2 , 1),;()();();(CI(h;z)為指示協方差， 為指示變差函數 );(zhI是拉格朗日乘數 )(z有多少個門限值z，就有多少個IK方程組 IK方程組：方程組： 97從IK方程組中解出 ), 2 , 1)(;(nzxi*(x;z) = 求出i(x;z)的估計值 nzxizx1);();(i*(x;z)實際上是I(x;z)的條件數學期望EI(x;z)|(n)的估計值。 EI(x;z)|(n) = 0PI(x;z)=0|(n)+1PI(x;z)=1|(n) = PI(x;z)=1|(n) = PZ(x） z)|(n) = Fx;z|(n)