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

時頻流形自適應稀疏重構的遙測振動信號特征增強方法

2022-04-04 18:55:38劉學孫翱李冬黃銳
振動工程學報 2022年1期

劉學 孫翱 李冬 黃銳

摘要:針對遙測振動信號非線性、非平穩性、瞬態沖擊性等特點,提出一種基于時頻流形自適應稀疏重構的遙測振動信號特征增強方法,對振動信號進行相空間重構提取其時頻流形;以時頻流形為基礎,采用KSVD算法自適應構建過完備字典,并從中找到最匹配的時頻原子,根據得到的原子與相空間展開信號的時頻分布,依次匹配計算獲得其重構的稀疏系數;利用稀疏系數和時頻原子對相空間中各維信號的時頻分布進行重構,通過時頻分布的逆運算和相空間還原得到特征增強信號。仿真和實測信號處理結果驗證了算法的有效性。

關鍵詞:信號處理;遙測振動信號;時頻流形;稀疏重構;特征增強

中圖分類號:TN911.7

文獻標志碼:A

文章編號:10044523( 2022)01-024609

DOI: 10.1638 5/j .cnki.issn.10044523.2022.01.027

引 言

遙測振動信號是由安裝在被試飛行器內的振動加速度或位移、溫度、壓力等傳感器采集的反映系統運行狀態的時間序列。受飛行器本身振動、飛行環境、電磁環境和傳輸條件等因素的影響,各結構部位的振動會相互影響、調制和疊加,振動的傳輸路徑復雜多變,采集的遙測振動信號往往夾雜著大量的高頻、低頻和沖擊噪聲、各階次的諧波分量,頻譜成分異常復雜,且相關性耦合程度高,呈現較強的非線性和非平穩性[1]。如何遙測振動信號在不改變其所反映的系統動力學特征情況下進行特征增強直接關系到飛行狀態分析的準確性。

信號的稀疏分解重構降噪一直是非線性信號降噪領域的研究熱點,傳統的特征增強分析往往從單一的時域或頻域構造一系列瞬態脈沖原子對信號進行稀疏分解,如Cui等[2]利用振動信號進行Gabor變換生成一組脈沖原子實現對信號的稀疏分析;Wang等[3]采用Morlet小波作為時域原子構造過完備字典,采用正交匹配追蹤( Orthogonal Matching Pur-suit,OMP)對旋轉機械信號進行稀疏分解;Fan等[4]采用Laplace小波變換結合譜分析技術在時域實現對軸承信號的稀疏表示。但是實際系統工況多變,沖擊過程復雜,如果只用單一的、確定性的、公式化的原子很難表征實際的信號特征,最近,為提高對信號的表征能力,很多改進的變公式原子模型相繼被提出,如連續雙Laplace小波原子[5]、雙邊非對稱Morlet小波原子[6-7]、變窗寬Gabor小波等[8],根據先驗信息不斷改變原子公式模型去逼近真實信號沖擊情況,取得了一定的效果,但是這種變公式模型存在固有缺陷。需要有先驗信息,但在實際情況下很難提前獲知沖擊過程,且實際信號瞬態特征波形大多為非公式化的,這就需要搜索優化更多的公式化原子模型參數,使得稀疏過程變得更加復雜。

另外,經He等[9-10]研究發現,傳統的時域稀疏方法對信號中的相位信息具有較差的匹配效果,導致恢復出的結果較原信號存在較大的誤差,由于時頻分布綜合了時間域與頻率域的聯合分布信息,可以有效描述動態系統中的非平穩特性。因此,He等[11]采用時頻域重構可以有效解決相位匹配的問題,并可獲取更高的重構精度,但該文獻依然采用公式化原子的時頻分布構造時頻原子對原信號的時頻分布進行稀疏匹配,不可避免還會受到公式化模型固有缺陷的限制,為了解決上述問題,本文提出一種基于時頻流形( Time-Frequency Manifold,TFM)自適應稀疏重構的信號特征增強方法,采用KSVD算法[12]從時頻流形中自適應構建過完備字典,并從中找到最匹配的時頻原子,避免對原子庫模型構造的依賴性,利用保持時頻分布的相位信息,有效提高了信號瞬態特征增強的準確性以及稀疏方法的通用性。

1 時頻流形

時頻流形是嵌入在非平穩信號時頻分布中的一種內在的非線性流形結構,采用流形學習從信號相空間重構分量的時頻分布(高維空間)提取到該信號嵌入低維空間的特征結構[9]。在流形學習過程中綜合了信號本身的非平穩性和非線性兩種信息,因此可有效地挖掘和表征信號的時頻模態,刻畫信號時頻分布特性。

1.1 相空間重構和時頻分布

采用相空間重構的方法將振動信號擴展到多維相空間,使混疊的各源信號在多維空間中重新展開,以確保滿足Takens定理要求和尺度空間的統一性。對于振動信號x(t)=[x1,x2,...XN],相點向量可通過下式重構得到:

1.2

LTSA時頻流形提取

局部切空間排列算法(Local Tangent SpaceAlignment,LTSA)基于數據流形滿足局部線性條件假設,即高維數據空間的局部區域和局部切空間之間存在著一個線性投影,同時全局低維空間的局部區域和局部切空間之間也存在著一個線性投影[15]。LTSA通過逼近每一樣本點的切空間來構建低維流形的局部幾何結構,觀測數據點在局部切空間的投影獲得局部低維坐標,由局部仿射變換而得到低維的全局坐標[16]。

采用上述LTSA算法提取時頻流形的主要思想:在高維的重構相空間中,相對于嵌入維數而言,系統有用信號吸引子所在的主流形是低維的,其分布只局限于相空間中某個低維的子空間內,而白噪聲則在相空間的所有維度中都有分布。根據信號和噪聲分布的不同,通過LTSA流形學習,在時間序列的嵌入相空間降維過程中,保留系統的主流形結構[9]。但LTSA算法的輸人為兩維,而相空間m維子序列的時頻分布矩陣均為三維矩陣,因此要把它們轉換成兩維矩陣,以時間軸為基準,將幅值矩陣Aj(k,廠)的每一列首尾相連,構造一個一維向量TjFFD,將m個一維向量組成兩維數據集XTFD為:

2 時頻流形自適應稀疏重構

稀疏重構就是使用過完備字典對信號進行稀疏表示,將信號的能量集中在字典的少量原子上,通過少量原子來表征信號的本質特征。從數學和信號的角度,大多數零分量或小幅度分量和少數非零大幅度分量揭示了信號的內在結構和本質屬性,由于遙測振動信號具有瞬態沖擊性,滿足稀疏分解條件,同時采用時頻流形可以有效提取信號的瞬態特征,也能消除帶內噪聲,但是受流形學習算法式(6)的限制,提取的時頻流形特征幅值損失嚴重,為了在去噪的同時還原信號的幅值信息,以時頻流形為基礎,自適應提取時頻原子,結合相位信息保持,通過學習到時頻原子對信號進行稀疏重構,恢復幅值信息,對信號特征進行增強。

2.1 基于圖像塊的稀疏域建模式中 μ為懲罰因子。式(9)屬于NP難求解問題,理論上需要轉化為1范數才能進行多項式解析求解,但依然求解困難,通常采用稀疏表達近似求解方法,主要分為兩大類,一是匹配追蹤算法,通過構造一系列具有稀疏表達能力的原子字典,通過相關匹配分析,學習出信號中的主要特征成分,采用的是貪婪迭代方法;另一種是基追蹤算法,通過稀疏約束在全局范圍內的極值問題,使用線性規劃方法來實現對信號稀疏分解,采用的是優化重構。本文采用正交匹配追蹤OMP算法進行近似求解。在對時頻流形分布T的每一個圖像塊建立稀疏域模型的基礎上,對T進行貝葉斯最大后驗估計,得到:

2.2 自適應時頻原子字典構造

為了克服公式化模型原子字典固有缺陷的限制,本節采用K-SVD算法自適應訓練時頻原子字典,基本思想是利用時頻流形瞬態特征提取與帶內噪聲去除的能力,采用DCT字典或從時頻流形分布中隨機抽取冗余時頻原子作為初始化的過完備字典D,使用OMP算法求解NP問題得到稀疏系數aij的近似優化解,然后根據稀疏系數aij使用K-SVD算法[12.17]依次對字典D中的每一個原子進行更新。算法流程如下:式中 Ei表示第Z次迭代后的殘差能量,E。表示原信號能量,通過實測遙測振動信號驗證一般K≤7,式(11)即可保證稀疏主成分被提取完畢。因此使用式(11)作為匹配追蹤算法的迭代終止條件,對每一個圖像塊Rij進行上述的迭代優化得到其對應的稀疏系數aij;

2.3 相空間信號時頻分布稀疏重構

根據上節訓練得到的基于時頻流形的過完備字典D,則模型求解問題即可等價于求解局部稀疏系數aij和重構時頻分布T,可分兩步求解,第一步為局部圖像塊稀疏匹配,令T=T,采用塊系數最小化方法計算aij的最優值,則式(10)可轉換為

采用OMP算法對式(13)進行求解,同2.2節一樣依然使用式(11)作為匹配追蹤算法的迭代終止條件;

3 時頻流形自適應稀疏重構特征增強方法

稀疏重構方法的降噪能力與原子選取有關,通常采用公式化的原子創建過完備字典,但在實際情況下信號瞬態特征波形大多為非公式化的,且很難提前獲知沖擊過程以及相位變化等先驗信息,這就導致匹配追蹤時對原子形態要求過高,因此,公式化原子字典很難在實際工況下對原信號進行較好的稀疏逼近。另外,如果直接在時頻域上做匹配追蹤稀疏分析,雖然可以保持信號的相位信息,使得恢復出的信號與原信號具有相同的波形特征,但受噪聲的影響(特別是帶內噪聲),匹配稀疏的特征表達能力將受到嚴重的削弱。

針對這些問題,本文將時頻流形與稀疏分析相結合,充分利用時頻流形良好的瞬態特征提取與帶內噪聲去除的能力,以及稀疏表示良好的特征表達能力,從時頻流形中自適應構建過完備字典,并從中找到最匹配的時頻原子,在保持原信號時頻分布相位信息的同時,采用時頻原子匹配稀疏的方法,對相空間中含噪信號時頻的分布進行稀疏表達,實現沖擊特征的提取和信號幅值的恢復。算法流程如圖1所示。算法分為基于時頻流形的白適應字典學習和相空間信號時頻分布的稀疏重構兩個部分。

(1)首先對采集的遙測振動信號依據指令時刻進行特征段順序選取,若相鄰指令時刻間隔較近(前一指令響應未結束,后一指令響應即開始),可按前一特征段最大幅值能量的10%進行截取;若重疊部分超過50%,則將相鄰兩特征段合并處理。對所選特征段進行預處理:根據《GJB2238A-2004》的規范進行零漂修正、趨勢項去除、野值剔除等;

(2)對特征段信號x(t)采用式(l)和(2)進行相空間重構,對重構相空間每一維子序列Pi(j=1,2,…,m),采用式(3)計算其時頻復數矩陣Si(k,f),將幅值矩陣A,(k,f)和相位矩陣θ,(k,f)進行分離;

(3)對由m維相空間子序列時頻分布的Aj(k,f)的幅值矩陣組成高維矩陣進行LTSA主流形提取,得到時頻流形分布T(k,f);

(4)采用2.2節的方法從T(k,f)白適應訓練過完備字典D;

(5)利用學習得到的過完備字典D,采用式(12)~(15)對m維相空間子序列時頻分布的A,(k,廠)進行稀疏重構,得到重構后的幅值矩陣T1(k,f),結合步驟(2)得到的相位矩陣θ,(k,f)對時頻復數矩陣進行重新合成,得到Sj(k,f);

(6)對Si(k,f)進行STFT逆變換更新相空間每一維子序列Px(j=1,2,…,m);

(7)對相空間矩陣進行還原得到去噪信號y(t)及采用式(3)計算其時頻分布。

式中N為信號長度,集合Ii為包含信號元素i的下標集,Ci為Ii中元素個數。

4 仿真與實測信號分析

采用仿真和實測信號對本文所提方法( Time-Frequency Manifold Adaptive Sparse Reconstruc-tion,TFMASR)同基于時域稀疏重構(Time Do-main Sparse Reconstruction,TDSR)、時頻域稀疏重構方法( Time-Frequency Domain Sparse Reconstruc-tion,TFDSR)進行性能對比驗證,為實現對瞬時脈沖特征的挖掘,時域稀疏重構采用雙尺度Gabor變換構造時域原子[8],時頻域稀疏重構采用時域原子的STFT變換得到的時頻分布作為初始時頻原子[11],各算法均采用OMP匹配追蹤算法進行稀疏重構,OMP算法參數有稀疏度、最大迭代次數和迭代終止條件。稀疏度可默認為字典的維數,即使得所有字典遍歷一次,最大迭代次數設置為100,迭代終止條件如式(13)所示,殘差能量的變化率閾值ε設置為0.005。經實驗驗證,一般不超過50次即可滿足迭代終止條件,跳出循環。結果采用信噪比( Signal to Noise Ratio,SNR)和衡量時頻域上瞬態脈沖的敏感特征特性一聯合時頻熵[11]( Joint Time-Frequency Entropy,JTFE)作為量化指標。算法的運行平臺:Inter Core i7-4790(主頻3.6 GHz)CPU,8GBDDR3內存,Matlab 2015b,Windows 7 64位專業版操作系統。

4.1 仿真分析

驅動頻率f=l kHz,采樣頻率f=1 kHz,信號長度N- 1024,A=[0.76,0.96,0.92.1.09]和θ=[ π/6,π/4,π/3,π/2]分別為給定的幅值和初始相位向量,阻尼系數ξ=0.01,r=[0.02,0.04,0.06,0.08]為脈沖起始時刻,η=0.02 s為脈沖持續時間,n(t)為加入-8 dB的高斯白噪聲,結果如表1和圖2~7所示。

圖2和3分別給出了仿真信號和加噪信號的時域波形、頻譜和時頻分布。通過對比發現,加噪信號中的4個瞬態沖擊成分被噪聲嚴重污染,導致瞬態特征被削弱。圖4給出了時域稀疏重構方法獲得的結果圖。從表1給出的輸出信噪比看僅為-5.19dB,提升不明顯,重構出的結果與原信號存在較大的誤差,這說明在初始相位不同的情況下,相位匹配在一定程度上會影響時域稀疏重構算法的性能,由于稀疏重構只是采用少數原子去對原信號進行稀疏表達,噪聲將嚴重降低(消弱)其匹配稀疏的特征表達能力,時域稀疏僅在時域內進行匹配追蹤,這就對時域原子的形態要求過高,很難對幅值和相位同時匹配以保證重構波形不失真,這也驗證了文獻[11]的結論。

圖5為時頻域稀疏重構方法的結果圖,可以看出重構的結果遠比時域重構的結果要好,輸出信噪比為-3.82,這說明時頻域稀疏重構方法采用時頻分析(STFT變換)對幅值和相位進行分離,通過相位保持降低了對原子庫模型構造的依賴性,與仿真信號的波形特點更為吻合,保證了瞬態特征提取的準確性,但時頻域稀疏重構方法將噪聲也恢復出來,瞬態特征增強效果不強;圖6為加噪信號的時頻流形分布圖,從圖中可以看出,帶內噪聲基本被去除,具有良好的瞬態特征提取能力,但受LTSA流形學習算法公式(7)的約束,提取的時頻流形特征對比原信號幅值損失嚴重,最高幅值僅為0.037,相差約2個數量級,這將嚴重影響后續的信號分析;

圖7給出了TFMASR方法的輸出結果圖,可以明顯可以看出,基于時頻流形稀疏重構方法具有高效準確的瞬態特征增強能力。表1給出輸出信噪比為6.21 dB,JTFE為最小的0.7655,反映出輸出信號的頻帶能量分布更為集中,且具有更少的帶內噪聲,這表明所提方法利用了時頻流形帶內噪聲去除能力,降低了噪聲對時頻原子構造學習的影響,通過相位保持對信號波形特點增強的準確性,極大地提高了瞬態沖擊特征的挖掘與增強能力。

4.2 實測信號驗證

為驗證所提方法的有效性,采用某型飛行器試驗任務中高頻振動傳感器采集的遙測振動信號進行處理驗證,采樣頻率為5.12 kHz,結果如圖8~12所示。由于SNR的計算需要無噪信號的先驗信息,因此其在仿真研究中可以作為較好的量化指標,但對實測信號卻無法使用。因此對于實測遙測振動信號,只采用JTFE對時頻能量分布情況進行量化。

圖8給出了實測遙測振動信號的時域波形、頻譜和時頻分布,驅動頻率大致在1600~2000 Hz之間,噪聲覆蓋整個時頻平面,信號瞬態特征被嚴重削弱,且這種瞬態特征呈現較強的非線性,波形特征也隨時變化。圖9給出了時域稀疏重構的結果,可以看出時域稀疏雖然可以去除一部分噪聲,但帶內噪聲去除得不夠明顯,JTFE=0.7633時,由于沒有對相位信息進行保持,同上一節的仿真結果一樣,重構出的信號波形較原信號有一定的差異。

圖10給出了時頻域稀疏重構的結果圖,采用時頻分析技術將幅值和相位分離,通過相位保持使得恢復出的信號同原信號的波形特點較為吻合,在一定程度上提升了信號瞬態特征的稀疏表達效果,僅當JTFE= 0.7697時,反而比時域稀疏重構方法大,這說明受噪聲的影響,特別是帶內噪聲,使得在時頻原子構造過程中添加了噪聲信息,在重構過程中難免將噪聲信息也一并恢復,使得降噪效果不佳;圖11給出了實測信號的時頻流形分布圖,可以看出其具有良好的帶內噪聲去除效果,但提取的時頻流形的幅值受LTSA算法數據中心化的影響有較大衰減,最高幅值僅為0.045,較原信號幅值28.24相差約3個數量級,這將嚴重影響后續信號分析的準確性;圖12為基于時頻流形自適應稀疏重構方法的輸出結果圖,不難發現,該方法具有最小的JTFE-0.7418,重構出的信號時頻分布具有最高的時頻能量聚集性,瞬態特征增強效果最好。這說明從時頻流形中白適應學習構造時頻字典充分利用時頻流形的帶內噪聲去除能力,降低了噪聲對構造原子質量的影響,因此,本文所提方法在波形特點保持、帶內噪聲去除、瞬態特征增強以及降低原子庫模型構造的依賴性等方面較時域稀疏重構和時頻域稀疏重構方法具有更好的效果。但從表1和2各算法運行時間對比來看,TFMASR方法的主要缺點是運行時間較長,1024點的時序仿真信號需要539.56 s的計算時間,經算法流程分析和實驗得出,算法的計算量主要集中在時頻流形的提取上,由于相空間信號時頻分布樣本數據量較大,導致后續的LTSA非線性流形學習效率降低,計算損耗增加;另外,TFMASR在字典白適應學習和相空間信號時頻分布需要2次的全局匹配追蹤搜索,也需要一定的時間損耗,為了提高算法的運行效率,可采取只對感興趣的分析頻段、時間段進行時頻流形提取和匹配追蹤搜索,降低時頻分布樣本數據量和提高匹配搜索的效率。

5 結 論

針對遙測振動信號非線性、強噪聲、瞬態沖擊性等特點,提出一種基于時頻流形自適應稀疏重構的信號特征增強方法,以時頻流形為基礎自適應構建過完備字典,避免對原子庫模型構造的依賴性;利用保持時頻分布的相位信息,有效提高了信號瞬態特征增強的準確性以及稀疏方法的通用性。仿真和實測信號實驗結果表明新方法在波形特點保持、帶內噪聲去除、瞬態特征增強以及降低原子庫模型構造的依賴性等方面優于時域稀疏重構和時頻域稀疏重構方法。

參考文獻:

[1] 劉學.虛部噪聲輔助LCD方法及其在遙測振動信號處理中的應用[J].振動與沖擊,2017,36(12):1-6.

Liu Xue. An image noise assisted local characteristicscale decomposition method and its application in telem-etry vibration signal processing[J].Journal of Vibrationand Shock, 2017, 36( 12): 1-6.

[2]

Cui L,Wang J,Lee S C.Matching pursuit of an adap-tive impulse dictionary for bearing fault diagnosis[Jl.Journal of Sound and Vibration, 2014, 333: 28402862.

[3]

Wang S,Cai G, Zhu Z,et al.Transient signal analysisbased on LevenbergMarquardt method for fault featureextraction of rotating machines[J].Mechanical Systemsand Signal Processing, 2015, 54: 1640.

[4]

Fan W, Cai G, Zhu Z K, et al.Sparse representation oftransients in wavelet basis and its application in gearboxfault feature extraction[J].Mechanical Systems and Sig-nal Processing, 2015, 56: 230245.

[5]Tang H, Chen J, Dong G. Sparse representation basedlatent components analysis for machinery weak fault de-tection [ J] . Mechanical Systems and Signal Processing ,2014, 46(2) :373388.

[6]Lin J, Qu L. Feature extraction based on Morlet wave-let and its application for mechanical fault diagnosis [Jl.Journal of Sound and Vibration, 2000, 234 (1) :135-148.

[7]Cui L, Wu N, Wang W, et al. Sensorbased vibrationsignal feature extraction using an improved compositedictionary matching pursuit algorithm [Jl. Sensors,2014, 14(9) : 16715-16739.

[8]Feng Z P, Liang M. Complex signal analysis for plane-tary gearbox fault diagnosis via shift invariant dictionarylearning[ J] . Measurement , 2016 . 90 : 382395.

[9]He Q, Liu Y, Long Q, et al. Timefrequency manifoldas a signature for machine health diagnosis [Jl. IEEETransactions on Instrumentation and Measurement,2012 , 61( 5) : 1218-1230.

[10]He Q. Timefrequency manifold for nonlinear feature ex-traction in machinery fault diagnosis [Jl. MechanicalSystems and Signal Processing , 2013 , 35( 2) : 200218.

[11]He Q B, Ding X X. Sparse representation based on lo-cal timefrequency template matching for bearing tran-sient fault feature extraction [Jl. Journal of Sound andVibration , 2016 , 370 : 424443.

[12]Ahron M, Elad M, Bruckstein A M. KSVD: an algo-rithm for designing overcomplete dictionaries for sparserepresentation[J]. IEEE Transactions on Instrumenta-tion and Measurement, 2006, 15( 11): 4311-4322.

[13]T akens F.Detecting strange attractors in turbulence[J].Lecture Notes in Mathematics, 1981, 898: 361-381.

[14]昌金虎,陸君安,陳士華.混沌時間序列分析及其應用[M].武漢:武漢大學出版社,2002.

Zhang Z, Zha H. Principal manifolds and nonlinear di-mensionality reduction via tangent space alignment[J].SIAM Joumal on Scientific Computing, 2004, 26 (1):313-338.

[16]楊慶,陳桂明,江良洲,等.帶標志點的LTSA算法及其在軸承故障診斷中的應用[J].振動工程學報,2012,25(6):732-738.

Yang Qing, Chen Guiming, Jiang Liangzhou, et al. Lo-cal tangent space alignment algorithm based on selectinglandmark points and its application in fault diagnosis forrolling bearing [Jl. Joumal of Vibration Engineering,2012,25(6):732-738.

[17]Zhu H, Wang X, Zhao Y, et al. Sparse representationbased on adaptive multiscale features for robust machin-ery fault diagnosis[Jl. Proceedings of the Institution ofMechanical Engineers, Part C: Journal of MechanicalEngineering Science, 2015, 229( 12): 23032313.

主站蜘蛛池模板: 中美日韩在线网免费毛片视频| vvvv98国产成人综合青青| 午夜影院a级片| 亚洲综合网在线观看| 国产成人免费视频精品一区二区 | 美女国内精品自产拍在线播放 | 97av视频在线观看| 亚洲无码高清一区二区| 小13箩利洗澡无码视频免费网站| 丰满的少妇人妻无码区| 欧美不卡二区| 国产成人免费手机在线观看视频| 日韩欧美国产成人| 色哟哟国产成人精品| 97久久人人超碰国产精品| 午夜欧美在线| 97成人在线观看| 青草视频在线观看国产| 国产裸舞福利在线视频合集| 国产乱人伦AV在线A| 波多野结衣一二三| 久久久久国产精品熟女影院| 日日碰狠狠添天天爽| 亚洲日本中文字幕乱码中文| 国产96在线 | 国产SUV精品一区二区6| 国产高清不卡视频| 亚洲一区免费看| 99精品一区二区免费视频| 久草视频中文| 日韩av资源在线| 波多野结衣无码AV在线| 亚洲国产系列| 免费在线成人网| 国产精品2| 国产成人无码久久久久毛片| 制服无码网站| 无码福利视频| 国产女人18水真多毛片18精品| 婷婷色狠狠干| 99在线观看精品视频| 亚洲成年人片| 热九九精品| 国产精品网拍在线| 精品自窥自偷在线看| 成人免费一级片| 丁香六月综合网| 国产69囗曝护士吞精在线视频| 亚洲人成网站在线观看播放不卡| 日韩亚洲综合在线| 国产麻豆精品在线观看| 免费国产黄线在线观看| 青青草原偷拍视频| 欧美亚洲综合免费精品高清在线观看 | 色婷婷天天综合在线| 一级看片免费视频| 亚洲日韩久久综合中文字幕| 成年看免费观看视频拍拍| 亚洲综合香蕉| 国产情精品嫩草影院88av| 亚洲色图欧美视频| 国产91蝌蚪窝| 国产成人h在线观看网站站| 国产麻豆va精品视频| 久草热视频在线| 国产乱肥老妇精品视频| 欧美日韩va| 天天综合网色| 亚洲国产在一区二区三区| 亚洲国产清纯| 日本在线视频免费| 亚洲国产在一区二区三区| 日本一区二区三区精品视频| 97人人模人人爽人人喊小说| 国产精品30p| 亚洲欧美日韩动漫| 亚洲第七页| 亚洲精品国产精品乱码不卞| 熟妇丰满人妻| 国产精品成| 人人艹人人爽| 欧美在线综合视频|