中國省域化肥面源污染時(shí)空格局演變與分組預(yù)測*

雷俊華，蘇時(shí)鵬，余文夢，孫小霞**

(1.福建農(nóng)林大學(xué)農(nóng)村區(qū)域競爭力研究中心 福州 350002；2.中國人民大學(xué)環(huán)境學(xué)院 北京 100872)

化肥既是重要的農(nóng)業(yè)生產(chǎn)要素，又是導(dǎo)致農(nóng)業(yè)面源污染的主要因素之一。中國化肥施用強(qiáng)度較大，利用效率總體較低[1]，對糧食增產(chǎn)的邊際貢獻(xiàn)率不斷下降[2]，單位耕地面積化肥施用量遠(yuǎn)超國際安全使用水平[3]，嚴(yán)重威脅到中國農(nóng)業(yè)經(jīng)濟(jì)的可持續(xù)發(fā)展[4-5]。過量的化肥投入不僅增加了農(nóng)業(yè)生產(chǎn)的經(jīng)濟(jì)成本，還帶來了嚴(yán)重的生態(tài)環(huán)境負(fù)擔(dān)[6]。探討農(nóng)業(yè)發(fā)展與化肥面源污染之間時(shí)空演化關(guān)系，可以更好地推動(dòng)化肥零增長行動(dòng)，促進(jìn)中國農(nóng)業(yè)高質(zhì)量綠色發(fā)展。

環(huán)境庫茲涅茨曲線(Environment Kuznets Curve，EKC)是經(jīng)濟(jì)增長與環(huán)境關(guān)系實(shí)證研究的重要理論。Grossman 等[7]提出環(huán)境與經(jīng)濟(jì)之間可能存在類似庫茲涅茨曲線的“倒 U 型”關(guān)系。Antle 等[8]和McConnell[9]較早探討了農(nóng)業(yè)面源污染與經(jīng)濟(jì)發(fā)展間的關(guān)系，認(rèn)為農(nóng)業(yè)面源污染也有EKC 的變化特征。相關(guān)研究表明，農(nóng)業(yè)面源污染EKC 可能呈現(xiàn)為“倒U 型”[10-12]、“U 型”[13-15]、“N 型”或“倒N 型”[16-18]以及線性增長[19-20]等多種演變趨勢。然而，以上研究存在兩方面的不足：1)用化肥投入量來表征農(nóng)業(yè)面源污染，忽視了農(nóng)業(yè)面源污染是化肥施用后，流失的養(yǎng)分進(jìn)入水體才構(gòu)成污染的事實(shí)[21]，從而無法反映化肥面源污染實(shí)際狀況。2)基于同質(zhì)性假設(shè)，只關(guān)注國家或地區(qū)的整體狀況，忽略了區(qū)域間異質(zhì)性的影響。各地區(qū)經(jīng)濟(jì)發(fā)展水平、地理環(huán)境、資源稟賦、勞動(dòng)力水平等差異都可能對當(dāng)?shù)剞r(nóng)業(yè)面源污染帶來不同影響[22]，使得基于同質(zhì)性假設(shè)的EKC 擬合結(jié)果不符合區(qū)域特質(zhì)。部分研究按照東、中、西部分組[23]、省份單獨(dú)檢驗(yàn)[24]、數(shù)理方法選取分類變量[22]等方式進(jìn)行分組EKC 擬合，試圖克服同質(zhì)性問題。然而，東、中、西部的劃分方法不能合理反映化肥面源污染的地區(qū)特質(zhì)[25]，而各省份分別檢驗(yàn)的方法則忽視了面源污染在區(qū)域間的流動(dòng)性和分散性[26]，以及省份間產(chǎn)業(yè)流動(dòng)和勞動(dòng)力轉(zhuǎn)移等因素導(dǎo)致的空間上的相關(guān)性，進(jìn)而影響估計(jì)結(jié)果的準(zhǔn)確性。數(shù)理方法分組雖有利于運(yùn)用數(shù)據(jù)的內(nèi)生特性優(yōu)勢，但指標(biāo)選取往往具有主觀性且假定一個(gè)區(qū)域的環(huán)境污染與相鄰區(qū)域在地理空間上不存在空間相關(guān)性[18]，從而可能使擬合結(jié)果與實(shí)際情況產(chǎn)生偏差。

本文利用全國1997—2018年的省際面板數(shù)據(jù)，首先對各省化肥面源污染排放強(qiáng)度進(jìn)行核算，更加準(zhǔn)確地反映各省歷年化肥面源污染狀況。在此基礎(chǔ)上，利用ArcGIS 空間統(tǒng)計(jì)模型對其進(jìn)行空間自相關(guān)分析和熱點(diǎn)分析，綜合時(shí)空格局演變分析結(jié)果對各省份進(jìn)行分組。此種分組方式既考慮了不同區(qū)域化肥面源污染的實(shí)際情況，又考慮了區(qū)域的空間相關(guān)性，更加符合各區(qū)域的實(shí)際?；诜纸M，擬合各組化肥面源污染排放強(qiáng)度與人均農(nóng)業(yè)生產(chǎn)總值之間EKC關(guān)系的時(shí)間路徑，考察空間異質(zhì)性條件下各組時(shí)間路徑的差異性，最后對時(shí)間路徑的未來演變趨勢進(jìn)行預(yù)測，試圖為化肥減量化和農(nóng)業(yè)面源污染防治工作提供一定參考。

1 研究方法與數(shù)據(jù)來源

1.1 研究方法

1.1.1 面源污染排放強(qiáng)度核算

當(dāng)前相關(guān)研究多采用單元調(diào)查評估方法[27-28]和清單分析法[29]進(jìn)行大尺度區(qū)域的化肥面源污染核算，然而此類方法沒有考慮土地類型、自然條件(如氣候)等因素的差異對化肥流失產(chǎn)生的異質(zhì)性影響，利用化肥施用強(qiáng)度和施用比例確定的化肥流失率也可能導(dǎo)致估計(jì)結(jié)果不夠準(zhǔn)確[30]。因此，本文在單元調(diào)查評估法和清單分析法的基礎(chǔ)上，結(jié)合流失系數(shù)法進(jìn)行化肥面源污染的核算。參考《第一次全國污染源普查—農(nóng)業(yè)污染源—肥料流失系數(shù)手冊》的6 大分區(qū)，綜合考慮各省氣候、地形、土壤等主要影響因素，將全國化肥面源污染源分為6 個(gè)區(qū)域類別，各類別的流失系數(shù)主要運(yùn)用“一污普”后發(fā)布的《農(nóng)業(yè)污染源肥料流失系數(shù)手冊》中各區(qū)域各模式流失系數(shù)的均值，并通過廣泛的文獻(xiàn)調(diào)研對各類區(qū)域的流失系數(shù)進(jìn)行確認(rèn)和校正，以使流失系數(shù)更符合各區(qū)域的實(shí)際情況，保證流失系數(shù)的科學(xué)性和準(zhǔn)確性，流失系數(shù)見表1。

本文核算的化肥面源污染單元為磷肥、氮肥與復(fù)合肥。污染物指標(biāo)采用總氮(TN)、硝氮(NO3--N )、氨氮(NH4+-N)、總磷(TP)和可溶性總磷(DTP)[31]。運(yùn)用化肥流失系數(shù)法，核算氮的總排放量和磷的總排放量，進(jìn)而核算化肥面源污染排放總量及強(qiáng)度。計(jì)算公式為：

表1 中國31 省(市、自治區(qū))化肥流失情況分區(qū)表 Table1 Fertilizer loss rate in different regions of China

式中：E為化肥面源污染總排放量；Eij為污染單元i產(chǎn)生并流失進(jìn)入水體的第j種污染物的量；Ci是污染單元i的指標(biāo)統(tǒng)計(jì)數(shù)(氮肥的統(tǒng)計(jì)數(shù)為N 的量，磷肥統(tǒng)計(jì)數(shù)為P2O5的量，復(fù)合肥按32%折算P2O5的量，按40%折算TN 的量[32])；ηij為污染單元i中第j種污染物的流失系數(shù)；A為各省農(nóng)作物播種面積[22，27]；I(kg·hm-2)為單位面積的污染排放量，即化肥面源污染排放強(qiáng)度。

1.1.2 空間統(tǒng)計(jì)模型

1)全局空間自相關(guān)?；拭嬖次廴九欧艔?qiáng)度在 全局空間上的相關(guān)程度用Global Moran’sI測度，其含義是空間相鄰省份的化肥面源污染排放強(qiáng)度在整體區(qū)域內(nèi)空間相關(guān)的總體趨勢[33]。在給定的顯著水平下，若Moran’sI值為負(fù)，則污染排放強(qiáng)度呈顯著的空間離散趨勢；若Moran’sI值為正，則污染排放強(qiáng)度呈顯著的空間集聚趨勢[34]。Global Moran’sI的計(jì)算公式為：

式中：n為省份個(gè)數(shù)，Wij為空間權(quán)重矩陣(兩省相鄰為1，不相鄰為0)，Xi和Xj為i省和j省的觀測值，為觀測值的平均值，S2為觀測值的方差。

2)局部空間自相關(guān)。局部自相關(guān)特征用熱點(diǎn)分析的Gi*指數(shù)進(jìn)行測算，用于識別區(qū)域內(nèi)的冷點(diǎn)區(qū)和熱點(diǎn)區(qū)。若Gi*指數(shù)顯著為正，則i省周圍污染排放強(qiáng)度相對較高，屬于熱點(diǎn)區(qū)；若Gi*指數(shù)顯著為負(fù)，則i省周圍污染排放強(qiáng)度相對較低，屬于冷點(diǎn)區(qū)；若Gi*指數(shù)不顯著，則i省不是熱(冷)點(diǎn)區(qū)。計(jì)算公式為：

式中：Wij(d)是根據(jù)距離規(guī)則定義的空間權(quán)重。

1.1.3 EKC 模型設(shè)定

EKC 模型通常先設(shè)為三次曲線，若不顯著，則再設(shè)為二次曲線，若二次曲線仍不顯著，則擬合成線性模型[35]。因此，EKC 曲線可能表現(xiàn)為“N 型” “倒N 型” “U 型” “倒U 型”或者直線型。為減少異方差，本文對各指標(biāo)進(jìn)行對數(shù)化處理，設(shè)定化肥面源污染排放的EKC 模型為：

式中：I表示化肥面源污染排放強(qiáng)度；G表示人均農(nóng)業(yè)產(chǎn)值，作為解釋變量，用來表征農(nóng)業(yè)經(jīng)濟(jì)發(fā)展水平；G(元·人–1)=農(nóng)業(yè)總產(chǎn)值(狹義)/鄉(xiāng)村總?cè)丝跀?shù)(狹義農(nóng)業(yè)總產(chǎn)值指種植業(yè)和其他農(nóng)業(yè)的主產(chǎn)品和副產(chǎn)品的總產(chǎn)值，主要包括糧、棉、油料、麻類、煙葉、蔬菜、藥材、瓜類和其他農(nóng)作物的種植，以及茶園、桑園、果園的生產(chǎn)經(jīng)營。狹義上的農(nóng)業(yè)是化肥面源污染的主要來源，因而此處用狹義的農(nóng)業(yè)總產(chǎn)值計(jì)算)。其他指標(biāo)是參考現(xiàn)有相關(guān)研究選取的控制變量：P表示作物種植結(jié)構(gòu)，用糧食播種面積與農(nóng)作物總播種面積的比值表征，作物種植結(jié)構(gòu)差異可能致使化肥施用量的差異；M表示農(nóng)業(yè)機(jī)械投入強(qiáng)度(kg·hm-2)，用單位農(nóng)作物播種面積的農(nóng)業(yè)機(jī)械投入強(qiáng)度表征，機(jī)械化水平提高可能會(huì)對化肥施用量產(chǎn)生正負(fù)兩面的影響；S表示產(chǎn)業(yè)結(jié)構(gòu)，用農(nóng)業(yè)產(chǎn)值占GDP 的比重表征

農(nóng)業(yè)產(chǎn)值所占比重大小影響到化肥施用量的大小，地區(qū)產(chǎn)業(yè)結(jié)構(gòu)會(huì)影響當(dāng)?shù)氐沫h(huán)境質(zhì)量；R表示農(nóng)村居民收入結(jié)構(gòu)，農(nóng)村居民收入結(jié)構(gòu)用農(nóng)村居民人均工資性收入與農(nóng)村居民人均純收入的比值表征，農(nóng)村居民人均工資性收入比重的高低可能影響化肥施用量的大??；N表示農(nóng)村居民人均純收入，農(nóng)村居民收入的高低會(huì)影響化肥的投入；U表示城鄉(xiāng)居民收入差距，用城鎮(zhèn)居民人均可支配收入與農(nóng)村居民人均純收入之比表征，收入差距可能是環(huán)境退化的重要因素，均衡的收入分配可能減少污染[36]；A表示農(nóng)業(yè)產(chǎn)業(yè)結(jié)構(gòu)，用狹義農(nóng)業(yè)總產(chǎn)值占農(nóng)林牧漁業(yè)總產(chǎn)值的比重表征，以反映種養(yǎng)結(jié)構(gòu)調(diào)整對化肥面源污染的影響；μ為隨機(jī)誤差項(xiàng)；i表示省份，t表示年份；回歸模型的估計(jì)系數(shù)βn(n=1～10)表示彈性系數(shù)。

1.2 數(shù)據(jù)來源

本文在我國化肥面源污染排放強(qiáng)度核算的基礎(chǔ)上，對其進(jìn)行時(shí)空格局演變分析，并以時(shí)空特征為分類條件，對空間異質(zhì)性條件下化肥面源污染EKC的時(shí)間路徑進(jìn)行擬合與預(yù)測。研究數(shù)據(jù)為1997—2018年31 省(市、自治區(qū)，港澳臺數(shù)據(jù)缺失，以下均稱省)共22年的面板數(shù)據(jù)，來源于《中國農(nóng)村統(tǒng)計(jì)年鑒》(1998—2019年)和《中國統(tǒng)計(jì)年鑒》(1998—2019年)，并用插值法補(bǔ)齊少量年份因年鑒統(tǒng)計(jì)口徑變化造成的缺失值。為實(shí)現(xiàn)數(shù)據(jù)的可比性，農(nóng)業(yè)總產(chǎn)值通過農(nóng)業(yè)產(chǎn)值指數(shù)折算為1997年不變價(jià)，農(nóng)村居民人均純收入和城鎮(zhèn)居民人均可支配收入運(yùn)用消費(fèi)者價(jià)格指數(shù)折算為1997年不變價(jià)。

2 化肥面源污染排放強(qiáng)度核算及其時(shí)空格局演變

2.1 化肥面源污染排放強(qiáng)度核算

根據(jù)公式(1)和(2)對我國1997—2018年各省化肥面源污染排放強(qiáng)度進(jìn)行核算，各分區(qū)核算結(jié)果如圖1所示。我國化肥面源污染排放強(qiáng)度省際差異較大。北京、江蘇、上海、浙江、湖北等經(jīng)濟(jì)發(fā)達(dá)省份排放強(qiáng)度較大，22年間大多數(shù)年份在5 kg·hm-2以上，超過其他省份的平均強(qiáng)度。青海、西藏、黑龍江等省份排放強(qiáng)度較小。各省22年間化肥面源污染排放強(qiáng)度的走勢各異，表現(xiàn)出不同升降趨勢，表明我國化肥面源污染排放強(qiáng)度時(shí)空差異較大。各區(qū)域內(nèi)各省份排放強(qiáng)度差異較小，因此高(低)強(qiáng)度區(qū)域可能呈現(xiàn)一定的集聚特征。所以有必要進(jìn)一步分析化肥面源污染排放強(qiáng)度的總體和局部時(shí)空格局演變。

2.2 省際化肥面源污染排放強(qiáng)度時(shí)空格局演變特征

2.2.1 總體時(shí)空格局演化特征

圖1 1997—2018年各分區(qū)化肥面源污染排放強(qiáng)度Fig.1 Emission intensities of non-point source pollution of chemical fertilizers in different regions from 1997 to 2018

表2 1997—2018年各省(市、自治區(qū))化肥面源污染排放強(qiáng)度Global Moran’sI值Table2 Global Moran’s I values of chemical fertilizer non-point source pollution emission intensity in various provinces (cities，autonomous regions) from 1997 to 2018

運(yùn)用ArcGIS 10.2 軟件計(jì)算1997—2018年各省域化肥面源污染排放強(qiáng)度的Global Moran’sI指數(shù)，結(jié)果如表2所示。1997—2018年間，我國省域化肥面源污染排放強(qiáng)度的Moran’sI指數(shù)均通過了顯著性檢驗(yàn)，且均大于0.307，Z得分均大于4.352，說明我國化肥污染排放強(qiáng)度呈正向空間自相關(guān)性，化肥面源污染在區(qū)域上一直處于較高的集聚狀態(tài)，區(qū)域間存在較強(qiáng)的相關(guān)性和依賴性。Moran’sI指數(shù)22年間呈現(xiàn)小幅波動(dòng)變化，1997—2006年總體較為平穩(wěn)，2006年開始呈現(xiàn)出波動(dòng)下降的趨勢，說明我國化肥面源污染排放強(qiáng)度空間格局總體較穩(wěn)定，空間自相關(guān)性呈現(xiàn)出先平穩(wěn)后下降的特征?？臻g自相關(guān)性降低表明區(qū)域與周邊地區(qū)變化趨勢的一致性相對減弱，應(yīng)當(dāng)注意區(qū)域間的異質(zhì)性特征和空間溢出效應(yīng)的影響。

2.2.2 局部時(shí)空格局演化特征

利用ArcGIS 10.2 軟件的空間統(tǒng)計(jì)工具計(jì)算化肥面源污染排放強(qiáng)度省域空間的Gi*指數(shù)。依據(jù)自然斷點(diǎn)法將Gi*值分為7 類：排放強(qiáng)度冷點(diǎn)(99%置信)、排放強(qiáng)度冷點(diǎn)(95%置信)、排放強(qiáng)度冷點(diǎn)(90%置信)、排放強(qiáng)度熱點(diǎn)(99%置信)、排放強(qiáng)度熱點(diǎn)(95%置信)、排放強(qiáng)度熱點(diǎn)(90%置信)和不顯著。部分年份空間格局如表3所示(限于篇幅，本文只展現(xiàn)時(shí)空格局變化較大的1997年、2003年、2007年、2009年、2014年和2018年的空間格局分布情況)。分析可知：

1)化肥面源污染排放強(qiáng)度空間格局總體較穩(wěn)定。熱點(diǎn)區(qū)主要集中分布在華中和華南地區(qū)，其中長江中下游地區(qū)尤其顯著，常年處于99%的置信熱點(diǎn)。冷點(diǎn)區(qū)主要集中在西部地區(qū)的青海、甘肅和西藏，東北的黑龍江也在2003年后長期呈現(xiàn)為冷點(diǎn)狀態(tài)。化肥面源污染熱(冷)點(diǎn)不顯著的區(qū)域主要集中在我國西南部以及東北和西北的部分省份。

2)化肥面源污染排放強(qiáng)度空間格局在各階段呈現(xiàn)出局部變動(dòng)，黃淮海地區(qū)呈現(xiàn)出熱點(diǎn)程度減弱的變化。1997—2002年，全國化肥面源污染空間格局較為穩(wěn)定，各省份熱(冷)點(diǎn)的屬性未發(fā)生大的變動(dòng)，只有少量省份在置信度上表現(xiàn)出一定的波動(dòng)。2003年和2004年，黑龍江在90%置信冷點(diǎn)和不顯著間波動(dòng)，2005年后則持續(xù)呈現(xiàn)為90%的置信冷點(diǎn)。2007年全國空間格局發(fā)生較大變化，京津冀地區(qū)首次從熱點(diǎn)區(qū)消失，山東的熱點(diǎn)顯著度也出現(xiàn)下降，黃淮海部分地區(qū)迎來化肥面源污染的拐點(diǎn)。2009年廣東首次進(jìn)入熱點(diǎn)區(qū)的版圖，至今雖有個(gè)別年份呈現(xiàn)為不顯著，但總體呈現(xiàn)為熱點(diǎn)狀態(tài)。2014年后，部分熱點(diǎn)省份發(fā)生波動(dòng)，如2014年，廣東由90%置信熱點(diǎn)變?yōu)椴伙@著，江西由99%置信熱點(diǎn)變?yōu)?5%置信熱點(diǎn)，山東由95%置信熱點(diǎn)變?yōu)?0%置信熱點(diǎn)；2017年，山東由90%置信熱點(diǎn)進(jìn)而變?yōu)椴伙@著，河南由長期的99%置信熱點(diǎn)變?yōu)?5%置信熱點(diǎn)。從全國冷熱點(diǎn)時(shí)空格局的變化來看，近年來對于化肥面源污染的治理取得一定成效，部分區(qū)域的熱點(diǎn)程度有所下降，尤其是黃淮海地區(qū)的效果較為明顯。近年來，高強(qiáng)度的減量化措施使我國化肥施用量短期減少，但目前農(nóng)業(yè)生產(chǎn)方式對于化肥的依賴性以及政策執(zhí)行的波動(dòng)性還是會(huì)使局部地區(qū)化肥面源污染狀況產(chǎn)生波動(dòng)。

表3 1997—2018年各省(市、自治區(qū))化肥面源污染排放強(qiáng)度空間格局分布Table3 Spatial pattern of chemical fertilizer non-point source pollution emission intensity in different provinces (cities，autonomous regions) from 1997 to 2018

對31 個(gè)省22年熱(冷)點(diǎn)發(fā)生的頻次和顯著水平進(jìn)行統(tǒng)計(jì)，據(jù)此將31 個(gè)省份分為熱點(diǎn)區(qū)、冷點(diǎn)區(qū)、非熱(冷)點(diǎn)區(qū)3 類區(qū)域，結(jié)果如表4所示。3 類區(qū)域化肥面源污染狀況的差異體現(xiàn)出了我國化肥面源污染在空間上的異質(zhì)性。

表4 化肥面源污染排放強(qiáng)度熱點(diǎn)分析的分組結(jié)果Table4 Grouping results of hot spot analysis of chemical fertilizer non-point source pollution emission intensity

3 化肥面源污染EKC 分組模擬與預(yù)測

3.1 面板單位根檢驗(yàn)和協(xié)整檢驗(yàn)

面板數(shù)據(jù)兼具截面數(shù)據(jù)和時(shí)間序列數(shù)據(jù)的特點(diǎn)，為避免偽回歸，在構(gòu)建面板回歸模型前需用單位根方法檢驗(yàn)各個(gè)變量的平穩(wěn)性和協(xié)整關(guān)系。常見的面板數(shù)據(jù)單位根檢驗(yàn)有LLC、IPS 和費(fèi)雪式檢驗(yàn)等方法，本文運(yùn)用以上方法對各組分別進(jìn)行檢驗(yàn)，結(jié)果表明各種檢驗(yàn)結(jié)論一致，限于篇幅僅 列出LLC 檢驗(yàn)結(jié)果。LLC 檢驗(yàn)結(jié)果(表5)顯示，多數(shù)變量的原始序列存在單位根。因此對原始序列進(jìn)行一階差分，結(jié)果顯示各變量在一階差分時(shí)為平穩(wěn)，序列為一階單整序列I(1)，因此可以進(jìn)行協(xié)整檢驗(yàn)。

選取Kao 檢驗(yàn)方法對各組回歸模型的二次曲線和三次曲線進(jìn)行協(xié)整檢驗(yàn)，結(jié)果如表6所示。各組ADF 值均通過了顯著性檢驗(yàn)，表明組內(nèi)各個(gè)變量間的關(guān)系長期穩(wěn)定，可構(gòu)建面板模型。

表5 化肥面源污染排放強(qiáng)度不同熱點(diǎn)區(qū)域的面板單位根(LLC)檢驗(yàn)結(jié)果Table5 Results of panel unit root test (LLC test) of different hotspot regions of fertilizer non-point source pollution emission intensity

續(xù)表5

表6 化肥面源污染排放強(qiáng)度不同熱點(diǎn)區(qū)域的協(xié)整檢驗(yàn)結(jié)果Table6 Results of panel cointegration test of different hotspot regions of fertilizer non-point source pollution emission intensity

3.2 組間異方差、組內(nèi)自相關(guān)和組間同期相關(guān)的檢驗(yàn)

由于3 類區(qū)域的分組數(shù)據(jù)均為長面板數(shù)據(jù)，因此需要對各組數(shù)據(jù)進(jìn)行組間異方差、組內(nèi)自相關(guān)和組間同期相關(guān)的檢驗(yàn)，并根據(jù)結(jié)果確定合適的計(jì)量方法。檢驗(yàn)結(jié)果(表7)顯示3 類區(qū)域的二次和三次模型均存在組間異方差和組間同期相關(guān)；而組內(nèi)自相關(guān)檢驗(yàn)表明，熱點(diǎn)區(qū)和非熱(冷)點(diǎn)區(qū)的二次模型和三次模型均存在組內(nèi)一階自相關(guān)，冷點(diǎn)區(qū)的二次模型和三次模型則不存在組內(nèi)自相關(guān)。

3.3 化肥面源污染EKC 空間分組擬合結(jié)果

為解決組間異方差、組內(nèi)自相關(guān)和組間同期相關(guān)問題，熱點(diǎn)區(qū)和非熱(冷)點(diǎn)區(qū)均采用全面FGLS 方法進(jìn)行估計(jì)，該方法同時(shí)考慮了這3 個(gè)因素，能夠保證估計(jì)的效率。對于不存在組內(nèi)自相關(guān)的冷點(diǎn)區(qū)，采用LSDV 方法進(jìn)行估計(jì)，并使用組間異方差和組間同期相關(guān)的穩(wěn)健標(biāo)準(zhǔn)誤差來保證估計(jì)結(jié)果的穩(wěn)健性。表8報(bào)告了各組的回歸估計(jì)結(jié)果。

回歸模型結(jié)果顯示，熱點(diǎn)區(qū)的三次曲線模型在10%顯著水平未通過檢驗(yàn)，二次曲線模型在1%顯著水平上通過檢驗(yàn)，故采用二次曲線形式；冷點(diǎn)區(qū)的三次曲線形式模型在10%顯著水平通過檢驗(yàn)，故采用三次曲線形式；非熱(冷)點(diǎn)區(qū)的三次曲線在5%顯著水平通過檢驗(yàn)，故采用三次曲線形式。通過lnG的各次項(xiàng)的正負(fù)，可以判斷熱點(diǎn)區(qū)的EKC 形狀為“倒U 型”，冷點(diǎn)區(qū)和非熱(冷)點(diǎn)區(qū)為均“倒N 型”。

1)熱點(diǎn)區(qū)。本組包含上海、江蘇、浙江、河南、安徽、湖北等10 個(gè)省份，這些省份主要是我國經(jīng)濟(jì)發(fā)達(dá)地區(qū)或農(nóng)業(yè)大省，化肥面源污染強(qiáng)度普遍較大。本組EKC 的時(shí)間路徑呈現(xiàn)為“倒U 型”曲線，說明伴隨農(nóng)業(yè)發(fā)展，化肥面源污染表現(xiàn)為“升—降”趨勢。該曲線拐點(diǎn)對應(yīng)的人均農(nóng)業(yè)生產(chǎn)總值G 分別為5 889.933 元。從熱點(diǎn)區(qū)10 省的情況看，江蘇(2011年)、山東(2016年)、福建(2017年)、河南(2017年)、河北(2017年)5 個(gè)省份已經(jīng)跨過了“倒U 型”曲線的拐點(diǎn)，進(jìn)入了“隨著農(nóng)業(yè)經(jīng)濟(jì)增長，化肥面源污染減弱”的下降階段。從時(shí)間上看，多數(shù)省份在2017年左右跨過拐點(diǎn)，其原因除經(jīng)濟(jì)不斷增長促使化肥面源污染到達(dá)拐點(diǎn)外，國家糧食安全戰(zhàn)略下糧食主產(chǎn)區(qū)近年糧食播種面積比重的回升、農(nóng)村居民人均純收入的提高和農(nóng)村居民對農(nóng)業(yè)收入依賴的降低都可能促進(jìn)了拐點(diǎn)的到來。此外，2015年以來，國家頒布的一系列化肥減量化政策，如《到2020年化肥使用量零增長行動(dòng)方案》以及近年中央一號文件對農(nóng)業(yè)面源污染重視程度的不斷加強(qiáng)也可能促進(jìn)了拐點(diǎn)的到來。安徽、湖南、江西、上海和浙江則均尚未到達(dá)拐點(diǎn)，仍處于“隨著農(nóng)業(yè)經(jīng)濟(jì)增長，化肥面源污染增強(qiáng)”的階段。安徽、湖南和江西作為我國重要的產(chǎn)糧區(qū)，對于化肥的依賴程度較高，難以在短期擺脫對化肥的依賴。上海和浙江作為糧食主銷區(qū)，較低的種植業(yè)比重反而致使化肥施用的集約化程度更低，發(fā)展種養(yǎng)循環(huán)的動(dòng)力更小，從而未充分發(fā)揮有機(jī)肥代替化肥等方面的作用。熱點(diǎn)區(qū)農(nóng)業(yè)機(jī)械投入強(qiáng)度與化肥面源污染程度呈正相關(guān)，表明科技和機(jī)械水平雖然在提高農(nóng)業(yè)產(chǎn)值方面發(fā)揮了一定作用，但是對于提高化肥利用效率和控制化肥流失率的作用尚未明顯體現(xiàn)，應(yīng)該充分發(fā)揮農(nóng)業(yè)新技術(shù)的作用，如生態(tài)循環(huán)農(nóng)業(yè)、配方施肥和精準(zhǔn)施肥等。

表7 各組組間異方差、組內(nèi)自相關(guān)、組間同期相關(guān)檢驗(yàn)結(jié)果Table7 Results of test for groupwise heteroskedasticity，autocorrelation within panel and contemporaneous correlation in each regions

表8 化肥面源污染EKC 分組擬合結(jié)果Table8 EKC group fitting results of chemical fertilizer non-point source pollution

2)冷點(diǎn)區(qū)。本組包含甘肅、青海、西藏、黑龍江4 個(gè)省份，主要分布在我國西部和東北。本組EKC 的時(shí)間路徑為“倒N 型”曲線，說明伴隨農(nóng)業(yè)發(fā)展，化肥面源污染表現(xiàn)為“降—升—降”趨勢。該曲線的拐點(diǎn)對應(yīng)的人均農(nóng)業(yè)生產(chǎn)總值G分別為 851.500 元和4 433.746 元。冷點(diǎn)區(qū)4 省均在1997年以前就跨過了“倒N 型”曲線的第1 個(gè)拐點(diǎn)，已經(jīng)跨過第2 個(gè)拐點(diǎn)的省份為甘肅(2017年)和黑龍江(2006年)，表明這兩個(gè)省份目前處于曲線下降階段，即“隨著農(nóng)業(yè)經(jīng)濟(jì)增長，化肥面源污染減弱”。青海和西藏則還未跨過第2 個(gè)拐點(diǎn)，仍處于曲線的上升階段，且表現(xiàn)出短期內(nèi)不會(huì)達(dá)到拐點(diǎn)的趨勢。對于青海和西藏而言，其本身較低的農(nóng)業(yè)經(jīng)濟(jì)增長水平促使其距離拐點(diǎn)較遠(yuǎn)。黑龍江農(nóng)業(yè)產(chǎn)值高但化肥面源污染排放強(qiáng)度卻一直較低，其原因在于東北黑土資源肥沃、供肥能力強(qiáng)，且多年來耕地質(zhì)量普遍不斷提高[37]，化肥流失率低，從而使黑龍江較早跨過拐點(diǎn)進(jìn)入曲線下降階段。

3)非熱(冷)點(diǎn)區(qū)。本組主要包括廣西、云南、四川、吉林、遼寧、陜西等17 個(gè)省份。本組EKC 的時(shí)間路徑為“倒N 型”曲線，說明伴隨農(nóng)業(yè)發(fā)展，化肥面源污染表現(xiàn)為“降—升—降”趨勢。該曲線兩個(gè)拐點(diǎn)對應(yīng)的人均農(nóng)業(yè)生產(chǎn)總值G分別為1 490.698元和9 405.031 元。從非熱(冷)點(diǎn)區(qū)17 省的情況看，多數(shù)省份在1997年以前就跨過第1 個(gè)拐點(diǎn)，西南和華北部分省份則在2000—2010年間跨過第1 個(gè)拐點(diǎn)，進(jìn)入了“隨著農(nóng)業(yè)經(jīng)濟(jì)增長，化肥面源污染增強(qiáng)”的上升階段。然而，當(dāng)前非熱(冷)點(diǎn)區(qū)跨過第2 個(gè)拐點(diǎn)的省份只有海南(2016年)和新疆(2016年)。表明當(dāng)前非熱(冷)點(diǎn)區(qū)的多數(shù)省份均處于曲線的上升階段，即“隨著農(nóng)業(yè)經(jīng)濟(jì)增長，化肥面源污染增強(qiáng)”。按其曲線趨勢，非熱(冷)點(diǎn)區(qū)各省的拐點(diǎn)相對熱點(diǎn)區(qū)較晚，且由于空間格局上近年空間自相關(guān)有所減弱，因此應(yīng)當(dāng)注意非熱(冷)點(diǎn)區(qū)和熱點(diǎn)區(qū)之間空間相關(guān)性的變化，討論其相互間的影響。

3.4 化肥面源污染EKC 時(shí)間路徑分組預(yù)測

綜合來看，3 類區(qū)域的曲線形態(tài)和拐點(diǎn)差異較大，證明了空間異質(zhì)性分析的必要性。截至2018年，熱點(diǎn)區(qū)仍有5 個(gè)省份尚未跨過“倒U 型”曲線的拐點(diǎn)，處于曲線的上升階段；冷點(diǎn)區(qū)4 個(gè)省份中還有兩個(gè)未跨過“倒N 型”曲線的第2 個(gè)拐點(diǎn)；非熱(冷)點(diǎn)區(qū)仍有15 個(gè)省份未跨過“倒N 型”的第2 個(gè)拐點(diǎn)。因此當(dāng)前我國大部分省份正處于“隨著農(nóng)業(yè)經(jīng)濟(jì)增長，化肥面源污染增強(qiáng)”的曲線上升階段。在此情況下，有必要預(yù)測這些省份EKC 曲線的時(shí)間路徑，以根據(jù)其階段特征和未來趨勢采取相應(yīng)的措施。本文運(yùn)用趨勢延伸法，依據(jù)各省22年的人均農(nóng)業(yè)產(chǎn)值年均增長率來預(yù)估下一拐點(diǎn)的出現(xiàn)年份，結(jié)果如表9所示。對于熱點(diǎn)區(qū)省份，除上海外，其余省份可能在近年迎來拐點(diǎn)。對于上海而言，其農(nóng)業(yè)產(chǎn)業(yè)比重較低，應(yīng)當(dāng)進(jìn)一步加強(qiáng)產(chǎn)業(yè)結(jié)構(gòu)的調(diào)整，通過區(qū)域間產(chǎn)業(yè)轉(zhuǎn)移等方式來尋求經(jīng)濟(jì)效益和生態(tài)效益的平衡，獲取更高的效益。對于冷點(diǎn)區(qū)，青海和西藏距離拐點(diǎn)的時(shí)間較長，且兩省生態(tài)脆弱，應(yīng)進(jìn)一步調(diào)整農(nóng)業(yè)結(jié)構(gòu)，推廣農(nóng)業(yè)新技術(shù)的應(yīng)用，利用其畜牧業(yè)優(yōu)勢，發(fā)展種養(yǎng)循環(huán)、生態(tài)農(nóng)業(yè)，使兩省的拐點(diǎn)提前到來，進(jìn)入“倒N 型”曲線的下降階段。非熱(冷)點(diǎn)區(qū)中除遼寧、吉林、內(nèi)蒙古等省份外，其他省份距離“倒N 型”曲線的第2 個(gè)拐點(diǎn)距離較遠(yuǎn)，將較長期處于曲線上升階段，應(yīng)當(dāng)尤其引起重視。

4 討論與結(jié)論

4.1 討論

近年來，化肥減量化和化肥面源污染防治成為我國農(nóng)業(yè)工作的重點(diǎn)和研究的熱點(diǎn)。隨著農(nóng)業(yè)面源污染對環(huán)境質(zhì)量的影響越來越顯著[38]，基于區(qū)域異質(zhì)性的化肥面源污染研究能使研究結(jié)果更加符合區(qū)域的實(shí)際情況，從而更加準(zhǔn)確地總結(jié)和預(yù)測我國化肥面源污染狀況。當(dāng)前有學(xué)者結(jié)合區(qū)域的實(shí)際情況對中國化肥面源污染環(huán)境風(fēng)險(xiǎn)時(shí)空變化進(jìn)行分析，發(fā)現(xiàn)各省差異較大，地域聯(lián)系較為明顯，呈現(xiàn)出聚集分布的特點(diǎn)[39]，但并未探討各區(qū)域差異的規(guī)律及原因。有學(xué)者對化肥施用及環(huán)境風(fēng)險(xiǎn)的時(shí)空變化特征進(jìn)行了模擬，得出了區(qū)域間差異化的結(jié)論，然而并未考慮經(jīng)濟(jì)等因素對環(huán)境風(fēng)險(xiǎn)的影響[40]。為研究區(qū)域異質(zhì)性，有學(xué)者通過面板門檻模型選取變量進(jìn)行區(qū)域分組，結(jié)合經(jīng)濟(jì)發(fā)展?fàn)顩r，對化肥面源污染狀況進(jìn)行異質(zhì)性條件下的EKC 曲線檢驗(yàn)[22]，然而數(shù)理方法選取的分組變量可能存在主觀性，使得驗(yàn)證結(jié)果不能準(zhǔn)確反映各區(qū)域的特質(zhì)，組間差異性也不明顯。本文在我國化肥面源污染時(shí)空格局特征和相關(guān)性分析的基礎(chǔ)上，根據(jù)面源污染的空間異質(zhì)性分組，對化肥面源污染與農(nóng)業(yè)經(jīng)濟(jì)發(fā)展之間的EKC 時(shí)間路徑進(jìn)行模擬與趨勢預(yù)測，使研究結(jié)果更加符合區(qū)域的客觀情況，并得出了區(qū)域間差異化的結(jié)論，其中一些現(xiàn)象值得進(jìn)一步討論。

表9 部分省(市、自治區(qū))化肥面源污染跨過下一拐點(diǎn)的時(shí)間預(yù)測Table9 Time forecast for some provinces (cities，autonomous regions) cross the next inflection point of non-source pollution of fertilizers

1)從化肥面源污染空間格局角度看，我國化肥面源污染狀況呈現(xiàn)出顯著的空間正自相關(guān)性，表明存在較強(qiáng)的空間集聚效應(yīng)。近年來，空間自相關(guān)性有所減弱，長期的集聚格局可能發(fā)生一定變動(dòng)，不同區(qū)域間的流動(dòng)性可能增強(qiáng)，環(huán)境污染的空間溢出效應(yīng)開始愈發(fā)顯現(xiàn)[41]。其原因可能是各區(qū)域間貿(mào)易[42-43]和產(chǎn)業(yè)結(jié)構(gòu)的調(diào)整和流動(dòng)[44-47]帶來的污染物的轉(zhuǎn)移，因此部分省份從熱點(diǎn)區(qū)消失的部分原因可能是通過產(chǎn)業(yè)轉(zhuǎn)移將污染物轉(zhuǎn)移到了非熱點(diǎn)區(qū)域。對于北京、上海等省市而言，產(chǎn)業(yè)上的轉(zhuǎn)移有可能提高區(qū)域的生產(chǎn)效率，降低化肥面源污染的總體程度，因此可以進(jìn)一步適當(dāng)推進(jìn)農(nóng)業(yè)產(chǎn)業(yè)的轉(zhuǎn)移，同時(shí)通過生態(tài)補(bǔ)償、排污權(quán)交易等方式促進(jìn)區(qū)域間的利益協(xié)調(diào)。從根本而言，則應(yīng)當(dāng)進(jìn)一步推廣農(nóng)業(yè)新技術(shù)的運(yùn)用，從根本上減少全國范圍內(nèi)的化肥面源污染。

2)從EKC 模型的分區(qū)模擬結(jié)果來看，各區(qū)域在曲線形狀和拐點(diǎn)上存在較大的差異。當(dāng)前大多省份均處于曲線的上升階段，但距離拐點(diǎn)的時(shí)間則差異較大。熱點(diǎn)區(qū)的多數(shù)省份距離拐點(diǎn)較近，這表明熱點(diǎn)區(qū)近年來化肥面源污染防控工作取得一定成效。但熱點(diǎn)區(qū)迎來拐點(diǎn)后，非熱(冷)點(diǎn)區(qū)和冷點(diǎn)區(qū)卻依然會(huì)在較長期內(nèi)處于曲線上升階段。其原因一方面可能是以上討論的污染物空間上的轉(zhuǎn)移，另一方面則是非熱(冷)點(diǎn)區(qū)的化肥面源污染防控工作亟待進(jìn)一步加強(qiáng)，應(yīng)提高農(nóng)業(yè)新技術(shù)的普及率，推廣有機(jī)肥和配方施肥。

3)從控制變量的驅(qū)動(dòng)作用來看，3 類區(qū)域農(nóng)業(yè)機(jī)械強(qiáng)度均正向影響化肥面源污染，這表明當(dāng)前農(nóng)業(yè)生產(chǎn)過程中，機(jī)械水平的提高在提高農(nóng)業(yè)生產(chǎn)效率的同時(shí)也增加了面源污染的風(fēng)險(xiǎn)，因此在機(jī)械生產(chǎn)過程中應(yīng)當(dāng)提高精準(zhǔn)施肥和配方施肥技術(shù)的運(yùn)用。農(nóng)村居民工資性收入的提高也使農(nóng)業(yè)生產(chǎn)對于農(nóng)民生計(jì)影響變小，可能致使化肥施用更加粗放，種養(yǎng)循環(huán)難度變大。此外，相同因素在不同區(qū)域可能產(chǎn)生不同的影響，該結(jié)論與崔鑫生等[48]“同一因素對不同國家的作用存在著顯著差異”的研究結(jié)論具有相似之處。糧食播種比例、產(chǎn)業(yè)結(jié)構(gòu)的變化、種養(yǎng)結(jié)構(gòu)的變化在不同區(qū)域會(huì)對化肥面源污染產(chǎn)生不同的影響。如畜牧業(yè)比重較大的冷點(diǎn)區(qū)，種植業(yè)比重的下降明顯地對化肥面源污染產(chǎn)生積極影響。

由于各區(qū)域在空間上的相關(guān)性，因此在化肥面源污染治理中既要根據(jù)各區(qū)域的特性施策，又要從整體上注重區(qū)域間的相互作用。對于熱點(diǎn)區(qū)大部分省份，應(yīng)當(dāng)進(jìn)一步大力推廣測土配方施肥，精準(zhǔn)施肥，研發(fā)推廣適用施肥設(shè)備，機(jī)械深施，提高化肥利用率。冷點(diǎn)區(qū)部分省份生態(tài)較為脆弱，應(yīng)實(shí)施保護(hù)性耕作，增施有機(jī)肥，推廣水肥一體化等高效節(jié)水灌溉技術(shù)。非熱(冷)點(diǎn)區(qū)應(yīng)當(dāng)注重種養(yǎng)循環(huán)，推進(jìn)有機(jī)肥的利用。從整體而言，在化肥面源污染治理中，應(yīng)當(dāng)注意區(qū)域間的相互影響，尤其是注重產(chǎn)業(yè)結(jié)構(gòu)調(diào)整可能帶來的污染轉(zhuǎn)移，因此應(yīng)當(dāng)根據(jù)各區(qū)域省份的實(shí)際情況，適當(dāng)調(diào)整產(chǎn)業(yè)結(jié)構(gòu)，協(xié)調(diào)好經(jīng)濟(jì)效益和生態(tài)效益的平衡，通過生態(tài)補(bǔ)償、排污權(quán)交易等手段保障區(qū)域間發(fā)展均衡。

4.2 結(jié)論

本文在1997—2018年中國31 個(gè)省(市、自治區(qū))化肥面源污染排放強(qiáng)度核算的基礎(chǔ)上，運(yùn)用ArcGIS空間自相關(guān)分析和熱點(diǎn)分析空間統(tǒng)計(jì)模型，分析了化肥面源污染排放強(qiáng)度時(shí)空格局演變特征，并根據(jù)該特征，將全國分為化肥面源污染排放強(qiáng)度的熱點(diǎn)區(qū)、冷點(diǎn)區(qū)和非熱(冷)點(diǎn)區(qū)3 組，在空間異質(zhì)性條件下模擬和預(yù)測各組化肥面源污染和人均農(nóng)業(yè)產(chǎn)值間的時(shí)間路徑關(guān)系，得出以下結(jié)論：

1)我國化肥面源污染排放強(qiáng)度總體較高，且省際差異較大。從時(shí)間上看，各省化肥面源污染排放強(qiáng)度變化趨勢各異，表現(xiàn)出不同的升降趨勢；從空間上看，空間自相關(guān)分析表明，我國化肥面源污染排放強(qiáng)度表現(xiàn)出顯著的空間正向自相關(guān)，呈現(xiàn)為集聚模式，集聚程度近年有所下降；熱點(diǎn)分析顯示，我國1997—2018年化肥面源污染排放強(qiáng)度時(shí)空格局演變較為穩(wěn)定，熱點(diǎn)區(qū)主要分布在華中地區(qū)和華南地區(qū)，黃淮海地區(qū)的熱點(diǎn)程度下降較為明顯，冷點(diǎn)區(qū)主要集中在西部地區(qū)和東北的黑龍江。各區(qū)域間表現(xiàn)出一定的空間溢出效應(yīng)。

2)化肥面源污染排放強(qiáng)度時(shí)空格局演變特征分組的EKC 模擬結(jié)果表明，我國化肥面源污染與農(nóng)業(yè)經(jīng)濟(jì)增長之間存在著顯著的非線性EKC 關(guān)系，并且各組EKC 時(shí)間路徑曲線形狀和拐點(diǎn)存在明顯差異。熱點(diǎn)區(qū)呈現(xiàn)為“倒U 型”曲線、冷點(diǎn)區(qū)和非熱(冷)點(diǎn)區(qū)呈現(xiàn)為“倒N 型”曲線。我國大多數(shù)省份正處在化肥面源污染狀況隨農(nóng)業(yè)經(jīng)濟(jì)發(fā)展而惡化的階段，且惡化的原因由于區(qū)域特質(zhì)的不同存在差異，亟需進(jìn)一步因時(shí)因地推進(jìn)各區(qū)域化肥面源污染防治，并注意區(qū)域間的相互影響，實(shí)現(xiàn)區(qū)域間的協(xié)同治理。

3)化肥面源污染EKC 的時(shí)間路徑預(yù)測表明，熱點(diǎn)區(qū)上海、冷點(diǎn)區(qū)的青海和西藏以及非熱(冷)點(diǎn)區(qū)的多數(shù)省份距離曲線拐點(diǎn)仍較遠(yuǎn)，需因時(shí)因地施策，促使拐點(diǎn)盡快到來。