亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        一種基于改進(jìn)EMD分解人車地震動信號識別算法

        2022-05-10 06:29:12鄒瑛珂賈云飛劉素蕓
        中國測試 2022年4期
        關(guān)鍵詞:模態(tài)特征信號

        鄒瑛珂, 賈云飛, 劉素蕓

        (1. 南京理工大學(xué)機(jī)械工程學(xué)院,江蘇 南京 210094; 2. 山西北方興安化學(xué)工業(yè)有限公司,山西 太原 030008)

        0 引 言

        近年來,隨著安全形勢的變化以及傳感器技術(shù)的不斷成熟,自動化預(yù)警、“無人值守”概念在安保領(lǐng)域和軍事領(lǐng)域越來越受重視。作為周界安防系統(tǒng)的一環(huán),對地面目標(biāo)的識別需求最為廣泛;因此,如何從傳感器所采集到的含噪音地面地震動信號中提取需要的信息來實(shí)現(xiàn)對人車目標(biāo)的識別就顯得尤為重要。其中,研究在復(fù)雜振動信號中提取有用的特征值的方法對解決地震動信號人車識別問題具有重要意義。針對此問題所開展的研究相對較少,目前比較主流的傳統(tǒng)方法是以短時過零點(diǎn)[1]、峭度[2]、小波變換、經(jīng)驗(yàn)?zāi)B(tài)分解為代表的方法;在國內(nèi),針對人車信號的專門研究比較少,類似問題的前沿方法包括由陶良小等提出的基于線性陣列包絡(luò)線偏移疊加的檢測方法[3],吳甄非等提出的HHT變換[4]來幫助識別人車信號。在國外,針對類似震動信號,Muralidharan等也提出了一種基于隨機(jī)森林和變分模態(tài)分解的震動信號特征提取方法[5]。劉自然提出一種基于改進(jìn)經(jīng)驗(yàn)小波變換的沖擊信號特征提取方法[6]。這些方法在平坦且震動傳導(dǎo)良好的水泥路面、實(shí)驗(yàn)室環(huán)境中表現(xiàn)出色。但在野外環(huán)境中,由于土質(zhì)、天氣、地下生物活動等因素影響,地震動信號常常包含著許多噪音,且往往需要探測較遠(yuǎn)距離的目標(biāo),這些因素造成整個信號信噪比偏低,且后兩者計算量偏大,不易進(jìn)行實(shí)際應(yīng)用。因此上述方法雖然在實(shí)驗(yàn)環(huán)境表現(xiàn)出色,但在野外松軟土質(zhì)環(huán)境中識別率較低?;诃h(huán)境和被測信號的特征兩個因素,本文使用希爾伯特變換提取包絡(luò)線來實(shí)現(xiàn)一次降噪,并將一次降噪后的信號利用改進(jìn)型的EMD分解方法來獲得具有高信噪比的IMF分量,然后通過互相關(guān)函數(shù)、統(tǒng)計學(xué)量等方法獲得特征值,最后投入隨機(jī)森林分類器中進(jìn)行分類從而實(shí)現(xiàn)較遠(yuǎn)距離人車信號的識別。

        1 算法原理

        1.1 希爾伯特變換

        希爾伯特變換是一種常用的提取信號包絡(luò)線的方法。該變換方法可以看做是一個正交濾波器,可將所有的正頻率分量移相–90°,對負(fù)頻率分量移相90°,從而能將一個實(shí)信號變換為一個復(fù)信號的虛部。通過求解該復(fù)信號的幅值即可求得原信號的包絡(luò)從而將其低頻分量解調(diào)出來[7]。

        一個實(shí)信號x(t)的希爾伯特變換定義為:

        將原信號作為實(shí)部,式(1)信號作為虛部,可得到一個解析復(fù)信號:

        得到該復(fù)信號后,通過求其幅值則可得到其復(fù)包絡(luò)信號:

        1.2 經(jīng)驗(yàn)?zāi)B(tài)分解

        經(jīng)驗(yàn)?zāi)B(tài)分解方法(EMD)是依據(jù)數(shù)據(jù)自身的時間尺度特征來進(jìn)行信號分解而無須預(yù)先設(shè)定任何基函數(shù)。正是這樣的特性,使EMD 方法在理論上可以應(yīng)用于任何類型的信號的分解, 因而在處理非平穩(wěn)及非線性數(shù)據(jù)上,具有非常明顯的優(yōu)勢和很高的信噪比。

        EMD分解步驟可歸類為如下幾步[8]:

        1)找出原始信號所有極值點(diǎn)。

        4)原始信號減均值包絡(luò)線,得到中間信號:

        5)判斷該中間信號是否滿足IMF的兩個條件:

        ①在整個數(shù)據(jù)段內(nèi),極值點(diǎn)的個數(shù)和過零點(diǎn)的個數(shù)必須相等或相差最多不能超過一個。

        ②在任意時刻,由局部極大值點(diǎn)形成的上包絡(luò)線和由局部極小值點(diǎn)形成的下包絡(luò)線的平均值為零,即上、下包絡(luò)線相對于時間軸局部對稱。

        如果滿足,該信號就是一個IMF分量IMFi(n),用原始信號減去該IMF,作為新的原始信號,返回第一步。

        重復(fù)步驟1)~5)k次,直到滿足IMF篩選停止的原則。IMF篩選停止原則為:

        EMD分解的優(yōu)點(diǎn)是克服了基函數(shù)無自適應(yīng)性的問題,對于一段信號可以不進(jìn)行事先的復(fù)雜分析或參數(shù)調(diào)整就可直接分解出所需的特征信號,且運(yùn)算量相對小波分析等算法來說更較小,這是其可運(yùn)用在下位機(jī)上的關(guān)鍵優(yōu)勢。但其算法本身的原理會導(dǎo)致其分解出來的IMF出現(xiàn)模態(tài)混疊現(xiàn)象甚至造成信號失真;并且在低信噪比的信號中分解出來的IMF信噪比也會較低,導(dǎo)致這種現(xiàn)象會被進(jìn)一步放大,從而影響后續(xù)識別的結(jié)果。雖然有人也提出了包括EEMD[9]等改進(jìn)方法,但針對此類問題的效果提升不大。因此需要對EMD分解進(jìn)行改進(jìn)。

        1.3 協(xié)方差疊加經(jīng)驗(yàn)?zāi)B(tài)分解

        針對低信噪比信號的EMD分解方法,本文提出了一種名為協(xié)方差疊加經(jīng)驗(yàn)?zāi)B(tài)分解(covariance superposition empirical mode decomposition, CSEMD)的改進(jìn)型EMD方法。其基本原理是對信號中具有高相關(guān)性的信號進(jìn)行加重,低相關(guān)性的信號(一般為噪音)進(jìn)行減重。在此基礎(chǔ)上進(jìn)行EMD分解,從而得到擁有較高信噪比的IMF分量。CSEMD分解可歸類為如下幾步:

        1)對信號x(n)進(jìn)行EMD分解。

        2)將EMD分解所得出的IMF信號IMFi(n)與原始信號x(n)進(jìn)行協(xié)方差運(yùn)算。所得到的協(xié)方差矩陣為:

        3)將所得協(xié)方差標(biāo)準(zhǔn)化后進(jìn)行排序,取協(xié)方差最大的IMF信號為IMFmax(n),取協(xié)方差最小的IMF信號為IMFmin(n)。二者相減并取其平均值,將值加至原始信號得到高信噪比信號x'(n)。其公式為:

        4)對x'(n)再次進(jìn)行EMD分解,可得高信噪比IMF。

        1.4 隨機(jī)森林算法

        隨機(jī)森林算法是機(jī)器學(xué)習(xí)、計算機(jī)視覺等領(lǐng)域內(nèi)應(yīng)用較為廣泛的一個分類器,也可以作為一種數(shù)據(jù)降維的手段[10]。通過集成學(xué)習(xí)的思想,將原本作為弱分類器的通過不同訓(xùn)練集訓(xùn)練出來多個決策樹算法集成到一起從而形成一個強(qiáng)分類器,實(shí)現(xiàn)復(fù)雜的分類。

        隨機(jī)森林在變量(列)的使用和數(shù)據(jù)(行)的使用上進(jìn)行隨機(jī)化,生成很多分類樹,每個樹都是一個獨(dú)立的判斷分支,互相之間彼此獨(dú)立。其中生成分類樹的算法本文中采用基尼指數(shù)[11]:

        當(dāng)在基于某屬性對一個新的對象進(jìn)行分類判別時,隨機(jī)森林中的每一棵樹都會給出自己的分類選擇,并由此進(jìn)行加權(quán)輸出。本算法相比其他機(jī)器算法優(yōu)點(diǎn)是不易過擬合,不用進(jìn)行相關(guān)參數(shù)的調(diào)整,其算法理解起來較容易,可以很好地移植到下位機(jī)上進(jìn)行獨(dú)立運(yùn)算,且判斷過程的運(yùn)算量較低,不用運(yùn)算復(fù)雜的矩陣。

        2 信號仿真實(shí)驗(yàn)

        下文將通過人員和車輛兩種仿真信號各自使用希爾伯特變換+CSEMD與希爾伯特變換+EMD算法進(jìn)行分解,并將所得的特征信號在信噪比和模態(tài)混疊現(xiàn)象方面進(jìn)行對比。

        2.1 人員仿真信號生成

        根據(jù)特點(diǎn)分析,可以將人員信號看作是一種類沖擊信號,其特征是不連續(xù)出現(xiàn)但是信號本身頻率較高。由于在野外環(huán)境下人員信號會受到土壤中某些高頻振動信號(如風(fēng)吹樹導(dǎo)致樹根震動、地下動物的活動)的干擾,因此在信號中會有能量較高的高頻震動信號和幅值較低的隨機(jī)信號。根據(jù)分析生成如下信號,長度為1024個數(shù)據(jù)點(diǎn),采樣頻率為1000 Hz。其中人員仿真信號方程為:

        式中:N——數(shù)據(jù)點(diǎn)位置;

        f——需要提取的特征信號;

        rand——絕對值介于0到0.25的高斯白噪音。

        其中,260

        圖1 人員仿真信號

        2.2 人員仿真信號實(shí)驗(yàn)

        利用希爾伯特變換提取的包絡(luò)線如圖2所示。

        圖2 人員仿真信號包絡(luò)線

        可以看出包絡(luò)線雖然提取出了沖擊信號,但由于毛刺過多會嚴(yán)重影響后續(xù)算法的準(zhǔn)確性。直接使用CSEMD提取人員腳步特征結(jié)果如圖3所示。

        圖3 包絡(luò)線CSEMD分解結(jié)果IMF3

        通過IMF3模態(tài)可以看出雖然CSEMD可以提取出腳步信號,但其信噪比較低,且有比較嚴(yán)重的模態(tài)混疊現(xiàn)象,不利于后續(xù)包括均值、方差、峭度等特征量的提取。因此必須對包絡(luò)線進(jìn)行平滑處理。由于人行走信號是一種類沖擊信號,幅值較大。因此進(jìn)行的平滑處理不可以過濾該類信號。同時觀察到包絡(luò)線中的毛刺主要是圍繞某個值所產(chǎn)生的上下小范圍跳變,是高頻低幅值噪音,使用滑動平均濾波效果好??紤]到二者性質(zhì),使用滑動平均濾波可使信號平滑的同時保留沖擊信號特征[12]。其濾波公式為:

        其中n為幀長。此處步長取1。

        通過滑動平均濾波平滑處理后的包絡(luò)線信號如圖4所示。

        圖4 人員仿真信號包絡(luò)線(平滑后)

        可以看出在保留了沖擊信號的同時過濾掉了絕大多數(shù)的毛刺信號,便于后續(xù)的算法分解。使用EMD和CSEMD提取人員腳步特征結(jié)果如圖5和圖6所示。

        圖5 平滑后包絡(luò)線EMD分解結(jié)果

        圖6 平滑后包絡(luò)線CSEMD分解結(jié)果IMF3

        對比圖3,可以看出對平滑后的包絡(luò)線進(jìn)行分解后獲得的腳步信號信噪比更高,更有利于對包括翹度、方差等特征量提取。但其中EMD分解提取出的第一個腳步信號有所失真。只有平滑后的包絡(luò)線通過CSEMD分解提取出的腳步信號在保證信號不失真的同時擁有更高信噪比,模態(tài)混疊現(xiàn)象影響較小,更有利于從中提取特征量從而完成識別。

        2.3 車輛仿真信號生成

        根據(jù)特點(diǎn)分析,車輛信號可以看做是一種連續(xù)低頻信號。且由于技術(shù)指標(biāo)要求車輛探測距離較遠(yuǎn),因此幅值較低。在野外環(huán)境下會受到土壤中某些高頻振動信號(如風(fēng)吹樹導(dǎo)致樹根震動、地下動物的活動)的干擾,因此在信號中會有能量較高的高頻震動信號和幅值較低的隨機(jī)信號。生成信號如圖7所示,長度為 900個數(shù)據(jù)點(diǎn),采樣頻率為 1000 Hz。由于車輛震動信號也與人行走震動信號一樣是自身產(chǎn)生的震動信號與地面中的高頻信號的疊加。因此車輛仿真信號方程為:

        圖7 車輛仿真信號

        2.4 車輛仿真信號實(shí)驗(yàn)

        參照人員信號,對其進(jìn)行希爾伯特變換,所得包絡(luò)線如圖8所示。

        圖8 車輛仿真信號包絡(luò)線

        直接使用CSEMD從中提取車輛腳步特征結(jié)果如圖9和圖10所示。

        圖9 車輛仿真信號包絡(luò)線CSEMD分解結(jié)果IMF3

        圖10 包絡(luò)線CSEMD分解結(jié)果IMF3瞬時頻率

        可以觀察到由于大量毛刺的影響,使用CSEMD進(jìn)行信號分解有比較嚴(yán)重的模態(tài)混疊。由于車輛信號屬于一種低頻信號,進(jìn)行滑動平均濾波處理對車輛信號的影響較小。因此同樣使用滑動平均濾波,所得圖像如圖11所示。

        圖11 車輛仿真信號包絡(luò)線(平滑后)

        使用EMD和CSEMD對平滑后的包絡(luò)線提取車輛特征結(jié)果如圖12和圖13所示。

        圖12 車輛仿真信號平滑后包絡(luò)線EMD分解結(jié)果IMF3

        圖13 車輛仿真信號平滑后包絡(luò)線CSEMD分解結(jié)果IMF3

        由于在時域上二者的對比不是很明顯,因此觀察二者的瞬時頻率結(jié)果如圖14和圖15所示。

        圖14 平滑后包絡(luò)線EMD分解結(jié)果IMF3瞬時頻率

        圖15 平滑后包絡(luò)線CSEMD分解結(jié)果IMF3瞬時頻率

        從瞬時頻率圖中與圖10相比可以看出,通過CSEMD提取出的車輛信號的IMF3模態(tài)可改善模態(tài)混疊。雖然對連續(xù)車輛信號的提取改善相較于EMD沒有腳步信號明顯,但對其進(jìn)行特征提取仍然是有利的。

        綜上所述,本文提出的希爾伯特+CSEMD算法對人員腳步信號的提取相對于傳統(tǒng)希爾伯特+EMD算法在提高信噪比、防止提取信號失真和抑制模態(tài)混疊方面有明顯正面促進(jìn)作用。在針對連續(xù)車輛信號方面也有一定改善。因此本算法對于提取人車信號是有效的。

        3 實(shí)際信號實(shí)驗(yàn)與結(jié)果分析

        3.1 人員與車輛信號樣本說明

        在晴天、低速風(fēng)(風(fēng)速不高于3級)、均質(zhì)土壤的環(huán)境中通過采集卡和VAS-100地震動傳感器獲得人車原始信號,采集頻率為 1000 Hz。地震動傳感器插入地表以下約20 cm處,測試車輛為某品牌SUV在距離傳感器100 m以20 km/h勻速向著傳感器方向行駛,測試人員為單人從距離傳感器50 m處向著傳感器方向正常勻速行走。車輛和人員的部分原始信號如圖16和圖17所示。

        圖16 車輛原始信號

        圖17 人員原始信號

        由于信噪比較低,進(jìn)行傅里葉變換后其簡單的時域特征量(包括平均數(shù),中位數(shù),方差,最大最小值等)和頻譜功率譜的各個特征量(包括平均數(shù),中位數(shù),方差,最大最小值等)都較為相近,因此無法單純從頻域上進(jìn)行識別。而運(yùn)用一些較為復(fù)雜的頻域特征提取算法如小波變換等計算量又過大,無法在下位機(jī)中實(shí)現(xiàn)。因此需要在時域或時頻域上尋找其他計算量相對較小又能夠提取其特征的方式來實(shí)現(xiàn)目標(biāo)功能。

        3.2 實(shí)際信號實(shí)驗(yàn)

        處理數(shù)據(jù)為上述實(shí)地采集的車輛人員地震動數(shù)據(jù),其中車輛數(shù)據(jù)為距離傳感器100 m處開始向著傳感器方向接近的過程,人員數(shù)據(jù)為距離傳感器50 m處向著傳感器方向接近的過程,采樣頻率F=1000 Hz,2048個點(diǎn)為一幀數(shù)據(jù)。使用本文提出的希爾伯特變換+CSEMD算法和希爾伯特變換+EMD算法進(jìn)行處理并比較結(jié)果。

        首先通過希爾伯特變換提取二者的包絡(luò)線相當(dāng)于對原信號進(jìn)行了解調(diào)。去掉包絡(luò)線偏置,對其進(jìn)行EMD分解,可得到EMD分量如圖18和圖19所示。

        圖18 人員信號EMD分解

        圖19 車輛信號EMD分解IMF1-5

        再將信號的IMF與原始信號進(jìn)行協(xié)方差運(yùn)算,篩選出協(xié)方差最大和最小的IMF進(jìn)行相減取平均,并加在原始信號中,再對其進(jìn)行EMD分解,得到的CSEMD分量如圖20和圖21所示。

        圖20 人員信號 CSEMD分解

        圖21 CSEMD車輛信號IMF1-5

        可以看出,EMD分解人車信號所得的IMF2與CSEMD的IMF1等效,且都包含了顯著的特征信號:人員信號對應(yīng)的IMF有效提取出了人員腳步行走較離散的特征信號,車輛信號對應(yīng)的IMF提取出了車輛運(yùn)動較為連續(xù)的特征信號。同時,車輛信號的幅值都提高了50%以上。

        但其模態(tài)混疊段的幅值均出現(xiàn)了下降。其中人員信號的模態(tài)混疊幅值降低了。其中人員信號模態(tài)混疊段幅值降低了80%,車輛信號模態(tài)混疊段幅值降低了50%。EMD的IMF2與CSEMD的IMF1如圖22~圖25所示。

        圖22 EMD人員信號IMF2模態(tài)混疊段

        圖23 CSEMD人員信號IMF1模態(tài)混疊段

        圖24 EMD車輛信號IMF2模態(tài)混疊段

        圖25 CSEMD車輛信號IMF1模態(tài)混疊段

        最后,可采用CSEMD分解得到的IMF1與IMF2進(jìn)行統(tǒng)計學(xué)方面的分析,從而得出特征指標(biāo)。筆者在本文中使用方差閾值、平均值、中位數(shù)、峭度、自動編碼器算法來提取特征。將2S信號分幀為20段,計算每段的特征值。方差超過閾值的點(diǎn)則視為一次有效腳步信號。將信號的腳步信號量作為一個特征量。同時,使用互相關(guān)函數(shù)對分幀信號進(jìn)行特征提取。若相關(guān)系數(shù)大于閾值,則將該幀視為腳步信號,同樣也可作為一個特征量。經(jīng)過大量數(shù)據(jù)比對和統(tǒng)計,車輛由于是連續(xù)信號且信號頻率不穩(wěn)定,其腳步特征量一般小于人行走信號。而人行走無論快慢,每秒步頻大約為3~4步,因此腳步信號量值較高且穩(wěn)定。

        3.3 實(shí)驗(yàn)結(jié)果與分析

        總結(jié)前面的步驟,整個算法流程圖如圖26所示。

        圖26 算法流程圖

        將本文算法所得的兩個特征量、過零點(diǎn)與峭度、直接EMD分解所得特征量投入隨機(jī)森林算法中進(jìn)行訓(xùn)練并測試,車輛訓(xùn)練樣本數(shù)量為100組,人員訓(xùn)練樣本數(shù)量為100組,車輛測試樣本數(shù)量為100組,人員測試樣本數(shù)量為100組,測試樣本和訓(xùn)練樣本相互獨(dú)立。其中樣本均為在陰天與小雨,風(fēng)力不高于6級,在擁有均勻硬土質(zhì)地的公園所采集。其中車輛信號為2噸家用小轎車以40 km/h的速度從距離傳感器150~10 m處勻速駛近時采集的各距離樣本,人員信號為體重為80 kg的單人以正常步行速度從距離50~5 m處勻速走進(jìn)時采集的各距離樣本。

        所得結(jié)果如表1所示。

        表1 本文算法與傳統(tǒng)算法識別結(jié)果比較

        由表1可知,在使用同一種分類算法的情況下,最簡單最傳統(tǒng)的零點(diǎn)+峭度算法從結(jié)果來看完全無法在野外環(huán)境中識別人員和車輛。兩個特征量在遠(yuǎn)場和近場環(huán)境中其結(jié)果完全不同甚至所得結(jié)果相反,造成這種現(xiàn)象的原因是遠(yuǎn)場信噪比較近場更低,因此在遠(yuǎn)場環(huán)境中所得到的特征量更偏向于是噪音本身的特征量。

        而使用EMD分解所得特征量識別率偏低,且識別成功的樣本多為近場樣本,實(shí)用性偏低。

        直接對遠(yuǎn)場信號進(jìn)行EMD分解會造成比較嚴(yán)重的模態(tài)混疊,從而影響后續(xù)的特征量提取,最后造成信號識別結(jié)果錯誤。同理,使用希爾伯特變換提取包絡(luò)線后再對包絡(luò)線進(jìn)行EMD分解雖然可以改善模態(tài)混疊的現(xiàn)象,但其仍然對識別率有較大影響。使用本文算法可以有效改善模態(tài)混疊現(xiàn)象,使識別成功率得到較大提升,在遠(yuǎn)場和近場環(huán)境下均有良好表現(xiàn)。且相對其他方法計算量并沒有太大提升。

        4 結(jié)束語

        本文針對在野外環(huán)境中對人車地震動信號進(jìn)行正確識別這一問題,提出了通過希爾伯特變換提取包絡(luò)線來實(shí)現(xiàn)降噪,并將降噪后的信號利用改進(jìn)型EMD分解方法——CSEMD分解來獲得具有特征量的IMF分量,再利用一些統(tǒng)計學(xué)方法提取特征量并將其送入隨機(jī)森林分類算法中進(jìn)行人車識別。實(shí)驗(yàn)結(jié)果表明:本文使用的識別算法相比于傳統(tǒng)算法結(jié)果更好,更可用于實(shí)際,且具有較好的抗噪性能。下一步工作準(zhǔn)備在此基礎(chǔ)上嘗試使用其他機(jī)器學(xué)習(xí)算法,進(jìn)一步優(yōu)化CSEMD分解算法效率并改進(jìn)包絡(luò)線提取算法,從而降低運(yùn)算量同時提升識別正確率;優(yōu)化算法準(zhǔn)確度,使其可以在惡劣天氣(如下雨大風(fēng)等情況)進(jìn)行檢測;同時將該算法應(yīng)用于諸如無人機(jī)飛行噪聲等其他目標(biāo)探測中,提高其泛用性。

        猜你喜歡
        模態(tài)特征信號
        信號
        鴨綠江(2021年35期)2021-04-19 12:24:18
        完形填空二則
        如何表達(dá)“特征”
        不忠誠的四個特征
        基于FPGA的多功能信號發(fā)生器的設(shè)計
        電子制作(2018年11期)2018-08-04 03:25:42
        抓住特征巧觀察
        基于LabVIEW的力加載信號采集與PID控制
        國內(nèi)多模態(tài)教學(xué)研究回顧與展望
        基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識別
        由單個模態(tài)構(gòu)造對稱簡支梁的抗彎剛度
        計算物理(2014年2期)2014-03-11 17:01:39
        夹得好湿真拔不出来了动态图| 激情五月婷婷六月俺也去| 最新日本女优中文字幕视频| 影音先锋久久久久av综合网成人| 国产肥熟女视频一区二区三区 | 日本一区二区在线播放| 亚洲中文字幕无码永久在线| 欧美精品高清在线xxxx| 亚洲一区二区精品在线| 欧美老妇多毛xxxxx极瑞视频| 国产在线精品一区二区不卡| 久久亚洲精彩无码天堂 | aa视频在线观看播放免费| 亚洲一区二区刺激的视频| 久久精品国产亚洲av无码娇色| 欧美成人三级一区二区在线观看| 看全色黄大黄大色免费久久| 开心五月婷婷激情综合网| 在线视频观看免费视频18| 久久久久中文字幕无码少妇| av在线不卡一区二区三区| 18禁在线永久免费观看| 免费人成视频在线观看网站| 久久久久久久国产精品电影| 亚洲一区二区综合精品| 国产成人精品999视频| 亚洲最大天堂无码精品区| 亚洲一区二区女优av| 亚洲av熟女少妇久久| 国产精品三级在线观看无码| 久久露脸国产精品WWW| 精品精品国产三级av在线| 玩弄丰满奶水的女邻居 | 美腿丝袜日韩在线观看| 国产乱人激情h在线观看| 国产欧美一区二区成人影院| 日韩精品久久伊人中文字幕| 丰满人妻被两个按摩师| 国产精品白浆一区二小说| 一区二区亚洲 av免费| 久久久久久久久无码精品亚洲日韩|