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

一類隨機Amensalism種群模型的動力學行為

2022-04-29 00:00:00艾曉輝,鄧文波,李鵬喆
哈爾濱理工大學學報 2022年3期

摘要:針對受環境因素影響的偏害共生動力學問題,依據隨機微分方程理論,采用引入白噪聲的方法,建立了隨機偏害共生動力學模型,并研究了隨機偏害共生種群模型的動力學行為。首先證明了隨機偏害共生動力學模型全局正解的存在性,其次利用伊藤公式和切比雪夫不等式證明了隨機偏害共生動力學模型解的隨機最終有界性,再次分析了隨機偏害共生動力學模型解的矩的漸近性質,最后通過數值模擬驗證了理論結果。

關鍵詞:隨機模型;全局正解的存在性;隨機最終有界性;矩的漸近性

DOI:10.15938/j.jhust.2022.03.019

中圖分類號: O211.63文獻標志碼: A文章編號: 1007-2683(2022)03-0140-07

Dynamic Behavior of a Stochastic Amensalism Population Model

AI Xiao-hui,DENG Wen-bo,LI Peng-zhe

(School of Science, Northeast Forestry University, Harbin 150040, China)

Abstract:Aiming at the problem of biased symbiosis dynamics affected by environmental factors, according to the theory of stochastic differential equation, the stochastic biased symbiosis dynamics model is established by introducing white noise, and the dynamic behavior of the stochastic biased symbiosis population model is studied. Firstly, the existence of global positive solution of the stochastic partial hazard symbiosis dynamics model is proved. Secondly, the stochastic ultimate boundedness of the solution of the stochastic partial hazard symbiosis dynamics model is proved by using the It formula and Chebyshev inequality. Thirdly, the asymptotic properties of the moment of the solution of the stochastic partial hazard symbiosis dynamics model are analyzed. Finally, the theoretical results are verified by numerical simulation.

Keywords:stochastic models; global positive solutions; stochastic ultimate boundedness; asymptotic property of moment

0引言

生物數學是生命科學與數學之間的邊緣學科[1],該學科的發展歷史悠久,在16世紀,我國科學家徐光啟[2]利用數學估算人口增長,這是最早利用數學解決種群問題的先例,在1798年,英國人口學家Malthus[2]提出了Malthus人口指數模型,并通過此模型來預測人口的增長規律,但在現實生活中資源和生存環境是有限的,人口的增長不可能是無限的,因此在1838年,Verhulst[2]提出了著名的Logistic模型,在20世紀40年代, Lotka[2] 和Volterra[2]提出了著名的Lotka-Volterra捕食者模型,對于該模型的研究有諸多成果[3-10],并且該模型至今仍是生物數學領域最經典和最重要的模型。

在20世紀,越來越多的學者通過數學對生物科學進行研究,并且以數學模型來解決生態學和生物學問題,使人們對生物發展規律的認識更加全面和深刻,也使得生物數學作為重要學科得到了長足發展[11],因此數學建模逐漸變為研究種群模型不可缺少的部分。從20世紀80年代開始,國內對于生物數學的研究越來越多,因此使我國的生物數學的研究水平不斷獲得提高。

然而在現實生活中,對于任何一個或多個種群的發展都會受到各種環境因素的影響,這些環境因素通常被稱為環境噪聲,它會對增長率、競爭系數等其他參數造成影響[12],因此如果將這些環境噪聲會忽略,會使研究的結果與現實情況會有一些偏差,因此為了研究結果更加準確,建立隨機生物種群是十分重要的。對于隨機生物種群的研究,首先是在確定性模型中加入環境噪聲干擾項,再通過隨機微分方程的相關知識進行研究[13]。這已經成為種群模型研究的熱點,并取得了一定的成果[14-20]。

1模型簡介

經典的非自治Logistic模型表示如下

dx(t)dt=x(t)[r(t)-a(t)x(t)](1)

其中:r(t)表示增長率;r(t)/a(t)表示時刻環境容量,a(t)和r(t)為連續有界函數。

在2016年,熊煥煥等[21]學者通過給典型的Logistic模型增添偏害作用項,創建了與Lotka-Volterra相似的偏害共生動力學模型

dN1dt=r1N11-N1P1-uN2P2

dN2dt=r2N21-N2P2(2)

其中N1≥0,N2≥0,r1gt;0,r2gt;0,P1gt;0,P2gt;0,ugt;0。

作者研究了該模型各奇點的穩定性,同時說明了每個奇點代表的生物意義,并且得出該系統存在唯一穩定點的結論。

但由于任何物種發展都會受到環境噪聲影響,我們在系統(2)中考慮環境中的白噪聲干擾,于是可以用

ri→ri+λiBi(t),i=1,2。

來代替系統(2)中的增長率ri(i=1,2),其中Bi(t)為標準白噪聲,λ2i為隨機干擾強度。從而建立出如下隨機偏害共生動力學模型

dN1=r1N11-N1P1-uN2P1dt+

λ1N11-N1P1-uN2P1dB1(t)

dN2=r2N21-N2P2dt+λ2N21-N2P2dB2(t)(3)

式中:N1≥0,N2≥0,r1gt;0,r2gt;0,P1gt;0,P2gt;0,ugt;0,λ1gt;0,λ2gt;0,B1(t)和B2(t)是獨立的標準Brownian運動。

2主要結果

2.1隨機偏害共生動力學模型的全局正解

接下來對系統(3)的全局正解存在性進行證明。

定理1對于任意給定的初值條件N(0)=(N1(0),N2(0))∈ 瘙 綆 2+,系統(3)具有唯一解N(t)=(N1(t),N2(t))在t≥0上定義,并且此解依概率1停留在 瘙 綆 2+中。

證明: 令

F(N(t))=(f1(N(t)),f2(N(t))),

G(N(t))=(g1(N(t)),g2(N(t))),

f1(N(t))=r1N1(t)1-N1(t)P1-uN2(t)P1,

g1(N(t))=λ1N1(t)1-N1(t)P1-uN2(t)P1,

f2(N(t))=r2N2(t)1-N2(t)P2,

g2(N(t))=λ2N2(t)1-N2(t)P2。

因為系統(3)的系數不滿足線性增長條件,但滿足局部Lipschitz條件。

對任意給定的初值N(0)=(N1(0),N2(0))∈ 瘙 綆 2+,在[0,τe)上存在唯一局部解N(t)=(N1(t),N2(t)),這里τe是爆破時間,只需證明τe=∞。

令P0充分大,使N1(0),N2(0)都在區間1P0,P0中,對每一個整數ngt;n0,定義停時

τP=inft∈[0,τe):Ni(t)1P,p,i=1,2

在這里我們設inf=+∞(代表空集),顯然當p→+∞時,τp單增。

令τ∞=limP→+∞τp,這里τ∞≤τea.s.,如果能證明τ∞=∞a.s.,那么τe=∞a.s.,并且N(t)∈ 瘙 綆 2+,tgt;0,即只需證明τ∞=∞a.s.,否則Tgt;0和ε∈(0,1),使得

P(τ∞≤T)gt;ε

因此P1gt;P0,使得對所有Pgt;P1,有

P(τp≤T)≥ε(4)

定義一個C2( 瘙 綆 2+)上的函數V: 瘙 綆 2+→ 瘙 綆 +如下:

V(Y)=(N1-1-0.5lnN1)+(N2-1-0.5lnN2)

這里當Nigt;0,i=1,2,有

Ni-1-0.5lnNi≥0,i=1,2

如果N(t)=(N1(t),N2(t))∈ 瘙 綆 2+,由It公式有

dV(N(t))=∑2i=1VNidNi+12VNiNi(dNi)2=

(0.5N-0.51-0.5N-11)r1N11-N1P1-uN2P1dt+

λ1N11-N1P1-uN2P1dB1(t)+0.5(-0.25N-1.51+

0.5N-21)λ21N211-N1P1-uN2P12dt+(0.5N-0.52-

0.5N-12)r2N21-N2P2dt+λ2N21-N2P2dB2(t)+

0.5(-0.25N-1.52+0.5N-22)λ22N221-N2P22dt=

(0.5N0.51-0.5)r11-N1P1-uN2P1dt+

λ11-N1P1-uN2P1dB1(t)+0.5(-0.25N0.51+

0.5)λ211-N1P1-uN2P12dt+(0.5N0.52-

0.5)r21-N2P2dt+λ21-N2P2dB2(t)+

0.5(-0.25N0.52+0.5)λ221-N2P22dt(5)

其中擴散算子

LV=(0.5N0.51-0.5)r11-N1P1-uN2P1+

0.5(-0.25N0.51+0.5)λ211-N1P1-uN2P12+

(0.5N0.52-0.5)r21-N2P2+0.5(-0.25N0.52+

0.5)λ221-N2P22

因此Mgt;0,使得

LV≤M。(6)

把式(6)代入(5)得

dV(N(t))≤M+(0.5N0.51-0.5)λ11-N1P1-

uN2P1dB1(t)+(0.5N0.52-0.5)λ21-N2P2dB2(t)

∫τpΛT0dV(N)≤∫τpΛT0Mdt+∫τpΛT0(0.5N0.51-

0.5)λ11-N1P1-uN2P1dB1(t)+

∫τpΛT0(0.5N0.52-0.5)λ21-N2P2dB2(t)(7)

在式(7)兩邊取均值,得

E[V(N(τpΛT))]≤V(N(0))+ME(τpΛT)≤

V(N(0))+MT

E[V(N1(τpΛT),N2(τpΛT))]≤

V(N1(0),N2(0))+MT(8)

記ΩP={τp≤T},那么由式(4)推出P(ΩP)≥ε,注意到對每一個ω∈Ωp,N(τp,ω)要么等于P要么等于1P,因此V(N(τp,ω))不小于

min{P-1-0.5lnP,1P-1+0.5lnP}

利用式(8)得

V(N(0))+MT≥E[1ΩP(ω)V(N(τp))]≥

εminP-1-0.5lnP,1P-1+0.5lnP

其中1ΩP為ΩP的指標函數。

令P→+∞就得到∞gt;V(N(0))+MT=∞,矛盾。因此一定有τ∞=∞a.s.。

2.2隨機偏害共生動力學模型解的最終有界性

本小結對系統(3)解的隨機最終有界性進行證明。

假設1mini=1,2ri-12(1-θ)λ2igt;0。

定義1系統(3)的解被稱為隨機最終有界的,如果對ε∈(0,1),χ=χ(ω)gt;0,使得對任意初值(k1,k2)∈ 瘙 綆 2+,系統(3)的解N(t)=(N1(t),N2(t))有以下性質:

lim supt→+∞{|N(t)|=N21(t)+N22(t)gt;χ}lt;ε

接下來給出一個與最終有界性的引理。

引理1對于r1,P1,u,r2,P2,λ1,λ2gt;0和θ∈(0,1),存在正常數M=M(θ),獨立于初值N(0)=(N1(0),N2(0)),使得在tgt;0上系統(3)的解N(t)=(N1(t),N2(t))有如下的性質:

lim supt→+∞E|N(t)|θ≤M

證明:定義

V(N(t))=∑2j=1Nθj,N(t)∈ 瘙 綆 2+

由It公式,我們有

dV(N(t))=∑2i=1VNidNi+12VNiNi(dNi)2=

θNθ-11r1N11-N1P1-uN2P1dt+

λ1N11-N1P1-uN2P1dB1(t)+

12θ(θ-1)Nθ-21λ21N211-N1P1-uN2P12dt+

θNθ-12r2N21-N2P2dt+λ2N21-N2P2dB2(t)+

12θ(θ-1)Nθ-22λ22N221-N2P22dt=

θNθ1r11-N1P1-uN2P1+

12θ(θ-1)Nθ1λ211-N1P1-uN2P12+

θNθ2r21-N2P2+12θ(θ-1)Nθ2λ221-N2P22dt+

θNθ1λ11-N1P1-uN2P1dB1(t)+

θNθ2λ21-N2P2dB2(t)(9)

定義:

LV=θNθ1r11-N1P1-uN2P1+

12θ(θ-1)Nθ1λ211-N1P1-uN2P12+

θNθ2r21-N2P2+12θ(θ-1)Nθ2λ221-N2P22

則由已知與假設2.2.1,M1gt;0使得

LV≤θr1Nθ1+12θ(θ-1)λ21Nθ1+θr2Nθ2+12θ(θ-1)λ22Nθ2=

θr1-12(1-θ)λ21Nθ1+θNθ2+θr2-12(1-θ)λ22Nθ2=

θ1θ+r1-12(1-θ)λ21Nθ1+θ1θ+r2-12(1-θ)λ22Nθ2-

(Nθ1+Nθ2)≤M1-V(N)

再由It公式,我們有

d[etV(N)]=etV(N)dt+etdV(N)=

et[V(N)dt+dV(N)]≤

etM1dt+θλ11-N1P1-uN2P1×Nθ1dB1(t)+

θλ21-N2P2Nθ2dB2(t)

對上式兩邊從0到t積分并同時取期望,得

E[V(N)]≤e-tV(N(0))+M1(1-e-t)

進而有

lim supt→+∞EV(N)≤M1

再由

|(N1,N2)|θ≤2θ2max{Nθ1,Nθ2}≤2θ2V(N1,N2)

lim supt→+∞E|(N1,N2)|θ≤2θ2lim supt→+∞EV(N1,N2)≤2θ2M1

令M=2θ2M1,則結論得證。

定理2系統(3)的解是隨機最終有界的。

證明: 由引理1,令θ=12,則M2gt;0,使得

lim supt→+∞E|(N1(t),N2(t))|12≤M2

于是對εgt;0,令χ=M22ε2,由Chebyshev不等式有

P(|(N1(t),N2(t))|gt;χ)≤E|(N1(t),N2(t))|12/χ

lim supt→+∞P(|(N1(t),N2(t))|gt;χ)≤ε

因此結論得證。

2.3隨機偏害共生動力學模型解的矩的漸近性質

接下來研究系統(3)解的矩的漸近性質。

定理3對任意給定的θ∈(0,1),Q=Q(θ),使得系統(3)的解N(t)=(N1(t),N2(t))滿足:

lim supt→+∞1t∫t0E[N2+θ1(s)+N2+θ2(s)]ds≤Q

證明:定義V: 瘙 綆 2+→ 瘙 綆 +為

V(N(t))=∑2j=1Nθj,N(t)∈ 瘙 綆 2+

則可由引理1得

LV(N(t))=θr11-N1P1-uN2P1Nθ1+

12θ(θ-1)λ211-N1P1-uN2P12Nθ1+

θr21-N2P2Nθ2+12θ(θ-1)λ221-N2P22Nθ2

進而

LV(N(t))+14θ(1-θ)λ21Nθ+21+14θ(1-θ)λ22Nθ+22≤

θr1Nθ1+12θ(θ-1)λ21Nθ1+θr2Nθ2+

12θ(θ-1)λ22Nθ2+14θ(1-θ)λ21Nθ+21+

14θ(1-θ)λ22Nθ+22=

θr1Nθ1+θr2Nθ2-14θ(1-θ)λ21Nθ+21-

14θ(1-θ)λ22Nθ+22≤M1

取λ2=minj=1,2{λ2j},有

LV(N(t))+14θ(1-θ)λ2Nθ+21+14θ(1-θ)λ2Nθ+22≤

LV(N(t))+14θ(1-θ)λ21Nθ+21+

14θ(1-θ)λ22Nθ+22≤M1(10)

其中M1是一個正數。

把式(10)代入式(9),并兩邊取期望,得

E[V(N)]+14θ(1-θ)λ2∫t0E[Nθ+21(s)+Nθ+22(s)]ds≤

V(N1(0),N2(0))+M1t

∫t0E[Nθ+21(s)+Nθ+22(s)]ds≤

4θ(1-θ)λ2[V(N1(0),N2(0))+M1t](11)

令Q=Q(θ)4θ(1-θ)λ2,在式(11)兩邊同除t,并令t→+∞,我們可以得到:

lim supt→+∞1t∫t0E[Nθ+21(s)+Nθ+22(s)]ds≤Q

因此結論得證。

3數值模擬

首先通過Milstein方法[22]來模擬上面關于隨機偏害共生動力模型 (3) 所得到的結果。

考慮系統 (3) 相應的離散化系統

N1(k+1)=N1(k)+r1N1(k)1-N1(k)P1-uN2(k)P1Δt+

λ1N1(k)1-N1(k)P1-uN2(k)P1Δtεk+

λ212N1(k)1-N1(k)P1-uN2(k)P11-2N1(k)P1-

uN2(k)P1(ε2k-1)Δt

N2(k+1)=N2(k)+r2N2(k)1-N2(k)P2Δt+

λ2N2(k)1-N2(k)P2Δtηk+

λ222N2(k)1-N2(k)P21-2N2(k)P2(η2k-1)Δt

其中εk和ηk(k=1,2,…,n)都是服從N(0,1)的Gaussian隨機變量。

使用上述的方法和R語言,選取適當的參數,我們給出隨機系統(3)的數值仿真。對系統(3)在初始條件N1(0)=N10,N2(0)=N20下進行數值模擬,選取N10=5000,N20=2000,r1=0.5,r2=0.3,u=2,P1=8000,P2=5000,λ1=0.5,λ2=0.5,并用R語言進行求解,得到圖1,模擬結果表明,該模型全局正解是存在的,即驗證了該模型全局正解的存在性。

對系統(3)在初始條件N1(0)=N10,N2(0)=N20下進行數值模擬,選取N10=7,N20=3,r1=0.5,r2=0.3,u=2,P1=8,P2=5,λ1=0.5,λ2=0.5,θ=0.6,M=3.8,并用R語言進行求解,得到圖2,圖中E=E|N(t)|θ,模擬結果表明在tgt;0下E的值小于等于M,因此lim supt→+∞E|N(t)|θ≤M成立,即引理1成立。

對系統(3)在初始條件N1(0)=N10,N2(0)=N20下進行數值模擬,選取N10=7,N20=4,r1=0.5,r2=0.3,u=1.65,P1=8,P2=5,λ1=0.5,λ2=0.5,χ=5.02,ε=0.8,并用R語言進行求解,得到圖3,圖中P=P{|N(t)|=N21(t)+N22(t)gt;χ},結果表明P在t→+∞下的值均小于ε,則lim supt→+∞P{|N(t)|=N21(t)+N22(t)gt;χ}lt;ε成立,即定義1成立。

對系統(3)的初始條件N1(0)=N10,N2(0)=N20下進行數值模擬,選取N10=7,N20=3,r1=0.5,r2=0.3,u=2,P1=8,P2=5,λ1=0.5,λ2=0.5,θ=0.6,Q=176.7,并用R語言進行求解,得到圖4,圖中N=E[N(θ+2)1(s)+N(θ+2)2(s)],模擬結果表明N在t∈[0,+∞)上的值均小于等于Q,則有lim supt→+∞1t∫t0E[N2+θ1(s)+N2+θ2(s)]ds≤Q成立,即定理3成立。

參 考 文 獻:

[1]徐克學. 試論生物數學的特點與展望[J].生物數學學報,1986(2):73.

XU Kexue. Some Viewpoint to the Characters and Prospects of Biomathematics[J]. Journal of Biological Mathematics, 1986(2):73.

[2]馬知恩. 種群生態學的數學建模與研究[M].合肥:安徽教育出版社,1996.

[3]HIROTSUGU M, NAOFUMI O, AKIRA S, et al. Statisticalmechanics of Population: the Lattice Lotka/Volterra Model[J]. Progress of Theoretical Physics,1992(6):1035.

[4]王暉,李冬梅,郭秀微,一類具有時滯的Lotka-Volterra捕食系統的持久性與全局穩定性[J].哈爾濱理工大學學報,2009,14(6):73.

WANG Hui, LI Dongmei, GUO Xiuwei. Persistence and Global Stability of Lotka-Volterra Predatorprey System with Time Delay[J].Journal of Harbin University of Science and Technology,2009,14(6):73.

[5]崔景安. 時滯Lotka-Volterra系統的持久性和周期解[J].數學學報,2004,47(3):511.

CUI Jingan. Permanence and Periodic Solution of Lotka-Volterra System with Time Delay[J]. ACTA MATHEMATICA SINICA,2004,47(3):511.

[6]王暉,李冬梅,郭秀微. 一類具有時滯的Lotka-Volterra捕食系統的持久性與全局穩定性[J].哈爾濱理工大學學報, 2009(6):73.

WANG Hui, LI Dongmei, GUO Xiuwei. Persistence and Global Stability of Lotka-Volterra Predatorprey System with Time Delay[J]. Journal of HarbinUniversity of Technology,2009(6):73.

[7]呂敬亮,王克. 隨機的改進Lotka-Volterra競爭模型[J].數學學報(中文版),2011,54(5):853.

LV Jingliang, WANG Ke. A Stochastic Modified Lotka-Volterra Competition Model[J]. ACTA MATHEMATICA SINICA,Chinese Series,2011,54(5):853.

[8]WAN Li, ZHOU Qinghua. Stochastic Lotka-Volterra Model with Infinite Delay[J]. Statistics amp; Probability Letters, 2009,79:698.

[9]BAHAR A, MAO Xuerong. Stochastic Delay Lotka-Volterra Model[J]. Journal of Mathematical Analysis amp; Applications,2004,292(2):364.

[10]BAO Jianhai, MAO Xuerong, YIN Geroge, et al. Competitive Lotka-Volterra Population Dynamics with Jumps[J]. Nonlinear Analysis,2011,74(17):6601.

[11]陳蘭蓀.數學生態學模型與研究方法[M].北京:科學出版社,1988.

[12]王克. 隨機生物數學模型[M].北京:科學出版社,2010.

[13]龔光魯. 隨機微分方程及其應用概要[M]. 北京:清華大學出版社, 2008.

[14]GARD TC. Stability for Multispecies Population Models in Random Environments[J]. Nonlinear Analysis: Theory, Methods amp; Applications, 1986, 10(12):1411.

[15]JIANG Daqing, SHI Ningzhong. A Note on Nonautonomous Logistic Equation with Random Perturbation[J]. Journal of Mathematical Analysis amp; Applications, 2005, 303(1):164.

[16]LIU Meng, WANG Ke. Survival Analysis of Stochastic Singlespecies Population Models in Polluted Environments[J]. Ecological Modelling,2009,220(9/10):1347.

[17]趙亞男,高海音. 具有隨機擾動的廣義\"食物有限\"種群模型正解的全局吸引性[J].吉林大學學報(理學版),2011,49(2):263.

ZHAO Yanan, GAO Haiyin. Global Attractivity of Positive Solutions to General “Food-Limited” Species Model with Random Perturbation[J].Journal of Jilin University (Science Edition),2011,49(2):263.

[18]孟笑瑩,鄧飛其,彭云建. 具有隨機擾動的食餌捕食系統的穩定性[J].系統工程與電子技術,2011,33(2):385.

MENG Xiaoying, DENG Feiqi, PENG Yunjian. Stability of a Prey-predator System with Random Perturbation[J].Systems Engineering and Electronics,2011,33(2):385.

[19]劉蒙. 隨機種群模型若干性質的研究[D].哈爾濱:哈爾濱工業大學,2012.

[20]趙愛民,趙曉丹,劉桂榮. 一類隨機互惠種群模型的漸近行為[J].山西大學學報(自然科學版),2019,42(2):281.

ZHAO Aimin, ZHAO Xiaodan, LIU Guirong. Asymptotic Properties of a Stochastic Mutualism Population Model[J]. Journal of Shanxi University(Nat.sci.Ed.),2019,42(2):281.

[21]熊煥煥,王斌斌,張海亮.魚群偏害共生動力學模型的穩定性分析[J].應用數學進展,2016,5(2):255.

XIONG Huanhuan,WANG Binbin,ZHANG Hailaing.Stability Analysis on the Dynamic Model of Fish Swarm Amensalism[J].Advances in Applied Mathematics,2016,5(2):255.

[22]HIGHAM D J. An Algorithmic Introduction Tonumerical Simulation of Stochastic Differential Equations[J].SIAM Rev,2001,43(3):525.

(編輯:溫澤宇)

主站蜘蛛池模板: 熟女日韩精品2区| 亚洲欧美成人在线视频| 日韩在线网址| 成人精品区| 欧美性精品| 免费可以看的无遮挡av无码 | 狠狠操夜夜爽| 精品一区二区三区四区五区| 国产主播福利在线观看| 高清无码手机在线观看| 九九热免费在线视频| 亚洲第一黄色网址| 真人高潮娇喘嗯啊在线观看| 亚洲成人播放| 欧美日在线观看| 被公侵犯人妻少妇一区二区三区| 日本国产精品| 2022精品国偷自产免费观看| 幺女国产一级毛片| 97在线视频免费观看| 日本国产精品| 国产白浆在线观看| 成人一级免费视频| 高清无码一本到东京热| 在线中文字幕网| 99精品免费欧美成人小视频| 性色一区| 午夜a视频| 任我操在线视频| 欧美色99| 午夜国产精品视频黄| 日韩美女福利视频| 免费国产黄线在线观看| jizz国产视频| 亚洲精品在线观看91| AV老司机AV天堂| 国产成人高清精品免费5388| 色久综合在线| 国产精品55夜色66夜色| 女人av社区男人的天堂| 日韩黄色大片免费看| 免费国产不卡午夜福在线观看| 国产xx在线观看| 一区二区在线视频免费观看| 99草精品视频| 亚洲欧美成人网| 欧美精品1区2区| 999精品视频在线| 91精品情国产情侣高潮对白蜜| 中文字幕资源站| 国产精品美女网站| 精品国产免费观看| 欧美一级视频免费| 97亚洲色综久久精品| 人妻21p大胆| 手机成人午夜在线视频| 日韩天堂网| 精品成人一区二区三区电影| 一级一毛片a级毛片| 亚洲女同欧美在线| 国产成人高清亚洲一区久久| 乱色熟女综合一区二区| 色综合久久久久8天国| 久久天天躁夜夜躁狠狠| 女人18毛片久久| 国产成人精品一区二区秒拍1o| 一级毛片不卡片免费观看| 欧美亚洲欧美| 成色7777精品在线| 综合成人国产| 亚洲一区二区三区国产精华液| 国产欧美在线观看一区| 99久久精品国产综合婷婷| 99九九成人免费视频精品| 亚洲欧美综合在线观看| 激情在线网| 欧美色图第一页| 中文字幕第4页| 91无码视频在线观看| 大乳丰满人妻中文字幕日本| 毛片在线播放网址| 国产视频入口|