近年来，随着计算机技术和人工智能的快速发展，机器学习算法在媒介生物传染病中的应用逐渐广泛，已成为媒介生物传染病预测预警的重要工具[3]。例如，神经网络、随机森林、支持向量机已较好的应用于疟疾、布鲁氏菌病、肾综合征出血热和登革热的发病率预测中[4-7]，取得良好效果。中国鼠疫自然疫源地分布广、类型多、结构复杂，这决定了我国鼠疫防控工作的艰巨性和长期性[8]。鼠疫预测预警可提前发现鼠疫流行区域和时间，从而及时采取干预措施，优化防控策略，是鼠疫防控工作的重要基础和关键环节[9]。目前我国发现的15种鼠疫自然疫源地中[10]，仅有少数学者针对特定疫源地对动物鼠疫流行进行预测研究，一方面是由于建模预测相关研究需要长期、连续的数据作为支撑，目前我国鼠疫监测预警体系还不够完善，获得连续鼠疫监测资料较为困难；二是由于动物鼠疫流行是一个复杂的过程，宿主动物生态行为、气候变化及人类活动均会影响动物鼠疫的流行。新疆地区复杂的生态地理环境为鼠疫研究提供了丰富的资源基础。准噶尔盆地大沙鼠鼠疫自然疫源地是我国新发现的疫源地，2005－2016年连续开展了大范围综合性鼠疫疫源性调查和研究，鼠疫流行病学和宿主媒介调查资料连续、完整，不仅适于复杂动物鼠疫流行预测研究，也为鼠疫预测预警体系的建设和发展开辟了新的思路和方向。本文收集了新疆准噶尔盆地大沙鼠鼠疫自然疫源地古尔班通古特沙漠（古尔班通古特沙漠大沙鼠鼠疫自然疫源地）2005－2016年动物鼠疫监测和气象资料，分别利用BP神经网络（back propagation neural network，BPNN）和多元线性逐步回归（multiple linear stepwise regression，MSR）算法建立该疫源地动物鼠疫流行预测模型，评价筛选出预测精度较高的模型，为该地区鼠疫防控提供理论依据和技术支持。
1 材料与方法
1.1 数据来源与处理
古尔班通古特沙漠大沙鼠鼠疫自然疫源地动物鼠疫监测资料来源于新疆疾病预防控制中心2005－2016年的鼠疫监测数据，数据包括该区域昌吉回族自治州的木垒哈萨克自治县、奇台县、吉木萨尔县、阜康市、昌吉市、呼图壁县、玛纳斯县和乌鲁木齐市米东区8个县（市、区），数据指标包括大沙鼠血清鼠疫耶尔森菌（Yersinia pestia，鼠疫菌）F1抗体阳性率（血清抗体阳性率）、染蚤率和蚤指数、捕获数、密度，以及大沙鼠洞群覆盖率、洞群密度和洞群栖息率，鼠疫监测数据调查和计算采用文献[11]中描述的方法进行。气象数据由国家气象数据网（https://data.cma.cn）下载，采用距离阜康市鼠疫监测点最近的气象台站（51377）站点数据，包括以月为尺度的降水量和日照时数。
1.2 统计分析
采用SPSS 25.0软件中的Wilcoxon秩和检验对气象分组的降水量和日照时数进行差异性分析。采用Spearman秩相关分析方法对大沙鼠密度、染蚤率等生物学监测因子和气象因子与血清抗体阳性率的相关性进行分析。
1.3 预测模型构建
模型构建前，对数据集进行划分。将古尔班通古特沙漠大沙鼠鼠疫自然疫源地2005－2014年数据作为训练集，用于模型训练，2015－2016年数据作为测试集，用来验证模型预测效果，输出变量（应变量）为血清抗体阳性率，输入变量（自变量）为经相关分析筛选出的与血清抗体阳性率显著相关的生物及气象因子。
1.3.1 MSR模型
采用SPSS 25.0软件中的逐步法对自变量进行筛选，考察引入模型的自变量与应变量间的回归关系是否仍具有统计学意义，判断其是否应继续保留在方程中。以P值最小且具有统计学意义为依据进行自变量的引入和剔除。具体标准为：引入自变量的显著性概率为P < 0.05，剔除自变量的显著性概率为P > 0.10。对回归方程进行F检验，若F > $ {F}_{0.01} $ ，说明建立的回归方程有统计学意义。
1.3.2 BPNN模型
BPNN的基本结构由输入层、隐藏层、输出层以及各层之间的节点连接构成，通过梯度下降技术、利用反向传播算法不断调整网络权重和阈值，从而最小化网络实际输出值与目标预测值之间的误差。在BPNN参数设置中，隐藏层激活函数采用非线性双曲正切S型传递函数（hyperbolic tangent sigmoid transfer function，Tansig），其能将输入数据映射到非线性输出空间，从而允许神经网络学习复杂的数据模式和关系。训练算法使用trainlm函数，训练次数为1 000次，学习速率设为0.01，训练目标最小误差设为0.001，输出层神经元传递函数采用线性函数（linear function）。此外，隐藏层神经元节点数的确定至关重要，若隐藏层神经元节点数过多，会增加计算量并可能导致过拟合；过少则神经网络可能无法充分提取样本信息，从而达不到好的预测效果。因此，本研究通过经验公式$ q=\sqrt[]{k+m}+\mathrm{a} $ 来确定隐藏层神经元节点数。其中，k为输入层神经元个数，即自变量数，m为输出层神经元个数，即应变量数，a是一个可调节的常数因子，范围为[1, 10]，以Spearman秩相关分析筛选出的与血清抗体阳性率相关的变量作为输入层，血清抗体阳性率作为输出层，利用MATLAB 2024a软件建立BPNN模型用于动物鼠疫流行预测。
1.4 模型精度评价
表1 新疆维吾尔自治区大沙鼠鼠疫自然疫源地古尔班通古特沙漠动物鼠疫流行强度分级Tab. 1 Animal plague epizootic intensities classification in the Gurbantunggut Desert at the natural foci of great gerbil plague in Xinjiang Uygur Autonomous Region |
| 大沙鼠血清鼠疫菌F1抗体阳性率（%） | 对应流行强度分级 |
| 0 | 0 |
| 0 < X≤4.10 | 1 |
| 4.10 < X≤9.40 | 2 |
| 9.40 < X≤16.90 | 3 |
| X > 16.90 | 4 |
注：X代表抗体阳性率。 |
其次，通过计算决定系数（coefficient of determination，R2）、均方根误差（root mean square error，RMSE）和平均绝对百分比误差（mean absolute percentage error，MAPE）对预测模型进行精度评价和误差分析，R2越大、RMSE和MAPE越小，预测效果越好。利用MATLAB 2024a软件计算上述误差指标，各指标计算公式如下：
式中，$ {y}_{i} $ 表示实际值，$ {\widehat{y}}_{i} $ 表示预测值，$ {\stackrel{-}{y}}_{i} $ 表示实际值的平均值，n表示样本数量。
2 结果
2.1 动物鼠疫监测结果
2005－2016年共捕获大沙鼠3 973只，染蚤鼠3 446只，梳检体蚤43 565只，平均染蚤率为（87.22±7.19）%；平均蚤指数10.58±3.25；检验大沙鼠血清4 069份，鼠疫菌F1抗体阳性血清390份，平均阳性率为（10.25±6.73）%；洞群覆盖率、洞群栖息率和洞群密度的平均值分别为（25.17±10.64）%、（75.64±15.54）%和（5.68±6.93）个/hm2，各年份的动物鼠疫监测结果见表 2。
表2 2005－2016年新疆维吾尔自治区大沙鼠鼠疫自然疫源地古尔班通古特沙漠动物鼠疫监测结果Tab. 2 Surveillance data of animal plague in the Gurbantunggut Desert at the natural foci of great gerbil plague in Xinjiang Uygur Autonomous Region, 2005-2016 |
| 年份 | 阳性率b（%） | 染蚤率（%） | 蚤指数 | 鼠密度（只/hm2） | 洞群覆盖率（%） | 洞群栖息率（%） | 洞群密度（个/hm2） |
| 2005 | 7.00 | 93.10 | 14.77 | - | 33.00 | 72.00 | 25.36 |
| 2006 | 13.40 | 87.20 | 6.30 | 8.52 | 24.00 | 77.00 | 12.08 |
| 2007 | - | - | - | 22.6 | 11.00 | - | - |
| 2008 | 1.60 | 74.00 | 11.41 | 8.35 | 6.00 | 37.00 | 6.78 |
| 2009 | 11.20 | 91.70 | 14.98 | 4.68 | 22.00 | 88.00 | 1.78 |
| 2010 | 2.80 | 82.60 | 10.58 | 1.78 | 20.00 | 67.00 | 1.39 |
| 2011 | 3.70 | 79.00 | 8.39 | 2.77 | 25.00 | 82.00 | 1.45 |
| 2012 | 8.50 | 95.40 | 12.35 | 6.63 | 22.00 | 97.00 | 1.24 |
| 2013 | 20.80 | 89.00 | 5.23 | 12.13 | 35.00 | 87.00 | 2.35 |
| 2014 | 20.40 | 97.00 | 9.86 | 16.39 | 49.00 | 86.00 | 3.28 |
| 2015 | 18.10 | 79.00 | 14.61 | 6.89 | 29.00 | 62.00 | 3.56 |
| 2016 | 5.30 | 91.40 | 7.95 | 8.46 | 26.00 | 77.00 | 3.17 |
| 平均值a | 10.25±6.73 | 87.22±7.19 | 10.58±3.25 | 9.02±5.80 | 25.17±10.64 | 75.64±15.54 | 5.68±6.93 |
注：a平均值用 |
2.2 气象资料分布
2005－2016年该区域平均降水量为225.6 mm，最大降水量为319.5 mm，最小为126.7 mm，中位数为224.9 mm；平均日照时数为2 590.7 h，最大日照时数为2 774.0 h，最小值为2 143.3 h，中位数为2 615.4 h，见表 3。气象分组的降水量和日照时数间均存在统计学差异（Z=-3.175，P=0.001；Z=-4.158，P < 0.001）。
表3 新疆维吾尔自治区大沙鼠鼠疫自然疫源地古尔班通古特沙漠气象资料统计Tab. 3 Statistics of meteorological data of the Gurbantunggut Desert at the natural foci of great gerbil plague in Xinjiang Uygur Autonomous Region |
| 变量 | 平均值 | 最小值 | 中位数 | 最大值 |
| 年降水量（mm） | 225.6 | 126.7 | 224.9 | 319.5 |
| 秋末冬春降水量（mm） | 79.5 | 43.0 | 78.6 | 128.4 |
| 春末夏秋降水量（mm） | 146.1 | 79.8 | 159.3 | 232.7 |
| 年日照时数（h） | 2 590.7 | 2 143.3 | 2 615.4 | 2 774.0 |
| 秋末冬春日照时数（h） | 773.6 | 471.5 | 771.1 | 900.3 |
| 春末夏秋日照时数（h） | 1 817.1 | 1 671.6 | 1 832.7 | 1 906.8 |
2.3 影响血清抗体阳性率相关因子筛选
2.3.1 鼠疫生物学监测因子
Spearman秩相关分析结果显示，大沙鼠密度和洞群覆盖率与血清抗体阳性率间存在正相关关系（ρ=0.648，P=0.043；ρ=0.706，P=0.015），洞群栖息率、洞群密度、染蚤率和蚤指数与血清抗体阳性率的相关关系无统计学意义（ρ=0.465，P=0.150；ρ=0.145，P=0.670；ρ=0.055，P=0.873；ρ=-0.182，P=0.593）。
2.3.2 气象因子
Spearman秩相关分析结果显示，滞后1年春末夏秋平均日照时数与血清抗体阳性率间存在正相关关系（ρ=0.645，P=0.032），滞后1年年均降水量与血清抗体阳性率存在负相关关系（ρ=-0.609，P=0.047），其余气象因子与血清抗体阳性率的相关关系无统计学意义（均P > 0.05）。见表 4。
表4 大沙鼠鼠疫自然疫源地古尔班通古特沙漠气象因子与大沙鼠血清鼠疫耶尔森菌抗体阳性率相关性Tab. 4 Correlation between meteorological factors and serum yersinia pestis antibody positivity rate in the Gurbantunggut Desert at the natural foci of great gerbil plague |
| 气象因子 | 相关系数 | P值 |
| 年均降水量 | 0.082 | 0.811 |
| 滞后1年年均降水 | -0.609 | 0.047 |
| 春末夏秋平均降水 | 0.291 | 0.385 |
| 滞后1年春末夏秋平均降水 | -0.278 | 0.408 |
| 秋末冬春平均降水 | -0.127 | 0.709 |
| 滞后1年秋末冬春平均降水 | -0.127 | 0.709 |
| 年均日照时数 | 0.087 | 0.800 |
| 滞后1年年均日照时数 | 0.533 | 0.091 |
| 春末夏秋平均日照时数 | 0.087 | 0.800 |
| 滞后1年春末夏秋平均日照时数 | 0.645 | 0.032 |
| 秋末冬春平均日照时数 | 0.000 | 1.000 |
| 滞后1年秋末冬春平均日照时数 | 0.427 | 0.190 |
2.4 预测模型构建
2.4.1 MSR模型
选择2005－2014年血清抗体阳性率作为应变量，大沙鼠密度、洞群覆盖率、滞后1年年均降水量和滞后1年春末夏秋平均日照时数作为自变量，利用逐步回归法进行建模，经过模型自动识别，最终余下洞群覆盖度1项在模型中，回归系数为0.390（t=2.802，P=0.023）。模型通过了F检验（F=7.849，P=0.023），R2=0.495，其他结果见表 5。模型公式为：$ Y=0.390X+0.003 $ ，式中Y为血清抗体阳性率，X为洞群覆盖率。
表5 大沙鼠鼠疫自然疫源地鼠疫生物学监测因子和气象因子与大沙鼠血清鼠疫耶尔森菌抗体阳性率多元线性逐步回归模型结果Tab. 5 Results of multiple linear stepwise regression model of biological and meteorological factors and serum Yersinia pestis antibody positivity rate at the natural foci of great gerbil plague |
| 参数 | 非标准化系数 | 标准化系数（β值） | t值 | P值 | R²值 | F检验 | ||
| B值 | 标准误 | F值 | P值 | |||||
| 常量 | 0.003 | 0.038 | - | 0.090 | 0.931 | 0.495 | 7.849 | 0.023 |
| 洞群覆盖率 | 0.390 | 0.139 | 0.704 | 2.802 | 0.023 | |||
注：- 表示值空缺。 |
2.4.2 BPNN模型
将大沙鼠密度、洞群覆盖率、滞后1年年均降水量和滞后1年春末夏秋平均日照时数进行归一化处理，作为输入层，输出层为血清抗体阳性率，对数据集进行网络训练。采用公式$ q=\sqrt[]{k+m}+\mathrm{a} $ 确定隐藏层神经元节点数，此处k=4，m=1，a为[1, 10]之间的常数，确定隐藏层神经元节点数区间为[3, 12]。通过10次试验结果可见，当隐藏层神经元节点数为5时，RMSE值最小，为0.029，训练结果最佳，见图 1。
选择2005－2014年数据为训练集进行训练，2015－2016年数据进行测试，最终选出训练集R2=0.790、测试集R2=0.892的模型为最佳模型。
2.5 2种模型预测精度评价
图2 大沙鼠鼠疫自然疫源地动物鼠疫2种模型拟合预测图注：BPNN反向传播神经网络模型；MSR多元线性逐步回归模型。 Fig. 2 Fitting plots of the two prediction models for animal plague at the natural foci of great gerbil plague |
表6 大沙鼠鼠疫自然疫源地动物鼠疫2种预测模型训练及测试结果Tab. 6 Training and test results of the two prediction models for animal plague at the natural foci of great gerbil plague |
| 样本 | 年份 | 实际阳性率（%） | 流行强度分级 | MSR预测阳性率（%） | 流行强度分级 | BPNN预测阳性率（%） | 流行强度分级 |
| 训练集 | 2005 | 7.00 | 2 | 13.20 | 3 | 14.00 | 3 |
| 2006 | 13.40 | 3 | 9.70 | 3 | 10.00 | 3 | |
| 2007 | 10.30 | 3 | 4.60 | 2 | 9.30 | 2 | |
| 2008 | 1.60 | 1 | 2.60 | 1 | 1.50 | 1 | |
| 2009 | 11.20 | 3 | 8.90 | 2 | 16.10 | 3 | |
| 2010 | 2.80 | 1 | 8.10 | 2 | 2.60 | 1 | |
| 2011 | 3.70 | 1 | 10.10 | 3 | 3.40 | 1 | |
| 2012 | 8.50 | 2 | 8.90 | 2 | 8.40 | 2 | |
| 2013 | 20.80 | 4 | 14.00 | 3 | 20.00 | 4 | |
| 2014 | 20.40 | 4 | 19.40 | 4 | 20.00 | 4 | |
| 测试集 | 2015 | 18.10 | 4 | 11.60 | 3 | 17.10 | 4 |
| 2016 | 5.30 | 2 | 10.40 | 3 | 8.10 | 2 |
注：阳性率为大沙鼠血清鼠疫耶尔森菌F1抗体阳性率；MSR为多元线性逐步回归；BPNN为反向传播神经网络模型。 |
比较2种预测模型的R2、RMSE和MAPE，BPNN在训练集和预测集上的R2分别为0.790和0.892，明显高于MSR模型的0.496和0.167，在数据拟合和预测方面有更高的精度。进一步看，BPNN在训练集和测试集上的RMSE和MAPE均明显小于MSR模型，预测精度较高。见表 7。
表7 大沙鼠鼠疫自然疫源地动物鼠疫2种模型拟合预测精度比较Tab. 7 Comparison of prediction accuracy of the two models for animal plague at the natural foci of great gerbil plague |
| 模型 | 训练集 | 测试集 | |||||
| RMSE | MAPE | R2值 | RMSE | MAPE | R2值 | ||
| MSR | 0.046 | 0.659 | 0.496 | 0.058 | 0.661 | 0.167 | |
| BPNN | 0.029 | 0.209 | 0.790 | 0.021 | 0.292 | 0.892 | |
注：BPNN反向传播神经网络模型；MSR多元线性逐步回归模型；R2决定系数；RMSE均方根误差；MAPE平均绝对百分比误差。 |
3 讨论
构建预测模型时，模型预测变量的选择至关重要，纳入过多预测变量会增加模型复杂性，易引发过拟合问题，进而降低模型预测精度和泛化能力[14]。因此，选择最具代表性和关键性的变量作为预测因子，不仅能够简化模型结构，还能保证预测效果、提高模型实用性和可操作性。筛选与动物鼠疫流行相关的影响因子，是建立动物鼠疫预测模型的基础和提升预测性能的关键。既往研究表明，动物鼠疫流行与生物和气象因子显著相关[15-16]。因此，本研究通过分析大沙鼠密度、洞群覆盖率、洞群密度等6个生物因子及年均降水量、秋末冬春平均降水量和年均日照时数等12个气象因子的相关性，共筛选出大沙鼠密度、洞群覆盖率、滞后1年春末夏秋平均日照时数和滞后1年年均降水量4个与动物血清鼠疫菌F1抗体阳性率显著相关的因子进行模型构建。其中，大沙鼠密度和洞群覆盖率可作为衡量疫源地宿主密度的指标，且与血清抗体阳性率存在正相关关系，这与杨永刚等[17]、周晓磊等[18]的研究结果基本一致。其原因可能是宿主密度变化影响了寄生蚤数量，间接影响动物鼠疫流行。研究表明，鼠密度与鼠体蚤指数呈正相关，而体蚤指数升高会增加对宿主的叮咬次数，从而提高了鼠疫菌在宿主间的传播概率[19-20]。此外，本研究发现气象因子对动物鼠疫流行的影响存在滞后性且不同气象因子对动物鼠疫流行的影响存在差异，如滞后1年春末夏秋平均日照时数与血清抗体阳性率存在正相关关系，而滞后1年年均降水量与血清抗体阳性率存在负相关关系。这种现象在以往研究中也有发现，孙静等[12]研究了降水对新疆北天山长尾黄鼠（Spermophilus undulatus）鼠疫流行的影响，发现当年和前1年年降水量、春末夏秋和秋末冬春降水量均与长尾黄鼠血清抗体阳性率呈负相关；尚猛等[21]研究发现日照时数的增加会促进宿主动物繁殖，从而增加鼠疫宿主动物数量，提高鼠疫流行概率。
建模方法选择是影响动物鼠疫流行预测精度的一个重要因素。以往研究中，学者们利用灰色预测和判别分析，对我国特定鼠疫疫源地动物鼠疫流行进行定量预测研究[22-23]，这类模型简单易用，但忽略了多种因素对鼠疫流行的综合影响。本研究基于生物和气象因子建立了MSR和BPNN 2种动物鼠疫流行预测模型。MSR分析中，当自变量间存在共线性问题时，采用“逐步法”可剔除冗余信息，以确保自变量之间不存在高度相关性，并可简化模型，提高模型的解释性和预测精度[24]。石杲等[25]采用逐步回归法建立达乌尔黄鼠（Spermophilus dauricus）鼠疫流行预测模型，发现模型对当年鼠疫流行拟合率为93.33%，次年预报准确率为71.43%。然而，该模型适合处理自变量与应变量间存在线性关系的情况，在处理非线性关系时，可能导致预测精度下降。BPNN是一种基于误差逆向传播算法进行训练的多层前馈网络，不仅能存储和学习大量输入输出数据，而且无需对变量之间的映射关系进行明确表述，具备出色的非线性处理能力[26]。本研究对古尔班通古特沙漠大沙鼠鼠疫自然疫源地动物鼠疫流行进行预测研究，发现BPNN要优于MSR模型。此结果与李华宇[13]预测达乌尔黄鼠鼠疫自然疫源地鼠疫流行模型研究的结论一致。
准噶尔盆地大沙鼠鼠疫自然疫源地动物鼠疫流行区域广泛，包括古尔班通古特沙漠流行区和准噶尔盆地西部流行区，各流行区由于地理地貌、生态结构和气候等因素的差异，动物鼠疫流行存在异质性[2]，因而本研究建立的BPNN模型仅适用于准噶尔盆地古尔班通古特沙漠区域动物鼠疫流行预测，尚未对外推至准噶尔盆地大沙鼠鼠疫自然疫源地的其他流行区域的可靠性进行验证和评估，后期有必要进行验证，或针对各个流行区域分别建立动物鼠疫预测模型，为该地区鼠疫防控提供全面的理论依据和技术支持。