未來(lái)發(fā)展情景選擇對(duì)農(nóng)業(yè)面源污染風(fēng)險(xiǎn)的影響分析

朱康文，陳玉成，2*，熊海靈，張晟，楊志敏，2，黃磊，2，雷波

（1.西南大學(xué)資源環(huán)境學(xué)院，重慶 400716；2.農(nóng)村清潔工程重慶市工程研究中心，重慶 400716；3.西南大學(xué)計(jì)算機(jī)與信息科學(xué)學(xué)院，重慶 400716；4.重慶市生態(tài)環(huán)境科學(xué)研究院，重慶 401147）

農(nóng)業(yè)面源污染（Agricultural non-point source pollution，ANSP）是農(nóng)業(yè)生產(chǎn)活動(dòng)引起的各種污染物在降雨、徑流等驅(qū)動(dòng)作用下從土壤向水體擴(kuò)散，引起受納水體污染的過(guò)程[1]。ANSP 的研究主要從20 世紀(jì)70 年代開(kāi)始，美國(guó)環(huán)保局為控制日益嚴(yán)重的水體環(huán)境污染，要求各州編制非點(diǎn)源污染控制方案[2]。我國(guó)在20世紀(jì)80年代才開(kāi)始關(guān)注ANSP，但并沒(méi)有采取相應(yīng)的控制措施[3]。ANSP 具有滯后性，其污染凸顯需要一個(gè)過(guò)程，20 世紀(jì)90 年代開(kāi)始我國(guó)農(nóng)村生態(tài)環(huán)境問(wèn)題開(kāi)始顯現(xiàn)。當(dāng)前國(guó)內(nèi)外ANSP 形勢(shì)依然嚴(yán)峻，中國(guó)第一次、第二次全國(guó)污染源普查公報(bào)顯示，農(nóng)業(yè)源總氮（TN）、總磷（TP）分別占水體污染排放總量的比例由2007 年的57.2%、67.3%變?yōu)?017 年的46.5%、67.2%。涪陵區(qū)2018 年統(tǒng)計(jì)年鑒數(shù)據(jù)顯示，化肥、農(nóng)藥施用強(qiáng)度分別為221.19、7.41 kg·hm-2。美國(guó)的調(diào)查和研究指出，15%的排水良好的農(nóng)業(yè)區(qū)其地下水的硝態(tài)氮濃度超過(guò)了飲用水質(zhì)標(biāo)準(zhǔn)[4]，這源于土壤中施用了高含量的氮。歐洲因農(nóng)業(yè)活動(dòng)輸入到北海河口的TN、TP 分別占入海通量的60%和25%[5]。因此為解決ANSP 問(wèn)題各國(guó)學(xué)者均開(kāi)展了大量研究，隨著農(nóng)學(xué)與地理學(xué)的深入融合，ANSP 研究不斷深入，如LULL 等[6]、TAHMASEBI 等[7]、CHENG 等[8]、龍?zhí)煊宓萚9]將地理學(xué)領(lǐng)域方法融入ANSP 研究中，這有利于準(zhǔn)確識(shí)別污染物的輸移過(guò)程、空間變化、風(fēng)險(xiǎn)趨勢(shì)等情況。地理學(xué)融入ANSP 分析側(cè)重于土地利用和景觀格局對(duì)ANSP 的影響研究及地理學(xué)分析方法的應(yīng)用，主要包括ANSP 隨土地利用變化的演變情況[1]，ANSP 對(duì)景觀格局優(yōu)化調(diào)整的響應(yīng)分析[10]，地理學(xué)方法對(duì)ANSP空間演變[9]、負(fù)荷強(qiáng)度[8]、輸移路徑[11]等方面的研究。

綜合已有研究和當(dāng)前ANSP 問(wèn)題在區(qū)縣層面防控的實(shí)際情況發(fā)現(xiàn)：當(dāng)前區(qū)域現(xiàn)狀是土地利用與ANSP 風(fēng)險(xiǎn)之間的關(guān)系研究較多，但對(duì)于未來(lái)土地利用變化對(duì)ANSP 風(fēng)險(xiǎn)演變影響的研究較少，掌握未來(lái)ANSP 風(fēng)險(xiǎn)演變有利于提升區(qū)縣層面的風(fēng)險(xiǎn)防控效果；同時(shí)ANSP 受當(dāng)?shù)亟?jīng)濟(jì)發(fā)展模式選擇影響較大，如自然發(fā)展和耕地保護(hù)兩種模式下地類結(jié)構(gòu)、肥藥施用水平均存在一定差異，因此探討不同模式下的ANSP 演變情況也是區(qū)縣層面需要開(kāi)展的研究；目前已有研究中較為成熟的思路是結(jié)合土地利用數(shù)據(jù)與輸出系數(shù)法進(jìn)行區(qū)域ANSP 演變分析[1]，結(jié)果展示方面以子流域?yàn)榻y(tǒng)計(jì)單元的效果最佳[9]，但是輸出系數(shù)法在測(cè)算時(shí)較少考慮子流域內(nèi)的自然環(huán)境差異，導(dǎo)致結(jié)果難以反映子流域內(nèi)部空間差異情況。因此通過(guò)優(yōu)化改進(jìn)模型達(dá)到區(qū)分子流域內(nèi)部空間異質(zhì)性的目的也是目前研究需要解決的問(wèn)題。

未來(lái)土地利用數(shù)據(jù)模擬方法主要包括CA-Markov模型[12]、神經(jīng)網(wǎng)絡(luò)算法[13]、系統(tǒng)動(dòng)力學(xué)方法[14]、CLUE-S模型[15-16]等，其中以 CLUE-S 和 CA-Markov 模型應(yīng)用較多。因此為豐富區(qū)縣尺度ANSP 風(fēng)險(xiǎn)演變研究工作，本研究利用CLUE-S[17-18]和Markov[16]模型在土地利用情景模擬方面的優(yōu)勢(shì)，以具有生態(tài)地位高、生態(tài)風(fēng)險(xiǎn)大、農(nóng)業(yè)經(jīng)濟(jì)比重高、ANSP風(fēng)險(xiǎn)大等特點(diǎn)的三峽庫(kù)區(qū)典型區(qū)縣——涪陵區(qū)為研究對(duì)象，進(jìn)行自然發(fā)展（ND）和耕地保護(hù)（CP）兩種發(fā)展模式下2020—2030年土地利用覆蓋變化模擬，并融入距水域距離、坡長(zhǎng)坡度因子等對(duì)輸出風(fēng)險(xiǎn)模型進(jìn)行改進(jìn)，以分析區(qū)域內(nèi)各子流域ANSP 風(fēng)險(xiǎn)演變趨勢(shì)，并進(jìn)行輸出風(fēng)險(xiǎn)響應(yīng)分析。研究路徑與方法可很好地支撐區(qū)域的ANSP未來(lái)發(fā)展模式選擇與風(fēng)險(xiǎn)調(diào)控等決策工作，改進(jìn)后的輸出風(fēng)險(xiǎn)模型可較好地區(qū)分子流域內(nèi)部空間異質(zhì)性，對(duì)區(qū)域ANSP 風(fēng)險(xiǎn)防控具有重要的數(shù)據(jù)支撐和借鑒意義。

1 研究方法

1.1 研究區(qū)域概況與數(shù)據(jù)來(lái)源

涪陵區(qū)位于中國(guó)重慶市中部，東經(jīng)106°56′~107°43′，北緯29°21′~30°01′。長(zhǎng)江和烏江在涪陵城區(qū)交匯，區(qū)域內(nèi)匯入長(zhǎng)江的一級(jí)支流有35 條，其中流域面積大于100 km2的河流有烏江、梨香溪、小溪、渠溪河等12 條，可劃分470 個(gè)子流域（圖1，子流域劃分在ArcGIS 軟件中進(jìn)行，詳細(xì)方法見(jiàn)文獻(xiàn)[19]），該區(qū)域在保障長(zhǎng)江中下游和三峽庫(kù)區(qū)生態(tài)安全方面具有重要地位[20]。此外涪陵區(qū)生態(tài)環(huán)境較好，包括大木山自然保護(hù)區(qū)、江東桫欏自然保護(hù)區(qū)、武陵山森林公園、烏江森林公園等自然保護(hù)地。涪陵區(qū)經(jīng)濟(jì)發(fā)展水平較高，2018年全區(qū)生產(chǎn)總值為1 076.13億元，第一、二、三產(chǎn)業(yè)比例為6.1∶57.0∶36.9。同時(shí)涪陵區(qū)農(nóng)業(yè)較為發(fā)達(dá)，2018 年涪陵區(qū)農(nóng)林牧漁業(yè)總產(chǎn)值為96.67 億元，全年糧食播種面積9.18萬(wàn)hm2，全年糧食總產(chǎn)量44.03萬(wàn)t。

以上表明涪陵區(qū)具有生態(tài)地位高、農(nóng)業(yè)占比高和ANSP 風(fēng)險(xiǎn)高等特點(diǎn)，其區(qū)域內(nèi)的ANSP 問(wèn)題備受學(xué)者關(guān)注，LIU 等[21]、梁增芳等[22]、WANG 等[11]在涪陵區(qū)部分小流域已開(kāi)展了ANSP 相關(guān)研究工作。因此涪陵區(qū)的ANSP 時(shí)空演變及未來(lái)發(fā)展趨勢(shì)研究具有重要的現(xiàn)實(shí)意義。

1.1.1 土地利用現(xiàn)狀分析

涪陵區(qū)2010—2015 年的土地利用數(shù)據(jù)反映出其地類比例變化較大，2010 年和2015 年耕地、林地、草地、水域、建設(shè)用地、其他用地面積比例分別為52.85%、41.94%、0.14%、3.13%、1.91%、0.03% 和36.96%、53.13%、1.20%、3.50%、5.21%、0（圖2）。采用轉(zhuǎn)移矩陣對(duì)各地類轉(zhuǎn)化情況分析發(fā)現(xiàn)，耕地、林地轉(zhuǎn)出較多。耕地主要存在林地轉(zhuǎn)入，比例為7.77%；林地主要存在耕地轉(zhuǎn)入，比例為29.21%；草地主要存在耕地、林地轉(zhuǎn)入，比例分別為28.76%、65.62%；水域主要存在耕地、林地轉(zhuǎn)入，比例分別為13.05%、7.31%；建設(shè)用地主要存在耕地、林地轉(zhuǎn)入，比例分別為52.08%、13.24%?？傮w上林地和水域空間格局變化較小，建設(shè)用地外擴(kuò)現(xiàn)象明顯，水田、旱地、建設(shè)用地主要分布在長(zhǎng)江、烏江沿岸，林地主要分布在中部黃草山和南部大木山、武陵山區(qū)域。

1.1.2 數(shù)據(jù)來(lái)源

數(shù)據(jù)類型主要分為面板數(shù)據(jù)、遙感數(shù)據(jù)和統(tǒng)計(jì)數(shù)據(jù)，遙感數(shù)據(jù)包括2010年和2015年的土地利用數(shù)據(jù)、土壤類型及質(zhì)地?cái)?shù)據(jù)、數(shù)字高程模型（DEM）地形數(shù)據(jù)和坡度數(shù)據(jù)。土地利用數(shù)據(jù)來(lái)源于中國(guó)生態(tài)環(huán)境遙感評(píng)估項(xiàng)目[23]，當(dāng)前此數(shù)據(jù)集僅發(fā)布到2015 年的數(shù)據(jù)。DEM 數(shù)據(jù)來(lái)源于資源環(huán)境數(shù)據(jù)云平臺(tái)（http：//www.resdc.cn/），坡度數(shù)據(jù)利用GIS 測(cè)算獲取，所有數(shù)據(jù)在GIS軟件中進(jìn)行重采樣后分辨率為30 m。

1.2 技術(shù)路線與研究方法

研究采用 CLUE-S 和 Markov 模型分 ND、CP 兩種發(fā)展模式進(jìn)行未來(lái)土地利用變化模擬，采用Kappa 系數(shù)進(jìn)行初步精度分析[24]。同時(shí)根據(jù)孫毅中等[25]采用的Pontius，通過(guò)統(tǒng)計(jì)模擬結(jié)果和實(shí)際數(shù)據(jù)之間差異區(qū)和共性區(qū)構(gòu)建的FoM 指標(biāo)再次進(jìn)行空間模擬精度驗(yàn)證，一般認(rèn)為FoM大于0.21時(shí)具有較高精度。在此基礎(chǔ)上開(kāi)展子流域劃定，結(jié)合改進(jìn)的輸出風(fēng)險(xiǎn)模型進(jìn)行ANSP 輸出風(fēng)險(xiǎn)演化分析，并分析ANSP 風(fēng)險(xiǎn)對(duì)輸出系數(shù)變化的響應(yīng)情況（圖3）。

1.2.1 基于CLUE-S和Markov模型的土地利用情景預(yù)測(cè)

研究采用 CLUE-S 和Markov 模型分 ND、CP 兩種發(fā)展模式進(jìn)行未來(lái)土地利用變化模擬。以2010 年為基礎(chǔ)，模擬2015年土地利用數(shù)據(jù)，并與2015年土地利用現(xiàn)狀數(shù)據(jù)進(jìn)行對(duì)比以分析模型在研究區(qū)的模擬精度，開(kāi)展ND、CP 兩種發(fā)展情景下2020、2025、2030 年研究區(qū)的土地利用模擬。

（1）CLUE-S模型及輸入要素介紹

CLUE-S 模型是 Peter-verburg 團(tuán)隊(duì)在 CLUE 模型的基礎(chǔ)上改進(jìn)的土地利用變化模擬模型[26-27]，CLUES 模型主要由需求模塊和空間分配模塊構(gòu)成[27-28]。需求模塊用于計(jì)算研究區(qū)域內(nèi)各土地利用類型在目標(biāo)年的需求量，需要通過(guò)外部模型或方法進(jìn)行計(jì)算，本研究將Markov 模型引入用來(lái)計(jì)算各地類需求量[12]?？臻g分配模塊輸入要素主要包括4 個(gè)方面：土地需求文件、各土地利用類型的轉(zhuǎn)換規(guī)則及轉(zhuǎn)換彈性系數(shù)、各土地利用類型與各驅(qū)動(dòng)因子之間的關(guān)系、限制區(qū)文件，具體模擬過(guò)程詳見(jiàn)文獻(xiàn)[17]。

（2）Markov模型測(cè)算土地利用需求

土地利用需求是指計(jì)算模擬期間內(nèi)每年各土地利用類型的需求量。Markov模型是狀態(tài)離散和時(shí)間離散的動(dòng)力學(xué)模型，是主要針對(duì)數(shù)量上的預(yù)測(cè)[29]。Markov模型由于具有長(zhǎng)期預(yù)測(cè)的能力而被廣泛應(yīng)用，例如HULST[30]、LIPPE 等[31]、劉春艷等[32]曾將 Markov 模型應(yīng)用在解決植被生態(tài)或生物多樣性功能變化預(yù)測(cè)方面等，表明Markov 模型在未來(lái)土地利用類型預(yù)測(cè)方面具有較好的效果。本研究在GIS、IDRISI（Markov 模塊）軟件的支持下，通過(guò)2010年和2015年土地覆蓋數(shù)據(jù)，結(jié)合Markov 預(yù)測(cè)的需求量以及涪陵區(qū)未來(lái)土地利用規(guī)劃、城鄉(xiāng)發(fā)展規(guī)劃等，得到CLUE-S 模型需要的 2020、2025、2030 年 ND 情景模式下土地利用類型面積需求文件。CP情景在ND 情景基礎(chǔ)上將耕地、草地、其他用地降速分別調(diào)減10%，林地、水域、建設(shè)用地增速分別調(diào)減10%以平衡面積變化。

（3）轉(zhuǎn)換系數(shù)

根據(jù)不同土地利用類型的歷史變化情況和未來(lái)實(shí)際土地的規(guī)劃情況設(shè)置各土地利用類型的穩(wěn)定程度，由模型參數(shù)ELAS（轉(zhuǎn)換系數(shù)）[33]定義：①一般情況下，不會(huì)轉(zhuǎn)換為其他土地利用類型的ELAS 設(shè)為1；②極易變化的土地利用類型的ELAS 設(shè)為0；③0

（4）土地利用類型與驅(qū)動(dòng)因子之間的關(guān)系方程構(gòu)建

二項(xiàng)Logistic 回歸是目前最常用的二項(xiàng)選擇模型，其因變量只有0和1兩個(gè)值，選取某一土地利用類型作為因變量，驅(qū)動(dòng)因子作為自變量[28]，回歸分析結(jié)果表示不同土地利用類型與不同驅(qū)動(dòng)因子之間關(guān)系，具體公式詳見(jiàn)文獻(xiàn)[17]。本研究選取離城鎮(zhèn)道路距離（sc1gr0）、離鄉(xiāng)村道路距離（sc1gr1）、離城市距離（sc1gr2）、離建制鎮(zhèn)距離（sc1gr3）、離村莊居民點(diǎn)距離（sc1gr4）、離水域距離（sc1gr5）、高程（sc1gr6）、坡度（sc1gr7）、起伏度（sc1gr8）等9個(gè)因子作為土地利用變化的驅(qū)動(dòng)因子。使用SPSS 軟件將6 種土地利用類型與9種驅(qū)動(dòng)因子進(jìn)行Logistic 回歸分析，得到表1中回歸方程的系數(shù)和常數(shù)，用Pontius 等提出的ROC（Receiver operating characteristic）方法對(duì)回歸結(jié)果進(jìn)行檢驗(yàn)[26，36]。ROC值均大于0.7，且大部分高于0.8，說(shuō)明土地利用類型與驅(qū)動(dòng)因子之間的相關(guān)性較好，滿足模擬的條件。

表1 Logistic回歸分析結(jié)果Table 1 Results of logistic regression analysis

（5）創(chuàng)建限制區(qū)數(shù)據(jù)

限制區(qū)可使模擬結(jié)果更符合研究區(qū)實(shí)際情況，一般將土地利用結(jié)構(gòu)長(zhǎng)期保持穩(wěn)定的區(qū)域作為限制區(qū)，包括森林公園、自然保護(hù)區(qū)等，其對(duì)于未來(lái)土地利用模擬有很好的修正作用[37]。已有研究較多考慮的限制因子主要包括國(guó)家公園、自然保護(hù)區(qū)等生態(tài)重點(diǎn)保護(hù)區(qū)，也有研究者將耕地優(yōu)化作為限制因子，如坡度30°以上的區(qū)域禁止轉(zhuǎn)化為耕地和建設(shè)用地[17]。張丁軒等[38]將基本農(nóng)田作為模擬的限制區(qū)域，WAIYASUSRI等[24]將野生動(dòng)物保護(hù)區(qū)作為一個(gè)有限的區(qū)域來(lái)模擬Huai Thap Salao 流域的土地利用。本研究以生態(tài)保護(hù)紅線作為限制區(qū)，生態(tài)保護(hù)紅線是繼“18 億畝耕地紅線”（1 畝=0.067 hm2）后另一條被提到國(guó)家層面的“生命線”，涪陵區(qū)內(nèi)涵蓋自然保護(hù)區(qū)、森林公園等，區(qū)域內(nèi)土地利用長(zhǎng)期保持穩(wěn)定，符合涪陵區(qū)實(shí)際情況。涪陵區(qū)限制區(qū)面積約占涪陵區(qū)總面積的8%，主要分布在西南部的大木山地區(qū)、江東街道中部、荔枝街道中部、荔枝和白濤交界處以及長(zhǎng)江沿岸等區(qū)域。

1.2.2 改進(jìn)的輸出風(fēng)險(xiǎn)模型

改進(jìn)后的輸出風(fēng)險(xiǎn)模型計(jì)算公式如下：

式中：j為污染物類型，本研究將N、P作為主要污染物類型進(jìn)行研究；i為研究區(qū)主要的土地利用類型，包括耕地、林地、草地、水域、建設(shè)用地、其他用地；Lj表示流域內(nèi)污染物j帶權(quán)重的平均輸出系數(shù)，kg·hm-2·a-1；Eij表示污染物j在第i種土地利用類型上的輸出系數(shù)，kg·hm-2·a-1；Ai表示流域第i種地類的面積，hm2；A為流域內(nèi)的總面積，hm2；Eh為林地的最大輸出系數(shù)，由于是帶權(quán)重計(jì)算，這里假設(shè)權(quán)重為1 時(shí)作為林地最大輸出系數(shù)，kg·hm-2·a-1；Pj為污染物j的輸出風(fēng)險(xiǎn)概率[39-40]。林地相比其他用地類型，在自然狀態(tài)下一般不產(chǎn)生人類活動(dòng)影響下的面源污染，因此以林地為參考，將林地的最大輸出系數(shù)作為評(píng)價(jià)閾值，當(dāng)某地類輸出系數(shù)高于閾值時(shí)證明該地類存在面源污染風(fēng)險(xiǎn)，Pj介于-1 到1 之間，小于0 的值賦值為0。由于輸出系數(shù)重點(diǎn)考慮輸入方面的壓力，對(duì)于輸移過(guò)程的難易程度未體現(xiàn)，考慮到增加的因子應(yīng)是在長(zhǎng)時(shí)間內(nèi)保持穩(wěn)定的因子。因此采用GWd（歸一化處理的距水域距離值）和GLS（歸一化處理的坡長(zhǎng)坡度值）作為校正因子以體現(xiàn)驅(qū)動(dòng)動(dòng)能相關(guān)因子的作用。R為優(yōu)化的輸出風(fēng)險(xiǎn)指數(shù)，R風(fēng)險(xiǎn)分級(jí)時(shí)以0.1 為間隔分為10 個(gè)等級(jí)，并以Ⅰ、Ⅱ、Ⅲ、…、Ⅹ表示，等級(jí)越高表示風(fēng)險(xiǎn)越高。

各土地利用類型的輸出系數(shù)是研究的關(guān)鍵，建設(shè)用地參考龍?zhí)煊宓萚9]在三峽庫(kù)區(qū)進(jìn)行的相關(guān)研究；林地、草地參考龍?zhí)煊宓萚9]、楊彥蘭等[41]在三峽庫(kù)區(qū)進(jìn)行的相關(guān)研究，取二者研究的均值；耕地參考唐肖陽(yáng)等[42]的計(jì)算方法，取目前不同區(qū)域相關(guān)研究的均值；水域和其他用地不計(jì)算污染排放（表2）。為了判斷肥藥施用變化對(duì)ANSP 風(fēng)險(xiǎn)的影響，將目前輸出系數(shù)中耕地輸出系數(shù)分別按±5%和±10%進(jìn)行測(cè)算。

表2 涪陵區(qū)土地利用類型的輸出系數(shù)（kg·hm-2·a-1）Table 2 Output coefficients of different land use types in Fuling District（kg·hm-2·a-1）

2 結(jié)果與分析

2.1 模擬精度與未來(lái)土地利用情景模擬結(jié)果

2.1.1 模擬2015年土地利用數(shù)據(jù)進(jìn)行精度驗(yàn)證

將2015年模擬結(jié)果與2015年土地利用現(xiàn)狀數(shù)據(jù)進(jìn)行對(duì)比，各地類面積總量上，耕地、林地、草地、水域、建設(shè)用地、其他用地誤差分別為8.87%、5.42%、5.12%、0.56%、7.24%、0.09%，總量模擬精度均高于90%；空間布局上，疊加2015 年土地利用模擬數(shù)據(jù)與土地利用現(xiàn)狀數(shù)據(jù)發(fā)現(xiàn)，Kappa系數(shù)為0.75，位于0.6~0.8 之間，屬于高度一致性。FoM 指標(biāo)結(jié)果為0.23，大于0.21，說(shuō)明具有較高精度。因此總量和空間布局精度兩方面結(jié)果表明模擬精度可以滿足研究需要。

2.1.2 兩種情景下土地利用變化情況

2020、2025、2030 年 ND 和 CP 兩種發(fā)展情景下涪陵區(qū)土地利用空間分布的模擬結(jié)果如圖4 所示。CP情景的耕地保有量高于ND 情景，CP情景的建設(shè)用地面積低于ND 情景?？傮w上各土地利用類型面積變化規(guī)律符合情景設(shè)定情況。CP 情景下城市發(fā)展更側(cè)重于侵占林地、草地，符合實(shí)際發(fā)展規(guī)律。兩種情景下各用地類型的空間分布基本一致：（1）耕地廣泛分布，但主要分布于長(zhǎng)江以北區(qū)域。（2）林地分布也比較廣泛，但是主要以黃草山、大木山、焦石鎮(zhèn)、江東街道等區(qū)域?yàn)橹鳌＃?）水域主要以長(zhǎng)江、烏江及其支流為主。（4）建設(shè)用地主要分布在涪陵城區(qū)、李渡街道、白濤街道等區(qū)域，這與涪陵區(qū)老城、新城、工業(yè)園區(qū)規(guī)劃布局一致。（5）草地、其他用地分布面積較少且分散。

2.2 2010—2015年輸出風(fēng)險(xiǎn)分布時(shí)空變化

2010—2015 年期間風(fēng)險(xiǎn)概率明顯降低，尤其是靠近大木山、武陵山區(qū)域的鎮(zhèn)街和子流域。風(fēng)險(xiǎn)概率等級(jí)變化分析發(fā)現(xiàn)，區(qū)域內(nèi)存在非常明顯的等級(jí)降低趨勢(shì)。從各鎮(zhèn)街不同等級(jí)TN輸出風(fēng)險(xiǎn)面積比例構(gòu)成來(lái)看（TP 趨勢(shì)一致），總體上各鎮(zhèn)街TN 的ANSP 風(fēng)險(xiǎn)等級(jí)集中在Ⅲ~Ⅶ級(jí)，Ⅲ級(jí)及以下等級(jí)占比較高的鎮(zhèn)街有大木鄉(xiāng)、武陵山鄉(xiāng)、藺市鎮(zhèn)和龍?zhí)舵?zhèn)，表明這些鎮(zhèn)街的低等級(jí)風(fēng)險(xiǎn)比例較高；Ⅶ級(jí)占比較高的鎮(zhèn)街有清溪鎮(zhèn)、江東街道、荔枝街道、藺市鎮(zhèn)、涪陵城市區(qū)等，表明這些鎮(zhèn)街的高風(fēng)險(xiǎn)等級(jí)比例較高；Ⅷ級(jí)及以上等級(jí)占比較高的鎮(zhèn)街有江東街道、南沱鎮(zhèn)、白濤街道、新妙鎮(zhèn)、青羊鎮(zhèn)、藺市鎮(zhèn)、明家鄉(xiāng)和龍橋街道，其他鎮(zhèn)街比例較低。2010 年和2015 年TN 均僅在白濤街道存在分別0.02%和0.01%的Ⅸ級(jí)風(fēng)險(xiǎn)分布；2010 年在南沱鎮(zhèn)、江東街道、白濤街道存在6%以上面積比例的Ⅷ級(jí)風(fēng)險(xiǎn)分布，2015 年僅有南沱鎮(zhèn)存在6%以上面積比例的Ⅷ級(jí)風(fēng)險(xiǎn)分布；2010 年在清溪鎮(zhèn)、江東街道存在40%以上面積比例的Ⅶ級(jí)風(fēng)險(xiǎn)分布，2015 年僅有江東街道存在30%以上面積比例分布。數(shù)據(jù)分析發(fā)現(xiàn)，風(fēng)險(xiǎn)降低的主要原因是區(qū)域存在較為明顯的林地面積上升、耕地面積減少的趨勢(shì)（圖5）。

2.3 2020、2025、2030年多情景下輸出風(fēng)險(xiǎn)演變分析

2020、2025、2030 年 ND、CP 情景模式下 TN 和 TP輸出風(fēng)險(xiǎn)概率總體上變化較小，風(fēng)險(xiǎn)概率處于降低的趨勢(shì)中，CP情景下風(fēng)險(xiǎn)概率總體高于ND情景（圖6）。TP 風(fēng)險(xiǎn)概率高于TN 風(fēng)險(xiǎn)但趨勢(shì)一致，因此以TP 風(fēng)險(xiǎn)概率作為各鎮(zhèn)街的風(fēng)險(xiǎn)概率等級(jí)進(jìn)行分析。從空間分布上可以發(fā)現(xiàn)，兩種模式下在荔枝街道、焦石鎮(zhèn)、白濤街道、大木鄉(xiāng)、武陵山鄉(xiāng)等區(qū)域的風(fēng)險(xiǎn)明顯低于其他鎮(zhèn)街，同時(shí)長(zhǎng)江以北、烏江以西區(qū)域風(fēng)險(xiǎn)明顯偏高，這與區(qū)域自然資源本底條件吻合。兩種發(fā)展模式下各風(fēng)險(xiǎn)概率等級(jí)的空間差異較小，但CP 情景的風(fēng)險(xiǎn)等級(jí)高于ND情景。

從2010—2030 年各鎮(zhèn)街Ⅶ級(jí)以上等級(jí)風(fēng)險(xiǎn)統(tǒng)計(jì)結(jié)果來(lái)看，2010—2015 年風(fēng)險(xiǎn)降幅明顯，2015 年后風(fēng)險(xiǎn)等級(jí)降低速率明顯放緩，例如ND 情景下，2010、2015、2020、2025、2030 年各年份清溪鎮(zhèn)Ⅶ以上等級(jí)風(fēng)險(xiǎn)比例分別為46.20%、29.69%、25.72%、22.20%、20.00%。圖7 明顯可看出江東街道、清溪鎮(zhèn)、涪陵城市區(qū)、藺市鎮(zhèn)等Ⅶ級(jí)以上等級(jí)風(fēng)險(xiǎn)較多，龍?zhí)舵?zhèn)、大木鄉(xiāng)、武陵山鄉(xiāng)、羅云鄉(xiāng)、焦石鎮(zhèn)等高風(fēng)險(xiǎn)分布比例較少。同時(shí)從ND 和CP 兩種情景來(lái)看，CP 情景下各鎮(zhèn)街高風(fēng)險(xiǎn)分布比例明顯高于ND情景。

2.4 輸出風(fēng)險(xiǎn)對(duì)輸出系數(shù)變化的響應(yīng)分析

輸出系數(shù)對(duì)區(qū)域輸出風(fēng)險(xiǎn)概率有一定影響，可在一定程度上反映區(qū)域施肥狀況。鑒于本文重點(diǎn)分析輸出風(fēng)險(xiǎn)對(duì)輸出系數(shù)變化的響應(yīng)關(guān)系，因此測(cè)算中按照改進(jìn)前的輸出風(fēng)險(xiǎn)模型進(jìn)行測(cè)算，以便于單純分析肥藥用量變化可能帶來(lái)的影響，更好地統(tǒng)計(jì)子流域的變化情況。假設(shè)耕地的TP 輸出系數(shù)分別按照±5%和±10%進(jìn)行調(diào)整，其他土地利用類型的輸出系數(shù)保持不變，對(duì)ND情景下TP輸出風(fēng)險(xiǎn)概率變化情況進(jìn)行分析。

結(jié)果表明ND 情景下，當(dāng)輸出系數(shù)分別上調(diào)5%和 10%時(shí)，2020、2025、2030 年的 TP 輸出風(fēng)險(xiǎn)概率在Ⅸ級(jí)及以上等級(jí)的子流域個(gè)數(shù)分別增加6、9 個(gè)，7、14個(gè)，4、11 個(gè)。當(dāng)輸出系數(shù)分別下調(diào)5%和10%時(shí)，2020、2025、2030年的TP輸出風(fēng)險(xiǎn)概率在Ⅸ級(jí)及以上等級(jí)的子流域個(gè)數(shù)分別減少4、15個(gè)，7、10個(gè)，3、8個(gè)。從圖8 分析可以發(fā)現(xiàn)，當(dāng)輸出系數(shù)上調(diào)時(shí)低風(fēng)險(xiǎn)等級(jí)堆積數(shù)降低，反之低風(fēng)險(xiǎn)等級(jí)堆積數(shù)上升，結(jié)果表明區(qū)域TP 輸出風(fēng)險(xiǎn)概率等級(jí)與耕地輸出系數(shù)密切相關(guān)，同時(shí)表明其與施肥量密切相關(guān)。本研究可以通過(guò)空間變化很好地判斷輸出系數(shù)變化對(duì)哪些具體子流域影響較大，識(shí)別出的子流域是區(qū)域土地利用優(yōu)化和化肥施用控制的關(guān)鍵區(qū)域，如第3、75、104、141、202、211、259、292、330、398、461 號(hào)等子流域?qū)敵鱿禂?shù)調(diào)減作出了波動(dòng)響應(yīng)，第79、118、199、288、315、333、381、400、454 號(hào)等子流域?qū)敵鱿禂?shù)調(diào)增作出了波動(dòng)響應(yīng)。這些子流域主要位于荔枝街道、馬武鎮(zhèn)、白濤街道、焦石鎮(zhèn)、增福鄉(xiāng)、江東街道、同樂(lè)鄉(xiāng)、龍?zhí)舵?zhèn)等區(qū)域，表明輸出系數(shù)調(diào)整（即農(nóng)藥化肥用量調(diào)整）后風(fēng)險(xiǎn)等級(jí)變化較快的區(qū)域主要位于長(zhǎng)江以南。

3 討論

3.1 子流域角度的測(cè)算利于展現(xiàn)區(qū)域風(fēng)險(xiǎn)實(shí)際

已有研究一般以行政單元或大流域?yàn)樽钚⊙芯繂卧_(kāi)展輸出系數(shù)測(cè)算。例如，SHEN等[43]在三峽庫(kù)區(qū)的研究以大流域?yàn)檠芯繂卧M(jìn)行風(fēng)險(xiǎn)測(cè)算，DURAND等[44]開(kāi)展了流域尺度的氮輸入變化引起的污染變化研究。本研究在衡量輸出風(fēng)險(xiǎn)概率時(shí)，選取最能反映區(qū)域地理和流域?qū)嶋H狀況的子流域范圍進(jìn)行分析，這種測(cè)算方式使測(cè)算結(jié)果更接近區(qū)域風(fēng)險(xiǎn)實(shí)際，便于研究者和管理者清晰地認(rèn)識(shí)視野范圍內(nèi)的各子流域的TN、TP輸出風(fēng)險(xiǎn)概率演變情況，有利于土地利用優(yōu)化調(diào)整或肥藥施用方面的管控。

3.2 土地利用演變是產(chǎn)出風(fēng)險(xiǎn)變化的主要驅(qū)動(dòng)力

從產(chǎn)出風(fēng)險(xiǎn)測(cè)度結(jié)果來(lái)看，涪陵大部分地區(qū)的風(fēng)險(xiǎn)相對(duì)較高，符合區(qū)域農(nóng)業(yè)發(fā)展的實(shí)際情況，這與CHEN等[45]、ZHU等[46]的研究基本一致。在不同時(shí)期，子流域產(chǎn)出風(fēng)險(xiǎn)的降低或增加與區(qū)域土地利用類型的構(gòu)成密切相關(guān)，輸出風(fēng)險(xiǎn)高的主要原因是耕地或建設(shè)用地比重大。根據(jù)土地利用數(shù)據(jù)統(tǒng)計(jì)，2010—2015年存在較多的耕地轉(zhuǎn)為林地，主要原因是該區(qū)坡度基本在20°左右，耕作條件不適宜，導(dǎo)致耕地向林地演變。同時(shí)，我國(guó)農(nóng)業(yè)人口外流現(xiàn)象普遍存在，導(dǎo)致山區(qū)出現(xiàn)耕地拋荒現(xiàn)象，這也是輸出風(fēng)險(xiǎn)整體降低的原因之一。在未來(lái)情景分析中，CP情景相對(duì)ND 情景耕地面積較多，這是CP情景產(chǎn)出風(fēng)險(xiǎn)普遍高于ND情景的主要原因之一。

3.3 改進(jìn)的輸出風(fēng)險(xiǎn)模型更能體現(xiàn)污染輸移過(guò)程差異

改進(jìn)的輸出風(fēng)險(xiǎn)模型引入的距離水域距離因子和坡度坡長(zhǎng)因子是影響中小尺度ANSP 從“源”到“匯”的重要輸移介質(zhì)。王琦等[47]在查干湖匯水區(qū)開(kāi)展的ANSP 風(fēng)險(xiǎn)識(shí)別中引入了距離水域距離和坡度等因子，反映出這些因子對(duì)ANSP 影響明顯。輸出風(fēng)險(xiǎn)模型測(cè)算原理是以統(tǒng)計(jì)單元為基礎(chǔ)，對(duì)于統(tǒng)計(jì)單元內(nèi)部差異未能體現(xiàn)，改進(jìn)的輸出風(fēng)險(xiǎn)模型可以很好地展現(xiàn)區(qū)域內(nèi)部的風(fēng)險(xiǎn)等級(jí)差異。

研究路徑與方法在未來(lái)可以針對(duì)識(shí)別出的鄉(xiāng)鎮(zhèn)、子流域進(jìn)行土地利用結(jié)構(gòu)優(yōu)化或肥藥施用水平管控，達(dá)到ANSP 風(fēng)險(xiǎn)防控的目的。由于本研究在輸出系數(shù)調(diào)整中未考慮作物類型，未來(lái)將進(jìn)行深入分析。

4 結(jié)論

（1）通過(guò)結(jié)合CLUE-S 與Markov模型的優(yōu)勢(shì)在涪陵區(qū)開(kāi)展自然發(fā)展和耕地保護(hù)兩種情景下2020、2025、2030 年土地利用空間模擬研究，取得了較好的效果。

（2）涪陵區(qū)2010—2015 年間農(nóng)業(yè)面源污染風(fēng)險(xiǎn)等級(jí)明顯降低，但對(duì)于南沱鎮(zhèn)、江東街道、白濤街道等風(fēng)險(xiǎn)較高的鎮(zhèn)街在管理中應(yīng)重點(diǎn)關(guān)注區(qū)域肥藥施用強(qiáng)度是否合理，在加強(qiáng)宣傳和培訓(xùn)提升肥藥利用率的同時(shí)降低流失率。

（3）未來(lái)的情景模擬有利于探尋區(qū)域不同發(fā)展情景下輸出風(fēng)險(xiǎn)概率與土地利用變化之間的響應(yīng)關(guān)系，以及各區(qū)域風(fēng)險(xiǎn)概率演變情況。針對(duì)識(shí)別出的未來(lái)高風(fēng)險(xiǎn)分布較多的區(qū)域應(yīng)加大肥藥施用強(qiáng)度管控和景觀格局優(yōu)化力度，有條件的地區(qū)可構(gòu)建林草緩沖帶或開(kāi)展示范區(qū)建設(shè)以提升農(nóng)業(yè)面源污染消減能力。

（4）子流域輸出風(fēng)險(xiǎn)等級(jí)對(duì)于耕地輸出系數(shù)變化（即化肥施用水平調(diào)整）具有很好的響應(yīng)關(guān)系。自然發(fā)展情景下部分子流域?qū)敵鱿禂?shù)調(diào)減作出了波動(dòng)響應(yīng)，這些子流域是風(fēng)險(xiǎn)防控的關(guān)鍵區(qū)域，在未來(lái)管理中應(yīng)深入調(diào)研分析，提出合理的肥藥施用水平以達(dá)到降低農(nóng)業(yè)面源污染風(fēng)險(xiǎn)的目的。