孫燕群，張守剛，陸墨原，3，張 艷，潘衍宇，王 沖，吳起新，姚美雪，李成國

1南京市疾病預(yù)防控制中心，南京醫(yī)科大學附屬南京疾病預(yù)防控制中心，江蘇 南京 210003；2軍事科學院軍事醫(yī)學研究院微生物流行病研究所，病原微生物生物安全國家重點實驗室，北京 100071；3東南大學公共衛(wèi)生學院，江蘇 南京 210009；4中國疾病預(yù)防控制中心傳染病預(yù)防控制所，北京 102206；5徐州醫(yī)科大學公共衛(wèi)生學院流行病與衛(wèi)生統(tǒng)計學系，江蘇 徐州 221004

蚊類（Mosquito）屬于昆蟲綱（Insecta）、雙翅目（Diptera）、長角亞目（Nematocera）、蚊科（Culicidae）。蚊蟲最重要的生態(tài)習性是刺叮吸血，從而傳播多種傳染?。?-3］。由于全球氣候、環(huán)境、交通、城市化的變化以及昆蟲自身抗藥性等因素，蚊蟲等病媒生物對人類的威脅不斷上升，重大病媒生物傳播疾病總共約占全球傳染性疾病負擔的17%，每年造成70多萬人死亡，其中以蚊蟲造成的疾病負擔最重［4］。蚊蟲種類多，分布廣，適應(yīng)能力強，易產(chǎn)生抗藥性，其侵害和控制已成為世界關(guān)注的焦點。

差分自回歸移動平均（differential autoregressive moving average，ARIMA）模型是由美國統(tǒng)計學家Box和英國統(tǒng)計學家Jenkins 于20 世紀70 年代初提出的時間序列分析、預(yù)測和控制的方法，又稱Box-Jenkins法，主要用于擬合具有平穩(wěn)性或者可以被轉(zhuǎn)換為平穩(wěn)序列的時間序列，結(jié)合了自回歸和移動平均的長處，具有不受數(shù)據(jù)類型束縛和適用性強的特點［5-6］，廣泛應(yīng)用在公共衛(wèi)生領(lǐng)域中的疫情預(yù)警預(yù)測等［7-9］，在病媒生物侵害密度領(lǐng)域應(yīng)用不多。本研究利用南京地區(qū)長期、連續(xù)、系統(tǒng)開展的蚊蟲侵害橫斷面調(diào)查數(shù)據(jù)，運用季節(jié)性差分自回歸移動平均（SARIMA）模型進行數(shù)據(jù)擬合，對本地區(qū)蚊蟲侵害數(shù)據(jù)開展預(yù)測研究，為進一步防控蚊媒病和開展愛國衛(wèi)生運動提供新的思路和方法。

1 資料和方法

1.1 資料

蚊蟲侵害數(shù)據(jù)來源于江蘇省疾病預(yù)防控制綜合業(yè)務(wù)集成平臺中的病媒生物監(jiān)測網(wǎng)絡(luò)直報系統(tǒng)，包括監(jiān)測點編號、監(jiān)測時間、蚊蟲數(shù)量等。

1.2 方法

1.2.1 研究變量

以蚊蟲密度（D）為主要研究對象，其計算公式為：D=Nm/（Nl·T），式中：Nm為蚊蟲數(shù)量，單位為只；Nl為燈的數(shù)量，單位為燈；T為誘蚊小時數(shù)，單位為h。

1.2.2 數(shù)據(jù)處理

以月為時間段，統(tǒng)計2015年1月—2019年12月南京地區(qū)南京蚊蟲密度，其中每年的1月、2月和12月為非蚊蟲活動時間，蚊蟲密度以0 填充。選擇2015 年1 月—2018 年12 月數(shù)據(jù)為訓練集，2019 年1—12月數(shù)據(jù)為驗證集。

1.2.3 模型預(yù)處理

純隨機性檢驗：純隨機性檢驗又稱白噪聲檢驗，是專門用來檢驗序列是否為純隨機的方法，使用Box.test 中的Box-Pierce 函數(shù)進行檢驗。單位根檢驗是對序列是否平穩(wěn)的檢驗，一般使用ADF檢驗進行判斷。

1.2.4 模型構(gòu)建

SARIMA 是ARIMA 中的一種，主要分為季節(jié)乘積和季節(jié)相加兩種，以季節(jié)乘積為多，主要處理序列的季節(jié)效應(yīng)、長期趨勢效應(yīng)和隨機波動之間存在復(fù)雜的交互影響關(guān)系。SARIMA 乘積模型公式為：，也可以簡化為SARIMA（p，d，q）×（P，D，Q）s，其中，p、d、q分別表示短期相關(guān)模型的自回歸、趨勢差分、移動平均的階數(shù)，P、D、Q 分別表示季節(jié)趨勢的自回歸、趨勢差分、移動平均的階數(shù)，s 表示周期時間間隔。使用Box.test 中的Ljung-Box 函數(shù)進行殘差白噪聲檢驗。經(jīng)過殘差診斷合格后可以對模型進行預(yù)測。

1.2.5 模型評價

選擇平均誤差（mean error，MSE）、均方根誤差（root mean squared error，RMSE）、平均絕對誤差（mean absolute error，MAE）、平均絕對比例誤差（mean absolute sacled error，MASE）和決定系數(shù)R2進行模型擬合效果的評價。

1.3 統(tǒng)計學方法

研究人員將蚊密度按照時間2015 年1 月—2019年12月進行整理，用Excel2013進行匯總，建立時間序列，利用R x64 4.0.3 軟件及下載的forecast、tseries包進行統(tǒng)計學分析，P＜0.05為差異有統(tǒng)計學意義。

2 結(jié)果

2.1 序列概況

繪制2015—2018 年蚊蟲侵害數(shù)據(jù)的時間序列分解圖（圖1），可以看出蚊蟲侵害數(shù)據(jù)時間序列季節(jié)效應(yīng)較明顯，而長期趨勢不明顯。

圖1 2015年1月—2018年12月時間序列分解圖

2.2 SARIMA模型構(gòu)建

2.2.1 模型準備階段

經(jīng)Box-Pierce 函數(shù)進行檢驗，序列為非白噪聲序列（P＜0.01），經(jīng)ADF 檢驗不平穩(wěn)（Dickey-Fuller=-2.275，P=0.466），后經(jīng)一階差分后ADF檢驗平穩(wěn)（Dickey-Fuller=-3.853，P＜0.05）。

2.2.2 SARIMA模型預(yù)測

參數(shù)判斷：分別對2 次差分繪制自相關(guān)函數(shù)（auto correlation function，ACF）圖和偏自相關(guān)函數(shù)（partial ACF，PACF）圖，差分后的ACF圖和PACF圖見圖2。1 階差分的自相關(guān)圖非截尾，q 考慮取0，偏自相關(guān)圖為2 階截尾，p 考慮可取1 和2；季節(jié)差分之后的自相關(guān)圖呈現(xiàn)顯著非零，Q 考慮取0，偏自相關(guān)圖在lag12 處有spike，P 考慮取1。利用可能的參數(shù)，分別建立ARIMA（1，1，0）（1，1，0）12或者ARIMA（2，1，0）（1，1，0）12模型，計算最小信息量準則（AIC）值，其中前者AIC=114.57，后者AIC=113.54，綜合考慮ARIMA（2，1，0）（1，1，0）12模型。

圖2 經(jīng)1階和1次季節(jié)差分后時間序列自相關(guān)圖和偏自相關(guān)圖

模型診斷：對ARIMA（2，1，0）（1，1，0）12進行Ljung-Box殘差檢測，可以認為該模型殘差序列不存在自相關(guān)，為白噪聲序列（χ2=0.079，P=0.778）。

模型預(yù)測：用forecast函數(shù)預(yù)測2019年1—12月的蚊蟲侵害密度，預(yù)測結(jié)果見圖3。

圖3 2015 年1 月—2018 年12 月實際時間序列及2019 年1—12月預(yù)測序列

2.2.3 SARIMA模型評價

SARIMA預(yù)測的2019年1—12月結(jié)果與實際監(jiān)測結(jié)果比較見表1。

表1 2019年1—12月時間序列預(yù)測與實際監(jiān)測比較

根據(jù)預(yù)測密度和實際密度，計算ME、RMSE、MAE、MASE和R2。計算如下：ME=0.029 074，RMSE=0.441 683，MAE=0.332 771，MASE=0.517 030，R2=0.907。結(jié)果表明預(yù)測精度較高，預(yù)測與監(jiān)測曲線擬合見圖4。

圖4 2019年1—12月時間序列預(yù)測和監(jiān)測結(jié)果擬合曲線

3 討論

蚊蟲侵害受到多種因素，如氣溫、濕度、氣壓等氣象因素以及植被、水網(wǎng)、土地利用模式等地理因素的影響。由于不斷進行的城市化進程和人員交通物流持續(xù)增強的社會化活動，各種社會性因素也在影響著蚊蟲侵害問題。人民對于日益增長的美好生活的要求，對于人居環(huán)境改善的迫切需求等，使得蚊蟲侵害顯然不僅僅是一個疾病預(yù)防的問題，而是逐漸演變成一個公眾關(guān)注的公共衛(wèi)生問題，因此要更加關(guān)注蚊蟲侵害的現(xiàn)狀和未來，更加注重蚊蟲控制的效果評價，更加注重蚊蟲和蚊媒病的相關(guān)性研究［10-13］。

蚊蟲侵害密度監(jiān)測是國內(nèi)開展較為成熟的橫斷面研究，方法較為系統(tǒng)和成熟，有長時間的工作積累和分析方法，但大部分都局限在橫斷面分析，如對蚊蟲的密度、種群、群落結(jié)構(gòu)和季節(jié)消長進行描述性研究，僅有少部分開展前瞻性的預(yù)測預(yù)警研究，方法較為單一。SARIMA 方法是ARIMA 模型中對于季節(jié)變化趨勢較為明顯的時間序列進行研究的方法，相對于簡單的AR 或者MA 模型，SARIMA 加入了季節(jié)性變化因素，使模型分析更為精準。有學者利用ARIMA 模型對三帶喙庫蚊或白紋伊蚊密度開展過類似研究，預(yù)測的擬合效果均較好，說明ARIMA 模型適用于開展蚊蟲侵害預(yù)測研究［14-16］。從本研究來看，SARIMA 對于本地區(qū)蚊蟲侵害的預(yù)測精度較高，決定系數(shù)R2達到了0.9 以上，預(yù)測的成功率大大提升，對提前進行蚊媒病防制動員和愛國衛(wèi)生運動具有前瞻性意義。

但是也要看到SARIMA 模型仍是一種傳統(tǒng)的線性時間序列模型，預(yù)測時僅僅考慮從歷史看未來，而沒有加入可能影響未來的各種因素（如上述所說的氣象環(huán)境、社會因素等），因此本研究預(yù)測的高峰值與實際高峰值有所差別，所以這是本研究之后要考慮開展的方向。如氣象因素對于蚊蟲密度也有滯后效應(yīng)，參照肖揚［13］的研究可以使用分布滯后非線性模型（DLNM）開展類似研究。對于時間序列中的非線性部分本文未涉及，參照國內(nèi)類似研究可以在ARIMA的基礎(chǔ)上加入支持向量機（SVM）或神經(jīng)網(wǎng)絡(luò)模型等人工智能技術(shù)［17-20］，使得預(yù)測擬合效果更加貼近真實值，這樣才能真正服務(wù)于現(xiàn)實公共衛(wèi)生工作需要。