長(zhǎng)三角平原區(qū)縣域土壤磷素流失風(fēng)險(xiǎn)及其空間不確定性的快速評(píng)估

陳 劍,瞿明凱,王 燕,萬(wàn)夢(mèng)雪,黃 標(biāo),趙永存

1 中國(guó)科學(xué)院南京土壤研究所土壤環(huán)境與污染修復(fù)重點(diǎn)實(shí)驗(yàn)室, 南京 210008 2 中國(guó)科學(xué)院大學(xué), 北京 100049

土壤磷素是作物生長(zhǎng)的主要營(yíng)養(yǎng)元素,同時(shí)也是水體富營(yíng)養(yǎng)化的關(guān)鍵限制因子[1- 5]。由于農(nóng)田土壤磷素具有較強(qiáng)的空間變異性,而且影響磷素流失的環(huán)境因子較多,因此對(duì)其流失進(jìn)行控制和治理難度較大,需耗費(fèi)大量的人力、物力和財(cái)力[6]。許多研究結(jié)果顯示,在土壤磷素向水體輸出的過(guò)程中,少數(shù)景觀單元輸出的磷素往往占了整個(gè)流域污染負(fù)荷的大部分,這些景觀單元?jiǎng)t被稱為土壤磷素流失的關(guān)鍵源區(qū)[7](Critical Source Areas, CSAs)。而對(duì)土壤磷素污染展開(kāi)風(fēng)險(xiǎn)評(píng)估和CSAs的劃定或識(shí)別,是高效治理污染的關(guān)鍵[8]。在眾多的評(píng)估方法中,磷指數(shù)法(Phosphorus Index,PI)因其簡(jiǎn)便實(shí)用,且評(píng)估結(jié)果具有較好可視性的特點(diǎn)而得到了廣泛的運(yùn)用。但傳統(tǒng)的PI在評(píng)估過(guò)程中,對(duì)各影響因子提前分等定級(jí),具有一定的人為主觀性。而且在縣域尺度下的多數(shù)評(píng)估未考慮研究區(qū)域特征,評(píng)估要求的數(shù)據(jù)量大,過(guò)程復(fù)雜。

對(duì)于現(xiàn)有的PI法評(píng)估,往往未考慮土壤磷素具有較強(qiáng)的空間變異性,且忽視了預(yù)測(cè)結(jié)果的空間不確定性。區(qū)域土壤磷素濃度水平的高低往往直接影響其流失風(fēng)險(xiǎn)程度[9- 14]。土壤全磷含量通常用于衡量區(qū)域土壤磷素的本底水平[15- 16]。地統(tǒng)計(jì)學(xué)模型通常被用于區(qū)域土壤屬性的空間預(yù)測(cè)。然而,基于有限樣本點(diǎn)預(yù)測(cè)得到的土壤全磷空間分布必然存在一定的空間不確定性[17]。傳統(tǒng)的地統(tǒng)計(jì)學(xué)插值模型只能得到唯一的預(yù)測(cè)結(jié)果,故不能反映這一不確定性特征。而且,區(qū)域土壤數(shù)據(jù)集可能含有一定數(shù)量的離群值,這些往往會(huì)對(duì)傳統(tǒng)地統(tǒng)計(jì)學(xué)模型預(yù)測(cè)精度產(chǎn)生一定影響。構(gòu)建穩(wěn)健不確定性模擬方法以去除離群值對(duì)模擬結(jié)果的影響,且能產(chǎn)生多個(gè)可能的空間模擬實(shí)現(xiàn)。這些模擬實(shí)現(xiàn)可以輸入土壤磷素流失風(fēng)險(xiǎn)評(píng)估模型,進(jìn)一步評(píng)估其流失風(fēng)險(xiǎn)及相應(yīng)的不確定性。除了區(qū)域土壤全磷含量,土壤磷素的流失風(fēng)險(xiǎn)還與其他環(huán)境因子密切相關(guān),如距受納水體距離、地表徑流潛力、土壤侵蝕潛力、磷肥施用量等[18-28]。

磷肥的施用是土壤磷素的一個(gè)重要外源輸入,縣域尺度下磷肥的施用量主要和土地利用類型相關(guān)；雖然磷素在不同性質(zhì)的下墊面遷移過(guò)程中衰減機(jī)制尚具有較大的差異[29],但磷素在遷移過(guò)程中不斷被截留和稀釋,距離受納水體較遠(yuǎn)的磷素源區(qū)比距離較近的源區(qū)貢獻(xiàn)量小,因此距離因子可在一定程度上表征土壤磷素在遷移過(guò)程中隨著距離衰減的規(guī)律；縣域尺度下地表徑流量主要由徑流系數(shù)決定,徑流系數(shù)與土地利用類型密切相關(guān)[30]。由于長(zhǎng)三角農(nóng)田水利設(shè)施完善,排水溝渠及擋土墻可有效防治土壤侵蝕,故本研究針對(duì)土壤侵蝕不予考慮。

本研究以長(zhǎng)三角典型縣域金壇為研究區(qū),具體目的在于：(1)采用經(jīng)穩(wěn)健處理的序貫高斯模型(RSGS)模擬金壇區(qū)土壤全磷可能的空間分布格局,并采用獨(dú)立驗(yàn)證點(diǎn),與傳統(tǒng)序貫高斯模擬模型(SGS)的模擬精度進(jìn)行對(duì)比,以驗(yàn)證RSGS在提高磷素空間模擬精度方面的有效性；(2)結(jié)合相關(guān)環(huán)境因子(即磷肥施用量、距受納水體距離、地表徑流潛力)構(gòu)建快速磷指數(shù)(RPI)評(píng)估模型；(3)評(píng)估長(zhǎng)三角典型縣域土壤磷素的流失風(fēng)險(xiǎn)及其空間不確定性,并識(shí)別土壤磷素流失的關(guān)鍵源區(qū)。本研究以期從“不確定性”角度為區(qū)域土壤磷素的調(diào)控,劃定更符合實(shí)際限制條件的管控區(qū)域。

1 材料與方法

1.1 研究區(qū)域概況

本研究選取長(zhǎng)江三角洲典型農(nóng)業(yè)區(qū)金壇區(qū)(31°33′42″—31°53′22″ N,119°17′45″—119°44′59″ E)作為研究區(qū)域(圖1)。其地處江蘇省南部,茅山東麓,是南京、上海和杭州三角地帶的中心。全區(qū)總面積約975平方公里,平均海拔13m。金壇區(qū)內(nèi)水網(wǎng)密布,水域面積達(dá)194.22km2,境內(nèi)有大小河流216條,總長(zhǎng)512km。東南部的洮湖(長(zhǎng)蕩湖)是江蘇省十大淡水湖之一。金壇區(qū)屬北亞熱帶季風(fēng)區(qū),雨熱同季,四季分明；雨量充沛,年平均降雨量為1063.5 mm；日照充足,年平均氣溫為15.3 ℃,無(wú)霜期長(zhǎng)達(dá)228d,年平均濕度為78%。區(qū)內(nèi)農(nóng)用地類型單一,主要為水田和旱地,農(nóng)田水利工程完善,排水溝渠及擋土墻可有效防治土壤侵蝕,但密集的溝渠亦是地表徑流磷素流失的快速通道,有利于土壤磷素的遷移。金壇區(qū)部分區(qū)域化肥施用量大,全年每畝農(nóng)田化肥施用量達(dá)112.5 kg,其中水稻全年每畝施肥量約是全國(guó)平均水平的3.7倍[31]。該區(qū)主要糧食作物為水稻和小麥,是太湖地區(qū)重要的糧食生產(chǎn)基地之一。

1.2 評(píng)價(jià)因子的來(lái)源與處理

(1)土壤全磷

金壇區(qū)內(nèi)共采集土壤樣點(diǎn)259個(gè)(圖1),采樣時(shí)間為2016年5月底(即小麥?zhǔn)崭詈笾了驹苑N前)。采集每個(gè)樣點(diǎn)周圍100 m2范圍內(nèi)4—5處表層(0—20 cm)土樣,均勻混合后縮分至1—2 kg用于實(shí)驗(yàn)室分析,每個(gè)樣本點(diǎn)采用GPS定位并記錄其詳細(xì)信息。所有樣本在室溫(20—22 ℃)風(fēng)干后去除雜物及碎石,磨細(xì),過(guò)100目篩,收集50 g土樣用于土壤全磷的測(cè)定。土壤全磷采用HClO4-H2SO4消煮,鉬銻抗比色法測(cè)定。具體分析方法參見(jiàn)魯如坤[32]。

圖1 金壇區(qū)土地利用類型及土壤采樣點(diǎn)分布圖Fig.1 Map of study area, land-use types, and soil sampling sites in Jintan District

(2)環(huán)境因子

用于評(píng)估長(zhǎng)三角典型縣域金壇區(qū)土壤磷素流失風(fēng)險(xiǎn)的三個(gè)重要環(huán)境因子為磷肥施用量因子、采樣點(diǎn)距受納水體距離因子及地表徑流潛力因子。磷肥施用量以2014年、2015年和2016年金壇區(qū)的統(tǒng)計(jì)年鑒為依據(jù),以三年磷肥施用量平均值計(jì)算,其中旱地磷肥施用量以水田的2.5倍計(jì)[33],其他土地利用類型均考慮為無(wú)肥區(qū),單位面積的磷肥施用量因子由金壇區(qū)磷肥施用總量和水田、旱地的面積比來(lái)獲取；采樣點(diǎn)距受納水體的距離則通過(guò)金壇區(qū)內(nèi)水系矢量數(shù)據(jù)(1∶25萬(wàn))采用近鄰分析方法得到；地表徑流潛力可采用徑流深R(cm)來(lái)度量,其計(jì)算公式如下：

R=p·α

(1)

式中,p表示年降雨量(cm)；α表示年徑流系數(shù)。

本研究中,考慮到縣域尺度內(nèi)降雨量區(qū)域差異相對(duì)較小,地表徑流主要受年徑流系數(shù)α的影響。因此,本研究根據(jù)金壇區(qū)的特點(diǎn)和相關(guān)研究來(lái)確定該區(qū)域不同土地利用類型下的平均徑流系數(shù)(見(jiàn)表1)。

表1 不同土地利用類型的地表徑流系數(shù)[34]

1.3 穩(wěn)健半方差估計(jì)識(shí)別空間離群值

在滿足本征假設(shè)下,傳統(tǒng)的半方差函數(shù)γ{h)的計(jì)算公式為[35]：

(2)

式中：γ*{h)為半方差函數(shù)γ{h)的估計(jì)值；h為樣本間距；N{h)為間隔距離為h的樣本觀測(cè)點(diǎn)對(duì)數(shù)；Z{xi+h)和Z{xi)分別為區(qū)域化變量Z{x)在位置xi+h和xi處的預(yù)測(cè)值。

傳統(tǒng)半方差函數(shù)容易受到樣本離群值的影響,為減小土壤全磷樣本離群值對(duì)其影響,本研究采用穩(wěn)健方法識(shí)別土壤全磷樣本數(shù)據(jù)可能的空間離群值。其中穩(wěn)健半方差主要有Cressie and Hawkins′s[36]估計(jì),Dowd′s[37]估計(jì)以及Genton′s[38]估計(jì)。而選擇最優(yōu)的穩(wěn)健方法是識(shí)別空間離群值的關(guān)鍵。我們通過(guò)普通克里格的交叉驗(yàn)證方法來(lái)判別傳統(tǒng)半方差和三種穩(wěn)健半方差預(yù)測(cè)的優(yōu)劣,其計(jì)算公式如下：

(3)

對(duì)于最優(yōu)的穩(wěn)健半方差普通克里格預(yù)測(cè),我們通過(guò)標(biāo)準(zhǔn)誤εs(xi)來(lái)識(shí)別可能的空間離群值,其計(jì)算公式如下：

(4)

1.4 序貫高斯模擬

序貫高斯模擬 (Sequential Gaussian Simulation,SGS) 是一種對(duì)連續(xù)變量進(jìn)行隨機(jī)模擬的地統(tǒng)計(jì)學(xué)方法。SGS根據(jù)現(xiàn)有樣本數(shù)據(jù)計(jì)算待模擬點(diǎn)值的條件概率分布,從該分布中隨機(jī)抽取一個(gè)值作為模擬實(shí)現(xiàn),每得到一個(gè)模擬值,就把它連同原始數(shù)據(jù)、模擬數(shù)據(jù)一起作為條件數(shù)據(jù)代入到下一點(diǎn)的模擬,該模型能產(chǎn)生多個(gè)可能的空間分布格局,有效避免普通克里格插值結(jié)果的平滑效應(yīng)[41]。為降低樣本離群值對(duì)土壤全磷SGS結(jié)果的影響,我們首先通過(guò)穩(wěn)健半方差函數(shù)來(lái)剔除空間離群值,再對(duì)剩下的樣本數(shù)據(jù)進(jìn)行SGS。對(duì)整個(gè)金壇區(qū)劃分成規(guī)則模擬格網(wǎng)(200 m × 200 m)后,SGS的基本步驟[17,42]為：

(1)確定一條訪問(wèn)每個(gè)網(wǎng)格節(jié)點(diǎn)一次的隨機(jī)路徑；

(2)對(duì)隨機(jī)路徑上的一個(gè)節(jié)點(diǎn)x0,利用樣本數(shù)據(jù)的方差函數(shù)和簡(jiǎn)單克里格方法求取該節(jié)點(diǎn)搜索半徑內(nèi)已知節(jié)點(diǎn)的高斯條件累積分布函數(shù)(ccdf)的參數(shù)(平均值和方差)；

(3)從求得的ccdf中隨機(jī)抽取一個(gè)值作為x0節(jié)點(diǎn)的模擬值,將x0節(jié)點(diǎn)處的模擬值加入到條件數(shù)據(jù)集中用于下一網(wǎng)格節(jié)點(diǎn)x1的模擬；

(4)沿著該隨機(jī)模擬路徑重復(fù)步驟(2)—(4)獲取整個(gè)路徑上所有網(wǎng)格節(jié)點(diǎn)的模擬結(jié)果,即可獲得一次隨機(jī)模擬實(shí)現(xiàn)。

為了產(chǎn)生500個(gè)等概率模擬實(shí)現(xiàn),上述步驟應(yīng)采用不同的隨機(jī)路徑來(lái)模擬?；谕寥廊椎?00個(gè)可能的空間模擬實(shí)現(xiàn),可計(jì)算對(duì)應(yīng)的500個(gè)可能的土壤磷素流失風(fēng)險(xiǎn)值的模擬實(shí)現(xiàn)。

1.5 模型評(píng)估標(biāo)準(zhǔn)

從259個(gè)土壤樣本點(diǎn)中隨機(jī)抽取20%樣本(52個(gè))作為模型精度驗(yàn)證數(shù)據(jù)集,用于對(duì)比兩種SGS模擬的E-type估計(jì)精度。分別計(jì)算驗(yàn)證點(diǎn)處土壤全磷的實(shí)測(cè)值與兩種隨機(jī)模擬方法產(chǎn)生的E-type估計(jì)之間的平均預(yù)測(cè)誤差(ME)、均方根誤差(RMSE)及Pearson相關(guān)系數(shù)(r),其中ME和RMSE計(jì)算方法如下：

(5)

(6)

式中：n為驗(yàn)證點(diǎn)數(shù)目；z{u) 和z*{u) 分別為驗(yàn)證點(diǎn)u處土壤全磷的實(shí)測(cè)值和預(yù)測(cè)值。較高的r和較低的ME、RMSE代表更高的預(yù)測(cè)精度。

RSGS較傳統(tǒng)SGS的E-type估計(jì)結(jié)果的相對(duì)提高指數(shù)(RI)為：

(7)

式中：RMSERSGS和RMSESGS分別為RSGS方法和傳統(tǒng)SGS方法的均方根誤差。若RI為正,則表示RSGS的E-type估計(jì)較傳統(tǒng)SGS的E-type估計(jì)具有更高的精度,反之,則RSGS的E-type估計(jì)預(yù)測(cè)精度低于SGS方法[43]。

1.6 土壤磷素流失風(fēng)險(xiǎn)評(píng)估模型及風(fēng)險(xiǎn)劃分

由于長(zhǎng)三角典型縣域金壇區(qū)具有農(nóng)田水利設(shè)施完善、多為平原等特點(diǎn),土壤侵蝕對(duì)土壤磷素流失的影響較弱,本研究對(duì)其不予考慮。我們基于傳統(tǒng)的磷指數(shù)法[6,14]構(gòu)建快速磷指數(shù)評(píng)估模型來(lái)計(jì)算研究區(qū)的土壤磷素流失風(fēng)險(xiǎn)指數(shù),其模型計(jì)算公式如下：

(8)

式中,RPI(u)表示土壤采樣點(diǎn)u處的快速磷指數(shù)值；P(u)表示土壤采樣點(diǎn)u處的土壤全磷濃度；F(u)為土壤采樣點(diǎn)u處的磷肥施用量,同一土地利用類型的磷肥施用量相同；D(u)為土壤采樣點(diǎn)u處到受納水體的最短距離；α(u)是空間位置u處的徑流系數(shù)；β1、β2、β3、β4分別為土壤全磷濃度、磷肥施用量、距受納水體距離、徑流系數(shù)的權(quán)重,參見(jiàn)Sharply、張淑榮[6- 7]等人的研究,本研究β1、β2、β3、β4分別取值為1、0.75、1、0.5。P(u)、F(u)、D(u)、α(u) 均歸一化處理,無(wú)量綱。

本研究采用自然間斷點(diǎn)分級(jí)法[44](Jenks)進(jìn)行土壤磷素流失風(fēng)險(xiǎn)的分級(jí),將金壇區(qū)分為為低風(fēng)險(xiǎn)區(qū)、較低風(fēng)險(xiǎn)區(qū)、較高風(fēng)險(xiǎn)區(qū)和高風(fēng)險(xiǎn)區(qū)。

1.7 不確定性評(píng)估

具體空間位置u處RPI值的超標(biāo)概率為在位置u的RPI(u)高于給定臨界閾值c的概率。該概率計(jì)算方法為：

(9)

式中,500為模擬的土壤磷素流失風(fēng)險(xiǎn)值的次數(shù)；RPI(u)是空間位置u處的RPI值；c是給定的臨界閾值,即以自然間斷點(diǎn)分級(jí)法劃定為高風(fēng)險(xiǎn)所對(duì)應(yīng)的值作為本研究的臨界閾值；n(u)是位置u上產(chǎn)生的模擬值高于給定閾值c的個(gè)數(shù)。

1.8 數(shù)據(jù)處理

常規(guī)統(tǒng)計(jì)分析在SPSS 24.0中完成。穩(wěn)健半方差擬合及離群值篩選采用R軟件。RSGS、近鄰分析及空間制圖通過(guò)ArcGIS 10.3完成。

2 結(jié)果與討論

2.1 土壤全磷的描述性統(tǒng)計(jì)

金壇區(qū)不同土地利用類型(即水田、旱地和其他用地)下土壤全磷的描述性統(tǒng)計(jì)量見(jiàn)表2。金壇區(qū)水田土壤全磷介于0.19—2.97 g/kg,其均值為0.74 g/kg；旱地土壤全磷介于0.27—0.98 g/kg,其均值為0.59 g/kg。金壇區(qū)表層(0—20 cm)土壤全磷平均值為0.67 g/kg,高于第二次土壤普查時(shí)耕層土壤全磷平均含量(0.50 g/kg),表明該區(qū)域在長(zhǎng)期磷肥的施用下土壤全磷在地表已產(chǎn)生一定程度的積累。變異系數(shù)的大小反映區(qū)域土壤屬性受地形、土地利用及不同耕作方式等因素影響的分布不均勻程度[45]。本研究中水田、旱地、其他用地的土壤全磷空間變異系數(shù)分別為52.24%、30.51%和43.24%,表明金壇區(qū)的土壤全磷變化受內(nèi)部因素和外部因素的共同作用,存在中等強(qiáng)度的變異性[46]。

表2 金壇區(qū)不同土地利用類型下土壤全磷描述性統(tǒng)計(jì)量

2.2 穩(wěn)健半方差模型參數(shù)及基于最優(yōu)穩(wěn)健半方差的空間離群值的篩選

土壤全磷經(jīng)對(duì)數(shù)轉(zhuǎn)換的四種半方差擬合參數(shù)見(jiàn)表3。與傳統(tǒng)的Matheron估計(jì)對(duì)比可知,三種穩(wěn)健估計(jì)的塊金值,基臺(tái)值均低于Matheron估計(jì),這可能是傳統(tǒng)實(shí)驗(yàn)半方差函數(shù)Matheron估計(jì)量γ*{h)是基于偏差的平方,對(duì)離群值較敏感,不具有穩(wěn)健性。所以當(dāng)結(jié)果中出現(xiàn)離群值時(shí),可能會(huì)高估實(shí)驗(yàn)半方差函數(shù)值。

四種半方差交叉檢驗(yàn)的統(tǒng)計(jì)量SSPE(xi)的中位數(shù)見(jiàn)表4。傳統(tǒng)的Matheron估計(jì)得到的SSPE(xi)的中位數(shù)明顯低于0.455。其中Dowd穩(wěn)健估計(jì),經(jīng)交叉檢驗(yàn)得到的SSPE(xi)的中位數(shù)為0.454,最接近目標(biāo)值0.455。因此我們基于Dowd估計(jì)進(jìn)行空間離群值的篩選。依據(jù)公式(4),我們共識(shí)別出9個(gè)土壤全磷的空間離群值,這些離群值不參與后續(xù)的SGS。

塊金系數(shù)C0/(C0+C)表示隨機(jī)因素部分引起的空間異質(zhì)性占系統(tǒng)總變異的比重,反映了區(qū)域土壤性質(zhì)空間相關(guān)性的強(qiáng)弱[47]。若C0/(C0+C)小于25%,則表現(xiàn)為強(qiáng)的空間相關(guān)性；在25%—75%之間,則存在中等強(qiáng)度的空間相關(guān)性；而大于75%,則表示空間相關(guān)性較弱?；谧顑?yōu)穩(wěn)健的Down估計(jì)所擬合模型的塊金系數(shù)為49%,屬中等強(qiáng)度空間自相關(guān)性,這是因?yàn)閰^(qū)域土壤養(yǎng)分分布是由結(jié)構(gòu)性因素(如氣候、母質(zhì)、地形、土壤類型等)和隨機(jī)性因素(如耕作措施、種植制度等)共同作用的結(jié)果。變程反映了變量空間自相關(guān)范圍的大小,當(dāng)觀測(cè)值之間的距離小于該值時(shí),說(shuō)明它們之間存在一定的相關(guān)關(guān)系,若大于該值,則說(shuō)明它們之間是相互獨(dú)立的。本研究中所擬合模型的變程為8.5 km,一定程度上可反應(yīng)在變程范圍內(nèi),土壤全磷表現(xiàn)出中等強(qiáng)度的空間自相關(guān)關(guān)系。本研究中土壤全磷的這種空間結(jié)構(gòu)特征可能與區(qū)域種植制度有關(guān),在相似氣候與地形的長(zhǎng)三角縣域尺度下,研究區(qū)內(nèi)相對(duì)一致的種植制度使得農(nóng)田管理和養(yǎng)分的投入長(zhǎng)期保持穩(wěn)定,因此土壤全磷的空間分布表現(xiàn)出一定的規(guī)律性。

表3 土壤全磷基于四種半方差估計(jì)的的變異函數(shù)參數(shù)

表4 基于四種半方差模型的交叉檢驗(yàn)的統(tǒng)計(jì)量SSPE(xi)的中位數(shù)

2.3 模擬精度對(duì)比

如表5所示,52個(gè)驗(yàn)證點(diǎn)處土壤全磷采用SGS的E-type估計(jì)的ME、RMSE以及相關(guān)系數(shù)r分別為0.06、0.28和0.68,而采用RSGS的E-type估計(jì)的ME、RMSE分別為0.03、0.21和0.80,RSGS的E-type估計(jì)較SGS的E-type估計(jì)的相對(duì)提高指數(shù)為25%。說(shuō)明RSGS較傳統(tǒng)的SGS有更高的模擬精度,因此后續(xù)的土壤全磷的空間格局、土壤磷素流失風(fēng)險(xiǎn)以及不確定性評(píng)估均采用RSGS的模擬結(jié)果。

表5 SGS的E-type和RSGS的E-type估計(jì)模擬土壤全磷的預(yù)測(cè)精度對(duì)比

2.4 土壤全磷的空間分布

采用RSGS對(duì)金壇區(qū)的土壤全磷含量進(jìn)行空間隨機(jī)模擬,其產(chǎn)生的第100(a)、300(b)、500(c)次模擬實(shí)現(xiàn)以及E-type估計(jì)圖如圖2所示。每個(gè)模擬實(shí)現(xiàn)代表土壤全磷一個(gè)可能的空間分布模式,這些模擬實(shí)現(xiàn)所對(duì)應(yīng)的E-type估計(jì)圖代表了最佳估計(jì),且該最佳估計(jì)具有一定的平滑效應(yīng)。由圖2可知,E-type預(yù)測(cè)圖在表現(xiàn)區(qū)域土壤全磷含量空間變化的總體趨勢(shì)上較為明顯,土壤全磷在研究區(qū)的北部、中部以及洮湖周邊含量較高,而西部的土壤全磷含量明顯低于全區(qū)的平均值。雖然各個(gè)模擬實(shí)現(xiàn)在總體上顯示出相似的空間分布模式,但在局部空間變化上卻是有差別的,這種局部的變化能夠反映出該區(qū)域土壤全磷的可能空間分布水平。

本研究中土壤全磷空間分布差異可能是由于區(qū)域土壤特性造成的。金壇區(qū)的中東部主要為水田,其土壤主要發(fā)育于石灰?guī)r,土壤質(zhì)地較粘,保肥能力較強(qiáng),因此其土壤磷素含量普遍較高。西部林地的土壤主要發(fā)育于石英砂巖、頁(yè)巖、角礫巖、玄武巖等母巖,多數(shù)為中壤土,且西部林地受到的人為干擾較少,其土壤全磷主要來(lái)源于母質(zhì),因此其土壤磷素含量較低。而對(duì)于西部的旱地區(qū),經(jīng)調(diào)查發(fā)現(xiàn)其肥料投入較水田多,但該地區(qū)土壤淋洗作用較強(qiáng),所以其土壤全磷的含量較中東部水田低。土壤全磷含量可表征區(qū)域土壤磷素的總儲(chǔ)量,較高的土壤全磷含量說(shuō)明影響土壤磷素流失的本底值較高,但難以完全反應(yīng)本研究區(qū)域土壤磷素流失的風(fēng)險(xiǎn)程度。磷肥的施用量、地表徑流以及距受納水體的距離也是重要的影響因子。

圖2 土壤全磷穩(wěn)健序貫高斯模擬(RSGS)的第100(a)、300(b)、500(c)次以及E-type(d)圖Fig.2 The 100th realization (a), 300th realization (b), 500th realization (c) and E-type estimate (d) of soil TP generated by RSGS

2.5 土壤磷素流失風(fēng)險(xiǎn)指數(shù)空間分布與綜合風(fēng)險(xiǎn)評(píng)估

金壇區(qū)土壤磷素流失風(fēng)險(xiǎn)指數(shù)的空間分布模式如圖3所示,土壤磷素流失風(fēng)險(xiǎn)較高的區(qū)域主要集中分布在洮湖周邊、金壇主城區(qū)的東部及北部,而金壇區(qū)西部的RPI值較低。圖3d與土壤全磷R(shí)SGS模擬生成的E-type圖(圖2d)對(duì)比可知,結(jié)合環(huán)境因子評(píng)估的土壤磷素流失高風(fēng)險(xiǎn)區(qū)與土壤全磷的高值區(qū)在空間分布上具有一定的相似性,在洮湖的周邊,土壤全磷的含量較高,且離受納水體距離較近,故RPI值較大,該區(qū)域土壤磷素必然存在較高的流失風(fēng)險(xiǎn)。同時(shí),結(jié)合環(huán)境因子的評(píng)估結(jié)果與土壤全磷的分布又存在明顯的空間差異,研究區(qū)中西部的旱地區(qū),土壤全磷的含量較低,但該區(qū)有較大河流,同時(shí)地表徑流系數(shù)較水田高,故該區(qū)域表現(xiàn)出較高的土壤磷素流失風(fēng)險(xiǎn)。自然間斷點(diǎn)分級(jí)法將RPI值劃分為四個(gè)等級(jí),區(qū)間為0.47—0.76、0.77—0.92、0.93—1.05、1.06—1.46分別對(duì)應(yīng)低、較低、較高、高四個(gè)風(fēng)險(xiǎn)等級(jí)(表6)。土壤磷素流失風(fēng)險(xiǎn)性“高”的區(qū)域即為金壇區(qū)土壤磷素流失的關(guān)鍵源區(qū)。由表6可知,金壇區(qū)土壤磷素流失較高風(fēng)險(xiǎn)區(qū)的面積占金壇區(qū)總面積的40.94%,而高風(fēng)險(xiǎn)區(qū)的面積占24.94%,RPI值在0.93—1.46(即較高和高風(fēng)險(xiǎn)區(qū))的面積占金壇區(qū)總面積65.88%,而低及較低風(fēng)險(xiǎn)區(qū)的面積占比僅為34.12%,由此可見(jiàn)整個(gè)研究區(qū)內(nèi)土壤磷素潛在流失風(fēng)險(xiǎn)較大。

金壇區(qū)土壤磷素流失風(fēng)險(xiǎn)等級(jí)見(jiàn)圖4。研究區(qū)內(nèi)的土壤磷素流失高風(fēng)險(xiǎn)區(qū)主要沿著河流呈現(xiàn)條帶狀及斑塊狀的分布特征,主要分布在研究區(qū)的北部及中東部。西部林地的流失風(fēng)險(xiǎn)較低,但西部的旱地由于磷肥投入量、地表徑流潛力較大,因此該區(qū)域土壤磷素也存在較高的潛在流失風(fēng)險(xiǎn),但主要分布在河流周邊,初步分析金壇區(qū)的土壤磷素流失風(fēng)險(xiǎn)分布格局,主要是由土壤全磷、磷肥施用量、距離因子以及地表徑流潛力共同作用的結(jié)果。

本研究提供的方法僅僅模擬土壤磷素流失的風(fēng)險(xiǎn)指數(shù),而土壤磷素的實(shí)際流失往往是一個(gè)復(fù)雜的自然生態(tài)過(guò)程?，F(xiàn)有的研究表明,林地的土壤磷素流失量較低,耕地土壤磷素流失量較高[48]。這主要是因?yàn)榱值赝寥懒姿氐耐庠摧斎脒h(yuǎn)遠(yuǎn)小于耕地,且林地較好的植被覆蓋能夠有效截留徑流中的磷素,在本研究的模擬結(jié)果中,金壇區(qū)西部林地的土壤磷素流失風(fēng)險(xiǎn)指數(shù)遠(yuǎn)遠(yuǎn)低于中東部的耕地。孫金華等[49]利用AnnAGNPS模型模擬太湖流域雪堰鎮(zhèn)平原水網(wǎng)圩區(qū)的土壤磷素流失發(fā)現(xiàn),土壤磷素的流失主要集中在河道兩側(cè),這與本研究的模擬結(jié)果一致。而降雨作為土壤磷素的遷移動(dòng)力,也顯著影響著土壤磷素的流失。而本研究針對(duì)降雨對(duì)土壤磷素流失的影響存在一定的局限性,即沒(méi)有考慮由降雨的季節(jié)性差異而導(dǎo)致的土壤磷素流失風(fēng)險(xiǎn)在時(shí)間尺度上的差異性。傳統(tǒng)的PI法用于評(píng)估區(qū)域土壤磷素流失風(fēng)險(xiǎn)已經(jīng)不斷得到應(yīng)用和改進(jìn),本研究從土壤全磷空間模擬及風(fēng)險(xiǎn)的不確定角度出發(fā),鑒于區(qū)域土壤數(shù)據(jù)集一般含有一定數(shù)量的離群值,這些往往會(huì)對(duì)傳統(tǒng)地統(tǒng)計(jì)學(xué)模型預(yù)測(cè)精度產(chǎn)生一定影響,故本研究采用最優(yōu)穩(wěn)健不確定性模擬方法,降低空間離群值對(duì)模擬結(jié)果的影響,且能產(chǎn)生多個(gè)可能的空間模擬實(shí)現(xiàn),這樣的結(jié)果輸入到RPI模型中,便可快速評(píng)估具有類似生態(tài)特征的土壤磷素流失風(fēng)險(xiǎn)及其不確定性。

圖3 土壤磷素流失風(fēng)險(xiǎn)快速指數(shù)值(RPI)的第100(a)、300(b)、500(c)次模擬實(shí)現(xiàn)以及E-type圖(d) Fig.3 Map of the 100th realization (a), 300th realization (b), 500th realization (c) and E-type estimate (d) of RPI

表6 金壇區(qū)快速磷指數(shù)值(RPI)的等級(jí)劃分與面積占比

Table 6 Classification and area ratio of RPI in Jintan District

土壤磷素流失風(fēng)險(xiǎn)指數(shù)值RPI0.47—0.760.77—0.920.93—1.051.06—1.46風(fēng)險(xiǎn)等級(jí)Risk grade 低較低較高高面積占比Area ratio8.04%26.08%40.94%24.94%

2.6 土壤磷素流失風(fēng)險(xiǎn)的不確定性評(píng)估

金壇區(qū)土壤磷素流失風(fēng)險(xiǎn)指數(shù)超過(guò)給定閾值的概率分布如圖5所示。RPI相對(duì)較高(P>0.5)的區(qū)域主要集中在中、東部以及西部沿河流周邊區(qū)域。金壇區(qū)土壤磷素流失風(fēng)險(xiǎn)指數(shù)的不確定性評(píng)價(jià)結(jié)果顯示,當(dāng)閾值概率劃定為0.50、0.75、0.85、0.95時(shí),超閾值概率的區(qū)域面積分別占金壇區(qū)總面積的16.71%、5.74%、2.84%、1.04%(圖6)。可以看出隨著閾值概率的增加,所劃定的區(qū)域面積不斷減小。當(dāng)閾值概率為0.85時(shí)(圖6c),研究區(qū)中主城區(qū)的西北部、洮湖的周邊以及西部旱地部分地區(qū)土壤磷素流失風(fēng)險(xiǎn)等級(jí)高具有很大的可能性。對(duì)于那些概率值接近50%的區(qū)域(圖6a),其評(píng)價(jià)結(jié)果就有很大的不確定性。所以這樣的概率分布能夠給管理部門的決策提供更多空間信息,而不是直接按某個(gè)閾值將整個(gè)地區(qū)機(jī)械地劃分為不同等級(jí)的污染區(qū)域[50]。閾值概率設(shè)定為0.85時(shí),圖6c中的劃定的區(qū)域就要重點(diǎn)實(shí)施管控措施。

圖4 金壇區(qū)土壤磷素流失風(fēng)險(xiǎn)等級(jí)圖Fig.4 Level of soil phosphorus loss risk in Jintan District

圖5 金壇區(qū)土壤磷素流失風(fēng)險(xiǎn)指數(shù)(RPI)超給定閾值(1.06)的概率分布圖Fig.5 Probability map of RPI exceeding the threshold (RPI≥1.06) in Jintan District

圖6 分別按0.50、0.75、0.85、0.95的臨界概率劃定的土壤磷素流失高風(fēng)險(xiǎn)區(qū)域Fig.6 The highest risk areas of soil phosphorus loss were defined according to the critical probability of 0.50(a), 0.75(b), 0.85(c) and 0.95(d), respectively

3 結(jié)論

本研究結(jié)合長(zhǎng)三角水網(wǎng)發(fā)達(dá)、多為平原、縣域景觀異質(zhì)性小、農(nóng)田水利工程完善等縣域特征,在傳統(tǒng)分析土壤磷素空間分布格局的基礎(chǔ)上,引入磷肥施用量、距受納水體距離、地表徑流潛力構(gòu)建了針對(duì)長(zhǎng)三角典型縣域土壤磷素流失風(fēng)險(xiǎn)的RPI評(píng)估模型；同時(shí),我們將金壇區(qū)土壤全磷可能的RSGS實(shí)現(xiàn)輸入該模型,進(jìn)而評(píng)估了流失風(fēng)險(xiǎn)指數(shù)可能的空間分布格局；最后我們依據(jù)不同的超標(biāo)概率閾值劃定了金壇土壤磷素流失風(fēng)險(xiǎn)的關(guān)鍵源區(qū)即高風(fēng)險(xiǎn)區(qū)。結(jié)果顯示土壤磷素流失的高風(fēng)險(xiǎn)區(qū)主要沿著河流呈現(xiàn)條帶狀及斑塊狀分布,集中分布在洮湖周邊、金壇主城區(qū)的東部及西部,較高及高風(fēng)險(xiǎn)區(qū)(快速磷指數(shù)值大于0.93)的面積占金壇區(qū)面積的65.88%。概率閾值分別設(shè)定為0.50、0.75、0.85、0.95時(shí),其超標(biāo)面積占金壇區(qū)總面積分別達(dá)到16.71%、5.74%、2.84%、1.04%。本研究提供的土壤磷素流失風(fēng)險(xiǎn)信息,有助于決策者從“不確定性”角度為區(qū)域土壤磷素的調(diào)控,劃定更符合實(shí)際限制條件的管控區(qū)域。