亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        基于平均曲率模態(tài)和最小二乘支持向量機的混凝土拱壩損傷識別方法研究

        2013-08-09 01:50:27劉明軍馬奕仁郭法旺
        長江科學院院報 2013年11期
        關鍵詞:拱壩曲率振型

        李 波,劉明軍,馬奕仁,曹 浩,郭法旺

        (1.長江科學院a.工程安全與災害防治研究所;b.水利部水工程安全與病害防治工程技術研究中心;c.國家大壩安全工程技術研究中心,武漢 430010;2.中國電力投資集團公司南方分公司,廣州 510130;3.長江工程監(jiān)理咨詢有限公司,武漢 430010;4.中國水電工程顧問集團貴陽勘測設計研究院,貴陽 550081)

        基于平均曲率模態(tài)和最小二乘支持向量機的混凝土拱壩損傷識別方法研究

        李 波1a,1b,1c,劉明軍2,馬奕仁3,曹 浩1a,1b,1c,郭法旺4

        (1.長江科學院a.工程安全與災害防治研究所;b.水利部水工程安全與病害防治工程技術研究中心;c.國家大壩安全工程技術研究中心,武漢 430010;2.中國電力投資集團公司南方分公司,廣州 510130;
        3.長江工程監(jiān)理咨詢有限公司,武漢 430010;4.中國水電工程顧問集團貴陽勘測設計研究院,貴陽 550081)

        受眾多外界因素的影響,混凝土拱壩結構損傷與模態(tài)信息之間表現出明顯的非線性特征,這使得傳統(tǒng)的模態(tài)分析很難精確識別結構的損傷程度。針對上述問題,提出一種基于平均曲率模態(tài)和最小二乘支持向量機的混凝土拱壩損傷識別方法。該方法在數值模擬的基礎上,首先利用平均曲率模態(tài)對混凝土拱壩損傷位置進行識別,然后利用最小二乘支持向量機建立平均曲率模態(tài)和損傷程度間的非線性關系,對混凝土拱壩損傷程度進行識別。工程實例分析表明,該方法能有效地識別混凝土拱壩同時發(fā)生多處不同程度損傷的位置及損傷程度。

        平均曲率模態(tài);最小二乘支持向量機;混凝土拱壩;損傷識別

        1 研究背景

        近些年來我國制定了西部大開發(fā)戰(zhàn)略和“西電東輸”戰(zhàn)略,大壩建設事業(yè)得到了前所未有的發(fā)展,尤其是錦屏(壩高305 m)、溪洛渡(壩高273 m)、小灣(壩高292m)、拉西瓦(壩高254m)等特高混凝土拱壩的建設,對國民經濟作出了重大的貢獻,產生巨大的經濟效益。由于不利因素的影響,混凝土拱壩在服役期不可避免出現各種損傷,如果在該范圍內沒有觀測點,現有的監(jiān)測位移場、應力場和溫度場變化的靜態(tài)監(jiān)測方法將不會監(jiān)測到損傷的出現。這些損傷若不能及時發(fā)現,任其自由發(fā)展,將導致混凝土拱壩出現裂縫,危及大壩的安全。因此,對混凝土拱壩進行損傷識別成為亟待研究解決的一個重要課題。

        損傷會使結構的剛度降低,進而導致其結構動力特性(如模態(tài)頻率和振型等)發(fā)生變化。近年來,國內外許多專家學者對結構損傷識別進行了大量研究,并提出一系列基于模態(tài)分析的損傷識別方法,主要有基于固有頻率、振型和曲率模態(tài)的損傷識別方法?;诠逃蓄l率的損傷識別方法主要通過損傷前后固有頻率的變化識別結構的損傷狀況,該方法理論基礎清晰,測試簡單方便,但是固有頻率對結構早期損傷有時并不十分敏感,同時固有頻率是結構的整體性能描述,很難識別損傷的確切位置。基于振型的損傷識別方法,通過分析損傷前后的振型變化情況來識別結構損傷,但是振型對局部損傷的位置和程度不敏感。曲率模態(tài)實質上是振型的二階導數,可以通過各階振型得到,大量研究表明,曲率模態(tài)比固有頻率和振型對結構局部損傷更為敏感[1-2]。

        混凝土拱壩受眾多外界因素的影響,其結構損傷與模態(tài)信息之間表現出明顯的非線性特征,這使得傳統(tǒng)的基于模態(tài)分析識別結構損傷的方法遇到了困難,特別是對損傷程度的精確識別更是困難。人工神經網絡的出現在一定程度上解決了上述問題[3-4],然而,人工神經網絡存在過擬合、收斂速度慢、易發(fā)散、易陷入局部極值等實際問題,這嚴重影響了它的實用性。支持向量機是在統(tǒng)計學習理論基礎上發(fā)展起來的一種新的學習方法,它采用結構風險最小化原則,具有很強的泛化能力,并克服了神經網絡中存在的過擬合、收斂速度慢、容易陷入局部極值等缺點,在損傷識別中有著很好的應用前景[5]。

        針對上述問題,本文將模態(tài)分析技術與最小二乘支持向量機有機結合,提出基于平均曲率模態(tài)和最小二乘支持向量機的混凝土拱壩損傷識別方法。該方法在數值模擬的基礎上,首先,利用平均曲率模態(tài)對局部損傷靈敏性高的優(yōu)點,對混凝土拱壩損傷位置進行識別;然后,進一步利用最小二乘支持向量機泛化能力強的優(yōu)點,建立平均曲率模態(tài)和損傷程度間的非線性關系,對混凝土拱壩損傷程度進行識別。

        2 平均曲率模態(tài)方法

        基于位移模態(tài)計算出的曲率模態(tài),對于結構幾何尺寸及工作性能的變化能夠產生較明顯的效果,以梁為例來討論曲率模態(tài)的相關理論依據,其結論可以適用于任何類型的線性結構[6-7]。梁振動的微分方程為

        式中:E為梁的彈性模量;I(x)為梁的截面抵抗矩;v(x,t)為橫向振動位移;a1為剛度比例因子;m(x),c(x)分別表示梁的質量和阻尼。

        根據模態(tài)理論,方程(1)的解可以表示為模態(tài)貢獻的疊加形式,即

        式中:j為模態(tài)階數;φj(x),qj(t)分別為位移模態(tài)振型和模態(tài)坐標。

        依據材料力學理論中彈性梁彎曲變形曲線曲率與位移的關系,任意截面x處梁彎曲振動曲線的曲率變化函數為

        式中:φ″j(x)為梁的j階曲率模態(tài);k(x,t)為曲率;ρ(x,t)為曲率半徑。

        由材料力學知,梁的彎曲靜力關系式為

        式中M為梁截面彎矩。

        由式(3)和式(4)可知,曲率模態(tài)隨結構剛度的變化而變化,即曲率模態(tài)對結構損傷敏感,而且曲率模態(tài)與位移模態(tài)是一一對應關系。

        在模態(tài)分析中,通過頻響函數可求解出結構的模態(tài)參數,然后在位移模態(tài)振型基礎上按中心差分法近似計算曲率模態(tài)振型。當各節(jié)點之間距離相同時,第i個節(jié)點的第j階曲率模態(tài)為

        式中:φ″i,j為第j階、第i個節(jié)點的曲率模態(tài);l為各節(jié)點間的距離。

        當節(jié)點間的距離不相同時,式(5)應作以下修正:

        式中l(wèi)i為第i個節(jié)點與第i+1個節(jié)點之間的距離。

        第1個節(jié)點與最后1個節(jié)點的曲率模態(tài)無法通過式(6)算出,而它們不能簡單設為0(只有當該點位于簡支端或自由端的端點時,才可以作此處理),可以通過下式進行估算(假設共有n個測點):

        由模態(tài)振型衍生出來的曲率模態(tài)、曲率模態(tài)比、曲率模態(tài)差、平均曲率模態(tài)均可作為損傷因子。以下將重點分析平均曲率模態(tài)損傷因子(CDF),用下式表示為

        式中:N為所考慮的模態(tài)總數;φ″uj為無損結構的曲率模態(tài);φ″dj為損傷結構的曲率模態(tài)。

        3 最小二乘支持向量機

        支持向量機克服了人工神經網絡出現的問題,同時提高了計算效率。最小二乘支持向量機與標準支持向量機的主要區(qū)別在于采用不同的優(yōu)化目標函數,并且用等式約束代替不等式約束,它是一種很有潛力的數據分類和回歸工具[8-10]。已知一組訓練集(x1,y1),…,(xl,yl),x∈Rn,y∈R對于非線性問題,可以通過非線性變換φ(·)將輸入向量映射到高維特征空間,轉化為類似的線性回歸問題。最小二乘支持向量機利用結構風險最小化原則構造了下面的最小化目標函數:

        其中,w為權向量;γ為正則化參數;ei為樣本訓練誤差。

        同時將支持向量機算法中的不等式約束轉化為等式約束,即

        式中ai為Lagrange乘子。

        根據優(yōu)化條件:

        可以得到式(11)的最優(yōu)條件:

        消去式(13)的w和e,可得

        式中:y=[y1,y2,…,yl]T,a=[a1,a2,…,al]T,Z=[φ(x1),φ(x2),…,φ(xl)]T,1=[1,…,1]T。

        解這個線性方程組求得b和a,最小二乘支持向量機回歸函數為

        為了避免高維特征空間中的“維數災難問題”,用輸入空間的一個核函數K(x,xi)等效高維空間的內積形式,可以解決高維計算問題。所要求的回歸函數為

        本文選擇徑向基核函數:采用徑向基核函數的最小二乘支持向量機的主要參數是正則化參數和核函數寬度σ,這2個參數在很大程度上決定了最小二乘支持向量機的學習和泛化能力。常用網格搜索法來確定這2個參數。為了克服該方法計算量較大的缺點,本節(jié)采用改進的網格搜索法來確定參數,即先用一個大的步長進行粗糙的搜索,確定最好組合后,再在這個組合的周圍用小的步長進行更精細的搜索,具體步驟如下:

        (1)將γ和σ分別取M和N個值,組成M×N個γσ的組合。這里使用的參數范圍是:γ∈[22,24,…,224],σ∈[22,24,…,210]。

        (2)對訓練集歸一化后,對步驟(1)中劃分的組合采用交叉驗證。為了提高訓練的速度,在不影響模型精度的情況下,這里采用6折的交叉驗證,求得最高的學習精度,得到最好參數組合。

        (3)對步驟(2)中確定的參數正負22范圍內,再以20.5為步長進行更精細的網格搜索,求得精度最高的參數組合。

        相比一般的網格搜索法,改進的網格搜索法不僅可以減少訓練量,而且可以提高模型的預測精度。

        4 混凝土拱壩損傷識別步驟

        混凝土拱壩損傷識別的具體步驟如下:

        (1)根據實際資料,建立混凝土拱壩的三維有限元模型。

        (2)使用MSC.Marc強大的動力分析求解功能,求解出混凝土拱壩損傷前后的各階固有頻率。將損傷前后的自振頻率進行比較,初步判斷混凝土拱壩是否存在損傷。

        (3)求解出混凝土拱壩的損傷前后的各階位移模態(tài),利用上述介紹的曲率模態(tài)計算方法,計算出各種工況下的平均曲率模態(tài),對混凝土拱壩的損傷位置進行識別。

        (4)在已知損傷位置的情況下,計算出混凝土拱壩在所有可能損傷程度下的平均曲率模態(tài)。

        (5)以損傷單元所屬節(jié)點的平均曲率模態(tài)為最小二乘支持向量機的輸入樣本,實際損傷量為輸出樣本。

        (6)結合改進的網格搜索法,利用最小二乘支持向量機進行訓練和預測,識別出混凝土拱壩的損傷程度。

        5 工程實例

        某同心圓變半徑的混凝土重力拱壩,自左向右有28個壩段,壩頂高程為126.3 m,最大壩高為76.3 m,壩頂弧長419 m,壩頂寬8 m,最大壩底寬53.5 m,有限元計算模型的范圍:上下游方向取2倍左右壩高(各約150 m),左右壩肩各取約150 m。單元采用六面體8節(jié)點等參單元。模型共劃分了9 471個單元,11 997個節(jié)點。如圖1所示。

        假定人工邊界范圍以內的壩基均是“無質量的彈簧”,在形成整個系統(tǒng)的特征矩陣時,壩基單元只考慮彈性,不考慮質量,以消除壩體在振動時形成的振動波的傳遞效應,避免人為的放大作用。壩體混凝土的泊松比和材料密度分別取為0.167和2 400 kg/m3,壩體和基巖的彈性模量分別為24.0 GPa和19.5 GPa。

        圖1 混凝土拱壩有限元模型Fig.1 Finite elem entm odel of concrete arch dam

        在壩基-壩體-庫水的耦合系統(tǒng)中,假定壩基是無質量的彈性體,庫水是不可壓縮的,庫水動水壓力的影響按附加質量法考慮,折算為單位加速度相應的上游壩面結點附加質量,在對應結點上附加質量元。

        5.1 損傷位置的識別

        為了更清楚地描述損傷出現的位置,圖1給出了大壩上游面關鍵單元和所屬關鍵節(jié)點示意圖,圖中7939和4707為損傷單元的編號,1~11為關鍵節(jié)點的編號。為識別混凝土壩損傷的位置,本文在壩體的不同位置模擬了6種不同的損傷,詳見表1。

        表1 單元損傷情況Table 1 Conditions of element damage

        在正常蓄水位下,使用MSC.Marc強大的動力分析求解功能進行模態(tài)分析,各種損傷情況下,混凝土拱壩前4階自振頻率變化如表2所示。

        表2 各種損傷情況下的自振頻率Table 2 Natural frequency in damage conditions

        從表2中可以看出,在1處損傷和2處損傷下,模態(tài)自振頻率值隨損傷程度的加深而呈下降趨勢,在理論上證實了剛度降低導致頻率降低這一結論。但是變化的程度很小,實際測量中由于環(huán)境變化、噪音干擾等因素,很難識別混凝土拱壩的局部損傷。因此,以下采用平均曲率模態(tài)作為損傷因子對混凝土拱壩進行損傷識別。

        圖2為單元7939損傷前后各關鍵節(jié)點順河向位移的前4階平均曲率模態(tài),圖3為單元7939和4707損傷前后各關鍵節(jié)點順河向位移的前4階平均曲率模態(tài)。

        圖2 1處損傷下平均曲率模態(tài)Fig.2 M ean curvaturemode in the presence of one damage

        圖3 2處損傷下平均曲率模態(tài)Fig.3 M ean curvaturemode in the presence of two damages

        由圖2和圖3可以看出,在1處損傷和2處損傷下,損傷單元所屬兩節(jié)點的平均曲率模態(tài)損傷因子明顯突變,并且隨著損傷程度的增加,突變的程度增加。因此,平均曲率模態(tài)損傷因子可以有效地估計混凝土拱壩局部損傷的位置,并且可以定性分析損傷的程度。

        5.2 損傷程度的識別

        本文針對第5.1節(jié)中2處不同損傷進行研究。經過損傷位置分析可知混凝土單元7939和4707出現損傷,其所有可能的損傷程度為10%,20%,30%,構造2處損傷單元所有可能的損傷程度組合,并計算相應單元所屬節(jié)點的平均曲率模態(tài),各損傷單元訓練樣本為32=9組。單元7939和4707的訓練樣本分別見表3和表4。結合改進的網格搜索法,利用最小二乘支持向量機進行訓練,得到的最優(yōu)參數組合(γ,σ)分別為(1 048 576,181.019 3)和(1 048 576,128),最終訓練結果分別見表3和表4。

        從表3和表4可看出,單元7939和4707的訓練的相對殘差絕對值很小,說明訓練結果非常好。

        在上述建立的最小二乘支持向量機模型基礎上,對單元7939和4707分別出現15%和25%的損傷程度進行預測。預測結果分別見表5和表6,從表中可以看出,預測的相對殘差絕對值都很小,說明

        表3 單元7939的訓練樣本和訓練結果Table 3 Training sam p les and results for element 7939

        表4 單元4707的訓練樣本和訓練結果Table 4 Training samp les and results for element 4707

        表5 單元7939的預測樣本和預測結果Table 5 Prediction sam ple and result for element 7939

        表6 單元4707的預測樣本和預測結果Table 6 Prediction sam p le and result for element 4707

        建立的最小二乘支持向量機模型能有效地定量識別混凝土拱壩的損傷程度。

        6 結 論

        (1)混凝土拱壩模態(tài)自振頻率隨損傷程度的加深而呈下降趨勢,但是變化的程度很小,實際測量中由于環(huán)境變化、噪音干擾等因素,很難識別混凝土拱壩的局部損傷。

        (2)混凝土拱壩損傷處的前4階的平均曲率模態(tài)會發(fā)生突變,且隨著損傷程度的增大,突變程度也在增加。前4階的平均曲率模態(tài)是一個對混凝土拱壩損傷比較敏感的參數,可用于對混凝土拱壩的多處損傷位置進行有效的識別。

        (3)最小二乘支持向量機具有很強的泛化能力,將模態(tài)分析技術與最小二乘支持向量機有機結合,能有效地識別混凝土拱壩同時發(fā)生多處不同程度損傷的損傷位置和損傷程度。

        [1] 劉義倫,時圣鵬,廖 偉.利用曲率模態(tài)識別橋梁損傷的研究[J].振動與沖擊,2011,30(8):77-81.(LIU Yi-lun,SHI Sheng-peng,LIAO Wei.Bridge Damage I-dentification Using Curvature Mode Shapes[J].Journal of Vibration and Shock,2011,30(8):77-81.(in Chinese))

        [2] 張 旭,葛繼平,李胡生,等.基于曲率模態(tài)的箱梁損傷識別方法研究[J].上海應用技術學院學報(自然科學版),2012,12(3):201-206.(ZHANG Xu,GE Jiping,LI Hu-sheng,et al.Damage Identification of Box Beam Based on Curvature Theory Method[J].Journal of Shanghai Institute of Technology(Natural Science),2012,12(3):201-206.(in Chinese))

        [3] 李小榮,郭永剛.基于遺傳算法優(yōu)化神經網絡權值的大壩結構損傷識別[J].震災防御技術,2008,3(2):189-196.(LIXiao-rong,GUO Yong-gang.Dam Damage Identification on the Basis of Optimizing Neural Network Weight by Genetic Algorithm[J].Technology for Earthquake Disaster Prevention,2008,3(2):189-196.(in Chinese))

        [4] 于 菲,刁延松,佟顯能,等.基于振型差值曲率與神經網絡的海洋平臺結構損傷識別研究[J].振動與沖擊,2011,30(10):183-187.(YU Fei,DIAO Yan-song,TONG Xian-neng,et al.Damage Identification of an Offshore Platform Based on Curvature ofModal Shape Difference and BP Neural Network[J].Journal of Vibration and Shock,2011,30(10):183-187.(in Chinese))

        [5] 練繼建,李松輝.基于支持向量機和模態(tài)參數識別的導墻結構損傷診斷研究[J].水利學報,2008,39(6):652-658.(LIAN Ji-jian,LISong-hui.Damage Diagnosis of Spillway Guide Wall Based on Support Vector Machine and Modal Parameter Identification[J].Journal of Hydraulic Engineering,2008,39(6):652-658.(in Chinese))

        [6] 馬立英,周 ?A,彭曉?。誓B(tài)及其在汽車后橋損傷識別中的應用[J].同濟大學學報(自然科學版),2011,39(8):1208-1211.(MA Li-ying,ZHOU Hong,PENG Xiao-jun.Application of Curvature Mode Shape to Damage Identification of Rear Axle[J].Journal of Tongji University(Natural Science),2011,39(8):1208-1211.(in Chinese))

        [7] 張 靜,靜 行,袁海慶.基于曲率模態(tài)變化率指標的結構損傷識別[J].華中科技大學學報(城市科學版),2010,27(2):82-86.(ZHANG Jing,JING Hang,YUAN Hai-qing.Structural Damage Detection Based on Change Rate of Curvature Mode[J].Journal of Huazhong University of Science and Technology(Urban Science),2010,27(2):82-86.(in Chinese))

        [8] 李 波,顧沖時,李智錄,等.基于偏最小二乘回歸和最小二乘支持向量機的大壩滲流監(jiān)控模型[J].水利學報,2008,39(12):1390-1394.(LIBo,GU Chongshi,LI Zhi-lu,et al.Monitoring Model for Dam Seepage Based on Partial Least Squares Regression and Least Square Support Vector Machine[J].Journal of Hydraulic Engineering,2008,39(12):1390-1394.(in Chinese))

        [9] 李 波,徐寶松,武金坤,等.基于最小二乘支持向量機的大壩力學參數反演[J].巖土工程學報,2008,30(11):1722-1725.(LIBo,XU Bao-song,WU Jin-kun,et al.Back Analysis of Dam Mechanical Parameters Based on Least Squares Support Vector Machine[J].Chinese Journal of Geotechnical Engineering,2008,30(11):1722-1725.(in Chinese))

        [10]GU CS,LIB.Back Analysis ofMechanical Parameters of Roller Compacted Concrete Dam[J].Science China:Technological Sciences,2010,53(3):848-853.

        (編輯:黃 玲)

        Damage Identification of Concrete Arch Dam Using Mean Curvature Mode and Least Squares Support Vector Machine

        LIBo1,LIU Ming-jun2,MA Yi-ren3,CAO Hao1,GUO Fa-wang4
        (1.Yangtze River Scientific Research Institute,Wuhan 430010,China;2.South Branch of China Power Investment Corporation,Guangzhou 510130,China;3.Changjiang Project Supervision&Consultant Company,Ltd.,Wuhan 430010,China;4.Hydrochina Guiyang Survey and Design Institute,Guiyang 550081,China)

        Affected bymany external factors,the relation between concrete arch dam damage andmodal information is apparently nonlinear,whichmakes it difficult to accurately identify the degree of structural damage by traditional modal analysis.Aimed at these problems,amethod of identifying the damage of concrete arch dam by usingmean curvaturemode and least squares support vectormachine is proposed.On the basis of numerical simulation,the damage location is firstly identified usingmean curvaturemode,then the non-linear relationship betweenmean curvaturemode and damage degree is established using least squares support vectormachine to identify the damage degree.Engineering example shows that by using thismethod,the location and degree of various damages occurring simultaneously can be effectively identified.

        mean curvaturemode;least squares support vectormachine;concrete arch dam;damage identification

        TV312

        A

        1001-5485(2013)11-0113-06

        10.3969/j.issn.1001-5485.2013.11.023

        2013-09-10

        “十二五”國家科技支撐計劃項目(2012BAK10B04)

        李 波(1980-),男,湖北天門人,博士,主要從事水工結構安全的研究,(電話)027-82926142(電子信箱)lb007403@163.com。

        猜你喜歡
        拱壩曲率振型
        大曲率沉管安裝關鍵技術研究
        關于模態(tài)綜合法的注記
        力學與實踐(2022年5期)2022-10-21 08:10:34
        一類雙曲平均曲率流的對稱與整體解
        Phytochemicals targeting NF-κB signaling:Potential anti-cancer interventions
        縱向激勵下大跨鋼桁拱橋高階振型效應分析
        半正迷向曲率的四維Shrinking Gradient Ricci Solitons
        塔腿加過渡段輸電塔動力特性分析
        特種結構(2019年2期)2019-08-19 10:05:52
        淺議高拱壩壩踵實測與計算應力差異原因
        水電站設計(2018年1期)2018-04-12 05:32:02
        砌石雙曲拱壩拱冠梁設計的探討和實踐
        結構振型幾何辨識及應用研究
        山西建筑(2015年14期)2015-06-05 09:37:07
        99久久亚洲国产高清观看| 久久精品无码av| www插插插无码免费视频网站| 亚洲免费一区二区三区视频| 国产一区二区三区探花| 国产高颜值女主播在线| 亚洲中文字幕在线观看| 欧洲中文字幕| 亚洲影院在线观看av| 国产一区亚洲二区三区| 精品少妇人妻av无码专区 | 樱桃视频影院在线播放| 少妇厨房愉情理伦片bd在线观看| 五月天婷婷一区二区三区久久 | 国产免费无码9191精品| 亚洲熟少妇一区二区三区| 中国美女a级毛片| 激情偷乱人伦小说视频在线| 亚洲色欲色欲欲www在线 | 人妻体内射精一区二区三区| 国产精品国产午夜免费看福利| 国产精品久久久久久久专区| 午夜一区二区三区福利视频| 乱人伦中文视频在线| 中字幕久久久人妻熟女| 精品少妇后入一区二区三区| 手机在线国产福利av| 天天躁日日躁狠狠躁av麻豆| 精品久久亚洲中文无码| 99在线无码精品秘 入口九色| 风韵人妻丰满熟妇老熟| 免费看美女被靠的网站| 亚洲精品永久在线观看| 麻豆av一区二区天堂| 国产亚洲人成在线观看| 毛多水多www偷窥小便| 精品亚洲少妇一区二区三区| 久久久人妻一区二区三区蜜桃d | 国产美女自慰在线观看| 久久露脸国产精品WWW| 中文字幕一区二区三区精品在线|