气象因素以多种交互作用的方式对人类的健康及生存造成一系列影响^{［1, 2, 3, 4］}。由于多种气象因素对人群健康均有影响，确定气象风险阈值，提前预警不良气象条件的健康风险，从而采取及时有效的应对措施，可降低不良气象事件对人群健康危害^{［5, 6］}。

近年来在国内开展了气象健康风险预警的相关研究^{［7, 8］}。大部分仅考虑气温，或仅对高温热浪进行健康风险预警^{［9, 10］}，很少研究综合考虑多种气象因素。然而，气象因素对人体健康的影响不只体现在气温，极端气象事件也不局限于高温热浪，为了综合评估气象因素的健康风险和预警，应考虑多种气象因素的暴露，选取综合多种气象因素指标作为气象健康风险的暴露指标。体感温度（apparent temperature，AT）是基于人体舒适度而开发的综合多种气象因素的生物气象指标^{［11, 12, 13］}，有研究表明该指标更能反映不良气象因素对人群健康的综合影响^［14］，该指标既能反映人体在不同气象条件下的舒适度，也能评价气象条件对人体健康的影响，是建立气象健康风险预警较为理想的指标。

目前国内外气象健康风险预警阈值的确定没有统一标准^{［15, 16］}，大部分研究以气象健康风险的百分比确定风险预警阈值，缺乏科学性。有研究采用分类回归决策树的方法分析与传染病发生相关的自变量预警阈值^{［17, 18, 19］}，但尚未有研究将该方法应用到气象健康风险预警阈值确定。

在既往相关的研究中，大部分以死亡数作为健康效应终点。然而死亡数这一指标仅考虑了死亡的数量，忽略了死亡时的年龄，而不良气象条件往往对老年人死亡的影响更大^{［2，20, 21］}，因此采用死亡数作为效应终点可能会高估不良气象条件对人群健康的影响。寿命损失年（years of life lost，YLL）同时考虑了死亡数以及不同死亡年龄的权重，能更准确地反映气象因素对人群健康的影响^［22］。

在本研究中，我们首先分析AT与YLL率的暴露-反应关系，基于该暴露-反应关系采用回归决策树模型确定各风险等级的预警阈值，为开展气象健康风险预警服务提供研究基础。

1. 研究地点：本研究选取全国364个区（县）作为研究点。选取标准是区（县）总人口数>20万，或年死亡率>4‰。地理分布见图1。

点击查看大图

图1

364个研究区（县）分布情况

点击查看大图

图1

364个研究区（县）分布情况

2. 数据收集

（1）死亡数据：云南、广东、湖南、浙江和吉林省研究区（县）的死亡数据从各省CDC获取，研究时间为2013年1月1日至2017年12月31日，其他省份的研究区（县）死亡数据从国家CDC获取，研究时间为2006年1月1日至2012年12月31日。获取的死亡数据包含年龄、性别、根本死因、死亡日期、国际疾病分类（ICD）编码等信息。根据ICD第10版（ICD-10）对死因分类，包括非意外死亡（A00~R99），循环系统疾病（I00~I99，CVD）和呼吸系统疾病（J00~J98，RESP）。

（2）气象数据：从中国气象科学数据共享服务平台（http：//data.cma.cn/）获取全国698个气象站点2006-2017年每日气象资料，包括平均温度（℃）、相对湿度（%）、风速（m/s）等。由于气象站点并未覆盖所有的研究区（县），所以运用气象空间插值方法Anusplin对每日平均温度与相对湿度进行插值。Anusplin的原理是基于薄板样条函数理论，引入经纬度作为自变量，海拔作为协变量进行气象要素空间插值^［23］。本研究利用全国698个气象站点的数据，插值形成分辨率为0.01°×0.01°的全国日均温度的栅格数据。通过10折交叉方法进行验证，日均温度的插值模型的R²=0.96，均方根误差（root mean square error，RMSE）=2.37 ℃。采用同样的方法获取每日相对湿度栅格数据，相对湿度插值模型的R²=0.81，RMSE=7.7%。随后从每日的全国栅格数据提取每个研究区（县）的每日平均温度和每日相对湿度。研究区（县）每日风速数据来源于距离最近的气象监测站。

（3）大气污染物数据：由于难以获得其他大气污染物在整个研究周期的数据，选取每日PM₁₀平均浓度作为分析AT与YLL关系的大气污染物控制变量。从全国城市空气质量实时发布平台获取2006-2017年国家污染物监测站PM₁₀平均浓度数据。由于大气污染物监测点并未覆盖所有的研究地点，运用随机森林模型进行预测每个研究区（县）每日PM₁₀平均浓度。利用每个大气污染物监测站点的经纬度、海拔、每日平均温度和相对湿度，以及每个监测站周围1 300 m缓冲区内的人口密度、道路密度、土地利用类型以及人均GDP和每日PM₁₀数据构建随机森林模型^{［24, 25］}。同时在模型中加入平滑时间函数以控制污染物浓度的长期趋势和季节变化。采用与气象数据插值模型相似的10折交叉验证方法，结果表明，该模型R²=0.78，RMSE=13.2 μg/m³。

（4）人口与地理信息数据：从南安普顿大学地理数据研究所获取2015年全国人口密度栅格数据 （https：//www.worldpop.org/）。从中国科学院资源环境数据云平台（http：//www.resdc.cn）获取道路密度栅格数据、土地利用类型栅格数据、全国县级矢量地图等地理信息数据以及全国GDP栅格数据。

3. 相关指标计算

（1）YLL率：利用各省2010年普查数据中人口信息和死亡数信息，按照WHO提供的方法计算各省寿命表。匹配死亡个例的性别与年龄，计算个体YLL，汇总每日总YLL。由于各研究点YLL值与人口数相关，为了合并各区（县）的效应并进行比较，将每日总YLL除以各区（县）的人口数，得到YLL率（每10万人口）作为健康结局指标。

（2）体感温度：

式中，T为气温（℃），e为水汽压（hPa），v为风速（m/s），RH为相对湿度（%）。

4. 统计学分析：采用R 3.6.0软件进行统计分析。采用dlnm R程序包构建DLNM模型，采用mvmeta R程序包进行多变量Meta分析，采用rpart R程序包构建回归决策树模型。

（1）AT与YLL率暴露反应关系：本研究利用2006-2016年各区（县）时间序列数据，采用两阶段分析方法分析AT与YLL率的关系。第一阶段，采用分布滞后非线性模型（distributed lag non-linear model，DLNM）连接高斯函数分析单个区（县）AT与YLL率暴露反应关系，模型公式：

式中，YLL_t为第td的YLL率；cb（AT，lag）为AT与滞后天数的“交叉基”函数，其中包含一个带有3个内部节点（AT范围的第10、50和90百分位数）的二次样条函数和一个包含滞后效应的三次样条函数^［26］，由于低温效应往往持续3周左右，为获取总体的气温滞后效应，滞后天数设定为21 d^{［3, 4］}；ns（time，6）为时间变量的自然立方样条函数，用于控制长期趋势和季节趋势，根据赤池信息量准则值（Akaike information criterion，AIC）最小化原则确定每年的df设为6/年；β为回归系数，DOW为控制“星期几效应”的分类变量。

第二阶段，采用多变量Meta分析合并所有区（县）的结果，分析AT与YLL率的累积暴露-反应关系^{［27, 28］}。根据累积暴露-反应曲线，估计最低YLL率对应的AT（minimum YLL apparent temperature，MYT）。以MYT为分界划分冷效应和热效应，AT对应的效应值称为归因YLL。对不同年龄组和死因人群进行分层分析。

（2）体感温度健康风险预警阈值确定：本研究构建以AT为自变量，以每日归因YLL率为因变量的回归树模型，分别确定冷效应和热效应的健康风险预警阈值。具体而言，回归树模型以每日归因YLL的均方误差（mean square error，MSE）最小化为原则，将AT划分为m个范围，并输出各个AT范围对应的日均归因YLL，输出预测值大则代表该AT范围的健康危害效应大，即风险等级高，划分AT的分节点为对应风险等级的阈值。为兼顾最终叶节点的数量以及模型误差的大小，通过预数据分析，本研究将回归树模型最终叶节点数确定为4，即将冷效应和热效应划分为4个风险等级，分别定义为高风险、中风险、低风险和舒适范围。

（3）体感温度健康风险预警阈值的验证：计算各个健康风险等级的超额死亡数评估风险等级划分的合理性。首先利用2006-2017年全国数据，采用两阶段分析方法拟合AT和非意外死亡数的暴露反应关系。在该暴露-反应关系中获取最低死亡风险体感温度（minimum mortality apparent temperature，MMT）以及AT导致的死亡风险（相对危险值RR）。随后，利用MMT对应的日均死亡数作为基线数据，结合死亡风险计算每个健康风险等级的日均超额死亡数。计算公式：

式中，E_d为日均超额死亡例数，N_MMT为MMT对应的日均死亡例数，RR_t为AT对应的死亡风险值。计算2006-2016年（即参与回归树模型的期间数据）的日均超额死亡数进行内部验证，计算2017年（即不参与回归树模型的数据）的日均超额死亡数进行外部验证，若日均超额死亡数随风险等级的升高而升高，则该健康风险预警阈值达到预警效果。

（4）敏感性分析：通过改变AT与YLL率DLNM模型滞后天数（14 d、21 d、28 d），对模型进行敏感性分析。

1. 基本情况：本研究总人群平均每日YLL率为23.6人年/10万，≥65岁人群平均每日YLL率为107.6人年/10万，0~64岁人群平均每日YLL率为14.2人年/10万。研究区（县）每日平均体感温度为15.7 ℃。见表1。

点击查看表格

表1

2006-2017年每日YLL率、气象条件与大气污染物基本情况（x±s）

表1

2006-2017年每日YLL率、气象条件与大气污染物基本情况（x±s）

变量	非意外死亡	循环系统疾病	呼吸系统疾病
YLL率（人年/10万）
总人群	23.6±15.9	8.1±8.2	2.4±4.1
≥65岁	107.6±74.4	49.9±50.5	18.2±45.5
0~64岁	14.2±14.7	3.7±6.8	0.7±2.6
气象因素
平均气温（℃）	16.0±9.8
相对湿度（%）	73.0±13.9
风速（m/s）	2.1±1.2
体感温度（℃）	15.7±12.6
PM10（μg/m3）	81.1±41.4

2. 暴露反应关系：AT与YLL率的暴露反应关系类似“U”形。在不同的死因人群中，AT对≥65岁人群的影响均远大于0~64岁人群，而循环系统疾病人群对AT敏感性也较呼吸系统疾病高（图2）。

点击查看大图

图2

不同人群体感温度与寿命损失年率累积暴露-反应关系

点击查看大图

图2

不同人群体感温度与寿命损失年率累积暴露-反应关系

3. 体感温度健康风险预警阈值：回归决策树模型经交叉验证得到预测误差均<0.01。≥65岁组的冷效应风险预警阈值均比0~64岁组高。在非意外死因和呼吸系统疾病死因中，0~64岁组热效应风险预警阈值比≥65岁人群稍低。见表2。

点击查看表格

表2

不同人群体感温度健康风险等级阈值及日均归因寿命损失年（YLL）

表2

不同人群体感温度健康风险等级阈值及日均归因寿命损失年（YLL）

效应等级	非意外死亡		循环系统疾病		呼吸系统疾病
	体感温度阈值（℃）	日均归因YLL	体感温度阈值（℃）	日均归因YLL	体感温度阈值（℃）	日均归因YLL
0~64岁
冷效应
	高风险	<-14.5	4.04	<-17.5	1.61	<-16.9	0.37
	中风险	<-3.7	2.12	<-7.5	0.88	<-7.6	0.22
	低风险	<6.7	0.89	<2.5	0.33	<0.1	0.09
	舒适	<26.9	0.12	<32.3	0.06	<9.6	0.02
热效应
	舒适	≥26.9	0.02	≥32.3	0.02	≥9.6	0.01
	低风险	≥33.6	0.22	≥34.6	0.10	≥16.7	0.05
	中风险	≥35.2	0.54	≥36.0	0.21	≥21.3	0.10
	高风险	≥37.2	1.16	≥37.7	0.43	≥29.5	0.14
≥65岁
冷效应
	高风险	<-9.1	57.35	<-9.9	31.24	<-11.6	9.63
	中风险	<4.2	27.44	<3.2	14.87	<-1.1	4.90
	低风险	<14.6	14.17	<13.7	7.70	<7.2	2.24
	舒适	<32.2	3.98	<32.2	2.14	<32.1	0.42
热效应
	舒适	≥32.2	0.50	≥32.2	0.29	≥32.1	0.12
	低风险	≥34.5	2.81	≥34.5	1.57	≥34.4	0.64
	中风险	≥35.9	6.28	≥35.9	3.49	≥35.8	1.41
	高风险	≥37.6	12.72	≥37.6	7.04	≥37.5	2.85

4. 体感温度健康风险预警阈值验证：内部数据和外部数据均显示，日均超额死亡数随着风险等级的增大而上升，内部及外部数据验证结果都符合预期设想（图3）。

点击查看大图

图3

不同人群各体感温度健康风险等级日均超额死亡数

点击查看大图

图3

不同人群各体感温度健康风险等级日均超额死亡数

5. 敏感性分析：改变DLNM模型中不同的滞后天数，对各健康风险等级的预警阈值的影响较小（表3），结果相对稳定。

点击查看表格

表3

不同滞后天数下体感温度健康风险预警阈值（℃）

表3

不同滞后天数下体感温度健康风险预警阈值（℃）

效应等级	滞后14 d	滞后21 d	滞后28 d
冷效应
高风险	<-10.0	<-12.0	<-12.3
中风险	<1.0	<-0.1	<-0.8
低风险	<9.8	<10.3	<9.2
热效应
低风险	≥34.2	≥34.3	≥33.7
中风险	≥35.7	≥35.8	≥35.3
高风险	≥37.5	≥37.5	≥37.2

在全球气候变化的大背景下，各气象因素相互作用对人群健康产生影响，为减少不良气象条件的健康危害，提前预警不良气象事件，做好应对措施，能有效降低气象健康风险。本研究基于AT与YLL的暴露-反应关系，采用回归树模型确定AT的健康风险预警阈值，为制定我国气象健康风险预警提供理论支撑。

为减少不良气象条件对人群健康的影响，国内不少学者开发了气象健康风险等级。在既往的研究中，仅将温度作为暴露指标确定气象风险阈值，如南京市、深圳市的高温健康风险早期预测系统的确定都只是考虑了气温^{［9, 10］}，将湿度作为混杂因素在模型中控制，而在本研究中，考虑不同气象因素的综合效应，将AT作为气象健康风险预警的暴露指标。AT是人体在某个环境感受下的冷暖程度转化为同一参考湿度水平和同一参考风速水平下的温度，该指标是基于热物理学模型，综合了湿度、风速对温度的修饰效应而得到的^{［11, 12, 13］}。曾韦霖等^［29］研究发现如果将不同的单一气象指标以独立变量引入模型，会忽略各个气象因素之间的交互作用，而降低模型的准确性。因此，将基于热物理模型以及人体舒适感，综合了多种单一气象指标的暴露指标纳入模型，能更准确地评价气象因素对人群健康的影响，有利于提高气象健康风险预警的准确性。

在既往气象健康风险预警的相关研究中，大部分研究采用回归方程法或广义加性模型评估气象因素的健康效应，忽略了气象因素的滞后效应^［8，30］。本研究采用DLNM模型拟合AT和YLL的暴露反应关系可以控制时间变量等混杂因素以及气象因素滞后效应的影响。本研究发现，过高或过低AT均会使YLL率增高，而且低AT的效应高于高AT的效应，这与以往关于气温与死亡的相关研究结果相似^{［3, 4，31, 32］}。

既往的相关研究中，在气象风险等级的划分上，大多数研究根据日最高温度与持续时间相结合或气象健康风险的百分比确定气象指标的健康风险阈值^{［7, 8, 9，33］}。如黄照等^［8］基于气温与死亡数据的暴露反应关系，将由气温导致的死亡效应（RR值）平均分为3个等级，进而确定各个气象健康风险等级的气温阈值。本研究采用回归决策树模型以AT的效应值（归因YLL）的MSE最小化为原则，划分AT范围，确定冷效应和热效应的风险预警阈值，利用该机器学习方法以实际数据为导向确定气象健康风险预警阈值，比起直接以人为确定的百分位数更客观并具有可靠性。

本研究针对不同人群确定气象健康风险预警阈值。由于不同人群对气象因素的敏感性和适应性不同，对较为脆弱的人群的预警标准应较低。比如：老年人冷效应的风险预警阈值较年轻人低，且相差5 ℃以上，这是因为相比年轻人，老年人对寒冷的脆弱性更大。然而，老年人热效应的风险预警阈值稍高于年轻人，这可能是由于年轻人高温暴露较老年人多，这提示年轻人在出现热效应风险预警时要做好相关防暑措施，尤其是户外工作者尽量避免在高温时段进行户外工作。针对不同人群提供针对性的气象健康风险预警预报，有利于对脆弱性人群采取更加及时的保护性措施，有利于提高预警系统的特异性，同时也节省了风险预警系统的运行成本。

本研究存在局限性。首先，本研究以区（县）的气象数据替代个体的气象暴露情况，可能存在一定的测量偏倚；其次，本研究仅控制大气污染物PM₁₀，未考虑其他大气污染物对气象因素与人群健康暴露反应关系的影响，但之前有研究指出空气污染对温度-死亡之间的关系影响不大^［34］；再次，既往研究表明，温度变异也是影响人群健康的一个重要气象因素^{［35, 36, 37］}，本研究在确定气象健康风险阈值时尚未考虑温度变异的影响，在以后的研究中要考虑温度变异的影响。

本研究根据体感温度与YLL率的暴露-反应关系，采用回归决策树模型确定不同人群的冷、热效应体感温度的风险预警阈值，经验证预警效果良好。本研究为我国建立气象健康风险预警系统提供科学依据，有利于精准降低气象健康危害。