环境科学研究  2018, Vol. 31 Issue (1): 25-33  DOI: 10.13198/j.issn.1001-6929.2017.03.47

引用本文  

吴钲, 谢旻, 高阳华, 等. 利用集合均方根卡尔曼滤波反演重庆地区SO2源排放[J]. 环境科学研究, 2018, 31(1): 25-33.
WU Zheng, XIE Min, GAO Yanghua, et al. Inversion of SO2 Emissions over Chongqing with Ensemble Square Root Kalman Filter[J]. Research of Environmental Sciences, 2018, 31(1): 25-33.

基金项目

国家科技支撑计划项目(2014BAC16B06);重庆市气象局开放式研究基金项目(KFJJ-201607)
Supported by National Key Technology Research and Development Program of China (No.2014BAC16B06); Open Research Fund of Chongqing Meteorological Service, China (No.KFJJ-201607)

责任作者

谢旻(1978-), 男, 湖北黄冈人, 副教授, 博士, 主要从事城市区域空气质量模拟研究, minxie@nju.edu.cn

作者简介

吴钲(1985-), 男, 江苏兴化人, 工程师, 博士, 主要从事空气质量预报及相关研究, wukgdqghg@163.com

文章历史

收稿日期:2017-05-16
修订日期:2017-09-19
利用集合均方根卡尔曼滤波反演重庆地区SO2源排放
吴钲1 , 谢旻2 , 高阳华1 , 芦华1 , 赵磊1 , 高松1     
1. 重庆市气象科学研究所, 重庆 401147;
2. 南京大学, 江苏 南京 210093
摘要:传统的"自下而上"清单方法估算的排放清单,其数据的准确性和时效性存在较大局限.基于集合均方根卡尔曼滤波的源清单反演方法,结合WRF-CMAQ(天气研究和预报模式-公共多尺度空气质量模型)被用于对以清华大学编制的2010年MEIC(中国多尺度排放清单模型)排放清单为基础制作的重庆地区SO2排放源进行反演试验以解决准确性和时效性问题,试验时间段为2014年10月15-31日,重庆主城17个环境空气质量国控监测点ρ(SO2)小时观测资料用于反演及检验.结果表明:该方法能够反演重庆地区SO2源排放量,随着反演次数增加,基于反演排放源预报的ρ(SO2)预报误差持续减小,反演4次后预报误差达到比较低的稳定的水平,其均方根误差均低于20 μg/m3. 5次反演后SO2源排放量用于2014年10月24-29日每天起始预报,其预报的站点、时间平均的均方根误差从100~400 μg/m3降至30 μg/m3以下.反演中应用局地化尺度减少集合取样误差影响,54与81 km两个局地化尺度反演结果对预报改善效果相当,表明主要影响重庆主城ρ(SO2)的源排放位于主城及周边地区,也说明内源排放对重庆主城ρ(SO2)起主要影响.反演后面源排放量主城区降幅约为30 kg/(d·km2),周边地区减少10~20 kg/(d·km2),主城区部分SO2点源排放量降幅约为25 kg/(d·km2),说明2010年MEIC排放清单高估了试验时段重庆地区的SO2排放.
关键词SO2    WRF-CMAQ    源排放反演    集合均方根卡尔曼滤波    
Inversion of SO2 Emissions over Chongqing with Ensemble Square Root Kalman Filter
WU Zheng1 , XIE Min2 , GAO Yanghua1 , LU Hua1 , ZHAO Lei1 , GAO Song1     
1. Chongqing Institute of Meteorological Sciences, Chongqing 401147, China;
2. Nanjing University, Nanjing 210093, China
Abstract: Emissions inventories based on the 'bottom-up' inventory method have many problems with respect to accuracy and time. To solve this problem, an inversion of the SO2 emission over Chongqing from October 15th to 31st, 2014 was performed with WRF-CMAQ coupled with an ensemble square root Kalman filter. The initial SO2 emissions estimates come from the 2010 Multi-resolution Emission Inventory for China developed by Tsinghua University. Hourly observations of SO2 concentrations at 17 state air quality monitoring stations over the downtown area of Chongqing were used for inversion and verification. The results show that the ensemble square root Kalman filter is suitable for inverting SO2 emissions over Chongqing. The forecast error decreased gradually as the number of inversion emissions increased; additionally, the root mean square errors of the forecast concentration of SO2 were all less than 20 μg/m3 after four inversions. Localization was employed to account for ensemble sampling error. The experiments using a localization distance of 54 and 81 km had similar results, demonstrating that the concentrations of SO2 in Chongqing are primarily affected by the SO2 emission inventories surrounding the area and that the inner sources are primarily responsible for the concentration of SO2 in the downtown area of Chongqing. The average root mean square error of several station forecasts, initialized at each day from October 24th to 29th, 2014, decreased from 100-400 μg/m3 to concentrations less than 30 μg/m3. The estimated SO2 emissions of the area decreased by approximately 30 kg/(d·km2) and 10-20 kg/(d·km2) in the downtown and surrounding areas of Chongqing, respectively; meanwhile, point emissions decreased by approximately 25 kg/(d·km2) in part of the downtown area. The regional SO2 emissions inventories from 2010 MEIC were overestimated over Chongqing during the experiment.
Keywords: SO2    WRF-CMAQ    inversion estimation    ensemble square root Kalman filter    

重庆市地处四川盆地,是我国空气污染严重的四大区域之一,也是我国重要的酸雨和SO2控制区[1].在20世纪八九十年代,重庆地区能源结构以及煤为主,煤的含硫量高,颗粒物和SO2的污染十分严重,1991年后,特别是1997年重庆市直辖之后,重庆市加大了大气污染防治力度,开展了一系列计划和行动,调整了能源结构,降低了煤的比例,SO2污染得到了有效治理[2].李九彬等[3]分析了2001—2011年重庆市空气质量特征,虽然SO2污染天数呈逐年减少趋势,但依然是重要污染物之一. SO2是硫酸盐气溶胶的重要前体物,对大气灰霾污染和大气酸沉降有重要贡献,对霾的预报有重要的影响.

公共多尺度空气质量模型CMAQ[4]被广泛用于赛事管控效果分析评估[5-7]、污染物传输-沉降研究[8-9]和重污染过程分析[10],有较好的应用效果.重庆市气象局运行的环境气象数值预报系统采用了WRF-CMAQ模式,人为源清单为清华大学编制的2010年MEIC排放清单(http://www.meicmodel.org),系统评估检验表明,ρ(SO2)预报值异常偏高,约为观测值的10倍.在环境气象数值预报中,源排放数据会持续进入模式影响预报结果,其准确性对环境预报有重要影响[11].不论是从实际重庆市治理措施还是从模式结果评估看,对源清单数据中的SO2源排放进行处理均十分必要.

源清单信息的获取大致有两类方法:①传统的“自下而上”地通过实地调查结合污染排放模型得到,污染源的分布和变化十分复杂,编制清单需要较长时间以及人力来清查以及建立模型估算区域内各类污染物排放情况,由于工业快速发展或政府减少SO2排放政策的影响,这种估算方法得到的SO2排放会很快失去时效性[12]. ②在现有清单基础上,利用模式预报结合观测,通过某种方法对现有清单进行反演.如Lee等[12]利用地面观测和卫星观测结合GEOS-Chem模式的伴随模式反演了全球SO2排放,Elbern等[13]利用4DVAR反演欧洲区域多种污染物排放,Bousquet等[14]用Bayes变分分析方法反演全球CO2的源和汇,Kopacz等[15]比较了Bayes变分分析和伴随模式反演方法对亚洲地区CO排放的反演效果,蔡旭晖等[16]用逆向轨迹方法反演了北京地区甲烷排放源,程兴宏等[17]用“Nudging”源同化方法反演了华北地区重霾过程的SO2和NOx源排放.除此之外,集合卡尔曼滤波在源反演也有较多应用,如朱江等[18]对集合卡尔曼平滑和集合卡尔曼滤波在污染源反演中应用进行了研究,Miyazaki等[19]应用集合转换卡尔曼滤波反演2005—2006年全球NOx排放,TANG等[20]用集合卡尔曼滤波反演了北京及周边地区的CO排放源.

集合卡尔曼滤波有多种表现形式[21-23],该研究应用集合均方根卡尔曼滤波[23]对重庆地区SO2源进行反演,该方法不需要对观测进行扰动.应用集合卡尔曼滤波时需要考虑集合样本数不足引起的取样误差影响,“局地化”为减小该类误差的有效方法[24-25],最优局地化尺度与集合样本数相关[24].大量研究[19-20, 26-29]中应用集合卡尔曼滤波研究不同的对象时采用不同的局地化尺度,具体从几公里至几千公里不等,说明最优局地化尺度与研究对象有关.

该研究利用集合均方根卡尔曼滤波方法结合WRF-CMAQ建立区域SO2源排放反演系统,将重庆主城环境空气质量国控监测点的观测资料用于源反演,分析评估其对WRF-CMAQ在重庆区域预报的影响,以期提高排放源清单的准确性、改进业务系统ρ(SO2)的预报,同时为评估SO2排放源的影响、如何有效减排提供技术方法.

1 材料方法 1.1 模式和资料

该研究中气象模式为WRF V3.5.1,环境气象模式为CMAQ V4.7.1. WRF模式为美国国家大气研究中心(NCAR)、美国环境预测中心(NCEP)、美国国家海洋和大气管理局(NOAA)、美国俄克拉荷马大学等科研机构联合开发的中尺度数值天气预报系统(http://www.wrf-model.org/index.php),该模式使用完全可压缩、非静力平衡的欧拉原始方程组,具体细节参考模式说明[30]. CMAQ模式为是美国环境保护局(US EPA)开发的第3代区域空气质量模式(https://www.cmascenter.org/cmaq),该模式将大气作为一个整体,对环境大气中的物理、化学过程进行综合考虑. WRF模式采用Lambert投影,标准纬度为20°N和40°N,标准经度为106.5°E,双层嵌套网格,水平分辨率分别为27、9 km,格点数分别为200×160、288×216,中心点经度分别为104.5°E、106.5°E,中心点纬度分别为34.5°N、29.8°N,垂直方向层数为51层,物理参数化方案采用MYJ边界层方案、BMJ积云参数化方案、WSM6微物理过程. CMAQ模式同样采用双重嵌套,网格分辨率跟WRF模式一致,网格区域小于WRF模式,分别为160×130、198×166,其余区域设置与WRF一致.垂直方向,化学模式采用15个sigma层,其中边界层内分为8层. CMAQ模式为单向嵌套,先进行27 km区域计算,得到每小时浓度场作为9 km区域的初始和边界条件,再进行9 km区域的预报. CMAQ的水平输送和垂直输送采用ppm方案,水平扩散采用multi-scale方案,垂直扩散采用eddy方案,气相化学机制采用cb05cl机制,气溶胶为AERO4方案.

自然源排放由全球自然源排放模式MEGAN[31]计算.人为源排放采用2010年MEIC排放清单. WRF驱动场为0.5°×0.5°的GFS资料(3 h间隔).观测资料来源于重庆市环境保护局网站(http://www.cepb.gov.cn/xxgk/hjzl/index.shtml),研究使用了重庆主城17个环境空气质量国控监测点2014年10月15—31日的ρ(SO2)小时观测资料(见图 1).

注:底图数据来源来自国家基础地理信息中心, 1:4 000 000基础地理信息数据.下同. 图 1 重庆主城17个环境空气质量国控监测点分布 Fig.1 Distribution of 17 state air quality monitoring stations in the downtown area of Chongqing
1.2 反演方法

与集合数值天气预报扰动生成不同,源排放不需要满足一定的连续性条件以及平衡条件,该试验的集合扰动由蒙特卡洛方法[32]生成,对所有源排放加入随机扰动,扰动服从期望值为0、标准差为原清单值0.15倍的高斯分布,扰动后值低于0的设为0.试验仅考虑SO2排放对ρ(SO2)的影响,只对SO2排放量进行反演订正.为避免模式预报日变化的误差对反演的影响,应用一天24次的小时观测值进行反演.由于观测点比较密集(见图 1),为减小计算量以及代表性误差影响,基于最优估计理论对观测进行“超级观测”处理[19].采用临近点插值方法得到观测对应的预报值,根据最优估计理论将所有位于同一格点的m个观测值yi(i=1, 2, 3, …, m)和n个成员对应的预报值xik(i=1, 2, …, m; k=1, 2, …, n)构成超级观测值以及对应预报值,假设不同时次、不同站点的观测误差相互独立,yi对应的观测标准差为ri,则超级观测值ynew、观测标准差rnew和对应的预报值xnewn满足下列关系:

$ \frac{1}{{{r_\rm{new}}^2}} = \sum\limits_{i = 1}^m {\frac{1}{{{r_i}^2}}} $ (1)
$ {y_\rm{new}} = \left( {\sum\limits_{i = 1}^m {\frac{{{y_i}}}{{{r_i}^2}}} } \right)/\left( {\sum\limits_{i = 1}^m {\frac{1}{{{r_i}^2}}} } \right) $ (2)
$ {\rm{ }}{x_\rm{new}}^k = \left( {\sum\limits_{i = 1}^m {\frac{{{x_i}^k}}{{{r_i}^2}}} } \right)/\left( {\sum\limits_{i = 1}^m {\frac{1}{{{r_i}^2}}} } \right) $ (3)

反演应用超级观测.唐晓等[28, 33]的研究中均假定观测误差为观测值的10%,试验阶段ρ(SO2)观测值低于50 μg/m3,为了便于处理,ρ(SO2)观测标准差均设为5 μg/m3.该研究应用集合均方根卡尔曼滤波进行反演试验,其计算公式:

$ {\mathit{\boldsymbol{\overline x}} ^{\rm{a}}} = {\mathit{\boldsymbol{\overline x}} ^{\rm{b}}} + \mathit{\boldsymbol{K}}({\mathit{\boldsymbol{y}}^{\rm{o}}} - \mathit{\boldsymbol{H}}{\mathit{\boldsymbol{\overline x}} ^{\rm{b}}}){\rm{ }} $ (4)
$ \mathit{\boldsymbol{x}}{' ^{\rm{a}}} = \mathit{\boldsymbol{x}}{' ^{\rm{b}}} - \mathit{\boldsymbol{\widetilde K}}\mathit{\boldsymbol{Hx}}{' ^{\rm{b}}} $ (5)
$ \mathit{\boldsymbol{K}} = {\mathit{\boldsymbol{P}}^{\rm{b}}}{\mathit{\boldsymbol{H}}^{\rm{T}}}{(\mathit{\boldsymbol{H}}{\mathit{\boldsymbol{P}}^{\rm{b}}}{\mathit{\boldsymbol{H}}^{\rm{T}}} + \mathit{\boldsymbol{R}})^{ - 1}} $ (6)

顺序同化独立观测时

$ \mathit{\boldsymbol{\widetilde K}} = \left( {1 + {\sqrt {\frac{\mathit{\boldsymbol{R}}}{{\mathit{\boldsymbol{H}}{\mathit{\boldsymbol{P}}^{\rm{b}}}{\mathit{\boldsymbol{H}}^{\rm{T}}} + \mathit{\boldsymbol{R}}}}}}} \right) ^{ - 1}\mathit{\boldsymbol{K}} $ (7)

式中,xaxb为反演前后的平均源排放量,xaxb为对应的源排放量扰动,yoynew组成的矢量,R为观测误差协方差矩阵,H为观测算子,${\mathit{\boldsymbol{P}}^{\rm{b}}} = \frac{1}{{n - 1}}\mathit{\boldsymbol{x}}{' ^{\rm{b}}}{(\mathit{\boldsymbol{x}}{' ^{\rm{b}}})^{\rm{T}}}$为背景误差协方差矩阵,K$\mathit{\boldsymbol{\widetilde K}}$为集合平均和集合扰动的卡尔曼增益矩阵.滤波在观测空间进行,直接应用观测对应的预报值作为Hxb计算,试验中观测算子H表示WRF-CMAQ基于源排放量计算出ρ(SO2)格点预报值后再插值到观测站点的整个计算过程.在观测相互独立情况下同时同化所有观测与顺序同化的结果相同[29].该研究采用顺序同化方案[34],每次反演一个观测,此时R以及HPb为标量,矩阵求逆计算变为倒数计算,计算量和实现难度降低.

1.3 试验设计

唐晓等[28]研究臭氧资料同化时发现集合样本数超过20即可达到比较好的效果,Miyazaki等[19]研究全球NO2集合反演时发现32个集合成员比较合适,笔者使用31个集合样本进行试验.反演采用与业务相似的流程(见图 2).首先对源清单进行满足高斯分布的随机扰动,生成30组源清单并用于第1天模式预报,结合未扰动清单预报组成31个成员的集合预报,通过最优估计将25~48 h的小时预报以及观测得到超级观测和对应预报,在集合均方根卡尔曼滤波中同化超级观测得到反演后源清单,将平均源清单作为新的未扰动源清单应用于第2天的集合生成、预报、反演中,重复上述步骤就可以循环更新源清单.该研究的源清单反演更新时间从2014年10月15日开始,模式起报时刻均为08:00(北京时间,下同).主要反演时间段内以重庆沙坪坝站为中心的1°×1°区域内国家基本站10 m高度平均纬向和经向风速约为0和0.18 m/s,基本为静风态,气象模式WRF预报10 m纬向和经向风速约为-0.43和0.21 m/s,预报值略高于观测值,预报东风略偏强,总体反演期间气象模式预报与实际差异不大,可以用于源反演.

注:ENSRF为集合均方根卡尔曼滤波. 图 2 计算流程示意 Fig.2 Schematic diagram of calculation flow

远距离的源排放理论上对本地污染物浓度没有影响,但受取样误差、模式误差影响,二者可能存在虚假相关,因此需要引入局地化方案将反演区域限定在观测周边一定范围以内.气象场为连续场,同一时间相邻区域有一定的相关性,应用局地化时需考虑随距离增加而减少的相关性特征[29].由于大气传输扩散作用[8, 35-37],污染物浓度与源排放的相关性没有明显的随距离减少的特点,应用随距离减少的局地化方案可能会订正不足,同时源排放量可以不连续,应用阶跃函数不会产生计算不稳定等问题.综上考虑,试验中采用如下局地化函数:

$ C\left( d \right) = \left\{ \begin{array}{l} 2, d \leqslant L\\ 0, d > L \end{array} \right. $ (8)

式中:d为反演格点与观测点之间距离,km; L为局地化尺度,km.最优局地化尺度与大气传输能力和物种生命史相关,大气平流传输越强,源排放影响范围越大,最优局地化尺度越大; 物种生命史越长,扩散和传输作用时间越长,最优局地化尺度越大. SO2冬季生命史约为(58±20)h[12],具有区域污染特征[17],最优尺度应该较大,但反演时段重庆地区风速较低,大气传输作用相对较小,最优尺度应该较小.为研究局地化尺度的影响,设计两组试验,局地化尺度分别采用54和81 km,对应27 km网格下的2倍和3倍格距以及9 km网格下的6倍和9倍格距.该研究中将两组试验分别称为L6试验(L=54 km)和L9试验(L=81 km),初始源预报称为控制试验.为了便于描述,对于对比、检验、源排放变化等均只给出9 km的结果,27 km网格区域与9 km结果类似.

2 反演结果 2.1 反演一次结果

集合均方根卡尔曼滤波反演源排放本质上是通过预报值与源排放量的相关性实现.位于格点(100,79)的ρ(SO2)预报值(其中一个观测所在的格点)与周边SO2面源和点源排放量的相关系数分布如图 3所示.根据统计理论,自由度为30时满足0.05和0.01显著水平的相关系数分别需要达到0.349和0.449.以0.01显著水平为标准,同一格点的面源排放量与ρ(SO2)预报值达到显著相关,其相关系数超过0.8〔见图 3(a)〕,其他格点的面源、点源排放量与该格点ρ(SO2)预报值的相关系数较低,说明同一格点的面源排放量对该格点ρ(SO2)预报值有主要影响.位于格点(100,79)周边部分格点面源和点源排放量与ρ(SO2)预报值的相关达到了0.05的显著水平,说明可能存在一定的短距离平流传输作用〔见图 3(b)〕.部分格点相关系数为负值,表示源排放量增加(减少)时ρ(SO2)减小(增加),与理论不符,这与取样误差有关.

注:黑色圆点为格点(100, 79). 图 3 2014年10月15日08:00起报的格点(100,79)SO2预报值与格点(90~110,69~89)SO2排放量的相关系数分布 Fig.3 Correlations between SO2 emissions in grid points (90-110, 69-89) and the forecast ρ(SO2) in grid point (100, 79) initialized at 08:00 October 15th, 2014

集合平均(未扰动)的清单被认为相对最合理,因此未扰动成员的预报被用于检验反演效果.为避免初始浓度误差以及气象模式起转磨合(spin-up)影响,只有模式16~64 h的ρ(SO2)预报结果被用于模式检验. 图 4为一次反演订正后L6、L9和控制试验2014年10月16日08:00起报站点平均的ρ(SO2)预报对比结果,平均误差从265.2 μg/m3降至145.3 μg/m3(L6试验)和144.7 μg/m3(L9试验),均方根误差从333.2 μg/m3降至186.7 μg/m3(L6试验)和185.8 μg/m3(L9试验),两种误差均为原误差的55%,误差明显降低但依然较大,其中16个站点预报误差减小,但鱼洞站预报误差增加(图略).

注:模式起报时间为10月16日08:00. 图 4 反演一次后各试验组站点平均的ρ(SO2)预报值和观测值时间序列 Fig.4 Average temporal variation of SO2 concentration forecast of each experiment and observation after one inversion

SO2源清单变化情况如图 5所示,反演后观测所在区域排放量绝大部分减少,面源减少范围相比点源更广.部分站点反演订正后排放量有所增加,这应该是鱼洞站预报变差原因,本质是取样误差引起.两种试验源清单变化在重庆主城区(观测站所在区域)差异不大,在局地化半径影响下离主城较远的区域反演结果有一定差别.

图 5 L6和L9试验反演一次后SO2面源和点源排放量变化情况 Fig.5 Innovations of area and point SO2 emissions of experiment L6 and L9 after one inversion
2.2 反演多次结果

图 6为源排放循环反演7次过程中站点和预报时间平均SO2平均误差和均方根误差变化结果.随着反演次数增加,预报误差持续减少,L6和L9试验的结果差异很小.反演一次后均方根误差约为原误差的55%,反演两次后约为原误差29%,反演4次后预报误差达到比较低的稳定的水平,均方根误差约为原误差的5%左右,从平均300 μg/m3减至低于20 μg/m3,平均误差约为原误差的3%,从平均280 μg/m3减至低于10 μg/m3.集合均方根卡尔曼滤波根据预报和观测的差来反演源排放,反演3次后预报值与观测没有较大的偏差,同时清单的不确定性在反演后有所减小,因此清单变化较小.需要说明的是,两次反演后所有站点的误差均降低(图略).

图 6 各试验组站点和预报时间平均的ρ(SO2)平均误差和均方根误差随反演次数变化情况 Fig.6 Variation in inversion times of the average mean errors and root mean square errors of the forecast time and stations of all experiments

图 7为5次反演后清单排放量相比初始清单变化结果,最大变化区域位于重庆主城,SO2面源排放量降幅约为30 kg/(d·km2),主城周边面源排放量减少10~20 kg/(d·km2),点源变化区域和量级小于面源变化,主城区部分点源排放量降幅约25 kg/(d·km2),说明原清单高估了试验时段重庆主城区SO2排放量.重庆主城附近两组试验清单变化相似,不同之处在于L9试验中在离主城较远的区域清单变化明显,特别是东部区域.反演5次后从预报误差看,L9试验略优于L6试验,但考虑到观测误差等影响,可以认为两组试验预报误差没有显著区别,L9试验比L6试验多反演的源对预报没有显著影响,54 km(L6)的局地化尺度对这一段时间的重庆主城区SO2反演更合适.

图 7 L6和L9试验反演5次后SO2面源和点源排放量相比原清单变化情况 Fig.7 Innovations of area and point SO2 emissions in experiment L6 and L9 after five inversions compared with the initial SO2 emissions estimates

将反演5次的SO2源排放用于10月24—29日每日起始的预报,其预报误差相比采用原始源排放的预报有显著减少,两种局地化尺度反演预报的误差接近,不同日期预报的均方根误差差异较小,26—27日控制试验预报平均误差较小,原因是模式预报低层南风异常偏强,原MEIC排放清单的SO2排放量在重庆南部区域排放低,南风增强有利于污染物排出,因而ρ(SO2)偏低与观测接近,误差较其他日期小(见图 8).除26日外,其他日期基于反演5次源预报的平均误差和均方根误差均低于原误差的8%,均方根误差从100~400 μg/m3降至30 μg/m3以下,平均误差绝对值从80~300 μg/m3降至8 μg/m3以下,这与源排放改进后模式对排放源的输送扩散等过程的影响有关[38].

图 8 2014年10月24—29日每日起报的采用初始源和L6、L9试验反演5次源的站点平均和预报时间平均的ρ(SO2)均方根误差 Fig.8 Average root mean square error of SO2 concentrations initialized each day from October 24th to 29 th, 2014 with the initial emissions and emissions after five inversions of experiments L6 and L9
3 讨论

TANG等[20]应用了集合卡尔曼滤波对北京CO排放进行反演,站点检验说明反演后预报结果优于反演前.该研究中基于集合均方根卡尔曼滤波的与系统业务运行相似的循环更新方案被设计用于反演重庆地区SO2源清单排放,随着反演订正次数增加,均方根误差与原误差的比例减少为约55%、29%,直至达到约5%,误差值低于20 μg/m3,反演5次后源清单应用于10月24—29日预报,其误差也有明显降低.影响大气成分的主要过程包括排放、输送扩散、化学变化以及沉降这些物理化学过程.已有模拟经验显示,排放源的不确定性对最终模拟结果的影响很大[11].经过多次源反演后,ρ(SO2)模拟误差显著减小的主要原因是集合均方根卡尔曼滤波对源的反演改进了排放,进而改进了ρ(SO2)预报结果.该结果同时说明CMAQ模式对SO2输送扩散、化学变化和沉降过程的描述比较理想.反演后SO2源排放量相比原始清单减少,这与政府治理调控、能源结构调整的结果[2-3]相一致,说明清华大学MEIC 2010年的清单高估了反演试验阶段的SO2源排放.

集合样本数限制引起取样误差是集合卡尔曼滤波应用存在的主要问题之一[28, 32],取样误差表现为远距离虚假相关,如试验中部分源与观测的相关性为负值(见表 1),表示源排放量增加(减少)时,预报浓度减小(增加),与实际情况不符.局地化是解决远距离虚假相关的有效方法[20, 24-25, 29],如Whitaker等[29]在NCEP全球预报系统中应用集合卡尔曼滤波同化资料,试验中局地化尺度为2 800 km;唐晓等[28]在区域臭氧同化试验中局地化尺度采用54 km;TANG等[20]在北京地区CO源反演中分别采用54和45 km用于同化和反演;Miyazaki等[19]研究表明,450 km是全球1月NO2源排放分析最优局地化尺度.笔者采用54和81 km局地化尺度分别进行试验,检验结果表明两组试验反演效果相似,在54 km局地化尺度影响范围内源清单变化相似,说明观测主要影响区域在54 km范围内,本地源排放对ρ(SO2)影响明显,在重庆地区本地源排放控制对本地SO2污染治理至关重要,这与李荔等[5]南京青奥会期间管控措施评估得出的内源控制重要性结论一致.反演时段内主要区域基本为静风,模拟风速略高于观测,存在一定的虚假平流输送,对离观测较远区域源的反演存在较大的不确定性,两组试验订正效果相似,在东部地区差异较大,预报东风比实际偏强,存在虚假平流,说明L9试验中东部地区可能存在虚假反演,54 km局地化尺度相对更合理.实况风速弱于模式预报,污染物更容易堆积,本地源影响更强,气象预报的误差不影响初始源清单高估该区域SO2排放量的结果.

受限于集合方案设计,该方法可以有效改善排放量较大的源,很难反演出虚假的源以及初始源没有但实际存在的排放.该试验中,重庆主城区及周边初始源SO2排放量较高,观测与预报误差较大,该方法效果明显.源反演的观测算子隐含在模式预报中,模式预报的性能即模式对气象场的描述能力、对源排放物种输送、扩散、沉降、相互作用等各种过程的描述能力对源反演有决定性影响,该试验没有充分考虑气象模式和和环境模式的模拟预报能力对反演效果的影响,需要进一步深入的研究.

4 结论

a) 基于集合均方根卡尔曼滤波设计的反演方法能够反演重庆地区SO2源排放量,试验结果表明,随着反演次数增加,预报误差持续减少,L6和L9试验的结果差异很小,反演4次后预报误差达到比较低的稳定的水平,均方根误差约为原误差的5%左右,从平均300 μg/m3减至低于20 μg/m3,平均误差约为原误差的3%,从平均280 μg/m3减至低于10 μg/m3.基于反演5次的源排放量预报均方根误差与原误差比例基本低于8%,均方根误差从100~400 μg/m3降至30 μg/m3以下,平均误差绝对值从80~300 μg/m3降至8 μg/m3以下,说明该方法能较好地反演源排放.

b) 反演后重庆地区SO2源排放量相比原清单明显减小,反演5次后重庆主城区SO2面源排放量降幅约为30 kg/(d·km2),周边地区面源排放量减少10~20 kg/(d·km2),主城区部分SO2点源排放量降幅约为25 kg/(d·km2),说明2010年MEIC排放清单高估了重庆地区试验阶段SO2排放量.

c) 局地化尺度试验表明,54和81 km局地化尺度反演的结果对预报改善效果类似,说明对重庆主城主要影响源排放区域主城及周边区域,54 km局地化尺度相对更适合,本地及近周边源排放量对ρ(SO2)有重要作用,SO2污染治理的关键是内源排放控制.

参考文献
[1]
徐鹏, 郝庆菊, 吉东生, 等. 重庆市北碚大气中PM2.5、NOx、SO2和O3浓度变化特征研究[J]. 环境科学学报, 2016, 36(5): 1539-1547.
XU Peng, HAO Qingju, JI Dongsheng, et al. Characteristics of atmospheric PM2.5, NOx, SO2 and O3 in Beibei District, Chongqing[J]. Acta Scientiae Circumstantiae, 2016, 36(5): 1539-1547. (0)
[2]
江文华, 刘德, 陈勇航, 等. 1980—2012年重庆地区霾日时空变化特征[J]. 干旱气象, 2015, 33(4): 602-606.
JIANG Wenhua, LIU De, CHEN Yonghang, et al. Temporal and spatial variations of haze days in Chongqing from 1980 to 2012[J]. Journal of Arid Meteorology, 2015, 33(4): 602-606. (0)
[3]
李九彬, 王建力. 2001—2011年重庆市空气质量特征分析[J]. 西南大学学报(自然科学版), 2013, 35(9): 145-153.
LI Jiubin, WANG Jianli. An analysis of air quality characteristics of Chongqing(2001-2011)[J]. Journal of Southwest University (Natural Science Edition), 2013, 35(9): 145-153. (0)
[4]
BINKOWSKI F S, ROSELLE S J. Models-3 community multiscale air quality (CMAQ) model aerosol component 1.model description[J]. Journal of Geophysical Research Atmospheres, 2003, 108(D6): 335-346. (0)
[5]
李荔, 刘倩, 李冰, 等. 南京青奥会期间管控措施空气质量改善效果评估[J]. 环境科学研究, 2016, 29(2): 175-182.
LI Li, LIU Qian, LI Bing, et al. Assessment of air quality benefits from control measures during Nanjing Youth Olympic Games[J]. Research of Environmental Sciences, 2016, 29(2): 175-182. (0)
[6]
陈焕盛, 吴其重, 王自发, 等. 华北火电厂脱硫对奥运期间区域空气质量的影响[J]. 环境科学学报, 2014, 34(3): 598-605.
CHEN Huansheng, WU Qizhong, WANG Zifa, et al. Impacts of power plant desulfurization on regional air quality over North China Plain during the Olympic Games[J]. Acta Scientiae Circumstantiae, 2014, 34(3): 598-605. (0)
[7]
李岩, 安兴琴, 左洪超, 等. 2008年北京奥运交通限制效果的模式研究[J]. 高原气象, 2010, 29(6): 1619-1626.
LI Yan, AN Xingqin, ZUO Hongchao, et al. Study of traffic restriction effects on air quality during Beijing Olympic Games in 2008 using CMAQ model[J]. Plateau Meteorology, 2010, 29(6): 1619-1626. (0)
[8]
张珏, 孟凡, 何友江, 等. 长江三角洲地区大气二英类污染物输送-沉降模拟研究[J]. 环境科学研究, 2011, 24(12): 1393-1402.
ZHANG Jue, MENG Fan, HE Youjiang, et al. Simulation of transport-deposition of dibenzo-p-dioxins and dibenzofurans in the atmosphere in the Yangtze River Delta region[J]. Research of Environmental Sciences, 2011, 24(12): 1393-1402. (0)
[9]
边敏娟, 钱骏, 刘志红, 等. 四川中东部地区2009年大气硫沉降模拟[J]. 环境科学研究, 2012, 25(10): 1115-1119.
BIAN Minjuan, QIAN Jun, LIU Zhihong, et al. Modeling study on atmospheric sulfur deposition in the eastern and central regions of Sichuan area in 2009[J]. Research of Environmental Sciences, 2012, 25(10): 1115-1119. (0)
[10]
余纬, 罗栩羽, 范绍佳, 等. 珠三角一次重空气污染过程特征分析及数值模拟[J]. 环境科学研究, 2011, 24(6): 645-653.
YU Wei, LUO Xuyu, FAN Shaojia, et al. Characteristics analysis and numerical simulation study of a severe air pollution episode over the Pearl River Delta[J]. Research of Environmental Sciences, 2011, 24(6): 645-653. (0)
[11]
唐孝炎, 张远航, 邵敏. 大气环境化学[M]. 2版.北京: 高等教育出版社, 2006, 447-453. (0)
[12]
LEE C, MARTIN R V, VAN DONKELAAR A, et al. SO2 emissions and lifetimes:estimates from inverse modeling using in situ and global, space-based (SCIAMACHY and OMI) observations[J]. Journal of Geophysical Research, 2011. DOI:10.1029/2010JD014758 (0)
[13]
ELBERN H, STRUNK A, SCHMIDT H, et al. Emission rate and chemical state estimation by 4-dimensional variational inversion[J]. Atmospheric Chemistry and Physics, 2007, 7(1): 3749-3769. (0)
[14]
BOUSQUET P, CIAIS P, PEYLIN P, et al. Inverse modeling of annual atmospheric CO2 sources and sinks 1.method and control inversion[J]. Journal of Geophysical Research Atmospheres, 1999, 104(D21): 26161-26178. DOI:10.1029/1999JD900342 (0)
[15]
KOPACZ M, JACOB D J, HENZE D K, et al. Comparison of adjoint and analytical Bayesian inversion methods for constraining Asian sources of carbon monoxide using satellite (MOPITT) measurements of CO columns[J]. Journal of Geophysical Research Atmospheres, 2009. DOI:10.1029/2007JD009264 (0)
[16]
蔡旭晖, 邵敏, 苏芳. 甲烷排放源逆向轨迹反演模式研究[J]. 环境科学, 2002, 23(5): 19-24.
CAI Xuhui, SHAO Min, SU Fang. A backward trajectory inversion model for methane emission over Beijing area[J]. Environmental Science, 2002, 23(5): 19-24. (0)
[17]
程兴宏, 徐祥德, 安兴琴, 等. 2013年1月华北地区重霾污染过程SO2和NOx的CMAQ源同化模拟研究[J]. 环境科学学报, 2016, 36(2): 638-648.
CHENG Xinghong, XU Xiangde, AN Xingqin, et al. Inverse modeling of SO2 and NOx emissions using an adaptive nudging scheme implemented in CMAQ model in North China during heavy haze episodes in January 2013[J]. Acta Scientiae Circumstantiae, 2016, 36(2): 638-648. (0)
[18]
朱江, 汪萍. 集合卡尔曼平滑和集合卡尔曼滤波在污染源反演中的应用[J]. 大气科学, 2006, 30(5): 871-882.
ZHU Jiang, WANG Ping. Ensemble Kalman smoother and ensemble Kalman filter approaches to the joint air quality state and emission estimation problem[J]. Chinese Journal of Atmospheric Sciences, 2006, 30(5): 871-882. (0)
[19]
MIYAZAKI K, ESKES H J, SUDO K. Global NOx emission estimates derived from an assimilation of OMI tropospheric NO2 columns[J]. Atmospheric Chemistry & Physics, 2012, 11(5): 31523-31583. (0)
[20]
TANG X, ZHU J, WANG Z F, et al. Inversion of CO emissions over Beijing and its surrounding areas with ensemble Kalman filter[J]. Atmospheric Environment, 2013, 81(4): 676-686. (0)
[21]
BISHOP C H, ETHERTON B J, MAJUMDAR S J. Adaptive sampling with the ensemble transform Kalman filter.part I:theoretical aspects[J]. Monthly Weather Review, 2001, 129(3): 420-436. DOI:10.1175/1520-0493(2001)129<0420:ASWTET>2.0.CO;2 (0)
[22]
ANDERSON J L. An ensemble adjustment Kalman filter for data assimilation[J]. Monthly Weather Review, 2003, 129(12): 2884-2903. (0)
[23]
WHITAKER J S, HAMILL T M. Ensemble data assimilation without perturbed observations[J]. Monthly Weather Review, 2002, 130(7): 1913-1924. DOI:10.1175/1520-0493(2002)130<1913:EDAWPO>2.0.CO;2 (0)
[24]
HAMILL T M, WHITAKER J S, SNYDER C. Distance-dependent filtering of background error covariance estimates in an ensemble Kalman filter[J]. Monthly Weather Review, 2001, 129(11): 2776-2790. DOI:10.1175/1520-0493(2001)129<2776:DDFOBE>2.0.CO;2 (0)
[25]
ANDERSON J L. Exploring the need for localization in ensemble data assimilation using a hierarchical ensemble filter[J]. Physica D:Nonlinear Phenomena, 2007, 230(1/2): 99-111. (0)
[26]
SEKIYAMA T T, TANAKA T Y, SHIMIZU A, et al. Data assimilation of CALIPSO aerosol observations[J]. Atmospheric Chemistry and Physics, 2009, 10(1): 39-49. (0)
[27]
SOBASH R A, STENSRUD D J. The impact of covariance localization for radar data on EnKF analyses of a developing MCS:observing system simulation experiments[J]. Monthly Weather Review, 2013, 141(11): 3691-3709. DOI:10.1175/MWR-D-12-00203.1 (0)
[28]
唐晓, 朱江, 王自发, 等. 基于集合卡尔曼滤波的区域臭氧资料同化试验[J]. 环境科学学报, 2013, 33(3): 796-805.
TANG Xiao, ZHU Jiang, WANG Zifa, et al. Regional ozone data assimilation experiment based on ensemble Kalman filter[J]. Acta Scientiae Circumstantiae, 2013, 33(3): 796-805. (0)
[29]
WHITAKER J S, HAMILL T M, WEI X, et al. Ensemble data assimilation with the NCEP global forecast system[J]. Monthly Weather Review, 2008, 136(2): 463-482. DOI:10.1175/2007MWR2018.1 (0)
[30]
SKAMAROCK W C, KLEMP J B, DUDHIA J, et al.A description of the advanced research WRF version 3[EB/OL].Boulder, Colorado, USA:National Center for Atmospheric Research, 2008[2008-06].http://www2.mmm.ucar.edu/wrf/users/docs/arw_v3.pdf. (0)
[31]
GUENTHER A, KARL T, HARLEY P, et al. Estimates of global terrestrial isoprene emissions using MEGAN (Model of Emissions of Gases and Aerosols from Nature)[J]. Atmospheric Chemistry & Physics, 2006, 6(11): 3181-3210. (0)
[32]
EVENSEN G. Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics[J]. Journal of Geophysical Research Oceans, 1994, 99(C5): 10143-10162. DOI:10.1029/94JC00572 (0)
[33]
崔应杰, 王自发, 朱江, 等. 空气质量数值模式预报中资料同化的初步研究[J]. 气候与环境研究, 2006, 11(5): 616-626.
CUI Yingjie, WANG Zifa, ZHU Jiang, et al. A preliminary study on data assimilation for numerical air quality model prediction[J]. Climatic and Environmental Research, 2006, 11(5): 616-626. (0)
[34]
HOUTEKAMER P L, MITCHELL H L. A sequential ensemble Kalman filter for atmospheric data assimilation[J]. Monthly Weather Review, 2001, 129(1): 123-137. DOI:10.1175/1520-0493(2001)129<0123:ASEKFF>2.0.CO;2 (0)
[35]
张玉洽, 杨迎春, 李杰, 等. 东南亚生物质燃烧对我国春季PM2.5质量浓度影响的数值模拟[J]. 环境科学研究, 2016, 29(7): 952-962.
ZHANG Yuqia, YANG Yingchun, LI Jie, et al. Modeling the impacts of biomass burning in Southeast Asia on PM2.5 over China in spring[J]. Research of Environmental Sciences, 2016, 29(7): 952-962. (0)
[36]
王芳, 陈东升, 程水源, 等. 基于气流轨迹聚类的大气污染输送影响[J]. 环境科学研究, 2009, 22(6): 637-642.
WANG Fang, CHEN Dongsheng, CHENG Shuiyuan, et al. Impacts of air pollutant transport based on air trajectory clustering[J]. Research of Environmental Sciences, 2009, 22(6): 637-642. (0)
[37]
程念亮, 李云婷, 张大伟, 等. 2014年10月北京市4次典型空气重污染过程成因分析[J]. 环境科学研究, 2015, 28(2): 163-170.
CHENG Nianliang, LI Yunting, ZHANG Dawei, et al. Analysis about the characteristics and formation mechanisms of serious pollution events in October 2014 in Beijing[J]. Research of Environmental Sciences, 2015, 28(2): 163-170. (0)
[38]
程兴宏. 空气质量模式"源同化"模型及排放源影响效应研究[D]. 北京: 中国科学院研究生院, 2008: 1-149. (0)