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

        ?

        向量自回歸和T2控制圖在變工況齒輪箱故障檢測中的應(yīng)用

        2023-09-13 12:52:42毛雨琳
        關(guān)鍵詞:故障診斷振動故障

        薛 成,李 鑫,2,焦 靜,毛雨琳

        (1.南京工程學(xué)院 汽車與軌道交通學(xué)院, 南京 211167;2.南京航空航天大學(xué) 民航學(xué)院, 南京 211106)

        0 引言

        齒輪箱以其結(jié)構(gòu)緊湊、傳動比高、承載能力大的獨特特點被廣泛地集成到各種旋轉(zhuǎn)機(jī)械中。由于工作環(huán)境惡劣,齒輪箱容易損壞,導(dǎo)致整個機(jī)械系統(tǒng)性能下降甚至停機(jī)。為了有效地提高機(jī)器運行的可靠性和降低維護(hù)成本,對齒輪箱進(jìn)行在線監(jiān)測和故障診斷是十分必要的[1-2]。

        基于振動信號的分析是最常見的齒輪箱故障診斷方法,但在面對多工況交替等情況時,其信號波動幅度大,各個工況的振幅不同,導(dǎo)致監(jiān)測診斷能力有所下降。然而,振動信號的特征階次在非平穩(wěn)工況下,不會隨著轉(zhuǎn)速的變化而發(fā)生改變。因此,應(yīng)用階次分析的方法可以有效進(jìn)行非平穩(wěn)工況下齒輪箱的故障診斷。趙登利等[3]推導(dǎo)并總結(jié)了齒輪箱特征階次,并通過真實風(fēng)場數(shù)據(jù)進(jìn)行了故障診斷分析與實驗驗證。王博等[4]根據(jù)齒輪轉(zhuǎn)速對故障頻率的影響,總結(jié)了齒輪箱各部件的故障階次特征表,實現(xiàn)了變轉(zhuǎn)速工況下的齒輪箱故障診斷。唐敏[5]也運用階次跟蹤技術(shù)將變速齒輪箱的時域非平穩(wěn)特征轉(zhuǎn)化為平穩(wěn)的瞬態(tài)特征,優(yōu)化了變速齒輪箱的狀態(tài)評估效率。杜麗君等[6]提出了一種時域稀疏調(diào)頻字典的齒輪復(fù)合故障診斷方法,從變速齒輪箱復(fù)雜的信號中準(zhǔn)確提取出故障信號分量,確保了后續(xù)信號的幅值重構(gòu)精度,完成了復(fù)合故障信號的分離診斷。張小剛等[7]提出了一種深度殘差網(wǎng)絡(luò)的方法,能夠更加有效地挖掘信號特征之間的信息,提高診斷精度。

        程龍歡等[8]則認(rèn)為,單源的振動信號不足以體現(xiàn)設(shè)備的運行狀態(tài),應(yīng)該將多個振動傳感器進(jìn)行組合使用,融合多個傳感器的振動信號,才能夠正確表明設(shè)備的狀態(tài)。為了解決單一故障特征診斷準(zhǔn)確性低的問題,張智禹等[9]提出了融合注意力機(jī)制的改進(jìn)深度置信網(wǎng)絡(luò)的變工況齒輪箱故障診斷方法,以融合時域、頻域、小波包時頻域特征作為故障診斷參數(shù),有效提高了故障診斷精確度和泛化能力。Kumar等[10]將多傳感器信息融合技術(shù)應(yīng)用在齒輪箱故障診斷中,使用振動加速度計、麥克風(fēng)和聲發(fā)射傳感器從齒輪箱中獲取不同加載工況和轉(zhuǎn)速變化下的機(jī)器狀態(tài)數(shù)據(jù),對采集到的振動、聲音和聲發(fā)射信號進(jìn)行小波特征提取,識別出突出特征,并采用機(jī)器學(xué)習(xí)算法以實現(xiàn)故障診斷過程的自動化。

        統(tǒng)計過程控制可以有效監(jiān)測數(shù)據(jù)的隨機(jī)波動和異常波動狀況。馬曉原[11]運用統(tǒng)計過程控制方法對數(shù)控機(jī)床加工單零件的能耗進(jìn)行控制,實現(xiàn)了單零件能耗的提前預(yù)警。梁坤等[12]基于統(tǒng)計過程方法,采用指數(shù)加權(quán)移動平均(exponential weighted moving average,EWMA)控制圖,提高了故障預(yù)警的準(zhǔn)確度。

        本研究將時間序列模型與統(tǒng)計過程控制方法相結(jié)合,進(jìn)行了變工況下齒輪箱的異常檢測。在時域同步平均算法(time synchronous averaging,TSA)預(yù)處理振動信號的基礎(chǔ)上,選取健康數(shù)據(jù)建立向量自回歸模型(vector autoregression,VAR),提取特征參數(shù),從而建立統(tǒng)計過程控制圖,監(jiān)測數(shù)據(jù)的波動情況,評估提取的故障特征值,確定齒輪箱在正常工作狀態(tài)下故障特征值的幅值,準(zhǔn)確預(yù)測、判斷齒輪箱異常發(fā)生點,及時地進(jìn)行故障預(yù)警。

        1 數(shù)據(jù)獲取

        1.1 數(shù)據(jù)來源

        試驗數(shù)據(jù)來源于如圖1所示的賓夕法尼亞大學(xué)的機(jī)械故障診斷試驗臺[13]。通過此試驗臺讀取轉(zhuǎn)速與扭矩數(shù)據(jù),可對齒輪箱壽命進(jìn)行可靠性驗證試驗。

        圖1 試驗臺設(shè)備圖

        試驗所用的齒輪箱為14號齒輪箱,該齒輪箱額定輸入轉(zhuǎn)速為1 750 r/min,最大額定輸出扭矩為555 in-lbs(1 in-lbs=0.113 N·m),最大額定輸入功率為4.66 hp(1 hp=745.7 W),包含一個齒數(shù)為70的驅(qū)動齒輪和一個齒數(shù)為21的小齒輪,齒輪嚙合頻率為613 Hz。齒輪箱振動數(shù)據(jù)來源于三軸加速度傳感器。測試的振動數(shù)據(jù)以離散的間隔采樣,每個樣本被保存為一個文件[14]。

        圖2是從A10、A11以及A12三軸加速度計中采集的齒輪箱原始的振動數(shù)據(jù)曲線,其中A10設(shè)置為軸向,A11與地面方向平行,而A12則與地面方向垂直[15]。觀察原始振動數(shù)據(jù),發(fā)現(xiàn)數(shù)據(jù)點的數(shù)量十分巨大,前96 h,輸出扭矩恒定在100%時,振幅穩(wěn)定,但工況開始變化以后,振幅變化不明顯,無法提取有用的信號來進(jìn)行異常檢測。

        圖2 整個試驗過程的原始振動數(shù)據(jù)曲線

        1.2 采樣過程

        齒輪箱全壽命周期試驗在不同工況時分別采樣一次,采樣頻率為20 kHz,每次采樣時間為10 s。試驗在2種工況下進(jìn)行:轉(zhuǎn)速恒定不變,齒輪箱在100%額定扭矩下運行96 h;然后以50%、100%、150%、200%、250%、300%扭矩和功率運行了19.3 h,直到出現(xiàn)故障,失效停機(jī)。本試驗在100%額定扭矩工況下,齒輪箱正常工作。在變工況狀態(tài)下,齒輪箱出現(xiàn)異常直至故障停機(jī),所以只需選取變工況下齒輪箱的振動信號進(jìn)行分析。TR#14齒輪箱在不同負(fù)載扭矩下共采樣152個文件。將故障停機(jī)的齒輪箱拆開進(jìn)行檢查,檢查發(fā)現(xiàn)小齒輪正常,驅(qū)動齒輪有5個完全損壞的齒和2個部分損壞的齒。

        2 數(shù)據(jù)預(yù)處理

        2.1 TSA

        TSA算法選取信號時應(yīng)按照齒輪的旋轉(zhuǎn)周期進(jìn)行,提取齒輪旋轉(zhuǎn)的旋轉(zhuǎn)頻率與嚙合頻率,然后將同步相加、平均[16]。隨著同步平均次數(shù)增加,噪聲干擾信號逐漸減少,直到不足以對故障診斷產(chǎn)生影響。假設(shè)每個振動數(shù)據(jù)文件中有n個數(shù)據(jù)點。研究齒輪旋轉(zhuǎn)一周的完整的采樣點數(shù)由以下公式求得:

        (1)

        式中:fs為取樣頻率,fm為研究齒輪的基本嚙合頻率,Nt為齒輪的齒數(shù),「·?表示對計算結(jié)果向上取整。

        需要平均的循環(huán)次數(shù)可以由下公式求得:

        (2)

        式中:N為每個數(shù)據(jù)文件中采樣點的總數(shù),?·」為對計算結(jié)果向下取整。

        時域信號中,振動數(shù)據(jù)中的TSA信號可由下式計算:

        (3)

        式中:K為齒輪每運行一個周期的采樣點數(shù),M0為單個采樣文件所包含的齒輪的轉(zhuǎn)數(shù)。

        將原始振動信號導(dǎo)入Matlab軟件中,利用TSA算法得到TSA處理曲線,如圖3所示。通過TSA信號的健康數(shù)據(jù)可以建立VAR模型來擬合數(shù)據(jù),再將原始數(shù)據(jù)與擬合數(shù)據(jù)相減得到殘差。

        圖3 變工況下A10、A11、A12 TSA信號

        2.2 分工況處理

        在對原始振動信號進(jìn)行TSA處理之后,由于不同工況的扭矩負(fù)載不同,因此需要對數(shù)據(jù)進(jìn)一步按不同的工況進(jìn)行分類。文件序號對應(yīng)的輸出扭矩值如圖4所示。

        圖4 齒輪箱輸出扭矩曲線

        將相同工況的文件中的振動數(shù)據(jù)組合,即可得到各個工況下的TSA振動數(shù)據(jù),其中比較典型的A12分工況處理后的TSA振動數(shù)據(jù)如圖5所示。

        圖5 分工況TSA處理后A12的振動數(shù)據(jù)曲線

        2.3 VAR模型及殘差計算

        VAR模型是一種常用的計量經(jīng)濟(jì)模型,常使用模型中所有當(dāng)期變量對所有變量的若干滯后變量進(jìn)行回歸。VAR模型可以在沒有預(yù)先限制的情況下,對聯(lián)合內(nèi)生變量進(jìn)行動態(tài)預(yù)測,實現(xiàn)對2個或多個相關(guān)聯(lián)的時間序列的預(yù)測。VAR模型描述在同一樣本期間內(nèi)的n個變量(內(nèi)生變量)可以作為它們過去值的線性函數(shù)[17]。一個VAR(p)模型可以寫成:

        yt=c+A1yt-1+A2yt-2+…+Apyt-p+εt

        (4)

        式中:c為常數(shù)向量,yt為n維內(nèi)生變量,Ai(i=1,2,…,p)為n×n矩陣,εt為誤差向量。

        選取經(jīng)TSA處理后健康狀態(tài)下的數(shù)據(jù)作為訓(xùn)練數(shù)據(jù),使用Eviews軟件建立VAR模型。因為齒輪箱大部分時間處于正常情況下運行,為減小計算量,只需選取一部分平穩(wěn)數(shù)據(jù)進(jìn)行建模。選取前5組數(shù)據(jù)TSA信號建立VAR模型。

        數(shù)據(jù)序列經(jīng)過檢驗后證明平穩(wěn)。有代表性的300%工況下建立的VAR模型的階數(shù)確定如表1所示,根據(jù)各個信息準(zhǔn)則的參數(shù),綜合考慮確定該VAR模型的滯后階數(shù)為47。

        表1 300%工況建立的VAR模型的階數(shù)

        通過單位圓檢驗對該模型進(jìn)行穩(wěn)定性判斷,如圖6所示,所有的點均處于單位圓內(nèi),因此判斷該VAR模型穩(wěn)定。

        圖6 300%工況建立的VAR模型穩(wěn)定性判斷

        模型的表達(dá)式如下:

        A10=1.334 969 002 31×A10(-1)-

        0.563 750 025 639×A10(-2) +…+

        0.001 687 210 172 87×

        A12(-47)-1.718 841 198 06

        A11=0.088 596 424 477×A10(-1)-

        0.163 006 450 186×A10(-2) +…-

        0.007 271 001 622 06×

        A12(-47)+6.807 399 845 83

        A12=-0.065 850 411 545×A10(-1) +

        0.061 530 633 529 9×A10(-2) +…+

        0.079 860 198 146 3×

        A12(-47)-1.704 486 071 05

        分析殘差需要選取全部的文件,通過VAR模型擬合的數(shù)據(jù)與原始數(shù)據(jù)進(jìn)行對比,整體分析判斷齒輪箱的健康狀態(tài),得到的殘差圖如圖7所示。

        圖7 殘差曲線

        2.4 特征值提取

        常用的時域信號的有量綱特征值有最大值、最小值、峰值、整流平均值、均值、方差、均方根值;無量綱特征值有值因子、脈沖因子、裕度因子、峭度因子、波形因子和偏度等。

        整流平均值可表示為信號絕對值積分的平均值:

        (5)

        均方根(root mean square,RMS)又稱有效值。將所有數(shù)據(jù)平方求和,求其均值,再開平方,就得到均方根值,其計算公式如下:

        (6)

        有量綱指標(biāo)盡管對信號特性非常敏感,但是在工作條件(例如負(fù)荷)變化時也會發(fā)生改變,而且容易受到外界的干擾,因此其性能不穩(wěn)定。相對來說,無量綱指標(biāo)可以消除干擾,因而在特征抽取方面得到了廣泛的應(yīng)用。

        峭度因子是四階中心矩和標(biāo)準(zhǔn)差的四次方的比值,也是對振動信號沖擊特性的反應(yīng)。偏度和峭度有一定的相關(guān)性,偏度因子是三階中心矩和標(biāo)準(zhǔn)差的三次方的比值,描述的是分布。

        為了確保提取的特征值能夠較為準(zhǔn)確地反映出齒輪箱的運行狀態(tài),選取了有量綱特征值整流平均值、均方根以及無量綱特征值偏度、峭度因子進(jìn)行比較。以齒輪箱變工況運行時A10的振動數(shù)據(jù)為例,提取了其VAR模型與TSA數(shù)據(jù)殘差的特征值,如圖8—11所示。

        圖8 整流平均值曲線

        齒輪箱在變工況狀態(tài)下運行,其振動數(shù)據(jù)應(yīng)具有一定的規(guī)律性,同一種工況時的振動數(shù)據(jù)應(yīng)大致相等或隨時間而略有增加,振動數(shù)據(jù)應(yīng)隨著扭矩增加或減小而相應(yīng)的增加或減小。根據(jù)上述判斷,觀察選取符合實際情況的特征值圖像。觀察圖8數(shù)據(jù),曲線具有規(guī)律性,隨時間穩(wěn)定上升,符合推斷。圖9曲線與圖8曲線特征相似,符合推斷。

        圖9 均方根曲線

        圖10曲線波動幅度基本相同,不符合推斷。圖11曲線復(fù)雜無規(guī)律,難以找到特征。

        圖10 偏度因子

        圖11 峭度因子

        因此,選用整流平均值以及均方根這2個特征指標(biāo)對數(shù)據(jù)進(jìn)行提取。

        殘差是實際數(shù)據(jù)與擬合數(shù)據(jù)的差值,采用時域分析將殘差數(shù)據(jù)中的故障特征提取可以進(jìn)一步實現(xiàn)對齒輪箱的故障診斷。通過特征值的波動情況可以大概判斷齒輪箱的異常情況。利用時域分析計算殘差的整流平均值(圖12)和均方根(圖13)指標(biāo)來初步判斷齒輪的故障發(fā)生點。

        圖12 殘差的整流平均值

        圖13 殘差的均方根值

        通過對圖12、圖13的觀察,發(fā)現(xiàn)齒輪箱在3個方向上的振動信號基本穩(wěn)定,但在曲線末端都有明顯的上升,曲線十分陡峭,可以確定有故障出現(xiàn),但在其他位置處是否產(chǎn)生故障無法判斷,接下來進(jìn)行統(tǒng)計過程分析,從而有效地找到異常點。

        3 統(tǒng)計過程控制

        控制圖是通過測量、記錄和評價工藝質(zhì)量特性,從而監(jiān)測工藝是否在控制中的一種圖表。它可以區(qū)別正常和不正常的波動,是一種重要的統(tǒng)計手段。在這張圖表上,有3條與橫向軸線平行的線:中心線(central line,CL)、上控制限(upper control limit,UCL)、下控制限(lower control limit,LCL)以及按時間順序提取的采樣統(tǒng)計數(shù)字的描點序列。UCL、CL和LCL共同稱為控制極限,一般將控制極限設(shè)置在±3個標(biāo)準(zhǔn)偏差處。CL是統(tǒng)計數(shù)據(jù)的平均值,其上、下控制邊界和CL之間有幾個標(biāo)準(zhǔn)偏差。如果在控制圖中的描點處在UCL和LCL以外,或者在UCL和LCL之間沒有任意排列,那么說明該過程是不正常的[17]。

        EWMA控制圖可以檢測過程均值中的小偏差,而不會受到低值和高值的影響。由圖3中TSA信號與圖7中殘差值可以明顯得出,A10振動信號振幅大于A11與A12,波動強(qiáng)于A11與A12。選取A10的VAR模型與TSA殘差的整流平均值與均方根數(shù)據(jù),繪制EWMA控制圖。在EWMA控制圖選項中使用默認(rèn)設(shè)定,即距離CL大于3個標(biāo)準(zhǔn)差為檢驗不合格。EWMA控制圖如圖14、圖15所示。

        圖14 整流平均值EWMA控制圖

        圖15 均方根EWMA控制圖

        圖14與圖12中的曲線形狀趨勢類似,在圖14中檢測到多個點不合格,A10振動信號TSA數(shù)據(jù)的整流平均值穩(wěn)定性較差,可能導(dǎo)致故障檢測結(jié)果準(zhǔn)確率嚴(yán)重下降。圖15與圖13中的曲線類似,在圖15中檢測出5個點不合格,且曲線末端的2個點確定發(fā)生故障。初步推斷均方根數(shù)據(jù)用于故障檢測效果更好,更準(zhǔn)確的異常檢測結(jié)果還需進(jìn)一步觀察三維數(shù)據(jù)融合后的多變量控制圖。

        MEWMA控制圖是EWMA控制圖的多變量形式,可以在EWMA控制圖中同時監(jiān)視2個或多個相關(guān)過程特征。MEWMA控制圖中,每個標(biāo)繪點都包括來自所有以前數(shù)據(jù)的加權(quán),可以檢測較小的過程偏移。

        MEWMA統(tǒng)計量的計算方法如下:

        (7)

        式中:Zi表示先前子組的權(quán)重:

        Zi=rXi+(1-r)Zi

        (8)

        加權(quán)的協(xié)方差矩陣為:

        (9)

        式中:

        (10)

        選取A10、A11和A12的VAR模型與TSA殘差的整流平均值與均方根,繪制三維變量的MEWMA控制圖,如圖16、圖17所示。

        圖16 整流平均值MEWMA控制圖

        圖17 均方根MEWMA控制圖

        圖16中在曲線的首尾兩端檢測出多個異常點,且超過控制限距離較大,與實際不符。圖17中雖然同樣檢測出多個異常點,但大部分異常點都接近控制限,這是由于MEWMA控制圖常用于檢測數(shù)據(jù)中較小的過程偏移,導(dǎo)致其敏感度較高,最終造成虛警率過高。

        廣義方差控制圖可以同時監(jiān)視2個或多個相關(guān)變量的過程變異性是否受控制,此控制圖是極差控制圖、標(biāo)準(zhǔn)差控制圖和移動極差控制圖的多變量形式。選取A10、A11和A12的VAR模型與TSA殘差的整流平均值與均方根,繪制三維變量的廣義方差控制圖,如圖18、圖19所示。

        圖18 整流平均值廣義方差控制圖

        圖19 均方根廣義方差控制圖

        圖18中在曲線首端檢測出2個異常點,然而齒輪箱在試驗開始階段正常運行,因此不符合實際。圖19中在曲線前段檢測出49、60 2個異常點,在曲線后段檢測出102、103、114、125 4個異常點,最后在曲線末端檢測出151、152 2個異常點。然而,齒輪箱在試驗中一直運行直至故障停機(jī),控制圖中不應(yīng)出現(xiàn)較多的異常點。因此,均方根廣義方差控制圖的虛警率略高。

        使用T2控制圖可以同時監(jiān)視2個或多個相關(guān)變量的過程位置是否受控制。此控制圖是均值控制圖和單值控制圖的多變量形式。

        當(dāng)數(shù)據(jù)位于子組中時,T2的計算方式如下:

        (11)

        式中:

        (12)

        (13)

        T2控制圖的控制上限為:

        (14)

        圖20是選取了A10、A11和A12的VAR模型與TSA殘差的整流平均值,繪制了三維變量的T2控制圖。該控制圖在齒輪箱運行早期檢測出多個異常點,不符合實際。

        圖20 整流平均值T2控制圖

        圖21是選取了A10、A11和A12的VAR模型與TSA殘差的均方根數(shù)據(jù),繪制了三維變量的T2控制圖。

        圖21 均方根T2控制圖

        均方根T2控制圖檢驗結(jié)果如表2所示。

        表2 均方根異常結(jié)果

        表2表明均方根T2控制圖檢測出了3個異常點,文件編號分別為60、151、152,結(jié)合圖21進(jìn)行分析,得出齒輪箱故障診斷結(jié)論:齒輪箱在文件編號為60時,出現(xiàn)微小故障,產(chǎn)生了早期微弱的故障信號。在文件編號為151、152時,齒輪箱振動幅度巨大,出現(xiàn)故障,失效停機(jī)。

        4 結(jié)論

        1) 以14號齒輪箱為研究對象,采用TSA算法去除振動信號中的背景噪聲。對振動文件進(jìn)行分工況處理后,選取部分健康數(shù)據(jù)建立VAR模型,得到整體數(shù)據(jù)的殘差,并提取殘差的整流平均值、均方根指標(biāo)。

        2) 以殘差均方根和整流平均值指標(biāo)建立統(tǒng)計過程控制曲線,比較MEWMA、廣義方差以及T2控制曲線,最終確定均方根T2控制圖檢測效果最佳,判斷齒輪箱異常位置為第60、151、152個文件。以多維數(shù)據(jù)進(jìn)行分析,可以更為準(zhǔn)確地判斷出異常點,不僅能找出故障出現(xiàn)的異常點,也能捕捉到早期微弱的故障信號。

        3) 相比于傳統(tǒng)的針對單一工況下的齒輪箱故障診斷方法,本研究對齒輪箱振動數(shù)據(jù)進(jìn)行了分工況處理,可以判斷齒輪箱全運行狀態(tài)下的健康狀態(tài),提升了應(yīng)用性。

        猜你喜歡
        故障診斷振動故障
        振動的思考
        振動與頻率
        故障一點通
        中立型Emden-Fowler微分方程的振動性
        奔馳R320車ABS、ESP故障燈異常點亮
        因果圖定性分析法及其在故障診斷中的應(yīng)用
        故障一點通
        江淮車故障3例
        UF6振動激發(fā)態(tài)分子的振動-振動馳豫
        計算物理(2014年2期)2014-03-11 17:01:44
        基于LCD和排列熵的滾動軸承故障診斷
        日日碰日日摸日日澡视频播放| 91极品尤物在线观看播放| 中文亚洲第一av一区二区| 国产精品亚洲三级一区二区三区| 妺妺窝人体色www聚色窝| 亚洲妓女综合网99| 色综合久久久久综合一本到桃花网| 国产中文色婷婷久久久精品| 台湾佬中文网站| 久久精品国内一区二区三区| 四虎国产精品免费久久麻豆| av天堂手机在线看片资源| 观看在线人视频| 无码人妻丰满熟妇啪啪7774| 精品国产1区2区3区AV| 大陆国产乱人伦| 亚洲av精品一区二区| 久久国产精品免费专区| 国产丝袜美女| 天堂网在线最新版www中文网| 亚洲av永久无码精品水牛影视| 日本本土精品午夜视频| 日本大骚b视频在线| 亚洲永久无码动态图| 午夜香蕉av一区二区三区| 产美女被爽到高潮免费a| 久久精品麻豆日日躁夜夜躁| 女的把腿张开男的猛戳出浆| 久久开心婷婷综合中文| 亚洲精品无码永久在线观看| 亚洲的天堂av无码| 美腿丝袜一区二区三区| 偷拍一区二区三区四区| 乱人伦中文无码视频在线观看 | 先锋五月婷婷丁香草草| 国产精品jizz观看| 久久精品国产精品亚洲婷婷| 日本添下边视频全过程| 成人做爰高潮尖叫声免费观看| 中国免费一级毛片| 色婷婷精品久久二区二区蜜臀av|