CN110503348B - 一种基于位置匹配的个体空气污染暴露模拟测量方法 - Google Patents

Abstract

本发明公开了一种基于位置匹配的个体空气污染暴露模拟测量方法，首先利用污染监测站点和个体移动监测设备获取的六种空气污染物浓度数据进行污染精细制图；其次，对于未携带移动监测设备的个体，基于其手机APP定位进行邻近设备搜索；最后，基于个体模拟污染物浓度数据，结合个体暴露参数数据，计算个体空气污染模拟暴露剂量。本发明是针对当前空气污染暴露测量工作主要围绕大范围人群、携带移动监测设备个体的情况下，提出的一种针对未携带移动监测设备的综合性高、可参考性强的个体暴露精细模拟测量方法，为个体采取正确的空气污染暴露防护措施提供了参考。

Description

一种基于位置匹配的个体空气污染暴露模拟测量方法

技术领域

本发明涉及环境暴露评估领域，特别涉及利用污染制图方法与暴露行为活动模式对未携带便携式测量仪器个体的空气污染暴露量进行模拟的一种基于位置匹配的个体空气污染暴露模拟测量方法。

背景技术

随着我国社会经济的快速发展和城市化、工业化的飞速推进，以细颗粒物和臭氧为代表的空气污染问题日益突出，频发的空气污染浓度超标事件严重影响公众的生活与健康。截止2018年，全国338个城市仍有64.2％的城市环境空气污染超标，平均超标天数比例为20.7％。浓度超标的空气污染严重危害暴露人群的健康，中国地区因室外空气污染死亡人数为110.8万人，造成的寿命年损失为2177.9万年。

空气污染对人体健康的影响始于人与空气的接触(即暴露)。定量评估空气污染暴露剂量，是评价健康风险的前提，也是人群健康防护的基础。随着人们对环境与健康保护意识的增强，人们日常活动类型及活动范围的多样性发展，人群的平均暴露水平已不能真实的反映人们在不同微环境下的暴露特征，因此，个体空气污染暴露测量成为一个重要研究方向。

个体空气污染暴露测量是对个体暴露于空气中有害因子的强度、频率、时间、暴露的途径和方式进行测量、估算或预测的过程，其准确性一方面取决于先进的、高时间分辨率的个体暴露监测装置，另一方面依赖于适当的个体暴露评价技术方法。目前携带便携式暴露监测装置的个体暴露测量技术已较为成熟，但针对未携带仪器个体的精细暴露测量研究较少。

发明内容

针对现有技术存在的上述不足，本发明的目的是提供一种能够实现准确评估未携带测量仪器个体的空气污染暴露量的个体空气污染暴露模拟测量方法。

为了解决上述问题，本发明采用了以下的技术方案。

一种基于位置匹配的个体空气污染暴露模拟测量方法，其特征在于，包括以下步骤：

(一)利用污染监测站点或个体移动监测设备实时采集得到的空气污染物浓度数据，并结合区域内土地利用度数据、人口数据以及气象数据进行污染程度制图，污染程度制图包括以下子步骤：

(1)确定特征变量：选取建设用地X₁₁、耕地X₁₂、林地X₁₃、水体X₁₄、裸地X₁₅、人口密度X₂、高速道路X₃₁、主要道路X₃₂、次要道路X₃₃、其他道路X₃₄、高程X₄、温度X₅₁、相对湿度X₅₂、气压X₅₃、风速X₅₄和降水量X₅₅作为模型构建用预测变量；

(2)建立不同空间尺度预测变量：

以污染监测站点或个体移动监测设备所在位置的人口密度X₂、高程X₄、温度X₅₁、相对湿度X₅₂、气压X₅₃、风速X₅₄和降水量X₅₅的数值来表征影响因子，得出7个影响因子；

以污染监测站点或个体移动监测设备所在位置为圆心，并分别以r₁、r₂、r₃...r_n为半径建立系列缓冲区，统计各个缓冲区内建设用地X₁₁、耕地X₁₂、林地X₁₃、水体X₁₄、裸地X₁₅、高速道路X₃₁、主要道路X₃₂、次要道路X₃₃、其他道路X₃₄的数值，得出9n个影响因子；

(3)建立多元回归模型：将步骤(2)得到的各影响因子与空气污染物最近一周的浓度均值进行双变量相关分析，筛选出与污染物浓度显著相关的影响因子；将筛选出的与污染物浓度显著相关影响因子与污染物浓度值进行逐步线性回归处理：

C＝a₀+a₁X₁+a₂X₂+a₃X₃+...+a_mX_m+μ (1)

式中，C为污染物浓度值，X₁、X₂、X₃…X_m为筛选出的与污染物浓度显著相关的影响因子，a为未知参数，μ为随机误差项，未知参数a和随机误差项μ使用最小二乘法确定；

(4)污染物浓度模拟：在研究区内均匀设置网格点，根据公式(1)计算出各个网格点处污染物浓度预测值；采用克里格插值法生成研究区内大气污染物浓度值的空间模拟分布图；

(二)对于未持移动监测设备的个体，基于其手机APP定位进行邻近设备搜索，依据其邻近设备的污染物浓度确定个体所在位置各污染物项目浓度值；

(三)根据步骤(一)和步骤(二)的结果，结合个体体重、活动时间以及活动期间呼吸速率等暴露参数，计算个体空气污染暴露模拟剂量，其中,时间分辨率为D，D≤1h。

进一步的，所述步骤(二)包括：

(1)以未持个体监测设备的个体为圆心，设定半径为200m的搜索圆；搜索圆的外接正方形的经纬度范围内的邻近设备；

基于公式——

其中，haversin运算与正弦运算、余弦运算的关系为：

haversine(θ)＝sin²(θ/2)＝(1-cos(θ))/2 (3)

得出：

外接正方形纬度范围的

所对应的弧度值为：

外接正方形经度范围的

所对应的弧度值为：

将Δλ,

转为角度，分别用dlat,dlon表示：

其中：

R为地球半径；

d为搜索范围半径，取200m；

为地球上任意两点的坐标；

则搜索圆的外接正方形四个顶点的经纬度坐标分别为：

左上顶点经纬度left_top:(lat+dlat,lon-dlon)

右上顶点经纬度right_top:(lat+dlat,lon+dlon)

左下顶点left_bottom:(lat-dlat,lon-dlon)

右下顶点right_bottom:(lat-dlat,lon+dlon)

其中，(lat,lon)为未持个体监测设备的个体所在位置的经纬度坐标；

(2)根据移动监测设备采集的志愿者所在的经纬度坐标M(lat_i,lon_i)，判断志愿者是否处于步骤(1)得出的外接正方形范围内：

当M(lat_i,lon_i)同时满足(lat-dlat)≤lat_i≤(lat+dlat)和(lon-dlon)≤lon_i≤(lon+dlon)时，则志愿者处于搜索圆外接正方形范围内，该志愿者记作第i个志愿者M_i；i为大于1的自然数；

(3)计算未持移动监测设备的个体所在位置到邻近设备搜索范围内各志愿者所在位置的欧氏距离，计算公式如下：

式中：d_i代表未持移动监测设备的个体所在位置到第i个志愿者M_i所在位置的欧氏距离；(x₀,y₀)是未持移动监测设备的个体所在位置的坐标，(x_i,y_i)是第i个志愿者M_i所在位置的坐标；

当距离d_i≤200m时，则第i个志愿者M_i处搜索圆范围内，并记录志愿者编号i与距离d_i；反之，第i个志愿者M_i不处于搜索圆范围内。

(4)当所述搜索圆范围内存在志愿者时，未持移动监测设备的个体所在位置污染物模拟浓度取该个体距搜索范围内所有志愿者距离的加权值：

式中：C_p代表未持移动监测设备的个体所在位置的污染物p的模拟浓度值；

λ_i为未持移动监测设备的个体到第i个志愿者M_i距离的权重系数：

其中，d_i为未持移动监测设备的个体所在位置到第i个志愿者M_i所在位置的欧氏距离；

C_{p_i}为搜索范围内第i个志愿者的移动监测设备的污染物p的浓度；

当搜索范围内不存在志愿者时，个体所在位置的污染物模拟浓度为污染精细制图结果对应网格的浓度值。

相比现有技术，本发明具有如下优点：

本发明是一种用于计算个体空气污染暴露模拟剂量的一种方法，是个体空气污染暴露风险评价的前提。相对于已有的暴露模拟测量方法只针对携带移动监测设备的个体，本方法通过搜索个体邻近设备进行个体所在位置污染物浓度模拟，结合个体暴露参数数据进行个体暴露模拟测量，克服了实地监测的困难，可在暴露风险评价工作中节省大量财力。同时，本方法不仅适用于携带移动监测设备的个体，更为未携带移动监测设备的个体提供了个暴露风险评价的可靠手段，扩大了可进行暴露风险评价的群体范围。

附图说明

图1为本发明的示意图；

图2为本发明搜索圆的外接正方形的经纬度确定方法示意图。

具体实施方式

下面结合附图和实例，对本发明做进一步详细说明。

本发明基于暴露时空行为活动模式与污染制图方法，运用地理信息系统和空间分析等手段，提供一种能够实现准确评估未携带测量仪器个体的空气污染暴露量的个体空气污染暴露模拟测量方法。

首先利用污染监测站点和个体移动监测设备获取的六种空气污染物(SO₂、NO₂、CO、O₃、PM₁₀、PM_2.5)浓度数据进行污染精细制图(空间分辨率≤200m)；其次，对于未携带移动监测设备的个体，基于其手机APP定位进行邻近设备搜索(搜索半径为200m)；最后，基于个体模拟污染物浓度数据，结合个体暴露参数数据，计算个体空气污染模拟暴露剂量(时间分辨率≤1h)。这是当前空气污染暴露测量工作主要围绕大范围人群、携带移动监测设备个体的情况下，发明的一种针对未携带移动监测设备的综合性高、可参考性强的个体暴露精细模拟测量方法，为个体采取正确的空气污染暴露防护措施提供参考。

如图1所示，包括以下步骤：

(一)利用污染监测站点或个体移动监测设备实时采集得到的空气污染物浓度数据，并结合区域内土地利用度数据、人口数据以及气象数据进行污染程度制图，污染程度制图包括以下子步骤：

(1)确定特征变量：

选取土地利用(建设用地X₁₁、耕地X₁₂、林地X₁₃、水体X₁₄、裸地X₁₅)、人口密度X₂、道路交通(高速道路X₃₁、主要道路X₃₂、次要道路X₃₃、其他道路X₃₄)、高程X₄、气象(温度X₅₁、相对湿度X₅₂、气压X₅₃、风速X₅₄和降水量X₅₅)共16类作为模型构建的预测变量；

(2)建立不同空间尺度预测变量：

例如，以污染监测站点或个体移动监测设备所在位置为中心建立系列缓冲区(0.1km,0.2km,0.3km,0.4km,0.5km,0.6km,0.7km,0.8km,0.9km,1km)，统计不同半径缓冲区内影响因子数值(如各监测点1km缓冲区内建设用地面积)。由于人口密度、高程和气象要素均以监测点的数值来表征，不涉及缓冲区分析，因此，针对X₁₁、X₁₂、X₁₃、X₁₄、X₁₅、X₃₁、X₃₂、X₃₃、X₃₄)，缓冲区分析后共得到90(9×10)个影响因子；再加上X₂、X₄、X₅₁、X₅₂、X₅₃、X₅₄、X₅₅，共计97(90+7)个影响因子。

(3)建立多元回归模型：

将单个影响因子与各污染物(SO₂、NO₂、CO、O₃、PM₁₀、PM_2.5)近一周的浓度均值进行双变量相关分析，筛选出与污染物浓度显著相关(置信度p＜0.05)的影响因子；为了降低后续多元线性回归中同类变量(如建设用地)不同缓冲半径下影响因子间的共线性，先选出同类预测变量中与污染物浓度相关性最高的影响因子(如1km缓冲区内建设用地面积)，然后去除与该因子相关性较高(相关系数r＞0.6)的其他同类因子(如0.8km缓冲区内建设用地面积)；最后将剩下的所有影响因子与污染物浓度值进行逐步线性回归，得到多元线性回归方程，即研究区的LUR模型，如公式(1)所示。

C＝a₀+a₁X₁+a₂X₂+a₃X₃+...+a_mX_m+μ (1)

式中,C为污染物浓度值，X为最终进入模型的地理要素特征变量，a为未知参数，μ为随机误差项，a和μ使用最小二乘法确定。

此外，本发明还提供了模型验证方法对模型精度进行评估。污染监测站点或个体移动监测设备共n个，将污染监测站点或个体移动监测设备分为90％的实验集和10％的验证集，先用实验集中污染监测站点或个体移动监测设备的各污染物浓度与LUR模型中的预测变量进行回归，利用回归方程预测验证集的各污染物浓度，并将该点的预测值和监测值比较；该过程共重复n次，得到区域LUR模型的模拟精度和均方根误差。

(4)污染物浓度模拟：

在研究区均匀布点(2km×2km)，并计算网格点相应地理变量相关数值，标准化处理后，代入公式(1)，计算出网格点处各污染物浓度预测值。采用克里格插值法生成研究区大气污染物浓度值的空间模拟分布图。

(二)对于未持移动监测设备的个体，基于其手机APP定位进行邻近设备搜索，依据其邻近设备的污染物浓度确定个体所在位置各污染物项目浓度值；

(1)对未持移动监测设备的个体进行邻近设备搜索的范围确定：

如图2所示，以未持个体监测设备的个体为圆心，设定半径为200m的搜索圆；首先搜索圆的外接正方形的经纬度范围内的邻近设备；

(2)根据移动监测设备采集的志愿者所在的经纬度坐标(lat_i,lng_i)，判断志愿者是否处于步骤(1)得出的外接正方形范围内：

基于公式——

其中，haversin运算与正弦运算、余弦运算的关系为：

haversine(θ)＝sin²(θ/2)＝(1-cos(θ))/2 (3)

得出：

如图2所示，首先求东西两侧的边界范围，将O、F两点的坐标带入公式(2)，即令

则外接正方形纬度范围的一半所对应的弧度值为：

然后求南北两侧的边界范围，将O、E两点的坐标带入公式(2)，即令λ₁＝λ₂，外接正方形纬度范围的

所对应的弧度值为：

外接正方形经度范围的

所对应的弧度值为：

将Δλ,

转为角度，分别用dlat,dlon表示：

其中：

R为地球半径；

d为搜索范围半径，取200m；

为地球上任意两点的坐标；

则搜索圆的外接正方形四个顶点的经纬度坐标分别为：

左上顶点A经纬度left_top:(lat+dlat,lon-dlon)

右上顶点B经纬度right_top:(lat+dlat,lon+dlon)

左下顶点Dleft_bottom:(lat-dlat,lon-dlon)

右下顶点Cright_bottom:(lat-dlat,lon+dlon)

其中，(lat,lon)为未持个体监测设备的个体所在位置的经纬度坐标；

(2)根据移动监测设备采集的志愿者所在的经纬度坐标M(lat_i,lon_i)，判断志愿者是否处于步骤(1)得出的外接正方形范围内：

当M(lat_i,lon_i)同时满足(lat-dlat)≤lat_i≤(lat+dlat)和(lon-dlon)≤lon_i≤(lon+dlon)时，则志愿者处于搜索圆外接正方形范围内，该志愿者记作M_i；i为志愿者设备编号，为大于1的自然数；

也即是说这里标记的志愿者M_i可能是M₁、M₂、M₃直至若干个志愿者，这些志愿者均处于搜索圆外接正方形范围内。

(3)计算未持移动监测设备的个体所在位置到邻近设备搜索圆外接正方形范围内各志愿者所在位置的欧氏距离，计算公式如下：

这里即是计算步骤(2)中判断找出的处于外接正方形范围内，并被标记为志愿者M_i的若干个志愿者与未持移动监测设备的个体之间的欧式距离。以实现进一步的志愿者筛选。

式中：d_i代表未持移动监测设备的个体所在位置到志愿者M_i所在位置的欧氏距离；(x₀,y₀)是未持移动监测设备的个体所在位置的坐标，(x_i,y_i)是志愿者M_i所在位置的坐标；

当距离d_i≤200m时，则志愿者M_i处搜索圆范围内，并记录志愿者编号i与距离d_i；反之，志愿者M_i不处于搜索圆范围内。也即是说这里最后找出的志愿者设备所检测的污染物浓度数据将用于未持移动监测设备的个体所在位置污染物模拟浓度计算。

(4)当所述搜索圆范围内存在志愿者时，未持移动监测设备的个体所在位置污染物模拟浓度取该个体距搜索范围内所有志愿者距离的加权值：

式中：C_p代表未持移动监测设备的个体所在位置的污染物p的模拟浓度值；

λ_i为未持移动监测设备的个体到第i个志愿者M_i距离的权重系数：

其中，d_i为未持移动监测设备的个体所在位置到第i个志愿者M_i所在位置的欧氏距离；

C_{p_i}为搜索范围内第i个志愿者的移动监测设备的污染物p的浓度；

当搜索范围内不存在志愿者时，个体所在位置的污染物模拟浓度为污染精细制图结果对应网格的浓度值。

现有邻近设备搜索方法是将计算个体所在位置的坐标与所有移动监测设备的距离，并与搜索半径进行比较，以此来筛选处于个体所在位置的200m范围内的移动监测设备，这种搜索方法时间复杂度高，搜索效率低。本发明提供的搜索方法首先确定个体所在位置的搜索圆的外接正方形的顶点坐标，然后通过各移动监测设备的坐标与其顶点坐标进行比较，确定处于搜索圆外接正方形范围内的移动监测设备，最后再计算处于外界正方形范围内的移动监测设备与个体所在位置的距离并与搜索半径进行比较。因此减少了需进行距离计算的设备数量，降低了时间复杂度，使得搜索效率得到提高。

(三)根据步骤(一)和步骤(二)的结果，结合个体体重、活动时间以及活动期间呼吸速率等暴露参数，计算个体空气污染暴露模拟剂量，其中,时间分辨率≤1h。

具体的，暴露剂量是用于描述人体经呼吸道接触空气污染物的量，主要计算各污染物的个体单位体重时均暴露剂量值，如下式所示。

式中：D_p为个体暴露于污染物P的单位体重时均暴露剂量；C_p为污染物P的浓度值；IR_i为个体呼吸速率；ET_i为各运动状态持续时间；BW为个体体重。六种空气污染物(SO₂、NO₂、CO、O₃、PM₁₀、PM_2.5)浓度值为步骤4中计算的个体所在位置的污染物模拟浓度值；个体体重、个体运动持续时间通过个体基本信息表获取；个体呼吸速率来自于《中国人群暴露参数手册》，通过性别、年龄、运动状态等条件进行查询。

最后需要说明的是，以上实施实例仅用以说明本发明的技术方案而非限制，尽管参照较佳实例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或等同替换，而不脱离本发明技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (1)

1.一种基于位置匹配的个体空气污染暴露模拟测量方法，其特征在于，包括以下步骤：

(一)利用污染监测站点或个体移动监测设备实时采集得到的空气污染物浓度数据，并结合区域内土地利用度数据、人口数据以及气象数据进行污染程度制图，污染程度制图包括以下子步骤：

(1)确定特征变量：选取建设用地X₁₁、耕地X₁₂、林地X₁₃、水体X₁₄、裸地X₁₅、人口密度X₂、高速道路X₃₁、主要道路X₃₂、次要道路X₃₃、其他道路X₃₄、高程X₄、温度X₅₁、相对湿度X₅₂、气压X₅₃、风速X₅₄和降水量X₅₅作为模型构建用预测变量；

(2)建立不同空间尺度预测变量：

以污染监测站点或个体移动监测设备所在位置的人口密度X₂、高程X₄、温度X₅₁、相对湿度X₅₂、气压X₅₃、风速X₅₄和降水量X₅₅的数值来表征影响因子，得出7个影响因子；

以污染监测站点或个体移动监测设备所在位置为圆心，并分别以r₁、r₂、r₃...r_n为半径建立系列缓冲区，统计各个缓冲区内建设用地X₁₁、耕地X₁₂、林地X₁₃、水体X₁₄、裸地X₁₅、高速道路X₃₁、主要道路X₃₂、次要道路X₃₃、其他道路X₃₄的数值，得出9n个影响因子；

(3)建立多元回归模型：将步骤(2)得到的各影响因子与空气污染物最近一周的浓度均值进行双变量相关分析，筛选出与污染物浓度显著相关的影响因子；将筛选出的与污染物浓度显著相关影响因子与污染物浓度值进行逐步线性回归处理：

式中，C为污染物浓度值，X₁、X₂、X₃…X_m为筛选出的与污染物浓度显著相关的影响因子a_i为未知参数，a_i为a₀、a₁、a₂、a₃...a_m；μ为随机误差项，未知参数a_i和随机误差项μ使用最小二乘法确定；

(4)污染物浓度模拟：在研究区内均匀设置网格点，根据公式(1)计算出各个网格点处污染物浓度预测值；采用克里格插值法生成研究区内大气污染物浓度值的空间模拟分布图；

(二)对于未持移动监测设备的个体，基于其手机APP定位进行邻近设备搜索，依据其邻近设备的污染物浓度确定个体所在位置各污染物项目浓度值；

(1)以未持个体监测设备的个体为圆心，设定半径为200m的搜索圆；搜索圆的外接正方形的经纬度范围内的邻近设备；

基于公式——

其中，haversin运算与正弦运算、余弦运算的关系为：

haversine(θ)＝sin²(θ/2)＝(1-cos(θ))/2 (3)

得出：

外接正方形纬度范围的

所对应的弧度值为：

外接正方形经度范围的

所对应的弧度值为：

将Δλ,

转为角度，分别用dlat,dlon表示：

其中：

R为地球半径；

d为搜索范围半径，取200m；

为地球上任意两点的坐标；

则搜索圆的外接正方形四个顶点的经纬度坐标分别为：

左上顶点经纬度left_top:(lat+dlat,lon-dlon)

右上顶点经纬度right_top:(lat+dlat,lon+dlon)

左下顶点left_bottom:(lat-dlat,lon-dlon)

右下顶点right_bottom:(lat-dlat,lon+dlon)

其中，(lat,lon)为未持个体监测设备的个体所在位置的经纬度坐标；

(2)根据移动监测设备采集的志愿者所在的经纬度坐标M(lat_i,lon_i)，判断志愿者是否处于步骤(1)得出的外接正方形范围内：

当M(lat_i,lon_i)同时满足(lat-dlat)≤lat_i≤(lat+dlat)和(lon-dlon)≤lon_i≤(lon+dlon)时，则志愿者处于搜索圆外接正方形范围内，该志愿者记作第i个志愿者M_i；i为大于1的自然数；

(3)计算未持移动监测设备的个体所在位置到邻近设备搜索范围内各志愿者所在位置的欧氏距离，计算公式如下：

式中：d_i代表未持移动监测设备的个体所在位置到第i个志愿者M_i所在位置的欧氏距离；(x₀,y₀)是未持移动监测设备的个体所在位置的坐标，(x_i,y_i)是第i个志愿者M_i所在位置的坐标；

当距离d_i≤200m时，则第i个志愿者M_i处搜索圆范围内，并记录志愿者编号i与距离d_i；反之，第i个志愿者M_i不处于搜索圆范围内；

(4)当所述搜索圆范围内存在志愿者时，未持移动监测设备的个体所在位置污染物模拟浓度取该个体距搜索范围内所有志愿者距离的加权值：

式中：C_p代表未持移动监测设备的个体所在位置的污染物p的模拟浓度值；

λ_i为未持移动监测设备的个体到第i个志愿者M_i距离的权重系数：

其中，d_i为未持移动监测设备的个体所在位置到第i个志愿者M_i所在位置的欧氏距离；

C_{p_i}为搜索范围内第i个志愿者的移动监测设备的污染物p的浓度；

当搜索范围内不存在志愿者时，个体所在位置的污染物模拟浓度为污染精细制图结果对应网格的浓度值；

(三)根据步骤(一)和步骤(二)的结果，结合个体体重、活动时间以及活动期间呼吸速率暴露参数，计算个体空气污染暴露模拟剂量。

Images