包頭市細顆粒物遙感監測混合線性模型

同麗嘎，李雪銘，張靖，3

(1.包頭師范學院 資源與環境學院，內蒙古 包頭 014030；2.遼寧師范大學 城市與環境學院，遼寧 大連 116029；3.大連民族大學 環境與資源學院，遼寧 大連 116600)

0 引言

隨著我國城市化進程的加快，建筑物逐漸取代城市及周邊的植被，城市建設、工業生產等帶來一系列的生態環境問題[1]。其中大氣污染日益嚴重，大氣顆粒物污染是影響我國城市空氣質量的首要因素之一[2]。其中的大氣細顆粒物(fine particulate matter 2.5，PM2.5，動力學直徑≤2.5 μm)可進入人體呼吸系統的深部及血液循環系統，長期暴露于PM2.5污染中將嚴重危害人體健康[3-5]。

針對PM2.5的暴露研究，如流行病學調查等多通過實測數據或地面空氣質量監測站完成。然而，結果卻受限于監測站點布設位置或站點分布密度的影響[6]。有別于站點定位監測技術，遙感手段有著大范圍同步監測的優點，應用遙感技術反演城市顆粒物濃度和分布等方面得到了很好的應用[7-9]。大氣中氣溶膠光學厚度(aerosol optical depth，AOD)是表征地表大氣中顆粒物濃度含量的指標之一[10]，未知區域地表PM2.5濃度可通過建立“PM2.5-AOD”二者關系來估算[9]。方法通常有3種：基于傳統的統計學方法(如線性模型)、基于變量的物理關聯以及基于環境因子的耦合建立模型(如高級統計模型)[11]。線性模型最早用于建立二者之間關系[12]，但未考慮變量間的機理聯系，精度較差[6]。物理訂正方法考慮建立大氣濕度，AOD垂直分布、粒徑分布等與PM2.5濃度和分布的物理聯系[13]。應用這類方法的研究有：利用站點測量數據(如激光雷達測量等)訂正AOD垂直分布特征[14-15]、結合大氣邊界層厚度和近地面相對濕度估算PM2.5濃度等[16]。這類方法估算精度高，但部分方法需要先驗的監測數據、且估算方法復雜。

相對于物理訂正法關注于大氣物理狀態的垂直變異，耦合模型方法則側重于研究“PM2.5-AOD”關系隨時空變異情況，即PM2.5的濃度與局地環境氣候條件、地表類型、季節、污染狀況有關[11]。在考慮二者關系空間變異方面，國內外學者建立：包含氣象條件[17]、土地利用[18]等因素的多元回歸模型；結合氣象因素與大氣邊界層高度等，利用神經網絡來構建AOD估算PM2.5的非線性模型[19]；綜合考慮氣象模式、土地參數、地面觀測數據建立的考慮局地信息的地理加權模型[20-21]，這些模型均有效減少了AOD和近地面顆粒物濃度間的不確定性。在PM2.5隨時間變異方面，考慮AOD隨時間/季節變化的問題，恰恰“PM2.5-AOD”關系的許多影響因子都是隨時間變化(如溫度、風速、濕度、PM2.5垂直濃度分布等)，因此通過引入時間因素消除引起兩者關系中的偏頗[22]。這類方法有支持向量機結合模糊粒化時間序列方法[23]、多元自適應樣條回歸[24]、混合效應(/線性)模型[6，25]等。但是，也存在著一定的局限性，需要使用高時間分辨率遙感與監測數據。

而對于在城市尺度上的城市人居環境領域研究中，如城市格局、土地利用布局與PM2.5污染的相互關系，及人口PM2.5污染暴露風險等方面，高時間分辨率的MODIS氣溶膠產品10 km的分辨率過于粗糙；Landsat或類似傳感器的衛星，雖然在空間分辨率能滿足要求，但時間分辨率稍顯不足。因此本文使用混合線性模型，針對Landsat衛星時間分辨率不足而引入季節因子，將季節PM2.5變異作為隨機效應，選擇包頭市建成區為研究區，嘗試以容易獲取的Landsat-8 OLI和MOD09反射率產品為數據源，通過6S模型建立查找表，采用深藍(deep blue)算法反演包頭市氣溶膠光學厚度(AOD)，隨后應用混合效應模型，獲得AOD與實測PM2.5關系模型，同時探討包頭市PM2.5的分布特征，以期為使用類似傳感器的相關工作提供方法借鑒，并為相關部門制定PM2.5防治措施提供一定參考。

1 數據與研究方法

1.1 研究區概況

包頭市為內蒙古自治區最大的工業城市，具有“草原鋼城”之稱，轄昆都侖、青山、東河、九原等四區，建成區面積223 km2(圖1)；屬溫帶大陸性氣候，年降水量300 mm，年均溫6.4 ℃。包頭市擁有眾多高耗能企業，城市以大氣顆粒物為主的環境污染問題突出。本文根據2014—2015年PM2.5監測數據測算，2014年、2015年PM2.5年均值分別為54.9 μg·m-3、50.5 μg·m-3，均超過國家環境污染二級標準，其PM2.5濃度均值季節排序為：冬季>秋季>春季>夏季。

圖1 研究區及環境監測站點分布圖

1.2 氣溶膠反演原理與算法

設大氣水平均一，地表為朗伯面，則大氣上界衛星觀測到的表觀反射率可以表示為：

(1)

式中：ρTOA是大氣頂部反射率；ρ0是大氣的路徑輻射項等效反射率；μs=cos θs，μv=cos θv，θs和θv分別為太陽天頂角與觀測天頂角；T(μs)和T(μv)分別是太陽輻射方向和觀測方向的大氣透過率，可合并表示為T；ρs為地表反射率；S為大氣下界的半球反射率。ρ0、S和T都是氣溶膠模式的函數，代表了大氣的狀態。AOD反演時，首先設定不同的大氣模型、氣溶膠特性、傳感器光譜特性等參數，利用輻射傳輸模型構造查找表[26]；其次依據衛星不同的觀測幾何條件，讀取查找表中不同光學厚度值對應ρ0、S、T的值或其內插值；然后，通過某波段與受氣溶膠影響較小波段的線性關系求解得到ρs；最后，將這些參數代入公式(1)，計算不同光學厚度下假定的表觀反射率ρTOA，與真實的ρTOA比較，差距最小所對應的AOD即為所求[27]。

因城市地區在缺少濃密植被、且冬季植被落葉的問題，而采用Hsu等人提出的適用于亮目標區域的深藍算法[28]，該方法在反演城市AOD方面得到了廣泛的應用[25，27，29]。

1.3 數據來源與處理

選用2014—2015年58期Landsat-8 OLI圖像，軌道號127-32與128-32，云覆蓋小于2%；按照包頭市建成區邊界裁切出研究區范圍。在所獲得的遙感影像中，春季13景、夏季10景、秋季17景、冬季18景(3月到5月為春季，6月到8月為夏季，9月到11月為秋季，12月到次年的2月為冬季)。根據王中挺等的研究結果[30]，將紅光波段反射率大于0.2作為云去除閾值。Landsat-8 OLI的表觀反射率(ρTOA)計算使用ENVI 5.3 軟件的Radiometric Calibration工具。藍波段地表反射率采用MODIS地表反射率產品MOD09A1的藍光反射率，其原理在亮像元地表，藍光波段的地表反射率較小，同期地表反射率不變的特點。重采樣300 m×300 m，根據呂春光等的研究[31]，采用公式(2)校正OLI和MODIS藍波段的光譜差異：

(2)

以0.1為步長至2.0將AOD設定為21個分級(0.000 1，0.1，0.2，……，1.9，2.0)，使用6S模型建立AOD反演查找表，由于缺少全球氣溶膠地基觀測網(AERONET)的氣溶膠觀測數據，將氣溶膠模式固定為大陸型；反射率值大于0.15的像元作為亮目標處理[8，27]，在獲得AOD后，對結果圖像進行10×10窗口平滑處理，重采樣為300 m，消除大氣的不穩定性與填充未參與計算的云或暗目標區的柵格。查找表構建與AOD反演在IDL語言編程下完成[32]。

氣象數據下載自“中國氣象數據網”的中國地面氣候資料日值數據集(V3.0) (http：∥data.cma.cn/data/cdcdetail/dataCode/SURF_CLI_CHN_MUL_DAY_V3.0.html)。2014—2015年包頭市6個站點PM2.5日均值數據(圖1)，來源于中國空氣質量在線監測分析平臺(http：∥www.aqistudy.cn/)。

1.4 分析模型

首先，以6個站點監測的PM2.5值為應變量，以AOD、地表溫度(LST，采用覃志豪的單窗算法[33]反演而來)、日平均氣溫(TEM)、日平均相對濕度(RH)和日平均風速(WIN)作為自變量，采用逐步回歸的方法，發現只有AOD、LST和WIN通過統計檢驗。因此，分別采用線性模型、多元線性模型和混合線性模型(公式(3)、公式(4)、公式(5))，描述“PM2.5-AOD”關系。

PM2.5i，j=a0+a1AODi，j

(3)

PM2.5i，j=a0+a1AODi，j+a2LSTi，j+a3WINi，j

(4)

PM2.5i，j=a0+a1AODi，j+a2LSTi，j+a3WINi，j+
b1，jAODi，j+b2，j+ei，j

(5)

式中：PM2.5i，j、WINi，j、AODi，j與LSTi，j分別為第i個監測站點在季節j(j=1～4，分別表示春、夏、秋、冬4個季節)的PM2.5站點、日均風速的實測值，及AOD、LST的反演值；α0～α3分別為(或混合效應模型中固定效應的)截距和回歸系數；β1，j和β2，j分別為混合效應模型中隨機效應的斜率和截距；對于公式(5)，固定效應是每個地區所有監測站的研究期間在影響因素下的“PM2.5-AOD”間平均關系，隨機效應代表區域二者關系的季節變化。線性模型和多元線性模型可以用決定系(R2)描述方程擬合程度的高低，而混合線性模型的擬合精度則通過比較3個模型的均方根誤差(RMSE)來判斷，即隨機選擇符合要求數據中的90%參與建模，剩余的10%作為驗證數據，去除受到云和暗目標區影響的AOD反演錯誤的點，最后獲得190個樣點。使用Excel軟件進行數據預處理，在SPSS 22.0中進行模型擬合。

2 結果分析

2.1 “PM2.5-AOD”關系擬合

二者關系擬合如表1所示，線性模型回歸的確定系數(R2)為0.28，即AOD的變化可以解釋28%的PM2.5監測值的變異，其均方根誤差(RMSE)較大，為26.53 μg·m-3；引入地表溫度和日平均風速后的多元線性模型R2有較大幅度的提高(0.41)，這些因子的引入對PM2.5監測值變異的解釋能力增強，RMSE也相應的減小到19.70 μg·m-3；混合線性模型的RMSE是在這3個模型中最小，為19.70 μg·m-3。

表1 線性模型、多元線性模型及混合線性模型擬合精度

2.2 包頭市PM2.5變化特征

包頭市PM2.5濃度空間差別明顯，形成以昆都侖區西部，經青山區南部、九原區西北部，至東河區東部為較高污染區，并向西北經九原區、青山區向昆都侖區污染減輕然后又加強的趨勢。從空間分布來看(圖2)，昆都區的西部(包鋼廠區)、中部及其與青山區和九原區交界處建筑密集區、九原區東部(鐵西)、東河區中東部等地區PM2.5濃度較高，昆都侖區(北沙梁)和青山區北部(一機、二機廠區)是PM2.5濃度次之，九原區南部和昆都侖區南部(大片農田)PM2.5濃度較低。

圖2 包頭市PM2.5濃度季節分布特征(單位：μg·m-3)

包頭市PM2.5污染水平較高，且季節變化明顯，即冬季>秋季>夏季>春季。PM2.5的分布特征從時間上來看(圖2和圖3)，包頭市PM2.5濃度冬季濃度水平最高(均值為83.19 μg·m-3)，秋季次之(均值為53.03 μg·m-3)，夏季和春季最低(均值分別為35.98和31.22 μg·m-3)；從各城區PM2.5污染來看，昆都侖區PM2.5的濃度最高(年均值為57.82 μg·m-3)，青山區緊隨其后(年均值為54.56 μg·m-3)，東河區和九原區PM2.5的濃度較低(年均值分別為53.12和52.74 μg·m-3)。

圖3 包頭市PM2.5濃度季節變化統計表

3 討論

3.1 各季節濃度變化分析

PM2.5作為我國大氣的首要污染物，其對人體健康、城市人居環境的影響已成為人們普遍關注的焦點。包頭市PM2.5濃度在時間分布上季均值表現為：冬季>秋季>夏季>春季。對比其他學者的研究，北京市[34-35]、河南省[36]、青島[37]、成都城區[38]等地區均表現類似的分布特征。首先，秋冬季大氣邊界層較低，不利于污染物擴散，加之北方開始采暖，是PM2.5濃度較高；而夏季，太陽輻射強烈，使大氣邊界層升高，加之降雨的清除作用，使得PM2.5濃度處于較低水平[38-39]。在本文中，經計算的PM2.5春季的濃度值低于夏季，這可能與當地的實際情況有些出入，該地區春季干旱、大風、沙塵等天氣也將增加PM2.5濃度。且使用2014—2015年PM2.5實測數據計算的結果是冬季(79.9 μg·m-3)>秋季(52.8 μg·m-3)>春季(43.9 μg·m-3)>夏季(34.7 μg·m-3)，出現這種的原因可能與所獲得的春季數據量少有關，或因遙感數據收集過程中沒有收集到沙塵天氣的影像，以后的研究可以通過獲取多源遙感數據以解決這個問題。

PM2.5濃度空間分布與土地利用因素有關[30]。PM2.5濃度空間分布上的不均勻基本與包頭市工業布局相關，PM2.5濃度高值區為：昆都侖區西部和南部分別有包鋼工業園區和化工區，青山區北部有一機、二機集團和北方重工集團，九原區鋼鐵稀土工業區，東河區河東工業區、磴口工業區和二道沙河小工業區等，以及一些建筑密集區域、棚戶區；低值區則分布與這些區的周邊農田、未建設用地的空地處。

深藍算法為Hsu為了研究沙漠、城市區域、雪覆蓋等地區而建立的[28]，這些亮地表通常在紅光和近紅波段反射率較強而藍光波段反射率較弱[40]，該方法在國內外許多反演AOD的研究中得到應用。但是在包頭市建成區的部分工礦用地和倉儲用地存在著藍色或白色屋頂，及一些水體區域，這些區域強烈的反射藍光，造成深藍算法的誤判，夸大AOD的濃度值，最終導致PM2.5的濃度回歸出現錯誤。

3.2 模型的不足與改進

本研究中混合線性模型的RMSE要低于多元線性模型和線性模型的擬合精度。采用線性模型研究“PM2.5-AOD”的關系，造成其擬合精度較低(通常R2在0.2左右)，近80%無法解釋的變異，這些變異可能是由PM2.5日平均值濃度、氣象因素、AOD的季節變異等因素造成的。首先，PM2.5的監測數據是某天的平均值，遙感數據是某一時刻的快照數據，二者之間存在著時間不匹配的問題。其次，AOD的季節分布可能會出現夏、秋季高于秋、冬季的現象，而在我國國內通常夏季的PM2.5濃度恰恰比較低，所以導致了這個原因。如李成林等研究鄭州市建成區內AOD月度數據也發現，大致在7月和8月AOD值較大，11月、12月、1月較小[41]；本文將張小娟中國2003—2012年AOD分布圖配準，而獲得的包頭地區的AOD值，AOD季節變化：為夏季>春季>冬季>秋季[42]。夏季出現高值的原因可以與高濕、高溫天氣或其促進人為排放污染物發生光化學反應等原因有關，進入秋季氣溫下降，二次氣溶膠轉化減弱導致秋季的AOD濃度減少。混合線性模型導入引起“PM2.5-AOD”的關系的氣象因素、AOD的季節變異等因素，而提高了二者關系的擬合精度。

在本研究中混合線性模型和多元線性模型的RMSE相差不大，并沒有達到Sorek-Hamer結果的精度[6]，采樣點數量較少是造成這種現象的主要原因，本文在抽樣過程中發現，隨著樣本數量的增加，二者的差距在逐漸增大，混合線性模型的擬合精度逐漸變好。

4 結束語

選擇包頭市建成區為研究區，引入季節因子，使用混合線性模型與深藍算法反演得到的包頭市AOD建立回歸關系，并獲取各季節包頭市PM2.5的分布特征。結論如下：與線性模型相比，該模型有較好的回歸精度，適用于監測站點較少與重復周期較長遙感數據擬合的情況；包頭市PM2.5濃度時空差別明顯。PM2.5的分布特征從時間上來看，包頭市PM2.5濃度冬季濃度水平最高，春季最低；從空間分布來看昆都侖區污染較高，青山區北部次之，東河區和九原區最低，形成以昆都侖區西部經青山區南部、九原區西北部，至東河區東部較高污染帶。