指示克立格估值計(jì)算的理論和方法

指示克立格估值計(jì)算的理論和方法

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