馬金寶，譚亞男，袁建偉，常曉敏

(1.北京中水利德科技發(fā)展有限公司，北京 100071；2.中國水利水電科學(xué)研究院，北京 100038；3.長江生態(tài)環(huán)保集團有限公司，湖北 武漢 430062)

1 研究背景

農(nóng)業(yè)面源污染由于其隨機性強、時空范圍大、影響周期長、監(jiān)測難度大等特殊性，已成為制約我國農(nóng)業(yè)經(jīng)濟可持續(xù)發(fā)展的重要因素之一[1-2]。已有研究表明，農(nóng)業(yè)面源污染將大量氮、磷、農(nóng)藥等污染物帶入水體，對區(qū)域水環(huán)境和人類健康產(chǎn)生嚴重危害，且已成為中國流域性水體污染的重要來源[3-5]。第二次全國污染源普查公報數(shù)據(jù)顯示，農(nóng)業(yè)面源COD、TN和TP排放量約占到全國排放量的50%、47%和67%，農(nóng)業(yè)源污染物的占比仍然很高，其中種植業(yè)與畜禽養(yǎng)殖業(yè)的貢獻占到農(nóng)業(yè)源的90%以上[6]。農(nóng)業(yè)面源污染已成為危害生態(tài)環(huán)境的一個重要因素。因此，農(nóng)業(yè)面源污染防控與治理是實現(xiàn)農(nóng)業(yè)綠色發(fā)展的本質(zhì)要求，而農(nóng)業(yè)面源污染排放負荷定量估算、污染主要來源及關(guān)鍵區(qū)域的識別是農(nóng)業(yè)面源污染防控與治理的重要步驟[7]，對區(qū)域農(nóng)業(yè)可持續(xù)發(fā)展具有重要意義[8]。目前，針對農(nóng)業(yè)面源污染的負荷估算、主要影響因素分析[9-10]、時空特征分布[11]等已開展了的大量的研究。農(nóng)業(yè)面源污染的影響因素眾多且復(fù)雜，在時間尺度上與降雨量有直接的正相關(guān)關(guān)系[12]，Lu等指出長江上游流域顆粒態(tài)磷輸出負荷在降雨集中的雨季最高[13]。在空間尺度上，不同景觀的分布格局對面源污染的形成具有不同的作用，從而影響面源污染負荷產(chǎn)生和輸移的空間分布特征。面源污染負荷定量估算是其時空分布特征研究的基礎(chǔ)，目前國內(nèi)外學(xué)者在農(nóng)業(yè)面源污染負荷估算方面應(yīng)用較多的有排污系數(shù)法[14-15]、輸出系數(shù)法[16]、經(jīng)驗系數(shù)法[17]、SWAT模型[18]等。但關(guān)于長時間序列農(nóng)業(yè)面源污染排放負荷定量估算及特征解析，主要污染源貢獻及污染區(qū)劃分等還需進一步深入研究。

沂河主要干流位于山東省臨沂市，沂河水系是臨沂市的重要水源地，反之，臨沂市水生態(tài)環(huán)境對沂河下游水環(huán)境有重要的影響。由于沂河流域水系比較發(fā)達，使得農(nóng)業(yè)面源中的化肥農(nóng)藥、生活污水、畜禽排泄等污染物很容易隨地表徑流流失到周邊水體。在生態(tài)環(huán)境保護和高質(zhì)量發(fā)展背景下，協(xié)同推進農(nóng)業(yè)生產(chǎn)與農(nóng)業(yè)面源污染防治是沂河流域發(fā)展農(nóng)業(yè)要面對的現(xiàn)實之一。目前關(guān)于臨沂市農(nóng)業(yè)面源污染定量估算及特征解析的研究相對較少，臨沂市農(nóng)業(yè)面源污染時空動態(tài)變化及污染來源定量解析尚不清晰。本文以山東省臨沂市為研究對象，采用GIS與排污系數(shù)相結(jié)合的方法，定量估算臨沂市不同行政單元2000—2019年TN、TP、COD的排放總量及排放強度。定性和定量的分析研究區(qū)農(nóng)業(yè)面源污染時空變化特征及污染來源貢獻，揭示不同類型農(nóng)業(yè)源污染物排放規(guī)律及其區(qū)域差異，根據(jù)污染物排放程度對研究進行分區(qū)。研究成果可為協(xié)調(diào)區(qū)域農(nóng)業(yè)發(fā)展，制定區(qū)域農(nóng)業(yè)面源污染防控策略提供科學(xué)依據(jù)。

2 研究區(qū)域概況及數(shù)據(jù)來源

2.1 研究區(qū)概況沂河是淮河流域沂沭泗水系中一大支流，發(fā)源于魯南山區(qū)，全長約333 km，流域面積11 820 km2，多年平均徑流量達到36.06億m3。其中臨沂市以上河道長約223 km，沂河區(qū)(臨沂段)面積為9383 km2，多年平均地表水資源量24.30億m3[19]。本研究以臨沂市為研究對象，轄3區(qū)9縣，總面積為17 185 km2，屬于魯東南低山丘陵區(qū)(圖1)。全市水系發(fā)育，河流眾多，呈脈狀分布，水資源利用程度較高，較容易發(fā)生污染現(xiàn)象，加之近年來經(jīng)濟發(fā)展迅速，從而給水資源和水生態(tài)環(huán)境帶來了嚴峻的挑戰(zhàn)。

圖1 研究區(qū)地理位置圖Fig.1 Location map of the study area

2.2 數(shù)據(jù)來源本研究涉及的農(nóng)作物播種面積、作物種植類型、農(nóng)作物產(chǎn)量、化肥農(nóng)藥施用量、畜禽養(yǎng)殖數(shù)量、農(nóng)村人口等數(shù)據(jù)來源于2000—2019年《臨沂市統(tǒng)計年鑒》。農(nóng)村生活、農(nóng)業(yè)種植、畜禽養(yǎng)殖三類污染源中污染物TN、TP、COD排放量的核算方法及相關(guān)系數(shù)參考中華人民共和國生態(tài)環(huán)境部2021年6月最新發(fā)布的《排放源統(tǒng)計調(diào)查產(chǎn)排污核算方法和系數(shù)手冊》(簡稱“手冊”)[20]。

3 研究方法

3.1 農(nóng)村生活源污染負荷量估算根據(jù)《手冊》中農(nóng)村生活源污染物核算方法[20]，對研究區(qū)農(nóng)村生活污水污染物排放量進行估算，農(nóng)村生活污水污染物排放量=污染物產(chǎn)生量(1-對生活的污水進行處理的行政村比例污染物綜合去除率)，其中污染物產(chǎn)生量=農(nóng)村常駐人口(萬人)污染物產(chǎn)污強度(g/人/d)365 d/100。山東省臨沂市農(nóng)村生活污水污染物產(chǎn)污強度、污染物綜合去除率等見下表1，對農(nóng)村生活污水進行處理的行政村比例來自《中國城鄉(xiāng)建設(shè)統(tǒng)計年鑒》。

表1 山東省臨沂市農(nóng)村生活污水污染物產(chǎn)污強度及污染物綜合去除率

3.2 種植業(yè)源污染負荷量估算種植業(yè)污染物排放量參考“手冊”核算方法，采用農(nóng)作物播種面積與污染物流失系數(shù)進行核算。山東省臨沂市種植業(yè)氮磷排放系數(shù)見表2。

表2 山東省臨沂市種植業(yè)氮磷排放系數(shù)

3.3 畜禽養(yǎng)殖業(yè)源污染負荷量估算畜禽養(yǎng)殖業(yè)污染物(COD、TN、TP)產(chǎn)生量和排放量均采用《手冊》中產(chǎn)排污核算法計算，畜禽養(yǎng)殖業(yè)污染物排放量等于各類畜禽養(yǎng)殖的污染物排放量之和，其中，第i類畜禽養(yǎng)殖污染物排放量，等于第i類畜禽養(yǎng)殖量與污染物排放系數(shù)的乘積。研究區(qū)畜禽面源污染排放系數(shù)見表3。

表3 山東省臨沂市畜禽面源污染排污系數(shù)

4 結(jié)果與分析

4.1 農(nóng)業(yè)面源污染主要影響因素分析圖2為研究區(qū)2000—2019年農(nóng)業(yè)及畜禽養(yǎng)殖業(yè)各要素時間變化圖，可看出，臨沂市農(nóng)藥使用量和農(nóng)用化肥施用總量年際變化較為一致，呈先增后減的趨勢，在2009年左右發(fā)生拐點變化。臨沂市2019年農(nóng)藥使用量為11254 t，比2000年減少了11.98%；氮肥、磷肥、鉀肥、復(fù)合肥折純后施用量為90304、24310、33012、170641 t，比2000年減少了12.26%。農(nóng)作物播種面積在2013年發(fā)生拐點，呈先增后減的趨勢，2019年農(nóng)作物播種面積比2000年和2009年分別下降了7.35%和9.06%。畜禽養(yǎng)殖按照《畜禽養(yǎng)殖污染物排放標準》(GB18596-2001)，將不同類型畜禽養(yǎng)殖量換算成豬的養(yǎng)殖量，研究區(qū)當量豬數(shù)量呈明顯增加趨勢，2019年比2000年增加了1.89倍，隨人們生活水平的提高，畜禽養(yǎng)殖業(yè)迅速發(fā)展，由此帶動了畜禽養(yǎng)殖量的增加。以往研究也表明，淮河流域人口集聚，畜禽養(yǎng)殖發(fā)展迅速，是造成當?shù)剞r(nóng)業(yè)面源污染的主要原因[21]。

從2019年研究區(qū)不同行政單元農(nóng)業(yè)及養(yǎng)殖業(yè)量占臨沂市的比例來看(圖3)，蘭陵縣、郯城縣和莒南縣在化肥施用和農(nóng)作物播種面積上占比最高，化肥施用量占比分別為18.30%、13.55%、11.54%；農(nóng)作物播種面積占比分別為17.29%、11.12%、11.04%。蘭陵縣、沂南區(qū)、平邑縣在農(nóng)藥使用上占比最高，分別為17.86%、12.49%、11.99%。沂南縣、沂水縣和費縣在畜禽養(yǎng)殖數(shù)量上占比最高，分別為20.03%、13.60%、12.03%。這也與后面不同類型農(nóng)業(yè)面源污染物排放量相一致，蘭陵縣TN和TP污染排放總量中來自種植業(yè)源的占比最高；沂南縣和沂水縣TN、TP、COD污染排放總量中來自畜禽養(yǎng)殖源的占比最高。

圖3 2019年研究區(qū)各行政單元農(nóng)業(yè)及養(yǎng)殖業(yè)量占臨沂市的比例Fig.3 Spatial distribution of agriculture and aquaculture in each administrative unit of the study area

4.2 農(nóng)業(yè)面源污染時間變化分析從圖4可看出，2019年TN和TP總排放量分別為154585.30 t/a，98363.34 t/a，比2000年增加了約9.00%，比2009年減少了20.60%。在TN和TP污染物排放量變化中，來自畜禽養(yǎng)殖業(yè)源污染排放量最大，且在2014年發(fā)生拐點變化，呈先增加后減少的趨勢；其次是來自農(nóng)村生活源污染的排放量，呈明顯減少趨勢；最后是來自種植業(yè)源污染的排放量，呈輕微減少趨勢。2019年COD總排放量為221414.35 t/a，比2000年和2009年分別增加了37.28%和7.21%。在COD排放量變化中，來自畜禽養(yǎng)殖業(yè)源的排放量最大且逐漸增加，其次是來自農(nóng)村生活源的排放量逐漸減少，COD總排放量呈先增加后平穩(wěn)變化的趨勢。整體來看，研究區(qū)的重點污染源為畜禽養(yǎng)殖業(yè)源。

圖4 研究區(qū)2000—2019年農(nóng)業(yè)面源污染物TN、TP、COD排放量變化Fig.4 Changes of TN，TP and COD emissions of agricultural non-point source pollutants in the study area from 2000 to 2019

由圖5可看出，COD是臨沂市農(nóng)業(yè)面源污染的主要污染物，其次是TN、最后是TP。2000—2019年TN、TP、COD排放貢獻率均值分別為35.55%、22.55%、41.90%。以2014年為拐點，2014年以前，TN、TP、COD排放貢獻率變化不大，2000—2014年TN、TP、COD排放貢獻率均值分別為36.59%、23.28%、40.12%。2014年以后，COD排放貢獻率逐漸增大，TN和TP排放貢獻率減小，2014—2019年TN、TP、COD排放貢獻率均值分別為34.50%、21.81%、43.68%。這可能受區(qū)域產(chǎn)業(yè)發(fā)展影響較大，近年來畜禽養(yǎng)殖的快速發(fā)展，使得研究區(qū)的COD排放明顯增加。因此，研究區(qū)應(yīng)重點關(guān)注污染物COD的減排和治理。已有研究表明，畜禽養(yǎng)殖業(yè)目前已成為我國環(huán)境污染的主要污染源之一，據(jù)統(tǒng)計，2010年規(guī)模戶養(yǎng)殖場COD排放量占農(nóng)業(yè)源總排放量的比例為95%，占全國COD排放總量的45%[22]。

圖5 研究區(qū)2000—2019年農(nóng)業(yè)面源污染物TN、TP、COD排放負荷貢獻率分布圖Fig.5 Distribution of contribution rate of agricultural non-point source pollutants TN，TP and COD emission load in the study area from 2000 to 2019

從表4可看出，研究區(qū)污染負荷排放貢獻率為畜禽養(yǎng)殖源>農(nóng)村生活源>種植業(yè)源，2019年臨沂市畜禽養(yǎng)殖業(yè)污染對TN、TP、COD的排放貢獻率分別為97.63%、99.78%、81.45%，農(nóng)村生活源為1.85%、0.20%、18.55%，種植業(yè)源TN為0.52%、TP為0.02%。臨沂市農(nóng)業(yè)面源污染主要來自畜禽養(yǎng)殖業(yè)源，這也與以往的研究較為一致[23]。根據(jù)山東省農(nóng)業(yè)畜禽養(yǎng)殖調(diào)查和基礎(chǔ)數(shù)據(jù)分析，來自畜禽養(yǎng)殖的排泄物對環(huán)境影響較大，目前山東省畜禽糞便無害化處理率不足40%，大量未經(jīng)處理的畜禽糞便直接排入農(nóng)田，對區(qū)域環(huán)境造成嚴重污染[24]。

表4 2019年臨沂市農(nóng)業(yè)面源污染物排放量及貢獻率

受區(qū)域農(nóng)業(yè)產(chǎn)業(yè)發(fā)展影響，研究區(qū)不同行政單元農(nóng)業(yè)面源污染主要污染源貢獻率有一定的差異(圖6)。在2019年TN和TP排放中，來自畜禽養(yǎng)殖業(yè)的排放對沂南縣和沂水縣的貢獻率最大，分別為23.25%和19.47%，來自種植業(yè)源的排放對蘭陵縣和郯城縣的貢獻率最大，分別為17.30%和11.13%。這也與前面農(nóng)業(yè)面源污染主要影響因素分析較為一致，化肥施用、農(nóng)作物播種面積、畜禽養(yǎng)殖量較大的區(qū)縣，產(chǎn)生的面源污染排放量也較大。

圖6 2019年研究區(qū)不同行政單元農(nóng)業(yè)面源污染排放貢獻分布圖Fig.6 Distribution map of agricultural no-point source pollution emission contribution of different administrative units in the study area in 2019

4.3 農(nóng)業(yè)面源污染空間變化分析從臨沂市不同行政單元農(nóng)業(yè)面源污染排放總量空間分布來看，各行政區(qū)污染物排放總量差異明顯。圖7為研究區(qū)TN、TP、COD排放強度分布圖，2019年研究區(qū)TN、TP、COD排放強度均值分別為8886.28 kg/km2、5663.59 kg/km2、12650.23 kg/km2。TN和TP排放強度比2000年減少了2.57%和2.47%，比2009年減少了20.71%和21.05%。COD排放強度比2000年增加了28.26%，比2009增加了6.13%。2000—2019年間TN和TP污染物排放強度呈先增加后減小的趨勢，COD排放強度呈逐漸增加的趨勢。

圖7 2000、2009、2019年研究區(qū)農(nóng)業(yè)面源污染TN、TP、COD排放強度空間分布Fig.7 Spatial distribution of TN，TP and COD emission intensity of agricultural non-point source pollution in the study area in 2000，2009 and 2019

綜合2000—2019年單位面積污染物排放強度GIS數(shù)據(jù)的自然間斷點[22]，將區(qū)域單位面積排放強度分為3個等級，重度污染區(qū)(TN>12000 kg/km2、TP>7000 kg/km2、COD>15000 kg/km2)、中度污染區(qū)(7000 kg/km2

從空間分布的時間變化來看，TN排放強度重度、中度、輕度的面積占臨沂市總面積的比例2000年分別為10.78%、49.19%、10.01%；2009年分別為50.15%、26.85%、22.98%；2019年分別為30.20%、22.14%、47.65%。TP排放強度重度、中度、輕度的面積占臨沂市總面積的比例2000年分別為20.83%、32.15%、47.00%；2009年分別為50.15%、26.85%、22.98%；2019年分別為30.20%、11.89%、57.90%?？煽闯?000—2019年間TN和TP排放強度重度區(qū)面積先增后減，中度區(qū)逐漸減少，輕度區(qū)先減后增。COD排放強度重度、中度、輕度的面積占臨沂市總面積的比例2000年分別為0.00%、20.69%、79.30%；2009年分別為15.96%、52.57%、31.45%；2019年分別為25.67%、44.48%、29.84%。COD排放強度重度區(qū)面積逐漸增加，中度區(qū)緩慢增加，輕度區(qū)逐漸減少。

整體來看，TN和TP排放強度在2009年最大，即重度污染區(qū)覆蓋面積最大，占總面積的比例達50.15%。隨著化肥減量、農(nóng)業(yè)綠色發(fā)展行動等的實施，面源污染治理取得較好效果，2019年TN和TP重度污染區(qū)覆蓋面積明顯減小，部分區(qū)縣降級為中度或輕度污染區(qū)。而COD排放強度重度區(qū)面積逐漸增加，是研究區(qū)重要的污染物。

從圖8可看出，2019年臨沂市TN、TP、COD的平均排放強度為8886.28 kg/km2、5663.59 kg/km2、12650.23 kg/km2。在全市12個區(qū)縣中，其中沂南縣、郯城縣、沂水縣和臨沭縣4個縣TN的排放強度超過臨沂市平均排放強度的百分比分別為127.58%、12.18%、37.69%、99.14%；TP的排放強度超過臨沂市平均排放強度的百分比分別為134.60%、12.86%、38.89%、106.45%。沂南縣、郯城縣、蘭陵縣、費縣、臨沭縣5個縣COD的排放強度超過臨沂市平均排放強度的百分比分別為90.58%、10.90%、5.42%、25.83%、43.52%。這些縣多位于畜禽養(yǎng)殖業(yè)和種植業(yè)較為發(fā)達區(qū)域，農(nóng)業(yè)面源污染排放量較大。因此，針對不同地區(qū)畜禽養(yǎng)殖污染的嚴重程度，因地制宜，實行差別化的污染治理至關(guān)重要，對畜禽養(yǎng)殖重度污染區(qū)，一方面要著力控制大型集約化畜禽養(yǎng)殖業(yè)排泄物的排放量，優(yōu)化畜禽養(yǎng)殖技術(shù)，加強生態(tài)養(yǎng)殖和清潔養(yǎng)殖治理技術(shù)；另一方面對散養(yǎng)畜禽養(yǎng)殖業(yè)應(yīng)促進畜禽養(yǎng)殖廢棄物的無害化處理和資源化利用。同時在空間布局上合理規(guī)劃畜牧養(yǎng)殖規(guī)模，防止由于養(yǎng)殖過于集中造成單體區(qū)域的環(huán)境污染。針對COD的處理，具體可采用固液分離處理技術(shù)降低畜禽排污水COD，包括重力沉降技術(shù)、蒸發(fā)技術(shù)、絮凝分離技術(shù)、篩分技術(shù)、壓濾技術(shù)及沉淀離心技術(shù)，其中絮凝分離技術(shù)方法簡單、成本低廉，應(yīng)用較為廣泛[25]。

圖8 2019年研究區(qū)農(nóng)業(yè)面源污染TN、TP、COD排放強度對比Fig.8 Comparison of TN，TP and COD emission intensity of agricultural non-point source pollution in the study area in 2019

5 結(jié)論

本研究基于排污系數(shù)法對臨沂市不同行政單元長時間序列的農(nóng)業(yè)種植、畜禽養(yǎng)殖、農(nóng)村生活3類污染源TN、TP、COD排放量進行了定量估算，系統(tǒng)分析了不同類型農(nóng)業(yè)面源污染排放時空特征及來源貢獻，得到以下結(jié)論：

(1)從時間尺度上看，2000—2019年間研究區(qū)農(nóng)業(yè)面源污染物TN和TP排放量呈“先增后減”，COD排放量呈“先增加后平穩(wěn)”的變化趨勢，其中2019年TN和TP總排放量比2009年減少了20.60%，COD總排放量比2009年增加了7.21%。

(2)從污染源和污染物結(jié)構(gòu)來看，臨沂市污染源負荷排放貢獻率排序為畜禽養(yǎng)殖源>農(nóng)村生活源>種植業(yè)源，畜禽養(yǎng)殖業(yè)污染對TN、TP、COD的排放貢獻率均高于80%；污染物排放量為COD>TN>TP，對應(yīng)的排放貢獻率均值分別為41.90%、35.55%、22.55%。因此，研究區(qū)應(yīng)重點加強畜禽養(yǎng)殖污染物的治理，尤其是COD的減排治理。

(3)從空間分布上看，TN和TP排放強度在2009年最大，即重度污染區(qū)覆蓋面積最大，占總面積的比例達50.15%，隨著化肥減量、農(nóng)業(yè)綠色發(fā)展行動等的實施，2019年TN和TP重度污染區(qū)覆蓋面積明顯減小，降為30.20%。而COD排放強度重度區(qū)面積逐漸增加，由2000年的15.96%增加到25.67%，是研究區(qū)的首要污染物。

(4)從污染物排放強度劃分結(jié)果來看，沂南縣、臨沭縣、沂水縣為TN和TP重度污染區(qū)，其中TN排放強度超過臨沂市平均排放強度的百分比對應(yīng)為127.58%、99.14%、37.69%，TP對應(yīng)為134.60%、106.45%、38.89%。沂南縣、臨沭縣、費縣為COD重度污染區(qū)，COD的排放強度超過臨沂市平均排放強度的百分比對應(yīng)為90.58%、43.52%、25.83%，以上幾個縣應(yīng)作為臨沂市農(nóng)業(yè)面源污染的優(yōu)先控制區(qū)進行重點治理。