多元指示克立格法在黃土丘陵溝壑區(qū)的應(yīng)用

多元指示克立格法在黃土丘陵溝壑區(qū)的應(yīng)用

干旱和半干旱地區(qū)耕地的土壤水分和養(yǎng)分是限制作物生長的主要因素之一。監(jiān)測和預(yù)測土壤水分和養(yǎng)分的時空分布對指導(dǎo)這些地區(qū)的節(jié)水灌溉、堿土壤改良和農(nóng)業(yè)可持續(xù)發(fā)展具有重要意義。由于“鹽隨水來,鹽隨水走”,土壤水分和鹽分的時空分布及其動態(tài)并不是相互獨立的過程,而是具有一定的協(xié)同相關(guān)性,即為協(xié)同區(qū)域化變量,因此,國內(nèi)外的學者們也應(yīng)用協(xié)同克立格法對土壤水分和鹽分的空間分布進行了估計。然而,無論是常用的普通克立格,還是單元指示克立格,或協(xié)同克立格,只能給出某一時期土壤水分或鹽分的空間分布等值線圖,卻不能同時從土壤的水分、鹽分含量及其它特性幾方面綜合評價土壤的質(zhì)量,如果再加上時域特征,分析起來就更為困難了,這就需要一種可以將這些時空變量整合在一起形成一個綜合變量,并考慮綜合變量空間結(jié)構(gòu)性的方法,多元指示克立格正是這樣一種方法。多元指示克立格法是利用多元指示變換(Multiple-VariableIndicatorTransform)整合各種指示變量(例如,水土資源的物理性質(zhì)、化學性質(zhì)、生物性質(zhì)及其它描述其質(zhì)量的變量都可以作為水土資源質(zhì)量指示變量),利用非參數(shù)地質(zhì)統(tǒng)計學——指示克立格方法(IndicatorKriging)估計滿足某一閾值集合(由各個指示變量的閾值組成)條件下綜合指示變量(CombinedIndicator)條件概率的一種新方法。該方法1993年由Smith等提出,此后,Halvorson等應(yīng)用該方法分析了灌木-草原生態(tài)系統(tǒng)中土壤特性的空間變異模型,但將該方法應(yīng)用于較大面積不同時期土壤水鹽特性的分析少有報道。本文應(yīng)用多元指示克立格法,結(jié)合單元指示克立格對河套灌區(qū)約55hm2農(nóng)田的兩個時期土壤水分和鹽分的空間變異性進行了綜合分析和評價,這是該法在中尺度農(nóng)田的土壤質(zhì)量評價中應(yīng)用的一個很好的例子,對應(yīng)用多元地質(zhì)統(tǒng)計學理論綜合各種影響因素及其空間分布,從而制定科學的農(nóng)業(yè)管理措施,具有一定的借鑒意義和參考價值。1基本原則和方法1.1傳統(tǒng)克立格估計算法土壤具有高度的空間異質(zhì)性,其特性的觀測數(shù)據(jù)中經(jīng)常會出現(xiàn)一些特異值,它們比所研究的全部觀測值的算術(shù)平均值或中位數(shù)值要高(或低)得多,但并非誤差所致,特異值在全部觀測值中雖只有少部分,但它對全部觀測值的統(tǒng)計結(jié)果影響很大,往往使觀測值偏離正態(tài)分布,造成變異函數(shù)的不穩(wěn)健,從而影響克立格估計結(jié)果。通常的解決辦法是通過非線性變換(如,平方根變換,對數(shù)變換,反正弦變換等)“平滑”這些特異值,對變換后的數(shù)據(jù)進行(地質(zhì))統(tǒng)計分析,最后再進行逆變換。但經(jīng)非線性變換后的數(shù)據(jù)不能直接應(yīng)用普通克立格等線性估計技術(shù),除非有辦法將克立格估計的數(shù)據(jù)轉(zhuǎn)換回原尺度,否則就會產(chǎn)生有偏估計。但目前只有對數(shù)變換可以應(yīng)用對數(shù)正態(tài)克立格加以轉(zhuǎn)換,其他變換沒有相應(yīng)的轉(zhuǎn)換方法。另一個途徑是剔除特異值或用估計鄰域和影響系數(shù)等方法處理特異值,但這些方法都在一定程度上破壞了變量空間變異的真實信息。第三種途徑就是應(yīng)用指示克立格法進行分析。指示克立格是條件普通克立格的非參數(shù)形式,因為它能抑制特異值對變異函數(shù)的影響,因此是處理有偏數(shù)據(jù)的有力工具。假設(shè)在區(qū)域D上采樣測定了某變量,若該變量在區(qū)域D上的閾值或邊界值為z,則在x∈D處可定義如下階梯函數(shù)(即指示函數(shù)):I(x?z)={?1Z(x)≥z?0Z(x)＜z(1)在給定的閾值Z的條件下,記Z(x)≥z(即I(x,z)=1)的概率為F(z)=Prob[I(x?z)=1](2)則指示函數(shù)I(x,z)的期望值為E{I(x?z)}=1×F(z)+0×[1?F(z)]=F(z)(3)式(3)表明x處Z(x)≥z(即I(x,z)=1)的概率等于指示變換值的平均值。當I(x+h,z)和I(x,z)為被矢量h分隔的兩個指示變換數(shù)據(jù)點時,則指示變異函數(shù)可定義為γI(h?z)=12E{[I(x+h?z)?I(x?z)]2}(4)可見用式(4)求得的指示變異函數(shù),可通過普通克立格估計得到的未知點的I(x,z)值,即為Z(x)≥z在該點出現(xiàn)的概率,即F?(Z)=∑α=1nλαI(xα)(5)式中F*為Z(x)≥z的概率F的克立格估計;n是在待估點x的鄰域中的有效點數(shù);xα是第α個采樣點的位置;λα是克立格權(quán)系數(shù)。用式(5)可以估計Z(x)≥z的條件概率圖。式(4)表明,指示變異函數(shù)是用指示函數(shù)計算得到的,計算中所有的指示變換值1(或0)對變異函數(shù)的貢獻都是一樣的,因此,不受特異值的影響。須注意的是,指示函數(shù)除式(1)外,還有其它形式,一般是出現(xiàn)預(yù)期結(jié)果變換為1,否則為0,選用何種形式,按具體問題而定。1.2指標閾值集的確定在水土資源的分析評價中,影響決策的指標往往不只一個,如何將眾多的指標整合為一個綜合指標,是非常重要的。多元指示變換是實現(xiàn)這一目標的有用的方法,其步驟為:①根據(jù)問題的要求(如,農(nóng)業(yè)生產(chǎn)管理或環(huán)境要求等),為參與分析的所有指標(即變量)獨立地確定一個閾值,每個閾值可是一個臨界值,也可以是一個范圍,所有指標的閾值構(gòu)成閾值集,閾值原則上是可任選的;②為每個指標確定指示函數(shù),并根據(jù)指示函數(shù)對相應(yīng)的指標進行二態(tài)指示變換;③進行指標整合,確定綜合指標,即當某一點上有m(m≤n,n為變量數(shù))個指標的指示變換值為1,則該點的綜合指示變換值為1,否則為0。經(jīng)多元指示變換后得到的由綜合指示變換值(即1或0)組成的數(shù)據(jù)集即為綜合指標,可以用作評價相應(yīng)采樣點上水土資源質(zhì)量優(yōu)劣的指標。然而,這一方法更大的價值在于可以利用綜合指標的指示變異函數(shù)(求法與單變量指示變異函數(shù)相同),通過普通克立格法估計未采樣點的值,進而繪制整個研究區(qū)域上所有或部分指標(即變量)同時滿足相應(yīng)標準(或要求)的概率圖。由于這種方法可以整合多個指標(即變量),得到綜合克立格概率分布圖,故稱之為多元指示克立格(MVIK)。多元指示克立格雖然只是單元指示克立格的簡單外延,卻具有更廣泛的用途。該方法本身并不定義水土資源的質(zhì)量,而是整合那些代表水土資源質(zhì)量的標準(即閾值),因而根據(jù)要解決的問題選擇合適的標準很關(guān)鍵。2研究區(qū)土壤水分、土壤樣品采集與分析試驗區(qū)設(shè)在黃河河套灌區(qū)長勝試驗站,該地區(qū)氣候干旱,地下水埋深淺,水平排水不暢,是我國鹽漬土發(fā)育的典型地區(qū)。試驗區(qū)是正在改良的鹽漬化土地,試驗區(qū)總面積為1200hm2,重點試驗區(qū)為規(guī)則的長方形,面積為55hm2,1998年秋澆前(9月中旬)與1999年夏灌前(4月下旬)沿東西、南北方向每隔100m用剖面法采集1m深度(分5層)的土樣,樣品分析成果包括土壤全鹽量(占干土重%)、土壤水分含量、EC值、pH值以及土壤8個主要離子含量。本文以重點試驗區(qū)秋澆前與夏灌前0～100cm內(nèi)土壤全鹽量及土壤水分含量的平均值為例,應(yīng)用指示克立格法分析土壤水鹽特性的空間分布,討論多元指示克立格法的適用性。3結(jié)果分析3.1夏灌前水分變異情況表1為土壤水分和鹽分的統(tǒng)計特征值,水分的變異性明顯小于鹽分。若按一般對變異系數(shù)值Cv的評估,當Cv≤0.1時,為弱變異性,0.1≤Cv≤1.0為中等變異性,則夏灌前水分屬弱變異性,秋澆前水分屬中等偏弱變異,鹽分都屬于中等變異性。水分的偏度都較接近于0,峰度也較接近于3,所以水分接近正態(tài)分布,但鹽分的峰度和偏度都很大程度偏離了3和0,即嚴重偏離正態(tài)分布。實際上,鹽分的頻率分布圖(圖1只給出了夏灌前鹽分的頻率分布圖)明顯向左偏倚,不符合正態(tài)分布,經(jīng)對數(shù)變換后,大致接近正態(tài)分布,但右側(cè)由于特異值存在仍有長尾現(xiàn)象。3.2指揮克立格分析(1)夏灌前土壤鹽分的行為統(tǒng)計特征(表1和圖1)表明,土壤鹽分不僅呈非正態(tài)分布,而且存在特異值,兩者都會影響變異函數(shù)的穩(wěn)健性,進而影響克立格估計結(jié)果;此外,盡管各種克立格方法(包括指示克立格和普通克立格等)能給出單個變量,即某一時期土壤鹽分的空間分布圖,但有時人們更感興趣的不是某一點處鹽分的具體含量,而是大于(或小于)的某一閾值,或處于某一范圍內(nèi)的鹽分含量在空間上的分布概率。指示克立格法不僅能削弱有偏分布和特異值對變異函數(shù)及估計結(jié)果的影響,而且能有效地估計區(qū)域化變量的分布概率,成為解決這類問題的有力工具。以夏灌前鹽分為例,進行單元指示克立格分析。其計算的第一步(也是比較關(guān)鍵的一步)是選擇閾值,對于夏灌前鹽分本文選0.3%,且置大于0.3%的含鹽量為1,否則為0。選0.3%為閾值的原因是:河套灌區(qū)鹽漬化土地夏灌前未播種,這時如知道哪些是鹽漬化發(fā)生的高概率區(qū),哪些是低概率區(qū),就可為種植作物種類的選擇提供依據(jù)(即鹽漬化發(fā)生的低概率區(qū)種植耐鹽性較差的玉米等作物,反之種植耐鹽性較好的小麥、向日葵等);而按一般鹽土的分級,1m深度內(nèi)土層全鹽量小于0.2%或0.3%為非鹽漬化土壤,否則為鹽漬化土壤,考慮本試驗區(qū)為正在改良的鹽漬化土壤,因而選較大的0.3%為閾值。圖2和圖3分別為夏灌前土壤鹽分數(shù)據(jù)的指示變異函數(shù)圖和指示克立格估計的鹽分含量大于0.3%的概率分布圖。圖4是夏灌前土壤鹽分的指示克立格估值分布圖。指示克立格估計值分布圖(圖4)實際是根據(jù)條件概率計算得到的鹽分含量期望值的分布圖,由于該圖的估計不受特異值和數(shù)據(jù)有偏分布影響,無需處理特異值,也無需對數(shù)據(jù)進行正態(tài)變換,能避免普通克立格估計中由于這些處理而引起的估計誤差。實際上,交叉驗證的結(jié)果(略)表明,用指示克立格法對夏灌前各實測點鹽分進行估值,所得誤差均方比普通(對數(shù)正態(tài))克立格估計(即先進行對數(shù)正態(tài)變換,再用影響系數(shù)法處理特異值,然后計算變異函數(shù),最后用對數(shù)正態(tài)克立格法估值)的誤差均方縮小27%。另一方面,由于普通克立格強烈的“平滑效應(yīng)”,使得所得鹽分分布圖(未給出)中左側(cè)的重度鹽漬化(全鹽量大于0.7%)未完全反映出來,而實際試驗區(qū)的西南角(圖4中的左上角)與鹽荒地相鄰,鹽漬化較為嚴重。說明指示克立格所得鹽分分布圖較普通克立格更能反映實際情況。但有一點須承認:指示克立格也有一定程度的平滑效應(yīng),只是不如普通克立格明顯。從圖3能清楚地看出夏灌前鹽分含量大于0.3%的高概率區(qū)域占很大一部分,對這些區(qū)域可選擇較耐鹽的作物種植,且在夏灌時(播種前)對鹽分含量大的高概率區(qū)(圖3中的深色部分)要適當增加灌水量以沖洗鹽分。對含鹽量超過0.3%的低概率區(qū)(圖3中的淺色部分)為爭取農(nóng)時可酌情考慮夏灌前提早播種。(2)夏灌前鹽漬化情況的調(diào)查與防治上面是只考察一個時期土壤鹽分的情形,如果要同時考察兩個或兩個以上時期土壤的鹽分,或同時考察土壤水鹽的分布時,就需要用到多元指示克立格。這里所謂的多元指示克立格是(單元)指示克立格的簡單外延,為的是將多個診斷指標(即變量)整合為一個指標,以便于分析和估值。下面以秋澆前和夏灌前兩個時期的鹽分為指標診斷試驗區(qū)鹽分的分布情況,并說明多元指示克立格具體計算步驟:①分別為每個指標選擇各自的閾值;本例兩個時期的土壤鹽分閾值都取0.4%,原因是:對于河套灌區(qū)夏灌前和秋澆前是兩個較特殊的時期,秋澆前反映的是經(jīng)整個作物生長期人工灌溉干預(yù)下土壤水分和鹽分的分布狀況,夏灌前非飽和帶中土壤水分和鹽分已經(jīng)完成了凍融季節(jié)的再分配過程,基本代表冬春自然過程影響下的分布狀態(tài),如果某區(qū)域兩個時期含鹽量都較大,則基本上可以認為該區(qū)域鹽漬化程度較嚴重;而按一般鹽土的分級,1m深度內(nèi)土層全鹽量大于0.4,則表明鹽漬化已達中度以上,中度鹽漬化以上的土壤中作物正常生長就會受到抑制;②為每個指標建立指示函數(shù),按指示函數(shù)實施指示變換;本例兩指標選相同的指示函數(shù)形式,即式(1);③將多個指標合并為一個綜合指標(或多個指示變換值合并為一個綜合指示變換值),即當某采樣點上m(m≤n,n表示變量數(shù),m依問題不同而定)個指示變換值如果為1,則綜合指標(或綜合指示變換值)為1,否則綜合指標為0;本例m=n=2,即采樣點上兩個時期的鹽分含量都大于0.4%時,綜合指標為1,否則為0;④計算綜合指標的指示變異函數(shù)(圖5(a));⑤進行普通克立格插值,得到綜合指標的概率圖(圖5(b))。圖5(b)中顏色較深的部分即為秋澆前和夏灌前兩個時期鹽分含量超過0.4%的高概率區(qū),即這些區(qū)域鹽漬化達中度以上的概率較大,需要采取適當?shù)拇胧┻M行鹽堿化防治。圖6為秋澆前水分(大于中值21.3%變換為1,否則為0)和鹽分含量(小于0.4%變?yōu)?,否則為0)的綜合指示變異函數(shù)圖與概率圖。圖6(b)的深色部分表示土壤水分含量較大(大于中值21.3%)、鹽分含量較小(小于0.4%)的高概率區(qū),所以,在這些區(qū)域秋季的儲水灌溉(目的是造墑壓鹽)灌水定額可以低于其它區(qū)域,從而節(jié)約灌水量,減小地下水位抬高的程度。制約作物生長的土壤特性不僅是水分和鹽分,還有各種養(yǎng)分含量、pH值等,多元指示克立格允許同時考慮這些因素,繪制綜合概率圖,可以為土壤或其它資源環(huán)境區(qū)域信息的合理化管理決策提供有力依據(jù),它為綜合、全面地評價土壤及其它資源質(zhì)量提供了新的方法。4多變量量化方法在水土資源特性的數(shù)據(jù)分析中,往往由于特異值的存在使數(shù)據(jù)偏離正態(tài)分布,而不能得到穩(wěn)健的變異函數(shù),影響克立格估計(產(chǎn)生有偏估計)或條件模擬結(jié)果的精度。指示克立格法是解決這類問題的有效工具。本例中,夏灌前土壤鹽分近似符合對數(shù)正態(tài)分布,但存在特異值,經(jīng)指示克立格估值,其結(jié)果的誤差均方比普通(對數(shù))克立格縮小27%,且指示克立格估計的分布圖比普通克立格更能反映真實情況。特別是當變量存在臨界值或某一合理取值范圍時,指示克立格估計的條件概率圖能給出各點取值的風險概率(如本文的夏灌前