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

        ?

        基于狀態(tài)集員估計的主動故障檢測

        2021-06-20 10:11:04王晶史雨茹周萌
        自動化學報 2021年5期
        關鍵詞:故障信號檢測

        王晶 史雨茹 周萌

        系統(tǒng)的操作安全以及產(chǎn)品的質量對于工業(yè)過程控制至關重要,而一旦系統(tǒng)發(fā)生故障不僅會影響產(chǎn)品質量,甚至會給整個系統(tǒng)帶來安全隱患,因此故障檢測及診斷具有重要的研究意義與研究價值[1].通常,系統(tǒng)根據(jù)是否發(fā)生故障及故障的大小可分為三個階段:第1 階段為正常階段,此時系統(tǒng)雖然受到未知干擾和測量噪聲的影響,但是被監(jiān)控系統(tǒng)的變量處在可接受的范圍內.當系統(tǒng)的被監(jiān)控變量的幅值或者特征發(fā)生微小變化,使得系統(tǒng)的某些量超出可接受的閾值,但是仍小于某一較大的閾值時,此時為第2 階段.這些微小的變化一般很難被觀測到,此階段為微小故障或者潛伏故障階段,該階段的變化不會大程度地影響系統(tǒng).但是微小故障會隨著時間的變化逐漸變大,當超過較大的這個閾值時,此時系統(tǒng)處在第3 階段,即顯著故障階段.該階段被監(jiān)控變量的幅值發(fā)生大幅度變化,很容易被監(jiān)控人員發(fā)現(xiàn)[2].系統(tǒng)處在第3 階段時,由于系統(tǒng)的故障特征十分明顯,往往會存在較為重大的安全隱患.因此,針對此階段的故障檢測與診斷引起了高度重視.在此階段通過系統(tǒng)的測量數(shù)據(jù)生成殘差,而無需設計輔助輸入信號激勵系統(tǒng)來判斷故障發(fā)生與否,此類故障診斷方法通常稱為被動故障診斷.被動故障診斷主要分為三類:基于模型的方法[3-5]、基于知識的方法[6-8] 以及基于數(shù)據(jù)驅動[9-11]的方法.

        當系統(tǒng)處在第2 階段,即系統(tǒng)發(fā)生微小故障時,由于其對系統(tǒng)的影響較小,往往未造成實質性傷害,此時,對微小故障進行及時有效地監(jiān)控,將對預防和減少災難的發(fā)生具有重要作用[12].但同時由于微小故障發(fā)生在初始階段,幅值較小,受未知干擾信號的影響,采用被動故障檢測及診斷方法往往不能準確地檢測出來.相較于被動故障診斷,主動故障診斷方法通過生成輔助輸入信號對原系統(tǒng)加以激勵,從而使得微小故障的影響得以顯現(xiàn),近年來受到越來越多研究者的關注[13-14].主動故障診斷方法根據(jù)未知干擾信號的處理方式可以分為確定性方法和概率性方法[15].確定性方法是將系統(tǒng)的干擾和噪聲等不確定的因素假設為具有已知上下界或者已知內部結構分布的信號,如控制器和檢測器參數(shù)設計[16-17]、設計輔助輸入信號[18]等方法.概率性方法是將系統(tǒng)不確定性假設為具有已知概率密度分布的變量,如結合貝葉斯方法的主動故障診斷[19]等方法.

        由于基于輔助輸入信號設計的主動故障診斷原理較為簡單,且不會改變系統(tǒng)的結構,所以受到研究者關注.傳統(tǒng)基于輔助輸入信號設計的主動故障診斷方法一般假設系統(tǒng)的干擾和測量噪聲服從某種概率分布,而在實際工業(yè)中的干擾和噪聲分布的先驗信息往往很難獲得.為此,研究者假設干擾和噪聲內部結構未知但是上下界已知,提出了集員估計理論[20].由于該假設符合實際工業(yè)過程干擾和噪聲的特性,得到研究者的青睞.集員估計理論期望得到與系統(tǒng)模型、測量輸出以及未知干擾相包容的狀態(tài)或參數(shù)的可行集,通常利用幾何體,如橢球、平行多面體、區(qū)間、全對稱多胞形(Zonotope)等來近似該可行集[21].文獻[22]針對于多模型系統(tǒng),利用不同模型輸出所在橢球交集為空,設計最小輔助輸入信號,對不同模型類型進行檢測分離.文獻[23]針對線性時不變系統(tǒng),利用0~N 時刻輸出所在的橢球,實現(xiàn)了離線故障檢測及輔助輸入信號設計.相比較于橢球,全對稱多胞形的結構以及運算方式簡單,保守性較小等優(yōu)點,近幾年得到大量研究.針對線性加性故障系統(tǒng),文獻[24]通過引入全對稱多胞形的思想,設計并比較了開環(huán)輸入信號和閉環(huán)輸入信號的大小.文獻[25]利用0~N 時刻正常和故障模型輸出所在的全對稱多胞形的交集為空,進行故障檢測,同時得到相應輔助輸入信號激勵系統(tǒng).

        主動故障檢測主要是利用設計得到的輔助輸入信號激勵系統(tǒng),實現(xiàn)微小故障檢測.因此在保證該輔助輸入信號能夠激勵故障的前提下,同時也應該對系統(tǒng)的影響較小.目前,基于集員估計的主動故障檢測的研究主要集中在通過分析系統(tǒng)的輸出所在的集合,來進行故障檢測和輔助輸入信號設計[26].但是由于系統(tǒng)輸出的測量值所在的集合中不僅包含了系統(tǒng)的過程干擾也包含了測量噪聲等不確定因素,其保守性較大.而系統(tǒng)的狀態(tài)所在的集合比輸出所在的集合少了下一時刻的測量噪聲,因此利用狀態(tài)所在集合交集為空設計輔助輸入信號可降低一部分保守性,使得設計得到的輔助輸入信號較小.為此,本文主要針對微小故障,研究基于狀態(tài)集員估計的輔助輸入信號設計方法.首先,利用全對稱多胞形卡爾曼濾波對系統(tǒng)狀態(tài)進行估計.然后利用該觀測器觀測得到系統(tǒng)狀態(tài)所在的集合,在保證輔助輸入信號最小的情況下,通過使正常和故障模型的狀態(tài)所在的全對稱多胞形的交集為空,得到優(yōu)化問題.將該優(yōu)化問題轉化為混合整數(shù)二次規(guī)劃(Mixed integer quadratic programming,MIQP),通過求解該MIQP 問題,得到最優(yōu)的輔助輸入信號.

        1 問題描述及相關概念

        1.1 問題描述

        考慮如下線性離散不確定系統(tǒng)

        其中,i ∈{0,1,···,q},xk∈Rnx,uk∈Rnu,yk∈Rny分別為系統(tǒng)在k 時刻的狀態(tài),輸入及測量輸出向量.q 代表發(fā)生故障的個數(shù).當i=0 時,代表系統(tǒng)正常,當i 為其他值時,則代表系統(tǒng)發(fā)生故障.A[i],B[i],C[i],E[i],F[i]分別為適當維度的參數(shù)矩陣.wwwk∈Rnw和Wk∈Rnv分別為系統(tǒng)的過程干擾及測量噪聲.不失一般性,本文假設系統(tǒng)的初始狀態(tài),過程干擾和測量噪聲分別滿足:x0∈X0,wwwk∈W,Wk∈V,其中,X0,W,V 分別為初始狀態(tài)、過程干擾以及測量噪聲所在的全對稱多胞形集合.

        通常,微小故障很難被觀測人員察覺,主動故障檢測利用系統(tǒng)的輸出數(shù)據(jù)產(chǎn)生決策及輔助輸入信號,將該輔助輸入信號注入給系統(tǒng),提高檢測質量.集員估計理論是主動故障檢測常用的一種方法,將系統(tǒng)的干擾及測量噪聲等不確定的因素用集合的形式表示成一種幾何形狀,如區(qū)間、橢球、多面體、全對稱多胞形等,然后設計最優(yōu)的輔助輸入信號,使得正常模型和故障模型的輸出所在的集合交集為空,即

        其中,Y[0]代表正常模型輸出所在的集合,Y[i]為故障模型輸出所在的集合.

        主動故障檢測的動態(tài)過程如圖1 所示.其中,圖1(a)是系統(tǒng)的輔助輸入信號為零時輸出的集合,當故障為微小故障時兩個集合相交,當實際輸出屬于相交部分時無法確認系統(tǒng)是正常還是故障.基于輸出集合的主動故障檢測的原理是通過設計輔助輸入信號,使得正常輸出集合與故障模型輸出集合相分離,如圖1(b)所示.

        圖1 在正常模型和故障模型下,加入輔助信號與未加輔助信號結果的對比Fig.1 Comparison of the results of adding auxiliary signals and unassisted signals under normal and fault models

        注意到已有的方法采用的輸出集合交集為空的方法,由式(1)可得輸出集合的表示

        其中,Y 代表輸出所在的集合,X 為狀態(tài)所在的集合,V 為測量噪聲所在的集合.

        由于CX 為線性變化,并不影響兩個交集的相交與否,如圖2(b)所示.但由于FV 的引入,使得原本相切的兩個集合相交,相交面積的大小與引入的測量噪聲V 的大小有關.若針對圖2(c)中利用輸出集合設計輔助輸入信號使得輸出集合交集為空,則引入的信號必將比利用狀態(tài)集合的交集為空的方法大,進而對原系統(tǒng)的影響大.

        圖2 在正常模型和故障模型下,加入輔助信號對狀態(tài)和輸出的影響Fig.2 The effect of adding an auxiliary signal on state and output under normal and fault models

        綜上分析,利用狀態(tài)交集為空設計輔助輸入信號的保守性較低.因此,本文主要針對系統(tǒng)(1)的狀態(tài)所在的集合交集為空設計最優(yōu)輔助輸入信號,用該最優(yōu)輔助輸入信號激勵系統(tǒng),提高故障檢測質量.

        1.2 基本概念

        目前主動故障檢測方法中常用的集合有橢球、全對稱多胞形等.相較于其他集合,全對稱多胞形的形式簡單、計算靈活性高,可以描述任意凸面體,并能夠適當?shù)目刂茋猍27]等.因此,本文主要研究利用全對稱多胞形對狀態(tài)集合進行描述,相關定義與基本性質描述如下.

        定義1[25].r 階全對稱多胞形Z 的表達式定義為

        其中,-1 ≤αi≤1,超立方體Br=[-1,+1]r,c ∈Rn為Z 的中心,向量g1,g2,···,gr為Z 的生成器,G={g1,g2,···,gr}∈Rn×r為Z 的生成器矩陣.

        全對稱多胞形滿足如下基本計算規(guī)則:

        Frobenius 半徑可以作為衡量全對稱多胞形大小的指標,一般由生成矩陣G 的Frobenius 范數(shù)表示:

        由式(5)可知,隨著時間的推移,全對稱多胞形的階數(shù)會不斷增加,最終無法計算,為此,需要對全對稱多胞形的生成矩陣的維數(shù)加以限制.

        引理1[28].全對稱多胞形可以由一個最小的盒子包含

        其中,□Z 代表包含集合Z 的最小盒子,rs(G)為對角矩陣,滿足

        引理2[29].考慮r 維的全對稱多胞形Z=c ⊕GBr?Rn以及整數(shù)s,且滿足n ≤s ≤r,定義ˉG為將G 的列向量按歐氏范數(shù)降序排列所得矩陣,則全對稱多胞形Z 可以包含一個s 維的全對稱多胞形,即

        引理3[25].已知全對稱多胞形Z=〈az+bz,Gz〉和Y=〈ay+by,Gy〉,則Z ∩Y=? 的充要條件為

        2 全對稱多胞形卡爾曼濾波器設計

        本文主要研究基于狀態(tài)集合交集為空設計輔助輸入信號的方法.由于狀態(tài)量無法實時獲取,本節(jié)設計一個全對稱多胞形卡爾曼濾波器對狀態(tài)集合進行估計,并利用全對稱多胞形對受干擾影響的狀態(tài)集合進行描述.

        假設系統(tǒng)(1)的初始狀態(tài),過程干擾和測量噪聲分別滿足:

        針對該系統(tǒng),設計全對稱多胞形卡爾曼濾波器形式如下:

        定理1.設初始時刻的狀態(tài)滿足已知k 時刻估計得到的狀態(tài)為則k+1 時刻的狀態(tài)為

        其中,

        證明.由式(15)可得:

        根據(jù)全對稱多胞形的性質以及引理2,化簡式(17)可得k+1 時刻估計得到的狀態(tài)所在的全對稱多胞形的中心及生成矩陣□

        定理2.依據(jù)卡爾曼濾波器設計方法,針對狀態(tài)全對稱多胞形(16),設計的增益矩陣為式(18)~(23),則觀測器(15)可以實時估計系統(tǒng)(1)的狀態(tài)所在的集合.

        證明.定義為評判全對稱多胞形大小的指標,當時得到最優(yōu)的觀測器增益具體證明過程詳見文獻[30],本文不再贅述. □

        3 最優(yōu)輔助輸入信號設計

        本節(jié)利用估計的狀態(tài)集合來設計輔助輸入信號,通過引入輔助輸入信號去激勵實際系統(tǒng),從而實現(xiàn)故障檢測.由于輔助輸入信號的加入很可能會影響系統(tǒng)的實際運行狀態(tài),需要設計的輔助輸入信號越小越好.本節(jié)將輔助輸入信號的設計轉化為求解混合整數(shù)二次規(guī)劃問題,得到最優(yōu)的輔助輸入信號.

        為了簡化后文計算過程的表達式,定義

        定理3.針對系統(tǒng)(1),設計輔助輸入信號uk,滿足

        則在輔助輸入信號uk影響下的正常模型與故障模型狀態(tài)集合的交集為空,即

        其中,△[0i]=B[i]-B[0],Zm(0,i)=M[0]⊕(-M[i]).

        證明.由定理1 以及式(27)整理得:

        根據(jù)引理3,需要滿足

        經(jīng)整理可得式(27). □

        為了使得所設計的輔助輸入信號對系統(tǒng)的影響最小,需所設計的輔助輸入信號最小,因此通過求解以下優(yōu)化問題,得到輔助輸入信號:

        其中,R 為半正定矩陣.依據(jù)文獻[25],由于優(yōu)化問題(31)是一個非凸優(yōu)化問題,不易得到最優(yōu)解,為了得到有效的求解方法,結合定理3,本文對優(yōu)化問題(31)進行重構,將式(31)轉化為混合整數(shù)二次規(guī)劃問題進行求解.

        其中

        其中

        當δ[0i]≤0 時,則不滿足△[0i]uk的條件,只有當δ[0i]>0 時,才能滿足該條件.

        因此對于每個uk,則可以將優(yōu)化問題(31)重構為線性規(guī)劃問題

        進一步,優(yōu)化問題(36)可以轉化成

        注意到式(39)為雙層優(yōu)化問題,不易得到最優(yōu)解,本文通過引入拉格朗日乘子使得式(39)中的優(yōu)化問題轉化為單優(yōu)化問題.

        定理4.對于系統(tǒng)(1),設計輔助輸入信號uk,使其滿足以下條件

        則引入輔助輸入信號uk后,正常模型和故障模型狀態(tài)所在的集合交集為空.其中,∈R2r為拉格朗日乘子,為二進制變量,向量111=[11 ··· 1]T.

        證明.設拉格朗日乘子分別為

        其中,

        構造拉格朗日函數(shù)為

        為了求得δ[0i]的局部最小值,利用文獻[31]中的定理3.3.1和3.4.1,通過對拉格朗日函數(shù)局部求導,可以分別得到

        化簡式(45)和式(46),可得:

        由于式(47)和式(48)是非凸函數(shù),所以不能作為優(yōu)化問題的約束函數(shù).因此,通過引入二進制變量則式(47)和式(48)可以轉化為

        結合式(52),式(35)可以化為

        可以得到

        同理可得

        因此,整理式(56)和式(57)可以得到以下邊界條件

        通過式(58)中的邊界條件,以及在式(51)中引入二進制,可以將式(47)和式(48)轉化為

        整理式(39),(49),(50),(59),可得優(yōu)化問題式(40). □

        利用正常模型和故障模型交集為空的思想設計得到輔助輸入信號以后,全對稱多胞形的故障檢測邏輯主要的目標是判斷實際系統(tǒng)估計得到的狀態(tài)是否發(fā)生故障,則判斷系統(tǒng)是否發(fā)生故障可以轉化成以下邏輯問題:

        其中,fresult=0 時,系統(tǒng)處于正常情況下,否則,系統(tǒng)發(fā)生故障.式(60)中的邏輯判斷條件根據(jù)全對稱多胞形的定義可以轉換為以下約束問題:

        4 仿真驗證

        為了驗證本文所提方法的正確性與有效性,采用文獻[25]中的永磁直流電動機的低頻線性模型進行仿真.該系統(tǒng)模型為

        其中,u,i,Ra,L,Ke,Kt,J1,fr分別為電樞電壓,電流,電阻,電感,轉矩常數(shù),反電動勢常數(shù),電機慣量和摩擦系數(shù),轉矩常數(shù)Kt=1.0005×Ke.u=uc+ua,其中,uc為系統(tǒng)的控制輸入,ua為系統(tǒng)的輔助輸入.電機處在正常情況下系統(tǒng)參數(shù)分別為Ra=1.2030 Ω,L=5.5840×10-3H,Ke=8.1876 ×10-2V·rad/s,J1=1.3528×10-4J·s2/rad,fr=2.3396×10-4J·s/rad.當電機電感發(fā)生變化時,即L=8.7548×10-3H 時,系統(tǒng)發(fā)生故障.根據(jù)電動機的機理過程,為將電動機的轉速維持在70.3 rad/s,給定系統(tǒng)的控制輸入uc=6 V.利用前向歐拉差分方法將電動機模型離散化,其中采樣時間為5 ms.由于離散過程會使得參數(shù)引入不確定性,因此離散化的狀態(tài)動態(tài)過程增加Ewwwk,其中

        假設系統(tǒng)的初始狀態(tài)、測量噪聲和過程干擾所在的全對稱多胞形分別滿足

        系統(tǒng)的初始輔助輸入信號給定為ua=0.通過求解最優(yōu)化問題(40),可得輔助輸入信號,在求解優(yōu)化問題(40)時,取最小分離閾值最大分離閾值

        給定實際系統(tǒng),假設實際系統(tǒng)前50 時刻發(fā)生故障的情況如下:

        為了驗證所提方法的有效性,本文與基于輸出集合的輔助輸入信號方法進行對比.

        圖3 是兩種方法得到的輔助輸入信號.其中實線為基于本文所提方法所得的輔助輸入信號,虛線為基于輸出集合交集所得到的輔助輸入信號.由圖3可以較為明顯地看出,在系統(tǒng)穩(wěn)定后本文提出的方法得到的輔助輸入信號的幅值相對較小.

        圖3 輔助輸入信號Fig.3 The auxiliary input signal

        圖4 給出了不同方法得到的輔助輸入信號對系統(tǒng)輸出影響的對比.其中實線表示系統(tǒng)未加輔助輸入信號時系統(tǒng)的狀態(tài)以及相應的輸出,虛線表示文中所提方法得到的輔助輸入信號對系統(tǒng)狀態(tài)和輸出影響,點畫線表示利用輸出集合得到的輔助輸入信號對系統(tǒng)狀態(tài)和輸出的影響.可以看出本文所提出的方法對原系統(tǒng)的影響較小.

        圖4 加入輔助信號和未加輔助信號對系統(tǒng)輸出的影響Fig.4 Effect of adding auxiliary signal and unassisted signal on system output

        圖5 分別為k=20,21,22,23 時刻以及k=40,41,42,43 時刻利用文中所提方法得到的輔助輸入信號,將該輔助輸入信號注入系統(tǒng)后,通過觀測器觀測得到的正常模型和故障模型的狀態(tài)所在的全對稱多胞形的分離情況.其中,三角符號以及方形符號分別為相應全對稱多胞形的中心,圓圈為實際系統(tǒng)觀測得到的狀態(tài).根據(jù)式(64)可得,系統(tǒng)在k=20時刻發(fā)生故障,由圖5 可以看出,實際系統(tǒng)在k=21 時刻落在故障全對稱多胞形中,可以判斷系統(tǒng)在k=21 處在故障的情況下.同理看得,實際系統(tǒng)在k=41 時刻落在正常全對稱多胞形中.可以判斷系統(tǒng)在k=41 處在正常的情況下.由于實際系統(tǒng)在k=20 時刻已經(jīng)發(fā)生故障,所以系統(tǒng)由正常到故障不會突變,都有個緩慢的過程,所以在k=20 時刻時未檢測出故障,故障檢測結果有一個時刻的延遲.圖6為利用全對稱多胞形故障檢測邏輯得到的故障檢測結果,可以直觀地看出利用文中方法對實際系統(tǒng)的故障檢測結果.

        圖5 狀態(tài)多胞形:通過狀態(tài)交集為空的方法Fig.5 Zonotope of the system state:the method of emptying the state intersection

        圖6 故障檢測結果Fig.6 The results of fault detection

        為了更好地展現(xiàn)出利用文中所提方法比基于輸出集合交集為空設計輔助輸入信號的保守性小,給出了k=20,21,22,23 相應的全對稱多胞形,見仿真圖7和圖8.圖7和圖8 為基于輸出集合交集為空的主動故障檢測方法所得到的結果.圖7 為系統(tǒng)輸出所在的全對稱多胞形的分離情況,圖8 為系統(tǒng)輸出所在全對稱多胞形剛好相切的情況下,對應的狀態(tài)所在全對稱多胞形的情況.可以看出正常模型和故障模型的輸出所在的全對稱多胞形相切時,此時狀態(tài)所在的全對稱多胞形已經(jīng)相離.因此,這說明從系統(tǒng)狀態(tài)的角度出發(fā),設計輔助輸入信號去激勵系統(tǒng)實現(xiàn)故障分離,該方法要比從系統(tǒng)輸出角度出發(fā),設計得到的輔助信號小,保守性小.

        圖7 輸出多胞形:通過輸出交集為空的方法Fig.7 Zonotope of the system output:the method of emptying the output intersection

        圖8 狀態(tài)多胞形:通過輸出交集為空的方法Fig.8 zonotope of the system state:the method of emptying the output intersection

        5 結論

        本文在系統(tǒng)過程干擾和測量噪聲未知但有界的前提下,設計了一種基于狀態(tài)估計的輔助輸入信號方法.首先利用全對稱多胞形卡爾曼濾波設計系統(tǒng)的狀態(tài)觀測器,對系統(tǒng)的狀態(tài)進行估計,得到實時狀態(tài)所在的集合.之后在對系統(tǒng)影響最小的前提下,利用正常模型的狀態(tài)所在的多胞形和故障模型的狀態(tài)所在的多胞形的交集為空,設計得到實時的最優(yōu)輔助輸入信號.將該輔助輸入信號注入系統(tǒng),充分激勵系統(tǒng),提高系統(tǒng)的故障檢測質量.最后通過分析數(shù)值模型驗證文章所提方法.本文所用到的設計狀態(tài)觀測器的方法是通過實時最小化狀態(tài)所在多胞形的F半徑,該方法雖然可以實時得到狀態(tài)觀測器增益,但在進行狀態(tài)估計時仍不能十分準確地估計狀態(tài),因此如何選擇狀態(tài)觀測器較為準確估計狀態(tài)仍需進一步進行研究.

        猜你喜歡
        故障信號檢測
        “不等式”檢測題
        “一元一次不等式”檢測題
        “一元一次不等式組”檢測題
        信號
        鴨綠江(2021年35期)2021-04-19 12:24:18
        完形填空二則
        故障一點通
        基于FPGA的多功能信號發(fā)生器的設計
        電子制作(2018年11期)2018-08-04 03:25:42
        奔馳R320車ABS、ESP故障燈異常點亮
        小波變換在PCB缺陷檢測中的應用
        基于LabVIEW的力加載信號采集與PID控制
        日韩av一区二区网址| 婷婷五月亚洲综合图区| 国产一区二区在线观看视频免费| 日本久久精品福利视频| av中文字幕潮喷人妻系列| 国产精品亚洲成在人线| 亚洲av日韩片在线观看| 日韩成人高清不卡av| 精品国产粉嫩内射白浆内射双马尾| 国产丝袜在线精品丝袜| 国产人成无码视频在线| 激情视频在线播放一区二区三区| 91成人自拍国语对白| 久激情内射婷内射蜜桃人妖| 国产午夜亚洲精品理论片不卡 | av大全亚洲一区二区三区| 特级无码毛片免费视频尤物| 大陆一级毛片免费播放| 久久蜜臀av一区三区| 丁香五月缴情在线| 无码精品人妻一区二区三区人妻斩 | 亚洲中文字幕无码中文字| 欧美性群另类交| 极品人妻少妇一区二区| 日本一区二区三区视频免费在线| 亚洲国产成人久久综合| 在线视频你懂的国产福利| 国产一区不卡视频在线| 乱码窝窝久久国产无人精品| 99re热视频这里只精品| 2021国产最新无码视频| 熟女少妇精品一区二区三区| 挺进邻居丰满少妇的身体| 精品国产午夜福利在线观看| 亚洲av乱码一区二区三区女同| 爆操丝袜美女在线观看| 无码人妻av一区二区三区蜜臀| 亚洲av乱码专区国产乱码| 情av一区二区三区在线观看| 国产精品女人呻吟在线观看| 69精品免费视频|