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

時頻數據驅動的典型復雜供輸機構健康狀態軟測量方法

2022-05-13 05:25:20張鋼梁偉閣佘博田福慶
兵工學報 2022年4期
關鍵詞:振動特征測量

張鋼, 梁偉閣, 佘博, 田福慶

(1.海軍工程大學 兵器工程學院, 湖北 武漢 430033; 2.海軍大連艦艇學院 導彈與艦炮系, 遼寧 大連 116018)

0 引言

復雜供輸機構是特種裝備的重要組成部分,其可靠性和穩定性直接影響特種裝備功能的發揮,如火炮供輸機構,直接影響彈藥的運動狀態[1],且其運動過程伴隨有猛烈的沖擊、振動等,力學環境十分復雜,對此復雜供輸機構進行狀態監測已經成為裝備健康管理領域關注的焦點問題[2]。

目前,復雜機械設備的健康狀態評估方法主要分為機理模型和數據驅動兩類。基于機理模型的健康狀態評估方法通過建立機械設備失效機理的數學模型[3-5]描述性能退化過程[6-8]。但是,由于類似火炮供輸機構的復雜機械設備,運動過程復雜、力學環境惡劣、性能退化機理復雜多變,難以建立準確的機理模型反映其健康狀態。

隨著傳感器技術的不斷發展,數據驅動的健康狀態評估方法成為研究熱點[9-11]。但是傳感器測得的振動信號數據往往只能反映機械設備某一方面的特征,并不能直接反映其健康狀態。湯健等[12]提出一種適用于旋轉機械設備狀態監測的軟測量技術,通過傳感器測量旋轉機械運行數據,并提取機械設備多方面的特征信息,構建這些特征信息的融合模型,從而得到能夠有效反映其健康狀態信息的軟測量模型。概括而言,軟測量技術主要包括特征提取[13]、建立軟測量模型[14]、模型評價[15]等部分。

目前的研究均是基于公共數據集,研究結果對復雜供輸機構的測量不具有普適性,潘宏俠等[16]采集火炮供輸機構振動信號,利用同步壓縮時頻變換分析得到反映不同運行狀態的二維時頻圖像,對圖像進行灰度化處理后,提取其紋理特征作為輸入,搭建基于模糊聚類的供輸機構軟測量模型,有效識別機構不同健康狀態。張航等[17]提出一種基于多分辨奇異值分解與灰色理論的火炮供輸機構軟測量模型,有效區分火炮供輸機構不同的健康狀態。席茂松等[18]對火炮供輸機構振動信號進行多尺度分解,提取分量信號信息熵作為特征,利用概率神經網絡搭建軟測量模型,有效識別供輸機構早期健康狀態變化情況。然而,上述研究重點在火炮供輸機構健康狀態識別,性能退化階段健康狀態的定量變化情況有待進一步研究。同時,在特征提取過程中需要依賴一定的專家經驗。

針對以上問題,本文擬采集供輸機構振動信號,提取振動信號時頻圖像,利用深度卷積網絡搭建軟測量模型,定量描述振動信號與供輸機構健康狀態關系,并在復雜供輸機構上進行試驗驗證,以期為供輸機構健康管理提供一定的理論參考。

1 軟測量理論介紹

1.1 連續小波變換

機械設備健康狀態退化過程中,振動加速度信號存在非平穩特征。時頻分析能夠有效表征非平穩信號時頻特征,可利用連續小波變換提取非平穩振動信號的特征,模型為

(1)

1.2 卷積神經網絡

典型的卷積神經網絡主要包括3層:卷積層、池化層和全連接層。首先,卷積神經網絡通過堆疊的卷積層和池化層提取原始圖片的二維特征;然后,將二維特征展平成一維特征向量;最后,將一維特征向量輸入到全連接層,得到輸出數值。

卷積層中,卷積核與前一層的特征圖進行卷積操作,然后將卷積結果輸入到激活函數,構建當前層的特征圖,即

(2)

利用展平層將通過堆疊卷積層和池化層構建的特征圖展平成一維特征向量。將該特征向量輸入到全連接層中,得到最終輸出值,即

(3)

1.3 基于變分推斷的模型不確定性估計理論

對于1.2節構建的卷積神經網絡,本文采用dropout正則化方法防止網絡訓練過擬合。假設網絡權重隨機變量ω服從高斯先驗分布,X表示給定的訓練數據,Y表示對應的輸出。通過貝葉斯理論可得隨機變量ω的后驗概率分布為

(4)

隨機變量ω包括深度學習網絡中所有可學習得到的參數,包括卷積層的內核和偏置項,全連接層的權重和偏置項。

(5)

式中:y*表示預測輸出值;x*表示網絡輸入值;q(ω) 為權重ω的取值概率;q*(ω)為使KL散度最小的數值。因此,不難得出,卷積神經網絡的優化目的就是使得KL散度最小,即相當于最大化證據下邊界。因此,卷積神經網絡的目標函數可以表示為

(6)

式中:xr和yr分別為第r個訓練數據和對應的輸出值。

(6)式右邊第1項可以采用蒙特卡洛積分法求解,可得

(7)

(8)

通過分析(7)式、(8)式可知,卷積神經網絡模型的輸出值yr隨采樣權重ω的變化,可以看做是隨機刪除權重矩陣中的某些行,這與深度學習網絡中的dropout過程是一致的。因此,深度學習網絡的模型的不確定性可以通過dropout算法進行計算。

網絡訓練過程中的目標函數為

(9)

目標函數Ldropout可以通過常見的優化器(如隨機梯度下降算法、Adam算法等)進行最小化,然后利用蒙特卡洛估計算法(也稱為蒙特卡洛dropout算法)進行預測。對于新的輸入樣本x*,蒙特卡洛dropout算法通過V折交叉驗證,可以得到預測結果的平均值和方差,模型為

(10)

2 復雜供輸機構健康狀態軟測量方法

2.1 復雜供輸機構健康狀態軟測量模型

復雜供輸機構振動加速度信號中蘊含豐富的健康狀態信息,但是其運行過程存在較大的沖擊、振動,導致采集的監測信號具有噪聲多、非線性、非平穩等特點,無法直接用于健康狀態評估。小波分析構建時頻圖像數據能夠有效反映振動信號中的瞬變特征。通過提取復雜供輸機構不同循環周期時頻圖,記錄機構不同運行階段的健康狀態。深度卷積網絡能夠有效提取圖像特征,通過搭建不同的軟測量模型分別實現復雜供輸機構健康狀態識別和性能退化程度定量評估。健康狀態識別模型首先將試驗數據分為正常運行階段監測數據與性能退化階段數據兩類,利用深度卷積網絡識別振動信號時頻圖特征,區分不同運行階段;性能退化程度定量評估模型是將試驗數據不同運行時刻的時頻圖作為輸入數據,相應的運行時刻在整個運行周期中的百分比作為標簽,利用深度卷積網絡提取不同運行時刻的時頻圖特征,建立時頻特征與運行時刻之間的映射關系。將測試數據時頻圖輸入至已訓練的網絡,即可輸出其對應運行時刻,進而實現健康狀態評估。文獻[20]表明,回歸預測過程中,深度學習網絡算法一個新點的輸出應大致與最接近的訓練點的輸出相同。因此,模型輸出結果即可認為是性能退化狀態軟測量值。在性能退化程度定量評估軟測量模型中,采用dropout技術隨機“丟棄”部分網絡層之間的連接,減緩過擬合現象對評估結果的影響,導致軟測量模型存在一定的不確定性,表現在評估結果中,即存在一定的置信區間。針對模型不確定性對軟測量結果的影響,利用變分推斷原理推導出軟測量結果不確定的定量表達式,有效表示了模型不確定與最終評估結果之間的定量關系。

時頻數據驅動的典型復雜供輸機構軟測量模型機構如圖1所示,評估健康狀態的步驟如下:

1)利用振動加速度傳感器采集供輸彈系統振動信號;

2)利用小波變換得到振動信號時頻圖;

3)將步驟2得到的時頻圖作為輸入信號,分別輸入到健康狀態識別軟測量模型和性能退化程度定量評估軟測量模型,對模型進行訓練;

4)將測試數據信號輸入到構建好的模型中,得到供輸彈系統健康狀態軟測量結果。

圖1 復雜供輸機構健康狀態軟測量模型Fig.1 Soft measurement model of complex feeding and ramming mechanism

2.2 基于多尺度卷積網絡的圖像特征提取

深度卷積網絡通常將最后一層池化層提取的特征輸入至全連接層進行回歸預測,該全局特征比上一層卷積層提取的特征更加穩定。但是,卷積層提取的特征經過池化層后,會丟失一些局部特征,導致特征提取不充分。針對該問題,本文提出一種多尺度深度卷積神經網絡結構,如圖2所示。多尺度卷積網絡是在卷積網絡結構基礎上,將最后一層池化層與上一層卷積層展平后相連接,組成一個特征向量,輸入到全連接層進行回歸預測。

圖2 多尺度卷積網絡結構圖Fig.2 Structure of multi-scale convolutional network

卷積網絡展平層的輸出為

(11)

經過改進后的卷積網絡,其多尺度展平層的輸出為

(12)

對比(11)式、(12)式可知,多尺度卷積網絡的多尺度展平層同時包含了池化層與卷積層提取的特征,增強了網絡特征提取能力,有利于反映滾動軸承早期性能退化狀態。

2.3 考慮模型不確定性的性能退化狀態定量評估

考慮深度學習模型不確定性的性能退化狀態定量評估方法流程如圖3所示。

3 供輸機構臺架試驗驗證

3.1 臺架試驗設計

某型復雜供輸機構試驗臺架如圖4所示。該試驗臺架包含動力裝置、試驗臺架體、控制裝置及擺動機構組成。試驗平臺采用移動式液壓站為動力,通過外能給復進機儲能并保持,在人工后坐油缸返回到位后,扳機釋放尾部鐵塊,尾部鐵塊在復進過程中為擺彈機儲能。

圖4 復雜供輸機構試驗臺架Fig.4 Test rig for complex feeding and ramming mechanism

尾部鐵塊在滾輪帶動下向后移動,擺動機構帶動擺臂由豎直方向擺動90°,使擺臂與地面平行。滾輪健康狀態好壞將直接影響尾部鐵塊運行能否到位,進而影響擺臂能否正常擺動到位。由于滾輪在每次擺臂過程中受到一定的載重力,且反復運用,磨損較大,其健康狀態將直接影響復雜供輸機構能否正常運行,一定程度上反映了復雜供輸機構的健康狀態。由于滾輪位于凹槽內,位置隱蔽,空間狹小,無法直接安裝傳感器(見圖5)。而擺動機滑板、壓板機在滾輪運動過程中處于滾輪正上方且距離滾輪較近,便于布置傳感器,因此在擺動機滑板、壓板機附近布置6個振動加速度傳感器,分別為測點1~測點6,如圖6所示。

圖5 供輸機構滾輪及滑槽Fig.5 Roller and spout for feeding and ramming mechanism

圖6 傳感器布置Fig.6 Sensor layout

傳感器類型為壓電集成電路加速度傳感器,采樣頻率為10 kHz,采用德國西門子公司生產的32通道LMS信號采集系統。本文試驗分為2組:第1組采用健康狀態的滾輪循環運行11組、每組循環包括20次動作;第2組換裝性能退化階段的滾輪循環運行9組,每組循環包括20次動作。

試驗后期,滾輪出現較大裂紋,為安全起見,停止試驗,認為滾輪已經達到壽命閾值。性能退化階段的滾輪和試驗后出現較大裂紋的滾輪如圖7所示。

圖7 性能退化階段及出現較大裂紋的滾輪Fig.7 Performance degradation stage and large cracks in the roller

3.2 臺架試驗數據采集

選取位于滾輪正上方測點1處的振動傳感器數據進行分析。健康狀態和性能退化狀態下測得的振動加速度信號分別如圖8(a)和圖8(b)所示。由圖8可知,供輸系統測得的振動加速度信號中有較大的沖擊振動。同時,僅從原始振動加速度信號來看,無法有效區分健康狀態和性能退化階段。

圖8 原始振動加速度信號Fig.8 Original vibration acceleration signal

供輸系統試驗平臺采用自動控制系統控制每個動作的運行時間,一個循環周期的時間大約為31.2 s。提取其中一個運行周期的振動加速度信號,如圖9所示。該循環過程主要包括4個運動過程:擺動機構儲能、扳機解脫、拉鉤解脫、擺動到位。通過圖9可知,該循環過程的振動加速度信號屬于非平穩信號,頻譜分析難以得到系統性能退化的特征。

圖9 一個循環的振動加速度信號Fig.9 A cyclic vibration acceleration signal

基于(1)式可得,供輸系統健康狀態和性能退化狀態的時頻圖如圖10所示。

圖10 一個循環的時頻圖Fig.10 Time-frequency map of a cycle

觀察圖10可知,供輸系統振動信號的能量幾乎分布在所有頻段,且中低頻段內的能量高于高頻段內的能量。在一個循環周期內,出現明顯的能量波動,信號的瞬變特征明顯。對比健康狀態和性能退化階段的時頻圖可知,性能退化階段的信號能量高于健康狀態的信號能量,同時二者瞬變過程也存在差異。

提取性能退化階段不同時期的時頻圖,如圖11所示。由圖11可知,性能退化階段的時頻圖變化并不明顯,無法直接判斷性能退化程度。因此,需要通過卷積網絡提取深層特征,定量表征性能退化程度。

圖11 性能退化階段不同時期的時頻圖Fig.11 Time-frequency maps of different degradation periods

健康狀態和性能退化階段的振動信號時頻圖共400組數據,將其隨機打亂后選擇其中的300組數據作為訓練數據,100組數據作為測試數據。

3.3 試驗結果及分析

由于學習過程中,網絡結構的增加、學習率的大小、卷積核寬度的增減、Dropout值的高低均可能造成模型過擬合。目前,還沒有統一的確定超參數的理論。本文綜合考慮試驗數據量大小,輸入的時頻圖大小,采用手動調整超參數的方法,確定使得驗證誤差相對較小的超參數值,如表1所示。

表1 網絡超參數

為了進一步分析不同網絡層學習得到的特征,將測試集中的時頻圖輸入到已訓練的網絡中,分別提取第1層和第2層卷積層學習到的圖像特征,如圖12所示。

圖12 不同卷積層提取的特征圖Fig.12 Feature graphs extracted by different convolutional layers

由圖12可知,隨著卷積層數的增加,網絡提取的特征愈加抽象。由第1層卷積層提取的特征圖2、特征圖6、特征圖7、特征圖9、特征圖10、特征圖11可知,性能退化階段的特征大約出現在一個循環的22 s附近,此時供輸系統進行擺彈動作,性能退化的滾輪在凹槽內滑動,與實際情況相符合,說明網絡學習到了性能退化階段滾輪運行的時間特征;通過第2層卷積層輸出的12個特征圖可知,能量較高的模塊大都集中在中低頻段,可知表征滾輪性能退化階段的頻率也大都集中在中低頻段。由此可見,網絡能夠有效學習供輸系統健康狀態與性能退化階段的頻譜特征。

為了驗證本文所提方法的有效性,與兩種常用的識別方法進行對比。識別準確率曲線如圖13所示。第1類信息熵結合神經網絡的識別方法是指提取振動信號時域信息熵作為輸入,搭建基于神經網絡的軟測量模型,輸出健康狀態識別結果(簡稱信息熵+神經網絡方法)。第2類方法是提取振動信號時域均方根(RMS)作為輸入,搭建基于神經網絡的軟測量模型,輸出健康狀態識別結果(簡稱RMS+神經網絡方法)。由圖13可知,本文所提方法在迭代訓練80次后識別準確率可達90%,信息熵+神經網絡方法在迭代訓練50次后,識別準確率基本穩定在70%,RMS+神經網絡方法識別準確率較低,僅有55%。

圖13 不同方法的識別準確率Fig.13 Recognition accuracies of different methods

當供輸系統進入性能退化階段時,本文提出一種基于深度卷積神經網絡的回歸預測模型,能夠定量測量性能退化狀態,進而為下一步的剩余壽命預測奠定基礎。本文搭建的用于回歸預測的深度卷積神經網絡是在分類深度卷積神經網絡的基礎上,將最后的二分類層改為回歸預測層,輸出節點數為1。供輸系統進入性能退化階段后是一個循序漸進的過程,即從性能退化開始一直到滾輪出現裂紋無法工作的過程是單調不可逆的。將供輸系統性能退化階段采集到的振動加速度信號轉換成時頻圖后作為深度卷積神經網絡的輸入,可表示為It={xt,yt},其中:xt為t時刻一個循環內尺寸為H×H的時頻圖,xt∈RN×N,H為整數,表示時頻的尺寸,本文采用的時頻圖尺寸為100×100;yt為t時刻循環對應的供輸系統性能退化程度,yt∈[0,1],yt=t/T,T表示總的運行時間。通過多次試驗可知訓練數據標簽yt與循環時刻t的關系并不影響最終的健康因子構建結果。因此,在此選擇線性模型構造訓練數據標簽。此處,供輸系統共運行了400個循環,則第200個循環對應的退化程度yt=200/400=0.5。將所有數據隨機打亂后分為訓練組與測試組,利用訓練集對深度卷積網絡進行訓練,將測試數據輸入到已訓練的深度卷積網絡,得到性能退化定量測量結果。將所有數據輸入到訓練好的網絡中,得到供輸系統性能退化曲線,如圖14所示。由圖14可知,網絡能夠根據設定標簽有效學習到滾輪性能退化過程。將測試數據輸入到訓練好的網絡中,得到性能退化程度預測值。

圖14 供輸系統性能退化曲線Fig.14 Performance degradation curve of feeding and ramming mechanism

利用預測結果與真實結果的相對誤差作為評估預測精度的準則,計算公式為

(13)

式中:Erη表示第η次循環得到的預測精度值;ActDgredη表示第η次循環實際的性能退化程度;Dgredη表示第η次循環軟測量模型預測的性能退化程度。

將本文所提方法的軟測量精度與其他3類方法進行對比分析,結果如圖15所示。RMS+神經網絡方法表示提取RMS作為特征,利用神經網絡預測性能退化程度;信息熵+神經網絡方法表示提取信息熵作為特征,利用神經網絡預測性能退化程度;時頻圖+神經網絡方法表示提取時頻圖作為特征,利用神經網絡預測性能退化程度。時頻圖+神經網絡方法來自文獻[21],其卷積網絡結與本文方法不同,沒有設置dropout層。由圖15可知,本文方法的相對誤差最小,即測量精度最高。同時,測量結果誤差范圍相對較小,表明該方法具有一定的穩定性。對比時頻圖+神經網絡方法可知,沒有采用dropout層的卷積網絡測量精度與本文所提方法接近,但是測量結果誤差范圍較大,同時,出現兩個異常點,可能的原因是網絡出現一定程度的過擬合。

圖15 不同方法的精度對比圖Fig.15 Accuracy comparison of different soft measurement methods

為了驗證dropout正則化項是否能夠緩解過擬合現象,將實驗數據隨機分為兩組,其中350組數據作為訓練集,50組數據作為驗證集,將dropout值分別設為0.1、0.2、0.3、0.4,計算訓練集、驗證集的損失函數值,結果如表2所示,當dropout取值為0時,訓練集損失函數值最小,擬合度最高,但是在驗證集上的損失函數值最大,說明過擬合現象最嚴重;隨著dropout值的增加,訓練集上的損失函數值先增加、后平穩變化,而驗證集上的損失函數值則先減小、后增加。訓練集和測試集損失函數差值先減小、后增加。當dropout值為0.2時,其在驗證集上的損失函數值最小,且訓練集與驗證集損失函數差值最小,說明添加dropout正則化項后一定程度上緩解了過擬合問題,且不同dropout值對過擬合問題的緩解程度存在差異。

表2 不同dropout值下的損失函數值

為了進一步評估不同dropout取值對模型不確定性的影響,以第49個測試數據為例,dropout取值分別為0.1、0.2、0.3、0.4時的概率分布如表3和圖16所示。圖16中,4組曲線對應的均值分別為0.602、0.610、0.511、0.562,這4組值是將第49個測試數據的時頻圖輸入到訓練好的模型中,模型對應的輸出值,即為軟測量結果。總的試驗循環次數為400次,因此對應的軟測量循環次數依次240.8、244.0、224.8、204.4,實際的循環次數應為245.2。由于循環次數應為整數,四舍五入后的軟測量結果分別為241、244、224、204。將軟測量結果代入(13)式中,即可計算得到對應的軟測量精度分別為1.63%、0.41%、8.57%、16.73%。由此可知,隨著dropout的取值增加,軟測量精度先是逐漸增加,之后逐漸減小。當dropout值由0.1變為0.2時,軟測量精度增加,同時數據分布方差變小,即模型不確定性降低。以上結果表明,當dropout值為0.2時,模型容量與數據最匹配,得到的軟測量精度最高,模型不確定性相對較小。由表3和圖16可知,隨著dropout值的增加,模型不確定性不斷增大,測量結果的方差增加,即結果不確定性增加。

表3 不同dropout值的測量結果

圖16 軟測量結果概率密度曲線Fig.16 Probability density curve of soft measured results

4 結論

為解決復雜力學環境下工作的供輸機構健康狀態評估問題,本文提出了一種基于時頻圖和深度卷積網絡的健康狀態軟測量方法。通過morlet小波變換提取振動信號時頻圖作為輸入特征,構建了基于深度卷積網絡的軟測量模型,并引入dropout正則化層,實現了供輸機構健康狀態識別和性能退化程度定量評估。利用典型供輸機構臺架試驗驗證了所提方法的有效性,與不同方法的對比分析表明,基于時頻圖和深度卷積網絡的軟測量方法能夠有效提高供輸機構健康狀態識別準確率和性能退化程度軟測量精度,降低測量結果的不確定性。本文研究方法對復雜供輸機構健康狀態評估提供了一種新的思路,也為供輸機構健康管理提供一定的理論參考。

猜你喜歡
振動特征測量
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
如何表達“特征”
把握四個“三” 測量變簡單
不忠誠的四個特征
當代陜西(2019年10期)2019-06-03 10:12:04
滑動摩擦力的測量和計算
滑動摩擦力的測量與計算
中立型Emden-Fowler微分方程的振動性
抓住特征巧觀察
測量
主站蜘蛛池模板: 久久香蕉国产线看精品| 伊人婷婷色香五月综合缴缴情| 国产嫖妓91东北老熟女久久一| 国产麻豆福利av在线播放| 日韩在线第三页| 成人亚洲天堂| 国产精品播放| 试看120秒男女啪啪免费| 国产91蝌蚪窝| 中文无码精品A∨在线观看不卡| 美女高潮全身流白浆福利区| 91口爆吞精国产对白第三集| 国产无码精品在线播放| 国产一级α片| 日韩区欧美国产区在线观看| 亚洲精品大秀视频| 日韩欧美国产精品| 色综合综合网| 中国一级特黄大片在线观看| 国产真实自在自线免费精品| 国产精品自拍露脸视频| 国产十八禁在线观看免费| 日韩精品成人网页视频在线| 精品三级在线| 91成人在线免费视频| 国产精品私拍在线爆乳| 国产91久久久久久| 四虎国产精品永久一区| 全部毛片免费看| 国产欧美视频在线观看| 国产激爽大片在线播放| 久久99久久无码毛片一区二区| 亚洲欧洲自拍拍偷午夜色| 国产精品亚洲а∨天堂免下载| 色综合a怡红院怡红院首页| julia中文字幕久久亚洲| 欧美在线精品怡红院| 思思热在线视频精品| 亚洲久悠悠色悠在线播放| 免费一级毛片不卡在线播放| 欧美亚洲一区二区三区在线| 人人妻人人澡人人爽欧美一区| 秋霞午夜国产精品成人片| 日本在线亚洲| 亚洲日韩精品无码专区| 亚洲第一成年网| 国产香蕉97碰碰视频VA碰碰看| 精品乱码久久久久久久| 国产区免费| 色婷婷综合激情视频免费看| 无码乱人伦一区二区亚洲一| 国内精品久久久久久久久久影视| 真人高潮娇喘嗯啊在线观看 | 久久精品无码中文字幕| 色香蕉影院| 九九免费观看全部免费视频| 国产xxxxx免费视频| 久久久久久国产精品mv| 亚洲大学生视频在线播放| 一级做a爰片久久毛片毛片| 国产香蕉在线| 欧美成人综合在线| 亚洲第一区欧美国产综合| 国内精品久久九九国产精品| 亚洲综合欧美在线一区在线播放| 精品国产成人国产在线| 久久久久人妻一区精品| 久久99精品久久久久久不卡| 亚洲视频在线青青| 亚洲天堂视频网| 欧美成人区| 无码国产偷倩在线播放老年人| 精品无码一区二区在线观看| 日韩毛片免费观看| 四虎永久免费在线| 国产午夜人做人免费视频中文| 日本三级欧美三级| 天天色综网| 色婷婷综合激情视频免费看| 色婷婷成人| 久视频免费精品6| 久久黄色免费电影|