花粉症是由致敏花粉导致的过敏性鼻炎, 近年来全球范围花粉症患病率和严重性均呈逐年增加态势^[1]。我国多中心流行病学研究显示, 花粉过敏源是北方引发过敏性鼻炎最主要的过敏原, 发病人口数达7500万^[2]。大量研究表明花粉症发病季节与气传花粉浓度高峰期吻合, 且患者就诊人次与致敏花粉浓度呈高度相关^[3—5]。北京作为我国北方地区超大城市之一, 常住人口数量大、人员居住密集, 伴随绿色生态建设推进带来植被状况改善, 致敏花粉引发的问题也日趋明显^[6—7]。

北京地区植物资源丰富, 花粉致敏植物种类繁多, 前人已开展大量研究工作探讨北京地区花粉致敏植物种类、分布及变化等特征, 结果显示北京地区花粉总量年变化呈双峰型, 春季高峰期以乔木植物为主, 秋季高峰期以草本植物居多^[8—11]。研究表明不同时段内的气象条件是影响不同种类花粉散粉期和散粉量的主要自然因素之一^[12—15], 基于气象条件对花粉浓度影响分析, 国内外大量学者开展了致敏植物生育期、花期以及不同时间尺度花粉浓度预测的研究工作^[16—18]。针对北京地区：张德山等^[19]应用正交筛选回归方法, 使用多气象因子建立了基于站点的北京地区未来1—4天花粉浓度预报模式；赵文芳等^[20]提出了基于粒子群优化算法和支持向量机、使用多种气象要素提升花粉未来24小时浓度预测精度的方法；卞萌等^[21]利用具有外部输入的非线性自回归神经网络模型, 结合气候因子及叶面积指数观测, 建立了北京市次日花粉浓度预测模型。

相较于不同时间尺度的花粉浓度预测工作, 鲜有关于花粉浓度高峰期出现时间的预测研究。本研究使用北京地区多年分类花粉浓度和气象要素观测数据, 基于花粉年浓度变化规律及年浓度高峰期内优势致敏花粉种类分析, 使用模糊逻辑相关方法探讨气象条件对北京地区主要气传致敏花粉年浓度峰值出现时间的影响, 结合北京地区高分辨率区域天气预报模式, 可对北京不同地区主要气传致敏花粉年浓度峰值出现日期进行预测, 为花粉过敏症预防和临床诊断工作提供参考。

1 花粉样本采集概况

本研究使用花粉浓度观测数据为北京地区12个花粉采样站点2012至2020年逐日花粉分类浓度观测数据。所有花粉采样站点均布设在所在区国家级地面观测站观测场内, 站点分布如图 1所示。

2012至2020年花粉采样期集中在每年3至10月, 各年、各站实际采样期长度如图 2所示。花粉样本使用Durham花粉采样器基于叶氏重力沉降法获取, 样本采样片每24小时更换一次, 获取后使用碱性复红染色液对采样片进行染色处理, 之后使用显微镜进行人工分类计数, 花粉浓度计数结果为单位面积花粉数量(粒/(500mm²))。

2 研究方法 2.1 花粉年浓度峰值对应日期确定方法

某类花粉年浓度峰值为该年花粉季内其逐日浓度观测值的最大值^[22], 浓度最大值出现日期即为该种类花粉年浓度峰值日期。为避免年花粉总量过少或同一花粉季存在多个浓度峰值等少数情况对建模影响, 当出现以下三种情况中的任意一种, 则认为某年度某类花粉浓度峰值存在明显不确定性, 对该年数据进行剔除处理：(1)年采样期内某种花粉浓度最大值小于50粒/(500mm²)；(2)年采样期内某种花粉存在多个数值接近的花粉浓度峰值, 并且出现日期相差超过1周；(3)某年花粉浓度峰值对应日期明显偏离2012至2020年平均值, 相差日期大于30d。

2.2 气象要素累积值计算方法

本研究使用气象数据为花粉采样站点对应的国家级地面观测站2012至2020年气象要素日值观测资料, 包含日平均气温、日平均地表温度、日降水量和日照时数。基于气象要素日值观测资料, 统计逐年、各站不同种类花粉年浓度峰值日期对应气象要素累积值, 包含有效积温(本研究所有种类花粉对应生物学下限温度均使用0℃)、累积地表温度(大于0℃日平均地表温度累积值)、累积降水量(日降水量累积值)和累积日照时数。

其中对于在6月1日之前其浓度达到峰值的花粉种类, 其气象要素累积值计算起始日期为同年1月1日；在6月1日之后浓度达到峰值的花粉种类, 其气象要素累积值计算起始日期为同年6月1日。花粉年浓度峰值日期对应有效积温、累积地表温度、累积降水量和累积日照时数为年浓度峰值日期当天、前三天以及后三天(共计7天)气象要素累积值的平均值。

2.3 模糊逻辑原理

20世纪60年代, Zadeh^[23]使用了一种基于模糊集和模糊逻辑的方法来处理语言学中的模糊信息。此后这一方法不断发展并逐步形成了模糊逻辑系统理论, 广泛应用于气象学、医学、计算机科学等多种领域。基于隶属度函数的概念, 使用模糊逻辑可以模仿人的不确定性判断和推理方式, 描述并区分模糊集合, 处理模糊关系, 从而对经典逻辑难以处理的非线性、不确定问题进行定量化描述^[24]。因此可应用模糊逻辑理论将界限不清晰的定性知识和经验转化为确定性定量关系。

2.4 模型建立、检验方法

本研究使用多种气象要素的模糊逻辑隶属度函数的组合, 描述气象条件对不同种类花粉年浓度峰值出现时间的非线性影响, 其中不同种类花粉对应有效积温、累积地表温度、累积降水量和累积日照时数的隶属度函数通过花粉年浓度峰值日期对应气象要素日值数据累积值的统计结果拟合得出。

基于逐日气象要素累积值, 使用公式P=A×B+(1-A×B)×(C×D)预测花粉年浓度峰值出现日期。其中P为某一日是当年某种花粉浓度达到峰值的概率, A、B、C、D分别表示4种气象要素隶属度函数, A、B为影响某种花粉浓度峰值日期的主要隶属度函数。使用训练集通过排列组合方式筛选出年浓度峰值日期计算值与观测值绝对误差天数最小的组合形式, 最终针对北京地区花粉年浓度高峰期内优势花粉种类建立相应的花粉年浓度峰值日期预测模型, 并使用测试集数据验证模型预测准确率。

模型预测准确率定义为测试集中年浓度峰值日期计算值与观测值误差小于或等于±7d的样本量占测试集样本总量的比例^[25]。并选择模型预测日期与观测日期的平均绝对误差(Mean absolute error, MAE)和均方根误差(Root mean square error, RMSE)等统计指标反映模型模拟误差大小和整体模拟效果。

3 结果与分析 3.1 北京地区花粉年浓度变化规律

基于2012至2020年逐日花粉分类浓度观测数据, 计算北京地区3至10月逐日花粉总浓度平均值, 其变化规律如图 3灰色部分所示。北京地区年花粉浓度存在两个明显高峰期：其一为3月上旬至5月中旬, 其二为8月中旬至9月中旬, 该规律与前人研究成果基本一致^[8—11]。其中北京地区3月上旬至5月中旬花粉浓度高峰期可进一步划分为3月中旬至4月上旬以及4月下旬至5月上旬两个高峰期。

3.2 北京地区花粉年浓度高峰期内优势花粉种类

分别对3月1日至4月15日、4月15日至5月15日、8月15日至9月15日北京地区各类花粉总浓度年平均值进行统计排序, 结果如图 4所示：北京地区3月上旬至4月中旬花粉浓度高峰期内柏科和杨柳科花粉浓度占比为57.4%, 4月中旬至5月中旬花粉浓度高峰期内优势花粉种类为豆科和松科植物花粉, 8月中旬至9月中旬花粉浓度高峰期内桑科、菊科蒿属和藜科植物花粉浓度占比总和为77.7%。同时可以看到, 北京地区3月上旬至5月中旬花粉浓度高峰期以木本植物为主, 8月中旬至9月中旬高峰期则以草本植物为主。

结合北京地区主要致敏花粉种类^[26—29], 确定本研究以柏科、杨柳科、松科植物花粉为北京地区3月上旬至5月中旬花粉浓度高峰期内优势花粉种类, 桑科、菊科蒿属、藜科植物花粉为8月中旬至9月中旬花粉浓度高峰期内优势花粉种类, 其花粉浓度年变化规律如图 3所示。

3.3 主要气传致敏花粉年浓度峰值及其对应日期年际、站间差异

对2012至2020年北京地区花粉年浓度高峰期内6种优势花粉种类(柏科、杨柳科、松科、桑科、菊科蒿属、藜科植物花粉)年浓度峰值存在不确定性的数据进行剔除后, 分别对6种花粉在不同采样站点、不同年份下花粉年浓度峰值及其对应日期差异进行统计, 结果如图 5所示。

柏科和杨柳科花粉年浓度峰值主要出现在3月中旬到4月中旬, 平均出现时间为3月底, 松科花粉年浓度峰值对应日期波动范围较大, 为4月中旬到5月下旬, 平均出现时间为5月上旬。相较于3月上旬至5月中旬花粉浓度高峰期内优势花粉种类, 8月中旬至9月中旬3种主要致敏花粉年浓度峰值日期波动范围较小, 集中在8月下旬到9月上旬出现, 平均出现时间为9月初。

对于同一种类花粉, 其年浓度峰值存在采样站点间和年际间显著差异, 年浓度峰值显著偏离平均值的情况普遍存在。而年浓度峰值对应日期的观测站点间和年际间差异较小, 即花粉年浓度峰值出现时间的可预测性高于年浓度最大值的可预测性, 因此花粉年浓度峰值对应日期的预测研究对于花粉症防治更具指导意义。

3.4 北京地区主要气传致敏花粉年浓度峰值日期预测模型

花粉浓度变化与致敏植物生育期、尤其是花期关系密切, 而气象条件是直接影响植物生长发育的主要环境因素^[30]。作物模型以光、温、水、土壤等条件作为环境驱动变量, 再现植物生长发育过程^[31]。基于作物模拟模型概念, 采用模糊逻辑原理, 可建立以气象要素为驱动的数值模拟模型对致敏植物生长发育的生物学过程进行动态模拟^[32], 从而实现气传致敏花粉浓度峰值出现日期的预测。针对上述北京地区花粉年浓度高峰期内优势花粉种类, 本文选择可表征光、温、水等环境因素变化的气象要素累积值来描述气象条件对致敏植物生长发育的动态影响, 基于多年各类花粉年浓度峰值日期对应当年多种气象要素累积值统计分析结果, 建立北京地区主要气传致敏花粉年浓度峰值日期预测模型。

3.4.1 花粉年浓度峰值日期与气象要素累积值的相关性分析

使用剔除年浓度峰值存在不确定性的数据, 基于北京地区6种主要气传致敏花粉(柏科、杨柳科、松科、桑科、菊科蒿属和藜科植物花粉)多站、多年的年浓度峰值对应日期及匹配的气象要素累积值数据, 构建原始数据集。

选择有效积温、累积地表温度、累积降水量和累积日照时数表征温度、水分、光照等气象要素对致敏植物生长发育影响^[33]。使用原始数据集采用Spearman相关系数分别计算不同花粉年浓度峰值日期对应日序与气象条件累积值之间的相关性, 结果如表 1所示。6种花粉年浓度峰值出现时间与不同气象要素累积值之间的相关性存在差异, 有效积温、累积地表温度与花粉浓度峰值出现时间之间有较强相关性, 而累积降水量和累积日照时数与花粉浓度峰值出现时间(尤其是桑科、菊科蒿属和藜科花粉)的相关性较弱, 使用传统回归模型难以有效表达多种气象要素变化对花粉浓度峰值出现日期的综合影响, 因此本文引入模糊逻辑原理, 使用不同气象要素隶属度函数组合的形式描述光、温、水等环境条件对花粉浓度峰值出现时间的非线性影响。

3.4.2 隶属度函数的建立

对6种花粉分别随机选取原始数据集中80%的样本作为模型训练集, 用于确定气象要素隶属度函数及其组合方式；剩余20%样本作为测试集, 用于判断模型预测准确率。基于模型训练集, 分别统计6种北京地区主要气传致敏花粉年浓度峰值日期对应有效积温、累积地表温度、累积降水量和累积日照时数的分布特征, 具体方法如下。

以柏科花粉的有效积温隶属度函数建立为例：经统计, 2012至2020年北京地区12个花粉采样站柏科花粉年浓度峰值日期对应有效积温变化范围为110℃d至320℃d, 各区间对应个例数量如图 6所示。将个例数最大值对应区间的隶属度设为1, 其余区间对应隶属度为本区间个例数与最大个例数的比值。最终基于划分区间与隶属度值计算结果, 使用分段函数拟合气象要素分布规律, 建立柏科花粉有效积温隶属度函数(图 6虚线所示), 其中分段函数拐点为个例数最大值对应区间左侧区间左端点及右侧区间右端点。

基于训练集数据统计得出北京地区6种主要气传致敏花粉年浓度峰值日期对应4种气象要素累积值分布规律(图 7、图 8)。依照上述隶属度函数拟合方式, 分别建立柏科、杨柳科、松科、桑科、菊科蒿属和藜科花粉对应有效积温、累积地表温度、累积日照时数以及累积降水量的隶属度函数, 函数具体表达式见附表1。

基于上述北京地区6种主要气传致敏花粉隶属度函数拟合结果, 使用模型训练集数据, 通过排列组合方式, 确定最优隶属度函数组合方式, 具体筛选方法如下：

根据公式P=A×B+(1-A×B)×(C×D), 设置A、B、C、D对应不同隶属度函数的全部组合方式, 分别计算6种花粉在其每年花粉季中, 逐日花粉浓度为该年浓度峰值的概率P, 选择其中最大值对应日期作为当年该种类花粉浓度达到峰值的日期。使用模型训练集计算每种花粉对应所有隶属度函数组合结果, 筛选出年浓度峰值日期计算值与观测值绝对误差天数最小的组合形式, 结果如表 2所示。

由隶属度函数组合方式筛选结果可以看到, 3月上旬至5月中旬花粉浓度高峰期内主要气传致敏花粉对应的木本植物生长对气温和日照时数条件敏感, 而8月中旬到9月中旬花粉浓度高峰期内优势花粉种类对应的草本植物生长对气温和降水条件更敏感。

3.4.3 模型预测效果检验

对柏科、杨柳科、松科、桑科、菊科蒿属以及藜科花粉年浓度峰值日期预测模型, 分别使用对应测试集验证模型预测准确率。检验结果为5次随机分组检验结果的平均值, 不同种类花粉对应模型检验结果如表 3所示。

柏科和桑科花粉对应模型预测准确率较高, 预测平均绝对误差值约4d, 杨柳科、菊科蒿属和藜科花粉对应模型预测准确率在80%左右, 平均绝对误差在5d左右, 松科花粉对应模型预测准确率较其他花粉对应模型较差, 预测平均误差以及预测结果波动性较大。整体看, 基于模糊逻辑隶属度函数建立的北京地区主要气传致敏花粉年浓度峰值日期预测模型可以较为准确地预测出不同种类花粉年浓度峰值出现时间, 从而为北京地区花粉症防治工作开展提供理论依据。

4 讨论 4.1 花粉分类精度对模型性能的影响

本研究建模使用的花粉分类原始数据, 除菊科蒿属外, 其余花粉采集样本均为分类到科数据。对于北京地区而言, 柏科、杨柳科、松科、桑科、藜科均包含多种不同属植物^[34], 其花期以及花粉年浓度峰值出现时间等均存在不同程度的差异。因此本文建立的花粉年浓度峰值日期预测模型性能在一定程度上受原始观测数据的花粉分类精度限制, 柏科、杨柳科和松科花粉对应模型预测准确率的差异, 可能与该科数据中所含属之间花期等因素的差异性有关, 北京地区松科花粉对应植物多样性及差异性较大可能是导致模型预测准确率较低的主要原因之一。

花粉分类精度对模型性能的影响还体现在不同属植物之间在生物学下限温度、花期影响因素等生物学特性的差异不确定性上。本研究在计算不同种类花粉年浓度峰值日期对应有效积温时, 生物学下限温度均采用0℃。前人研究表明不同下限温度设定对植物花期模拟有影响^[35], 对于生物学下限温度等生长发育适宜气象条件明确的致敏植物, 在建立对应气象条件隶属度函数的过程中应考虑其生物学特性, 以增强模型机理性。本研究中同一科不同属植物间生物学下限温度的差异对有效积温隶属度函数具体形式的影响、进而对模型性能的影响有待进一步研究。同时本研究计算气象要素累积值的起始时间为固定值, 研究表明不同植物在各生育期内对不同气象要素的敏感性存在差异^[36], 对于种类确定植物应选取对花期影响最显著气象要素的最显著影响时段, 建立对应气象要素的隶属度函数。由于原始数据分类精度限制, 本研究建模过程中不同花粉对应气象要素累积值计算起始时间的确定有待进一步研究, 其对模型性能的影响仍存在不确定性。

4.2 前期气候条件和临近天气条件对花粉年浓度峰值日期预测的影响

本文基于作物模型概念, 使用模糊逻辑原理建立了北京地区主要气传致敏花粉年浓度峰值日期预测模型, 模型以气象要素累积值为驱动, 可逐日动态判断花粉季内每日不同种类花粉浓度达到峰值的可能性。前人研究结果表明, 前期气候条件对植物花期等物候有重要影响^[37—38], 同时天气条件对当日花粉浓度变化也存在显著影响^[13—14]。因此在应用本研究建立模型预测花粉年浓度峰值日期时应考虑前期气候条件与临近天气变化对模型预测结果进行合理修正。

基于本研究建模原始数据集, 对柏科、杨柳科和松科花粉年浓度峰值出现日期对应日序与前一年冬季(前一年12月至当年2月)气象条件的相关性进行分析。结果表明当年花粉浓度峰值出现时间与前一年冬季平均气温和平均地表温度存在负相关关系, 柏科花粉对应负相关系数大于0.7(表 4)。研究表明气象条件对花粉浓度的影响存在滞后性^[39], 在实际应用模型预测当年花粉浓度峰值出现日期时, 应考虑前期气候条件对模型预测结果进行合理修正。

基于逐日花粉浓度观测数据, 分析逐日气象要素的变化对当日花粉浓度的影响, 结果如图 9所示。花粉日浓度高值集中出现在5—20℃气温区间内, 当气温低于5℃或高于20℃时, 花粉日浓度出现高值的可能性较小。在较低气压和较高水汽压(高湿)的情况下, 花粉日浓度出现高值的概率较小, 同时降水对于花粉日浓度有明显抑制作用。当日平均风速过小或大于4 m/s时, 花粉日浓度出现高值的可能性较小, 在偏西风情况下(西北或西南), 花粉浓度出现高值的可能性较大。上述结论与其他研究结果相似^[12—14], 即某日花粉浓度受当日气温、气压、湿度、降水和风向风速等气象要素影响明显, 在特定天气条件下, 花粉日浓度出现高值的可能性较大。在使用本研究建立模型预测花粉年浓度峰值日期时, 除考虑前期气候条件对花粉年浓度峰值日期影响外, 还应考虑预测日期前后具体天气情况, 对预测结果进行修正。

5 结论

本文基于2012至2020年北京地区12个花粉采样站点花粉季逐日分类浓度观测数据, 对北京地区花粉年浓度变化规律、年浓度高峰期内优势花粉种类及其浓度峰值出现日期特征等进行分析, 并基于作物模型概念和模糊逻辑原理, 以花粉年浓度峰值日期对应气象要素累积值为核心, 针对北京地区主要气传致敏花粉分别建立了年浓度峰值日期预测模型。本研究主要结论如下：

(1) 北京地区花粉浓度年变化规律呈现双峰型, 总浓度高峰期分别出现在3月上旬至5月中旬以及8月中旬至9月中旬, 其中3月上旬至5月中旬花粉浓度高峰期可进一步划分为3月中旬至4月上旬以及4月下旬至5月上旬两个高峰期。

(2) 北京地区3月上旬至5月中旬花粉浓度高峰期内优势致敏花粉种类为柏科、杨柳科和松科, 其中3月中旬至4月上旬致敏花粉以柏科、杨柳科为主, 4月下旬至5月上旬致敏花粉以松科为主；而8月中旬至9月中旬花粉浓度高峰期内优势致敏花粉种类为桑科、菊科蒿属和藜科植物花粉。

(3) 基于作物模型概念和模糊逻辑原理, 使用花粉年浓度峰值日期及其对应气象要素累积值(有效积温、累积地表温度、累积降水量、累积日照时数), 建立北京地区主要气传致敏花粉年浓度峰值日期预测模型。经检验, 柏科、杨柳科、桑科、菊科蒿属和藜科花粉年浓度峰值日期预测模型预测准确率分别为87.8%、80.0%、86.7%、78.8%和81.8%, 松科花粉模型预测准率较低为64.4%。使用基于模糊逻辑隶属度函数、以气象要素累积值作为驱动的花粉年浓度峰值日期预测模型, 综合花粉季前期气候条件以及花粉季内天气变化, 以北京区域高分辨率数值预报模式作为输入, 可实现北京地区主要气传致敏花粉年浓度峰值日期分区预测, 为花粉症防治提供理论参考。