中國(guó)霧霾污染的空間計(jì)量分析

林楚海

（暨南大學(xué) 經(jīng)濟(jì)學(xué)院，廣州 510000）

0 引言

霧霾，即霧和霾的組合詞。霧是由大量懸浮在近地面空氣中的微小水滴或冰晶組成的氣溶膠系統(tǒng)，而霾是指由于灰塵、硫酸、硝酸、有機(jī)碳?xì)浠衔锏攘Ｗ討腋≡诳諝庵?，而使水平能?jiàn)度小于1000米一種現(xiàn)象。霧霾天氣主要由空氣中的二氧化硫、氮氧化物和可吸入顆粒物這三項(xiàng)引起，它們與霧氣結(jié)合在一起，會(huì)讓天空瞬間變得陰沉灰暗，也對(duì)人體和社會(huì)帶來(lái)各方面的影響。2014年，中國(guó)31個(gè)省份中，PM10年平均濃度達(dá)到國(guó)家一級(jí)標(biāo)準(zhǔn)的只有海南省，而SO2達(dá)到國(guó)家一級(jí)標(biāo)準(zhǔn)的也只有7個(gè)省份，從整體上看，霧霾污染依然不容樂(lè)觀。再加上霧霾的源頭多種多樣，比如汽車(chē)尾氣、工業(yè)排放、建筑揚(yáng)塵、垃圾焚燒等，因此，如何在眾多因素中快速找到造成霧霾的主導(dǎo)因素并采取有效的治理措施，成為中國(guó)政府和居民急需解決的難題。基于此，本文將選取中國(guó)31個(gè)省份2014年三種霧霾主要組分——PM10、SO2、NO2的年平均濃度作為目標(biāo)變量，并利用LASSO高維變量選擇方法從22個(gè)社會(huì)影響因子中篩選出非零變量進(jìn)行空間回歸分析，從而揭示中國(guó)內(nèi)地霧霾污染的真正來(lái)源和空間影響因素。

1 研究方法

1.1 空間計(jì)量分析

1.1.1 空間權(quán)重矩陣

空間權(quán)重矩陣即地理單元間的空間結(jié)構(gòu)關(guān)系的一種表達(dá)形式，描述了各單元間的鄰近關(guān)系?？臻g權(quán)重矩陣W一般表示為一個(gè)n維矩陣，其基本形式可表示如下：

其中，wij表示區(qū)域i和j的鄰近關(guān)系，一般分為基于鄰近概念的空間權(quán)重矩陣和基于距離的空間權(quán)重矩陣。其定義如下：

1.1.2 全局空間自相關(guān)

全局空間自相關(guān)描述的是研究范圍內(nèi)所有單元之間的整體空間關(guān)系。一般采用Moran’s I統(tǒng)計(jì)量來(lái)描述整個(gè)研究區(qū)域上所有空間單元之間的平均關(guān)聯(lián)程度及其顯著性。其計(jì)算公式如下：

其中，xi表示地區(qū)i的觀測(cè)值，，n為地區(qū)數(shù)，wij為空間權(quán)重矩陣。Moran’s I統(tǒng)計(jì)量的取值范圍為[-1，1]。I越接近1，表示區(qū)域間空間正相關(guān)的程度越強(qiáng)；當(dāng)I越接近-1時(shí)，表示區(qū)域間空間負(fù)相關(guān)的程度越強(qiáng)；當(dāng)I接近0時(shí)，則表示區(qū)域間不存在空間自相關(guān)性。

1.1.3 局域空間自相關(guān)

局域空間自相關(guān)統(tǒng)計(jì)量就是用來(lái)識(shí)別隨空間位置不同而存在的不同空間關(guān)聯(lián)模式，從而觀察空間局部不平穩(wěn)性。本文的局域空間自相關(guān)分析將采用局部Moran’s I指數(shù)繪制的集聚地圖和顯著性地圖。

局部Moran’s I指數(shù)的計(jì)算公式如下：

其中，Ii表示i地區(qū)與其相鄰區(qū)域之間的空間相關(guān)程度，其他符號(hào)的含義同式（4），這里不再贅述。當(dāng)Ii＞0時(shí)，表示第i個(gè)地區(qū)某種屬性值與其周?chē)貐^(qū)的屬性值呈空間正相關(guān)，表現(xiàn)為高高集聚或低低集聚；當(dāng)Ii＜0時(shí)表示負(fù)相關(guān)，表現(xiàn)為高低集聚或低高集聚。

1.1.4 空間回歸分析

（1）空間滯后模型：假定空間變量的空間依賴(lài)性?xún)H通過(guò)因變量產(chǎn)生，即某一空間對(duì)象上的因變量不僅與同一對(duì)象上的自變量有關(guān)，還與相鄰對(duì)象的因變量有關(guān)，其形式為：

其中，y為n×1階截面數(shù)據(jù)被解釋變量，W為空間權(quán)重矩陣，Wy即為被解釋變量的空間滯后向量，其系數(shù)ρ為空間自回歸系數(shù)，反映了相鄰區(qū)域?qū)Ρ緟^(qū)域的作用，當(dāng)ρ顯著時(shí)，說(shuō)明因變量間存在明顯的空間依賴(lài)性，當(dāng)ρ＞0時(shí)，表明相鄰空間區(qū)域的空間相互作用表現(xiàn)為“溢出效應(yīng)”;反之，則表現(xiàn)為“負(fù)效應(yīng)”或“擴(kuò)散效應(yīng)”。X即為n×k外生解釋變量矩陣，其中k為解釋變量的個(gè)數(shù)。ε為n×1階隨機(jī)誤差向量。

（2）空間誤差模型：在某些情況下，空間依賴(lài)可能是某些對(duì)局域內(nèi)多個(gè)單元均有影響的遺漏變量造成的，從而造成區(qū)域“一損俱損，一榮俱榮”的效果，這時(shí)候需要考慮誤差項(xiàng)的空間自回歸過(guò)程。其形式為：

其中，λ是誤差項(xiàng)的空間自回歸系數(shù)，Wε是誤差項(xiàng)的空間滯后向量，μ為不相關(guān)的、均值為0且同方差的誤差項(xiàng)。若λ顯著，則說(shuō)明了鄰近區(qū)域的關(guān)于因變量的誤差沖擊對(duì)本地區(qū)的觀測(cè)值具有顯著性影響。

（3）空間通用模型：空間通用模型綜合了空間滯后模型和空間誤差模型，其形式為：

其中，空間權(quán)重矩陣W1與W2可以相同，也可以不同，其他符號(hào)含義與上述兩個(gè)模型相同，這里不再贅述。

（4）空間相關(guān)性檢驗(yàn)與模型選擇：一般可通過(guò)兩個(gè)拉格朗日乘數(shù)LMERR、LMLAG及其穩(wěn)健的R-LMERR、R-LMLAG等統(tǒng)計(jì)量來(lái)實(shí)現(xiàn)。當(dāng)LMLAG比LMERR顯著，且R-LMLAG顯著而R-LMERR不顯著時(shí)，則選擇空間滯后模型更為合適；相反地，如果LMERR比LMLAG顯著，且R-LMERR顯著而R-LMLAG不顯著時(shí)，則選擇空間誤差模型更為合適。其他可供參考的模型比較準(zhǔn)則還有對(duì)數(shù)似然值、赤池信息準(zhǔn)則（AIC）和施瓦茨準(zhǔn)則（SC）等。當(dāng)對(duì)數(shù)似然值越大，AIC和SC值越小，則模型擬合效果越好。

1.2 LASSO高維變量選擇

設(shè)有p個(gè)自變量 x1，x2，…，xp和因變量y，并建立如下的線(xiàn)性回歸模型：

其中，α為常數(shù)項(xiàng)，β1，β2，…，βp為自變量的回歸系數(shù)，ε為隨機(jī)擾動(dòng)項(xiàng)。設(shè)(xi1，xi2，…，xip;yi)，i=1，2，…，n 為n組觀測(cè)值，假定數(shù)據(jù)已經(jīng)經(jīng)過(guò)中心化和標(biāo)準(zhǔn)化，即，…，βp)T,則 β 的最小二乘估計(jì)為：

1996年，Tibshirani提出了LASSO變量選擇方法，其基本思想是在最小二乘的基礎(chǔ)上施加l1懲罰項(xiàng)，將一些變量的系數(shù)恰好壓縮為零，從而實(shí)現(xiàn)變量選擇。其表達(dá)式如下：

其中，s≥0是一個(gè)懲罰參數(shù)，控制著壓縮的程度，s越小，壓縮程度越強(qiáng)，會(huì)有更多的系數(shù)被壓縮至零。記最小二乘估計(jì)得到的參數(shù)估計(jì)為實(shí)現(xiàn)壓縮，即有些系數(shù)會(huì)變?yōu)?，這些等于0的變量就會(huì)被刪除，從而達(dá)到變量選擇的目的。其等價(jià)形式如下：

若設(shè)計(jì)陣正交，即XTX=I，LASSO的參數(shù)估計(jì)為：

2 實(shí)證分析

2.1 數(shù)據(jù)說(shuō)明

本文選取中國(guó)31個(gè)省份2014年的PM10、SO2、NO2的年平均濃度（微克/立方米）作為被解釋變量，同時(shí)從產(chǎn)業(yè)結(jié)構(gòu)、人口因素、固定資產(chǎn)投資和房地產(chǎn)、能源因素、城市建設(shè)五個(gè)大方面選擇22個(gè)指標(biāo)作為自變量，數(shù)據(jù)均來(lái)源于各省2015年的《統(tǒng)計(jì)年鑒》和《環(huán)境統(tǒng)計(jì)公報(bào)》，如下頁(yè)表1所示。

采用ArcGis 10.0軟件計(jì)算空間權(quán)重矩陣，相鄰規(guī)則采用Rook相鄰。另外，由于海南島遠(yuǎn)離陸地，與其他省份沒(méi)有邊界相鄰，為了避免出現(xiàn)“孤島”現(xiàn)象，將鄰居個(gè)數(shù)最低值設(shè)為1，從而使得每個(gè)省份至少都有一個(gè)鄰居相鄰。

表1 指標(biāo)含義與說(shuō)明

2.2 全局空間自相關(guān)分析

通過(guò)Geoda軟件進(jìn)行全局空間自相關(guān)分析，得到31個(gè)省份2014年的PM10、SO2、NO2年平均濃度對(duì)應(yīng)的全局Moran’s I指數(shù)，如表2所示。結(jié)果表明，PM10、SO2、NO2年平均濃度對(duì)應(yīng)的全局Moran’s I指數(shù)分別為0.5377、0.3501、0.5044，并且對(duì)應(yīng)的P值均小于0.05，說(shuō)明31個(gè)省份之間的PM10、SO2、NO2的空間自相關(guān)性非常顯著，且呈正的空間自相關(guān)，而非完全隨機(jī)分布。

表2 全局Moran’s I指數(shù)及相關(guān)統(tǒng)計(jì)指標(biāo)

2.3 局域空間自相關(guān)分析

在現(xiàn)實(shí)中，各地區(qū)霧霾組分受到的空間影響的模式完全一致一般是很少見(jiàn)的，為了進(jìn)一步探究各省份與相鄰省份之間的空間影響關(guān)系，本文繪制了各省份PM10、SO2、NO2的顯著性地圖和聚類(lèi)地圖（圖略）。在0.05的顯著性水平下，華北地區(qū)的北京、天津、山西、河北以及華東地區(qū)的山東、華中地區(qū)的河南這六個(gè)省份的PM10年平均濃度具有顯著性的空間正自相關(guān)，呈“高高集聚”的分布形態(tài)；而華南地區(qū)的廣東、廣西和西南地區(qū)的云南、貴州的PM10年平均濃度也具有顯著性的正自相關(guān)，但呈“低低集聚”的分布形態(tài)。其他省份的PM10年平均濃度的空間自相關(guān)性則不顯著，呈隨機(jī)分布的形態(tài)?？梢园l(fā)現(xiàn)，中國(guó)PM10的高濃度聚集區(qū)主要分布在華北地區(qū)以及華東、華中地區(qū)的個(gè)別省份，這幾個(gè)地區(qū)是中國(guó)人口較為密集和工業(yè)生產(chǎn)較為頻繁的地區(qū)，工業(yè)氣體廢棄物的排放和眾多人口冬季的采暖燃煤可能是造成這幾個(gè)地區(qū)PM10濃度高的主要原因。

對(duì)于SO2來(lái)說(shuō)，山東、河北、河南這三個(gè)省份的SO2年平均濃度具有顯著性的空間正相關(guān)關(guān)系，呈現(xiàn)“高高集聚”的空間布局。而北京和內(nèi)蒙古雖然都與高SO2的河北相鄰，但卻均呈現(xiàn)出“低高集聚”，對(duì)于北京來(lái)說(shuō)，這可能與政府近幾年的污染治理息息相關(guān)，SO2是與燃煤直接相關(guān)的大氣污染物。自1998年以來(lái)，在燃煤鍋爐清潔能源改造上，北京城六區(qū)累計(jì)完成約5.13萬(wàn)蒸噸燃煤鍋爐改造，基本實(shí)現(xiàn)了城六區(qū)無(wú)燃煤鍋爐，因此，北京的SO2年平均濃度才能在高SO2的相鄰區(qū)域下保持較低的濃度，也從側(cè)面反映出工業(yè)污染治理對(duì)環(huán)境帶來(lái)的巨大效應(yīng)。而對(duì)于內(nèi)蒙古來(lái)說(shuō)，其SO2年平均濃度之所以能保持較低的濃度，可能與其地理優(yōu)勢(shì)有關(guān)，內(nèi)蒙古北有大青山、東南部有蠻漢山、西南及南部則為平原，且由東北向西南逐漸傾斜，因此，形成了阻礙SO2擴(kuò)散的低氣壓帶，對(duì)河北的SO2的空間傳輸形成了較強(qiáng)的阻礙作用。而廣東屬于沿海省份，海陸風(fēng)對(duì)霧霾組分的擴(kuò)散和稀釋起了很大的作用，再加上自身污染沒(méi)有北方省份嚴(yán)重，因此，其SO2年平均濃度依然呈現(xiàn)“低低集聚”，即與周?chē)》菹嗨?，SO2年平均濃度都相對(duì)較低。

另外，北京、天津、河北三省的NO2年平均濃度呈現(xiàn)“高高集聚”的分布形態(tài)，說(shuō)明這三個(gè)省份也是NO2的高濃度地區(qū)，而廣東的NO2年平均濃度則相對(duì)較低，呈“低低集聚”分布。

2.4 基于LASSO方法的變量選擇

由全局空間自相關(guān)分析和局域空間自相關(guān)分析可以發(fā)現(xiàn)，中國(guó)各省份之間的霧霾組分確實(shí)存在空間效應(yīng)，因此，需要對(duì)其進(jìn)行空間回歸分析，以探討哪種回歸模型更適合于三種霧霾組分。在建模之前，先使用LASSO方法對(duì)變量進(jìn)行篩選，以使得模型更加簡(jiǎn)潔和更有解釋性。為了消除各變量量綱的影響，先對(duì)各變量取對(duì)數(shù)，其中，X1(第二產(chǎn)業(yè)增加值占GDP比重)、X3(城鎮(zhèn)人口比重)、X11(建成區(qū)綠化覆蓋率)、X15(生活垃圾無(wú)害化處理率)、X16(森林覆蓋率)由于是百分比的形式，因此不進(jìn)行對(duì)數(shù)化處理。然后，對(duì)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化、中心化處理。最后通過(guò)LASSO方法得到各組分的非零變量如下頁(yè)表3所示。

2.5 空間回歸分析

接下來(lái)，利用R語(yǔ)言的spdep安裝包，根據(jù)LASSO篩選出的非零變量對(duì)三種霧霾組分分別進(jìn)行空間回歸分析，以探討最優(yōu)的擬合模型。

2.5.1 PM10的空間回歸分析與檢驗(yàn)

建立PM10的回歸模型，并進(jìn)行診斷性檢驗(yàn)，結(jié)果如表4所示。結(jié)果表明，在lnPM10的普通最小二乘估計(jì)中，X1、lnX10、X16的回歸系數(shù)對(duì)應(yīng)的P值均小于0.05，均通過(guò)顯著性檢驗(yàn)。從拉格朗日乘數(shù)檢驗(yàn)結(jié)果來(lái)看，LMLAG對(duì)應(yīng)的P值為0.00214，遠(yuǎn)小于LMERR對(duì)應(yīng)的P值0.12300，而且R-LMLAR的P值為0.00194＜0.05，通過(guò)了顯著性檢驗(yàn)，而R-LMERR的P值為0.10980，在5%的顯著性水平下不顯著。因此，應(yīng)該選擇空間滯后模型較為合適。從空間滯后模型的擬合效果來(lái)看，擬合優(yōu)度為74.669%，比最小二乘估計(jì)的擬合優(yōu)度（60.470%）高，而且其對(duì)數(shù)似然值(-23.75857)比最小二乘估計(jì)的（-29.09351）大，AIC和SC也比最小二乘估計(jì)的小，因此，采用空間滯后模型是最優(yōu)的，其表達(dá)式如下：

表3 LASSO方法篩選出的非零變量

從模型的回歸系數(shù)來(lái)看，各變量依然對(duì)PM10具有顯著性的影響，并且其空間回歸系數(shù)ρ=0.582，對(duì)應(yīng)的P值為0.00001，也具有很強(qiáng)的顯著性，說(shuō)明中國(guó)各省份之間的PM10濃度具有很強(qiáng)的空間溢出效應(yīng)，即當(dāng)空間滯后值每增加1%時(shí)，本地的PM10濃度將平均上升0.582%。從社會(huì)影響因子來(lái)看，當(dāng)?shù)诙a(chǎn)業(yè)占GDP比重每上升1%時(shí)，PM10的年平均濃度將平均上升0.262%，這說(shuō)明第二產(chǎn)業(yè)的發(fā)展對(duì)PM10濃度帶來(lái)了重要的影響。第二產(chǎn)業(yè)主要分為工業(yè)和建筑業(yè)，對(duì)于工業(yè)來(lái)說(shuō)，其對(duì)煤炭、原油等化石燃料的燃燒是該行業(yè)產(chǎn)生PM10的主要原因。僅2013年，我國(guó)工業(yè)煤炭消費(fèi)總量占煤炭消費(fèi)總量的比重就已經(jīng)達(dá)到94.99%。因此，第二產(chǎn)業(yè)工業(yè)越發(fā)達(dá)的地區(qū)，其PM10將可能越高。而建筑業(yè)由于對(duì)鋼筋、水泥的需求比較大，因此，也間接地增加了化石燃料的燃燒，另外，建筑施工過(guò)程引起的揚(yáng)塵也是PM10的主要來(lái)源。從天然氣的供氣量（X10）來(lái)看，當(dāng)天然氣供氣總量每增加1%時(shí)，PM10的年平均濃度將平均上升0.296%，這結(jié)果可能與天然氣作為一種清潔能源相悖，但從目前我國(guó)天然氣的用途來(lái)看，其用途主要分為城市燃?xì)?、化工原料、天然氣發(fā)電、工業(yè)燃料四大塊，并且化工原料和工業(yè)燃料的用氣量占絕大部分。2013年，化工原料和工業(yè)燃料的天然氣用氣量占天然氣用氣總量的41%。因此，天然氣供氣總量很大一部分還是用于工業(yè)發(fā)展，天然氣供氣總量大的地區(qū)，其工業(yè)的比重也可能比較高，因此，工業(yè)污染對(duì)PM10的年平均濃度的影響可能也比較大。而森林覆蓋率（X16）與PM10的年平均濃度則呈負(fù)相關(guān)作用，當(dāng)森林覆蓋率每增加1%時(shí)，PM10的年平均濃度將平均下降0.308%，這主要是由于森林對(duì)PM10具有滯留、吸附和攔截功能，因此，森林覆蓋率越大，PM10的年平均濃度將越低。2.5.2 SO2、NO2的空間回歸分析與檢驗(yàn)

表4 PM10的回歸分析與診斷性檢驗(yàn)

同理，對(duì)SO2、NO2進(jìn)行空間回歸分析，得到下頁(yè)表5、表6。結(jié)果表明，SO2的空間回歸模型更適合空間滯后模型，其表達(dá)式如式（15）所示；而NO2則沒(méi)有通過(guò)拉格朗日乘數(shù)檢驗(yàn)，其回歸模型使用普通最小二乘估計(jì)即可，將不顯著變量剔除后，得到修正后的模型結(jié)果如式（16）、下頁(yè)表7所示。從SO2模型中各變量的回歸系數(shù)和顯著性檢驗(yàn)來(lái)看，在5%的顯著性水平下，各變量均通過(guò)了顯著性檢驗(yàn)。在空間區(qū)域影響方面，空間回歸系數(shù)ρ對(duì)應(yīng)的P值為0.00055，說(shuō)明我國(guó)各省份之間的SO2也存在空間溢出效應(yīng)，并且當(dāng)空間滯后值每增加1%時(shí)，本地的SO2年平均濃度將平均增加0.491%。從社會(huì)發(fā)展影響來(lái)看，SO2年平均濃度與第二產(chǎn)業(yè)增加值占GDP的比重呈正相關(guān)，當(dāng)比重每增加1%時(shí)，SO2年平均濃度將平均上升0.429%，這說(shuō)明第二產(chǎn)業(yè)增加值所占地區(qū)經(jīng)濟(jì)的比重越大，其SO2年平均濃度將可能越高，這主要還是由于第二產(chǎn)業(yè)中的工業(yè)廢氣的排放巨大。2014年，全國(guó)廢氣中SO2排放量為1974.4萬(wàn)噸，其中，工業(yè)SO2排放量為1740.4萬(wàn)噸，所占比重達(dá)到88.15%。因此，治理工業(yè)廢棄項(xiàng)目對(duì)霧霾治理至關(guān)重要。從治理廢棄項(xiàng)目完成投資(X18)的回歸系數(shù)可以發(fā)現(xiàn)，當(dāng)治理廢棄項(xiàng)目完成投資每增加1%時(shí)，SO2年平均濃度將平均下降0.270%。另外，森林也有降低SO2濃度的作用，當(dāng)森林覆蓋率每增加1%時(shí)，SO2年平均濃度將平均下降0.233%。這主要是由于植物對(duì)SO2有吸收作用，可對(duì)空氣進(jìn)行凈化。而對(duì)于NO2，X3（城鎮(zhèn)人口比重）和X21（民用汽車(chē)擁有量）與NO2年平均濃度呈正相關(guān)，城鎮(zhèn)人口比重每增加1%，NO2年平均濃度將平均增加0.365%，民用汽車(chē)擁有量每增加1%，NO2年平均濃度將平均增加0.503%。從這兩個(gè)指標(biāo)可以看出，我國(guó)各省份大氣中的NO2濃度主要源自城市中汽車(chē)尾氣的排放，當(dāng)城市人口比重越大，所擁有的民用汽車(chē)數(shù)量將越多，其大氣中的NO2濃度將越高。此外，與PM10和SO2相似，提高森林覆蓋率將是降低NO2濃度一個(gè)有效的措施。

表5 SO2的回歸分析與診斷性檢驗(yàn)

表6 NO2的回歸分析與診斷性檢驗(yàn)

表7 NO2修正后的回歸模型

3 結(jié)論

本文通過(guò)空間計(jì)量分析和LASSO變量選擇方法，對(duì)中國(guó)31個(gè)省份造成霧霾天氣的三種主要污染物——PM10、SO2和NO2進(jìn)行研究分析。結(jié)果發(fā)現(xiàn)：（1）華北地區(qū)及其周邊霧霾天氣最嚴(yán)重，呈“高高集聚”的空間分布。（2）第二產(chǎn)業(yè)是空氣中PM10和SO2的主要來(lái)源，且空間效應(yīng)顯著。第二產(chǎn)業(yè)中煤炭等化石燃料的燃燒和廢棄項(xiàng)目的排放是PM10和SO2的主要來(lái)源，通過(guò)治理工業(yè)廢棄項(xiàng)目完成投資的增加有助于SO2排放的減少。另外，各省份的PM10和SO2確實(shí)存在空間效應(yīng)，采用空間滯后模型最合適。（3）汽車(chē)尾氣排放是NO2的主要來(lái)源。隨著城市化進(jìn)程的加快，人們的生活水平不斷提高，越來(lái)越多的人擁有自己的私人汽車(chē)，而隨著汽車(chē)數(shù)量的急劇增加，汽車(chē)的尾氣排放成為NO2的主要來(lái)源，并且其回歸模型采用最小二乘估計(jì)比較合適。