陳 劍,瞿明凱,王 燕,萬夢雪,黃 標,趙永存
1 中國科學院南京土壤研究所土壤環境與污染修復重點實驗室, 南京 210008 2 中國科學院大學, 北京 100049
土壤磷素是作物生長的主要營養元素,同時也是水體富營養化的關鍵限制因子[1- 5]。由于農田土壤磷素具有較強的空間變異性,而且影響磷素流失的環境因子較多,因此對其流失進行控制和治理難度較大,需耗費大量的人力、物力和財力[6]。許多研究結果顯示,在土壤磷素向水體輸出的過程中,少數景觀單元輸出的磷素往往占了整個流域污染負荷的大部分,這些景觀單元則被稱為土壤磷素流失的關鍵源區[7](Critical Source Areas, CSAs)。而對土壤磷素污染展開風險評估和CSAs的劃定或識別,是高效治理污染的關鍵[8]。在眾多的評估方法中,磷指數法(Phosphorus Index,PI)因其簡便實用,且評估結果具有較好可視性的特點而得到了廣泛的運用。但傳統的PI在評估過程中,對各影響因子提前分等定級,具有一定的人為主觀性。而且在縣域尺度下的多數評估未考慮研究區域特征,評估要求的數據量大,過程復雜。
對于現有的PI法評估,往往未考慮土壤磷素具有較強的空間變異性,且忽視了預測結果的空間不確定性。區域土壤磷素濃度水平的高低往往直接影響其流失風險程度[9- 14]。土壤全磷含量通常用于衡量區域土壤磷素的本底水平[15- 16]。地統計學模型通常被用于區域土壤屬性的空間預測。然而,基于有限樣本點預測得到的土壤全磷空間分布必然存在一定的空間不確定性[17]。傳統的地統計學插值模型只能得到唯一的預測結果,故不能反映這一不確定性特征。而且,區域土壤數據集可能含有一定數量的離群值,這些往往會對傳統地統計學模型預測精度產生一定影響。構建穩健不確定性模擬方法以去除離群值對模擬結果的影響,且能產生多個可能的空間模擬實現。這些模擬實現可以輸入土壤磷素流失風險評估模型,進一步評估其流失風險及相應的不確定性。除了區域土壤全磷含量,土壤磷素的流失風險還與其他環境因子密切相關,如距受納水體距離、地表徑流潛力、土壤侵蝕潛力、磷肥施用量等[18-28]。
磷肥的施用是土壤磷素的一個重要外源輸入,縣域尺度下磷肥的施用量主要和土地利用類型相關;雖然磷素在不同性質的下墊面遷移過程中衰減機制尚具有較大的差異[29],但磷素在遷移過程中不斷被截留和稀釋,距離受納水體較遠的磷素源區比距離較近的源區貢獻量小,因此距離因子可在一定程度上表征土壤磷素在遷移過程中隨著距離衰減的規律;縣域尺度下地表徑流量主要由徑流系數決定,徑流系數與土地利用類型密切相關[30]。由于長三角農田水利設施完善,排水溝渠及擋土墻可有效防治土壤侵蝕,故本研究針對土壤侵蝕不予考慮。
本研究以長三角典型縣域金壇為研究區,具體目的在于:(1)采用經穩健處理的序貫高斯模型(RSGS)模擬金壇區土壤全磷可能的空間分布格局,并采用獨立驗證點,與傳統序貫高斯模擬模型(SGS)的模擬精度進行對比,以驗證RSGS在提高磷素空間模擬精度方面的有效性;(2)結合相關環境因子(即磷肥施用量、距受納水體距離、地表徑流潛力)構建快速磷指數(RPI)評估模型;(3)評估長三角典型縣域土壤磷素的流失風險及其空間不確定性,并識別土壤磷素流失的關鍵源區。本研究以期從“不確定性”角度為區域土壤磷素的調控,劃定更符合實際限制條件的管控區域。
本研究選取長江三角洲典型農業區金壇區(31°33′42″—31°53′22″ N,119°17′45″—119°44′59″ E)作為研究區域(圖1)。其地處江蘇省南部,茅山東麓,是南京、上海和杭州三角地帶的中心。全區總面積約975平方公里,平均海拔13m。金壇區內水網密布,水域面積達194.22km2,境內有大小河流216條,總長512km。東南部的洮湖(長蕩湖)是江蘇省十大淡水湖之一。金壇區屬北亞熱帶季風區,雨熱同季,四季分明;雨量充沛,年平均降雨量為1063.5 mm;日照充足,年平均氣溫為15.3 ℃,無霜期長達228d,年平均濕度為78%。區內農用地類型單一,主要為水田和旱地,農田水利工程完善,排水溝渠及擋土墻可有效防治土壤侵蝕,但密集的溝渠亦是地表徑流磷素流失的快速通道,有利于土壤磷素的遷移。金壇區部分區域化肥施用量大,全年每畝農田化肥施用量達112.5 kg,其中水稻全年每畝施肥量約是全國平均水平的3.7倍[31]。該區主要糧食作物為水稻和小麥,是太湖地區重要的糧食生產基地之一。
(1)土壤全磷
金壇區內共采集土壤樣點259個(圖1),采樣時間為2016年5月底(即小麥收割后至水稻栽種前)。采集每個樣點周圍100 m2范圍內4—5處表層(0—20 cm)土樣,均勻混合后縮分至1—2 kg用于實驗室分析,每個樣本點采用GPS定位并記錄其詳細信息。所有樣本在室溫(20—22 ℃)風干后去除雜物及碎石,磨細,過100目篩,收集50 g土樣用于土壤全磷的測定。土壤全磷采用HClO4-H2SO4消煮,鉬銻抗比色法測定。具體分析方法參見魯如坤[32]。

圖1 金壇區土地利用類型及土壤采樣點分布圖Fig.1 Map of study area, land-use types, and soil sampling sites in Jintan District
(2)環境因子
用于評估長三角典型縣域金壇區土壤磷素流失風險的三個重要環境因子為磷肥施用量因子、采樣點距受納水體距離因子及地表徑流潛力因子。磷肥施用量以2014年、2015年和2016年金壇區的統計年鑒為依據,以三年磷肥施用量平均值計算,其中旱地磷肥施用量以水田的2.5倍計[33],其他土地利用類型均考慮為無肥區,單位面積的磷肥施用量因子由金壇區磷肥施用總量和水田、旱地的面積比來獲取;采樣點距受納水體的距離則通過金壇區內水系矢量數據(1∶25萬)采用近鄰分析方法得到;地表徑流潛力可采用徑流深R(cm)來度量,其計算公式如下:
R=p·α
(1)
式中,p表示年降雨量(cm);α表示年徑流系數。
本研究中,考慮到縣域尺度內降雨量區域差異相對較小,地表徑流主要受年徑流系數α的影響。因此,本研究根據金壇區的特點和相關研究來確定該區域不同土地利用類型下的平均徑流系數(見表1)。

表1 不同土地利用類型的地表徑流系數[34]
在滿足本征假設下,傳統的半方差函數γ{h)的計算公式為[35]:
(2)
式中:γ*{h)為半方差函數γ{h)的估計值;h為樣本間距;N{h)為間隔距離為h的樣本觀測點對數;Z{xi+h)和Z{xi)分別為區域化變量Z{x)在位置xi+h和xi處的預測值。
傳統半方差函數容易受到樣本離群值的影響,為減小土壤全磷樣本離群值對其影響,本研究采用穩健方法識別土壤全磷樣本數據可能的空間離群值。其中穩健半方差主要有Cressie and Hawkins′s[36]估計,Dowd′s[37]估計以及Genton′s[38]估計。而選擇最優的穩健方法是識別空間離群值的關鍵。我們通過普通克里格的交叉驗證方法來判別傳統半方差和三種穩健半方差預測的優劣,其計算公式如下:
(3)

對于最優的穩健半方差普通克里格預測,我們通過標準誤εs(xi)來識別可能的空間離群值,其計算公式如下:
(4)

序貫高斯模擬 (Sequential Gaussian Simulation,SGS) 是一種對連續變量進行隨機模擬的地統計學方法。SGS根據現有樣本數據計算待模擬點值的條件概率分布,從該分布中隨機抽取一個值作為模擬實現,每得到一個模擬值,就把它連同原始數據、模擬數據一起作為條件數據代入到下一點的模擬,該模型能產生多個可能的空間分布格局,有效避免普通克里格插值結果的平滑效應[41]。為降低樣本離群值對土壤全磷SGS結果的影響,我們首先通過穩健半方差函數來剔除空間離群值,再對剩下的樣本數據進行SGS。對整個金壇區劃分成規則模擬格網(200 m × 200 m)后,SGS的基本步驟[17,42]為:
(1)確定一條訪問每個網格節點一次的隨機路徑;
(2)對隨機路徑上的一個節點x0,利用樣本數據的方差函數和簡單克里格方法求取該節點搜索半徑內已知節點的高斯條件累積分布函數(ccdf)的參數(平均值和方差);
(3)從求得的ccdf中隨機抽取一個值作為x0節點的模擬值,將x0節點處的模擬值加入到條件數據集中用于下一網格節點x1的模擬;
(4)沿著該隨機模擬路徑重復步驟(2)—(4)獲取整個路徑上所有網格節點的模擬結果,即可獲得一次隨機模擬實現。
為了產生500個等概率模擬實現,上述步驟應采用不同的隨機路徑來模擬。基于土壤全磷的500個可能的空間模擬實現,可計算對應的500個可能的土壤磷素流失風險值的模擬實現。
從259個土壤樣本點中隨機抽取20%樣本(52個)作為模型精度驗證數據集,用于對比兩種SGS模擬的E-type估計精度。分別計算驗證點處土壤全磷的實測值與兩種隨機模擬方法產生的E-type估計之間的平均預測誤差(ME)、均方根誤差(RMSE)及Pearson相關系數(r),其中ME和RMSE計算方法如下:
(5)
(6)
式中:n為驗證點數目;z{u) 和z*{u) 分別為驗證點u處土壤全磷的實測值和預測值。較高的r和較低的ME、RMSE代表更高的預測精度。
RSGS較傳統SGS的E-type估計結果的相對提高指數(RI)為:
(7)
式中:RMSERSGS和RMSESGS分別為RSGS方法和傳統SGS方法的均方根誤差。若RI為正,則表示RSGS的E-type估計較傳統SGS的E-type估計具有更高的精度,反之,則RSGS的E-type估計預測精度低于SGS方法[43]。
由于長三角典型縣域金壇區具有農田水利設施完善、多為平原等特點,土壤侵蝕對土壤磷素流失的影響較弱,本研究對其不予考慮。我們基于傳統的磷指數法[6,14]構建快速磷指數評估模型來計算研究區的土壤磷素流失風險指數,其模型計算公式如下:
(8)
式中,RPI(u)表示土壤采樣點u處的快速磷指數值;P(u)表示土壤采樣點u處的土壤全磷濃度;F(u)為土壤采樣點u處的磷肥施用量,同一土地利用類型的磷肥施用量相同;D(u)為土壤采樣點u處到受納水體的最短距離;α(u)是空間位置u處的徑流系數;β1、β2、β3、β4分別為土壤全磷濃度、磷肥施用量、距受納水體距離、徑流系數的權重,參見Sharply、張淑榮[6- 7]等人的研究,本研究β1、β2、β3、β4分別取值為1、0.75、1、0.5。P(u)、F(u)、D(u)、α(u) 均歸一化處理,無量綱。
本研究采用自然間斷點分級法[44](Jenks)進行土壤磷素流失風險的分級,將金壇區分為為低風險區、較低風險區、較高風險區和高風險區。
具體空間位置u處RPI值的超標概率為在位置u的RPI(u)高于給定臨界閾值c的概率。該概率計算方法為:
(9)
式中,500為模擬的土壤磷素流失風險值的次數;RPI(u)是空間位置u處的RPI值;c是給定的臨界閾值,即以自然間斷點分級法劃定為高風險所對應的值作為本研究的臨界閾值;n(u)是位置u上產生的模擬值高于給定閾值c的個數。
常規統計分析在SPSS 24.0中完成。穩健半方差擬合及離群值篩選采用R軟件。RSGS、近鄰分析及空間制圖通過ArcGIS 10.3完成。
金壇區不同土地利用類型(即水田、旱地和其他用地)下土壤全磷的描述性統計量見表2。金壇區水田土壤全磷介于0.19—2.97 g/kg,其均值為0.74 g/kg;旱地土壤全磷介于0.27—0.98 g/kg,其均值為0.59 g/kg。金壇區表層(0—20 cm)土壤全磷平均值為0.67 g/kg,高于第二次土壤普查時耕層土壤全磷平均含量(0.50 g/kg),表明該區域在長期磷肥的施用下土壤全磷在地表已產生一定程度的積累。變異系數的大小反映區域土壤屬性受地形、土地利用及不同耕作方式等因素影響的分布不均勻程度[45]。本研究中水田、旱地、其他用地的土壤全磷空間變異系數分別為52.24%、30.51%和43.24%,表明金壇區的土壤全磷變化受內部因素和外部因素的共同作用,存在中等強度的變異性[46]。

表2 金壇區不同土地利用類型下土壤全磷描述性統計量
土壤全磷經對數轉換的四種半方差擬合參數見表3。與傳統的Matheron估計對比可知,三種穩健估計的塊金值,基臺值均低于Matheron估計,這可能是傳統實驗半方差函數Matheron估計量γ*{h)是基于偏差的平方,對離群值較敏感,不具有穩健性。所以當結果中出現離群值時,可能會高估實驗半方差函數值。
四種半方差交叉檢驗的統計量SSPE(xi)的中位數見表4。傳統的Matheron估計得到的SSPE(xi)的中位數明顯低于0.455。其中Dowd穩健估計,經交叉檢驗得到的SSPE(xi)的中位數為0.454,最接近目標值0.455。因此我們基于Dowd估計進行空間離群值的篩選。依據公式(4),我們共識別出9個土壤全磷的空間離群值,這些離群值不參與后續的SGS。
塊金系數C0/(C0+C)表示隨機因素部分引起的空間異質性占系統總變異的比重,反映了區域土壤性質空間相關性的強弱[47]。若C0/(C0+C)小于25%,則表現為強的空間相關性;在25%—75%之間,則存在中等強度的空間相關性;而大于75%,則表示空間相關性較弱。基于最優穩健的Down估計所擬合模型的塊金系數為49%,屬中等強度空間自相關性,這是因為區域土壤養分分布是由結構性因素(如氣候、母質、地形、土壤類型等)和隨機性因素(如耕作措施、種植制度等)共同作用的結果。變程反映了變量空間自相關范圍的大小,當觀測值之間的距離小于該值時,說明它們之間存在一定的相關關系,若大于該值,則說明它們之間是相互獨立的。本研究中所擬合模型的變程為8.5 km,一定程度上可反應在變程范圍內,土壤全磷表現出中等強度的空間自相關關系。本研究中土壤全磷的這種空間結構特征可能與區域種植制度有關,在相似氣候與地形的長三角縣域尺度下,研究區內相對一致的種植制度使得農田管理和養分的投入長期保持穩定,因此土壤全磷的空間分布表現出一定的規律性。

表3 土壤全磷基于四種半方差估計的的變異函數參數

表4 基于四種半方差模型的交叉檢驗的統計量SSPE(xi)的中位數
如表5所示,52個驗證點處土壤全磷采用SGS的E-type估計的ME、RMSE以及相關系數r分別為0.06、0.28和0.68,而采用RSGS的E-type估計的ME、RMSE分別為0.03、0.21和0.80,RSGS的E-type估計較SGS的E-type估計的相對提高指數為25%。說明RSGS較傳統的SGS有更高的模擬精度,因此后續的土壤全磷的空間格局、土壤磷素流失風險以及不確定性評估均采用RSGS的模擬結果。

表5 SGS的E-type和RSGS的E-type估計模擬土壤全磷的預測精度對比
采用RSGS對金壇區的土壤全磷含量進行空間隨機模擬,其產生的第100(a)、300(b)、500(c)次模擬實現以及E-type估計圖如圖2所示。每個模擬實現代表土壤全磷一個可能的空間分布模式,這些模擬實現所對應的E-type估計圖代表了最佳估計,且該最佳估計具有一定的平滑效應。由圖2可知,E-type預測圖在表現區域土壤全磷含量空間變化的總體趨勢上較為明顯,土壤全磷在研究區的北部、中部以及洮湖周邊含量較高,而西部的土壤全磷含量明顯低于全區的平均值。雖然各個模擬實現在總體上顯示出相似的空間分布模式,但在局部空間變化上卻是有差別的,這種局部的變化能夠反映出該區域土壤全磷的可能空間分布水平。
本研究中土壤全磷空間分布差異可能是由于區域土壤特性造成的。金壇區的中東部主要為水田,其土壤主要發育于石灰巖,土壤質地較粘,保肥能力較強,因此其土壤磷素含量普遍較高。西部林地的土壤主要發育于石英砂巖、頁巖、角礫巖、玄武巖等母巖,多數為中壤土,且西部林地受到的人為干擾較少,其土壤全磷主要來源于母質,因此其土壤磷素含量較低。而對于西部的旱地區,經調查發現其肥料投入較水田多,但該地區土壤淋洗作用較強,所以其土壤全磷的含量較中東部水田低。土壤全磷含量可表征區域土壤磷素的總儲量,較高的土壤全磷含量說明影響土壤磷素流失的本底值較高,但難以完全反應本研究區域土壤磷素流失的風險程度。磷肥的施用量、地表徑流以及距受納水體的距離也是重要的影響因子。

圖2 土壤全磷穩健序貫高斯模擬(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
金壇區土壤磷素流失風險指數的空間分布模式如圖3所示,土壤磷素流失風險較高的區域主要集中分布在洮湖周邊、金壇主城區的東部及北部,而金壇區西部的RPI值較低。圖3d與土壤全磷RSGS模擬生成的E-type圖(圖2d)對比可知,結合環境因子評估的土壤磷素流失高風險區與土壤全磷的高值區在空間分布上具有一定的相似性,在洮湖的周邊,土壤全磷的含量較高,且離受納水體距離較近,故RPI值較大,該區域土壤磷素必然存在較高的流失風險。同時,結合環境因子的評估結果與土壤全磷的分布又存在明顯的空間差異,研究區中西部的旱地區,土壤全磷的含量較低,但該區有較大河流,同時地表徑流系數較水田高,故該區域表現出較高的土壤磷素流失風險。自然間斷點分級法將RPI值劃分為四個等級,區間為0.47—0.76、0.77—0.92、0.93—1.05、1.06—1.46分別對應低、較低、較高、高四個風險等級(表6)。土壤磷素流失風險性“高”的區域即為金壇區土壤磷素流失的關鍵源區。由表6可知,金壇區土壤磷素流失較高風險區的面積占金壇區總面積的40.94%,而高風險區的面積占24.94%,RPI值在0.93—1.46(即較高和高風險區)的面積占金壇區總面積65.88%,而低及較低風險區的面積占比僅為34.12%,由此可見整個研究區內土壤磷素潛在流失風險較大。
金壇區土壤磷素流失風險等級見圖4。研究區內的土壤磷素流失高風險區主要沿著河流呈現條帶狀及斑塊狀的分布特征,主要分布在研究區的北部及中東部。西部林地的流失風險較低,但西部的旱地由于磷肥投入量、地表徑流潛力較大,因此該區域土壤磷素也存在較高的潛在流失風險,但主要分布在河流周邊,初步分析金壇區的土壤磷素流失風險分布格局,主要是由土壤全磷、磷肥施用量、距離因子以及地表徑流潛力共同作用的結果。
本研究提供的方法僅僅模擬土壤磷素流失的風險指數,而土壤磷素的實際流失往往是一個復雜的自然生態過程。現有的研究表明,林地的土壤磷素流失量較低,耕地土壤磷素流失量較高[48]。這主要是因為林地土壤磷素的外源輸入遠遠小于耕地,且林地較好的植被覆蓋能夠有效截留徑流中的磷素,在本研究的模擬結果中,金壇區西部林地的土壤磷素流失風險指數遠遠低于中東部的耕地。孫金華等[49]利用AnnAGNPS模型模擬太湖流域雪堰鎮平原水網圩區的土壤磷素流失發現,土壤磷素的流失主要集中在河道兩側,這與本研究的模擬結果一致。而降雨作為土壤磷素的遷移動力,也顯著影響著土壤磷素的流失。而本研究針對降雨對土壤磷素流失的影響存在一定的局限性,即沒有考慮由降雨的季節性差異而導致的土壤磷素流失風險在時間尺度上的差異性。傳統的PI法用于評估區域土壤磷素流失風險已經不斷得到應用和改進,本研究從土壤全磷空間模擬及風險的不確定角度出發,鑒于區域土壤數據集一般含有一定數量的離群值,這些往往會對傳統地統計學模型預測精度產生一定影響,故本研究采用最優穩健不確定性模擬方法,降低空間離群值對模擬結果的影響,且能產生多個可能的空間模擬實現,這樣的結果輸入到RPI模型中,便可快速評估具有類似生態特征的土壤磷素流失風險及其不確定性。

圖3 土壤磷素流失風險快速指數值(RPI)的第100(a)、300(b)、500(c)次模擬實現以及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 金壇區快速磷指數值(RPI)的等級劃分與面積占比
Table 6 Classification and area ratio of RPI in Jintan District

土壤磷素流失風險指數值RPI0.47—0.760.77—0.920.93—1.051.06—1.46風險等級Risk grade 低較低較高高面積占比Area ratio8.04%26.08%40.94%24.94%
金壇區土壤磷素流失風險指數超過給定閾值的概率分布如圖5所示。RPI相對較高(P>0.5)的區域主要集中在中、東部以及西部沿河流周邊區域。金壇區土壤磷素流失風險指數的不確定性評價結果顯示,當閾值概率劃定為0.50、0.75、0.85、0.95時,超閾值概率的區域面積分別占金壇區總面積的16.71%、5.74%、2.84%、1.04%(圖6)。可以看出隨著閾值概率的增加,所劃定的區域面積不斷減小。當閾值概率為0.85時(圖6c),研究區中主城區的西北部、洮湖的周邊以及西部旱地部分地區土壤磷素流失風險等級高具有很大的可能性。對于那些概率值接近50%的區域(圖6a),其評價結果就有很大的不確定性。所以這樣的概率分布能夠給管理部門的決策提供更多空間信息,而不是直接按某個閾值將整個地區機械地劃分為不同等級的污染區域[50]。閾值概率設定為0.85時,圖6c中的劃定的區域就要重點實施管控措施。

圖4 金壇區土壤磷素流失風險等級圖Fig.4 Level of soil phosphorus loss risk in Jintan District

圖5 金壇區土壤磷素流失風險指數(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的臨界概率劃定的土壤磷素流失高風險區域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
本研究結合長三角水網發達、多為平原、縣域景觀異質性小、農田水利工程完善等縣域特征,在傳統分析土壤磷素空間分布格局的基礎上,引入磷肥施用量、距受納水體距離、地表徑流潛力構建了針對長三角典型縣域土壤磷素流失風險的RPI評估模型;同時,我們將金壇區土壤全磷可能的RSGS實現輸入該模型,進而評估了流失風險指數可能的空間分布格局;最后我們依據不同的超標概率閾值劃定了金壇土壤磷素流失風險的關鍵源區即高風險區。結果顯示土壤磷素流失的高風險區主要沿著河流呈現條帶狀及斑塊狀分布,集中分布在洮湖周邊、金壇主城區的東部及西部,較高及高風險區(快速磷指數值大于0.93)的面積占金壇區面積的65.88%。概率閾值分別設定為0.50、0.75、0.85、0.95時,其超標面積占金壇區總面積分別達到16.71%、5.74%、2.84%、1.04%。本研究提供的土壤磷素流失風險信息,有助于決策者從“不確定性”角度為區域土壤磷素的調控,劃定更符合實際限制條件的管控區域。