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

基于NLM-CEEMDAN 和樣本熵的水電機組振動信號去噪

2023-07-04 00:36:38章芳情王成城
中國農村水利水電 2023年6期
關鍵詞:模態振動信號

章芳情,袁 方,賀 玉,王成城,郭 江

(1. 武漢大學動力與機械學院,湖北 武漢 430072; 2. 武漢理工大學自動化學院,湖北 武漢 430070;3. 機械工業儀器儀表綜合技術經濟研究所,北京 100050)

0 引 言

水電機組是實現能量轉換的高安全裝備,在國家實現雙碳戰略中發揮著重要作用。隨著經濟社會的發展,對水電機組安全性、可靠性和穩定性的要求越來越高,但水電機組運行環境相對比較復雜,這對其安全穩定運行造成了很大威脅,如何對機組運行狀態進行有效監測與診斷是目前行業研究熱點[1]。據相關研究表明,水電機組近80%的故障均可以從設備的振動中體現出來[2]。因此,通過對水電機組的振動信號進行在線監測診斷,能及時發現裝備的異常狀況,從而保障水電機組的安全穩定運行[3],為后續故障診斷工作的開展提供可靠依據。

水電機組的振動信號是機組運行狀態評價的重要指標之一[4],其本身就是一種非線性、非平穩性的信號,再加上信號采集時一般會受到設備運行所產生的噪聲影響,使得采樣信號難以直接真實反映設備運行狀態。因此,如何將設備振動信號去噪,進而獲取真實信息十分必要。傳統的傅里葉變換比較適合用于平穩規則的線性信號分析[5],小波分析雖然可以進行相應的非平穩信號分析,但是其參數設置缺乏自適應性[6]。經驗模態分解(Empirical Mode Decomposition ,EMD)作為一種較新的時頻分析算法,通過將信號分解成若干個固有模態分量(Intrinsic Mode Function ,IMF),可在時頻上對信號進行全面分析,但其分解過程中容易產生模態混疊、端點效應等問題[7]。集合經驗模態分解(Ensemble Empirical Mode Decomposition,EEMD)是EMD的一種改進方法[8],它雖然可以部分抑制EMD分解所產生的模態混疊和端點效應,但該方法本身也容易帶來輔助噪聲殘留等新問題,影響去噪效果[9,10]。為此,Yeh 和Huang 等[11]提出了互補集合經驗模態分解(Complementary Ensemble Empirical Mode Decomposition,CEEMD)方法,有效抑制了殘留輔助噪聲,但信號分解不完備、計算效率低等問題未能很好解決。可變模態分解(Variational Mode Decomposition,VMD)能很好地避免類EMD 方法中的端點效應問題[12],但需要事先設定分解的K值導致無法自適應分解,難以滿足水電機組振動信號在線自動監測分析的需求[13]。近年來,利用自適應噪聲完備集合經驗模態分解(Complementary Ensemble Empirical Mode Decomposition with Adaptive Noise, CEEMDAN)方法[14]對非平穩振動信號的降噪已經在球磨機筒體振動信號去噪[15]、腦電信號噪聲濾除[16]和軸承故障診斷[17]等領域有了初步應用,并取得了較好的進展,但在水電機組振動信號處理方面還有待推廣。

為彌補CEEMDAN方法對低信噪比信號去噪效果不理想的缺點以及傳統分量檢測指標近似熵依賴數據長度和計算效率低的不足,本文在已有的研究基礎上,提出了一種基于NLMCEEMDAN 和樣本熵的水電機組振動信號去噪新方法。首先對原始信號進行非局部均值(Non-Local Mean,NLM)預處理降噪[18],然后采用CEEMDAN 方法分解預處理信號得到一系列IMF分量,同時計算各個分量的樣本熵,并根據樣本熵值大小將分量劃分為高頻含噪分量、信噪混合分量和低頻有效分量。最后進一步通過小波閾值去噪方法濾出信噪混合分量中的噪聲成分,連同高頻含噪分量一起從原始振動信號中濾除,從而完成水電機組振動信號的去噪。以水電機組擺度信號監測分析為例進行應用,通過擺度仿真信號和擺度實測信號進行驗證,并與小波閾值去噪、CEEMDAN 濾波去噪作對比試驗,結果表明所提出的NLM-CEEMDAN 去噪方法能有效提取水電機組振動信號中的有效信息,是一種優異的信號處理方法。

1 基本原理

1.1 非局部均值濾波

非局部均值濾波算法是利用圖像中存在眾多相似結構的特性,通過對這些相似結構進行加權平均操作來達到圖像去噪的目的,由此NLM算法被顯著用于二維圖像去噪[19]。然而這些相似特性也存在于一維信號中,并且NLM 算法已被成功應用到一維的滾動軸承振動信號處理中[20],因此NLM算法也能被用于處理水電機組振動信號。

本文利用NLM 算法對水電機組振動信號進行預處理降噪,以方便后續信號特征提取。NLM算法原理如下:

假設實際帶噪聲的水電機組振動采集信號y(t)為真實振動信號x(t)和外部干擾噪聲n(t)的疊加,即:

NLM 算法通過計算y(t)中的全部相似塊的加權平均值K(t)來估算真實信號x(t),即:

式中:D(t)表示以t為中心的搜索范圍內所有點的集合;Z(t)為歸一化因子,表示全部搜索塊相似度之和,其計算公式為:

式中:ω(t,s)表示權重,指以t和s為中心的2個搜索塊之間的相似度,且須滿足0 ≤ω(t,s)/Z(t) ≤1 和∑tω(t,s) = 1 的基本條件,有:

式中:λ為濾波器帶寬參數,影響著去噪信號平滑度;Δ 為以t為中心的搜索塊,取K為Δ 區域長度的一半,參數K影響著NLM算法的計算量和計算時間;LΔ為以s為中心的鄰域塊,有LΔ=2P+ 1,參數P影響著算法運行過程中所發現的相似結構塊的數量。

參數λ,K,P為NLM算法的決定性參數,很大程度上影響著算法對一維信號的去噪效果,但這些參數的設置仍較多的依賴于人為經驗。

1.2 自適應噪聲完備集合經驗模態分解

EEMD 方法能部分抑制EMD 方法的模態混疊,但由于殘存輔助噪音未加隔離,使得噪聲會從高頻向低頻轉移,影響后續信號分析處理[21]。針對上述缺陷,CEEMDAN方法解決措施為:①在待分解信號中自適應地加入EMD 分解后含輔助噪音的IMF分量,解決了EEMD 分解過程中存在的完備性缺失和重構誤差過大的問題[22];②得到EMD 分解后的第一階IMF分量后就進行總體平均計算,有效解決了噪聲從高頻向低頻轉移的問題;③能減少很多EEMD 分解過程中的對于分析意義不大的幅值很小的低頻分量,具有更好的模態分解結果。

CEEMDAN分解算法原理步驟如下:

(1)在原信號X(t)里添加滿足N(0,1)的白噪音εωi(t),ε表示噪聲標準差,并將得到的信號進行第一階段的EMD 分解,得到1個第一階段的IMF分量。通過對X(t)加N種不同的白噪音分別進行第一階段的EMD 分解,得到N個固有模態分量IMFi1(t),再對這些固有模態分量集合進行總體平均,得到最終的第一階固有模態分量IMF1(t),即:

(2)計算原信號X(t)減去IMF1(t)的剩余量。

(3)定義Ek表示給定信號經EMD 分解后獲得的第k個模態分量。將剩余量X1(t)加入ε1E1(ωi[t])后繼續EMD 分解直至得到滿足IMF2(t)條件的分量,之后對所得到的N個分量作整體平均,得到原信號的第二階IMF分量,即:

(4)設最終分解次數為M,k= 2,3,4,…,M時,計算第k個余量:

(5)向Xk(t)中添加輔助噪音εkEk[ωi(t)]后進行第一階段EMD 分解獲取第一階的固有模態分量集合,計算出原信號的第k+ 1階固有模態分量:

(6)重復執行步驟(4)、(5),當剩余量的極值點數小于或等于2 時停止上述迭代步驟。最終的剩余量即殘差分量,用R(t)來表示。

(7)經CEEMDAN分解后的原信號表示成:

1.3 樣本熵原理

樣本熵(Sample Entropy, SampEn)是由Richman 等人提出的一種通過度量原始信號中產生新模式的概率大小來衡量時間序列復雜程度的方法[23],當時間序列越復雜,含有噪聲分量越多時,樣本熵的值就越大[24]。對于具有非線性特征的水電機組振動信號具有較好的特征提取作用。作為近似熵的改進算法,樣本熵的計算減少了對數據長度的依賴程度,并且具有更好的一致性,算法更為簡單,計算速度更快。

由樣本熵計算時間序列的算法步驟如下:

(1)將定長為N的時間序列X(i)按順序組成m維矢量,即:

式中:i= 1,2,…,N-m+ 1。

(2)計算矢量X(i)和其他向量X(j)間的距離:

式中:j= 1,2,…,N-m+ 1;k= 0,1,…,m- 1,i≠j。

計算其均值可得:

(4)將維數增加到m+ 1,同時重復步驟(2)、(3),可得:

(5)當時間序列X(i)的長度N為有限值時,可將樣本熵定義為:

由式(16)可知,參數m,r,N的取值對于樣本熵值的計算準確性相當重要。根據文獻[25]以及多次試驗,發現預先取m=2,r=(0.1~0.25)δ(δ為原始數據標準差)時,能使樣本熵的統計特性得到最有效的表達,結果更能真實的反映信號特征。

2 基于NLM-CEEMDAN 和樣本熵的水電機組振動信號去噪

由CEEMDAN 理論可知,通過CEEMDAN 分解所獲得的若干個固有模態分量的局部頻域是從高頻到低頻變化的,而原信號所含的噪音主要集中于前幾個高頻分量中。引入樣本熵算法對IMF分量集合進行計算,IMF分量的樣本熵值越大,說明該分量越不規則,含噪聲越多,反之則反。

基于NLM-CEEMDAN 和樣本熵的振動信號去噪方法的流程如圖1所示,具體步驟為:

圖1 基于NLM-CEEMDAN和樣本熵的信號去噪流程Fig.1 Process of signal denoising method based on NLM-CEEMDAN

(1)先利用NLM 算法對水電機組振動信號X(t)預處理降噪,后將預處理信號進行CEEMDAN 分解獲得一系列的固有模態分量IMFk(t)和一個趨勢分量R(t);

(2)經多次試驗,設置樣本熵算法的模式維數m= 2,相似容限系數r= 0.2δ,并分別計算各模態分量和趨勢分量的樣本熵值;

(3)根據樣本熵值的大小,將分量集合劃分為高頻噪聲分量、信噪混合分量和低頻有效分量三部分,并且進一步利用小波閾值去噪方法從信噪混合分量中濾出噪聲成分。

(4)最后從原始信號中濾除高頻噪聲分量和信噪混合分量中的噪聲成分即可完成水電機組振動信號的去噪。

3 仿真分析

3.1 仿真信號的構造

水電機組的擺度信號是機組的重要監測指標,為驗證基于NLM-CEEMDAN 和樣本熵的振動信號去噪方法的有效性,本文選擇水電機組擺度信號進行仿真分析。水電機組擺度主要受到機械激振和水力激振的影響,機械激振一般以中頻(1、2、3倍轉頻)為主,水力激振一般以低頻(0.2~0.45 倍轉頻)為主,故構造仿真模擬信號如下[26]:

式中:A1~A6分別為20、4.5、2.55、1.5、0.4、0.3 μm;f1~f6分別為1.25、2×1.25、3×1.25、4×1.25、0.2×1.25、0.3×1.25 Hz。采樣頻率取1 000 Hz,可得原始未加噪仿真信號如圖2 所示。對原始未加噪仿真信號疊加一個信噪比為20 dB 的高斯白噪音,得到的仿真加噪信號如圖3所示。

圖2 仿真未加噪信號Fig.2 Simulation without noise signal

圖3 仿真加噪信號Fig.3 Simulation noise signal

3.2 仿真信號的CEEMDAN分解

為對仿真加噪信號進行NLM 預處理降噪,在結合文獻[27,28]所提的NLM 參數設置方法的基礎上,進行多次試驗,發現設置NLM 參數λ= 0.3σ(σ為仿真加噪信號的標準差),K= 20,P= 12 時能較好的濾除信號中的背景噪聲,同時很好地保存信號中的有效信息,方便后續對預處理信號的處理。

對預處理信號進行CEEMDAN 分解時,需要對加入原信號的噪聲幅值系數ε、執行CEEMDAN 的總次數M和最大允許迭代次數MaxIter進行合理的設置。加入噪聲幅值與理想信號的標準差εn服從如下規律[29]:

ε值越小,則εn越小,即分解精度越高,但當ε過小時,可能不足以引起信號局部極值點的變化,從而使加入噪聲以改變信號的局部時間跨度失去意義[30]。而M和MaxIter越大則信號分解越完備,但當M和MaxIter過大時會導致效率過低,耗時增加,不利于信號實時監測分析。

經過多次試驗,設置參數ε=0.2,M=1 000,MaxIter=1 000 后對經NLM 預處理降噪后的信號進行CEEMDAN 分解,可以得到11個IMF分量和1個趨勢分量R,m 如圖4所示。計算這些分量的樣本熵,結果如表1所示。

表1 仿真分析時各分量的樣本熵Tab.1 Sample entropy of each component for simulation analysis

圖4 仿真信號的CEEMDAN分解結果Fig.4 CEEMDAN decomposition of the simulation signals

為了更方便看出仿真分析時各分量的樣本熵值的變化趨勢,作其趨勢分布如圖5所示。

圖5 仿真信號各分量樣本熵值Fig.5 Sample entropy of each component for simulation analysis

從圖5中能明顯看出IMF1、IMF2、IMF33個分量的樣本熵值遠大于其他分量,而中間的IMF4、IMF5、IMF6分量的樣本熵值則是大小居中,剩余分量的樣本熵值基本都很小。由此可將前3個分量IMF1、IMF2、IMF3劃分為高頻噪聲分量,中間3 個分量IMF4、IMF5、IMF6劃分為信噪混合分量,剩余分量則劃分為低頻有效分量。通對信噪混合分量進行小波閾值去噪,篩選出其中的噪聲成分,之后在仿真加噪信號中濾除高頻噪聲分量和篩選出來的噪聲成分,完成水電機組擺度仿真信號的去噪。

3.3 分量重構與濾波去噪方法的對比

分解分量重構是目前常用的獲取去噪信號的方法,該方法通過類EMD 算法將原始信號分解,將所得到的IMFS分量直接或者濾波降噪后重構來獲取去噪信號。但該方法易受類EMD算法分解過程中存在的模態混疊以及端點效應等問題的影響,使得去噪重構信號中仍存在一定量的噪聲。針對上述問題,本文提出一種新的濾波去噪方法獲取去噪信號,不同于以往利用類EMD 算法獲取有效信號成分,新方法是通過類EMD 算法將原始信號中的噪聲部分篩選出來,再從原始信號中濾除噪聲部分,進而可獲得噪聲成分更少的去噪信號。

為便于對不同方法的去噪性能進行定量對比,下面定義了信噪比SNR和均方根誤差RMSE[31],且信噪比越大,均方根誤差越小,去噪效果越好。

(1)信噪比SNR。

(2)均方根誤差RMSE。

式中:M為采樣點數;Xi為原始未加噪仿真信號;Yi為去噪后的信號。

下面利用CEEMDAN分解對分解分量重構和濾波去噪兩種方法進行性能對比,限于篇幅,此處僅給出去噪效果性能比較,結果如表2所示。

表2 兩種方法的去噪性能Tab.2 Denoising performance of two methods

從表2可以發現分解分量重構法的去噪性能指標明顯差于濾波去噪法,說明濾波去噪的方法所得的去噪信號更為接近原始未加噪仿真信號,其去噪效果更好。

3.4 不同去噪方法的對比分析

為了驗證本文所提出的NLM-CEEMDAN 和樣本熵去噪方法在水電機組擺度信號去噪中的有效性,分別采用小波閾值去噪方法、CEEMDAN 去噪方法和本文提出的NLM-CEEMDAN 去噪方法對仿真加噪信號進行去噪處理,并通過去噪波形圖以及去噪性能指標來對不同方法的去噪效果進行對比分析。

小波閾值去噪過程中的參數選取會對原始信號的去噪效果產生較大的影響。本文根據文獻[15]中的小波閾值去噪參數的確定方法,將仿真加噪信號的小波閾值去噪參數設置為:sym9 小波基函數、分解層數為6、heursure 閾值準則和硬閾值函數。去噪后的波形圖如圖6~8 所示,去噪效果性能指標如表3所示。

表3 不同方法的去噪性能Tab.3 Denoising performance of different methods

圖6 小波閾值去噪信號Fig.6 De-noised simulation signals by wavelet threshold

將圖6 中小波閾值去噪信號與圖2 中的仿真未加噪信號相比較,發現小波閾值去噪已基本能濾除仿真加噪信號中的大部分噪聲,但去噪信號波形中的毛刺現象較為明顯,去噪效果不太理想。

將圖7和圖8中的波形同圖2波形相比,可以發現兩種去噪信號都基本很好地實現了仿真加噪信號的去噪。但CEEMDAN濾波去噪信號波形曲線過于光滑,說明在舍棄分量的同時也丟失了部分有效信息。而NLM-CEEMDAN 去噪則是把CEEMDAN 濾波去噪方法所舍棄的信噪混合分量進行降噪處理,從而在保持波形光滑的同時,使得原始信號中的一些有效特征也得到相應保留。

圖7 CEEMDAN去噪信號Fig.7 De-noised simulation signals by CEEMDAN

圖8 NLM-CEEMDAN去噪信號Fig.8 De-noised simulation signals by combined method

根據去噪性能評判標準,對比表3 中的數據可知,NLMCEEMDAN 去噪方法的去噪效果最好,相較于小波閾值去噪方法和CEEMDAN 濾波去噪方法,信噪比SNR分別提高3.877 和0.722 5 dB,均方根誤差RMSE分別降低0.139 9和0.021 4。

綜上,通過對不同去噪方法的去噪信號的波形特征和去噪性能評價指標進行對比分析,可以發現NLM-CEEMDAN 聯合去噪的方法是明顯優于其他兩種方法。

4 實例分析

4.1 實例信號的采集

為了驗證本文所提出的NLM-CEEMDAN 和樣本熵的去噪方法在水電機組振動信號去噪過程中的有效性,進一步選取水電機組擺度監測實例進行分析,

水電機組擺度信號采集過程中的機組轉速為250 r/min,試驗的采樣頻率fs為2 048 Hz,截取采樣數據中的6 000 個點進行分析,以確保信號處理分析中的特征參數能全面真實地反映實際工況。

4.2 實例分析

水電機組上導擺度實測信號波形如圖9 所示,其中包含了大量毛刺。先對該實測信號進行NLM 預處理降噪,再對預處理信號進行CEEMDAN 分解,得到11 個IMF分量1 個趨勢分量。計算每個分量的樣本熵值,結果如表4 所示。同樣根據各分量樣本熵值大小將前3 個分量IMF1、IMF2、IMF3劃分為高頻噪聲分量,中間3 個分量IMF4、IMF5、IMF6劃分為信噪混合分量,剩余分量則劃分為有效信息分量。由小波閾值去噪算法篩選出信噪混合分量中的噪聲成分,最后從上導擺度實測信號里濾除高頻噪聲分量和篩選出的噪聲成分,完成水電機組上導擺度實測信號的去噪。從圖10 中可以看出上導擺度實測信號中的噪聲得到了有效的濾除。

表4 擺度信號各分量的樣本熵Tab.4 Sample entropy of each component for throw signals

圖9 上導擺度實測信號Fig.9 Real signals of upper guide bearing throw

圖10 去噪后的上導擺度信號Fig.10 De-noised upper guide bearing throw signals

由于無法從實測信號得到原始理想的擺度信號值,無法直接計算SNR和RMSE,下面定義了信號去噪前后的噪聲抑制比(Noise rejection ratio, NRR),來表征去噪后有效信號的突出程度,該值越大,則去噪后的有效信號越突出[32]。

3 種去噪方法處理擺度實測信號后的噪聲抑制比如表5 所示。從表5中可以看出NLM-CEEMDAN 去噪的NRR值最高,濾波效果最好,小波閾值去噪次之,CEEMDAN去噪效果最差。

表5 不同方法去噪后的擺度信號噪聲抑制比Tab.5 The noise rejection ratio of de-noised swing signals by different methods

水電機組下導、水導擺度實測信號波形分別如圖11、13 所示。對這兩種信號NLM 預處理降噪后進行CEEMDAN 分解,都分別得到11 個IMF分量和1 個趨勢分量。計算每個分量的樣本熵值,計算結果如表4 所示。根據樣本熵值,進行分量的劃分,依據同樣的方式完成水電機組下導、水導擺度實測信號的去噪,分別如圖12、14所示。

圖11 下導擺度實測信號Fig.11 Real signals of lower guide bearing throw

圖12 去噪后的下導擺度信號Fig.12 De-noised lower guide bearing throw signals

圖13 水導擺度實測信號Fig.13 Real signals of turbine guide bearing throw

圖14 去噪后的水導擺度信號Fig.14 De-noised turbine guide bearing throw signals

從圖中可以看出,基于NLM-CEEMDAN 和樣本熵的方法能夠有效地去除水電機組下導和水導擺度實測信號中含有的大量背景噪聲。同樣計算下導、水導擺度信號的3 種不同去噪方法的噪聲抑制比。從表5 中可以看出,NLM-CEEMDAN 去噪的方法去噪效果最好。

5 結 論

本文提出了一種基于NLM-CEEMDAN 和樣本熵的水電機組振動信號去噪方法。通過該方法分別對水電機組擺度仿真信號和3種實測信號進行去噪分析,并與小波閾值去噪、CEEMDAN 去噪的方法進行對比試驗,比較波形差異和去噪性能指標大小后,得到以下結論。

(1)相較于傳統的原始信號分解分量重構的方法,本文提出的從原始信號中直接濾除噪聲成分的方法使去噪后的振動信號去噪效果得到極大的提高,更接近于真實振動信號。

(2)基于NLM-CEEMDAN 和樣本熵的振動信號去噪方法對于仿真信號及實測信號的去噪效果都優于傳統方法,能有效濾除噪聲成分,為設備振動信號在線監測提供有力幫助。

(3)本文的研究成果也可廣泛應用到核電、化工等行業的其他高安全裝備振動信號去噪。

猜你喜歡
模態振動信號
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
基于FPGA的多功能信號發生器的設計
電子制作(2018年11期)2018-08-04 03:25:42
中立型Emden-Fowler微分方程的振動性
基于LabVIEW的力加載信號采集與PID控制
國內多模態教學研究回顧與展望
基于HHT和Prony算法的電力系統低頻振蕩模態識別
UF6振動激發態分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
主站蜘蛛池模板: 国产剧情无码视频在线观看| 国产乱子精品一区二区在线观看| 国产欧美日韩另类| 在线国产91| 香蕉99国内自产自拍视频| 中文字幕 日韩 欧美| 中文字幕欧美日韩高清| 狠狠做深爱婷婷综合一区| 麻豆精品在线视频| 高清国产va日韩亚洲免费午夜电影| 少妇精品网站| 国产精品高清国产三级囯产AV| 毛片免费高清免费| 欧美日韩免费在线视频| 狠狠色噜噜狠狠狠狠奇米777 | 再看日本中文字幕在线观看| 国产美女在线观看| 8090午夜无码专区| 美女一级免费毛片| 国产成人亚洲毛片| 国产嫩草在线观看| 亚洲第一区在线| 久久精品丝袜高跟鞋| 亚洲乱伦视频| 亚洲国产天堂久久综合| 亚洲国产亚综合在线区| 99精品福利视频| 2022国产91精品久久久久久| 亚洲a级在线观看| 久久动漫精品| 日本成人在线不卡视频| 亚洲一区二区三区麻豆| 国产乱人伦精品一区二区| 成人午夜网址| 伊人久久婷婷五月综合97色| 久久国产成人精品国产成人亚洲| 亚洲欧美另类日本| 中文字幕在线观| 呦女精品网站| 成人免费黄色小视频| 久久久久无码国产精品不卡| 又粗又大又爽又紧免费视频| 亚洲欧洲日产国产无码AV| 一区二区自拍| 538精品在线观看| 日本午夜三级| 99青青青精品视频在线| 她的性爱视频| 国产香蕉97碰碰视频VA碰碰看| 无码网站免费观看| 日韩欧美成人高清在线观看| 国产在线自乱拍播放| 2021亚洲精品不卡a| 久草青青在线视频| 日韩欧美综合在线制服| 亚洲毛片网站| 高清久久精品亚洲日韩Av| 亚洲美女久久| 欧美日本视频在线观看| 18禁影院亚洲专区| 久久国产精品无码hdav| 亚洲人成影院在线观看| 欧美日韩国产在线人| 在线精品视频成人网| 中国一级特黄视频| 久久精品国产一区二区小说| 高清免费毛片| 亚洲精品自产拍在线观看APP| 看av免费毛片手机播放| 在线观看免费AV网| 国产在线一区视频| 无码人中文字幕| AV老司机AV天堂| 亚洲成aⅴ人片在线影院八| 91精品在线视频观看| 中文字幕欧美日韩高清| 国产精品成人一区二区| 精品福利国产| a毛片基地免费大全| 在线色国产| 欧美一级高清视频在线播放| 天天躁狠狠躁|