指示克立格法在區(qū)域淺層地下水硝酸鹽含量分析中的應(yīng)用

指示克立格法在區(qū)域淺層地下水硝酸鹽含量分析中的應(yīng)用

通常，在分析土壤資源特征時(shí)會(huì)出現(xiàn)一些特定值，這些特定值不是由分析錯(cuò)誤或采樣方法引起的，而是由研究對(duì)象的客觀部分引起的。它只占所有數(shù)據(jù)的最后一部分，但其作用不容忽視。因此，地質(zhì)學(xué)家和統(tǒng)計(jì)學(xué)家進(jìn)行了大量有效的研究，并生成了非參數(shù)統(tǒng)計(jì)方法。也被稱(chēng)為無(wú)分布法，并描述了一種不確定的統(tǒng)計(jì)方法。在不排除重要和真實(shí)數(shù)據(jù)的情況下，處理各種現(xiàn)象，提供非知識(shí)的估計(jì)值和概率，并根據(jù)一定的風(fēng)險(xiǎn)條件給出不知道的估計(jì)值和概率。侯景如用這種方法估算了礦漿石的儲(chǔ)量。然而，關(guān)于該區(qū)域地下水環(huán)境評(píng)價(jià)的文獻(xiàn)很少。在這項(xiàng)工作中，我們使用克立格法研究了區(qū)域平坦地下水中氮含量的空間分布，并提供了超過(guò)國(guó)家規(guī)定的最大允許值的概率分布圖。1基本原則1.1閾值的二元化設(shè)在區(qū)域D上取樣并測(cè)定了某觀測(cè)項(xiàng)目,若區(qū)域D上的閾值或邊界值為Z,則在D上的每個(gè)樣品點(diǎn)x∈D上定義1個(gè)Z的如下階梯函數(shù):i(x;Ζ)={1當(dāng)x點(diǎn)上的觀測(cè)值Ζ(x)≤Ζ0當(dāng)x點(diǎn)上的觀測(cè)值Ζ(x)>Ζ(1)i(x;Z)={10當(dāng)x點(diǎn)上的觀測(cè)值Z(x)≤Z當(dāng)x點(diǎn)上的觀測(cè)值Z(x)>Z(1)在D上的任一區(qū)域A∈D內(nèi),低于閾值Z的觀測(cè)值Z(x)所占區(qū)域A的比例表示為:?(A;Ζ)=1A∫Ai(x;Ζ)dx∈(2)?(A;Z)=1A∫Ai(x;Z)dx∈(2)?(A;Z)是關(guān)于Z(x)和Z的二元函數(shù),即小于閾值Z的全部i(x;Z)(x∈A)的平均值.1.2rob/nb運(yùn)行參數(shù)選取在給定的閾值Z的條件下,隨機(jī)函數(shù)I(x;Z)服從二項(xiàng)分布,其期望值是:E{Ι(x;Ζ)}=1?Ρrob{Ζ(x)≤Ζ}+0?Ρrob{Ζ(x)>Ζ}=Ρrob{Ζ(x)≤Ζ}=F(Ζ)(3)E{I(x;Z)}=1?Prob{Z(x)≤Z}+0?Prob{Z(x)>Z}=Prob{Z(x)≤Z}=F(Z)(3)式中的F(Z)為觀測(cè)值Z(x)的分布函數(shù)在Z(閾值)處的值.當(dāng)I(x+h;Z)和I(x;Z)為被矢量h分隔的兩個(gè)隨機(jī)變量時(shí),則指示半變異函數(shù)可定義如下:γ1(h;Ζ)=12E{[Ι(x+h;Ζ)-Ι(x;Ζ)]2}(4)γ1(h;Z)=12E{[I(x+h;Z)?I(x;Z)]2}(4)1.3指示克立格方程和空間分布的指示克立格估算1.3.1指示函數(shù)空間的建立設(shè)在區(qū)域D上有N個(gè)有效數(shù)據(jù),在D上的一個(gè)域A(A∈D)內(nèi)有n個(gè)有效數(shù)據(jù),{Z(xα),xα∈A,α=1,2,…n},在給定閾值之后,得到樣品的指示函數(shù)空間;{i(xα;Z),α=1,2,…,n},則?(A;Z)的估計(jì)值可表示如下:?*(A;Ζ)=n∑α=1λα(Ζ)?i(xα;Ζ)(5)??(A;Z)=∑α=1nλα(Z)?i(xα;Z)(5)1.3.2拉格朗日乘子法估計(jì)方差如同普通克立格法一樣,為了求得?*(A;Z)必須在無(wú)偏和估計(jì)方差極小的條件下求式(5)中的權(quán)系數(shù)λα(Z)(α=1,2,…,n).根據(jù)式(2)、式(3),偏差的期望可寫(xiě)成:E{?*(A;Ζ)-?(A;Ζ)}=n∑α=1λα?E{Ι(xα;Ζ)}-F(Ζ)=[n∑α=1λα-1]?F(Ζ)=0(6)E{??(A;Z)??(A;Z)}=∑α=1nλα?E{I(xα;Z)}?F(Z)=[∑α=1nλα?1]?F(Z)=0(6)當(dāng)n∑α=1λα=1∑α=1nλα=1時(shí),滿(mǎn)足上式,所以,無(wú)偏條件是:n∑α=1λα=1.與普通克立格法相似,也可以給出指示克立格法的估計(jì)方差σ2EΙ:σ2EΙ=E{[?*(A;Ζ)-?(A;Ζ)]2}=2n∑α=1ˉγi(xα,A；Ζ)-ˉγi(A,A;Ζ)-2n∑α=1n∑β=1λαλβˉγi(xα,xβ;Ζ)(7)要使估計(jì)方差在無(wú)偏條件下變?yōu)樽钚?這是一個(gè)求條件極值的問(wèn)題,同樣,采用標(biāo)準(zhǔn)拉格朗日乘子法,最后得到對(duì)于n+1個(gè)未知數(shù)(n個(gè)λ和1個(gè)μ)的n+1個(gè)方程的指示克立格方程組:{n∑β=1λβˉγi(xα,xβ；Ζ)+μ=ˉγi(xα,A;Ζ)n∑α=1λα=1(α=1,2,?,n)(8)式(7)和式(8)中:μ為拉格朗日乘子;ˉγi(xα,xβ;Z)為在給定的閾值Z條件下,矢量h的兩個(gè)端點(diǎn)分別在信息域xα、xβ內(nèi)的所有對(duì)點(diǎn)的平均指示半變異函數(shù)值;ˉγi(A,A;Z)為在給定閾值Z條件下,矢量h的兩個(gè)端點(diǎn)在待估域A內(nèi)所有對(duì)點(diǎn)的平均指示半變異函數(shù)值;ˉγi(xα,A;Z)為在給定閾值Z條件下,矢量h的一個(gè)端點(diǎn)在信息域xα內(nèi),另一端點(diǎn)在待估域A內(nèi)所有對(duì)點(diǎn)的平均指示半變異函數(shù)值.1.3.3點(diǎn)x的估計(jì)值ai(xα;Ζ)=Ρrob{Ζ(xα)≤Ζ|Ζ(xα)=Ζα}(9)這時(shí),某點(diǎn)x的估計(jì)值i*(x;Z)可表示為:i*(x;Ζ)=Ρrob{Ζ(x)≤Ζ|周?chē)鷶?shù)據(jù)}(10)而某待估域A的?(A;Z)的估值?*(A;Z)可表示如下:?*(A;Ζ)=1A∫AΡrob{Ζ(x)≤Ζ}|周?chē)鷶?shù)據(jù)}dx(11)1.4普通克立格法A平均值的指示克立格估計(jì)指示克立格法是應(yīng)用某種克立格法(如普通克立格法)求得?(A;Z)的線(xiàn)性估算值?*(A;Z),最后得到待估域A的平均估計(jì)值[Z(x)]*.[Ζ(x)]*=L∑l=1[Ιl(x)]*?[Ζ(x)|x∈Ζl]*(12)1.5半變異函數(shù)的計(jì)算(1)對(duì)原始數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析:包括樣品分布特征的研究,繪制樣品累積分布曲線(xiàn)等;(2)將原始數(shù)據(jù)變換成指示值:將原始數(shù)據(jù)變換成順序統(tǒng)計(jì)量,做累積直方圖,以10%,20%,…,90%9個(gè)分位數(shù)(根據(jù)實(shí)際情況而定)對(duì)應(yīng)的樣品觀測(cè)值{Z(x10),Z(x20),…,Z(x90)}為確定的閾值Zl=(l=1,2,…,9),求出每一樣品對(duì)應(yīng)于Zl(l=1,2,…,9)的指示值并列成表格.這時(shí),指示值的定義為:i(x;Ζl)={1?當(dāng)Ζ(xα)≤Ζl(l=1,2,?,9)0；當(dāng)Ζ(xα)>Ζl(α=1,2,?,n)(3)計(jì)算每個(gè)閾值下的指示半變異函數(shù)γ*l(h;Zl)(l=1,2,…,9).為了減少工作量,在條件允許時(shí),可用中位數(shù)(50%分位數(shù))指示半變異函數(shù)進(jìn)行研究,或只選那些有明顯差異的指示半變異函數(shù);(4)進(jìn)行結(jié)構(gòu)分析,給出不同閾值下的指示半變異函數(shù)模型;(5)求待估域的平均估計(jì)值和條件概率.2實(shí)驗(yàn)計(jì)劃2.1地下水樣的采集研究區(qū)設(shè)在河北省曲周縣境內(nèi),面積為667km2,位于漳河沖積扇、漳河滏陽(yáng)河沖積平原和黃河沖積平原的交匯處.經(jīng)度36°35′43″～36°57′56″,緯度114°50′22.3″～115°13′27.4″,屬暖溫帶、半濕潤(rùn)、大陸性季風(fēng)氣候區(qū).年平均氣溫為13.1℃,平均降水量為556.2mm,降水主要集中在7～9月份,占全年降水量的2/3,年平均蒸發(fā)量是年平均降水量的3.3倍,地下水埋深一般在4m以下.作物種植制度主要為:小麥-玉米,小麥-大豆,小麥-棉花.小麥生育期內(nèi)主要依賴(lài)地下水灌溉,林地、荒地均很少,有部分菜地.在該縣境內(nèi)農(nóng)田灌溉的淺層水井里隨機(jī)、但盡可能均一地提取地下水樣,共進(jìn)行139個(gè)點(diǎn)的采樣,其地理位置(經(jīng)度、緯度)用GPS定位,采樣時(shí)間為1999年5月.2.2水位含鹽量、地下水硝酸鹽含量(1)地下水位埋深:采取水樣時(shí)用井繩標(biāo)記地下水水面與地面的距離,用刻度尺量得水位.(2)地下水含鹽量:使用DDS-11AT電導(dǎo)率儀測(cè)定.(3)地下水硝酸鹽含量:取樣后加1ml500ppmPMA抑制微生物活性,在冰箱中保存,用自動(dòng)流動(dòng)分析儀測(cè)定.3結(jié)果分析3.1不服從正態(tài)分布從表1、圖1及圖2可以看出:硝酸鹽含量既不服從正態(tài)分布也不服從對(duì)數(shù)正態(tài)分布,其變異系數(shù)為2.25,屬于高變異強(qiáng)度,并且含有特異值.因此,適合用指示克立格法進(jìn)行分析.3.2閾值的選取和說(shuō)明由于地下水中硝酸鹽含量超過(guò)了一定的濃度就會(huì)威脅人體健康,世界衛(wèi)生組織(WHO)頒布的飲用水質(zhì)標(biāo)準(zhǔn)規(guī)定硝酸鹽含量的最大允許濃度為10mg/l.我國(guó)規(guī)定為20mg/l,因而取20mg/l作為閾值.根據(jù)圖1試驗(yàn)累積頻率分布圖,可知閾值20mg/l所對(duì)應(yīng)的百分?jǐn)?shù)為0.97,中位數(shù)0.50所對(duì)應(yīng)的濃度值為0.3mg/l.故選取0.3mg/l和20mg/l作為閾值.求出每個(gè)樣品值所對(duì)應(yīng)的閾值的指示值并列成表格,以便計(jì)算各個(gè)閾值的半變異函數(shù).3.3層地下水中硝酸鹽含量的半變異函數(shù)通過(guò)計(jì)算硝酸鹽原始數(shù)據(jù)不同方向(0°30°60°90°4個(gè)方向)的半變異函數(shù),得到了如圖3的結(jié)果.從圖3可以看出,曲周縣淺層地下水中硝酸鹽含量的分布在不同的方向上的差異不明顯(不存在帶狀和幾何各向異性),故可以認(rèn)為其分布是各向同性的.根據(jù)式(4)可以得到不同閾值下的指示半變異函數(shù)值,用不同的模型來(lái)進(jìn)行擬合,最后得到如圖5的半變異函數(shù)模型.圖5(a)的半變異函數(shù)模型為純塊金形式,表示在所取的采樣密度下,各個(gè)點(diǎn)之間是相互獨(dú)立的.圖5(b)表示用球狀模型擬合了閾值為20mg/l的半變異函數(shù),并給出了理論模型的3個(gè)參數(shù).3.4曲周縣淺層地下水中硝酸鹽含量分布根據(jù)式(12)可以得到不同待估域的指示克立格的平均估計(jì)值(見(jiàn)圖5).根據(jù)式(11)可以得到超過(guò)最大允許濃度20mg/l的條件概率(見(jiàn)圖6).從圖5可以看出,地下水硝酸鹽含量高的地方主要集中在曲周縣的中西部和北部的一小塊區(qū)域,這主要是因?yàn)橹形鞑刻幱诳h城周?chē)粠?這里容納了大量富含氮素的生活廢水,種植蔬菜的面積居全縣之首,并且此區(qū)域也是流經(jīng)曲周縣境內(nèi)幾大河流(滏陽(yáng)河,漳河)匯集之處,而這些河流的水質(zhì)已遭到明顯的污染,所有這些因素都會(huì)造成地下水中硝酸鹽含量升高.北部的一小塊區(qū)域?yàn)榍芸h的第二大鎮(zhèn)——河南疃鎮(zhèn),其生活污水,蔬菜種植面積居全縣次位,故也造成了其硝酸鹽含量較高.從圖6可以看出,其濃度超過(guò)最大允許值20mg/l的概率達(dá)到95%的地方,也主要是這兩個(gè)區(qū)域.因此,就目前而言,應(yīng)在這兩個(gè)區(qū)域采取有效管理措施(如節(jié)水灌溉、科學(xué)施用氮肥等),對(duì)淺層地下水質(zhì)量進(jìn)行保護(hù).將圖5和圖6與用傳統(tǒng)的普通克立格法得到的區(qū)域硝酸鹽含量分布圖對(duì)比分析,可知,在普通克立格法得到的曲周縣淺層地下水硝酸鹽含量分布圖中,河南疃鎮(zhèn)受污染的程度沒(méi)有反映出來(lái),這說(shuō)明了普通克立格法具有明顯的平滑效應(yīng),而采用指示克立格法更能準(zhǔn)確地反映實(shí)際情況,且得到的信息更為豐富.4非參數(shù)地統(tǒng)計(jì)學(xué)方法在水土資源特性的數(shù)據(jù)分析中經(jīng)常會(huì)出現(xiàn)特異值,這些值不可任意地舍棄,因?yàn)樗鼈兯鸬淖饔猛遣豢珊鲆暤?普通克