陳優(yōu)芳，周一敏，趙昕奕

(1. 北京大學(xué)城市與環(huán)境學(xué)院，北京100871；2. 地表過程分析與模擬教育部重點(diǎn)實驗室，北京100871)

1 引 言

氣溶膠細(xì)顆粒物(空氣動力學(xué)直徑小于2.5 μm，PM2.5) 在霾的形成和發(fā)展中扮演主要角色[1]，由此造成的空氣污染導(dǎo)致了許多嚴(yán)重的環(huán)境問題，如城市能見度下降、影響交通安全、糧食產(chǎn)量下降等[2-4]。更嚴(yán)重的是，由于體積很小，它們可滲入人體支氣管和肺部，進(jìn)而引起嚴(yán)重的心血管疾病、呼吸道疾病甚至肺癌[5]。最新的研究已經(jīng)證實PM2.5會影響人體發(fā)病率和死亡率水平，對大眾健康有重要影響[6]，并在所有健康風(fēng)險影響因素中被列為第九大影響因素[7]。由于工業(yè)的快速發(fā)展、發(fā)達(dá)的交通網(wǎng)絡(luò)和密度較大的人口，華北地區(qū)是世界上大氣氣溶膠濃度最高的地區(qū)之一[8-9]。自1980年代以來，PM2.5污染天數(shù)顯著增加[10-12]，更重要的是，長期的重污染事件越來越頻繁。根據(jù)中華人民共和國環(huán)境保護(hù)部(MEP)發(fā)布的2013年度74個主要城市的年度空氣質(zhì)量數(shù)據(jù)，京津冀地區(qū)所有城市的PM2.5年平均濃度均超標(biāo)，其中7 個城市排名前10 位。一些城市的空氣污染天數(shù)占全年總天數(shù)的40%以上。

重污染事件給人類健康、氣候、食品生產(chǎn)、交通安全等帶來了一系列不利影響。因此，研究PM2.5的時空變化特征對進(jìn)一步預(yù)測PM2.5濃度和防治PM2.5污染具有重要意義?；诘孛嬲军c(diǎn)的觀測可提供準(zhǔn)確、及時的PM2.5濃度。但是，因為運(yùn)營成本和技術(shù)要求較高，地面站點(diǎn)的空間和時間覆蓋范圍有限，無法捕捉到高時空分辨率的PM2.5濃度。而衛(wèi)星遙感可提供長時間序列、大范圍的觀測數(shù)據(jù)而被廣泛應(yīng)用于各領(lǐng)域研究中[13-16]。在PM2.5濃度研究中，衛(wèi)星遙感可提供各種尺度的氣溶膠光學(xué)厚度(AOD) 產(chǎn)品，從而反演出未覆蓋地面PM2.5觀測區(qū)域的PM2.5濃度[17-21]。AOD 是氣溶膠最基本的光學(xué)特性之一，定義為介質(zhì)的消光系數(shù)在垂直方向上的積分，可用來描述大氣氣溶膠對光的衰減作用，推算大氣氣溶膠的含量、評估大氣污染的程度與確定氣溶膠的氣候效應(yīng)[22-23]。衛(wèi)星觀測AOD 始于1970s，目前已發(fā)展出包括MODIS、MISR、GEOS、POLDER、OMI、SeaWiFS、CALIOP等在內(nèi)的多數(shù)據(jù)來源、多種分辨率的衛(wèi)星AOD 數(shù)據(jù)，許多研究使用這些數(shù)據(jù)對氣溶膠濃度進(jìn)行反演[24]。在眾多衛(wèi)星AOD 產(chǎn)品中，中分辨率成像光譜儀(MODIS)AOD 產(chǎn)品由于其高分辨率和準(zhǔn)確度而被廣泛應(yīng)用于地面PM2.5濃度的反演[25]。反演效果被許多研究認(rèn)可，通過與地面觀測系統(tǒng)的比較，一般認(rèn)為MODIS AOD 的系統(tǒng)誤差為±0.05，反演算法帶來的誤差約為15%，達(dá)到設(shè)計目標(biāo)[9,26-27]。焦利民等[28]發(fā)現(xiàn)京津冀地區(qū)AOD 與PM2.5濃度具有很高的空間相關(guān)性，二者年均值相關(guān)系數(shù)約為0.6，但相關(guān)性存在時空差異，秦瑋等[29]對南京地區(qū)、范萌等[30]對珠三角地區(qū)的研究也有類似結(jié)論。可見，利用MODIS_AOD 對某區(qū)域PM2.5演變情況進(jìn)行研究是可行的。

反演方法主要包括簡單線性模型，土地利用回歸模型，地理加權(quán)回歸(GWR)和具有附加變量的復(fù)雜模型等[31-33]。這些方法中PM2.5與AOD 之間的關(guān)系是固定的，但是由于PM2.5濃度受排放條件、氣象條件等影響，而這些影響因素具有時空變化，因此PM2.5-AOD 關(guān)系隨著時空變化而變化[34]。因此，我們在反演時使用了混合效應(yīng)模型(Mixed effect Model)?；旌闲?yīng)模型是包含固定效應(yīng)和隨機(jī)效應(yīng)的統(tǒng)計模型。這些模型在對相同統(tǒng)計單位進(jìn)行重復(fù)測量的情況下特別有用。在1950年代，亨德森提供了固定效應(yīng)的最佳線性無偏估計(BLUE) 和隨機(jī)效應(yīng)的最佳線性無偏預(yù)測(BLUP)[35]。從那時起，混合建模已成為統(tǒng)計研究的一種主要手段。Lee 等[36]通過開發(fā)線性混合效應(yīng)(LME)模型將這種方法引入到PM2.5反演研究中。以前的研究表明，LME 模型的表現(xiàn)優(yōu)于其他統(tǒng)計模型，如地理加權(quán)模型[36]和線性回歸模型[37]。

目前，現(xiàn)有的研究主要是結(jié)合特定的污染過程中氣象要素和環(huán)流對PM2.5濃度的影響進(jìn)行研究，即對某次短時間內(nèi)的霾事件過程進(jìn)行研究[38-39]，但在長時間尺度的研究上，對于華北地區(qū)PM2.5濃度的時空分布特征研究較少，難以為PM2.5污染的預(yù)測和防治提供更有效的信息。因此，本研究利用了2013—2015年MODIS 的AOD 3 km 分辨率產(chǎn)品(MYDO4_3k)和環(huán)保部公布的2014年4月—2015年3月各個城市的PM2.5濃度數(shù)據(jù)對華北地區(qū)長時間序列的PM2.5濃度的時空分布特征進(jìn)行了分析。

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

2.1 研究區(qū)概況

本研究的研究區(qū)域為華北地區(qū)及其周邊小部分地區(qū)，包括北京、天津、河北、山東、陜西、河南、內(nèi)蒙古、遼寧等八個省區(qū)市，位于113～123 °E，34～43 °N 之間(圖1)。整個研究區(qū)域的地形呈西部高、東南低的分布，是大氣氣溶膠分布的重要成因之一。由于豐富的礦產(chǎn)資源，華北平原是我國的重要工業(yè)基地，其產(chǎn)業(yè)結(jié)構(gòu)也是形成空氣污染的重要原因之一。

圖1 研究區(qū)及測站城市示意圖

2.2 研究數(shù)據(jù)

2.2.1 衛(wèi)星遙感數(shù)據(jù)

衛(wèi)星遙感數(shù)據(jù)為NASA 發(fā)布的2003年1月—2015年3月的MODIS 傳感器MYD04(http://ladsweb.nascom.nasa.gov/data/search.html) 的 3 km 分辨率AOD 逐日數(shù)據(jù)產(chǎn)品。該衛(wèi)星每日經(jīng)過華北地區(qū)的時間為北京時間13 時，即該衛(wèi)星每日對華北地區(qū)的觀測時間為每日13 時。MYD04_3k數(shù)據(jù)包括了氣溶膠的光學(xué)性質(zhì)(例如光學(xué)厚度和粒徑分布)、質(zhì)量濃度等特征，具有高時空分辨率的特點(diǎn)。在MODIS 系列5 數(shù)據(jù)之前，氣溶膠數(shù)據(jù)一般使用MYD04 產(chǎn)品，但是其空間分辨率只有10 km，導(dǎo)致目前空氣污染的相關(guān)研究精度不夠。因此，在系列6 產(chǎn)品中，NASA 使用暗目標(biāo)氣溶膠算法創(chuàng)建了3 km 氣溶膠產(chǎn)品，MYD04_3k 與MYD04 兩個產(chǎn)品的算法一致，但是MYD04_3k空間分辨率更高。本研究將使用與地面觀測數(shù)據(jù)相匹配的2014年4—12月的數(shù)據(jù)進(jìn)行模型擬合，并用2015年1—3月的衛(wèi)星數(shù)據(jù)對模型效果進(jìn)行檢驗。

2.2.2 地面觀測數(shù)據(jù)

地面觀測數(shù)據(jù)為環(huán)保部公布的各地(城市)實時空氣質(zhì)量PM2.5濃度數(shù)據(jù)(http://air.ePMap.org/)。本研究使用了上述八個省區(qū)市中地處113～123°E，34～43 °N 范圍內(nèi)的41 個城市總計215 個站點(diǎn)、2014年4月—2015年3月的逐小時PM2.5數(shù)據(jù)。為了匹配衛(wèi)星觀測時間，取11：00—15：00 的小時數(shù)據(jù)的平均數(shù)據(jù)作為日值數(shù)據(jù)。根據(jù)《中國國家環(huán)境空氣質(zhì)量標(biāo)準(zhǔn)(GB3095-2012)》，本研究中的PM2.5污染指的是PM2.5質(zhì)量濃度大于75 μg/m3時的濃度，而污染日數(shù)比則是指在研究時段內(nèi)的PM2.5質(zhì)量濃度值大于75 μg/m3的天數(shù)總和與有PM2.5質(zhì)量濃度數(shù)據(jù)的天數(shù)總和的比值。

2.3 研究方法

2.3.1 聚類分析

MYD04_3k 的AOD 數(shù)據(jù)為日值數(shù)據(jù)，觀測時間為北京時間每天13 時左右。由于衛(wèi)星觀測的原因，MYD04_3k 的AOD 數(shù)據(jù)存在大量缺測值并且缺失情況每日不同，所以需要對缺失的觀測結(jié)果進(jìn)行補(bǔ)充?？紤]到氣溶膠在大氣中的空間一致性和連續(xù)性，本研究通過空間聚類分析將研究區(qū)域分為幾個不同區(qū)域，并且使用每一類中的有觀測數(shù)據(jù)的點(diǎn)的AOD 大小的平均值作為該區(qū)域無觀測點(diǎn)的替代；若是一個區(qū)域內(nèi)所有點(diǎn)在當(dāng)天均無觀測值，則視為該區(qū)域在當(dāng)天沒有觀測結(jié)果。

2.3.2 Mixed Effect Model 一元線性回歸算法反演PM2.5濃度

對華北地區(qū)進(jìn)行聚類分區(qū)后，本研究使用Mixed Effect Model 一元線性回歸算法，通過MYD04_3k 的AOD 數(shù)據(jù)反演了2003—2014年的華北地區(qū)各個小區(qū)域PM2.5濃度值。

首先，利用地面觀測的PM2.5濃度數(shù)據(jù)校準(zhǔn)MYD04_3k 的AOD 數(shù)據(jù)。AOD 數(shù)據(jù)取該城市中當(dāng)天所有點(diǎn)的AOD 值的平均值作為該城市當(dāng)天的AOD 值。對每個城市逐日的PM2.5濃度數(shù)據(jù)和AOD 數(shù)據(jù)進(jìn)行一元線性回歸。本研究既不是對每個城市分別進(jìn)行線性回歸，也不是對所有城市的所有天數(shù)統(tǒng)一進(jìn)行線性回歸，而是參考Lee 等[36]和Xie 等[40]的方法利用Mixed Effect Model 進(jìn)行線性回歸。Mixed effect Model 主要是在一元線性回歸中的固定參數(shù)的基礎(chǔ)上，對不同站點(diǎn)和不同時間的數(shù)據(jù)添加隨機(jī)參數(shù)，從而更精確地模擬PM2.5的濃度。本研究在反演過程中先將反演日期內(nèi)研究區(qū)中地面觀測數(shù)據(jù)和MODIS AOD 數(shù)據(jù)都有觀測數(shù)據(jù)的天數(shù)挑選出來，而后在SPSS 中通過分析-混合模型-線性的工具對數(shù)據(jù)進(jìn)行混合效應(yīng)模型擬合，最終獲得各個城市的參數(shù)。在Mixed Effect Model 中：

其中β 和α 是修正后的斜率和截距，εi為不同城市的修正值，對于不同城市來說，ε 不同，且∑εi=0。

在獲得了各個城市的一元線性回歸參數(shù)后，對ε 在整個研究區(qū)域進(jìn)行數(shù)值插值，并利用插值的結(jié)果反演出整個研究區(qū)域的2003—2014年的PM2.5濃度數(shù)值大小，結(jié)合聚類分析的結(jié)果最終得到了華北地區(qū)2003—2014年的PM2.5濃度。

2.3.3 EOF 分解

本文使用經(jīng)驗正交函數(shù)(EOF)對反演的PM2.5濃度時間和空間特征分析。EOF，也稱為特征向量分析或主成分分析(PCA)，是通過矩陣數(shù)據(jù)中的主要成分對數(shù)據(jù)結(jié)構(gòu)進(jìn)行分析的方法。 Lorenz[41]在1950年代首次將其引入氣象學(xué)和氣候?qū)W研究，現(xiàn)在廣泛應(yīng)用于地球科學(xué)和水電學(xué)等領(lǐng)域。EOF 分析能夠?qū)r變信號或數(shù)據(jù)集分解為空間函數(shù)部分和時間函數(shù)部分。空間函數(shù)部分不隨時間變化并且可總結(jié)數(shù)據(jù)的地理分布特征，而僅依賴于時間變化的時間函數(shù)部分由變量的線性組合組成，因此可展現(xiàn)數(shù)據(jù)的時間變化。前幾個分解成分占原始信號總方差的很大一部分，可看作原始信號的主要信息集中在這些成分上。該方法快速收斂，并且容易集中大量數(shù)據(jù)信息。因此，我們通過EOF分解可了解PM2.5的時間序列和空間模式。

EOF 分解過程主要如下。

(1) 將PM2.5反演數(shù)據(jù)預(yù)處理成距平形式，得到數(shù)據(jù)矩陣Xm×n，即m 個變量、n 次觀測。

(2) 計算X 矩陣與其轉(zhuǎn)置矩陣XT的交叉積，得到協(xié)方差矩陣，

(3) 計算方陣C 的特征根(λ1，征根λm)和特征向量Vm×m，二者滿足：

其中E 是m×m 維對角陣，即，

一般將特征根從大到小順序排列，即λ1＞λ2＞……＞λm。因為數(shù)據(jù)X 是PM2.5濃度，所以λ 應(yīng)該大于或者等于0。每個非0 的特征根對應(yīng)一列特征向量值，也稱為EOF。如λ1對應(yīng)的特征向量值稱為第一個EOF 模態(tài)，也就是V 的第一列，它反映的是PM2.5的空間分布特征。

(4) 計算主成分。將EOF 投影到原始資料矩陣上，就得到所有空間特征向量對應(yīng)的時間系數(shù)(即主成分)，即，

其中PC 中每行數(shù)據(jù)就是對應(yīng)每個特征向量的時間系數(shù)，也就是PM2.5的時間變化特征。

3 研究結(jié)果

3.1 對研究區(qū)進(jìn)行空間聚類分析

影響AOD 反演PM2.5效果的重要因素是AOD 數(shù)據(jù)的缺失。由于衛(wèi)星觀測數(shù)據(jù)大量缺失，必須對缺失的數(shù)據(jù)進(jìn)行補(bǔ)充，而補(bǔ)充數(shù)據(jù)的前提則是對研究區(qū)進(jìn)行空間聚類?？臻g聚類分析結(jié)果如圖2，表1 對9 個區(qū)域所在位置進(jìn)行了具體分析，結(jié)合EOF 分析的結(jié)果可看到，區(qū)域7 為污染情況最嚴(yán)重的區(qū)域，區(qū)域5、區(qū)域8 和區(qū)域6 的污染情況僅次于區(qū)域7；區(qū)域4 為污染情況快速減少的區(qū)域；區(qū)域9 是污染較少的區(qū)域；區(qū)域1、區(qū)域2、區(qū)域3 則是污染最少甚至可說幾乎沒有污染的區(qū)域。9 個小區(qū)域內(nèi)的各個觀察點(diǎn)PM2.5濃度的相關(guān)系數(shù)顯著高于區(qū)域外的各個點(diǎn)，且在該區(qū)域與其他區(qū)域的分界處相關(guān)系數(shù)快速下降，說明這幾個區(qū)域內(nèi)部一致性較高，區(qū)域聚類效果較好。

圖2 空間聚類分析結(jié)果行政區(qū)切割

表1 9 個區(qū)域所在位置

3.2 分區(qū)反演濃度

為了評估模型反演的擬合效果，計算了反演PM2.5濃度和地面觀測PM2.5濃度之間的確定系數(shù)(R2)，平均誤差(MPE)和均方根誤差(RMSE)。

R2如圖3 所示，大多數(shù)城市通過0.01 顯著性檢驗(P＜0.01)，威海、青島和德州通過0.10(P＜0.10)的顯著性檢驗，秦皇島、菏澤、大同和葫蘆島通過0.05 顯著性檢驗(P＜0.05)，但棗莊和鄭州未通過0.10 顯著性檢驗(P＞0.10)。 因此，我們使用了通過0.10 顯著性檢驗的39 個城市的數(shù)據(jù)來構(gòu)建混合效應(yīng)模型。總體而言，北方城市的相關(guān)系數(shù)略高于南方城市，尤其是河北南部、河南北部和山東南部城市的相關(guān)系數(shù)。這些城市通常是PM2.5污染嚴(yán)重的地區(qū)。也就是說，對于PM2.5污染嚴(yán)重的地區(qū)，模型擬合結(jié)果較差，而污染較少的地區(qū)則效果較好。

圖3 混合效應(yīng)模型中每個城市的擬合系數(shù)(R2)點(diǎn)的顏色越深，R2 越大。

最后對相關(guān)系數(shù)通過0.10 的顯著性檢驗的39個城市進(jìn)行Mixed Effect 一元線性回歸，使用2014年4—12月的地 面PM2.5和MODIS 衛(wèi)星AOD 作為模型的訓(xùn)練樣本，2015年1—3月的數(shù)據(jù)作為檢測樣本，對比發(fā)現(xiàn)經(jīng)過Mixed Effect Model 處理后，整個區(qū)域的地面觀測與衛(wèi)星觀測的數(shù)據(jù)的相關(guān)性顯著提高(擬合優(yōu)度提高50%左右)。2014年4—12月共有1 973 個數(shù)據(jù)，MPE 和RMSE 分別為21.6 和30.1。此外我們利用2015年1—3月8 個城市共141 個數(shù)據(jù)檢驗了模型效果，MPE 和RMSE 分別為25.6 和38.5(圖4)。

圖5 為各個城市經(jīng)過Mixed Effect 一元線性回歸算法后得到的截距的空間插值結(jié)果?？煽吹浇鼐嘣谡麄€研究區(qū)內(nèi)呈西南高東北低的趨勢，特別是河北南部和山東西部地區(qū)的截距明顯高于其他地區(qū)，表明在衛(wèi)星觀測到相同AOD 的情況下，河北南部和山東西部地區(qū)的PM2.5濃度要高于其他地區(qū)，這可能是由于不同地區(qū)的污染物化學(xué)類型不同造成的。

圖6 為對研究區(qū)反演后的PM2.5濃度進(jìn)行經(jīng)驗正交(EOF)分析得到的第一空間模態(tài)，其對整個PM2.5濃度場的方差貢獻(xiàn)率為73%，且第二～第六空間模態(tài)的方差貢獻(xiàn)率相當(dāng)并且都較少，故認(rèn)為第一空間模態(tài)可代表整個研究區(qū)PM2.5濃度的空間分布情況：河北南部、河南北部和山東西部地區(qū)為PM2.5污染(下稱污染)空間高值區(qū)域，從其向周圍擴(kuò)散一直到太行山脈和山東半島，整個區(qū)域的污染空間分布較均勻，梯度較小；其東部與山東半島交界的地區(qū)污染梯度較大，隨著經(jīng)度增加快速減少。河北中部到北京南部、天津北部，污染的等值線與地形海拔高度的等高線幾乎重合，在此地區(qū)污染情況隨海拔高度的增加迅速減少。在40 °N以北地區(qū)，污染情況總體較少，特別是張家口及其周邊地區(qū)幾乎沒有污染。

從整個EOF 分析的結(jié)果來看，地形對污染的影響非常明顯，整個山脈將污染分為兩個部分，山脈的東南平原地區(qū)為污染較嚴(yán)重的地區(qū)，而山脈西北則污染較少。其原因主要有兩方面，一是平原地區(qū)為主要的工業(yè)產(chǎn)地，特別是污染較大的重工業(yè)；二是由于地形的阻擋，平原地區(qū)產(chǎn)生的污染無法擴(kuò)散到山脈的西北地區(qū)，故而山脈西北地區(qū)的污染較少。

圖4 通過0.10 顯著性檢驗的39 個城市混合效應(yīng)模型模擬效果 a. 2014年4—12月的反演PM2.5 濃度和實測PM2.5 濃度擬合情況；b. 2015年1—3月反演PM2.5 濃度和實測PM2.5 濃度擬合情況，灰線表示擬合的線性回歸線。

圖5 研究區(qū)各個城市在Mixed Effect 算法中截距空間插值結(jié)果

對數(shù)據(jù)進(jìn)行反演后，利用空間聚類的結(jié)果對9個研究分區(qū)的衛(wèi)星觀測數(shù)據(jù)進(jìn)行補(bǔ)充，使用每一類中有觀測數(shù)據(jù)的格點(diǎn)的反演PM2.5濃度的平均值作為該區(qū)域無觀測數(shù)據(jù)格點(diǎn)的替代濃度；若是一個區(qū)域內(nèi)所有點(diǎn)在當(dāng)天均無觀測值，則視為該區(qū)域在當(dāng)天沒有觀測結(jié)果。在經(jīng)過數(shù)據(jù)補(bǔ)充之后，在研究的總共4 808日中，研究區(qū)內(nèi)9 個區(qū)域數(shù)據(jù)齊全的日數(shù)共有2 017 天，8 個區(qū)域有數(shù)據(jù)的有740 天。3.3 研究區(qū)PM2.5濃度在2003—2014年的時空分布

圖6 研究區(qū)EOF 第一空間模態(tài)

在補(bǔ)充了部分缺失數(shù)據(jù)后，通過修正過的衛(wèi)星反演的PM2.5濃度得到華北地區(qū)PM2.5濃度在2003—2014年的變化情況。用每個區(qū)域中的平均值代表這個區(qū)域的污染情況得到不同區(qū)域污染日數(shù)百分比的年際變化(圖7)，通過該圖也可看出研究區(qū)污染情況的時空分布特征。從空間分布看，總體情況與EOF 第一主成分的空間格局類似，太行山脈是污染強(qiáng)弱的明顯的分界線，山脈東南部的污染顯著高于西北部；整個研究區(qū)北部和山東半島地區(qū)污染情況較少；河北南部和山東西部地區(qū)是污染最嚴(yán)重的區(qū)域，其周圍區(qū)域污染程度僅次于這一區(qū)域，南北方向上的污染程度稍高于東西方向。從區(qū)域來看，北方4 個區(qū)域(區(qū)域1～區(qū)域4)污染情況相似，幾乎無污染，其中區(qū)域1 最好，區(qū)域4 最差；區(qū)域7 為污染情況最嚴(yán)重的地區(qū)，區(qū)域5 和8 次之，區(qū)域6 再次之。區(qū)域與區(qū)域間的相關(guān)系數(shù)也可以反映出這種空間聯(lián)系性，從表2 可發(fā)現(xiàn)，區(qū)域1、區(qū)域2、區(qū)域3 和區(qū)域4 之間都有較高的相關(guān)系數(shù)，區(qū)域5、區(qū)域6 和區(qū)域7 之間的相關(guān)系數(shù)較高(下劃線數(shù)據(jù))，區(qū)域6、區(qū)域7 和區(qū)域8之間相關(guān)系數(shù)較高(下劃線數(shù)據(jù))；區(qū)域9 污染較少，但與其他8個區(qū)域的相關(guān)性均較低，可能與其獨(dú)特的地理位置有關(guān)。

圖7 整個研究區(qū)和9 個區(qū)域污染日數(shù)百分比的年際變化 藍(lán)色系曲線為污染較少的區(qū)域，與左側(cè)坐標(biāo)相對應(yīng)；紅色系曲線為污染較多的區(qū)域，與右側(cè)坐標(biāo)相對應(yīng)；黑色曲線為整個研究的污染日數(shù)比，與左側(cè)坐標(biāo)相對應(yīng)。

從時間分布來看，2004年和2005年兩年為污染情況較少的年份，從2006年開始，污染程度急劇加強(qiáng)，并且在隨后兩年保持高值，但是2009年稍稍緩和，可能與保障2008年奧運(yùn)會順利召開所采取的大氣污染治理措施有關(guān)。之后污染再次加劇，不過2013年后，污染呈現(xiàn)下降趨勢。

表2 不同區(qū)域污染情況的相關(guān)系數(shù) 均通過0.01 顯著性檢驗。

4 討論與展望

本研究使用了MODIS 的AOD 數(shù)據(jù)和PM2.5的地面觀測數(shù)據(jù)進(jìn)行統(tǒng)計分析。此前，大部分研究主要使用單站或者多個站點(diǎn)的觀測結(jié)果進(jìn)行研究，導(dǎo)致研究結(jié)果在空間和時間上的局限性。本研究使用了上述兩套數(shù)據(jù)，可以提供大范圍和長時間的氣溶膠觀測數(shù)據(jù)，大大提高了研究的可靠性。但是由于衛(wèi)星遙感的AOD 有數(shù)據(jù)缺失的問題，本研究中采用了聚類分區(qū)，然后用區(qū)域平均的結(jié)果對區(qū)域中的缺失數(shù)據(jù)進(jìn)行填補(bǔ)的方法，引入了誤差。雖然兩者的相關(guān)性較高，但需要探討是否有其他更好的方法可對缺失數(shù)據(jù)進(jìn)行填補(bǔ)。

此外，從長期或平均的狀態(tài)來說，城市污染程度取決于城市的能源結(jié)構(gòu)、交通和工業(yè)排放污染物的多少[42-43]，但從短期或?qū)崟r狀態(tài)而言，卻主要與當(dāng)?shù)?、?dāng)時的氣象條件有關(guān)。大尺度環(huán)流形勢和局地氣象條件不僅可影響污染物的生成、積累和清除，同時也是影響區(qū)域輸送的重要外部條件[44-45]。在華北地區(qū)，污染物濃度的變化很大程度上取決于局地的氣象條件，如風(fēng)速、邊界層高度、相對濕度、降水等[46-47]，因而在未來需要著重對氣象要素如何影響PM2.5濃度時空分布特征進(jìn)行研究。

5 結(jié) 論

本研究利用了衛(wèi)星遙感氣溶膠光學(xué)厚度資料、地面觀測的PM2.5濃度數(shù)據(jù)，采用了混合效應(yīng)模型、聚類分析、EOF 分解等統(tǒng)計分析方法，研究了2003—2014年華北及周邊地區(qū)PM2.5濃度的時空分布特征。

(1) 衛(wèi)星遙感的AOD 與地面觀測PM2.5濃度值有較高的相關(guān)系數(shù)，可利用MODIS 衛(wèi)星遙感的AOD 數(shù)據(jù)對地面觀測的PM2.5濃度進(jìn)行反演。

(2) 在華北和周邊地區(qū)，PM2.5濃度呈現(xiàn)出明顯的空間分布特征，40 °N 以北的區(qū)域和山東半島地區(qū)污染情況較少，太行山脈是污染強(qiáng)弱的明顯的分界線，山脈東南部的污染顯著高于西北部，且在地勢變化的地方出現(xiàn)明顯的突變；河北南部、河南北部和山東西北小部分區(qū)域是污染最嚴(yán)重的地區(qū)，其周圍區(qū)域污染程度僅次于這一區(qū)域，南北方向的污染程度稍高于東西方向。

(3) 從時間分布來看，2004 和2005年兩年為污染情況較少的年份，從2006年開始，污染程度急劇加強(qiáng)，2009年稍稍緩和，但從2010年開始又逐漸加重，2013年后呈現(xiàn)下降趨勢。河北南部和山東西北小部分區(qū)域的污染情況與北京南部、天津和河北中部的污染情況相關(guān)性較高。