基于指示克里格法的土壤重金屬pb空間分布特征分析

基于指示克里格法的土壤重金屬pb空間分布特征分析

隨著城市化進程和工農(nóng)業(yè)的快速發(fā)展，土壤中積累了大量的重金屬元素，不僅會造成土壤和耕地的污染，還會通過鏈式石柱進入人體，從而損害人類的健康。因此,土壤中重金屬的含量及其空間分布、環(huán)境質(zhì)量、風險評價等研究都受到了人們的廣泛關(guān)注,再加上近年來GIS和地統(tǒng)計理論被引入這一研究中,推動了研究向前發(fā)展,目前這已成為土壤科學的熱點問題之一。然而,人們在土壤重金屬時空變異研究中發(fā)現(xiàn),由于受來源、形態(tài)及控制因素等影響,區(qū)域土壤重金屬往往存在特異值,這些特異值既非分析誤差所致,也非采樣方法等人為誤差所致,卻往往使觀測值偏離正態(tài)分布,影響變異函數(shù)的穩(wěn)健性。目前解決途徑主要有三種:一種是非線性變換(如平方根變換,對數(shù)變換,反正弦變換等),但目前只有對數(shù)變換可以應(yīng)用對數(shù)正態(tài)克里格估值并轉(zhuǎn)換回原尺度,其他變換沒有相應(yīng)的轉(zhuǎn)換方法;另一個途徑是剔除或用估計鄰域和影響系數(shù)等方法處理特異值,但會改變實驗數(shù)據(jù)的原始空間結(jié)構(gòu);第三種途徑是非參數(shù)統(tǒng)計方法或稱無分布方法,常見的如指示克里格法,是條件普通克里格的非參數(shù)形式,可有效抑制特異值對變異函數(shù)穩(wěn)健性的影響。指示克立格法自提出以來,在礦產(chǎn)、石油、水文等領(lǐng)域有廣泛應(yīng)用。目前,使用指示克立格法對水土資源條件概率空間分布的研究較多,而對區(qū)域土壤重金屬含量空間分布估值研究鮮有報道。為此,本文以國家重要糧食生產(chǎn)基地的魯西北禹城市為例,采用指示克里格法對該地區(qū)土壤重金屬Pb含量分布進行估值研究,揭示土壤中Pb含量小于特定閾值的條件概率分布及其含量空間分布特征,并對普通克里格和指示克里格的估值效果進行交叉驗證,旨在為該地區(qū)土壤環(huán)境質(zhì)量評價和重金屬污染防控提供科學依據(jù)。1材料和方法1.1鹽堿土現(xiàn)狀分析禹城市地處山東省西北部,徒駭河中游,地理位置介于E116°22′~116°45′,N36°40′~37°12′,屬暖溫帶半濕潤季風氣候,多年平均氣溫13.1℃,降雨量616mm,蒸發(fā)量1379.5mm,蒸降比2.24,降雨季節(jié)性分布不均,其中6~9月的夏季平均為470.7mm,占全年的76.4%。禹城市屬典型的沖積平原地質(zhì)區(qū),海拔高度介于19.2~27.3m,微地貌比較復(fù)雜,現(xiàn)有耕地5.05×104ha,占全市總面積的51.01%,主要有潮土和潮化鹽土兩大土類,潮土占總面積的97.77%。歷史上,禹城市為黃淮海平原典型的鹽堿土分布區(qū),土壤鹽堿和瘠薄嚴重制約農(nóng)業(yè)生產(chǎn)發(fā)展,經(jīng)過引黃灌溉和多年的鹽堿土整治改良,鹽堿地得到綜合治理,耕地地力得到培育。目前,禹城市已成為國家重要的商品糧生產(chǎn)基地,耕地種植強度和農(nóng)田集約化生產(chǎn)水平大幅提升,耕地環(huán)境質(zhì)量安全越來越受重視。1.2采樣點布設(shè)和樣品采集于2008年5月下旬小麥收獲前后采集土壤樣點。土壤采樣方案為規(guī)則分層抽樣法,采樣過程分為3步:首先,以1:5萬禹城市地形圖、行政區(qū)劃圖、土壤類型圖、地形地貌類型圖等為底圖,通過數(shù)字化獲得研究區(qū)域的邊界圖、土壤圖、地形地貌圖等矢量圖件,在研究區(qū)域邊界圖上進行1.5km×1.5km網(wǎng)格布點;其次,根據(jù)預(yù)設(shè)布點進行野外采樣,實際采樣過程中對于落在村莊、河流、道路等地的樣點,在附近500m內(nèi)調(diào)整,并記錄實際采樣點坐標,共采集396個樣點;最后,在土地利用類型較為復(fù)雜的城郊結(jié)合帶按照500m×500m網(wǎng)格加密采集285個樣點。本研究共計采集有效采樣點681個。每個樣點均按照10m×10m正方形4個頂點和中心點分別采樣,各取表層(0~20cm)土壤約1kg,現(xiàn)場均勻混合后用四分法從中選取1kg土樣作為代表該點的混合樣品。為防止樣品污染,在采樣、樣品保存和樣品處理過程中,避免與金屬器皿直接接觸。采集的土樣帶回實驗室內(nèi)自然風干,磨碎、過100目尼龍網(wǎng)篩后,進行重金屬Pb全量含量的測定,測定方法為石墨爐原子吸收分光光度法,土壤樣品的化學分析按照國家標準(GB15618-1995)執(zhí)行完成。1.3估計累積分布函數(shù)指示克里格方法(IndicatorKriging,IK)是一種非參數(shù)估計方法,其主要目的是對隨機函數(shù)Z(u)在非取樣點的不確定性作出估計。它首先要對數(shù)據(jù)作指示變換,然后估計出待估點處的累積分布函數(shù)(cdf),以此為基礎(chǔ)可解決各種估計和模擬問題。指示克里格方法的具體計算過程與步驟可參考文獻[15-17].本文采用GS+軟件計算各閾值條件下的指示變異函數(shù),計算時以采樣點間的平均間距作為步長X、Y方向滯后距,以22.5°作為方向容許偏差,分別計算0°、45°、90°和135°方向指示變異函數(shù)以及各向同性下的指示變異函數(shù)。利用ArcGIS9.3軟件繪制不同閾值條件下土壤Pb含量的概率分布,并進行空間運算,形成土壤Pb含量的指示克里格空間分布圖。2結(jié)果與分析2.1研究區(qū)表層土壤pb含量分布情況研究區(qū)土壤Pb含量的統(tǒng)計特征值如表1所示?？梢钥闯?土壤Pb含量介于1.18~159.43mg·kg-1,平均值為36.33mg·kg-1,表明研究區(qū)土壤Pb含量整體不高,其平均值略高于自然背景值(≤35mg·kg-1),但不同部位土壤Pb含量差異顯著。從圖1顯示的頻率分布與累積頻數(shù)分布圖來看,土壤Pb含量分布均明顯呈現(xiàn)向左偏倚的特點,在左側(cè)的低值表現(xiàn)出高頻數(shù)分布,右側(cè)由于特異值的存在出現(xiàn)長尾現(xiàn)象。單樣本K-S正態(tài)檢驗和對數(shù)正態(tài)檢驗(p<0.05,2-tailed)顯示,研究區(qū)土壤Pb含量及其對數(shù)轉(zhuǎn)換值均不服從正態(tài)分布(括號內(nèi)為對數(shù)轉(zhuǎn)換值K-S檢驗顯著性水平)。合理采樣數(shù)計算結(jié)果表明,在99%置信水平5%相對誤差條件下,土壤Pb的合理采樣數(shù)目為372,因此研究區(qū)土壤Pb采樣點密度已滿足該置信水平和相對誤差條件的合理采樣數(shù)量要求,即該樣點布設(shè)模式和樣本量能較客觀地反映研究區(qū)土壤Pb的空間變異特征。按照國家GB15618-1995土壤環(huán)境質(zhì)量標準分級可知,土壤全量Pb的44.05%采樣點屬于一類;參照二類土壤執(zhí)行標準,即能保證農(nóng)業(yè)生產(chǎn)和維護人體健康為標準,根據(jù)研究區(qū)土壤調(diào)查,該地區(qū)土壤pH值在7.52~9.15之間,所以選擇國家標準的第二級pH>7.5的閾值為標準進行評價,其余55.95%樣點屬于二類。研究區(qū)表層土壤Pb含量未超過國家標準的限值(≤350mg·kg-1)。因此,研究區(qū)域總體土壤環(huán)境質(zhì)量狀況良好。2.2閾值區(qū)域化測量的方差分析及其擬合參數(shù)取土壤Pb含量0.1~0.9分位數(shù)共9個值作為閾值T,分別為18.64、24.23、28.40、33.24、37.96、41.32、44.48、47.51和51.19mg·kg-1,計算小于各閾值條件的指示變異函數(shù)。計算時X、Y方向滯后距均為500m,分別計算0°、45°、90°和135°方向指示變異函數(shù)以及各向同性下的指示變異函數(shù),并將所得到的指示變異函數(shù)采用球狀模型進行擬合。限于篇幅,圖2僅列出了部分閾值的指示變異函數(shù)及其擬合曲線,其中圖2(a)、圖2(b)、圖2(c)和圖2(d)的閾值分別為0.1、0.3、0.5和0.7分位數(shù),即閾值T分別為18.65、28.40、37.96和44.48mg·kg-1。從各方向指示變異函數(shù)對比可以看出,盡管研究區(qū)域內(nèi)土壤Pb含量在各閾值條件下的概率分布其變化規(guī)律趨于一致,但在不同的方向上存在差異,即區(qū)域化變量表現(xiàn)出一定的各向異性,尤其是在NE45°和NE135°方向上。各閾值條件下的指示變異函數(shù)及其擬合參數(shù)如表2所示。塊金值C0通常表示由實驗誤差和小于實驗取樣尺度引起的變異;基臺值C通常表示系統(tǒng)內(nèi)總的變異;塊金值/基臺值C0/C表示隨機部分引起的空間異質(zhì)性占系統(tǒng)總變異的比例,它反映了土壤屬性的空間依賴性,常被用作研究變量空間相關(guān)的分類依據(jù);變程A反映區(qū)域化變量的相關(guān)程度,變程越大,相關(guān)性越好,越小則相關(guān)性越差。由表2可知,不同的閾值其指示變異函數(shù)擬合參數(shù)差異較大,閾值為0.5分位數(shù)條件下,塊金值C0和基臺值C達到最大;C0/C在49.51%~80.85%之間,隨著閾值的不斷升高,C0/C呈逐漸增加趨勢,其空間依賴性降低;閾值為0.1~0.4分位數(shù)條件下,指示變異函數(shù)呈中等變異性,閾值大于0.5分位數(shù)其指示變異函數(shù)均呈強變異性,空間結(jié)構(gòu)性減弱,這一點可從表2中擬合殘差RSS和決定系數(shù)r2的對比可知。2.3土壤pb的空間估值利用表2中的指示變異函數(shù)及其參數(shù)進行克里格估值,獲得的土壤Pb含量小于各閾值T的條件概率空間分布如圖3所示,其中圖3a~圖3i分別顯示閾值為0.1~0.9分位數(shù)的條件概率分布。為了便于根據(jù)條件累積分布函數(shù)利用指示克里格法對研究區(qū)域各部位的土壤Pb含量進行估算,須將研究區(qū)域各閾值的條件概率分布劃分為規(guī)則的網(wǎng)格單元。綜合考慮研究區(qū)土壤采樣方案與實際指導(dǎo)性,將各條件概率分布圖轉(zhuǎn)換為500m×500m柵格數(shù)據(jù),整個研究區(qū)域共有3922個柵格參與計算。根據(jù)表2中各設(shè)定閾值和圖3的條件累積概率分布圖,利用式(6)計算各網(wǎng)格單元土壤Pb含量指示克里格估計值,獲得的研究區(qū)土壤Pb含量空間分布見圖4(圖4b)。圖4分別顯示了研究區(qū)土壤Pb含量分別采用普通克里格(OK)和指示克里格(IK)估值形成的空間分布圖。從空間分布特點來看,普通克里格法和指示克里格法分別形成的空間分布具有相似的格局,在局部平滑區(qū)和斑塊部位具有一致性;從分布規(guī)律來看,研究區(qū)土壤Pb含量高值區(qū)主要集中在各鄉(xiāng)鎮(zhèn)駐地附近的農(nóng)田,如圖4中,安仁鎮(zhèn)、倫鎮(zhèn)、張莊鎮(zhèn)等鄉(xiāng)鎮(zhèn)附近的土壤Pb含量均呈明顯的斑塊狀分布;從整個研究區(qū)域統(tǒng)計結(jié)果來看,由普通克里格法獲得的土壤Pb空間分布圖中(圖4a),Pb含量符合土壤環(huán)境質(zhì)量二級標準值(35~350mg·kg-1)的分布面積占研究區(qū)49.06%,該區(qū)域平均值為42.37mg·kg-1,而指示克里格法獲得的分布圖中(圖4b),二級標準分布面積占62.57%,其平均值為47.66mg·kg-1,普通克里格法的平滑效果強于指示克里格法,這是由于在土壤Pb含量空間結(jié)構(gòu)分析和空間估值過程中,盡管普通克里格和指示克里格都限制了特異值的影響范圍,但指示克里格仍然以概率的方式考慮了特異值的存在?？傮w上,研究區(qū)域土壤環(huán)境質(zhì)量良好,但對于土壤Pb含量相對較高的城郊結(jié)合部、工礦企業(yè)集中區(qū)和鄉(xiāng)鎮(zhèn)駐地附近,應(yīng)重視土壤環(huán)境質(zhì)量的監(jiān)測監(jiān)控,加強土壤重金屬含量的時間變異性研究。由表3列出了指示克里格和普通克里格對研究區(qū)土壤Pb含量空間估值的交叉驗證結(jié)果。可以看出,指示克里格法的平均誤差ME(MeanError)、均方根誤差RMSE(RootMeanSquareError)和平均標準誤ASE(AverageStandardEerror)均要小于普通克里格法,指示克里格法的均方標準誤RMSSE(Rootmean-SquareStandardizedError)較普通克里格法更接近1,且指示克里格估計值與實測值間的相關(guān)系數(shù)亦明顯高于普通克里格法,這表明指示克里格對土壤Pb含量空間估值的精度優(yōu)于普通克里格,指示克里格法較普通克里格法更適用于原始數(shù)據(jù)呈非正態(tài)分布且特異值對估值效果有較大影響的土壤屬性。因此,指示克里格法在該區(qū)域尺度范圍內(nèi)進行土壤Pb含量空間估值是可行的。3土壤pb空間分布增長中的應(yīng)用參數(shù)地質(zhì)統(tǒng)計學針對不同的土壤屬性提出了不同的研究方法。在平穩(wěn)條件下,可使用簡單克立格和普通克立格;在土壤屬性具漂移或趨勢變化特點時,可采用泛克立格或趨勢克立格;若各土壤屬性之間具有一定相關(guān)性時,則可應(yīng)用各種協(xié)克立格。參數(shù)地質(zhì)統(tǒng)計學往往對樣點均勻分布和樣本量充足有一定要求,且數(shù)據(jù)須服從正態(tài)分布。有些土壤屬性由于受區(qū)域內(nèi)多種因素的控制與影響,其分布常出現(xiàn)拖尾現(xiàn)象,如本文的土壤Pb含量,這種分布特征會影響數(shù)據(jù)的整體分布規(guī)律及估值的精度,采用普通的克立格法往往穩(wěn)健性較差。目前常采用的解決方法包括特異值剔除、數(shù)據(jù)平滑、進行對數(shù)轉(zhuǎn)換或某些線性轉(zhuǎn)換等,所有這些方法都是根據(jù)數(shù)據(jù)的分布特點進行變換,使數(shù)據(jù)分布特征滿足估值方法要求。然而在實際應(yīng)用中,被研究的土壤屬性常存在樣本數(shù)量有限且分布不均的現(xiàn)象,特異值的存在使得變異函數(shù)不穩(wěn)健,數(shù)據(jù)的分布假設(shè)無法檢驗等一系列問題使得參數(shù)地質(zhì)統(tǒng)計學的應(yīng)用受到一定的限制。為此Journel于1983年提出了指示克立格法,它不要求數(shù)據(jù)服從正態(tài)分布,且不受特異值的影響,是處理有偏數(shù)據(jù)的有力工具。本文正是針對土壤Pb含量的特異值和非正態(tài)特點,采用指示克里格法對區(qū)域土壤Pb含量空間分布進行估值,并與普通克里格法的估值精度進行對比。張新宇等(2006)對阿舍勒銅礦體進行了指示克里格可視化儲量估算,表明基于指示克里格法的儲量計算能夠取得較理想的預(yù)測結(jié)果;尚曉穎等(2007)利用普通克里格和指示克里格對水環(huán)境模擬參數(shù)的空間分布進行估值,結(jié)果顯示指示克里格較普通克里格更好地描述水環(huán)境模擬參數(shù)隨機場的空間分布規(guī)律;趙玉杰等(2009)探討了普通克里格法及指示克里格法在農(nóng)產(chǎn)品產(chǎn)地禁產(chǎn)區(qū)篩選中應(yīng)用的可行性,結(jié)果表明指示克里格法的誤判率顯著低于普通克里格法,指示克里格法較普通克里格法更適合于研究區(qū)水稻禁產(chǎn)區(qū)的篩選;這與本文的研究結(jié)果也是一致的,因而指示克里格法為呈非正態(tài)分布的水土屬性的空間分布估值提供了一種新的可行性方法。此外,基于均值和中位值的穩(wěn)健變異函數(shù)估算亦是處理有偏分布數(shù)據(jù)的有效手段,它可以在不必剔除或修改特異值的條件下進行分析,常用的方法有克瑞西-豪金斯法、中位調(diào)節(jié)法、中位數(shù)中位偏差法等,目前已在水文水資源、土壤屬性空間變異研究中得到一定程度應(yīng)用。在下一階段研究中,作者將探討非參數(shù)地質(zhì)統(tǒng)計法(指示克里格)與變異函數(shù)穩(wěn)健估計法應(yīng)用效果的對比。4土壤環(huán)境質(zhì)量1)研究區(qū)表層土壤Pb含量總體較低,平均值為36.33mg·kg-1,但不同部位土壤Pb含量差異較大,其頻數(shù)分