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

沖擊加速度校準中微分過程及峰值濾波的問題討論

2017-11-30 06:01:20梁志國
振動與沖擊 2017年22期

梁志國

(北京長城計量測試技術研究所 計量與校準技術重點實驗室,北京 100095)

沖擊加速度校準中微分過程及峰值濾波的問題討論

梁志國

(北京長城計量測試技術研究所 計量與校準技術重點實驗室,北京 100095)

針對于激光干涉法沖擊加速度校準中,微分運算之前的濾波問題,進行了實驗研究,提出了一種基于沖擊速度曲線階躍上升時間參數為基礎的,局域平均濾波差分方式,實現了沖擊速度的獲取。從而避免了ISO相關標準中基于沖擊加速度波形脈寬調控濾波器截止頻率的不便之處,并對參數的取值范圍進行了實驗研究。通過與ISO標準方法的對比實驗,驗證了該數據處理方法可獲得的峰值更高些,可以高10%~20%左右。在不同傳感器上所作的實驗驗證了該方法的有效性與可行性,為沖擊加速度的計算提供了另外一種有效實用的數據處理方式。

計量學; 沖擊加速度; 微分; 校準

多普勒效應是有相對運動的物體間信號傳遞時體現的一種自然物理現象,自從被發現以來,獲得了人們的廣泛研究和利用,特別是在測量、測試領域,具有不可替代的地位與作用,在沖擊計量和校準方面,其作用尤其突出。當今世界上最高水平的沖擊測量與校準裝置,均是借助于穩頻激光技術和激光多普勒原理實現的[1-5]。其中,文獻[1]和文獻[3]所限定的加速度量程均為102~105m/s2,不確定度為2%;文獻[5]限定的加速度量程為102~105m/s2,不確定度沒有具體規定;文獻[2]是還在開展的工作,尚未形成正式文件,其本質則為文獻[1]所述量程的拓展,使用標準加速度方式,已經將Hopkinson棒激勵的沖擊加速度量程拓展到2 000 km/s2峰值,但不確定度降為6%(速度比較法)和10%(加速度比較法)。其數據濾波及處理方式與文獻[1]相比沒有特別變化。

其中,干涉測量的基本思想,主要是借助于激光干涉條紋或光柵衍射條紋,在與光電探測器存在相對運動條件下,運動狀態被光電探測器接收,產生實時頻率隨運動速度成正比的多普勒頻移,從而導致光電探測器輸出的信號是一種包含實時多普勒頻移的調頻信號,該多普勒頻移表述的運動速度波形,可以使用FM(Frequency Modulation)信號解調方法解調獲得。有多種FM解調方法可供使用,例如,滑動平均方法[6-8],希爾波特變換方法等。

由該多普勒頻移可以獲知運動物體與光電探測器之間相對運動的速度波形,將該速度對時間進行微分運算,即可獲知沖擊加速度波形。由此可見,在沖擊加速度波形的光電測量與校準中,微分運算是一個不可避免的環節,若使用希爾波特變換法進行波形解調時,則需要用到先后兩次微分運算。

與正弦振動不同,沖擊實際上是一種沒有精確模型的運動過程,盡管通常稱其為半正弦波形,但實際上它與半正弦波形間的差異可能非常巨大。尤其在高沖擊和超高沖擊過程中,由于受沖擊體的形狀、材料的非均勻、塑性形變、測量噪聲、傳感器件的非線性、測量的離散性等眾多因素的影響和作用,使得人們很難真正確認每一次沖擊的真實波形到底是哪種樣子的,直接導致了簡單的相鄰數值的差分運算無法獲得真實有效的加速度波形。必須進行濾波、加窗或采取其它措施進行處理,而濾波器及濾波參數的選取,以及加窗方式和窗函數的選取,一直困擾著沖擊計量和校準行業,也極大制約了沖擊計量準確度的提高。因為每變化一下濾波參數,都會導致沖擊加速度峰值的變動,從而給出不同的計量結果。針對這一問題,國際標準化組織ISO專門在其所發布的標準里明確規定了濾波器的種類、階次、以及截止頻率的選取方式。但仍然存在一些實際問題,例如,濾波參數需要用到微分后才能獲得的加速度脈沖寬度T,導致出現倒果為因的操作過程;并且需要選取4階Butterworth濾波器進行濾波,對于近方形沖擊加速度波形,其截止頻率不大于16/T;對于近高斯形沖擊加速度波形,其截止頻率不大于5/T。需要有極豐富的實踐經驗才能合理選取參數值。另外,在超高沖擊條件下,所獲得的沖擊波形通常即非高斯形狀,也非近方波形狀,而往往體現出子母峰等復雜形狀,更加難以確定合適的濾波參數。

本文后續內容,將主要針對這些問題以及其在微分過程中的不便展開討論,并試圖給出一種新的方法和簡潔判據,用以處理和應對沖擊加速度計量校準中的微分問題。

1 沖擊加速度校準的一般過程

典型的沖擊加速度校準裝置主要由Hopkinson沖擊機、差動式光柵激光干涉儀、數字示波器,以及計算機等設備組成,如圖1所示。

(1) 使用Hopkinson沖擊機發射彈體,利用Hopkinson棒在校準端面上產生波形良好的沖擊加速度過程,直接加載到光柵和被校準的加速度計上。

(2) 使用差動式激光干涉儀,利用光電探測器檢測光柵運動產生的激光多普勒信號x(t)。

(3) 使用數字示波器對激光多普勒信號x(t)以及經過送信號適調儀調理后的被校加速度計輸出信號進行數據采集。

(1)

式中,h(n)為Hilbert變換器的單位沖擊響應

(2)

x(n)的瞬時幅度

(3)

x(n)的瞬時相位

(4)

(5) 對式(4)所示的實時相位波形通過疊加kπ(k=0, 1, 2,…)進行相位展開,獲得連續相位波形。

(6) 對相位展開后的連續相位波形進行數字濾波;濾波器為4階Butterworth濾波器,截止頻率上限為16/T(對于近方波脈沖)~5/T(對于近高斯形狀脈沖),T為沖擊加速度脈寬。

(7) 對上述濾波后的連續相位波形進行微分運算,解算出多普勒頻移信號Δf(i)

(5)

式中,Δt為采樣序列的采樣間隔。

使用差動式衍射激光干涉儀進行沖擊測量,如圖1所示。其運動速度波形信號v(i)

(6)

式中:d為光柵的柵距;p、q分別為差動激光干涉儀所使用的光柵衍射級數。

(8) 對式(6)所述的沖擊速度波形進行數字濾波;濾波器為4階Butterworth濾波器,截止頻率上限為16/T(對于近方波脈沖)~5/T(對于近高斯形狀脈沖),T為沖擊加速度脈寬。

(9) 將上述濾波后的速度波形微分得到沖擊加速度a(i)

(7)

圖1 沖擊加速度校準裝置框圖Fig.1 Primary shock calibration unit of acceleration

2 微分過程的問題及討論

前已說明,由于量化誤差、隨機噪聲誤差、以及其它許多因素的影響,人們獲得的光電探測器輸出的光電信號并不光滑,而經過希爾波特變換處理及相位展開后的連續相位波形依然如此。直接進行相鄰點差分運算無法獲得有效可用的微分波形。因此,ISO標準針對該問題給出的解決方式是使用4階Butterworth濾波器進行濾波。并規定,若最終沖擊加速度波形的脈寬為T,則該濾波器的截止頻率上限為16/T~5/T。當沖擊加速度波形為近似方波的情況下,濾波器截止頻率上限選16/T,而當沖擊加速度波形為近似高斯脈沖形狀時,濾波器截止頻率上限選5/T。這一原則同時適用于從相位波形到頻率波形的微分運算,以及從速度波形到加速度波形的微分運算。該方式在實際執行中主要面臨以下幾個問題:

(1) 人們在進行速度信號微分以獲取加速度波形時,尚未知曉加速度波形的真正形狀,更加無從知曉其“脈寬”等參數。因而,ISO的上述推薦方法在邏輯順序上存在倒果為因問題。這導致人們在使用該方法時,往往先用種種手段“預估計”一下沖擊加速度的波形,獲得其脈寬的預估計值,以此為微分前的濾波器預置濾波參數?;蛘吒纱嘀苯邮褂脹_擊加速度傳感器的輸出波形的脈寬作為微分前的濾波器預置濾波參數。

無論哪一種方式,都會帶來額外的負擔及運算成本,并給濾波參數的選擇帶來不確定性因素。

(2) 用16/T~5/T原則選取Butterworth濾波器截止頻率,在低量程沖擊加速度校準時不會出現大問題,此時,沖擊加速度波形往往接近半正弦的高斯形狀。在高沖擊及超高沖擊加速度校準時,其加速度波形的不對稱程度加大,甚至出現子母峰(母峰上長出一個細小的子峰)等現象,此時濾波器的截止頻率的變化對任何微分結果的加速度峰值都會產生較大影響,波動可達10%左右。而如何選取濾波器參數在此時并無明確結論。

實際上,沖擊速度波形是一個近似階躍波形,以階躍上升時間作為參數來限定微分運算前的濾波參數將更加合理可行,且避免了倒果為因的邏輯問題。另外,每執行一次積分運算,都會使信號的有效頻譜寬度變窄,而每執行一次微分運算也會使信號的頻譜變寬。因而,以微分后的信號波形脈寬來調整和限定微分前的濾波器截止頻率本身也存在因果關聯不夠明確的問題。而以同一種濾波參數來分別處理由相位鑒頻和由速度計算加速度兩處頻譜截然不同的信號的做法同樣顯得有些粗放。在不同環節上,理應使用不同特征參數的濾波器。

(3) Butterworth濾波器以及其它眾多濾波器,其特性通常以頻域特征參數體現和表述,其局域的時域特征表現并不明顯,且與濾波參數缺乏一一對應關系。而沖擊加速度計量校準主要關注的是其峰值、脈寬、上升時間等時域特征參量。針對每一個具體波形,人們通常很難明確濾波器參數變化與導致的其波形參數變化的對應關系。因此,導致人們在能不使用濾波器的情況下,盡量不去使用濾波器,以減低測量結果的不確定性,提高測量準確度。

3 局域平均法實現微分的具體過程

為了避免ISO標準推薦的濾波方法中,濾波器截止頻率參數選取中的倒果為因問題,本文提出一種直接時域局部平均方式進行微分運算,用于實現鑒頻微分和求加速度微分,以提高測量效率,降低結果的波動。其過程如下:

(1) 執行上述第“1”節中的(1)~(5),獲得連續相位波形;其中,使用式(1)進行希爾波特變換運算時,卷積區間為[-M, +M],M應當使卷積窗口寬于載頻信號的周期,且M≥100。

(2) 若圖1中的差動式激光干涉儀為外差式,則須對其上述相位展開后的相位序列φ1,φ2, …,φn進行線性趨勢剔除工作,使用端基直線法獲得參數G0和φ0

(8)

φ0=φ1-G0×1

(9)

計算獲得剔除線性趨勢后的相位序列φi

φi=φi-G0·i-φ0(i=1, …,n)

(10)

若圖1中的差動式激光干涉儀為零差式,則不必剔除線性趨勢。有φi=φi,(i=1,…,n)。

(3) 選取與希爾波特變換卷積區間相同的差分點數2M+1,進行滑動平均方式多點差分運算,獲得瞬時頻率值fi(i=1,…,n)

(11)

(4) 按式(6)計算沖擊速度序列v1,v2,…,vn;它是一個典型的階躍波形。

(5) 用眾數法計算沖擊速度波形序列的上升時間tr。

(12)

為了判斷L取值是否最優,采取如下判定方式:

將平均值窗口內的速度值序列vi-L,vi-L+1, …,vi+L-1,vi+L,做端基直線擬合,得

v(k)=k·g+v0

(13)

(14)

以局域平均窗口中的誤差均值Δv(t)曲線作為微分中的誤差指針,可見,當窗口由大到小時,Δv(t)曲線在加速度峰值附近的曲線值呈下降趨勢,且有明顯的局部極小值,當該極小值下降不明顯,且局部極小值特征微弱時,即可判定窗口選擇值已經接近最佳。

4 實驗驗證

在如圖1所示的沖擊加速度校準裝置上進行校準實驗。其中:

數字示波器型號為TDS544A,帶寬500 MHz,8 Bits A/D,雙通道采樣工作方式,垂直電壓設置為0.1 V/div,通道采樣速率500 MSps,采集數據點數n=231 001。該裝置里差動激光干涉儀所使用的光柵柵距d=1/300 mm,式(6)所用衍射級數分別為p=1,q=-1。

圖2為使用圖1所示校準裝置對某沖擊加速度傳感器,按照上述第“1”節中的(1)~(3)進行校準實驗獲得的沖擊過程中激光多普勒干涉信號。

圖2 沖擊過程中的激光多普勒干涉信號Fig.2 Laser doppler interference signal of impulse

經過第“1”節中的(4)~(9)后獲得的沖擊加速度曲線波形,如圖3所示。加速度峰值為1 048 301.679 8 m/s2,此時,Butterworth濾波器截止頻率選取10/T,T為脈沖的脈寬。

圖3 沖擊過程產生的加速度波形曲線Fig.3 Acceleration waveform of impulse

將圖2所示測量曲線按照第“3”節(1)~(4)所述過程和方法進行處理,獲得沖擊速度曲線波形,如圖4所示。從圖4可得其沖擊速度波形的上升時間為:tr=92.16 μs(46 080個采樣間隔)。

圖4 沖擊過程產生的速度波形曲線(tr=92.16 μs)Fig.4 Velocity waveform of impulse

選取不同寬度的局域平均窗口,按照第“3”節中的(6)所述計算加速度波形曲線,如圖5~圖12所示。它們的沖擊峰值最大值和最小值相差18%。由此可見,平均窗口的變化對加速度峰值的影響是不可忽視的。對于比較不同傳感器特性來說,應該使用相同的平均值窗口執行運算操作才更有意義和價值。

從表面上看,平均窗口長的圖5濾除了更多的細節特征,從而使得沖擊曲線更加光滑,而平均窗口越短的圖則保留了更多的細節特征,圖8~圖12中甚至出現了子母峰現象。但是,由于物理機理掌握尚不夠深入,很難判定圖5和圖11哪一個更加接近真實的沖擊加速度波形。這一直是沖擊計量校準中最難以解決的問題之一。

盡管如此,從圖12(a)和圖5(a)相比,可以看出,僅僅由于濾波參數的選擇,峰值變化可以達到20%之多!

而從圖12(a)和圖3相比可知,本文所述方法與ISO標準所列方法可以有近10%的峰值差異。由于大多數濾波器對于沖擊峰值的作用都是將峰值幅度降低,因此,有理由相信,本文方法獲得的子母峰現象是真實的,且所獲得峰值應該比圖3所示用ISO標準獲得的量值更接近真實峰值。

圖5 沖擊加速度波形曲線和局部速度偏差均值Δv(t)(1/2上升時間窗口)Fig.5 Acceleration waveform of impulse and local error Δv(t) of velocity of impulse(1/2 rise time)

圖6 沖擊加速度波形曲線和局部速度偏差均值Δv(t)(1/3上升時間窗口)Fig.6 Acceleration waveform of impulse and local error Δv(t) of velocity of impulse(1/3 rise time)

圖7 沖擊加速度波形曲線和局部速度偏差均值Δv(t)(1/4上升時間窗口)Fig.7 Acceleration waveform of impulse and local error Δv(t) of velocity of impulse(1/4 rise time)

圖8 沖擊加速度波形曲線和局部速度偏差均值Δv(t)(1/6上升時間窗口)Fig.8 Acceleration waveform of impulse and local error Δv(t) of velocity of impulse(1/6 rise time)

圖9 沖擊加速度波形曲線和局部速度偏差均值Δv(t)(1/7上升時間窗口)Fig.9 Acceleration waveform of impulse and local error Δv(t) of velocity of impulse(1/7 rise time)

圖10 沖擊加速度波形曲線和局部速度偏差均值Δv(t)(1/8上升時間窗口)Fig.10 Acceleration waveform of impulse and local error Δv(t) of velocity of impulse(1/8 rise time)

圖11 沖擊加速度波形曲線和局部速度偏差均值Δv(t)(1/9上升時間窗口)Fig.11 Acceleration waveform of impulse and local error Δv(t) of velocity of impulse (1/9 rise time)

圖12 沖擊加速度波形曲線和局部速度偏差均值Δv(t)(1/12上升時間窗口)Fig.12 Acceleration waveform of impulse and local error Δv(t) of velocity of impulse (1/12 rise time)

本文中,為了判定濾波窗口的優劣,特別引入窗口內局部速度偏差均值曲線Δv(t),從圖5(b)~圖12(b)可知,隨著窗口的變小,加速度峰值附近的局部小區間內的Δv(t)呈下降趨勢,且存在局部極小值特征,若極小值底部比較尖銳,則表明平均值窗口寬度過寬,應減小窗口寬度,直至出現比較平緩的局部極小值特征,或者直到某一寬度以后,局部極小值便穩定在一個范圍內波動,局部極小值特征不再明顯,這時,可以認定,再減小濾波窗口已經無意義了,此時的窗口接近最佳。

令橫坐標為上升時間與濾波窗口寬度之比r,則有圖13所示的Δv(t)-r關系曲線圖。

圖13 Δv(t)-r關系曲線圖Fig.13 Δv(t)-r curve

從圖13可知,隨著r值的增大,Δv(t)總體呈下降趨勢,但下降速度到一定程度后趨于穩定而不再下降,本實驗中,r=8以后,局部速度偏差Δv(t)變化趨穩,此時,加速度變化也趨于緩慢??烧J為獲得的加速度復現性較好,量值比較準確可靠。

通常,認為曲線急速下降至較為平緩的分界拐點,是濾波窗口的特殊取值點。本實驗中,r=4為分界拐點,此時窗口長度為上升時間的1/4。

比較平均值窗口取上升時間1/4的圖7與圖3所示曲線峰值[9],兩者相差僅為1.44%,是一個非常小的值,可見,此時,本文方法的分界拐點與ISO推薦的方法具有良好的符合性。

5 討 論

通過上述過程可見:

(1) 使用本文所述方法進行沖擊加速度計算,每一步微分運算所用到的參數均是在當時就可以從波形參數中直接獲取的,未使用需要后續處理才能獲得的參數作依據,沒有ISO標準中存在的倒果為因的邏輯問題。

(2) 直接使用具有階躍特征的沖擊速度波形的上升時間參數確定局域平均值窗口的寬度,不必判定其沖擊加速度波形是否介于方波和高斯形狀波形的特點,具有統一性和一致性,更加有利于執行微分操作。

因為無論其最終的加速度波形的形狀如何,其沖擊速度波形都可以看作一個典型的階躍響應波形,具有明確的階躍響應波形特征。

(3) 直接使用平均值方式進行濾波,相比于Butterworth濾波器等其它形式濾波器而言,本身屬于時域局部濾波方式,其時域濾波特性簡潔直觀,特征明顯,濾波的影響與不足分析評價容易,尤其是它給微分后的加速度峰值造成的影響可以簡單估計和定量給出,是本文方法相較于其它方法最突出的優點。

(4) 使用局域平均方式執行微分運算,用速度曲線微分計算加速度時,給后續的沖擊加速度峰值帶來的影響,直接與所選取的平均值窗口長度內最大斜率部分階躍曲線段的線性度相關聯。該線性度已知后,即可以簡單判定本條沖擊曲線峰值加速度的波動范圍。其加速度峰值波動偏差的上界為最大非線性偏差除以平均值窗口的時間寬度,而波動偏差的下界則是窗口內相鄰點最大偏差除以采樣間隔Δt。

(5) 通過變化平均值窗口的長度,在速度曲線最大斜率附近可以按照速度偏差值Δv(t)下降趨勢進行尋優,找到令速度偏差Δv(t)下降趨緩的拐點時的窗口長度L,以2L~3L作為最佳平均值濾波窗口長度,從而能夠實現一種最優濾波窗口的選取。

若不進行尋優計算,通??梢赃x取窗口長度為上升時間的1/3~1/12。

(6) 通過與ISO標準推薦方法相比較可見,本文所提出的局域平均濾波微分方法,與ISO推薦的方法具有總體上的相符合性,但所獲得的值比ISO標準方法要大些(本文實驗要大10%左右),最大可能大20%左右,體現出其濾波效應更小。其子母峰形式的波形特性細節,展現出鋁墊層(實際上是一個圓柱形活動墊塊)在超高沖擊過程中的微小傾斜及不均勻造成的效果,波形更加接近真實曲線。在不同傳感器上所作的校準實驗結果,證明了本文所述方法的切實可行性,可以用于沖擊加速度計量校準中的微分過程處理。并且,與ISO推薦的方法相比,本文所述方法可以更好地實現自適應濾波和校準,并能給出更加明確的不確定度評定結果。

6 結 論

綜上所述可見,本文主要討論了沖擊加速度計量校準中的微分過程及濾波的實現問題,針對ISO相關標準中推薦方法存在的問題,提出了一種基于階躍上升時間為關聯參數的局域平均方式,同時實現數據濾波和數值微分,具有過程簡潔,不確定度明確的優點,并且具有局部尋優的條件判據,同時,避免了原方法中倒果為因的邏輯問題、以及需要判定是否峰值符合高斯形狀或矩形脈沖形狀的問題。對于具有子母峰形狀的超高沖擊波形依然適用。理論與實踐均證明了所述方法的有效性與可行性。

[ 1 ] Methods for the calibration of vibration and shock transducers-Part 13: primary shock calibration using laser interferometry: ISO 16063-13: 2001 [S]. [S.l.]: International Organization for Standardization, 2001.

[ 2 ] Methods for the calibration of vibration and shock transducers-Part 22: shock calibration by comparison to a reference transducer: ISO 16063-22: 2005/DAM 1 (DS) [S]. [S.l.]: [s.n.], 2005.

[ 3 ] 振動與沖擊傳感器校準方法—第13部分:激光干涉法沖擊絕對校準: GB/T 20485.13—2007 [S].北京:中國標準出版社,2008.

[ 4 ] 沖擊力法沖擊加速度校準裝置: JJG 791—1992 [S].北京:中國計量出版社,1993.

[ 5 ] 沖擊加速度計(絕對法)校準規范: JJF 1153—2006 [S].北京:中國計量出版社,2007.

[ 6 ] 壓電加速度計檢定規程: JJG 233—2008 [S].北京:中國計量出版社,2008.

[ 7 ] 梁志國,李新良,連大鴻.激光干涉法一次沖擊加速度校準[J].電子測量與儀器學報,2006,20(1): 68-72.

LIANG Zhiguo, LI Xinliang, LIAN Dahong. Primary shock calibration using laser interferometry[J]. Journal of Electronic Measurement and Instrument, 2006, 20(10): 68-72.

[ 8 ] 梁志國,孫璟宇.調頻信號的數字化解調[J].測試技術學報,2005,19(2): 190-194.

LIANG Zhiguo, SUN Jingyu. A precise digitizing demodulation method of the FM signal[J]. Journal of Test and Measurement Technology, 2005, 19(2): 190-194.

[ 9 ] 梁志國,李新良,朱振宇.一種基于殘周期正弦擬合的沖擊峰值計算方法[J].振動與沖擊,2015,34(1): 49-52.

LIANG Zhiguo, LI Xinliang, ZHU Zhenyu. Novel calculation method for peak of impulse waveforms based on partial sinusoidal curve-fitting[J]. Journal of Vibration and Shock, 2015, 34(1): 49-52.

Adiscussionondifferentialprocessingandpeakfilteringinimpulseaccelerationcalibration

LIANG Zhiguo

(National Key Laboratory of Science and Technology on Metrology amp; Calibration, Changcheng Institute of Metrology and Measurement, Beijing 100095, China)

Aiming at the differential processing in impulse acceleration by using laser interference methods, some problems of pre-filtering were discussed. Some experiments were performed. And a novel method of pre-filtering and differential processing was introduced. It is based on the rise time of impulse velocity waveform. Through local average and differential calculation, the acceleration waveform from the velocity waveforms could be obtained. The discomfort of determining filtering parameter based on acceleration impulse width in the ISO standard were avoided, the range of filtering-parameter were studied by experiment. Compared with the ISO method, the proposed method can get a higher peak value. It may be 10%~20% higher. By the experiment in different transducers, the validity and feasibility were proved, and it must be a novel processing method of impulse acceleration signals by using laser interference methods.

metrology; impulse acceleration; differential; calibration

2016-06-17 修改稿收到日期: 2016-08-19

梁志國 男,博士,研究員,1962年生

TB973

A

10.13465/j.cnki.jvs.2017.22.012

主站蜘蛛池模板: 国产精品亚洲片在线va| 一级香蕉视频在线观看| 国产精品刺激对白在线| 国产精品视频3p| 亚洲天堂网在线观看视频| 欧美精品伊人久久| 亚洲第七页| 久久99国产精品成人欧美| 国产一级在线观看www色| 伊人激情久久综合中文字幕| 无码日韩视频| 亚洲婷婷在线视频| 日本欧美视频在线观看| 91欧美在线| 亚洲一区无码在线| 国产成人一区二区| 亚洲黄色成人| 亚洲国产成人久久精品软件| 欧美人与牲动交a欧美精品| 久久久久亚洲精品成人网 | 亚洲国产亚综合在线区| 精品国产电影久久九九| 亚洲国产欧美目韩成人综合| 国产特级毛片| 99在线视频免费观看| AV熟女乱| 91久草视频| 国产成人免费视频精品一区二区 | 国产a网站| 亚洲一本大道在线| 好吊色国产欧美日韩免费观看| 国产手机在线观看| 免费一级成人毛片| 国产精品亚洲а∨天堂免下载| 日韩毛片免费观看| 国产va免费精品观看| 美女无遮挡拍拍拍免费视频| 97se亚洲| 国产精品白浆无码流出在线看| 亚洲精品视频免费看| 国产第四页| 久久黄色视频影| 高清视频一区| 国产精品区网红主播在线观看| 国产乱子伦视频在线播放| 免费国产高清视频| 亚洲国产综合精品一区| 日本亚洲欧美在线| 狠狠躁天天躁夜夜躁婷婷| 五月婷婷伊人网| 中文无码伦av中文字幕| 在线观看亚洲天堂| 免费观看男人免费桶女人视频| 美女无遮挡免费网站| 国产AV无码专区亚洲精品网站| 中文无码毛片又爽又刺激| 人禽伦免费交视频网页播放| 亚洲精品男人天堂| 欧美色99| 成人午夜精品一级毛片| 欧美成人二区| 一本大道香蕉高清久久| 国产精品亚洲日韩AⅤ在线观看| 区国产精品搜索视频| 激情亚洲天堂| 国产成在线观看免费视频| 国产精品无码AⅤ在线观看播放| 亚洲欧美综合另类图片小说区| 国产成人精品在线| 久久精品嫩草研究院| 亚洲精品无码成人片在线观看| 久久一本精品久久久ー99| 国产在线高清一级毛片| 91 九色视频丝袜| 国内嫩模私拍精品视频| 99精品视频播放| 国产成人综合在线观看| 国产一区二区三区在线观看免费| 午夜性爽视频男人的天堂| 亚洲一道AV无码午夜福利| 在线观看国产黄色| 久久女人网|