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

基于信息熵與譜有限元法的傳感器優化布置

2016-07-26 08:15:00張加培朱宏平
振動與沖擊 2016年2期

張加培, 尹 濤, 朱宏平, 丁 蘭

(1.武漢大學 土木建筑工程學院,武漢 430072;2.華中科技大學 土木工程與力學學院,武漢 430074)

?

基于信息熵與譜有限元法的傳感器優化布置

張加培1,2, 尹濤1, 朱宏平2, 丁蘭2

(1.武漢大學 土木建筑工程學院,武漢430072;2.華中科技大學 土木工程與力學學院,武漢430074)

摘要:為從測量數據中獲得盡可能多信息,減少待識別模型參數的不確定性,提出面向結構模型參數識別的傳感器優化布置方法。為避免用靜態形函數傳統有限元方法建模對結構動力特性及傳感器優化布置影響,采用高精確動力學法即譜有限元法對結構進行動力學建模。以結構模型參數識別結果的不確定性最小作為傳感器優化布置準則,該不確定性程度通過信息熵標量指標量化,用貝葉斯統計系統識別法進行識別。采用整數編碼遺傳算法在所有可能的傳感器配置組合中極小化信息熵指標,獲得給定數目的傳感器最優布置位置。通過彈性地基帶彈性接頭的周期管梁模型數值仿真及模型試驗驗證所提方法。

關鍵詞:結構健康監測;傳感器優化布置;信息熵;譜有限元法;遺傳算法;周期結構

隨土木工程的迅猛發展,結構形式日趨復雜,結構健康監測尤其基于振動測量數據的結構健康監測成為重要研究方向[1-2]。該監測方法基于觀測的振動數據,而振動數據由傳感器采集獲得,即傳感器個數及其位置決定數據的數量及質量,會直接影響結構健康監測成功與否。顯然,測量的信息量越大結構健康監測會越準確,但受測試條件、實際環境及成本所限,不可能無限布置傳感器。因此對給定數目傳感器,如何布置才能獲得盡可能多的信息成為傳感器優化布置解決的問題[3]。

目前傳感器優化布置方法均基于傳統有限元法,但其所用與頻率無關的靜態形函數求解動力問題尤其頻率較高問題時不夠精確,會影響傳感器優化布置結果[4]。而譜有限元法基于快速傅里葉變換,采用精確波動方程解作為形函數,是傳統有限元方法在頻域內的表達。理論上,譜有限元法為高精度結構動力學分析方法,用較少數量自由度獲得頻域的精確解[5]。故本文用譜有限元法對模型進行動力響應分析,并以此為基礎進行傳感器優化布置研究。

解決傳感器最優布置問題,主要有優化布置準則及優化方法兩種。目前大多所用優化準則有識別誤差最小準則、模態應變能準則、模型縮減法準則、插值擬合準則及可控度可觀察度準則等[6],但對結構健康監測意義不大。Papadimitriou等[7]提出以結構模型參數識別結果不確定性最小為優化布置準則的傳感器優化布置方法,實質為通過貝葉斯統計系統識別法計算模型參數的不確定性[8],用信息熵度量參數估計不確定性大小,信息熵最小時即為參數識別不確定性最小,傳感器位置最優。信息熵方法不僅可用于比較不同數目傳感器優化布置方案的優劣,且適用于線性、非線性模型識別。尹濤[9]將文獻[7]方法推廣到分布參數結構,基于連續坐標結構體系進行傳感器優化布置問題研究,實現傳感器與激勵器優化布置。

在傳感器所有可能配置組合中快速搜索滿足優化準則最優需優化算法實現。理論上該優化問題可用窮舉法得到答案。而當結構有大量自由度時仍采用窮舉法在所有可能組合中最小化信息熵,效率較低甚至不可行。遺傳算法則非常適用此離散優化,故本文基于MATLAB遺傳算法工具箱的整數編碼遺傳算法,針對傳感器優化布置問題對遺傳算法目標函數進行適當約束[10],解決傳感器最優布置問題。

圖1 基于信息熵及譜有限元法的傳感器優化布置算法流程Fig.1 Procedure of optimal sensor configuration by spectral finite element method and information entropy

本文推廣基于信息熵的傳感優化布置方法并與譜有限元法結合,形成基于譜有限元法的信息熵指標,采用遺傳算法極小化該信息熵指標,進行傳感器布置優化問題求解。算法流程見圖1。用實驗室建立的彈性地基帶彈性接頭的周期管梁模型數值仿真及實驗數據研究,對本文所提方法進行驗證。

1理論背景

1.1貝葉斯統計系統識別法

模型參數不確定性可由貝葉斯統計系統識別方法計算獲得。設不同模型參數a決定不同模型M,一旦a確定模型M隨之確定。對特定模型M(a),定義q(n;a)∈RNd為模型在所有Nd個自由度上tn=nΔt時刻輸出,其中Δt為采樣間隔。設在所有Nd個自由度上僅布置N0個傳感器,則傳感器位置可用向量δ∈RNd表示。當第i個自由度上有傳感器時向量元素δi=1,否則為0。因此共有N0個自由度可觀察到,系統在tn時刻N0個自由度上輸出向量y(n)∈RN0可表示為

y(n)=S0q(n;a)+S0e(n;a)

(1)

式中:e(n;a)為由模型誤差、測量噪聲產生的預測誤差;S0為選擇矩陣,每行只有1個等于1的元素,其余全為0。S0中1元素位置由傳感器位置向量δ決定。

據貝葉斯統計系統識別理論,待識別參數值及相應不確定性可通過測量動態數據D識別。當J(a)最小時結構模型參數a即為最優值,表示為

(2)

式中:‖·‖為2的范數;N為采樣時間點總數。

參數a的不確定性可通過建立概率模型用概率密度函數量化。設模型為線性且預測誤差不確定性滿足標準正態分布,則a的概率密度函數可由漸近逼近[11]表示為

(3)

式中:c為標準化常數;π(a)為a的先驗分布。

1.2信息熵

(4)

(5)

為計算方便,式(5)進一步簡化為

(6)

(7)

對兩種不同傳感器位置向量δ及δ0,H及H0分別為其信息熵,其中δ0可在傳感器數量、位置上均與δ不同。由式(4)獲得H-H0即信息熵的改變量為

(8)

因此關于傳感器位置向量δ及δ0兩種不同分布,參數不確定性比值可表示為

(9)

式中:Na為模型參數個數。

由式(9)看出,參數不確定性比值僅依賴于H-H0的改變量及a的數量,該比值可用于度量兩種不同傳感器配置的參數不確定性改變量。

1.3譜有限元法建模

用彈性地基帶彈性接頭的周期管梁模型研究傳感器優化布置問題。該模型圓管段用梁單元模擬,管間接頭僅考慮彎曲變形用轉動彈簧kθ模擬。利用彈簧將眾多管梁聯結形成彈性地基的周期性梁-彈簧模型,見圖2。

圖2 彈性地基周期管梁模型Fig.2 Periodic pipe-beam model on elastic foundation

考慮Winkler彈性地基,梁的自由彎曲振動方程[12]可表示為

(10)

式中:w(x,t)為橫向位移;EI為彎曲剛度;A為梁橫截面積;ρ為質量密度;Kf為彈性地基剛度系數。

橫向彎矩、剪力分別表示為

(11)

對式(10)兩邊進行快速傅里葉變換[13],頻域內梁的平衡方程為

W″″-kFW=0

(12)

式中:W(x,ω)為頻域內橫向位移;kF為彎曲波數,即

(13)

對長度為L的管梁段,平衡方程通解表達式為

W(x,ω)=a1e-ikFx+a2e-kFx+a3eikFx+a4ekFx(14)

梁單元譜節點位移表達式為

(15)

梁單元譜節點荷載表達式為

(16)

將節點譜位移向量與譜荷載向量寫成矩陣相乘形式,即

SB(ω)d=fc(ω)

(17)

式中:SB(ω)為梁的譜單元剛度矩陣,即

(18)

式中:

頻域內接頭平衡方程為

(20)

式中:kθ為接頭彎曲剛度。

接頭彈簧單元譜剛度矩陣可表示為

(21)

形成圓管梁與彈簧譜單元剛度陣后,采用與傳統有限元法相同步驟,在頻域內將所有單元譜單元剛度陣組裝成整體剛度陣,將頻域內荷載向量組裝成整體向量,求解方程組即可獲得頻域內節點位移,再通過逆傅里葉變換即可獲得各節點位移時間歷程。

2數值仿真

基于MATLAB遺傳算法工具箱中整數編碼遺傳算法求解本文傳感器優化布置問題。以整數表示傳感器布置,如在9個待測自由度上布置3個傳感器,變量[1,5,8]表示第1、5、8自由度位置布置傳感器。為避免出現兩個及以上數目傳感器布置在相同位置,本文在遺傳算法的目標函數中增設一不等式約束進行限制。此外,由于遺傳算法工具箱中優化函數默認尋求目標函數最小值,當det(Q)取最大值時即表示傳感器布置位置最優,因此用-det(Q)作為遺傳算法目標函數。遺傳算法主要控制參數中再生參數取2,交叉率取0.8,遷移率取0.2。

對彈性地基上通過環間接頭連接的周期性管梁模型進行數值仿真研究,模型見圖3。共100段梁,每段梁長均為0.5 m,橫截面為圓環,模型具體幾何、材料參數見表1。考慮本模型彎曲變形顯著的結構特性,接頭處僅考慮彎曲剛度。僅將管梁接頭彎曲剛度作為不確定結構參數進行識別,研究傳感器最優布置問題。該模型由99個接頭彈簧連接而成,共計99個待識別參數。用系數向量a修正接頭彈簧彎曲剛度,即ki=aikθ(i=1,2,…,99)。此外,該模型主要有豎向位移及轉角兩種變形,但考慮轉角測量較困難,僅考慮豎向加速度響應測量。由圖3可知,在所有接頭處共有99個豎向自由度可供布置傳感器。

圖3周期管梁模型
Fig.3 Periodic pipe-beam model

表1 幾何及材料參數

基于譜有限元法對模型進行動力學建模,考慮實驗時外荷載施加方便,僅在模型左端作用一個豎向沖擊荷載。經計算,本模型關心信號最高頻率為500 Hz,據采樣頻率定理,采樣頻率取1 280 Hz,采樣時長為3.2 s,即采樣點數N=4 096。

由式(9)知,s/s0可用于度量兩種不同傳感器配置的不確定性程度改變量,值越小表明參數識別不確性越小。用本文遺傳算法進行傳感器布置優化計算,獲得不同數目傳感器最優布置下信息熵s/s0值及具體位置,結果見表2。其中參考值s0對應所有99個自由度均布置傳感器情況。由表2看出,在給定傳感器數目下,經優化布置的傳感器較任意布置能更有效減小模型參數的不確定性,且隨優化布置傳感器數目增加信息熵隨之減小。因此,減小模型參數識別結果的不確定性及從測量數據中獲得盡可能多的信息,除單純增加傳感器數量外,需對傳感器進行優化布置,以最大限度節省測試成本。本文所提傳感器優化布置方法可實現此目的。

表2 信息熵s/s0及最優位置

3實驗研究

為進一步驗證本文傳感器優化方法的正確性,在實驗室搭建彈性地基的周期管梁實驗模型及測試系統,見圖4。該模型幾何及材料特性參數與數值仿真相同(表1)。受模型土箱尺寸限制,僅制作5段圓管(圖4(a)),由4個環間彈性接頭連接各段圓管。與數值仿真類似,僅將各接頭彎曲剛度作為結構待識別不確定性參數,即本模型共有4個待識別參數。因豎向位移自由度可測,即所有接頭共4個自由度供傳感器布置。

實驗所用儀器設備為,4個加速度傳感器(圖4中(b))分別置于4個彈性接頭附近,用于測量模型第1~4號自由度豎向加速度響應,由DH5922信號調理采集(圖4中(c)),再由裝DHDAS動態信號采集分析系統的筆記本電腦記錄(圖4中(d)),并將數據轉存為MATLAB格式用于后續分析計算。本實驗模型用沖擊力錘在結構左端頂部施加一豎向沖擊激勵,采樣頻率、時長及點數均與數值仿真保持一致。

圖4 實驗系統Fig.4 Experimental system

因僅有4個位置供加速度傳感器布置,傳感器配置組合較少,盡管采用窮舉法可輕易獲得所有可能組合的結果,但仍用遺傳算法對窮舉法結果對比驗證。分3種工況,分別對應傳感器數目1、2、3,所有可能組合的信息熵s/s0值見表3,其中,參考值s0對應所有4個自由度均布置傳感器情況,各測點自由度加速度時程響應見圖5。

表3 信息熵s/s0

圖5 加速度響應時程曲線Fig.5The measured acceleration time history

由表3看出,改變傳感器數目、位置時信息熵隨之變化。據本文傳感器優化布置準則,信息熵s/s0越小,參數識別不確定性越小,傳感器位置更優。據表3結果,傳感器數目為2時信息熵s/s0值均小于傳感器數目為1時;傳感器數目為3時信息熵s/s0值均小于傳感器數目為2時。對本模型而言,傳感器布置越多,所得信息亦越豐富,參數識別不確定性越小,識別結果越準確。對僅布置單個傳感器情況,s/s0最大值為3.27,對應傳感器安裝在4號自由度,即最差位置。s/s0最小值為2.2,對應傳感器置于2號自由度,即傳感器最優位置。對2個傳感器情況,最優位置為1、2號自由度,而對3個傳感器情況,最優位置為1、2、3號自由度。窮舉法與遺傳算法結果一致,證明本文改進遺傳算法的正確性。對本算例而言,傳感器數目為N+1的最優位置始終含傳感器數目為N的最優位置。

為進一步驗證本文傳感器優化布置結果,給出第2、4號自由度加速度響應幅值譜,見圖6。由圖6看出,2號自由度幅值譜出現的峰值數目、大小均較4號幅值譜明顯,表明2號自由度頻域成分更豐富。而對系統參數識別而言,表示由2號自由度獲得數據所含信息量更大,即2號較4號更適合布置傳感器。表3中2號自由度信息熵為2.2,4號自由度信息熵為3.27。由于信息熵越小傳感器位置越優,故2號優于4號,因此實驗與理論計算結果吻合。

圖6 2、4號自由度幅值譜Fig.6 Amplitude spectrum at DOFs 2 and 4

圖7 2、3號自由度幅值譜Fig.7 Amplitude spectrum at DOFs 2 and 3

2、3號自由度加速度響應幅值譜見圖7。由圖7看出,盡管兩者頻譜曲線無太大區別,但2號自由度頻譜曲線所含信息較3號豐富(表現為低頻段),此與表3的理論計算吻合。事實上,2號自由度信息熵比值為2.2,3號自由度信息熵比值為2.24。雖然2號自由度信息熵比值較3號小,但兩者數值較接近,由此表明圖7中兩者頻譜圖無太大區別原因。盡管如此,3號自由度仍可作為次優位置供傳感器布置。

研究結果表明,模型實驗與理論計算結果相符,進一步驗證本文的基于信息熵與譜有限元法傳感器優化布置方法的正確性。

4結論

(1) 提出基于信息熵指標與譜有限元法結合的傳感器優化布置方法,并對基于信息熵的傳感器優化布置方法進行推廣,與譜有限元法結合可避免用靜態形函數的傳統有限元方法建模對結構動力特性及傳感器優化布置影響,形成基于譜有限元法信息熵指標,用以定量評估結構模型參數識別結果的不確定性,并用貝葉斯統計系統識別理論識別該不確定性,將傳感器優化布置問題轉化為數值優化問題,極小化信息熵指標,進行傳感器布置優化問題求解。并通過實驗室彈性地基帶彈性接頭的周期管梁模型數值仿真、實驗數據研究,對本文所提方法進行驗證。

(2) 通過本文方法優化布置傳感器,可有效減少模型參數識別結果不確定性。在給定傳感器數目情況下,本文方法可有效增加傳感器采集數據中所含與結構模型參數識別有關的信息量,保證待識別模型參數不確定性最小,利于結構健康監測工作開展。本文方法亦可方便地對不同數目、不同位置傳感器配置下參數識別結果不確定性程度進行比較,以平衡實際應用中傳感器數量(或測試成本)與獲取信息量間矛盾,實用性明顯。雖僅以彈性地基的周期接頭管梁模型進行理論、實驗驗證,但對其它可用數值建模的一般工程結構均可適用。結構模型規模、復雜程度不同,計算量會有較大差別,即數值模型越復雜待布置傳感器位置越多,計算量會越大。

參 考 文 獻

[1] 何浩祥,閆維明,張愛林. 面向結構健康監測的傳感器數量及位置優化研究[J]. 振動與沖擊,2008,27(9):131-134.

HE Hao-xiang,YAN Wei-ming, ZHANG Ai-lin. Optimization of number and placement of sensors for structural health monitoring[J]. Journal of Vibration and Shock,2008,27(9):131-134.

[2] 何旭輝,陳政清,黃方林,等. 南京長江大橋安全監測和狀態評估的初步研究[J]. 振動與沖擊,2003,22(1):77-80.

HE Xu-hui, CHEN Zheng-qing, HUANG Fang-lin, et al. Preliminary studies on safety monitoring and state assessment for Nanjing Yangtse River Bridge[J]. Journal of Vibration and Shock, 2003, 22(1):77-80.

[3] Meo M, Zumpano G. Optimal sensor placement on a large scale civil structure[C]//Proceedings of SPIE-The International Society for Optical Engineering, v5394, Health Monitoring and Smart Nondestructive Evaluation of Structural and Biological Systems III, 2004:108-117.

[4] 黃民水,朱宏平,李煒明. 基于改進遺傳算法的橋梁結構傳感器優化布置[J]. 振動與沖擊,2008,27(3):82-86.

HUANG Min-shui, ZHU Hong-ping, LI Wei-ming. Optimal sensor placement on bridge structure based on genetic algorithm[J]. Journal of Vibration and Shock, 2008, 27(3):82-86.

[5] LeeU, Kim J, Andrew Y T L. The spectral element method in structural dynamics[J]. The Shock and Vibration, 2000, 32(6): 451-465.

[6] 黃民水,朱宏平,宋金強. 傳感器優化布置在橋梁結構模態參數測試中的應用[J]. 公路交通科技,2008,25(2):85-88.HUANG Min-shui, ZHU Hong-ping, SONG Jin-qiang. Application of optimal sensor placement in modal parameters test of bridge structure[J]. Journal of Highway and Transportation Research and Development, 2008, 25(2):85-88.[ 7] Papadimitriou C, Beck J L, Au S. Entropy-based optimal sensor location for structural model updating[J]. Journal of Vibration and Control, 2000, 6(5):781-800.

[8] Beck J L, Katafygiotis L S. Updating models and their uncertainties. I: bayesian statistical framework[J]. Journal of Engineering Mechanics, 1998, 124(4):455-461.

[9] 尹濤. 一種基于信息熵的分布參數結構傳感器/激勵器優化布置方法[J]. 振動與沖擊,2014,33(22):51-57.

YIN Tao. A probabilistic approach for optimal sensor/actuator configuration of distributed-parameter systems based on information entropy[J]. Journal of Vibration and Shock, 2014, 33(22):51-57.

[10] 伊廷華,李宏男,顧明. 基于MATLAB平臺的傳感器優化布置工具箱的開發及應用[J]. 土木工程學報,2010,43(12):87-93.YI Ting-hua, LI Hong-nan, GU Ming. Development of MATLAB based optimal sensor placement toolbox and its appliction[J]. Civil Engineering Journal, 2010, 43(12):87-93.

[11] Chow H M, Lam H F, Yin T, et al. Optimal sensor configuration of a typical transmission tower for the purpose of structural model updating[J]. Structural Control and Health Monitoring, 2011, 18(3): 305-320.

[12] 克拉夫R,彭津J,著.王光遠,譯.結構動力學[M].北京:高等教育出版社, 2006.

[13] Cho J, Go H, Lee U. Dynamic response of the spectral elment model by using the FFT[J]. Key Engineering Materials, 2007, 345/346:845-848.

[14] 尹濤,余嶺,朱宏平. 一種基于模型修正的結構損傷識別方法[J]. 振動與沖擊, 2007, 26(6):59-62.

YIN Tao, YU Ling, ZHU Hong-ping. Model updating based approach for structural damage idenfitication[J]. Journal of Vibration and Shock, 2007, 26(6):59-62.

[15] 雷英杰,張善文,李續武,等. MATLAB遺傳算法工具箱及應用[M].西安:西安電子科技大學出版社,2005.

[16] 尹濤,朱宏平,余嶺. 運用改進的遺傳算法進行框架結構損傷檢測[J]. 振動工程學報,2006,19(4):525-531.

YIN Tao, ZHU Hong-ping, YU Ling.Application study of an improved genetic algorithm for frame structure damage detection[J]. Journal of Vibration Engineering, 2006,19(4):525-531.

基金項目:國家自然科學基金資助項目(51208390);國家重點基礎研究發展 (973)計劃資助項目(2011CB013800);湖北省自然科學基金資助項目(2011CDB265);中央高校基本科研業務專項經費資助項目(271198, 273766)

收稿日期:2014-08-19修改稿收到日期:2015-01-13

通信作者尹濤 男,博士,副教授,1979年生

中圖分類號:O211;O321

文獻標志碼:A

DOI:10.13465/j.cnki.jvs.2016.02.013

Optimal sensor configuration based on spectral finite element method and information entropy

ZHANG Jia-pei1,2, YIN Tao1, ZHU Hong-ping2, DING Lan2

(1. School of Civil and Architectural Engineering, Wuhan University, Wuhan 430072, China;2. School of Civil Engineering and Mechanics, Huazhong University of Science and Technology, Wuhan 430074, China)

Abstract:In structural health monitoring (SHM) based on vibration data, the quantity and quality of the measured data, i.e., the number of sensors and the corresponding locations are very important for the success of SHM utilizing measured dynamic responses. In order to extract the most information from the measured data and reduce the uncertainties of the identified model parameters, a method of optimal sensor configuration for structural model parameters identification was presented. In order to avoid the influence of modeling error induced by traditional finite element method based on static shape function on the results of structural dynamic characteristics and optimal sensor placement, the spectral finite element method, being a dynamic modeling method with high-accuracy, was employed to model the target structure in the proposed method. In addition, the minimum of the uncertainties in model parameter estimates was taken as the optimality criterion for placing sensors, and the information entropy measure was used to quantify these uncertainties which were calculated by the Bayesian statistical identification method. The minimal information entropy measure was drawn from a set of possible sensor configurations to optimally locate a given number of sensors by using an integer-coded genetic algorithm. Both numerical simulation and laboratory experiment were carried out on a periodic pipe-beam model with flexible joints on elastic foundations to verify the effectiveness of the proposed method.

Key words:structural health monitoring; optimal sensor placement; information entropy; spectral finite element method; genetic algorithm; periodic structure

第一作者 張加培 男,碩士生,1989年12月生

郵箱:tyin@whu.edu.cn

主站蜘蛛池模板: 久久久久久午夜精品| 久久婷婷国产综合尤物精品| 亚洲精品成人片在线观看| 色综合久久88| 精品视频一区在线观看| 日韩国产精品无码一区二区三区| 国产理论一区| 青青青国产视频手机| 欧美精品在线看| 全免费a级毛片免费看不卡| 午夜精品久久久久久久无码软件| 91久久大香线蕉| 一级毛片在线播放| 成人免费黄色小视频| 狠狠综合久久| 欧美三級片黃色三級片黃色1| 精品国产免费第一区二区三区日韩| 五月婷婷综合色| 午夜毛片免费观看视频 | 特黄日韩免费一区二区三区| 国产香蕉97碰碰视频VA碰碰看| 亚洲成人高清在线观看| 性欧美在线| 午夜小视频在线| 国产av一码二码三码无码 | 国产精品欧美在线观看| 不卡网亚洲无码| 亚洲欧洲天堂色AV| 成人字幕网视频在线观看| 人妻精品久久久无码区色视| 丁香综合在线| 人妻丰满熟妇αv无码| 国产欧美性爱网| 日韩成人在线视频| 国产aⅴ无码专区亚洲av综合网 | 91成人精品视频| v天堂中文在线| 国产精品亚洲一区二区三区z | 日韩福利视频导航| 极品性荡少妇一区二区色欲 | 综合色区亚洲熟妇在线| 波多野结衣AV无码久久一区| 免费观看三级毛片| 亚洲av无码专区久久蜜芽| 亚洲日产2021三区在线| 福利在线免费视频| 精品久久蜜桃| 国产综合日韩另类一区二区| 欧美成人一区午夜福利在线| 亚洲天堂久久| 成人无码一区二区三区视频在线观看| 亚洲国产精品人久久电影| 国产xxxxx免费视频| 国产精品福利导航| 色偷偷av男人的天堂不卡| 亚洲中久无码永久在线观看软件| 久久久久国产一级毛片高清板| 亚洲第一在线播放| 婷婷色中文网| 天天婬欲婬香婬色婬视频播放| 国产欧美日韩91| 久久久精品无码一区二区三区| 在线人成精品免费视频| 无码国产偷倩在线播放老年人| 欧美在线视频不卡| 国内老司机精品视频在线播出| 国产白浆在线| 国产爽歪歪免费视频在线观看| 67194在线午夜亚洲| 国产在线观看第二页| 欧美中文字幕在线二区| 波多野结衣AV无码久久一区| 国内熟女少妇一线天| swag国产精品| 亚洲色图综合在线| 高清不卡一区二区三区香蕉| 日韩乱码免费一区二区三区| 狠狠色综合久久狠狠色综合| 欧美伦理一区| 精品福利网| 手机成人午夜在线视频| 97国产在线播放|