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

        ?

        基于稀疏分解的軸承聲陣列信號特征提取

        2018-08-25 07:29:22郭瑩瑩趙學智上官文斌張春良
        振動、測試與診斷 2018年4期
        關(guān)鍵詞:裂紋故障信號

        郭瑩瑩, 趙學智, 上官文斌, 張春良

        (1.華南理工大學機械與汽車工程學院 廣州,510641) (2.廣州大學機械與電氣工程學院 廣州,510006)

        引 言

        滾動軸承是機械設(shè)備中最重要的部件之一,其狀態(tài)的好壞直接關(guān)系到設(shè)備的正常運轉(zhuǎn),因此對滾動軸承的狀態(tài)監(jiān)測和故障診斷一直是非常重要的研究課題[1]。利用振動傳感器檢測軸承振動,并通過對振動信號的特征提取來識別軸承的故障狀況是一種應(yīng)用最為廣泛的方法,其中的特征提取方法主要有小波及小波包分析[2-3]、經(jīng)驗?zāi)J椒纸鈁4]、稀疏分解[5-6]、獨立分量分析[7]、譜峭度[8]和形態(tài)濾波[9]等方法。

        雖然基于振動傳感器的軸承故障診斷技術(shù)取得了廣泛應(yīng)用,但是在高溫、高腐蝕及運動等難于接近的場合下,安置振動傳感器非常困難,此時不適于采用振動傳感器來監(jiān)測軸承的運動狀態(tài)。由于傳聲器可以在較遠距離采集到聲音信息,且聲音是振動產(chǎn)生的,其頻率為振動頻率,因此在這些難于接近的場合下采用傳聲器來檢測軸承的振動是一種不錯的選擇。鑒于傳聲器的非接觸、可實現(xiàn)較遠距離檢測等特點,研究人員已經(jīng)較早開展了聲學故障診斷研究。美國鐵道協(xié)會及所屬研究機構(gòu)開發(fā)出了一套軸承道旁聲學監(jiān)測系統(tǒng),采取傳聲器對列車軸承的運動狀態(tài)進行監(jiān)測[10]。Deblauwe[11]和Hald[12]基于發(fā)動機噪聲完成了制動器噪聲源分析、發(fā)動機的撞擊啟動及撞擊噪聲分析。Lu等[13]將近場聲全息和灰度共生矩陣結(jié)合用于軸承故障診斷等。這些研究都證明了利用聲信號進行故障診斷是可行的。

        利用聲信號對軸承進行狀態(tài)監(jiān)測和故障診斷遇到的最大問題是噪聲干擾。現(xiàn)場采集到的軸承故障聲信號信噪比差,成分復雜,在介質(zhì)中會迅速衰減。針對這一問題,研究者們采用了不同的方法來解決。Routray等[14]用小波方法和盲源分離的方法消除噪聲。Amarnath等[15]用經(jīng)驗?zāi)B(tài)分解(empirical mode decomposition,簡稱EMD)方法處理滾動軸承故障的聲信號,提取振動特征。周俊等[16]將時域盲解卷積、形態(tài)濾波和頻域壓縮感知重構(gòu)的稀疏分量分析相結(jié)合應(yīng)用于軸承復合故障聲學診斷。這些方法對消除軸承聲信號中的噪聲都取得了較好的效果。

        筆者結(jié)合全息技術(shù)與稀疏分解算法的優(yōu)點,提出一種基于稀疏分解的軸承聲陣列信號特征提取方法。首先,利用近場聲全息聲源識別技術(shù)判斷機械設(shè)備的主要噪聲源;然后,根據(jù)軸承故障聲信號的頻率特征和冗余字典的稀疏特性建立由coif4小波字典和局部余弦字典組成的冗余字典來稀疏逼近原信號,并采用拉格朗日基追蹤算法去噪。對于基追蹤的基不穩(wěn)定問題,利用拉格朗日基追蹤方法直接求解稀疏去噪會更有效[17]。與以往方法相比[14-16],本研究方法沒有EMD方法存在的模態(tài)混淆和端點效應(yīng),而文獻[14]中的盲源分離處理實際非線性信號的能力欠佳,且本稀疏分解算法中的基追蹤所挑選的向量為線性無關(guān),比匹配追蹤(貪婪算法)更適用于分解含有相近時頻結(jié)構(gòu)的信號。筆者對模擬的含噪軸承聲信號進行了分析,搭建了軸承聲陣列信號故障診斷實驗平臺。利用兩種傳聲器陣列識別設(shè)備主要噪聲源后再對“熱點”所對應(yīng)信號進行稀疏分解,提取清晰的故障特征頻率,準確識別了軸承故障。

        1 軸承振動的近場聲全息聲源識別原理

        近場聲全息聲源識別技術(shù)的最大特點是由傳聲器陣列信號采集設(shè)備在近場區(qū)域記錄被測對象發(fā)出的聲波,再用分析處理算法對聲源的輻射聲壓直接進行分析成像或重建,其中傳聲器測量面為全息面。在自由與半自由聲場中,由點聲源聲壓級關(guān)系可知,每當測量距離增大一倍,聲壓級就會減少6dB,且當全息面距離聲源大于一個波長時,傳聲器只能采集到聲源的宏觀信息(傳播波)而丟失了聲源的細節(jié)信息(倏逝波,其幅值隨距離的增加而指數(shù)衰減)。根據(jù)波速(v=340m/s)等于波長λ與頻率f的乘積,當軸承故障特征頻率在f=20~500Hz范圍內(nèi)時,聲波長λ=0.68~17m。因此,可以采用傳聲器陣列在近場(測量距離小于聲波波長d?λ)區(qū)域記錄并顯示軸承的全息數(shù)據(jù)。

        軸承元件受迫振動時輻射的聲波隨時間和空間快速改變,利用全息面有效聲壓場可以獲得聲源個數(shù)、位置及強度等信息,并通過對比正常與故障軸承工況下的設(shè)備聲壓圖,識別出主要噪聲源位置,進而判斷關(guān)鍵部位軸承是否異常。

        2 稀疏分解原理

        (1)

        其中:D為信號空間RN中的過完備冗余字典,D=D1,D2,,Dk;αk為對應(yīng)于每個字典的分解表示系數(shù)。

        2.1 拉格朗日基追蹤算法

        在冗余字典中,利用基追蹤算法選取多個字典向量計算信號的最優(yōu)逼近使誤差最小,這是一個非確定性多項式難題(non-deterministic polynomial-hard,簡稱NP-hard問題),所以未必得到最優(yōu)的逼近。式(1)的解可通過極小化l0拉格朗日函數(shù)獲得。l0偽范數(shù)是非凸的,此極小化問題難以求解,且當字典個數(shù)增加時計算將更加復雜。因此,為了逼近信號x,將關(guān)于分解系數(shù)αk的l0范數(shù)優(yōu)化求解問題歸結(jié)為式(2)中ak的l1范數(shù)優(yōu)化問題

        (2)

        式(2)為凸函數(shù),通過其拉格朗日追蹤的松弛公式求解

        (3)

        其中:T為拉格朗日乘子,即閾值。

        對于信號xk,可得

        (4)

        由式(3),(4)可知,式(2)中計算系數(shù){a1,,ak}可轉(zhuǎn)換為計算{x1,,xk}的問題

        (5)

        式(5)中{x1,,xk}的求解算法步驟為:

        1) 選合適的字典Dk構(gòu)造冗余字典D;

        2) 設(shè)定閾值T,初始化迭代次數(shù)n;

        3) 初始化各信號xk=0;

        2.2 冗余字典的構(gòu)建

        圖1為在半消聲實驗室內(nèi)采集的軸承內(nèi)圈(SKF6205-2RSJEM,裂紋尺寸約為0.1 mm×15 mm×0.28mm,轉(zhuǎn)速為1 800r/min)故障聲信號的時域波形圖,其故障沖擊間隔理論值為0.006 2s。圖中的故障聲信號包含沖擊信號、諧振信號和噪聲,故不宜選用單一的小字典稀疏表示,需選擇能稀疏表示各有效信號的字典來構(gòu)建冗余字典(Dk的個數(shù)大于k)。

        圖1 軸承內(nèi)圈故障聲信號的時域波形圖Fig.1 Acoustic signal time-domain waveform of bearing inner ring fault

        常用的冗余字典有Dirac字典、Fourier字典、小波字典和局部余弦字典等。一組最優(yōu)基是一組與信號時頻結(jié)構(gòu)匹配最好的時頻鋪疊。根據(jù)圖1信號特點,為了減少稀疏分解中最佳基逼近的復雜度,可選擇時頻標準正交字典,如小波字典和局部余弦字典。

        小波字典由一個積分為零的母小波Ψ構(gòu)造,通過伸縮函數(shù)的尺度參數(shù)s和平移時域u后可得

        (6)

        由式(6)可知,小波字典將頻率軸分成不同長度的區(qū)間,能為奇異點提供更稀疏的表示,適合逼近圖1中的沖擊信號。

        常用的小波函數(shù)有Haar小波、Daubechies小波、Meyer小波、symN小波族(變量N表示小波函數(shù)的消失矩階數(shù),一般N=2,3,)和coifN小波族等。小波字典與信號的相似性越高,字典的稀疏性越好,因沖擊信號具有快速衰減的特點,需選擇具有指數(shù)性和正交性(便于重構(gòu)),且存在緊支集的小波基函數(shù)。其中,sym8小波和coif4小波為正交小波,具有較長的支集長度和較大的消失矩,且對稱性較好,重構(gòu)時不易失真。sym8小波消失矩比coif4小波的大,信號分解后高頻分量更少,低頻分量更多,使運算量增加,沖擊信號峰值被削弱的更多。因此,coif4小波字典更適宜匹配軸承故障聲信號,其尺度函數(shù)和小波函數(shù)如圖2所示。

        圖2 Coif4小波的尺度函數(shù)和小波函數(shù)Fig.2 Scaling function and wavelet function of coif4

        (7)

        由式(7)可知,局部余弦字典將時間軸分成不同長度的區(qū)間,且Heisenberg時頻盒沿頻率平移,其時頻分辨率與諧波結(jié)構(gòu)相適應(yīng),適合逼近諧振信號。

        綜上所述,通過分析軸承內(nèi)圈故障聲信號時頻特征和冗余字典稀疏特性,選用coif4小波字典和局部余弦字典構(gòu)建冗余字典。

        3 仿真分析

        利用模擬的聲信號驗證本研究方法可以在復雜的噪聲背景下有效提取軸承聲信號的故障特征。以內(nèi)圈故障聲信號為例,其模型為

        x(t)=r(t)+s(t)+h(t)

        (8)

        其中:x(t)為復合聲信號;r(t)為噪聲信號;s(t)為沖擊信號;h(t)為諧振信號。

        u(t-mT)exp-β(t-mT)

        (9)

        其中:s(t)由M個幅值為B、衰減系數(shù)為β、高頻共振頻率為fn的沖擊成分組成;T為沖擊間隔;f=1/T為沖擊故障頻率;u(t)為單位階躍函數(shù);fr為轉(zhuǎn)頻。

        (10)

        其中:h(t)由K個幅值為A、頻率為fi的諧波構(gòu)成。

        為了使仿真更接近實際,取采樣頻率fs=10kHz,采樣點數(shù)N=99 000,加入高斯白噪聲, 使信噪比SNR=-15,并將表1中的參數(shù)代入式(8)~(10),得到仿真合成聲信號時域圖及Hilbert包絡(luò)頻譜圖,如圖3所示??梢钥闯觯瑳_擊特征已被強噪聲完全淹沒,普通的Hilbert包絡(luò)解調(diào)方法無法提取仿真信號的沖擊故障特征。

        利用本研究方法對圖3所示的仿真合成信號進行稀疏分解并包絡(luò)解調(diào),提取到的沖擊信號及其Hilbert頻譜圖如圖4所示??梢钥闯?,噪聲已基本被分離,時域圖中有明顯的周期性沖擊特征,其Hilbert頻譜圖中轉(zhuǎn)頻31.3Hz及二倍轉(zhuǎn)頻60.1Hz、內(nèi)圈故障特征頻率164.1Hz及其2~5次諧波和邊頻處峰值(接近理論值)比較凸顯。仿真實例證明,本研究方法可在復雜的噪聲背景下準確提取故障特征。

        表1 仿真信號各參數(shù)值

        圖3 仿真合成聲信號時域圖及Hilbert頻譜圖Fig.3 Time-domain diagram and Hilbert spectrum of acoustic synthetic signal

        圖4 稀疏分解提取到的沖擊信號及Hilbert頻譜圖Fig.4 Impact signal extracted by sparse decomposition and its Hilbert spectrum

        4 實例分析

        4.1 實驗平臺

        為了驗證基于稀疏分解的軸承聲陣列信號特征提取方法的有效性,在實驗室內(nèi)搭建軸承聲陣列信號故障診斷實驗平臺。平臺四周堆疊高為2m的吸聲尖劈形成簡易半消聲室,實驗數(shù)據(jù)在凌晨0~6點采集,環(huán)境噪聲較小。如圖5所示,該平臺主要由軸承實驗臺、傳聲器陣列以及聲頻信號采集處理系統(tǒng)組成。圖5(a)所示,工作原理為:a.利用數(shù)據(jù)采集箱(PXI1033)和采集卡(Express Card-8306)采集MP201電容型傳聲器(直徑為0.012 7m,頻率響應(yīng)范圍為6.3~20kHz)陣列測得的本特利RK4軸承實驗臺聲信號數(shù)據(jù);b.對比分析傳聲器陣列全息面三維有效聲壓場圖,根據(jù)“熱點”判斷關(guān)鍵部位軸承工作狀態(tài)是否異常;c.利用上述特征提取方法對“熱點”聲信號進行故障特征提取,進而判斷故障類型。

        圖5(c)和(d)中,傳聲器陣列架為7×7平面網(wǎng)格(0.1m×0.1m)結(jié)構(gòu),置于軸承實驗臺正前方(陣列傳聲器端面與轉(zhuǎn)軸軸心線的垂直距離為18.2cm)。考慮到軸承實驗臺為長方形結(jié)構(gòu),可適當減少縱向傳聲器個數(shù)。因此,實驗擬采用6×6和4×6兩種傳聲器陣列。其中,6×6和4×6陣列中的第4和第19通道傳聲器正對軸承實驗臺的電機,而第22和第7通道正對實驗軸承所在軸承座。

        圖5 軸承聲陣列信號故障診斷實驗平臺Fig.5 Microphone array signal fault diagnosis experiment platform of bearing

        軸承處于空載狀態(tài),轉(zhuǎn)速分別取1 800r/min和1 500r/min。圖6為線切割加工的外圈和內(nèi)圈裂紋故障。根據(jù)表2中的軸承結(jié)構(gòu)參數(shù),計算出外圈和內(nèi)圈的故障頻率分別為3.584 8fr和5.415 2fr(fr為轉(zhuǎn)頻)。當軸承轉(zhuǎn)速為1 800r/min(fr=30Hz)時,外圈和內(nèi)圈的故障頻率理論值分別為107.54Hz和162.46Hz。當轉(zhuǎn)速為1 500r/min(fr=25Hz)時,外圈和內(nèi)圈的故障頻率理論值分別為89.62和135.38Hz。

        圖6 故障軸承Fig.6 Failure bearings

        參數(shù)數(shù)值外徑/mm52內(nèi)徑/mm25寬度/mm15滾動體接觸角α/(°)0外圈裂紋寬度1/mm0.5內(nèi)圈裂紋寬度1/mm0.5內(nèi)圈裂紋寬度2/mm0.1裂紋深度/mm0.28裂紋位置滾道寬度方向

        4.2 實驗結(jié)果及分析

        實驗中,信號采樣頻率為10kHz,每個通道數(shù)據(jù)長度為98 304。圖7為傳聲器陣列采集的正常、內(nèi)圈和外圈裂紋(0.1或0.5mm)故障軸承在不同轉(zhuǎn)速下的全息面有效聲壓場及投影圖。其中,圖7(a)~(e)采用6×6傳聲器陣列,圖7(f)和(g)采用4×6傳聲器陣列。經(jīng)對比分析可知

        1) 圖7(a)中,只有電機附近的坐標(-0.2,0)有“熱點”,且正常軸承有效聲壓低于故障軸承。

        2) 圖7(b)和(c)中的“熱點”均出現(xiàn)在坐標(0,0)附近,即靠近軸承實驗臺中心位置,且實驗軸承端聲壓較高。這是因為外圈裂紋引起的異常振動通過軸承座和軸傳遞,其產(chǎn)生的聲音與電機噪聲疊加而形成此熱點。

        3) 圖7(d)中存在兩個“熱點”,一個在坐標(0,-0.1)處,其分布和成因與圖7(b)的基本一致;另一個在坐標(0.1,0.1)處,即對應(yīng)故障軸承,其聲壓略低于第1個“熱點”。原因是電機和故障軸承均為噪聲源,而第2個“熱點”主要受故障軸承的影響。類似的,圖(e)中也存在兩個“熱點”,其性質(zhì)與圖7(d)中的相同。

        圖7 不同工況下全息面有效聲壓場與投影圖Fig.7 Effective sound pressure fields and projections of holographic planes under different working conditions

        4) 圖7(f)和(g)中,當采用4×6陣列傳聲器測量內(nèi)圈裂紋(0.1mm)故障軸承聲場時,在坐標(0.25, 0.15)和(0.1, 0.1)處出現(xiàn)兩個“熱點”,與圖7(d)和(e)的熱點位置基本對應(yīng)。

        “熱點”意味著主要噪聲源位置,據(jù)此可初步判斷關(guān)鍵部位軸承是否異常,為進一步識別軸承故障類型提供依據(jù)。

        按噪聲功率信噪比SNR=-15的方式在信號整個時域范圍內(nèi)添加高斯白噪聲,模擬外界干擾以形成信噪比很差的實際工況。利用Hilbert頻譜分析和筆者提出的方法處理故障軸承不同工況下所對應(yīng)的“熱點”聲信號,并將稀疏分解前后的時域圖與Hilbert頻譜圖進行對比分析。

        圖8和9分別為外圈裂紋(0.5mm)故障軸承在轉(zhuǎn)速為1 800和1 500r/min時的聲信號處理結(jié)果。圖8(a)為半消聲室內(nèi)測得轉(zhuǎn)速為1 800r/min時外圈軸承故障聲信號的時域圖及Hilbert頻譜圖,可以看到明顯的周期性沖擊,其Hilbert頻譜圖呈現(xiàn)出清晰的外圈故障特征頻率107.1Hz及2~4次諧波。對兩種轉(zhuǎn)速下的加噪信號進行處理,稀疏分解前后的時域圖與Hilbert頻譜圖如圖8(b),(c)和9(a),(b)所示。圖8(b)和9(a)時域圖中噪聲較大,完全掩蓋了故障特征,因受到雜頻和噪聲干擾,即使從其Hilbert頻譜圖也只能依稀識別出外圈故障特征頻率(分別為107.1和89.19Hz,接近相應(yīng)的理論值)。圖8(c)和9(b)為稀疏分解處理得到的沖擊信號??梢钥吹?,在信噪比很差的情況下,稀疏分解提取到了明顯的周期性沖擊,噪聲已基本被去除,其Hilbert頻譜圖呈現(xiàn)出清晰的外圈故障特征頻率(分別為107.1和89.21Hz)及2~5次諧波,故障頻率峰值也更加凸顯。

        圖8 外圈0.5mm裂紋故障分析(1 800r/min)Fig.8 Analysis of 0.5mm outrace crack fault (1 800r/min)

        圖9 外圈為0.5mm裂紋故障分析(1 500r/min)Fig.9 Analysis of 0.5mm outrace crack fault (1 500r/min)

        圖10 內(nèi)圈為0.5mm裂紋故障分析(1 800r/min)Fig.10 Analysis of 0.5mm inner race crack fault (1 800r/min)

        圖11 內(nèi)圈為0.5mm裂紋故障分析(1 500r/min)Fig.11 Analysis of 0.5mm inner race crack fault (1 500r/min)

        圖12 內(nèi)圈為0.1mm裂紋故障分析(1 800r/min)Fig.12 Analysis of 0.1mm inner race crack fault (1 800r/min)

        圖13 內(nèi)圈為0.1mm裂紋故障分析(1 500r/min)Fig.13 Analysis of 0.1mm inner race crack fault (1 500r/min)

        圖10和11分別為內(nèi)圈裂紋(0.5mm)故障軸承在轉(zhuǎn)速為1 800和1 500r/min時的加噪聲信號處理結(jié)果。由圖10(a)和11(a)可知,從時域圖和Hilbert頻譜圖均無法判斷軸承故障類型。從圖10(b)和11(b)稀疏分解處理得到的沖擊信號時域圖中可看到明顯的周期性沖擊,噪聲已基本被去除,并且其Hilbert頻譜圖中相應(yīng)的二倍轉(zhuǎn)頻59.81和49.83Hz以及故障特征頻率161.7和134.8Hz的峰值也更加凸顯。

        圖12和13為內(nèi)圈裂紋(0.1mm)故障軸承在轉(zhuǎn)速分別為1 800和1 500r/min時的加噪聲信號處理結(jié)果。由圖12(a)和13(a)可知,時域圖和Hilbert頻譜圖仍然無法準確判斷軸承故障類型。從圖12(b)和13(b)稀疏分解處理得到的沖擊信號時域圖中可看到明顯的周期性沖擊,噪聲已基本被去除,并且其Hilbert頻譜圖中相應(yīng)的二倍轉(zhuǎn)頻59.92和49.74Hz以及故障特征頻率161.9和134.5Hz的峰值也更加凸顯。

        實例分析可見,在不同的故障類型和轉(zhuǎn)速下,傳統(tǒng)的Hilbert包絡(luò)譜分析只能勉強甚至完全不能提取出故障特征,而稀疏分解方法在強噪聲干擾下有效提取了軸承聲陣列信號的故障特征,識別出了相應(yīng)的故障,驗證了該方法的可行性。

        5 結(jié)束語

        利用傳聲器陣列對軸承滾道缺陷故障進行檢測,結(jié)合全息技術(shù)與稀疏分解算法的優(yōu)點,提出一種基于稀疏分解的軸承聲陣列信號特征提取方法。利用仿真分析證明了此方法在強噪聲背景下提取軸承故障聲信號特征的有效性。

        在半消聲室內(nèi)搭建了軸承聲陣列信號故障診斷實驗平臺,利用正常、內(nèi)外圈故障軸承進行實驗研究,對比分析了不同工況下得到的全息面三維有效聲壓場、提取的沖擊信號時域圖和Hilbert頻譜圖。結(jié)果表明,與傳統(tǒng)的Hilbert包絡(luò)譜分析方法相比,本研究方法消除了聲信號中的噪聲干擾,提取到了清晰的沖擊特征,其Hilbert頻譜圖中故障特征頻率峰清晰,可準確識別軸承故障類型。

        猜你喜歡
        裂紋故障信號
        裂紋長度對焊接接頭裂紋擴展驅(qū)動力的影響
        信號
        鴨綠江(2021年35期)2021-04-19 12:24:18
        完形填空二則
        故障一點通
        Epidermal growth factor receptor rs17337023 polymorphism in hypertensive gestational diabetic women: A pilot study
        基于FPGA的多功能信號發(fā)生器的設(shè)計
        電子制作(2018年11期)2018-08-04 03:25:42
        奔馳R320車ABS、ESP故障燈異常點亮
        基于LabVIEW的力加載信號采集與PID控制
        故障一點通
        江淮車故障3例
        国产乱理伦片在线观看| 精品女同一区二区三区| 日韩夜夜高潮夜夜爽无码| 男人添女人下部高潮全视频| 国产精品嫩草影院午夜| 亚洲精品久久久中文字| 亚洲中文字幕日韩综合| 亚洲av无码国产精品永久一区| 99久久免费精品高清特色大片| 日韩精人妻无码一区二区三区| 亚洲国产一区二区av| 中文字幕无线码一区二区| 久久无码人妻精品一区二区三区 | 欧美日韩性视频| 在线日韩中文字幕乱码视频| 国产一区二区三区三区四区精品| 亚洲av无码专区首页| 亚洲色成人网一二三区| 日韩亚洲午夜精品一区二区三区| 日本女优在线一区二区三区| 日韩毛片免费无码无毒视频观看| 精品人妻无码中文字幕在线| 日韩一区二区中文字幕视频| 97久久婷婷五月综合色d啪蜜芽| 日本熟妇人妻xxxxx视频| 久久久久久久尹人综合网亚洲 | 久久精品国产亚洲av四虎| 国产精彩视频| 亚洲一区二区日韩精品| 欧美精品国产综合久久| 草草网站影院白丝内射| 极品美女销魂一区二区三| 日韩精品熟女中文字幕| 日韩成人大屁股内射喷水| 午夜国产精品久久久久| 国产一区二区三免费视频| 无码丰满熟妇一区二区| 久久av无码精品人妻糸列| 视频在线亚洲视频在线| 奇米影视7777久久精品| 国产3p视频|