贝叶斯层次时空模型在甘肃省AOD 时空分布中的应用

吕建康, 孙建国, 谢 甫

(1.兰州交通大学 测绘与地理信息学院,甘肃 兰州 730070;2.地理国情监测技术应用国家地方联合研究中心,甘肃 兰州730070;3.甘肃省地理国情监测工程实验室,甘肃 兰州730070)

0 引言

作为地气系统的重要组成成分之一,大气气溶胶具有结构复杂、生命周期短和时空变化剧烈等特点[1]。 气溶胶通过吸收散射太阳辐射、改变云物理特性和大气化学过程,对地-气辐射能量收支平衡、全球气候变化以及空气质量产生重要影响[2]。 同时,大气气溶胶对人体健康有很大影响,其中的可吸入颗粒物会通过肺部血管进入血液,增加致癌风险[3]。 气 溶 胶 光 学 厚 度(Aerosol Optical Depth,AOD)是遥感监测气溶胶的一个重要光学特性参数,该参数描述了气溶胶对光的衰减作用,在一定程度上可以指示空气污染程度[4]。 因此,全面、准确地探测和分析气溶胶时空分布特征是应对气候变化和防治大气污染的重要基础工作,具有重要的理论价值和实践意义。

目前,对大气气溶胶的监测手段主要有地基监测和卫星遥感监测。 其中,地基监测精度较高,但地面站点较少、位置分布不均且多为单点观测,无法进行连续的、大范围的污染物浓度变化监测;卫星遥感监测具有空间上和时间上连续性强、覆盖范围广等优势,在气溶胶监测和分析工作中发挥重要作用[5-6]。 为满足不同区域气溶胶研究的需要,不同传感器平台均发布相应的气溶胶产品,如中分辨率成像光谱仪(Moderate Resolution Imaging Spectroradiometer,MODIS)、多角度成像光谱辐射计(Multiangle Imaging Spectro-radiometer,MISR)、超高分辨率扫描辐射计(Advanced Very High Resolution Radiometer,AVHRR) 和臭氧监测仪(Ozone Monitoring Instrument,OMI)等。 MODIS 气溶胶产品已经积累超过15 年的全球气溶胶观测数据,为不同区域、长时间序列的气溶胶研究提供了良好的数据基础[7],一直受学者们的青睐。 赵仕伟等[8]将MODIS C6 AOD 融合数据与AERONET 站点数据进行验证,确定了MODIS C6 在中国西北地区的适用性,分析得到了西北不同气候区域AOD 的时空变化特点及形成原因。 张瑞芳等[9]研究发现,西北地区AOD 月均值变化起伏较大,主要因为沙尘天气频繁发生,AOD高值区处于塔克拉玛干沙漠和陕西关中地区。 胡俊等[10]和李霞等[11]分别利用MYD04_L2_C006 数据与CE318 自动跟踪太阳光度计数据对乌鲁木齐的四季AOD 均值进行分析,得出AOD 值在春季最高,夏冬两季次之,秋季最低,春季AOD 值高原因是风沙天气的影响。 林泓锦等[12]分析得出内蒙古东北部和中西部AOD 水平较高,主要受人为源气溶胶粒子影响,内蒙古中部主要以沙尘气溶胶为主,全年AOD 峰值在4 月和6 月出现。 目前针对AOD 时空特征的研究主要采用Mann-Kendall 趋势检验、Sen's slope 变化趋势分析等方法,只能定性或定量表达空间分布在时间序列上的整体特征,难以揭示时间动态和空间分布的耦合过程。 贝叶斯时空模型(Bayesian Spatio-Temporal Model,BSTM)将时空过程分解为空间格局、时间趋势和局部变化3 个部分,精确地分析AOD 的时空特征,很好地弥补了传统方法的不足[13]。 为了能够更深层次地挖掘研究对象的时空规律,BSTM 已经出现了许多拓展模型[14],贝叶斯层次时空模型(Bayesian Hierarchical Spatio-Temporal Model,BHM)是拓展模型中的一种,是贝叶斯层次模型和时空交互模型的结合,该模型通过充分利用先验信息,考虑时空交互项,可以更加准确地分析AOD 时空分布。 BHM 在流行病学、统计学[15-17]等方面已经得到广泛应用,但在AOD 时空特征分析中的应用还不多见。

本文以甘肃省为研究区,利用BHM 从时间、空间以及局部变化3 个方面对甘肃省AOD 展开研究,以期为应对环境污染和气候变化的决策提出合理依据。

1 研究区概况与数据源

1.1 研究区概况

甘肃省位于北纬32°11′~42°57′,东经92°13′~108°46′,中国西北部,地处黄土高原、青藏高原和内蒙古高原三大高原和西北干旱区、青藏高寒区、东部季风区三大自然区域的交汇地带,地形复杂、山脉纵横、海拔相差悬殊,地势由东南向西北倾斜[18]。 气候类型多样,从南向北包括了亚热带季风气候、温带季风气候、温带大陆性(干旱)气候和高原高寒气候等四大气候类型,年平均气温0~15 ℃,全省各地年降水量在36. 6~734. 9 mm,大致从东南向西北递减[19]。 全省土地总面积达4. 3×105km2,南北宽530 km,东西长1 655 km,包括14 个地级市(自治州),共86 个县(市、区)。 从自然地理角度,甘肃省可大体分为河西走廊(酒泉、嘉峪关、金昌、张掖和武威市)、陇中高原(兰州、白银、临夏回族自治州和定西市)、陇南山地(陇南和天水市)、陇东高原(庆阳和平凉)和甘南高原5 个自然资源区[20]。

1.2 数据源与处理

MCD19A2 数据是美国国家航空航天局(National Aeronautics and Space Administration,NASA)发布的Terra & Aqua MAIAC Land Aerosol Optical Depth Daily 1 km V006 产品,空间分辨率为1 km,波长为550 nm。 基于谷歌地球引擎(Google Earth Engine,GEE)对数据进行去云、重投影、镶嵌和裁剪等操作,在此基础上通过均值合成法生成月平均值(剔除无效值和缺失值),再通过月平均AOD 数据计算出季节平均AOD 值和年平均AOD 值。 在研究季节变化时,为避免某一年份出现季节间变化的偶然性,季节划分为春季(3—5 月)、夏季(6—8 月)、秋季(9—11月)和冬季(12—次年2 月)。

2 方法

BHM 包括数据模型、过程模型和参数模型3 个基本子模型[21],本文使用BHM 从空间格局、时间趋势和局部变化3 个方面以县域为空间子单元对AOD 值进行分析[22]。 假设在研究区内,研究对象服从正态分布,将模型精简为:

式中,yit为i地区t年AOD 对应的数值;θyit为对估计出yit的均值;σ2为AOD 的方差;α为研究时间段内该区域AOD 的恒定系数;si为研究时间范围内i地区AOD 相对变化的程度参数;b0it+vt表示AOD变化的整体趋势,b0it表示时间趋势中的线性趋势,vt表示时间趋势中的非线性趋势;b1it表示空间子单元在研究时间段内的AOD 局部变化趋势;b1i为局部偏离整体趋势的参数;εit为高斯噪声随机变量。si与b1i的先验分布由条件自回归模型(Conditional Autoregressive Models,CAR)建模确定[23],b0的先验分布服从均匀分布,vt的先验分布服从N(0,σv2),误差项εit的先验分布服从 N (σε2 )。 参照Gelman[24]的研究,模型中涉及随机变量的均方差(σv,σε)的先验分布服从正值半高斯分布。

为了评价BHM 对研究对象时空变化过程的解释程度,引入方差成分系数(VPC)进行评估[25],计算式为:

3 结果与分析

基于BHM,获得研究区内AOD 的整体时间趋势、空间格局和局部变化。 BHM 模型估算结果的VPC 均达到92.4%以上,说明对时空变化解释程度较好。

3.1 AOD 时间变化特征

3.1.1 AOD 总体格局

2001—2020 年甘肃省AOD 二十年平均分布如图1 所示。 分析可得,AOD 分布高的地区主要集中于河西北部、陇中和陇东等地,河西南部、甘南和陇南等地AOD 分布则相对较低。 河西西北部和陇中北部为显著的高值中心或高值聚集区,AOD 高值区人口密度大,汽车尾气排放量大,并且产业集中导致工业生产中能源燃烧排放大量尾气,人为产生气溶胶较多[26]。 陇中和陇东地区为黄土高原地带,河西北部荒漠带土地沙漠化严重,2 个地区AOD 水平也相对较高。 河西南部以及甘南、陇南地区地形多为高原、山地,人口相对较少、工厂少,AOD 水平较低。

图1 2001—2020 年甘肃省AOD 多年平均分布Fig.1 Multi-year average distribution of AOD in Gansu province from 2001 to 2020

3.1.2 AOD 年际变化趋势

依据贝叶斯层次时空模型中exp(b0t+υt)的后验中值表示AOD 整体年际变化趋势如图2(a)所示,甘肃省AOD 年均值变化趋势如图2(b)所示。图2(a)与图2(b)变化趋势基本一致,进一步说明了模型的适用性。 以2010 年和2013 年为转折点分为3 个阶段。 第1 个阶段为2001—2010 年,这个时间段内趋势系数有微小的上下起伏,2001—2004 年小幅度上升,2004 年之后逐渐下降后又上升,这一阶段整体呈现上升的趋势。 第2 阶段为2010—2013 年,这个阶段变化显著,先明显下降后再大幅提高,整体呈现增长趋势,其中2010—2012 年,趋势系数大幅度下降,2012—2013 年趋势系数骤增。2010—2013 年,甘肃省天气状况发生了明显的变化,2012 年甘肃省大风天气与前几年比相对较少,降雨较多,干旱影响造成的影响小,且年沙尘暴日数和扬沙日数均为近52 年最少,这是2012 年趋势系数最低的主要原因。 2013 年3 月出现连续高温天气,近54 年同期最高,甘肃省大面积出现长时间沙尘天气,全省81 个气象站中有73 个气象站监测到沙尘天气,导致2013 年趋势系数大幅度上升并达到20 年中最高值。 第3 个阶段为2013—2020 年,趋势系数先快速后缓慢降低,这与国家出台生态环境保护政策有重要关系,2020 年趋势系数为0.380,处于较低水平,主要是受全球新冠肺炎疫情的影响,多家工厂和企业停产,排放到大气中的污染物大幅减少,空气质量得到有效改善。

图2 甘肃省AOD 年际变化趋势贝叶斯时空模型估计结果及2001—2020 年AOD 年际变化趋势Fig.2 Estimation results of Bayesian spatio-temporal model of AOD interannual variation trend and interannual variation trend of AOD from 2001 to 2020 in Gansu province

3.1.3 AOD 季节变化趋势

对多年四季AOD 求均值进行AOD 季节变化趋势分析,如图3 所示,保证了准确性和可靠性。 由于春季是西北地区沙尘天气高发期,西北地区不同区域年均AOD 在春季显著升高,AOD 季节变化整体呈现出春季整体高,夏季、秋季局部中心高,冬季水平低的变化特征。 其中,春季甘肃省AOD 水平较高面积大,高值中心处AOD 水平极高,达到四季中最高水平;夏季和秋季AOD 水平较春季有所降低,冬季AOD 水平整体较低。

图3 2001—2020 年甘肃省多年平均四季AOD 值分布图Fig.3 Distribution map of annual average AOD values in Gansu Province from 2001 to 2020

3.2 AOD 空间变化特征

3.2.1 AOD 空间格局

AOD 分布(Si)如图4(a)所示,存在明显差异。从总体上看,空间格局为“东高西低、中部过渡”。AOD 高于均值的地区主要分布在河西西部(酒泉)、陇中、陇南和甘南东部地区,低于均值的地区主要分布在河西中部区域以及陇东地区。

为了进一步明确AOD 空间格局(Si)的冷热区域聚集情况,利用exp(Si)>1 的后验概率,将其划分为热-温-冷点区域,分别代表不同县域AOD 变化程度与研究区域AOD 变化程度高低水平的差异。

基于Richardson 提出的分类准则[27],若p(exp(Si)>1|data)>0.8,划分为热点区域;若p(exp(Si)>1|data)<0.2,划分为冷点区域;若0.21|data)<0.8,则划分为温点区域。 冷热区域聚集情况如图4(b)所示。 AOD 冷、热点区域在空间分布上有显著差别,冷点区域集中分布在河西地区(除酒泉市)、陇东和甘南西部;热点区域主要分布在陇中、陇南和甘南一带,例如兰州、天水、白银、临夏、定西、陇南和甘南东部地区,此外,酒泉和平凉市也存在少数热点区域,AOD 高热地区主要为城市群集中地,大面积的工业生产和频繁的人类活动会导致AOD 的大面积增长;温点区域集中分布在酒泉、嘉峪关、金昌也有少量分布。

图4 AOD 空间分布格局及冷热区域聚集情况Fig.4 Spatial distribution pattern of AOD and clustering situation in hot and cold regions

3.2.2 AOD 局部变化趋势

局部区域甚至单个的空间统计单元的局部变化趋势与总体变化趋势不同,局部空间区域上的变化趋势(即局部变化趋势)更能体现出时空变化过程中稳态结构之外的时空变化规律[28]。 以(b0+b1i)的后验中值表示2001—2020 年AOD 的局部区域变化情况,在考虑空间维度上的影响的同时,也考虑了时间过程中的影响。 若(b0+b1i)的值为负数,表明该空间统计单元上AOD 呈下降趋势,b1i<0(>0),说明局部区域变化强度高于(低于)总体趋势;若(b0+b1i)的值为正数,表明该空间统计单元上AOD 呈上升趋势,b1i>0(<0),说明局部区域变化强度高于(低于)总体趋势。

局部变化(b0+b1i)和局部变化偏离总体变化的强度(b1i)分别如图5(a)和图5(b)所示。

图5 局部变化和局部变化偏离总体变化的强度Fig.5 Local change and the intensity of local change deviation from overall change

由图5(a)可以看出,AOD 分布情况均有上升,整体来看,白银、兰州、临夏回族自治州、甘南藏族自治州和酒泉上升明显,其他地区变化较小。 由图5(b)可以看出,河西中部(除酒泉、嘉峪关)、陇东和甘南部分地区的变化强度低于整体趋势,陇中、甘南以及河西西部(酒泉、嘉峪关)的增长速度高于整体趋势。 具体而言,张掖、金昌、武威、庆阳、平凉以及甘南大部分地区变化强度低于整体趋势,呈现弱增长趋势;酒泉、嘉峪关、白银和临夏回族自治州变化强度高于总体趋势且增长速度较明显,其余地区的增长速度则较为平稳。

为了更好地描述局部变化冷热区域聚集情况,以b1i>0 的后验概率对局部AOD 变化进行分类。 若p(b1i>0| data)>0. 8,划分为热点区域;若p(b1i>0|data)<0.2,划分为冷点区域;若0. 20|data)<0.8,则划分为温点区域,结果如图6 所示。 局部变化冷点区域主要分布在金昌、武威、天水和陇东部分地区;热点区域集中分布在陇中和甘南地区;温点区域主要分布在河西、陇南一带,此外,陇中、陇东以及甘南地区也存在少数温点区域。

图6 局部变化热点分布Fig.6 Distribution of local change hot spots

3.2.3 AOD 空间格局和局部变化联合分析

空间格局总体变化趋势只是时空变化的一个概括性描述,局部区域可以得到空间统计单元上的趋势变化,将空间格局与局部变化联合分析可以更清楚地描述该区域AOD 未来发展趋势,因此将图5(b)和图6进行叠加分析可知,陇中和甘南部分地区不仅是AOD 空间格局的热点区域,而且是局部变化热点区域,表明这2 个地区AOD 水平一直较高,并且未来有逐渐升高的趋势。 天水和平凉部分地区属于AOD 空间格局的热点区域,同时处于局部变化冷点区域,即在研究期内,该地区的AOD 水平较高,但未来增长速度有所减缓,AOD 水平逐渐降低。 武威西部、庆阳大部分地区处于AOD 空间格局的冷点区域,同时处于局部变化冷点区域,表明该地区的AOD 水平较低且在未来一段时间内AOD 值持续降低。 酒泉的大部分地区处于AOD 空间格局和局部变化中都处于温点区域,说明该地区未来AOD 水平逐渐升高。 总的来讲,甘肃省未来在陇中、甘南和酒泉的大部分地区AOD 水平存在升高的趋势,其余地区AOD 水平趋于稳定。

4 结束语

本 文 采 用 GEE 提 供 的 2001—2020 年 的MCD19A2 数据,基于BHM 对甘肃省AOD 进行时空演变分析,得出以下结论:

① 从年际尺度来看,2001—2020 年AOD 趋势系数维持在0. 3~0. 5,2001—2010 年AOD 呈现上升趋势,2010—2013 年年际变化趋势显著,呈“V”形,2013—2020 年趋势系数逐渐降低,AOD 呈现逐年降低趋势。

② 从季节尺度上看,AOD 季节变化呈现出春季整体高,夏季、秋季局部中心高,冬季水平低的变化特征。 春季和冬季AOD 高值中心明显,主要分布在敦煌、兰州和白银地区。 四季AOD 最低值均在河西走廊南部祁连山地区。

③ 从空间格局上看,甘肃省东南和西北地区AOD 分布水平较高,其他地区相对较低。 冷、热点区域的空间分布差异明显,热点区域集中分布在甘肃省陇中、甘南和陇南地区。

④ 局部AOD 分布水平存在明显的区域差异,整体呈现以陇中和河西西部地区为中心向四周降低的趋势,河西西部、陇中以及甘南地区的变化强度高于整体变化,低于整体变化的地区主要分布在河西东部、陇南和陇东一带,局部变化热点区域主要分布在陇中和甘南的部分地区。