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

基于嶺估計(jì)和L曲線的損傷識(shí)別技術(shù)研究

2016-01-18 03:39:19郭惠勇,羅樂,盛懋
振動(dòng)與沖擊 2015年4期

第一作者 郭惠勇 男,博士,副教授,1971年生

基于嶺估計(jì)和L曲線的損傷識(shí)別技術(shù)研究

郭惠勇1, 羅樂2, 盛懋1, 李正良1

(1. 重慶大學(xué) 土木工程學(xué)院, 重慶400045; 2. 成都城電電力工程設(shè)計(jì)有限公司,成都610041)

摘要:為了解決結(jié)構(gòu)損傷方程求解的可靠性問題,提出了基于嶺估計(jì)和L曲線的損傷方程求解方法。首先描述了基于模態(tài)應(yīng)變能靈敏度的結(jié)構(gòu)損傷方程,然后考慮測(cè)量噪聲等因素易造成病態(tài)的損傷方程問題,提出了采用嶺估計(jì)求解結(jié)構(gòu)的損傷方程,并利用L曲線確定最優(yōu)的嶺參數(shù),最后建立了損傷估計(jì)值的修正方法。數(shù)值仿真結(jié)果表明,基于修正的嶺估計(jì)和L曲線方法可以較為精確地識(shí)別出損傷的位置和程度,而且其定量識(shí)別結(jié)果明顯好于基本的嶺估計(jì)和L曲線方法以及一般的截?cái)嗥娈愔捣纸夥ā?/p>

關(guān)鍵詞:損傷識(shí)別;損傷方程;嶺估計(jì);L曲線;位移模態(tài)

基金項(xiàng)目:國家自然科學(xué)基金(51468058);教育部留學(xué)回國人員科研啟動(dòng)基金(20121707)

收稿日期:2014-02-08修改稿收到日期:2014-04-03

中圖分類號(hào):TB123文獻(xiàn)標(biāo)志碼: A

Structural damage identification method based on ridge estimation and L-curve

GUOHui-yong1,2,LUOLe2,SHENGMao1,LIZheng-liang1(1. School of Civil Engineering, Chongqing University, Chongqing 400045, China;2. Chengdu Chengdian Electric Power Engineering Design Co., Let., Chengdu 610041, China)

Abstract:In order to solve the reliability problem of structural damage equation, a damage detection method based on ridge estimation and L-curve was presented. Firstly, a structural damage equation based on modal strain energy sensitivity was described. Then, considering measurement noise easily making the damage equation be ill, the ridge estimation was utilized to solve the damage equation and L-curve is utilized to identify optimal ridge parameters. Finally, an improved ridge estimation method was proposed to enhance the identification precision. Simulation results demonstrated that the proposed ridge estimation and L-curve methods can detect structural damage locations and level with good accuracy, their quantitative detection results are obviously better than those of the basic ridge estimation and L-curve methods and the truncated singular value decomposition method.

Key words: damage identification; damage equation; ridge estimation; L-curve; displacement modal shape

工程結(jié)構(gòu)在投入使用后,由于環(huán)境的作用,常會(huì)產(chǎn)生各種損傷。結(jié)構(gòu)的嚴(yán)重?fù)p傷會(huì)引起建筑物的倒塌或者失效,從而造成巨大的經(jīng)濟(jì)和人員生命的損失。因此,對(duì)結(jié)構(gòu)損傷的識(shí)別研究一直是國際上的研究熱點(diǎn)[1-3]。利用損傷前后動(dòng)力特性的改變,建立相應(yīng)的損傷方程,從而求解出損傷系數(shù)是較為有效的方法[4-5]。但是,在工程實(shí)際中,測(cè)試數(shù)據(jù)不可避免的會(huì)受到測(cè)量噪聲、環(huán)境因素等的干擾,常會(huì)使損傷方程產(chǎn)生病態(tài)問題,從而給直接求解帶來困難,造成損傷的解極不穩(wěn)定,與真實(shí)值相差較遠(yuǎn)。這樣,損傷方程組的求解和精度成為問題的癥結(jié)。采用截?cái)嗥娈愔捣纸馑惴?TSVD)能夠減少方程的參數(shù)空間,使損傷方程組變的數(shù)值穩(wěn)定[6-7]。但是,該方法需要人為的選擇截?cái)嗾`差或者截?cái)嗨剑煌慕財(cái)嗾`差或截?cái)嗨綍?huì)產(chǎn)生不同的識(shí)別結(jié)果,故選擇恰當(dāng)?shù)慕財(cái)嗾`差或者截?cái)嗨經(jīng)Q定了方程求解的精度。而結(jié)構(gòu)的改變,甚至同一結(jié)構(gòu)單元?jiǎng)澐植煌紩?huì)影響截?cái)嗾`差或者截?cái)嗨降倪x取。嶺估計(jì)(RE)法[8]可以較好的處理病態(tài)問題,而L曲線方法[9]可以較好地處理優(yōu)化參數(shù)的選取問題,故本文擬采用嶺估計(jì)和L曲線法相結(jié)合,來解決損傷方程的病態(tài)問題以及優(yōu)化參數(shù)的選取問題。

1基于模態(tài)應(yīng)變能靈敏度的損傷方程

基于模態(tài)應(yīng)變能的損傷識(shí)別方法是較為有效的一種方法[10],可以根據(jù)模態(tài)應(yīng)變能的靈敏度建立相應(yīng)的損傷方程。

第j個(gè)單元和第r階模態(tài)的模態(tài)應(yīng)變能

(1)

式中:Erj是模態(tài)應(yīng)變能;Φr為第r階質(zhì)量歸一化位移模態(tài);Kj為第j個(gè)單元的剛度矩陣。則第j個(gè)單元和第r階模態(tài)的模態(tài)應(yīng)變能對(duì)任意設(shè)計(jì)參數(shù)p的敏感性為[6]

(2)

式中:λr是第r階特征值; K和M是整體剛度矩陣和質(zhì)量矩陣;p為設(shè)計(jì)參數(shù)。

當(dāng)結(jié)構(gòu)損傷時(shí),其第j個(gè)單元損傷后的剛度矩陣

Kjd=Kj-cjKj

(3)

SΔC=ΔR

(4)

式中:S是系數(shù)矩陣;ΔR是殘余向量;ΔC是損傷系數(shù)向量。其組成具體如下:

(5)

(6)

通過求解該方程可以同時(shí)識(shí)別結(jié)構(gòu)的損傷位置和程度。但是,由于結(jié)構(gòu)損傷常常出現(xiàn)在結(jié)構(gòu)的局部某些位置或單元,而大部分結(jié)構(gòu)單元未發(fā)生損傷,因此,在系數(shù)矩陣中常會(huì)出現(xiàn)幾乎線性相關(guān)的列,特別是模態(tài)參數(shù)在測(cè)試時(shí)會(huì)出現(xiàn)噪聲的干擾,故方程組的求解和精度往往成為問題的關(guān)鍵癥結(jié)。采用截?cái)嗥娈愔捣纸夥梢暂^為有效的處理該問題,但是,由于在現(xiàn)實(shí)中,損傷的位置和程度是未知的,截?cái)嗥娈愔捣纸夥ㄟx擇不同的截?cái)嗾`差或者截?cái)嗨剑糜?jì)算結(jié)果往往差別很大,從而對(duì)損傷識(shí)別結(jié)果造成較大的干擾[7]。故本文采取嶺估計(jì)和L曲線法相結(jié)合,來求解損傷方程,以獲取更具有可靠性的計(jì)算結(jié)果。

2嶺估計(jì)和L曲線法及其改進(jìn)

2.1損傷的嶺估計(jì)

可以采用嶺估計(jì)法求解損傷方程,以避免因測(cè)量噪聲等引起的損傷方程病態(tài)問題。嶺估計(jì)方法是一種在均方誤差意義下優(yōu)于最小二乘估計(jì)的有偏估計(jì)方法,是從Tikhonov正則化原理推導(dǎo)出來,該方法可以克服矩陣方程的病態(tài)性。根據(jù)Tikhonov正則化原理,對(duì)應(yīng)于損傷方程嶺估計(jì)的估計(jì)準(zhǔn)則為:

(7)

(8)

式中:I為單位陣。αI項(xiàng)的參與,可以有效的抑制方程的病態(tài)性,從而可以得到可靠的損傷估計(jì)值。嶺估計(jì)的核心是如何確定嶺參數(shù),這里選擇采用L曲線法來選取恰當(dāng)?shù)膸X參數(shù)。

2.2L曲線法

(9)

兩邊取對(duì)數(shù),可得

(10)

(11)

2.3估計(jì)值的修正

由于嶺估計(jì)是有偏估計(jì),而L曲線法求解的也是近似最優(yōu)解,故損傷的估計(jì)值需要修正。由嶺估計(jì)的估計(jì)準(zhǔn)則方程式(7)可知,該準(zhǔn)則是用來確定損傷系數(shù)的2-范數(shù)項(xiàng)最小,這樣,該嶺估計(jì)的解會(huì)造成近似線性的偏差,故采用簡單的線性修正方式來提高識(shí)別精度。另外,損傷系數(shù)值應(yīng)大于0,當(dāng)損傷系數(shù)估計(jì)值小于0時(shí),則認(rèn)為是無損傷,而該項(xiàng)本身是由于嶺估計(jì)的有偏估計(jì)造成的,故可根據(jù)所有無損單元的估計(jì)值平均偏移量來確定線性修正的幅值。這樣,就可以對(duì)可能損傷單元的損傷系數(shù)估計(jì)值進(jìn)行修正,修正公式如下

(12)

3數(shù)值算例

考慮如圖1所示的一個(gè)Euler-Bernoulli簡支梁。梁總長為6 m,劃分為20個(gè)單元,假設(shè)未損傷結(jié)構(gòu)的幾何參數(shù)和物理參數(shù)如下:梁的橫截面積A為0.005 m2,慣性矩I=1.67 m4,彈性模量E=32 GPa,密度ρ=2 500kg/m3,基準(zhǔn)參數(shù)為0.25。假設(shè)有三種損傷情況,第一種情況,在單元3和12發(fā)生損傷,剛度分別降低15%和20%;第二種情況,在單元2、10、14發(fā)生損傷,剛度分別降低25%,15%和20%;第三種情況,在單元5、11、19發(fā)生損傷,剛度分別降低20%,20%和30%。分別采用截?cái)嗥娈愔捣纸?(TSVD)法和嶺估計(jì)(RE)法求解損傷方程。由于低階模態(tài)比高階模態(tài)更能反映結(jié)構(gòu)的損傷,其抗噪能力更強(qiáng),而且更易被測(cè)量和獲取[11]。因此,本文僅采用第1階位移模態(tài)進(jìn)行結(jié)構(gòu)的損傷識(shí)別研究。嶺參數(shù)的確定主要是通過對(duì)式(11)求最大值κmax,κmax所對(duì)應(yīng)的點(diǎn)就是L曲線的所求的點(diǎn)。這樣就定位了L曲線上曲率最大的點(diǎn),這個(gè)點(diǎn)所對(duì)應(yīng)的α值就是所求的嶺參數(shù)。

圖1Euler-Bernoulli梁
Fig.1 Euler-Bernoulli beam

3.1例1

當(dāng)單元3和單元12發(fā)生損傷時(shí),其損傷系數(shù)的計(jì)算結(jié)果如圖2所示。采用了四種方案,即截?cái)嗾`差0.01的TSVD法、截?cái)嗾`差0.05的TSVD法、基本RE和L曲線法,修正RE和L曲線法。其中L曲線法求得的嶺參數(shù)為25.616 3。從圖2中可以發(fā)現(xiàn),截?cái)嗾`差0.01的TSVD法的損傷定位結(jié)果不理想,易對(duì)單元1和單元20產(chǎn)生誤判,其余三種方法定位識(shí)別較好,而從定量角度分析,可以較為明顯地發(fā)現(xiàn)修正RE和L曲線法的定量識(shí)別結(jié)果最好。將單元3和單元12的損傷系數(shù)單獨(dú)提取分析,結(jié)果如表1所示。從表1可以觀察到,修正RE和L曲線法的平均誤差最小,識(shí)別結(jié)果最好。

圖2 單元3、單元12損傷時(shí)四種損傷識(shí)別方案的識(shí)別結(jié)果 Fig.2 Identification results of four damage methods when multiple damages occur in the 3rd and 12th elements

識(shí)別方案單元3單元12計(jì)算結(jié)果計(jì)算結(jié)果總的平均誤差/%TSVD法(截?cái)嗾`差0.01)0.13520.18927.6TSVD法(截?cái)嗾`差0.05)0.16330.21498.2基本RE法+L曲線0.13540.19526.1修正RE法+L曲線0.14710.20081.2

圖3 單元3、單元12損傷時(shí)四種損傷識(shí)別方案的 識(shí)別結(jié)果(考慮3%的測(cè)量噪聲) Fig.3 Identification results of four damage methods when multiple damages occur in the 3rd and 12th elements(considering 3% noise)

當(dāng)考慮位移模態(tài)有3%的測(cè)量噪聲干擾時(shí),則四種損傷識(shí)別方案的一次計(jì)算結(jié)果如圖3所示。首先從損傷定位角度分析,截?cái)嗾`差0.01的TSVD法已經(jīng)喪失了損傷定位功能,截?cái)嗾`差0.05的TSVD法則具有一定的損傷定位效果,但是易對(duì)單元2和單元19產(chǎn)生誤判,而兩種RE和L曲線法的損傷定位效果較好。然后從損傷定量角度分析,從圖3可知,本建議的修正RE和L曲線法相對(duì)較好。故修正RE和L曲線法不僅具有較好的定位能力,具有一定的定量分析能力。由于測(cè)量噪聲具有一定的隨機(jī)性,則每次考慮隨機(jī)噪聲的計(jì)算結(jié)果都會(huì)產(chǎn)生一些小的變化,故未列出定量結(jié)果的平均誤差表。

3.2例2

當(dāng)單元2、單元10和單元14發(fā)生損傷時(shí),其損傷系數(shù)的計(jì)算結(jié)果如圖4所示。L曲線法求得的嶺參數(shù)為469.226 6。從圖4中可以發(fā)現(xiàn),截?cái)嗾`差0.01的TSVD法的損傷定位結(jié)果不理想,易對(duì)單元1、單元19、單元20產(chǎn)生誤判,其余三種方法的定位識(shí)別結(jié)果較好,再從定量角度分析,可以較為明顯地發(fā)現(xiàn)修正RE和L曲線法的定量識(shí)別結(jié)果最好。將單元2、單元10和單元14的損傷系數(shù)單獨(dú)提取分析結(jié)果見表2。從表2可知,修正RE和L曲線法的平均誤差最小,識(shí)別結(jié)果最好。

圖4 單元2、單元10和單元14損傷時(shí) 四種損傷識(shí)別方案的識(shí)別結(jié)果 Fig.4 Identification results of four damage methods when multiple damages occur in the 2nd, 10th and 14th elements

識(shí)別方案單元2單元10單元14計(jì)算結(jié)果計(jì)算結(jié)果計(jì)算結(jié)果總的平均誤差/%TSVD法(截?cái)嗾`差0.01)0.28610.15590.21348.4TSVD法(截?cái)嗾`差0.05)0.26290.12700.18659.3基本RE法+L曲線0.24260.12530.18499.0修正RE法+L曲線0.24360.14260.19393.5

圖5 單元2、單元10和單元14損傷時(shí)四種損傷識(shí)別 方案的識(shí)別結(jié)果(考慮3%的測(cè)量噪聲) Fig.5 Identification results of four damage methods when multiple damages occur in the 2nd, 10th and 14th elements(considering 3% noise)

當(dāng)考慮位移模態(tài)有3%的測(cè)量噪聲干擾時(shí),計(jì)算結(jié)果如圖5所示。首先從損傷定位角度分析,截?cái)嗾`差0.01的TSVD法已經(jīng)喪失了損傷定位功能,截?cái)嗾`差0.05的TSVD法則具有一定的損傷定位效果,但是易對(duì)單元3、單元18和單元19產(chǎn)生誤判,而兩種RE和L曲線法的損傷定位效果較好。然后從損傷定量角度分析,從圖3可知,本建議的修正RE和L曲線法相對(duì)較好。由于測(cè)量噪聲具有一定的隨機(jī)性,故未列出定量結(jié)果的平均誤差表。

3.3例3

當(dāng)單元5、單元11和單元19發(fā)生損傷時(shí),其損傷系數(shù)的計(jì)算結(jié)果如圖6所示。L曲線法求得的嶺參數(shù)為500.58。從圖6中依然可以發(fā)現(xiàn),截?cái)嗾`差0.01的TSVD法的損傷定位結(jié)果仍不理想,易對(duì)單元1、2和20產(chǎn)生誤判,其余三種方法的定位識(shí)別結(jié)果較好,再從定量角度分析,可以較為明顯地發(fā)現(xiàn)修正RE和L曲線法的定量識(shí)別結(jié)果最好。將單元5、單元11和單元19的損傷系數(shù)單獨(dú)提取分析結(jié)果見表3。從表3可知,本建議的修正RE和L曲線法的平均誤差最小,識(shí)別結(jié)果最好。

圖6 單元5、單元11和單元19損傷時(shí) 四種損傷識(shí)別方案的識(shí)別結(jié)果 Fig.6 Identification results of four damage methods when multiple damages occur in the 5th, 11th and 19th elements

當(dāng)考慮位移模態(tài)有3%的測(cè)量噪聲干擾時(shí),計(jì)算結(jié)果見圖7。首先從損傷定位角度分析,截?cái)嗾`差0.01的TSVD法已經(jīng)基本喪失了損傷定位功能,截?cái)嗾`差0.05的TSVD法則具有一定的損傷定位效果,但是易對(duì)單元2、3和18產(chǎn)生誤判,而基本和修正的RE和L曲線法的損傷定位效果較好。然后從損傷定量角度分析,從圖7可知,本建議的修正RE和L曲線法相對(duì)較好。由于測(cè)量噪聲具有一定的隨機(jī)性,故未列出定量結(jié)果的平均誤差表。

表3 例3的損傷定量識(shí)別結(jié)果

圖7 單元5、單元11和單元19損傷時(shí)四種損傷識(shí)別 方案的識(shí)別結(jié)果(考慮3%的測(cè)量噪聲) Fig.7 Identification results of four damage methods when multiple damages occur in the 5th, 11th and 19th elements(considering 3% noise)

從以上三個(gè)例子可以觀察到,當(dāng)僅采用第一階位移模態(tài)進(jìn)行損傷識(shí)別時(shí),截?cái)嗾`差為0.01時(shí),TSVD法的識(shí)別效果較差,截?cái)嗾`差為0.05時(shí),TSVD法的識(shí)別效果相對(duì)較好,故TSVD法的關(guān)鍵是選取恰當(dāng)?shù)慕財(cái)嗾`差值。由于結(jié)構(gòu)損傷是未知的,故截?cái)嗾`差的選擇也應(yīng)存在一個(gè)優(yōu)化問題。當(dāng)采用基本嶺估計(jì)和L曲線法,其定位和定量結(jié)果相對(duì)更好,這是因?yàn)長曲線可以選擇近似最優(yōu)的嶺參數(shù),從而解決了參數(shù)優(yōu)化問題。修正嶺估計(jì)和L曲線法的識(shí)別效果最好,這是由于該方法采用了未損單元的有偏估計(jì)均值對(duì)損傷系數(shù)估計(jì)值進(jìn)行了修正和改進(jìn)。

4結(jié)論

提出了基于L曲線和嶺估計(jì)的方法來求解損傷方程。該方法首先建立基于模態(tài)應(yīng)變能的損傷方程,然后利用嶺估計(jì)和L曲線法相結(jié)合來求解該方程。嶺估計(jì)可以較好的處理測(cè)量噪聲等引起的方程病態(tài)問題,L曲線法可以選擇近似最優(yōu)的嶺參數(shù),使嶺估計(jì)計(jì)算結(jié)果更為可靠。最后利用未損單元的有偏估計(jì)均值對(duì)嶺估計(jì)方法進(jìn)行了修正和改進(jìn)。數(shù)值分析結(jié)果表明,基本嶺估計(jì)和L曲線法可以較可靠地求解出損傷方程,并較好地識(shí)別出結(jié)構(gòu)的損傷,而修正的嶺估計(jì)和L曲線法可以更可靠地進(jìn)行進(jìn)行損傷的定性和定量識(shí)別,該方法優(yōu)于基本嶺估計(jì)和L曲線法和未經(jīng)過優(yōu)化的截?cái)嗥娈愔捣纸夥ā?/p>

參 考 文 獻(xiàn)

[1] Salawu O S. Detection of structural damage through changes in frequencies: a review[J]. Engineering Structures, 1997, 19(9): 718-723.

[2] Guo H Y, Li Z L. Structural damage identification based on Bayesian theory and improved immune genetic algorithm[J]. Expert Systems with Applications, 2012,39: 6426-6434.

[3] Seyedpoor S M. A two stage method for structural damage detection using a modal strain energy based index and particle swarm optimization[J]. International Journal of Non-Linear Mechanics, 2012, 47: 1-8.

[4] Shi Z Y, Law S S, Zhang L M. Structural damage localization from modal strain energy change[J].Journal of Sound and Vibration, 1998, 218(5): 825-844.

[5] 郭惠勇,李正良. 基于測(cè)量位移和頻率的結(jié)構(gòu)損傷二次識(shí)別方法[J].振動(dòng)與沖擊,2007,26(4): 94-96.

GUO Hui-yong, LI Zheng-liang. Two-stage structural damage identification method based on measured displacements and frequencies[J]. Journal of Vibration and Shock, 2007,26(4): 94-96.

[6] 顏王吉. 單元模態(tài)應(yīng)變能靈敏度及其在結(jié)構(gòu)損傷識(shí)別中的應(yīng)用[D]. 長沙:中南大學(xué),2008:38-45

[7] 張立濤,李兆霞,費(fèi)慶國. 基于加速度時(shí)域信息的結(jié)構(gòu)損傷識(shí)別方法研究[J].振動(dòng)與沖擊,2007,26(9): 138-141.

ZHANG Li-tao, LI Zhao-xia, FEI Qing-guo. Study on structural damage identification using acceleration data in time domain[J]. Journal of Vibration and Shock, 2007,26(9): 138-141.

[8] Hu Hong-chang. Ridge estimation of a semi-parametric regression mode[J]. Journal of Computational and Applied Mathematics, 2005, 176(1): 215-222.

[9] Hansen P C. Analysis of discrete ill-posed problems by means of the L-curve[J]. SIAM Review, 1992, 34(4): 561-580.

[10] 馮坤,楊永春,鄭軼. 基于軸向模態(tài)應(yīng)變能比法的三維桁架結(jié)構(gòu)損傷定位方法[J].振動(dòng)與沖擊,2013,32(12): 169-173.

FENG Kun, YANG Yong-chun, ZHENG Yi. Damage detection method of a three-dimensional truss structure based on axial modal strain energy ratio[J]. Journal of Vibration and Shock, 2013,32(12): 169-173.

[11] Guan H, Karbhari V M. Improved damage detection method based on element modal strain index using sparse measurement[J]. Journal of Sound and Vibration, 2008, 309(3-5): 465-494.

主站蜘蛛池模板: 自慰高潮喷白浆在线观看| 久久青草视频| 美女亚洲一区| 午夜激情婷婷| 久久a级片| 亚洲国产一区在线观看| av尤物免费在线观看| 亚洲天堂网视频| 日本人妻一区二区三区不卡影院| 亚洲国产高清精品线久久| 97se亚洲综合| 国产亚洲美日韩AV中文字幕无码成人| 欧美国产在线看| 国产区在线观看视频| 2021国产精品自拍| 亚洲小视频网站| 精品国产一区二区三区在线观看| 亚洲精品视频网| 国产肉感大码AV无码| 精品偷拍一区二区| 欧美日一级片| 99国产精品国产| 亚洲精品天堂自在久久77| av大片在线无码免费| 国产亚洲精| av大片在线无码免费| 日本高清有码人妻| 色综合天天综合| 日韩天堂网| 久久久久久尹人网香蕉| 亚洲综合极品香蕉久久网| 国产成人乱码一区二区三区在线| 成人综合网址| 国产一区二区丝袜高跟鞋| 狠狠色噜噜狠狠狠狠色综合久| 亚洲AⅤ无码日韩AV无码网站| 国产超薄肉色丝袜网站| 99在线小视频| 亚洲精品无码久久毛片波多野吉| 国产91av在线| 精品久久国产综合精麻豆| 97久久人人超碰国产精品| 无码在线激情片| 日本在线国产| 四虎永久免费在线| 毛片大全免费观看| 亚洲视频黄| 免费人成黄页在线观看国产| 无码电影在线观看| 国产精品99一区不卡| 精品無碼一區在線觀看 | 在线色国产| 国产一区三区二区中文在线| 亚洲第一综合天堂另类专| 精品人妻无码中字系列| 91精品国产福利| 国产h视频免费观看| 99er精品视频| a级毛片一区二区免费视频| 亚洲香蕉伊综合在人在线| 91福利在线观看视频| 国产日韩欧美一区二区三区在线| 亚洲欧美日本国产综合在线| 爆乳熟妇一区二区三区| 亚洲成A人V欧美综合天堂| 国产三级韩国三级理| 亚洲精品成人7777在线观看| 国产欧美专区在线观看| 97色婷婷成人综合在线观看| 国产在线一区视频| 搞黄网站免费观看| 国产精品视频导航| 国产在线观看一区精品| 91午夜福利在线观看| 永久天堂网Av| 国产91九色在线播放| 久久a毛片| 亚洲AⅤ无码日韩AV无码网站| 久久一级电影| 99视频在线看| 在线无码av一区二区三区| 国产精品三级av及在线观看|