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

沖擊信號零漂修正的沖擊響應譜互相關系數分析

2016-09-18 02:45:47李海廣潘宏俠任海鋒
振動與沖擊 2016年16期
關鍵詞:信號

李海廣 , 潘宏俠, 任海鋒

(1.中北大學 機電工程學院,太原 030051; 2.內蒙古自治區白云鄂博礦多金屬資源綜合利用重點實驗室,包頭 014010)

?

沖擊信號零漂修正的沖擊響應譜互相關系數分析

李海廣1,2, 潘宏俠1, 任海鋒1

(1.中北大學 機電工程學院,太原030051; 2.內蒙古自治區白云鄂博礦多金屬資源綜合利用重點實驗室,包頭014010)

沖擊振動信號通常因為沖擊能量強、傳感器特性等因素的影響,存在基線漂移現象,進而影響到后續的數據分析處理。針對這一問題,首先通過計算沖擊信號的沖擊響應譜,提出正負沖擊響應譜互相關系數的定義,利用互相關系數開展檢測基線漂移程度、檢驗修正效果等工作;其次,提出利用響應譜互相關系數為重構條件的自適應經驗模式分解(EMD)漂移修正新方法,并與其他方法進行了比較。在對某機槍自動機動作過程中的沖擊振動基線漂移信號分析和驗證,數據結果表明,沖擊響應譜互相關系數較傳統的響應譜低頻斜率等響應譜特征值,更能準確有效的表征基線漂移程度、以及修正效果的評判。

基線漂移;沖擊響應譜; 互相關系數;經驗模式分解

沖擊振動信號在測量過程中往往出現基線漂移現象,影響到測量精度,甚至數據的獲取。基線漂移產生的主要原因是測量傳感器在高沖擊下產生的失效特性[1-2],如壓電式加速度傳感器的零漂現象[3-4]。對于基線漂移的數據最理想的處理是棄之不用[5],然而往往受到測量環境、測量設備等條件限制,測量數據具有唯一性[6],促使學者研究一些修正方法去除基線漂移,如采用小波分解的方法[5,7],EMD方法[8]等。然而對于沖擊信號基線漂移的研究內容中,還有兩個基本問題沒有解決,如何度量基線漂移程度、如何比較修正方法。對于這兩個問題的研究中,沖擊響應譜是一種很好的工具,沖擊響應譜的最大響應譜、速度譜及正負響應譜都涉及到沖擊信號的對稱性,已經有學者研究表明在雙對數坐標下最大響應譜低頻直線段斜率應滿足在6 dB/oct與12 dB/oct之間,低于6 dB/oct則存在基線漂移[9],正負譜在某些頻率的差值應小于6 dB/oct,否則存在基線漂移[10]。筆者在自動武器的自動機動作沖擊振動信號處理過程中,同樣遇到基線漂移的問題,經過研究比較,本文提出正負沖擊響應譜互相關系數的定義,利用互相關系數可以量度沖擊信號的基線漂移程度;提出利用互相關系數為重構條件的EMD自適應漂移修正方法,同樣可以利用互相關系數比較不同的漂移修正方法的修正效果。試驗數據表明響應譜互相關系數可以有效的量度偏移程度及比較修正方法的效果。

1 理論基礎

1.1沖擊響應譜

沖擊響應譜是基于沖擊加速度歷史數據的計算函數,主要思想是將沖擊激勵施加在一個標準的單自由度的質量彈簧阻尼振動系統,如圖1所示,每個單元由質量為mi的質塊,剛度為ki的彈簧,和一個阻尼為ci的阻尼器組成的單自由度系統,每個單元系統阻尼比ξ是相等的,沖擊加速度數據作為輸入信號,作用在系統基礎上,得到的每個質塊的加速度響應,并將最大加速度響應Gi作為質塊自然頻率ωi的函數構成沖擊響應譜。

圖1 沖擊響應譜原理示意圖Fig.1 Shock response spectrum model

將加速度信號作為基礎激勵源,質塊的運動方程:

(1)

式中:x為質塊的絕對位移,y為基座的位移,定義z=x-y為相對位移,式(1)轉化為

(2)

求解式(2),可以得到沖擊響應峰值與自然頻率之間的函數,即該系統的沖擊響應譜。目前應用較廣泛的沖擊響應譜算法是改進的遞歸數字濾波法[11],該方法于由SMALLWOOD[12]提出來。

根據沖擊作用時間,響應譜可以分為初始響應譜、殘余響應譜;同時根據響應參變量矢量方向,響應譜可以分為正譜和負譜,正負譜可以反映出沖擊系統的對稱性,每時刻正負譜的絕對峰值譜作為系統的最大響應譜,圖2為各響應譜的關系圖。

圖2 不同沖擊響應譜關系圖Fig.2 Links between the different definitions of SRS

圖3 不同頻率的正弦沖擊激勵信號Fig.3 Different shape on symmetric impact signal

1.2基線漂移的沖擊響應譜特征提取

沖激響應譜能夠反映沖擊信號對稱性關系,圖3(a)~(f)分別為頻率為50 Hz、100 Hz、200 Hz,400 Hz、600 Hz、800 Hz正弦沖擊信號,采樣頻率為10 000 Hz,采樣時長為0.03 s,其中沖擊信號時間間隔為0.01 s,形成不同對稱程度的沖擊信號。

分別計算各沖擊信號的沖擊響應譜,其中沖擊響應譜計算參數為:自然頻率初始值10 Hz,最大自然頻率的選擇采集頻率的1/8、倍頻為21/6,即自然頻率序列為10,10.595,11.225,11.892,…,100×21/6(N-1),這里N=85;阻尼比ξ= 0.05,響應譜的計算采用遞歸數字濾波法。

圖4為圖3正弦沖擊信號的最大沖擊響應譜,從圖中可以發現,不同對稱形式的沖擊信號,所對應最大沖擊響應譜不同。一般最大沖擊響應譜可以劃分為三個區域[13]:脈沖區域(impulse domain),響應譜幅值小于沖擊幅值區域;靜態區域(static domain) ,響應譜幅值在高頻區域趨向于沖擊幅值的區域;中間區域(intermediate),響應譜峰值是沖擊幅值動態放大,放大系數顯著的依賴于沖擊的形狀和系統的阻尼。圖4中各響應譜線在對數坐標下的直線段位于脈沖區域,該區域的斜率在6 dB/oct與12 dB/oct之間[9],比較小的斜率是由于傳感器缺陷產生基線漂移造成的,同時PIERSOL[10]認為由于背景噪音造成斜率小于6 dB/oct的沖擊信號是無效的。如圖所示,半波正弦函數譜線直線段的斜率明顯小于其他對稱正弦沖擊信號譜線。

圖4 不同沖擊激勵對應的最大響應譜Fig.4 Max shock response spectrum of different impact vibration signals

正負響應譜同樣可以反映出沖擊信號的對稱性,圖5為圖3沖擊激勵信號的正負響應譜,圖中半波正弦信號正負譜線明顯不同,但隨著沖擊信號對稱性增強,正負譜線趨于重合。為進一步定量衡量譜線重合度,本文定義正負響應譜線的互相關系數。

(3)

式中:Cov(SRS+,SRS-)為正負響應譜互相關函數,Var(SRS+)、Var(SRS-)為正負響應譜方差,分別計算系列沖擊激勵的正負譜線互相關系數,如圖6(a)所示,互相關系數除半坡正弦沖擊較低外,其余沖擊信號逐漸增大趨于一穩定值,同時計算響應譜低頻段斜率,圖6(b)中,系列沖擊信號響應譜斜率逐漸減小并趨于一穩定值,說明對于本文采用的系列正弦沖擊激勵信號,沖擊響應譜互相關系數同樣具有指示激勵信號對稱性的作用,可以反映基線漂移程度,而且相關系數是取值0~1之間的相對值,較正負譜具體差值方法[10]更有利于比較,同時由于隨著信號對稱性增強,相關系數趨于穩定值,可以根據信號特性設置相關系數有效閥值,相關系數在閥值之上可以認為信號基線漂移程度符合要求。

圖5 不同沖擊激勵對應正負響應譜Fig.5 Maximax shock response spectrum of different impact vibration signals

圖6 不同沖擊激勵沖擊響應譜特征提取Fig.6 Shock response spectrum features extraction from different impact vibration signals

1.3基線去除方法

目前,去基線漂移的濾波算法已經比較多,比如基于平滑濾波器的濾波算法,五點三次算法、形態濾波算法,基線擬合算法等,還有基于信號分解,將低頻分量分離消除基線漂移的濾波算法,如小波變換、EMD等,本文分別列舉五點三次算法、小波算法和EMD算法基本原理。

(1) 五點三次平滑算法

五點三次平滑算法是等間距數值基礎上的一種局部平滑的方法,在五點的局部數據窗口內,采用三次多項式擬合,確定加權平均系數,利用數據窗口對所有數據點進行平滑處理后,即平滑1次。在進行多次平滑后,可以有效消除數據的毛刺,得到數據光滑的趨勢項,如果將原數據減去趨勢項后,即可得到信號的高頻成份,所以五點三次平滑算法可以用于基線漂移的去除。

(2) 小波算法

利用小波變換多尺度多分辨率的特點,將信號進行多尺度小波分解。由于基線漂移的主要成分為緩變趨勢分量,在小波分解中會直接顯現于較大的尺度下,只要在重構過程中將這一尺度下的分量直接去除,便可以實現基線漂移的修正。

(3) EMD分解去基線

EMD分解[14]是由黃鍔博士提出的一種非線性信號分析方法,其原理為根據時間尺度特征對于信號進行分解,將復雜信號分解為有限個本構函數(IMF),本構函數包含了相應時間尺度局部信號特征,每個IMF為滿足以下兩個條件的局部信號:① 函數在整個時間范圍中,局部極值點和過零點的數目相同,或者是最多相差一個;② 在任意時間點,局部最大值包絡(上包絡線)以及局部最小值的包絡(下包絡線) 平均值是零。

EMD分解結果為:原信號x(t) 可表示為級若干IMF分量和一個趨勢項之和,如式(4)表示。

(4)

每個IMF為一個單一頻率范圍。EMD分解本質上具備高通濾波特征,可以將IMF由高頻到低頻依次的分解出來,而基線漂移是疊加在信號上的低頻分量,因此將低頻IMF去除后重構信號即可完成基線漂移的消除。

2 試驗數據分析

2.1自動機動作沖擊信號的獲取

(1) 自動機沖擊振動試驗裝置

本文試驗采用某高射機槍自動機,該自動機為導氣式自動機,其工作原理為利用彈藥在槍管內的高溫高壓燃氣推動活塞,由自動機構完成彈藥的擊發、擊發復位、開鎖、開閂、抽筒等循環過程。考慮到自動機主要動作構件及傳感器安裝條件因素,測點布置如圖7所示,測點1為槍尾上方,測點2為機匣左側沖擊振動信號由三向壓電加速度傳感器獲得,經放大后由便攜式DASP采集系統記錄,采集頻率為51.2 kHz,生成2個測點位置3個不同方向的6組沖擊加速度信號,數據模式根據射擊模式不同,可分為單發、三連發、五連發數據,相關傳感器參數可由表1獲得。

圖7 試驗機槍及測點分布Fig.7 Gun test and measurement point distribution

測點位置靈敏度/(pC·ms-2)XYZ型號測點1槍尾上方23.021.819.0CA-YD-116測點2機匣左側18.221.820.8CA-YD-116

(2) 沖擊振動信號的基線漂移

圖8為單發射擊模式下的沖擊振動時域圖,由振動信號可以看出,信號具有很強的沖擊特性,能量在撞擊后出現峰值,之后迅速地衰減,而每一次撞擊正是完成自動機動作過程中的一環節。沖擊振動信號包含了自動機工作時正常高頻振動信號,但也存在低頻的漂移干擾信號,如圖中1z,2z方向,產生的原理主要有,傳感器在劇烈沖擊振動中,進入非線性工作區,以及傳感器固定裝置產生松弛等原因。低頻漂移并未攜帶有用信息,同時會給后續的信號特征提取帶來不便,因此有必要去除基線漂移信號。

圖8 兩測點不同方向的沖擊信號Fig.8 Signals at each measurement point from different directions

(3) 基線漂移的沖擊響應譜分析

觀察圖8中的沖擊振動信號可以發現:測點1的z方向沖擊信號基線漂移比較嚴重,其他信號也存在不同程度的基線漂移。根據上文分析,可以通過計算沖擊響應譜提取響應譜特征來定量的比較,計算相應響應譜如圖9所示,最大響應譜低頻段斜率為4.821 6 dB/oct,正負譜線顯示出明顯的差異,其互相關系數為0.270 32。

圖9 基線漂移的沖擊響應譜分析Fig.9 Shock response spectrum analysis for baseline shift

計算并比較所有測量信號的沖擊響應譜特征值,圖10為計算結果,其中,根據PIERSO的分析,測點1x、2x方向的響應譜低頻斜率大于6 dB/oct,沖擊信號有效,其它方向均小于6 dB/oct,信號存在基線漂移,但從圖中可以發現,1z方向的低頻斜率并不是最小,表明低頻斜率特征值并不能用于漂移程度的比較;相反,正負譜線互相關系數可以量化漂移程度的不同,各測點相關系數依次為0.983 246、0.991 672、0.270 32、0.989 849、0.990 441、0.936 797,如圖10(b)所示,其中1z方向0.270 3明顯低于其它方向,而且除1z測點外,2z方向相關系數0.936 797也低于其它測量值,根據觀察可以設置沖擊信號有效閥值為0.990 00,同時相關系數的比較關系完全符合圖8的觀察結果。

圖10 測量沖擊信號基線漂移響應譜特征值比較Fig.10 Comparison chart of shock response spectrum features for impact signals

2.2基線漂移的修正效果比較

(1) 五點三次平滑漂移修正

利用五點三次平滑算法對測點1的z方向沖擊信號平滑1 000次,得到趨勢項后,計算原信號與趨勢項的差值得到漂移修正信號,并計算沖擊響應譜,如圖11所示,其中最大響應譜低頻段斜率為6 dB/oct,正負譜線互相關系數為0.998 861,而且從圖中可以看出正負譜線幾乎完全重合,說明趨勢項已經很大程度上被消除。

圖11 五點三次平滑漂移修正的沖擊響應譜分析Fig.11 Shock response spectrum analysis for baseline correction using five-spot triple smoothing algorithm

(2) 基于小波漂移修正

采用離散Meyer小波,將測點1的z方向沖擊信號,進行10尺度小波分解,提取第10尺度上的低頻小波分解系數,即趨勢項系數,將系數置零后重構沖擊信號,并計算沖擊響應譜,如圖12所示。其中最大響應譜低頻段斜率為6.045 dB/oct,正負譜線互相關系數為0.881 497,而且從圖中可以看出正負譜線部分重合,說明趨勢項并沒有被完全消除。

圖12 小波漂移修正的沖擊響應譜分析Fig.12 Shock response spectrum analysis for baseline correction using wavelet algorithm

(3) 基于EMD漂移修正

根據本文1.3節的分析,運用EMD分解方法對測點1的z方向基線漂移信號分解,根據RILLING等[15]的收斂條件,共提取15個IMF和一個趨勢項,如圖13所示。

圖13 沖擊振動信號EMD分解后的所有IMFFig.13 EMD decomposition of a impact vibration signal

去除趨勢項和低頻IMF可以有效的消除基線漂移,然而如何選擇低頻IMF的范圍成為新問題,本文提出一種自適應的選擇方法,具體為:從消除趨勢項開始,每消除一項IMF,計算重構沖擊信號的正負沖擊響應譜,比較正負譜線的互相關系數,最大互相關系數所對應的重構信號為消除基線漂移的輸出信號,整個過程如圖14所示。

圖14 正負相應譜最大互相關系數選擇重構信號Fig.14 Reconstruct signal based on shock response spectrum correlation coefficients

圖15(a)、(b)分別為去除趨勢項、趨勢項至第五個低頻IMF沖擊響應譜分析,對比圖9原漂移信號及其沖擊響應譜,可以看出重構信號的基線漂移得到了明顯的改善,計算相應信號的最大響應譜低頻斜率分別為-0.097 8、2.204;正負譜線互相關系數為0.418 9、 0.902 1。

圖15 去除趨勢項、趨勢項至第五個低頻IMF沖擊響應譜分析Fig.15 Shock response spectrum analysis for trend item extraction and trend item to No.5 IMF extraction

根據上文的自適應方法,從去除趨勢項開始,依次計算去除每一個IMF后重構信號沖擊響應譜的低頻斜率和正負譜互相關系數,如圖16(a)、(b)所示,隨著低頻分量的去除,響應譜低頻斜率趨于12附近,而互相關系數逐漸增加,在去除到第11個IMF時達到最大值0.998 500,其后的相關系數分別為0.997 820,0.996 890、0.994 196、0.990 139、0.990 059逐漸降低,這是因為重構信號高頻分量有著趨向于高頻噪音的趨勢,響應譜趨于高頻噪音的響應譜形式。圖16(c)~(e)為相關系數最大值對應的輸出重構信號、最大響應譜、正負響應譜,如圖所示,重構信號消除了基線漂移保留了信號的高頻成份。

圖16 去除低頻IMF后重構信號的響應譜分析Fig.16 Shock response spectrum analysis for reconstruct signal after low frequency IMF extraction

(4) 五點三次、小波、EMD漂移修正比較

根據上文的分析,五點三次、基于EMD的基線漂移修正算法均可以消除基線漂移,而小波修正算法部分消除了基線漂移,通過修正信號的正負沖擊響應譜互相關系數比較,EMD算法的相關系數為0.998 400,五點三次算法的0.998 861,小波的相關系數為0.881 497,說明EMD算法的基線漂移修正效果基本和五點三次算法一致,強于小波修正方法,不同的是EMD算法在消除了低頻漂移分量的同時對信號進行了分解,更有助于信號的分析及特征的提取。

3 結 論

沖擊振動信號,受到傳感器特性等因素的影響,沖擊信號往往存在不同程度的基線漂移現象。本文提出通過計算沖擊振動信號的沖擊響應譜互相系數,開展信號漂移程度的量度及消除漂移成份效果評價的新方法。在本文開展的試驗中,響應譜互相關系數較最大沖擊響應譜的低頻段斜率方法,更有效的對于基線漂移程度進行了定量描述;利用響應譜互相關系數,有效的選擇EMD基線修正算法中去除本構函數的最優值,同時利用互相關系數比較了五點三算法、小波算法和EMD算法基線漂移的去除效果。試驗數據證明,該方法可以有效的開展沖擊信號基線漂移的程度判定和消除方法比較等相關工作。

[1] 趙小龍,馬鐵華,范錦彪. 高g值加速度計在高沖擊下的失效特性的研究[J]. 傳感技術學報,2012,(12):1668-1672.

ZHAO Xiaolong, MA Tiehua, FAN Jinbiao. The study of failure characteristic of highgacceleration in high shock[J]. Journal of Transduction Technology,2012,(12):1668-1672.

[2] 張國偉,馮順山,俞為民. 高沖擊過載加速度傳感器零漂分析[J]. 華北工學院學報, 2004(1): 64-67.

ZHANG Guowei, FENG Shunshan, YU Weimin.Analysis of zero drift of the acceleration sensor in high-impulsion and high overloading[J]. Journal of North China Institute of Technology, 2004(1): 64-67.

[3] 夏偉強,馬鐵華,范錦彪,等. 壓電式加速度傳感器在高沖擊環境下的零漂分析[J]. 傳感技術學報, 2007(7): 1522-1527.

XIA Weiqiang, MA Tiehua, FAN Jinbiao, et al.Analysis of zero drift of the piezoelectric acceleration sensor in high impact testing[J]. Journal of Transduction Technology, 2007(7): 1522-1527.

[4] 宋穎,杜彥良,孫寶臣. 壓電傳感技術在輪軌力實時監測中的應用探討[J]. 振動與沖擊, 2010,29(1): 228-232.

SONG Ying, DU Yanliang, SUN Baochen. Travel rate analysis for a non-harmonically vibrating conveyor[J].Journal of Vibration and Shock, 2010,29(1): 228-232.

[5] 陳建軍,王樹樂. 沖擊加速度零漂的一種改進小波修正方法[J]. 水雷戰與艦船防護, 2012(4):18-23.

CHEN Jianjun, WANG Shule. An improved wavelet correction for zero shifted accelerometer data[J]. Mine Warfare & Ship Self-Defence, 2012(4):18-23.

[6] 袁宏杰,姜同敏. 實測爆炸分離沖擊數據的分析和處理[J]. 固體火箭技術, 2006(1): 72-74.

YUAN Hongjie, JIANG Tongmin.Analysis and treatment of measured pyrotechnic shock data[J]. Solid Rocket Technology 2006(1): 72-74.

[7] EDWARDS T S. An improved wavelet correction for zero shifted accelerometer data[J]. Shock and Vibration, 2003, 10(3): 159-167.

[8] 胡燦陽,陳清軍. 基于EMD和最小二乘法的基線飄移研究[J]. 振動與沖擊, 2010, 29(3): 162-167.

HU Canyang, CHEN Qingjun.Research on baseline drift using least-square and EMD[J]. Journal of Vibration and Shock, 2010, 29(3): 162-167.

[9] WRIGHT C. Effective data validation methodology for pyrotechnic shock testing[J]. Journal of the IEST,2010,53(1): 9-30.

[10] PIERSOL A. Guidelines for dynamic data acquisition and analysis[J]. Journal of the IES, 1992,35(5): 21-26.

[11] 華師韓,田恒春. 沖擊響應譜計算相關參數選擇的研究[J]. 遙測遙控, 2005(6): 50-55.

HUA Shihan, TIAN Hengchun. Consideration of the parameters choice for shock response spectrum evaluation [J]. Telemetry and Telecontrol, 2005(6): 50-55.

[12] SMALLWOOD D O. Improved recursive formula for calculating shock response spectra[C]//Shock and Vibration Symposium. San Diego, CA, 1980.

[13] LALANNE C.‘Properties of shock response spectra’ in mechanical shock[M].Hoboken, NJ:John Wiley & Sons, Ltd, 2014:103-173.

[14] HUANG N E,SHEN Z,LONG S R,et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J]. Proceedings of the Royal Society a—Mathematical Physical and Engineering Sciences, 1998, 454(1971): 903-995.

[15] RILLING G,FLANDRIN P,GONCALVES P. On empirical mode decomposition and its algorithms[C]. NSIP-03, Grado (Ⅰ),2003:8-11.

Baseline correction of impact signals using the cross-correlation coefficient of shock response spectrum

LI Haiguang1,2, PAN Hongxia1,REN Haifeng1

(1. College of Mechatronic Engineering, North University, Tai Yuan 030051, China;2 . Inner Mongolia Key Laboratory for Utilization of Bayan Obo Multi-Metallic Resources:Elected State Key Laboratory, Baotou 003310, China)

The baseline shift components in the impact vibration signals are associated with the high impact energy and sensor characteristics and they influence data analysis. This paper presented a new method for detecting the shifted baseline level and correction of the baseline shift. The proposed method was based on the cross-correlation coefficient between the positive shock response spectrum and the negative shock response spectrum. The new method for the baseline correction of impact vibration signal which represents the reconstruction of signal was realized by using the empirical mode decomposition(EMD) based on the cross-correlation coefficient. The validity and precision of the proposed method was demonstrated with analysis and validation of the signals generated in the operation process of an automation gun .

baseline shift; shock response spectrum; cross-correlation coefficient; empirical mode decomposition

國家自然科學基金資助項目(51175480)

2015-10-30修改稿收到日期:2016-02-17

李海廣 男,博士生,副教授,1975年12月生

潘宏俠 男,教授,1950年10月生

TJ06

A

10.13465/j.cnki.jvs.2016.16.035

猜你喜歡
信號
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
7個信號,警惕寶寶要感冒
媽媽寶寶(2019年10期)2019-10-26 02:45:34
孩子停止長個的信號
《鐵道通信信號》訂閱單
基于FPGA的多功能信號發生器的設計
電子制作(2018年11期)2018-08-04 03:25:42
基于Arduino的聯鎖信號控制接口研究
《鐵道通信信號》訂閱單
基于LabVIEW的力加載信號采集與PID控制
Kisspeptin/GPR54信號通路促使性早熟形成的作用觀察
主站蜘蛛池模板: 日韩久草视频| 国产成人无码Av在线播放无广告| 亚洲国产精品不卡在线| 国产成人福利在线| 米奇精品一区二区三区| 毛片网站观看| 亚欧乱色视频网站大全| 亚洲综合久久成人AV| 无码电影在线观看| 欧美在线一二区| 国产流白浆视频| 亚洲精品无码日韩国产不卡| 国内熟女少妇一线天| 国产91精品久久| 成人噜噜噜视频在线观看| 亚洲最大综合网| 国产凹凸视频在线观看| 国产精品久线在线观看| 国内精品免费| 欧美日韩北条麻妃一区二区| 好吊妞欧美视频免费| 欧美色综合网站| 午夜毛片免费看| 99中文字幕亚洲一区二区| 欧美丝袜高跟鞋一区二区| 99免费视频观看| 国产精品熟女亚洲AV麻豆| 自拍偷拍欧美日韩| 国产成人精品男人的天堂| 一级爆乳无码av| 国产丰满大乳无码免费播放| 青青草一区| 亚洲欧美天堂网| 国产成人夜色91| 欧美一级99在线观看国产| 免费Aⅴ片在线观看蜜芽Tⅴ| 亚洲欧美另类中文字幕| 国产精品欧美激情| 成人亚洲视频| 色偷偷一区二区三区| 国产成人精品三级| 国产内射一区亚洲| 成人欧美日韩| 国产一区免费在线观看| 亚洲一区二区精品无码久久久| 97精品国产高清久久久久蜜芽| 亚洲AⅤ永久无码精品毛片| 成人永久免费A∨一级在线播放| 国产91丝袜在线播放动漫| 国产精品无码AV片在线观看播放| 国产人成乱码视频免费观看| 97国产精品视频自在拍| 58av国产精品| 真实国产精品vr专区| 亚洲色成人www在线观看| 亚洲AV免费一区二区三区| 99热最新在线| 一级毛片基地| 国产亚洲精品自在线| 国产成人久视频免费| 六月婷婷综合| 毛片久久网站小视频| 日韩毛片基地| 亚洲人成日本在线观看| 精品三级在线| 欧美成人日韩| 国产在线精品99一区不卡| 久久人体视频| 亚洲一区二区在线无码| 国产sm重味一区二区三区| 高潮爽到爆的喷水女主播视频 | 亚洲欧美日韩中文字幕在线| 第一区免费在线观看| 毛片一级在线| 国产亚洲欧美在线专区| 8090成人午夜精品| 国产中文在线亚洲精品官网| 911亚洲精品| 亚洲Aⅴ无码专区在线观看q| 国产成人一区二区| 在线观看无码av五月花| 国产h视频免费观看|