內(nèi)蒙古河套灌區(qū)土壤水鹽空間變異性的ik法閾值選擇

內(nèi)蒙古河套灌區(qū)土壤水鹽空間變異性的ik法閾值選擇

指導(dǎo)思想的不規(guī)則性（縮寫ik）法是條庫(kù)爾提出的一種非參數(shù)估計(jì)方法。史密斯和霍夫最初用于土壤和生態(tài)科學(xué)。IK法與其他插值估計(jì)方法比較有明顯優(yōu)點(diǎn):它不依賴空間現(xiàn)象的平穩(wěn)性;也不要求區(qū)域化變量服從某種分布假設(shè),由于對(duì)隨機(jī)變量經(jīng)過特有的指示變換后,使隨機(jī)函數(shù)的貢獻(xiàn)都一樣,可自然消除特異值(Outlier)的影響;可對(duì)非取樣點(diǎn)處的不確定性進(jìn)行條件概率分布函數(shù)(ConditionalCumulativeDistributionFunction簡(jiǎn)寫CCDF)估計(jì),尤其是當(dāng)人們感興趣的不是某些點(diǎn)的具體估值,而是大于或小于某閾值(CutoffValue)或某范圍內(nèi)物質(zhì)品位在空間的分布概率時(shí),可簡(jiǎn)捷獲得在一定風(fēng)險(xiǎn)條件下未知量分布概率的估值,符合地學(xué)科學(xué)家用人工智能思維處理空間問題的思路。IK法自創(chuàng)立以來(lái),在國(guó)外廣泛用于地質(zhì)、土壤、水文(Neuman)等領(lǐng)域,成為一種直接估計(jì)局部區(qū)域品位分布的數(shù)學(xué)地質(zhì)新方法。20世紀(jì)90年代開始在我國(guó)地學(xué)及土壤學(xué)科應(yīng)用,也有一些報(bào)導(dǎo):如地質(zhì)礦產(chǎn)方面侯景儒、土壤水文方面李保國(guó)等,水庫(kù)工程地質(zhì)方面吳蓉等,農(nóng)業(yè)水土工程方面劉全明等、徐英等均有一定成效,但多為局部估計(jì)的初步試探,而對(duì)IK法的運(yùn)用規(guī)則及區(qū)域性IK法閾值及其函數(shù)關(guān)系的系統(tǒng)研究尚少見。本文以黃河大型灌區(qū)(河套)土壤水鹽監(jiān)測(cè)實(shí)驗(yàn)的水土資源空間變異為案例,重點(diǎn)分析IK方程三項(xiàng)主要結(jié)構(gòu)因子:指示閾值、指示變異函數(shù)和條件概率分布函數(shù)的變化趨勢(shì)、相互關(guān)系和取值規(guī)則,為IK法的應(yīng)用提供參考。1ik平均估計(jì)的確定設(shè)在某實(shí)驗(yàn)區(qū)D上取樣并測(cè)定其品位觀測(cè)值,若該區(qū)域上的邊界門限(即閾值)為Z,則D內(nèi)每一個(gè)樣品點(diǎn)x∈D上定義一個(gè)Z的階梯函數(shù),即指示函數(shù):i(x;Ζ)={1當(dāng)x點(diǎn)上的觀測(cè)值Ζ(x)≤Ζ0當(dāng)x點(diǎn)上的觀測(cè)值Ζ(x)>Ζ(1)D內(nèi)任一區(qū)域A∈D,低于閾值Z的觀測(cè)值Z(x)所占A區(qū)的比例表示為:?(A;Ζ)=1A∫Ai(x;Ζ)dx∈(2)式中?(A;Z)是關(guān)于Z(x)和Z的二元函數(shù),即小于閾值Z的全部i(x;Z),(x∈A)的平均值。指示函數(shù)i(x;Z)在給定閾值Z的條件下服從二項(xiàng)分布,期望值是:E{Ι(x;Ζ)}=1?Ρrob{Ζ(x)≤Ζ}+0?Ρrob{Ζ(x)>Ζ}=Ρrob{Ζ(x)≤Ζ}=F(Ζ)(3)式中F(Z)稱觀測(cè)值Z(x)的分布函數(shù)在閾值Z處的數(shù)值。當(dāng)I(x+h;Z)和I(x;Z)是被矢量h分隔的兩個(gè)隨機(jī)變量,即兩個(gè)指示變換的數(shù)據(jù)時(shí),則指示變異函數(shù)(IndicatorVariogrom)定義為:γΙ(h;Ζ)=12E{[Ι(x+h;Ζ)-Ι(x;Ζ)]2}(4)用指示變異函數(shù)可描述IK的空間變異程度,并可獲得條件概率分布函數(shù)(CCDF)的估計(jì)值:F*[z|(n)]=prob*[Ζ≤z|(n)]=1nn∑α=1i(Ζα;Ζ)?∈[0?1](5)這個(gè)估計(jì)能被接受的前提是:各樣本觀測(cè)值Zα=zα是獨(dú)立的,它們每一個(gè)對(duì)Z的作用是相同的,否則應(yīng)采用不等權(quán)估計(jì):F*[z|(n)]=prob*[Ζ≤z|(n)]=n∑α=1λα(z)?i(za；z)α=1,?,n(6)可以看出指示值i(zα;z)是用Zα-zα單獨(dú)估計(jì)Z時(shí)的ccdf:F*[z|(n)=zα]=prob*[Ζ≤z|Ζ=Ζα]=i(zα;z)α=1,?,n(7)故式(6)可以寫成:F*[z|(n)]=prob*[Ζ≤z|(n)]=n∑α=1λα(z)?F*(z|Ζ=Ζα)(8)上式說明:用多個(gè)觀測(cè)值對(duì)CCDF做出的估計(jì)等于用單個(gè)觀測(cè)值進(jìn)行估計(jì)時(shí)所得CCDF的加權(quán)平均。指示函數(shù)值的方差S2(Z),即指示函數(shù)的基臺(tái)值C0,可從γ*I(h;Z)的圖形求得,故分布函數(shù)F*(z|(n))可從解以下二次方程得到F2(z|(n))-F(z|(n))-C0=0(9)同普通Kriging(OK)法一樣,再運(yùn)用線性無(wú)偏最優(yōu)插值和估計(jì)方差最小的OK要求,結(jié)合指示函數(shù)的有關(guān)條件建立類似普通克立格(OK)方程的IK方程組式(10),可獲得待估域A的IK平均估計(jì)。{n∑β=1λβˉγi(xα,xβ;Ζ)+μ=ˉγi(xα,A;Ζ)n∑α=1λα=1(α=1,2,?,n)(10)用式(10)可求出這些權(quán)系數(shù)λα及未知點(diǎn)的i*(x;Z)估值,即Z(x)≤Z在該點(diǎn)出現(xiàn)的概率和整個(gè)CCDF的估值,可不必剔除特異值而獲得具有一定風(fēng)險(xiǎn)條件下的未知Z(x)的估計(jì)量及空間分布。2案例分析2.1土樣及測(cè)定方法研究區(qū)域設(shè)在黃河河套灌區(qū)內(nèi)的沙壕渠試驗(yàn)區(qū)。試驗(yàn)區(qū)采樣面積為4km2,用手持GPS接收機(jī)定位。以間隔h=200m(采樣點(diǎn)數(shù)為121個(gè),中等田間尺度)的網(wǎng)狀進(jìn)行采樣,采樣分兩層取土,深度為0～20、20～40cm。土樣的測(cè)定分析項(xiàng)目為含水率θ和電導(dǎo)率EC、pH值。其中EC的測(cè)定采用水土比5∶1浸提,電導(dǎo)法測(cè)定水溶性鹽總量。采樣點(diǎn)布置見圖1。2.2數(shù)據(jù)的比較及預(yù)處理IK法關(guān)鍵是閾值的合理確定,目前還沒有從理論上研究水鹽的閾值合理選擇問題。閾值的選擇主要是根據(jù)研究目標(biāo)的專業(yè)門限要求、估計(jì)誤差和允許風(fēng)險(xiǎn)大小進(jìn)行選擇,在缺乏經(jīng)驗(yàn)時(shí),有人建議用中位數(shù)作為初步估計(jì)參考值,經(jīng)作者初步實(shí)踐認(rèn)為:中位數(shù)并非最佳閾值,估計(jì)結(jié)果具有平滑性,且接近OK法,也可能會(huì)產(chǎn)生IK法與OK法無(wú)明顯區(qū)別的誤導(dǎo),為此應(yīng)進(jìn)行不同方案的比較及探索研究。表1水分?jǐn)?shù)據(jù)的偏態(tài)系數(shù)接近于零,表明數(shù)據(jù)分布接近于正態(tài)分布的左右對(duì)稱,變異系數(shù)小反映水分變異性不大;表2鹽分?jǐn)?shù)據(jù)的偏態(tài)系數(shù)均大于零,表明數(shù)據(jù)分布呈右偏態(tài),變異系數(shù)大反映鹽分變異性較大;水鹽數(shù)據(jù)的峰度系數(shù)均大于零,表明數(shù)據(jù)分布比正態(tài)分布高聳而狹窄。按照均值的3倍標(biāo)準(zhǔn)差鄰域判斷,發(fā)現(xiàn)鹽分?jǐn)?shù)據(jù)有明顯特異值存在,0～20cm土層水分有特異值存在,適宜用IK法估計(jì)。對(duì)于大、中尺度樣本應(yīng)根據(jù)不同作物、不同土壤和氣候及地理?xiàng)l件,統(tǒng)一選定合理閾值,如水分閾值可選用田間持水量、凋萎系數(shù)、適宜含水量等;鹽分閾值可采用FAO推薦的作物全生育期平均土壤濃度臨界值、作物能容忍的土壤最大鹽分值等。本文綜合考慮數(shù)據(jù)的分位數(shù)特點(diǎn),水分初步采用20%、22%、25%、27%4個(gè)閾值;鹽分初步采用0.2、0.4、0.6、0.8、1.0、1.2mS/cm6個(gè)閾值進(jìn)行參數(shù)的分析與比較,旨在探索選擇合理水鹽閾值的原則。下文僅列出0～40cm土層土壤的水鹽計(jì)算成果。3閾值與函數(shù)關(guān)系的研究3.1水鹽閾值對(duì)鹽脅迫下自適應(yīng)的影響水鹽指示條件概率的計(jì)算采用ARCGIS9.2軟件的地質(zhì)統(tǒng)計(jì)學(xué)模塊完成。表3中水分預(yù)測(cè)平均條件概率隨著閾值的增加由99.2%減少至18.2%,反映了超過大閾值的概率逐漸減少;RMSE結(jié)果表明隨著閾值的增加,均方根誤差由9.5%增加至22.4%、51.1%,后又減小至39.0%,有先增加、后減小的變化趨勢(shì);表3中鹽分預(yù)測(cè)平均條件概率隨著閾值的增加由82.0%減少至1.8%,反映了超過大閾值的概率逐漸減少;RMSE結(jié)果表明隨著閾值的增加,均方根誤差由38.4%增加至52.2%,后減小至9.4%,也有相似的先增加、后減小的趨勢(shì)。從圖2可以直觀地看出水鹽閾值與指示概率均值的關(guān)系:水分預(yù)測(cè)平均概率隨著閾值的增加而減少,減小的速度不同,超過22%閾值對(duì)應(yīng)的指示概率均值減小速度加大。鹽分預(yù)測(cè)平均概率隨著閾值的增加而減少,減小的速度在0.6mS/cm閾值處變慢。水分預(yù)測(cè)概率均方根誤差隨著閾值的增加,先增加后減少,在25%～26%(接近中位數(shù)25.5%)區(qū)間概率預(yù)測(cè)誤差最大,接近50%,這將會(huì)降低風(fēng)險(xiǎn)評(píng)價(jià)與IK法估計(jì)的可靠性,在20%閾值處概率預(yù)測(cè)誤差最小,可見單從提高概率預(yù)測(cè)精度考慮,水分指示閾值應(yīng)該盡量選擇低于中位數(shù)的偏小值。鹽分預(yù)測(cè)概率均方根誤差隨著閾值的增加,先增加后減少,在0.4mS/cm(接近中位數(shù)0.38mS/cm)處最大超過了50%,也將會(huì)降低風(fēng)險(xiǎn)評(píng)價(jià)與IK法估計(jì)的精度,可見單從提高概率預(yù)測(cè)精度考慮,鹽分指示閾值應(yīng)該盡量選擇超過中位數(shù)的偏大值。圖3—6為不同閾值下土壤水鹽指示概率預(yù)測(cè)圖。分析后發(fā)現(xiàn)高閾值的大概率區(qū)均包含在低閾值的大概率區(qū)內(nèi),其分布位置隨閾值的增加而逐漸減少和集中,可用于不同閾值目標(biāo)下的概率與風(fēng)險(xiǎn)評(píng)價(jià)。分析水分閾值25%、27%的指示概率預(yù)測(cè)圖,含水率超過閾值25%、27%的概率在試驗(yàn)區(qū)中部及左上角較大,含水率超過27%的概率最大值為50%,有4個(gè)明顯大概率區(qū)均包含在超過25%的大概率區(qū)域內(nèi)。分析鹽分閾值0.4、0.6mS/cm相應(yīng)的指示概率預(yù)測(cè)圖,可以了解鹽分積聚的可能性,如電導(dǎo)率超過閾值0.6mS/cm的概率最大為36%,分布在試驗(yàn)區(qū)北部?jī)蓚?cè)及西南角,且均包含在超過閾值0.4mS/cm的大概率區(qū)域內(nèi)。3.2中性關(guān)系變異函數(shù)性變異性函數(shù)的檢測(cè)由表3可見:當(dāng)水分閾值由小逐漸增大時(shí),指示變異函數(shù)值有一種由小到大、再到小的趨勢(shì),轉(zhuǎn)折點(diǎn)在25%處,可見接近中位數(shù)閾值的指示變異函數(shù)值最大,其指示變異函數(shù)結(jié)構(gòu)性最差,所反映的指示概率的空間變化性最大,由此得到的指示概率預(yù)測(cè)誤差最大,這與上文RMSE在閾值接近中位數(shù)處最大的結(jié)果相一致。鹽分指示閾值與指示變異函數(shù)的關(guān)系同水分,指示變異函數(shù)值轉(zhuǎn)折點(diǎn)在0.4mS/cm處。3.3閾值對(duì)結(jié)構(gòu)性分析結(jié)果的影響從表5分析水分指示變異函數(shù)的結(jié)構(gòu)性,閾值22%的指示變異函數(shù)結(jié)構(gòu)性表現(xiàn)為純隨機(jī),因?yàn)槠銫0/(C0+C)=1;結(jié)構(gòu)性較強(qiáng)對(duì)應(yīng)的閾值為20%、25%、27%,隨著閾值增大,結(jié)構(gòu)性有由強(qiáng)到弱,再由弱到強(qiáng)的變化趨勢(shì),為非降函數(shù),其他兩層指示變異函數(shù)的結(jié)構(gòu)性分析具有相同的結(jié)論。從表6分析鹽分指示變異函數(shù)的結(jié)構(gòu)性,無(wú)純隨機(jī)現(xiàn)象,隨著閾值增大,結(jié)構(gòu)性有由強(qiáng)到弱,再由弱到強(qiáng)的變化趨勢(shì),0.4～0.6mS/cm區(qū)間內(nèi)閾值對(duì)應(yīng)的指示變異函數(shù)結(jié)構(gòu)性最差,其指標(biāo)值接近于1。4與其他方法插值對(duì)比本文選用含水率25%、電導(dǎo)率0.6mS/cm分別作為水鹽的指示閾值進(jìn)行IK法區(qū)域水鹽估計(jì),首先借助于IK方程組求解?(A;Z)的線性估計(jì)量?*(A;Z);然后計(jì)算區(qū)域A平均值的估計(jì)值:Z**(A;Ζe)=n∑e=1[?*(A;Ze)×Z*(A;Ze)],最后與OK法、傳統(tǒng)反距離加權(quán)平均法(簡(jiǎn)寫IDW法)的插值結(jié)果做了對(duì)比。圖7—8為土壤水鹽IK法插值圖,可以反映試驗(yàn)區(qū)的水鹽分布條件。水分分布比較連續(xù),有3處明顯濕潤(rùn)區(qū);鹽分在西南、東北及試區(qū)中部(對(duì)角線方向)有明顯條帶狀積聚。從表7可以分析3種方法的水分預(yù)測(cè)均值與實(shí)測(cè)值基本相同,在121個(gè)采樣點(diǎn)上的預(yù)測(cè)平均誤差為1.6%,說明了IK法同樣具有精確的插值功能。IK法的估計(jì)方差、變異系數(shù)最大,但卻小于采樣點(diǎn)對(duì)應(yīng)參數(shù),說明IK法對(duì)水分普通克立格插值平滑效應(yīng)有一定程度的消弱。從表8可以分析3種方法的鹽分預(yù)測(cè)均值與實(shí)測(cè)值基本相同,在121個(gè)采樣點(diǎn)上的預(yù)測(cè)平均誤差為0.230mS/cm,說明了IK法同樣具有精確的插值功能。IK法的估計(jì)方差、變異系數(shù)大于OK法,但卻小于采樣點(diǎn)對(duì)應(yīng)參數(shù),說明IK法對(duì)鹽分普通克立格插值平滑效應(yīng)有一定程度的消弱。5加指示概率估計(jì)的優(yōu)缺點(diǎn)(1)IK法的難度在于閾值的合理選擇,不僅涉及研究目標(biāo)和專業(yè)門限要求,也與研究允許風(fēng)險(xiǎn)概率及變異程度有關(guān),目前尚無(wú)理論分析方法,多數(shù)情況仍然需要試算探求。(2)經(jīng)過模型與參數(shù)分析,發(fā)現(xiàn)水鹽指示Kriging平均概率與指示閾值、指示變異函數(shù)的關(guān)系為:隨著閾值的增加指示Kriging平均概率減少;當(dāng)閾值由小逐漸增大時(shí),指示變異函數(shù)值有一種由小到大,再到小的趨勢(shì)。(3)作者認(rèn)為選擇中位數(shù)作閾值并非最佳估計(jì),接近中位數(shù)閾值的指示變異函數(shù)值最大,其指示變異函數(shù)結(jié)構(gòu)性最差,所反映的指示概率的空間變化性最大,由此得到的指示概率預(yù)測(cè)誤差最大。本