孙二昌 , 麻金继 *, 吴文涵 , 杨光 , 郭金雨

( 1 安徽师范大学地理与旅游学院, 安徽 芜湖 241002;2 资源环境与地理信息工程安徽省工程技术研究中心, 安徽 芜湖 241002 )

0 引 言

PM2.5污染危害人体健康[1]、降低大气能见度[2]以及对气候变化[3,4]有着重要影响。PM2.5污染的准确模拟不仅可以为公众提供及时的预报预警, 而且能在大气污染防控中起到重要作用。

大气化学模式基于天气动力学与大气化学基础, 能够重现大气污染物在三维空间中的反应与传输过程,是 PM2.5数值模拟的重要手段[5]。然而, 由于 PM2.5污染模拟的复杂性以及模式背景场 (初始场)、污染排放源等的不确定性[6-9], PM2.5污染模拟存在着偏差, 特别是重污染期间模拟偏差更大。资料同化是结合观测资料对模式背景场 "真实" 状态的一种估计方法, 是改善大气污染物预报精度的重要途径[10-12]。以往PM2.5污染模拟的同化资料多为地面观测资料或极轨卫星气溶胶光学厚度 (AOD) 资料。例如: Feng等[13]、Werner 等[14]通过地面观测资料与中分辨率成像光谱仪 (MODIS) AOD资料的同化实验证明了观测资料同化可减少背景场气溶胶的不确定性, 进而可改善 PM2.5浓度的预报精度。

新一代地球静止轨道卫星Himawari-8 发射于2014 年12 月7 日, 搭载了先进光学传感器AHI (Advanced Himawari Imager), 观测区域覆盖了我国大部分地区, 与极轨卫星相比拥有更高的观测频次, 为大气污染物模拟同化提供了更加丰富的观测资料[15]。Zhang等[16]、韦海宁等[17]、赵垒等[18]使用AERONET站点观测数据对中国区域Himawari-8 AOD 产品进行验证, 结果表明Himawari-8 AOD 产品在中国区域有较好的精度。Hong等[19]使用初始场优化同化方法评价了Himawari-8 AOD 产品估算的近地面 PM2.5资料同化对 PM2.5污染预报的影响, 结果表明由AOD 估算得到的 PM2.5资料同化能够改进 PM2.5污染模拟的精度。地球静止轨道卫星AOD资料同化研究目前处于发展阶段[20,21], AOD资料的同化方法需要进一步探讨。

本文使用Himawari-8 AOD产品与气象再分析数据, 基于地理加权回归模型 (GWR) 估算了近地表 PM2.5浓度, 利用WRF-Chem大气化学模式与GSI三维变分同化系统对 PM2.5污染模拟进行同化实验, 进一步探讨了高时间分辨率、高空间覆盖率的地球静止轨道卫星AOD资料在 PM2.5污染模拟同化中的应用。本文研究成果能为卫星AOD资料同化相关应用提供技术支撑。

1 数据与方法

1.1 研究区域与数据

2019年12月7―10日, 华北中南部、黄淮、江淮北部、汾渭平原、江汉和四川盆地出现大范围 PM2.5污染,其中河南、山东西部和河南南部均达到重污染程度。本研究基于WRF-Chem大气化学模式模拟了7―10日我国中东部地区 PM2.5污染, 模拟区域范围为 100° E～130° E, 20° N～40° N。

地面观测 PM2.5数据来自中国环境监测总站的全国城市空气质量实时发布系统。监测站连续自动监测系统有严格的质控技术规范与日常运行维护要求, 数据质量可靠。Himawari-8 AOD产品来自Level 3产品,覆盖时间12月7―10日, 空间分辨率为0.05°, 由Level 2产品经过严格的质量控制后合成。在使用观测数据前, 对数据进行了极值检查与异常值检查, 剔除了无效值与异常值。

1.2 PM2.5 资料同化

基于地理回归加权模型 (GWR), 使用气象再分析场数据 [包括大气边界层高度 (PBLH)、大气相对湿度(RH)、风场 (WINDS)] 与Himawari-8 AOD产品构建了近地面 PM2.5浓度估算模型[22], 从而进行近地表 PM2.5浓度的估算。

PM2.5污染模拟通过WRF-Chem模式实现。WRF-Chem模式考虑了大气污染物与气象双向反馈作用, 具有 "在线耦合" 的特点[23]。其中, 气溶胶模块选用GOCART[24], 气溶胶组成包括未分类的细粒子气溶胶 P2.5、硫酸盐 (Sulfate, 有效粒径为0.242 μm)、黑碳气溶胶 (BC1、BC2, 有效粒径为0.036 μm)、亲水性和疏水性有机质气溶胶 (OC1、OC2, 有效粒径为0.087 μm) 以及五档不同粒径的沙尘气溶胶 (D1: 0.2～2.0 μm; D2: 2.0～3.6 μm; D3: 3.6～6.0 μm; D4: 6.0～12.0 μm; D5: 12.0～20.0 μm) 与四档不同粒径的海盐气溶胶 (SS1: 0.2～1.0 μm; SS2: 1.0～3.0 μm; SS3: 3.0～10.0 μm; SS4: 10.0～20.0 μm)。PM2.5的浓度cPM2.5为

式中ρd为干空气密度。

PM2.5观测资料 (包括由AOD估算的PM2.5资料) 的同化使用GSI三维变分同化系统。三维变分同化方法可归结为求解函数极小值问题[25], 代价函数J的一般形式为

式中X表示分析场,Xb表示背景场,B-1表示背景场误差协方差矩阵的逆,Y表示观测数据,H表示观测算子,R-1表示观测场误差协方差矩阵的逆。

在GSI同化系统中, 地面观测 PM2.5资料同化的观测算子见式 (1)。地面观测不包含 PM2.5污染的垂直廓线信息, 所以仅同化近地面污染物浓度。地面 PM2.5观测资料的观测误差由站点自身观测误差以及插值距离决定[26], 计算公式为

式中εm为仪器误差,α为可调节参数, Δx为模拟的网格分辨率,Lrepr为观测影响的半径。

在GSI三维变分同化系统中, PM2.5观测浓度上限一般设为100 μg/m3, 考虑到此次污染过程 PM2.5观测浓度远大于此, 所以调整 PM2.5观测浓度上限至400 μg/m3。

背景场误差通过使用国家气象中心 (NMC) 方法统计获得。首先模拟2019年12月7―10日每日03:00和06:00 (UTC时间) 的48 h预报和24 h预报, 然后使用GEN_BE程序[27]计算得到。NMC方法[28]的计算公式为

式中Xb为背景场,Xt为真实场,X48h、X24h分别为对同一时刻提前48 h、24 h预报。

1.3 实验设计

为探究地球静止轨道卫星Himawari-8 AOD产品估算的近地面 PM2.5资料同化对 PM2.5污染模拟的影响,设计了四组 PM2.5污染模拟实验 (参见表1)。第一组实验为控制实验 (CONTROL), 用于观测资料同化效果的对照, 其余三组实验分别为AOD估算资料同化实验 (AODDA)、地面观测资料同化实验 (SITEDA)、同时同化地面观测资料与AOD估算资料实验 (BOTHDA)。在所有的实验中, WRF-Chem模式的模拟区域与物理化学方案设置相同, 模拟时间为2019年12月7―10日。

表 1 同化实验方案Table 1 Assimilation experiments scheme

在PM2.5污染模拟实验中, 相关参数设置为: 投影为Lambert, 区域中心坐标为 (110° E, 32° S), 网格尺度为27 km × 27 km, 水平网格个数为108 (东西方向) × 98 (南北方向), 垂直网格个数为40, 侧边界条件设置为理想大气成分。Spin-up 时间设置为前5 天。化学参数: 化学机制为RACM (Regional atmospheric chemistry mechanism) 方案, 光解作用方案为TUV。物理参数: 微物理方案为WMS5, 边界层方案为YSU, 积云对流参数化方案为Grell3D, 长波辐射方案为RTTM, 短波辐射方案为Goddard, 陆面过程方案为Noah。污染物排放清单数据为中国多尺度排放清单模型MEIC排放清单数据 (http://www.meicmodel.org/)。气象资料为美国国家大气中心 (NCAR) 提供的时间分辨率为6 h、空间分辨率为0.25° × 0.25°的全球天气再分析资料。

在资料同化实验中, 采用的是循环同化方法。同化时间窗口为3 h。2019年12月7日03:00同化后, 向后预报3 h作为7日06:00的同化背景场; 7日6时同化后, 向后预报20 h做为8日03:00的同化背景场, 一直循环同化预报至2019年12月10日结束。

1.4 同化效果评价方法

为定量评价资料同化对 PM2.5污染模拟的影响, 选取了绝对平均误差 (MAE,EMA)、均方根误差 (RMSE,ERMS) 以及相关系数r三个指标用于同化效果的评价, 指标的计算公式分别为

2 结果与讨论

2.1 近地面 PM2.5 估算结果

图1 为2019 年12 月7―10 日 Himawari-8 AOD 产品估算的近地表 PM2.5浓度空间分布情况。由图可知,AOD估算的近地面 PM2.5具有空间覆盖范围广、聚集性高的特点, 可为 PM2.5模拟同化提供更为丰富的 PM2.5空间分布信息。从图1 中各分图的比较中可以看出, 估算的 PM2.5浓度具有明显的时空差异, 7日华北地区、华东地区北部、中南地区北部是本次污染中 PM2.5的高值区, 随着时间的推移, 8―10日 PM2.5污染向北推移。这一现象与徐冉等[29]的分析结果一致。为进一步证实 PM2.5浓度估算结果的可靠性, 将地面观测 PM2.5数据与估算结果进行比较。图2为2019年12月7―10日估算的近地面 PM2.5浓度与地面观测 PM2.5数据的比较结果。由图可知, 7―10 日估算的 PM2.5相关系数r均值为0.85, MAE 与RMSE 均值分别为16.36 μg/m3、27.07 μg/m3, 与同类研究的估算水平相当[30], 表明估算结果是可靠的。

图1 Himawari-8 AOD产品估算的近地表PM2.5浓度空间分布。(a) 12月7日; (b) 12月8日; (c) 12月9日; (d) 12月10日Fig. 1 The spatial distribution of PM2.5 concentration near the surface estimated by Himawari-8 AOD. (a) December 7;(b) December 8; (c) December 9; (d) December 10

图2 Himawari-8 AOD产品估算的近地表PM2.5浓度与地面观测PM2.5浓度的比较。(a) 12月7日; (b) 12月8日;(c) 12月9日; (d) 12月10日Fig. 2 Comparison of the near-surface PM2.5 concentration estimated by Himawari-8 AOD product and the surface-observed PM2.5 concentration. (a) December 7; (b) December 8; (c) December 9; (d) December 10

2.2 背景场PM2.5的改进效果

为探究不同实验组资料同化对PM2.5背景场的改进作用, 分别将实验组背景场PM2.5浓度与地面观测PM2.5浓度进行了对比。从图3 (a) 中可以看出, CONTROL实验中背景场PM2.5浓度与地面观测PM2.5浓度的线性回归方程为y= 0.38x+ 41.12,r为0.43, MA、RMSE 分别为38.01 μg/m3和55.61 μg/m3。图3 (b)中,AODDA 实验中的背景场PM2.5浓度与地面观测PM2.5浓度的线性回归方程为y= 0.57x+ 33.56,r为0.67,MAE、RMSE分别 27.79 μg/m3、42.29 μg/m3。AOD估算的PM2.5资料同化后,r较CONTROL实验提升了0.24,MAE、RMSE 分别减少了10.22 μg/m3和13.32 μg/m3。图3 (c)中, SITEDA 实验背景场PM2.5浓度与地面观测PM2.5浓度的回归方程为y= 0.71x+ 19.95,r为0.86, MAE、RMSE 分别为17.36 μg/m3、27.93 μg/m3。图3 (d)中, BOTHDA 实验背景场PM2.5浓度与地面观测PM2.5浓度的回归方程为y= 0.71x+ 20.80,r为0.87, MAE、RMSE分别为17.39 μg/m3、27.17 μg/m3。综上可知, SITEDA与BOTHDA实验对背景场PM2.5浓度的改进效果相当, 好于AODDA实验的改进效果, 这可能与近地面PM2.5的估算精度以及同化站点与验证站点的空间位置有关, 因为当估算的PM2.5资料混入地面观测PM2.5资料后, 同化结果的评价指标有了较大的改进。

图3 背景场PM2.5浓度与地面观测PM2.5浓度散点图。(a) CONTROL实验; (b) AODDA实验;(c) SITEDA实验; (d) BOTHDA实验Fig. 3 Scatter plot of PM2.5 concentration in background field and PM2.5 concentration observed on the ground. (a) CONTROL experiment; (b) AODDA experiment; (c) SITEDA experiment; (d) BOTHDA experiment

表2统计了各组同化实验在不同时刻同化效果的评价指标。由表可知, AODDA实验中对背景场PM2.5浓度MAE改进的均值为10 μg/m3, RMSE改进的均值为13 μg/m3,r提升区间0.11～0.39, 极大值出现在T6时刻(9日06:00); SITEDA 实验中, MAE、RMSE改进的均值分别为21 μg/m3、27 μg/m3,r改进区间为0.29～0.53, 提升极大值也出现在T6时刻; BOTHDA实验中出现了相同的现象。可以认为, 资料同化对背景场PM2.5的改进效果存在时间上的差异, 随模拟时间增加改进效果呈增加趋势, T6时刻为提升效果极大时刻。分析原因可能为: 在7―10日的这次PM2.5污染中, 徐冉等[29]指出自7日开始发展, 8日凌晨西北地区东部以及中南地区北部等地的PM2.5快速累积至重度污染, 并于8―10日发展成为全国性的大范围PM2.5污染天气, 至10日中午前后,受新一轮冷空气的影响, PM2.5污染开始消散。在这一污染发展过程中, PM2.5污染逐渐加重, PM2.5以及相关前体物质开始大范围累积, 在重污染条件下相关部门往往会对部分污染排放进行人为管控, 导致了排放清单数据的偏差, 致使PM2.5模拟精度下降。这一点可从CONTROL实验的MAE、RMSE增加以及r指标减小的趋势看出。因此可认为, 在PM2.5污染过程的不同发展阶段, 同化对背景场改进效果总体上表现为随模拟偏差的增大而增大。

表2 2019年12月7―10日背景场PM2.5同化评价指标统计与评价指标统计差值Table 2 The difference between the PM2.5 assimilation inspection index and the verification index in the background field on December 7―10， 2019

2.3 预报场PM2.5改进的空间差异

在各组PM2.5污染模拟实验中, 每日07:00至次日02:00, 即下一次同化时刻前为PM2.5预报时段。为评价同化实验中预报场PM2.5的改进效果, 将实验中PM2.5模拟浓度与地面观测PM2.5浓度进行逐站点比较。图4给出了对比实验中PM2.5逐站点的MAE、RMSE 以及r统计情况。依据图4 (a)、(b), PM2.5浓度模拟结果的偏差存在空间差异。其中, 中南地区南部、华东地区中南部以及西北地区中部PM2.5浓度模拟结果的MAE、RMSE相对较小, 偏差均小于50 μg/m3, 对于偏差低值聚集区而言, 区域内部估计指标差异较小; 西南地区、中南地区北部以及华北地区模拟偏差较大且出现 "扎堆" 现象, MAE、RMSE 检验指标均在70～170 μg/m3范围;PM2.5模拟值与观测值r较好的站点主要出现在中南地区、西南地区以及西北地区东部,r检验指标大于0.6;华东地区东北部以及西北地区中部的相关性较差。分析CONTROAL实验中PM2.5模拟偏差出现空间差异的原因, 可能与PM2.5复杂的反应过程及前体污染物、污染排放源、气象条件等因素的不确定性有关。

以CONTROL实验中模拟偏差较大的西南地区、中南地区北部以及华东地区北部为例, 分析资料同化对PM2.5预报的改进作用。从图4 (d)、(e)中可以发现, AOD估算的PM2.5资料同化对PM2.5预报的偏差改进区间在1～6 μg/m3。通过比较图4 (g)、(h) 与图4 (j)、(k) 可以发现, SITEDA与BOTHDA实验中PM2.5偏差的改进效果相近。依据图4 (j)、(k)可以发现, 同时同化AOD估算的PM2.5资料与地面观测的PM2.5资料对PM2.5的MAE与RMSE偏差改进区间分别为1～17 μg/m3和 1～18 μg/m3, RMSE改进高于3 μg/m3的站点较AODDA实验中的数量增加了一倍左右。

综合比较图4 (d)、(e)、(g)、(h)、(j)、(k) 可以得出, 同化实验在CONTROL实验中PM2.5模拟偏差较大的地方对模拟偏差MAE、RMSE进效果会更好。具体表现在, CONTROL实验中PM2.5浓度模拟偏差较大的西南地区、中南地区北部以及华北地区是同化改进效果最好的区域。同样地, 中南地区南部、华东地区中南部以及西北地区中部在CONTROL实验中PM2.5模拟偏差MAE、RMSE相对较小, 同化对这些区域PM2.5模拟偏差的改进效果也较弱。同化对PM2.5预报的相关性改进在空间分布也呈现类似的特点。由图4 (f)、(i)、(l) 可知,CONTROAL实验中华东地区东北部是PM2.5模拟相关性较差的地区, 也是同化后r指标改进最好的地区。

图4 MAE、RMSE、r 指标空间分布。(a)―(c) CONTROL实验; (d)―(f) AODDA实验; (g)―(i) SITEDA实验;(j)―(l) BOTHDA实验Fig. 4 The spatial distribution of MAE, RMSE and r indicators. (a)―(c) CONTROL experiment; (d)―(f) AODDA experiment;(g)―(i) SITEDA experiment; (j)―(l) BOTHDA experiment

将同化对PM2.5污染模拟的改进效果与其他学者[31]的研究结果进行了比较。在RMSE、r指标上, 改进效果相近。但在西南地区, 其他学者同化结果表明RMSE 改进的数值在10～30 μg/m3范围, 而在本研究中BOTHDA实验的提升区间在4～17 μg/m3。综上所述, 观测资料同化对PM2.5预报场的改进作用积极, 三组同化实验改进的空间分布上较为接近。

2.4 预报场PM2.5改进的持续性

图5为同化实验预报改善效果的时序比较结果。其中横轴为预报时长, 0～20 h表示12月7日06:00同化后至第12月8日03:00同化前的预报时段, 20～40 h表示12月8日06:00同化后至12月9日03:00同化前的预报时长, 以此类推。从图5中各项评价指标曲线的变化趋势可以看出, 各组同化实验同化后的模拟偏差比参照实验的偏差更低, 相关性也有所提升, 改进时长可持续至两次同化间隔的20 h。通过比较AODDA、SITEDA以及BOTHDA曲线的变化可以发现, BOTHDA与SITEDA实验中对PM2.5预报的改进能够持续相对更长的时间。AODDA、SITEDA 以及BOTHDA 实验MAE 改进的极大值分别为 2.9、5.1、5.0 μg/m3, RMSE 极大值分别为 4.1、7.1、7.1 μg/m3。从预报改进的均值来看, 三者的MAE分别为 1、1.8、2.1 μg/m3, RMSE分别为1.5、2.7、3.1 μg/m3。

图 5 同化实验预报改善效果的时序比较。(a) MAE; (b) RMSE; (c) rFig. 5 Time series comparison of assimilation experiment forecast improvement effect. (a) MAE; (b) RMSE; (c) r

从PM2.5预报的不同时段来看, 预报效果改进幅度随模拟时长的增加而增加, 12月9―10日同化改进效果要好于7―8日的改进。随着PM2.5模拟时长的增加, 预报模式积累的误差也随之增加, 表现为参照实验中的偏差指标曲线呈上升趋势, 相关性指标数值降低, 但同化改进效果却随着误差积累增加而增加。可以认为,观测资料同化在PM2.5模拟偏差较大的时段对PM2.5预报的改进效果更好。

2.5 同化对垂直方向PM2.5的影响

PM2.5浓度不仅受到水平方向传输影响, 而且受到垂直方向上的影响。Jiang等[32]指出, 垂直方向的传输是影响同化效果的重要因素。为探究近地面观测资料同化对垂直方向PM2.5浓度的影响, 对不同实验组垂直方向的PM2.5浓度进行了统计分析。图6 (a) 描述了参照实验中模拟的PM2.5垂直方向整层均值随时间变化情况, 横轴的含义与图5中的相同, 纵轴为垂直层高度。从图6 (a) 中可以看出, 整层PM2.5模拟均值在垂直方向随高度增加而减少, 模拟的PM2.5主要分布在高度1 km 以下, 浓度均值在15～55 μg/m3, 高于1 km 且小于3 km 的PM2.5浓度均值在1～10 μg/m3, 高于3 km 的PM2.5浓度均值小于1 μg/m3。图6 (b)、(c)、(d) 分别为AODDA、SITEDA 以及BOTHDA 实验PM2.5模拟的垂直方向整层PM2.5均值与CONTROL 实验的差值。从近地面PM2.5资料同化对垂直方向PM2.5浓度模拟结果的影响来看, 观测资料同化对模拟的近地面PM2.5浓度影响最大, 在5～7 μg/m3范围; 对近地表至垂直高度2 km范围而言, 影响小于2 μg/m3; 高于2 km范围影响可以忽略不计。从同化对垂直方向模拟的PM2.5浓度影响的程度来看, AODDA与BOTHDA实验对PM2.5垂直方向的影响范围要大于SITEDA 站点资料同化。总结来看, 近地表PM2.5观测资料的同化对不同垂直高度上的PM2.5浓度的分布存在一定的影响, 不同垂直高度的PM2.5存在相互影响, 同化不同垂直高度上的PM2.5对近地表PM2.5模拟改进是有益的。

图 6 同化前后PM2.5垂直方向整层均值的时序比较。(a) 参照实验; (b)参照实验与AODDA的差值; (c) 参照实验与SITEDA差值; (d) 参照实验与BOTHDA差值Fig. 6 Temporal series comparison of the vertical whole layer mean of PM2.5 before and after assimilation. (a) Reference experiment; (b) difference between reference experiment and AODA; (c) difference between reference experiment and SITEDA; (d) difference between reference experiment and BOTHDA

3 结 论

使用三维变分同化系统 (GSI) 同化了由Himawari-8 AOD产品估算的近地面PM2.5资料, 设计控制实验与同化实验并结合地面观测PM2.5数据对同化效果进行了评价。主要结论: 1) 估算的PM2.5资料同化对PM2.5污染模拟的改进作用明显, 对背景场PM2.5的绝对平均误差 (MAE)、均方根误差 (RMSE) 分别改进了10.22 μg/m3、13.32 μg/m3,r较CONTROL实验提升了0.24; 对预报场PM2.5的MAE、RMSE改进区间在1～6 μg/m3;对预报场持续改进方面, MAE 与RMSE 改进的极大值分别为2.9 μg/m3和4.1 μg/m3, 改进均值分别为1 μg/m3和1.5 μg/m3。2) 资料同化对PM2.5模拟改进效果存在时空差异, 在西南地区、中南地区北部以及华东地区北部等CONTROL实验中模拟效果表现较差的地区改进效果更加显著; 9―10日模拟偏差较大的时段是改善效果最明显的时段。较其他研究中使用AOD观测算子直接同化AOD资料的方法, 同化AOD估算的PM2.5资料避免了AOD资料与硫酸盐、黑碳气溶胶、有机质气溶胶等PM2.5组成物质之间复杂的观测算子推演过程, 使得AOD资料的同化变得简单可行。

研究还发现, 影响同化AOD估算的PM2.5资料对PM2.5污染模拟改进效果的因素包括近地面的PM2.5的估算精度与PM2.5的垂直廓线信息。一方面可以使用深度学习、物理模型加深度学习的方式提高AOD对近地面PM2.5浓度的估算精度; 另一方面可以同化PM2.5的垂直维信息, 以改进PM2.5的模拟效果。