一种城市颗粒物时空管控方法

1.本发明涉及一种城市颗粒物时空管控方法，属于大气环境监管领域，涉及城市颗粒物管控方法。


背景技术：

[0002][0003]
目前中国的区域环境管控工作主要针对各发达城市群，划分结果较为粗略，其空间维度主要针对省、市级别，时间尺度多数以年、月、日为单位，同时也忽略了空气污染的时空异质性。因此，如何有效把握城市各区域颗粒物的变化规律，并提出针对性的管控方法是保障城市空气质量的关键。


技术实现要素：

[0004]
为解决上述问题，本发明提出一种城市颗粒物时空管控方法。
[0005]
本发明提供了一种新的管控框架，对城市进行中尺度(时间单位为日，空间单位为千米)的管控区域划分，分析各管控区域颗粒物波动幅值与各地理密度值的关系，得到城市各区域各时间段相关性最大的排放源类型，促进城市颗粒物时空精细化管控。
[0006]
本发明首先利用普通克里金插值法分析得到城市污染时空网格分布(1km
×
1km单位网格)，并以日均浓度值和日间浓度变化幅值为特征展开聚类分析，划分中尺度的颗粒物管控区域，利用单日浓度值(最高10％的浓度均值)标记各污染管控等级，分析各管控区域与当地地理因子核密度值的时空地理加权回归系数，得到颗粒物波动与各类排放源的关系，可为环保部门提供较为精准、科学的管控参考对象，推动城市颗粒物时空精细化管控，从而达到提升城市空气质量的目的。
[0007]
本发明的一种城市颗粒物时空管控方法，具体步骤如下：
[0008]
1、区域网格划分处理。
[0009]
将城市区域投影至平面坐标系内，根据城市的边界或需管控的区域划定最大矩形框，并展开中分辨率的网格划分，按照最小单元格分为m行和n列。
[0010]
2、时空污染分布分析。
[0011]
获取城市各乡镇监测站的地理信息和颗粒物浓度数据，计算每日各监测点位的均值，并模拟样本的变异函数，运用普通克里金插值法计算得到日均空间污染分布数据。
[0012]
3、中尺度管控区域划分。
[0013]
利用日均空间污染分布网格数据，进行日间均值浓度差值计算，以代表第m天的第i行j列单元格的浓度均值，代表第m-1天的第i行j列单元格的浓度均值，代表第m天与第m-1天的第i 行j列单元格的浓度变化幅值：
[0014][0015]
以上述网格的浓度变化幅值和日均浓度值为特征值，利用层次聚类算法进行区域
划分，其中聚类的距离函数如下所示：
[0016][0017]
式中代表两个不同单元格(p1和p2)在第m天的距离，和表示两单元格在第m天的浓度均值，和表示两单元格在第m天的浓度变化幅值。
[0018]
在不考虑气象条件等短暂性影响外，采用浓度最大值作为评价标准，一定程度上避免监测误差和短暂性意外污染源的影响。
[0019]
因此，分析得到各区域划分结果后，需利用相关等级评判该区域的污染风险。由于分析的区域范围为中尺度，除去外界污染扩散影响和季节性变化等重要因素，为充分表达区域性污染程度，将污染最大均值通过归一化处理后，划分5个等级，等级计算方法如下：
[0020][0021]
上式中，ci表示第i个区域的浓度均值，c
min
表示各管控区的最小浓度均值，c
max
表示各管控区的最大浓度均值。
[0022]
得到中尺度管控区域划分结果后，分别对各管控区展开潜在颗粒物时空排放源相关性分析，因此重新划分区域，以中尺度管控区域东南西北端划定矩形框，按照单元格划分矩形网格。
[0023]
4、核密度估计。
[0024]
各类poi以及道路分布情况与颗粒物浓度变化相关性较大。为分析各区域典型poi和道路的聚集性，可估计不同核密度值来描述各类的分布情况。对于独立分布的对象(道路和poi)，核函数密度公式可表达为：
[0025][0026]
式中的n代表统计对象的总数；k表示核密度计算中使用的对象数量的核函数，此处选用高斯核函数；h代表带宽，也可表示搜索半径，根据各类场所对周边颗粒物浓度的影响范围。
[0027]
5、时空自相关性分析。
[0028]
区域性污染扩散是空气污染的重要特征之一，因此颗粒物浓度及其影响因素通常具有空间相关性。莫兰指数可用于确定空间相邻区域的颗粒物浓度及其影响因素是否相似，其计算方式如下：
[0029][0030]
上式中，n为管控区域内的总网格数，xi和xj表示单元格i和单元格j的颗粒物浓度值，为所有空间单元格颗粒物浓度均值，ω
ij
为各单位之间的权重系数。
[0031]
利用莫兰指数分析颗粒物浓度在0～23小时的空间自相关情况，观察其时空异质性特点，说明引入时空地理加权回归模型的重要性。
[0032]
6、时空地理加权回归分析；
[0033]
大气污染变化影响因素可分为两部分，外部污染扩散和本地活动排放因素。由于研究范围程度较小，而冬季的污染浓度变化普遍较大，直接从各时空维度浓度分析会存在较多干扰因素。因此，为分析微尺度下各类典型poi和道路类型对当地颗粒物的影响程度，以研究区域内所有监测站的小时浓度均值作为区域整体浓度值，计算各监测点位的小时浓度均值与区域整体浓度值的差值，得到各监测点位所在单元格的浓度变化幅值，该幅值能充分体现监测点位的污染变化规律。
[0034]
gtwr模型(时空地理加权回归模型)通过构建权重矩阵来解释参数估计的时空非平稳特性，该权重矩阵基于观测值和所有其他观测值之间的距离确定相关性，其计算方法如下：
[0035][0036]
式中，yi表示第i个位置的预测值，即第i个监测点位的小时浓度均值与区域整体浓度值的差值；x
ik
为第i个监测点位归一化后的各类核密度值；(ui,vi,ti)表示第i个时空参数，包括网格的行数和列数以及0～23时；β0(ui,vi,ti)表示第i个时空维度下的模型回归参数；βk(ui,vi,ti)表示第i个时空维度下的模型与第k个维度的回归参数；εi表示第i个时空维度的随机性误差。
[0037]
最后，计算得到的βk(ui,vi,ti)可表示各类输入因子的影响大小，此处可作为各类poi或道路等级对当地颗粒物的影响程度，并在时空维度进行加权，得到各时间段对应典型因子的影响系数，最终参考该影响系数可得出精准的时空管控方案。
[0038]
优选地，步骤1所述的最小单元格尺寸为1km
×
1km。
[0039]
优选地，步骤3所述的浓度最大值为10％。
[0040]
优选地，步骤3所述的单元格尺寸为200
×
200m。
[0041]
本发明首先利用普通克里金插值法分析得到城市各区域时空污染分布情况，并以污染浓度均值和波动幅值为特征，利用层次聚类法划分管控区域，并以各管控区内的浓度均值(最大10％)为指标，进行管控区等级标记。将各中尺度管控区进一步精细划分为200
×
200m 的单元格，计算各类poi和道路的核密度分布，并利用莫兰指数验证颗粒物的时空异质性特点，采用时空地理加权模型拟合颗粒物波动值与各类核密度归一化值之间的关系，判断不同区域不同时间段的各类场所对颗粒物变化的影响程度，为环保部门提供参考性管控对象，最大程度保障市民活动和空气质量的和谐发展。
[0042]
本发明优点：从中尺度到微尺度研究管控区域的划分范围和管控等级以及主要排放源类型，可有效提高环保部门的管控效率，针对重要管控时间段，在保障城市空气质量，精细管控各类场所。
附图说明
[0043]
图1是本发明方法的流程图。
[0044]
图2是区域网格划分图。
[0045]
图3是春季管控区域划分结果。
[0046]
图4是夏季管控区域划分结果。
[0047]
图5是秋季管控区域划分结果。
[0048]
图6是冬季管控区域划分结果。
[0049]
图7是精细化网格划分图。
[0050]
图8(a)～图8(d)是典型对象的核密度分布图，其中，图8(a) 是景区核密度分布图，图8(b)是医疗服务核密度分布图，图8(c) 是一级道路核密度分布图，图8(d)是二级道路核密度分布图。
[0051]
图9是莫兰指数和z值变化图。
[0052]
图10是各因子的影响系数变化图。
具体实施方式
[0053]
下面结合附图，进一步说明本发明的技术方案。
[0054]
城市监测站数据为杭州市2020年68个乡镇监测站的逐小时监测数据，包括监测站地理信息、气象数据、pm
2.5
和pm
10
逐小时数据；选用的杭州区县行政区地图为2020年版本；poi和道路数据也均为 2020年版本，其中poi数据主要包括餐饮服务、景区、购物服务、金融保险、科教文化、汽车服务、生活服务、政府机构、住宅、道路附属设施、交通设施服务以及工厂，道路数据包括国道、省道、高速、城市一、二、三、四级道路。另外，本发明选择pm
2.5
为颗粒物研究对象并展开说明。
[0055]
如图1所示，一种城市颗粒物时空管控方法，具体步骤如下：
[0056]
1、区域网格划分处理。
[0057]
如图2所示，根据2020年杭州市行政区划分地图映射至平面坐标系(wgs_1984_utm_zone_51n)中，对杭州市主城区进行1km
ꢀ×
1km单元网格划分，经度方向共58个单元格，纬度方向共24个单元格。其中，为选取萧山(大江东区域)，不可避免地涉及海宁市部分领土，但该区域地网格面积约为35km2，仅占总研究面积的2.5％，可忽略不计。
[0058]
2、时空污染分布分析。
[0059]
获取城市68个乡镇监测站的地理数据和颗粒物浓度数据，计算每日各监测点位的均值，根据季节特征分类，对四批数据模拟样本的变异函数，并运用普通克里金插值法计算日均空间污染分布数据。
[0060]
3、中尺度管控区域划分。
[0061]
利用日均空间污染分布网格数据，进行相邻日间均值浓度差值计算，以代表第m天的第i行j列单元格的浓度均值，代表第 m-1天的第i行j列单元格的浓度均值，代表第m天与第m-1天的第i行第j列单元格的浓度变化幅值：
[0062][0063]
以上述单元格的浓度变化幅值和单元格日均浓度均值为特征值，利用凝聚层次聚类算法进行区域划分，其中聚类的距离函数如下所示：
[0064][0065]
式中代表两个不同的单元格(p1和p2)第m天的距离大小，和表示两单元格第m天的浓度均值，和表示两单元格第m天的浓度变化幅值。
[0066]
为贴近大气污染管控工作的实际需求，管控区域不可能零散而孤立存在，希望划
分结果的每个区域达到最小行政区面积，根据分层聚类算法结果，采用聚类最小的单元格数量不小于15(区域面积至少 15km2)的约束条件，合并较小子区域，达到方便统一管控目的。
[0067]
以各管控区域内最大10％的浓度数据均值为评价标准，在通过归一化后划分为5个等级，计算方法如下：
[0068][0069]
上式中，ci表示第i个区域的浓度均值，c
min
表示各区域的最小浓度均值，c
max
表示各区域的最大浓度均值。
[0070]
最终，根据季节类型分别划分了4个不同管控区域版本，如图3、图4、图5和图6所示。
[0071]
q值检验法被广泛应用于衡量生态现象的空间异质性程度。该指数可以估计集群之间的差异和集群内的相似性，q值统计方法是评估 pm
2.5
污染时空聚类有效性的可靠指标，该方法的计算方程如下：
[0072][0073]
式中，n表示分类的总数；nh表示第h个分类；为第h个分类的浓度方差值；σ2表示所有单元格的pm
2.5
浓度的方差。q值在0 到1之间，该值越大表示污染区域划分效率越高。
[0074]
根据上述结果将各季节下各区域内包含的各监测站分为不同类，以3个月为一个季节周期计算各区域的q值，计算过程中，考虑到涉及方差的计算，除去无监测站或仅有一个监测站的分析区域。具体结果如下表1所示。
[0075]
表1
[0076][0077]
为评价聚类结果的有效性，参考中国pm
2.5
污染划分的相关研究,q值阈值设置为0.27。而上述提出方法q统计量高达0.38，表明季节性浓度变化特征聚类方法是正确识别联合污染控制区域的有效方法。选取主城区管控区域进一步划分，具体如图7所示。
[0078]
4、核密度估计。
[0079]
各类poi以及道路分布情况与颗粒物浓度变化相关性较大。为分析各区域典型poi和道路的聚集性，可估计不同核密度值来描述各类的分布情况。对于独立分布的对象(道路和poi)，核函数密度公式可表达为：
[0080][0081]
式中的n代表统计对象的总数；k表示核密度计算中使用的对象数量的核函数，此处选用高斯核函数；h代表带宽，也可表示搜索半径，根据各类场所对周边pm
2.5
浓度的影响范围，搜索半径具体如表 2和表3所示。
[0082]
表2
[0083][0084]
表3
[0085][0086][0087]
最后，分析得到的核密度计算分布情况。如图8所示，列举了4 类较为典型的核密度分布情况。
[0088]
5、时空自相关性分析。
[0089]
莫兰指数可用于确定空间相邻区域的pm
2.5
浓度及其影响因素是否相似，其计算方式如下：
[0090][0091]
上式中，n为管控区域内的总网格数，xi和xj表示单元格i和单元格j的颗粒物浓度值，为所有空间单元格颗粒物浓度均值，ω
ij
为各单位之间的权重系数。
[0092]
另外，莫兰指数的显著性通常使用标准化z分数进行测试，其计算如下：
[0093][0094]
当z值为正且显著时，各pm
2.5
浓度之间存在正的空间自相关。相反，当z为负且显著，pm
2.5
浓度呈现负的空间自相关。
[0095]
利用莫兰指数和z值分析颗粒物浓度在0～23小时的空间自相关情况，观察其时空
异质性特点，强调引入时空地理加权回归模型的重要性。
[0096]
最终，时空异质分析结果如图9所示。主城区的莫兰指数与z 值在0至23小时内的波动较大。其中，莫兰指数的范围为0.05至0.45，呈双峰变化趋势，6时和16时的莫兰指数值为高峰值，13时和21时处于低谷值，而z值作为莫兰指数的显著性评估值，其变化趋势与莫兰指数相似，因此寻找各类poi、道路等级与pm
2.5
关系性对污染管控工作具有较大意义。
[0097]
6、时空地理加权回归分析
[0098]
为分析微尺度下各类典型poi和道路类型对当地颗粒物的影响程度，以研究区域内所有监测站的小时浓度均值作为区域整体浓度值，计算各监测点位的小时浓度均值与区域整体浓度值的差值，得到各监测点位的浓度变化幅值，该幅值能充分体现各网格的污染变化规律。若幅值为正值，表示该点位的污染排放强度明显；相反，则表示当地的污染净化能力较强，体现对pm
2.5
吸附、净化能力。
[0099]
gtwr模型(时空地理加权回归模型)作为gwr(地理加权回归模型)的一种补充方法，通过构建权重矩阵来解释参数估计的时空非平稳特性，该权重矩阵基于观测值和所有其他观测值之间的距离确定相关性，其计算方法如下：
[0100][0101]
式中，yi表示第i个监测点位的小时浓度均值与区域整体浓度值的差值；x
ik
为第i个监测点位归一化后的各类核密度值；(ui,vi,ti)表示第i个时空参数，分别对应网格的行数、列数以及0～23时；β0(ui,vi,ti) 表示第i个时空维度下的模型回归参数；βk(ui,vi,ti)表示第i个时空维度下的模型与第k个维度的回归参数；εi表示第i个时空维度的随机性误差。
[0102]
最后，计算得到的βk(ui,vi,ti)可表示各类输入因子的影响大小，此处可作为各类poi或道路等级对当地颗粒物的影响程度，并在时空维度为进行加权。当回归参数为正值时，表示该类型为颗粒物排放源，相反则表示该类型对颗粒物具有净化效果。
[0103]
利用主城区2020年冬季的监测站小时浓度变化幅值特征进行分析，各模型的具体分析结果如表4所示。gtwr模型在aicc和r2都优于普通最小二乘法和时空地理加权模型，并且r2高达0.84，表示输入的影响因子与pm
2.5
变化幅值存在较大相关型，能较好拟合各参数之间的关系。
[0104]
表4
[0105][0106]
如图10所示，各因子的影响程度在时间维度上都呈现出明显的变化规律。从道路等级角度分析，一级道路总体影响程度最小，在5 时达到的峰值(均小于4)；二级道路影响程度呈现凹形趋势，15时处于低谷值，23时处于峰值，高达7.9，表示该道路的pm
2.5
排放强度较高；三级道路的影响变化程度趋于平缓，总体处于1.6至4.4之间；而四级道路的呈现双谷变化趋势，变化程度最大，0至6时不断减小，6至11时呈现上升趋势，至15时又有所下降，然后剧烈上升至峰值(9.82)。总体而言，各级道路的影响程度主要与各道路的交通流量相关，
都呈现pm
2.5
排放特性，特别是二级和四级道路最为明显，在时间维度上变化较大。
[0107]
从四类典型poi角度出发，影响较大的为医疗服务和景区。医疗服务类型中主要分为各类医院、诊所和药店，该类场所的影响程度呈现双峰变化趋势，于9时和19时达到峰值(7.8和6.6)，且变化幅度较小，pm
2.5
的排放强度较均匀；景区类的影响程度均为负值，表示景区内的植被覆盖率对pm
2.5
有较大净化能力，夜间的净化水平较为稳定，而15时至16时的净化能力最小，该时段的变化与该场所内的市民活动强度相关。另外，餐饮服务和住宅服务(包括小区、酒店和饭店等)影响程度较小。
[0108]
根据pm
2.5
与各类潜在污染源的相关性分析结果，主城区的pm
2.5
污染存在较明显管控方法。
[0109]
交通源为主城区pm
2.5
的重点排放源。冬季的上下班高峰期，车辆发动机的不完全燃烧问题，导致颗粒物排放占比变高，而气温较低，又不利于空气扩散，在17时至23时的pm
2.5
聚集性污染问题越来越突出，因此针对主城区的pm
2.5
管控措施不仅仅局限于夜间，傍晚高峰期才为管控的“黄金期”。由于二级道路和四级道路相关性变化明显，而一级道路和三级道路的相关性相对稳定，因此可展开交通管控措施，控制二级道路和四级道路的车流量，引导部分车流量至一级道路和三级道路或引西湖景区西北部，如图8(a)，充分发挥景区对pm
2.5
明显的净化能力，降低燃油车尾气排放的影响程度。
[0110]
医疗服务场所也可作为pm
2.5
主要管控对象，由于医疗服务场所 (主要为各类大、中型医院)附近的市民活动水平较高，停车场容量和道路宽度有限，相关管控部门应保证该场所周围的交通畅通、合理的市民活动强度，并可考虑提高其周围绿化覆盖率。