指示克立格估值計算的理論和方法

指示克立格估值計算的理論和方法

近十年來，隨著人口的增加和經(jīng)濟的發(fā)展，許多地區(qū)的地下水開發(fā)利用程度日益嚴重，導致了一系列區(qū)域地質(zhì)環(huán)境問題。地下水數(shù)值模擬作為一種常規(guī)研究工具,規(guī)模和范圍也日趨擴大。此類區(qū)域性模型應用中的一個突出問題是勘探資料特別是含水層的滲透性參數(shù)資料分布不均或缺乏,在通常的地下水數(shù)值模擬中只能用簡單的參數(shù)分區(qū)來描述滲透性的非均質(zhì)性。由于分區(qū)過大而忽略分區(qū)內(nèi)部參數(shù)的差異性,導致模擬誤差的出現(xiàn)。尋求一種盡可能利用有限的勘探資料,對未知區(qū)域內(nèi)的含水層參數(shù)進行合理估值的方法,是目前大區(qū)域地下水流模擬中的關(guān)鍵性問題。地質(zhì)統(tǒng)計學中的克立格方法是線性無偏最優(yōu)的估值方法。20世紀70年代末之后,克立格方法在水文地質(zhì)學領(lǐng)域得到了廣泛應用。Delhomme(1979)用克立格法估計法國下諾曼底含水層的導水系數(shù)并繪制了導水系數(shù)的最優(yōu)等值線圖;Dagan等(1982)用相同方法繪制了意大利威尼斯市地下水頭分布的等水頭線并以此作為研究地面沉降的依據(jù);Neuman和Cliftun(1982)用克立格法給出了美國亞利桑拉Avra河谷含水層滲透系數(shù)的估計;Ahmed等(1987)也結(jié)合實例對使用克立格法求滲透系數(shù)與其他幾種方法進行了比較。由于自然界大量的數(shù)據(jù)樣本不服從正態(tài)分布,用高斯克立格法進行估值比較困難,為此Journel于1983年創(chuàng)立了指示克立格法,它不要求數(shù)據(jù)總體服從正態(tài)分布,且不受特異值的影響。指示克立格法自創(chuàng)立以來,在礦產(chǎn)(Fytas和Chaouai,1990)、土壤(Goovaerts,Webster和Dubois,1997)、土地管理(傅佩紅,2004)、石油(Data-Gupta等,1999和Hohn等,1994)等領(lǐng)域有廣泛應用。然而,使用指示克立格法對第四紀含水層滲透性空間分布進行估值國內(nèi)外還鮮有報道。本文以華北某第四系地層為例,以指示變異函數(shù)為結(jié)構(gòu)分析工具,運用指示克立格法對該區(qū)域含水層滲透性的空間分布進行三維估計,結(jié)果表明指示克立格法能很好地描述第四系含水層滲透性的空間分布規(guī)律。1拉格朗日函數(shù)中的克立格估計在研究含水層滲透性時,把含水層的滲透系數(shù)作為區(qū)域化變量,記作Z(u),它同時具有隨機性和結(jié)構(gòu)性的特點。建立變量Z(u)的模型時,設滲透系數(shù)在區(qū)域范圍內(nèi)的變化為二階平穩(wěn)過程,即其均值函數(shù)和協(xié)方差函數(shù)存在且平穩(wěn),則所有取樣點所形成的分布可以代替總體分布,而區(qū)域內(nèi)任一點處變量取值的概率分布服從這一總體。作所有取樣點的累積頻率曲線,將累積頻率曲線由小到大截成K+1段,則有K個截斷值(閾值)Zk,k=1,2,…,K,通常以其十分位數(shù)作為閾值(如果是類型變量,則以其類型作為閾值)。為求出位置u處的滲透系數(shù)Zu小于等于閾值Zk的概率,引進二值指示函數(shù),其定義為Ι(u?Ζk)={1當位置u處的滲透系數(shù)Ζu小于等于閾值Ζk時0其他(1)I(u?Zk)=?????1當位置u處的滲透系數(shù)Zu小于等于閾值Zk時0其他(1)指示變異函數(shù)可描述為γΙ(h)=1Ν(h)Ν(h)∑α=1[i(uα+h,Ζk)-i(uα,Ζk)]2(2)γI(h)=1N(h)∑α=1N(h)[i(uα+h,Zk)?i(uα,Zk)]2(2)它與協(xié)方差函數(shù)有如下關(guān)系:γI(h)=CI(0)-CI(h),CI(0)=Var[I(u,z)]為方差。在u處出現(xiàn)滲透系數(shù)Zu小于等于閾值Zk的數(shù)學期望是:E{Ι(u,Ζk)}=1?Ρrob{Ζ(u)≤Ζk}+0?Ρrob{Ζ(u)>Ζk}=Ρrob{Ζ(u)≤Ζk}=F(Ζk)(3)E{I(u,Zk)}=1?Prob{Z(u)≤Zk}+0?Prob{Z(u)>Zk}=Prob{Z(u)≤Zk}=F(Zk)(3)上式表明在u處出現(xiàn)滲透系數(shù)Zu小于等于閾值Zk的數(shù)學期望等于它出現(xiàn)的概率。所以,指示克立格法估計量能夠用于估計特征出現(xiàn)的概率,這個估計值可通過實測滲透系數(shù)的線性函數(shù)表示:F(u,Ζk)*=E{Ι(u,Ζk)ㄧ(n)}*=Ρrob*{Ζ(u)≤Ζkㄧ(n)}=n∑α=1λα(u,Ζk)Ι(uα,Ζk)(4)F(u,Zk)?=E{I(u,Zk)ㄧ(n)}?=Prob?{Z(u)≤Zkㄧ(n)}=∑α=1nλα(u,Zk)I(uα,Zk)(4)式中F(u,Zk)*為滲透系數(shù)在位置u處真實概率F(u,Zk)的克立格估計;uα是第α個實測滲透系數(shù)的位置,α=1,2,…,n;λα(u,Zk)是相應的截斷值Zk的簡單克立格法權(quán)重。要確定每個權(quán),必須使估計方差達到最小,即:min{E[F(u,Ζk)*-F(u,Ζk)]2}min{E[F(u,Zk)??F(u,Zk)]2}且受無偏性條件n∑α=1λα(u,Ζk)=1的限制。構(gòu)造拉格朗日函數(shù)L(λ1,…,λn,μ),μ為拉格朗日乘數(shù)。用L對每一個未知的權(quán)系數(shù)和拉格朗日乘數(shù)求偏導并使之為0,得正規(guī)方程組:{n∑β=1λβ(u,Ζk)CΙ(uβ-uα,Ζk)+μ=CΙ(u-uα,Ζk),α=1?2???nn∑α=1λα(u,Ζk)=1(5)解此方程組可求得諸λα(u,Zk),代入(4)式即可得到概率F(u,Zk)的克立格估值,它就是待估點處閾值為Zk的條件累積分布概率。這里的CI(uα-uβ,Z)為指示協(xié)方差函數(shù)。因滲透系數(shù)為二階平穩(wěn),指示協(xié)方差函數(shù)和變異函數(shù)與其位置無關(guān),只與位置的增量h有關(guān),故有:CΙ(uα-uβ,Ζ)=Cov{Ι(uα,Ζ),Ι(uβ,Ζ)}=CΙ(h),|α-β|=h(6)相應的最小估計方差為σ2(u)=C(0)-n∑β=1λβ(u,Ζk)CΙ(u-uβ,Ζk)(7)指示克立格法過程就是重復一系列K個截斷值Zk,k=1,2,…,K的過程,對每一個指示函數(shù)都同樣地使用上述的指示克立格算法。通過(4)式可以建立條件累積分布函數(shù)代表未采樣值Z(u)的不確定性概率模型。估計累計分布函數(shù)(cdf)實質(zhì)上是在每個估計位置或單元上綜合所有的指示估計。圖1是用多重指示克立格法估計出的條件cdf的示意圖。應指出的是條件cdf的估計需要K個指示函數(shù)的K個變差函數(shù)模型。通過對區(qū)域化變量Z的多個截斷值Zk處[i(u,Z)]*分別估值,得到待估點處的各滲透系數(shù)截斷值所對應的條件累積分布概率[Prob{Z(u)≤Zk|(n)}]*。然后采用E型估計可最終得到待估點處的滲透系數(shù)估計值Z(u),即:[Ζ(u)]*=∫+∞-∞ΖdF(u,Ζ|(n))≈Κ+1∑k=1Ζk[F(u,Ζk)*-F(u,Ζk-1)](8)對每個待估點或單元分別采用以上步驟,就可得到所有單元上的滲透系數(shù)克立格估值。2研究區(qū)地質(zhì)背景及含水層滲透率華北平原位于燕山山脈以南,淮河以北,太行山及豫西山地以東,東臨渤海,由灤河、海河、黃河沖積而成,是一個地勢坦蕩、第四系深厚的全國最遼闊的沖洪積平原。華北平原自山前向濱海分為沖洪積扇區(qū)、洪泛平原區(qū)、濱海平原區(qū),地層結(jié)構(gòu)分區(qū)明顯受物源區(qū)控制,平面上展示近山前是滾動組分占優(yōu)勢的沉積物,稍離山前地帶即以跳躍組分為主的沉積物,再遠離山前地帶是跳躍、懸移物質(zhì)組成的沉積物,其間以泥質(zhì)為絕對優(yōu)勢的沉積物在洼地沉積下來。沉積物的粒度級配變化隨離物源區(qū)漸遠,水動力減弱而變細。表現(xiàn)在滲透性上由山前向濱海逐漸變低。垂向上,第四系地層結(jié)構(gòu)與大型河流沖積扇發(fā)育歷史、構(gòu)造、氣候背景等息息相關(guān)。平原沉積物主要來自河流搬運,出山河流可以有平行、匯聚、分散、曲線、雜亂、旋轉(zhuǎn)、倒流等多種形式,造成河流擺動非常復雜,因而擺動形成的砂體在時空展布上的規(guī)律也非常復雜。本次選擇華北平原中部鉆孔控制程度較高的區(qū)域為研究區(qū),其西面緊靠山前、東面臨近濱海,范圍為250km×200km,該范圍內(nèi)共有各類鉆孔192個。為研究含水層的滲透性,將含水層的滲透率作為區(qū)域化變量?？紤]到鉆孔中水文地質(zhì)鉆孔較少,含水層滲透率的直接數(shù)據(jù)不多,這里以鉆孔巖性為依據(jù),參考當?shù)貛r性滲透性的經(jīng)驗值,并結(jié)合水文地質(zhì)鉆孔的數(shù)據(jù)給出各巖性的滲透系數(shù)。因克立格估值要求估計區(qū)域規(guī)則,且最好是區(qū)域化變量的最大變化方向與坐標軸方向一致,因此將坐標系統(tǒng)作一變換:先旋轉(zhuǎn)坐標系使X軸與區(qū)域化變量的最大變化方向一致,然后平移坐標軸使區(qū)域的左下角與原點重合。為簡化計算,X、Y軸的坐標值分別縮小1000倍,以km為單位。Z軸方向上,為消除地形起伏的影響,將高程轉(zhuǎn)化為埋深。此次研究埋深為200m。鉆孔揭露的不同巖層厚薄差異懸殊給變異函數(shù)的搜索造成困難,因此在鉆孔方向上以10m為單位統(tǒng)一承載,以巖層厚度為權(quán)進行調(diào)和平均計算各承載上的滲透系數(shù)。3數(shù)據(jù)總體特征為右偏態(tài)分布,基于微織構(gòu)在進行坐標轉(zhuǎn)換和滲透系數(shù)數(shù)據(jù)整理后,首先進行傳統(tǒng)統(tǒng)計學分析,作出滲透系數(shù)的直方圖,考察其是否服從正態(tài)分布或?qū)?shù)正態(tài)分布,并求出各分位數(shù)。如果服從正態(tài)分布或?qū)?shù)正態(tài)分布,則可以進行高斯克立格估值,否則可用指示克立格進行估值,因為指示克立格法是非參數(shù)估值方法,它不要求研究的數(shù)據(jù)總體服從何種分布。其直方圖如圖2所示。本次分析共有3592個滲透系數(shù)數(shù)據(jù),最大值為150m/d,最小值為0.001m/d,均值為1.969m/d,中位數(shù)為0.031m/d,數(shù)據(jù)總體呈右偏態(tài)分布,極端大值將算術(shù)平均數(shù)拉向極端大值一方,極端大值對均值貢獻大。離散系數(shù)為5.701,說明變異值的離散程度極大,平均數(shù)代表性極差。數(shù)據(jù)總體不服從正態(tài)分布或?qū)?shù)正態(tài)分布。4實驗變異函數(shù)擬合取中位數(shù)0.03,和1.0、10.0作為閾值,分X、Y、Z三個方向計算各閾值的指示變異函數(shù)。滯后距X、Y方向為10,Z方向為20。由于鉆孔位置分布不規(guī)則,在搜索數(shù)據(jù)對時以滯后距的一半作為距離容限,以22.5°作為方向容限,計算實驗變異函數(shù)。將所得到的實驗變異函數(shù)以球狀模型進行擬合,即:γ(h)=C0+γΙ(h)=C0+CΙ(3h2a-h32a3)(9)其中:C0為塊金常數(shù),CI為拱高,C0+CI為基臺值,a為變程。限于篇幅只給出域值為1.0實驗變異函數(shù)擬合圖(圖3)(實線為實驗變異函數(shù),虛線為擬合的理論變異函數(shù))。三個閾值指示變異函數(shù)擬合的參數(shù)見表1。參數(shù)a(變程)反映區(qū)域化變量的相關(guān)程度,變程越大,滲透系數(shù)的相關(guān)性越好,越小則相關(guān)性越差。塊金常數(shù)(C0)反映可測范圍內(nèi)原點處的性狀,塊金常數(shù)越大,滲透系數(shù)即使在很短矩離內(nèi)其差異性也很大。由圖表可以看出,不同閾值不同方向其參數(shù)各不相同。X軸方向變程約30km,Y軸方向約為20km,Z軸方向為30m,即X軸變程大于Y軸變程,當然遠大于Z軸方向,反映出X軸方向滲透性的連續(xù)性較Y軸方向好,且遠遠大于Z軸方向。各個閾值各個方向的塊金常數(shù)都較大,表明滲透性在相鄰很短的距離內(nèi)其變異性就很大。5滲透系數(shù)估計值的估計克立格估值為線性最優(yōu)無偏估計,即估計值與真值的均值相同、方差最小,指示克立格法同樣也遵循這一原則。將研究區(qū)250km×200km×200km的三維區(qū)域劃分為50×40×20的網(wǎng)格,每個網(wǎng)格大小為5km×5km×10km,共有40000個矩形格子。對每個格子進行估值。根據(jù)各閾值指示變異函數(shù)擬合參數(shù),由公式(9)得到各距離h處的γ(h)和CI(h),再代入公式(3)可得到[i(u,Z)]*的普通克立格估值,此為各閾值在每個格子上的條件累積分布概率[Prob{Z(u)≤Zk|(n)}]*。最后利用公式(7)進行E型估計即可得到各個格子的滲透系數(shù)估計值Z(u)。圖4是利用斯坦福大學Gslib軟件制作的研究區(qū)滲透性分布的水平和垂向兩個方向的切面,它們分別是Z軸方向的第20個切面、Y軸方向的第20個切面。由于Z軸方向坐標單位比X軸、Y軸方向小3個數(shù)量級,垂直方向的切面圖所反映的結(jié)構(gòu)特征不甚明顯。從切面圖可以看出,水平方向上,近山前地區(qū)地層滲透性最大,約為80m/d,向濱海區(qū)逐漸過渡變小,約為0.1m/d,在東部有一滲透性異常小的區(qū)域約為0.01m/d,是為以泥質(zhì)占絕對優(yōu)勢的沉積物在洼地沉積;在垂直方向上,滲透性變化不明顯,淺部比深部略好,這與該區(qū)域的地層結(jié)構(gòu)水平分布及其物性組成特征是一致的?？肆⒏穹椒ú粌H能夠給出待估點處的估計值,同時也能夠給出估計方差,它是度量估計精度的一個很好的尺度。圖5為滲透性估計方差等值線圖,對估計方差進行標準化,使其在0~1范圍之間,圖中圓點為鉆孔的位置。由圖中可以看出,有鉆孔控制的地方,其估計方差就小,估計精度高;鉆孔密度小的地方,則估計方差相對就大。局部地方如圖的西部估計方差較大,達到0.8以上,說明該區(qū)域地層結(jié)構(gòu)復雜,滲透性變異性大,工程控制程度低。為使估計精度滿足要求,應該在該區(qū)域增加工程控制。區(qū)域化變量的空間變異結(jié)構(gòu)分析是否正確、結(jié)構(gòu)模型是否適用于本區(qū)等需要經(jīng)過交叉驗證,只有經(jīng)過驗證證明其對本區(qū)滲透性的指示克立格法估計確實可靠且適用時才能采納。以滲透性的空間結(jié)構(gòu)模型參數(shù)為基礎,對該區(qū)已知樣點進行指示克立格法交叉驗證,作差值頻數(shù)分布直方圖(圖6)。由圖6可以看出,差值基本服從均值為0的正態(tài)分布,變化范圍在(-14.56,8.67)內(nèi),均值為0.0014;差值在(-1,1)區(qū)間內(nèi)的約占總數(shù)的67.2%,標準差為1.89。上述統(tǒng)計結(jié)果表明,用該空間變異結(jié)構(gòu)模型來對研究區(qū)滲透性進行指示克立格法估計,其結(jié)果是可信的。6含水層滲