一种登革热传播风险预测方法及受影响人口大小确定方法

技术领域

本发明涉及传染病传播风险分析技术领域，尤其涉及一种基于SIR动力学模型和集合校正卡尔曼滤波耦合设计的登革热传播风险预测方法及受影响人口大小确定方法。

背景技术

登革热是近年来世界上广泛流行的热带传染病之一。有研究估计每年约有3.9亿人感染登革热，约75％的登革热负担集中在东南亚和西太平洋地区。传染病的流行级别与地理传播范围息息相关。自从20世纪中期以来，全球化促进了人口的快速流动，特别是交通科技的进步带来了短期旅行的增长，使得全球人口流动的速度、数量以及涉及的地域都在飞速增长，加剧登革热的传播和流行。虽然已有一些研究致力于了解蚊媒传播疾病在当前和未来的疾病风险地理分布以帮助指导干预和保护公众健康，但是在日益加快的全球化背景下，对登革热在人群中引起暴发及不同规模的驱动因素的了解尚不清楚。因此，亟需针对登革热流行高风险地区的时空传播动力学机制开展深入研究，建立准确的登革热实时预测模型。这将有助于加强登革热监测、推动防控方案实施和减少全球人口登革热疾病负担。

我国华南地区是登革热的主要流行区，并且在2014年经历了近10年来最严重的登革热疫情。近年来，我国登革热流行主要有以下特点：(1)登革热每年流行规模存在明显差异，这种差别不仅体现在时间尺度上(以广州市历年登革热流行为例，峰值低时可少于10例，高时可达8000例)，也体现在地理空间上(在同时期，珠三角地区登革热病例数远远超过非珠三角地区病例数)；(2)在流行高峰时期前后的几周内病例数消长速度快，而在一年的大部分时间病例发生频率低，使得登革热监测数据呈现明显的零膨胀(zero－inflated)特征。所谓零膨胀数据是指实际研究中所收集的计数资料中观察事件发生数含有大量的零值。由于在单位时间内包含过量零计数数据的随机事件，这种数据超出了Poisson分布等一般离散分布的预测能力[。我们在前期研究中对广东省登革热流行的季节性和周期性进行综合评估，发现当地登革热病例时间序列呈明显的零膨胀特征。这种特殊的离散分布数据给模型预测方法带来了挑战。

登革热给我国造成严重的经济损失和社会负担。因此，亟需针对我国登革热流行高风险地区的时空传播动力学机制开展深入研究，并在此基础上采用数学建模技术建立准确的登革热实时预测模型。因此，为了更好地预防和控制登革热在大规模人群的流行和暴发，建立准确的登革热预警模型是当务之急。

已有相关研究采用概率型统计预测方法(主要包括季节性自回归滑动平均模型模型、广义可加模型和人工神经网络等)和基于传染病传播动力学特性的仓室模型(也称为传染病动力学模型，如SIR等)构建登革热预测模型。上述两大类方法存在一些明显的缺点。下面分别阐述以上两类预测方法存在的缺点。

第一，传统的概率型统计模型存在的主要问题包括：1、依赖于大量的输入数据，可能不适合小样本的研究；2、对于发生频率较低的疾病，模型的基本假设可能无法满足；3、有些模型误差需要满足独立性条件。然而，传染病时间序列是高度相关的。直接使用这类模型对环境－病原体－宿主之间的关联进行分析可能会导致有偏估计；4、这类方法往往不考虑传染病在人群中传播的机制和特性，忽视了相关重要信息。导致该方法存在相应缺陷的主要原因包括：1、这类模型的预测性能往往受到模型参数的影响，而对参数进行优化和调整以寻找最优参数组合是一大方法学挑战；2、这类模型对历史数据依赖性很高，且往往不考虑传染病传播过程相关因素(如人群易感性、传播媒介的性质等)的影响，导致模型预测能力有待进一步提升；3、由于我国登革热发病人数时间序列具有明显的零膨胀特征，导致上述研究方法均不能很好地拟合这种特殊的离散分布数据。第二，经典的动力学模型存在的主要问题包括：1、模型中参数的确定对模型预测结果影响较大，即具有明显的参数初始值敏感性；2、模型预测误差的非线性增长导致部分动力学模型生成的传染病疫情预报不准确。这可能归结于模型没有准确地对非线性误差的结构进行量化和调整。这也提示，为了抵消这些由非线性模型动态产生的结构性误差，需要建立一种模型误差动态校正技术。

发明内容

本发明实施例所要解决的技术问题在于，提供一种登革热传播风险预测方法。可实现动力学模型系统状态的实时更新，实现模型误差动态校正，从而优化模型预测。

为了解决上述技术问题，本发明实施例提供了一种登革热传播风险预测方法，包括以下步骤：

S1：获取历史蚊媒信息数据，转换成蚊媒信息的自然出生率时间序列；获取历史气温数据，转换为随温度变化的登革热病毒经蚊子叮咬的转播至人的概率时间序列；

S2：利用所述自然出生率时间序列、传播至人的概率时间序列构建蚊媒介导的登革热传播SIR动力学模型；

S3：由所述登革热传播SIR动力学模型预测未来一个单位时间内的新增感染病例数，通过EAKF算法与新增感染病例的观测值同化，并校正动力学模型中的状态参数；

S4：在观测时间连续的情况下，反复步骤S3，达到SIR－EAKF的耦合；

S5：所述SIR－EAKF的框架在给定初始状态参数，输入登革热新增感染病例观测流数据下，对每个观测时间点同化观测数据，并校正所述登革热传播SIR动力学模型状态参数集合，返回后验状态参数集合，以及由所述后验状态参数集合向前预测，得到未来的先验状态参数集合；

S6：由所述先验状态参数集合计算未来的集合均值、集合区间以及集合预测峰值，峰值时间和总感染病例数。

其中，所述S1中，所述获取历史蚊媒信息数据为获取诱蚊诱卵指数，并采用傅里叶变换方法进行平滑。

其中，所述S1中，所述获取历史气温数据为获取每日平均温度时间序列。

其中，所述S2中，所述登革热传播SIR动力学模型纳入状态参数包括：易感蚊子数量S

其中，所述登革热传播SIR动力学模型采用龙格－库塔法估计微分方程。

其中，所述登革热传播SIR动力学模型的状态参数需完成参数寻优，所述参数寻优采用随机搜索法或经验判断法。

其中，所述登革热传播SIR动力学模型在所述SIR－EAKF的框架中提供预测手段，预测未来一段时间内的状态变化趋势，并作为系统过程耦合EAKF。

其中，所述EAKF算法中对于未观测的状态参数，采用线性最优算法进行更新校正。

其中，所述S4中，所述初始状态参数采用蒙特卡洛模拟法，基于各自参数的分布和范围，模拟各集合的状态参数随机组合；所述观测流数据，包括新增感染病例数和其相应的观测方差。

本发明实施例还提供了一种受疫情影响的人口大小自确定方法，由历史数据计算当前受疫情影响的人口大小，并使用上述的SIR－EAKF的固定参数，重新运行所述SIR－EAKF对已获得的观测流数据的同化更新。

实施本发明实施例，具有如下有益效果：本发明构建SIR－EAKF耦合算法，实现整个预测系统状态的实时更新，校正由非线性模型动态产生的结构性误差，降低模型中参数初始值的敏感性，提升预测结果的准确性和鲁棒性。通过EAKF算法和蚊媒介导SIR模型的耦合设计，实现模型和参数随观测数据流的同化而更新，获取状态参数的精确估计值，从而，该耦合算法能够对登革热流行和传播的动态变化趋势进行精确地预测，避免了传统统计模型中对历史数据的高度依赖以及对零膨胀时间序列的失效的问题，解决了经典动力学模型中参数初始值敏感性和模型预测误差非线性增长的问题。实现了模型随观测时间变化的状态参数更新，从而更准确的进行传染病疫情预报，并返回动力学模型中其他系统状态的变化预测和流行病学参数的估计。

附图说明

图1为根据本发明实施例的自然出生率时间序列集合构造流程图；

图2为根据本发明实施例的登革热病毒经蚊叮咬传播概率时间序列集合构造流程图；

图3为根据本发明实施例的SIR－EAKF耦合框架流程图；

图4为根据本发明实施例的受疫情人口影响大小自确定嵌合SIR－EAKF框架流程图；

图5为根据本发明实施例的登革热疫情传播实时预测的实际操作使用流程示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作进一步地详细描述。

本发明实施例的一种登革热传播风险预测方法，是基于SIR动力学模型和集合校正卡尔曼滤波的耦合算法应用于登革热病毒传播的预测中

，包括以下步骤：

1)、获取历史的蚊媒信息数据，转换成蚊媒信息的自然出生率时间序列。获取历史气温数据，转换为随温度变化的登革热病毒经蚊子叮咬转播至人的概率的时间序列；

2)、利用所述历史蚊媒信息，气象温度信息构建蚊媒介导的登革热传播SIR动力学模型；

3)、由所述蚊媒介导的登革热传播SIR模型预测未来一个单位时间内的新增感染病例数，通过EAKF算法与新增感染病例的观测值同化，并校正动力学模型中其他的状态参数。在观测时间连续的情况下，反复所述步骤，实现SIR－EAKF的耦合。

4)、所述SIR－EAKF框架，给定初始状态参数，输入登革热新增感染病例观测流数据，即可在每个观测时间点同化观测数据，并校正SIR模型状态参数集合，返回后验的状态参数集合，以及由所述后验状态参数集合向前预测，得到未来一段时间内的先验状态参数集合。

5)、提供一种受疫情影响的人口大小自确定方法，由历史数据计算当前受疫情影响的人口大小，并作为SIR－EAKF的固定参数，重新运行SIR－EAKF对已获得的观测流数据的同化更新。

6)、由所述先验状态参数集合，计算未来的集合均值，集合区间，以及集合预测峰值，峰值时间和总感染病例数。

本发明构建SIR－EAKF耦合算法，实现整个预测系统状态的实时更新，校正由非线性模型动态产生的结构性误差，降低模型中参数初始值的敏感性，提升预测结果的准确性和鲁棒性。通过EAKF算法和蚊媒介导SIR模型的耦合设计，实现模型和参数随观测数据流的同化而更新，获取状态参数的精确估计值，从而，该耦合算法能够对登革热流行和传播的动态变化趋势进行精确地预测，避免了传统统计模型中对历史数据的高度依赖以及对零膨胀时间序列的失效的问题，解决了经典动力学模型中参数初始值敏感性和模型预测误差非线性增长的问题。实现了模型随观测时间变化的状态参数更新，从而更准确的进行传染病疫情预报，并返回动力学模型中其他系统状态的变化预测和流行病学参数的估计。

在上述实施方式中，对于实施例构建登革热病毒疫情传播实时预测的SIR－EAKF耦合算法包括以下步骤：

A、针对实施例中获取的蚊媒信息数据为诱蚊诱卵指数，通过图1的流程，将其转换为自然出生率时间序列集合，具体实施如下：

a、首先采用傅里叶变换方法对诱蚊诱卵指数进行平滑；

b、通过将一个季度的诱蚊诱卵指数除以一阶导数为极小值的情况下的诱蚊诱卵指数，得到标准化的自然增长比例序列；

c、自然状态下，假定蚊子的产卵周期为16天，则自然出生率时间序列为上述除法得到的标准化自然增长比例序列再开16次方根再减1；

d、计算历史数据中自然出生率时间序列的季度序列均值及标准差，假设自然出生率时间序列的误差均服从正态分布，通过蒙特卡洛模拟法在该分布中随机抽取自定义集合数量的次数，模拟出相同次数的自然出生率时间序列集合。

B、针对实施例中获取的气温数据为每日平均温度时间序列，通过图2的流程，将其转换为登革热病毒经蚊叮咬传播概率时间序列集合，具体实施如下：

a、每日平均温度与登革热病毒经蚊叮咬传播概率的函数关系符合：当温度为18℃时，登革热病毒经蚊叮咬传播概率为0％；当温度为23℃时，登革热病毒经蚊叮咬传播概率为13.3％；当温度为28℃时，登革热病毒经蚊叮咬传播概率为40％；当温度为32℃时，登革热病毒经蚊叮咬传播概率为66.7％。假定温度低于18℃时，登革热病毒经蚊叮咬传播概率为0％；温度高于32℃时，登革热病毒经蚊叮咬传播概率为66.7％。以线性方式连接各节点，估计不同温度下登革热病毒经蚊叮咬传播概率；

b、针对获得的登革热病毒经蚊叮咬传播概率时间序列，计算经平滑处理后的历史数据中登革热病毒经蚊叮咬传播概率时间序列的季度序列均值及标准差；

c、假设登革热病毒经蚊叮咬传播概率时间序列的误差均服从正态分布，通过蒙特卡洛模拟法在该分布中随机抽取自定义集合数量的次数，模拟出相同次数的登革热病毒经蚊叮咬传播概率时间序列集合。

C、针对登革热病毒传播的动力关系，构建蚊媒介导的登革热传播SIR动力学模型，其中模型纳入状态参数包括：易感蚊子数量S

D、由参数寻优法采用随机搜索法或经验判断法，设置所述蚊媒介导的登革热传播SIR动力学模型的初始状态参数。

E、设置集合数量ens{1,2,...j,...,m}，所述初始状态参数，采用蒙特卡洛模拟法，基于各自参数的分布和范围，模拟各集合的状态参数随机组合。

F、构建蚊媒介导的登革热传播SIR动力学模型集合，给定t-1时刻的系统状态和参数PB{X

G、登革热新增感染病例数为观测流obs{n}，由经验法则给出观测误差流obs_var{n}＝(obs{n}^2+100)/25。

H、由t-1时刻蚊媒介导的动力学模型后验集合PB{X

I、通过校正卡尔曼滤波计算后验方差

J、其他系统状态和参数的增益在线性最优的假设下计算为

K、针对本实施例，上述状态参数用于EAKF同化的状态是新增感染病例NewI

L、最终实现t时刻整个系统状态参数的同化更新

本发明实施例还提供了一种受疫情影响的人口大小自确定方法，由历史新发感染病例数的时间序列季度中各时间点对总感染病例数的平均占比，估计当前观测中的受疫情影响的人口大小。

通过估计t时刻的受疫情影响的人口大小，作为新的受疫情影响的人口大小固定参数输入SIR－EAKF框架，重新运行{1,2,...,t}时间段模型的同化更新，实现SIR－EAKF嵌合受疫情影响的人口大小自确定方法，如图4所示。

所述受疫情影响的人口大小自确定方法，设置从第5周开始截至平均峰值时间22周运行该方法。

所述蚊媒介导的登革热传播SIR动力学模型采用龙格－库塔法估计微分方程。

所述蚊媒介导的登革热传播SIR动力学模型在SIR－EAKF框架中提供预测手段，预测未来一段时间内的状态变化趋势，并作为系统过程耦合EAKF，完成如图3所示的SIR－EAKF耦合框架预测。

所述SIR－EAKF框架从每年的第20周开始投入使用，至第二年的第19周截止，完成对该季度的预测，以降低零膨胀时间序列对EAKF同化更新的影响。

在实际操作中，步骤D的初始状态参数默认应用于中国广东省广州市的登革热预警中，已经通过参数寻优法确定，应用于中国疫情预测也应具有普适性，可不作变动。

根据学习任务的不同，可选择是否启用受疫情人口影响大小自确定方法。

综上所述，本发明提出一种能够提高预测效率和预测准确性的基于SIR动力学模型和集合校正卡尔曼滤波耦合算法的登革热传播风险预测方法。基于SIR动力学模型和集合校正卡尔曼滤波的耦合算法，实现整个预测系统状态的实时更新，校正由非线性模型动态产生的结构性误差，降低模型中参数初始值的敏感性，提升预测结果的准确性和鲁棒性。为适应登革热的动力学传播及季节性变化，构建登革热传播的SIR动力学模型，可以反映人、蚊系统状态的变化，以及准确的反映登革热的季节性特征。本方法可以快速地同化登革热观测流数据并更新系统状态参数，并准确的预测未来一段时间登革热传播的发展态势，十分适用于目前登革热预报的应用场景。

本发明基于R软件实现代码编程，并完成函数封装。如图5所示，仅需输入初始参数及范围和观测流数据即可运行本发明，返回各观测时间点的未来的预测状态参数集合和后验的状态参数集合。

以上所揭露的仅为本发明一种较佳实施例而已，当然不能以此来限定本发明之权利范围，因此依本发明权利要求所作的等同变化，仍属本发明所涵盖的范围。