地質(zhì)統(tǒng)計(jì)分析方法2012

1、地質(zhì)統(tǒng)計(jì)分析方法Geostatistic Analysis 經(jīng)典統(tǒng)計(jì)學(xué)在應(yīng)用于地質(zhì)變量時(shí)存在不少缺陷，這是帶有普遍性和根本性的。主要有三：在統(tǒng)計(jì)樣品品位的頻率和做直方圖時(shí)均不考慮樣品的空間分布；經(jīng)典概率統(tǒng)計(jì)學(xué)的研究對象必須是純隨機(jī)變量（獨(dú)立性），許多地質(zhì)變量并不是純隨機(jī)變量（空間相關(guān)性），而是既有隨機(jī)性有又結(jié)構(gòu)性（在空間分布上有某種程度的相關(guān)性或連續(xù)性）的變量。經(jīng)典概率統(tǒng)計(jì)學(xué)所研究的變量原則上都是可以無限次重復(fù)試驗(yàn)或大量觀測的，但地質(zhì)變量不行。因?yàn)橐坏┰诘V體某處取一樣品后，嚴(yán)格來說，就不可能在同一地方再次取道樣品了。同一批品位數(shù)據(jù)的空間變化性（它們的平均品位和方差都一樣） 為了解決在地質(zhì)變量具

2、有隨機(jī)性和結(jié)構(gòu)性的條件下仍能使用統(tǒng)計(jì)方法的問題，20世紀(jì)40年代末出現(xiàn)了變差函數(shù)，或稱變差圖(Variogram)，它能夠同時(shí)描述地質(zhì)變量的隨機(jī)性和結(jié)構(gòu)性變化。地質(zhì)統(tǒng)計(jì)學(xué)的誕生 從1951年起，南非的礦山地質(zhì)工程師D.G.克立格和統(tǒng)計(jì)學(xué)家H.S.西舍爾等人根據(jù)他們對南非金礦多年來工作的經(jīng)驗(yàn)，提出了根據(jù)樣品空間位置不同和樣品間相關(guān)程度的不同，對每個(gè)樣品品位賦予一定的權(quán)，進(jìn)行滑動(dòng)加權(quán)平均，來估計(jì)中心塊段平均品位的方法，這就是克立格法（也稱“克立金”kriging）。地質(zhì)統(tǒng)計(jì)學(xué)（Geostatistics）是針對像礦產(chǎn)、資源、生物群落、地貌等有著特定的地域分布特征而發(fā)展的一門新的統(tǒng)計(jì)學(xué)分支。 其理論

3、基礎(chǔ)是由法國著名統(tǒng)計(jì)學(xué)家G. Matheron于20世紀(jì)60年代在大量理論研究的基礎(chǔ)上逐漸形成的。由于最先在地學(xué)領(lǐng)域應(yīng)用，故又稱地統(tǒng)計(jì)學(xué)。 目前，地質(zhì)統(tǒng)計(jì)學(xué)已經(jīng)被廣泛用于地理學(xué)、生態(tài)學(xué)、環(huán)境科學(xué)、土壤學(xué)等諸多領(lǐng)域的研究中。特別是GIS的發(fā)展帶來的空間數(shù)據(jù)的極大豐富，越來越多的科學(xué)家求助于地統(tǒng)計(jì)學(xué)來分析空間數(shù)據(jù)。 地質(zhì)統(tǒng)計(jì)學(xué)與經(jīng)典統(tǒng)計(jì)學(xué)的共同之處在于：它們都是在大量采樣的基礎(chǔ)上，通過對樣本屬性值的頻率分布或均值、方差關(guān)系及其相應(yīng)規(guī)則的分析，確定其空間分布格局與相關(guān)關(guān)系。地質(zhì)統(tǒng)計(jì)學(xué)區(qū)別于經(jīng)典統(tǒng)計(jì)學(xué)的最大特點(diǎn)：地統(tǒng)計(jì)學(xué)既考慮到樣本值的大小，又重視樣本空間位置及樣本間的距離，彌補(bǔ)了經(jīng)典統(tǒng)計(jì)學(xué)忽略空間方

4、位的缺陷。 地質(zhì)統(tǒng)計(jì)學(xué)是以區(qū)域化變量理論為基礎(chǔ)，以變異函數(shù)為主要工具，研究那些在空間分布上既有隨機(jī)性又有結(jié)構(gòu)性，或具有空間相關(guān)和依賴性的自然現(xiàn)象的科學(xué)。 克立格法是地質(zhì)統(tǒng)計(jì)分析的主要方法之一。Contents 區(qū)域化變量 協(xié)方差函數(shù) 變異函數(shù) 克立格法簡介 當(dāng)一個(gè)變量具有空間分布特征時(shí)，就稱之為區(qū)域化變量（Regionalized Variable）。一、區(qū)域化變量 1. 定義 區(qū)域化變量，亦稱區(qū)域化隨機(jī)變量，G. Matheron（1963）將它定義為以空間點(diǎn)x的三個(gè)直角坐標(biāo)為自變量的隨機(jī)場 。這種變量常常反映某種空間現(xiàn)象的特征，用區(qū)域化變量來描述的現(xiàn)象稱之為區(qū)域化現(xiàn)象-礦產(chǎn)、地質(zhì)、海洋、土

5、壤、氣象、水文、生態(tài)、溫度、濃度等領(lǐng)域都具有某種空間屬性。 具有局部的、隨機(jī)的、異常的性質(zhì)； 具有一般的或平均的結(jié)構(gòu)性質(zhì)。 具有空間的局限性、不同程度的連續(xù)性和不同程度的各向異性等特征。 空間協(xié)方差函數(shù)和變異函數(shù)能反映區(qū)域化變量的隨機(jī)性和結(jié)構(gòu)性。 一、區(qū)域化變量 2. 區(qū)域化變量的性質(zhì)變量在點(diǎn)X與偏離空間距離為h的點(diǎn)Xh處的隨機(jī)量Z(X)與Z(X+h)具有某種程度的自相關(guān)，而且這種自相關(guān)性依賴于兩點(diǎn)間的距離h與變量特征。 區(qū)域化隨機(jī)變量之間的差異，定義為區(qū)域化變量Z(x)的自協(xié)方差函數(shù)。 它是 在空間點(diǎn)x和x+h處的兩個(gè)隨機(jī)變量的二階混合中心矩：二、協(xié)方差函數(shù) 1.協(xié)方差函數(shù)的概念二階混合中心

6、矩是關(guān)于隨機(jī)向量兩個(gè)分量之間關(guān)聯(lián)程度的一個(gè)特征設(shè) 為一個(gè)n 維隨機(jī)向量, 稱 為xi 和xj 的二階混合中心矩; 其中E 為數(shù)學(xué)期望符號.X=（x1，x2，x3xn）二、協(xié)方差函數(shù)2.協(xié)方差函數(shù)的計(jì)算公式 式中：h為兩樣本點(diǎn)空間分隔距離或距離滯后； 為 在空間位置 處的實(shí)測值； 是 在 處距離偏離h的實(shí)測值； 是分隔距離為h時(shí)的樣本點(diǎn)對總數(shù)； 和 分別為 和 的樣本平均數(shù)。若 = =m（常數(shù)），則上式可以改寫為：式中：m為樣本平均數(shù)，可由一般算術(shù)平均數(shù)公式求得，即： 1.變異函數(shù)的概念 變異函數(shù)(Variograms)，又稱變差函數(shù)、變異矩，是地質(zhì)統(tǒng)計(jì)分析所特有的基本工具。三、變異函數(shù) 在一維

7、條件下變異函數(shù)定義為：當(dāng)空間點(diǎn)Z在一維x軸上變化時(shí)，區(qū)域化變量Z(x)在點(diǎn)x和x+h處的值Z(x)與Z(x+h)差的方差的一半為變量Z(x)在x軸方向上的變異函數(shù)，記為(h)，即：在二階平穩(wěn)假設(shè)條件下，對任意的h有因此，公式可以改寫為 從上式可知，變異函數(shù)依賴于兩個(gè)自變量x和h，當(dāng)變異函數(shù) 僅僅依賴于距離h而與位置x無關(guān)時(shí)， 可改寫成 ，即： 平穩(wěn)假設(shè)是指區(qū)域化變量的分布不因位置而改變, 對于二階平穩(wěn), 滿足以下兩個(gè)條件: 整個(gè)研究區(qū)內(nèi)區(qū)域化變量Z(x)存在且不取決于x; 整個(gè)研究區(qū)內(nèi)空間變量的協(xié)方差存在且相同, 即:C (x ) = EZ (x ) - m Z (x + h) - m = E

8、 Z (x ) Z (x + h) - m2 (1)(0)=0，即在h=0處，變異函數(shù)為0； (2)(h)=(-h)，即(h)關(guān)于直線h=0是對稱的，它是一個(gè)偶函數(shù)； (3)(h)0，(h)表示的方差只能大于或等于0；三、變異函數(shù) 2.變異函數(shù)的性質(zhì)設(shè)Z(x)是區(qū)域化變量，在滿足二階平穩(wěn)假設(shè)條件下，變異函數(shù)式具有如下性質(zhì)：(4)|h|時(shí)，(h)c(0)，或() =c(0)，即當(dāng)空間距離增大時(shí)，變異函數(shù)接近先驗(yàn)方差： 3.變異函數(shù)的計(jì)算公式 設(shè)Z(x)為一區(qū)域化隨機(jī)變量，并滿足二階平穩(wěn)假設(shè)，Z(x)是系統(tǒng)某屬性Z在空間位置x處的值，h為兩樣本點(diǎn)空間分隔距離， Z(xi) 和Z(xi+h) 分別是

9、區(qū)域化變量Z(x)在空間位置xi和xi+h處的實(shí)測值i=1,2,N(h)，那么，變異函數(shù)(h)的離散計(jì)算公式為：三、變異函數(shù) 這樣，對不同的空間分隔距離h，可計(jì)算出相應(yīng)的c(h)和(h)值。如果分別以h為橫坐標(biāo)，c(h)或(h) 為縱坐標(biāo)，畫出協(xié)方差函數(shù)和變異函數(shù)曲線圖，就可直接展示區(qū)域化變量Z(x)的空間變異特點(diǎn)?？梢姡儺惡瘮?shù)能同時(shí)描述區(qū)域化變量的隨機(jī)性和結(jié)構(gòu)性，從而在數(shù)學(xué)上對區(qū)域化變量進(jìn)行嚴(yán)格分析，是空間變異規(guī)律分析和空間結(jié)構(gòu)分析的有效工具。 例如，假設(shè)某地區(qū)降水量Z(x)（單位：mm）是二維區(qū)域化隨機(jī)變量，滿足二階平穩(wěn)假設(shè)，其觀測值的空間正方形網(wǎng)格數(shù)據(jù)如圖所示（點(diǎn)與點(diǎn)之間的距離為h=1

10、km）。試計(jì)算其南北方向及西北-東南方向的變異函數(shù)。圖9.1空間正方形網(wǎng)格降水量數(shù)據(jù)如果沒有缺失值，可直接對正方形網(wǎng)格數(shù)據(jù)結(jié)構(gòu)計(jì)算變異函數(shù)；在有缺失值的情況下，也可以計(jì)算變異函數(shù)。從圖可以看出，空間上有些點(diǎn)，由于某種原因沒有采集到，可采用“跳過”缺失點(diǎn)位置的方法進(jìn)行計(jì)算。圖9.2缺失值情況下樣本數(shù)對的組成和計(jì)算過程 為缺失值 首先計(jì)算南北方向上的變異函數(shù)值，由變異函數(shù)的計(jì)算公式可得： =385/72=5.35同樣計(jì)算出最后，得到南北方向和西北東南上的變異函數(shù)計(jì)算結(jié)果見下表。同樣可以計(jì)算東西方向上的變異函數(shù)。 南北方向西北東南方向h12345h1.412.824.245.657.07N(h) 3

11、6 27 21 13 5 N(h) 32211382(h)5.359.2617.5525.6922.90(h)7.0612.9530.8558.1350.004.變異函數(shù)的參數(shù) 三、變異函數(shù) 這些參數(shù)決定了變異函數(shù)的形狀，反映了自然現(xiàn)象空間分布的結(jié)構(gòu)或空間相關(guān)的類型，同時(shí)還能給出這種空間相關(guān)的范圍。變異函數(shù)有四個(gè)非常重要的參數(shù)，即: 基臺值(Sill); 變程(Range); 塊金值(Nugget); 分維數(shù)(Fractal Dimension)。 基臺值意味著在對應(yīng)（或大于）距離的樣點(diǎn)之間沒有空間相關(guān)性，因?yàn)榉讲畈辉匐S距離變化。 (1)基臺值(sill) 當(dāng)變異函數(shù)隨著間隔距離h的增大，從非

12、零值達(dá)到一個(gè)相對穩(wěn)定的常數(shù)時(shí)，該常數(shù)稱為基臺值C0+C，它是系統(tǒng)或系統(tǒng)屬性中最大的變異。a(2)變程(Range) 曲線從較低的方差值升高，到一定的間隔值時(shí)到達(dá)基臺值，這一間隔稱為變程。 a 在變程內(nèi)，樣點(diǎn)越接近，兩點(diǎn)之間相似性、即空間上的相關(guān)性越強(qiáng)。很明顯，如果某點(diǎn)與已知點(diǎn)距離大于變程，那么該點(diǎn)數(shù)據(jù)不能用于數(shù)據(jù)內(nèi)插（或外推），因?yàn)榭臻g上的自相關(guān)性不復(fù)存在。 變程是變異函數(shù)最重要的參數(shù)，它描述了該間隔內(nèi)樣點(diǎn)的空間相關(guān)特征。(3)塊金(nugget ) 理論方差函數(shù)曲線不穿過原點(diǎn)，而是存在一個(gè)最小的方差值。 a 當(dāng)不計(jì)分析誤差時(shí)，一個(gè)樣品自身的品位誤差應(yīng)等于0，當(dāng)兩個(gè)樣品即使相距很近，品位間仍存

13、在差異時(shí)，這種現(xiàn)象稱為“塊金效應(yīng)”，C0反映了區(qū)域化變量內(nèi)部隨機(jī)性的可能程度。 塊金是在間隔距離小于采樣間距時(shí)的測量誤差或空間變異，或者是二者的和。測量誤差是由儀器的內(nèi)在誤差引起的，空間變異是自然現(xiàn)象在一定空間范圍內(nèi)的變化。小于采樣間距的微觀尺度上空間變異是塊金的一部分。 上述三個(gè)參數(shù)可從變異函數(shù)曲線圖直接得到，或通過估計(jì)曲線回歸參數(shù)得到。 塊金與基臺的比值（C0/(C+C0)，基底效應(yīng)）可以用來說明空間的變異特征，該值越大，說明空間變異更多的是隨機(jī)成分引起的，否則，則是由特定的地理過程或多個(gè)過程綜合引起的。空間相關(guān)性的強(qiáng)弱，可用C/(C+C0)表示，該值越高，表明空間相關(guān)性越強(qiáng)。 （4）分維

14、數(shù) 用于表示變異函數(shù)的特性，由變異函數(shù)(h)和間隔距離h之間的關(guān)系確定： 分維數(shù)D的大小，表示變異函數(shù)曲線的曲率，可以作為隨機(jī)變異的量度。5.變異函數(shù)的功能 通過“變程”反映變量的影響范圍；在原點(diǎn)處的性狀可反映變量的空間連續(xù)性；不同方向上變差圖可反映區(qū)域化變量的各向異性；基臺值的大小可反映變量在該方向上變化幅度的大小；塊金常數(shù)C0的大小可反映區(qū)域化變量的隨機(jī)性。三、變異函數(shù) 6.變異函數(shù)的理論模型 地質(zhì)統(tǒng)計(jì)學(xué)將變異函數(shù)理論模型分為三大類：有基臺值模型，包括球狀模型、指數(shù)模型、高斯模型、線性有基臺值模型和純塊金效應(yīng)模型；無基臺值模型，包括冪函數(shù)模型、線性無基臺值模型、拋物線模型；孔穴效應(yīng)模型。三

15、、變異函數(shù) （1）純塊金效應(yīng)模型。其一般公式為：下面有代表性地介紹幾種常見的變異函數(shù)理論模型。式中：c00，為先驗(yàn)方差。該模型相當(dāng)于區(qū)域化變量為隨機(jī)分布，樣本點(diǎn)間的協(xié)方差函數(shù)對于所有距離h均等于0，變量的空間相關(guān)性不存在。 （2）球狀模型(Spherical model) 其一般公式為： 式中：c0為塊金（效應(yīng)）常數(shù)，c為拱高，c0+c為基臺值，a為變程。當(dāng)c0=0，c=1時(shí)，稱為標(biāo)準(zhǔn)球狀模型。 球狀模型是地質(zhì)統(tǒng)計(jì)分析中應(yīng)用最廣泛的理論模型，許多區(qū)域化變量的理論模型都可以用該模型去擬合。式中：c0和c意義與前相同，但a不是變程。當(dāng)h=3a時(shí)， ，即 ，從而指數(shù)模型的變程a約為3a。當(dāng)c0=0，

16、c=1時(shí)，稱為標(biāo)準(zhǔn)指數(shù)模型。（3）指數(shù)模型(Exponential model) 其一般公式為：（4）高斯模型（Gaussian model） 其一般公式為：式中：c0和c意義與前相同，a也不是變程。當(dāng)時(shí)， ，即 ，因此高斯模型的變程a約為 a 。當(dāng)c0=0，c=1時(shí)，稱為標(biāo)準(zhǔn)高斯函數(shù)模型。球狀模型、指數(shù)模型和高斯模型的比較 球狀模型的變程最小，指數(shù)模型的變程最大，高斯模型的變程介于二者之間。球狀模型和指數(shù)模型過原點(diǎn)存在切線，高斯模型則沒有。 （5）冪函數(shù)模型。其一般公式為： 式中：為冪指數(shù)。當(dāng)變化時(shí)，這種模型可以反映在原點(diǎn)附近的各種性狀。（6）對數(shù)模型。其一般公式為： 顯然，當(dāng) ，這與變異函

17、數(shù)的性質(zhì)(h) 0不符。因此，對數(shù)模型不能描述點(diǎn)支撐上的區(qū)域化變量的結(jié)構(gòu)。（7）線性有基臺值模型。其一般公式為： 式中該模型的變程為a，基臺值為 。 （8）線性無基臺值模型。其一般公式為 從式中可以看出，該模型沒有基臺值，也沒有變程。 在空間分析中，根據(jù)數(shù)據(jù)的空間自相關(guān)性和研究對象的先驗(yàn)知識，選擇要使用的模型。 選定了理論模型后，通常是用最小二乘法計(jì)算方程的各個(gè)參數(shù)，并用最大似然法（ML）來選擇擬合效果最好的模型。此外，在擬合中可以疊加使用上面的模型，模型的參數(shù)也可以進(jìn)行調(diào)節(jié)，以達(dá)到最大程度地?cái)M合實(shí)驗(yàn)數(shù)據(jù)。 例：某地區(qū)降水量是一個(gè)區(qū)域化變量，其變異函數(shù) 的實(shí)測值及距離h的關(guān)系見下表，試用回歸分

18、析方法建立其球狀變異函數(shù)模型。從上面的介紹和討論，我們知道，球狀變異函數(shù)的一般形式為：當(dāng) 時(shí)，有：如果記: 計(jì)算可知，上式的顯著性檢驗(yàn)參數(shù): F=114.054，F(xiàn)0.05(2,9)=4.26，R2=0.962可見模型的擬合效果是很好的。(9-1)則可以得到線性模型：根據(jù)表中的數(shù)據(jù)，對上式進(jìn)行最小二乘擬合，得到：(9-2)比較(9-2)式與(9-1)式，并做簡單計(jì)算可知： c0=2.048，c=1.154，a=8.353所以，球狀變異函數(shù)模型為：(9-3)1.克立格法概述 克立格（Kriging）法，又稱空間局部估計(jì)或空間局部插值法，是地質(zhì)統(tǒng)計(jì)學(xué)的主要內(nèi)容之一。四、克立格法簡介 克立格法根據(jù)待

19、估樣本點(diǎn)（或塊段）有限鄰域內(nèi)若干已測定的樣本點(diǎn)數(shù)據(jù)，對待估樣本點(diǎn)值進(jìn)行估計(jì)，其優(yōu)點(diǎn)是考慮了樣本點(diǎn)的形狀、大小和空間相互位置關(guān)系，及與待估樣本點(diǎn)的相互空間位置關(guān)系?？肆⒏穹ㄊ墙⒃谧儺惡瘮?shù)理論及結(jié)構(gòu)分析基礎(chǔ)之上的，它是在有限區(qū)域內(nèi)對區(qū)域化變量的取值進(jìn)行無偏最優(yōu)估計(jì)的一種方法。1.克立格法概述 四、克立格法簡介 （1）適用條件 變異函數(shù)和相關(guān)分析的結(jié)果表明區(qū)域化變量存在空間相關(guān)性。 其實(shí)質(zhì)是利用區(qū)域化變量的原始數(shù)據(jù)和變異函數(shù)的結(jié)構(gòu)特點(diǎn)，對未采樣點(diǎn)的區(qū)域化變量的取值進(jìn)行線性無偏、最優(yōu)估計(jì)。1.克立格法概述 四、克立格法簡介 （2）克立格法的類型 克里格插值（Kriging Interpolation

20、），是根據(jù)變異函數(shù)模型而發(fā)展起來的一系列地質(zhì)統(tǒng)計(jì)的空間插值方法，包括：普通克里格法（Ordinary Kriging）泛克里格法（Universal Kriging）指示克里格法（Indicator Kriging）析取克里格法（Disjunctive Kriging）協(xié)同克里格法（Cokriging）等。2.克立格估計(jì)量 四、克立格法簡介 對于研究區(qū)內(nèi)任一點(diǎn)x的測量值Z(x)，其估計(jì)值 的估算公式為： 估計(jì)量 是實(shí)際值Z(x)的克立格估計(jì)值，其中i 為權(quán)重系數(shù)，表示各空間樣本點(diǎn)xi 處的觀測值Z(xi)對估計(jì)值 的貢獻(xiàn)程度。 估計(jì)量 的好壞取決于怎樣計(jì)算或選擇權(quán)重系數(shù)i ，問題的關(guān)鍵在于求各

21、點(diǎn)的權(quán)重系數(shù)。3.普通克立格法 四、克立格法簡介 設(shè)Z(x)為區(qū)域化變量，滿足二階平穩(wěn)假設(shè)和本征假設(shè)，其數(shù)學(xué)期望為m，協(xié)方差函數(shù)c(h)及變異函數(shù)(h)存在。（1）克立格估計(jì)量假設(shè)在待估計(jì)點(diǎn)x0的鄰域內(nèi)共有n個(gè)實(shí)測點(diǎn)，即x1,x2,xn，其樣本值為Z(xi)。那么，普通克里格法的插值公式為：3.普通克立格法 四、克立格法簡介 （2）確定權(quán)重系數(shù) 權(quán)重系數(shù)的求取必須滿足兩個(gè)條件：無偏性。即偏差的數(shù)學(xué)期望為零；最優(yōu)性。即使估計(jì)值 和實(shí)際值Z(x)之差的平方和最小。使用協(xié)方差函數(shù)表達(dá)，它可以進(jìn)一步寫為： (9-4) 為使估計(jì)方差最小，根據(jù)拉格朗日乘數(shù)原理，令： (9-5) 求F對i和的偏導(dǎo)數(shù)，并令其為0，得克立格方程組： (9-6) 整理后得： (9-7) 解此線性方程組，求出權(quán)重系數(shù)i和拉格朗日乘數(shù)，代入公式(9