包頭市細(xì)顆粒物遙感監(jiān)測(cè)混合線性模型

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

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

0 引言

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

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

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

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

1 數(shù)據(jù)與研究方法

1.1 研究區(qū)概況

包頭市為內(nèi)蒙古自治區(qū)最大的工業(yè)城市，具有“草原鋼城”之稱，轄昆都侖、青山、東河、九原等四區(qū)，建成區(qū)面積223 km2(圖1)；屬溫帶大陸性氣候，年降水量300 mm，年均溫6.4 ℃。包頭市擁有眾多高耗能企業(yè)，城市以大氣顆粒物為主的環(huán)境污染問(wèn)題突出。本文根據(jù)2014—2015年P(guān)M2.5監(jiān)測(cè)數(shù)據(jù)測(cè)算，2014年、2015年P(guān)M2.5年均值分別為54.9 μg·m-3、50.5 μg·m-3，均超過(guò)國(guó)家環(huán)境污染二級(jí)標(biāo)準(zhǔn)，其PM2.5濃度均值季節(jié)排序?yàn)椋憾?秋季>春季>夏季。

圖1 研究區(qū)及環(huán)境監(jiān)測(cè)站點(diǎn)分布圖

1.2 氣溶膠反演原理與算法

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

(1)

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

因城市地區(qū)在缺少濃密植被、且冬季植被落葉的問(wèn)題，而采用Hsu等人提出的適用于亮目標(biāo)區(qū)域的深藍(lán)算法[28]，該方法在反演城市AOD方面得到了廣泛的應(yīng)用[25，27，29]。

1.3 數(shù)據(jù)來(lái)源與處理

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

(2)

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

氣象數(shù)據(jù)下載自“中國(guó)氣象數(shù)據(jù)網(wǎng)”的中國(guó)地面氣候資料日值數(shù)據(jù)集(V3.0) (http：∥data.cma.cn/data/cdcdetail/dataCode/SURF_CLI_CHN_MUL_DAY_V3.0.html)。2014—2015年包頭市6個(gè)站點(diǎn)PM2.5日均值數(shù)據(jù)(圖1)，來(lái)源于中國(guó)空氣質(zhì)量在線監(jiān)測(cè)分析平臺(tái)(http：∥www.aqistudy.cn/)。

1.4 分析模型

首先，以6個(gè)站點(diǎn)監(jiān)測(cè)的PM2.5值為應(yīng)變量，以AOD、地表溫度(LST，采用覃志豪的單窗算法[33]反演而來(lái))、日平均氣溫(TEM)、日平均相對(duì)濕度(RH)和日平均風(fēng)速(WIN)作為自變量，采用逐步回歸的方法，發(fā)現(xiàn)只有AOD、LST和WIN通過(guò)統(tǒng)計(jì)檢驗(yàn)。因此，分別采用線性模型、多元線性模型和混合線性模型(公式(3)、公式(4)、公式(5))，描述“PM2.5-AOD”關(guān)系。

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

2 結(jié)果分析

2.1 “PM2.5-AOD”關(guān)系擬合

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

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

2.2 包頭市PM2.5變化特征

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

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

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

圖3 包頭市PM2.5濃度季節(jié)變化統(tǒng)計(jì)表

3 討論

3.1 各季節(jié)濃度變化分析

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

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

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

3.2 模型的不足與改進(jìn)

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

在本研究中混合線性模型和多元線性模型的RMSE相差不大，并沒(méi)有達(dá)到Sorek-Hamer結(jié)果的精度[6]，采樣點(diǎn)數(shù)量較少是造成這種現(xiàn)象的主要原因，本文在抽樣過(guò)程中發(fā)現(xiàn)，隨著樣本數(shù)量的增加，二者的差距在逐漸增大，混合線性模型的擬合精度逐漸變好。

4 結(jié)束語(yǔ)

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