雷俊華，蘇時鵬，余文夢，孫小霞**

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

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

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

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

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

1.1 研究方法

1.1.1 面源污染排放強度核算

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

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

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

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

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

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

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

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

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

1.1.3 EKC 模型設(shè)定

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

式中：I表示化肥面源污染排放強度；G表示人均農(nóng)業(yè)產(chǎn)值，作為解釋變量，用來表征農(nóng)業(yè)經(jīng)濟發(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)值計算)。其他指標(biāo)是參考現(xiàn)有相關(guān)研究選取的控制變量：P表示作物種植結(jié)構(gòu)，用糧食播種面積與農(nóng)作物總播種面積的比值表征，作物種植結(jié)構(gòu)差異可能致使化肥施用量的差異；M表示農(nóng)業(yè)機械投入強度(kg·hm-2)，用單位農(nóng)作物播種面積的農(nóng)業(yè)機械投入強度表征，機械化水平提高可能會對化肥施用量產(chǎn)生正負兩面的影響；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)會影響當(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)村居民收入的高低會影響化肥的投入；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)整對化肥面源污染的影響；μ為隨機誤差項；i表示省份，t表示年份；回歸模型的估計系數(shù)βn(n=1～10)表示彈性系數(shù)。

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

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

2 化肥面源污染排放強度核算及其時空格局演變

2.1 化肥面源污染排放強度核算

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

2.2 省際化肥面源污染排放強度時空格局演變特征

2.2.1 總體時空格局演化特征

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

表2 1997—2018年各省(市、自治區(qū))化肥面源污染排放強度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

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

2.2.2 局部時空格局演化特征

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

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

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

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

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

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

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

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

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

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

表5 化肥面源污染排放強度不同熱點區(qū)域的面板單位根(LLC)檢驗結(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 化肥面源污染排放強度不同熱點區(qū)域的協(xié)整檢驗結(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)的檢驗

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

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

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

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

1)熱點區(qū)。本組包含上海、江蘇、浙江、河南、安徽、湖北等10 個省份，這些省份主要是我國經(jīng)濟發(fā)達地區(qū)或農(nóng)業(yè)大省，化肥面源污染強度普遍較大。本組EKC 的時間路徑呈現(xiàn)為“倒U 型”曲線，說明伴隨農(nóng)業(yè)發(fā)展，化肥面源污染表現(xiàn)為“升—降”趨勢。該曲線拐點對應(yīng)的人均農(nóng)業(yè)生產(chǎn)總值G 分別為5 889.933 元。從熱點區(qū)10 省的情況看，江蘇(2011年)、山東(2016年)、福建(2017年)、河南(2017年)、河北(2017年)5 個省份已經(jīng)跨過了“倒U 型”曲線的拐點，進入了“隨著農(nóng)業(yè)經(jīng)濟增長，化肥面源污染減弱”的下降階段。從時間上看，多數(shù)省份在2017年左右跨過拐點，其原因除經(jīng)濟不斷增長促使化肥面源污染到達拐點外，國家糧食安全戰(zhàn)略下糧食主產(chǎn)區(qū)近年糧食播種面積比重的回升、農(nóng)村居民人均純收入的提高和農(nóng)村居民對農(nóng)業(yè)收入依賴的降低都可能促進了拐點的到來。此外，2015年以來，國家頒布的一系列化肥減量化政策，如《到2020年化肥使用量零增長行動方案》以及近年中央一號文件對農(nóng)業(yè)面源污染重視程度的不斷加強也可能促進了拐點的到來。安徽、湖南、江西、上海和浙江則均尚未到達拐點，仍處于“隨著農(nóng)業(yè)經(jīng)濟增長，化肥面源污染增強”的階段。安徽、湖南和江西作為我國重要的產(chǎn)糧區(qū)，對于化肥的依賴程度較高，難以在短期擺脫對化肥的依賴。上海和浙江作為糧食主銷區(qū)，較低的種植業(yè)比重反而致使化肥施用的集約化程度更低，發(fā)展種養(yǎng)循環(huán)的動力更小，從而未充分發(fā)揮有機肥代替化肥等方面的作用。熱點區(qū)農(nóng)業(yè)機械投入強度與化肥面源污染程度呈正相關(guān)，表明科技和機械水平雖然在提高農(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)檢驗結(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)冷點區(qū)。本組包含甘肅、青海、西藏、黑龍江4 個省份，主要分布在我國西部和東北。本組EKC 的時間路徑為“倒N 型”曲線，說明伴隨農(nóng)業(yè)發(fā)展，化肥面源污染表現(xiàn)為“降—升—降”趨勢。該曲線的拐點對應(yīng)的人均農(nóng)業(yè)生產(chǎn)總值G分別為 851.500 元和4 433.746 元。冷點區(qū)4 省均在1997年以前就跨過了“倒N 型”曲線的第1 個拐點，已經(jīng)跨過第2 個拐點的省份為甘肅(2017年)和黑龍江(2006年)，表明這兩個省份目前處于曲線下降階段，即“隨著農(nóng)業(yè)經(jīng)濟增長，化肥面源污染減弱”。青海和西藏則還未跨過第2 個拐點，仍處于曲線的上升階段，且表現(xiàn)出短期內(nèi)不會達到拐點的趨勢。對于青海和西藏而言，其本身較低的農(nóng)業(yè)經(jīng)濟增長水平促使其距離拐點較遠。黑龍江農(nóng)業(yè)產(chǎn)值高但化肥面源污染排放強度卻一直較低，其原因在于東北黑土資源肥沃、供肥能力強，且多年來耕地質(zhì)量普遍不斷提高[37]，化肥流失率低，從而使黑龍江較早跨過拐點進入曲線下降階段。

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

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

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

4 討論與結(jié)論

4.1 討論

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

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

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

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

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

4.2 結(jié)論

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

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

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

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