1 引言(Introduction)

气溶胶是由固相和液相小质粒共同组成的多相态混合体系, 存在于不断变化的大气中, 由于地面下垫面的组成复杂, 气溶胶的分布呈现出显著的多变性和非均匀性(宗雪梅等, 2005).气溶胶在季节更替中颗粒大小、颗粒类型也发生改变, 是全球气候变化数值模拟中不确定的因子.在2007年第四次IPPC评估报告中明确指出气溶胶辐射强迫的不确定性原因是缺乏对气溶胶光学特性和时空分布特征认识(Solomon et al., 2007).要全面掌握和模拟全球气候变化, 必须研究气溶胶的光学特性和时空变化特征.环渤海地区位于我国东北部沿海地区, 经济发展迅速、自然环境优越、人口相对密集, 是全国经济发展的核心区域之一, 污染也较为严重, 然而对该区域气溶胶时空分布研究相对较少, 对环渤海气溶胶的光学特征的研究具有重要价值.

在气溶胶光学特性的研究中, 地基遥感观测研究被认为是较为精确的, 常作为卫星遥感的验证数据(Holben et al., 2001 ).卫星遥感观测空间范围覆盖广, 时间序列连续性好, 大气气溶胶光学厚度(Aerosol Optical Depth, AOD)是研究气溶胶的基本光学参数, 能够反映整个大气层气溶胶的消光作用并评估空气污染程度(刘海知等, 2019).针对卫星遥感获取气溶胶参数以及研究气溶胶分布特征这一重要科学前沿问题, 国内外学者取得了广泛的研究成果.张军华、毛节泰等利用日本静止气象卫星(GMS)遥感中国25个湖面上空的气溶胶光学厚度数据, 对比地面观测的数据, 取得有价值的结论(张军华等, 2003).齐玉磊等对比了2006—2009年中分辨率成像光谱仪(MODIS)和多角度成像光谱仪(MISR)传感器与我国北方4个地面AERONET站点观测的气溶胶光学厚度, 发现北方主要有两个高值区, 分别是塔克拉玛干沙漠和华北南部以及华东北部地区(齐玉磊等, 2013).刘莹等利用现代研究与应用再分析数据集(MERRA2), 通过趋势分析法分析了1990—2017年中国地区的气溶胶光学厚度数据, 发现气溶胶光学厚度年均值整体呈波动上升趋势, 季节尺度上整体呈春夏秋依次递减, 冬季有所提升的趋势, 地形和气候环境是其主要影响因素(刘莹等, 2019).

本文利用最新版本MODIS/Aqua C6.1气溶胶产品数据, 验证并优化出适合环渤海区域的新融合数据, 分析环渤海地区近10年气溶胶光学厚度分布状况, 以及季节与月变化特征, 以期为气候变化数值准确模拟, 应对全球变化提出合理依据.

2 材料与方法(Materials and methods) 2.1 研究数据与预处理

研究采用的卫星数据是NASA发布的MODIS/Aqua C6.1气溶胶日产品数据集(https://ladsweb.modaps.eosdis.nasa.gov/), 所有质量等级为3级的数据, 波段为550 nm, 时间为2009年1月—2018年12月, 空间分辨率为10 km×10 km.图像是栅格格式, 利用MATLAB软件将数据按照经纬度重采样在新建立的标准经纬度网格中, 网格为34°~44°N, 113°~126°E, 网格大小为0.1°×0.1°, MODIS数据空间分辨率较低, 为避免出现条带, 将月均和季节平均数据重采样在0.25°×0.25°网格.利用程序把原始数据×10^-3转换成气溶胶光学厚度的实际值, 对像元值小于0的无意义值进行剔除.日产品数据主要做验证, 把日产品数据处理之后生成月数据和季节数据以及年数据, 对该区域进行分析.用于验证的地基数据来源于全球自动观测网(https://aeronet.gsfc.nasa.gov/), 选择有质量保证的2级数据产品, 选取有代表性的3个站点：北京、香河和兴隆, 用于验证的实测数据站点位于天津科技大学滨海校区图书馆楼顶, 用到仪器为Microtops Ⅱ手持太阳光度计, 测量时间为2018年3月—2018年11月, 所有站点详细信息见表 1.在波段对应上, 利用地基和实测440 nm和675 nm波段AOD数据插值出550 nm波段上的AOD数据, 该计算过程直接在Excel数据表格中完成；在空间匹配上, 选择监测站点周围30 km缓冲区内的MODIS AOD平均值进行验证, 即每个验证站点选择标准网格0.3°×0.3°的数据；时间匹配上, 根据卫星过境时间, 地基及实测数据选择13:00—14:00的数据, 作为真实数据.

用来生成新融合数据的归一化植被指数(NDVI)数据来自NASA官网, MODIS/Aqua的MYD13C2, 空间分辨率为0.05°×0.05°, 栅格格式存储, 按照经纬度重采样在新建立的标准经纬度网格中.

2.2 研究区域

渤海为中国内海, 总面积达7.7万km², 环渤海地区包括：与渤海直接相接的河北省、辽宁省、山东省、天津市以及未直接靠海的北京市等邻近省份区域(图 1中灰色区域为研究区域, 白色为海洋), 研究区范围为34°~44°N, 113°~125°E, 典型的季风气候, 夏季炎热多雨, 冬季寒冷干燥.

2.3 研究方法

Sen′s slope趋势分析法是通过计算数据序列处理后的中值来反映数据变化的趋势度(王佃来等, 2013), 具体公式如式(1)所示.

(1)

式中, x_i和x_j表示数据在i时刻和j时刻的值, 本文是第i年和第j年的AOD值, 当β>0时, AOD在时间序列上呈上升的趋势, 反之呈下降的趋势.

Sen′s slope分析可以很好的反映变化趋势, 但对变化趋势的显著性缺乏判断能力.MannOR＝Kendall为一种非参数检验法, 用来检验趋势变化的显著程度, 两种方法结合增强抗噪性, 提高准确性.MannOR＝Kendall方法做出两种假设, 原假设：时间序列上由n个样本随机独立排列组成, 备用假设：在时间序列上样本呈递增或递减排列组成, 检验统计量由S表示：

(2)

(3)

根据n的大小差异采用不同的显著性检验方式, 当n < 10时, 采用双边趋势检验, 给定显著性水平α下, 拒绝原假设接受备用假设, 反之则反.如果S>0, 认为序列存在上升趋势；S=0无明显趋势；S < 0表明序列存在下降趋势.

当n≥10时, 使用检验统计量Z值进行显著性检验, 用如下公式表示：

(4)

其中, , n是时间序列中数据个数, m是重复出现数据组的个数, t_i是第i组重复出现数据组中的重复数据个数, 在已知置信水平的条件下, 从正态分布表确认临界值Z_1－α/2, |Z|≤Z_1－α/2则接受原假设, 说明趋势变化不显著, |Z|>Z_1－α/2则拒绝原假设, 说明趋势变化显著.

本研究中时间从2009到2018, 即n=10, 采用Z值检验, 取α=0.05, 查表Z_1－α/2=1.96.分为4类趋势：数据序列呈显著上升(β>0, |Z|>1.96)；数据序列呈显著下降(β < 0, |Z|>1.96)；数据序列呈上升趋势但不显著(β>0, |Z|≤1.96)；数据序列呈下降趋势但不显著(β < 0, |Z|≤1.96).

3 结果与分析(Results and analysis) 3.1 数据验证与新融合数据生成

暗像元算法(Dark Target, DT)又称浓密植被算法, 针对低反射率的表面, 反演效果最佳, 对裸露的亮区域上空AOD反演精度偏低.深蓝算法(Depth Blue, DB)弥补了暗像元算法的适用区域的局限性, 但对于地表反射率随时间变化较快的过渡区反演精度显著下降.卫星数据的验证结果如表 2所示, 北京站作为城市区的代表, 相较于DT算法产品, DB算法产品在相关系数R(变量间相互关系的密切程度, 越接近1关系越密切)、均方根误差(Root Mean Squaye Error, RMSE, 气溶胶属性改变造成随机误差的大小, 接近0最优)以及相对平均偏差(Relative Mean Bias, RMB, 大于小于1表征验证值的高估和低估)都优于DT算法产品.香河站作为郊区的代表, DT算法的相关系数高于DB算法, 但高估数据出现更多.兴隆作为林区的代表, DT算法RMB优于DB算法, 在相关系数和RMSE方面DB算法更优.滨海新区站点作为近海区代表, DB算法在相关系数优于DT算法, RMSE和RMB较DT算法差.

环渤海地区林地和裸地的过度地区占据大部分面积, 地表反射率随时间和空间变化明显, 单一算法无法达到反演精度.融合算法提供了很好的思路：当NDVI>0.3时选择DT AOD; 当NDVI < 0.2时选择DB AOD; 当0.2≤NDVI≤0.3时, 如果DT和DB的AOD产品质量等级QA_DT=3或QA_DB≥2, Combined AOD选二者中最优, 若两者QA均为3, 则取二者平均值, 如果QA_DT < 3且QA_DB < 2, Combined AOD为无效值(Levy et al., 2013).

在环渤海地区选择地基数据较多的三站点对融合算法进行评价, 在NDVI < 0.2覆盖区域, 北京和香河站点DT算法有效匹配率(光度计与卫星同天同站点有数据)都很低, 只有11%且高估严重, 认为融合算法对NDVI < 0.2时数据选取DB数据是合理的；在当0.2≤NDVI≤0.3时, 从RMB来看, 三站点分别为1.14、1.13、1.34, DT算法显著高估, 且有效匹配率很低, 分别只有35%、24%、11%, 融合的数据有效匹配率, 相关系数并没有提高, 甚至有所降低；NDVI大于0.3时, 融合算法只选取了DT算法数据集, DT算法拟合的相关系数更高, 三站点DB算法数据有效匹配率很高, 直接舍弃DB算法的融合在环渤海地区并不是很适用.

为找出DT算法在环渤海地区的适用范围, 对NDVI>0.3的数据进行细分为0.3 < NDVI≤0.4、0.4 < NDVI≤0.5、0.5 < NDVI≤0.6、NDVI>0.6进行验证.验证结果表明, 在北京和香河从低的NDVI暗像元算法已经在有效匹配率很高, 但在兴隆站点在NDVI>0.5时才有好的有效匹配率和RMB.对此提出优化后的融合算法：对环渤海区域10年的NDVI做平均, 分为NDVI高值区(大于0.5)格点和NDVI低值区(小于等于0.5)格点, 对于高值区, 月均NDVI＞0.5时取DT和DB的均值, 月均NDVI≤0.5时, 只取DB的数据, 没有值为NAN；对于低值区格点, 月均NDVI＞0.3时取DT和DB的均值, 月均NDVI≤0.3取DB的数据, 没有值为NAN.

新融合数据验证如表 5所示, 相对于原始3种算法, 在保证相关性高的情况下, 提高了数据的匹配率, RMB更接近1, RMSE更小.

3.2 环渤海AOD时空分布特征 3.2.1 月均AOD分布变化特征

将10年AOD月均空间分布图如图 2所示, 从11月开始, 内蒙区域AOD值开始增加, 到2月达到最高, 内蒙区域AOD值冬春季的上升, 和偶尔几次自然产生的沙尘气溶胶有直接关系, 同时东南区域也在这个时段增至一个峰值, 达0.8左右, 环渤海地区形成两个高值区, 东南平原地区高值区的出现主要是受到人类活动持续产生污染颗粒物的影响, 而11月—1月, 东北辽宁、吉林区域AOD呈不断下降趋势, 1月降到0.1左右, 为该区域全年最低, 这个时期大尺度、高强度的天气系统作用, 稀释了这个区域的气溶胶.2—5月, 区域整体AOD呈下降趋势, 尤其在环渤海西北的内蒙及河北变化幅度最大, 低值区区域不断扩大延伸, 从0.7一直降到0.1~0.2左右.6月又一次出现高值月, 本月高值的范围也是各个月份中最大, 包括山东半岛在内的全部东南区域, 以及东北辽宁、吉林也达到全月份的最高值, 平均有0.5~0.6左右, 这个时段升温明显, 气温高强辐射有助于光化学反应进行生成二次气溶胶, 降雨增加, 气溶胶的吸湿增长明显.6—9月, 区域AOD第二次整体下降, 从东南大面积高值到只有河北南方的小区域高值, 到9月降到最低, 整个区域平均只有0.25左右；9—11月区域AOD变化不明显, 整体有微弱的上升.

2009—2018年研究区域全部陆地地区的月AOD时间变化如图 3所示, 1—12月总体呈先升高后降低的趋势.2014年、2015年2月月均AOD达0.7以上, 显著高于其他年份, 频繁出现重度雾霾天气, 是研究区东南区域空间分布出现高值的原因之一.月平均AOD时间变化显示6、7月在各年都处于高值水平, 近年来有下降的区域, 2018年6月只有0.36.9月份在各个年份都处于低值水平, 最高只有0.3.

3.2.2 季节AOD分布变化分析

季节分析中, 3—5月为春季, 6—8月为夏季, 9—11月为秋季, 12月—次年2月为冬季.结果如图 4所示, 4个季节空间分布总体表现为东南高西北低, 区域整体夏季平均为0.433, 为四个季节最高, 高值以河北东南为中心, 向北逐步降低, 到内蒙西部最低为0.2左右, 夏季高温高湿高辐射, “气OR＝粒”转化过程大大缩短, 水溶性气溶胶吸湿膨胀, AOD值快速增高(Kang et al., 2016), 夏季环渤海近海周围AOD较高, 原因可能受季风影响, 沿海区湿度更高, 气溶胶吸湿增长, 相对于内陆的干气溶胶质量浓度, AOD值较大；春季平均为0.404, 仅次于夏季, 冬季局地沙尘的滞留, 以及季节空气干燥, 北方风大将地表浮尘卷起；冬季平均为0.388, 东南高值区连成一片, 内蒙高原出现中心为0.6的高值区域, 向外进行扩散, 这些高值的出现除受到内蒙古高原的风沙、扬尘重要影响外, 稳定的大气层不利于污染颗粒物的扩散, 很容易生成雾霾也是重要原因；秋季平均显著低于其他3个季节, 与夏季相比, 西北部低值区地区明显扩大, 向北延伸至辽宁西部, 向南到山西区域, 这部分区域在0.1左右, 整个区域平均只有0.275, 空气层结渐趋稳定, 混合层高度降低, 温度降低不易于空气污染物的扩散(谢志英等, 2015).

3.2.3 年际AOD分布变化分析

结合图 5与表 6可知, β值在整个研究区域99.8%表现为负值, 环渤海99.8%区域近10年AOD呈降低趋势, 北京奥运会前后对环境的治理以及“生态文明”建设方针的提出起到重要作用(Wang et al., 2012；刘浩等, 2015), 23.24%的区域AOD显著下降, 环渤海东北角及西北的张家口、承德区域下降显著, 这可能与张家口承办2022年冬季奥运会部分项目有很大联系, 研究区域东北小面积和西部也有一部分下降显著, 环境保护和改善效果显著.微弱下降的区域占据大部分面积, 达76.56%, 只有辽宁中部地区和内蒙古小面积区域AOD呈上升趋势, 但上升幅度不高.另有研究显示, 2008—2015年中国北部每年平均下降0.005~0.01, 而在工业地区可以达到每年0.02的下降速度(Sun et al., 2019), 近年来空气质量有所提高, 空气污染有所改善, 微弱上升区域的存在说明针对易出现沙尘暴的区域, 预防和治理措施需不断加强, 针对林区, 偷伐树木、森林火灾的情况还在发生, 生态文明建设需持续推进.

3.2.4 环渤海地区AOD空间分布

10年平均AOD的空间分布显示见图 6, 整体上AOD呈东南高西北和东北低的趋势.高值区主要分布在以河北省中部和东南部平原地区为中心, 中心区域AOD值达0.7左右, 向北延伸至北京东南、天津的平原区域, 向南延伸山东西部、北部黄河冲积形成的平原, 河南、安徽、江苏北部, 这些区域也在0.6~0.7之间.这些区域平均海拔低, 空气流动性较差, 气溶胶难以向外扩散, 工业发达, 人口密集, 人类活动频繁, 气溶胶排放量大, 造成了高的AOD值(He et al., 2016).地形的客观存在限制了气溶胶在水平方向的扩散(李成才等, 2003), 东南的气溶胶很难通过水平运输到达北部的高原山地, 低值区以河北西北和山西东北部、内蒙西部交界, 以及内蒙古高原东部和辽宁东部最低, 只有0.2左右, 其余北部区域AOD在0.3~0.4左右.这些区域人口不集中, 污染物排放就相对较少, 平均海拔高, 风速高, 有助于气溶胶的快速释放, 地形以山地和丘陵为主, 森林草地覆盖率高, 有利于空气的净化.韦晶等(2015)指出AOD＜0.3时, 空气质量良好无污染, 而在0.3~0.8之间, 空气质量一般, 当AOD＞0.8时, 空气质量明显较差, 有一定的污染, 甚至可能有雾霾.环渤海AOD低值区最低区域空气质量良好, 整个区域空气质量一般的区域占据大部分面积, 大面积空气质量的实时监控和环境改善措施的精准实施尤为重要, 遥感观测起到关键作用.

4 结论(Conclusions)

1) 经实测和地基AOD数据检验, MODIS C6.1产品数据在环渤海地区并不完全适用, 优化出新融合数据：对环渤海区域10年NDVI做平均, 分成NDVI高值区(＞0.5)格点和NDVI低值区(≤0.5)格点, 对于高值区, 月均NDVI＞0.5时取DT和DB的均值, 月均NDVI≤0.5时, 只取DB的数据；对于低值区格点, 月均NDVI＞0.3时, 取DT和DB的均值, 月均NDVI≤0.3时, 只取DB的数据.

2) 月均AOD的研究中：2月和6月成为全年变化的两个高值月, 2014年和2015年2月AOD达0.7以上, 每年6、7月均表现为高值, 9月为低值；季节分析中, 区域多年平均夏季表现为最高值(0.433), 秋季表现为最低值(0.275).Sen+MannOR＝Kendall趋势分析法分析年际AOD变化特征, 环渤海99.8%区域AOD呈降低趋势, 微弱下降的区域占据大部分面积, 但在内蒙古高原和辽宁中部仍有微弱上升区域, 生态文明建设和环境保护政策取得巨大的成果, 但防风固沙、退耕还林相关措施仍然需要进一步推进.

3) 环渤海地区气溶胶光学厚度整体呈东南高西北低的趋势.AOD高值(0.6~0.8)集中在除山东半岛和山东中部山地丘陵外的东南区域, 低值(0.2~0.3)分布在河北西北部高原地区、内蒙古高原, 北京的西北山地, 以及辽宁西部和东部的山地和丘陵区域等区域.