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

旋轉機械運行穩定性劣化的1.5維譜特征提取方法

2016-01-28 03:32:02蔣章雷徐小力
中國機械工程 2015年23期
關鍵詞:特征提取特征故障

蔣章雷 徐小力

北京信息科技大學現代測控技術教育部重點實驗室,北京,100192

?

旋轉機械運行穩定性劣化的1.5維譜特征提取方法

蔣章雷徐小力

北京信息科技大學現代測控技術教育部重點實驗室,北京,100192

摘要:風電機組是一種典型的大型旋轉機械,其運行狀態具有變工況、非平穩特點,運行中工況和負載等非故障因素的變化會造成信號能量產生變化,因此傳統的基于能量的振動級值及功率譜方法難以有效實現運行穩定性劣化特征的提取。針對該情況,研究了高階累積量運行穩定性劣化特征提取方法,提出了基于敏感性、趨勢性、差異性、一致性判斷特征提取方法的趨勢預測適用性的方法,通過轉子實驗臺多種劣化類型下不同劣化程度狀態的實驗研究,檢驗了1.5維譜方法對狀態劣化的敏感性、趨勢性、差異性、一致性,評估了該方法作為特征提取手段的性能,解決了風電機組傳動系統運行穩定性劣化的狀態診斷、劣化趨勢預測中特征提取方法的選擇缺少理論依據的問題。

關鍵詞:傳動系統;穩定性劣化;1.5維譜;特征提取

0引言

隨著大型風電機組的建設和運行,風電機組故障數量逐年增加,由此造成的經濟損失越來越嚴重。風電機組傳動系統中的傳動軸、齒輪箱等是風電機組故障多發部件[1]。從一系列公開發表的資料和文獻來看,傳動系統的狀態診斷及趨勢預測研究日益成為國內外科學研究中一個熱點課題[2-3],而特征提取方法是研究的難點[4]。

在信號處理領域,如果信號或噪聲被假定為服從高斯分布,那么僅用基于二階累積量的功率譜和相關函數等方法便可提取特征信息,這就是功率譜方法成為最基本和最有用的信號處理工具的原因。但對于非線性系統,即使輸入是高斯信號,其輸出也是非高斯信號。由于功率譜等價[5-6]現象的存在,此時應用基于二階累積量的方法進行特征提取就會丟失信號的相位信息且無法避免背景噪聲的干擾。

高階累積量是解決非線性、非平穩、非高斯、非最小相位、非因果信號的主要手段,它具有許多優良特性。理論上高階累積量可以完全抑制高斯有色噪聲的影響。建立在二階累積量基礎上的方法抑制了信號的相位信息,不能全面描述信號特性,而建立在高階累積量基礎上的雙譜和三譜則保留了非高斯信號的幅度和相位信息。因此,利用高階累積量提取不同于二階累積量的信號特征具有重要的研究價值[7-8]。

1高階累積量的基本理論

1.1高階累積量的定義

令單個隨機變量x具有概率密度f(x),定義其第一特征函數如下:

(1)

定義單個隨機變量x的累積量生成函數(第二特征函數)如下:

Ψ(ω)=lnφ(ω)

(2)

定義單個隨機變量x的k階累積量為第二特征函數的k階導數在原點的值:

(3)

將單個隨機變量的高階累積量定義進行推廣,得到隨機向量的高階矩與高階累積量的定義如下[5]。

令x=(x1,x2,…,xk)T是一隨機向量,且ω=(ω1,ω2,…,ωk),則隨機向量的第一特征函數為φ(ω1,ω2,…,ωk)=E(ej(ω1x1+…+ωkxk)),對其求r=v1+v2+…+vk次偏導數,則有

(4)

若令ω1=…=ωk=0,則由式(4)得到隨即向量x的r階矩的定義:

mv1…vk=

(5)

參照隨機變量的第二特征函數定義得到隨機向量的第二特征函數如下:

Ψ(ω1,ω2,…,ωk)=lnφ(ω1,ω2,…,ωk)

(6)

則隨機向量的r階累積量定義為

(7)

特別地,取v1=v2=…=vk=1,得到k階矩和k階累積量,記作:mk=m1,…,1=mom(x1,…xk),ck=c1,…,1=cum(x1,…xk)。

設數列{x(n)}為零均值k階平穩隨機過程,則其k階矩定義為

mkx(τ1,τ2,…,τk-1)=

mom(x(n),x(n+τ1),…,x(n+τn-1))

(8)

其k階累積量定義為

ckx(τ1,τ2,…,τk-1)=

cum(x(n),x(n+τ1),…,x(n+τn-1))

(9)

1.2高階累積量譜的定義

高階累積量的使用可避免高斯有色觀測噪聲的影響;理論上,高斯過程大于二階的累積量譜為零,而非零累積量譜提供了一種非高斯性度量,工程上通常用高階累積量譜作為分析非高斯隨機過程的主要數學工具[5,9]。

若式(9)中高階累積量ckx(τ1,τ2,…,τk-1)是絕對可和的,即

(10)

則k階累積量譜定義為k階累積量的k-1維Fourier變換,即

Skx(ω1,ω2,…,ωk-1)=

(11)

根據式(9),三階累積量定義為

c3x(τ1,τ2)=cum(x(n),x(n+τ1),x(n+τ2))

(12)

根據式(11),三階累積量譜定義為

(13)

習慣上,高階累積量譜常簡稱作高階譜或多譜,最常用的高階譜是三階譜和四階譜,我們稱三階譜為雙譜,稱四階譜為三譜。

特別地,S2x(ω)稱為數列{x(n)}的功率譜。

21.5維譜特征提取方法

高階累積量及其高階累積量譜是多維函數,維數的增加將導致對實際信號做高階累積量估計時的計算復雜性大大提高;同時對于三維以上的數據體難以用圖形進行表示,因此必須降低高階累積量的維數。Mendel[10]提出利用高階累積量的一維片段及其一維傅里葉變換作為從高階累積量中提取有用信息的方法,即對高階累積量做切片處理,這主要包括1.5維譜以及四階累積量對角切片譜方法。以下給出 1.5維譜的定義。

考慮一個非高斯平穩隨機過程x(n),它的三階累積量為c3x(τ1,τ2)=cum(x(n),x(n+τ1),x(n+τ2)),定義三階累積量一維切片為

(14)

式(14)稱為三階累積量一維對角切片[11]。

定義三階累積量一維對角切片譜為

(15)

式(15)為三階累積量一維對角切片的一維傅里葉變換,定義為三階累積量對角切片譜,或稱為1.5維譜。1.5維譜方法僅使用了三階累積量的對角線切片,而且只進行了一維傅里葉變換,因此可以大大減小計算量。

31.5維譜特征提取方法的特征提取性能實驗研究

本文對風電機組傳動系統的多種劣化類型下不同劣化程度狀態進行特征提取研究,以驗證高階累積量特征提取方法的特征提取性能。為進行風場風電機組難以實施的各種劣化實驗,通過轉子實驗臺獲取風電機組傳動系統相關振動信息。

3.1實驗設備及實驗數據的采集

實驗數據采集現場如圖1所示,轉子轉速設定為1700r/min,以模擬額定風速下某型風電機組齒輪箱輸出軸轉速;數據采集系統選用丹麥B&K公司的PULSE系統,采集儀型號為3050-B-060,為六通道采集儀;通過加速度振動傳感器采集振動信號,傳感器型號為INV9832A,為三向加速度振動傳感器;采樣頻率為10 240Hz,每組振動數據包含2048個采樣點。

圖1 轉子實驗臺數據采集

(1)轉子正常、轉子不平衡運行狀態。轉子不平衡是風電機組中較常見的一種故障狀態[12],以下設置轉子多個不平衡運行狀態,進行模擬實驗。轉子上固定有旋轉圓盤,圓盤沿圓周均布著16個圓孔,通過在圓孔中擰入不同數量的配重螺釘實現轉子不平衡運行狀態下不同劣化程度的模擬:將兩個、四個、六個配重螺釘擰入旋轉圓盤,設定此時轉子運行狀態為輕度不平衡狀態、中度不平衡、重度不平衡。給出每種狀態下各自旋入配重螺釘的總質量與旋轉圓盤質量的百分比,以此作為不平衡劣化程度的度量標準,如表1所示。通過實驗得到轉子正常運行、不平衡劣化類型下不同劣化程度的振動數據。

表1 平衡度列表

(2)基礎松動狀態。松動故障會使得機械系統出現偏心而產生不平衡力,當不平衡力超過重力時,機器就會被不平衡力周期性地抬起,導致機械系統出現更嚴重的故障狀態[13]。松動轉子實驗臺上軸承座支撐架與底座的緊固螺釘,以減弱支撐結構彈性剛度,模擬基礎松動狀態。將軸承座支撐架的一個緊固螺釘旋松一圈代表輕度基礎松動狀態;將緊固螺釘旋松兩圈代表中度基礎松動狀態。采集得到基礎松動狀態下不同劣化程度狀態的振動數據。

(3)動靜碰摩狀態。動靜碰摩除了引起轉子、定子之間的間隙增大,軸承支承磨損之外,更重要的是會引發轉軸不穩定的反進動,導致接觸面上產生很大的法向力和摩擦力,轉子的反向渦動不收斂,反向渦動轉速迅速增大,最終引發旋轉機械的嚴重損壞[13]。在轉子實驗臺中利用一個頂桿來進行動靜碰摩實驗[14],將動靜件碰摩定性地分成4個階段:剛開始碰摩、早期尖銳觸碰、中期半尖銳碰摩和晚期平鈍型碰摩,動靜碰摩原理圖見圖2a。結合本文的研究,選擇在轉子實驗臺上固定一個可以旋入碰摩塊的支撐架,碰摩塊與轉軸之間的間隙可通過螺紋調整,如圖2b所示,其中支撐架為可以旋入碰摩塊的支撐架。設置動靜碰摩不同劣化程度的兩個狀態,將剛開始碰摩到早期尖銳型觸碰作為輕度動靜碰摩狀態,將早期尖銳型觸碰到中期半尖銳型碰摩作為中度動靜碰摩狀態,采集轉子系統振動信號。

(4)耦合故障。在轉子系統實際運行中,兩種或兩種以上形式的故障有可能同時出現,即出現耦合故障,耦合故障的存在會增加轉子系統故障診斷、趨勢預測的復雜性。設置轉子輕度不平衡與輕度動靜碰摩狀態為耦合故障輕度狀態;中度不平衡與中度動靜碰摩狀態為耦合故障中度狀態。同理,設置轉子不平衡與基礎松動的耦合故障的輕度與中度狀態、動靜碰摩與基礎松動的耦合故障的輕度與中度狀態。采集每種狀態下的振動信號。

(a)動靜碰摩原理圖

(b)動靜碰摩實驗圖圖2 動靜碰摩實驗

3.2特征提取方法的特征提取性能分析

應用1.5維譜方法對傳動系統多種劣化類型不同劣化程度數據進行特征提取處理,顯示特征提取結果,得到其劣化特征譜圖。以下通過特征譜圖分析特征提取方法的劣化特征提取性能。

3.2.1單一劣化類型的特征提取

對轉子正常運行狀態、不平衡運行狀態的振動數據進行基于1.5維譜特征提取操作,得到轉子相關運行狀態的特征譜圖,如圖3~圖6所示。

圖3 轉子正常運行狀態的時域波形圖以及特征譜圖

由圖3~圖6可看出,隨著轉子系統不平衡狀態的逐漸劣化,時域波形圖給出的趨勢性信息并不明顯,加速度幅值隨狀態劣化沒有呈現明顯規律性。對特征譜圖的具體分析如下:

圖4 轉子輕度不平衡狀態的時域波形圖以及特征譜圖

圖5 轉子中度不平衡狀態的時域波形圖以及特征譜圖

圖6 轉子重度不平衡狀態的時域波形圖以及特征譜圖

(1)轉子系統正常運行狀態下的特征譜圖中最大峰值位于224Hz處,數值為1941.5V2,系統在轉頻的高次諧波附近存在峰值,說明系統存在非線性耦合現象。系統出現輕度不平衡后,最大峰值數值急劇增大到4519V2。特征譜圖中最大峰值在系統由正常狀態劣化到輕度不平衡狀態時變化十分顯著,這表明1.5維譜特征提取方法對于系統狀態劣化具有敏感性,這有利于揭示系統早期劣化進程,適合進行運行穩定性劣化趨勢的預測研究。

(2)在轉子不平衡狀態特征譜圖的最大峰值數值明顯變化的同時,譜峰的分布形態并未出現明顯變化,表明同一種劣化類型下的劣化特征具有一致性,這對于風電機組傳動系統的狀態診斷研究提供了良好基礎。

(3)隨著轉子系統不平衡狀態的劣化,最大峰值數值在輕度不平衡狀態時急劇增大到4519V2,位于224Hz處;中度不平衡狀態時達到5070.5V2,出現在208Hz處;重度不平衡狀態時達到5504.7V2,出現在232Hz處。轉子系統不平衡狀態下不同劣化程度的最大峰值數值存在明顯的趨勢性,這給傳動系統運行穩定性狀態劣化趨勢研究提供了良好的基礎。

(4)轉子系統正常狀態下頻率成份比較單一,劣化到不平衡狀態后,譜圖中低頻部分與高頻部分的峰值數量明顯增多,出現了多譜峰現象。正常運行狀態下的特征譜圖與不平衡運行狀態下的特征譜圖具有明顯的差異性。

多譜峰現象的出現表明系統具有了更強的非高斯性和非平穩性,峰值由低頻部分向高頻部分的擴散表明不平衡旋轉引起的沖擊頻率及其各階倍頻通過非線性耦合產生了新的頻率成分。

對基礎松動劣化狀態數據、動靜碰摩劣化狀態數據進行特征提取,結果表明,系統出現輕度劣化后,特征譜圖中峰值數值變化明顯,特征提取方法對系統由正常運行狀態遷移到輕度劣化狀態具有敏感性;同種劣化類型下不同劣化程度的特征譜圖具有一致性;系統由正常運行狀態逐漸劣化的過程表現出趨勢性。同時,轉子正常運行狀態特征譜圖中譜峰分布形態與不同劣化類型下特征譜圖中譜峰分布形態具有比較明顯的差異性。

3.2.2耦合劣化類型的特征提取

利用1.5維譜特征提取方法對轉子輕度不平衡耦合輕度基礎松動狀態、轉子中度不平衡耦合中度基礎松動狀態下的數據進行特征提取操作,得到耦合故障輕度狀態、中度狀態下的特征譜圖,如圖7所示。

由轉子系統正常運行狀態下的特征譜圖及圖7可看出:

(a)輕度狀態

(b)中度狀態圖7 耦合故障輕度、中度狀態(不平衡與基礎松動)的特征譜圖

(1)隨著轉子系統由正常運行狀態劣化到耦合故障輕度狀態,特征譜圖中最大峰值數值由1941.5V2(位于224Hz處)減小到1762.1V2(最大峰值出現在208Hz處),這表明該特征提取方法對于系統耦合劣化狀態的敏感性低于單一劣化狀態。

(2)特征譜圖中譜峰的分布形態不盡相同。在轉子系統正常運行狀態下的特征譜圖中只有一個明顯波峰;出現耦合故障中度劣化后,特征譜圖中出現了明顯的多個波峰:784Hz處的數值為979.3V2的波峰、264Hz處的數值為849V2的波峰,以及656Hz處的數值為511.9V2的波峰。這表明對于耦合故障類型,該特征提取方法提取的劣化特征的一致性低于單一劣化類型。

(3)特征譜圖中最大峰值數值隨著劣化程度的增加呈現不斷減小的趨勢,這表明,耦合故障狀態下不同劣化程度的最大峰值數值存在一定的趨勢性,但這種趨勢性與單一劣化類型下所表現出的趨勢性相反。

分析不平衡耦合動靜碰摩狀態、動靜碰摩耦合基礎松動狀態的特征譜圖,結果表明:當耦合故障狀態中某個單一劣化類型占主要成分時,該特征提取方法對于狀態劣化的敏感性、一致性有所提高。

3.3特征提取方法的趨勢預測適用性討論

針對多種劣化類型下不同劣化程度狀態數據,基于1.5維譜方法進行特征提取分析,得到其對狀態劣化的差異性、一致性、敏感性、趨勢性信息,如表2所示。

表2 1.5維譜的劣化特征提取性能分析

以下從特征提取方法對于狀態劣化的敏感性、趨勢性、差異性、一致性角度討論其作為穩定性劣化趨勢預測的適用性。

(1)對于單一劣化類型,1.5維譜特征提取方法具有良好的差異性、一致性,因此適合進行單一劣化類型的診斷研究;1.5維譜特征提取方法也具有良好的敏感性、趨勢性,因此也適合進行穩定性劣化趨勢預測研究。

(2)對于耦合劣化類型,1.5維譜特征提取方法具有良好的差異性,但不具有一致性,不能用于耦合劣化類型的診斷研究;其具有良好的敏感性,但不具有趨勢性,不能用于穩定性劣化趨勢預測研究。

4結束語

針對風電機組變工況、非平穩運行狀態,研究利用1.5維譜特征提取方法分離劣化特征與變工況等非劣化特征,提出了基于敏感性、趨勢性、差異性、一致性判斷特征提取方法的趨勢預測適用性的方法,并通過轉子試驗臺驗證了該方法的有效性,為風電機組傳動系統運行穩定性劣化的狀態診斷、劣化趨勢預測中特征提取方法的選擇提供理論依據。

參考文獻:

[1]JohanR,MargaretaBL.SurveyofFailuresinWindPowerSystemswithFocusonSwedishWindPowerPlantsduring1997-2005[J].IEEETransactionsonEnergyConversion,2007,22(1):167-173.

[2]于傳洋,徐寶杰,吳國新,等.一種風電機組的故障預測算法[J].北京信息科技大學學報,2013,28(4):35-38.

YuChuanyang,XuBaojie,WuGuoxin,etal.AFaultPredictionAlgorithmofWindTurbineGenerator[J].JournalofBeijingInformationScienceandTechnologyUniversity, 2013,28(4): 35-38.

[3]郭學鵬,王吉芳.無線傳感網絡在風機齒輪箱狀態監測中的應用[J].北京信息科技大學學報,2013,28(3): 19-22.

GuoXuepeng,WangJifang.ApplicationResearchofWirelessSensorNetworksintheGearboxofWindTurbineGeneratorMonitoring[J].JournalofBeijingInformationScienceandTechnologyUniversity, 2013,28(3): 19-22.

[4]王國彪,何正嘉,陳雪峰,等.機械故障診斷基礎研究“何去何從”[J].機械工程學報,2013,49(1): 63-72.

WangGuobiao,HeZhengjia,ChenXuefeng.etal.BasicResearchonMachineryFaultDiagnosis-whatisthePrescription[J].JournalofMechanicalEngineering, 2013,49(1): 63-72.

[5]張賢達.時間序列分析——高階統計量方法[M].北京:清華大學出版社,1996.

[6]BoxGE,JenkinsGM.TimeSeriesAnalysis:ForecastingandControl[M].SanFrancisko:Holden-Day, 1970.

[7]NikiasL,RaghuveerMR.BispectrumEstimation:aDigitalSignalProeessingFramework[J].Proc.ofIEEE,1987,75(7):869-891.

[8]黃方林,顧松年.結構故障高階譜診斷的一種方法[J].振動工程學報,1991(1): 33-40.

HuangFanglin,GuSongnian.AMethodofStructuralFaultDiagnosisbyHigherOrderSpectrum[J].JournalofVibrationEngineering, 1991(1): 33-40.

[9]邱天爽,張旭秀,李小兵,等.統計信號處理—非高斯信號處理及其應用[M].北京:電子工業出版社,2004.

[10]MendelJM.TutorialonHigher-orderStatistics(Spectra)inSignalProcessingandSystemTheoryTheoreticalResultsandSomeApplication[J].ProceedingsoftheIEEE,1991,79(3):278-305.

[11]NikiasCL,MendelJM.SignalProcessingwithHigher-orderSpectra[J].IEEESignalProcessingMagazine,1993,10(3):10-37.

[12]VillaLF,RenonesA,PeranJR,etal.StatisticalFaultDiagnosisBasedonVibrationAnalysisforGearTest-benchunderNon-stationaryConditionsofSpeedandLoad[J].MechanicalSystemsAndSignalProcessing,2012,29:436-446.

[13]聞邦椿,武新華,丁千,等.故障旋轉機械非線性動力學的理論與試驗[M].北京:科學出版社,2004.

[14]胡蔦慶,張雨,劉耀宗,等. 轉子系統動靜件間尖銳碰摩時的振動特征試驗研究[J].中國機械工程,2002,13(9): 777-780.

HuYuanqing,ZhangYu,LiuYaozong,etal.ExperimentResearchonVibrationCharacteristicsofSharpRup-impactbetweenRotorandStatorofaRotorSystem[J].ChinaMechanicalNgineering,2002,13(9): 777-780.

(編輯袁興玲)

Feature Extraction Method of 1.5-dimensional Faced to

Running Stability Deterioration of Rotating Machinery

Jiang ZhangleiXu Xiaoli

Key Laboratory of Modern Measurement & Control Technology,Ministry of Education,

Beijing Information Science and Technology University,Beijing,100192

Abstract:Wind turbines were typical large rotating machinery with variable conditions and non-stationary states, and vibration signal energy was changed because of non-fault factor such as variable conditions and loads. Fault diagnosis method based on vibration signal energy such as vibration level and power spectrum method could not extract running stability deterioration features from non-stationary signals under the influences of variable conditions and loads. 1.5-dimensional spectrum feature extraction method was studied based on higher-order cumulant method,and an evaluation method was proposed by using sensitivity property, trend property, difference property, and consistency property. Experimental data of varying degrees of deterioration under various types of deterioration were carried out to validate the proposed method, and the results show that feature extracting method can separate the deterioration characteristics from non-degradation characteristics, and the evaluation method can provide theoretical basis for selecting feature extraction method for fault diagnosis and trend prediction of running stability deterioration.

Key words:transmission system; stability deterioration; 1.5-dimensional spectrum; feature extraction

基金項目:國家自然科學基金資助項目(51275052);北京市自然科學基金資助重點項目(3131002);北京市教委科研計劃資助重點項目(KZ201311232036);現代測控技術教育部重點實驗室開放課題(KF20141123202)

收稿日期:2014-10-20修回日期:2015-09-30

中圖分類號:TH17DOI:10.3969/j.issn.1004-132X.2015.23.015

作者簡介:蔣章雷,男,1983年生。北京信息科技大學現代測控技術教育部重點實驗室助理研究員、博士。主要研究方向為機電裝備運行狀態監控檢測技術。獲北京市科學技術三等獎1項。發表論文10余篇。徐小力,男,1951年生。北京信息科技大學現代測控技術教育部重點實驗室教授、博士研究生導師。

猜你喜歡
特征提取特征故障
故障一點通
如何表達“特征”
基于Gazebo仿真環境的ORB特征提取與比對的研究
電子制作(2019年15期)2019-08-27 01:12:00
不忠誠的四個特征
當代陜西(2019年10期)2019-06-03 10:12:04
抓住特征巧觀察
一種基于LBP 特征提取和稀疏表示的肝病識別算法
奔馳R320車ABS、ESP故障燈異常點亮
故障一點通
江淮車故障3例
基于MED和循環域解調的多故障特征提取
主站蜘蛛池模板: 全午夜免费一级毛片| 伊人精品成人久久综合| 在线免费a视频| 99热国产这里只有精品无卡顿" | 激情无码视频在线看| 国产一在线观看| 91探花在线观看国产最新| 天天做天天爱夜夜爽毛片毛片| 亚洲伊人久久精品影院| 中文字幕调教一区二区视频| 亚洲男人天堂2020| 尤物亚洲最大AV无码网站| 亚洲天堂精品在线观看| 亚洲伊人久久精品影院| 极品尤物av美乳在线观看| 国产国语一级毛片| 91精品国产91久无码网站| 国产视频入口| 国产精品免费福利久久播放| 国产精品久久久久久久伊一| 99精品这里只有精品高清视频| 国产亚洲日韩av在线| 色婷婷成人| 色婷婷国产精品视频| 久久99热这里只有精品免费看| 亚洲国产中文欧美在线人成大黄瓜 | 女人毛片a级大学毛片免费| 久久婷婷国产综合尤物精品| 99视频全部免费| 91www在线观看| 精品国产中文一级毛片在线看| 久久99蜜桃精品久久久久小说| 亚洲黄色网站视频| 国产在线观看一区精品| 青青草国产在线视频| 中国国产高清免费AV片| 无码福利视频| 伊人丁香五月天久久综合| 亚洲国产日韩欧美在线| 69视频国产| 成人免费一区二区三区| 熟妇无码人妻| 亚洲v日韩v欧美在线观看| 在线观看网站国产| 久久国产成人精品国产成人亚洲 | 日本午夜三级| 国内精品一区二区在线观看| 精品无码一区二区在线观看| 欧美人人干| 99青青青精品视频在线| 亚洲国产精品VA在线看黑人| 成年网址网站在线观看| 国产精品一老牛影视频| 国产中文在线亚洲精品官网| 成人午夜免费观看| 中文字幕自拍偷拍| 九九热这里只有国产精品| 国产一区二区三区免费观看 | 97青草最新免费精品视频| 在线观看国产网址你懂的| 最新日本中文字幕| 国产主播喷水| 国产高清在线观看| 91美女在线| 国产人妖视频一区在线观看| 日韩天堂网| 国产日韩欧美中文| 热热久久狠狠偷偷色男同 | 少妇精品久久久一区二区三区| 精品国产中文一级毛片在线看 | 国产成人凹凸视频在线| 亚洲区第一页| 亚洲精品免费网站| 国产一级毛片yw| 播五月综合| 九色国产在线| 久久综合色天堂av| 乱人伦99久久| 中文字幕久久波多野结衣 | 91啦中文字幕| 狠狠做深爱婷婷综合一区| 国产在线精品网址你懂的|