999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

基于參考信號頻域半盲提取的機械故障特征聲學診斷

2015-04-29 00:44:03羿澤光潘楠劉鳳
河北科技大學學報 2015年4期

羿澤光 潘楠 劉鳳

摘要:針對生產現場機械設備零部件結構復雜、設備運行時背景噪聲干擾嚴重等造成的監測診斷難題,以及傳統盲信號處理算法在機械聲信號處理方面的局限性,提出一種基于參考信號約束頻域半盲提取的機械故障聲學診斷算法。詳細介紹了該算法的關鍵技術:以頻域盲解卷積算法為基礎,使用利于全局尋優的人工魚群算法,構建適用于機械故障特征的改進多尺度形態學濾波器,以最大程度削弱背景噪聲干擾;結合機械設備零部件結構參數構建參考信號,通過單元參考信號約束頻域半盲提取算法,對降噪后的信號逐段進行復數盲分離;利用改進KL距離,解決復分量間次序不確定性問題,最終實現機械故障特征信號的提取與分離。實際聲場環境中的滾動軸承故障聲學診斷實驗驗證了該算法的有效性。

關鍵詞:算法理論;參考信號約束;頻域半盲提取;人工魚群算法;聲學診斷

中圖分類號:TH1653文獻標志碼:A

收稿日期:2014-11-02;修回日期:2015-01-01;責任編輯:張士瑩

基金項目:國家自然科學基金(51305186);昆明理工大學校級人培基金(KK3201301026)

作者簡介:羿澤光(1989—),男,河北唐山人,碩士研究生,主要從事狀態監測與故障診斷方面的研究。

通訊作者:潘楠博士。E-mail:15808867407@163.com

羿澤光,潘楠,劉鳳.基于參考信號頻域半盲提取的機械故障特征聲學診斷[J].河北科技大學學報,2015,36(4):351-358.

YI Zeguang,PAN Nan,LIU Feng. Acoustic diagnosis of mechanical fault feature based on reference signal frequency domain semi-blind extraction[J].Journal of Hebei University of Science and Technology,2015,36(4):351-358.Acoustic diagnosis of mechanical fault feature based on reference signal

frequency domain semi-blind extraction

YI Zeguang, PAN Nan, LIU Feng

(Faculty of Mechanical & Electrical Engineering, Kunming University of Science and Technology, Kunming, Yunnan 650500, China)

Abstract:Aiming at fault diagnosis problems caused by complex machinery parts, serious background noises and the application limitations of traditional blind signal processing algorithm to the mechanical acoustic signal processing, a failure acoustic diagnosis based on reference signal frequency domain semi-blind extraction is proposed. Key technologies are introduced: Based on frequency-domain blind deconvolution algorithm, the artificial fish swarm algorithm which is good for global optimization is used to construct improved multi-scale morphological filters which is applicable to mechanical failure in order to weaken the background noises; combining the structural parameters of parts to build a reference signal, complex components blind separation is carried out on the signals after noise reduction paragraph by paragraph by reference signal unit semi-blind extraction algorithm; then the improved KL-distance of complex independent components is employed as distance measure to resolve the permutation, and finally the mechanical fault characteristic signals are extracted and separated. The actual acoustic diagnosis of rolling bearing fault in sound field environment results proves the effectiveness of this algorithm.

Keywords:algorithm theory; reference signal constraints; frequency-domain semi-blind extraction; artificial fish swarm algorithm; acoustic diagnosis

機械設備運轉時,其振動信號中蘊含著很多狀態信息,尤其是當設備出現故障時,振動特征信號內會產生明顯的沖擊成分,由此衍生的聲信號特性將隨之發生改變[1]。機械設備故障診斷能否順利進行很大程度上依賴于能否從繁多復雜的機械狀態信號中提取足夠數量且能夠真實而客觀地反映診斷對象工況的信息[2]。聲信號測試具有非接觸測量、測試方式簡便、在線測試和無附著物影響等諸多優點,尤其適用于工業現場狀態監測與故障診斷[3]。然而實際聲場環境復雜多變,待識別信號(故障源信號)常與各種干擾信號或噪聲相互混雜,無法被有效辨識[4]。因此,為了準確提取機械故障特征,需要將這些噪聲抑制或排除,隨后根據故障征兆識別故障原因[5]。

河北科技大學學報2015年第4期羿澤光,等:基于參考信號頻域半盲提取的機械故障特征聲學診斷 近年來,能夠在幾乎沒有任何先驗知識的情況下,從混合信號中恢復或估計出源信號的盲信號處理技術為機械故障信號的提取提供了一個有力的解決手段[6-7]。然而,傳統的盲分離算法在應用到機械聲信號處理時往往無法滿足現實情況,不能有效進行機械故障特征的識別與提取,而盲解卷積算法則更適用于實際工業聲場環境[8-11]。因此,本文嘗試以頻域盲解卷積算法為基礎,基于全局尋優能力較強的人工魚群算法,構建針對故障信號特征的多尺度形態學濾波器,最大程度削弱背景噪聲干擾;結合機械零部件結構參數構建參考信號,進而通過單元參考信號約束半盲提取算法,逐段進行復分量盲分離;利用改進KL距離解決復分量間次序不確定性問題,最終實現實際聲場環境中的故障特征提取。

1頻域盲解卷積基本模型

基于文獻\[3\]的論述,可將頻域盲解卷積算法的基本步驟整理如下。

通過加窗短時傅里葉變換(STFT)將拾取到的時域觀測信號轉換到頻域中;此時,時域卷積混合就被轉換成為對應頻域上各頻段的瞬時混合,從而引入復數盲分離算法對源信號進行估計。針對子信號輸出次序不確定性問題,將估計的復數信號分頻段進行重排序,通過加窗逆短時傅里葉變換(ISTFT)將信號轉換回時域,最終得到估計信號。

機械系統產生的信號多種多樣,加之設備所處環境中的背景噪聲干擾,使得頻域盲解卷積算法應用于實際復雜聲場中針對機械故障特征的提取時存在很多問題[12-16]:

1)機械聲場較為復雜且具有多源干擾,故障特征聲易淹沒在較強的高斯噪聲、復雜周期信號及其他非平穩信號之中,無法被有效識別;

2)聲學傳感器(傳聲器)不同于振動傳感器,其安裝位置與故障源距離較遠,聲信號在傳輸過程中容易因長卷積而發生衰減,以致算法求解的復雜度不斷增加;

3)使用頻域盲解卷積算法提取故障特征信號時,循環部分卷積誤差和次序不確定性等固有問題均會影響分離性能;

4)傳統頻域盲信號處理算法往往無法直接適用于機械故障聲振信號的提取。

2參考信號約束頻域半盲提取故障聲信號

針對前文所述問題,通過對原有頻域盲解卷積算法的各項步驟進行改造優化,降低背景噪聲及凸顯特征頻率范圍,將其應用到現實聲場故障特征提取之中。

2.1構造改進多尺度形態濾波器

2.1.1形態濾波簡介

近年來,能夠提取信號細節和抑制干擾噪聲形態濾波也被逐漸應用于機械聲振信號背景降噪[17]。其濾波器結構都是基于差值濾波器或者形態開-閉(OC)和閉-開(CO)平均組合,且多為單一結構元素[18]。然而,實際工業聲場中往往不止存在某一種干擾噪聲,且信號中的噪聲通常是隨機的,需要使用多尺度的不同結構元素構建形態濾波器,以避免較嚴重的濾波器輸出統計偏移。

2.1.2人工魚群算法優化結構元素

形態濾波的結構元素構成包括形狀、幅值和尺寸等多個方面。實驗結果表明,半圓形結構元素能夠較好地濾除隨機噪聲,三角形結構元素對于脈沖噪聲的濾波效果較好。為了避免目標尋優過早地進入局部收斂,特采用基于集群體智能思想的人工魚群算法進行全局尋優,將經過均值化處理后的觀測信號中相鄰峰值間隔的最大值和最小值通過極值尋優,確定結構元素的長度;隨后,根據信號峰值的最大值和最小值確定高度范圍;最后,將相應的結構元素尺寸分別代入半圓形和三角形結構公式,計算出對應的結構元素集合。整個算法流程如圖1所示。

圖1人工魚群算法尋優流程

Fig.1Flow of improved multi-scale morphological filtering algorithm

圖1中:N為人工魚群大小;{Xi}為人工魚個體的狀態位置;Yi=f(Xi)為第i條人工魚當前所在位置的實物濃度;Visual為人工魚的感知距離;Step為人工魚移動的最大步長;delta為擁擠度;n為當前覓食行為次數;gen為迭代次數;maxgen為最大迭代次數。

2.1.3改進多尺度形態濾波

閉運算可以較完整地提取信號中的正脈沖,平均組合形態濾波器則可以同時抑制信號中的正負脈沖噪聲;但由于在形態聯級中使用相同尺寸的結構元素,容易導致濾波器輸出統計偏移較為嚴重[17],因此需構建多尺寸、多結構聯級而成的結構元素閉-開組合形態濾波器。式(1)中:和Θ分別表示膨脹和腐蝕運算。

y1(n)=(fgΘgΘg)(n),

y2(n)=(fgΘgΘgg)(n),

y(n)=[y1(n)+y2(n)]/2。 (1)

整體改進多尺度形態濾波算法流程如下:

1) 對觀測信號x(t)進行均值化處理,同時構建初始化多尺度元λg;

2) 計算觀測信號的局部極大值和極小值,確定結構元素的高度和長度集合;

3) 將計算的高度和長度代入三角形和半圓形結構元素公式,構造結構元素集合;

4) 代入式(1)得到y(n)組合濾波器集合;

5) 采用人工魚群算法對濾波器集合進行尋優,使用y(n)對觀測信號x(t)進行濾波處理,直至得到最終去噪信號。具體流程見圖2。

圖2改進多尺度形態濾波算法流程

Fig.2Flow of improved multi-scale morphological

filtering algorithm

2.2參考信號約束半盲提取算法

為了減少算法復雜度,常針對機械設備故障診斷過程中重點關注的包含零部件故障特征頻率的低頻包絡周期信號,通過匹配關鍵零部件結構參數,基于先驗知識構造出理論上的參考信號。以此參考信號為約束,將觀測信號與構造參考信號進行相似性度量,進而提取有限數量的估計信號,此即稱為半盲提取。僅利用單個參考信號提取估計信號的情況稱之為單元參考信號約束半盲提取,即轉換為式(2)的約束問題[20]:

J(y)≈ρ[E{G(y)}-E{G(r)}]2,

g(w)=ε(y,r)-ξ≤0。

約束條件:

g(w)≤0,h(w)=E{y2}-1=0。 (2)

式中:J(y)為ICA算法的負熵目標函數;ε(y,r)表示觀測信號與參考信號的相似性測度;ξ為區分估計信號與其他信號的閾值,其選取為漸變適應過程,在保證約束目標函數有解的前提下同時避免陷入局部最優。

將式(2)寫為無約束拉格朗日函數:

L(w,μ)=J(y)-12γmax2[μ+rg(w),0]-μ2-λh(w)-12γ‖h(w)‖2。(3)

式中:μ和λ為拉格朗日乘子;γ為尺度懲罰參數。通過對式(3)的極值進行求解,即可得到牛頓迭代算法:

wk+1=wk-η R-1xxL″wk/δ(wk),(4)

L′wk=ρE{xG′y(y)}-0.5μE{xg″y(wk)}-λE{xy},

δ(wk)=ρE{xG″y2(y)}-0.5μE{xg″y2(wk)}-λ。(5)

式中:η為學習率;Rxx為混合信號的協方差矩陣;G′y,G″y2和g″y2(wk)分別表示Gy與gy(wk)的一階導數和二階導數;μ和λ由式(6)求得:

μk+1=max{0,μk+γg(wk)},

λk+1=λk+λh(w)。 (6)

當信號經過白化與中心化處理之后,R-1xx=1,則式(4)可寫為如下形式:

wk+1=wk-ηL′wk/δ(wk)。(7)

2.3解決次序不確定性問題

距離是衡量2個信號相似性的測度,距離越大,說明2個信號相似性越小;距離越小,則說明2個信號相似度越大。距離互參數法的基本原理是通過比較距離,將每個頻率段的輸出調整至同一通道對應著同一源信號。

高階統計量信息較二階統計量如相關系數等更能清晰地反映信號的統計特性,改進KL距離可被用以更好地描述相鄰2個頻率段內2個復值信號對應概率密度之間的距離:

KL(pj(ω),pi(ω+1))=∑Mm=1pj(ω,m)·logpj(ω,m)pi(ω+1,m)。 (8)

式中:pj(ω)=[pj(ω,1),pj(ω,2),…,pj(ω,M)];pj(ω,m)=pj,fre(ω,m)·|Yuj(ω,m)|2∑Mr=1|Yuj(ω,r)|2;pj,fre(ω,m)為Yuj(ω,m)在[Yuj(ω,1),Yuj(ω,2),…,Yuj(ω,M)]中出現的概率,其能更為準確地反映信號的復數信息。

圖3相似度聚類效果對比

Fig.3Comparison of similarity clustering effects

為了驗證改進KL距離在復分量相似測度計算上的優勢,特進行如下仿真:分別利用改進KL距離及常被用于計算譜聚類的余弦測度(相似度),對18組復分量進行相似度聚類計算。從聚類散點圖(見圖3)可以看出,利用改進KL距離求解的距離聚類散點線性程度較佳,且其運算速度較基于余弦測度的算法快近50%。

2.4算法整體步驟

綜合前文所述,現將參考信號約束頻域半盲提取算法整體執行步驟歸納如下:

1) 初始化零部件結構參數,根據結構參數計算特征頻率,構造參考信號r(t),隨后對觀測信號x(t)進行中心化處理;

2) 利用改進多尺度形態濾波器對x(t)進行噪聲抑制,得到濾波后的信號(t);

3) 使用STFT將信號(t)及r(t)轉換到頻域中,得到其在頻域的表現形式X(ω,t),R(ω,t);

4) 利用單元參考信號約束半盲提取算法,使用復分量間改進KL距離進行相似度測量,對X(ω,t)進行逐段復數盲提取,得到與參考信號R(ω,t)相似度最高的各頻段復值估計信號Y(ω,t);

5) 通過ISTFT將信號轉換回時域,最終得到時域估計信號Y(t),求解估計信號的包絡解調譜,對特征頻率進行分析,進行故障判斷。

3實驗驗證

為驗證本文提出算法的實際效能,于現實聲場環境中構建實驗裝置,以軸承故障特征信號為約束進行半盲提取實驗。本次實驗使用故障類型未知的NU205滾動軸承作為診斷對象,其相關物理參數如表1所示。

表1故障滾動軸承參數

Tab.1Parameters of the rolling bearing with fault

型號節圓直徑/mm滾動體直徑/mm滾動體個數接觸角度/(°)NU 205397.5120

將軸承外圈固定,內圈隨軸旋轉。使用NI CRIO 9082及相關控制模塊,配合NI-9234四通道采集模塊和3個聲望MPA416 1/4英寸TEDS傳聲器(頻率響應范圍為20 Hz~20 kHz)組成嵌入式采集裝置。利用LabVIEW開發平臺編寫測試程序進行信號采集,由于軸承故障頻率多集中在低頻段,故采樣頻率設置為8 192 Hz,算法由Matlab進行編寫。工況及軸承故障特征頻率如表2所示。

表2實驗工況及軸承故障特征頻率

Tab.2Operation condition and bearing failure frequencies in experiment

轉速/(r·min-1)旋轉頻率/Hz外圈特征頻率/Hz內圈特征頻率/Hz80013.33 64.6195.38

圖4試驗臺及傳聲器位置關系圖

Fig.4Position chart of test-rig and microphones

所有傳聲器均距地面高度1 m擺放并指向軸承座方向,其中傳聲器1與傳聲器2互呈90°,傳聲器3垂直試驗臺擺放。傳聲器探頭距離試驗臺邊緣的直線距離分別為0.64,1.81,1.5 m,并未采用“抵近測量”原則及吸音板等輔助材料,以期盡量貼近現實工況環境,試驗臺及傳聲器的位置關系如圖4所示。

開啟試驗臺驅動電機后觸發采集程序進行聲信號數據采集,采集完畢后通過Matlab從穩態運行情況下的原始數據中截取8 912點數據(1 s)進行分析。圖5為傳聲器拾取到3路信號的時域波形及幅值譜。由圖5可見,3路信號互相混雜在一起,唯有電機轉頻信號譜線(13 Hz)隱約可見,故障沖擊信號已被干擾噪聲完全淹沒。

圖5觀測信號

Fig.5Measurement signals

與本文提出算法的提取效果進行對比,首先使用文獻\[11\]提出的基于FCM模糊聚類的時域盲解卷積算法對觀測信號進行處理,步長為2,聚類個數設置為3,算法提取結果如圖6所示。經該算法運算,得到3個估計信號,但從圖6 a)中未能發現明顯沖擊成分,在對應包絡譜(圖6 b))中第2個分離信號雖可以發現13 Hz處譜線,但其他故障特征頻率并沒有得以體現。隨后利用本文提出的參考信號約束頻域半盲提取算法對圖6基于FCM模糊聚類的時域盲解卷積算法分離結果

Fig.6Signals separated by the time-domain blind deconvolution algorithm based on FCM

觀測聲信號進行處理,算法步驟如2.4節所述,其中STFT分幀長度設置為512,加漢寧窗,窗長512,窗口移動長度64。分別選用三角形和半圓形結構元素,利用人工魚群算法構建多尺度形態濾波,對觀測信號進行背景降噪。算法提取結果如圖7所示。圖7 a)為時域波形,利用光標對圖7 b)第1個分離信號的包絡譜數據進行讀數,可明顯發現65,130,195 Hz 3條譜線(分辨率f=1 Hz),通過對比表2中軸承故障特征頻率可知,其與外圈特征大致吻合。由此可判斷故障所在,拆解軸承后發現軸承外圈出現裂紋故障(如圖8所示),與聲學分析診斷結果相吻合。

圖7參考信號約束頻域半盲提取算法分離結果

Fig.7Signals separated by reference signal frequency domain semi-blind extraction method

圖8軸承外圈裂紋故障示意圖

Fig.8Schematic diagram of bearing outer crack fault4結語

針對目前機械聲信號盲提取所面臨的多種問題,提出了一種基于參考信號約束頻域半盲提取的機械故障特征聲學診斷算法,詳細介紹了人工魚群算法、參考信號構建、半盲提取算法、改進KL距離等關鍵技術的原理與應用,滾動軸承故障聲信號半盲提取實驗證明了該算法的有效性和實用性。

參考文獻/References:

[1]王華民,陳霞,安鋼,等.基于高階累積量的齒輪箱故障診斷研究[J].機械強度,2004,26(3):247-249.

WANG Huamin,CHEN Xia,AN Gang,et al.Fault dlagnosis of gearbox based on higher order cumulants [J].Journal of Mechanical Strength,2004,26(3):247-249.

[2]何正嘉,陳進,王太勇,等.機械故障診斷理論及應用[M].北京: 高等教育出版社,2010.

HE Zhengjia,CHEN Jin,WANG Taiyong,et al.Mechanical Fault Diagnosis Theory and Application [M].Beijing:Higher Education Press,2010.

[3]潘楠,伍星,遲毅林,等.基于頻域盲解卷積的機械設備狀態檢測與故障診斷[J].振動與沖擊,2012,31(12):34-41.

PAN Nan,WU Xing,CHI Yilin,et al.Mechanical equipment status detection and fault diagnosis based on frequency domain blind convolution[J].Journal of Vibration and Shock,2012,31(12):34-41.

[4]潘楠,伍星,遲毅林,等.欠定盲解卷積用于滾動軸承復合故障聲學診斷[J].振動、測試與診斷,2013,33(2):284-289.

PAN Nan,WU Xing,CHI Yilin,et al.Underdetermined blind convolution is used for rolling bearings composite fault acoustic diagnosis[J].Journal of Vibration, Measurement & Diagnosis,2013,33(2): 284-289.

[5]李豫川,伍星,遲毅林,等.基于形態濾波和稀疏分量分析的滾動軸承故障盲分離[J].振動與沖擊,2011,30(12):170-174.

LI Yuchuan,WU Xing,CHI Yilin,et al.Based on the morphological filter and sparse component analysis of the rolling bearing fault separation [J].Journal of Vibration and Shock,2011,30(12):170-174.

[6]NIKOLAOU N G, ANTONIADIS I A.Application of morphological operators as envelope extractors for impulsive-type periodic signals [J]. Mechanical Systems and Signal Processing, 2003, 17(6): 1147-1162.

[7] GELLE G, COLAS M, DELAUNAY G. Blind source separation applied to rotating machin- es monitoring by acoustical and vibrations analysis[J]. Mechanical Systems and Signal Processing, 2000, 14(3):427- 442.

[8]ROUTRAY A D, DASH N. Robust preprocessing: Denoising and whitening in the context of blind source separation of instantaneous mixtures[C]∥2007 5th IEEE International Conference.[S.l.]:[s.n.],2007: 377-380.

[9]THEODOR D P. Blind separation of vibration signals and source change detection:Application to machine monitoring[J]. Applied Mathematical Modelling,2010,34:3408-3421.

[10]KOLDOVSKY Z, TICHAVSKY P. Blind instantaneous noisy mixture separation with best interference-plus-noise rejection[J]. Lecture Notes in Computer Science,2007,4666:730-737.

[11]王宇, 伍星, 遲毅林,等.基于盲解卷積和聚類的機械弱沖擊聲信號提取[J]. 振動工程學報, 2009,22(6):620-624.

WANG Yu,WU Xing,CHI Yilin,et al.Weak transient impulse signal extraction based on blind deconvolution and cluster in acoustical machine diagnosis [J].Journal of Vibration Engineering, 2009, 22(6):620-624.

[12]PAN N, WU X, CHI Y L, et al. Application of frequency-domain blind deconvolution in mechanical fault detection[J]. Applied Mecha-nics and Materials, 2011,134: 2128-2132.

[13]DOUGLAS S C. Fixed-point fast ICA algorithms for the blind separation of complex-valued signal mixtures[C]//Conference Record of the Thirty-Ninth Asilomar Conference on Signals, Systems and Computers.[S.l.]:[s.n.], 2005: 1320-1325.

[14]LAN X, ZHANG Y G, KASSAM S A. Robust source separation using ranks[C]∥Proceedings of the Tenth IEEE Workshop on Statistical Signal and Array Processing.[S.l.]:[s.n.], 2000:324-328.

[15]REJU V G. Partial separation method for solving permutation problem in frequency domain blind source separation of speech signals[J]. Neurocomputing, 2000, 71(10/11/12): 2098-2112.

[16]MAZUR R, MERTINS A. An approach for solving the permutation problem of convolutive blind source separation based on statistical signal models[J]. IEEE Transactions on Audio, Speech and Language Processing, 2009, 17(1): 734-739.

[17]沈路,周曉軍,張文斌,等.廣義數學形態濾波器的旋轉機械振動信號降噪[J].振動與沖擊,2009,28(9):70-73.

SHEN Lu,ZHOU Xiaojun,ZHANG Wenbin,et al.Denoising for vibration signals of a rotating machinery based on generalized mathematical morphological filter[J].Journal of Vibration and Shock,2009,28(9):70-73.

[18]沈路. 數學形態學在機械故障診斷中的應用[D]. 杭州:浙江大學,2010.

SHEN Lu, Study on Application of Morphology in Machinery Fault Diagnosis[D].Hangzhou: Zhejiang University, 2010.

[19]潘楠,羿澤光.基于頻域盲信號處理的汽車制動異響定位方法研究[J].河北科技大學學報,2014,35(5):410-416.

PAN Nan,YI Zeguang. Study of brake abnormal sound localization method based on frequency-domain blind signal processing[J]. Journal of Hebei University of Science and Technology,2014,35(5):410-416.

[20]許海翔,叢豐裕,史習智.基于復信息量的頻域多通道盲解卷積算法[J].數據采集與處理,2008, 23(3):278-283.

主站蜘蛛池模板: 最新国产高清在线| 精品欧美一区二区三区在线| 白浆免费视频国产精品视频| 色婷婷视频在线| 国产美女精品一区二区| 色国产视频| 久久久久国产精品熟女影院| 在线不卡免费视频| 91高清在线视频| 日韩黄色大片免费看| 免费国产好深啊好涨好硬视频| 国产日本欧美亚洲精品视| 五月天综合网亚洲综合天堂网| 日韩精品一区二区三区swag| 国产呦精品一区二区三区网站| 欧美成人午夜视频| 亚洲欧美日韩中文字幕在线一区| 国产69精品久久久久孕妇大杂乱| 国产在线91在线电影| 日韩中文欧美| 91热爆在线| AⅤ色综合久久天堂AV色综合| 伊人国产无码高清视频| 日韩无码黄色网站| 综合亚洲网| 热99re99首页精品亚洲五月天| 97视频免费在线观看| 97se亚洲综合在线| 很黄的网站在线观看| 狠狠色香婷婷久久亚洲精品| 热久久这里是精品6免费观看| 亚洲精品视频免费看| 久久窝窝国产精品午夜看片| 5555国产在线观看| 青青青视频91在线 | 色欲国产一区二区日韩欧美| 亚洲国产黄色| 九色综合伊人久久富二代| 草草线在成年免费视频2| 欧美69视频在线| 无码av免费不卡在线观看| jijzzizz老师出水喷水喷出| www中文字幕在线观看| 天天色天天操综合网| 欧美日韩国产综合视频在线观看 | 国产三级韩国三级理| 中文天堂在线视频| 精品一区二区三区中文字幕| 久久久久亚洲AV成人网站软件| 99热国产这里只有精品无卡顿"| 国产免费精彩视频| 日本欧美一二三区色视频| 久久久久国产一级毛片高清板| 亚洲精品第1页| 亚洲综合18p| 亚洲精品成人福利在线电影| 青青青国产在线播放| 成人亚洲国产| 亚洲欧美另类视频| 四虎精品国产AV二区| 中国丰满人妻无码束缚啪啪| 一级片免费网站| 欧美亚洲国产视频| 无码一区18禁| 中文字幕欧美成人免费| 欧美日本激情| 激情六月丁香婷婷| 天堂中文在线资源| 国产精品久久久久无码网站| 91麻豆国产视频| 久草中文网| 国产在线八区| 国产av剧情无码精品色午夜| 99精品这里只有精品高清视频| 欧美日韩导航| 成人一区在线| 亚洲国产成人麻豆精品| 97se亚洲综合在线韩国专区福利| 久久精品国产999大香线焦| 97狠狠操| 91高清在线视频| 精品超清无码视频在线观看|