京津冀气溶胶光学厚度时空演变及气象解释

文韶鑫， 刘海新， 高叶鹏， 钱以临

(河北工程大学矿业与测绘工程学院， 邯郸 056038)

大气气溶胶是大气中粒径介于10-3～102μm的分子团、固态和液态微粒的总称，具有稳定性高、沉降速度小等特点[1]。气溶胶一方面会影响气候变化，例如影响温度垂直廓线进而影响地球辐射平衡、通过成为云凝结核影响云的光学特性及寿命等；另一方面，气溶胶会降低大气能见度[2-4]。气溶胶中的细颗粒物甚至会危及人类的身体健康，直接影响到人类的生产生活。由于气溶胶处在不断变化之中，理化性质也较为复杂，平均寿命十分短暂，目前对气溶胶的时空分布状况的掌握不够全面，因此亟需增进对气溶胶时空分布及变化趋势的了解。

相关研究表明，气溶胶含量与气溶胶光学厚度(aerosol optical depth,AOD)有极强的相关性，AOD作为大气中的气溶胶含量的表征亦被学术界所认可，并采用遥感手段对AOD展开遥感监测[5-7]。王海林等[8]用AERONET站点做验证，发现C006版本的MOD04数据精度有所提高，且在京津冀地区有良好的适用性。Jia等[9]表明MCD19A2数据反演精度高，可靠性强，在京津冀地区具有良好的适用性。Aklesso等[10]分析了2005—2015年非洲几内亚湾地区AOD变化趋势，发现550 nm时该地区AOD呈下降趋势，而470～660 nm时AOD呈上升趋势。郑小波等[11]根据2000—2010年MODIS气溶胶产品，发现胡焕庸线可作为中国AOD分布的分界线，AOD呈东高西低分布，且东部AOD呈增加趋势。孙忠保等[12]发现1980—2017年中国AOD时空演变趋势大致可分为三个阶段，2008年是AOD由升到降的转折点。通过对COVID-19期间AOD的时空变化进行分析，赵庆志等[13]发现了AOD的周末效应，且不同季节周末效应的差异较大。以上研究使人们对气溶胶时空分布及演变规律有了一定的认识，但针对京津冀地区长时间序列的时空演变特征研究还十分有限。

影响AOD时空变化的因素有很多，大致上可以分为两类，一类是影响气溶胶生成的因素，另一类是影响气溶胶传输扩散的气象因素[12]。目前研究认为人为因素是区域AOD变化的主要驱动力，但气象因子对AOD的驱动作用仍不容忽视[14]。有学者指出，气温、降水、气压、风速等都会对AOD产生一定的影响，且影响具有时空异质性[12,15-16]。陈裕迪等[17]发现气温、湿度以及降水都会对长三角地区AOD产生影响，且气温的影响最为显著，并对此现象做出了详尽解释。景悦等[18]利用地理探测器发现了京津冀地区不同年份AOD产生的主导因子，但未能从长时间序列分析。因此，现利用2001—2019年气溶胶融合数据分析京津冀地区的时空演变特征，并采用地理探测器对AOD的演变进行气象解释，以期增进对京津冀地区气溶胶时空分布及其演变规律的了解。

1 研究区域与数据来源

1.1 研究区域

图1 研究区概况图Fig.1 Overview map of the study area

京津冀地区位于中国华北地区，介于36°05′N～42°40′N，113°27′N～119°50′E，具体如图1所示。京津冀地区是中国重工业基地之一，具有重要的战略地位[19]。由于工业及城市的快速发展，京津冀地区的空气质量在逐渐下降，对人们的生产生活造成了一定的影响[1,3,11]。

1.2 研究数据

1.2.1 MODIS数据

目前使用的气溶胶数据多为MODIS MOD04产品，其深蓝产品(DB数据)覆盖范围广，但分辨率较低，不适合监测局部地区的AOD分布；而暗目标产品(DT数据)覆盖范围有限，难以覆盖全地区。MCD19A2是MODIS官方2017年发布的AOD产品，是NASA官方发布的最新C006版本气溶胶数据，使用MAIAC(多角度大气校正算法)生产的气溶胶日数据产品，具有更高的空间分辨率，能对更为复杂的地表进行AOD的反演[20-21]。由于MCD19A2冬季数据缺失严重，而缺失斑块多呈菱形，导致月合成及年合成数据存在斑块效应，需要MOD04数据作为补充。DB数据具有更高的反演覆盖面积及更广泛的适用性，所以采用 C006版本的MOD04 DB数据和MCD19A2数据作为融合数据。

在数据融合之前先对MODIS数据进行预处理，即影像的提取、投影转换、裁剪与镶嵌，并对MCD19A2数据进行日均值处理。为消除云雾所带来的影响，在融合前还需要将AOD较高(AOD>2)的区域进行掩膜处理。数据融合具体过程如下：对每日MCD19A2数据和MOD04数据进行逐像元扫描比较，若某一像元MCD19A2存在有效值，则采用MCD19A2数据；若不存在有效值，则采用MOD04数据，经过处理得到每日的融合数据。经过月平均处理，得到月AOD融合数据和年AOD融合数据。经对比，融合后的月数据和年数据斑块效应有所减轻。

1.2.2 AERONET数据

AERONET监控网是美国航空航天局(NASA)和法国国家科学研究中心(CNRS)联合创建的全球气溶胶地面监测网络，主要采用CE-318太阳光度计对大气进行实时监测[20]。AERONET的AOD反演精度高，可达0.01～0.02[1]，因此AERONET站点数据是可靠的地面验证数据，常用于验证卫星AOD产品的精度。AERONET全球共有800多个地面站点，其中位于京津冀地区的有4个站点，分别是北京观测站点(39.877°N，116.381°E)、兴隆观测站点(40.396°N，117.578°E)、香河观测站点(39.754°N，116.962°E)和石家庄观测站点(38.000°N，114.550°E)，4个站点的分布如图1所示。本文选用2001—2019年4个站点的月观测数据对卫星气溶胶数据进行精度验证。

由于MODIS数据和AERONET数据的中心波长和时空尺度不同，在进行数据验证前，先对AERONET站点数据进行波段匹配和时空匹配。波段匹配采用Angstrom关系式用AERONET站点波长为440 nm和675 nm数据通过插值得到与MODIS数据波长(550 nm)一致的数据，Angstrom关系式为

(1)

式(1)中：α为Angstrom波长指数；τ为气溶胶光学厚度；λ1为长波波长；λ2为短波波长。

由于云覆盖导致部分数据缺失，需要保证卫星有效数据和地面站点数据的日期一致性，因此本文需要对卫星数据进行采样。同时，由于地面站点监测30 km范围内的AOD，本文选择以对应站点为中心30 km×30 km范围内卫星数据的平均值与站点数据加以比对。

1.2.3 气象数据

气象数据来自中国科学院资源环境科学与数据中心(www.resdc.cn)，为月均值站点数据，经过克里金插值和数据修正得到月气象数据。月气象数据经过均值处理得到季际及年际气象数据。克里金插值是无偏最优估计法，具有较高的精度，适合具有强自相关性的气象数据[21]。受气象站点测量仪器的限制，针对气象数据存在缺失值的问题，研究中取相邻两天数据的均值替代。

2 研究方法

2.1 EOF时空分解

EOF时空分解是对矩阵结构加以分析，提取主要数据特征量的一种方法，该方法在1950年首次被Lorenz引入气象研究中，目前在地学中得到广泛的应用[4]。影像数据集可看作是若干模态及其时间系数的线性组合，EOF可以用少量的模态来呈现AOD的空间分布，并用对应的时间系数图反映AOD的时间分布。若用数据集的距平形式，则模态表示AOD空间变化的分布。EOF分解的数学表达式为

Xmn=VmmZmn

(2)

式(2)中：Xmn为数据矩阵；Vmm为模态矩阵；Zmn为时间系数矩阵；m为采样点个数；n为时间个数。

利用正交函数分析法对2001—2019年均AOD数据集进行时空分解，由于观测点的个数远大于时间个数，需要时空转换来降低计算机的运算量。具体实现方法如下。

(1)对数据集预处理，将数据集处理成距平形式。同时将数据集整合为二维矩阵并去除影像背景值，减小误差的产生。最终得到矩阵Xmn。

(2)计算X的转置矩阵XT与X的交叉积Cnn，即Cnn=XTX。

(3)计算矩阵C的特征向量Vm和特征根λ。

(4)通过时空转换得到X的特征向量V。

(6)对结果进行north检验，得到相互独立的模态。

2.2 线性回归分析法

回归分析是基于像元统计的常用方法，利用回归分析可以探寻AOD时间序列变化趋势的空间分布情况。通过构建年均AOD数据集的拟合函数，进而根据斜率判断AOD变化趋势，斜率k的计算公式为

(3)

式(3)中：t为时间序列长度；i为年份；AODi则表示第i年的AOD值。

根据斜率将AOD分为逐年增大(k>0)和逐年减小(k<0)两类。同时根据相关系数，并依据相关系数查找表，判断变化趋势的显著性：若r>r0.01，则为极显著；若r0.05

2.3 地理探测器

空间分异现象是地学基本现象之一，而地理探测器是探索空间分异性和揭示驱动因子驱动力的统计学工具[22]，通过计算驱动因子与待探测因子的空间相似性来衡量驱动力的大小，结果用q表示，q的计算公式为

(4)

式(4)中：k为驱动因子的类型量；Nk为第k类型驱动因子的样本计数；N为驱动因子总体样本计数；σ2为驱动因子总方差。

地理探测器的独特优势是可以分析两个驱动因子之间的交互驱动作用，其交互作用从弱到强依次为非线性减弱、单因子非线性减弱、双因子增强、独立和非线性增强。本文采用地理探测器探索气象因子对AOD空间分异的驱动作用，由于气象数据是连续变量，需要离散化成为类型量，在综合对比自然断点法、ISODATA以及几何间隔法后，选用ISODATA进行离散化[23]。然后，用10 km×10 km的网格对数据进行采样，使用地理探测器对采样结果进行探测。

3 结果与分析

3.1 融合精度验证

基于AREONET地面站点数据对2001—2019年AOD月数据进行精度验证，结果如图2所示。从图2可以看出，MOD04 DB数据、MCD19A2数据和融合数据与AERONET站点的相关系数r分别为0.71、0.76和0.80，证实：基于本文提出的融合方法，融合后的数据与地面站点的相关性有了稍微提高。两侧的虚线是NASA官方给出的误差范围，为(±0.05±0.15)AOD，EE%为期望误差通过率，结果表明融合数据的期望误差通过率为59%，相比MOD04 DB数据和MCD19A2数据有所增加[24]。多个角度对比说明，融合数据的误差精度要优于MODIS DB数据和MCD19A2数据。

图2 2001—2019年卫星AOD月数据精度验证Fig.2 Accuracy verification of monthly satellitic AOD data from 2001 to 2019

图3展示了研究区AOD年均值统计结果，从图3可以看出，2001—2013年AOD在平均值0.47附近上下波动，其中2003年AOD的异常升高可能与气象有关，年均大气湿度为63%，达到2001—2019年最高值，高出平均值6%，较高的大气湿度使气溶胶粒子聚集变大，AOD上升；同时，平均风速也较其他年份低，加剧了气溶胶粒子的聚集[12]。2008年由于举办北京奥运会，北京及周边地区的环境得到有效改善，AOD显著下降。2014—2019年AOD呈下降趋势，其中2016年AOD显著下降，这可能与2016年河北出台大气污染防治条例有关[5]，政府的大气污染治理取得的显著效果在AOD的迅速下降上得到了体现。

图3 2001—2019年均AOD时间变化Fig.3 Temporal changes of annual average AOD from 2001 to 2009

3.2 AOD年际时空演变特征

AOD的年际空间分布及变化情况如图4所示，根据自然断点法，将AOD≤0.45的地区作为低值区,而AOD>0.56的地区作为高值区。从整体看，AOD低值区主要位于京津冀西北部，包括承德市以及张家口市；而AOD高值区主要位于京津冀东南部，包括了河北省南部的大部分地区。结合图1进一步分析发现，AOD和高程呈明显的负相关性，AOD低值区主要位于山地及高原，包括坝上高原、太行山脉以及燕山山脉，高程较高；AOD高值区主要位于华北平原，高程较低。一方面，气溶胶具有易于扩散的特点，而太行山脉阻挡了气溶胶的扩散；另一方面，人类活动与地势的高低有明显的关系，华北平原地区工业发达，工业区分布较为密集，而工业废气的排放会增加AOD。

图4 2001—2019年AOD空间分布变化Fig.4 Changes in the spatial distribution of AOD from 2001 to 2019

为深入探讨研究区AOD的空间分布模式，本文对2001—2019年京津冀地区年均AOD数据集进行EOF分解，得到了能表征京津冀地区AOD空间分布特征的主要特征向量。第一以及第二模态的贡献度为75%，并通过了north检验，可以代表AOD的空间分布模式。EOF分解后的第一模态和第二模态如图5所示，第一模态的贡献度为68.7%，是主要的空间分布模式。该模态下的所有EOF载荷值符号一致，为全局一致型，表示AOD整体上升或降低。EOF载荷值由东南向西北逐渐降低，表明变化程度从西北地区到东南地区逐渐增大。第二模态的贡献度为7.5%，是对第一模态的补充。第二模态南部地区和北部地区EOF载荷值符号相反，为南北分异型，表示变化方向相反。EOF分解结果表明，京津冀地区AOD的空间分布模式以全局一致型为主，南北分异型为辅。

图5 京津冀地区特征向量Fig.5 Characteristic vector of Beijing-Tianjin-Hebei region

以往研究根据行政区对研究区域进行区域划分，但气溶胶具有空间一致性和连续性，聚类分析更贴合气溶胶特性[4]。基于迭代自组织数据分析算法(ISODATA)对研究区2001—2019年AOD年均结果进行分类，结果如图6(a)所示，从图6(a)中可以看出，研究区被分为明显的3个区域。分别对3个区域进行逐年统计，结果如图6(b)所示，由图6(b)可知，不同区域间的变化情况基本一致。结合图1进一步分析发现，AOD区域划分和海拔高度有密切联系。其中，区域1海拔较高，主要为山地地区，人口密度小，植被覆盖较高，AOD相对较低且变化幅度小。区域2位于山脚前，是区域1和区域3的过渡区，而区域3为平原地区，工业发达，人口密集，AOD较京津冀地区整体平均值高，且变化幅度较大。

图6 ISODATA分区及各分区AOD时间变化Fig.6 ISODATA partition and AOD temporal change of each partition

3.3 AOD季际时空演变特征

基于，依据研究区的气候特点将3—5月划分为春季，6—8月为夏季，9—11月为秋季， 12月—次年2月为冬季[1]。基于上述季节划分及2001—2019年月际AOD时空分布数据进行统计，结果如图7所示，其中，春季、夏季、春季和冬季的AOD分别为(0.46±0.07)、(0.59±0.08)、(0.40±0.04)、(0.37±0.05)。从图7中可以发现，夏季AOD较其他季节高，可能与夏季风带来海洋性气溶胶有关[17,25-26]。春季和夏季AOD标准差较大，说明AOD变化较大，由此可以判定，年际AOD的变化更多来自春季和夏季。

图7 2001—2019年AOD季节变化Fig.7 AOD seasonal changes from 2001 to 2019

为了进一步了解季节变化趋势，采用一元线性回归法，获得研究区AOD不同季节的变化趋势，结果如图8所示。同时，对各季节不同趋势所占比例进行统计，结果如表1所示。由表1可以得出，趋势变化比例最高的是减少不显著，春、夏、秋、冬和全年的比例分别为84.23%、76.83%、59.06%、46.58%和82.39%，而呈显著及极显著增加趋势的比例很低，基本为0。由此可以看出，AOD呈现减少的趋势。春季和夏季京津冀地区主要呈现减少趋势，而秋季和冬季增加不显著趋势的比例明显比春季和夏季高，和减少不显著趋势比例大致相等。结合图8进行分析发现，春季呈增加不显著趋势的地区主要位于冀南地区，而呈显著及极显著减少趋势的地区零星散布在研究区西北部。夏季呈显著减少及极显著趋势的地区明显比其他季节多，零星散布在研究区各地，以北部地区居多；呈增加不显著趋势的地区位于东北角。秋季呈显著及极显著增加趋势的区域在西北部山区，呈增加不显著趋势的区域分布在研究区各地，以北部地区为多。冬季呈现增加不显著趋势的区域主要位于研究区东部，而呈现减少不显著趋势的区域主要位于研究区西部。

图8 AOD四季变化趋势分布Fig.8 Seasonal trend distribution of AOD

表1 AOD不同季节变化趋势所占比例Table 1 The proportion of AOD trends in different seasons

3.4 AOD时空演变的气象解释

为了探究不同气象因子对AOD演变的驱动作用，本文采用地理探测器，在前人的研究基础上选取气温、降水、风速、气压和湿度作为气象因子，并考虑了影响的时空异质性，对不同区域、不同季节的时空演变做出气象解释[12-16]。

京津冀地区年均AOD的地理探测结果如表2所示。表2中对角线表示单个气象因子的解释力(q)，研究中所采用的5个气象因子中湿度、气压和气温的解释力都超过了0.6，解释力最大的是气温，其主要原因为，一方面气溶胶会对地表太阳能量产生散射作用，使气温升高；另一方面气温升高会促进气溶胶的生成，因此气温的解释力最强[12,17]。另外，空气湿度变大，亲水性气溶胶的粒径会增大，AOD同样升高[27]。而大气气压升高，会增加大气厚度导致AOD变大。平均风速的解释力较弱，因为气流对AOD的作用还要受到风向的影响，若是低气压中心，大气向中心流动，气溶胶易集聚，AOD会升高；而在高气压中心，大气的流动会将气溶胶带走，AOD会降低[17]。因而风速和气压呈非线性增强。单一气象因子对AOD的影响有限，更多的是通过多个气象因子的共同作用对AOD产生影响[24]。其中气温与湿度是主导交互因子，解释力高达0.9，说明两者的共同作用对AOD产生更强的影响，温度升高使得气溶胶粒子数增加，而湿度增大会使气溶胶粒径增大，粒子对地表辐射的作用更强。

对图6(a)所划分的区域分别进行地理探测，各区域气象因子的解释力如表3所示。湿度的解释力从区域3到区域1逐渐减小，可能由于来自海洋的海盐气溶胶受到山脉的阻挡，其对AOD的影响逐渐减弱[17]。而气压的解释力从区域3到区域1逐渐增强。风速解释力最强的是区域2，同时风速也是区域2的主导因子，可能是因为风受到山体的阻挡，在此处易形成环流，气溶胶在此处聚集，所以风速越快，AOD越高。

对2001—2019年季均AOD进行地理探测，结果如表4所示。从表4中可以看出，大气湿度不同季节驱动力变化幅度较大，由于气温与湿度的共同作用较强，可能气温对湿度的驱动力会产生较强的影响。降水对气溶胶有沉降作用，在冬季尤为明显，这是因为冬季降水少，每次的沉降作用明显。

表2 年平均AOD气象因子交互解释力及相互作用Table 2 Interaction explanatory power and interaction of annual average AOD meteorological factors

表3 分区气象因子解释力Table 3 The explanatory power of regional meteorological factors

表4 不同季节气象因子解释力Table 4 The explanatory power of different seasonal meteorological factors

4 结论

根据2001—2019年融合数据，采用EOF分析法、一元线性趋势法以及地理探测器等方法，从年、季两个角度分析了京津冀地区的时空演变特征，并探讨了气象要素对AOD演变的影响，对京津冀地区大气污染的防治具有借鉴意义，得出如下结论。

(1)MCD19A2月合成数据和年合成数据会出现斑块效应，可能是由于MCD19A2采用的多角度大气校正算法在遇到云遮盖时会呈菱形形状反演。为了缓解斑块效应，本文根据MCD19A2数据和MOD04 DB数据，采用MCD19A2优先法及平均值合成法生成了AOD融合月数据。同时对MCD19A2数据、DB数据和融合数据进行地面站点验证，发现融合数据的误差精度比MCD19A2数据和MODIS DB数据高。数据融合缓解了斑块效应带来的影响，但完全消除斑块效应有待进一步的研究。

(2)年际时空演变分析表明，京津冀地区AOD在2001—2013年围绕平均值0.47上下波动，2013年以后快速下降，并且空间分布和高程呈明显的负相关性。EOF分解结果显示京津冀地区的年际AOD空间分布模式以全局一致型为主，南北分异型为辅。季节时空演变分析表明，京津冀地区夏季AOD最高，春季和夏季变化程度最大；AOD整体呈现减少趋势，四季大部分地区呈现不显著下降趋势，且秋季和冬季部分地区呈不显著增加趋势。

(3)根据先前的研究，本文选取了湿度、降水、气压，气温、风速这5个气象因子来探究京津冀地区AOD变化的驱动力，在筛选的5个气象因子当中，气温的解释力最强。多个气象因子之间会相互影响，对AOD变化产生更高的影响力。风速和气压呈非线性增强，而气温和湿度是主导交互因子。湿度从区域1到区域3湿度解释力逐渐增强，而气压解释力逐渐减弱，区域2风速解释力比其他区域强。春季和秋季湿度比夏季和冬季的解释力强，而冬季降水的解释力比其他季节强。由于选取的因子有限，可能有更多的气象因子会对AOD的变化产生驱动力，需要进一步的探究。