[0001] 本发明涉及传染病风险监测技术领域，更具体的说是涉及一种基于蚊媒诱卵器指数的蚊媒传染病传播风险分析方法。

[0002] 目前，蚊媒传染病，是由蚊子作为蚊媒传播的传染性疾病，如疟疾、黄热病、寨卡病毒病、登革热等。蚊媒传染病是以“人‑蚊‑人”的传播模式进行传播的，即，通过雌性蚊子叮咬人类，将病原体从一个宿主传播到另一个宿主。随着城市化和国际旅行的增加，蚊媒传染病的传播风险显著增加，传播范围和频率不断扩大，且传播风险受多种因素影响，包括病原体特性、蚊媒生态学特征、气候条件和社会经济因素等。

[0003] 但是，现有技术中，往往采用病例数量反应传播风险，而只有当人们感染并被确诊后才会记录为病例，这往往意味着传播已经发生了一段时间，因此病例数量往往滞后于实际传播情况，仅能作为间接指标。

[0004] 因此，如何准确直观的分析传播风险实现实时监测是本领域技术人员亟需解决的问题。

[0041] 下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

[0042] 本发明实施例公开了一种基于蚊媒诱卵器指数的蚊媒传染病传播风险分析方法，包括以下步骤：

[0043] S1：获取病例数据和序列间隔分布数据；其中，序列间隔分布数据包括人类序列间隔分布和蚊媒序列间隔分布。

[0044] S2：通过计算人类序列间隔分布和蚊媒序列间隔分布的联合分布，得到风险分布。

[0045] S3：根据病例数据和风险分布，计算各个时刻的有效再生数，实现传播风险的实时监测分析。

[0046] 其中，有效再生数（ ），又称时变再生数，是指在特定时间 （例如：第 天），考虑实际的人群免疫水平、防控措施及社会行为变化等因素后，一个典型病例直接造成新感染者的平均数，用于描述传染病在特定时间点 的实时传播风险。当 ＞1时，传染病将会在人群中迅速传播，造成传播病情暴发；当 ＜1时，传染病传播趋向终止，传播病情逐渐解除。该指标是一个随时间变化的动态值，一系列不同时间点 上的有效再生数值按照时间顺序排列后构成了一个时间序列，揭示了传染病的实际传播动态。因此，使用有效再生数的时间序列分析蚊媒传染病的传播风险变化。

[0047] 在一种实施例中，有效再生数表示为 时刻的新病例数与截至 时刻累积的总传播风险之比，可通过更新方程来估算，并通过贝叶斯方法来估算其95%可信区间，具体公式如下：

[0048] （1）

[0049] 其中， 是 时刻的新病例数，加权平均值 是截至 时刻累积的总传播风险，n表示可成功传播的最大时间跨度（即序代间隔的最大时间跨度），权重 由序代间隔分布确定，该分布描述了感染期间传播风险的分布。

[0050] 在本实施例中，蚊媒传染病以蚊子作为病原的传播媒介，携带病原体的蚊子通过叮咬人类宿主，将病原传播给人类宿主的疾病（有传播蚊媒，即：“人‑蚊‑人”）。我们将感染的序代间隔（SI）分布 建模为感染人类和感染蚊媒的联合序代间隔分布， 可以通过感染人类的序代间隔分布 和感染蚊媒的序代间隔分布 “卷积”推导而得，具体公式如下：

[0051] （2）

[0052] 其中，k表示感染人类和感染蚊媒的序代间隔之和的随机变量；表示感染人类的序代间隔的随机变量。公式（2）中 和 的具体解释如下：将人类的序代间隔分布 设定为某种固定的概率分布，例如伽马分布。由于蚊媒在感染后可以持续传播病原体（给人类）直到死亡，蚊媒的序代间隔分布 需要同时考虑已经感染的情况下，已感染的蚊媒序代间隔分布 和蚊媒感染时间点的分布，通过已感染的蚊媒序代间隔分布和蚊媒感染时间点分布的“卷积”推导得到 。其中， 通常被设定为某种固定
的概率分布，变量 表示已感染的蚊媒的序代间隔随机变量。

[0053] 将蚊媒感染时间点（记为 ）建模为均匀分布，范围为雌性蚊媒寿命范围，从而得到蚊媒感染时间点的分布（用 ）。 是在考虑序代间隔为 的情况下，对于每个可能出现的蚊子感染时间点 ，计算蚊子生命剩余时间中对应的感染时间点分布的值。通过对不同的 与 进行（从0到 ）的积分，可以得到完整的蚊媒的序
代间隔分布。如下：

[0054] 。

[0055] 特别注意，是在卷积积分内部的被积分变量，当卷积积分完成后，所得到的结果与 无关。

[0056] 在另一实施例中，在得到有效再生数序列后，还可以基于有效再生数序列结合蚊媒诱卵器指数构建回归模型，并实现对有效再生值基线的估算。

[0057] 有效再生数基线值（ ），反映了传播风险的季节性波动，这种季节性由当地蚊子丰度和蚊子活动情况驱动，并间接受到其他影响蚊子种群的因素的影响。具体的，蚊媒传染病的传播媒介是蚊虫。蚊虫生长受到温度、降雨量、湿度等气象因素的影响。高温高湿降雨量大都会给蚊虫创造一个良好的栖息环境。我国广州地区作为亚热带季风气候，夏季高温、高湿、多雨。其夏季的气候特征为登革热传播媒介的繁殖提供过了有利条件。有效再生数基线由当地蚊子活动决定，当地蚊子的繁殖、存活及叮咬频率等活动状况构成了登革热病毒传播的基础条件，有效再生数基线可以量化这些影响，反映在自然状态下（无特殊防控或其他重大干扰因素时）蚊子活动所具有的病毒传播能力。而有效再生数是一个时刻水平。通过对有效再生数与其相应基线值的比较，可以判断蚊媒传染病传播风险的变化源于季节性波动或是其他因素影响。两者的比较过程是：若有效再生数＞有效再生数基线值，则表明在该时间点病毒传播处于相对活跃状态，可能存在其他促进传播的因素；反之，则说明传播相对受到抑制，可能得益于防控措施或其他不利传播的条件。因此，步骤还包括：

[0058] S4：获取蚊媒诱卵器指数数据，并根据所述蚊媒诱卵器指数数据和所述有效在生数构建第一回归模型。

[0059] S5：通过持续监测蚊媒诱卵器指数的变化，结合所述第一回归模型，实现早期预警。

[0060] S6：确认有效再生基线值，并将其与有效再生值进行比较，根据比较结果进行单侧检验。

[0061] 在本实施例中，通过当对前整体情势确认的第一回归模型进行最大似然估计，确认有效再生基线值作为有效再生数的稳定状态，结合实际的有效再生数值，判断该实际有效再生数相较稳定状态的显著性，当显著性达到单侧检验的显著性要求时，判断其为激增点。

[0062] 在本实施例中，蚊虫指数（例如：月度蚊媒诱卵器指数）反映某地区蚊子活跃度和种群丰度的总体状况。为估算有效再生数基线值，基于（有 天滞后的）月度蚊媒诱卵器指数的指数函数构建与有效再生数的回归模型，具体公式如下：

[0063] （3）

[0064] 其中， 是截距项， 是回归系数， 是时刻 的误差项。 是 的连续形式，即， ，其中z是被积分的时间变化量，并且积分后z消失。MOI( )表示时刻 的月度蚊媒诱卵器指数。根据蚊子从幼虫变为成虫的周期，设置滞后项x。 表示为：

[0065] （4）

[0066] 和 表示回归模型的参数点估计值，采用经典最大似然函数法（或最小二乘法）计算得到。假设在特定时间 ， 大于（或者小于） ，进行单侧假设性检验。

[0067] 为了进一步实施上述技术方案，单侧检验步骤包括：

[0068] 设置轮廓伪似然框架进行局部采样，根据局部采样数据构建第二回归模型，得到局部回归系数；由第一回归模型确认总体回归系数，计算局部回归系数相较总体回归系数的显著性；当显著性满足预设阈值时，根据当前轮廓伪似然框架输出相应的风险激增时点。

[0069] 具体的，传播风险激增时点，反映了在某一时点某种突发事件或条件的出现，导致风险在短时间内迅速扩散。进一步分析传播风险激增的驱动因素，如病原体的遗传变异、气象条件的改变、人口迁移模式等因素，可为探索传染病情爆发的潜在机制提供依据，有助于预测和应对未来的传染病情爆发。在现实世界中，我们无法通过传染病监测的方式来获得传播风险激增时点的准确时间，因此需要借助统计模型工具来计算这一时点。可以在这个轮廓（伪）似然框架中，通过构建一个与研究问题相关的似然函数（这里是伪似然函数），并寻找使该似然函数达到最大值的参数值，从而估计未知参数，这里的未知参数即为传播风险激增时点。

[0070] 通过比较公式（3）的总体回归与仅使用部分数据的局部回归之间的异质性，可以估算传播风险激增时点。总体回归是利用研究期间内所有数据进行拟合，用 表示总体回归的β估计值。局部回归仅限于使用1年周期内的数据，即时间窗口长度为期1年的数据，表示局部回归的 估计值。通过单侧t检验比较 和 ，得到每个 对应的显著性水平 值，从而确定最大（伪）似然估计值，这是传播风险激增时点出现可能性的测量值，这个过程即为确定伪似然框架。传播风险增高时，局部回归系数显著高于总体回归系数，即 ＞ ，假设“ ＞ ”。

[0071]

[0072] 因此，可以基于轮廓（伪）似然框架估算传播风险激增时点 的最大（伪）似然估计值及其95%置信区间。其中，95%置信区间是一种基于样本数据对总体参数（登革热病毒株入侵时间）进行区间估计的方法。它表示在多次重复抽样和参数估计的过程中，有95%的情况下所得到的区间会包含真实的传播风险激增时点t0参数值。对于本研究估计的新病毒株入侵时间不是单一的点时间估计（如最大似然估计值2014年4月），这考虑到了由于抽样误差和数据不确定性导致的估计不准确。简单来说，轮廓（伪）似然框架是寻找这个传播风险激增时点的过程，而最大似然估计值是我们想要确定的传播风险激增时点。

[0073] 为了证明本发明的有益效果，本发明采用实际数据进行了验证，具体验证过程如下：

[0074] 首先，根据监测数据可获得某地区多年内登革热每月病例数和蚊虫指数（每月蚊媒诱卵器指数，MOI）的时间序列。在数据分析前，对数据预处理，使用季节性差分自回归滑动平均模型（SARIMA模型）对缺失或不完整的历史数据进行插补，并对所有数据以“每日”为单位进行缩放处理，采用局部加权回归平滑法（LOESS）生成平滑曲线。

[0075] 其次，根据登革热病毒感染后的流行病学特征，将人类序代间隔分布设定为平均值为8天、标准差为1天的伽马分布（见图2中的部分a）；将蚊媒序代间隔分布设定为平均值为7天、标准差为1.5天的累积伽马分布（见图2中的部分b）。基于雌蚊寿命范围为5到8周，将蚊媒感染时间点的概率分布 设定为范围5到8周的均等分布。利用公式（2）进一步确定联合序代间隔分布（见图2中的部分c）。图2中，横轴为登革热感染后的时间。内源潜伏期和外源潜伏期分别用靛蓝和绿色表示。

[0076] 此外，利用公式（1）基于病例数估计有效再生数的时间序列，结果如图3所示。然后，利用公式（3）拟合有效再生数和月度蚊媒诱卵器指数，用于回归拟合的月度蚊媒诱卵器指数滞后15天，使用最小二乘法估计参数值，从而获得效再生数基线值，结果见图4中的部分a。其中，有效再生数估计值为红色点状曲线、95%可信区间为粉色竖条状、有效再生值基线为黑色曲线。

[0077] 比较有效再生数与其相对应的基线值，假设在特定时间 ，有效再生数大于（/小于）有效再生基线值，进行单侧假设检验，结果见图4中的部分b。

[0078] 其中，有效再生数估计值为红色点状曲线、95%可信区间为粉色竖条状、有效再生值基线为黑色曲线。

[0079] “ 大于（红色竖条）或小于（蓝色竖条）回归拟合基线 ”的单侧检验的显著性水平（ 值）（b）。

[0080] 图4中的部分c中，彩色曲线分别为2012‑2015年每年的局部回归拟合结果，值分别为“如果传播风险激增时点为2012‑2015年每年1月”的假设检验的显著性水平。部分d中表示的“2011‑2016年局部回归的回归系数（即斜率）（部分c中的彩色曲线）大于2008‑2018年整体回归的回归系数（a和c中的黑色曲线）”的单侧检验的显著性水平。

[0081] 最后，基于轮廓（伪）似然框架估算传播风险峰值出现的具体时间点。利用公式（3）分别以总体回归与局部回归方式拟合有效再生数和月度蚊媒诱卵器指数，估计总体回归系数 和局部回归系数 ，结果见图4c。并比较总体回归系数与局部回归系数，假设“＞ ”并进行单侧假设检验，结果见图4d。根据假设检验的显著水平（即：值），构建轮廓（伪）似然框架以估算传播风险激增时点及其95%置信区间，结果表明传播风险激增时点的最大似然估计值为2014年4月（前后2个月），见图5。

[0082] 先前基因研究发现在2014年大规模登革热传染病情中分离出了高传染性的登革热新病原株，系统发育分析也表明新病原株是在2014年传染病爆发前从国外输入的。这些研究结果与本技术估算的传播风险激增时点相一致，证明本发明能够准确确定传播风险激增时点，准确反映疫情的发展趋势，及时预测和响应疫情爆发。

[0083] 本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的装置而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。

[0084] 对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。