999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

基于辟謠機制的時滯謠言傳播模型的動力學分析*

2020-02-18 03:17:36朱霖河李玲
物理學報 2020年2期
關鍵詞:模型系統

朱霖河 李玲

(江蘇大學理學院,鎮江 212013)

在社交網絡謠言傳播模型中,考慮到辟謠機制和時滯效應對網絡謠言傳播的影響,建立基于辟謠機制和時滯效應的SIR謠言傳播模型.利用再生矩陣譜半徑方法得到R0;根據二次函數圖像特征給出謠言盛行平衡點存在的條件;通過特征值理論和Routh-Hurwitz判據確定無謠言平衡點和謠言盛行平衡點的局部穩定性以及發生Hopf分支的條件;數值仿真結果表明政府和媒體發布的辟謠信息會加快謠言收斂的速度和降低謠言傳播者的最大密度.

1 引 言

隨著移動通信網絡技術的快速發展,數據傳輸和消息傳送的速率飛速上升,使得社交網絡在我們的日常生活中扮演著愈加重要的角色,社交網絡的多元化、靈活性、快捷性,使其成為最流行的信息發布和傳播方式[1].然而由于社交網絡平臺的低門檻的信息傳播特征,使人為發布謠言與傳播謠言的可能性大大增加.控制網絡謠言的迅猛傳播需要更好地揭示網絡謠言傳播的內在規律,因此建立數學模型來分析網絡謠言傳播的動力學特征迫在眉睫.

近年來,有很多基于數學模型研究網絡謠言傳播的控制方法[2?11],如 2011年Zhao等[2]結合遺忘因素對社交網絡謠言傳播的影響,提出了具有遺忘機制的S IR 模型.2012 年,顧亦然和夏玲玲[3]考慮到人群中還具有潛伏期的謠言傳播者(接受謠言但并未立即傳播的人群),建立了在線社交網絡的謠言傳播 S EIR 模型,并且給出了基于熟人免疫控制策略的謠言傳播控制方法.2015年,趙洪涌等[5]綜合時空延遲效應和媒體覆蓋對社交網絡謠言傳播的影響,構建了具有媒體影響的 I SM 模型,仿真結果表明:時間滯后會增加模型的收斂時長,媒體報道可以抑制謠言的快速傳播.2016年Liu等[6]論證了免疫機制具有暫時性,建立了具有暫時免疫機制的 S IR 模型.2018 年,Ma 等[7]考慮到網絡節點間具有相互作用,將獨立的謠言傳播者引入SIR社交網絡謠言傳播模型中,仿真結果表明獨立傳播者出現的可能性與謠言傳播范圍呈正相關;吳曉等[8]認為用戶遺忘率和文化程度以及消息緩存可以影響網絡謠言的動態傳播過程,提出了全新的URBD網絡謠言傳播動力學模型,通過平均場理論研究模型的傳播規律;王飛雪和李芳[9]注意到由于個體行為的異質性導致謠言在不同節點之間傳播概率不同,將網絡節點自身的影響引入傳統S IR模型中建立新的謠言傳播模型,通過仿真實驗發現與原有SIR模型相比,改進后的模型能夠更有效地控制網絡謠言的傳播.2019 年,張菊平等[10]考慮遺忘因素對謠言傳播的影響,構建了具有真實信息傳播者的SITR社交網絡謠言傳播模型,理論方面給出了系統穩定性的基本判別條件,數值模擬結果表明真實信息傳播者的初值越高,系統收斂時間越短,謠言傳播者所能達到峰值越低.

在社交網絡中,用戶可以和朋友進行一對一的交流,也可以完成一對多的交流:比如在群聊中傳播信息,或通過更新其動態傳播信息[12?15].因此可以將社交網絡中信息傳播分為兩類:群體傳播和一對一傳播.事實上,謠言也是通過上述兩個渠道進行擴散的.Jia等[16]考慮傳播謠言的兩種方式:群體傳播和一對一傳播,創新性地建立了如下的SIR模型:

在此模型中,用戶節點分為三種狀態:易感狀態S、感染狀態I、恢復狀態R.易感節點未知謠言、感染節點傳播謠言、恢復節點接受謠言但不傳播謠言;感染節點通過群體傳播向易感節點傳播謠言使其轉變為一個感染節點,概率為 λ1,或使其轉變為恢復狀態,概率為 θ.感染節點通過一對一傳播向易感節點傳播謠言使其轉變為一個感染節點,概率為 λ2,或使其轉變為恢復狀態,概率為 θ.γ1表示感染節點選擇群體傳播的概率,γ2表示感染節點選擇一對一傳播的概率,α 表示恢復率,β 表示暫時免疫(恢復的節點轉變為易感節點)的概率,〈k〉表示網絡的平均度,Λ 表示出生率,u 表示死亡率.作者詳細研究了系統(1)式的動力學特征,分析了無謠言平衡點與謠言盛行平衡點的局部漸近穩定性和全局穩定性,給出了簡單易行的判別準則.然而,事實上,在謠言持續傳播的過程中,為了降低謠言對社會穩定造成的影響,政府或媒體部門將采取合理的辟謠措施干預謠言傳播.此外,由于環境影響辟謠過程可能存在時間滯后現象,進而影響謠言干預效果.因此,本文基于上述事實將進一步改進模型(1)式,提出新的具有辟謠機制的謠言傳播模型.

本文首先引入時延辟謠機制建立新的謠言傳播模型;其次研究模型平衡點的存在性、局部漸近穩定性及發生Hopf分支[17?20]的條件;接下來通過數值模擬驗證前面的理論結果的可靠性,最后總結上述研究內容.

2 模型建立

眾所周知,政府和媒體發布的辟謠信息可以抑制謠言在社交網絡中的快速傳播[21?23],如冉茂潔等[22]綜合個體興趣度差異和辟謠機制在網絡謠言傳播過程中的影響,提出了新的 I WSR 謠言傳播模型,數值仿真結果表明提升個體辨識能力和加強辟謠力度可以降低網絡謠言傳播速度;王筱莉等[23]考慮到政府辟謠策略對謠言傳播的約束,將辟謠機制引入非均勻網絡謠言傳播模型中得到新的謠言傳播模型,論證了政府辟謠后謠言傳播者的峰值將會下降的重要事實.可以看到,模型(1)式中考慮線性遏制機制 αI 來刻畫政府和媒體發布的辟謠信息是為了簡化問題,但是在現實生活中,政府和媒體對謠言傳播的約束能力是有限的,所以應進一步考慮選取合適的辟謠函數,以刻畫有限的約束能力,如圖1 中所示.飽和辟謠函數[24?26]h(I)=rI/(1+αI)由Capasso和Serio引入流行病模型,其中 α 表示感染者擁擠或變化對易感者個體的影響,rI用于衡量疾病的感染能力,1 /(1+αI) 代表治療延時的感染者的反作用,r I/(1+αI) 為當易感人群增加時,衡量行為變化引起的抑制效應.由于網絡謠言傳播和傳染病傳播在傳播機理上十分相似,所以在網絡謠言傳播中,引入辟謠函數h(I)=rI/(1+αI)來刻畫政府和媒體對謠言傳播的約束能力.

圖1 α 對辟謠函數 h (I)=rI/(1+αI) 的影響Fig.1.The influence of α on the rejection mechanism h(I)=rI/(1+αI).

此外,政府和媒體發布辟謠信息到謠言傳播者接受信息并恢復為免疫者是有時間滯后的[27?29],比如謠言傳播者未能及時收到辟謠信息或接受信息后未能立刻做出反應,因此還需在原有辟謠函數中加入時滯,得到新的辟謠函數用以刻畫這種滯后作用.此外,在網絡謠言傳播模型中,對于同一個謠言而言,接受謠言并且不傳播謠言的用戶已經對該謠言進行分析并且做出了自己的主觀判斷;從記憶的角度來說,該類用戶對此謠言已經進行了信息加工并存儲,所以接受謠言且不傳播謠言的用戶通過“遺忘”這種“暫時免疫機制”轉化為未知謠言的用戶的可能性相對較小,在理論上是可以將其舍去的.基于上述分析并結合文獻[16]的模型,得到以下改進的 S IR 社交網絡謠言傳播模型:

3 平衡點的存在性與后向分支分析

眾所周知,基本再生數在網絡謠言傳播的動力學分析中起著重要的作用.一般來說,基本再生數是用來衡量謠言傳播的感染能力的,為判斷謠言在社交網絡中是流行還是消亡提供了重要依據.顯然,系統(2)式總是存在一個無謠言平衡點E0=(Λ/u,0,0),根據Driessche[30]和 Al-Darabsah[31]計算基本再生數的方法,可以得到基本再生數為接下來將討論在何種條件下,系統(2)式將存在謠言盛行平衡點.

令系統(2)右端為零,可以得到

因此,可以得到

由(4)式有

當 I=0 時,有 S=Λ/u,R=0,即為無謠言平衡點E0=(Λ/u,0,0).將(5)式中S的表達式代入(8)式可得

等式兩邊同時乘以[1+αI][(m1+m2)I+u],得到

其中

因為α>0,u>0,m1>0,m2>0,可得 a>0.下面將對剩下的系數進行分類討論進而得到方程(9)的解的情況.

2) 當c=0 時,方程 (9)等價于

方程 (10)的解為I1=0或I2=-b/a.若 b<0,則有I2>0,也就是說方程(10)存在一個正解,這意味著系統(2)式存在一個謠言盛行平衡點記為;相反的若b≥0,則有I2<0,那么方程 (10)沒有正解,表明系統(2)式不存在謠言盛行平衡點.

3)當c>0 時,令

若b≥0,方程 (9)無正解,表明系統 (2)沒有謠言盛行平衡點.

若b<0,則有以下三種情況:

i)Δ>0,方程 (9)有兩個正解,即系統 (2)式存在兩個謠言盛行平衡點,記為

ii) Δ=0,方程 (9) 有一個正解,即系統 (2) 式存在一個謠言盛行平衡點;

iii) Δ<0,方程 (9)無正解,這表明系統 (2)式沒有謠言盛行平衡點.

由 c=m1Λ(1/R0-1),可得如下等價關系:

因此,有以下結論.

定理11)若R0>1,系統(2)式存在一個謠言盛行平衡點并且

2)若R0=1且b ≥0,系統 (2)式不存在謠言盛行平衡點;若R0=1且b<0,系統 (2) 式存在唯一的謠言盛行平衡點,并且

3)若R0<1且b ≥0,系統 (2)不存在謠言盛行平衡點;若R0<1且b<0,則有:

iii) Δ<0,系統(2)式不存在謠言盛行平衡點.

由以上分析可知,當R0<1,b<0,Δ>0 時,系統(2)式存在兩個謠言盛行的平衡點,這意味著謠言不會被消除,仍在社交網絡上傳播.因此,為了控制謠言的傳播,當后向分支發生時,將計算一個更小的閾值.

定理 2若R0<1,b<0,Δ>0 成立,則系統(2)式將會在R0=1 時發生后向分支.

證明:從定理1可以看出,當參數滿足R0<1,b<0且 Δ>0 時,系統(2)式存在兩個謠言盛行平衡點,換言之,當 R0落在區間 (R*,1) 時,一個后向分支將會在R0=1處出 現.同時,R*的值由Δ=0確定,因此為了得到一個較小的閾值 R*,令Δ=0,進而可以得到臨界的基本再生數為

同時有

關于上面定義的臨界基本再生數,可以總結得到如下定理.

定理31)若 R*<R0<1 系統(2)式存在兩個謠言盛行平衡點,換言之系統(2)式將會在R0=1處發生后向分支;

2)若R0>1,系統(2)式存在唯一的謠言盛行平衡點,換言之,謠言將會繼續在社交網絡中傳播;

3)若R0<R*<1,系統(2)式不存在謠言盛行平衡點,意味著謠言最終會消亡.

4 平衡點的穩定性分析

4.1 無病平衡點E0的穩定性分析

系統(2)式在無謠言平衡點E0的雅可比矩陣為

因此,系統在平衡點E0處的特征方程為

顯然,λ=-u 是方程 (11)的二重根,接下來只需要考慮的特征根的正負情況.

當 τ>0 時,令

若 F (0)<0,即R0>1 時,方程 (12)至少有一個正根,因此,E0是不穩定的.

證 :當時,E0是局部漸近穩定的;

假設 λ=iω(ω>0) 是方程 (12)的解,將其代入方程(12)可得

分離實部虛部并使等式兩邊平方相加我們可以得到

則有以下兩種情況:

從上面的分析中可知:當R0>1 時,E0是不穩定的,因此舍棄情況 2);又因為 m1,Λ,u≥ 0,所以可以得知 u>r,因此當是局部漸近穩定的.

若 (u-m1Λ/u)2<r2,方程 (13)存在正解ω0>0,使得 λ=iω0為方程 (12) 的解,又

從而,系統在 τ=τ0附近發生Hopf分支現象.

定理 41) 當 0<Rc<1 且R0<1 時,無謠言平衡點E0對所有的 τ ≥0 都是局部漸近穩定的;

2) 當 Rc>1 且R0<1 時,無謠言平衡點E0對所有的 0 ≤ τ<τ0都是局部漸近穩定的;當τ>τ0無謠言平衡點E0不穩定,換言之,系統 (2) 在τ=τ0附近發生Hopf分支;

3) 當R0>1 時,無謠言平衡點E0是不穩定的.

4.2 謠言盛行平衡點 的穩定性分析

直接計算可得其特征方程為

化解得

其中

其中

當 τ=0 時,方程 (19)等價于

其中

分離實部虛部得

等式兩邊平方再相加得

引理 11)若f2≥0,方程 (23)無正實根.

2)若f2<0,方程 (23)總會有一個正實根.

證明根據可得

由此可得m42≥m32,那 么f1=u2+m42-m32≥0.

從上面的分析中可知,當f2≥0 時,有f1≥0,因此,根據二次方程圖像特征可知方程(23)沒有正根,這意味著不存在 ω>0滿足 λ=iω是方程(19)的一個解,也就是說方程(19)不存在正根,那么謠言盛行平衡點E1*是局部漸近穩定的.相反地,若f2<0,易知無論 f1取何值時,方程 (23) 至少存在一個正根,也就是說方程(19)有一對純虛根λ=±iω0*.

直接計算可得

從而有

由 (21)式,有

因此,

綜上所述,可以得到如下結論.

定理51)若R0>1,m4>0成立,則當τ=0時謠言盛行平衡點是局部漸近穩定的.

2) 若R0>1,m4>0,f2≥0 成立,則謠言盛行平衡點對所有的都是局部漸近穩定的.

3) 若R0>1,m4>0,f2<0成立,則當時,謠言盛行平衡點是局部漸近穩定的;當時,謠言盛行平衡點不穩定,即謠言盛行平衡點在附近發生Hopf分支.

5 數值模擬

例 1在系統 (2)式中,選取參數 λ1=0.5,λ2=0.8,γ1=0.4,γ2=0.5,θ=0.1,u=0.25,λ1=0.5,λ2=0.8,γ1=0.4,γ2=0.5,θ=0.1,u=0.25,Λ=0.25,α=5,r=0.15,〈k〉=8,通過簡單計算可以得到基本再生數R0=0.625<1,Rc=2.5>1.0,τ0=10.4723,因此,系統 (2) 式的無謠言平衡點E0=(1,0,0) 在τ=τ0處發生 Hopf分支,數值模

擬結果如圖2所示.

從圖2(a) 可以看出:隨著時間的推移,S (t) 趨向于 1,I (t)和R (t) 趨向于零,這意味著散播謠言的用戶的密度最終將轉變為零且不再發生波動,即謠言最終消除.同時,從圖2(b)可以看出:隨著時間推移,傳播謠言個體的密度處于持續波動狀態,這意味著謠言在網絡中以周期形式持續傳播.該現象表明時滯的擴大會導致平衡點穩定性發生變化,換言之,如果政府延遲執行有效的辟謠機制,將會導致謠言在網絡中呈現周期性爆發現象,以致于加大了網絡謠言傳播的控制難度,極大的危害了網絡信息安全.因此,通過分析時滯導致的Hopf分支現象,可以有效地反映政府辟謠機制的執行功效,從而有利地指導政府開展及時合理的辟謠政策.

例2本例主要研究恢復率r對傳播謠言用戶密度的影響.首先,考慮R0<1 的情況,在系統 (2)式中,不妨取參數 λ1=0.5,λ2=0.8,γ1=0.6,γ2=0.5,θ=0.1,u=0.25,Λ=0.25,α=5,〈k〉=8,r 的值為 0.05,0.15,0.25,0.35,得出I(t) 的變化趨勢如圖3(a),其次考慮R0>1 的情況,令 λ1=0.7,λ2=0.8,γ1=0.8,γ2=0.9,θ=0.1,u=0.25,Λ=0.25,α=5,〈k〉=8,r 的值為 0.05,0.15,0.25,0.35,得出 I (t) 的變化趨勢如圖3(b).

從圖3(a) 可以看出,隨著時間的推移,I (t) 逐漸接近零,并且隨著恢復率r的增大,I (t) 趨于穩定的時間越來越短,即恢復率r越大,謠言消亡的速度越快.換言之,恢復率越高,傳播謠言用戶的密度下降得越快,趨于穩定的時間越短.同時,從圖3(b)可以看出,隨著恢復率r的值增大,I (t) 的峰值變小,也即恢復率越高,那么傳播謠言用戶的密度的峰值越低.

例3本例主要考察抑制效應α對傳播謠言用戶密度的影響.首先,考慮R0<1 的情況:在系統 (2) 中,選取參數 λ1=0.5,λ2=0.8,γ1=0.4,γ2=0.5,θ=0.1,u=0.2,Λ=0.2,〈k〉=8,分別取 α為 0,3,5,8,得出I(t) 的變化趨勢如圖4(a),其次考慮R0>1 的情況,令 λ1=0.7,λ2=0.8,γ1=0.8,γ2=0.9,θ=0.1,u=0.4,Λ=0.4,〈k〉=8,同樣地,取α為 0,3,5,8,得出I(t) 的變化趨勢如圖4(b).

圖2 (a) τ=10.3 ;(b) τ=10.5Fig.2.(a) τ=10.3 ;(b) τ=10.5.

圖3 (a)R0<1 ;(b) R0>1Fig.3.(a)R0<1 ;(b)R0>1.

圖4 (a)R0<1 ;(b) R0>1Fig.4.(a)R0<1 ;(b)R0>1.

圖5 出生率 Λ 與恢復率 r 對基本再生數 R 0 的影響Fig.5.Influence of birth rate Λ and recover rate r on the basic productive number R 0.

從圖4(a) 可以看出:隨著時間的推移,I (t) 趨向于零,并且隨著抑制效應α的增大,I (t) 趨于穩定的時間越來越長,也就是抑制效應α的值越高,謠言消亡的速度越慢.換言之,抑制效應影響程度越高,那么傳播謠言用戶的密度下降得越慢,趨于穩定的時間越長.同時,從圖4(b)可以看出,隨著抑制效應α的值增大,I (t) 的峰值越高,這表明當R0>1時,如果抑制效應α越高,那么傳播謠言用戶的密度的峰值越高.

例4本例主要研究出生率Λ和恢復率r對基本再生數 R0的影響,選取系統 (2)式中參數 λ1=0.5,λ2=0.8,γ1=0.7,γ2=0.6,θ=0.1,u=0.35,α=2,〈k〉=8.觀察 R0的變化趨勢,如圖5 所示.圖5表示在不同的出生率 Λ 、恢復率r下R0的取值情況.從圖5(a)中可以發現:隨著恢復率r不斷增大,基本再生數 R0的值不斷減小;當恢復率r逐漸增大并增大到一定程度時,R0小于1,意味著謠言最終消亡.相反地,隨著出生率 Λ 的增大,基本再生數R0的值不斷增大,也就是說出生率 Λ 與基本再生數 R0呈正相關關系.圖5(b)給出了R0<1與R0>1 時對應的出生率 Λ 與恢復率r的取值范圍,明確了謠言傳播的可知措施,如增加恢復率r等.

例 5在本例中,選取參數 λ1=0.5,λ2=0.5,γ1=0.4,γ2=0.2,θ=0.1,u=0.22,Λ=0.8,α=2,r=0.55,〈k〉=8,通過簡單計算可以得到基本再生數R0=0.9917,R*=0.9734,即 R*<R0<1,滿足定理3條件,那么系統(2)式存在兩個謠言盛行平衡點,即系統(2)式將會在R0=1 時發生后向分支,如圖6所示.

圖6 后向分支Fig.6.Backward branch.

6 結 論

綜合考慮辟謠機制和時滯效應對社交網絡謠言傳播機制的影響,建立了基于辟謠機制的時滯謠言傳播模型.利用再生矩陣譜半徑方法求出模型的基本再生數 R0;根據二次函數圖像特征給出謠言盛行平衡點存在的條件:若R0>1,系統(2)式總會存在謠言盛行平衡點;若R0=1 且 b ≥0,系統(2)式沒有謠言盛行平衡點;若R0=1 且 b<0,系統(2)式存在唯一的謠言盛行平衡點;若R0<1且b≥0,系統(2)式沒有謠言盛行平衡點;若R0<1 且 b<0,Δ=0,系統 (2)式存在兩個謠言盛行平衡點與;若R0<1 且 b<0,Δ<0,系統(2)式存在謠言盛行平衡點;若R0<1 且b<0,Δ<0,系統 (2)式沒有謠言盛行平衡點.通過特征值理論和Routh-Hurwitz判據確立了無謠言平衡點和謠言盛行平衡點的局部漸近穩定性:對無謠言平衡點 E 0,當 0<Rc<1 且R0<1 時,E 0 對于所有的 τ ≥0 都是局部漸近穩定的;當R0>1 時,E0不穩定;當 Rc>1 且R0<1 時,E 0 在 τ=τ0時發 生 Hopf分 支.對 謠 言 盛 行 平 衡 點 E1*,若R0>1,m4>0,f2≥0,謠言盛行平衡點對所有 的 τ ≥0 都 是 局 部 漸 近 穩 定 的;若R0>1,m4>0,f2<0,謠言盛行平衡點在 τ=時發生 Hopf分支.最后,通過數值模擬,本文驗證了理論分析的可靠性并探討了個別參數對模型的影響:恢復率 r 的值越大,當R0<1 時,I (t) 趨于穩定的時間越短,速度越快,當R0>1 時,I (t) 所能達到的峰值越小,謠言傳播的范圍更小,更容易控制;抑制效應 α 越大,當R0<1 時,I (t) 趨于穩定的時間越長,速度越慢,當R0>1 時,I (t) 所能達到的峰值越高;恢復率r不斷增大,基本再生數 R0的值不斷減小;當恢復率r逐漸增大并增大到一定程度時,R0小于 1,意味著網絡謠言最終消亡.相反地,隨著出生率 Λ 的增大,基本再生數R0的值不斷增大,也就是說出生率 Λ 與基本再生數R0呈正相關關系;并給出了R0<1 與R0>1 時對應的出生率Λ與恢復率r的取值范圍;最后驗證了發生后向分支的條件.當前,本文主要分析了辟謠機制和時滯效應對社交網絡謠言傳播的影響,在未來的工作中將進一步考慮個體行為差異和空間異質性對網絡謠言傳播造成的影響[32],完善相關的網絡謠言傳播模型.

猜你喜歡
模型系統
一半模型
Smartflower POP 一體式光伏系統
工業設計(2022年8期)2022-09-09 07:43:20
WJ-700無人機系統
ZC系列無人機遙感系統
北京測繪(2020年12期)2020-12-29 01:33:58
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
基于PowerPC+FPGA顯示系統
半沸制皂系統(下)
連通與提升系統的最后一塊拼圖 Audiolab 傲立 M-DAC mini
3D打印中的模型分割與打包
主站蜘蛛池模板: 97视频精品全国在线观看| 一级黄色欧美| 国产精品区网红主播在线观看| 精品少妇人妻av无码久久| 亚洲福利网址| 日本伊人色综合网| 色天堂无毒不卡| 99久久精品国产自免费| 亚洲无码精彩视频在线观看| 国产成人1024精品下载| 黄色福利在线| 亚洲九九视频| 日韩不卡高清视频| 在线视频一区二区三区不卡| 国产网友愉拍精品视频| 黄色三级毛片网站| 久久久国产精品无码专区| 国产成人a毛片在线| 国产91av在线| 国内毛片视频| 在线观看亚洲人成网站| 九色免费视频| 国产激爽大片在线播放| 色综合久久久久8天国| 亚洲午夜国产精品无卡| 久久免费看片| 青青青伊人色综合久久| 久久国产香蕉| 色综合久久无码网| 91久久精品日日躁夜夜躁欧美| 99青青青精品视频在线| 呦系列视频一区二区三区| 九九久久精品国产av片囯产区| 日本一区二区不卡视频| 婷婷六月激情综合一区| 粉嫩国产白浆在线观看| 特级精品毛片免费观看| 在线日韩一区二区| 亚洲欧美激情小说另类| 色屁屁一区二区三区视频国产| 欧美日一级片| 亚洲区视频在线观看| 18黑白丝水手服自慰喷水网站| 国产又色又刺激高潮免费看| 国产区福利小视频在线观看尤物| 国模视频一区二区| 原味小视频在线www国产| 免费观看无遮挡www的小视频| 亚洲人成网站色7777| 亚洲国产在一区二区三区| 欧美成人一区午夜福利在线| 91www在线观看| 青青操视频免费观看| 精品国产免费观看| 97超爽成人免费视频在线播放| 日韩精品成人在线| 在线视频一区二区三区不卡| 99无码中文字幕视频| 欧美在线网| 91毛片网| 第一区免费在线观看| 亚洲第一视频网| 丝袜久久剧情精品国产| 欧美国产综合视频| 六月婷婷精品视频在线观看| 婷婷亚洲综合五月天在线| 熟妇丰满人妻av无码区| 国产一区二区三区夜色| 一本大道视频精品人妻| 久久中文电影| 永久免费无码日韩视频| 毛片视频网址| 激情六月丁香婷婷四房播| 欧美日在线观看| 永久毛片在线播| 特级毛片免费视频| 欧美伦理一区| 午夜少妇精品视频小电影| 少妇精品久久久一区二区三区| 日韩小视频网站hq| 国产成人精品亚洲日本对白优播| 久久96热在精品国产高清|