2017, Vol. 34 
随着经济的发展和城市化的快速推进，污染物排放水平迅速增加，中国大多数城市都面临着严重的复合型大气污染问题。空气污染已成为主要的经济和社会问题[1]，影响着人们的健康[2-3]。在中国环境监测总站发布的《2013年上半年74城市空气质量状况报告》中，仅4个城市PM2.5浓度达到或优于国家二级标准，达标率为5.4%。PM2.5平均日均值超标率为35.4%，成为影响空气质量的首要污染物。
近年来，越来越多的研究者开始关注以大气科学领域为基础建立的空气质量模型。其中，以“一个大气”为核心思想的第三代模型Models-3/CMAQ[4]是现有的一个既先进又成熟的空气质量模型，已经在国内外取得广泛的应用。Wang等[5]采用CMAQ研究2013年1月河北南部暴发的重度雾霾事件，发现需要建立更大范围的联防联控政策才能有效地缓解当地的空气污染。Zhang等[6]预测北卡罗纳州2018年的PM2.5和O3浓度，发现政府的减排措施能显著降低PM2.5的浓度水平。
但是，由于空气质量模拟系统是由气象模型、源清单处理模型和大气化学模型等构成的复杂系统，需要准备大量的输入文件，在实际应用中运算时间长，并且由于模型只能在Linux平台上运行，对使用人员要求较高。为解决这些问题，应用在空气质量领域的响应曲面模型RSM (response surface model)[7]应运而生。严茹莎等[8]利用RSM研究长三角地区O3的敏感性特征，并提出上海市O3污染的最优控制方案。劳苑雯等[9]基于所建立的RSM，对美国8个典型城市的PM2.5污染状况进行控制情景分析。本文基于CMAQ模拟得到对应的污染物浓度，统计、归纳出各种污染物排放控制因子与污染物浓度之间的函数关系，建立空气质量响应曲面模型RSM，分析不同污染物的PM2.5响应值以及评估杭州市2020年的政府减排方案的控制效果。
1 响应曲面模型 1.1 模型简介RSM是一个基于数学统计理论的纯解析模型，通过一定数量的多维实验，借助数学统计手段，便可归纳并建立起用以表征某一响应量与一系列独立控制因素之间的响应关系。RSM的概念表达式[10]如下所示：
| $ \eta = f\left( {{x_1},{x_2}, \cdots ,{x_n}} \right) + \varepsilon , $ |
式中：η表示响应变量；f表示响应关系；x1, x2, …, xn等表示独立的控制因素；ε表示统计误差。
一般地，建立一个RSM主要分3步：首先，确定独立变量及其变化水平；其次，选择实验设计方案，预测响应变量的模型方程；最后，求解该模型方程并检验模型的可靠性。3步之中的每一步都很重要，独立变量的选取决定了模型的输入，实验方案设计的统计样本是否具有代表性、数量是否合适、分布是否合理很大程度上决定了最后归纳得到的模型方程的适用范围与准确性，模型的可靠性检验说明最终得到模型的可信度。尽管RSM的建立过程是复杂繁琐的，但一旦模型建立以后，RSM便是一种高度经济、简洁的响应变量预测模型，既可通过输入独立变量的值来预测响应值，也可评价多个输入变量的影响及其交互作用。
1.2 空气质量RSM针对大气污染实际决策应用中需建立高效、快速的“污染物排放-空气质量”响应模型这一问题，USEPA借鉴蒙特卡洛 (Monte Carlo) 大规模采样实验的原理[11-12]，直接通过CMAQ的上百次模拟结果，采用多维克里金插值的非线性拟合方法，归纳得到O3和PM对10~20种分部门分区域排放源的RSM作为控制决策的支持工具。自此之后，RSM也受到了空气质量模拟领域的研究者们的广泛关注。在RSM成功建模后，它本质上就是一个可用于预测空气质量的“简化模型”。通过输入污染物排放控制因子组合，便可快速得到该因子组合对应的控制情景下的污染物浓度。RSM的建立方法很直接，并不需要涉及空气质量模型内部的复杂机制，因此它适用于任何一种空气质量模型，可以对任何污染物和排放源的响应情况进行分析。
2 空气质量响应曲面模型的构建 2.1 控制因子选取对于复杂的非线性响应问题，一般来说都会有大量的参数对其产生影响。由于不可能将所有的影响因素都考虑在内，选取其中具有明显作用的控制因子就显得十分重要。控制因子的数量决定了采样空间的维度和需要的样本数，是RSM建模的关键，直接关系到建模的成败。
控制因子的选取取决于污染物前体排放的类型和污染来源，应从决策目的出发，将其中涉及的来自不同地区、不同部门的相关污染物作为控制因子。本文选取杭州市作为研究对象，研究不同部门的相关污染物排放对杭州市PM2.5污染的影响。因此，控制因子的选取应基于PM2.5污染前体物种类以及排放源类别，选取来自不同部门的各种污染物中相对效能较大的因子。由于PM2.5来源中一次微粒与污染物浓度的响应关系基本为线性，而二次气溶胶组分的非线性特征最为明显。依据吴小芳[13]的研究，本文选取OC、BC和PM10等3个作为线性控制因子，电厂源NOX、工业源NOX、交通源NOX、电厂源SO2、交通源SO2、NH3和VOC等7个作为非线性控制因子。
对于参数水平的设置，约定参数水平等于目标年排放量与基准年 (2010年) 排放量的比值。参数水平为1，表示目标年与基准年排放量相同，不做控制；0为最小值，代表没有可控源污染物排放；可控源污染物如果在目标年排放量增大，则其参数水平将大于1。表 1所示即为选取的控制因子及相应的参数水平汇总。
|
表 1 响应曲面模型的控制因子 Table 1 Control factors in the RSM model |
RSM的建立是在大量实验结果基础上的统计归纳。通常来讲，实验点并不能随便选取，而且模型的准确性与选取的实验点的数量也并非简单的线性关系。所以，在确定控制因子之后、建立RSM之前，需要合理设计参数空间中实验点的分布情况。选取实验点的方法被称为实验设计。一般而言，最希望得到的实验设计方案是在保证统计模型准确性的前提下所需实验点越少越好，这一点在大气模拟领域显得更为重要。因为大气模拟的一次实验就是进行一次控制情景下的空气质量模拟，受限于CMAQ计算时间的限制，进行大量的实验势必带来超高规模的计算负荷。在当前有限的计算能力制约下，如何利用有限的样本建立可靠的响应模型，如何设计实验方案以解决基于微小样本量的海量连续控制情景的预测问题，是RSM建模的关键所在。
对于表 1中的OC，BC和PM10这3个线性控制因子，由于响应关系是线性的，所以只需对每个控制因子在0和1两种控制情景下分别进行CMAQ模拟，即可获得求解该线性响应关系所需的实验样本。对于NOX_PP，NOX_IN，NOX_TR，SO2_PP，SO2_TR，NH3和VOC等7个非线性控制因子，当它们的参数水平变化范围为0~2时，假设参数连续变化的间隔为0.2，此时对应7维变量各自10个变化可能的控制情景组合的全排列高达107。通过遍历107种控制情景的CMAQ模拟获得求解7维变量非线性响应关系的实验样本的做法显然不可取。所以，高效的采样显得至关重要。本文根据拉丁超立方采样方法 (LHS，Latin hypercube sampling)[14]和计算仿真实验[15]，生成建立RSM模型所需要的实验矩阵 (表 2)，包含1个基本算例 (所有控制因子的参数水平均设置为1) 和150个控制算例的实验设计方案。
|
|
表 2 实验矩阵 Table 2 Experimental matrix |
杭州市PM2.5污染在冬季比较严重、夏季情况较好、春秋两季受外界影响较多，相对不稳定。鉴于此，选取比较有代表性且更多是受本地因素影响的冬、夏两季分别进行RSM建模。
将实验方案设计中的每一个样本所代表的控制情景通过CMAQ进行模拟，将模拟结果导入到由华南理工大学大气环境与污染控制重点实验室研发的RSM-VAT (response surface model-visualization and analysis tool) 区域大气污染控制可视化辅助决策工具，对模拟结果进行插值及非线性拟合，运算构建响应曲面模型，确定模型方程及其参数。在RSM建模部分，软件采用的是基于最大似然估计法的实验最佳线性无偏预测多维克里金插值方法 (multidimensional Kriging approach) 对模拟结果进行非线性拟合。
2.4 可靠性验证构建的RSM的误差主要来源于：一是CMAQ模拟本身具有不确定性。例如因为对污染过程认识不足而造成的物理建模偏差，数值计算过程中的近似处理误差，模拟输入数据包括排放清单和气象数据的不确定性等。二是在利用各控制情景下的CMAQ模拟数据进行RSM插值、统计、归纳、拟合时也必然会存在误差。所以，为了提高RSM模型的仿真性能，为空气质量决策提供科学依据，在反映连续变化的控制因子与污染物浓度的实时响应关系的RSM建立以后，需要对模型的可靠性进行验证。在此，采用标准误差NME (normalized mean error) 这一参数量化评估RSM的不确定性，计算公式如下：
| $ {\text{NME = }}\sum\limits_{i = 1}^N {\left| {{C_{ri}} - {C_{ci}}} \right|} /\sum\limits_{i = 1}^N {{C_{ci}}} , $ |
式中：N表示每一种控制情景中的网格数量；Cri表示在第i个网格中RSM的预测值；Cci对应第i个网格中CMAQ的模拟值。下面采用两种方法对RSM结果进行可靠性分析和验证。
1) 留一法交叉验证
留一法交叉验证是指在逐一剔除样本的同时，利用剔除之后的所有剩余样本建立RSM，然后应用该RSM预测在剔除样本的控制情景下的PM2.5浓度，最后比较该预测值与所剔除的CMAQ“真实值”。该方法的优点是验证可在拟合样本范围内进行，有多少个样本，就可以做多少次的验证，无需花时间准备额外的CMAQ模拟情景，但其缺点是无法验证样本范围外RSM系统的稳定性。
图 1为采用留一法交叉验证冬季和夏季150个样本得到的标准误差折线图。冬季RSM的平均标准误差为13.5%，夏季RSM的平均标准误差为10.3%, 这可以表明本研究的RSM预测结果基本可用。
|
Download:
|
|
图 1 冬季和夏季标准误差折线图 Fig. 1 Standard error line charts for winter and summer |
|
2) 外部验证
外部验证是指利用所有样本拟合RSM，然后通过设置样本范围以外的控制情景进行CMAQ模拟，比较同一情景下RSM预测值与CMAQ“真实值”对整个RSM系统进行检验。受限于CMAQ模拟计算时长，外部验证额外追加的控制情景不宜过多，通常可以设定一些靠近边界的数值以验证RSM靠近边界处的模拟性能，或者设置位于曲面边界以外，但处于政策之中的数值来评估其曲面的外边缘效果。根据以上原则，本文设置了所有控制因子逐一变化的情景 (表 3)。相比留一法交叉验证，该方法能评价超出实验样本范围时RSM预测特定情景的可靠性。
|
|
表 3 外部验证追加实验设计 Table 3 Additional experimental design for external validation |
运用CMAQ模型，对表 3中的每个控制情景进行模拟。图 2分别为10个追加的控制情景下冬季和夏季的标准平均误差折线图。由图可知，冬季RSM的平均标准误差为11.4%，最大标准误差为27.4%，夏季RSM的平均标准误差为9.8%，最大标准误差为20.9%。RSM外部验证的平均误差的量级与采用留一法交叉验证差不多，说明研究建立的RSM系统对位于曲面边界处的控制水平时的PM2.5预测也具有一定的可靠性。
|
Download:
|
|
图 2 冬季和夏季误差折线图 Fig. 2 Standard error line charts for winter and summer |
|
首先，定义PM2.5响应值为PM2.5浓度的变化率与某一污染物排放的变化率的比值。
在10%，50%和90%这3种控制水平下，依据之前建立的RSM模型，得到10种污染物的PM2.5响应值，如图 3所示。总的来说，PM2.5对OC、BC以及PM10的响应值在3种不同控制水平下几乎是相同的，说明OC、BC以及PM10排放对PM2.5的贡献是线性的；而在不同的排放水平下，PM2.5对其他7种污染物的响应值有较大差别，甚至会出现逆作用 (响应值随控制水平增大而减小) 和反作用 (响应值为负)，进一步说明其他7种污染物对PM2.5的贡献是非线性的；此外，同一种污染物在不同季节的PM2.5响应值也有比较大的差异。
在10种污染物排放中，PM10的PM2.5响应值最高，在冬季为0.18，在夏季为0.36。OC和BC的PM2.5响应值在冬季相对比较大，而在夏季相对比较小，说明杭州市空气中PM2.5的来源中，冬季一次污染物的比例大于夏季。冬季时，PM2.5对NOX和SO2的响应值随着控制水平的增大而增大，而NH3和VOC则是在控制水平为50%的时候达到最大，分别为0.14和0.08。这主要是由于NOX、SO2以及NH3作为二次无机气溶胶的前体物，VOC作为二次有机气溶胶的前体物，其在生成二次气溶胶的过程中具有很强的非线性特征。夏季时，VOC的影响增大，PM2.5响应值平均达到0.15以上，而其他几种污染物都出现了逆作用的现象。PM2.5对NH3的响应值随着控制水平的增加而减小，SO2在低控制水平下出现反作用，NOX在高控制水平下出现反作用，这说明正好存在一个拐点使得NOX和SO2在生成气溶胶的过程中有竞争关系。
|
Download:
|
|
图 3 冬季和夏季污染物的PM2.5响应值 Fig. 3 Response values of pollutant to PM2.5 in winter and summer |
|
表 4列出杭州市2020年控制因子的参数水平。表中基准年排放是指运行CMAQ建立RSM时所采用的2010年杭州市原始排放数据，控制情景Ⅰ是指基于目前的污染物减排政策和现有的执法力度的情况下的污染物排放水平，控制情景Ⅱ是指以牺牲一定的经济增长为代价严格控制污染物排放和加强执法力度的情况下的污染物排放水平。具体的排放数据预测方法参考文献[16]。归一化排放水平是指将目标年排放量比上基准年排放量的比率。
|
|
表 4 2020年两种控制情景下各污染物相比于2010年的变化率 Table 4 Pollutant change rate of 2020 to 2010 under two control scenarios |
利用杭州市PM2.5响应曲面模型，计算得到各种情景下2020年杭州市PM2.5浓度预测值，然后与2010年的数据比较，得到PM2.5的浓度变化率，也即2020年相对于2010年PM2.5浓度的变化量与2010年相应浓度的比值。在冬季情况下，控制相比于2010年的变化率情景的PM2.5减少19.90%，控制情景Ⅱ减少39.34%；夏季，控制情景Ⅰ的PM2.5减少30.73%，控制情景Ⅱ减少40.99%。这说明：控制情景Ⅱ的减排效果明显高于控制情景Ⅰ，夏季的减排效果好于冬季。
为考察不同污染物控制的效果，还设计了只将某一种污染物进行控制、其他污染物保持基准年排放水平的10种对比情景。此时，将各对比情景与参考情景进行比较，其差异便可认为是由于该污染物变化而带来的影响。每种排放源的控制效果如图 4所示。
|
Download:
|
|
图 4 冬季和夏季两种情景下的PM2.5浓度变化率 Fig. 4 Change rates of PM2.5 concentrations in winter and summer under two scenarios |
|
总的来讲，在冬季，OC、BC和PM10的控制对减少PM2.5浓度有比较大的影响，这说明冬季的PM2.5来自于一次排放的比例较大。同时可以看到，当控制政策收紧，NOX与NH3的影响就显现出来了。而在夏季，PM10的影响也比较大，但小于冬季；而电厂NOX和工业NOX的影响明显高于其他。而对于SO2，由于涉及大气环境从富氨到贫氨的转变，因而会出现负效应。NH3与VOC的影响普遍较小，也有可能是因为相关政策不够，其控制措施实施后本身排放量的变化不大。
4 结论本文基于CMAQ模拟得到对应的污染物浓度，借助响应曲面方法，统计、归纳出各种污染物排放控制因子与污染物浓度之间的函数关系，建立空气质量响应曲面模型RSM，并通过留一法交叉验证和外部验证对模型进行可靠性验证。在使用RSM分析不同污染物的PM2.5响应值以及评估杭州市未来减排情景 (Ⅰ、Ⅱ) 的控制效果时，得出如下结论：
1) OC、BC和PM10在不同的污染物控制水平下的PM2.5响应值几乎保持不变，体现出这3种污染物对PM2.5的线性响应特征。NOX、SO2、NH3和VOC的PM2.5响应值因不同的污染物控制水平而异，体现出明显的非线性响应特征，甚至会出现逆作用和反作用。
2) 在现有政策不变 (Ⅰ) 和控制政策收紧 (Ⅱ) 两种不同的减排情景下，杭州市2020年PM2.5浓度相对于2010年的减少量分别为：19.90%(Ⅰ-冬季)，39.34%(Ⅱ-冬季)，30.37%(Ⅰ-夏季)，40.99%(Ⅱ-夏季)，控制情景Ⅱ的减排效果明显高于控制情景Ⅰ，夏季的减排效果好于冬季。在冬季，OC、BC和PM10的控制对减少PM2.5有比较大的影响；当控制政策收紧时，NOX与NH3的影响才能明显显现。在夏季，PM10的影响也比较大，电厂NOX和工业NOX的影响明显高于其他，SO2会出现负效应。
| [1] | Sheehan P, Cheng E, English A, et al. China's response to the air pollution shock[J]. Nature Climate Change, 2014, 4(5):306–309. DOI:10.1038/nclimate2197 |
| [2] | Chen R, Zhao Z, Kan H. Heavy smog and hospital visits in Beijing, China[J]. American Journal of Respiratory and Critical Care Medicine, 2013, 188(9):1170–1171. DOI:10.1164/rccm.201304-0678LE |
| [3] | Pope Ⅲ C A, Burnett R T, Thun M J, et al. Lung cancer, cardiopulmonary mortality, and long-term exposure to fine particulate air pollution[J]. Jama, 2002, 287(9):1132–1141. DOI:10.1001/jama.287.9.1132 |
| [4] | Byun D, Schere K L. Review of the governing equations, computational algorithms, and other components of the Models-3 community multiscale air quality (CMAQ) modeling system[J]. Applied Mechanics Reviews, 2006, 59(2):51–77. DOI:10.1115/1.2128636 |
| [5] | Wang L T, Wei Z, Yang J, et al. The 2013 severe haze over southern Hebei, China:model evaluation, source apportionment, and policy implications[J]. Atmos Chem Phys, 2014, 14(6):3151–3173. DOI:10.5194/acp-14-3151-2014 |
| [6] | Zhang Y, Liu X H, Olsen K M, et al. Responses of future air quality to emission controls over North Carolina. Part Ⅱ:Analyses of future-year predictions and their policy implications[J]. Atmospheric Environment, 2010, 44(23):2767–2779. DOI:10.1016/j.atmosenv.2010.03.022 |
| [7] | Khayet M, Seman M N A, Hilal N. Response surface modeling and optimization of composite nanofiltration modified membranes[J]. Journal of Membrane Science, 2010, 349(1):113–122. |
| [8] | 严茹莎, 李莉, 安静宇, 等. 夏季长三角地区臭氧非线性响应曲面模型的建立及应用[J]. 环境科学学报, 2016, 36(4):1383–1392. |
| [9] | 劳苑雯, 朱云, JangC, 等. 基于响应面模型的区域大气污染控制辅助决策工具研发[J]. 环境科学学报, 2012, 32(8):1913–1922. |
| [10] | Baş D, Boyacı i H. Modeling and optimization I:Usability of response surface methodology[J]. Journal of Food Engineering, 2007, 78(3):836–845. DOI:10.1016/j.jfoodeng.2005.11.024 |
| [11] | Hanna S R, Chang J C, Fernau M E. Monte Carlo estimates of uncertainties in predictions by a photochemical grid model (UAM-IV) due to uncertainties in input variables[J]. Atmospheric Environment, 1998, 32(21):3619–3628. DOI:10.1016/S1352-2310(97)00419-6 |
| [12] | Moore G E, Londergan R J. Sampled Monte Carlo uncertainty analysis for photochemical grid models[J]. Atmospheric Environment, 2001, 35(28):4863–4876. DOI:10.1016/S1352-2310(01)00260-6 |
| [13] | 吴小芳.基于区域多尺度空气质量模型对杭州市大气污染的模式研究[D].杭州:浙江大学, 2014. |
| [14] | Iman R L. Latin hypercube sampling. Encyclopedia of quantitative risk analysis and assessment[M]. New Jersey: John Wiley & Sons, 2008: 408-411. |
| [15] | 邢佳.大气污染排放与环境效应的非线性响应关系研究[D].北京:清华大学, 2011. |
| [16] | 魏巍.中国人为源挥发性有机化合物的排放现状及未来趋势[D].北京:清华大学, 2009. |