蚊虫是多种传染病的传播媒介，如流行性乙型脑炎（乙脑）、疟疾、登革热、丝虫病、黄热病和寨卡病毒病等。近几年来蚊媒传染病在我国发病率持续上升，且已成为病死率较高的媒介传染病^[1]。连续系统的蚊密度监测是有效控制病媒生物的基础和依据。差分自回归移动平均（autoregressive integrated moving average，ARIMA）模型是根据一组随时间变化而互相关联的动态数据而形成的一种分析方法，用于描述任何齐次非平稳时间序列^[2]，是时间序列分析方法中最常用的1种。本研究收集唐山市2010-2018年的蚊虫密度和季节消长调查数据，将2010-2017年唐山市的蚊密度监测数据引入ARIMA模型进行建模，预测2018年的蚊密度，对蚊密度的动态趋势进行预测预警，有效预防蚊媒疾病的传播，为唐山市未来的防蚊工作提供科学依据。

汇总2010-2018年唐山市各监测点蚊虫数量和季节消长数据。唐山市2010-2018年蚊密度监测时间为每年的5-10月，本研究中1-4月和11-12月的蚊密度设定为0只/（灯·h）。

依据《全国病媒生物监测方案（试行）》^[3]中的诱蚊灯法，选择唐山市城区居民区、公园、医院和郊区农村民房及牲畜棚，统一使用“功夫小帅”牌诱蚊灯（武汉吉星环保科技有限公司生产），远离干扰光源，避风，挂灯位置距离地面1.5 m，从日落后20 min开始接通电源，连续诱捕6 h，次日鉴定蚊种并计数，每月上、下旬各监测1次。

应用Excel 2007软件录入2010-2018年蚊密度监测数据，用SPSS 22.0软件进行数据分析，构成比的比较采用χ²检验，检验水准α=0.05。利用ARIMA模型进行建模与预测。

从时间序列本身的相关性出发，将各种数学表达式纳入模型，并努力创建时间序列的回归方程以达到预测未来值的目的。ARIMA（p，d，q）×（P，D，Q）s模型中p、d、q分别表示自回归的阶、差分的阶和移动平均的阶，P、D、Q则是指季节性的阶数，s表示季节性周期。

① 序列平稳化:对非平稳序列进行差分处理转化为平稳序列；②模型识别和创建:基于自相关和部分相关函数建立可能的模型；③参数估计:模型参数的估计和假设检验；④参数诊断:判断所建立的模型是否合适，模型的残差序列为白噪声序列，各参数显著，可用于预测；⑤预测应用:利用2010-2017年各月总蚊密度数据所建立的模型，预测2018年各月的总蚊密度，并根据2018年实际监测密度评估预测效果。

2010-2017年共捕获成蚊22 335只，其中淡色库蚊（Culex pipiens pallens）占捕获总数的98.02%，为优势蚊种；中华按蚊（Anopheles sinensis）占1.69%，白纹伊蚊（Aedes albopictus）占0.23%，三带喙库蚊（Cx. tritaeniorhynchus）占0.02%；其他蚊种占0.04%，见表 1。

表 1 河北省唐山市2010-2017年不同生境蚊虫数量及其构成 Table 1 Numbers and composition ratios of mosquitoes in different habitats of Tangshan, Hebei province, China, 2010-2017

表选项

在不同监测生境中，蚊种构成排在首位的均为淡色库蚊。不同生境蚊种构成比差异有统计学意义（χ²=249.177，P < 0.001），对不同生境蚊种构成进行两两比较，除居民区和公园蚊种构成比差异无统计学意义（P＞0.05）外，其他生境两两比较差异均有统计学意义（P < 0.05），见表 1。

2010-2017年共布放诱蚊灯870盏，成蚊总密度为2.14只/（灯·h），在5类生境中，牲畜棚中蚊密度最高为3.65只/（灯·h），居民区蚊密度最低为1.54只/（灯·h）；从不同年份看，2014年蚊密度最高为3.00只/（灯·h），2016年蚊密度最低为1.41只/（灯·h），见表 2。

表 2 2010-2017年河北省唐山市不同年份各生境蚊虫数量（只）及密度〔只/（灯·h） Table 2 Numbers and densities (individuals/lamp·hour) of mosquitoes in different years in different habitats of Tangshan, Hebei province, China, 2010-2017

表选项

2010-2017年总蚊密度以7月最高，为4.71只/（灯·h），各年份蚊密度均为单峰曲线，除2012、2013和2017年高峰出现在8月外，其余年份蚊密度高峰均出现在7月，见表 3。淡色库蚊在7月达到高峰，密度为4.70只/（灯·h）；中华按蚊在8月达到高峰，密度为0.18只/（灯·h）。在不同生镜中，蚊密度也呈单峰曲线，高峰均出现在7月，居民区、公园、医院、农村民房和牲畜棚7月密度分别为3.19、3.06、4.04、4.16和9.09只/（灯·h），见图 1。

表 3 2010-2017年河北省唐山市不同月份蚊密度〔只/（灯·h）〕 Table 3 Mosquito densities (individuals/lamp·hour) in Tangshan, Hebei province, China, different months of 2010-2017

表选项

图 1 2010-2017年河北省唐山市不同生境蚊密度季节消长情况 Figure 1 Seasonal fluctuations in mosquito density in different habitats of Tangshan, Hebei province, China, 2010-2017

图选项

将2010-2017年唐山市蚊密度绘制成时间序列图，各年蚊密度均呈明显的周期性。季节周期为12个月，该序列为非平稳序列，见图 2。为使平稳性满足ARIMA模型的需要，对数据进行1次差分和1次季节性差分，随机序列近似平稳，见图 3。

图 2 2010-2017年河北省唐山市蚊密度时间序列结果 Figure 2 Time series plot of mosquito density in Tangshan, Hebei province, China, 2010-2017

图选项

图 3 一阶非季节性差分和一阶季节性差分（12个月为1个周期）后的蚊密度序列结果 Figure 3 Time series plot of mosquito density after first-order non-seasonal difference and first-order seasonal difference (12 months were defined as 1 cycle)

图选项

原始序列经过一阶非季节差分和一阶季节差分，得出ARIMA（p，1，q）×（P，1，Q）12模型，12表示12个月为1个周期，对差分后平稳的时间序列做自相关（ACF）和偏自相关（PACF）图（图 4）。

图 4 原序列经一阶非季节差分和一阶季节差分后的自相关和偏自相关 Figure 4 Autocorrelation and partial autocorrelation of the original series after first-order non-seasonal difference and first-order seasonal difference

图选项

根据差分后的ACF和PACF图，p、q、P和Q从低到高选择，根据模型的参数及平稳R²和标准化贝叶斯信息准则（BIC）值进行分析，选择ARIMA（0，1，1）×（2，1，0）₁₂和ARIMA（1，1，0）×（2，1，0）₁₂为备选模型，依据平稳R²越大，BIC值越小拟合越好，最终选择ARIMA（0，1，1）×（2，1，0）₁₂作为最佳模型，见表 4。

表 4 ARIMA模型拟合结果比较 Table 4 Comparison of the fitting results of ARIMA models

表选项

采用Ljung-Box检验方法得出残差是白噪声序列（Q=20.654，P=0.148），表明残差是随机误差，残差ACF和PACF均落入95%可信区间（95%CI）范围内，见图 5。模型的参数评估结果见表 5，最终得出的模型数学表达式:（1+0.016B¹²+0.454B²⁴）（1-B）（1-B¹²）Y_t=（1- 0.585B）e_t。

注:ACF.自相关；PACF.偏自相关。 图 5 ARIMA（0，1，1）×（2，1，0）12模型残差序列的自相关和偏自相关 Figure 5 Autocorrelation and partial autocorrelation of the residual series of the ARIMA (0, 1, 1)×(2, 1, 0)12 model

图选项

表 5 ARIMA（0，1，1）×（2，1，0）₁₂模型参数评估 Table 5 Parameter estimation of the ARIMA (0, 1, 1)× (2, 1, 0)₁₂ model

表选项

根据建立的模型预测2018年的总蚊密度，5-10月的预测值分别为0.09、2.51、3.45、4.61、1.05和0.22只/（灯·h），实际的蚊密度分别为0.07、2.26、3.65、4.56、0.76和0.27只/（灯·h），二者的平均相对误差为19.44%，＜20%为合格模型，可用于预测。在预测过程中实际值均落入预测值95%CI范围内，说明蚊密度正常，用此模型预测的蚊密度季节消长趋势与实际值基本吻合，表明此模型能够较好地拟合实际的蚊密度，可用于中短期蚊密度预测，见图 6、表 6。

图 6 2010-2018年河北省唐山市蚊密度实际值与模型拟合值及预测值分布 Figure 6 Distributions of actual mosquito densities and the fitted and predicted values of the model in Tangshan, Hebei province, China, 2010-2018

图选项

表 6 河北省唐山市2018年5-10月蚊密度〔只/（灯·h）〕预测值与实际值 Table 6 Predicted and actual mosquito densities (individuals/ lamp·hour) in Tangshan, Hebei province, China, during May to October, 2018

表选项

2010-2017年在唐山市区共捕获成蚊22 335只，其中淡色库蚊为该地区的优势蚊种，中华按蚊、白纹伊蚊和三带喙库蚊的数量相对较少，中华按蚊是疟疾的主要传播媒介^[5]；白纹伊蚊可以传播登革热、寨卡病毒病和基孔肯雅热；三带喙库蚊可以传播乙脑^[6]。当前，唐山市主要的蚊媒传染病有乙脑和疟疾，做好蚊虫的监测是防治蚊媒传染病最有效的措施。

蚊虫的孳生和繁殖受到温度和降雨等气候的影响，唐山市地处渤海湾中心地带，位于河北省东部，属于温带季风气候，夏天炎热多雨，降雨主要集中在7-8月，导致积水增多，蚊密度迅速达到高峰，2010-2017年唐山市月平均蚊密度7月达到最高峰，各年蚊密度高峰均出现在7、8两月，由此也证明蚊密度与气温和降雨呈一定的相关关系^[7-8]。

时间序列是通过利用研究对象历史数据之间随时间不断变化而得出的规律预测未来值的一种统计方法^[9-10]。ARIMA模型考虑了季节因素对时间序列平稳性的影响，对短期趋势预测准确率较高^[11]。目前，ARIMA模型已应用于各种传染病发病率的预测^[12-14]，本研究将ARIMA模型应用于唐山市蚊密度分布特征预测的研究中，从2010-2017年，ARIMA模型用于拟合唐山市各月的蚊虫密度，最后选择ARIMA（0，1，1）×（2，1，0）₁₂作为最佳模型，并对2018年的蚊密度进行预测，预测值与实际值拟合较好，整体动态趋势和季节规律基本吻合，模型预测效果较为理想。该方法可以应用于举办大型活动前对蚊虫及蚊媒传染病的风险评估中。在举办活动前应用ARIMA模型对蚊密度进行预测，将其作为预警值，在活动举办期间，结合实际监测的蚊密度可更直观迅速的进行预警，当实际蚊密度达到或超过预警值时，应立即预警，并迅速组织应急队伍及时有效地开展灭蚊工作，并根据蚊密度与蚊媒疾病的关系，及时提醒相关医疗单位为蚊媒疾病的暴发流行做好准备，从而科学有效地保障大型活动时蚊媒疾病的防治。

本研究ARIMA模型可较好地预测未来的蚊虫密度，将拟合的数学模型应用于媒介蚊虫监测网络信息系统^[15]，并根据不同地区成蚊密度与蚊媒疾病之间发生的联系，乘以一定的系数可以用作网络监测信息系统中蚊密度的预警值^[16]，结合不同地区蚊密度上报的实际数值和动态趋势迅速准确地进行预报预警，从而可以第一时间监测到疑似疫情信息。

综上所述，今后可以采用本研究拟合的ARIMA（0，1，1）×（2，1，0）₁₂模型对唐山市的蚊密度进行早期预测，并结合实际监测情况及时进行预警，为唐山市蚊媒传染病的防控提供科学依据。