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

面向不平衡數據的云模型旋轉機械故障識別方法

2022-12-01 10:25:00趙榮珍
振動與沖擊 2022年22期
關鍵詞:模式識別分類特征

趙 楠,趙榮珍

(蘭州理工大學 機電工程學院,蘭州 730050)

旋轉機械故障診斷技術對于流程工業的安全可靠運行發揮著重要作用[1]。隨著該類裝置必須向優質高效運行目標邁進,其監測信息的規模也在擴大,這使得旋轉機械故障監測系統獲取海量數據,從而推動旋轉機械故障診斷領域進入“大數據”時代[2]。面對積累的工業大數據資源,必須使用合理、高效的智能決策技術將其開發利用,這對于快速發展旋轉機械智能制造至關重要[3-4]。

故障模式識別是旋轉機械故障診斷的本質問題,其目的在于獲得高識別精度[5]。為此,多種分類器應運而生,典型的分類算法包括神經網絡(neural network)、支持向量機(support vector machine,SVM)、K最近鄰分類算法(K-nearest neighbor,KNN)等。這些分類算法對于平衡數據能夠取得較好的識別效果,然而對于不平衡數據中少數類樣本的識別精度不高[6]。其原因在于這些算法受到數據不平衡的影響,在特征學習過程中偏向于多數類樣本學習,而對少數類樣本普遍存在著學習不足的情況,導致分類器存在對于少數類樣本識別精度較低的缺陷[7]。在實際情況下的機械設備監測中,設備故障的發生屬于少數類,檢測的數據大多數為正常工作數據,這種情況容易使得分類器將故障狀態數據誤判為正常狀態數據,進而錯過機械設備的最佳維護時間,造成難以估計的后果和損失[8]。

就不平衡數據的分類算法而言,目前較為典型的方法包括集成學習法、代價敏感法、單類學習法等。其中,集成學習法根據不平衡數據構建多個具有差異性的分類器,并按照一定方式將分類器識別結果整合以提高整體分類器對于不平衡數據的識別精度。但它存在著多個分類器訓練時間長、基分類器類型和數量選擇較為困難的缺陷。代價敏感法則以錯分代價這一概念為基礎,著重關注錯誤代價較高類別所對應的樣本,并使得分類錯誤的總代價最低,從而優化分類算法。其缺陷在于難以真實準確地估計錯分代價,并且不適用于少數類樣本數量過少的情況。而單類學習法僅利用單一類別的訓練數據進行訓練,能夠減少訓練時間,對于處理少數類樣本數量極少的情況具有一定優勢。然而該方法在訓練少數類樣本時容易發生過擬合現象,導致分類方法泛化能力偏低。為尋找一種泛化能力強、識別精度高的不平衡數據模式識別方法,本文擬將云理論中的云模型引入對于旋轉機械不平衡故障數據的處理過程中。

云模型具有優良的數學性質,能夠實現數據的可視化效果,因此被廣泛應用于各個領域,模式識別是云模型應用的重要方面之一[9-13]。云模型能夠實現定性概念與定量值間的相互轉換,在以不平衡數據訓練分類器的過程中,通過逆向云發生器(backward cloud generator,BCG)將多數類樣本與少數類樣本轉換為期望、熵、超熵3個數字特征,再通過正向云發生器(forward cloud generator,FCG)將數字特征繪制成相同云滴點數的云圖,有效地減少了不平衡數據對于分類過程的影響,解決對不平衡數據分類時少數類樣本識別精度較低的問題。

基于上述緣由,故本文擬對利用云模型改善旋轉機械故障數據集質量問題進行研究,將借助于云圖概念改進云分類器以解決不平衡故障數據集的分類與識別問題。欲為工業大數據的智能決策技術實現,提供數據運算的基礎理論依據。

1 云模型理論介紹

云模型理論是李德毅院士以概率論與模糊數學理論兩者為基礎所提出的研究定性概念與其定量表示的認知計算模型,其賦予樣本點隨機確定度以統一刻畫概念中的模糊性、隨機性以及關聯性[14]。生成云模型最為關鍵的兩個步驟是構建正向、逆向云發生器。其中:逆向云發生器計算樣本數據的3個數字特征,既期望(Ex)、熵(En)、超熵(He);正向云發生器依據逆向云發生器輸出的3個數字特征繪制云圖。兩者相輔相成。

云模型按照云滴概率分布的不同,可分為正態云、高斯云、冪律云等。正態分布是概率分布中最基本、最重要、最為廣泛應用的模型。因此,正態隸屬函數也成為模糊理論中最常使用的隸屬函數,本文所使用的正態云模型正是以正態分布和正態隸屬函數為基礎構建的。

1.1 云模型的數字特征

云模型通過逆向云發生器與正向云發生器實現定量數據與定性概念之間的相互轉換,以期望、熵、超熵這3個數字特征描述一個不確定性概念,這3個數字特征在云圖中的表示如圖1所示,所代表的含義如下:

期望Ex——期望是云滴在論域空間的中心值,也是最能代表定性概念的點。

熵En——熵值既論域空間中云滴可被定性概念接受的取值范圍,是定性概念不確定性的度量,能夠反應云滴的離散程度。

超熵He——超熵是熵值的熵,既熵值的不確定性度量,代表云滴的凝聚程度,超熵值越小則云滴的離散程度也越小。

圖1 云圖中的數字特征Fig.1 Digital features in cloud chart

1.2 正向云發生器

正向云發生器能夠完成定性概念到定量值的轉換,它能夠通過3個數字特征生成云圖中的云滴。構建正態云模型時所使用的二階正向正態云發生器算法FCG(Ex,En,He,n)如圖2所示,具體步驟如下:

輸入(Ex,En,He,n)3個數字特征及生成的云滴個數n

輸出n個云滴xi及其所對應的隸屬度μ(xi)

算法步驟:

步驟1生成一個正態隨機數yi=RN(En,He),yi以En為期望、以He2為方差;

步驟4根據隸屬度μ(xi)以及正態隨機數xi在數域中生成一個云滴;

步驟5循環步驟1~步驟4,直至生成n個云滴(i=1∶n)。

圖2 正向云發生器Fig.2 Forward cloud generator

1.3 逆向云發生器

與圖2所示的正向云發生器相反,逆向云發生器能夠完成定量值到定性概念的轉換,它將精確數值轉換為3個數字特征所表示的定性概念。本文選擇較為常用的無需確定隸屬的逆向云分類器算法BCG(xi),如圖3所示,具體步驟如下:

輸入樣本點xi(i=1,2,…,n)

輸出反映定性概念的3個數字特征Ex,En,He

算法步驟:

圖3 逆向云發生器Fig.3 Backward cloud generator

對于云模型在模式識別中的應用,多數方法是將由正向云發生器改進出的隸屬度公式與逆向云發生器相結合,由此構造云分類器對數據進行處理分析[16-18]。這種分類方法沒有很好的利用所繪制的云圖,得到的識別結果與云圖并無直接關系。為此,本文提出一種改進的云分類器,依據云圖間的位置關系判斷出測試數據所屬類別。將這種改進的云分類器應用于旋轉機械故障診斷的數據不平衡情況,提高模式識別精度的同時,實現識別結果的可視化。

2 設計的不平衡故障診斷方法與流程

融合了云模型的數據挖掘方法考慮到模糊性與隨機性,其結果相較于多種數據挖掘方法更加合理,有利于數據挖掘的智能化[19]。正態云模型通過逆向云發生器與正向云發生器二者結合的方式,將不平衡數據中的多數類樣本與少數類樣本繪制成云滴點數相同的云圖,能夠有效的降低不平衡數據對于模式識別結果的影響。本文借助云模型對不平衡數據分類的優勢,以繪制云圖為基礎提出一種改進的云分類器。

2.1 改進的云分類器

構建云分類器是云模型理論應用于模式識別領域的關鍵環節。本節依據各類云圖間的位置關系提出改進云分類器。首先,通過逆向與正向云發生器繪制各類數據云圖;然后,比較測試數據云圖與各類訓練數據云圖之間中心線的距離,以對測試數據進行分類;最后,選擇距離最小的訓練數據云圖對應類別作為測試數據的判別結果,完成模式識別的同時實現特征數據的可視化表達。以3個類別的訓練數據為例,測試數據通過改進云分類器進行分類的可視化識別結果表達,如圖4所示。其中:Ex1,Ex2,Ex3對應3種訓練數據云圖的中心線;Ex對應測試數據云圖的中心線。由圖4可知,測試數據云圖與第二類訓練數據云圖的中心線距離最短,因此將測試數據歸為第二類別。

圖4 改進云分類器的可視化結果表達Fig.4 Improve the visualization result expression of cloud classifier

由逆向云發生器及正向云發生器的過程推導可知,對于某一特征下的單個待測樣本,其期望值為特征值本身,熵值與超熵值為零,使得該樣本能夠計算得出云圖中心線對應于x軸的坐標值(即期望Ex),而無法求解出云滴在云圖中對應于y軸的坐標值(即隸屬度μ(xi)),不能繪制出單個待測樣本的云圖,進而模式識別也難以進行。為解決該問題,依據云圖中位置關系提出類別判別公式D=|Ex-Ext|,將待測樣本的期望值減去各個訓練樣本的期望值,絕對值最小的差值所對應訓練樣本狀態即為待測樣本狀態。通過上述判別準則,以測試數據繪制云圖、以判別公式對待測樣本進行分類,完成的改進云分類器構建流程如圖5所示,具體應用步驟如下:

輸入單一特征下的t組訓練數據與測試數據、云滴個數n,其中t為數據的類別個數

輸出識別結果d

算法步驟:

步驟1生成數字特征,將某一特征下各個類別的訓練數據與該特征下的測試數據作為逆向云發生器的輸入,得到t+1組數字特征(Ex,En,He);

步驟2繪制云圖,將生成的t+1組數字特征作為正向云發生器的輸入,分別得到由各組數字特征衍生出的n個云滴(xi,μi)(i=1∶n),將各組云滴繪制成所對應的訓練數據及測試數據云圖;

步驟3計算距離,通過公式D=|Exj-Ext|(j=1∶j),其中j為測試數據所包含的待測樣本數量。計算出待測樣本所對應的期望值Exj與各個狀態下訓練數據云圖所對應的期望值Ext之間差值的絕對值,即待測樣本期望值與各個訓練數據云圖中心線之間的距離;

步驟4模式識別,取t個距離中最小值Dmin所對應的期望值Exj=d(t∈1∶t),則判定第d個訓練數據云圖所對應的類別即為待測樣本的識別結果。即將待測樣本歸類于云圖中心線間距離最近的訓練數據所屬類別。

2.2 建立的模式識別方法流程

改進的云分類器能夠使識別結果與所繪制的云圖對應,實現識別結果的可視化。但該分類器在分類的過程中只利用數據集中的單一特征,不利于獲得較高識別精度,同時造成了數據資源浪費。為解決這一問題,將本文所提方法與集成學習方法相結合,以n個特征訓練對應個數的云分類器,形成n組識別標簽,再將這些標簽通過相對多數投票法進行整合,選取標簽數量最多的類別作為識別結果。

圖5 改進云分類器流程圖Fig.5 Improved cloud classifier flow chart

通過集成學習的方式將多個云分類器進行組合,使得分類方法能夠處理多個特征,但在分類過程中,并不是特征數量越多越好。過多的特征數量將會導致分類器訓練時間過長、分類模型復雜,甚至會降低模式識別精度[20]。因此,本研究采用ReliefF算法對初始特征集中的特征進行選擇,以選取出特征數量適當并有利于實施分類的低維特征集進行模式識別。構建低維特征集的過程中,通過ReliefF算法計算出初始特征集中各個特征的權重,并依據權重的降序排序對特征進行重新排序,以選取出權重較大的前m個優質特征構成低維特征集。ReliefF算法計算特征權重的過程中,首先從初始特征集中隨機選取一個樣本Ri(i=1,2,…,m),從Ri的同類樣本子集中尋找Ri的k個近鄰樣本,記為Hj(j=1,2,…,k),從Ri的不同類樣本子集中尋找Ri的k個近鄰樣本,記為Mj(j=1,2,…,k);其次將選好的各樣本代入特征A的權重w(A)迭代公式(式(1));最后對上述步驟進行m次重復,以完成權重公式的m次迭代,得到初始特征集中各個特征的權重。

(1)

式中:抽樣次數m為公式的迭代次數;diff(A,Ri,Hj)為樣本Ri與樣本Hj在特征A上的距離;class(Ri)為樣本Ri所屬類別;C≠class(Ri)為類別C與樣本Ri不為同類別;P(C)為C類別樣本占總樣本數量的比例。

為能夠達到較高的識別精度,同時減少分類算法運行時間、降低分類模型復雜度,通過網格搜索法選取最優特征數量m以及ReliefF算法中近鄰樣本數k。

將改進云分類器與上述集成學習算法、ReliefF算法相結合,所構建的不平衡數據模式識別方法流程圖如圖6所示,具體步驟如下:

步驟1采集旋轉機械的原始振動信號并進行預處理,構造初始特征集;

步驟2通過ReliefF算法對初始特征集進行特征選取,構建低維特征集;

步驟3依據低維特征集劃分出不平衡訓練集和測試集;

步驟4將單一特征下的各狀態訓練數據及待測樣本代入改進的云分類器,構建出低維特征集特征個數的云分類器,并繪制各狀態云圖;

步驟5將各個云分類器的識別結果通過相對多數投票法進行結合,輸出故障識別結果。

圖6 建立方法的流程圖Fig.6 Flow chart of the establishment method

3 試驗結果與分析

3.1 試驗數據及相關參數的選擇情況

為驗證本文所提不平衡數據模式識別方法的有效性,使用無錫厚德自動化儀表有限公司所提供的綜合故障模擬試驗臺,模擬滾動軸承運行狀態并對不同狀態下的振動信號進行采集,通過本方法對所采集的數據進行模式識別,使用的試驗臺如圖7所示。

試驗采用NSK6308型號軸承,通過加速度傳感器采集軸承在正常、內圈故障、外圈故障、滾動體故障、保持架故障共5種狀態下的振動信號。試驗設置采樣頻率為8 000 Hz,采樣轉速分別為2 600 r/min,2 800 r/min,3 000 r/min的運行狀態。對所采集的振動信號進行小波消噪,截取同一轉速下的各個狀態振動信號80組,共計400組樣本,樣本長度為1 024點。提取該樣本的時域特征6個、頻域特征10個、時頻域特征4個,共計20個特征,具體特征組成情況如表1所示。試驗臺設置有5個信號采集通道,故所構成的初始特征集特征維數為100。為便于后續計算、利于云圖可視化表達,將初始特征集進行歸一化,歸一化區間為[-1,1]。

圖7 綜合故障模擬試驗臺Fig.7 Comprehensive failure simulation test bench

表1 故障特征組成情況Tab.1 Fault feature composition

通過RelifF算法計算初始特征集中各個特征的權重,以3 000 r/min運行狀態下的初始特征集為例得到的權重狀態,如圖8所示。圖8中,特征序號1~20代表通道1振動信號的表1中所述20個特征,依此類推。將此權重降序排序,取權重較大的前m個權重所對應的特征構成低維特征集。采用網格搜索法,選取ReliefF算法中的k近鄰樣本數為7、低維特征集的特征數為29。3 000 r/min運行狀態下低維特征集的特征組成情況,如表2所示,特征排序按照其權重由大到小進行排列。以表2中序號5特征為例,其表示低維特征集中的第5個特征為通道1下的頻譜一階重心(表1中的序號10特征)。將此低維特征集按照5∶5劃分為平衡訓練集與測試集,兩者樣本數均為200。取訓練集中的正常樣本40個、內圈故障樣本20個、外圈故障樣本20個、滾動體故障樣本10個、保持架故障樣本10個,完成不平衡訓練集的構建。為便于后續試驗的表達,將本文提出的依據云圖間距離的不平衡數據模式識別方法命名為以云圖間距離為基礎的云模型模式識別方法(distance cloud model,DCT)。

圖8 特征權重圖Fig.8 Feature weight histogram

表2 低維特征集的特征組成情況Tab.2 Feature composition of low-dimensional feature set

3.2 本方法的識別性能驗證

獲得高的識別精度是模式識別研究的目的所在,為驗證本文所提方法對于不平衡數據的故障識別能力,選用3.1節所述的3 000 r/min運行狀態下的數據進行驗證,按照該節所述方法構造維數為29的低維特征集、劃分不平衡訓練集以及測試集。對于不平衡數據的各類訓練數據、各類測試數據的云圖繪制,如圖9所示。圖9中的測試數據云圖均能得到正確分類;DCT方法對于測試集的識別結果如表3所示。該結果證明了DCT方法對于不平衡數據的多數類樣本、少數類樣本均能實施正確分類,驗證了本文方法的可行性。為證明DCT方法與其他分類方法對于不平衡數據的分類具有一定優勢,將本方法與BP(back propagation)神經網絡、SVM、以及云模型和集成極限學習機相結合的滾動軸承故障診斷方法(ensemble-extreme learning machine,E-ELM),這3種故障識別算法進行識別精度對比。BP、SVM分類方法均選用之前所述的3 000 r/min運行狀態下的低維特征集所劃分的不平衡訓練集以及測試集,將不平衡訓練集代入各個分類方法進行訓練、測試集進行測試,其中SVM的核函數使用高斯徑向基函數(radial basis function,RBF),BP神經網絡的結構為29-100-5。E-ELM方法按照文獻[21]所述步驟,提取3 000 r/min 運行狀態振動信號的云特征構建初始特征集,特征集劃分訓練、測試集的比例和不平衡訓練集的構建方法均與本文方法相一致,將ELM分類器的隱層神經元個數設置為25,其余參數按照趙榮珍等研究中的參數選擇部分進行設置。本文所提方法以及各分類算法對于平衡數據、不平衡數據的識別結果見表3。

圖9 3 000 r/min下的各類數據云圖Fig.9 Various data cloud images under 3 000 r/min

表3 各分類方法識別精度Tab.3 Identification accuracy of classification methods

由表3分析可知,各算法對于平衡數據識別精度較高,BP神經網絡、SVM、E-ELM的各狀態平均識別精度分別為98.45%,97%,98.5%,而本文所提DCT方法的平均識別精度為100%,具有一定的優勢。對于不平衡數據,BP神經網絡、SVM、E-ELM的識別精度明顯下降,對于訓練樣本數量較少的保持架狀態識別精度下降至95.5%,77.5%,90%,平均識別精度也下降至97.3%,90%,96.5%,而本文所提DCT方法的平均識別精度保持不變。其原因在于,BP神經網絡對樣本的依賴性強,少數類樣本的識別精度會明顯低于多數類樣本識別精度;SVM在對不平衡數據實施分類時,分類面會向少數類樣本方向偏移,導致部分少數類樣本被誤判為多數類樣本,降低了少數類樣本的識別精度;E-ELM中所使用的ELM分類器與BP神經網絡都基于前饋神經網絡的架構之下,故ELM同樣對樣本有著較強的依賴性,少數類樣本識別精度相對較低。而DCT方法通過逆向、正向云發生器將多數類、少數類訓練樣本繪制成云滴點數相同的云圖,并依據云圖間距離對測試樣本進行分類,一定程度減少了不平衡數據對于分類的影響。因此,相較于其他分類算法,本文方法在不平衡數據故障識別中具有一定優勢。

3.3 不同工況下DCT方法的識別效果

對于某一轉速的軸承不平衡故障數據,所提方法有著較高的識別精度,然而旋轉機械的軸承在實際運轉情況下轉速不唯一,甚至會出現速度波動的干擾。因此,需要對本文方法對于不同工況的適用性以及抗速度波動干擾能力進行驗證。為驗證本文所提的DCT方法在不同轉速下的識別能力,分別采用3.1節介紹的2 600 r/min,2 800 r/min,3 000 r/min這3種工況下的軸承故障數據進行識別。為驗證本文方法的抗速度波動干擾能力,將兩種轉速下的軸承故障初始特征集按照5∶5進行混合,2 600 r/min,2 800 r/min初始特征集合并成維數相等的2 600~2 800 r/min初始特征集,2 800 r/min,3 000 r/min初始特征集合并成維數相等的2 800~3 000 r/min初始特征集,并通過DCT方法對兩個混合轉速的初始特征集進行識別。將2 600 r/min,2 800 r/min,3 000 r/min,2 600~2 800 r/min,2 800~3 000 r/min 故障初始特征集的轉速命名為轉速A、B、C、D、E,不平衡數據構建以及試驗參數選擇按照3.1節所述方法進行。DCT方法對于不同工況的不平衡數據識別結果,如表4所示。

表4 DCT方法對于不同工況下的不平衡數據識別結果Tab.4 DCT method for the identification results of unbalanced data under different working conditions

由表4可以看出,DCT方法對于不同工況下的不平衡數據仍有較高模式識別精度并有良好的抗速度波動干擾能力,故本文所提方法在軸承故障不平衡數據的模式識別領域具有較好的適用性。

3.4 DCT方法的泛化性能驗證

DCT方法在本實驗室所測的軸承數據集上對樣本分布不平衡的少數類樣本能夠取得較高的識別精度,而該方法是否適用于其他數據集需要進一步驗證。為驗證本文所提方法的泛化性能,采用美國凱斯西儲大學電氣工程實驗室的滾動軸承數據進行試驗驗證。所使用數據采集于機械系統驅動端的加速度傳感器,負載為1 hp、采樣頻率為12 kHz,故障損傷直徑為0.533 4 mm。數據集中包含了軸承運轉的正常狀態以及3種故障狀態:內圈故障、外圈故障、滾動體故障。對各個狀態振動信號進行小波消噪,并截取各類振動信號80組,共計320組樣本,樣本長度為2 048點。按照3.1節所述方法構建特征集,選取權重較大的前10個特征構成低維特征集,不平衡訓練集組成情況為:正常狀態樣本40個、內圈故障樣本20個、外圈故障樣本20個、滾動體故障樣本20個。由DCT方法得到的識別結果如圖10所示。

圖10 不平衡數據集識別結果圖Fig.10 Unbalanced data set identification result graph

由圖10分析得出,對于此不平衡數據集,內圈狀態的識別精度為97.5%,其中一個樣本被誤分為正常狀態,其余狀態的識別精度均能達到100%。由此可見,本文所提的DCT不平衡數據模式識別方法具有較好的泛化性。

4 結 論

針對不平衡數據分類時少樣本數據識別精度較低的問題,提出了一種基于云模型的不平衡數據模式識別方法。該方法以RelifF算法對特征集進行二次選擇,精簡分類器結構、減少計算量;通過云模型中的逆向、正向云發生器將多數類樣本以及少數類樣本繪制成云滴點數相同的各類數據云圖,依據測試數據云圖與訓練數據云圖間的距離判別測試數據類型,有效避免了不平衡數據對于識別過程的影響;結合集成學習方法,將各個特征的云分類器進行整合,解決了云分類器僅處理單一特征的問題,提高了識別精度。

研究結果表明,本文所建立的方法相較于其他模式識別方法對旋轉機械軸承的不平衡數據有著較高識別精度,同時本方法適用于不同工況數據以及公開數據集,具有一定的泛化性,為旋轉機械的智能故障識別方法提供了一種新的參考方案。

猜你喜歡
模式識別分類特征
分類算一算
如何表達“特征”
不忠誠的四個特征
當代陜西(2019年10期)2019-06-03 10:12:04
分類討論求坐標
數據分析中的分類討論
教你一招:數的分類
抓住特征巧觀察
淺談模式識別在圖像識別中的應用
電子測試(2017年23期)2017-04-04 05:06:50
第四屆亞洲模式識別會議
第3屆亞洲模式識別會議
主站蜘蛛池模板: 欧美日本视频在线观看| 伊人网址在线| AV熟女乱| 999福利激情视频| 中文字幕久久亚洲一区| 国产精品成| 国产精品自在自线免费观看| 欧美性猛交一区二区三区| 国产在线精彩视频二区| 欧美福利在线播放| 中文字幕久久波多野结衣 | 人人91人人澡人人妻人人爽| 丁香婷婷综合激情| 精品人妻无码中字系列| 国产精品999在线| 午夜欧美理论2019理论| 人妻无码AⅤ中文字| 亚洲AV免费一区二区三区| 国产成人在线无码免费视频| 久久国产精品电影| 伊人大杳蕉中文无码| 狠狠色综合久久狠狠色综合| 91丨九色丨首页在线播放| 欧洲高清无码在线| 99热这里只有精品国产99| 日本欧美一二三区色视频| aa级毛片毛片免费观看久| 国产毛片基地| h视频在线播放| 亚洲欧美在线综合一区二区三区 | 波多野结衣在线se| 国产综合色在线视频播放线视| 成人一级黄色毛片| 午夜国产大片免费观看| 亚洲中文字幕97久久精品少妇| 亚亚洲乱码一二三四区| 国产欧美亚洲精品第3页在线| 女人av社区男人的天堂| 亚亚洲乱码一二三四区| 丁香五月亚洲综合在线| 婷婷综合亚洲| 亚洲美女AV免费一区| 国产成人1024精品| 色噜噜在线观看| 免费在线成人网| 国产色婷婷| 国产精品成人AⅤ在线一二三四 | 亚洲视频免| 色偷偷av男人的天堂不卡| 国内精自线i品一区202| 国产一线在线| 欧美国产日本高清不卡| 亚洲精品国产日韩无码AV永久免费网 | 国产真实乱子伦视频播放| 欧美亚洲香蕉| 国产亚洲一区二区三区在线| 国产 在线视频无码| 成年人免费国产视频| 亚洲第一黄色网址| 人妻少妇久久久久久97人妻| 伊人久久久久久久久久| 欧美α片免费观看| 亚洲美女久久| 亚洲最新地址| 国产一二视频| 免费a在线观看播放| 99ri精品视频在线观看播放| 亚洲人成影院在线观看| 国产欧美又粗又猛又爽老| 欧美日韩午夜| 国产精品成人一区二区不卡| 九色在线视频导航91| 国产无遮挡猛进猛出免费软件| 成人免费网站久久久| 国产精品永久免费嫩草研究院| 日韩在线网址| 亚洲嫩模喷白浆| 91久久国产综合精品| 国产无码制服丝袜| 亚洲男人的天堂久久香蕉| 精品国产美女福到在线不卡f| 亚洲人成网站18禁动漫无码 |