石家庄市地处太行山中部山坳地带, 地势西高东低, 全年静风频率高, 大气环境容量很小, 对污染物的消化和释放能力很弱, 加之汽车尾气和燃煤污染物排放量大, 多年来大气污染已成为该地区最突出的环境问题。2013年, 石家庄一直处于全国空气质量较差城市的前列, 大气污染防治形势十分严峻。

1.2.1 大气颗粒物地面监测数据

目前, 石家庄全市有人民会堂、世纪公园、西南高教、化工学校、职工医院、西北水源、高新区、封龙山等8个国控空气质量自动监测点。除封龙山外, 其余站点均分布在城区内, 站点分布如图1所示^[11]。监测数据通过河北省空气质量自动发布系统, 实时对外发布。石家庄灰霾天气主要发生在秋冬季节, 本次研究收集了2013年11月10日— 12月31日城区内7个站点的PM₁₀每小时质量浓度数据。参考中国气象科学数据共享服务网提供的地面气候资料日值数据, 该时间段内, 研究区也没有降水等异常天气发生。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 石家庄市区空气质量监测站点分布Fig.1 Distribution of air quality monitoring stations in Shijiazhuang

1.2.2 卫星影像资料

与SPOT, TM等卫星影像相比, MODIS传感器通道多, 动态性好, 数据获取方便, 是目前大气环境应用最为广泛的卫星遥感数据。MODIS传感器分别搭载在EOS-TERRA和EOS-AQUA这两颗太阳同步极轨卫星上, TERRA卫星每日地方时上午过境, AQUA卫星每日地方时下午过境。由于TERRA-MODIS扫描观测资料在第3和第5通道上经常存在条带噪声, 影响地物光谱特征的分析和特征指数的构建精度^[12]。因此本研究采用经辐射定标和定位处理的AQUA-MODIS L1B(MYD02 HKM)标准产品数据, 空间分辨率500 m。

利用MODIS L1B数据快视图文件采用目视解译的方法进行初步判断, 筛选出研究区无云覆盖的影像。在ENVI中使用MODIS Conversion Toolkit插件工具将卫星影像统一配准到UTM投影50带、WGS84坐标系下, 同时校正图像上的“ 蝴蝶结” 效应(bowtie现象), 提取MYD02产品中的反射率数据。进一步使用阈值法剔除云像元, 也可以参考MODIS云产品MYD35数据进行云掩模, 以保证表观反射率信息的准确。根据地面空气质量自动监测站的地理坐标, 在配准后的反射率图像上, 精确定位地面监测站所在像元。为了减小定位误差产生的影响, 以监测站为中心, 计算1.5 km× 1.5 km区域内像元表观反射率平均值, 作为对应空气质量监测站的光谱特征值。

对各监测站点PM₁₀质量浓度小于60 μ g/m³所有日期的相应影像表观反射率计算平均值, 代表受大气颗粒物污染影响较小状态下的监测站下垫面波谱特征(图2)。由于MODIS传感器空间分辨率的

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 MODIS1-7波段各监测站点地物波谱曲线Fig.2 Spectral characters of land cover types in monitoring stations

限制, 遥感图像中存在着大量的混合像元, 这在城市地区尤为明显。考虑石家庄市各监测站点的实际分布位置, 图2所示波谱特征总体表现为建设用地和植被的混合光谱。其中, 西北水源站点位于城乡交错地区(图1), 与分布在建筑密集区的其他站点波谱相比, 西北水源站点波谱中植被“ 红边” 特征较为明显, 即近红外(MODIS2)、可见光红波段(MODIS1)反射率差异更为显著。由于未经大气校正, 受大气影响最严重的蓝光波段(MODIS3)反射率也存在一定差异。

以职工医院和西北水源站点为例, 图3展示了重污染天(PM₁₀质量浓度高于400 μ g/m³)和较低PM₁₀质量浓度天(低于60 μ g/m³)的光谱变化特征。可以看出, 随着PM₁₀质量浓度的增加, 各站点上空反射出去的太阳辐射值在各通道均有不同程度的增加。其中, MODIS1— 4通道的表观反射率增加幅度整体较MODIS5— 7波段大。原因可能是该地区气溶胶粒子类型对表观反射率的影响以散射作用为主, 而当波长大于气溶胶粒子直径后, 使得相应散射效应减弱。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 地物光谱特征随PM10质量浓度的变化Fig.3 Spectral characteristics and their changes with PM10

2.2.1 差值植被指数

已有研究表明, AVHRR, TM和SPOT差值植被指数(DVI)与近地面PM₁₀, PM_5.0质量浓度之间存在一定的相关关系^{[7, 8, 9]}。在陆地遥感观测研究中, 植被指数主要用于反映绿色植物的生长状况, 一般利用近红外通道和可见光红光通道探测数据组合而成。在影像分类研究中, 一般忽略大气影响, 采用未经大气校正的像元灰度值或表观反射率建立植被指数。而在植被生物量、叶面积指数等定量遥感研究中, 往往需要对遥感数据进行大气校正后计算植被指数^[13]。本文研究植被指数定量指示大气颗粒物污染状况的适用性, 因此采用未经大气校正的MODIS影像波段2(近红外波段, 841~876 nm)与波段1(可见光红波段, 620~670 nm)反射率之差计算DVI, 即

DVI=B₂-B₁ , (1)

式中, B₂和B₁分别为MODIS影像波段2和波段1的亮度值。

2.2.2 归一化灰霾指数

Zha等^[10]通过总结南京市霾天和无霾天2景MODIS影像波段1、波段2和波段4(545~565 nm)上典型地物的光谱变化特征, 提出了归一化灰霾指数(NDHI)。在这个研究基础上, 本文尝试提取了石家庄市各监测站点的NDHI, 即

NDHI=(B₁-B₄)/(B₁+B₄) , (2)

式中, B₄为MODIS影像波段4的亮度值。

2.2.3 建筑指数

除排放源和气象因素外, 下垫面性质也直接影响大气颗粒污染物的分布。相关研究发现, 建设用地对近地表颗粒物浓度的影响比绿地明显^[14]。查勇等^[15]基于TM影像典型地物波谱特征, 提出了反映城镇用地信息的归一化建筑指数(normalized difference built-up index, NDBI)。参照TM数据波段设置(1 550~1 750 nm和760~900 nm), 本文选择MODIS影像波段6(1 628~1 652 nm)和MODIS波段2(841~876 nm)构建NDBI, 并参考差值植被指数形式, 提出差值建筑指数(difference built-up index, DBI), 一并用于分析它们与大气颗粒物质量浓度之间的相关关系, 即

NDBI=( B₆-B₂)/(B₆+B₂) , (3)

DBI=B₆-B₂ , (4)

式中, B₆为MODIS影像波段6的亮度值。

根据筛选出的MODIS影像的每日过境时间, 选取地面空气质量自动监测数据中最接近时刻, 即时间相隔在半小时内的PM₁₀颗粒物质量浓度有效数据, 进行统计分析。由于西南高教站点在2013年11月整月及十二月上旬均无有效观测数据发布, 样本数较少, 为了保证各站点间观测时间跨度的统一, 本文最终只讨论人民会堂、世纪公园、化工学校、职工医院、西北水源和高新区等6个城区站点共计190个样本数据的应用。所选站点颗粒物质量浓度的统计特征描述见表1。

表1

Tab.1

表1(Tab.1)

表1 颗粒物质量浓度特征统计量 Tab.1 Statistics for the mass concentration of particulate matter 监测站点 统计量 最小值/ (μ g· m-3) 最大值/ (μ g· m-3) 平均值/ (μ g· m-3) 标准差 样本数 化工学校 31 646 169 139 32 职工医院 48 459 173 117 34 人民会堂 13 621 193 168 33 世纪公园 39 727 233 199 29 高新区 52 418 150 94 31 西北水源 31 463 180 121 31

表1 颗粒物质量浓度特征统计量 Tab.1 Statistics for the mass concentration of particulate matter

由表1可知, 各监测站点均接近或远远超出了环境空气质量标准中的浓度限值, 空气质量较差^[16]。

利用SPSS统计软件分别计算变量大气颗粒物污染指数和变量颗粒物质量浓度(PM₁₀)之间的Pearson相关系数(表2), 从而获得指示大气颗粒物污染程度的最佳污染指数。

表2

Tab.2

表2(Tab.2)

表2 大气颗粒物污染指数与PM10之间的相关关系 Tab.2 Correlation between the spectral indices and PM10 监测站点 统计量 DVI NDHI NDBI DBI 相关系数 sig值 相关系数 sig值 相关系数 sig值 相关系数 sig值 化工学校 -0.296 0.100 0.358 0.045 -0.668 0.000 -0.788 0.000 职工医院 -0.348 0.044 0.642 0.000 -0.664 0.000 -0.704 0.000 人民会堂 -0.511 0.002 0.442 0.010 -0.702 0.000 -0.797 0.000 高新区 -0.248 0.178 0.001 0.995 -0.630 0.000 -0.657 0.000 世纪公园 -0.501 0.006 0.497 0.006 -0.724 0.000 -0.887 0.000 西北水源 -0.515 0.003 0.250 0.175 -0.104 0.576 -0.128 0.494 全部站点 -0.315 0.000 0.308 0.000 -0.599 0.000 -0.704 0.000

表2 大气颗粒物污染指数与PM10之间的相关关系 Tab.2 Correlation between the spectral indices and PM10

从表2中可以看出, 除NDHI外, DVI, NDBI和DBI与大气颗粒物质量浓度均呈负相关关系。对于DVI, 西北水源站点相关性最高, 且通过了0.01的显著性水平; 而其他站点相关性较低, 特别是化工学校、职工医院和高新区站点, 与PM₁₀均未通过0.01的显著性水平检验。相反, 对于NDBI和DBI, 西北水源站点相关性却是最低的, 甚至未通过0.05的显著性检验。结合图1和图2, 西北水源站点分布位置和下垫面波谱特征的独特性, 是造成以上相关性结果的主要原因。就本研究区来说, 除职工医院一个站点外, NDHI与大气颗粒物质量浓度之间的相关性并不明显。Zha等^[10]的研究指出, 相对于其他地物, 水体像元的NDHI随PM₁₀质量浓度变化更为敏感, 而北方城市市区内很少有大面积水域存在, 一定程度上限制了NDHI的应用。Zha等的研究利用典型地物NDHI和整个城市平均PM₁₀进行相关分析, 这也可能是造成与本文相关性结果不一致

的原因。整体来看, 除西北水源站点外, NDBI和DBI在其他站点均与PM₁₀保持了较为稳定的显著相关性, 而且DBI相关性更好。

通过上文相关性分析, 发现土地覆被类型对大气颗粒物污染指数应用效果有一定的影响。而城市地区土地覆被类型复杂多样, 加之中分辨率卫星遥感图像中大量混合像元的存在, 因此, 本文考虑综合多种大气颗粒物污染指数, 采用逐步回归分析方法构建大气颗粒物质量浓度的估算模型。从190个样本数据中随机选取160个数据建立估算模型, 剩余30个数据作为独立样本对估算模型进行精度验证。

PM₁₀逐步回归分析模型中自变量加入顺序为DBI, NDHI, NDBI和DVI, 模型决定系数R²分别为0.620, 0.692, 0.733和0.747, 以上模型Sig值均< 0.001, 均通过了显著性检验。最终确定PM₁₀的最优遥感指数估算模型为

PM₁₀=218.178-17 432.9DBI+1 858.361NDHI+1 730.713NDBI-2 025.427DVI 。 (5)

以上模型能够解释PM₁₀74.7%的变异, 160个建模样本数据的拟合平均相对误差为18.2%。进一步利用剩余30个检验样本数据对模型进行估算评价, 平均相对误差为8.8%。图4显示了建模样本和检验样本大气颗粒物质量浓度实测值和估算值之间的散点图。整体来看, 大气颗粒物污染指数具有指示PM₁₀质量浓度变化的能力。基于表观反射率的大气颗粒物污染指数反映的是地表到大气层顶整层大气的变化状况, 与近地面颗粒物浓度的相关关系受风速、混合层高度、相对湿度等气象因素以及颗粒物形状、组成等的综合影响, 有可能会产生个别较大的估算误差。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 大气颗粒物质量浓度实测值与估算值Fig.4 Plots of measured against the estimated PM10

以2013年11月14日为例, 基于以上建立的估算模型, 从MODIS遥感影像中计算获得了PM₁₀质量浓度的空间分布, 如图5所示。化工学校、职工医院、人民会堂、高新区、世纪公园和西北水源等站点的估算相对误差分别为5.1%, 14.5%, 8.4%, 27.7%, 15.3%和35.8%。从图中可以看出, 城区西南和中部地区估算获得的PM₁₀变化趋势较为连续, 而北部和东部地区空间分异较为明显。在逐步回归分析中, 决定系数越大, 自变量对因变量的解释程度越高, 自变量引起的变异占总变异的百分比高, 相关的密切程度越大。所确定的最优遥感指数估算模型中仅单因子差值建筑指数(DBI)就能够解释PM₁₀ 62.0%的变异。对照图1中城区土地利用类型分布来看, 建设用地也主要集中在城区西南和中部地区。综合以上相关性分析和估算模型误差检验结果, 本文所建估算模型更适用于中分辨影像上城市建筑密集区的应用。具体应用中, 可以参考城市土地利用图或者采用阈值法对非建筑用地的部分(如水域、耕地等)进行掩模处理, 以获取较为准确的大气颗粒物浓度分布信息。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 PM10质量浓度估算值空间分布Fig.5 Distribution of the estimated PM10