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

基于HBA-ICEEMDAN和HWPE的行星齒輪箱故障診斷*

2023-08-31 03:23:54陳愛午王紅衛
機電工程 2023年8期
關鍵詞:特征提取優化故障

陳愛午,王紅衛

(1.江蘇聯合職業技術學院 泰興分院,江蘇 泰興 225400;2.東南大學 機械工程學院,江蘇 南京 210096)

0 引 言

行星齒輪箱目前被廣泛應用于風力發電領域。由于工作環境惡劣,齒輪箱輪齒的關鍵部位易發生點蝕、斷齒等單點故障或復合故障,若故障持續演化將誘發嚴重的安全事故。因此,準確可靠地檢測行星齒輪箱的故障具有重大意義[1,2]。

目前,研究人員一般基于振動信號來對行星齒輪箱進行故障診斷和健康監測[3,4]。李昌林等人[5]提出了集成經驗模態分解(ensemble empirical mode decomposition,EEMD),利用EEMD引入白噪聲,對振動信號進行分解,以改變信號的極值特征;但EEMD依然存在模態混疊現象。TORRES M E等人[6]提出了自適應噪聲完備經驗模態分解(complete ensemble empirical mode decomposition with adaptive noise,CEE-MDAN),在信號分解過程中引入自適應的白噪聲使重構信號非常精確;但CEEMDAN分解的分量中依然存在噪聲和偽分量。為此,COLOMINAS M A等人[7]對CEEMDAN進行了優化,提出了改進自適應噪聲完備經驗模態分解(ICEEMDAN),改善了IMF分量的噪聲殘留問題。顧云青等人[8]將ICEEMDAN與排列熵相結合,用于提取滾動軸承的故障特征,結果證明了ICEEMDAN的優越性;但ICEEMDAN方法的參數均人為設置,盲目設置參數無法發揮算法的最佳性能。為實現EEMD參數的自適應設置目的,CHEN WEI-jia等人[9]利用人工蜂群算法對EEMD的參數進行了優化,實現了噪聲幅值系數自適應設置的目的;但其未對總體平均次數進行優化。王海鋒等人[10]采用細菌覓食算法對EEMD的噪聲幅值系數和總體平均次數進行了協同優化,取得了較好的優化結果;但EEMD的性能不佳。

針對ICEEMDAN的參數優化問題,筆者引入蜜獾算法(HBA)對ICEEMDAN的參數進行優化搜索,以實現參數自適應最優化設置的目的,并利用HBA-ICEEMDAN剔除行星齒輪箱振動信號中的噪聲。

完成行星齒輪箱振動信號的降噪重構后,需要進行特征提取。以樣本熵和排列熵為代表的特征提取指標目前被廣泛用于行星齒輪箱的故障診斷中[11]。但上述指標均只是從單一尺度來表征信號的復雜度,難以全面地描述信號的故障特性。

為此,謝棕等人[12]采用多尺度排列熵(multiscale permutation entropy,MPE)提取齒輪信號的故障特征,取得了較高的識別準確率;但MPE無法提取信號的高頻信息。黃海濱[13]基于層次分析和排列熵,提出了層次排列熵(hierarchical permutation entropy,HPE)方法,并用于滾動軸承的故障診斷,結果表明,相較于MPE,HPE對特征提取得更加全面和充分;但HPE忽略了信號的幅值信息,特征存在遺漏。

針對上述問題,筆者提出一種基于HBA-ICEEM-DAN、層次加權排列熵(HWPE)的行星齒輪箱故障特征提取方法。

該方法利用HBA-ICEEMDAN的自適應信號分解優勢,結合HWPE的特性(能夠準確反映齒輪箱信號動態特性),完成信號的故障特征提取工作;最后,采用支持向量機(經灰狼算法優化的)對行星齒輪箱的故障進行診斷。

為驗證行星齒輪箱故障診斷方法的有效性和優越性,筆者設置多組對比實驗,進行多個維度的綜合分析,包括信號分解算法對比、熵值方法對比和分類器對比實驗。

1 算法原理

1.1 蜜獾算法優化的改進自適應噪聲完備經驗模態分解

1.1.1 改進自適應噪聲完備經驗模態分解

ICEEMDAN方法中加入的是第h個IMF分量(白高斯噪聲經過EMD分解后獲得),是一種特殊的噪聲Eh[wi]。ICEEMDAN的計算過程如下:

1)向原始信號X中添加白噪聲E1[wi]:

X(i)=X+β0·E1[wi]

(1)

式中:wi為被加入的第i個白噪聲;

2)利用ICEEMDAN方法對原始信號進行分解,得到第1個IMF分量:

(2)

3)求解第2個IMF分量:

(3)

4)以此類推,能夠得到第h個IMF分量:

(4)

1.1.2 蜜獾算法

通過模擬蜜獾挖掘和覓食蜂蜜的行為,HASHIM F A[14]提出了蜜獾算法,原理如下:

1)種群初始化

隨機初始化蜜獾個體的坐標:

xi=lbi+r1·(ubi-lbi)

(5)

式中:xi為第i個蜜獾個體的坐標;ubi,lbi為尋優空間的上邊界和下邊界;r1為[0,1]內的隨機變量;

2)嗅覺系數定義

嗅覺系數Ii與蜜獾和蜂蜜之間的距離di以及蜂蜜的聚集程度S相關,可寫為:

(6)

S=(xi-xi+1)2

(7)

di=xprey-xi

(8)

式中:r2為[0,1]內的隨機變量;

3)密度因子更新

密度因子a的大小隨著迭代次數的增加而不斷衰減,以保證從探索到采集的平滑過渡,可寫為:

a=C·exp(-t/tmax)

(9)

式中:t為當前迭代次數;tmax為總的迭代次數,C為不小于1的常數(默認為2);

4)個體坐標更新

蜜獾坐標的更新包括挖掘階段和采蜜階段。在挖掘部分,蜜獾在蜂巢周圍搜尋,其坐標的更新基于嗅覺系數Ii、蜜獾和蜂蜜之間的距離di和密度因子a。挖掘部分的數學形式如下:

xnovel=xprey+F·β·Ii·xprey+F·r3·

a·di·|cos(2π·r4)·[1-cos(2π·r5)]|

(10)

式中:r3,r4,r5為[0,1]內的隨機變量;β≥1(默認為6)為蜜獾得到蜂蜜的能力;xprey為蜂蜜的坐標,也即最佳個體的坐標;F為方向控制系數,能夠增強蜜獾探索的全局性能。

F的算式如下:

(11)

式中:r6為[0,1]內的隨機變量。

在采蜜部分,蜜獾追尋向導鳥直接往蜂巢移動,蜜獾根據距離信息di在xprey周圍進行探索。在該部分,探索受密度因子a的影響,探索的數學形式為:

xnovel=xprey+F·r7·a·di

(12)

式中:r7為[0,1]內的隨機變量。

1.1.3 參數優化流程

由于ICEEMDAN方法的分解效果取決于白噪聲幅值權重(Nstd)和噪聲添加次數(NE),因此,筆者采用HBA對ICEEMDAN的2個參數進行優化,適應度函數采用包絡熵,當適應度值越小,則代表分解的效果越好;通過優化和更新,來確定最終的最佳參數組合(Nstd,NE)。

優化的步驟如下:

1)HBA算法的種群初始化,設置HBA的迭代次數和種群規模,并設置ICEEMDAN算法的參數優化范圍;

2)利用ICEEMDAN分解齒輪箱振動信號,并計算各個IMF分量的包絡熵,以包絡熵的最小值為適應度函數;

3)判斷優化是否達到算法的終止條件,若是,則繼續下一步;若否,則更新種群位置,并返回第2)步;

4)保存最優的ICEEMDAN參數組合,并將其代入至ICEEMDAN中,為ICEEMDAN構建HBA算法;

5)利用HBA-ICEEMDAN方法分解齒輪箱振動信號,得到最佳的IMF分量。

1.2 層次加權排列熵

加權排列熵僅從單一尺度來描述信號的動態特性,描述得不夠全面,為此,陳柯宇等人[15]提出了多尺度加權排列熵(multi為scale weighted permutation entropy,MWPE),實現了信號多尺度計算的目的。但粗?;幚肀举|上是一個低通濾波器,其過濾了信號的高頻特征信息,因此,筆者采用層次處理方式對時間序列進行分割,結合加權排列熵,提出了層次加權排列熵[16]。

HWPE能夠同時提取信號的低頻和高頻節點信息,對信號進行全面分析,計算流程如下:

1)對于信號{u(i),i=1,…,N},定義平均算子Q0和差分算子Q1,公式如下:

(13)

(14)

其中,N=2n,n是正整數。算子Q0和Q1的長度為2n-1。

利用算子Q0和Q1,則原始信號可以重構為:

u={Q0(u)j+Q1(u)j,Q0(u)j-

Q1(u)j}j=0,1,2,…,2n-1

(15)

當j=0或j=1,定義矩陣Qj算子如下所示:

Qj(u)=

(16)

2)生成一個n維向量[δ1,δ2,…,δn]∈{0,1},則定義整數γ為:

(17)

其中,正整數γ對應于向量[δ1,δ2,…,δn];

3)在向量[δ1,δ2,…,δn]∈{0,1}的基礎上,定義信號u(i)的每一層分解的節點序列如下:

uk,e=Qδn·Qδn-1·…·Qδ1(u)

(18)

式中:k為層次分解中的第k層。

當k=2時,原始信號的兩層層次分解的示意圖如圖1所示。

圖1 信號u(i)的層次分解示意圖

4)計算每個層次分量的加權排列熵,得到2k個層次分量的加權排列熵值,即HWPE定義為:

HWPE(u,k,γ,m,t)=WPE(uk,γ,m,t)

(19)

綜上所述,Q0和Q1算子是信號的低頻和高頻成分,與Haar小波的低通和高通濾波器原理一致;

在圖1中最左側的分解節點u1,0和u2,0的加權排列熵值分布對應于多尺度計算中各尺度的加權排列熵值,也就是分解節點uk,0對應于多尺度計算中各尺度的加權排列熵值。

上述分析說明,多尺度加權排列熵僅分解了信號低頻部分的故障信息,遺漏了高頻成分的故障信息,而層次加權排列熵在提取信號高頻信息的同時,還提取了低頻成分的故障信息。

在實際工作中,行星齒輪箱振動信號的高頻部分也包含了關鍵的損傷信息,單一低頻成分的信息無法全面表征齒輪箱損傷的固有特性,這證明了對信號進行層次分析的必要性。

1.3 灰狼算法優化支持向量機

1.3.1 灰狼算法

針對支持向量機算法的懲罰系數c和核函數g需要優化的問題,筆者采用灰狼算法進行優化[17]。GWO的種群結構為金字塔形,最上方為α狼,第2層為β狼,第3層為δ狼,最下層ω狼是整個種群的基礎。狼群狩獵包含3個階段,分別是追捕階段、圍獵階段和攻擊階段。

根據算法生成GWO數學模型,設定t為迭代次數,XP為獵物坐標矩陣,X為種群坐標矩陣,則種群與獵物之間的距離d如下所示:

d=|C·Xp(t)-X(t)|

(20)

種群不斷更換坐標,如下所示:

X(t+1)=Xp(t)-A·d

(21)

式中:A,C為系數向量。

A,C可寫為:

(22)

式中:r1,r2為范圍0到1內的隨機數;a為收斂系數。

a可寫為:

(23)

式中:Tmax為允許的迭代次數。

狼群抓捕獵物坐標變換如下所示:

(24)

(25)

式中:dα,dβ和dδ為對應α,β和δ狼與獵物的距離;A1,A2,A3和C1,C2,C3為對應于α,β,δ的參數向量;Xα(t),Xβ(t),Xδ(t)為t時刻獵物的詳細坐標;X1,X2,X3為對應種群的坐標;Xα,Xβ,Xδ為對應獵物的坐標。

1.3.2 灰狼算法優化支持向量機流程

GWO算法優化支持向量機的流程如下:

1)設置狼種群的規模、允許的迭代次數、優化參數c和g的允許優化范圍;

2)隨機生成灰狼種群,各灰狼種群的個體坐標包含兩個維度,分別是c和g;

3)根據初始給定的參數c和g,對支持向量機進行訓練,個體的適應度值由訓練樣本的診斷準確率來表征;

4)計算每條灰狼的適應度,依據適應度值的大小將灰狼劃分為α,β,δ和ω,4個等級,對種群中的每個個體進行坐標更新;

5)若迭代次數超出允許次數,則結束優化,輸出最佳的組合參數;否則,繼續執行步驟4),直至獲得最優的參數。

2 行星齒輪箱故障診斷方法

筆者采用HBA對ICEEMDAN的參數進行優化,并利用具有最佳參數的ICEEMDAN對行星齒輪箱的振動信號進行分解,得到若干個IMF分量,篩選出具有較大相關系數的分量,進行重構;隨后,利用HWPE提取重構信號的故障特征;最后,利用GWO-SVM進行故障識別和分類。

基于HBA-ICEEMDAN-HWPE-GWO-SVM的行星齒輪箱故障診斷方法算法的具體流程如下:

1)利用HBA-ICEEMDAN對行星齒輪箱的振動信號進行分解,得到若干個IMF分量;

2)根據相關系數篩選出系數大于0.2的分量,進行重構;

3)利用HWPE提取重構信號的熵值,生成齒輪箱的故障特征樣本;

4)利用GWO對支持向量機的參數進行優化,并對特征樣本進行訓練和測試,輸出識別結果,完成故障類型的診斷工作。

基于HBA-ICEEMDAN-HWPE-GWO-SVM的行星齒輪箱故障診斷方法的流程圖如圖2所示。

圖2 基于HBA-ICEEMDAN-HWPE-GWO-SVM的行星齒輪箱故障診斷方法的流程圖

3 故障診斷實驗及分析

3.1 實驗數據

為了驗證基于HBA-ICEEMDAN-HWPE-GWO-SVM的行星齒輪箱故障診斷方法的有效性,筆者利用行星齒輪箱的振動數據進行實驗?;赒PZZ-II旋轉機械故障模擬實驗臺進行實驗數據的采集,利用該平臺能夠模擬行星齒輪箱的多種故障形式。

QPZZ-II旋轉機械故障模擬實驗平臺[18]如圖3所示。

圖3 QPZZ-II型旋轉機械故障實驗臺

齒輪箱實驗臺的結構配置為:小齒輪齒數為55,大齒輪齒數為75,模數為2 mm。實驗模擬了齒輪箱在轉速為1 500 r/min下的6種損傷類型。

樣本的詳細信息如表1所示。

表1 齒輪箱的損傷類型

表1中,6種損傷類型的振動信號由布置在齒輪箱上的加速度傳感器進行采集而得到。傳感器布置在輸入軸電機側軸承的Y方向;信號的采樣頻率為5 120 Hz。

筆者選擇2 048個數據點為一組樣本,每種工況選擇45組樣本,則共能夠獲得270組樣本;將樣本依據6∶4的比例進行分割,得到訓練樣本和測試樣本。

齒輪箱6種損傷振動信號的波形如圖4所示。

圖4 齒輪箱振動信號波形

從圖4可以發現:僅通過波形無法直接判斷齒輪箱的故障類型,需要采用更加智能的方法進行故障的識別。

3.2 特征提取

3.2.1 HBA-ICEEMDAN分解與信號重構

首先,筆者利用HBA優化ICEEMDAN的2個關鍵參數。其中,噪聲幅值權重的優化范圍設置為[0.15,0.6],噪聲添加次數的優化范圍為[50,600],而篩選迭代次數設置為100次。

為了驗證HBA算法在優化中的有效性,筆者利用粒子群算法(particle swarm optimization,PSO)、遺傳算法(genetic algorithm,GA)和鯨魚優化算法(whale optimization algorithm,WOA)進行對比;每種算法的種群規模和迭代次數均為10和30,以局部包絡熵最小為適應度函數。

4種算法的優化迭代曲線如圖5所示。

圖5 4種算法的適應度迭代曲線

從圖5可以發現:在迭代3次后,HBA就達到了最優,而且其最優解優于另外3種算法,證明了HBA在優化過程中的高效和性能。GA最終收斂于一個較大的適應度值,這證明了GA的全局優化性能較差,易陷入局部最優。WOA在迭代過程中未實現收斂目的,證明局部尋優能力較差,無法收斂。

經過優化后,4種算法優化得到的ICEEMDAN最佳參數如表2所示。

表2 ICEEMDAN最佳參數

筆者將經過HBA優化得到的參數代入至ICEEMDAN中,對齒輪箱振動信號進行分解。

以健康狀態信號的分解為例,其分解結果如圖6所示。

圖6 健康樣本的HBA-ICEEMDAN分解波形

隨后,為了去除信號中的噪聲,增強故障特征對齒輪箱工況的敏感性,筆者利用相關系數評估各IMF分量與原始分量的緊密度。

各分量與原始信號的相關程度如圖7所示。

圖7 IMF分量的相關系數

從圖7可以發現:隨著IMF分量階數的增加,相關系數也隨之減小,這證明低階IMF分量包含了主要的信息。而根據相關系數的相對關系來看,當階數大于4后,IMF分量的相關系數也比較小,因此,筆者選擇前4階IMF進行重構。

3.2.2 HWPE特征提取

行星齒輪箱的振動信號具有復雜的動態特性,為了準確地表征系統的復雜性,筆者采用HWPE對重構信號進行分析,提取故障特征。HWPE的參數設置為:嵌入維數m=5,時間延遲d=1,分解層數k=3。

首先,為了充分評估不同損傷部位的振動信號差異,筆者利用加權排列熵對6種振動信號進行分析,結果如圖8所示。

圖8 不同工況振動信號的加權排列熵

從圖8可以發現:不同工況樣本的加權排列熵分布得較為混亂,不具有明顯的差異性;但是可以發現點磨樣本的熵值普遍較小,這表明信號的復雜度可能較小。這是因為當齒輪正常運轉時,振動信號由嚙合沖擊和環境噪聲等相互作用而產生,信號比較不規則;而齒輪發生點磨故障后,信號中出現了周期性較為顯著的嚙合成分,復雜性降低,因此熵值較小。同時,根據該分析也可證明,對信號進行多尺度分析是有必要性的。

隨后,筆者對經HBA-ICEEMDAN分解和重構后的信號進行HWPE分析,得到6種齒輪箱工況在8個尺度下的熵值,如圖9所示。

圖9 振動信號的HWPE均值

由圖9可知:在各個特征尺度上,熵值特征具有一定的差異;但隨著尺度的增加,這種差異逐漸減小。

為了增加特征的區分度,提高故障識別的區分度,筆者選擇前4個尺度的特征來構造故障特征樣本。

3.2.3 模型的診斷

筆者將提取的故障特征按照6∶4的比例隨機地分割為訓練樣本和測試樣本,即從每組工況45個樣本中,隨機抽取27個作為訓練樣本,18個作為測試樣本(總共有135個訓練樣本,90個測試樣本),利用GWO-SVM對其進行故障識別。

測試樣本的GWO-SVM分類結果如圖10所示。

圖10 模型的混淆矩陣

從圖10可以發現:有5類齒輪箱樣本均被準確識別,只有第2類樣本(點蝕故障)被錯誤地識別為第1類樣本(健康),準確率為88.89%;在整個測試樣本中,被準確識別的樣本為106個,被錯誤識別的樣本為2個,總的識別準確率為98.15%。

3.3 模型比較

3.3.1 分類器的優化算法對比

為了研究不同參數尋優算法對模型準確率分類時間的影響,以突出GWO-SVM分類器的有效性,筆者使用遺傳算法(GA)、粒子群算法(PSO)、人工魚群算法(artificial fish school algorithm,AFSA)和蝙蝠算法(bat algorithm,BA)優化SVM的參數,并進行測試樣本的識別。

筆者記錄了各個分類器的識別準確率和運行時間,實驗中,懲罰系數和核函數的優化范圍設置為[0.01,100]。

不同分類器的診斷結果如表3所示。

表3 不同分類器的診斷結果

由表3可以發現:各個分類器都能取得超過90%的識別準確率,準確率最低的是BA-SVM,為92.59%;GA-SVM的分類時間最長,需要1.97 s;GWO-SVM的分類準確率最高,為98.15%,分類時間也僅需1.16 s。

根據上述結果可以證明,GWO-SVM在齒輪箱故障識別中具有一定的優越性。

3.3.2 信號分解算法對比

為了證明HBA-ICEEMDAN方法的有效性和優越性,筆者利用GA、PSO、WOA對ICEEMDAN的參數進行優化,并利用HWPE進行特征提取,采用GWO-SVM進行故障識別。此外,為了驗證信號去噪的重要性,筆者將原信號輸入到故障診斷模型中進行分類,并將結果進行對比。

降噪效果對比結果如圖11所示。

圖11 降噪效果對比

從圖11可以發現:直接利用原信號進行故障識別的表現最差,這是因為原信號中包含較多的噪聲,直接進行特征提取會受到較多的干擾;而就不同信號分解算法而言,HBA-ICEEMDAN的準確率最高,這是因為HBA算法具有較好的全局和局部尋優能力,優化所得到的參數能夠使IMF分量質量更高。

基于上述分析,HBA-ICEEMDAN方法能夠很好地消除信號中的噪聲分量,突出故障特性,以便于后續的特征提取和模式識別。

3.3.3 不同熵值指標對比

為了對比不同熵值指標在特征提取中的有效性和表現,以驗證HWPE方法的有效性,筆者分別采用層次排列熵(HPE)、多尺度排列熵(MPE)和多尺度加權排列熵(MWPE)來提取故障特征,結果如圖12所示。

圖12 齒輪箱振動信號的熵均值

從圖12可以發現:各個方法在大多數特征尺度上都能夠區分不同的故障類型;但健康和點蝕故障的熵曲線存在重疊和交織現象,區分效果需要進一步評估。

為了更準確地判斷HWPE方法的有效性,筆者分別將HWPE、HPE、MPE和MWPE提取的故障特征輸入至GWO-SVM分類器中,進行故障的識別,結果如圖13所示。

圖13 10次分類的平均診斷準確率

從圖13可以發現:利用HWPE進行特征提取所獲得的準確率最高,其次是HPE,而MPE方法的準確率最低,這證明了利用HWPE進行特征提取的優越性(造成這種現象的原因是因為HWPE和HPE充分提取了重構信號的高頻特征,而MWPE和MPE僅提取了重構信號的低頻特征,特征提取存在遺漏,因此無法全面地描述齒輪箱的故障特性)。

基于上述分析可知,利用HWPE進行特征提取是有效的,且優于其他對比方法。

4 結束語

針對行星齒輪箱的故障特征提取和故障診斷問題,筆者提出了一種結合HBA-ICEEMDAN、HWPE和GWO-SVM的行星齒輪箱故障診斷模型(方法)。

首先,筆者采用HBA算法優化了ICEEMDAN的2個關鍵參數,以剔除振動信號中的噪聲;然后,計算重構信號的HWPE熵值,生成故障特征樣本;最后,利用GWO優化支持向量機的參數,并對特征樣本進行了故障識別。

經過齒輪箱數據分析,得出了以下結論:

1)利用HBA對ICEEMDAN的參數進行優化是有效的,能夠實現參數自適應設置的目的,且HBA的優化效果優于GA、PSO和WOA的結果;經過HBA-ICEEMDAN算法去噪和HWPE特征提取的故障特征能夠準確地反映樣本的故障特性;

2)在HBA-ICEEMDAN-HWPE特征提取的基礎上,GWO-SVM模型的診斷準確率達到了98.15%;對比其他分類器,GWO-SVM的診斷準確率和效率最高;

3)和其他熵值指標相比,基于HWPE的特征提取方法具有更高的準確率,明顯優于HPE、MPE和MWPE。

雖然筆者所提出的故障診斷方法能較好地對齒輪箱進行故障識別,但在參數優化方面仍然存在不足,特別是ICEEMDAN的迭代次數仍需人為設置。

因此,在后續研究中,筆者將進一步對上述故障診斷方法進行優化,以實現參數自適應設置的目的。

猜你喜歡
特征提取優化故障
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
故障一點通
基于Gazebo仿真環境的ORB特征提取與比對的研究
電子制作(2019年15期)2019-08-27 01:12:00
一種基于LBP 特征提取和稀疏表示的肝病識別算法
奔馳R320車ABS、ESP故障燈異常點亮
故障一點通
江淮車故障3例
主站蜘蛛池模板: 99久久国产自偷自偷免费一区| 日本午夜精品一本在线观看| 人妻中文久热无码丝袜| 精品人妻系列无码专区久久| 国产欧美日本在线观看| 日本三区视频| 伊在人亞洲香蕉精品區| 麻豆a级片| 在线精品亚洲国产| 91色在线视频| 国产综合亚洲欧洲区精品无码| 无码久看视频| 幺女国产一级毛片| 无码久看视频| 五月综合色婷婷| 自慰网址在线观看| 亚洲日韩图片专区第1页| 国产午夜小视频| 国产自在线播放| 国产精品3p视频| 国产手机在线观看| 成人一级黄色毛片| 2021无码专区人妻系列日韩| 麻豆国产在线观看一区二区| 国产三级国产精品国产普男人| 亚洲av综合网| 无遮挡一级毛片呦女视频| 亚洲天堂精品在线| 97精品伊人久久大香线蕉| 大陆国产精品视频| 美女内射视频WWW网站午夜| 97国产精品视频自在拍| 国产农村精品一级毛片视频| 国产凹凸视频在线观看| 久久福利网| 欧美日韩中文国产va另类| 欧美一道本| 免费在线成人网| 中文字幕无码电影| 一级香蕉视频在线观看| 国产系列在线| 亚洲天堂福利视频| 四虎免费视频网站| 亚洲天堂福利视频| 超碰色了色| 女同久久精品国产99国| 小蝌蚪亚洲精品国产| 直接黄91麻豆网站| 免费高清a毛片| 日本高清视频在线www色| 婷婷亚洲视频| 999在线免费视频| 国产日本一线在线观看免费| 久久这里只精品国产99热8| 日韩无码一二三区| 欧美国产日产一区二区| 中国国产A一级毛片| 国产日韩av在线播放| 亚洲日本一本dvd高清| 中文字幕自拍偷拍| 亚洲综合香蕉| 亚洲一级毛片在线播放| 欧美日韩亚洲国产主播第一区| 自偷自拍三级全三级视频| 久久 午夜福利 张柏芝| 天天综合亚洲| 国产精品v欧美| 国产一级一级毛片永久| 免费无码一区二区| 国产一区二区三区在线观看免费| 国产午夜福利亚洲第一| 一级毛片基地| 成人无码一区二区三区视频在线观看 | 四虎永久在线精品影院| 精品日韩亚洲欧美高清a | 2021国产精品自产拍在线| 国产成人午夜福利免费无码r| 色偷偷一区| 国产00高中生在线播放| 日韩一二三区视频精品| 99re视频在线| 国产精品白浆无码流出在线看|