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

基于DSP-LightGBM算法的風電機組齒輪箱高速軸齒輪健康狀態監測研究

2024-12-31 00:00:00吳智泉于海瀛王松吳文韜左康會
太陽能 2024年11期

摘 要:齒輪箱高速軸齒輪作為齒輪箱最重要的部件之一,其健康狀態不僅對風電機組的安全運行與維護影響重大,還對風電機組的發電量有直接影響。基于甘肅省某風電場中風電機組的狀態監測系統(CMS)采集的大量風電機組齒輪箱高頻振動信號,利用數字信號處理(DSP)方法與LightGBM算法構建了可監測風電機組齒輪箱高速軸齒輪健康狀態的模型,利用DSP方法中的快速傅里葉變換(FFT)方法對信號進行樣本擴增;同時從時域與頻域的角度進行特征構造,篩選并剔除冗余特征,進而構建了基于Dsp-LightGBM算法的風電機組齒輪箱高速軸齒輪健康狀態監測分類模型,并利用CMS采集的風電機組實際的齒輪箱高速軸齒輪高頻振動信號對所構建模型進行效果驗證。驗證結果顯示:所構建模型預測效果良好,能夠高效檢測出風電機組齒輪箱高速軸齒輪的狀態。

關鍵詞:風電機組;齒輪箱;數字信號處理;LightGBM;快速傅里葉變換;狀態監測

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

收稿日期:2023-12-08

通信作者:吳文韜(1994—),男,碩士,主要從事風電領域數智化系統方面的研究。512283725@qq.com

0" 引言

風能作為一種環保、綠色、可再生的清潔能源,在全球節能減排進程中起到了越來越重要的作用。風電的發展不僅為中國解決了眾多能源問題,也成為中國清潔能源發展的一個重要應用方向。齒輪傳動系統是風電機組的核心組成部分,而齒輪箱高速軸齒輪作為齒輪傳動系統中的重要部件之一,當其處于異常狀態時,極易被損壞。利用振動傳感器對齒輪箱關鍵部件(齒輪、軸承等)進行振動信號采集、分析和監測的風電機組狀態監測系統(condition monitor system,CMS)已逐步成為風電場的標配,但其狀態監測效果仍有明顯的提升空間。隨著中國風電裝機容量的不斷增長,開發高效、準確的齒輪箱高速軸齒輪健康狀態監測技術[1-5],具有重要意義。

目前,基于數字信號處理(DSP)方法的健康狀態監測技術已經成為國內外公認的主流技術,其中,快速傅里葉變換(FFT)是一個非常重要的工具,可實現頻域分析和信號統計特征提取等[6];而近年來基于離散小波變換(DWT)、經驗模態分解(EMD)等的非穩態信號處理方法也獲得了廣泛關注[7]。但是這些方法對于狀態監測及故障診斷[8],仍無法做到有效識別;且DSP方法在特征選擇和狀態量化評估方面仍需改進[9]。基于此,本文以甘肅省某風電場中的風電機組為例,選取CMS采集的大量齒輪箱高速軸齒輪高頻振動信號,通過DSP方法進行準確的齒輪振動信號特征提取,再利用LightGBM算法實現其直觀狀態評估,從而構建可進行風電機組齒輪箱高速軸齒輪健康狀態監測的模型雛形[6](下文簡稱為“狀態監測模型”),并通過模型訓練,最終確定模型形式,最后進行效果驗證。

1" FFT介紹

DSP方法主要是利用計算機或專用處理設備,以數值計算的方法對信號進行采集、變換、綜合、估值與識別等加工處理,以達到提取信號統計特征和便于實際應用的目的。在DSP方法發展過程中,作為其技術分支的FFT的出現使信號分析的運算速度提高了幾百倍,將信號的處理技術推向了一個新階段,隨之其也形成了一套完整的理論和實現方法。因此,本文采用FFT對CMS采集的齒輪箱高速軸齒輪高頻振動信號進行分析處理。

傅里葉變換是將信號從時域變換到頻域的一種變換形式,是信號處理領域中一種重要的分析工具;而FFT是實現離散傅里葉變換(DFT)的一種高效方式。

對于有限長數字信號x(n),其經過DFT后的頻域信號X(k)可表示為:

(1)

式中:WN為旋轉因子,具有對稱性和周期性;N為x(n)的長度;n為x(n)的離散時間點索引,表示時域信號的采樣點序號,用于確定時域信號在時間上的離散位置;k為X(k)的頻率索引,表示頻域中的不同頻率成分。

旋轉因子可表示為:

WN=e" " " " " " " " " " " " " " " " " " " " " " " " " " " " " (2)

式中:j為虛數單位,滿足j2=-1;e為自然常數。

FFT利用旋轉因子的對稱性和周期性特性,通過分解DFT矩陣,將長序列DFT的計算分解為多個短序列DFT的計算,簡化了計算過程,降低了運算時間,顯著提高了頻譜計算的效率,使頻域分析成為可能;且FFT計算結果與直接DFT計算結果精度相同。

2" LightGBM算法原理

XGBoost和pGBRT等算法可能并不適用于處理大量數據和眾多特征[10],部分原因是因為這些算法需要遍歷每個特征的所有數據樣本,以估算所有可能出現的分割點的信息增益,這可能會導致計算效率低下。因此,LightGBM算法[11]對XGBoost算法進行了優化,其是一種基于梯度提升決策樹(GBDT)的機器學習算法,囊括了基于梯度的單側采樣(GOSS)算法、互斥特征捆綁(EFB)算法及直方圖(Histogram)算法。LightGBM算法中各算法的特征及運用方式如表1所示。

GOSS算法的計算過程為:

1)輸入訓練數據、迭代次數、大梯度數據的采樣率a、小梯度數據的采樣率b,選擇弱學習器的類型。

2) GOSS算法先將所有樣本的梯度數據絕對值進行從大到小的排序,然后選取具有最大梯度的前a個樣本作為集合A;剩下的(1–a)個樣本(稱為集合Ac),從中隨機抽取b個樣本(稱為集合B),則其中隨機抽取的小梯度樣本的占比為b·(1–a);通過這種方式,在保留部分重要樣本的同時,也引入了一定比例的相對不那么重要但具有隨機性的樣本,且在利用小梯度樣本計算信息增益時,需乘以1個常數系數(1–a)/b,這樣就只需要更多地關注訓練不足的樣本,而無需較多改變原始數據分布。

Histogram算法的計算步驟為:

1)為每個特征創建1個直方圖,每個直方圖需體現兩類統計值,分別為每個bin中樣本的梯度之和、每個bin中的樣本數量。

2)對于每個特征,遍歷所有的樣本,將上述兩類統計值累積到樣本所屬的bin中。

3)對于某個葉子節點,遍歷所有的bin,分別以當前bin作為分割點,累加其左邊的bin至當前bin中樣本的梯度和及樣本數量中,并將父節點上的總梯度和與總樣本數量相減,得到右邊所有bin中樣本的梯度和及樣本數量,并以此計算該節點的信息增益。在遍歷過程中,取最大的增益值,以最大信息增益值對應的那個時刻的特征和bin的特征值分別作為分裂節點的特征和分裂特征值取值。

4)針對所有的葉子節點,重復上述操作,遍歷所有的特征,找到信息增益值最大時所對應的特征及當前葉子節點的劃分值,以此來分裂該葉子節點。

3" 狀態監測模型構建時的特征分析

3.1" 時域特征

對CMS采集的設備高頻振動信號進行分析、統計,求出相應的時域特征值。不同時域特征的計算方法如表2所示。

3.2" 頻域特征

頻域特征的提取步驟為:

1) 倍頻對齊。提取頻率特性時,必須對基于FFT的幅值頻譜圖的橫軸進行均勻的倍頻變換,即需進行倍頻對齊。因為經過數據分析發現,故障頻率所具有的特殊性一般與通過CMS收集的數據中的轉速數據有關。倍頻對齊的具體做法是:①將轉速值除以60,從而實現分鐘到秒的變換,得到以赫茲為單位的轉頻;②將幅值頻譜圖中的頻率橫坐標除以轉頻,得到倍頻;③針對額定轉速為1200 rpm (對應轉頻為20 Hz)、采樣頻率為25600 Hz(依據Nyquist采樣定理,對應的最大分析頻率為12800 Hz)的工況,估算出合適的最大倍頻和倍頻步長,并為所有樣本設定統一的以倍頻為橫坐標的取值列表;④將原幅值頻譜圖橫軸的頻率轉換為倍頻。

2) 離散化(即“區間劃分”)。幅值頻譜圖橫坐標轉換為倍頻之后,可進行倍頻區間劃分,設定相鄰區間之間有一定的取值覆蓋,從而保證特征頻率可被充分檢測。倍頻區間的劃分方式如表3所示。

3)頻域特征提取。對每個倍頻區間進行特征提取,提取方法有兩種,分別為頻域描述方法和頻域峰值方法。頻域描述方法是依次從每個倍頻區間提取表2中的時域特征。頻域峰值方法是指提取倍頻區間內相對于鄰近值的最高點,可以提取多個較高的峰值點及其橫坐標值。

以某時段數據的60~70 Hz倍頻區間為例,依據頻域峰值方法提取了30個特征峰值,如圖1所示。需要說明的是,圖中深紅色為提取出的特征峰值。

由圖1可知:采用頻域峰值方法得到的特征峰值均為此倍頻區間內的較高值。

3.3" 特征選擇

由于在構建狀態監測模型時無法預知哪些特征對于監測結果的影響程度較大,只有選擇全部的特征來訓練模型,利用模型訓練完1輪后,再剔除無用或冗余的特征。這不僅是為了實現對數據的降維和優化[6],也是為了避免模型過擬合,提高模型的泛化能力,因此,需要先對所有的特征進行重要性排名。

如果信號的特征數量過多,則可以使用決策樹算法,基于特征分裂過程進行特征選擇。本文基于基尼系數[11],根據打亂特征值,利用得到的測試數據集進行測試,通過測試值和真實目標值計算損失函數,并評估損失函數的值因特征值的隨機排序而升高了多少,只留下能對訓練數據起到分類或回歸作用的特征,剔除冗余或不重要的特征,以提高模型的學習效率和準確度。

對基于頻域描述方法提取的所有特征進行統計,并列舉重要性排名前39位的特征,如表4所示;其重要性如圖2所示。

在模型訓練時,將重要性低的頻域特征刪除,一方面可提高運算效率,另一方面可使模型更快捷地捕捉對應特征。

從重要性排名前39位的特征中選取重要性排名前10位的特征,其重要性如表5所示。

4" 實證分析

4.1" 信號獲取

本文以甘肅省某風電場的50臺風電機組為例,選取CMS采集的這些風電機組2022年5—11月的齒輪箱信號,這些信號包含高頻振動信號和溫度信息。在50臺風電機組中,共有7臺風電機組在此期間發生過齒輪箱高速軸齒輪斷齒等故障。CMS采集齒輪箱高頻振動信號時的采樣頻率為25600 Hz,風力發電機額定轉速(等同于齒輪箱高速軸轉速)為1200 rpm,以此信息作為本文研究的風電機組的機理信息。為了兼顧樣本數量和樣本特征的一致性,對樣本進行工況篩選,選取額定轉速為1200±50 rpm的信號;同時,為了兼顧樣本的泛化能力,按比例選取了轉速在800~1050 rpm之間的信號,保證了樣本信號處于較為統一的轉速水平。

4.2" 樣本擴增

總體來說,CSM實際采集到的風電機組信號中,大部件的異常信號占比較小[9]。雖然本文使用的信號中有7臺風電機組在所選取的時間段發生過齒輪箱高速軸齒輪斷齒等故障,但在此期間均僅發生過1次,因此故障樣本信號占比較小。雖然在特征工程階段會對信號進行大量的變換和統計等特征構造、特征選擇工作,以方便通過機器學習算法更加快捷地捕捉到問題本質,但機器學習算法對信號樣本數量有較高的要求,因此,本文對樣本數量進行擴增。

樣本擴增是指通過合理(符合風電機組的機械工程學原理)的方式對原始信號進行分割、變換、組合、加噪等處理,構造出更多的樣本供機器學習算法學習,以提高模型預測的準確性和泛化能力。根據本文使用的CMS采集的信號特點,采用時長分割和組合均值兩種方式進行樣本擴增。

4.2.1" 時長分割

時長分割是指將原始文件按照時長平均分割成更小的文件,以達到擴增樣本數量的目的。例如:本文使用的CMS采集的信號原始文件的時長是4 s,可以平均分割為4個1 s的文件,也可以分割為兩個2 s的文件。分割成的小文件與原始文件之間會存在一定的差異性,在一定程度上提高了樣本庫的豐富性、多樣性;但隨著分割數量的增加,分割成的小文件包含的信息量會隨之減小,這將在一定程度上增加學習器的學習難度。因此,采用時長分割方式時,分割的樣本數量需折中考慮。

4.2.2" 組合均值

組合均值是指將兩個或兩個以上的原始信號文件,以均值的方式組合,從而形成新的樣本。由于樣本的時域、頻域信號形式之間存在差異,因此時域組合與頻域組合方法需分別討論。

1) 時域組合。由于原始文件內容為時域波形信號,不同信號文件的波峰、波谷、相位等的變化不一致,因此不能簡單地將兩個原始文件的時域波形信號求平均值,可將時域波形信號依據相應的統計指標計算出來之后,再將這些值求平均值。此處的統計方法主要是指各類時域特征,這意味著在時域組合之前,需要完成時域特征的提取工作。樣本的兩個原始文件通過時域組合方式得到的新樣本的特征值如表6所示。

2) 頻域組合。頻域組合可基于FFT的幅值頻譜信號進行,具體步驟為:①倍頻對齊。倍頻對齊的目的,一方面是為了特征構造,在一定程度上緩解轉速差異造成的特征頻率差異;另一方面是為了樣本擴增,為有一定轉速差異的樣本組合和樣本擴增打下基礎。②幅值求均值。通過頻域組合方式得到的新樣本的特征值如表7所示。由于以倍頻(fft_x)為橫坐標的樣本信號已經對齊,只需對以幅值(fft_y)為縱坐標的樣本信號進行均值處理即可;頻域組合求均值后,可按照正常的特征工程等步驟進行后續處理。

4.3" 模型訓練

下文基于狀態監測模型,進行模型訓練。

4.3.1" 數據劃分

數據劃分是指將數據集劃分為訓練集、驗證集和測試集的過程,其結果會直接影響模型的泛化能力與運算速度等性能。因此,為了避免數據劃分不均衡對模型預測結果產生的不利影響,本文數據劃分以風電機組機位號、故障時間等變量為維度,對數據進行分層抽樣。

同時,為了進一步提升狀態監測模型的泛化能力,本文將進行10折交叉驗證,步驟為:先將數據集隨機劃分為10份大小相同的互斥子集,輪流將1份作為測試集,剩下的9份作為訓練集,進行試驗。

4.3.2" 評價指標的計算

分類問題的本質是將預測目標分為不同的類別或標簽。設定以故障概率0.5作為預測結果的截斷點,將故障概率大于等于0.5的樣本標記為“預測為1”;將故障概率小于0.5的樣本標記為“預測為0”。在分類問題中,通常會出現4種不同的預測結果,分別為:真陽性(true positive,TP),即真實為1,預測為1;假陰性(1 negative,FN),即真實為1,預測為0;假陽性(1 positive,FP),即真實為0,預測為1;真陰性(true negative,TN),即真實為0,預測為0。對落入這些結果的樣本數量進行統計,并用于下述評價指標的計算,對分類模型的準確性進行多維度評價。

1) ROC_AUC。ROC_AUC是指受試者工作特征(receiver operating characteristic,ROC)曲線下方包圍的面積。ROC曲線是通過一系列不同的截斷點得到的真正率TPR和假正率FPR形成的坐標點連接繪制而成。真正率和假正率的計算式分別為:

(16)

式中:BTP為真陽性樣本數量;BFN為假陰性樣本數量。

(17)

式中:BFP為假陽性樣本數量;BTN為真陰性樣本數量。

2) 正確率的計算方法為:

(18)

3) 精確率的計算方法為:

(19)

4) 召回率的計算方法為:

(20)

5) F1分數的計算方法為:

(21)

4.3.3" 模型訓練結果

狀態監測模型訓練得到的評價指標對應的值如表8所示。

由表8可知:ROC_AUC值基本達到了0.95及以上,說明本文模型有較好的預測能力。通過觀察以0.5為故障概率截斷點時對應的各指標取值,可以預估在調整故障概率截斷點之后某個評價指標可達到的最佳狀態。表8也給出了第n(1≤n≤10)折交叉驗證時,正確率達到最佳值時對應的故障概率截斷點取值。

4.4" 模型的超參數設置

經過10折交叉驗證后,綜合狀態監測模型的ROC_AUC、正確率、召回率等多項評價指標,根據多項評價指標均較好的結果,確定了狀態監測模型的超參數,主要超參數設定如表9所示。

由表9可知:狀態監測模型經過訓練后,最終形成的模型為基于DSP-LightGBM算法的風電機組齒輪箱高速軸齒輪的狀態監測分類模型(以下簡稱為“本文模型”)。

5" 模型驗證

從該風電場的50臺風電機組中選出1臺在2018年6—11月發生過齒輪箱高速軸齒輪斷齒故障的風電機組,利用當時CMS采集的其齒輪箱高頻振動信號進行本文模型預測結果驗證,得到的故障概率預測結果如圖3所示。其中:2018年9月1日的風電場定檢報告上顯示齒輪箱發生了高速軸齒輪斷齒故障,并在2018年9月23日完成修復,該風電機組恢復正常運行;由于9月2—22日風電機組停機,無有效監測數據。

由圖3可知:對比齒輪箱高速軸齒輪斷齒故障修復前后的數據,本文模型的故障概率有明顯下降,說明本文模型的預測結果與風電機組實際情況相吻合,模型預測效果良好。預測的故障概率值可以作為齒輪箱高速軸齒輪斷齒故障的評價指標,當該指標較高時,代表故障概率高、齒輪箱健康狀況不佳;當該指標偏低時,代表故障概率低、齒輪箱健康狀況良好。

后期經過大量數據驗證發現:在實際應用中,可將本文模型的故障概率閾值設定為0.8,若本文模型的故障概率大于閾值,則模型告警,提示齒輪箱高速軸齒輪存在異常;否則,模型顯示正常。

6" 結論

本文基于甘肅省某風電場中風電機組的CMS采集的大量風電機組齒輪箱高頻振動信號,利用DSP方法與LightGBM算法構建了可監測風電機組齒輪箱高速軸齒輪健康狀態的模型,利用DSP方法中的FFT方法,對信號進行了樣本擴增;同時,從時域與頻域的角度進行特征構造,篩選并剔除冗余特征,進而構建了基于LightGBM算法的風電機組齒輪箱高速軸齒輪健康狀態監測分類模型,并利用CMS采集的風電機組實際的齒輪箱高速軸齒輪高頻振動信號對所構建模型的預測效果進行了驗證。驗證結果顯示:所構建模型預測效果良好,能夠高效檢測出風電機組齒輪箱高速軸齒輪的狀態。

目前該模型還存在一定的不足和局限性,主要原因在于故障案例、故障種類、樣本數量都相對偏少,后續可通過搜集更多的故障案例、故障種類、樣本數量來對模型進行調整和優化。

[參考文獻]

[1] 尹詩,侯國蓮,于曉東,等. 基于SCADA數據的風電機組齒輪箱狀態監測方法[J]. 太陽能學報,2021,42(1):324-332.

[2] REN J,YU Z P,GAO G L,et al. A CNN-LSTM-LightGBM based short-term wind power prediction method based on attention mechanism[J]. Energy reports,2022,8:437-443.

[3] 湯婷. 基于信息融合的風電齒輪箱軸承故障檢測方法研究[D]. 蘭州:蘭州理工大學,2021.

[4] 陳庭記,桂帆,郭政,等. 基于不同運行狀態的風電機組齒輪箱故障率預測模型[J]. 太陽能學報,2023,44(4):45-51.

[5] 黃宏臣,郭四洲,王子彥,等. 加速度包絡解調方法在風力發電機滾動軸承早期故障診斷應用研究[J]. 機械設計與制造,2022(3):251-253,257.

[6] 周慶梅,劉冰冰,胡號朋,等. 基于重構特征和寬度學習的海上風電機組齒輪箱高速軸故障預警[J]. 船舶工程,2022,44(增刊2):32-36.

[7] 彭媛寧. 基于音頻特征的風電葉片損傷監測技術研究[D]. 南京:南京航空航天大學,2018.

[8] 張鈺,陳珺,王曉峰,等. 隨機森林在滾動軸承故障診斷中的應用[J]. 計算機工程與應用,2018,54(6):100-104,114.

[9] 陳亞楠,胡凱凱,陳剛,等. 基于機器學習的風電機組齒輪箱故障預警[J]. 控制與信息技術,2021(5):108-112.

[10] 呂明珠. 風電軸承狀態監測與智能維護策略研究[J]. 電氣開關,2021,59(3):35-40.

[11] 陳維剛,張會林. 基于RF-LightGBM算法在風機葉片開裂故障預測中的應用[J]. 電子測量技術,2020,43(1):162-168.

[12] 蘇元浩,孟良,許同樂,等. 不平衡數據集下優化WGAN的風電機組齒輪箱故障診斷方法[J]. 太陽能學報,2022,43(11):148-155.

Research on Health Status Monitoring of High-Speed Shaft Gears in Wind Turbine Gearboxes Based on

DSP-LightGBM Algorithm

Wu Zhiquan1,Yu Haiying2,Wang Song1,Wu Wentao3,Zuo Kanghui4

(1. SPIC Yunnan International Power Investment Co.,Ltd,Kunming" 650228,China;

2. Harbin Energy Innovation Digital Technology Co.,Ltd,Harbin 150000,China;

3. Yunnan Power Investment Green Energy Technology Co.,Ltd,Kunming 650228,China;

4. New Energy Development Branch,SPIC Yunnan International Power Investment Co.,Ltd,Kunming 650228,China)

Abstract:As one of the most important components of the gearbox,health status of the high-speed shaft gear not only has a significant impact on the safe operation and maintenance of the wind turbine,but also has a direct impact on the power generation capacity of the wind turbine. This paper is based on a large number of high-frequency vibration signals of wind turbine gearboxes collected by the state monitoring system (CMS) of wind turbines in a wind farm in Gansu Province. A model that can monitor the health status of high-speed shaft gears in wind turbine gearboxes is constructed using digital signal processing (DSP) methods and LightGBM algorithm. The signal is simple amplified using the fast Fourier transform (FFT) method in DSP methods. Simultaneously constructing features from both time and frequency domains,filtering and removing redundant features,a classification model for the health status monitoring of high-speed shaft gears in wind turbine gearboxes based on DSP-LightGBM algorithm is constructed. The actual high-frequency vibration signals of high-speed shaft gears in wind turbine gearboxes collected by CMS are used to verify the effectiveness of the model. The verification results show that the constructed model has good prediction performance and can efficiently detect the high-speed shaft gear state of the wind turbine gearbox.

Keywords:wind turbines;gearbox;DSP;LightGBM;FFT;status monitoring

主站蜘蛛池模板: 亚洲日韩AV无码一区二区三区人| 国产欧美性爱网| 国产一区二区精品高清在线观看| 91精品啪在线观看国产91九色| 国产在线视频自拍| 久久9966精品国产免费| 欧美三级视频网站| 婷婷丁香在线观看| 日本道综合一本久久久88| 亚洲三级电影在线播放| 国产欧美视频在线| 欧美成人午夜在线全部免费| 国产精品手机视频| 亚洲AV无码乱码在线观看代蜜桃 | 无码人妻热线精品视频| 71pao成人国产永久免费视频 | 香蕉伊思人视频| 亚洲成人黄色网址| 一级在线毛片| 国产精品欧美在线观看| 欧美国产菊爆免费观看 | 中文无码精品A∨在线观看不卡 | 色网在线视频| 亚洲中久无码永久在线观看软件| 美女国内精品自产拍在线播放| 亚洲无码日韩一区| 亚洲无码高清视频在线观看| 一级毛片免费不卡在线| 一区二区偷拍美女撒尿视频| 国产在线精品99一区不卡| 成人在线观看不卡| 一本一道波多野结衣av黑人在线| 午夜电影在线观看国产1区| 亚洲欧美国产五月天综合| 四虎在线高清无码| hezyo加勒比一区二区三区| 国产成人亚洲综合A∨在线播放| 国产后式a一视频| 亚洲一级毛片| 无码日韩精品91超碰| 欧美日韩亚洲综合在线观看| 亚洲一区二区日韩欧美gif| 国产成人精品优优av| 国产麻豆福利av在线播放| 亚洲天堂精品视频| 国产亚洲欧美另类一区二区| 免费又黄又爽又猛大片午夜| 国产va在线观看免费| 亚洲无码91视频| 国产精品视频第一专区| 中文字幕亚洲综久久2021| 欧美伦理一区| 欧美另类第一页| 国产成人在线无码免费视频| a毛片免费在线观看| 无码一区二区波多野结衣播放搜索 | 国产激爽大片高清在线观看| 一区二区影院| 成人午夜亚洲影视在线观看| 国产日韩欧美精品区性色| 国产在线精品香蕉麻豆| 日本免费一区视频| 欧美日韩91| 青草视频久久| 亚洲无码高清一区二区| 全免费a级毛片免费看不卡| 国产人人射| 成人精品区| 亚洲天堂区| 国模极品一区二区三区| 亚洲二区视频| 国产国模一区二区三区四区| 91精品福利自产拍在线观看| 美女被操91视频| 国产裸舞福利在线视频合集| 亚欧乱色视频网站大全| 亚洲无线一二三四区男男| 一级黄色欧美| 国产免费久久精品99re不卡| 好紧太爽了视频免费无码| 国产福利不卡视频| 中文毛片无遮挡播放免费|