包头市细颗粒物遥感监测混合线性模型

同丽嘎，李雪铭，张靖，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浓度冬季浓度水平最高，春季最低；从空间分布来看昆都仑区污染较高，青山区北部次之，东河区和九原区最低，形成以昆都仑区西部经青山区南部、九原区西北部，至东河区东部较高污染带。