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

基于改進SEIR模型的COVID-19疫情度量分析

2023-03-11 05:01:48馬力文張亞琪
計算機仿真 2023年1期
關鍵詞:疫情模型

馬力文,周 穎,張亞琪

(1. 南京郵電大學自動化人工智能學院,江蘇 南京210023;2. 南京郵電大學江蘇先進材料協作研究中心,江蘇 南京 210023)

1 引言

隨著疫情席卷全球,破解COVID-19疫情的發展趨勢和傳播模式迫在眉睫:比如病毒的潛伏期多長,潛伏期是否具有傳染性,傳染強度有多大?根據目前的報道可知病毒的傳染源是攜帶COVID-19病毒的患者(包括潛伏期患者),主要的傳播途徑是飛沫傳播,人群普遍易感[1-3]。

SEIR模型是基于復雜網絡的經典傳染病模型。主要的復雜網絡傳染病模型包括SI模型,SIR模型和SEIR模型。SI模型將人群分為易感人群和感染人群。易感人群在單位時間內以一定的傳染概率受到感染轉化為感染人群。該模型并沒有考慮感染者被治愈或死亡的情況,與現實生活中的實際情況不相符。SIR模型增加了移出狀態(R),表征感染者治愈或死亡的狀態變化,感染人群在單位時間內以一定的治愈或死亡概率由感染狀態轉移到移出狀態。該模型已經可以表征很多傳染病的發展過程,但根據報道可知,潛伏期的COVID-19病毒患者同樣具有傳染性。由此提出具有潛伏狀態(E)的SEIR模型,潛伏期的患者在單位時間內以概率γ轉移至感染狀態[4-9]。基于SEIR模型的現有研究表明,對攜帶COVID-19病毒的患者進行有效隔離并及時進行病原追蹤可以有效遏制疫情大范圍擴散[10]。

COVID-19病毒的潛伏期在3-14天之間,傳染強度不如SARS病毒,但傳染范圍遠超SARS病毒[11-13]。文獻 [14]對國內疫情的省間傳播做了具體分析,文獻[15]建立了國內新冠肺炎傳播的智能體模型,文獻[16]采用周期性SEIR模型對傳染動力學進行建模。已有文獻主要是對國內新冠肺炎疫情的研究,且僅使用的SEIR模型仿真了四類人群的演變歷程,缺少對病毒潛伏期和傳染強度的仿真。當前新冠肺炎疫情的高風險地區已經從國內轉移到國外,各國疫情、民情、國情特征各不相同,COVID-19病毒的具體表征也各不相同[17],只分析國內的疫情信息,是遠遠不夠的。潛伏期長度和傳染強度是表征傳染病的重要指標,會隨著疫情發展而發生變化,想要高效防疫就必須掌握這兩個參數的變化。

本文選擇了中國(武漢),澳大利亞,印度,法國,韓國,美國幾個國家從2019年12月開始到2020年7月間的相關數據,分別進行了疫情特征描述和防疫效果度量。特征描述分為兩個部分:使用加入彈性變量的中位數模型求解平均潛伏期;基于SEIR模型狀態間的轉化關系求得病毒的基本傳染數。使用相應的回溯參數演繹各國防疫效果,具體分為三個參數判定指標:控制傳染源,切斷傳播途徑,保護易感人群。仿真結果認為中國、韓國、意大利、法國的防疫成果較好;澳大利亞的防疫成果一般;美國、印度的防疫成果較差。

2 COVID-19疫情特征描述

2.1 基于中位數模型的潛伏期計算

潛伏期的人(包括無癥狀感染者)具有傳染性且潛伏期是隨機的[18]。通過對感染者的跟蹤和回溯得到COVID-19的潛伏期長度,有較大的不確定性。建立估計潛伏期的模型并合理估計平均潛伏時間,有利于政策制定者更合理地安排防疫舉措。

2.1.1 模型假設

1)病毒的源頭為共同來源

2)統計的疫情期只有一個高峰

3)病例集中在最長和最短潛伏期之間。

2.1.2 模型的建立

傳染病潛伏期長度計算的方法有兩種:幾何均數法和中位數法。從官方新聞中得知,此次疫情的潛伏期以天為單位,是離散型計數數據,且潛伏期少有成指數模式增長的情況發生[2]。因此,本文采用中位數法。

由于病毒的傳播尚未結束,因而從世界衛生組織(WHO)網站以三天為一個步長,收集意大利、法國、澳大利亞、韓國、印度和中國從2020年1月24日到2020年7月13日的總確診人數和單日新增確診人數的數據。選取總確診數的16%,50%,84%作為中位數法選取數據的標志比例,將得到的首部,中部,尾部的總確診人數代入到模型中,估算出該國家的病毒潛伏期。

由于影響病例增長的因素錯綜復雜,因而在中位數模型中引入彈性變量,使得模型的適應性更強。平均潛伏期估算公式如下

(1)

由于各國對疫情的統計信息準確性不同,直接帶入中位數公式難以得到合適的估算值,因此引入彈性變量。彈性變量j為標志比例個數。

2.1.3 各國潛伏期估算結果

統計篩選WHO官網上公布的各國疫情數據,帶入式(1)中計算得到的潛伏期估算結果如下表所示。

表1 澳大利亞潛伏期估算結果

表2 意大利潛伏期估算結果

表3 法國潛伏期估算結果

表4 印度潛伏期估算結果

表5 韓國潛伏期估算結果

表6 中國潛伏期估算結果

由上表可知,各國的平均潛伏期分布在3.93-7.43天之間,總平均潛伏期為5.489天。本輪絕大多數COVID-19疫情患者感染病毒到發病只需經歷2-3天[19],而潛伏期長度的計算結果受各國確診病例統計數據的影響,所以潛伏期平均時間相對較短的國家發現潛伏病患的能力更強。由上述表格可以看出,中國在病例最多的情況下,病毒潛伏期估算值最小,說明中國的病毒檢測能力最強。法國,韓國和澳大利亞醫療基礎較強,潛伏期的估算值大小處于中間位置,病毒檢測能力較強。意大利和印度在疫情期間,確診病例增速先慢再突然爆發,且潛伏期的估算值較大,說明疫情初期,兩國的病毒檢測能力沒有跟上,大量得不到確診的病例引起了感染人數爆炸式的增長。

2.2 基于SEIR模型的病毒傳染強度仿真

從官方信息可知,攜帶COVID-19的潛伏者同樣具有感染健康人群的能力,因此,選用SEIR模型對人群結構進行建模。將總樣本分為易感者,潛伏者,感染者,移出者四部分。病毒的感染強度,政府的管控力度和民眾自我防范意識等因素都會影響四類人群的轉化率。本文通過計算各國不同時期四類人群轉化率的大小來仿真病毒傳染強度隨時間變化的規律,并使用基本傳染數R0作為衡量傳染強度大小的指標。

其中S為易感狀態的人群,為統計地區可能接觸到傳染源的全部人口;E為處于潛伏期的人群,表示已經被感染但沒有表現出感染癥狀的群體,具有傳染能力;I為感染人群,具有傳染能力;R為移出狀態人群,包括治愈,死亡,有效隔離的人,不具有病毒傳染能力且移出后不會再次被感染。四種類型人群的轉化方向如圖1所示。

記S(t),E(t),I(t),R(t)分別為t時刻的易感人群數,潛伏人群數,感染人群數,移出人群數[9]。N(t)為t時刻群體總人數,有

圖1 SEIR模型狀態轉化

S(t)+E(t)+I(t)+R(t)≡N(t)

(2)

假設S為易感人群總數,N為總人口數,易感人群在一定時間內被傳染的概率為β,比例為S/N,此時共有I(t)感染個體,易感人群的變化如下所示[20]。

(3)

假設感染者在一定時間內可以接觸k個人,接觸者被感染的概率為b,有

β=kb

(4)

相應地,潛伏個體的數目按照如下變化率增加,并且整體以單位時間概率γ1轉化為感染個體[20-21]。

(5)

潛伏期人群確診轉化為感染人群,同時感染人群治愈或死亡,以概率γ2轉化為移出人群[21]

(6)

相應的,感染群體以概率γ2由向移出群體轉化

(7)

綜上可得

(8)

基本傳染數R0為

(9)

其中,C為恢復期,C=14。

由于中國人口總數龐大,湖北地區與湖北以外地區疫情差距極大,為減少結果偏差,中國地區以武漢為研究對象。本次疫情的傳染周期一般為7天[22],故本文中將7天作為一組求b、R0的值,按照疫情爆發前期(1月25日至1月30日)、中期(2月25日至3月3日)、后期(3月8日至3月14日)選取武漢市衛生健康委員會公示的數據。2019年12月8日的確診病例為1人,為第0天,武漢人口為1300萬,確診病例的死亡率為3%,COVID-19的恢復期C為14天,總確診人數函數為I(t),有

I(0)=1

(10)

在疫情早期,由于防控意識淡薄,本題中假定感染者每天平均接觸人數k=5,系統中的全部人口都處于易感狀態。統計篩選武漢市衛生健康委公示的疫情數據,得到武漢市各疫情期總確診人數如下表所示。

表7 武漢市疫情前期的總確診人數

表8 武漢市疫情中期的總確診人數

表9 武漢市疫情晚期的總確診人數

表10 武漢傳染強度參數

經過計算得到表10中的數據,可以看出武漢地區病毒的傳染概率和傳染強度隨時間推移逐漸減弱。使用高斯煙羽模型[23]對病毒傳染能力進行可視化仿真

(11)

其中,C為空間點(x,y,z)的病毒傳染強度;σy和σz分別表示橫向和縱向上病毒傳染能力的擴散參數;μ為攜帶傳染源的患者平均移動距離;q為傳染源強度,x,y,z為空間中某點到點源的距離。得到可視化結果如下圖所示。

圖2 二維病毒傳染能力仿真

高斯煙羽模型中得到的仿真圖中,曲線上的數值為可傳染人數。可以看出距離傳染源30米的范圍內病原體的傳染能力很強,隨著社交距離的增加,病毒的傳染能力逐漸減弱。將病毒攜帶者和易感者有效隔離,可以大大削弱病毒的感染能力。

圖3 三維病毒傳染能力仿真

使用SEIR模型繼續求出各國病毒傳染強度變化:

表11 韓國傳染強度參數

表12 美國傳染強度參數

表13 意大利傳染強度參數

表14 法國傳染強度參數

表15 澳大利亞傳染強度參數

表16 印度傳染強度參數

表17 不同時期各國家病毒傳染強度

基本傳染數R0值越小說明疫情控制得越好,當R0遞減甚至小于1后,病毒的傳播就會逐漸消失。如圖4所示,各國R0的大小都在逐漸變小。說明從疫情爆發至今,各國政府的管控力度不斷增大,人民的防范意識不斷增強,病毒的傳染強度在人為控制下逐漸減弱了。

圖4 各國病毒傳染強度隨時間變化

中國在疫情初期,就對武漢封城,并在全國停工停產[19],很好地阻斷了病毒的傳播,所以中國的R0在各國中處于一個很低的水平。意大利早期的R0很大,但在政府決定封城后,中期的R0迅速變小,疫情發展逐漸向好。澳大利亞的R0一直處于很低的水平得益于廣大的領土和極少的人口,感染者可接觸的平均人數很低。而美國,在早期疫情爆發后,仍沒有有效的防疫舉措,可以看到雖然傳染強度在減弱,但防疫形式依然嚴峻。

3 防疫效果度量指標選取

以中國武漢為例,確定行之有效的防疫舉措。圖5中武漢已經封城,但前七天未對城內居民的行動進行過多限制。此時N=1100萬,C=14,K=5,R0=3.108。

由仿真結果可知,感染人群的曲線在第66天左右開始遞增,表明中國武漢的疫情在2 月6 日左右開始集中爆發,感染人群的曲線在第108天達到極大值,此時疫情的傳播到達高峰(3月末,四月初),潛伏期人群數在第150 天降為0,表明疫情的傳播接近尾聲(5月上旬)。

圖5 未采取措施的四類人群變化趨勢

圖6中,政府采取了一定的管控手段和人們自發的避免感染,會減小病毒的基本傳染數。此時N=1100萬,C=14,K=5,R0=2.6166。

疫情在第60天左右開始集中爆發(1月末),第95天左右達到高峰(3月中旬),第140天后接近尾聲(4月下旬)。這也說明早期武漢政府披露的感染人數可能比實際少得多,這使得初期的基本傳染數較大。引起這些問題的原因有:檢測跟不上,醫療設備缺乏,無法滿足檢測和治療的需求。

圖6 采取一定管控的手段后四類人群變化趨勢

圖7中,政府對武漢的交通進行了限制,并關閉了各種不必要的公共場所。此時保持傳染基本數不變,感染者每天平均接觸到的人數k減少。此時N=1100萬,C=14,K=1,R0=3.108。

疫情在第50天左右開始集中爆發(1 月25 日左右),第84天左右達到高峰(3月上旬),第120天后接近尾聲(4月中旬)。該結果表明,管控措施對疫情防治作用很明顯。

圖7 交通管制后四類人群變化趨勢

隨著時間推移,人們的防范意識和政府管控加強,基本傳染系數R0會減少。此時N=1100萬,C=14,K=1,R0=2.6166。

圖8 政府強化管制后四類人群變化趨勢

疫情在第48天左右開始集中爆發(1 月23 日左右),在第75天左右達到高峰(2月末),在第125天后接近尾聲(4月上旬)。與武漢實際情況比較相符,說明政府當前的管控力度非常強。

前面的分析驗證了中國疫情防控取得了很大成效;具體舉措何以分為三個方面:控制感染源,切斷傳播途徑,保護易感人群。

從控制感染源的角度來講,中國很早成功檢測并向世界公布了新冠病毒的基因序列,并大批量生產檢驗試紙。在疫情重點區域進行了大規模核酸檢測,提出了“應收盡收,不漏一人”的口號[23]。中國強大的科研能力和政府執行力使國家很快就控制住了傳染源。

從切斷傳播途徑的角度來看,在疫情初期,黨中央就迅速下令,對武漢進行封城,全國范圍內停運市內交通,并關閉人員密集型服務場所[23]。全國人民自發的在家中隔離,減少外出活動。這些措施都極大地阻斷了病毒從武漢向全國范圍內擴散,也極大地降低了潛伏者感染健康者的可能性。

從保護易感人群的角度來說,中國的方艙醫院對感染者進行了有效隔離,阻斷了二次感染[20],極大地降低了感染總人數,使疫情快速到達拐點,加快了疫情結束的腳步。

4 各國COVID-19防控評估

疫情數據非常繁雜,難以清晰明了地分析數據之間的關系和度量各國防疫效果。本節根據控制傳染源,切斷傳播途徑,保護易感人群三個指標,使用相應的回溯參數對各國疫情趨勢進行演繹,并將演繹結果與實際進行比較得出直觀的度量結論。

假設確診總人數為M,包括潛伏期人數和感染期人數,其中潛伏期人數為m,潛伏人群占總確診人數的比例為

(12)

其中,P由平均潛伏期內感染比例求得:

(13)

其中n為潛伏期仍具備感染能力(即未被有效隔離)的人數,N為總人數,T為從感染到確診的時間和。

由上面兩式可以回溯求得感染人數的變化過程。接下來將使用回溯參數結合控制傳染源,切斷傳播途徑,保護易感人群三個手段來演繹各國理想的防疫趨勢,并和實際趨勢進行比較來評價各國的防疫成果。使用表18中的數據進行回溯傳播建模。

表18 各國回溯參數數據

中國防疫工作的效果顯著,但圖9理想回溯的圖像和實際圖像并未完全重合。這是因為在疫情初期,人們對病毒的意識不足,政府的檢測能力較弱,并且存在大量的潛伏期病人仍在自由活動,因而在這一時期實際的感染人數應多于確診人數。這使得總確診人數略高于理想模型,疫情高峰、結束期的到來也相應延后。

控制感染源:美國每天可以確診數萬的感染病例,說明其具有強大的檢驗能力。切斷感染源:未停工停產,民眾不戴口罩且可以參加聚集性活動。保護易感人群:未建立方艙醫院,無法有效隔離易感人群。

圖9 中國(武漢)回溯仿真

從圖10中可以看出,易感人群的下降曲線推后了大概25天,這說明美國政府并沒有有效地切斷傳染源,實際的感染人數仍呈井噴式爆發,遠超每日新增的確診人數。直接導致疫情高峰相較于預測模型向后推遲了30天,結束日期則推遲了近50天。實際移出群體數量也大量減少,說明由于防疫失利,治愈率下降,死亡率升高,且大量民眾再確診前就已經死亡。如果美國仍不采取有效的防控措施,直至2020年年末,美國都不會迎來疫情拐點。

控制感染源:韓國政府的病毒檢測能力很強。切斷傳染源:未停工停產,但民眾普遍帶口罩。保護易感人群:因為較好的醫療基礎,韓國醫院對已經確診的人進行了較為有效的隔離。

圖10 美國回溯仿真

從圖11可以看出,韓國疫情的高峰期出現在3月15日前后,比預測大概推遲了15天左右,5月初疫情趨于穩定,實際趨勢與預測相差不大。這說明韓國政府的防控工作做的較好,民眾普遍帶口罩也極大地減弱了病毒的傳染,但未停工停產還是造成了部分地區性傳染事件的發生。韓國接下來需要重點防控小范圍區域傳染。

圖11 韓國回溯仿真

控制感染源:印度人口眾多,數量龐大的貧困人口得不到及時的檢測。切斷傳染源:占人口絕大多數的貧困地區是醫療和防控的真空區,且人口密度極大,印度政府無法真正地切斷傳染源。保護易感人群:印度基礎設施和醫療儲備很差,完全做不到有效隔離。

從圖12中可以看出,印度國內疫情爆發明顯滯后于其他國家,這說明印度政府的檢測速度遠遠落后于病毒的傳染速度,印度政府根本沒有辦法控制疫情的發展。占比超過90%的貧困人口預示印度的防疫措施將形同虛設,在接下來的時間,印度將一直處于疫情高峰期,直至政府發放免費的疫苗。

控制感染源:意大利政府在接受了中國政府的大量援助后,獲得了大量的檢驗試紙,檢驗能力較強。切斷感染源:停工停產,禁止了民眾的非必須外出活動,要求全民戴口罩。保護易感人群:學習了中國的方艙醫院,有效的阻斷了病毒的二次傳染。

圖12 印度回溯仿真

從圖13可以看出,易感人群的下降比預期滯后了一個月,意大利在一月出現輸入病例后,并沒有采取有效措施,例如馬拉松之類的大型活動就是病毒大范圍傳播的溫床,直至進入3月疫情爆發才開始提倡居家隔離。移出率比預測低了很多,這說明意大利出現了醫療擊穿的情況,死亡率高達6.6%,且意大利老齡化嚴重,超過80歲的人群死亡率超過14.4%[24]。意大利在得到中國的援助后,迅速封城,并建立了方艙醫院,有效的遏制了疫情的進一步惡化,在疫情開始后的110天(大概四月中旬)迎來了疫情的拐點。整體防疫趨勢向好。

控制感染源:作為老牌資本主義強國,法國政府的病毒檢驗能力很好。切斷傳染源:未停工停產,但鼓勵民眾戴口罩。保護易感人群:對感染的人群進行了有效隔離。

圖13 意大利回溯仿真

從圖14可以看出,疫情早期,法國仍然對申根國家開放邊境,輸入了大量流動的病毒攜帶者,疫情快速增長期出現在意大利宣布封國后的一周。但法國國力和醫療基礎實力雄厚,并沒有出現類似于意大利的死亡率整體的防疫效果不錯。在疫情發生的第100天左右(四月上旬),比意大利更早的迎來了疫情拐點。

控制感染源:澳大利亞政府病毒檢測能力較好。切斷感染源:未停工停產,民眾戴口罩的意識一般,但地廣人稀,有自身得天獨厚的優勢。保護易感人群:對感染的人群進行了有效隔離。

圖14 法國回溯仿真

從圖15可以看出,易感人群的下降時間,疫情拐點的到來,疫情結束時間都明顯滯后于預測模型,澳大利亞政府的疫情防控工作做得很一般,但是得益于自身較好的基礎醫療積累和地廣人稀的天然優勢,新冠肺炎并沒有爆炸式傳播。澳大利亞的移出人群比預期低,是治愈率較低。如果該國政府不提高防疫效率,來自美國的輸入病例就會引爆澳大利亞疫情爆發的火藥桶。

綜上所述,控制傳染源,切斷傳播途徑,保護易感人群是控制疫情的有效手段。檢測能力受國力影響很大,老牌強國的檢測能力都不錯,而印度受制于龐大的貧困人口,心有余而力不足。由于經濟原因,絕大多數國家無法像中國,意大利一樣,大面積停工停產,但減少不必要聚集活動和戴口罩是切斷傳播途徑的有效方法。保護易感人群可以極大地減緩感染人數的增加,在醫院隔離資源飽和的情況下,方艙醫院是有效隔離的最佳選擇。疫情防控不僅僅是政府工作,民眾配合也十分重要[25],中國民眾自發居家隔離,美國民眾不戴口罩且外出聚集,這都極大的影響了疫情發展。各國的實際情況同樣說明越及時采取符合國情的有效措施,防控效果越好。

圖15 澳大利亞回溯仿真

5 結束語

本文分析了病毒傳染能力的兩個指標:平均潛伏期和基本傳染數。各國平均潛伏期長度為5.489天,基本傳染數在3.93-7.43之間。使用改進SEIR模型確定了度量防疫效果的三個方向:控制傳染源,切斷傳播途徑,保護易感人群。使用相應的回溯參數對理想防疫趨勢進行演繹,通過和實際趨勢的比較對各國防疫效果進行客觀評價。本文認為在第一波疫情高峰中中國、韓國、意大利、法國的防疫成果較好;澳大利亞的防疫成果一般;美國、印度的防疫成果較差。

猜你喜歡
疫情模型
一半模型
戰疫情
重要模型『一線三等角』
抗疫情 顯擔當
人大建設(2020年5期)2020-09-25 08:56:22
疫情中的我
疫情當前 警察不退
北極光(2020年1期)2020-07-24 09:04:04
重尾非線性自回歸模型自加權M-估計的漸近分布
待疫情散去 春暖花開
文苑(2020年4期)2020-05-30 12:35:48
疫情期在家帶娃日常……
37°女人(2020年5期)2020-05-11 05:58:52
3D打印中的模型分割與打包
主站蜘蛛池模板: av在线无码浏览| 四虎国产在线观看| 热久久这里是精品6免费观看| 极品尤物av美乳在线观看| 日本a级免费| 国产一区三区二区中文在线| 国产视频大全| 无码AV动漫| 日韩免费成人| 亚洲乱码视频| 日韩精品专区免费无码aⅴ | 重口调教一区二区视频| 精品视频一区在线观看| 成人日韩精品| 欧美伦理一区| 亚洲国产系列| 四虎影视库国产精品一区| a级高清毛片| 亚洲成人网在线观看| 综合久久五月天| 伊人福利视频| 亚洲精品国产日韩无码AV永久免费网| 亚洲男人的天堂在线| 日韩中文字幕免费在线观看| 国产真实乱子伦视频播放| 中文字幕伦视频| 青青草欧美| 男女猛烈无遮挡午夜视频| 免费在线一区| 亚洲成年网站在线观看| 亚洲成av人无码综合在线观看| 99爱视频精品免视看| 国产乱人伦AV在线A| 欧美va亚洲va香蕉在线| 国产av一码二码三码无码| 国产白浆视频| 国产地址二永久伊甸园| 亚洲九九视频| a毛片免费在线观看| 亚洲中文字幕av无码区| 亚洲欧洲综合| 大陆精大陆国产国语精品1024| 一级毛片无毒不卡直接观看| 国产精品分类视频分类一区| 99热国产这里只有精品无卡顿"| 欧美区国产区| 四虎永久免费网站| 亚洲视频免| 国产精品自在自线免费观看| 高清不卡毛片| 在线观看国产精品第一区免费 | 亚洲码在线中文在线观看| 看看一级毛片| 国产爽妇精品| 亚洲精品爱草草视频在线| 国产综合色在线视频播放线视| 亚洲一级毛片免费观看| 99在线国产| 国产精品久久久久久久伊一| 精品一区二区三区视频免费观看| jizz亚洲高清在线观看| 亚洲日韩精品综合在线一区二区| 91视频国产高清| 国产国产人成免费视频77777| 欧美另类第一页| 国产成人1024精品| 久夜色精品国产噜噜| 国产成人综合在线观看| 色婷婷综合激情视频免费看 | a国产精品| 欧美伦理一区| 精品少妇人妻无码久久| 五月激激激综合网色播免费| 午夜毛片免费看| 国产最新无码专区在线| 婷婷午夜天| 国产鲁鲁视频在线观看| 在线看AV天堂| 国产99精品久久| 国产成人精品高清不卡在线| 一级做a爰片久久毛片毛片| 91热爆在线|