亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        色噪聲激勵(lì)的時(shí)滯非線性系統(tǒng)瞬態(tài)響應(yīng)研究

        2013-08-07 11:32:41戚魯媛高維廷
        關(guān)鍵詞:概率密度將式時(shí)滯

        戚魯媛,徐 偉,高維廷

        QI Luyuan1,XU Wei1,GAO Weiting2

        1.西北工業(yè)大學(xué) 理學(xué)院應(yīng)用數(shù)學(xué)系,西安 710129

        2.西北工業(yè)大學(xué),電子信息學(xué)院,西安 710129

        ◎博士論壇◎

        色噪聲激勵(lì)的時(shí)滯非線性系統(tǒng)瞬態(tài)響應(yīng)研究

        戚魯媛1,徐 偉1,高維廷2

        QI Luyuan1,XU Wei1,GAO Weiting2

        1.西北工業(yè)大學(xué) 理學(xué)院應(yīng)用數(shù)學(xué)系,西安 710129

        2.西北工業(yè)大學(xué),電子信息學(xué)院,西安 710129

        建立了色噪聲與時(shí)滯聯(lián)合作用的非線性系統(tǒng)模型,提出求解其瞬態(tài)概率密度的高效近似算法。利用等價(jià)變換將時(shí)滯系統(tǒng)簡(jiǎn)化為非時(shí)滯系統(tǒng);通過線性化方法和隨機(jī)平均原理得到原系統(tǒng)振幅過程的平均It?隨機(jī)微分方程和相應(yīng)的Fokker-Planck-Kolmogorov(FPK)方程?;谕嘶€性系統(tǒng)導(dǎo)出一組正交基,在該基空間內(nèi)進(jìn)行Galerkin變分得到近似瞬態(tài)概率密度。將該方法應(yīng)用到受色噪聲激勵(lì)的雙時(shí)滯Duffing-Van Der Pol振子得到理論解,采用蒙特卡羅模擬(MCS)驗(yàn)證理論解的正確性。分析了色噪聲參數(shù)和時(shí)滯參數(shù)對(duì)瞬態(tài)響應(yīng)的影響。研究結(jié)果表明:所提理論方法可有效求解受色噪聲激勵(lì)的時(shí)滯非線性系統(tǒng)的瞬態(tài)概率密度;算法求解效率高于MCS;色噪聲和時(shí)滯均明顯影響了系統(tǒng)瞬態(tài)響應(yīng)。

        時(shí)滯;色噪聲;FPK方程;瞬態(tài)概率密度;Galerkin算法

        1 引言

        不確定因素廣泛存在于圖像處理、自動(dòng)化、移動(dòng)通信、地球物理等研究領(lǐng)域[1-3]。理論上對(duì)不確定因素的描述主要分為兩類:一是將系統(tǒng)描述為不確定參數(shù)系統(tǒng);二是利用各種噪聲對(duì)系統(tǒng)進(jìn)行擾動(dòng)。研究證明幾乎所有的機(jī)械(結(jié)構(gòu))系統(tǒng)都具有一定程度的非線性性態(tài),非線性隨機(jī)系統(tǒng)的響應(yīng)預(yù)測(cè)已然成為理論研究和工程實(shí)踐的熱點(diǎn)問題[4]。眾多學(xué)者經(jīng)過長(zhǎng)期努力,發(fā)展了一系列預(yù)測(cè)系統(tǒng)響應(yīng)的理論方法,其中Fokker-Planck-Kolmogorov(FPK)方程方法是擴(kuò)散過程理論的主要方法。通過求解FPK方程得到穩(wěn)(瞬)態(tài)轉(zhuǎn)移概率密度,可用以分析系統(tǒng)控制[5]、信息熵[6]等實(shí)際問題。由于FPK方程的復(fù)雜性,其瞬態(tài)精確解只能對(duì)少數(shù)一階系統(tǒng)得到[7],一般意義的非線性系統(tǒng)的瞬態(tài)解只能結(jié)合理論分析與數(shù)值計(jì)算進(jìn)行近似求解。

        隨機(jī)平均原理是擴(kuò)散理論與確定性非線性振動(dòng)理論相結(jié)合的產(chǎn)物,是研究隨機(jī)非線性系統(tǒng)響應(yīng)問題的重要理論方法之一。利用隨機(jī)平均原理可有效避免非白噪聲帶來的FPK方程擴(kuò)維現(xiàn)象,并可降低FPK方程維數(shù),簡(jiǎn)化了理論分析和數(shù)值計(jì)算。Galerkin法是一種求解微分方程的變分法,1968年,Bhandari和Sherrer首次將Galerkin法應(yīng)用到求解FPK方程中,得到了穩(wěn)態(tài)概率密度[8]。繼而Wen應(yīng)用其求解了FPK方程的瞬態(tài)概率密度[9]。2007年,Spanos結(jié)合隨機(jī)平均理論和Galerkin法研究了白噪聲激勵(lì)的非線性系統(tǒng)的瞬態(tài)概率密度,且求解過程不牽扯攝動(dòng)原理,適用于含有較大非線性參數(shù)的系統(tǒng)[10]。

        時(shí)滯現(xiàn)象廣泛存在于計(jì)算機(jī)、控制等工程實(shí)際領(lǐng)域[11-13]。研究時(shí)滯作用對(duì)瞬態(tài)概率密度的影響具有重要的理論及應(yīng)用價(jià)值,并已經(jīng)引起相關(guān)學(xué)者的注意。文獻(xiàn)[14]將隨機(jī)平均原理和Galerkin法兩者結(jié)合分析了白噪聲作用的非線性時(shí)滯系統(tǒng)的瞬態(tài)概率密度。正如很多學(xué)者指出的,實(shí)際中存在的噪聲是具有一定相關(guān)時(shí)間的有色噪聲,故而研究色噪聲對(duì)瞬態(tài)概率密度的影響是研究此領(lǐng)域不可缺少的一部分??偨Y(jié)前人工作發(fā)現(xiàn),對(duì)色噪聲激勵(lì)的非線性系統(tǒng)瞬態(tài)概率密度的研究工作較少[15],研究時(shí)滯與色噪聲共同作用的非線性系統(tǒng)的瞬態(tài)概率密度仍處于空白階段。

        綜上所述,本文將等價(jià)線性化原理、隨機(jī)平均理論和Galerkin算法三者結(jié)合,求解色噪聲與時(shí)滯聯(lián)合作用的非線性振子的瞬態(tài)概率密度;系統(tǒng)地分析了該算法的有效性,并研究了色噪聲和時(shí)滯對(duì)系統(tǒng)瞬態(tài)概率密度的影響。

        2 模型構(gòu)造

        建立色噪聲激勵(lì)的時(shí)滯非線性系統(tǒng)模型如下:

        ζk(t)為零均值,強(qiáng)度為hk,譜密度函數(shù)為Sk(τ)的相互獨(dú)立的有色噪聲。假設(shè)存在一個(gè)δ小量且H和F是δ階小量,hk為δ1/2階小量。

        3 理論求解過程

        3.1 隨機(jī)平均原理

        對(duì)式(1)應(yīng)用隨機(jī)平均原理,引入如下變換:

        其中,A(t)、Φ(t)、Θ(t)均為隨機(jī)過程;A(t)和Θ(t)是慢變過程,Φ(t)是快變過程。在小時(shí)滯假設(shè)下將 Xτi、X˙τi近似展開為如下形式:

        利用式(4)、(5)可將 Xτi和 X˙τi轉(zhuǎn)化為非時(shí)滯情形,故而系統(tǒng)(1)等價(jià)為如下非時(shí)滯系統(tǒng):

        由等價(jià)線性化原理,將式(6)變形為如下形式:

        由式(3)可得:

        由式(3)、(7)、(9)可得:

        對(duì)式(10)應(yīng)用隨機(jī)平均原理[16],得到A(t)的It?隨機(jī)微分方程為:

        其中,W(t)為單位維納過程,式(11)對(duì)應(yīng)的FPK方程為:

        初值問題式(13)、(14)的理論精確解難以得到,下面通過Galerkin算法進(jìn)行近似求解。

        3.2 Galerkin算法原理

        首先考慮式(1)的退化線性方程,如下:

        在式(13)中令ξe(a)=0得與式(15)對(duì)應(yīng)的FPK方程:

        利用變量分離法[17],得式(16)的特征值與特征函數(shù):

        其中,λn表示特征值,An表示正交基函數(shù),σs表示系統(tǒng)(15)中X(t)的平穩(wěn)方差,Ln(˙)表示n階Laguerre多項(xiàng)式。

        方程(13)的解可近似寫成如下形式:

        其中,cn(t)是待求解項(xiàng),p0(a,t)是式(16)的解。表示 p(a,t)相對(duì)于 p0(a,t)的偏移程度,此偏移由非線性項(xiàng)和時(shí)滯項(xiàng)引起。實(shí)際計(jì)算中,級(jí)數(shù)求和形式需進(jìn)行適當(dāng)截?cái)?n=0,1,…,N)。

        將式(20)代入式(13),并結(jié)合FPK方程特征函數(shù)定義,得如下誤差參量:

        以Aj(a)A0(a)為權(quán)函數(shù),基于Galerkin方法,令誤差參量滿足下列條件:

        根據(jù)Laguerre多項(xiàng)式的性質(zhì),推得如下化簡(jiǎn)關(guān)系式:

        將式(23)代入式(24)并按式(25)進(jìn)行化簡(jiǎn),得到關(guān)于cj(t)的方程組,如下:

        由式(14)和式(20)可得式(26)的初始條件為:cj(t)=0。將式(26)的解代入式(20)即得瞬態(tài)響應(yīng)概率密度。

        4 算例分析

        選如下復(fù)雜色模型為例進(jìn)行分析。此模型稱為白噪聲二階濾過過程[18],其數(shù)學(xué)表達(dá)為:其中,常數(shù)?k和 βk是噪聲參數(shù),影響噪聲帶寬,?k和 βk增大代表帶寬變寬。Wk(t)是強(qiáng)度為2Dk的高斯白噪聲,此過程有如下譜密度:

        算例中取非線性項(xiàng)和時(shí)滯項(xiàng)為如下形式:

        其中,ω0和ξ0與(1)中相同,a1、a2、β0和ε均為常數(shù)。算例中取一個(gè)色噪聲激勵(lì)的情況(h1=1)。

        4.1 理論分析

        將式(29)、(30)代入式(1)即得色噪聲激勵(lì)的雙時(shí)滯Duffing-Van Der Pol振子。

        將式(29)、(30)代入式(8)并結(jié)合式(4)、(5)可得:

        將式(32)代入式(13)得與式(31)相應(yīng)的FPK方程為:

        由式(17)~(19)得特征值和基底函數(shù)為:

        將式(32)代入式(23)得誤差參量為:

        由式(18)并結(jié)合Laguerre多項(xiàng)式性質(zhì)得化簡(jiǎn)關(guān)系式:

        由式(32)代入式(26)并利用式(39)、(40)化簡(jiǎn),最終得到關(guān)于cj的方程組如下:

        式(41)中,cj(t)下標(biāo)為負(fù)時(shí)認(rèn)為不存在;式(43)中δj,n為Kronecker Delta符號(hào)。

        4.2 數(shù)值分析

        本文4.1節(jié)已舉例對(duì)第3章建立的算法進(jìn)行了理論分析,為說明色噪聲參數(shù)和時(shí)滯參數(shù)對(duì) p(a,t)的影響,本節(jié)選不同參數(shù)值進(jìn)行數(shù)值計(jì)算,同時(shí)對(duì)式(31)進(jìn)行MCS,將理論解和MCS解對(duì)比來驗(yàn)證算法可行性。

        圖1給出了方程(41)的解(N=8)。系統(tǒng)參數(shù)取值為:ω0=1,ξ0=0.02,ε=5,β0=0.5,α1=0.01,α2=0,τ1=τ2=0.5;噪聲參數(shù)取值為 D1=1,?1=3,β1=0.5。從圖中可直觀看出,隨時(shí)間增加,cj(t)將趨向于一系列穩(wěn)定值,穩(wěn)定時(shí)間約為100 s。

        圖1 cj(t)的時(shí)間演變情況

        圖2顯示了不同時(shí)刻振幅概率密度函數(shù) p(a,t) (N=20),其參數(shù)取值與圖1中的參數(shù)相同。采用Euler算法進(jìn)行MCS??梢钥闯?,由本文建立的算法得到的理論解和MCS解吻合程度好,這證明了本文建立的算法是有效的。由于系統(tǒng)在100 s后達(dá)平穩(wěn)狀態(tài),故該方法在一定程度上刻畫了穩(wěn)態(tài)響應(yīng),為求解穩(wěn)態(tài)概率密度提供了新途徑。

        圖3為振幅概率密度 p(a,t)(N=20),?1=3.5,其余參數(shù)與圖2中的參數(shù)相同。可以看出,理論近似解和MCS解擬合程度非常好。比較圖2和圖3易知,?1對(duì) p(a,t)產(chǎn)生了顯著影響,隨著?1增大,p(a,t)有一個(gè)向左的偏移,在較小幅值處達(dá)到峰值。從概率論角度解釋為:系統(tǒng)將以更大的概率在更小的幅值運(yùn)動(dòng),增大?1可降低系統(tǒng)響應(yīng)程度。

        圖2 p(a,t)(ω0=1,ξ0=0.02,ε=5,β0=0.5,α1=0.01,α2=0,τ1=τ2=0.5,D1=1,=3,β1=0.5)

        圖3 p(a,t)(ω0=1,ξ0=0.02,ε=5,β0=0.5,α1=0.01,α2=0,τ1=τ2=0.5,D1=1,=3.5,β1=0.5)

        圖4為令 β1=1.0其余參數(shù)與圖2中相同時(shí),p(a,t)圖像。比較圖2和圖4發(fā)現(xiàn),參數(shù)β1對(duì) p(a,t)的影響是顯著的。增大 β1使 p(a,t)向左偏移,在較小幅值處達(dá)到峰值。從概率角度可得如下結(jié)論:增大β1可降低系統(tǒng)響應(yīng)程度。

        圖4 p(a,t)(ω0=1,ξ0=0.02,ε=5,β0=0.5,α1=0.01,α2=0,τ1=τ2=0.5,D1=1,=3,β1=1)

        圖5~圖7分析了時(shí)滯參數(shù)對(duì) p(a,t)的影響。參數(shù)取值為:ω0=1,ξ0=0.01,ε=2.5,β0=0.5,α1=0.01,α2=0,D1=1,?1=3,β1=0.5。圖5中τ1=τ2=0.5,圖6中τ1=τ2=1.5,圖7中τ1=τ2=4.7。

        圖5~圖7表明以下兩點(diǎn):(1)理論解和MCS解吻合相當(dāng)一致,證明了本文算法的有效性;(2)對(duì)比圖5~圖7發(fā)現(xiàn),時(shí)滯對(duì) p(a,t)產(chǎn)生了影響,影響較為顯著。

        圖5 p(a,t)(ω0=1,ξ0=0.01,ε=2.5,β0=0.5,α1=0.01,α2=0,D1=1,=3,β1=0.5,τ1=τ2=0.5)

        圖6 p(a,t)(ω0=1,ξ0=0.01,ε=2.5,β0=0.5,α1=0.01,α2=0,D1=1,=3,β1=0.5,τ1=τ2=1.5)

        圖7 p(a,t)(ω0=1,ξ0=0.01,ε=2.5,β0=0.5,α1=0.01,α2=0,D1=1,=3,β1=0.5,τ1=τ2=4.7)

        綜合圖2~圖6分析,色噪聲和時(shí)滯均對(duì)系統(tǒng)瞬態(tài)概率密度 p(a,t)產(chǎn)生了影響,故在研究非線性系統(tǒng)瞬態(tài)概率密度時(shí)考慮色噪聲和時(shí)滯的聯(lián)合作用是非常必要的。數(shù)值計(jì)算過程中發(fā)現(xiàn),本文算法不但可以有效求解此類方程,且求解效率相對(duì)MCS是極高的。以N=15為例,利用本文算法所需計(jì)算時(shí)間與MCS相比縮減1/12;當(dāng)N=25時(shí),本文算法耗時(shí)為MCS的1/3。

        5 結(jié)論

        建立了色噪聲激勵(lì)的時(shí)滯非線性系統(tǒng)模型,并對(duì)此模型建立了求解瞬態(tài)響應(yīng)的相關(guān)算法。該算法包括三方面:

        (1)將時(shí)滯方程簡(jiǎn)化為非時(shí)滯方程;

        (2)利用等價(jià)線性化和隨機(jī)平均法得到系統(tǒng)振幅的平均It?隨機(jī)微分方程和平均FPK方程;

        (3)由FPK方程本征函數(shù)法得到正交基,進(jìn)行Galerkin變分并求解相應(yīng)的常微分方程組得到與時(shí)間相關(guān)的系數(shù),進(jìn)而得到瞬態(tài)概率密度函數(shù)。

        以白噪聲二階過濾過程為色噪聲模型,以Duffing-Van Der Pol振子為算例實(shí)現(xiàn)本文建立的算法。利用MCS對(duì)原系統(tǒng)進(jìn)行模擬,模擬解和理論解對(duì)比發(fā)現(xiàn)擬合程度非常好,證明本文的近似算法是有效的。通過對(duì)色噪聲參數(shù)和時(shí)滯參數(shù)的研究發(fā)現(xiàn)色噪聲和時(shí)滯對(duì)瞬態(tài)概率密度具有顯著影響,研究瞬態(tài)概率密度時(shí)考慮色噪聲和時(shí)滯的作用是非常必要的。

        本文算法可通過觀察與時(shí)間相關(guān)的系數(shù)的解,直觀看出系統(tǒng)是否具有穩(wěn)態(tài)解,并可估計(jì)系統(tǒng)到達(dá)穩(wěn)態(tài)的時(shí)間,使得利用MCS計(jì)算穩(wěn)態(tài)響應(yīng)時(shí)能大大減少盲目運(yùn)算時(shí)間,克服了使用MCS計(jì)算穩(wěn)態(tài)響應(yīng)時(shí)耗費(fèi)時(shí)間的問題。該算法是研究隨機(jī)非線性振動(dòng)系統(tǒng)瞬態(tài)概率密度的一種簡(jiǎn)便方法,對(duì)理論研究及工程應(yīng)用具有重要意義。

        [1]Mohamed W,Hamza A B.Medical image registration using stochastic optimization[J].Opt and Lasers Eng,2010,48:1213-1223.

        [2]Liu J,Liu X Z,Xie W C,et al.Stochastic consensus seeking with communication delays[J].Automatic,2011,47:2689-2696.

        [3]Repperger D W,F(xiàn)arris K A.Stochastic resonance—a nonlinear control theory interpretation[J].Int J Syst Sci,2010,41:897-907.

        [4]Lin Y K,Cai G Q.Probabilistic structural dynamics[M].New York:Mc Graw Hill,2004.

        [5]Zhu W Q,Deng M L.Optimal bounded control for minimizing the response of quasi-integraable Hamiltonian systems[J].Non-Lin Mech,2004,39:1535-1546.

        [6]Guo P R,Xu W,Liu D.Upper bound for the time derivative of entropy for a stochastic dynamical system with double singularitied driven by non-Gaussian noise[J].Chinese Phys B,2010,19:1-6.

        [7]朱位秋.隨機(jī)振動(dòng)[M].北京:科學(xué)出版社,1992:223-308.

        [8]Bhandari R G,Sherrer R E.Random vibrations in discrete nonlinear dynamic systems[J].J Mech Eng Sci,1968,10:168-174.

        [9]Wen Y K.Approximate methods for nonlinear random vibration[J].J Eng Mech Div,1975,101:389-401.

        [10]Spanos P D,Sofi A,Di Paola M.Nonstationary response envelope probability densities of nonlinear oscillators[J].J of Appl Mech,2007,74:315-324.

        [11]Hetel L,Daafouz J,Jungers M.Delay-dependent sampled-data control based on delay estimates[J].Systemsand Control Letters,2011,60:146-150.

        [12]Maccari A.Time delay control for two van der pol oscillators[J].J Comput Nonlin Dyn,2011,6:11-16.

        [13]Luan B L,Keiji K,Naoyuki H.Design and experimental verification of multiple delay feedback control for time-delay nonlinear oscillators[J].Nonlinear Dynam,2012,67:1407-1418.

        [14]Jin X L,Huang Z L.Nonstationary probability densities of strongly nonlinear single-degree-of-freedom oscillators with time delay[J].Nonlinear Dynam,2010,59:195-206.

        [15]Qi L Y,Wu W,Gu X D.Nonstationary probability densities of a class of nonlinear system excited by external colored noise[J].Sci China-Phys Mech Astron,2012,55:477-482.

        [16]Stratonovich R L.Topics in the theory of random noise,Vols 2[M].New York:Gordon and Breach,1986:170-305.

        [17]Gardiner C W.Handbook of stochastic methods for physics,chemistry and the natural science[M].Berlin:Springer-Velarge,1983:154-170.

        [18]Feng C S,Wu Y J,Zhu W Q.Response of duffing system with delayed feedback control under combined harmonic and real noise excitations[J].Commun Nonlinear Sci Numer Simulat,2009,14:2542-2550.

        1.Department of Applied Mathematics,Northwestern Polytechnical University,Xi'an 710129,China

        2.School of Electronic Information,Northwestern Polytechnical University,Xi'an 710129,China

        An effective approach to calculate the transient probability densities of the multi-delayed nonlinear system driven by colored noise excitations is developed.The system with time delay is simplified to an equivalent system without time delay.The linearization technique and the stochastic averaging method are adopted to obtain the averaged It? stochastic differential equation and the corresponding FPK(Fokker-Planck-Kolmogorov)equation for the amplitude process.A set of orthogonal base functions is obtained from the degenerated linear system.The Galerkin method is applied in the orthogonal base space to obtain the approximate probability densities.The proposed procedure is applied to the Duffing-Van Der Pol oscillator with two time delays and an external colored noise.The reliability of the theoretical results is verified by MCS(Monte Carlo simulation).Effects of the colored noise and the time delay are also discussed.The results show that the proposed method is effective at studying the transient probability densities of the time-delayed nonlinear system driven by colored noises;the theoretical calculation efficiency is higher than that of MCS;both of the colored noises and time delay affect the transient responses.

        time delay;colored noise;FPK(Fokker-Planck-Kolmogorov)equation;transient probability density;Galerkin method

        A

        TP301.5

        10.3778/j.issn.1002-8331.1210-0061

        QI Luyuan,XU Wei,GAO Weiting.Study on transient probability densities of delayed nonlinear system excited by colored noises.Computer Engineering and Applications,2013,49(7):1-5.

        國(guó)家自然科學(xué)基金(No.11172233,No.10932009,No.61171155);陜西省自然科學(xué)基金(No.2012JM8010);西北工業(yè)大學(xué)博士論文創(chuàng)新基金(No.CX201215)。

        戚魯媛(1986—),女,博士研究生,研究領(lǐng)域?yàn)閼?yīng)用數(shù)學(xué),非線性隨機(jī)振動(dòng);徐偉(1957—),男,博士,教授,博士生導(dǎo)師,研究方向?yàn)閼?yīng)用數(shù)學(xué);高維廷(1984—),男,博士研究生,研究方向?yàn)橥ㄐ判盘?hào)處理。E-mail:qiluyuan@gmail.com

        2012-10-09

        2012-12-18

        1002-8331(2013)07-0001-05

        CNKI出版日期:2012-12-26 http://www.cnki.net/kcms/detail/11.2127.TP.20121226.1120.001.html

        猜你喜歡
        概率密度將式時(shí)滯
        AKNS方程的三線性型及周期孤立波解
        連續(xù)型隨機(jī)變量函數(shù)的概率密度公式
        帶有時(shí)滯項(xiàng)的復(fù)Ginzburg-Landau方程的拉回吸引子
        因子von Neumann代數(shù)上非線性*-Lie導(dǎo)子的刻畫
        單自由度系統(tǒng)
        Hunt過程在Girsanov變換下的轉(zhuǎn)移概率密度的表示公式
        阻尼系統(tǒng)的特征
        隨機(jī)變量線性組合的分布的一個(gè)算法
        隨機(jī)結(jié)構(gòu)-TMD優(yōu)化設(shè)計(jì)與概率密度演化研究
        一階非線性時(shí)滯微分方程正周期解的存在性
        丰满人妻被猛烈进入中文字幕护士| 人妻丰满av无码中文字幕| 四虎国产精品永久在线| 大肉大捧一进一出视频出来呀| 色综合久久久久久久久五月| 亚洲中文字幕日产喷水| 国成成人av一区二区三区| 国产在线欧美日韩一区二区| 91免费国产高清在线| 国产男女做爰猛烈视频网站| 中文字幕有码手机视频| 日韩精品中文一区二区三区在线 | 色偷偷av男人的天堂| 暖暖视频在线观看免费| 亚洲Av午夜精品a区| 日本福利视频免费久久久| 青青草视频是针对华人| 精品综合久久久久久888蜜芽| 夜夜未满十八勿进的爽爽影院| 色综合久久久久综合999| 久久久99精品国产片| 国产精品女同一区二区免| 在线一区二区三区国产精品| 少妇高潮流白浆在线观看| 亚洲男人第一av网站| 国产目拍亚洲精品区一区| 久久综合精品国产丝袜长腿| 一本色道久久爱88av| 亚洲国产综合精品 在线 一区| 26uuu欧美日本在线播放| 国产二区中文字幕在线观看| 亚洲综合天堂av网站在线观看| 久久综合九色综合欧美狠狠| 欧美日韩国产在线成人网| 亚洲av无吗国产精品| 最美女人体内射精一区二区| 日韩高清在线观看永久| 国产一起色一起爱| 国产无卡视频在线观看| 少妇人妻中文字幕hd| 99久久国产综合精品麻豆|