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

        ?

        改進(jìn)CCM算法檢測(cè)外部擾動(dòng)下系統(tǒng)變量間的時(shí)滯和因果關(guān)系

        2016-12-14 09:00:12羅磊程非凡邱彤趙勁松
        化工學(xué)報(bào) 2016年12期
        關(guān)鍵詞:流形時(shí)滯擾動(dòng)

        羅磊,程非凡,邱彤,趙勁松

        ?

        改進(jìn)CCM算法檢測(cè)外部擾動(dòng)下系統(tǒng)變量間的時(shí)滯和因果關(guān)系

        羅磊,程非凡,邱彤,趙勁松

        (清華大學(xué)化學(xué)工程系,北京 100084)

        在化工過(guò)程中,可以通過(guò)分析變量間的相互作用和時(shí)滯關(guān)系,推理故障的傳播路徑和網(wǎng)絡(luò),指出故障的根原因。這對(duì)提高過(guò)程安全性,增強(qiáng)經(jīng)濟(jì)效益具有重要意義,是研究的熱點(diǎn)和難點(diǎn)。常見(jiàn)的互相關(guān)函數(shù)和傳遞熵等方法,由于只適用于線性或弱非線性系統(tǒng),或計(jì)算量較大,往往無(wú)法準(zhǔn)確地獲得變量間的時(shí)滯信息和作用強(qiáng)度,在實(shí)際應(yīng)用中存在不足。近年來(lái),在生態(tài)領(lǐng)域研究中提出的交叉收斂映射(CCM)算法,被認(rèn)為是一種適用于非線性耦合過(guò)程因果分析的方法,可適用于耦合變量間時(shí)滯關(guān)系的檢測(cè)。但對(duì)于帶有外部擾動(dòng)的化工過(guò)程,CCM無(wú)法根據(jù)隨時(shí)受到擾動(dòng)的過(guò)程數(shù)據(jù)構(gòu)造出穩(wěn)定的嵌入流形,導(dǎo)致了時(shí)滯和因果分析失敗。而基于CCM進(jìn)行改進(jìn)的擾動(dòng)過(guò)濾交叉映射(DFCM)方法,通過(guò)分析外部擾動(dòng)對(duì)系統(tǒng)的影響,預(yù)先篩選擾動(dòng)數(shù)據(jù),再將過(guò)濾后的數(shù)據(jù)代入交叉映射的計(jì)算中。算例表明,這種處理能有效地避免擾動(dòng)下嵌入流形不穩(wěn)定的問(wèn)題,適用于處于擾動(dòng)下的化工過(guò)程,并得到良好的時(shí)滯和因果關(guān)系分析效果。

        CCM;因果分析;時(shí)滯檢測(cè);算法;過(guò)程系統(tǒng);安全

        引 言

        一直以來(lái)針對(duì)大型復(fù)雜系統(tǒng)的因果關(guān)系識(shí)別和時(shí)滯檢測(cè)都較為困難,尤其是當(dāng)系統(tǒng)具有非線性特征時(shí)。目前已有人從不同角度出發(fā)提出了若干方法用于解決過(guò)程系統(tǒng)中變量間的因果關(guān)系識(shí)別和時(shí)滯分析。最簡(jiǎn)單的是互相關(guān)函數(shù)法(CCF)[1],這一方法通過(guò)計(jì)算在不同的時(shí)滯假設(shè)點(diǎn)上的互相關(guān)函數(shù)值,選取其絕對(duì)值的最大值對(duì)應(yīng)的時(shí)滯作為時(shí)滯檢測(cè)的結(jié)果。盡管該方法計(jì)算簡(jiǎn)單,也適用于某些實(shí)際過(guò)程,但它對(duì)非線性系統(tǒng)適應(yīng)性較差。另一種被廣泛使用的格蘭杰因果檢驗(yàn),則往往因?yàn)闊o(wú)法滿足待檢測(cè)變量間相互獨(dú)立的條件,不能適用于耦合系統(tǒng)的分析。Schreiber等[2-10]基于信息熵的概念提出了傳遞熵(transfer entropy)的方法,這一方法也被廣泛地應(yīng)用于化工、大氣科學(xué)、神經(jīng)科學(xué)等領(lǐng)域。但由于其基于概率計(jì)算的理論,對(duì)概率密度函數(shù)分布的計(jì)算必須足夠精確以確保結(jié)果的準(zhǔn)確性,這就帶來(lái)了較大的計(jì)算量。另外,有時(shí)傳遞熵并不能準(zhǔn)確地計(jì)算出變量間的時(shí)滯關(guān)系,所以無(wú)法獲得類似于符號(hào)有向圖(signed directed graph)[11-17]的結(jié)果以清晰地總結(jié)和分析結(jié)論。總之,上述的方法或者基于概率理論,計(jì)算量較大;或者只能適用于線性或者是弱非線性系統(tǒng),難以應(yīng)用于對(duì)非線性系統(tǒng)的分析。

        近年來(lái)提出的交叉收斂映射算法(convergent cross mapping,CCM)[18],利用非線性耦合系統(tǒng)中變量間微分同胚映射的特點(diǎn),將原有的因果關(guān)系識(shí)別轉(zhuǎn)化為對(duì)嵌入流形相互預(yù)測(cè)效果的對(duì)比。Ye等[19]的工作表明CCM算法適用于某些自振系統(tǒng)中變量因果關(guān)系的檢測(cè)和時(shí)滯分析。但是,化工過(guò)程中處于常見(jiàn)的外部擾動(dòng)時(shí),Ye等的方法也許不再適用。本文對(duì)CCM算法進(jìn)行改進(jìn),提出DFCM方法,解決了該問(wèn)題。

        1 CCM算法

        CCM算法由Sugihara等[18]在2012年提出,相比于互相關(guān)函數(shù)和傳遞熵,它能更清楚地把握耦合系統(tǒng)中變量間的相互作用,得到更準(zhǔn)確的結(jié)果。根據(jù)Takens嵌入定理,在一動(dòng)態(tài)系統(tǒng)中,相互關(guān)聯(lián)兩變量其對(duì)應(yīng)的影子流形MM都與原系統(tǒng)的流形微分同胚。這里的影子流形MM都是通過(guò)對(duì)原始數(shù)據(jù)和進(jìn)行時(shí)序嵌入重構(gòu)得到,其元素分別用xy表示。這樣,流形MM之間便形成了連續(xù)單射的聯(lián)系。對(duì)于M上選定的點(diǎn)x,1,x,2,…, x,k,在M上也一定能找到對(duì)應(yīng)的y,1,y,2,…, y,k,若所有的x(=1,…,)均收斂于特定的x,0,那么相應(yīng)地,y(=1,…,)也會(huì)收斂于x,0在M上所對(duì)應(yīng)的y,0。對(duì)于M上每個(gè)點(diǎn)x,CCM都尋找其相鄰點(diǎn)及其在流形M上的對(duì)應(yīng)點(diǎn),并對(duì)y,i進(jìn)行預(yù)測(cè),記作?,i。這樣變量對(duì)的預(yù)測(cè)能力ρ→y即是預(yù)測(cè)值{?,i}和真實(shí)值{y,i}之間的相關(guān)系數(shù)。隨著樣本總數(shù)的增加,預(yù)測(cè)值序列{?,i}逐漸收斂,ρ→y也將會(huì)收斂到一定值,這樣便可使用變量預(yù)測(cè),反之亦然。然而,如果變量和之間不存在耦合關(guān)系,由于信息傳遞的不對(duì)稱將會(huì)導(dǎo)致作為結(jié)果的某一變量對(duì)另一作為原因的變量預(yù)測(cè)效果的下降。圖1解釋了這一原理。

        對(duì)于維動(dòng)態(tài)系統(tǒng)

        d/d=(,) (1)

        其中和分別對(duì)應(yīng)變量和參數(shù)向量,=[1,2,…,x]T,=[1,2,…,u]T。對(duì)于中的兩變量和,其對(duì)應(yīng)的影子流形構(gòu)造如下

        是嵌入維數(shù),是延遲時(shí)間,這樣根據(jù)前面所述,便有

        若在M上找到一系列最接近M,0的點(diǎn){M,}={M,1,M,2,…,M,+1},那么M,0的值便可由{M,}進(jìn)行估計(jì)

        如果樣本數(shù)足夠大以保證{M,}足夠接近M,0,那么最終M預(yù)測(cè)M的能力ρ→y將會(huì)收斂到0和1之間。

        Ye等[19]將CCM的應(yīng)用拓展到對(duì)變量間時(shí)滯關(guān)系的確定中。根據(jù)他的想法,如果變量對(duì)的作用存在時(shí)滯Δ,即()=((?Δ)),那么使用CCM將對(duì)進(jìn)行估計(jì)最強(qiáng)的作用應(yīng)出現(xiàn)在時(shí)滯為+Δ的位置,這里正號(hào)表示落后于。圖2對(duì)Ye等的思想進(jìn)行了簡(jiǎn)要的解釋。

        經(jīng)過(guò)計(jì)算發(fā)現(xiàn),Ye等的方法在對(duì)自振系統(tǒng)進(jìn)行計(jì)算時(shí)能取得良好的效果,但是,如果系統(tǒng)變量振蕩來(lái)源于外部輸入的擾動(dòng),Ye等的方法并沒(méi)有考慮到這部分變化對(duì)CCM計(jì)算結(jié)果的影響,計(jì)算結(jié)果并不理想。圖3列出了將Ye等的方法應(yīng)用于3個(gè)不同的非線性時(shí)滯系統(tǒng)的計(jì)算結(jié)果,例1是一自振系統(tǒng),例2和例3中都加入了參數(shù)的擾動(dòng),正是該擾動(dòng)導(dǎo)致本應(yīng)穩(wěn)定的系統(tǒng)出現(xiàn)了波動(dòng)。從最后的結(jié)果可以看出,Ye等的方法能準(zhǔn)確地計(jì)算出例1中變量之間的時(shí)滯,但卻在例2和例3中無(wú)法得到準(zhǔn)確的結(jié)果。

        blue line—maps; red line—maps

        圖3中具體算例如下。

        例1

        x+1=x(3.78?3.78x?0.07y3)

        y+1=y(3.77?3.77y?0.08x5)

        例2

        x+1=0.3x+(1+0.3u?y?25)x/5

        y+1=0.9y+0.02(x?10?5.5)

        u= rand(0,1),=1, 2, 3, …

        例3

        x+1=x+D[0.4(1+0.1u)?0.2x+0.3y?80]

        y+1=y+D[0.2x?50?0.3y?0.4(1+0.1u)y]

        u= rand(0,1),=1, 2, 3, …

        2 擾動(dòng)過(guò)濾交叉映射(DFCM)法

        盡管CCM最初被應(yīng)用于生態(tài)學(xué)領(lǐng)域,但其中微分同胚的特點(diǎn)也常見(jiàn)于化工過(guò)程數(shù)據(jù)的重構(gòu)流形中,因此CCM算法的思想也適用于化工過(guò)程,但算例2和例3中參數(shù)的頻繁改變使得MM不再穩(wěn)定,這導(dǎo)致CCM算法無(wú)法識(shí)別到正確的時(shí)滯并對(duì)變量間的相互作用作出計(jì)算。以一個(gè)二元系統(tǒng)為例,如圖4所示,如果參數(shù)被設(shè)定為不同值,該系統(tǒng)變量和將會(huì)收斂到不同的穩(wěn)態(tài)值。同時(shí),由于穩(wěn)態(tài)點(diǎn)的不斷變化,相軌跡也會(huì)出現(xiàn)相互交叉,這將會(huì)導(dǎo)致原CCM方法計(jì)算準(zhǔn)確度的下降。

        例如圖4中紅色和藍(lán)色軌跡的交點(diǎn)對(duì)應(yīng)于圖5中曲線的突變,這說(shuō)明,即使在原流形或者M中發(fā)現(xiàn)兩相鄰點(diǎn),由于相軌跡的交叉等原因也并不能保證它們?cè)?i>M中也滿足相鄰條件。還有一種可能是它們對(duì)應(yīng)于不同的參數(shù),只是恰好在圖中相鄰而已。這樣,CCM尋找到的相鄰點(diǎn)在其之后的計(jì)算中是沒(méi)有意義的。

        根據(jù)以上對(duì)CCM算法的分析,對(duì)其進(jìn)行改進(jìn),希望它能對(duì)處于外部擾動(dòng)下的系統(tǒng)進(jìn)行計(jì)算,從而提出擾動(dòng)過(guò)濾交叉映射(disturbance filtrated cross-mapping,DFCM)法。該方法基于以下幾點(diǎn)假設(shè):

        (1)系統(tǒng)對(duì)于每個(gè)有且僅有一個(gè)穩(wěn)態(tài),的擾動(dòng)不會(huì)導(dǎo)致多穩(wěn)態(tài)現(xiàn)象的發(fā)生;

        (2)穩(wěn)態(tài)點(diǎn)的變化是連續(xù)的;

        (3)如果參數(shù)和0接近,那么對(duì)應(yīng)的從相同初值出發(fā)的相軌跡也會(huì)收斂到0所對(duì)應(yīng)的相軌跡上。

        圖6對(duì)該過(guò)程進(jìn)行了闡述。圖中從左至右參數(shù)的變化逐漸增大,不同下對(duì)應(yīng)的影子流形與原流形的差別也逐漸增大。這樣,通過(guò)限定擾動(dòng)的變化便可獲得穩(wěn)定的流形。

        blue—original manifold; red—new manifold

        圖7是DFCM計(jì)算的大致流程。在DFCM中,首先選定合適的嵌入維數(shù)和時(shí)間延遲,將過(guò)程數(shù)據(jù)根據(jù)進(jìn)行過(guò)濾篩選,然后再將篩選后的數(shù)據(jù)代入交叉收斂映射的計(jì)算。

        3 算例

        3.1 DFCM對(duì)例2和例3的計(jì)算

        圖8顯示了將DFCM用于之前的例2和例3后的計(jì)算結(jié)果,對(duì)于例2,DFCM的計(jì)算表明最強(qiáng)的作用出現(xiàn)在變量比領(lǐng)先Δ=12的位置,這一結(jié)果與圖中給出的兩變量的輸出結(jié)果相符合,注意為了方便比較,兩輸出的橫坐標(biāo)已按照DFCM的結(jié)果進(jìn)行了相應(yīng)平移。對(duì)于例3,DFCM不僅在比領(lǐng)先Δ1=50和Δ2=-80發(fā)現(xiàn)了時(shí)滯點(diǎn),還發(fā)現(xiàn)兩變量的最強(qiáng)作用出現(xiàn)的位置是在Δ=0的位置,這與實(shí)際的輸出相符合(圖中兩變量呈反相關(guān))。

        blue line—maps; red line—maps

        3.2 DFCM對(duì)四罐耦合系統(tǒng)的計(jì)算

        圖9是一個(gè)四罐耦合模型[20],模型輸入為泵的輸入電壓v,其對(duì)應(yīng)的流量為kv,γ表示流量分流比。為了維持整個(gè)系統(tǒng)的穩(wěn)定,分別在罐1和罐3加入了液位控制器。

        在該模型中假設(shè)液體從上游裝置流往下游需要經(jīng)過(guò)一段不可忽略的時(shí)間,用藍(lán)色數(shù)字表示。從圖9中很容易看出罐3和罐4應(yīng)該比其他兩個(gè)罐更容易受到2的擾動(dòng)的影響,它們應(yīng)該處于擾動(dòng)傳播路徑的上游位置。模型的輸出由圖10所示。

        從圖11的計(jì)算結(jié)果可以看出,3始終領(lǐng)先于其他變量,說(shuō)明它最接近擾動(dòng)的根源。之后便是4、1和2。這一計(jì)算結(jié)果與實(shí)際輸出和從流程結(jié)構(gòu)出發(fā)的分析相符。

        blue line—maps; red line—maps

        4 結(jié) 論

        CCM算法是一種新提出的可適用于線性和非線性系統(tǒng)因果分析和時(shí)滯檢測(cè)算法,但本文發(fā)現(xiàn)其對(duì)某些受到外部擾動(dòng)的系統(tǒng)并不能取得好的計(jì)算效果,原因是在外部擾動(dòng)存在的情況下無(wú)法構(gòu)造原方法中的穩(wěn)定的流形。針對(duì)這一問(wèn)題本文對(duì)原方法進(jìn)行改進(jìn),提出了擾動(dòng)過(guò)濾交叉映射(DFCM)方法。算例證明,該方法適用于處于外部擾動(dòng)下的化工過(guò)程中耦合變量之間時(shí)滯關(guān)系的檢測(cè)。未來(lái)將會(huì)進(jìn)一步完善相關(guān)工作,期待這一方法能在化工和其他領(lǐng)域的因果分析和時(shí)滯檢測(cè)中有進(jìn)一步的應(yīng)用。

        [1] BAUER M, THORNHILL N F. A practical method for identifying the propagation path of plant-wide disturbances [J]. J. Process Control, 2008, 18: 707-719.

        [2] BAUER M, COX J W, CAVENESS M H,. Finding the direction of disturbance propagation in a chemical process using transfer entropy [J]. IEEE Trans. Control Syst. Technol., 2007, 15: 12-21.

        [3] VICENTE R, WIBRAL M, LINDNER M,. Transfer entropy — a model-free measure of effective connectivity for the neurosciences [J]. J. Comput. Neurosci., 2011, 30: 45-67.

        [4] 謝中凱, 劉國(guó)華, 吳志根. 基于傳遞熵的梁結(jié)構(gòu)損傷動(dòng)力識(shí)別[J]. 浙江大學(xué)學(xué)報(bào)(工學(xué)版), 2012, 46(10): 1880-1886. XIE Z K, LIU G H, WU Z G. Dynamic damage identification for beam structures based on transfer entropy [J]. J. Zhejiang Univ. (Engineering Sci.), 2012, 46(10): 1880-1886.

        [5] KUEHNERT C, BEYERER J. Data-driven methods for the detection of causal structures in process technology [J]. Machines, 2014, 2: 255-274.

        [6] SHU Y, ZHAO J. Data-driven causal inference based on a modified transfer entropy [J]. Comput. Chem. Eng., 2013, 57: 173-180.

        [7] 張志森, 龔志強(qiáng), 支蓉. 利用傳遞熵對(duì)Lorenz系統(tǒng)和Walker環(huán)流信息傳輸方向的分析[J]. 物理學(xué)報(bào), 2013, 62(12): 548-557.ZHANG Z S, GONG Z Q, ZHI R. Analysis of the direction of information transfer of Lorenz system and Walker circulation with transfer entropy [J]. Acta Physica Sinica, 2013, 62(12): 548-557.

        [8] 賀丁, 趙勁松. 基于Hopfield網(wǎng)絡(luò)的時(shí)滯分析故障診斷策略[J]. 化工學(xué)報(bào), 2013, 64(2): 633-640. HE D, ZHAO J S. Fault strategy of time delay analysis based on Hopfield network [J]. CIESC Journal, 2013, 64(2): 633-640.

        [9] 賀丁, 舒逸聃, 趙勁松. 基于改進(jìn)的時(shí)滯分析算法的化工過(guò)程故障定位[J]. 化工學(xué)報(bào), 2012, 63(10): 3165-3172. HE D, SHU Y D, ZHAO J S. An improved time-delay algorithm for chemical processes fault identification [J]. CIESC Journal, 2012, 63(10): 3165-3172.

        [10] DUAN P, YANG F, SHAH L S,. Transfer zero-entropy and its application for capturing cause and effect relationship between variables [J]. IEEE Transactions on Control Systems Technology, 2015, 23(3): 855-867.

        [11] YANG F, XIAO D. Signed directed graph modeling of industrial processes and their validation by data-based methods[C]//2010 Conference on Control and Fault Tolerant Systems. Nice, France: IEEE, 2010: 387-392.

        [12] 王杭州, 陳丙珍, 何小榮, 等. 基于開(kāi)源組件的SDG推理平臺(tái)[J]. 化工學(xué)報(bào), 2010, 61(7): 1829-1836.WANG H Z, CHEN B Z, HE X R,. Open source signed digraph inference framework [J]. CIESC Journal, 2010, 61(7): 1829-1836.

        [13] 盧秉南, 張貝克, 馬昕, 等. 基于SDG模型的控制系統(tǒng)故障診斷方法[J]. 化工學(xué)報(bào), 2009, 60(9): 2243-2251. LU B N, ZHANG B K, MA X,. Fault diagnosis of control system based on signed directed graph model [J]. CIESC Journal, 2009, 60(9): 2243-2251.

        [14] YANG F, XIAO D Y. Progress in root cause and fault propagation analysis and large-scale industrial processes [J]. Journal of Control Science and Engineering, 2012: 478373.

        [15] DUAN P, CHEN T W, SHAH L S,. Methods for root cause diagnosis of plant-wide oscillations [J]. AIChE Journal, 2014, 60(6): 2019-2034.

        [16] HE B, CHEN T, YANG X H. Root cause analysis in multivariate statistical process monitoring: integrating reconstruction-based multivariate contribution analysis with fuzzy-signed directed graphs [J]. ComputersChemical Engineering, 2014, 64(7): 167-177.

        [17] WAN Y M, YANG F, Lü N,Statistical root cause analysis of novel faults based on digraph models [J]. Chemical Engineering Research and Design, 2013, 91: 87-99.

        [18] SUGIHARA G, MAY R, YE H,. Detecting causality in complex ecosystems [J]. Science, 2012, 338: 496-500.

        [19] YE H, DEYLE E R, GILATTANZ L J,. Distinguishing time-delayed causal interactions using convergent cross mapping [J]. Sci. Rep., 2015, 5:14750.

        [20] SADATI N, RAHMANI M, SAIF M,. Two-level robust optimal control of large scale nolinear systems [J]. IEEE Systems Journal, 2015, 9: 242-251.

        An improved convergent cross mapping algorithm for causality identification and time delay analysis between systemic variables under external disturbance

        LUO Lei, CHENG Feifan, QIU Tong, ZHAO Jinsong

        (Department of Chemical Engineering, Tsinghua University, Beijing 100084, China)

        In chemical processes, fault propagation pathways and root cause identification could be discovered though analysis of interactions and time delay relationships among different process variables. Because of its importance in improving process safety and operation profit, fault discovery has been a popular and challenging research topic. Common methods such as correlation and entropy transfer functions, which usually cannot get accurate time delay and interaction strength between variables by limited applicability for linear and weak nonlinear systems or high computation demand, have experienced many disadvantages in actual application. Recently, a new convergent cross mapping (CCM) algorithm in ecology has been considered suitable for causality analysis and time delay identification for nonlinear coupling process variables. However, CCM fails to find application in externally disturbed chemical processes because it cannot establish stable embedded flow from process data. An improved CCM, disturbance filtered cross mapping method (DFCM), overcame many challenges of creating stable embedded flow by analyzing external disturbance, filtering disturbed process data, and applying filtrated data to CCM calculation. Case studies showed good results of time delays and causality analysis, thus DFCM could be applied to chemical processes under external disturbance.

        CCM; causality analysis; time-delay detection; algorithm; process system; safety

        date: 2016-09-18.

        Prof. QIU Tong, qiutong@tsinghua.edu.cn

        10.11949/j.issn.0438-1157.20161300

        TP 277

        A

        0438—1157(2016)12—5122—09

        國(guó)家自然科學(xué)基金項(xiàng)目(U1462206)。

        supported by the National Natural Science Foundation of China (U1462206).

        2016-09-18收到初稿,2016-09-25收到修改稿。

        聯(lián)系人:邱彤。第一作者:羅磊(1992—),男,碩士研究生。

        猜你喜歡
        流形時(shí)滯擾動(dòng)
        Bernoulli泛函上典則酉對(duì)合的擾動(dòng)
        帶有時(shí)滯項(xiàng)的復(fù)Ginzburg-Landau方程的拉回吸引子
        緊流形上的Schr?dinger算子的譜間隙估計(jì)
        (h)性質(zhì)及其擾動(dòng)
        迷向表示分為6個(gè)不可約直和的旗流形上不變愛(ài)因斯坦度量
        Nearly Kaehler流形S3×S3上的切觸拉格朗日子流形
        小噪聲擾動(dòng)的二維擴(kuò)散的極大似然估計(jì)
        用于光伏MPPT中的模糊控制占空比擾動(dòng)法
        基于多故障流形的旋轉(zhuǎn)機(jī)械故障診斷
        一階非線性時(shí)滯微分方程正周期解的存在性
        国产内射视频在线免费观看| 最新国产精品亚洲二区| 国内精品久久久久国产盗摄| 日本大胆人体亚裔一区二区| 亚洲丰满熟女乱一区二区三区| 亚洲中国精品精华液| 激情综合色综合久久综合| 成全视频高清免费| 久久精品综合国产二区| 国产av精选一区二区| 无码人妻久久一区二区三区免费丨| 色拍自拍亚洲综合图区| 国产精品后入内射日本在线观看 | 日产无人区一线二线三线乱码蘑菇 | 亚洲第一成人网站| 国产一区二区三区视频大全| 日本一区二区三区丰满熟女| 天天躁夜夜躁狠狠是什么心态| 国产成人无码a区在线观看视频| 亚洲午夜无码久久yy6080| 在线视频一区二区三区中文字幕| 少妇又色又爽又高潮在线看| 人妻少妇出轨中文字幕| www插插插无码免费视频网站| 一区二区三区国产高潮| av毛片亚洲高清一区二区| 日本精品久久久久中文字幕| 亚洲精品乱码久久久久久蜜桃图片 | 久久777国产线看观看精品| 欧洲熟妇色xxxx欧美老妇多毛| 久久精品国产99精品国偷| 成年人男女啪啪网站视频| 成年人干逼视频水好多| 女女女女女裸体处开bbb| 九九热在线视频观看这里只有精品| 国产男女猛烈无遮挡免费视频网址| 丝袜美腿丝袜美腿丝袜美腿丝袜| 无码人妻一区二区三区兔费| 亚洲男人天堂| 按摩师玩弄少妇到高潮hd| 亚洲精品视频中文字幕|