2. 中国科学院大气物理研究所, 大气边界层和大气化学国家重点实验室, 北京 100029;
3. 中国科学院大学, 北京 100049;
4. 中国环境监测总站, 北京100012
2. State Key Laboratory of Atmospheric Boundary Layer Physics and Atmospheric Chemistry, Institute of Atmospheric Physics, Chinese Academy of Science, Beijing 100029;
3. University of Chinese Academy of Sciences, Beijing 100049;
4. China National Environmental Monitoring Centre, Beijing 100012
当前我国区域空气质量远未达标,大气重污染事件频发,京津冀区域尤为严重(王跃思等,2014).根据环保部发布的2013年中国环境状况公报显示,京津冀区域重度及以上污染天数占全年的平均比例为20.7%,显著高于同期长江三角洲的5.9%及珠江三角洲的0.3%(中华人民共和国环境保护部,2014).重污染天气会直接影响居民健康,并带来巨大的经济损失.依据空气质量流行病研究模型估算,2013年1月重污染过程大气污染浓度状况将导致北京地区增加早逝201例、呼吸病系统疾病1056例,相关的健康经济损失4.89亿元(谢元博等,2014).提前预告重污染天气过程的时段越长,就能争取更长的时间采取相应的响应措施,减少相应损失.
为有效地应对重污染天气过程,2013年国务院印发了《大气污染防治行动计划》,明确要求建立重污染天气监测预警体系,制定完善应急预案,将重污染天气应急响应纳入地方人民政府突发事件应急管理体系.同年5月,环保部印发了《城市大气重污染应急预案编制指南》,以提高城市大气重污染预防预警和应急响应能力.同年9月,中国气象局和环境部联合发布了《京津冀及周边地区重污染天气监测预警方案》(以下简称《方案》),依此对区域的重污染天气过程进行预警.
国际上空气质量预报的方法有两种,一种是基于统计方法建立的污染物浓度与气象参数间的统计预报模式,另一种则是建立在科学的理论和假设基础上,用数值方法描述大气中污染物的传输、扩散、化学反应及清除过程的数值预报模式(王自发等,2006).与统计预报相比,数值预报的优势在于:它可以深入分析污染过程演变特征,解析污染物的来源和去向,以及模拟预测应急控制效果等.对模式的预报能力进行评估,是实现预报预警的首要任务.目前,区域重污染天气过程预报评估方法相对于污染物浓度或者AQI的预报评估方法(陈焕盛等,2013;胡鸣等,2015)较少.因此,本文拟建立一种适用于京津冀区域重污染天气过程预报的评估方法,利用2013年和2014年日均值观测资料及NAQPMS系统AQI预报产品,评估模式系统提前1~3 d对重污染天气过程的预报能力.
造成京津冀区域重污染天气的原因有很多,污染源排放不确定性与未知性是造成京津冀区域重污染天气的内在原因,气象条件则是外在原因(王跃思等,2014;吴兑等,2012).但在模式系统中,排放相对稳定的情况下,气象条件往往是决定性因素.在实际污染预报中,业务部门往往会结合未来一周气象预报产品对系统的原始结果加以订正,从而提高预报的准确率.因此,本文在评估模式系统对重污染天气过程预报能力的基础上,从气象方面探讨导致预报产生偏差的原因,以期提高业务预报的能力,为重污染天气应急控制和响应提供支持.
2 资料与方法(Materials and methods) 2.1 空气质量模式系统嵌套网格空气质量预报模式(NAQPMS)是由中国科学院大气物理研究所自主开发的、基于“一个大气”理念设计的第三代空气质量模式.模式全面考虑了空气污染物在大气中的平流、扩散、干湿沉降及化学转化等过程(王自发等,2006).NAQPMS模式系统参数化方案详细介绍见文献(Li et al.,2012; 朱莉莉等,2015).
采用WRF模式为空气质量预报模式提供的主要气象参数包括气温、湿度、降水、风向、风速和云量等.
初始边界条件采用美国NCEP全球预报分析资料GFS,其时间间隔为24 h,空间分辨率为0.5°×0.5°.
排放清单采用亚洲区域大气污染物排放清单REAS2.1(Kurokawa et al.,2013),该清单空间分辨率为0.25°×0.25°.基准年为2007年,并根据GDP增长情况简单更新到预报年.
预报系统的区域设置覆盖了中国东部地区,模式区域的中心经纬度为25°N、115°E,网格数为120×162,水平网格格距为45 km.采用适合我国中低纬度特点的LAMBERT投影方式.垂直方向上采用地形追随高度坐标的形式,不等距地分为20层,边界层内约有7~8层,近地层中心高度约47 m,模式层顶的高度为20 km.
模式系统每日预报未来一周京津冀区域13个城市空气质量状况,每日预报以前一天北京时间20:00为起始时间,预测未来一周污染物浓度及AQI.本文主要分析此模式系统第2~7 d的预报结果.
2.2 观测数据目前,我国以空气质量指数(AQI)评价每日空气质量状况,它将PM2.5、PM10、SO2、NO2、CO和O3等6项污染物用统一的标准评价.按照AQI指数可将空气质量分为优、良、轻度污染、中度污染、重度污染及严重污染6个等级(环境空气质量标准GB 3095—2012).本文以AQI指数大于200作为重污染的统计标准.
2013年和2014年京津冀区域每日AQI资料来自于中国环境监测总站,逐日天气图来自于韩国气象厅.根据地理位置将所有城市分为河北北部(张家口、承德和秦皇岛)、北京、天津、河北中部(廊坊、唐山、沧州和保定)及河北南部(石家庄、衡水、邢台和邯郸).
2.3 污染过程评价标准 2.3.1 重污染天气过程标准依据《方案》,当预报出京津冀区域3个及以上城市AQI连续3 d大于200时,启动重污染天气预警.因此,本文定义京津冀区域3个及以上城市AQI连续3 d大于200为一个重污染天气过程.
2.3.2 预报能力评价标准本文主要从起始时间、持续时间和空间范围3方面综合评价重污染天气过程的预报能力.京津冀区域至少3个城市AQI大于200的第1 d为重污染天气过程的起始时间,最后一天为重污染天气过程的终止时刻,起始时间与终止时刻间包含的天数为持续时间.所有连续3 d及以上AQI大于200的城市都被认为处在该重污染天气过程空间范围内.
以实际发生重污染天气过程为参照,当预报的区域污染过程达到启动重污染预警标准且预报污染过程空间范围达实际的60%时,若预报起始时间与实际一致,则评定该次重污染天气过程“预报成功”;若起始时间较实际提前,评定为“早报”;若起始时间较实际滞后(滞后不能超过2 d),评定为“晚报”;而当预报的城市少于3个、持续时间少于3 d或污染过程空间范围低于实际的60%时,判为“漏报”.
参照《方案》,本文将AQI等于200作为标准Ⅰ;同时,通过对两年AQI预报结果统计发现,模式预报AQI达到150以上,实际出现重污染天气过程的概率较大,因此,将AQI等于150作为一个预警临界值,即标准Ⅱ.基于标准Ⅰ和标准Ⅱ对观测重污染天气过程预报效果进行评估.
预报成功、早报和晚报都能起到警示作用,本文将预报成功、早报和晚报归为“报出”一类.用预报准确率P(式(1))来评价2013年和2014年模式系统提前1 d、2 d和3 d对重污染天气过程预报能力.
|
(1) |
式中,Cb、Cs分别为报出重污染天气过程次数和实际重污染天气过程总次数.
3 结果与讨论(Results and discussion) 3.1 京津冀区域重污染分类与概况参考李令军等(2012)的分类方法,本文将京津冀区域的重污染天气过程分为静稳型、沙尘型和特殊型污染.静稳型污染指由于出现持续不利于扩散的静稳条件导致污染物大范围积累,使得污染物达到重污染水平.沙尘型污染指重污染天气过程有1 d及以上受外来沙尘或者局地扬尘造成的颗粒物重污染.特殊型污染是指重污染天气过程中有1 d及以上受到特定季节秸秆燃烧及春节期间鞭炮燃放等特殊原因的影响造成的重污染.
结合MODIS火点分布(NASA,https://firms.modaps.eosdis.nasa.gov/firemap/)、污染期间6项污染物浓度变化、天气形势、边界层垂直结构和地面气象要素对重污染过程分析统计表明:京津冀区域2013年和2014年共发生30次重污染天气过程(表 1),秋、冬季发生频次最多,达24次,其次是春、夏季,分别为4次和2次.这些重污染天气过程的首要污染物多为PM2.5或PM10.重污染天气过程平均持续时间为6 d,最严重一次发生于2013年1月,持续28 d AQI>200,污染覆盖除张家口和承德外的整个区域.河北北部、北京、天津、河北中部及河北南部受污染影响频繁程度呈递增的趋势.这些重污染天气过程主要为静稳型污染,达24次,2次沙尘型污染(2013年3月5—9日、2013年3月15—18日)和4次特殊型污染(秸秆燃烧:2013年6月15—17日、2013年10月3—6日;鞭炮燃放:2014年1月27日—2月2日、2013年2月9—15日).京津冀区域静稳型污染主要发生在秋冬季,秋季的静稳型污染一般持续3~5 d,影响范围大.冬季静稳型污染平均持续7 d,持续时间和影响范围均大于秋季,沙尘型污染主要发生在春季,而特殊型污染主要发生在秸秆燃烧季(如秋季)及春节期间.
| 表 1 2013年和2014年重污染天气过程信息表 Table 1 Information of heavy pollution in 2013 and 2104 |
随着重污染天气越来越受到关注,重污染天气预报也成为预报部门的核心业务,评估模式系统对重污染天气的预报效果就显得尤为重要,它不仅使我们对预报能力有大致了解,对于预报结果更为有效应用、预报准确率的提高及模式系统自身的改进和完善也具有重要的意义.因此,本文通过对比分析提前1、2和3 d对30次重污染天气过程的预报能力,发现以下几个特点.
在起始时间方面,图 1a给出了系统提前1、2和3 d的预报情况.提前1 d的预报效果好于提前2 d和3 d,重污染天气过程预报成功的次数分别为13、11次和11次.用标准Ⅰ进行评估,提前1、2和3 d重污染天气过程预报准确率均可达到57%以上.用标准Ⅱ进行评估,预报准确率均可达到70%以上,图 1b给出了不同季节的预报效果,不同季节的准确率分别为75%、0、90%和60%,秋季预报效果较其他季节好.如表 2所示,24次静稳型污染过程,提前1、2和3 d均可达18次以上,6次沙尘型和特殊型污染则最多可报出3次.
![]() |
| 图 1 2013年和2014年提前1、2和3 d预报结果(a)与不同季节预报效果(b) (O为实际污染的次数;标准I-1、标准I-2、标准I-3和标准Ⅱ-1、标准Ⅱ-2、标准Ⅱ-3含义同表 1) Fig. 1 Forecasting ability of one-day, two-day and three-day ahead (a) and in different seasons (b) |
| 表 2 不同类型重污染天气过程预报效果 Table 2 Forecasting ability on different types of heavy pollution |
在持续时间方面,提前1、2和3 d预报的重污染天气过程一般持续时间为3~4 d,最长可达6 d.对于冬季持续时间较长的重污染天气过程,系统可能预报出重污染天气过程的开始或中间时段,这时需要密切关注该时间段每日的预报结果,综合发布未来几日的预警信息.以2013年1月4日开始持续时间长达28 d的重污染天气过程为例,虽未能正确预报出此次过程的起始时间,但系统陆续报出该时段内7—11日、13—15日、18—23日和26—31日的重污染天气过程,考虑模式中静态的排放清单难以反映京津冀区域冬季复杂的污染排放可能是造成预报持续时间不连续的主要原因之一.
在空间范围预报方面,系统预报出的空间范围通常大于实际发生的重污染天气过程影响的空间范围,对秋冬季由于静稳条件导致的大范围的重污染天气过程,预报的空间范围通常可达实际的80%以上.如图 2所示,对河北中南部的预报效果较好,特别是保定、石家庄和邯郸,实际受到26、30和24次重污染天气过程影响,报出次数分别为20、22和18次.秦皇岛、北京、邢台则有漏报趋势,秦皇岛3次污染过程均未报出,北京17次污染过程漏报8次,邢台28次污染过程漏报了12次.REAS2排放清单对京津冀区域不同城市污染物排放估计与实际偏差可能是导致这一现象的主要原因之一.
![]() |
| 图 2 2013和2014年京津冀区域各城市受重污染天气过程影响况与预报效果对比 Fig. 2 Comparison of heavy pollution in Beijing-Tianjin-Hebei area between observation and forecast in 2013 and 2014 |
研究表明,气象条件包含天气形势和气象要素,天气形势决定了气象要素的分布、变化情况及大气扩散的能力(缪育聪等,2015).当华北地区出现严重污染时,地面通常主要受低压(孟燕军等,2002;谢付莹等,2010)、弱高压或华北地形槽(高书然等,1982)类系统控制.当华北地区500 hPa上空为偏西或西南气流控制,850 hPa上空受暖平流(西南气流)影响或暖舌影响时,华北地区比较容易发生重污染天气(谢付莹等,2010).500 hPa天气形势的变化对低层天气系统的维持与发展具有很大的影响.而850 hPa则大致相当于近地面边界层顶的高度,受局地地形影响较近地面小.因此,本文利用韩国气象厅8:00和20:00逐日的500 hPa和850 hPa天气图与对应的WRF模式预报结果进行对比,探讨天气形势对重污染天气过程早报、晚报及漏报的可能原因.
3.3.1 早报标准Ⅱ评估提前1、2和3 d共统计出11次早报现象,以2014年10月23—25日影响京津冀区域典型的静稳型污染过程为例,结合不同高度层的实况天气图来分析导致系统早报可能的原因.在此次过程中,22日,京津冀区域500 hPa高空以经向环流为主,受冷槽控制;850 hPa高度上,风场上由西北气流逐渐转为偏西气流,从温度场上看,锋区已经东移,暖脊尚未影响到本区;地面气压场受冷中心控制,从不同高度天气系统的配置来看,并不利于污染物的累积,故该日京津冀区域空气质量并未达到重度污染级别.23—25日,500 hPa高空,冷槽位于贝加尔湖以北,并缓慢东移,本区受冷槽前部偏西气流影响,以纬向环流为主;850 hPa高压中心位于黄淮地区,本区受高压顶部偏西气流影响,以辐散环流为主,并且有暖中心配合;在地面,冷高压迅速变性东移,而内蒙古地区出现一低压中心,该低压中心不断向东北及西南扩展,本区即受该低压中心向南伸展的低压辐合带控制.从不同高度的天气系统配置来看,23—25日,高空无明显冷空气活动,中低层层结日趋稳定,地面辐合形势明显,扩散条件不利,出现有利于污染物累积的典型天气形势,从而导致23—25日出现一次大范围静稳型的重污染天气过程.
模式系统于2014年10月20日(提前3 d)预报出上述污染过程,但预报结果较观测早1 d,持续时间4 d,空间范围达观测的60%.对比22日实况与提前3 d的WRF预报结果,22日WRF预报出500 hPa高空和地面天气系统位置和强度与观测较一致,而850 hPa则预报该区域已受暖脊控制(图 3).23—25日,不同高度天气系统位置和强度与观测较一致.可见,22日暖脊相对于实况提前影响京津冀区域可能是导致早报1 d的原因.
![]() |
| 图 3 2014年10月22日08时850 hPa实况(a)、预报(b)天气形势 Fig. 3 Actual (a) and forecasted (b) synoptic situation at 850 hPa at 08:00, 22 October 2014 |
在统计出的7次晚报现象中,挑选2014年1月4—7日的一次覆盖保定、邯郸、衡水、石家庄及邢台等城市的重污染天气过程为例,来分析导致系统晚报可能的原因.在该次过程中,4日,500 hPa高空位于槽后,受西北气流影响;850 hPa高空,风场上由偏西气流逐渐转为偏南气流,地面受低压系统控制,偏南风影响.5—7日,500 hPa高空随着冷槽东移,本区受偏西气流影响;中低层受偏南气流影响,层结日趋稳定.简而言之,4—7日大气扩散条件较差,污染物易在近地面累积.
对比观测与预报结果发现,提前1 d预报成功,提前2 d和3 d预报为晚报.对比提前1、2、3 d预报结果,4—7日,预报的500 hPa天气系统的位置和强度与实况较为一致,且5—7日预报的850 hPa情况也与实况也较为接近,但4日(图 4),850 hPa高空提前1 d预报与实况一致,而提前2 d和3 d预报出该区域位于槽后,受西北气流影响.可见4日,脊相对于实况滞后影响该区域可能是导致晚报1 d的原因.
![]() |
| 图 4 2014年1月4日08时850 hPa实况(a)及提前1 d(b)、2 d(c)和3 d(d)预报的天气形势 Fig. 4 Actual (a), and one-day (b), two-day (c), and three-day (d) ahead forecast of synoptic situation at 850 hPa at 08:00, 4 January 2014 |
通过对早报和晚报个例与实况对比,发现从气象方面能解释部分静稳型污染早报和晚报的原因.早报可能的主要原因为:当850 hPa天气形势不利于污染物扩散时,若WRF对500 hPa槽后西北气流预报较实况弱,可能会导致早报(2013年10月27—29日提前1 d预报结果);WRF对850 hPa暖脊或西南气流提前预报,该区域提前受逆温层的影响,抑制了污染物的垂直扩散,从而可能会导致早报(2014年10月23—25日提前3 d预报结果).相反,若WRF对500 hPa槽后西北预报较实况强(2014年1月10—19日提前3 d预报结果);850 hPa暖脊或西南气流滞后影响该区域,则可能会导致晚报(2014年1月4—7日提前3 d预报结果).
WRF预报出不利于污染扩散的静稳型天气系统的提前和滞后会导致重污染天气过程早报或晚报.因此,业务预报时要特别关注秋冬季该类型污染,结合模式系统中WRF及欧洲预报中心等不同气象模式的预报结果,对NAQPMS模式系统的原始预报结果进行修订.
3.3.3 漏报漏报分为两种,一种是对静稳型污染漏报,一种是对沙尘型或特殊型污染的漏报.模式系统中模式参数化方案和污染源清单都是既定的,对于突发的沙尘型和特殊型污染的漏报现象难以从天气形势预报上寻查原因. 例如,2013年6月15—17日的一次污染过程,15—17日天气形势总体不利于污染物的扩散,从MODIS卫星火点分布图(NASA)可看出,15—17日,石家庄、衡水、邢台及南边的河南省有大量的火点分布,故局地的生物质燃烧,使得石家庄、衡水和邢台污染物的排放量突增,加重了污染级别是导致该次污染的可能原因.因此,业务预报时要特别关注秸秆燃烧季节、春节期间烟花爆竹排放增加或重大活动减排信息,结合天气形势对原始系统预报结果进行人为的客观订正.对比静稳型污染漏报个例实况与预报结果之间的差异发现,500 hPa和850 hPa预报的天气系统的位置和强度都与观测较为相似.
对于突发的沙尘型和特殊型污染的漏报,预报时需密切地关注季节排放源的变化,修订原始结果和预报结果,提高预报的准确率;而对于漏报的静稳型污染,未来通过加强卫星资料应用(朱彬等,2010)和排放源方面(程兴宏等,2008;谢博文等,2013)的工作来提高该类预报效果.
4 结论(Conclusions)本文发展了一种适应于京津冀区域重污染天气过程预报评估的新方法,该方法分别从起始时间、持续时间及空间范围3个方面评估了2013年和2014年京津冀区域30次重污染天气的预报结果,分析NAQPMS模式系统提前1、2和3 d对重污染天气过程的预报能力,主要得出以下结论:
1)NAQPMS模式系统能预报出大部分的重污染天气过程.在起始时间方面,提前3 d的重污染天气过程预报准确率可达到57%;在持续时间方面,预报的重污染天气过程一般持续时间大致为3~4 d,最长可达到6 d;在空间范围方面,系统对于大范围的重污染天气过程的空间范围预报效果较好,对河北中南部的预报效果最好,对秦皇岛、北京和邢台则有漏报现象.
2)业务预报经验发现,预报系统对污染过程预报有低估的趋势,为了更好地利用预报结果,将AQI=150作为标准Ⅱ对预报能力进行评估,预报准确率将会提高.提前3 d的重污染天气过程预报准确率提高至70%以上,静稳型污染预报准确率从58%提高至80%.
3)对静稳型污染预报最易出现早报和晚报现象,从气象上可以解释部分早报和晚报的原因.造成该类污染早报可能主要有以下两个原因:当850 hPa天气形势不利于污染物扩散时,WRF对500 hPa高空槽后西北气流预报较实况弱;或WRF对850 hPa暖脊或西南气流提前预报;相反,则可能导致晚报.
| [${referVo.labelOrder}] | 陈焕盛, 王自发, 吴其重, 等.2013. 空气质量多模式系统在广州应用及对PM10预报效果评估[J]. 气候与环境研究 , 2013, 18 (4) : 427–435. |
| [${referVo.labelOrder}] | 程兴宏.2008.空气质量模式"源同化"模型及排放源影响效应研究[D].北京:中国科学院研究生院.11-41 |
| [${referVo.labelOrder}] | 高书然, 李郁竹.1982. 空气污染的天气形势分析和预报[J]. 气象 , 1982, 8 (1) : 33–34. |
| [${referVo.labelOrder}] | 胡鸣, 赵倩彪, 伏晴艳.2015. 上海环境空气质量预报考核评分方法研究和应用[J]. 中国环境监测 , 2015, 31 (4) : 54–57. |
| [${referVo.labelOrder}] | Kurokawa J, Ohara T, Morikawa T, et al. 2013. Emissions of air pollutants and greenhouse gases over Asian regions during 2000-2008:Regional Emission inventory in Asia (REAS) version 2[J]. Atmospheric Chemistry and Physics Discussions , 13 (4) : 10049–10123. DOI:10.5194/acpd-13-10049-2013 |
| [${referVo.labelOrder}] | 李令军, 王英, 李金香, 等.2012. 2000-2010北京大气重污染研究[J]. 中国环境科学 , 2012, 32 (1) : 23–30. |
| [${referVo.labelOrder}] | Li J, Wang Z F, Zhuang G, et al. 2012. Mixing of Asian mineral dust with anthropogenic pollutants over East Asia:a model case study of a super-duststorm in March 2010[J]. Atmos Chem Phys , 12 : 7591–7607. DOI:10.5194/acp-12-7591-2012 |
| [${referVo.labelOrder}] | 缪育聪, 郑亦佳, 王姝, 等.2015. 京津冀地区霾成因机制研究进展与展望[J]. 气候与环境研究 , 2015, 20 (3) : 356–368. |
| [${referVo.labelOrder}] | 孟燕军, 程丛兰.2002. 影响北京大气污染物变化的地面天气形势分析[J]. 气象 , 2002, 28 (4) : 42–47. |
| [${referVo.labelOrder}] | 吴兑.2012. 近十年中国灰霾天气研究综述[J]. 环境科学学报 , 2012, 32 (2) : 257–269. |
| [${referVo.labelOrder}] | 王跃思, 张军科, 王莉莉, 等.2014. 京津冀区域大气霾污染研究意义、现状及展望[J]. 地球科学进展 , 2014, 29 (3) : 388–396. |
| [${referVo.labelOrder}] | 王自发, 谢付莹, 王喜全, 等.2006. 嵌套网格空气质量预报模式系统的发展与应用[J]. 大气科学 , 2006, 30 (5) : 778–790. |
| [${referVo.labelOrder}] | 谢博文, 杨少英, 李永秀.2013. 2种气溶胶排放源对模拟中国地区气溶胶分布的影响[J]. 环境科学与技术 , 2013, 36 (8) : 187–193. |
| [${referVo.labelOrder}] | 谢付莹, 王自发, 王喜全, 等.2010. 2008年奥运会期间北京地区PM10污染天气形势和气象条件特征研究[J]. 气候与环境研究 , 2010, 15 (5) : 584–594. |
| [${referVo.labelOrder}] | 谢元博, 陈娟, 李巍.2014. 雾霾重污染期间北京居民对高浓度PM2.5持续暴露的健康风险及其损害价值评估[J]. 环境科学 , 2014, 35 (1) : 1–8. |
| [${referVo.labelOrder}] | 朱彬, 苏继锋, 韩志伟, 等.2010. 秸秆焚烧导致南京及周边地区一次严重空气污染过程的分析[J]. 中国环境科学 , 2010, 30 (5) : 585–592. |
| [${referVo.labelOrder}] | 中华人民共和国环境保护部.2014.中国环境状况公报[OL].北京:中华人名共和国环境保护部.2014-06-05.http://jcs.mep.gov.cn/hjzl/zkgb/2013zkgb/201406/t20140605_276521.htm?COLLCC=3178359310#38; |
| [${referVo.labelOrder}] | 朱莉莉, 晏平仲, 王自发, 等.2015. 江苏省级区域空气质量数值预报模式效果评估[J]. 中国环境监测 , 2015, 31 (2) : 17–23. |
2016, Vol. 36





