非洲豬瘟舍內(nèi)傳播元胞自動機(jī)仿真研究

楊婷婷， 薄興野， 田晶， 邵家源， 劉振宇*

(1.山西農(nóng)業(yè)大學(xué)信息科學(xué)與工程學(xué)院， 太谷 030801； 2.山西農(nóng)業(yè)大學(xué)農(nóng)業(yè)工程學(xué)院， 太谷 030801)

非洲豬瘟病毒在豬只之間的傳播，主要與畜禽舍環(huán)境的溫度和濕度等相關(guān)因素有關(guān)。爆發(fā)地大多分布于北方[1]。每年年初是非洲豬瘟爆發(fā)的高峰期，爆發(fā)特點(diǎn)具有聚集性[2]。畜禽舍內(nèi)豬只直接接觸是其傳播的主要途徑，具體包括豬只之間進(jìn)行口鼻接觸、共同飲食、空氣傳播等，空氣傳播一般不超過2 m，豬只間保持安全距離，可減少感染[3]。急性非洲豬瘟通常會集中發(fā)生，且致死率極高，從患病到死亡的時間非常的短，往往不超過2 d，同時還未研發(fā)出可治療的藥物和有效的疫苗[4]。

針對傳染病的研究，目前大多都是使用單一的傳染病模型來進(jìn)行分析。唐雄[5]從SIR(susceptible-infective-removed)模型的動力學(xué)行為作為研究，得出了傳染病在傳播過程中得以控制的條件。常見的傳染病模型，建立在微分方程之上，雖然計算方便，但存在很大的缺陷[6-10]。針對傳染病在畜禽舍內(nèi)真實(shí)傳播過程中，一些隨機(jī)發(fā)生的事件，不能進(jìn)行很好的處理。復(fù)雜的計算與現(xiàn)實(shí)畜禽舍內(nèi)，豬只的接觸產(chǎn)生的隨機(jī)事件，不能很好地結(jié)合。因此單純使用傳染病模型模擬傳染病的傳播過程存在明顯不足，忽略了很多局部存在的個體特性[11]。

隨著技術(shù)的發(fā)展，大量的仿真與數(shù)學(xué)計算為傳染病(如非洲豬瘟)的模擬仿真提供了有利條件。目前，元胞自動機(jī)已被廣泛應(yīng)用于人群擴(kuò)散[12]、魚群行為預(yù)測[13]、線路規(guī)劃[14]、新冠病毒傳播模型的建立[15]，驗(yàn)證了元胞自動機(jī)在行為移動預(yù)測方面的真實(shí)性。

非洲豬瘟自傳入中國以來，對養(yǎng)豬業(yè)和食品安全等相關(guān)行業(yè)，產(chǎn)生了嚴(yán)重的不良影響[16]。因此，研究非洲豬瘟的傳染規(guī)律，進(jìn)行仿真研究，當(dāng)非洲豬瘟在養(yǎng)殖環(huán)境中爆發(fā)之后，能夠根據(jù)其預(yù)測結(jié)果，對病毒進(jìn)行有效的防控?，F(xiàn)使用元胞自動機(jī)，考慮非洲豬瘟傳染病的隨機(jī)性，將元胞自動機(jī)與傳統(tǒng)SI(susceptible-infected)模型相結(jié)合，對非洲豬瘟傳染病的傳播情況進(jìn)行仿真分析，以期為相關(guān)傳染病的仿真研究提供理論指導(dǎo)。

1 研究方法

1.1 非洲豬瘟傳染模型選擇

非洲豬瘟在密閉豬舍內(nèi)的傳播特征為，豬只個體數(shù)量少，傳染目標(biāo)明確，與預(yù)期的效果相符。非洲豬瘟傳染病的特性：發(fā)病時間短，患病到死亡時間極端，而且暫無有效疫苗的研發(fā)，因此不存在治愈者。傳染病模型包括：SI、SIR、SEIR(susceptible-exposed-infective-removed)等，結(jié)合非洲豬瘟傳染病的特性，因?yàn)榉侵挢i瘟發(fā)病時間極端，所以不考慮E潛伏豬只的可能，即接觸過感染豬只但暫無傳染性的豬，而是全部都用易感豬代替，而SEIR多用于惡性傳染病，且潛伏期對于整個傳染周期較長。而對于致死率相對高的的傳染病，大多感染后就會死亡，不存在治愈的情況，R治愈豬只假設(shè)為不存在。因此針對非洲豬瘟選用SI模型，相比較更為合適。該模型是一個單向模型，隨著時間的推移，易感豬只數(shù)會逐漸向感染豬只數(shù)輸入，最終易感豬只數(shù)都會降為 0，所有均變成患病豬。這也是SI模型的局限性，因此該模型不適用于傳播周期長和有疫苗研發(fā)成功的的傳染病，但適合與非洲豬瘟傳染病的研究。因此，在其理論研究基礎(chǔ)上，基于動力學(xué)理論，建立非洲豬瘟傳播的動力學(xué)SI模型，并計算模型的解。通過非洲豬瘟傳染病模型的構(gòu)建與仿真，發(fā)現(xiàn)對控制和管理非洲豬瘟的傳播起了重要作用。

1.2 模型的構(gòu)建

對SI模型建立的基本假設(shè)。首先，假設(shè)已知的數(shù)據(jù)具有真實(shí)性，對于豬只的出生、死亡等情況不進(jìn)行考慮，設(shè)定該豬舍內(nèi)的豬只總數(shù)為Ntotal。其次，將該模型簡單的分為易感豬只S和感染豬I，其模型傳播機(jī)理如圖1所示，SI模型的感染規(guī)則[17]。最后，假設(shè)豬只不可流動，即無豬只流入和流出該豬舍，因?yàn)椴恍枰紤]了感染畜禽治愈后獲得免疫能力的情況和接種疫苗后的治愈行為。

圖1 SI模型Fig.1 SI model

圖1中，S為易感類，該類豬只沒有免疫能力，當(dāng)符合感染情況，過度接觸時，有一定的概率染上該豬瘟；I為感染類，該類豬只以及為患病豬，具有一定的概率把該病毒傳染給易感豬只；β為易感豬只被患病豬感染的概率。

假設(shè)總畜禽數(shù)不變，比例系數(shù)不變。在初始狀態(tài)時，設(shè)置一個或者多個為感染豬只。在單位時間內(nèi)，當(dāng)易感豬只和感染豬只發(fā)生接觸行為，易感豬只就會以概率為β的可能性轉(zhuǎn)變成感染豬只。t時刻易感豬只占總數(shù)的比例為S(t)，感染豬只占總數(shù)的比例為I(t)，兩者的關(guān)系可表示為

S(t)+I(t)=1

(1)

該模型是一個單向模型，隨著時間的推移，易感豬只數(shù)會逐漸向感染豬只數(shù)輸入，最終易感豬只數(shù)都會降為0，所有均變成患病豬。

1.3 異質(zhì)性

在現(xiàn)實(shí)生活中，豬只感染的概率與環(huán)境、群體數(shù)量、季節(jié)變化等相關(guān)因素之間存在一定的聯(lián)系。因此需要對不同的元胞，設(shè)置參數(shù)不同的感染系數(shù)，從而使模擬的畜禽舍傳染過程考慮到了群體的異質(zhì)性[18]。根據(jù)不同豬只個體體質(zhì)的考慮，設(shè)計抵抗能力測算值Ri,j；對患病后豬只的運(yùn)動情況的分析，設(shè)計不同的傳染能力fi,j，以上即為個體間存在異質(zhì)性[19]。

典型的SI模型的微分方程在計算的過程中，并沒有考慮豬只的個體異質(zhì)性，每個豬只都有相同的概率感染非洲豬瘟病毒。結(jié)合豬只具有異質(zhì)性這一特點(diǎn)，針對不同因素進(jìn)行全面的考慮，引入距離影響因子[式(2)]，進(jìn)而推導(dǎo)計算出豬只間通過有效接觸而引發(fā)的傳染概率Pi,j[式(5)]。

(2)

(3)

(4)

式中:Dci,jci′j′為距離影響因子；i和j為豬只移動前的坐標(biāo);i′和j′為豬只移動后的坐標(biāo);ci,j為感染概率計算豬只即中心元胞；ci′,j′為豬只8個方向的任一鄰居豬只；Nci,j為中心豬只的摩爾鄰居內(nèi)所有鄰居豬只的集合。

綜上所述，定義為中心元胞的豬只ci,j的有效接觸感染概率為Pi,j，Pi,j值越大，說明越容易感染非洲豬瘟。

(5)

式(5)中：fi,j為個體的感染能力；Ri,j為個體的抵抗能力。

從式(5)可以看出，在具有傳染病的豬舍內(nèi)，不同的豬只感染病毒的概率與距離、具體病毒的傳染力以及個體本身的抵抗能力3個值相關(guān)。

1.4 演化規(guī)則

演化規(guī)則為元胞自動機(jī)進(jìn)行演化分析的核心[20]，狀態(tài)轉(zhuǎn)移函數(shù)f可表示為

(6)

式(6)中：

系統(tǒng)時間記錄的基本單位為天數(shù)，記錄元胞此刻的存在狀態(tài)，在演化過程中，關(guān)于具體元胞在系統(tǒng)內(nèi)的移動性描述也是一個模型成功與否的關(guān)鍵因素，結(jié)合對豬舍內(nèi)豬只活動范圍的觀察，在忽略外界干預(yù)的情況下，豬舍內(nèi)豬只的移動是隨機(jī)的，因此對豬只移動規(guī)則進(jìn)行設(shè)計，因此本模型采用的是隨機(jī)移動步長的方式。

1.4.1 豬只移動與感知范圍

元胞關(guān)系為Moore型鄰元關(guān)系，圖2所示的摩爾鄰居即為豬舍內(nèi)的豬只移動范圍。

圖2 元胞(i,j)的摩爾鄰居Fig.2 Moore neighbor of cell (i, j)

系統(tǒng)將所模擬的豬舍進(jìn)行網(wǎng)格劃分，即劃分為M×N個相同大小的網(wǎng)格，網(wǎng)格元胞拼接即為模擬的豬舍的二維空間。設(shè)定每個元胞都屬于，設(shè)施(即畜禽舍的墻面和圍欄)、個體(豬只個體本身)、空白區(qū)域(豬舍可活動區(qū)域)，這3種屬性特征的其中一種。因此，對于任意位置(xi,yi)元胞表示的豬只，元胞狀態(tài)集合Scell的表現(xiàn)形式為

Scell=(xi,yi,Ui,j,ui,j,vi,j)

(7)

式(7)中：Ui,j、ui,j、vi,j分別為豬只在舍內(nèi)、向x軸移動時的速度、y軸移動時的速度，m/s。

1.4.2 模型移動行為假設(shè)

(1)避開豬舍內(nèi)設(shè)施移動。除空白活動區(qū)域外，墻體圍欄等設(shè)施以及豬只個體，都為不可活動的網(wǎng)格，因此豬只發(fā)現(xiàn)到會與墻體等發(fā)生碰撞的時候，會自覺地轉(zhuǎn)換移動方向。

(2)空白區(qū)域隨意移動行為。該模型不考慮豬只的行為的其他故意干擾等影響因素，因此豬只在豬舍內(nèi)的運(yùn)動為隨機(jī)活動，摩爾鄰居的各個方向都有一定的概率，而且速度等因素也不受控制

(3)優(yōu)先向目視范圍移動。一般豬比較習(xí)慣向自己目光所視方向移動，因此豬具有相當(dāng)大的概率向可以看到的范圍為移動，而非背后。

1.4.3 豬只移動規(guī)則

豬只在豬舍內(nèi)的移動規(guī)則，通過計算向不同的方向的權(quán)重值，進(jìn)而得到移動的概率，確定Agent下一步移動的網(wǎng)格[21]。具體的移動方向如圖3所示。

Agent向K方向移動的概率計算公式為

(8)

1～8表示方向圖3 Agent運(yùn)動方向Fig.3 Agent in the direction of motion

豬只不同方向網(wǎng)格的引力概率的計算公式為

(9)

(1)避開豬舍內(nèi)設(shè)施移動的表達(dá)式為

(10)

式(10)中：Kb為放置有設(shè)施的方向。

(2)空白區(qū)域隨意移動行為，模型中表示為豬只擁有向各個方向運(yùn)動的概率，可表示為

(11)

式(11)中：w為可能的移動方向；αran為w方向的吸引權(quán)重。

(3)優(yōu)先向目視范圍移動，豬具有相當(dāng)大的概率向可以看到的范圍為移動，其表達(dá)式為

(12)

式(12)中：Ksee為豬只的前方可視范圍的網(wǎng)格方向；αsee為可視范圍的網(wǎng)格吸引力權(quán)重值。

1.5 算法過程描述

步驟1使生成的豬只元胞，離散的處于代表豬舍的元胞空間內(nèi)。首先設(shè)定所有的豬只元胞狀態(tài)為0。然后隨機(jī)在這些元胞中，設(shè)定一定數(shù)量的感染豬只，這些元胞的狀態(tài)為1，初始設(shè)定數(shù)目不易過大，否則不能很好的模擬非洲豬瘟在畜禽舍內(nèi)的傳播過程。

步驟2設(shè)定元胞在空間中的最大移動步長為dmax。從系統(tǒng)模擬的迭代天數(shù)d=0時刻開始，每更新一次，都需要對系統(tǒng)中所有元胞的狀態(tài)進(jìn)行掃描與統(tǒng)計。同時，隨機(jī)選擇一定的元胞模擬其在畜禽舍內(nèi)的具體移動，標(biāo)記選擇的元胞坐標(biāo)為(i,j)，然后根據(jù)移動規(guī)則和各個方向的吸引權(quán)重因子，計算向周圍各個方向移動的最大值，得出移動的具體方向，即移動距離(di,dj)進(jìn)行移動，交換后坐標(biāo)為(i,j)和(i+di,j+dj)的兩個元胞，表示元胞發(fā)生了一次移動。

步驟3獲取每個時刻豬只元胞的狀態(tài)，然后判斷其是否為易感元胞，如果為易感元胞，再判斷其周圍是否存在感染元胞，并計算該易感元胞的具體感染概率值，再與規(guī)定的感染閾值X進(jìn)行比較，如果X

步驟4如果步驟3中的被選元胞不是易感元胞，判斷是否為感染元胞。如果是，則不進(jìn)行任何的改變，該元胞的狀態(tài)值保持不變，并且一直具有傳染能力。

步驟5對元胞空間內(nèi)所有的元胞進(jìn)行掃描，隨著系統(tǒng)的演化，不斷地更新每個元胞的狀態(tài)后，判斷元胞自動機(jī)是否達(dá)到特定的狀態(tài)，即穩(wěn)定狀態(tài)，當(dāng)所有的元胞都變?yōu)楦腥驹?，仿真過程正式結(jié)束，否則跳轉(zhuǎn)到步驟2重新選擇元胞。

2 元胞自動機(jī)仿真設(shè)計

2.1 物理模型的設(shè)計

按照豬場設(shè)計相關(guān)規(guī)范，來設(shè)計豬舍內(nèi)非洲豬瘟的傳播模型。實(shí)驗(yàn)地點(diǎn)為晉中市太谷區(qū)范村鎮(zhèn)北里村，東方希望養(yǎng)豬場，研究豬舍為育肥豬舍。模型的具體尺寸為：豬場長20 m、寬5 m，由10個大小相同的隔間組成，單個隔間的室寬為2 m、長為 4 m。不同隔間內(nèi)豬只的具體數(shù)量為6～10頭。豬舍的3D幾何模型如圖4所示。

圖4 豬場幾何模型Fig.4 Geometric model of pig farm

仿真研究嚴(yán)格按照豬場物理模型進(jìn)行模擬，通過對實(shí)地豬場的豬只飼養(yǎng)數(shù)量、活動情況、接觸數(shù)量等進(jìn)行實(shí)地的調(diào)查與試驗(yàn)進(jìn)行設(shè)計。

2.2 數(shù)據(jù)假設(shè)

不同傳染強(qiáng)度的非洲豬瘟傳染病，具體的元胞機(jī)設(shè)計也有所不同。首先，針對不同系統(tǒng)區(qū)域，設(shè)定不同元胞個數(shù)，來模擬養(yǎng)殖密度不同的豬舍隔間，方便觀察非洲豬瘟在豬舍內(nèi)的不同密度的豬舍的感染情況。其次，根據(jù)分析要求，系統(tǒng)初始患病豬只數(shù)為隨機(jī)生成的極少個數(shù)。然后演化系統(tǒng)內(nèi)的整體環(huán)境為豬只個體位于該離散系統(tǒng)的不同位置，并且按照設(shè)定的離散時間進(jìn)行輪次，在規(guī)定范圍內(nèi)按照設(shè)定的移動規(guī)則進(jìn)行移動。具體設(shè)計用四元組表示豬舍系統(tǒng)，即

A=(Ld,Scell,Nmoore,F)

(13)

式(13)中：A為模擬的整個豬舍系統(tǒng)；Ld為二維元胞空間網(wǎng)格的規(guī)模，設(shè)定為100×100的二維圖,假設(shè)的元胞個數(shù)為10 000個。但是并不是所有的元胞都表示為豬舍內(nèi)的豬只個體；Scell為元胞狀態(tài)集合；Nmoore為豬只的摩爾鄰居；F為豬只的感染和移動規(guī)則，即豬只狀態(tài)轉(zhuǎn)換函數(shù)。

假設(shè)該系統(tǒng)為不受其他外界因素影響的封閉的畜禽舍，不設(shè)置治愈率，感染后即為患病豬。當(dāng)健康豬只在系統(tǒng)移動過程中到達(dá)了新的位置后，遇到了另一只患病畜豬，不同的豬只傳染病的傳染性不同，根據(jù)豬只個體的異質(zhì)性，計算不同的感染概率。每只畜禽都有概率被傳染，或者去傳染臨近的其他鄰居元胞。

3 結(jié)果與討論

通過對比，分析元胞自動機(jī)對非洲豬瘟傳播的仿真模擬數(shù)據(jù)與原始SI模型傳播數(shù)據(jù)的擬合性。圖5為使用元胞自動機(jī)對非洲豬瘟自由傳播情況下的仿真結(jié)果。其中，可視圖的x、y坐標(biāo)大小分別設(shè)置為100、100，模擬整體豬場，小隔間的x、y大小分別為40、40,數(shù)量為10個，走廊大小設(shè)置為20、100，長度與實(shí)現(xiàn)選擇的豬場大小的長度比為 10∶100，寬度為5∶100，與實(shí)際豬場構(gòu)造擬合度較高；小隔間內(nèi)自動生產(chǎn)的豬只數(shù)目為2～11個，符合實(shí)際豬場平均隔間豬只數(shù)2～10只的數(shù)量。通過觀察可以發(fā)現(xiàn)，豬只每迭代一次的移動范圍在隔間內(nèi)，與現(xiàn)實(shí)豬只的分欄飼養(yǎng)相符。結(jié)果表明：元胞自動機(jī)仿真研究與現(xiàn)場豬舍較吻合，符合試驗(yàn)要求。

非洲豬瘟爆發(fā)初期，即模擬的迭代天數(shù)為1 d時[圖5(a)]，豬舍內(nèi)感染非洲豬瘟的豬只為4頭，感染率約為5%，此結(jié)果與SI模型模擬的第1～7.8天的感染數(shù)據(jù)(圖6)一致。

D為模擬的迭代天數(shù)，d圖5 自由傳播情況仿真結(jié)果Fig.5 Simulation results for the free propagation case

當(dāng)模擬的迭代天數(shù)為10 d時[圖5(b)]，此時的感染豬只為21頭，當(dāng)模擬的迭代天數(shù)為15 d時[圖5(c)]，此時的感染豬只高達(dá)69頭。在模擬的天數(shù)為10～15 d的時期增長率最高，與SI模型I=Ntotal/2時(圖6)，即感染豬只數(shù)目到達(dá)半左右時的，增長速率基本一致，都在較短的時間實(shí)現(xiàn)了大量的傳播。

當(dāng)模擬的迭代天數(shù)為50 d時[圖5(d)]，豬舍內(nèi)非洲豬瘟的感染率接近99%，只有一頭未染病豬只，其余基本全部感染，但相比模擬的迭代天數(shù)為15 d時的感染情況，增長趨于緩慢，此結(jié)果與SI模型(圖6)，第20天后數(shù)據(jù)趨于平穩(wěn)的結(jié)果一致。

通過對比兩種模型的仿真結(jié)果(圖7)可知，當(dāng)易感豬只全部轉(zhuǎn)化為感染豬只時SI模型d=13.8,元胞自動機(jī)迭代為d=14.3。擬合度為0.96。隨著時間的推移，感染豬只增長速率為：先快速上升再緩慢下降。SI模型有兩個主要結(jié)論:①指數(shù)增長率r=βNtotal，其中，Ntotal為豬舍內(nèi)養(yǎng)殖總數(shù)，即增長率與養(yǎng)殖總數(shù)成正比；②在I=Ntotal/2 時，I增加得最快，之后趨于平穩(wěn)。

圖6 SI模型傳播數(shù)據(jù)Fig.6 SI model propagation data

圖7 SI模型與CA模型傳播數(shù)據(jù)對比擬合Fig.7 SI model compared the propagation data with the CA model

分析可知，元胞自動機(jī)對非洲豬瘟傳播的仿真模擬與SI模型的數(shù)據(jù)能夠?qū)崿F(xiàn)較好的擬合，并且精確度優(yōu)于傳統(tǒng)的SI模型，可以更直觀的看到具體的傳播情況。

傳統(tǒng)的SI模型，只是對于非洲豬瘟爆發(fā)后感染豬只數(shù)量進(jìn)行預(yù)測計算，使用元胞自動機(jī)可以更好的對演化過程進(jìn)行可視化的顯示，通過迭代預(yù)測豬只個體在豬舍內(nèi)的活動，以及模擬非洲豬瘟隨著豬只的移動的傳播情況，對傳播過程與結(jié)果進(jìn)行預(yù)測。

4 結(jié)論

通過設(shè)計在封閉畜禽舍系統(tǒng)中，豬只感染非洲豬瘟傳染性疾病后傳播問題的仿真分析，得出如下結(jié)論。

(1)環(huán)境因素和豬個體體質(zhì)的差異，使每只豬被感染的概率是不同的，因此引入了距離影響因子。通過計算可知，易感豬只感染的概率與個體的感染能力、個體的抵抗能力及距離因子有關(guān)。

(2)考慮豬舍內(nèi)豬只的移動規(guī)則，建立的模型能夠模擬，避開豬舍內(nèi)設(shè)施移動、空白區(qū)域隨意移動行為和目視范圍移動3種行為模式，模擬結(jié)果能夠體現(xiàn)豬的運(yùn)動特征具有隨機(jī)性。

(3)進(jìn)行元胞自動機(jī)仿真試驗(yàn)，仿真結(jié)果與SI模型及實(shí)際情況下的傳播數(shù)據(jù)，進(jìn)行對比擬合，擬合度為0.96。結(jié)果表明：元胞自動機(jī)對非洲豬瘟的仿真結(jié)果與實(shí)際情況相符，且模擬精度更高。

(4)與傳統(tǒng)的傳染病SI模型相比，使用元胞自動機(jī)具有更好的預(yù)測性，具體包括可以更好的觀察非洲豬瘟在豬舍的具體傳播情況，以及對非洲豬瘟傳染病的演化過程可視化?？筛鶕?jù)仿真結(jié)果，更好的分析預(yù)測的傳播數(shù)據(jù)，進(jìn)而選擇適當(dāng)?shù)目刂撇呗詫?shí)現(xiàn)非洲豬瘟的有效控制。