多元线性回归方法对北京地区PM2.5预报的改进应用

2019-05-15潘锦秀晏平仲李云婷刘保献王占山

中国环境监测 2019年2期

潘锦秀，晏平仲，孙峰，李云婷，刘保献，王占山，董瑞

1.北京市环境保护监测中心,大气颗粒物监测技术北京市重点实验室,北京 100048 2.中国科学院大气物理研究所,大气边界层和大气化学国家重点实验室,北京 100029

近年来，在《北京市2013—2017年清洁行动计划》[1]和京津冀及周边地区大气污染联防联控政策的积极推动下，北京市PM2.5浓度逐年降低，从2013年89 μg/m3降至2017年58 μg/m3[2]，空气质量在一定程度上得以改善。但是，秋冬季仍然是空气重污染频发季节，以2016年为例，38 d PM2.5重污染中有22 d发生在采暖季，38 d 重污染对全年PM2.5平均浓度贡献为23 μg/m3(占31.5%)[3]。北京市的空气质量重污染过程形成受到气象、排放和地理条件等诸多因素影响[4-6]，重污染过程的发生有时无法避免，要达到年均值逐年下降的目标，重污染过程中污染物的“削峰降速”作用就显得尤为重要。北京市PM2.5浓度变化剧烈，在极短的时间内容易出现浓度快速上升或下降的巨大变化，在2017年新修订的《北京市空气重污染应急预案》中再次引入了“小时”概念[7]，这意味着对重污染的实际业务预报需要精确到小时级别，对于其起始时间、高值持续时段及污染的级别等均要做出相应的判断，为政府相关部门争取足够的时间采取减排措施，降低污染程度。

日常的空气质量预报往往是在业务化多模式空气质量预报系统[8]结果的基础上加以人工订正，从而对未来的空气质量进行预报，但是目前空气质量模式多样，预报结果具有多样性，应用时往往难以抉择，这就需要开发一个优于各个模式结果的预报产品来为我们服务。集合预报技术主要基于复杂的三维环境空气质量数值模式，通过构建产生多个具有差异的预报样本，利用多元回归、神经网络等数学方法产生最优确定性预报结果，并且可提供污染发生概率预报，同时定性定量表征区域污染程度和污染规律，为环境空气质量预报预警和污染控制决策支持提供更为丰富的预报信息[9-11]。王茜等[12]采用学习型线性回归方法对上海市PM2.5预报效果进行改进，明显提高了PM2.5浓度的预报效果。吴剑斌等[13]在臭氧业务预报中引入一种最优化集合方法，提高了臭氧预报的准确率。

目前，北京市多模式空气质量预报系统对PM2.5预报效果主要依赖不同模式成员的预报效果，并没有将观测信息有效应用于空气质量预报的优化，因此，该研究首先评估了多模式空气质量预报系统对2016年北京市PM2.5的预报技巧，并利用多元线性回归方法将一年的日均PM2.5浓度值及2016年12月16—22日(红色预警期间)的小时PM2.5浓度值与三模式的结果进行集成，得到优于各模式预报特点的集合预报方法，从而提高多模式集合预报系统对日常业务化空气质量预报及重污染天气预报预警的支撑力度。

1 研究方法

数值模式数据来自中国科学院大气物理研究所研发的多模式预报系统，主要由CMAQ、CAMx和NAQPMS 3个模式组成。不同模式具有相同的气象场输入、区域设置和排放源，该多模式空气质量预报系统的相关详细设置见参考文献[14]。模式系统每日预报未来7 d京津冀区域13个地市空气质量状况，每日预报以前一日北京时间20:00为起始时间，预测未来7 d污染物浓度。该研究模式相关数据均来自模式系统第二层嵌套预报结果，提取每日预报的第5～28 h模式数据，从而获得2016年全年的三模式PM2.5小时预报数据。

1.1 数据资料来源

观测资料来自于北京市环境保护监测中心，包括北京2016年12个国控点位数据PM2.5的逐小时观测浓度，北京市国控点位的具体信息如表1所示。

表1 北京市国控点位信息Table 1 The information of air monitoring stations in Beijing

1.2 多元线性回归方法

多元线性回归方法(REG)的主要思路为根据起报日之前一定历史时期的模式预报效果，总结出各模式预报值与实测值之间的函数关系。黄思等[15]利用多元线性回归方法改进了2010年北京PM10的预报效果，在此基础上，该研究利用该方法改进日均PM2.5的业务化预报效果，重点研究了该方法对小时PM2.5预报的改进效果，以期支撑重污染应急措施的启动和解除，多元线性回归方法的具体公式见参考文献[15]。在利用多元线性回归方程集成时有一个关键参数，即训练时长，它表征着需要多长时间的历史观测资料才能更好地构建回归模式，不同的训练时长订正的预报效果不同。为了寻求最优的训练时长，该研究利用北京2016年全年国控点位平均PM2.5日均数据和三模式预报数据，统计了训练时长为1～70 d 的情况下，预报与观测的相关系数(r)和均方根误差(RMSE)的变化，如图1所示。

图1 各点位平均误差随训练时长的变化Fig.1 Average error of all sites at different lengths of training time

由图1可以看出，当训练时长小于13 d时，多元线性回归的r值低于0.65，RMSE高于50 μg/m3，这表示在训练时长低于13 d时，多元线性回归方法并不能有效提高PM2.5的预报技巧。随着训练时长的增加，r不断增加，RMSE不断减小，这表示随着训练时长的增加，多元线性回归方法对预报的改进效果显著，当训练时长达到64 d后，r和RMSE的变化处于较平稳的变化中，敏感度下降，增加训练时长对预报技巧的影响较小，因此，该研究将训练时长设置为64 d，对北京市PM2.5的预报结果进行集合预报。

1.3 评估方法

选取2016年1—12月作为评估时段，并划分为春季(3—5月)、夏季(6—8月)、秋季(9—11月)和冬季(1—2月和12月)，主要从时间序列分析和统计参数2个方面评估模式系统预报能力。

时间序列分析主要考察各模式对污染物浓度的整体变化趋势的预测能力，并在气象、环境的模式研究领域通用的统计指标中，选取r、RMSE、平均偏差(MB)和标准化平均偏差(NMB)等共4项指标对PM2.5预报效果进行评估，具体指标定义和计算公式见参考文献[16]。其中，r反映了观测值与预报值随时间变化趋势的相似程度，其值越大表明相关性越好；RMSE表示的是一组数据偏离平均值的离散程度，其值越小表示预报效果越好；MB表示预报值与观测值的差值，MB>0表示预报相对观测高估，MB<0则表示预报相对于观测低估，RMSE和MB的单位均为μg/m3；NMB克服了在观测浓度差异较大时难以直接比较的缺点，可以用来表征预报高估或者低估的程度。

2 结果与讨论

2.1 日均值预报效果评估

为了解多模式集合预报系统中不同模式的预报能力和特点，首先评估了不同季节多模式系统中3个模式的预报能力。其次，为了充分利用3个模式的预报结果，并将历史的观测数据纳入预报中，结合多元线性回归方法集成3个模式的预报结果，并对其效果进行评估。

图2～图4为3个模式和多元线性回归方法在不同季节对日均PM2.5浓度值预报的时间序列图，表2为参数统计。挑选定陵点位代表北部地区，万寿西宫点位代表城区，国控平均值代表全市基本状况对预报结果进行评估。3个模式对PM2.5的预报随时间的变化趋势基本一致，说明气象因素是污染物日际间变化的主因。就春季而言，相关系数方面，CMAQ和NAQPMS的预报效果略好于CAMx，这可能主要与模式内部的化学机制有关系。对国控平均和城区的预报效果较北部地区好，具体表现为r在国控平均和城区基本达到0.8以上，RMSE为40 μg/m3左右，各模式预报的MB基本为负值，表示预报值对实况有低估的趋势，低估主要表现在对重污染峰值的低估，各城区点位MB值为-25～-8 μg/m3，NMB为-30%～-15%，处于较好水平，表明低估现象不是很明显。而对北部点位的预报结果评估发现，r为0.7以上，略低于城区；RMSE为40～60 μg/m3；MB为-40～-20 μg/m3，NMB为-60%～-20%，表明对北部地区PM2.5的变化趋势预报和观测基本一致，但是有较明显的低估现象，其低估趋势较城区更明显。

就夏季而言，CMAQ、CAMx和NAQPMS对不同点位的预报也存在差异，r值较春季低，但也处于0.65以上的较好水平，RMSE低于春季，处于25～32 μg/m3，表明其离散程度较春季好，MB和NMB不同模式有正有负，低估和高估现象均不严重，NMB处于-20%～26%，对夏季污染过程峰值浓度的预报也较好。针对北部点位，绝大部分r值大于0.7，RMSE为20 μg/m3左右，对峰值浓度有一定的低估现象，但是不明显。

就秋季而言，国控平均和城区r值绝大部分大于0.7，RMSE处于35～60 μg/m3，MB和NMB显示了模式对秋季的预报存在高估现象，其中CMAQ和CAMx的高估较NAQPMS明显，但是三模式对污染的峰值浓度以及清除时间和清除后浓度把握准确，准确预测污染累积和清除过程的时间及峰、谷值浓度，对重污染过程预报有良好指导作用，对日常业务预报有较强的指导意义。对北部点位的预报趋势较好，但对峰值浓度有明显的低估现象，这可能是由于排放源对该地区的估计量不够所致。

就冬季而言，各模式对城区的预报趋势一致，能较好地反映PM2.5的日均变化，但是由于2016年1—2月的污染水平较轻，模式系统中的排放估计量过高而导致MB、NMB及RMSE均处于较高水平，高估现象很明显，而对污染水平较高的12月而言，无论从污染起始时间、结束时间及污染程度上均把握较好。冬季各模式对北部点位存在低估趋势，主要表现在1—2月峰值等预报较好，但是对于12月的污染过程存在明显的低估现象，这跟排放源的空间分布有关。

图2 不同季节国控平均PM2.5时间序列图Fig.2 The series of observed and simulated PM2.5 concentrarion of averaged sitesin different seasons

图3 不同季节定陵PM2.5日均值时间序列图Fig.3 The series of observed and simulated PM2.5concentrarion of Dingling sitein different seasons

图4 不同季节万寿西宫PM2.5时间序列图Fig.4 The series of observed and simulated PM2.5concentrarion of WSXG in different seasons

季节和预报方法国控平均定陵万寿西宫rRMSE/(μg/m3)MB/(μg/m3)NMBrRMSE/(μg/m3)MB/(μg/m3)NMBrRMSE/(μg/m3)MB/(μg/m3)NMB春夏秋冬CMAQ0.8346.4 -25.6 -0.360.7657.7 -34.6 -0.540.8642.7 -21.9 -0.30CAMx0.8143.8 -21.0 -0.290.7455.5 -31.2 -0.480.8638.7 -16.2 -0.22NAQPMS0.8741.4 -21.5 -0.300.7951.3 -27.8 -0.430.8936.1 -14.0 -0.19REG0.8130.6 -2.3 -0.040.7535.2 -2.8 -0.050.8733.0 -8.1 -0.12CMAQ0.7224.2 -9.1 -0.160.7625.8 -15.1 -0.310.6627.2 -9.8 -0.16CAMx0.7123.2 0.6 0.010.7423.1 -8.0 -0.170.6725.7 1.2 0.02NAQPMS0.7027.8 10.7 0.180.7126.0 8.5 0.180.6532.3 14.0 0.23REG0.8021.8 9.0 0.150.7922.6 9.4 0.190.7524.5 8.9 0.14CMAQ0.7541.3 -3.2 -0.040.7840.6 -19.2 -0.300.7046.1 1.8 0.02CAMx0.7740.1 7.8 0.100.7836.8 -11.5 -0.180.7150.2 17.0 0.20NAQPMS0.8136.3 1.8 0.020.6841.7 -9.5 -0.150.7942.8 15.8 0.19REG0.7839.1 -2.7 -0.030.7338.0 -4.4 -0.070.7542.8 0.4 0.01CMAQ0.7859.4 26.2 0.320.7649.3 -11.1 -0.180.7674.1 36.2 0.38CAMx0.7668.6 40.1 0.490.7647.5 -4.7 -0.080.7589.4 59.1 0.62NAQPMS0.8059.7 28.5 0.350.8147.5 -9.3 -0.150.73102.2 68.3 0.72REG0.7370.6 9.3 0.070.7241.5 0.3 00.7977.9 2.2 0.01

总体而言，各模式均能较好预测出PM2.5变化趋势及峰值出现时间，大部分点位预报的r值为0.6～0.9，NMB为-0.6～0.6，但是不同季节、不同地区还存在比较明显的高估或低估趋势，基于多模式预报和观测数据的多元线性回归方法能较好地解决这一问题。图2～图4还显示了多元线性回归方法集成的PM2.5随时间变化趋势，结果表明集成后夏季r在CMAQ和NAQPMS的基础上仍有所提高，而对于春季和秋季则改进效果不明显，冬季对万寿西宫的改进效果好于定陵，r和RMSE均有较高的提升。利用多元线性回归方法集成后对2016年春季国控平均低估现象改进明显，三模式平均MB由-23 μg/m3改善至-2.3 μg/m3，NMB由-0.32改善至-0.04，RMSE由44 μg/m3降至31 μg/m3；对2016年冬季国控平均高估也有较好的改进，各模式平均MB由32 μg/m3降至9 μg/m3，NMB由0.39降至0.07；对万寿西宫的改进尤为明显，三模式平均MB由55 μg/m3降至2 μg/m3，NMB由0.54降至0.01。

2.2 对2016年12月红色预警的改进效果

根据不同模式对PM2.5日均值预报效果的评估发现，在日常业务预报中，模式对于重污染过程的PM2.5日均值预报存在明显的低估或高估现象，多元线性回归方法虽然对模式预报的PM2.5日均值整体高估或者低估有一定的改进，但是针对重污染过程的起始浓度、峰值浓度并未有明显的改进效果，这关系到重污染预警的启动与解除。因此，利用多元线性回归方法对2016年12月北京红色预警期间不同模式的小时预报结果进行订正，来探讨其改善效果，此时选择的最优训练时长为58 h。

2016年12月16—22日，北京市发生一次连续5 d的严重污染过程。此次重污染过程持续时间长，共有122个小时浓度处于150 μg/m3重度污染级别限值以上，PM2.5峰值浓度达到434 μg/m3。图5为不同预报方法红色预警期间PM2.5随时间的变化，各预报方法的统计指标见图6～图8。

图5 不同预报方法红色预警期间PM2.5随时间的变化Fig.5 The variation of PM2.5of different forecast methods during red alert period of air pollution

图6 不同预报方法红色预警期间国控平均PM2.5散点图Fig.6 Scatter plots of the observed and simulated PM2.5 concentration for averaged sites of different forecast methods

图7 不同预报方法红色预警期间定陵PM2.5散点图Fig.7 Scatter plots of the observed and simulated PM2.5 concentration for Dingling site of different forecast methods

图8 不同预报方法红色预警期间万寿西宫PM2.5散点图Fig.8 Scatter plots of the observed and simulated PM2.5 concentration for WSXG of different forecast methods

研究结果表明，利用多元线性回归方法集合后，PM2.5小时浓度随时间的变化与实况更为吻合，改进后国控平均预报与观测PM2.5的r值由低于0.6提升至0.73，预报值与观测值的r提高了0.13，RMSE降低了20～30 μg/m3,低估现象也得到了一定的缓解，NMB由平均-0.15变为-0.13。对万寿西宫的改进效果不明显，r值低于三模式平均值，RMSE降低了4 μg/m3,低估现象缓解不明显。在定陵的效果最为明显，不同模式均大大低估了该地区的污染过程，利用多元线性回归方法将不同模式预报结果和观测结果进行集成后，对峰值浓度有较好的调整，12月21日06：00—22日06：00，观测的最大值为493 μg/m3，CMAQ、CAMx和NAQPMS预报的对应值分别为134、160、169 μg/m3，3个模式的平均值相对观测值低了339 μg/m3，多元线性回归方法集合后的值为420 μg/m3，接近最大值。对比定陵与万寿西宫的预报结果发现，各预报方法对定陵的清除时间把握较好，而对万寿西宫的预报则较实况提前3～4 h，主要是由于模式预报冷空气到达市区的估计较实况偏早。利用多元线性回归方法改进后，国控平均和万寿西宫的清除时间较实况有所滞后，这也说明，集成的预报成员对污染过程的把握越好，改进后的效果就越佳，因此，提高多模式系统成员的预报能力是提高重污染预报准确率的关键。

当预报值处于观测值的0.5～2倍范围内则为合理[17]，散点图中采用FAC表示预报值落在观测值的0.5～2倍范围内的比例，其值越大表明预报值落在观测值的0.5～2倍范围内的点越多，预报效果越好。由图6～图8可见，散点图的横纵坐标分别为观测和预报的PM2.5红色预警期间的小时浓度，当散点落在k=0.5下方时表明模式对实况低估了，而当散点落在k=2的上方时表明模式对实况高估了。

对国控平均而言(图6)，不同预报方法在红色预警期间观测与预报PM2.5的FAC处于73%～79%之间，对高浓度PM2.5预报有低估现象，这主要与模式系统中排放源的分布有关；而在PM2.5浓度较低时容易高估，这主要是由于红色预警期间，受有利气象条件的扰动，PM2.5浓度短时处于优良水平，而模式对这种弱气象系统反应不明显使预报的PM2.5浓度处于较高浓度水平。多元线性回归方法订正后FAC达到82%，对高浓度PM2.5预报的低估现象有一定程度的缓解。图7显示定陵不同预报方法的FAC处于49%～65%之间，处于k=0.5线下的散点较多，低估现象较国控平均和城区点位明显，而多元线性回归方法改进后，FAC可达到93%，低估现象基本得到缓解。对万寿西宫改进效果进行评估时(图8)发现，不同预报方法FAC处于72%～80%，改进后的FAC较CMAQ(72%)和NAPQMS(75%)高，达到77%，但低于CAMx(80%)。这主要是由于多元线性回归方法的改进效果直接受单模式预报效果的影响，模式系统对万寿西宫的预报效果为CAMx好于CMAQ和NAPQMS，从而导致集成效果略低于CAMx。