余秋菊，武以敏，高鳳偉，欒詩(shī)晴

宿州學(xué)院統(tǒng)計(jì)調(diào)查咨詢服務(wù)中心，安徽宿州，234000

山東省近二十年來不斷擴(kuò)大規(guī)?？偭?，優(yōu)化產(chǎn)業(yè)結(jié)構(gòu)，其中農(nóng)業(yè)發(fā)展取得了輝煌成就，山東省已成為農(nóng)業(yè)生產(chǎn)大省之一。2018年，山東省農(nóng)業(yè)增加值約占農(nóng)林牧漁業(yè)總增加值的55.143%。隨著農(nóng)業(yè)生產(chǎn)規(guī)模的擴(kuò)大以及農(nóng)藥、化肥、地膜等各種生產(chǎn)要素的投入，山東省農(nóng)業(yè)面源污染逐漸加重，主要產(chǎn)生水體富營(yíng)養(yǎng)化、土壤次生鹽漬化和耕地通透性變差等污染現(xiàn)象。中共十八大和十九大會(huì)議中分別提出了建設(shè)“資源節(jié)約型、環(huán)境友好型”社會(huì)和“綠水青山就是金山銀山”理念，強(qiáng)調(diào)了人與自然和諧共生，堅(jiān)定走生態(tài)文明發(fā)展道路。2021年，國(guó)家相關(guān)部門發(fā)布了《“十四五”全國(guó)農(nóng)業(yè)綠色發(fā)展規(guī)劃》，其中將加強(qiáng)農(nóng)業(yè)面源污染防治作為總體要求之一。在此背景下，積極探討山東省農(nóng)業(yè)面源污染的來源以及影響因素具有重要的理論意義及實(shí)際意義。

目前，國(guó)內(nèi)外學(xué)者圍繞生態(tài)文明建設(shè)及綠色發(fā)展進(jìn)行了研究，且環(huán)境污染一直是研究的熱門話題。如屈文波等[1]從非正式環(huán)境規(guī)制減排角度出發(fā)，利用動(dòng)態(tài)空間面板模型探究公眾參與對(duì)環(huán)境污染的影響，研究表明公眾參與度的提高能夠有效減少環(huán)境污染物的排放，且環(huán)境污染具有空間溢出效應(yīng)；施震凱等[2]利用面板數(shù)據(jù)模型研究進(jìn)口復(fù)雜度對(duì)霧霾污染的影響效應(yīng)，實(shí)證發(fā)現(xiàn)進(jìn)口復(fù)雜度的提升能夠促進(jìn)產(chǎn)業(yè)結(jié)構(gòu)升級(jí)及技術(shù)創(chuàng)新等，進(jìn)而對(duì)霧霾污染具有正向積極作用；姚成勝等[3]運(yùn)用面板數(shù)據(jù)模型對(duì)全國(guó)31個(gè)省份的環(huán)境污染影響因素進(jìn)行研究，研究表明大部分地區(qū)的工業(yè)集聚水平對(duì)環(huán)境污染產(chǎn)生了抑制作用，并且土地城鎮(zhèn)化水平、工業(yè)集聚水平對(duì)全國(guó)的環(huán)境污染均具有正向促進(jìn)作用。梳理并歸納文獻(xiàn)[4-6]可知，關(guān)于綠色發(fā)展水平的研究主要體現(xiàn)在綠色發(fā)展效率的測(cè)算及影響因素分析，研究表明經(jīng)濟(jì)發(fā)展水平、產(chǎn)業(yè)結(jié)構(gòu)、科技進(jìn)步、外商直接投資等因素對(duì)綠色發(fā)展具有顯著的正向促進(jìn)作用。同時(shí)注意到，國(guó)內(nèi)外大多數(shù)學(xué)者對(duì)高質(zhì)量經(jīng)濟(jì)發(fā)展的研究?jī)H僅從工業(yè)污染、霧霾污染、綠色發(fā)展等視角進(jìn)行研究，而對(duì)農(nóng)業(yè)面源污染進(jìn)行量化分析的研究偏少。

在此背景下，本文選取山東省17個(gè)地市作為研究對(duì)象，利用清單分析法和固定效應(yīng)模型探究山東省17個(gè)地市的農(nóng)業(yè)面源污染排放量及其影響因素，并運(yùn)用Moran’I指數(shù)、LISA散點(diǎn)圖、集聚圖等分析山東省17個(gè)地市農(nóng)業(yè)面源污染是否具有空間溢出效應(yīng)。

1 研究方法

本文研究目的是運(yùn)用清單分析法測(cè)算山東省17個(gè)地市的農(nóng)業(yè)面源污染排放量，且根據(jù)Moran’I指數(shù)、LISA散點(diǎn)圖、集聚圖等分析山東省17個(gè)地市的農(nóng)業(yè)面源污染排放的空間效應(yīng)，并利用面板數(shù)據(jù)模型探究農(nóng)業(yè)面源污染的影響因素。

1.1 農(nóng)業(yè)面源污染排放量的測(cè)算

查閱文獻(xiàn)[7]，利用清單分析法測(cè)算農(nóng)業(yè)面源污染物全氮(TN)、全磷(TP)、化學(xué)需氧量(COD)。在現(xiàn)有研究基礎(chǔ)上進(jìn)行擴(kuò)展，得到農(nóng)業(yè)面源污染產(chǎn)污清單表，如表1所示。

表1 農(nóng)業(yè)面源污染產(chǎn)污清單

文中選取上表中五類污染源，其污染源產(chǎn)污系數(shù)分別來源于：肥料施用污染源中使用的系數(shù)參考梁流濤和張佳卓的文獻(xiàn)[8-9]、畜禽養(yǎng)殖中使用的系數(shù)參考《全國(guó)污染普查畜禽養(yǎng)殖業(yè)產(chǎn)污系數(shù)與排污系數(shù)手冊(cè)》、水產(chǎn)養(yǎng)殖中使用的排污系數(shù)參考《水產(chǎn)養(yǎng)殖業(yè)污染源產(chǎn)排污系數(shù)手冊(cè)》、農(nóng)作物污染系數(shù)參考《全國(guó)種植業(yè)污染源普查排污系數(shù)測(cè)算實(shí)施方案》，其污染排放量的計(jì)算公式如下：

其中，E表示農(nóng)業(yè)面源污染排放量；EUi表示污染單元i的統(tǒng)計(jì)量；Pi為污染單元i的產(chǎn)生系數(shù)；Ci為污染單元i的流失系數(shù);PEi為農(nóng)業(yè)面源污染產(chǎn)生量。通過上面的污染排放量計(jì)算公式，可分別測(cè)算2009—2018年山東省及17個(gè)地市的農(nóng)業(yè)面源污染物全氮、全磷、化學(xué)需氧量排放量。

1.2 空間效應(yīng)檢驗(yàn)

運(yùn)用全局Moran′s I指數(shù)檢驗(yàn)山東省農(nóng)業(yè)面源污染的整體空間相關(guān)性，用局部Moran′s I指數(shù)檢驗(yàn)各區(qū)域與周邊區(qū)域的空間相關(guān)性，用LISA統(tǒng)計(jì)量檢驗(yàn)空間集聚現(xiàn)象。其全局Moran′s I指數(shù)和局部Moran′s I指數(shù)的計(jì)算公式分別如下[10]：

1.3 面板數(shù)據(jù)模型

面板數(shù)據(jù)同時(shí)含有橫截面和時(shí)間序列的數(shù)據(jù)，是對(duì)一組固定調(diào)查對(duì)象的多次觀測(cè)得到的數(shù)據(jù)，即由橫截面上個(gè)體在不同時(shí)間的重復(fù)觀測(cè)而形成的數(shù)據(jù)。在建立面板數(shù)據(jù)模型需要確定模型類別：根據(jù)對(duì)截距項(xiàng)和解釋變量系數(shù)的不同假設(shè)，可以將面板數(shù)據(jù)回歸模型具體劃分為混合回歸模型、固定效應(yīng)模型和隨機(jī)效應(yīng)模型三種，其中固定效應(yīng)模型公式如下[10]：

yit=xitβ+zi′δ+ri+ui+εit

其中,yit,xit分別表示因變量和自變量在橫截面i和時(shí)間t上的數(shù)值，β,δ分別是系數(shù)項(xiàng),zi是不隨時(shí)間變化的個(gè)體特征，ui+εit是復(fù)合擾動(dòng)項(xiàng)，i=1,2,…N表示截面?zhèn)€體數(shù)，t=1,2,…N表示對(duì)每個(gè)截面的觀察時(shí)點(diǎn)數(shù)。

2 實(shí)證分析

2.1 時(shí)序變化分析

由于數(shù)據(jù)可獲取性，文中選取2009—2018年面源污染指標(biāo)數(shù)據(jù)進(jìn)行分析。根據(jù)清單分析法可得2009—2018年山東省農(nóng)業(yè)面源污染物全氮、全磷、化學(xué)需氧量的排放量，由于年份較多，本文選取2018年的農(nóng)業(yè)面源污染排放量進(jìn)行分析。如表2所示。

表2 2018年山東省農(nóng)業(yè)面源污染物排放量

由表2可知，2018年山東省農(nóng)業(yè)面源污染物排放量較多的是TN和COD，其分別為555 660.9 t、203 107.3 t，兩者分別占總污染排放量的66.35%、24.25%；而TP排放量為78 661.47 t，僅占總污染排放量的9.39%。并且從各類污染源看，由畜禽養(yǎng)殖產(chǎn)生的污染物為577 298.6 t，其占總污染排放量的68.94%；而由農(nóng)作物產(chǎn)生的污染物占總污染排放量的17.33%，這兩種污染源是農(nóng)業(yè)面源污染產(chǎn)生的主要部分。同理觀察2009—2018年山東省農(nóng)業(yè)面源污染物排放量也可發(fā)現(xiàn)，各年山東省畜禽養(yǎng)殖和農(nóng)作物產(chǎn)生的污染物也居于前位。

由圖1可知，2009—2018年山東省農(nóng)業(yè)面源污染物排放量最多的是TN，最少的是TP。此外，2009—2018年山東省農(nóng)業(yè)面源污染物TN、TP和COD處于較小波動(dòng)狀態(tài)，其中TP長(zhǎng)期內(nèi)呈平穩(wěn)趨勢(shì)，TN和COD從2014—2018年整體上呈緩慢下降趨勢(shì)。

圖1 2009—2018年山東省農(nóng)業(yè)面源污染物排放量

根據(jù)圖2可知，將2009年作為比較基期，2010—2018年山東省農(nóng)業(yè)面源污染物排放量增長(zhǎng)率變化趨勢(shì)較大的是2018年，且三種污染物排放量增長(zhǎng)率達(dá)到最大的負(fù)增長(zhǎng)。此外，2010—2014年山東省農(nóng)業(yè)面源污染物COD排放量增長(zhǎng)率均為正增長(zhǎng)，而2015—2018年三種農(nóng)業(yè)面源污染物排放量增長(zhǎng)率均處于負(fù)增長(zhǎng)。這表明2010—2014年山東省農(nóng)業(yè)種植和水產(chǎn)養(yǎng)殖過程中產(chǎn)生了大量的污染物COD，但從2015—2018年山東省逐漸意識(shí)到生態(tài)環(huán)境的重要性，注重綠色發(fā)展之路，進(jìn)而農(nóng)業(yè)面源污染物排放量逐步較弱。

圖2 2010—2018年污染物排放量增長(zhǎng)率變化趨勢(shì)圖

2.2 空間效應(yīng)分析

由于涉及年份較多，本文利用2018年數(shù)據(jù)對(duì)山東省17個(gè)地市農(nóng)業(yè)面源污染排放量進(jìn)行空間效應(yīng)分析。首先利用全局Moran’s I指數(shù)檢驗(yàn)山東省農(nóng)業(yè)面源污染排放量是否具有空間相關(guān)性，檢驗(yàn)結(jié)果如圖3所示。

圖3 2018年山東省各地市Moran’s I散點(diǎn)圖

由圖3可知，2018年山東省17個(gè)地市的Moran′s I指數(shù)為0.289，大于0。另外，在5%的顯著性水平下，p值小于0.05。由此表明，2018年山東省17個(gè)地市的農(nóng)業(yè)面源污染排放總量具有顯著的空間正向相關(guān)性，也即是山東省整體農(nóng)業(yè)面源污染具有空間溢出性。同理，若分別對(duì)TN、TP、COD污染排放量進(jìn)行全局Moran′s I指數(shù)檢驗(yàn)，可得TN、TP、COD污染排放量的全局Moran′s I指數(shù)在5%的顯著性水平下，三者均顯著。由此表明，2018年山東省各類農(nóng)業(yè)面源污染物也分別具有空間溢出性。下面利用局部空間相關(guān)性檢驗(yàn)?zāi)车貐^(qū)與周邊地區(qū)是否具有空間相關(guān)性，輸出結(jié)果如圖4所示。

圖4 2018年山東省各地市LISA顯著性和集聚圖

根據(jù)LISA顯著圖4可知，在1%的顯著性水平下，2018年山東省德州市和聊城市的農(nóng)業(yè)面源污染總量通過了檢驗(yàn)。在5%的顯著性水平下，2018年山東省濟(jì)南市和濰坊市農(nóng)業(yè)面源污染總量通過了檢驗(yàn)。由LISA集聚圖發(fā)現(xiàn)，2018年山東省的濟(jì)南市處于H-L象限，濰坊市處于L-L象限，德州市和聊城市均處于L-H象限。此外，從LISA和集聚圖發(fā)現(xiàn)，數(shù)十年山東省農(nóng)業(yè)面源污染聚集情況沒有產(chǎn)生較大變化。

2.3 農(nóng)業(yè)面源污染影響因素分析

選取2009—2018年山東省各地市的年末總?cè)丝跀?shù)(PE)、化肥使用量(HS)、農(nóng)藥使用量(NS)、地膜(DM)、灌溉面積(GG)、柴油使用量(NC)作為解釋變量；用上文測(cè)算出的2009—2018年山東省各地市的TP、TN、COD排放總量分別作為被解釋變量，用Y表示。所有數(shù)據(jù)均來源于2010—2019年《山東省統(tǒng)計(jì)年鑒》、各地市《統(tǒng)計(jì)公報(bào)》等。

首先利用F檢驗(yàn)和Hausman檢驗(yàn)確定適合的面板數(shù)據(jù)模型，得到的結(jié)果如下所示：

根據(jù)表3可知，以TN為被解釋變量的模型中，P值小于0.05，則表明應(yīng)接受備擇假設(shè)，建立固定效應(yīng)模型。

表3 污染物全氮隨機(jī)效應(yīng)模型的Hausman檢驗(yàn)

此外，通過F檢驗(yàn)發(fā)現(xiàn)(表4)，P值小于0.05，則表明應(yīng)拒絕原假設(shè)，建立固定效應(yīng)模型，這也與Hausman檢驗(yàn)結(jié)果相一致。同理，通過對(duì)TP和COD模型的一系列檢驗(yàn)發(fā)現(xiàn)，以TP和COD為被解釋變量的模型也應(yīng)建立固定效應(yīng)模型。

表4 污染物全氮固定效應(yīng)模型的F檢驗(yàn)

其次，文中分別以TN、TP、COD為被解釋變量建立固定效應(yīng)模型，整理結(jié)果如表5所示。

表5 固定效應(yīng)模型回歸結(jié)果

從表5可知，以TN、TP、COD為被解釋變量的三個(gè)固定效應(yīng)模型中，模型的F統(tǒng)計(jì)量都顯著超過相應(yīng)臨界值水平，這表明三個(gè)固定效應(yīng)模型的整體顯著性水平較好。觀察以TN被解釋變量的固定效應(yīng)模型發(fā)現(xiàn)，HS、NS、GG的回歸系數(shù)在1%的顯著性水平下顯著，且HS、NS的系數(shù)為正值，GG的系數(shù)為負(fù)值，這表明化肥使用量和農(nóng)藥使用量對(duì)TN污染物的排放具有顯著的正向作用，農(nóng)田灌溉面積對(duì)TN污染物的排放具有顯著的負(fù)向作用。DM的回歸系數(shù)在10%的顯著性水平下顯著，且DM的系數(shù)為正值，這表明地膜使用量對(duì)TN污染物的排放具有顯著的正向作用。同理，觀察以TP、COD為被解釋變量的固定效應(yīng)模型發(fā)現(xiàn)，PE、HS、NS的回歸系數(shù)在1%的顯著性水平下顯著，且年末總?cè)丝跀?shù)、化肥使用量對(duì)TP污染物的排放具有顯著的正向作用。PE、HS、GG的回歸系數(shù)在1%的顯著性水平下顯著，且年末總?cè)丝跀?shù)、化肥使用量、農(nóng)田灌溉面積對(duì)COD污染物的排放具有顯著的正向作用，這也與胡鞍鋼[11]、葛繼紅[12]等人的研究結(jié)果相一致。從回歸系數(shù)絕對(duì)值發(fā)現(xiàn)，三個(gè)模型中柴油使用量的系數(shù)較小，也表明這個(gè)變量對(duì)三種污染物的排放具有較小的影響。

3 結(jié) 語

文中主要利用清單分析法測(cè)算2009—2018年山東省各地市的農(nóng)業(yè)面源污染排放量，并對(duì)各地市的農(nóng)業(yè)面源污染排放量進(jìn)行時(shí)空效應(yīng)分析，以及運(yùn)用固定效應(yīng)模型分析農(nóng)業(yè)面源污染排放的影響因素，得到的結(jié)論有：第一，2009—2018年山東省整體農(nóng)業(yè)面源污染中TN和COD的排放量較多，TP最少；且2015—2018年山東省農(nóng)業(yè)面源污染物排放量增長(zhǎng)率均處于負(fù)增長(zhǎng)。第二，2009—2018年山東省農(nóng)業(yè)面源污染TN、TP、COD的排放量具有空間溢出性，且德州市、聊城市以及濰坊市的農(nóng)業(yè)污染排放量相對(duì)較高。第三，人口數(shù)、化肥使用量以及地膜等變量對(duì)TN、TP、COD的排放具有正向作用，也即是表明TN、TP、COD的排放量會(huì)隨著人口數(shù)、化肥和地膜使用量的增加而增加。因此，筆者提出以下建議：第一，政府應(yīng)積極提倡綠色發(fā)展，鼓勵(lì)居民和企業(yè)投身于生態(tài)文明建設(shè)，加強(qiáng)可持續(xù)發(fā)展理念。第二，在農(nóng)業(yè)發(fā)展進(jìn)程中，合理控制并規(guī)劃農(nóng)藥、地膜、化肥等污染環(huán)境的生產(chǎn)要素投入，積極引進(jìn)先進(jìn)農(nóng)業(yè)技術(shù)水平，提高生產(chǎn)效率和能源使用效率。第三，各地市政府要重點(diǎn)關(guān)注農(nóng)業(yè)面源污染的空間效應(yīng)，努力減少污染源，共同治理環(huán)境污染，協(xié)調(diào)發(fā)展。