陽成 王建波 許小可 杜占瑋
摘 要:流行病的傳播會對整個人類社會構成巨大威脅,因此迅速識別傳播源并及時采取控制措施至關重要。然而,由于流行病傳播過程具有多樣性、信息不確定性等因素,使得快速準確識別傳播源成為一項挑戰(zhàn)。結合反向感染算法、復合種群網絡模型以及馬爾可夫鏈理論,提出了一個在復合種群網絡中識別傳播源的新算法。該算法首先利用馬爾可夫鏈來初步估計子種群被感染的時間,被感染子種群根據感染時間獲得自己的身份信息,然后遍歷所有獲得感染子種群身份信息的子種群,將收集到的感染子種群身份信息傳播給其所有鄰居,最后根據獲得所有感染子種群身份信息的時間順序推斷出復合種群網絡的傳播源。在真實的航空網和人造復合種群網絡上進行大量仿真實驗,發(fā)現(xiàn)無論在已知全部感染快照還是部分感染快照的情況下,該算法與其他傳播溯源算法相比,識別傳播源的準確性都有顯著提升。該算法非常適合用于航空網這類復合種群網絡,對現(xiàn)實世界中的流行病傳播溯源和控制也具有參考意義。
關鍵詞:復合種群網絡; 傳播溯源算法; 反向感染; 計算機仿真
中圖分類號:TP391.9?? 文獻標志碼:A
文章編號:1001-3695(2023)09-019-2681-07
doi:10.19734/j.issn.1001-3695.2023.02.0034
RIA-based tracing algorithm for propagation of metapopulation networks
Yang Cheng1, Wang Jianbo1,2, Xu Xiaoke3,4, Du Zhanwei2
(1.School of Computer Science, Southwest Petroleum University, Chengdu 610500, China; 2.School of Public Health, University of Hong Kong,Hong Kong 999077, China; 3.Center for Computational Communication Research, Beijing Normal University, Zhuhai Guangdong 519085, China; 4.School of Journalism & Communication, Beijing Normal University, Beijing 100875, China)
Abstract:The spread of epidemics poses a significant threat to the entire human community. Therefore, it is critical to identify the sources of transmission quickly and take timely control measures. However, the diversity of epidemic transmission processes and information uncertainty makes it challenging to identify the sources of transmission quickly and accurately. This paper proposed a new algorithm for identifying transmission sources in a metapopulation network by combining the reverse infection algorithm and Markov chain theory. The algorithm firstly used a Markov chain to initially estimate the time when a subpopulation was infected, and the infected subpopulation obtained its own identity information based on the infection time. Then, it traversed all subpopulations that obtained the identity information of the infected subpopulation and spreaded the collected identity information of the infected subpopulation to all its neighbors. Finally, the spreading source of the metapopulation network could be inferred based on the temporal order in which all the identity information of the infected subpopulation was obtained. Simulation experiments conducted on real airports networks and artificial networks show that the accuracy of this algorithm is significantly improved compared to other algorithms, regardless of whether all or partial of the infection snapshots are known. This algorithm is well-suited for metapopulation networks such as aviation networks and is also useful for real-world epidemic transmission tracing and control.
Key words:metapopulation network; propagation traceability algorithm; reverse infection algorithm; computer simulation
0 引言
非典型性肺炎[1] 、甲型H1N1流感[2]、新冠肺炎[3]等呼吸道傳染病具有高致死率和傳染性,嚴重威脅人類健康,成為了制約人類社會發(fā)展的阻礙之一[4]。在互聯(lián)網上,電腦病毒[5,6]的存在也給人們帶來了許多負面的影響。2003年,出現(xiàn)了一種被稱為“2003蠕蟲王病毒”的病毒,當電腦感染該蠕蟲病毒后,網絡帶寬被大量占用,導致網絡癱瘓。這種病毒在亞洲、美洲、澳大利亞等地迅速傳播,造成了全球性的網絡災害。2006年,一種被稱為“熊貓燒香”的病毒在互聯(lián)網上傳播,上百萬的用戶受到影響。同時,隨著人們交流的便捷,謠言[7~9]也很容易流行開來,影響了人們的生活。傳播溯源[10,11]旨在基于感染快照、感染者流動信息等來識別傳播源頭的節(jié)點(例如,電腦病毒、謠言和流行病的傳播[12]),被廣泛地應用在流行病源頭推測、互聯(lián)網病毒源識別、社交網絡中的謠言源追蹤等方面。
傳播溯源問題包括流行病在人群和個體網絡中的時空傳播溯源,計算機網絡上的病毒溯源,以及社交網絡上的謠言溯源。對于傳播溯源問題的研究挑戰(zhàn)主要來自以下幾個方面:
a)網絡結構復雜多樣。從小世界網絡[13]、無標度網絡[14,15]再到多層網絡[16,17]、高階網絡[18,19],流行病或信息傳播的載體越來越復雜。在多尺度視角下,單個個體作為節(jié)點組成的接觸網絡可以傳播流行?。欢鄠€個體組成的城市作為節(jié)點組成的航空網也可以傳播流行病。網絡結構的復雜性,使得傳播溯源往往需要對不同類型的網絡設計不同的識別算法。
b)流行病類型多樣。在流行病傳播過程中,個體可以處于最常見的四個狀態(tài),即易感態(tài)(susceptible,記為S)、潛伏態(tài)(exposed,記為E)、感染態(tài)(infected,記為I)以及恢復態(tài)(re-covered,記為R)。生活中多種流行病可以用其組合的傳播模型描述,例如:易感—感染(susceptible-infected,SI)、易感—感染—易感(susceptible-infected-susceptible,SIS)[20]、易感—感染—恢復(susceptible-infected-recovered,SIR)、易感—潛伏—感染—恢復(susceptible-exposed-infected-recovered,SEIR)等模型。
c)已知信息不確定。在流行病傳播過程中,經過時間t,得到了包含每個節(jié)點的感染狀態(tài)的快照。這里的傳播時間t在某些情況下為已知量[21,22],有時候,傳播時間t也是未知量[23~25];同時,感染快照中的狀態(tài)信息有時候包含了所有的節(jié)點[26],有時候只覆蓋了部分節(jié)點[27~29]。
面對如此情況,國內外學者提出了各種方法來克服這些挑戰(zhàn)。2010年,Shah等人[30]針對正則樹上的 SI 傳播模型提出一個謠言中心性(rumor centrality)指標,并據此提出了第一個解決溯源問題的算法。在該算法中,會構建一個正則樹,某個節(jié)點是傳播源的概率,與流行病從該節(jié)點出發(fā)感染其他節(jié)點的所有可能順序的計數(shù)成正比。因此,依據文獻[30]定義的傳播中心性,得分最大的節(jié)點作為傳染源的估計。2016年,Zhu等人[31]提出,依據流行病傳播路徑和傳播源頭的關系,用最可能的傳播路徑對應的源點來推斷傳播源。文獻[31]借鑒文獻[32]中離心率(eccentricity)的思想,定義感染離心率(infection eccentricity)為到所有被感染節(jié)點距離的最大值,并類比于Jordan中心(Jordan center)的概念,定義Jordan感染中心(Jordan infection center)作為感染離心率最小的節(jié)點,并將之稱為反向感染算法(reverse infection algorithm)。Antulov-Fantulin等人[33]考慮在已知傳播過程中T時刻所有節(jié)點的狀態(tài)的情況下,提出了基于蒙特卡羅模擬的定位算法。算法非常容易理解,即通過直接模擬流行病傳播的方法來估計傳播源。對每個可能的傳播源,重復做多次獨立的、傳播時間為T的蒙特卡羅模擬傳播,然后通過計算T時刻所有節(jié)點狀態(tài)來估計傳播源。2017年,Xu等人[34]通過貝葉斯估計法提出了一種算法,在具有稀疏觀察者的復雜網絡中定位流行病源并找到初始時間。2019年,Xu等人[35]基于網絡中有限的觀察者,利用了網絡結構和擴散動力學之間的相關性,提出了一種算法,用以識別推斷傳播源。2021年,Das等人[36]結合多種中心性算法,提出了一種組合網絡中心性方法(CNCA)來推測傳播源。
然而,現(xiàn)有方法可能在某一類溯源檢測上具備一定的準確性,但如果擴展到其他的網絡上時,會出現(xiàn)準確率的下降。比如,謠言中心性算法[30]能夠在多種網絡上取得一定的效果,但其準確率卻沒有特別高;反向感染算法[31]為了最小化源頭到其他節(jié)點的最大距離而忽略了流行病傳播的動態(tài)性;蒙特卡羅模擬的算法[33]需要了解傳播時間?;谏鲜龇治?,本文提出了新算法,該算法首先利用馬爾可夫鏈來初步估計子種群被感染的時間,然后所有感染子種群基于感染時間將自己的身份信息傳播給其所有鄰居,最先獲得所有感染子種群身份信息的子種群推斷為傳播源。通過在多個網絡上的實驗以及多個方法的對比,本文算法的準確率較高,并且在不完全快照的情況下也能夠識別部分傳播源。本文的創(chuàng)新性主要有三點:a)基于反向感染算法,結合復合種群網絡模型以及馬爾可夫鏈理論,提出了一個在復合種群網絡中識別傳播源的新算法;b)利用馬爾可夫鏈說明了子種群感染時間與感染人口狀態(tài)轉移的正相關性,并給出了Pearson相關系數(shù)和Spearman系數(shù)的實驗結果;c)在真實航空網上驗證了本文算法的正確性,證明了本文算法在全部快照和部分快照的情況下,都具有推斷傳播源的優(yōu)良性能。
1 復合種群網絡上的流行病傳播溯源問題
1.1 復合種群網絡
在現(xiàn)實生活中,個體處于諸如社區(qū)和城市這樣的社會單位中。在復合種群網絡模型中,對于社區(qū)和城市這樣的社會單位,定義為子種群,不同的子種群通過交通路線的相互連接,從而形成人口流動。流行病的傳播依賴于人與人之間的接觸,在子種群內部,通過感染態(tài)個體與其他個體的接觸,使得流行病在子種群內部進行傳播;通過子種群之間的人口流動,使得流行病在不同的子種群之間擴散。
對于復合種群網絡,考慮使用一個有向圖G={V,E,W}來表示,其中V是節(jié)點的集合,表示復合種群網絡中的子種群,E是邊的集合,表示不同子種群之間的聯(lián)系,W是節(jié)點與節(jié)點之間遷移率的集合。每一個節(jié)點可以是一個社區(qū)、城市、甚至國家。每一個節(jié)點內部有多個個體,各個個體之間能夠進行隨機接觸,并發(fā)生相互反應(按SIR模型等)。在復合種群網絡中,對于不同的子種群,如果子種群i、j之間存在連邊eij∈E, 那么兩個種群i、j之間的人口就能夠依據遷移率wij∈W進行流動,否則不能進行人口流動。在每一個子種群內部,不同個體之間能夠隨機接觸。圖1給出了復合種群網絡的示意圖,不同的節(jié)點表示不同的子種群,當兩個子種群之間存在連邊時,它們的人口能夠進行流動。在子種群內部,不同的個體之間是能夠自由移動并相互接觸。
1.2 流行病傳播的SIR模型
在流行病模型中,每個個體可以處于最常見的四個狀態(tài):易感(S)、潛伏(暴露)(E)、感染(I)和恢復(R)。在本文中,使用SIR模型。最初,除了傳播源子種群v*中的部分個體處于感染態(tài),其他所有的個體都處于易感態(tài)。本文設置一個時間間隔(時間步),子種群內的每個個體在每個時間步開始時改變它們的狀態(tài)。在每個時間步的開始,每個感染態(tài)個體以感染率β感染它的易感態(tài)鄰居,使得易感態(tài)個體的狀態(tài)能夠轉換為感染態(tài);每個感染態(tài)個體以恢復率γ恢復至恢復態(tài)。由于易感態(tài)個體與恢復態(tài)個體均未表現(xiàn)出癥狀,所以這兩個狀態(tài)的個體在觀測時是不能區(qū)分的。個體a在時間t時刻的狀態(tài)用Pa(t)表示,它的取值如下:
Pa(t)=1? t時刻a個體處于感染態(tài)0? t時刻a個體處于易感態(tài)或恢復態(tài)(1)
在SIR模型中,將s(i,r)表示種群內的易感態(tài)(感染態(tài),恢復態(tài))的個體的比例。一個種群內,SIR模型中的狀態(tài)轉移方程可以由以下微分方程組表示:
dsdt=-βsi? didt=βsi-γi? dsdt=γi(2)
上述方程中忽略了網絡結構,這些微分方程可以用于個體自由移動的網絡的狀態(tài)粗略近似。在復合種群網絡模型中,在一個種群內,例如在子種群i,設置所有個體均可以自由移動并隨機接觸,使用Si(Ii,Ri)表示其易感態(tài)(感染態(tài),恢復態(tài))的人口數(shù),那么對于子種群i,可以使用方程式(3)來表示種群內部人口狀態(tài)的轉移。同時,在SIR模型中,對于不同狀態(tài)的人口,還具備以下性質:易感態(tài)個體可能被他/她的受感染的鄰居感染,感染態(tài)個體可能恢復,并且恢復的個體不能被再次感染。
dSidt=-βSiIiNi? dIidt=βSiIiNi-γIi? dRidt=γIi(3)
對于被感染的子種群c∈V,當其中的傳播源個體1被感染時,發(fā)生的流行病傳播過程可以用圖2表示。
在圖2中,圖(a)是不同個體之間的流行病傳播示意圖,個體1是傳播源。個體1在t=0時刻被感染,并分別在t=1和t=2時刻感染個體2和3。在t=3時刻獲取快照時,獲得的快照如圖(b)所示。對于子種群c,在時間t中的快照信息表示為Xc(t)={Pi(t)}。在t=3時刻,能夠知道個體1~7是否處于感染態(tài),本文將這個信息定義為Xc(3)。由于方程式(1),將無法區(qū)分恢復態(tài)的個體和易感態(tài)的個體,所以此時擁有的信息是
Xc(3)={0,1,0,0,1,1,1}
即某一個子種群c的快照用數(shù)學表示為
Xc(t)={Pi = 1, i∈VIc,Pj=0,j∈Vc\VIc}(4)
其中:VIc表示子種群c內的所有感染個體的集合;Vc表示子種群c內所有個體的集合。上式表示在t時刻獲得的子種群c內部的快照信息,其中,處于感染態(tài)的個體使用1表示,處于易感態(tài)或恢復態(tài)的個體用0表示。
1.3 傳播溯源問題描述
在流行病的傳播中,由于子種群v∈V內部的易感態(tài)個體和恢復態(tài)個體是不可區(qū)分的,所以也不能區(qū)分易感態(tài)的個體(子種群)和恢復態(tài)的個體(子種群)。所以,在時間t獲取的快照為Y(t)={Yv(t),v∈V},其中:
Yv=1? 如果Xv(t)不全為00? 如果Xv(t)全為0(5)
即當某個子種群內所有個體處于易感態(tài)或恢復態(tài)時,對應子種群的信息為0,其余情況取1。復合種群網絡傳播源檢測問題是在給定復合種群網絡圖G和感染快照Y的情況下識別v*,其中t是未知參數(shù)。表1列出了所有符號的描述。
2 基于樣本路徑的傳播溯源
在流行病傳播中,傳播源經常處于感染快照的中心。反向感染算法[31]是尋找感染快照的Jordan感染中心。在該算法中,首先定義網絡中兩個體i、 j的距離為d(i,j),其物理意義為兩個節(jié)點之間最短距離,數(shù)值上的取值為最短距離的值。依據圖論中的偏心率的定義,頂點v的偏心率e(v)是v和圖中任何其他頂點之間的最大距離。圖的Jordan中心是具有最小離心率的節(jié)點。按照類似的術語,定義感染偏心率e(v):給定感染快照Y和網絡圖G,節(jié)點v與所有感染節(jié)點距離的最大值作為感染偏心率。同樣地,定義圖的Jordan感染中心為具有最小感染偏心率的節(jié)點。d(v1,v2)是v1與v2之間的距離,例如,在圖3中,d(v1,v2)=1。感染網絡的感染偏心率e(i)是某個節(jié)點i與其他所有感染態(tài)節(jié)點j∈VI之間的距離的最大值。也就是說,感染偏心率e(v)表示節(jié)點v與其他感染節(jié)點的最大距離,即e(v)=arg maxj∈VId(v,j)。在反向感染算法中,推斷源是圖的Jordan感染中心,即具有最小感染偏心率e(v)的節(jié)點:
v=argminv∈V e(v)=argminv∈V argmaxj∈VI d(v,j)(6)
在圖3中,節(jié)點v3、v6、v8和v9被感染。v1、v2、v3、v4的感染偏心率分別為2、3、4、3,Jordan感染中心為v1。該方法依據感染快照的路徑,Zhu等人[31]命名為最佳樣本路徑推測,其對應的算法稱之為反向感染算法(reverse infection algorithm,RIA),對應的傳播源是具有最小感染偏心率的節(jié)點。
3 基于反向感染的復合種群網絡傳播溯源
反向感染算法在一般網絡特別是樹型網絡上的誤差較小,該算法認為網絡中的每一個節(jié)點代表一個個體,它們的權重相同。然而,隨著交通越來越便利,流行病往往會很快傳播到大尺度的空間范圍,不同個體之間的連接關系(連邊)也常常會發(fā)生變化。因此,現(xiàn)實中的網絡結構需要適應性地調整,流行病傳播溯源常在復合種群網絡(即一個節(jié)點代表一個城市/社區(qū))上進行考慮。在復合種群網絡中,直接使用反向感染算法常常會出現(xiàn)誤差,因為在復合種群網絡中,不同的節(jié)點代表著不同的子種群而不是某一個個體,所以不同的子種群的感染狀態(tài)并不完全相同,即很少出現(xiàn)不同的子種群之間處于感染態(tài)的人數(shù)、處于感染態(tài)人數(shù)的比例等信息完全一致,所以不同的感染子種群之間的權重應當不同。針對上述反向感染算法對復合種群網絡中的流行病傳播溯源的不完全適應問題,結合復合種群網絡,提出了本文算法,可以解決復合種群網絡中的流行病溯源問題。在復合種群網絡中,不同子種群的權重應當與該子種群內部感染狀態(tài)有關。首先,先分析不同子種群的人口狀態(tài)轉移方程。在復合種群網絡的 SIR 模型中,對子種群i而言,其感染態(tài)人口的變化按照下述方程進行變化:
dIdt=βSiIiNi-γIj+∑j∈N(i)wjiIj-∑j∈N(i)wijIi
dIdt=∑Sk=0(βIiNi)k(1-βIiNi)1-k+∑j∈N(i)∑Ijn=1n!∏nm=1wmjiNjm!-∑j∈N(i)∑Iin=1n!∏nm=1wmijNim?。?)
其中:S、I、N 分別代表易感態(tài)人口數(shù)、感染態(tài)人口數(shù)、總人數(shù)。對于子種群i而言,人口狀態(tài)的轉移包括自身人口狀態(tài)轉移以及人口的遷移。其中對于其所有鄰居的遷移對它的影響可以表示為
P=Ni∑j∈N(i)wij-∑j∈N(i)Niwji(8)
而兩鄰居種群i、 j之間的遷移人數(shù)為
Qij=NiwijQji=NjwjiQij=Qji(9)
根據人口流動分析,兩城市之間的人口變化很小,所以式(8)約為0,即P≈0,即忽略人口的遷移現(xiàn)象。此時,不同時刻的人口狀態(tài)與上一個時刻相關,人口狀態(tài)轉移方程用馬爾可夫鏈表示為
Si(t+1)=Si(t)-βSi(t)Ii(t)Ni(t)Ii(t+1)=Ii(t)+βSi(t)Ii(t)Ni(t)-γIi(t)Ri(t+1)=Ri(t)+γIi(t)(10)
通過式(10)可以看到,對于某一個子種群而言,所處的時間不同,其內部不同狀態(tài)的人口規(guī)模不同。忽略遷移影響較小的人口變化,可以依據式(10)來初步估計不同子種群的感染時間tv,來估計復合種群網絡模型下不同子種群的權重。
對于相鄰的子種群,如果在t時刻i種群被感染,t+1時刻其鄰居種群j被感染,顯然,他們的感染時間差與距離差相同,即dt(vi,vj)=d(vi,vj)=tvi-tvj=1。如果在t時刻i種群被感染,t+m(m>1)時刻其鄰居種群j被感染,顯然,他們的感染時間差大于1,即tvi-tvj>1,那么兩者之間的感染時間差表示為dt(vi,vj)=tvi-tvj>d(vi,vj)=1。如果在t時刻i子種群被感染,t+m(m≥1)時刻其鄰居子種群j被感染,t+n(n≥1)時刻其鄰居子種群k被感染,顯然,此時dt(vi,vj)≥1,dt(vi,vk)≥1,但另外兩個不相鄰的子種群j與k之間的距離可能會大于時間之差,即tvj-tvk<2,d(vj,vk)≥2,而在傳播過程中,顯然由子種群j傳播流行病信息至子種群k的時間不會小于距離,即連個不相鄰的子種群j與k的感染時間差表示為dt(vi,vj)=max(d(vi,vj),tvj-tvk)。統(tǒng)一表示,兩個不同的子種群i、 j之間的時間差為
dt(vi,vj)=max(d(vi,vj),tvi-tvj)(11)
所以,在復合種群網絡中,子種群i的感染偏心率表示為dt(vi)=argmin(max(dt(vi,vj)))。那么,推斷的傳播源可以表示為
v=argminvi∈V,vj∈VI(max(dt(vi,vj)))=argminvi∈V dt(vi)(12)
依據上述思想,本文得到了算法1(reverse infection-based-detection source-algorithm,RIDSA)。算法1的時間復雜度為VI+c×dI×VI,空間復雜度為VI+V×VI,其中,VI是被感染子種群的數(shù)量,V是所有子種群的數(shù)量,c是復合種群網絡G的平均度,dI是感染快照Y的直徑。在本算法中,首先需要估算被感染個體的馬爾可夫鏈時間,此時,時間復雜度為VI;其次,在傳播溯源階段,被感染子種群的信息需要傳遞給其鄰居,耗時c×VI,在推斷出傳播源時,信息已傳遍整個網絡,耗時c×VI×dI,總時間復雜度為VI+c×dI×VI。在傳播溯源過程中,本算法需要保存的信息為感染時間矩陣,占空間VI,還需要為每一個子種群儲存接收到的感染信息,占空間V×VI,總空間復雜度為VI+V×VI。
算法1 基于反向感染的復合種群網絡傳播溯源算法(RIDSA)
輸入:感染快照Y;傳播網絡G。
輸出:推斷的傳播源v。
1 for v∈Y do
2? 通過方程式(10)計算馬爾可夫鏈估計的時間MTv
3? 添加估計的時間到矩陣MT(v,MTv)
4 根據MTV對時間矩陣MT進行排序
5 初始化時間t
6 while true do
7? for v∈Y do
8?? if 子種群 v not in G
9??? if 子種群v估計的感染時間MTv≥t do
10???? 添加子種群v到G
11? else
12?? 反向感染子種群v的鄰居
13?? for vi∈N(v) do
14???? 將子種群v添加到子種群vi的信息集合Infvi
15? t=t+1
16? if Infv=Y do
17?? 推斷的傳播源v=v并退出循環(huán)
18? else
19?? 重復 7~19行
20 返回推斷的傳播源v
4 計算機仿真驗證
本章評估了本文算法在不同網絡上的性能,包括人造網絡和真實世界網絡。
4.1 基于復合種群網絡的蒙特卡羅模擬隨機流行病過程方法
在模擬流行病的傳播中,本研究采用蒙特卡羅方法來進行實驗??紤]到實驗的偶然性偏差,本文在不同的傳播網絡中篩選出人口和度數(shù)不同的子種群設計傳播實驗,每一個子種群進行10次的傳播模擬,使用得到的實驗數(shù)據進行傳播溯源。
本仿真的模擬過程主要包括種群內部的傳播以及不同種群之間的擴散兩個方面。
a)傳播過程,即子種群內部的感染過程。在每一個子種群內部,易感態(tài)個體、感染態(tài)個體和恢復態(tài)個體是均勻混合的,每一個個體都有幾率接觸到其他個體。在一個時間步上,根據SIR傳播模型,人數(shù)的變化的期望值是一個確定的值,但考慮到現(xiàn)實生活中的隨機變化,本文采用二項分布來模擬人數(shù)變化的隨機性。需要注意的是,人數(shù)變化的期望值應當是SIR模型下的人數(shù)變化值,即本文采用的模擬過程使得人數(shù)在期望值附近波動。
b)擴散過程,即子種群之間的人口遷移,相當于現(xiàn)實生活中的人員流動。每座城市或者社區(qū)之間,會存在一定的人員流動,且每一個時間步流動的人口數(shù)量可能會有所不同。此處,本文采用多項分布進行模擬,對每一種人口均設置流動性,而子種群的總的流動人數(shù)的期望符合本文給定的遷移率。此處需要注意的是,每個子種群的人數(shù)不同,遷移率也不同。
4.2 馬爾可夫鏈下的傳播時間相關性分析
為了驗證方程式(10)對感染時間先后順序的估計,本文設計了在六個不同的人造網絡(artificial network,AN)上進行模擬實驗,網絡信息如表2所示。
在現(xiàn)實生活中,不同子種群之間的感染率與恢復率會受到氣溫、醫(yī)療等影響而不同,因此本實驗通過設計不同感染率與恢復率來模擬流行病在不同地區(qū)的傳播情況。在不同的人造網絡中,從感染開始時,記錄不同網絡每天的感染狀況與真實時間,根據真實的傳播時間與通過方程式(10)得到的馬爾可夫鏈下估計的時間的數(shù)組,通過方程式(13)來計算Pearson相關系數(shù)和Spearman系數(shù),分析不同感染率與恢復率下的相關性,實驗結果如圖4所示。
r=N∑xi×yi-∑xi∑yiN∑x2i-(∑xi)2×N∑y2i-(∑yi)2r=∑i (xi-)(yi-)∑i (xi-)2∑i (yi-)2(13)
通過圖4可以看到,SIR模型下,當感染率與恢復率發(fā)生變化時,真實傳播時間與估計時間的相關性得分很高,均在0.80以上??梢哉f明,通過馬爾可夫鏈估計時間能夠代表不同子種群感染的先后順序。
4.3 真實航空網上和人造網絡完整快照的傳播溯源結果
本節(jié)使用真實航空數(shù)據構建的美國航空網絡(American airlines network,AAN)和四種常用的人造網絡(規(guī)則網絡RG、ER網絡、WS網絡和BA網絡)進行傳播過程仿真。美國航空網絡中,知道網絡拓撲,包括子人口規(guī)模和旅行人數(shù)和旅客的城市人口。人造網絡中,本文依據美國航空網絡中的遷移率分布等信息,參考文獻[37]構造出了對應的人造網絡。對于傳播源子種群的篩選,本實驗中選擇了低于平均度的子種群占40%,高于平均度的子種群占60%;同時,低于人口平均數(shù)和高于人口平均數(shù)的子種群均約占50%。考慮到不同算法是溯源條件不同,選擇與本文算法與溯源條件相同的算法如謠言中心性算法[30]、RIA算法[31]、STC算法[35]、CNCA算法[36]、有效距離溯源[38]、K-center算法[39]進行比較,通過不同算法的實驗結果來分析本文算法的準確率。
在復合種群網絡中,設置感染概率β=0.3、恢復概率γ=0.2進行流行病的傳播仿真。持續(xù)時間t是依據傳播范圍進行選擇的整數(shù)。對于整個復合種群網絡,通過對感染規(guī)模的估計,當所有子種群均出現(xiàn)感染者時,停止感染仿真,并獲取這個時候的快照。依據此時的快照數(shù)據,對不同算法的溯源準確率和誤差距離進行統(tǒng)計,得到的實驗結果如下。
為了直觀地表示本文算法的準確率,給出了圖5和6的實驗圖。
在圖5(a)、圖6(a)中,本文算法的準確率在80%以上,具備較高的準確率;而其他對比算法的準確率較低。常用的統(tǒng)計方式還有溯源誤差距離,本文統(tǒng)計了真實傳播源頭和推斷源頭之間的距離,實驗結果如圖5(b)、圖6(b)所示。在圖5(b)、圖6(b)中,本文算法的誤差距離在0.2跳以下,即真實的傳播源與推斷的傳播源之間的平均差距為0.2。STC算法利用了感染時間與路徑的相關性,誤差距離在0.7~0.9。而原始的反向感染算法、謠言中心性算法、K-center算法、有效距離的溯源算法的誤差距離在1.0~1.6,誤差大于本文算法。在復合種群網絡上的實驗結果表明,本文算法能夠減小溯源誤差距離。同時,考慮到同一個子種群在網絡的位置相同,那么它的溯源結果很有可能具有收斂的性質,本文統(tǒng)計了相同子種群的實驗結果,發(fā)現(xiàn)90%左右的子種群在多次模擬中的結果收斂,10%左右的子種群實驗結果不收斂。通過分析發(fā)現(xiàn),這些子種群的度值遠大于平均度,從而在本實驗中偶然性的影響下導致了實驗結果不收斂。
在準確率的測試中,本文算法取得了更好的準確率。分析主要原因如下:本文算法考慮了子種群內部的感染信息與未感染子種群的信息;其次,STC算法的表現(xiàn)優(yōu)于其他的算法,是因為其考慮了節(jié)點感染先后順序和有效距離路徑的相關性;再次,其他算法幾乎只從傳播路徑的角度分析了傳播源,而沒有考慮復合種群網絡上流行病傳播的更多信息,比如節(jié)點感染時間與節(jié)點感染先后順序。RC、RIA會考慮了感染節(jié)點與傳播路徑的關系;CNCA雖然考慮了多種中心性算法,但這些中心性算法之間本身就存在一定的耦合性;ED、K-center提出了設置有效距離的方法,但是對子種群本身的感染狀態(tài)等信息未加以利用。從結果上來看,本文算法可以利用復合種群網絡上的更多信息,從而能更準確地推斷真正的傳播源。當復合種群網絡的信息減少,比如,當所有的子種群的人數(shù)進行坍縮至1人時,復合種群網絡退化為普通的傳播網絡,此時,本文算法也會坍縮至RIA算法,從而會使實驗準確率下降。值得一提的是,當將復合種群網絡看做普通網絡時,其余算法的傳播溯源的結果仍能夠保持一定的準確率,而當復合種群網絡快照進行減小時,本文算法的準確率也會下降。復合種群網絡上的實驗結果表明,本文算法能夠提高溯源準確率。接著,本文統(tǒng)計了不同算法在不同誤差距離的分布情況,得到的結果如圖7和8所示。
可以看到,本文算法90%以上的誤差為0,即找到了傳播源頭,而未找到傳播源頭時,誤差也在1跳內。STC算法的誤差距離在0~2,也具備較高的準確率。而其他算法的誤差距離多在1~2,誤差大于本文算法。RC、RIA、CNCA算法利用感染節(jié)點在感染快照中的位置來進行傳播溯源,卻沒有考慮到復合種群網絡中,不同節(jié)點之間的距離與節(jié)點之間的遷移率存在關系;ED、K-center算法通過設置有效距離來進行傳播溯源,考慮了遷移率的因素;STC算法使用節(jié)點感染順序與有效距離之間的關系來進行傳播溯源,準確率高于其他的算法;本文算法考慮了復合種群網絡中的子種群內部的信息,結合子種群在感染快照中的位置來進行傳播溯源,相較于其他算法,結果有明顯的提升。在復合種群網絡中,不同的子種群內部感染人數(shù)和感染比例不同,那么處于子種群被感染的時間越久,那么該子種群是感染源的概率越大;該子種群所占的權重也應當越大。這是復合種群網絡的特性,也是現(xiàn)實生活中個體的現(xiàn)實屬性,是如今的傳播溯源算法需要考慮的問題。
4.4 真實航空網上部分快照的傳播溯源結果
對于信息獲取不完全的情況,通過隨機減少快照的比例來進行實驗。假設獲取快照的比例為20%、30%、50%、70%、100%時,本文算法的準確率和誤差距離的實驗結果如圖9所示。
在圖9中,可以看到,當能夠獲取20%的感染快照時,本文算法已經能夠達到20%左右的溯源準確率,且誤差距離在1.0附近。在不同的復合種群網絡中,當快照規(guī)模達到50%時,本文算法溯源準確率已經達到50%;當獲取完整快照時,幾乎已經能夠準確找到所有感染源頭。圖9(a)(b)顯示了本文算法在獲取不同快照信息時的準確率和誤差距離。隨著獲取信息的增加,準確率在逐漸增加,誤差距離在逐漸減小。在不要求更高的準確率時,部分快照的信息已經能夠用來溯源;如果要求更高的準確率,本文算法需要更準確的快照。在不同的流行病的感染率、恢復率可能會不同。為了驗證本文算法的魯棒性,本節(jié)在美國航空網絡(AAN)設置不同的感染率、恢復率進行流行病傳播仿真,使用仿真結果進行溯源實驗,得到的結果如圖10所示。
圖10(a)(b)分別表示在AAN中的不同感染率、恢復的溯源結果。可以發(fā)現(xiàn),不論是在固定感染率并設置不同的恢復率,還是固定恢復率并設置不同的感染率的實驗中,本文算法的準確率均在80%以上,準確率的分布相近。本文算法在不同的流行病參數(shù)下均表現(xiàn)出較高的準確率,表明其具有良好的擴展性。
4.5 2022年上海市新冠疫情的傳播溯源結果
2020年3月份左右,上海市爆發(fā)了新冠疫情。本節(jié)收集了上海市16個區(qū)的感染數(shù)據。依據2022年4月5日至8日的數(shù)據,使用重力模型[40]進行來模擬人口流動,當人口變化保持穩(wěn)定時,構造出的上海市的網絡如圖11所示。使用本文算法進行傳播溯源。不同區(qū)的得分如圖12所示。
在圖12中可以看到,從2022年4月5日至8日,浦東區(qū)的得分均為最高,從而推斷浦東區(qū)為傳播源頭。
5 結束語
本文研究了流行病在復合種群網絡中的傳播,通過分析不同子種群的感染規(guī)模及狀態(tài)來進行傳播源檢測的問題。結合RIA算法以及流行病在復合種群網絡上的傳播特性,本文考慮了復合種群網絡與普通網絡的區(qū)別、流行病傳播的馬爾可夫鏈、時間變化等因素,提出了一種適合復合種群網絡的傳播溯源算法。在幾個不同的網絡上評估了本文算法的性能。仿真結果表明,該算法在傳播源檢測問題上表現(xiàn)良好。同時,當實驗只能獲取部分快照時,本文算法也表現(xiàn)出了良好的準確率。本文算法考慮了子種群內部的感染信息與未感染子種群的信息,以及考慮了感染快照的路徑信息。其次,STC算法的表現(xiàn)優(yōu)于其他的算法,是因為其考慮了節(jié)點感染先后順序和有效距離路徑的相關性。再次,其他算法幾乎只從傳播路徑的角度分析了傳播源,而沒有考慮復合種群網絡上流行病傳播的更多信息,比如節(jié)點感染時間與節(jié)點感染先后順序。值得注意的是,其他算法幾乎沒有依賴于網絡的特殊性,所以這些方法在許多網絡上均具備一定的魯棒性,而本文算法則更加依賴復合種群網絡上的快照的實驗數(shù)據,這將直接與本文算法的準確率相關聯(lián)。作為進一步的研究,有三個方向值得探索。a)通過真實傳播時間與馬爾可夫鏈下估計的時間的相關性分析,本文發(fā)現(xiàn)在復合種群網絡模型中,給定感染率與恢復率,流行病的傳播先后順序能夠進行估計。如果能夠準確地估計不同子種群的傳播時間,或許能夠更加準確地識別傳播源。b)到目前為止,本文已經得出了SIR模型下的單一信息源的方法,將它擴展到其他模型下的多信息源檢測,對它們進行擴展是很有意義的。c)當快照足夠完整時,如果能夠準確估計傳播時間,那么新的方法或許能夠不依賴于傳播路徑,即能夠在時變網絡等模型下進行傳播溯源。
參考文獻:
[1]Bell D M. Public health interventions and SARS spread, 2003[J]. Emerging Infectious Diseases, 2004,10(11): 1900-1906.
[2]Neumann G, Noda T, Kawaoka Y. Emergence and pandemic potential of swine-origin H1N1 influenza virus[J]. Nature, 2009,459(7249): 931-939.
[3]Wu Fan, Zhao Su, Yu Bin, et al. A new coronavirus associated with human respiratory disease in China[J]. Nature, 2020, 579(7798): 265-269.
[4]李翔, 李聰, 王建波, 復雜網絡傳播理論——流行的隱秩序[M]. 北京: 高等教育出版社, 2020:10. (Li Xiang, Li Cong, Wang Jianbo. Theory of spreading on complex networks: hidden rules of epidemics [M]. Beijing: Higher Education Press, 2020: 10.)
[5]Shah D, Zaman T. Detecting sources of computer viruses in networks: theory and experiment[C]//Proc of ACM SIGMETRICS International Conference on Measurement and Modeling of Computer Systems. New York: ACM Press, 2010: 203-214.
[6]Wang Pu, González M C, Hidalgo C A, et al. Understanding the spreading patterns of mobile phone viruses[J]. Science, 2009,324(5930): 1071-1076.
[7]Bellet A, Guerraoui R, Hendrikx H. Who started this rumor? Quantifying the natural differential privacy of gossip protocols[C]//Proc of the 34th International Symposium on Distributed Computing. [S.l.]: Schloss Dagstuhl-Leibniz-Zentrum für Informatik, 2020: 8:1-8:18.
[8]Grinberg N, Joseph K, Friedland L, et al. Fake news on Twitter du-ring the 2016 U.S. presidential election[J]. Science, 2019,363(6425): 374-378.
[9]Cator E, Van M P. Second-order mean-field susceptible-infected-susceptible epidemic threshold[J]. Physical Review E, 2012,85(5): 056111.
[10]Jiang Jiaojiao,Wen Shui,Yu Shui,et al. Identifying propagation sources in networks:state-of-the-art and comparative studies[J]. IEEE Communications Surveys & Tutorials, 2017,19(1):465-481.
[11]Waniek M, Holme P, Farrahi K, et al. Trading contact tracing efficiency for finding patient zero[J]. Scientific Reports, 2022,12(1): 22582.
[12]陳鈺書, 劉影, 唐明. 具有非馬爾可夫旅途感染的流行病傳播模型研究[J]. 計算機應用研究, 2023,40(6): 1739-1749. (Chen Yushu, Liu Ying, Tang Ming. Study on epidemic transmission model with non-Markovian journey infection[J]. Application Research of Computers, 2023,40(6): 1739-1749.)
[13]Barabási A L, Albert R. Emergence of scaling in random networks[J]. Science, 1999,286(5439): 509-512.
[14]Pastor-Satorras R, Vespignani A. Epidemic spreading in scale-free networks[J]. Physical Review Letters, 2001,86(14): 3200-3203.
[15]Newman M E J. The structure and function of complex networks[J]. SIAM Review, 2003,45(2): 167-256.
[16]Belyi A, Bojic I, Sobolevsky S, et al. Global multi-layer network of human mobility[J]. International Journal of Geographical Information Science, 2017,31(7): 1381-1402.
[17]沈力峰, 王建波, 杜占瑋, 等. 基于社團結構和活躍性驅動的雙層網絡傳播動力學研究[J]. 物理學報, 2023,72(6): 356-364. (Shen Lifeng, Wang Jianbo, Du Zhanwei, et al. A study on propagation dynamics of two-layer networks based on association structure and activeness drive[J]. Acta Physica Sinica, 2023,72(6): 356-364.)
[18]Sartori F, Turchetto M, Bellingeri M, et al. A comparison of node vaccination strategies to halt SIR epidemic spreading in real-world complex networks[J]. Scientific Reports, 2022,12(1): 21355.
[19]Armbruster A, Holzer M, Roselli N, et al. Epidemic spreading on complex networks as front propagation into an unstable state[J]. Bulletin of Mathematical Biology, 2022,85(1): 4.
[20]周海平, 賴兵兵, 劉妮. 一類考慮病毒發(fā)生變異的SIS疾病傳播模型[J]. 計算機應用研究, 2014,31(9): 2773-2775. (Zhou Haiping, Lai Bingbing, Liu Ni. A class of SIS disease transmission models considering virus occurrence mutation[J].Application Research of Computers, 2014,31(9): 2773-2775.)
[21]Pinto P C, Thiran P, Vetterli M. Locating the source of diffusion in large-scale networks[J]. Physical Review Letters, 2012,109(6): 068702.
[22]Lokhov A Y, Mézard M, Ohta H, et al. Inferring the origin of an epidemic with a dynamic message-passing algorithm[J]. Physical Review E, 2014,90(1): 012801.
[23]Luo Wuqiong, Tay W P, Leng Mei. Identifying infection sources and regions in large networks[J]. IEEE Trans on Signal Processing, 2013,61(11): 2850-2865.
[24]Luo Wuqiong, Tay W P. Identifying multiple infection sources in a network[C]//Proc of Conference Record of the 46th Asilomar Confe-rence on Signals, Systems and Computers. Piscataway, NJ: IEEE Press, 2012: 1483-1489.
[25]Luo Wuqiong, Tay W P, Leng Mei. Identifying infection sources and regions in large networks[J]. IEEE Trans on Signal Processing, 2013, 61(11): 2850-2865.
[26]Dong Wenxiang, Zhang Wenyi, Tan C W. Rooting out the rumor culprit from suspects[C]//Proc of IEEE International Symposium on Information Theory. Piscataway, NJ: IEEE Press, 2013: 2671-2675.
[27]Karamchandani N, Franceschetti M. Rumor source detection under pro-babilistic sampling[C]//Proc of IEEE International Symposium on Information Theory. Piscataway, NJ: IEEE Press, 2013: 2184-2188.
[28]Louni A, Subbalakshmi K P. A two-stage algorithm to estimate the source of information diffusion in social media networks[C]//Proc of IEEE Conference on Computer Communications Workshops. Pisca-taway, NJ: IEEE Press, 2014: 329-333.
[29]Luo Wuqiong, Tay W P, Leng Mei. How to identify an infection source with limited observations[J]. IEEE Journal of Selected Topics in Signal Processing, 2014,8(4): 586-597.
[30]Shah D, Zaman T. Rumors in a network: whos the culprit?[J]. IEEE Trans on Information Theory, 2011,57(8): 5163-5181.
[31]Zhu Kai, Ying Lei. Information source detection in the SIR model: a sample-path-based approach[J]. IEEE/ACM Trans on Networking, 2016,24(1): 408-421.(下轉第2693頁)
收稿日期:2023-02-07;修回日期:2023-04-06? 基金項目:國家自然科學基金面上項目(62173065)
作者簡介:陽成(1996-),男,四川綿陽人,碩士,主要研究方向為復雜網絡、網絡傳播動力學;王建波(1980-),男(通信作者),四川成都人,講師,碩導,博士,主要研究方向為復雜網絡、網絡傳播動力學(phyjbw@gmail.com);許小可(1979-),男,遼寧人,教授,博導,博士,主要研究方向為社交網絡大數(shù)據、數(shù)據可視化、計算傳播學;杜占瑋(1988-),男,河北人,副研究員/研究助理教授,博導,博士,主要研究方向為計算流行病學、機器學習、數(shù)據科學.