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

基于聯邦擴展卡爾曼濾波的結構損傷識別方法

2017-11-30 06:56:49王路丹宋固全徐昌宏
振動與沖擊 2017年21期
關鍵詞:卡爾曼濾波故障信號

張 純, 王路丹, 宋固全, 徐昌宏, 廖 群

(南昌大學 建筑工程學院, 南昌 330031)

基于聯邦擴展卡爾曼濾波的結構損傷識別方法

張 純, 王路丹, 宋固全, 徐昌宏, 廖 群

(南昌大學 建筑工程學院, 南昌 330031)

當結構損傷與傳感器故障同時存在時,兩者的相互影響會導致損傷識別結果的劣化;為此,提出了一種基于聯邦擴展卡爾曼濾波的結構損傷識別算法。利用分散化濾波計算量小、濾波精度高的優點,聯邦擴展卡爾曼濾波方法能根據正常的傳感器信號準確識別結構的損傷位置與程度,具有良好的魯棒性;同時利用聯邦濾波容錯性好的特點,能實現對故障信號的自動檢測和剔除,并將剩余的正常子系統進行組合,以繼續提供準確的損傷識別結果。梁式結構的數值算例及實驗分析驗證了該算法的有效性及對故障傳感器信號的檢測隔離能力。

聯邦擴展卡爾曼濾波; 損傷識別; 傳感器故障; 殘差檢驗法; 分散化濾波

在結構健康監測領域中,基于結構振動信號的系統辨識和損傷識別技術是當前國內外研究的熱點。目前,時域內基于模型的損傷識別方法[1-3]主要包括最小二乘估計、擴展卡爾曼濾波、序貫非線性最小二乘方法等。其中,擴展卡爾曼濾波作為一種可用于非線性系統結構參數識別的實時遞推算法,近年來被逐步應用于工程結構的參數反演與損傷識別[4-5],可實現測量信息不完備[6]、系統參數時變[7]等情況下的結構損傷狀態反演。然而,這些研究工作均是基于傳感器不存在故障、測量數據完全正確(或含白噪聲干擾)的前提下進行,相應的結構損傷識別算法不具有傳感器故障信號診斷排除能力。

在實際結構健康監測系統的長期使用過程中,傳感器常出現信號漂移、性能蛻化、故障或失效的問題。為保證傳感器能準確獲取測量信息,研究者提出了眾多傳感器故障檢測方法[8],但絕大部分基于數學模型的診斷方法是以健康的結構模型為基礎的。當結構損傷和傳感器故障同時出現時,兩者之間的相互作用會嚴重影響各自的診斷效果。特別是在反問題計算不適定性的影響下,基于故障或不夠準確測量信號的損傷識別結果會迅速劣化,失去診斷能力。因此,為保證結構健康監測系統的工作性能,迫切需要發展具有自主故障信號診斷和排除能力的損傷識別算法。

隨著多傳感器數據融合算法的發展,由Carlson[9]提出的聯邦卡爾曼濾波器因設計方案靈活,計算量小、容錯性好的優點在組合導航領域得到了廣泛應用。作為一種分散式濾波方法,聯邦濾波算法不僅具有非常高的識別精度,而且能夠對各測量子系統進行故障檢測,通過分離故障子系統或調整信息分配因子來抑制故障信號,從而保證系統狀態的正確跟蹤[10]。

本文結合聯邦濾波方法的優點,提出了一種適用于結構損傷識別的聯邦擴展卡爾曼濾波方法(Federal Extended Kalman Filter method, FEKF),不僅具有良好的計算效率、識別精度和魯棒性,還能及時判別和隔離傳感器故障信號,保證結構損傷識別結果的精度和穩定性。梁式結構的損傷識別數值算例及實驗分析驗證了本文方法的有效性以及對故障信號隔離檢測能力。

1 基本理論

1.1結構動力學方程

經過有限元離散后,線性動力系統的運動方程可以寫為:

(1)

(2)

由于實際結構的自由振動響應可較為方便地通過脈沖激勵或自然激勵信號得到,因此,文中主要考慮結構的自由振動情況。于是,由式(2)可解得模態坐標的理論解:

pn(t)=e-ξnωntpn0cos(ωdnt)+

(3)

1.2聯邦擴展卡爾曼濾波算法

聯邦卡爾曼濾波是一種以最小方差為準則的遞推估計算法,在對結構參數進行辨識時,可視為一種遞推型的模型修正方法。將聯邦卡爾曼濾波應用于結構損傷識別時,需要把系統的狀態向量進行擴展,通過引入損傷參數來跟蹤單元彈性模量的變化規律,從而識別結構損傷。因此,將聯邦擴展卡爾曼濾波的系統狀態向量定義為:

(4)

根據聯邦濾波器理論框架[10],它是一種兩層結構的分散化濾波器,由若干個子濾波器和一個主濾波器組成。使用時,采用信息分配原理把整個系統中的信息分配到每個局部子濾波器中,利用“方差上界”技術使得每個子濾波器獨立工作,然后根據協方差陣把子濾波器的信息進行最優融合估計,得到系統的全局估計值。其結構圖如圖1所示。

圖1 聯邦卡爾曼濾波結構圖

為解決結構損傷識別問題,本文方法將所有測量傳感器進行分組,構成不同的子系統,并采用獨立的子濾波器進行參數識別。每個子系統均可形成如下形式的損傷參數估計問題:

(5)

zi(xi,t)=h(θ,xi,t)+vi

(6)

式中:i代表第i個子系統;相應的觀測噪聲vi為零均值、協方差矩陣為Ri的高斯白噪聲;zi是子系統的觀測向量。函數h(θ,xi,t)為子系統測點處的結構振動響應的估計值,根據式(3),其理論解可表示為

(7)

式中:φn(x)為第n階的模態向量,x對應觀測點的空間坐標。在各子濾波器中,由式(5)、(6)組成的非線性參數識別問題按照傳統擴展卡爾曼方法進行遞推;再將其識別結果送入主濾波器中進行融合計算,以得到狀態向量的最優估計,同時為各子濾波器提供更新后的狀態估計值。

根據以上分析,本文建立了基于聯邦擴展卡爾曼濾波的結構損傷識別遞推算法,其具體流程如下:

(1) 信息分配過程:將主濾波器第k時刻的全局狀態估計θg(k)及估計誤差協方差矩陣Pg(k)反饋到子濾波器中,以重置子濾波器估計值。具體的信息分配原則為

(8)

θi(k)=θg(k)

(9)

(2) 信息的時間更新:時間更新過程在各子濾波器中獨立進行,其中第i個子濾波器的時間更新算法為:

(10)

(11)

(3) 信息的測量更新:測量更新過程也只在各子濾波器中獨立進行

(12)

(13)

Ri(k+1))-1

(14)

其中因反問題計算的不適定性,式(14)中的矩陣求逆采用了正則化方法[11]。Ki(k+1)為增益矩陣,I為單位矩陣,Hi(k+1)為非線性函數h(θ,xi,t)的雅克比矩陣:

Hi(k+1)=?h(θi(k+1),xi,tk+1)/?θi(k+1)

(15)

當結構質量不發生變化時,(15)式計算中所需的固有頻率和模態關于損傷參數的靈敏度?ωn/?αi、?φn/?αi可用基準有限元模型計算[12]。

(4) 信息融合:將各個子濾波器k+1時刻通過時間更新以及測量更新得到的狀態向量送入主濾波器中進行融合處理,得到k+1時刻的最優全局狀態估計及估計誤差協方差矩陣:

(16)

(17)

方程組式(8)~(17)構成了聯邦擴展卡爾曼濾波遞推過程中的一步,重復以上4個步驟即可完成系統狀態向量的最優估計,識別出結構損傷。

1.3基于聯邦卡爾曼濾波的傳感器故障診斷

隨著傳感器工作性能的蛻化,測量信號逐步偏離結構的真實動力響應。此時的信號誤差與卡爾曼濾波模型中設定的高斯白噪聲有明顯差異,將導致顯著的損傷識別誤差。由于聯邦擴展卡爾曼濾波是一種分散式濾波,有良好的容錯性,故障傳感器的影響首先會集中體現在相應的子濾波器中。因此,結合經典的殘差χ2檢驗法[10],本文方法可以在識別結構損傷的同時,判斷測量信號質量,有效診斷、排除故障傳感器,保證識別結果的穩定性。

由各子濾波器第k時刻的一步狀態預測值和觀測值構造出每個子濾波器的新息:

(18)

如果k時刻以前整個系統沒有出現故障,各子濾波器新息均為零均值高斯白噪聲;當某個傳感器發生故障時,相應子濾波器新息將不再滿足零均值白噪聲條件。因此,根據新息的均值變化即可進行故障檢測。

對第i個子濾波器的新息作二元假定:

(1)H0:無故障,

E(ri(k))=0

(2)H1:有故障,

E(ri(k))=μi

E[(ri(k)-μi)(ri(k)-μi)T]=Si(k)

其中Si(k)為ri(k)的方差,根據極大似然比檢驗原理,可以構建如下故障檢測函數:

(19)

(1) 若λz(k)>TD,判定有故障

(2) 若λz(k)≤TD,判定無故障

其中TD為預先設定的門限,可根據子系統測點數和預先設定的誤警率,通過查詢χ2分布表得到。于是,利用殘差χ2檢驗法可以檢測出各個子系統是否存在故障,然后將各傳感器進行不同的子系統組合以確定故障傳感器;進而隔離故障信號,避免故障信號對識別結果的干擾。其具體流程如圖2所示。

2 數值算例分析

數值算例采用矩形截面的簡支梁結構,梁長10 m,截面尺寸為0.25 m×0.5 m,材料彈性模量和密度分別為210 GPa、8 200 kg/m3。簡支梁結構被等分為10個梁單元,在第2~10節點上設置傳感器,以測量節點豎向位移。同時,根據傳感器數量與布置,本文采用了3個子系統,每個子系統包含3個測點,具體的單元劃分及測點分配情況如圖3所示。

施加任意的初始條件后,將各測點的自由振動信號作為觀測值(采樣頻率設為1 000 Hz),同時按如下公式加入高斯白噪聲模擬觀測誤差:

zwgn(t)=z(t)+v(t)

(20)

圖2 故障傳感器診斷與隔離

圖3 簡支梁模型與各子系統測點信號分配

在本文算例中,結構自由振動響應以低階模態成分為主;理論上在狀態變量(4)中保留更多的模態坐標有利于提高計算精度,但也會增加計算負擔。因此,本文算法采用前3階模態坐標來構建狀態變量。

2.1單一損傷識別結果

考慮單元5損傷參數為0.3,其它單元無損傷的工況。對于結構響應數據中無噪聲、信噪比分別為30 dB、25 dB、20 dB和15 dB的工況,簡支梁結構的損傷識別結果如圖4所示。盡管輸入信號的噪聲不斷增大,本文方法始終保持很高的損傷識別精度。即使信噪比達到15 dB時,損傷單元彈性模量識別結果的最大相對誤差為1.35%。

圖5給出了信噪比15 dB時單元損傷參數的識別曲線,初始階段識別結果波動較大,隨后損傷參數快速收斂到實際損傷值,體現出良好的收斂性。

圖6給出了測點2的真實信號、加噪信號、濾波信號的結果對比圖(為清晰起見,0.25~0.40 s局部時間段的信號如圖6(b)所示)。圖中,本文算法估計出的測點位移響應信號與真實信號幾乎完全重合,濾波效果良好。表1列出了單損工況下狀態變量(4)中模態坐標初始值、模態坐標對時間導數的初始值、模態阻尼比系數的識別結果。

圖4 不同噪聲水平下結構單一損傷的識別結果

圖5 信噪比15 dB時損傷單元參數收斂曲線

參數無噪30dB25dB20dB15dB真實值識別值識別值識別值識別值識別值p101.51.501.501.501.491.51p200.80.800.800.800.810.79p30-0.5-0.50-0.50-0.50-0.50-0.5p101514.9414.9514.9814.8813.9p204545.0244.8143.3945.7645.4p30110108.58112.65113.62109.01101.07ξ10.0250.02500.02500.02500.02460.0251ξ20.0250.02500.02510.02510.02570.0257ξ30.0250.02470.02490.02490.02530.0247

2.2多處損傷結果識別

設第2、5、8號單元的損傷參數分別為0.4、0.3、0.2,其它單元無損傷。從圖7結果可以看出,對于結構多處損傷工況,FEKF方法依然能準確識別損傷的位置與程度。所設置的3處損傷單元彈性模量識別結果的最大相對誤差為4.43%,并且隨著噪聲級別增大,損傷識別結果始終保持很好的精度,體現出算法良好的魯棒性。從圖8及表2給出的結果,可以看出與單一損傷工況的結果類似,即多損工況下測點信號也能跟蹤準確,損傷參數也能較快收斂到真值,狀態變量的所有元素均能正確識別。

(a)

(b)

圖7 不同噪聲水平下結構多處損傷的識別結果

利用表1、表2結果,進一步分析不同噪聲級別的影響,可以發現狀態向量(4)中各元素識別結果的精度與其對結構響應的貢獻密切相關。與損傷參數、模態坐標初始值相比,模態坐標時間導數的初始值對動力響應的貢獻會由于阻尼的作用迅速衰減;因此,隨著噪聲級別的提高,模態坐標時間導數的識別結果會首先劣化;如果噪聲級別繼續增大,損傷參數、模態坐標初始值的識別誤差也會不斷增加。總體上,算法在抑制噪聲干擾方面具有良好的性能。

圖8 信噪比15 dB時損傷單元參數收斂曲線

參數無噪30dB25dB20dB15dB真實值識別值識別值識別值識別值識別值p101.51.501.501.501.501.50p200.80.800.800.790.810.81p30-0.5-0.50-0.49-0.48-0.51-0.5p101514.9614.8015.1714.9015.8p204544.9345.6547.6444.8252.6p30110116.84115.94112.77111.30103.40ξ10.0250.02500.02510.02500.02490.0251ξ20.0250.02500.02480.02490.02510.0257ξ30.0250.02520.02490.02480.02640.0248

2.3傳感器發生故障的結構損傷識別結果

上述計算結果表明,聯邦擴展卡爾曼濾波算法能夠準確識別結構的損傷位置與程度,且具有良好的魯棒性和收斂性。為進一步驗證本文方法的故障檢測與隔離能力,采用2.2節簡支梁結構多損傷工況(信噪比為20 dB)進行分析。受篇幅限制,本文僅考慮突變故障。假設測點5處的傳感器發生突變故障(即故障信號與真實信號有一個突變的差值),導致測點5位移信號從0.05 s開始出現一個大小為0.01 m的突變故障,如圖9所示。此時,如果不對故障信號進行檢測和排除,圖10所示結果表明,損傷識別結果會因測點5的突變故障信號干擾而出現明顯的劣化,從而導致損傷誤判。

利用聯邦擴展卡爾曼濾波方法檢測故障信號時,由于每個子系統使用3個測點信息,維數為3,誤警率設為0.5%,查χ2分布表可得門限值TD=12.838。根據圖2所示的流程,采用殘差χ2檢驗法進行檢驗:

(1) 由圖11(a)所示結果可知,第一個子濾波器的故障檢測函數值遠大于門限值,表明第一個子系統中存在故障信號;而在圖11(b)中第二個子濾波器的故障檢測函數值均小于門限值,表明信號正常;第三個子濾波器的現象與第二個子濾波器類似,未單獨列出。由此,可判定故障信號來自于第一個子系統,且第二、三子系統中的傳感器工作正常,這與算例中的設定是一致的。

(2) 在對各子系統是否存在故障傳感器進行判斷后,可交換各子系統中的傳感器,并根據子系統故障檢測函數圖的變化,進一步確認故障傳感器編號。即將第一個子系統(含故障傳感器)的測點2、5、8依次與第二個子系統的測點3進行交換。再利用分別重組的信號重新檢驗,計算結果顯示,當交換測點5與測點3時,第二個子系統出現故障信號,而第一個子系統信號正常,故認為測點5處的傳感器為故障傳感器。

圖9 測點5的故障信號(20 dB)

圖10 含故障信號的聯邦擴展卡爾曼濾波識別結果

在完成故障信號診斷后,即可通過對故障信號的隔離,有效保證損傷識別結果的正確性。圖12為剔除故障信號后的聯邦擴展卡爾曼濾波損傷識別結果。對比圖10和圖12可知,聯邦擴展卡爾曼濾波方法能結合殘差檢驗法實現對故障傳感器信號的檢測,通過故障信號隔離以及觀測信號的冗余性,保證了結構損傷的正確性與穩定性,有效避免了故障信號導致的損傷誤判。

(a) 第一個子濾波器

(b) 第二個子濾波器

圖12 排除傳感器故障后的FEKF損傷識別結果

3 懸臂梁實驗分析

實驗懸臂梁采用長75 cm,厚0.4 cm,寬4.0 cm的矩形截面鋼梁,材料的彈性模量及密度分別為167.64 GPa、7149.2 kg/m3。相應的單元劃分、加速度傳感器布置和測點分配情況如圖13(a)所示。

結構損傷通過在懸臂梁5號單元開孔(Ф=2 cm)模擬。在懸臂梁有限元模型中調整5號損傷單元的彈性模量折減量,使得有限元模型與實際結構固有頻率盡量一致,從而可確定實驗懸臂梁在5號單元的損傷參數為0.2。

(a) 懸臂梁有限元模型與各子系統測點信號分配

(b) 實驗懸臂梁

結構各測點的自由振動信號由力錘激勵后測得,采樣頻率為1 000 Hz。為模擬傳感器故障,人為在測點5處振動傳感器的測量信號中加入了0.002 m的突變(從0.1 s開始)。

利用本文提出的FEKF算法進行損傷識別。當存在傳感器故障時,圖14清晰地指示出故障傳感器所在的子濾波器,進一步通過傳感器交換檢測確定故障傳感器并進行信號隔離后,本文方法的損傷識別結果如圖15所示。為進行對比,圖中也列出了無故障傳感器時的損傷識別結果。由圖示結果可見,無論是否存在故障傳感器,本文提出的聯邦擴展卡爾曼濾波算法均能通過自主的故障診斷與信號處理,有效識別損傷,具有良好的識別精度與魯棒性。

(a) 第一個子濾波器

(b) 第二個子濾波器

圖15 懸臂梁損傷識別結果

4 結 論

本文提出了一種基于聯邦擴展卡爾曼濾波的結構損傷識別方法,能利用結構自由振動響應準確識別結構的損傷位置與程度。根據梁式結構的數值模擬結果及實驗分析可以得到以下結論:

(1) 在不同的損傷工況下,聯邦擴展卡爾曼濾波方法均能準確地識別出結構的損傷位置和損傷程度,且具有很好的魯棒性。識別精度受噪聲影響較小,即使信噪比增大到15 dB,也能保證較高的識別精度。

(2) 當結構損傷和傳感器故障同時存在時,聯邦擴展卡爾曼濾波算法能結合殘差χ2檢驗法準確識別出故障信號來源,并通過對故障信號的自動剔除,有效保證損傷識別結果的正確性與穩定性。

[1] 雷鷹, 江永強. 輸入輸出信息有限觀測下的結構損傷診斷[J]. 振動、測試與診斷, 2012, 32(5): 736-740.

LEI Ying, JIANG Yongqiang. Structural damage detection technique with limited input and output measurement signals [J]. Journal of Vibration, Measurement & Diagnosis, 2012, 32(5): 736-740.

[2] CALABRESE A, SERINO G, STRANO S, et al. An extended Kalman Filter procedure for damage detection of base-isolated structures [P]. 2014 IEEE Workshop on Environmental Energy and Structural Monitoring Systems (EESMS), Naples, Italy. IEEE, USA, 2014.

[3] YANG J N, HUANG H. Sequential non-linear least-square estimation for damage identification of structures with unknown inputs and unknown outputs[J]. International Journal of Non-linear Mechanics, 2007, 42(5): 789-801.

[4] 樊素英, 李忠獻. 基于廣義卡爾曼濾波的橋梁結構物理參數識別[J]. 計算力學學報, 2007, 24(4): 472-476.

FAN Suying, LI Zhongxian. Identification on physical parameters of bridge structures based on extended Kalman filter[J]. Chinese Journal of Computational Mechanics, 2007, 24(4): 472-476.

[5] 趙博宇, 丁勇, 吳斌. 基于擴展卡爾曼估計算法的地震模擬振動臺模型識別 [J]. 振動與沖擊, 2014, 33(12): 145-150.

ZHAO Boyu,DING Yong,WU Bin. Identification of shaking table model for seismic simulation based on an extended Kalman estimator[J]. Journal of Vibration and Shock, 2014, 33(12): 145-150.

[6] 謝獻忠, 易偉建. 結構物理參數時域識別的子結構方法研究[J]. 工程力學, 2005, 22(5): 94-98.

XIE Xianzhong, YI Weijian. A substructure method for parameter estimation in time domain[J]. Engineering Mechanics, 2005, 22(5): 94-98.

[7] 穆騰飛, 周麗. 輸入未知條件下基于自適應廣義卡爾曼濾波的結構損傷識別[J]. 振動工程學報, 2014, 27(6): 827-834.

MU Tengfei, ZHOU Li. Structural damage identification using adaptive extended Kalman filter with unknown inputs[J]. Journal of Vibration Engineering, 2014, 27(6): 827-834.

[8] 周東華, 胡艷艷. 動態系統的故障診斷技術[J]. 自動化學報, 2009, 35(6):748-758.

ZHOU Donghua, HU Yanyan. Fault diagnosis techniques for dynamic systems[J]. Acta Automatica Sinica, 2009, 35(6):748-758.

[9] CARLSON N A, BERARDUCCI M P. Federated kalman filter simulation result[J]. Journal of the Institute of Navigation, 1994, 41(3): 297-321.

[10] 卞鴻巍, 李安, 覃方君. 現代信息融合技術在組合導航中的應用[M]. 北京: 國防工業出版社, 2010.

[11] 戴霖. 基于擴展卡爾曼濾波的結構損傷識別方法[D]. 南昌: 南昌大學, 2014.

[12] WEBER B, PAULTRE P, PROULX J. Consistent regularization of nonlinear model updating for damage identification[J]. Mechanical Systems and Signal Processing, 2009, 23(6): 1965-1985.

Structuraldamageidentificationbasedonthefederalextendedkalmanfilter

ZHANGChun,WANGLudan,SONGGuquan,XUChanghong,LIAOQun

(School of Civil Engineering and Architecture, Nanchang University, Nanchang 330031, China)

The coexistence of structure damages and sensor faults will deteriorate identified results evidently, so an algorithm for the identification of structure damages based on the Federated Extended Kalman Filter method (FEKF) was proposed by using free vibration signals. The presented method can identify the location and extent of damages accurately, and shows good robustness when the sensors work normally. Combined with the residual chi-square test, the FEKF also can eliminate the effects of fault sensors by the automatic detection and removal of the fault sensor signals. Numerical simulations and experiments show that the FEKF can ensure the accuracy and stability of the damage identification results and detect fault signals effectively.

federal extended Kalman filter; damage identification; sensor fault; residual chi-square test; decentralized filtering

TU311

A

10.13465/j.cnki.jvs.2017.21.027

國家自然科學基金(51268045;51469016);教育部高等學校博士點基金(20123601120011);水沙科學與水利水電工程國家重點實驗室開放研究基金(sklhse-2014-C-03);廣東省水利科技創新基金(2014-08)

2015-12-23 修改稿收到日期:2016-07-09

張 純 男,博士,教 授,1976年10月

宋固全 男,博士,教 授,1964年2月

猜你喜歡
卡爾曼濾波故障信號
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
故障一點通
基于FPGA的多功能信號發生器的設計
電子制作(2018年11期)2018-08-04 03:25:42
基于遞推更新卡爾曼濾波的磁偶極子目標跟蹤
奔馳R320車ABS、ESP故障燈異常點亮
基于LabVIEW的力加載信號采集與PID控制
基于模糊卡爾曼濾波算法的動力電池SOC估計
電源技術(2016年9期)2016-02-27 09:05:39
基于擴展卡爾曼濾波的PMSM無位置傳感器控制
電源技術(2015年1期)2015-08-22 11:16:28
故障一點通
主站蜘蛛池模板: 亚洲国产成人综合精品2020 | 国产成人综合网在线观看| av天堂最新版在线| 成人精品午夜福利在线播放| 国产精品自在自线免费观看| 97久久精品人人做人人爽| 亚洲国产精品国自产拍A| 麻豆国产精品一二三在线观看| 欧洲成人免费视频| 国产jizzjizz视频| 国产91麻豆视频| 男人的天堂久久精品激情| 精品小视频在线观看| 日韩久草视频| 青青草国产精品久久久久| 99精品热视频这里只有精品7| 91精品福利自产拍在线观看| 亚洲综合久久成人AV| 国产免费怡红院视频| 一级毛片免费观看久| 亚洲成a人片77777在线播放| 青青极品在线| 亚洲综合色婷婷| 日韩精品少妇无码受不了| 欧美激情伊人| 亚洲美女久久| 亚洲三级色| 久久大香伊蕉在人线观看热2| 日本精品αv中文字幕| 亚洲永久视频| 亚洲欧美在线综合图区| 免费国产高清精品一区在线| 伊伊人成亚洲综合人网7777| 青青青国产视频手机| 欧美特级AAAAAA视频免费观看| 无码福利视频| 国产精品微拍| 亚洲最大福利网站| 久久一色本道亚洲| 日本精品影院| 麻豆精品在线视频| 熟女成人国产精品视频| 波多野衣结在线精品二区| 亚洲国产日韩在线观看| 亚洲AV无码一区二区三区牲色| 成人精品区| 亚洲人成电影在线播放| 免费看a级毛片| 亚洲一区二区视频在线观看| 丰满少妇αⅴ无码区| 国产精品高清国产三级囯产AV| 一本大道无码高清| 最新痴汉在线无码AV| jizz国产视频| 亚洲小视频网站| 免费中文字幕在在线不卡| 精品乱码久久久久久久| 欧美亚洲日韩中文| 日本午夜精品一本在线观看 | 亚洲大尺码专区影院| 尤物特级无码毛片免费| 亚洲精品卡2卡3卡4卡5卡区| 波多野结衣在线一区二区| 最新加勒比隔壁人妻| 欧美h在线观看| 国产在线观看一区精品| 国产精品一区二区不卡的视频| 波多野结衣亚洲一区| 漂亮人妻被中出中文字幕久久| 亚洲无码高清视频在线观看| 亚洲精品视频在线观看视频| 国产特一级毛片| 全部毛片免费看| 久青草网站| 国产在线精品99一区不卡| 欧美成人怡春院在线激情| 久久中文字幕2021精品| 亚洲乱码在线视频| a毛片在线| 久热99这里只有精品视频6| 亚洲男人的天堂久久香蕉| v天堂中文在线|