基于融合MODIS AOD 和LME 模型的京津冀地區PM2.5 季節濃度模擬

邵絲露，孫艷玲，馬振興，陳 莉，毛 健，高 爽，張 輝，付宏臣，景 悅

（天津師范大學地理與環境科學學院，天津 300387）

經濟的迅猛發展和城市化進程的加快在改善人類生活質量的同時造成了能源的大量消耗，導致環境污染日益嚴重，其中近年來頻發的霾備受關注.PM2.5又稱細顆粒物，其空氣動力學直徑≤2.5 μm，是霾的重要組分之一.有研究證實PM2.5不僅可以影響空氣質量，還會造成光的散射和吸收，從而降低能見度[1]，過度吸入會造成心、腦血管等方面的疾病，甚至威脅生命[2-3].因此，了解PM2.5濃度分布對治霾措施的選擇以及心、腦血管疾病的病因分析均具有非常重要的作用.

時間序列的PM2.5濃度可由地面監測站點提供，但目前監測站點并非均勻分布，無法得到空間連續的監測值，對有限的站點數據進行空間插值又會產生較大誤差.因此，利用衛星遙感（覆蓋范圍廣、分辨率高）獲取空間連續的氣溶膠光學厚度（aerosol optical depth，AOD），進而反演PM2.5濃度成為國內外學者常用的方法.AOD 是氣溶膠的消光系數在垂直方向的積分，可用于表征大氣的渾濁程度. AOD 產品包括MODIS AOD、MISR AOD、MAIAC AOD、風云系列極軌衛星產品和環境一號衛星產品等，其中MODIS AOD 是應用最為普遍的產品之一.利用AOD 反演PM2.5濃度的方法主要為建立AOD-PM2.5關系模型，該模型經歷了由簡單線性回歸到精度較高的土地利用回歸以及地理加權回歸等模型的轉變. 早期研究假定AOD-PM2.5的關系是恒定不變的，而Lee 等[4]首次提出的線性混合效應模型則解釋了其隨時間變化的情況，研究模擬所得美國東北部PM2.5濃度的可決系數R2為0.97，極大地提高了擬合精度，證明該模型對AOD-PM2.5具有較高的擬合能力. 隨后Xie 等[5]、Ma 等[6-7]、Zheng 等[8]和Yang 等[9]分別利用改進的線性混合效應模型模擬中國不同地區的PM2.5濃度（R2=0.71 ～0.85），均得到較高的擬合精度.

MODIS AOD 產品包括Terra 衛星的MOD04_L2和MOD04_3K 產品以及Aqua 衛星的MYD04_L2 和MYD04_3K 產品共四類.每一類AOD 產品由于云層、冰雪覆蓋等原因都會或多或少出現缺失，而數據量的不足會影響PM2.5濃度的模擬.因此，對多種產品進行融合，即使用多星融合數據可以相互彌補缺失數據，提高模擬準確度.Zheng 等[8]將MOD04_L2 和MYD04_L2兩類數據建立線性回歸來預測缺失的AOD，反演得到精度較高的京津冀、長江三角洲和珠江三角洲的PM2.5濃度（R2=0.77 ～0.80）.計羽西等[10]利用融合4 種氣溶膠產品所得AOD 日均值對湖北省PM2.5濃度進行模擬，同樣得到較好的擬合結果（R2=0.85 ～0.89）.

京津冀是霾較為嚴重的地區，現有研究大多通過地理加權回歸模型[11-12]和土地利用回歸模型等[13-14]模擬該地區PM2.5濃度，均獲得較好的擬合結果. 近些年，有研究將線性混合效應模型應用于京津冀地區的分析，郝靜等[15]和景悅等[16]均利用該模型模擬了京津冀地區PM2.5濃度，但由于采用單顆星AOD，數據缺失較多，研究僅能分冷暖兩季進行.本文同樣以京津冀為研究區域，對MODIS AOD 的四類產品進行融合，增加了AOD 的覆蓋范圍，使得數據量可以按季節展開研究.因此，本研究基于融合的AOD 數據和PM2.5站點觀測數據，按季節分別建立線性混合效應模型并估算PM2.5濃度、模擬不同季節PM2.5濃度的空間分布，以期為PM2.5濃度分布相關研究提供新的數據處理思路，也為流行病學的病源追溯研究提供更加全面的數據基礎.

1 數據來源與方法

1.1 研究區域概況

京津冀包括北京市、天津市以及河北省的保定、唐山、廊坊、石家莊、邯鄲、秦皇島、張家口、承德、滄州、邢臺和衡水等11 個地級市，是我國重要的工業中心之一，也是政治中心和文化中心.該地區屬暖溫帶大陸性季風氣候，冬季盛行西北風，夏季盛行東南風，四季分明，夏季炎熱多雨，冬季寒冷少雪，春多風沙，秋高氣爽.地勢由西北向東南傾斜，西北部為壩上高原、山地和丘陵，東南部為廣闊平原，該地區的季風氣候和地形皆不利于大氣污染物的擴散，導致大氣污染嚴重.

1.2 數據來源與處理

1.2.1 MODIS AOD 數據

用于融合的原始AOD 是由美國地球觀測系統（EOS）中Terra 和Aqua 衛星搭載的中分辨率成像光譜儀（MODIS）提供的產品，該產品與大氣顆粒物濃度有較好的相關性且對外免費. 上午軌道衛星Terra 過境時間約為地方時10 ∶30，下午軌道衛星Aqua 過境時間約為地方時13 ∶30，二者在時間上互補. MODIS AOD 產品包括10 km 和3 km 空間分辨率，分別由深藍算法和暗像元算法反演得到，深藍算法對城市和沙漠等明亮地區反演效果較好，暗像元算法在植被濃密區反演效果較好，二者融合可以彌補云層和冰雪等覆蓋造成的數據缺失.

本研究采用從NASA 網站（https：//ladsweb.modaps.eosdis.nasa.gov/search/）下載的C6.1 版本2016 年1 月1日—2016 年12 月31 日的MOD04_L2、MOD04_3K、MYD04_L2 和MYD04_3K 四類數據.首先利用ENVI5.3對其進行預處理，包括幾何校正和拼接，然后對拼接后的數據進行融合.較高的空間分辨率有助于提高反演結果的質量[14]，因此將空間分辨率為10 km 的MOD04_L2 和MYD04_L2 的數據分別重采樣為3 km 分辨率的數據，并對重采樣后的數據建立線性回歸來填補缺失，而后對填補后的兩類數據取均值得到融合后的3 km AOD（深藍算法）.MOD04_3 K 和MYD04_3K 同樣建立線性回歸填補缺失并取均值得到融合后的3 km AOD（暗像元算法）.最后基于重采樣后分辨率為3 km的土地覆蓋類型數據為不同土地覆蓋類型匹配相應的融合AOD 數據，得到空間分辨率為3 km 的優化后的日均AOD，用于PM2.5濃度的反演.

1.2.2 PM2.5數據

PM2.5監測數據來源于中國環境監測總站，時間分辨率為1 h.由于融合AOD 為日均值，PM2.5數據也由每小時值計算得到日均值，即2016 年1 月1 日—2016 年12 月31 日日均PM2.5數據.京津冀地區監測站點分布如圖1 所示.由圖1 可以看出，京津冀共有國家空氣質量監測站點80 個，其中北京市12 個，天津市15 個，河北省53 個，分別位于保定、滄州、承德、邯鄲、衡水、廊坊、秦皇島、石家莊、唐山、邢臺和張家口.

圖1 京津冀地區國家空氣質量監測點分布示意圖Fig.1 Distribution diagram of national air quality monitoring points in Beijing-Tianjin-Hebei region

1.3 模型擬合與驗證

1.3.1 線性混合效應模型

AOD 與PM2.5的關系受相對濕度、溫度、邊界層高度和地形等因素的影響，其中相對濕度、溫度和邊界層高度等皆逐日變化，因此AOD 與PM2.5的關系也呈現逐日變化的特點，而簡單的模型只能解釋AOD-PM2.5間的固定效應，無法表示逐日隨機效應.已有研究表明線性混合效應模型可以較好地闡釋這種逐日變化的隨機效應[4-10]，因此本研究選擇利用線性混合效應模型模擬京津冀地區的PM2.5濃度，模型表達式為

式（1）中：PPM2.5和AODij分別為監測站點i 在第j 天的顆粒物濃度和AOD 值；α 和β 分別為固定截距和固定斜率，其中β 代表AOD 對PM2.5濃度的固定影響，且不隨時間發生變化；μj和γj分別為隨機截距和隨機斜率，其中γj代表第j 天AOD 對PM2.5濃度的隨機影響，反映了時間上的差異；εj為第j 天的隨機殘差；Σ 為逐日變化的隨機效應對應的方差-協方差矩陣.此模型的建立基于R 語言完成.

1.3.2 模型驗證

為檢驗模型擬合是否存在過度擬合，采用留一交叉驗證方法進行檢驗.在AOD-PM2.5匹配數據對構成的樣本集中選取每一個樣本輪流作為測試集，其余全部樣本作為訓練集訓練模型，用測試集檢驗模型，不斷重復訓練模型、檢驗模型的過程，直至所有樣本均做過一次測試集，最終取所有交叉驗證結果的均值作為模型擬合的評價結果.

線性混合效應模型的擬合準確度通過PM2.5監測值與預測值間的可決系數R2和均方根誤差（RMSE）來判斷，其中

式（2）中：yi為第i 天PM2.5濃度的監測值為第i 天PM2.5濃度的預測值；n 為樣本數.

2 結果與分析

2.1 描述性統計

融合AOD 數據與PM2.5匹配過程中，剔除當天數據對少于2 對的數據，AOD-PM2.5數據最終全年可匹配天數共317 d，春夏秋冬各79 d、76 d、80 d 和82 d.與使用單顆星的AOD 數據相比[16]，融合AOD 與PM2.5的可匹配天數和樣本量明顯增多，尤其是AOD 數據嚴重缺乏的冬季.單顆星的AOD 數據由于缺失較多，建立模型時時間序列只能具體到冷暖季，而融合后的AOD 彌補了大量缺失的數據，因此本研究具有充足的樣本量，可以將所有數據分為四季，分別建立線性混合效應模型模擬不同季節PM2.5濃度分布.

AOD-PM2.5匹配數據的描述性統計結果如表1 所示.由表1 可知，2016 年PM2.5濃度變化范圍在3.000 ～688.690 μg/m3，年均值約為67.188 μg/m3，超過年均濃度一級標準，但未超過二級標準.四季AOD-PM2.5樣本量均在2000 對以上，可以建立線性混合效應模型.PM2.5濃度均值冬季最高，為89.131 μg/m3，秋季次之，夏季最低，為42.287 μg/m3，此季節差異可能與冬季集中供暖有關.融合AOD 年均值約為0.241，全年AOD的數值變化范圍為0.002～3.232.四季均值春季最高，為0.232，秋冬次之，夏季最低，為0.189.

表1 AOD-PM2.5 匹配數據的描述性統計Tab.1 Descriptive statistical of AOD-PM2.5 matching data

匹配京津冀地區日均PM2.5濃度和日均AOD，結果如圖2 所示. 由圖2 可以看出，日均AOD 與日均PM2.5濃度的變化趨勢基本一致，尤其是夏季.PM2.5濃度有237 d 超過日均濃度的一級標準（35 μg/m3），108 d超過日均濃度的二級標準（75 μg/m3），說明2016 年京津冀地區PM2.5濃度嚴重超標，霾污染問題嚴峻.

圖2 京津冀地區的日均PM2.5 和日均AOD 變化曲線圖Fig.2 Daily average change graph of PM2.5 and AOD of Beijing-Tianjin-Hebei region

2.2 模型擬合與驗證結果

本研究基于融合AOD 數據和PM2.5數據分別建立一元線性回歸模型和線性混合效應（LME）模型，二者所得結果如表2 所示，2 種模型固定截距和固定斜率皆有統計學意義（p ＜0.0001）.

表2 不同模型擬合結果Tab.2 Fitting results of different models

通過運行R 語言程序，得到PM2.5濃度各季節線性混合效應模型每日的隨機斜率和隨機截距.PM2.5濃度春季線性混合效應模型的隨機斜率和隨機截距變化范圍分別為-171.074 ～532.461 和-49.054 ～94.474，標準差分別為103.658 和26.234；夏季的隨機斜率和隨機截距變化范圍分別為-199.674 ～82.032 和-39.368 ～42.330，標準差分別為57.052 和21.903；秋季的隨機斜率和隨機截距變化范圍分別為-246.985 ～444.934 和-68.509～150.466，標準差分別為134.600 和35.511；冬季的隨機斜率和隨機截距變化范圍分別為-576.406 ～1151.724 和-66.886 ～458.905，標準差分別為302.570和73.309.由表2 可知，四季一元線性回歸模型擬合的R2分別為0.105、0.127、0.198 和0.108，而線性混合效應模型擬合的R2分別為0.738、0.668、0.644 和0.760.對于簡單一元線性回歸模型而言，PM2.5濃度較低的季節的擬合效果優于PM2.5濃度較高的季節的擬合效果，而線性混合效應模型對PM2.5濃度較高的季節擬合效果顯著.總體上，線性混合效應模型對AOD 和PM2.5濃度的擬合效果明顯優于簡單一元線性回歸模型.

為判斷模型是否存在過度擬合，采用留一交叉驗證進行檢驗，結果如圖3 所示，其中圖3（a）、圖3（c）、圖3（e）和圖3（g）分別為春夏秋冬四季的PM2.5濃度實測值和線性混合效應模型預測值的散點圖，圖3（b）、圖3（d）、圖3（f）和圖3（h）分別為四季的留一交叉驗證結果.對比圖3 結果可知，各季節驗證后的散點圖較驗證前更為分散，表明各季節的模型均存在過度擬合.

圖3 四季線性混合效應模型結果和留一交叉驗證結果Fig.3 Results of seasons linear mixed effect model and leave-one-out cross-validation

由圖3 可以看出，進行留一交叉驗證后，R2均有所下降，而RMSE 均有所上升，說明模型僅存在輕微過度擬合.春季和冬季R2的擬合效果較好，RMSE 冬季最高，這可能與冬季PM2.5濃度過高以及AOD 數據的融合有關.與其他應用該模型的多變量研究相比，本研究只使用單一變量，即融合后的AOD，擬合結果與他人研究[5-9]（R2=0.71 ～0.85）比較接近，且過擬合程度較低，說明使用高分辨率的融合AOD 數據比增加變量更容易提高模型性能.

對于京津冀整體而言，基于監測站點實測所得PM2.5四季平均濃度分別為57.197、42.287、72.745 和89.131 μg/m3，線性混合效應模型預測所得PM2.5四季平均濃度分別為57.277、42.487、72.740 和88.952 μg/m3.季節平均濃度實測值與預測值非常接近，說明模型比較適用于京津冀地區的顆粒物濃度模擬.

京津冀地區共有13 個城市，各城市不同季節PM2.5的實際濃度平均值和預測濃度平均值如圖4 所示.

圖4 不同季節各城市的PM2.5 平均監測濃度和平均模擬濃度Fig.4 PM2.5 average monitored concentration and average simulated concentration of each city in different seasons

由圖4 可以看出，大多城市實測值與預測值比較接近，證明模型擬合較好.各季節實際濃度平均值最低的城市包括張家口市、承德市和秦皇島市，這些城市人口密度小，工業占總產業比重較低，以旅游業為主，PM2.5排放量較低；此外，這些城市植被覆蓋度較高、對PM2.5的吸附作用強也是原因之一.衡水市、保定市和石家莊市各季節實際濃度平均值較高，與人口密度大、工業企業數量多、尤其是高能耗企業和污染企業較多等因素有關，上述因素均導致PM2.5排放量較高.春夏季承德市、衡水市和張家口市實測值與預測值差異顯著，秋冬季保定市、承德市、石家莊市和張家口市實際值與預測值差異顯著，導致該現象的原因一方面是以上城市均為PM2.5實際濃度最低或最高的城市，而過低或過高的顆粒物濃度在一定程度上會降低線性混合效應模型擬合的準確度；另一方面，以上城市均為城市面積較大而環境監測站點數量較少的城市，通過有限的站點數據模擬區域濃度誤差較大.

各城市不同季節PM2.5濃度的擬合結果如表3 所示.

表3 不同季節各城市PM2.5 擬合結果對比Tab.3 Comparison of PM2.5 fitting results of different cities in different seasons

由表3 可知，各城市春季PM2.5實測值與線性混合效應模型預測值的R2為0.406 ～0.903，春季RMSE 為14.817 ～30.097 μg/m3，平均值分別為0.718 μg/m3和20.294 μg/m3；夏季R2為0.397～0.825，RMSE 為11.666～21.334μg/m3，平均值分別為0.666μg/m3和15.290μg/m3；秋季R2為0.394～0.819，RMSE 為21.569～36.071 μg/m3，平均值分別為0.677 μg/m3和29.256 μg/m3；冬季R2為0.433～0.875，RMSE 為27.484 μg/m3～57.922 μg/m3，平均值分別為0.720 μg/m3和39.389 μg/m3.總體上，各城市均為春季和冬季的擬合效果較好，說明線性混合效應模型應用于PM2.5濃度模擬時，在濃度較高的季節表現更優. 由于監測站點數量相對較多、建模數據缺失較少以及大氣污染略重等原因，廊坊市、北京市和天津市在模擬不同季節PM2.5濃度時結果均較為理想.除張家口市以外，其余城市不同季節PM2.5濃度擬合的R2均大于0.5，大部分城市擬合效果比較理想.由于冬季AOD 數據通過融合方式補充較多，故各城市冬季RMSE 較高，絕大部分城市RMSE 在30 μg/m3以上.

圖5 2016 年京津冀地區PM2.5 季節濃度模擬Fig.5 Seasonal concentration simulation of PM2.5 of Beijing-Tianjin-Hebei region in 2016

2.3 PM2.5 季節濃度分布模擬

基于線性混合效應模型得到各個季節京津冀地區的逐日PM2.5濃度分布，利用ArcGIS10.2 對各季節的逐日PM2.5濃度分別進行加和平均，得到京津冀地區PM2.5各季節平均濃度分布如圖5 所示，圖5（d）中，由于冬季受冰雪覆蓋影響，北部山區一些地區數據缺失.

由圖5 可以看出，PM2.5季節濃度總體上冬季最高，春秋次之，夏季最低.春季PM2.5濃度高值區主要集中在京津冀中部和南部，邯鄲、唐山以及保定東南部PM2.5濃度較高，分別為63.922、62.969 和62.397 μg/m3；承德和張家口PM2.5濃度較低，分別為49.794 μg/m3和50.735 μg/m3.夏季受以東南風為主的夏季風影響，顆粒物向西北擴散，加之供暖期已過，污染物排放量減少，故PM2.5濃度明顯下降且各城市差異較小，衡水、保定和廊坊PM2.5濃度較高，分別為47.336、45.821 和45.795 μg/m3；秦皇島和天津PM2.5濃度較低，分別為33.117 μg/m3和39.305 μg/m3. 秋季各城市PM2.5濃度普遍上升，高值集中在京津冀中南部地區，保定、石家莊和邢臺PM2.5濃度較高，分別為84.988、80.600 和80.367 μg/m3；承德和秦皇島PM2.5濃度較低，分別為55.675 μg/m3和56.282 μg/m3.冬季進入供暖期，污染物排放量增多，PM2.5濃度進一步提高，受冬季西北季風的影響，顆粒物集中分布于東南部城市，西北部城市PM2.5濃度顯著下降，邯鄲、保定、滄州、衡水和石家莊一帶PM2.5濃度較高，分別為107.096、104.656、101.908、101.667 和101.055 μg/m3；秦皇島和承德PM2.5濃度較低，分別為62.633 μg/m3和65.765 μg/m3.圖4 和圖5均表明京津冀地區的PM2.5濃度存在明顯的季節差異以及地區差異.本文對京津冀地區PM2.5濃度分布的模擬結果與孟曉艷等[17]、溫佳薇等[18]對該區域PM2.5濃度變化特征的研究結果相近，說明本研究所得PM2.5濃度分布的模擬結果合理.

此外，由圖5 還可知，四季PM2.5濃度均為北部低、南部高.太行山和燕山的阻擋使南部的PM2.5難以向北擴散，且張家口和承德兩市以旅游業為主，污染物排放少，因此張家口市和承德市大部分區域PM2.5濃度較低，這一現象與圖4 所得結果相吻合. 渤海沿岸PM2.5濃度也較低，這一方面是由于沿海風力較大，有利于沿海地區的污染物擴散；另一方面沿海氣象條件復雜，所用數據和模型應用于沿海地區可能存在誤差. PM2.5濃度較高的地區主要為石家莊、保定、廊坊和唐山等城市，這些城市均以第二產業為主，污染物排放量較高.同時，私家車保有量逐年增加、冬季大量燃煤供暖、建筑揚塵以及二次顆粒物等都是導致PM2.5濃度高居不下的原因.

3 討論

基于融合AOD 數據，本文利用線性混合效應模型模擬的京津冀地區2016 年PM2.5季節平均濃度的結果（各季節交叉驗證的R2分別為0.733、0.658、0.633、0.756）與其他使用同種模型的國內研究[6-10，15-16，19-20]結果相近，但低于Lee 等[4]的研究結果，這可能是因為本研究沒有對數據進行篩選處理以及研究區域PM2.5濃度過高、AOD 數據融合產生誤差等原因所致.與應用廣義相加模型、BP 神經網絡、結合MODIS AOD 的LUR 時空模型等模型或方法的研究[11，14，21-23]相比，本研究所得結果與較高級統計模型結果相近，但仍低于一些應用改進的高級模型的研究[12，24]，如混合時空地理加權回歸模型和地理加權回歸空間降尺度方法等.雖然使用改進的融合AOD 數據提高了模擬精度，但本研究只使用AOD 一個變量模擬PM2.5濃度，而應用高級模型的研究基本都使用土地利用、氣象因素、人口密度和道路等多變量進行模擬，結果證明應用較多變量在一定程度上可以提高模型擬合精度，這也是本研究今后的改進方向之一.此外，也可以考慮對AOD 數據進行濕度訂正和垂直訂正，以降低數據給模型帶來的誤差，從而獲得更優的研究結果.

4 結論

本研究對AOD 四種產品的數據進行融合，基于線性混合效應模型對京津冀地區全年4 個季節分別建立模型，估算并模擬PM2.5濃度的分布情況，得到以下結論：

（1）京津冀地區PM2.5濃度春季交叉驗證后R2為0.733，RMSE 為21.040 μg/m3；夏季R2為0.658，RMSE 為16.201μg/m3；秋季R2為0.633，RMSE 為33.717μg/m3；冬季R2為0.756，RMSE 為38.402 μg/m3.和驗證前相比，各季節的模型僅存在輕微過度擬合，且擬合結果較為理想.

（2）除承德、衡水、張家口和石家莊外，各城市四季實際值與預測值接近，且各城市均為春冬季的擬合效果較好，這表明線性混合效應模型應用于PM2.5濃度模擬時，在濃度較高的季節表現更優.

（3）PM2.5四季濃度的分布情況大致為北部低、南部高，其中張家口市和承德市PM2.5濃度較低，而石家莊、保定、廊坊和唐山等城市PM2.5濃度較高；季節濃度總體上冬季最高，春秋次之，夏季最低.以上結果說明在氣候、地形、經濟和人口等多種因素的綜合影響下，京津冀地區的PM2.5濃度存在季節差異和地區差異.