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

采用參數化解調的變轉速下柱塞泵故障診斷方法

2021-10-11 02:34:10徐孜潮群高浩寒陶建峰劉成良孟成文
西安交通大學學報 2021年10期
關鍵詞:信號方法模型

徐孜,潮群,高浩寒,陶建峰,2,劉成良,2,孟成文

(1.上海交通大學機械與動力工程學院,200240,上海;2.上海交通大學機械系統與振動國家重點實驗室,200240,上海)

柱塞泵作為液壓系統的動力源和關鍵元件,具有功率密度大、效率高、變量調節方便等優點,因此在航空航天、機器人、工程機械等領域得到了廣泛應用[1]。在航空航天領域,電靜液作動器(Electro-Hydrostatic Actuator,EHA)取消了集中式液壓能源系統,提高了飛機的安全性、可維護性和工作效率,是多電飛行器的關鍵部件之一[2]。作為EHA核心元件之一的軸向柱塞泵,需要在10 000 r/min范圍內變速工作[3],高轉速下柱塞泵容易出現空化現象,引發氣蝕,破壞元件,產生強烈振動,造成泵出口流量降低[4],導致效率下降、壓力波動、設備故障,甚至造成災難性后果。因此,在變轉速工況下,對柱塞泵空化故障進行及時診斷,可以為液壓系統的維護帶來極大的方便,節省大量的時間和成本,提高液壓系統的可靠性與安全性。

常用的故障診斷方法有3類:基于模型的方法、基于知識的方法和數據驅動的方法[5]?;谀P偷姆椒ㄒ蠼⑤^為準確的數學模型,由于柱塞泵結構復雜且內部多物理場耦合,因此很難建立準確的數學模型?;谥R的方法需要大量專家經驗知識作為基礎,但是泵的動態復雜性使得專家經驗知識難以獲得。數據驅動的方法由于僅需要歷史數據受到越來越多的關注,人工神經網絡[6]、支持向量機[7]等方法均在泵空化診斷上有成功應用。近年來,深度學習迅速發展,因其強大的自動特征提取能力獲得研究者的青睞。Xiao等將小波包分解的特征輸入長短期記憶網絡(LSTM)[8],可準確識別柱塞泵的健康水平,相比支持向量機有更高的準確率。Chao等將三通道振動信號轉換成RGB圖像,輸入二維卷積神經網絡(CNN)實現柱塞泵空化等級識別[9],具有良好的抗噪性能。

然而,目前多數柱塞泵故障診斷的研究都基于恒定工況的假設,對變轉速條件下柱塞泵的故障特征缺乏深入了解,限制了在實際中的應用。變轉速工況下信號是非平穩的,其統計特征隨時間變化,使得常用的信號處理方法如傅里葉變換、帶通濾波等技術都不再適用。基于階次跟蹤的方法是公認有效的變轉速下信號處理方法,其核心思想是根據轉速曲線在角域內對信號進行重采樣,將時變信號轉換為角域內的周期平穩信號,使得恒定工況下的信號處理、故障診斷技術可以延用。傳統的硬件式階次跟蹤通過硬件直接獲取角域信號,成本高且有安裝空間要求。計算式階次跟蹤通過插值算法進行重采樣,簡化了硬件設備,在變轉速故障診斷上取得了一定的效果[10],但是存在3個問題:一是仍然依賴轉速測量設備,鑒于航空領域對于減小重量、體積的極端要求,額外的硬件安裝限制了該方法的應用;二是角域重采樣過程中的角加速度恒定假設與插值算法會帶來誤差,尤其在轉速波動大的時候誤差無法控制;三是故障信號通常含有多個分量,但階次跟蹤無法分離故障分量[11]。針對第一個問題,研究者提出無轉速計式頻率估計方法[12],基于時頻變換提取瞬時頻率曲線。Wang等對軸承振動信號進行短時傅里葉變換,使用峰度搜索從時頻圖中提取頻率曲線[13],再做角域重采樣,提取出軸承故障特征。針對第二個問題,Olhede等提出廣義解調算法[14],利用頻率曲線構造解調因子,解調變換將信號在時頻域內轉換成平行與時間軸的直線,無需重采樣即可實現時變信號的平穩化。Ma等將自適應線調頻模態分解算法與廣義解調相結合[15],增強了算法在變轉速滾動軸承故障診斷中的自適應能力和抗噪聲能力。針對第三個問題,傳統信號分解方法如經驗模態分解存在模態混疊等缺陷,研究者提出迭代廣義解調算法。Feng等利用迭代廣義解調精確估計其振幅包絡和瞬時頻率[16],從而構造無交叉項干擾的精細時頻分辨率的時頻表示,在非平穩條件下行星齒輪箱故障振動信號分析中取得較好的應用成果。變轉速下的故障診斷方法,瞬時頻率曲線的準確估計仍是其核心,目前基于時頻分析的方法非常依賴于時頻變換的精度,且時頻變換過程中需計算很大的參數矩陣,導致無法處理數據量大的信號[17]。因此,Yang等提出參數化解調方法[18],基于頻譜集中性指標和優化算法實現瞬時頻率曲線的準確估計,結合迭代解調變換,實現時變多分量信號的分解。

針對變轉速工況下柱塞泵故障診斷問題,本文在前人研究的基礎上做出改進,提出一種基于參數化解調結合1DCNN-LSTM網絡的柱塞泵空化故障程度識別方法,結合了參數化解調的非平穩信號處理能力與深度學習的特征提取能力。該方法利用參數化解調對變轉速下泵出口壓力信號自適應地提取重點分量并重構,重構信號輸入CNN-LSTM網絡提取局域特征并學習長期時序信息。本文使用Pumplinx軟件建立計算流體動力學(CFD)仿真模型,進行仿真實驗來獲取變轉速下柱塞泵出口壓力信號,仿真實驗驗證了參數化解調提取信號分量的有效性;采集實測信號進行驗證,表明該方法在強噪聲環境下仍然能以較高的精度對柱塞泵空化程度進行識別,且該方法有效抑制過擬合,具有良好的泛化性能。

1 診斷方法原理

1.1 參數化解調

變轉速工況下設備信號一般是非平穩多分量的,其統計特征隨時間變化。參數化解調算法是一種處理時變信號的方法,能夠分離提取各信號分量。其方法是通過參數估計構建解調因子,將時變信號轉換為頻率恒定的平穩信號,在使用帶通濾波過濾后反解調獲得分離出的信號分量,多次迭代,即可獲得多個信號分量。

本文用多項式信號模型來描述非平穩信號,多項式相位信號在旋轉機械、雷達探測等領域有廣泛應用[19],其相位是時間的多項式函數,定義為

(1)

式中:M為信號分量個數,可以事先通過時頻分析如STFT來判斷;xm(t)為觀測到的第m個信號分量;ε(t)為復雜噪聲;Am為第m個分量的振幅;P為多項式相位的階數;am,k(k=0,…,P)為第m個分量的相位系數。

使用參數化解調方法,對非平穩多分量信號進行重點分量提取并重構,算法總共分為4個步驟,流程如圖1所示。

圖1 基于參數化解調的信號分解重構方法流程Fig.1 Flow chart of signal decomposition and reconstruction method based on parameterized demodulation

(2)

(a)原信號FFT

(b)解調后信號FFT圖2 解調變換對頻譜的影響Fig.2 Effect of demodulation transform on frequency domain

因此,構建頻譜集中度(SCM)來量化解調后分量的集中程度,作為參數估計的指標,其定義為

DSCM=∑(|Zd(f)|q)

(3)

式中:Zd(f)是式(2)中xd(t)的傅里葉變換;q可以是任何正整數,更大的q在一定程度上能夠增益目標信號分量、抑制噪聲,從而減小優化算法陷入局部最小值的可能性,但是更大的值并不一定會顯著改善搜索結果,本文中取q=3。

(4)

在處理含有多個分量的信號時,能量高的分量將優先被提取。

1.1.2 解調變換 利用第1步中估計得的第i個分量的相位系數,按式(2)對信號進行解調變換。目標分量將被變換為一個平穩分量,集中在某一特定頻率上,在時頻圖上,該分量是一條與時間軸平行的直線。而其他分量會分布在較寬的帶寬中,噪聲則均勻地分散在所有頻段。

1.1.3 對解調信號進行帶通濾波 將解調信號頻譜峰值所在頻率作為濾波器的中心頻率

(5)

濾波后獲得的信號為

e(t)≈Aiexp(j(ai,0+ai,1t))+e(t)

(6)

式中:e(t)為噪聲和小部分不需要的分量,其能量相比目標分量很低;由于分量的系數被正確地估計了,Δai.k非常接近0。

1.1.4 恢復過濾后的分量 使用估計的系數來恢復過濾后的分量

(7)

從原信號中去除該分量得到信號殘差,按下式計算

(8)

迭代地執行上述4個步驟,直到獲得所有目標分量。

1.2 1DCNN-LSTM網絡

采用深度學習中較為常用的1D-CNN和LSTM網絡來搭建分類器。1D-CNN網絡能有效地從連續數據中提取特征,被廣泛應用于早期診斷、異常檢測等領域[20]。LSTM網絡同樣是處理時間序列的強大工具,能夠學習長期依賴關系。

故障診斷模型結構如圖3所示。原始數據經參數化解調處理后輸入到1DCNN-LSTM網絡中,該網絡主要包括卷積層、池化層、LSTM層、全連接層和輸出層。

卷積層利用權值共享的卷積核對數據的局部區域進行卷積運算,實現對輸入信號的局部特征提取,卷積器的運算公式為

yl+1,k=g(bl,k+Wl,k*Xl)

(9)

圖3 故障診斷模型結構Fig.3 Structure of the fault diagnosis model

式中:yl+1,k為第l+1層第k個卷積核卷積后的輸出;g為激活函數,此處采用線性整流函數(ReLU),增強網絡的非線性表達能力;Wl,k為第l層第k個卷積核的權重矩陣,bl,k為其偏置項;Xl為第l層的輸入;*為卷積運算。

池化層又稱下采樣層,可以減少數據處理量同時保留充足的有效信息,抑制模型過擬合。本模型采用最大池化方法,將池化區域內的最大值作為池化后對應神經元的值

(10)

式中:pl+1,k(m)為池化后第l+1層第k個特征圖第m個神經元的值;yl,k(i)為第l層第k個特征圖第i個神經元的值,i∈[(m-1)L,mL],L為池化核尺寸。

LSTM是一種特殊的循環神經網絡(RNN),能夠學習長期依賴關系,又解決了傳統RNN在面對大數據量時的梯度消失或梯度爆炸的問題[21]。圖4是單個LSTM單元的結構,主要由遺忘門、輸入門、輸出門組成,遺忘門決定是否保留前一單元的狀態,輸入門根據前一時刻狀態和此刻輸入來更新該單元的狀態,輸出門控制單元的輸出。

LSTM單元的更新如公式(11)所示

(11)

?—向量積;?—向量加;σ—Sigmoid函數;tanh—雙曲正切函數;ft—遺忘門;it—輸入門;ot—輸出門;xt—單元輸入;ht—單元輸出;ct—狀態向量;候選狀態向量t。圖4 LSTM單元內部結構Fig.4 Internal structure of LSTM unit

式中:Wf、Wi、Wo分別為遺忘門、輸入門、輸出門的權重矩陣;bf、bi、bo分別為遺忘門、輸入門、輸出門的偏置項。

全連接層與輸出層位于模型的最后,將學習到的特征映射到樣本標簽空間,輸出分類結果。

為提高診斷模型泛化能力,引入隨機丟棄機制,在池化層之后插入dropout層,dropout層隨機將一部分神經元激活值置0,減少過擬合現象。

訓練模型時,優化器采用Adamt[22]來更新網絡參數;損失函數采用稀疏多分類交叉熵,衡量真實標簽與模型預測標簽的相似性,其計算公式為

(12)

2 故障診斷流程

2.1 數據預處理

(1)使用參數化解調處理信號,按1.1節中的方法,提取泵出口壓力信號中能量高的前3個分量并重構成一維壓力信號。

(2)采用滑動窗口的形式對信號進行切片,要求切片后一個片段數據中至少包含一個轉動周期的信息。單個片段中含數據點1 024個,相鄰片段重合點數512個。對每個樣本數據進行標準化,公式如下

(13)

(3)將數據集按照80%和20%的比例劃分為測試集與訓練集。

2.2 1DCNN-LSTM模型建立

按照圖3的結構建立深度學習模型,通過一維卷積層提取數據局部特征,再輸入LSTM層學習長期時序信息,添加池化層與dropout層抑制過擬合提高泛化能力。表1給出了深度學習模型主要網絡層的特征及參數。

表1 深度學習模型主要網絡特征及參數Table 1 The main network layer’s characteristics and parameters of the deep learning model

3 仿真驗證

為了檢驗參數化解調能否有效處理多分量時變信號以應用于泵空化診斷,本節建立CFD仿真模型,通過仿真獲取變轉速工況下柱塞泵發生空化時出口動態壓力信號進行驗證,也為后續空化模擬實驗提供參照。采用Pumplinx仿真軟件對柱塞泵進行CFD仿真,仿真流程如下。

(1)使用三維建模軟件CREO提取被測泵的流體域模型導入Pumplinx,如圖5a所示,通過添加油膜來模擬柱塞泵中配流副、柱塞副、滑靴副的間隙,其中配流副和滑靴副油膜為5 μm,柱塞副油膜為10 μm。

(a)流體域

(b)網格劃分圖5 CFD模型流體域與網格Fig.5 Fluid domain and generated mesh for the CFD model

(2)對模型進行網格劃分,生成的網格模型中包含292 707個網格單元,如圖5b所示。

(3)添加空化模型。仿真中采用的空化模型在Singhal[23]提出的全空化模型進行了擴展,添加基于道爾頓-亨利定律的平衡溶解氣體模型[24],該模型能夠較好地模擬軸向柱塞泵流體特性[25]。

(4)設置仿真參數。表2列出了仿真模型參數,壓力為絕對壓力,轉速在6 000~10 000 r/min的范圍內變化。

表2 仿真模型參數Table 2 Parameters for the simulation model

由于柱塞泵的實際工作環境常常比較復雜,工作條件惡劣,泵的信號難以避免被環境噪聲所污染。針對此種情況,將高斯白噪聲疊加在原始信號上,構造噪聲信號,信噪比(SNR,用符號RSN表示)為0 dB。信噪比表征信號的污染程度,信噪比越低,污染越嚴重,其計算公式為

(14)

式中:Psignal和Pnoise分別為信號功率及噪聲功率。

對加噪聲信號按1.1節中的流程進行參數化解調處理,提取能量最高的3個分量并重構。圖6a~6c展示了原信號、加噪信號與重構信號的時域波形圖,可見處理后的重構信號比加噪信號更平滑,與原信號具有更高的相似性。圖6d~6f展示了原信號、加噪信號與重構信號的STFT時頻圖,原信號中能量較高(顏色較淺)的3個非平穩分量均被提取到重構信號中,且加噪信號中遍布所有頻帶的噪聲均在處理后被去除,有效分量十分清晰。

(a)原信號時域波形

(b)加噪信號時域波形

(c)重構信號時域波形

(d)原信號時頻圖

(e)加噪信號時頻圖

(f)重構信號時頻圖圖6 參數化解調處理前后的壓力信號Fig.6 Pressure signals before and after parameterized demodulation

仿真結果表明,在噪聲環境下,參數化解調方法能夠成功地提取出非平穩信號的分量,對噪聲有顯著的抑制效果。

4 實驗驗證

4.1 實驗方案與數據采集

為了監測軸向柱塞泵在變轉速工況下的空化情況,搭建了柱塞泵故障模擬試驗臺,系統由主工作回路(圖7a)與冷卻過濾回路(圖7b)組成,通過油箱連接,實物圖如圖7c所示。主工作回路中,通過離心泵增壓,結合壓力傳感器調節被測柱塞泵進口壓力,電機通過聯軸器帶動被測泵在不同的轉速下運行;電磁比例溢流閥調節出口壓力,模擬負載。冷卻過濾回路中熱交換器將油溫穩定在40 ℃左右,以減少油溫波動對柱塞泵工作狀態的影響。

(a)工作回路原理圖

(b)冷卻過濾回路原理圖

(c)試驗臺實物圖圖7 柱塞泵故障試驗臺實物圖Fig.7 Test rig for piston pump fault testing

實驗臺安裝了9個主要的傳感器:壓力傳感器PS1-PS3測量泵入口、出口、回油壓力,采樣率為102 400 Hz;流量計FS1-FS2測量泵出口、回油口流量,采樣率為20 Hz;轉速計RS1測量被測泵的轉速,采樣率為102 400 Hz;溫度傳感器TS1-TS3測量泵入口、出口、回油口的油液溫度,采樣率為20 Hz。本文主要使用出口壓力信號對柱塞泵空化故障進行診斷,由壓力傳感器PS2采集。

系統運行時,傳感器連接到NI數據采集卡,通過以太網傳輸至計算機上存儲并進行進一步分析。使用測控軟件Labview定義被測泵入口壓力、出口負載、轉速來配置不同的工況,并通過PLC執行。實驗工況設定泵出口壓力為表壓5 MPa,入口壓力為表壓0.1~0.3 MPa。入口壓力和油液溶解氣體體積分數是影響空化的主要參數[25],實驗中維持油液溫度穩定,減少其對溶解氣體體積分數的影響,通過調節泵入口壓力來模擬空化故障的發生;柱塞泵在3 000~10 000 r/min的范圍內運行,轉速變化如圖8所示,變速段為有效數據。

圖8 柱塞泵轉速曲線Fig.8 Speed curve of piston pump

實驗中,由于泵殼體不透明,無法通過觀察泵內部氣泡來判斷空化發生情況。考慮到空化會引起泵出口體積流量下降,且流量下降程度與空化嚴重程度正相關[26],本文使用容積效率損失Δη來間接表征柱塞泵的空化程度,容積效率損失越大,空化程度越高。容積效率損失的計算公式如下

(15)

式中:qo為泵出口流量;qr為泵回油流量;Q為泵理論吸油流量;n為轉速;V為泵的排量。由于流量計采樣率較低,要求出口、回油流量qo、qr實現同步采樣,非采樣點處的Δη通過相鄰采樣點的Δη來估計范圍,若相鄰采樣點的Δη相差過大,則該段數據舍棄不用。

如表3所示,按容積效率損失分為3個空化等級[26]:無空化0%~1.5%,輕微空化2%~3%,中等空化7%~8%。在實驗中,進一步降低入口壓力時,泵將發生嚴重空化,會對被測泵造成不可逆的破壞,在實際應用中,希望能在空化程度較輕的時候就能分辨出來,以便及時采取措施抑制空化發生。

表3 不同工況下的柱塞泵空化等級劃分Table 3 Cavitation degrees of piston pump under different working conditions

按第2節的流程處理數據、搭建診斷模型并進行診斷。數據集按照80%訓練集、20%測試集的比例隨機劃分,劃分情況如表4所示。

表4 訓練集與測試集在不同空化等級下的樣本數量Table 4 Sample sizes of training set and testing set under different cavitation degrees

4.2 實驗結果分析

4.2.1 參數化解調前處理的影響 一組原始信號經過了參數化解調前處理,另一組信號未經解調重構,分別將兩組數據作為輸入,用結構、參數完全相同的模型進行訓練。圖9給出了2種方法在測試集上的分類結果混淆矩陣,每列為各空化等級的預測樣本占總樣本之比,每行為各空化等級的真實樣本占總樣本之比,對角的數值越接近1,說明分類準確度越高。從圖9看出,采用參數化解調方法能夠很好地區分3種不同空化等級,而不采用該方法在無空化和輕微空化的分類上表現不佳,有較明顯的誤分類。為了進一步評價該方法的分類效果,采用準確率A、精確度P、回收率R及F1值這4個指標,其計算公式分別為

(16)

(17)

(18)

(19)

(a)參數化解調+CNN-LSTM

(b)CNN-LSTM圖9 使用兩種方法的診斷結果混淆矩陣Fig.9 Confusion matrixes of diagnosis results using two methods

式中:TP、FP、FN、TN是測試集中不同樣本的數量,TP是真正例,FP是假正例,FN是假負例,TN是真負例。準確率A是預測正確的樣本占總樣本的比例;精確度P是衡量實際正例占預測正例的比例;回收率R指實際正例被預測為正的概率;F1是精確度與回收率的調和平均,綜合評價診斷模型的性能。

表5對比了是/否使用參數化解調在訓練集和測試集上的分類性能評價指標,可見參數化解調方法在測試集上的各項指標都比較高,其分類準確率達到95.4%,比未解調的方法準確率高了6.3%,說明該方法在變轉速工況下能有效地提高柱塞泵空化故障等級的識別能力。兩種方法在訓練集上的準確率很相近,但是未解調的方法在測試集上的準確率下降得更多,說明對變工況下非平穩信號進行重點分量提取并重構,能夠有效抑制過擬合,提高模型的泛化能力。

表5 兩種方法的診斷性能對比Table 5 A comparison between the diagnostic performances of two methods

4.2.2 在噪聲數據集上的影響 在實際應用中,柱塞泵工作環境惡劣,信號不可避免會被噪聲污染。為了評估診斷模型在噪聲環境下的變轉速工況空化等級識別能力,構造一組含噪聲的信號,作為模型測試集。構造方法是將不同信噪比的高斯白噪聲疊加到原始壓力信號上,信噪比范圍是0~10 dB。實驗方法是,按照上一節的方法對原始壓力信號進行參數化解調預處理后訓練深度學習模型,然后對含噪聲信號同樣進行參數化解調預處理后進行測試。圖10展示了是否進行參數化解調對噪聲環境下分類能力影響的對比。在中等空化程度下,兩種方法的診斷能力都很高,參數化解調處理對于噪聲環境下的診斷能力沒有明顯的提高;在無空化與輕微空化條件下,不經參數化解調處理,隨著噪聲增強,診斷準確率顯著下降,在4 dB環境下無空化識別準確率為67.6%,輕微空化識別準確率僅為43.4%,無法滿足診斷需求,而參數化解調方法在噪聲環境下有較強的魯棒性,0 dB的強噪聲下準確率維持在90%以上,說明該方法能在空化初生時就實現有效識別,便于及時采取措施抑制空化。

(a)無空化

(b)輕微空化

(c)中等空化圖10 兩種方法在噪聲數據集上的診斷性能對比Fig.10 A comparison between the diagnostic performances of two methods on noise dataset

5 結 論

針對變轉速工況下軸向柱塞泵故障診斷問題,本文提出一種基于參數化解調和深度學習的診斷方法。柱塞泵出口動態壓力信號經過參數化解調處理,提取重點分量并重構,輸入到CNN-LSTM網絡訓練,用于空化故障等級的識別。此外,將不同信噪比的高斯白噪聲疊加在原始壓力信號上構造噪聲信號,經由參數化解調處理后作為測試集對訓練好的CNN-LSTM模型進行評估。實驗結果表明,本文方法能夠自動提取信號的重點分量,準確地對柱塞泵空化嚴重等級進行判斷,準確率達到95.4%,相比未經參數化解調的方法,本文方法能抑制過擬合,提高泛化能力,將分類準確率提高了6.5%。在噪聲環境下,本文方法具有較強的魯棒性,0 dB的強噪聲下準確率能維持在90%以上,而不采用參數化解調的方法,其準確率隨噪聲增強顯著下降。

猜你喜歡
信號方法模型
一半模型
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
重要模型『一線三等角』
完形填空二則
重尾非線性自回歸模型自加權M-估計的漸近分布
基于FPGA的多功能信號發生器的設計
電子制作(2018年11期)2018-08-04 03:25:42
3D打印中的模型分割與打包
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
基于LabVIEW的力加載信號采集與PID控制
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
主站蜘蛛池模板: 亚洲香蕉久久| 波多野结衣中文字幕一区二区| 丁香五月激情图片| 国产激情无码一区二区三区免费| 99久久国产综合精品2020| 国产免费怡红院视频| 国产成人精品高清在线| 婷婷亚洲最大| 日韩精品一区二区三区大桥未久| 国产福利微拍精品一区二区| 午夜视频www| 久久精品日日躁夜夜躁欧美| 国产流白浆视频| 精品精品国产高清A毛片| A级全黄试看30分钟小视频| 国产无码在线调教| 666精品国产精品亚洲| 日本久久网站| 国产丝袜一区二区三区视频免下载| 国产午夜福利亚洲第一| 国产成人一区| 麻豆精品久久久久久久99蜜桃| 农村乱人伦一区二区| 国产精品一区二区国产主播| 亚卅精品无码久久毛片乌克兰| 国产成人凹凸视频在线| 99热这里只有精品国产99| 国产精鲁鲁网在线视频| 成人亚洲国产| 免费一极毛片| 爆操波多野结衣| 女人18毛片久久| 欧美不卡二区| 亚洲综合天堂网| 97无码免费人妻超级碰碰碰| 久久6免费视频| 国产va免费精品| 无码网站免费观看| 国产va在线| 亚洲欧洲自拍拍偷午夜色无码| 91午夜福利在线观看精品| 亚洲一区二区约美女探花| 在线毛片免费| 日韩精品视频久久| 91精品福利自产拍在线观看| 中字无码av在线电影| 亚洲日韩高清无码| 国产网友愉拍精品| 青青草a国产免费观看| 国产精品白浆在线播放| 久久精品日日躁夜夜躁欧美| 亚洲综合九九| 毛片网站在线播放| 视频二区国产精品职场同事| 毛片网站在线看| 国产精品极品美女自在线网站| 日韩小视频网站hq| 久久久久免费看成人影片| av一区二区无码在线| 亚洲手机在线| 亚洲欧美色中文字幕| 欧美激情福利| 精品国产电影久久九九| 精品国产自| 五月婷婷精品| 欧美啪啪一区| av在线人妻熟妇| 亚洲国产第一区二区香蕉| 亚洲美女操| 日韩中文欧美| 久久综合九九亚洲一区| 97超级碰碰碰碰精品| 久青草免费视频| 在线无码av一区二区三区| 精品国产aⅴ一区二区三区| 一本一道波多野结衣一区二区| 国模极品一区二区三区| 黄色网页在线观看| 天天色综合4| 国产成人精品免费视频大全五级| 国产区免费精品视频| 国产乱码精品一区二区三区中文 |