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

分譜處理算法在基于超聲導波損傷識別中的應用

2012-02-12 11:40:12苗曉婷李富才
振動與沖擊 2012年1期
關鍵詞:信號檢測

苗曉婷,李富才,孟 光

(上海交通大學 機械系統與振動國家重點實驗室,上海 200240)

近20年內,超聲導波(如蘭姆波)已被廣泛應用于無損健康檢測領域[1-4]。基于超聲導波的損傷識別方法通常是把試件無缺陷時所采集的波信號作為基準,并把當前所采集的檢測波信號與基準波信號進行對比,從而提取出由損傷散射的波信號的特征信息,如飛行時間(ToF)等[5-6],達到損傷識別和定位等目的。由于超聲導波不可避免的頻散特性以及工作環境中的噪聲嚴重降低了檢測波信號的信噪比(SNR),很難依據原始的檢測波信號評估出損傷散射的波信號的ToF。因此信號處理技術對基于超聲導波的損傷識別起到至關重要的作用,決定了ToF的精確評估以及三角定位算法的有效性。至今已有多種信號處理技術被廣泛應用于提高信號的 SNR,比如 Hilbert-Huang變換(HHT)[7-8]和小波變換(WT)[9-10]。通過適當的降噪處理,可以直接把檢測波信號與基準波信號相減以得到損傷散射的波信號,并基于能量分析方法評估出損傷散射的波信號的ToF[11-12]。然而在實際的工作環境中,檢測波信號可能被不確定的噪聲干擾,如果噪聲能量分布在目標波信號的頻率帶內,此時通過HHT和WT仍然無法有效地排除噪聲,噪聲能量將淹沒損傷散射的波信號的能量,導致這種基于能量分析方法的失效。

由于分譜處理(SSP)算法具有較好的抵抗隨機噪聲干擾的能力,它已被廣泛應用于超聲波的無損探傷和雷達定位[13-16]。當目標信號在所采集的信號中出現時,通過SSP算法所重建的一系列時間域信號將呈現出一致的幅值信息[17]。利用適當的特征提取算法,在各個采樣時間點上提取這一系列重建的時間域信號的幅值并衡量它們的一致性,最終可以推斷那些呈現出較好的幅值一致性的采樣時間點對應于目標信號的出現。由于所衡量的幅值一致性不受噪聲的影響,即使所采集的信號的SNR較低,SSP算法仍然可以正確地評估目標信號的ToF。

本研究分別在無噪聲和有噪聲干擾檢測波信號的情況下利用SSP算法評估鋁板中的切縫缺陷所散射的基礎階對稱(S0)模式的ToF。結合傳感網絡中各條激勵-傳感路徑所評估的ToF,最終利用三角定位算法在傳感網絡所包圍的檢測區域內識別并定位切縫缺陷。

1 基本原理

1.1 基于ToF的三角定位算法

在被用于損傷識別的波信號的特征信息中,ToF是比較直觀的特征信息之一。利用三角定位算法,結合傳感網絡中各條傳感路徑所評估的損傷散射波信號的ToF,可以在傳感網絡所包圍的檢測區域內識別并定位缺陷[18-19]。薄板結構中的超聲導波——蘭姆波包括兩種模式,即對稱模式與反對稱模式,并且不同的模式各自以不同的速度在板內傳播。對板結構而言,在任何激勵頻率下至少會有兩種模式的蘭姆波被同時激勵出來并在結構中傳播。通過調整激勵的中心頻率,可以使結構中傳播的模式降低到最少,即基礎階的對稱(S0)和反對稱模式(A0)[20]。在這種情況下,由于 S0模式具有最快的傳播速度,本研究通過評估損傷散射的S0模式的ToF來定位缺陷以避免波信號中其它波模式的干擾。由損傷散射引起的 S0模式的 ToF可以表示為:

其中:

TA-D-S是激發的S0模式從激勵點傳播到缺陷,再從缺陷傳播到感應點所需要的飛行時間。TA-S是激發的S0模式從激勵點直接傳播到感應點所需要的飛行時間。LA-D是激勵點(xA,yA)與缺陷的中心位置(xD,yD)之間的距離。LD-S是缺陷的中心位置與感應點(xS,yS)之間的距離。LA-S是激勵點與感應點之間的距離。vg是S0模式的群速度。方程式(1)的解繪制出一個以激勵點和感應點為焦點的橢圓,這個橢圓便是缺陷中心位置的坐標軌跡。結合傳感網絡中各條傳感路徑所評估的ToF,可以繪制出多條坐標軌跡,它們的交匯區域便是缺陷可能出現的區域。

1.2 分譜處理(SSP)算法

一般而言,一個原始信號r(t)包括目標信號s(t)和噪聲n(t),可以表示為:

其中,T是信號的采樣時間。信噪比(SNR)被用于衡量信號被噪聲污染的程度,定義為[21]:

目標信號被噪聲污染得越嚴重,SNR的值就越小。SSP算法包括兩步:① 分譜;② 特征提取,算法流程如圖1所示。

圖1 SSP算法的流程圖Fig.1 Schematic flow chart of SSP algorithm

1.2.1 分譜

(1)快速傅里葉變換(FFT)

利用FFT,可以獲得原始信號r(t)的頻率 -能量譜:

在本研究中,通過FFT所獲得的能量譜的能量最大值(能量峰值)的一半所對應的起始頻率和截止頻率分別被用于指定信號的主要頻率帶,如圖2所示。

圖2 用于切分主要頻率帶的一組高斯帶通濾波器的頻率譜分布Fig.2 Frequency spectrum distribution of a series of Gaussian band-pass filters cutting the main frequency band

(2)帶通濾波器

利用一組帶通濾波器把所指定的主要頻率帶切分為一系列窄的頻率帶(Xi(ω),i=1~M),如圖2所示。SSP算法的有效性與以下參數的合理選取有關:① 帶通濾波器函數;② 帶通濾波器的總數;③ 這一組帶通濾波器的頻率譜分布形式。本研究采用了總數為M的一組高斯帶通濾波器,相鄰兩個帶通濾波器的中心頻率之間的頻率間隔為Δf。它們各自具有不同的中心頻率(fi,i=1~M),相同的帶寬(2× Δf)。M個高斯帶通濾波器把主要頻率帶切分為M個相互重疊的頻率帶以盡可能地減小信號成分的遺漏[22]。M和 Δf分別定義為:

其中T是信號的采樣時間,即信號總的時間長度;B為主要頻率帶的帶寬。

(3)傅里葉逆變換(IFFT)

所切分的M個窄的頻率帶(Xi(ω),i=1~M)各自進行傅里葉逆變換,重建了M個時間域信號:

本研究把xi(t),(i=1~M)的瞬時幅值的變化范圍定義為信號r(t)的瞬時幅值變化度(IAVD(t)),表示為:

其中

如果xi(t),(i=1~M)的瞬時幅值在較小的范圍內變化,則對應較小的IAVD(t)值;相反,則對應較大的IAVD(t)值。

1.2.2 特征提取

基于SSP算法的理論,當目標信號在所采集的信號中出現時,xi(t),(i=1~M),呈現出比較一致的瞬時幅值,對應較小的IAVD(t)值[14]。在以前的研究中,用于衡量瞬時相位一致性的特征提取方法(如極性閾值法和幅值最小值法)對SSP算法中所選用的帶通濾波器的帶寬以及頻率譜分布形式比較敏感,這種敏感性在一定程度上降低了SSP算法的有效性[17]。本文所提出的特征提取方法是通過比較基準波信號的IAVD(IAVDbenchmark(t))和檢測波信號的IAVD(IAVDdanage(t))來評估損傷散射的S0模式的ToF。當激勵的S0模式或由損傷散射的S0模式在波信號中出現時,波信號的IAVD值較小。如果在某些采樣時間點,IAVDbenchmark(t)的值較大,而IAVDdanage(t)的值較小,則意味著此時S0模式沒有出現在基準波信號中卻出現在檢測波信號中,于是可以推斷損傷散射的S0模式出現在這些采樣時間點,并由此評估出損傷散射的S0模式的ToF。特征提取方法可具體描述如下:

當激發的S0模式到達時(TA-S),IAVDbenchmark(t)與IAVDdanage(t)之間的差異(Δ(TA-S))被設定為閾值。在TA-S以后,如果 Δ(t)的值比 Δ(TA-S)大,則設定 SSP 算法在t時的輸出(y(t))等于檢測波信號rdanage(t)的幅值,否則設定y(t)的值為0。第一個值不為0的y(t)所對應的采樣時間點被評估為S0模式從激勵點傳播到缺陷,再從缺陷傳播到感應點所需要的飛行時間TA-D-S。把所評估的TA-D-S值代入方程式(1)中,則可以評估出缺陷散射的S0模式的飛行時間TD,并繪制出表示缺陷中心位置的坐標軌跡。

2 實驗評估

2.1 實驗裝置

實驗中利用一個基于VXI平臺的信號發生和采集系統,以20.48 MHz的采樣率采集波信號。這個系統包括:信號發生器(Agilent E1441)、信號放大器(Piezo-Sys EPA-104)、信號解調器(Agilent E3242A)和信號數字離散器(Agilent E1437A)。漢寧窗調制的、中心頻率為300 kHz的、5周正弦調幅脈沖(峰峰值為60V)被用做激勵信號激發板中的蘭姆波。四個直徑為6.9 mm、厚度為 0.5 mm的圓形壓電應變片(PZT PI PIC151,PRYY-0929)被表面固定在鋁板(600 mm ×600 mm ×2.0 mm)上,組成一個包圍了400 mm ×400 mm大小的檢測區域的傳感網絡,如圖3(a)所示。利用一個坐標系來定位檢測區域,它的原點設在P3,X軸由P3指向P4,Y軸由P3指向P1。選用了的六條激勵-感應波信號的傳感路徑,各自采用“一發一收”的工作方式激勵和采集蘭姆波信號,覆蓋了檢測區域(圖3(b))。實驗中首先對無缺陷的鋁板進行檢測,并把所采集的波信號作為基本波信號。然后用一個厚度為0.6 mm的鋸片在先前的鋁板上切割出一條切縫缺陷,表格1詳細列出了切縫缺陷的位置和尺寸。對帶有切縫缺陷的鋁板進行檢測,并把所采集的波信號作為檢測波信號。各個信號的采集時間均設定為200 μs。為了模擬嘈雜的工作環境,幅值大小不同的白噪聲被分別加載到檢測波信號上,相應地生成了具有不同信噪比的檢測波信號。

2.2 實驗結果與分析

在無缺陷的鋁板中所測量到的S0模式在最短的傳播路徑(400 mm)上所需要的飛行時間TA-S為77 μs(第1 578個采樣點),由此推算S0模式的傳播群速度vg為5 195 m/s。下文舉例展示對傳感路徑P3-P4所采集的波信號進行信號處理的過程。所采集的波信號進行統一的幅值正則化。

圖3 (a)帶有切縫缺陷的鋁板與傳感網絡的示意圖以及(b)所選用的六條傳感路徑Fig.3(a)Schematic diagram of the notched aluminum plate with a sensor network and(b)the six selected sensing paths

表1 鋁板中切縫缺陷的位置和尺寸Tab.1 Notch in the aluminum plate

2.2.1 無噪聲干擾時的損傷識別

圖4 傳感路徑P3-P4所采集的基準波信號與無噪聲干擾時的檢測波信號Fig.4 The captured wave signals by sensing path P3 - P4 from benchmark and damage case without noise

圖5 基準波信號與無噪聲干擾時的檢測波信號的頻率-能量譜Fig.5 Frequency-energy spectrum of the captured wave signals from benchmark and damage case without noise

無噪聲干擾時,傳感路徑P3-P4所采集的基準波信號與檢測波信號的幅值和相位之間存在一定的差異(圖4),采樣起始處的信號為電磁干擾。利用FFT,得到它們的頻率-能量譜(圖5)。對基準波信號而言,當頻率為330 kHz時對應能量的最大值(0.13),能量峰值的一半所對應的起始頻率為295 kHz、截止頻率為365 kHz。因此所指定的主要頻率帶的帶寬為70 kHz(從295 kHz到365 kHz),即:

對檢測波信號進行分譜的過程中,指定與基準波信號相同的主要頻率帶(從295 kHz到365 kHz)。實驗中所設定的采樣時間為200 μs,即:

基于式(6)和式(7)可得:

根據以上所設定的參數,本研究采用一組總數為15、相鄰兩個帶通濾波器的中心頻率之間的頻率間隔為5 kHz、各自帶寬為10 kHz的高斯帶通濾波器,把基準波信號以及檢測波信號的主要頻率帶切分成15個相互重疊的頻率帶(Xi(ω),i=1~15)。然后各自進行IFFT,重建了15個時間域的信號(xi(t),i=1~15)。圖6舉例展示了由基準波信號所重建的時間域信號x6(t),x7(t)和x8(t)。根據方程式(9)分別求得IAVDbenchmark(t)與IAVDdanage(t),如圖7所示。根據方程式(10)分別求得Δ(t)和SSP算法的輸出y(t),如圖8所示;最終所評估的TA-D-S對應于第1 910個采樣時間點(圖8(b))。

圖6 由基準波信號所重建的時間域信號Fig.6 The reconstructed time- domain signals(a)x6(t),(b)x7(t)and(c)x8(t)of the captured wave signal from benchmark

圖7 所求得的基準波信號與檢測波信號(無噪聲干擾)的IAVDFig.7 IAVDbenchmark(t)and IAVDdamage(t)(without noise)

圖8 無噪聲干擾時的(a)Δ(t)和(b)SSP算法的輸出y(t)Fig.8(a)Δ(t)and(b)the output of SSP algorithm(y(t))without noise

結合所選用的六條傳感路徑各自評估的TA-D-S,利用1.1節中所陳述的三角定位算法繪制出的缺陷中心位置的坐標軌跡的交匯區域如圖9所示。把交匯區域中距離坐標軌跡的距離和最小的坐標識別為缺陷的中心位置坐標,所識別的中心位置以‘+’標記。表2中列出了所識別的中心位置坐標以及相對于實際中心位置的距離。

表2 在無噪聲干擾和有噪聲干擾(SNR=1 dB)時所識別的結果Tab.2 Identification results without and with noise(SNR=1 dB)

圖9 無噪聲干擾時所識別的缺陷位置Fig.9 Identified result without noise

2.2.2 有噪聲干擾時的損傷識別

當檢測波信號被幅值大小不同的白噪聲干擾時,相應地,檢測波信號的SNR分別為13 dB,8 dB,5 dB和1 dB,檢測波信號被噪聲污染得越嚴重,SNR的值就越小,如圖10所示。利用SSP算法對SNR=1 dB的檢測波信號進行處理。首先得到檢測波信號的頻率-能量譜(圖11),相比于無噪聲干擾時的檢測波信號的頻率-能量譜(圖5)可以看出噪聲能量嚴重干擾了檢測波信號的能量,很大程度地改變了頻率-能量的分布。對檢測波信號進行如2.2.1節中所陳述的相同的信號處理,可以求得SNR=1 dB的檢測波信號的IAVD,IAVDdamage(SNR=1dB)(圖12),與圖7比較可知在無噪聲干擾時和SNR=1 dB時,所得到的檢測波信號的IAVD基本一致。根據它與IAVDdenchmark(t)的差(圖13(a))以及方程式(10)得到有噪聲干擾時的SSP的輸出y(t)(圖 13(b)),最終所評估的TA-D-S對應于第1 930個采樣時間點。可知在無噪聲干擾時和有強噪聲干擾時,利用SSP算法所評估的TA-D-S基本一致。同樣結合所選用的六條傳感路徑各自評估的TA-D-S,繪制出的缺陷中心位置的坐標軌跡的交匯區域如圖14所示。所識別的中心位置坐標以及相對于實際中心位置的距離如表2中所示。

圖10 傳感路徑P3-P4所采集的具有不同SNR的檢測波信號(a)SNR=13 dB,(b)SNR=8 dB,(c)SNR=5 dB和(d)SNR=1 dBFig.10 The captured wave signals by sensing path P3-P4 from damage case with noise(a)SNR=13 dB,(b)SNR=8 dB,(c)SNR=5 dB and(d)SNR=1 dB

圖11 基準波信號與有噪聲干擾時的檢測波信號(SNR=1 dB)的頻率-能量譜Fig.11 Frequency-energy spectrum of the captured wave signals from benchmark and damage case with noise(SNR=1 dB)

圖12 所求得的基準波信號與檢測波信號(SNR=1 dB)的IAVDFig.12 IAVDbenchmark(t)and IAVDdamage(SNR=1dB)(t)

圖13 有噪聲干擾時的(a)Δ(t)和(b)SSP算法的輸出y(t)Fig13(a)Δ(t)and(b)the output of SSP algorithm(y(t))with noise(SNR=1 dB)

圖14 有噪聲干擾時(SNR=1 dB)所識別的缺陷位置Fig.14 Identified result with noise(SNR=1 dB)

實驗的損傷識別結果表明無論檢測波信號是否被強噪聲干擾,本研究所提出的SSP算法都可以精確地評估出損傷散射的S0模式ToF。所識別的缺陷中心位置的誤差均小于33 mm(檢測區域邊長的8.25%),在允許的誤差范圍內。

3 結論

本研究提出一個SSP算法,通過比較基準波信號的IAVD與檢測波信號的IAVD來評估損傷散射的S0模式的ToF。分別在理想(無噪聲)和嘈雜(有噪聲)的工作環境中檢測帶有切縫缺陷的鋁板,利用SSP求得檢測波信號的IAVD。并且通過適當的特征提取方法,與基準波信號的IAVD比較可以精確地評估出損傷散射的S0模式的ToF。結合傳感網絡中各條傳感路徑所評估的ToF,利用三角定位算法在傳感網絡所包圍的檢測區域內成功地識別出缺陷的位置。實驗結果表明,無論檢測波信號是否被強噪聲干擾,利用本研究所提出的SSP算法所求得的信號的IAVD基本保持一致,進而保證在較低信噪比的情況下仍然能夠精確地評估損傷散射的S0模式的ToF,實現有效的損傷識別與定位。本研究所證實的SSP算法強大的抵抗噪聲干擾的能力,表明它在基于超聲導波的損傷識別領域里具有很好的實用性。

[1]彭海闊,孟 光,李富才.一種基于壓電晶片陣列的板結構損傷識別方法[J].振動與沖擊,2009,28(9):56-59.

[2] Yang Y J,Cascante G,Polak M A.Depth detection of surface-breaking cracks in concrete plates using fundamental Lamb modes[J].NDT & E International,2009,42(6):501-512.

[3]Lu Y,Ye L,Su Z,et al.Quantitative assessment of throughthickness crack size based on Lamb wave scattering in aluminium plates[J]. NDT & E International,2008,41(1):59-68.

[4]張偉偉,王志華,宏 馬.含缺陷管道超聲導波檢測信號的相關性分析[J].暨南大學學報:自然科學與醫學版,2009,30(3):269-272.

[5] Wilcox P,Lowe M,Cawley P.The effect of dispersion on long-range inspection using ultrasonic guided waves[J].NDT& E International,2001,34(1):1-9.

[6] Li F,Meng G,Ye L,et al.Dispersion analysis of Lamb waves and damage detection for aluminum structures using ridge in the time-scale domain[J].Measurement Science and Technology,2009,20,095704:1 -10.

[7]Quek S T,Tua P S,Wang Q.Detecting anomalies in beams and plate based on the Hilbert-Huang transform of real signals[J].Smart Materials and Structures,2003,12(3):447-460.

[8] Zhang H,Fan S,Lu D.Application of Hilbert-Huang transform to arrival time extraction of multi-mode lamb waves[J].Journal of Vibration,Measurement and Diagnosis,2008,28(3):216-219.

[9]Li F,Meng G,Kageyama K,et al.Optimal mother Wavelet selection for lamb wave analyses[J].Journal of Intelligent MaterialSystems and Structures, 2009, 20(10):1147-1161.

[10]Ip K H,Tse P W,Tam H Y.Extraction of patch-induced Lamb waves using a wavelet transform[J].Smart Materials and Structures,2004,13(4):861-872.

[11] Jennifer E M.Detection,localization and characterization of damage in plates with an in situ array of spatially distributed ultrasonic sensors[J].Smart Materials and Structures,2008,17(3):1-15.

[12]Michaels J E,Michaels T E.Guided wave signal processing and image fusion for in situ damage localization in plates[J].Wave Motion,2007,44(6):482-492.

[13] Gustafsson M G.Nonlinear clutter suppression using split spectrum processing and optimaldetection[J]. IEEE Transactions on Ultrasonics,Ferroelectrics,and Frequency Control,1996,43(1):109 -124.

[14]Gustafsson M G,Stepinski T.Split spectrum algorithms rely on instantaneous phase information-a geometrical approach[J].IEEE Transactions on Ultrasonics,Ferroelectrics,and Frequency Control,1993,40(6):659 -665.

[15] BenammarA, DraiR, Guessoum A. Detection of delamination defects in CFRP materials using ultrasonic signal processing[J].Ultrasonics,2008,48(8):731 -738.

[16] Roig I B,Domínguez L V.New insights and extensions of split-spectrum algorithms from an optimum distributed detection perspective[C].IEEE Ultrasonics Symposium,2005,1797-1800.

[17] Ericsson L,Stepinski T.Cut spectrum processing:a novel signal processing algorithm for ultrasonic flaw detection[J].NDT & E International,1992,25(2):59-64.

[18]Diamanti K,Soutis C,Hodgkinson J M.Lamb waves for the non-destructive inspection of monolithic and sandwich composite beams[J].Composites Part A:Applied Science and Manufacturing,2005,36(2 SPEC.ISS.):189 -195.

[19]Kessler S S,Spearing S M,Soutis C.Damage detection in composite materials using Lamb wave methods[J].Smart Materials and Structures,2002,11(2):269 -278.

[20] Santoni G B,Yu L,Xu B,et al.Lamb wave-mode tuning of piezoelectric wafer active sensors for structuralhealth monitoring[J]. JournalofVibration and Acoustics,Transactions of the ASME,2007,129(6):752-762.

[21] Karpur P,Canelones O J.Split spectrum processing:a new filtering approach for improved signal-to-noise ratio enhancement of ultrasonic signals[J].Ultrasonics,1992,30(6):351-357.

[22] Karpur P,Shankar P M,Rose J L,et al.Split spectrum processing:determination of the available bandwidth for spectral splitting[J].Ultrasonics,1988,26(4):204 -209.

猜你喜歡
信號檢測
“不等式”檢測題
“一元一次不等式”檢測題
“一元一次不等式組”檢測題
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
“幾何圖形”檢測題
“角”檢測題
完形填空二則
孩子停止長個的信號
小波變換在PCB缺陷檢測中的應用
基于LabVIEW的力加載信號采集與PID控制
主站蜘蛛池模板: 成人免费一级片| 久久国产免费观看| 一本久道久综合久久鬼色| 亚洲欧美综合在线观看| 四虎精品国产永久在线观看| 亚洲国产看片基地久久1024| 白丝美女办公室高潮喷水视频| 国产精品99久久久久久董美香| 热99精品视频| 亚洲一区网站| 亚洲天堂网在线播放| 国产成人精品高清在线| 无码电影在线观看| A级全黄试看30分钟小视频| 欧美日韩中文国产va另类| 国产在线视频福利资源站| 久久国产精品无码hdav| 国产免费观看av大片的网站| 日韩中文无码av超清| 免费激情网站| 色网在线视频| 中国成人在线视频| 一级香蕉人体视频| 日韩二区三区无| 午夜三级在线| 欧美人在线一区二区三区| 国产精品区视频中文字幕| 91成人精品视频| 欧美日韩成人在线观看| 免费观看成人久久网免费观看| 中文字幕亚洲无线码一区女同| 伊人色综合久久天天| 97国产成人无码精品久久久| 99久久精品国产麻豆婷婷| 国产麻豆91网在线看| 玖玖精品视频在线观看| 国产午夜人做人免费视频中文 | 精品人妻无码中字系列| 91精品情国产情侣高潮对白蜜| 性视频一区| 国产h视频免费观看| 亚洲爱婷婷色69堂| 四虎综合网| 国内精品视频在线| 欧美色香蕉| 亚洲乱码精品久久久久..| 国产一在线| 亚洲欧美日韩成人高清在线一区| 性色一区| 国产一区在线观看无码| 国产粉嫩粉嫩的18在线播放91| 中日韩欧亚无码视频| 国产在线91在线电影| 欧洲极品无码一区二区三区| 国产精品久久自在自2021| 亚洲欧美一区在线| 亚洲第一成年人网站| 91久久偷偷做嫩草影院电| 亚洲综合色吧| 色老头综合网| 91久久国产热精品免费| 国产精品妖精视频| 亚洲无码视频图片| 免费国产高清视频| 亚洲精品片911| 动漫精品啪啪一区二区三区| 无码精品一区二区久久久| 亚洲国产中文在线二区三区免| 被公侵犯人妻少妇一区二区三区| 99精品在线视频观看| 午夜福利视频一区| 东京热高清无码精品| 乱人伦中文视频在线观看免费| JIZZ亚洲国产| 四虎影视库国产精品一区| 麻豆国产原创视频在线播放| 久久精品国产精品一区二区| 污网站在线观看视频| 91亚洲视频下载| 韩国v欧美v亚洲v日本v| 国产欧美精品午夜在线播放| 国产无码制服丝袜|