高 翔，肖建華，王洪斌

（東北農業(yè)大學黑龍江省動物普通疾病防治重點實驗室， 哈爾濱 150030）

0 引 言

在養(yǎng)殖業(yè)中，家禽養(yǎng)殖受傳染性疾病影響較大。與發(fā)達國家相比，中國家禽養(yǎng)殖在疾病防控方面存在的一個問題是，從業(yè)者往往更依賴疫苗與抗生素進行傳染病的防控，忽視養(yǎng)殖管理、綜合防制的重要性。經驗證明，忽視日常管理，僅依靠抗生素與疫苗無法完全有效地控制家禽傳染性疾病的傳播流行[1]。疫苗雖然可以阻止動物發(fā)病，卻不能阻止在養(yǎng)殖場與養(yǎng)殖家禽種群內已存在的病毒的擴散與傳播[2]。同時，由于病原微生物不斷的變異，在面對新發(fā)傳染病時，疫苗的研發(fā)會有一個滯后期。在這種情況下，評估養(yǎng)殖場內外存在的影響疫病傳播的風險因素，針對性加強管理，提高動物抗病能力，降低動物接觸傳播病原的風險，是一種有效的防控措施。同時，抗生素殘留產生的耐藥問題已嚴重影響對相關疾病的治療?？股仉S著動物排泄物進入環(huán)境，在土壤和水體中慢慢富集，也最終對人類健康造成危害。對中國河流湖泊水樣進行抽樣調查后發(fā)現(xiàn)，在采集到的水體樣本中檢出的抗生素濃度普遍較高[3]。

綜上所述，國內外流行病學研究專家越來越重視利用現(xiàn)代信息技術加強動物養(yǎng)殖健康管理[4-5]。同時，傳染性疾病的傳播流行往往受到環(huán)境地理等因素，如氣象、植被、海拔等的影響[6]。為了及早發(fā)現(xiàn)傳染性疾病暴發(fā)流行的“苗頭”，國內外的學者利用各種疾病與風險因素數據，結合數學評估預測模型，開展對動物傳染性疾病風險因素評估與早期預警的研究，在實際應用中表現(xiàn)出了良好的性能[7-8]。

禽巴氏桿菌病(Avian pasteuretlosis)，俗稱禽霍亂，是一種在家禽中傳播的急性敗血性呼吸系統(tǒng)傳染病。長久以來，禽霍亂是中國《獸醫(yī)公報》中收錄發(fā)病數據最多的家禽傳染性疾病。由于該病發(fā)病快，致死率高，嚴重威脅著中國家禽養(yǎng)殖業(yè)的健康發(fā)展，中國將此病列入家禽傳染病控制二類名單[9-10]。通過對過往發(fā)病數據分析，發(fā)現(xiàn)禽霍亂在中國集中分布于南方高溫高濕氣候省份，并表現(xiàn)出季節(jié)性發(fā)病的特點，因此，氣象因素是影響禽霍亂發(fā)病的重要風險因素。

基于以上背景，該研究提出利用數學評估模型評估氣象數據對禽霍亂發(fā)病影響，利用網絡地理信息（Web GIS）技術建立禽霍亂監(jiān)測預警系統(tǒng)的思路。通過本文中的研究與系統(tǒng)開發(fā)，期望可以支持家禽養(yǎng)殖場疫病防控管理決策，降低經濟損失，為數學評估模型在各類動物疫病防控工作中的應用研究提供參考。

1 系統(tǒng)設計

1.1 系統(tǒng)結構設計

系統(tǒng)的建設采用 B/S架構，不僅能夠滿足用戶對傳統(tǒng)信息管理系統(tǒng)的需求，還可以進行數據的輸入、輸出、編輯、更改、查詢、計算等功能。而且，在地理信息技術的支持下，還可以進行數據的電子地圖可視化操作與空間分析功能。用戶在客戶端填寫的各類信息通過互聯(lián)網絡，傳輸到管理部門的后臺數據庫服務器，進行運算分析，得到的結果再通過網絡返回用戶。系統(tǒng)總體結構可以被抽象為用戶層、功能層、技術層及數據層 4個層次。

圖1 系統(tǒng)架構圖Fig.1 Architecture diagram of system

1.2 系統(tǒng)功能設計

系統(tǒng)功能大體上劃分為地圖基本操作、疫情數據上報、疫情數據查詢、疫情空間分析和風險預測 5大功能模塊。具體目標包括：1）地圖基本操作：包括地圖的放大、縮小與平移，鷹眼視圖等功能。2）疫情數據上報：疫情上報信息包括發(fā)病日期、發(fā)病地點、動物種類、動物來源、發(fā)病數量、報告人等。利用系統(tǒng)，使用者點擊地圖上緊急疫情發(fā)生的地點，彈出數據錄入對話框，錄入相關信息。信息錄入后，在地圖上生成特定標示的紅色點。數據錄入還包括對歷史疫情信息的收集整理，包括發(fā)病數、死亡數、捕殺數等。對于歷史疫情數據，還可以形成統(tǒng)計報表和統(tǒng)計圖形。3）疫情數據查詢：信息查詢功能包括圖查屬性功能與屬性查圖功能。圖查屬性功能允許使用者以感興趣的某一點為中心，通過設置查詢半徑，在生成的查詢區(qū)范圍內的所有信息都已列表的形式被顯示。這方便了使用者及時掌握疫情點周圍的信息；屬性查圖通過鍵入感興趣點的信息，查找其在地圖上的位置。這方便了使用者快速查找到所感興趣的疫情點。4）疫情空間分析：包括緩沖分析與路徑分析。緩沖分析為三重緩沖分析，對應滿足了疫情管理中需要劃分疫區(qū)、受影響區(qū)以及封鎖區(qū)的要求；路徑分析通過在地圖上設置起點及終點，由系統(tǒng)自動給出兩點間最短的路徑信息。方便使用者在應對緊急疫情時進行物質調配和人員的調動。5）風險預測：利用得到的各風險因素與疾病發(fā)病相關的風險系數，結合數據庫中的風險因素數據，計算疾病發(fā)生的風險值并生成風險地圖。地圖以顏色的深淺判斷疾病風險的大小。

1.3 系統(tǒng)數據庫設計

系統(tǒng)數據庫包括屬性數據庫、空間數據庫以及模型庫。數據庫支持系統(tǒng)實現(xiàn)對數據的管理、計算功能，是系統(tǒng)實現(xiàn)空間分析與預測的基礎。本研究采用的是ArcGIS 10.0內置數據庫來管理空間數據，以矢量地圖格式，按空間、地物屬性分類的依據分層存儲，地圖坐標為經緯度坐標。

圖2 監(jiān)測預警系統(tǒng)細節(jié)功能圖Fig.2 Function of monitoring early warning system

屬性數據庫主要包括養(yǎng)殖狀況數據以及歷史疫情數據等。養(yǎng)殖情況數據包括全國各省、市、自治區(qū)每年年底家禽存欄數、出欄數、各省的土地面積等數據，為顯示疫情信息等功能提供數據支持。疫情數據包括各省月家禽疫情數據，包括發(fā)病數、死亡數、捕殺數等數據，是建立評估預測模型的主要數據來源。

預測模型是整個決策支持系統(tǒng)的核心功能部分。模型庫主要是利用已有的專家知識和專家經驗，采用零膨脹負二項式回歸模型(zero-inflated negative binomial，ZINB)方法建立相關預警模型，通過編程技術建立系統(tǒng)模型。模型可調用數據庫中的數據及參數值，并將結果返回動態(tài)數據庫。

2 預測模型的建立

2.1 研究數據

該研究中，2006年1月至2016年12月中國各省禽霍亂每月發(fā)病數據來源于中國獸醫(yī)局網站下載的《獸醫(yī)公報》。將2006年1月至2015年12月間的數據資料用于預測研究，2016年1月至2016年12月的數據資料作為模型預測效果的檢驗數據。養(yǎng)殖數據來自于由國家統(tǒng)計局網站下載的2006―2016年《中國統(tǒng)計年鑒》。氣象數據來自于國家氣象局信息中心。

2.2 時間與季節(jié)性分析

該研究中，首先通過時間與季節(jié)性分析排除之前月份發(fā)病對本月發(fā)病的影響，以及季節(jié)因子對疾病發(fā)病的影響，為氣象因素風險作用的分析做準備。

時間序列分析方法通過判斷時間上前后兩點間數據的相關關系，探尋數據相關事件的發(fā)展規(guī)律。時間序列分析包括對數據隨機性與平穩(wěn)性的分析。

時間序列分析的一個前提是要求研究數據為非平穩(wěn)性序列。檢驗時間序列是否為平穩(wěn)性序列，目前常用的是增廣 ADF檢驗(augmented dickey-fuller test)，簡記為ADF檢驗[11]。ADF檢驗是單位根檢驗(unit root test)的一種，原假設為序列具有非平穩(wěn)性，顯著水平取值為0.01。當P＜0.01拒絕原假設。ADF檢驗任何一部分結果中出現(xiàn)P＞0.01的情況，都表明原假設成立。其計算公式如下所示：

式中μ為常數均值；β為回歸系數；εt為誤差項。

進行時間序列分析的另一個前提是作為研究對象的數據序列必須存在時間上的關聯(lián)性，純隨機序列沒有進行時間序列分析的價值[12]。理論上，純隨機序列數據的自相關系數（autocorrelation coefficients，AC）為 0。自相關系數可根據以下公式求得：

式中 n為序列觀測值的數目；ρk為樣本自相關系數；k(0＜k＜n)為延遲階數；xt為t時刻序列測量值。

該研究中，首先利用ADF方法判斷禽霍亂發(fā)病數據序列的平穩(wěn)性。然后，利用自相關與偏自相關系數判斷發(fā)病數據本身存在的自相關性。相關分析由SAS 9.2分析軟件完成。

該研究通過使用季節(jié)分解模型對2006年1月至2015年12月中國家禽傳染性疾病月發(fā)病數據進行分析，描述其季節(jié)分布特征與長期趨勢。通過季節(jié)分解模型，可以將原始發(fā)病數據分解得到季節(jié)因子序列（seasonal adjustment factors，SAF），通過SAF序列，可以判斷發(fā)病數據的季節(jié)性，當序列峰值大于 1時，發(fā)病數據具有季節(jié)性[13]。對于時間序列ARIMA（n, d, m），其季節(jié)分解模型的一般形式為：

式中e為時間序列ARIMA（n, d, m）的殘差；d為時間序列由非平穩(wěn)序列轉變?yōu)槠椒€(wěn)序列的d階差分；a為時間序列ARIMA（n, d, m）的白噪音，n與m分別為時間序列自相關系數最佳的階層和偏自相關系數最佳的階數，

2.3 風險分析

通過對禽霍亂發(fā)病數據時空分布特點的研究，發(fā)現(xiàn)禽霍亂發(fā)病會受到季節(jié)因子作用的影響。同時，上月發(fā)病會影響本月疾病的發(fā)生。因此，季節(jié)因子與上月發(fā)病數作為協(xié)變量納入模型分析。當地的家禽養(yǎng)殖量作為啞變量控制因養(yǎng)殖規(guī)模不同而對發(fā)病數據造成的影響。

為了檢驗氣象因子相互之間是否存在共線性。研究利用 Spearman 等級相關分析方法對氣象因子間的相關關系進行檢驗。Spearman 等級相關系數取值區(qū)間為0～1，一般認為相關系數大于 0.7時，兩變量間存在共線性，不能同時納入多變量的模型分析[14]。

對于存在共線性的 2個氣象因素，研究采用采用赤池信息量信息準則（akaike information criterion，AIC）對納入了不同變量的模型進行評估。根據相關準則，AIC的取值越小，模型的擬合預測結果越好[15-16]。

式中K是所擬合模型中參數的數量；L是對數似然值；n是觀測值數目。

本文中的研究使用互相關(cross-correlation)分析方法研究氣象因素對全國家禽傳染性疾病的滯后作用，其計算公式為[17]：

式中 Cov(yt,Xt-k)為延遲 k階的協(xié)方差，Var(yt)，Var(Xt-k)分別是相應變量的方差。

在傳染性疾病的研究中，由于疾病空間集聚以及季節(jié)性分布的特點，發(fā)病數的方差往往明顯偏離發(fā)病均數，不能滿足等離散的條件。對這些類型的傳染病數據進行分析，大大超過了傳統(tǒng)計數模型的分析能力。同時，在對發(fā)病數據此類計數資料進行研究的過程中，時常還會遇到在單位觀測時間段內，事件的發(fā)生次數為0的現(xiàn)象。傳統(tǒng)的計數模型對此類含有大量“0”的事件數據都不能有效的進行擬合。因此，Lambert于1992年提出了ZINB模型[18]。ZINB模型是一種將“零膨脹”事件的發(fā)生過程分作了零計數過程以及泊松計數過程 2個部分的混合概率分析方法。利用零計數過程可以解釋事件發(fā)展過程中出現(xiàn)大量“0”的原因，有效地彌補了傳統(tǒng)計數分析模型分析含有大量“0”的事件數據不足的問題，使得分析結果更加準確[19]。

ZINB模型形式如下：

式中φ為ZI參數；α為散度參數，E(Y)=(1-φ)/λ。

2.4 模型評估

該研究釆用組內相關系數（intraclass correlation coefficient，ICC）分析方法評估模型預測發(fā)病數與實際疾病發(fā)病數之間的一致性[20]。本文采用的ICC分析方法是經Shrout和Fleiss對傳統(tǒng)的組內相關系數分析方法理論作出改進，將其擴展到了評價重復量測數據之間相關性應用領域后的方法，組內相關系數被定義為被測量者間變異占總變異的比例[21]。改進后，ICC分析方法被逐漸應用于經濟學與醫(yī)學統(tǒng)計分析領域，組內相關系數的取值范圍在0～1，一般認為ICC大于0.75時預測值與實際值擬合良好[22]。其計算公式為：

式中MS處理為真實值與預測間的均方；MS區(qū)組代表不同月間病例數的均方；MS誤差為誤差的均方；n被觀察對象的例數；k為處理組數。

3 預測結果與預測系統(tǒng)實現(xiàn)

3.1 時間序列與季節(jié)性分析結果

由ADF檢驗結果可見，各檢驗方法均有P值大于0.01，不能拒絕無效假設，可認為序列為非平穩(wěn)性時間序列。

表1 禽霍亂發(fā)病數序列單位根檢驗結果Table 1 Results of ADF tests of Avian pasteuretlosis sequence

從自相關系數可以看出，除 1～4階自相關系數外，延遲5階及以后的自相關系數均在2倍的標準誤以內，沒有統(tǒng)計學意義。同時，由表2可以看出延遲1階的自相關系數及偏自相關系數對比其他階系數均明顯偏大，這說明了中國2006年1月至2015年12月禽霍亂病例具有延遲一階的自相關，即當月的禽霍亂發(fā)病數受上個月病例數的影響。

表2 禽霍亂發(fā)病數序列自相關與偏自相關系數Table 2 Autocorrelation coefficients(AC) and partial autocorrelation coefficients(PAC) of Avian pasteuretlosis sequence

通過禽霍亂各月份發(fā)病數與季節(jié)因子SAF的曲線可以看出，2006年10月至2015年12月中國禽霍亂發(fā)病多出現(xiàn)在秋冬季節(jié)（圖 3），具有明顯的季節(jié)性（季節(jié)因子峰值超過1）[13]。

圖3 禽霍亂各月份發(fā)病數Fig.3 Number of Avian pasteurellosis cases in every month

3.2 氣象因素風險分析

該研究將氣象因素的當月數據和之前 1～3月的數據作為自變量，將全國家禽傳染性疾病每月發(fā)病數作為響應變量，納入分析。模型納入前 1個月的家禽傳染性疾病月發(fā)病數用于控制疾病自相關；納入季節(jié)因子控制禽病疫情發(fā)生的季節(jié)性；納入各地區(qū)的禽養(yǎng)殖數量總數抵消養(yǎng)殖數量不同造成的差異。

圖4 禽霍亂的季節(jié)因子分析Fig.4 Seasonal factor analysis of Avian pasteurellosis

由于Spearman等級相關分析顯示，發(fā)現(xiàn)中國當月的降水量與月平均濕度以及月平均溫度，當月的降水量與滯后1個月的月最高溫度呈共線性，不能同時納入模型。因此，禽霍亂多因素研究有2個模型：I滯后本月的平均相對濕度+本月的平均溫度+本月的最高氣溫+本月的平均風速+上月的平均氣溫；II本月的總降水量+本月的最高氣溫+本月的平均風速+上月的平均氣溫。根據AIC準則，選擇模型I作為最優(yōu)模型。

利用回歸模型分析得到氣象因素對禽霍亂發(fā)病影響作用結果如表 1示。在計數過程中，本月最高氣溫、本月相對濕度與本月禽霍亂發(fā)生呈正相關，本月和上月的月平均氣溫也與其呈正相關；季節(jié)性因子與本月平均風速與禽霍亂發(fā)生呈負相關。在零計數過程中，上個月禽霍亂病例對零計數的發(fā)生概率呈負相關，即上個月的禽霍亂病例會降低本月不發(fā)生禽霍亂的概率。

表3 ZINB回歸模型結果Table 3 Results of multivariate analysis by using zero-inflated negative binomial model

3.3 模型評估結果

利用組內相關系數方法對預測值與中國2016年禽霍亂月實際發(fā)病數進行對比，結果顯示 ICC取值分別為0.802，說明模型預測值與實際病例數據有較好的一致性，預測值與觀測值之間的R2取值為0.7853。

圖5 2016年禽霍亂的月發(fā)病數預測值與觀測值擬合圖Fig.5 Predicted number of outbreaks versus observed number of monthly-reported outbreaks of Avia, 2016

3.4 系統(tǒng)實現(xiàn)

系統(tǒng)基于ArcGIS 10.0與Visual Studio.NET 2008開發(fā)平臺，以 C#. NET語言為基礎進行開發(fā)，數據庫采用Microsoft SQL Server 2005。開發(fā)完成的系統(tǒng)軟件可在Windows操作系統(tǒng)下運行，系統(tǒng)界面如圖6所示。系統(tǒng)要實現(xiàn)的主要目標包括：①數據管理：通過疫情上報的方式對家禽歷史疫情信息包括發(fā)病數、死亡數、捕殺數等進行收集；通過系統(tǒng)自身建立的環(huán)境數據專題庫對影響家禽傳染性疾病發(fā)病的氣象因素進行整理；②疾病預警：利用ZINB模型分析得到的各氣象風險因素與疾病發(fā)病相關的風險系數，計算疾病發(fā)生的風險值并生成風險地圖；③應急管理：通過緩沖分析、路徑分析等多種空間分析功能，方便使用者在應對緊急疫情時制定物質調配和人員調動等應急管理策略；④用戶可隨時查詢任意時間段內的監(jiān)測記錄，并生成相對應的報告圖片。系統(tǒng)部分功能如圖7所示。

4 討論

傳統(tǒng)研究中，利用邏輯回歸分析或相關性分析模型揭示氣象因素對傳染性疾病發(fā)生傳播的影響作用，在進行此類研究時，需要預先假設氣象因素與傳染病發(fā)病之間滿足線性或對數線性相關的關系。然而，在實際中，氣象因素與疾病發(fā)病之間的關系往往表現(xiàn)非線性特點。例如，在一定的范圍內，病原微生物的致病能力隨環(huán)境溫度的上升而增加。但當氣溫超過一定范圍時，環(huán)境溫度升高反而抑制病毒的感染能力[23]。此外，對于傳染性疾病發(fā)生發(fā)展，氣象因素還存在滯后性影響。例如，當氣象變化超過生物體適宜的生存范圍時，機體的免疫能力會受損下降。但這種變化并不是立即發(fā)生的，生物體的內環(huán)境具有一定的“穩(wěn)態(tài)”，對于外界條件的變化往往表現(xiàn)一定的抵抗力[24]。ZINB模型是一種改進后的泊松分析方法，它不僅能夠很好地解釋氣象因素對傳染性疾病的非線性影響，結合互相關分析，對氣象因子的滯后性作出評估。而且，對于疫情數據集中因為某些月份沒有疾病發(fā)生而存在的大量“零數據”， ZINB模型也能有效消除其對分析結果的影響，大幅提髙預測結果的準確度。因此，ZINB模型被越來越多的應用于氣象對疾病傳播流行影響的分析中[25]。

ZINB模型分析顯示禽霍亂發(fā)生與氣溫、濕度呈正相關，與風速呈負相關。氣象因素與氣候變化對傳染病發(fā)生的影響表現(xiàn)在各個方面，包括影響宿主的抵抗力和易感性等[26]。一般而言，機體內環(huán)境對于外界環(huán)境的變化存在一定的抵抗力，維持著一定的“穩(wěn)態(tài)”。但是當動物長期處于不適宜的外界環(huán)境中時，內環(huán)境的“穩(wěn)態(tài)”會受到破壞，引起免疫力的下降。Garren等[27]于1956年的研究表明，長期處于高溫、高濕的環(huán)境下會導致機體的免疫能力受到抑制。解剖后發(fā)現(xiàn)，長期在高溫、高濕環(huán)境下飼養(yǎng)的雞的胸腺、法氏囊等免疫器官發(fā)生壞死萎縮、淋巴液相較于正常的雞明顯減少，抗病能力顯著下降。

除了高溫、高濕的環(huán)境影響外，家禽養(yǎng)殖場內及其周邊環(huán)境的空氣質量也會造成家禽免疫力的下降。這是由于健康的動物體內或體外也會帶有一些致病微生物，這些微生物隨唾液或脫落的皮屑在養(yǎng)殖場內或養(yǎng)殖場周邊的空氣中飄蕩，經過長時間的風化，與空氣中的水、顆粒塵埃等進一步融合形成氣溶膠。在氣溶膠內的微生物往往可以在很長時間內保持感染能力[28]。同時，家禽養(yǎng)殖場內及周邊的空氣中，硫化氫與氨氣的濃度往往很高。這些氣體進入呼吸道，會引起呼吸道免疫系統(tǒng)的纖毛-黏液系統(tǒng)的功能受損，降低呼吸道的抗感染能力[29]。

綜上所述，高氣溫與濕度引起禽霍亂發(fā)病數的增加可能是因為高的氣溫與濕度會造成動物體免疫機能的損傷與下降。由于機體對外界環(huán)境的變動具有一定的抵抗力，因此氣溫對禽霍亂發(fā)病數的影響也表現(xiàn)出一定的滯后作用。較高的風速可能在一定程度上抵抗高溫、高濕的影響，同時，風速較大時可能會降低養(yǎng)殖場環(huán)境空氣中氨氣、粉塵及有害氣體的濃度。因此，風速對禽霍亂的發(fā)生呈現(xiàn)一定的負相關作用。

5 結論

本文通過將季節(jié)性分析、時間序列分析方法和零膨脹負二項式回歸模型技術相結合，對中國禽霍亂發(fā)病的分布特征進行了描述，對氣象因素對中國禽霍亂發(fā)病的影響作用進行了評估。結果表明：1）禽霍亂發(fā)生與氣溫、濕度呈正相關，與風速呈負相關。同時，氣象因素影響的滯后作用可能反映了動物機體對外界環(huán)境變動具有一定的抵抗能力；2）基于零膨脹負二項式回歸模型的分析結果，利用氣象數據對禽霍亂的發(fā)病作出預測，與實際發(fā)病數比較，準確率可達 80.2%，該模型的分析結果能夠較好地反映相關氣象因子對中國禽霍亂發(fā)病的影響作用。

家禽傳染性疾病的傳播過程復雜，受多種因素影響，同時，中國家禽養(yǎng)殖模式多種多樣，養(yǎng)殖水平高低不一。條件所限，本文只選擇對氣象因素對中國禽霍亂發(fā)病的影響進行了分析。如何將經濟、交通乃至候鳥遷徙等更多風險因素納入分析，進一步提高預測的準確度，是下一步研究工作的重點。