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

基于神經網絡和子空間的非線性系統載荷識別*

2022-11-04 08:37:04龐志雅馬志賽丁千
振動、測試與診斷 2022年5期
關鍵詞:模型系統

龐志雅,馬志賽,丁千

(1.天津大學機械工程學院天津,300350)

(2.天津市非線性動力學與控制重點實驗室 天津,300350)

引言

非線性系統廣泛存在于現實世界中,例如物體之間的接觸和滑動會引起摩擦非線性[1],加工超差、裝配誤差和使用磨損等因素會造成結構連接部位出現間隙非線性[2-3]等。隨著工程應用領域的不斷拓展和設計需求的日益增長,非線性系統的動力學反問題受到了越來越多的關注。結構動力學反問題包括系統辨識和載荷識別2類基本問題,但在處理和解決實際工程問題時,系統辨識和載荷識別常常是相互耦合的。

通過系統辨識手段描述非線性系統的動態特性和行為是解決其動力學問題的關鍵,非線性系統辨識已逐漸成為結構動力學領域的研究熱點之一。文獻[4-5]回顧了近幾十年來結構動力學中非線性系統辨識的研究進展,并對線性化、時域、頻域、時頻分析、非線性模態分析、黑箱建模和模型更新等參數估計方法進行了綜述,其中時域法由于無需將數據轉換至頻域,憑借可有效避免泄露和混頻等優點得到了迅速發展。文獻[6-7]提出了一種將非線性項視作標稱線性系統的內反饋力,進而利用線性子空間理論對非線性系統進行辨識的時域方法,即時域非線性子空間辨識方法,該方法具有良好的調節和計算能力,可以同時辨識多個非線性項,并且克服了由于信號相關性增加導致非線性項和標稱線性系統頻響函數在共振區附近估計精度較低的缺點。No?l等[8]采用時域非線性子空間方法對含強非線性特征的衛星結構進行了辨識,得到其標稱線性系統的模態參數和非線性項,并在辨識精度和計算效率等方面與頻域非線性子空間方法進行了比較。Sun等[9]采用時域非線性子空間辨識法對全動舵面進行了辨識,指出該方法由于利用幾何框架和數值線性代數分解,具有一定的魯棒性和高性能。Zhang等[10]將非線性系統分為標稱線性系統和局部非線性部分,分2步估計標稱線性系統的頻響函數和非線性部分的系數,提高了時域非線性子空間方法的精度與可靠性。

雖然子空間方法具有優良的性能,但對于非線性系統來說,非線性描述函數的選定仍然是辨識過程中的難點之一。例如,現有的多項式擬合法[11]、三線性函數法[12]等都會使辨識過程變得繁瑣,且具有一定的局限性。此外,不同的激勵水平往往會激發出具有明顯差異的非線性恢復力,使得每一次實驗都需要重復繁瑣的選定過程,這也在一定程度上造成了人力與時間的耗費。

在進行系統辨識時,時域非線性子空間辨識方法通常需要已知響應和外載荷數據,但在某些情況下,外載荷是未知或難以測量的。準確了解結構系統所承受的真實外載荷特性,對于結構響應預測、設計優化及故障診斷等都具有重要意義[13-14]。因此,載荷識別問題受到了許多學者的關注,并已發展出多種辨識方法[15-18]。楊智春等[19]對動載荷的識別方法進行了述評,并介紹了神經網絡在載荷識別中的應用。程光等[20]將一種改進的BP神經網絡方法用于雙跨轉子系統的載荷識別,通過拾取實驗信號進行預測證明了方法的有效性。然而,上述這些方法大多集中于處理與線性系統有關的載荷識別問題,但在實際工程中非線性系統更為普遍[21],且由于非線性恢復力的存在,其載荷識別問題往往更為困難。

神經網絡可以通過自適應訓練過程來學習變量之間的關系,進而代替傳統的機理建模[22]。經過大量數據訓練后的神經網絡模型可以等效代替響應-非線性恢復力映射關系,只需已知非線性位置的響應即可獲得非線性恢復力,可有效避免重復繁瑣的辨識過程。

因此,針對外載荷難以測量的情況,筆者旨在結合神經網絡和子空間法對非線性系統載荷識別方法進行研究,并開展相應的數值與實驗驗證。

1 時域非線性子空間辨識方法

1.1 狀態空間模型

對于一個包含h自由度的非線性時不變系統,其動力學方程可以表示為

其中:M,Cv和K分別為質量、阻尼和剛度矩陣;下標v表示黏性;z(t),(t),(t)分別為連續時間t時刻的位移向量、速度向量和加速度向量;f(t)為載荷向量;非線性項由p個分量的總和構成;μj為非線性參數;gj(t)為非線性描述函數,代表非線性類型;Lnj為非線性元素的位置向量,通常Lnj中元素的取值為1,0或-1。

引入狀態向量x(t)=[z(t),z˙(t)]T,則式(1)的狀態空間模型形式可以表示為

其中:u(t)=[f(t),-g1(t),…,-gp(t)]T為擴展輸入向量;下標c表示連續時間。

1.2 時域非線性子空間辨識

在實際測量中,采樣都是基于離散時間點的,故可以將式(2)描述的連續時間狀態空間模型轉換為離散時間狀態空間模型

其中:k為時間序列;A,B,C和D分別為離散時間狀態空間模型中的系統矩陣、輸入矩陣、輸出矩陣和反饋矩陣。

離散時間的系統矩陣和輸入矩陣與連續時間的系統矩陣和輸入矩陣的關系為

子空間方法的思想在于:通過未知的確定性系統生成測量數據總長度為s的輸入向量和輸出向量,利用線性代數工具辨識離散時間狀態空間模型中的狀態空間矩陣,一旦狀態已知,辨識問題就成為了未知系統矩陣中的線性最小二乘問題[23]。實現狀態空間矩陣的辨識首先需要確定擴展可觀測矩陣Γi。通過對式(3)進行迭代,可得到

根據式(5),輸入-輸出矩陣方程可以寫為

其中:下標p表示過去;f表示未來;塊行數i由用戶定義,應大于辨識系統的階數n。

類似于未來的輸入和輸出,未來狀態序列表示為Xf,狀態序列定義為

式(6)中的擴展可觀測矩陣Γi為

式(6)中的擴展下三角Toeplitz矩陣Hi為

將輸入、輸出數據收集到塊Hankel矩陣中

當列數l滿足l=s-2i+1時,意味著給定長度的數據可被充分利用。

根據式(10)和式(11),可將過去的輸入、輸出數據收集到塊Hankel矩陣Wp=def[Up/Yp]。Yf的行空間沿著Uf的行空間在Wp的行空間上的斜投影可定義為?i

其中:符號‘⊥’表示正交補;‘?’表示廣義逆[23]。

對加權斜投影進行奇異值分解

其中:W1和W2為用戶定義的權重矩陣,W1需滿秩,W2需遵循rank(Wp)=rank(WpW2)[23]。

當無噪聲存在時,斜投影只有n個非零奇異值;當噪聲存在時,則有n個以上非零奇異值。通過觀察奇異值可以確定模型階數n,進而確定U1和S1。可觀測矩陣則通過得到。

系統矩陣A通過確定,輸出矩陣C由Γi的前No行確定,No為輸出的數量;輸入矩陣B和反饋矩陣D則采用最小二乘求解以下線性方程組確定

在完成非線性系統的狀態空間矩陣的估計后,可進一步計算擴展頻響函數[6],即

其中:HL=(K+iωCv-ω2M)-1為非線性系統中的標稱線性頻響函數;下標L表示標稱線性;E表示擴展。

由式(15)可辨識得到標稱頻響函數及非線性參數。

2 基于神經網絡的載荷識別

2.1 神經網絡模型的訓練

為解決時域非線性子空間辨識過程中非線性描述函數選定困難造成的辨識過程重復繁瑣、計算耗時長等問題,本節旨在訓練一個從非線性位置響應信號映射到非線性恢復力的神經網絡模型,用以代替非線性描述函數的選定過程和時域非線性子空間辨識過程,如圖1所示。

圖1 神經網絡模型等效代替辨識過程的示意圖Fig.1 Schematic diagram of identification process equivalently represented by neural network models

針對非線性系統,基于大量數據對神經網絡模型進行訓練,網絡的輸入為非線性位置的響應,期望輸出為時域非線性子空間辨識中采用三線性函數法重構的非線性恢復力,使用所需的網絡層數、神經元個數、激活函數和訓練算法等參數構建網絡模型[19],然后進行迭代訓練,直到實際輸出與期望輸出之間的誤差達到可接受的范圍。訓練數據越多且標定質量越好,訓練后的神經網絡模型表征能力越強,越能以高精度逼近復雜的輸入-輸出關系。若對結果精度有更高的要求,可以采取增加隱含層或神經元的數量,獲得更大的訓練數據集,嘗試不同的訓練算法等方式。

需要注意的是增加訓練數據有助于提高神經網絡模型的逼近精度,但付出的代價是網絡模型可能會更加復雜,所需要的訓練時間也更長。因此,在用神經網絡模型等效描述響應-非線性恢復力映射關系時,通常需要在逼近精度和計算效率之間做適當折中。

2.2 載荷識別

根據非線性項視作反饋力的思想,式(1)可以寫為

其中:fsum(t)為非線性恢復力fnl(t)與外載荷f(t)的合力,此時系統可看作是一個受合力的線性系統;下標nl表示非線性。

與式(6)同理,載荷識別模型可以寫為

根據零初始條件x(0)=0,式(17)可以重寫為

一旦獲得標稱線性系統狀態空間矩陣AL,BL,CL和DL,就可以確定矩陣H,對式(18)進行最小二乘法估計,即可得到輸入合力fsum(t)。

根據式(16)可以計算作用在系統上的外載荷為

由基于子空間法建立的載荷識別模型可知,非線性恢復力是識別非線性系統外載荷的關鍵。因此,利用訓練后的神經網絡模型預測非線性恢復力,可使非線性系統載荷識別得以實現。需要注意的是,本研究識別的載荷為時域的時間歷程數據。

3 數值算例

間隙作為一種常見的非線性現象,廣泛存在于工程結構中。間隙的存在會改變工程結構的動力學特性,進而導致動力學響應預測困難,甚至引發工程結構出現精度低、壽命短等問題。在系統辨識中,間隙相關非線性參數的確定是非常關鍵的,也是消除和控制間隙的必要條件[24]。對于間隙非線性系統來說,需要辨識的非線性參數包括間隙值和接觸剛度。由于間隙值很小的誤差就會導致接觸剛度值出現較大的誤差[24],確定間隙值成為辨識間隙接觸剛度的前提。本節以3自由度間隙非線性結構為例,說明所提方法的過程和性能。

3.1 3自由度間隙非線性結構

考慮一個3自由度間隙非線性結構,如圖2所示,其動力學方程為

圖2 3自由度間隙非線性結構的示意圖Fig.2 Schematic diagram of a three degrees-of-freedom structure with clearance nonlinearity

其中:M,Cv和K分別為3自由度結構的質量、阻尼和剛度矩陣;zh(t),(t)和(t)分別為各自由度的位移、速度和加速度;fnl(t)為間隙非線性恢復力;fh(t)為各自由度的外載荷。

非線性恢復力可表示為

其中:kc為間隙接觸剛度;dc為間隙值;z2為非線性位置的位移響應。

式(21)表示的間隙非線性的恢復力曲線如圖3所示。

圖3 間隙非線性的恢復力曲線Fig.3 Restoring force curve of clearance nonlinearity

該3自由間隙非線性結構系統的基本參數如表1所示。

表1 數值算例的系統參數Tab.1 System parameters of the numerical example

3.2 時域非線性子空間的系統辨識

對該結構的每個自由度各施加一個激勵,選取的激勵均為零均值高斯白噪聲,即在任意兩個不同時刻的采樣信號均是統計獨立的。利用Runge-Kutta方法數值仿真得到振動位移響應。仿真過程中,采樣頻率設置為1 kHz,采樣時長為10 s,生成樣本數據長度為104。以1組激勵信號的均方根(root mean square,簡稱RMS)分別為39.445 4,40.178 3和39.839 8 N的仿真數據為例對辨識過程進行詳細說明。

對于間隙非線性系統來說,間隙值的確定是辨識間隙接觸剛度的前提。目前,已經發展了很多辨識間隙的方法,包括三線性函數法[6]、概率密度函數導數圖法[24]和恢復力曲面法[25]等。筆者使用三線性函數法對間隙值進行辨識。根據仿真響應數據,將正位移范圍[0,max(z(t))]劃分為11個等距的間隔[dj-1,dj]。只有當測量的正位移超過用戶定義的值時,每個非線性描述函數才有非零貢獻。參考式(21),11個非線性恢復力可被定義為

其中:dj為用戶定義的間隙值;fnlj和μj分別為第j個非線性項的非線性恢復力和非線性參數。

利用時域非線性子空間辨識方法,可以分別得到每個非線性項中非線性恢復力的貢獻和接觸剛度的貢獻,貢獻較大且集中的區域即為間隙值附近。區間[0,max(z(t))]內的非線性恢復力和非線性參數分布如圖4和圖5所示。

圖4 區間[0,max(z(t))]內的非線性恢復力分布Fig.4 Distribution of nonlinear restoring force in[0,max(z(t))]

圖5 區間[0,max(z(t))]內的非線性參數分布Fig.5 Distribution of nonlinear parameters in[0,max(z(t))]

為了得到更精確的間隙值,可以進一步縮小選定間隙范圍為[0.8 mm,1.2 mm],等距劃分為4份,如圖6和圖7所示。從圖中可以看出,間隙值主要集中在1 mm附近。圖8給出了重構的非線性恢復力特性曲線與真實非線性恢復力特性曲線的對比。同時,該方法也可以辨識結構的狀態空間矩陣,進而得到標稱線性頻響函數,如圖9所示。

圖6 區間[0.8 mm,1.2 mm]內的非線性恢復力分布Fig.6 Distribution of nonlinear restoring force in[0.8 mm,1.2 mm]

圖7 區間[0.8 mm,1.2 mm]內的非線性參數分布Fig.7 Distribution of nonlinear parameters in[0.8 mm,1.2 mm]

圖8 重構的非線性恢復力特性曲線與真實曲線Fig.8 The reconstructed and true nonlinear restoring force characteristic curve

圖9 辨識得到的標稱線性頻響函數H22曲線與真實曲線Fig.9 The identified and true nominal linear frequency response function H22 curve

結果表明,時域非線性子空間辨識方法可以獲取準確的非線性參數值,然而過程較為繁瑣,且每一次實驗都需要重復上述辨識過程。

3.3 響應-非線性恢復力映射關系的神經網絡模型訓練

使用25組不同激勵水平的載荷作用在結構上,利用時域非線性子空間辨識方法和三線性函數法重構對應的響應-非線性恢復力映射關系。用以訓練的神經網絡模型示意圖如圖10所示。

圖10 響應-非線性恢復力映射關系的神經網絡模型示意圖Fig.10 Schematic diagram of the neural network models for nonlinear restoring force-response mapping relations

圖10顯示的神經網絡模型含有2個隱含層,每個隱含層有10個神經元,使用Levenberg-Marquardt訓練算法,圖中w表示權重參數,b表示偏置參數。由于訓練過程中使用的數據較少,且響應-非線性恢復力映射關系較為簡單,在網絡模型結構選取時,通過對多種不同層數、單元數及激活函數進行預實驗以確定最優模型。從結果和計算效率的角度綜合考慮,選擇具有2個隱含層、每層10個神經元的神經網絡模型進行訓練即可達到較高精度。另外,反向傳播訓練算法Levenberg-Marquardt有利于梯度形式較簡單的激活函數,且激活函數的梯度計算涉及到優化過程,結合間隙非線性特征,筆者選取Positive linear函數作為激活函數,其中輸入為25組不同激勵水平下的非線性位置響應,輸出為時域非線性子空間辨識方法和三線性函數法重構的非線性恢復力。將輸入-輸出數據隨機分成3個集合,其中70%用于訓練,15%用于驗證,15%用于完全獨立的神經網絡模型泛化測試。使用與上節系統辨識所用的相同仿真數據,通過訓練后的神經網絡模型進行預測,發現預測值與真實值吻合良好,如圖11所示。

圖11 神經網絡模型預測的非線性恢復力特性曲線與真實曲線Fig.11 Predicted by neural network models and true nonlinear restoring force characteristic curve

結果表明,訓練后的神經網絡模型可以等效代替響應-非線性恢復力映射關系,即不論激勵水平大小如何,只需已知非線性位置的響應即可獲取非線性恢復力。

3.4 載荷識別

利用3.3節神經網絡模型預測的非線性恢復力和測量的振動響應可以對非線性系統的外載荷進行識別。同樣以激勵信號RMS分別為39.445 4,40.178 3和39.839 8 N的仿真數據為例對識別過程進行詳細說明。根據式(18)表示的基于子空間方法建立的載荷識別模型計算各自由度的合力,其中非線性位置的識別合力與真實合力對比結果如圖12所示。圖13為訓練后的神經網絡模型預測的非線性恢復力-時間關系。由式(19)可知,3自由度外載荷計算公式為

圖12 非線性位置識別的合力與真實合力Fig.12 Identified and true resultant force at the nonlinear position

為了量化非線性系統載荷識別結果的精度,引入平均絕對誤差評價指標(mean absolute error,簡稱MAE)

其中:q為數據總數量;f^為識別值;f為真實值。

結合圖12和圖13,可識別出非線性位置的外載荷與真實載荷,如圖14所示,其MAE值為3.243 2。對于1,3自由度來說,由于沒有非線性恢復力的存在,識別出的合力即為外載荷,以第1自由度載荷識別結果為例,如圖15所示,其MAE值為1.274 4。結果表明,識別出的外載荷曲線與真實載荷曲線吻合良好,驗證了所提載荷識別方法的有效性和可行性。

圖13 神經網絡模型預測的非線性恢復力與真實非線性恢復力Fig.13 Predicted by neural network models and true nonlinear restoring force

圖14 非線性位置處識別的外載荷與真實載荷Fig.14 Identified and true external force at the nonlinear position

圖15 第1自由度處識別的外載荷與真實載荷Fig.15 Identified and true external force at the first degree of freedom

4 實驗驗證

本節設計并搭建了1個3自由度間隙非線性實驗結構,用于對所提出的載荷識別方法進行實驗驗證。圖16和圖17分別為3自由度間隙非線性實驗結構的原理圖和裝置圖。該結構由3個鋁質質量塊、4根薄梁以及鋼質底座組成。由于鋁塊的質量相對較大,薄梁可以看作無質量彈性梁,質量塊的水平位移只與豎直薄梁的抗彎剛度有關,因而所研究結構可視為3自由度振動系統。由于實際條件限制,無法在每個自由度施加激勵,為識別非線性位置的外載荷,激振器作用在第2層。調節間隙大小的螺栓也設置在結構第2層。當激勵足夠大時,第2層鋁塊撞擊到具有固定間隙的螺栓,從而產生非線性行為。激振器上布有力傳感器用以測量實驗中真實外載荷,同時利用3個電渦流傳感器測量各層的位移響應,如圖17所示。進行不同激勵水平的實驗,每次采集數據時長為10 s,采樣頻率為1 kHz。

圖16 3自由度間隙非線性實驗結構的原理圖Fig.16 Schematic diagram of a three degrees-of-freedom experimental structure with clearance nonlinearity

圖17 3自由度間隙非線性實驗結構的裝置圖Fig.17 Set-up of the three degrees-of-freedom experimental structure with clearance nonlinearity

利用時域非線性子空間辨識方法重構的25組不同激勵水平下的非線性恢復力數據以及測量的非線性位置位移響應數據,對圖10所示的神經網絡模型進行訓練。使用RMS為4.190 0 N的白噪聲隨機激勵作用在實驗結構上,測量非線性位置響應,使用訓練后的神經網絡模型對非線性恢復力進行預測,發現預測結果與時域非線性子空間辨識結果幾乎一致,如圖18所示。圖19給出了標稱線性系統的頻響函數辨識結果。

圖18 神經網絡模型預測的非線性恢復力特性曲線與重構曲線Fig.18 Predicted by neural network models and reconstructed nonlinear restoring force characteristic curve

圖19 辨識得到的標稱線性頻響函數H22曲線Fig.19 The identified nominal linear frequency response function H22 curve

基于等效代替響應-非線性恢復力映射關系的神經網絡模型和實測的振動位移響應數據對外載荷進行識別。圖20為基于式(17)中載荷識別模型獲取的非線性位置合力識別結果。圖21為訓練后的神經網絡模型預測的非線性恢復力。結合圖20和圖21,根據式(23)可識別系統外載荷,如圖22所示,其中藍實線代表非線性位置的識別載荷,紅虛線代表實際測量的外載荷,MAE值為1.879 0。兩者趨勢基本一致,吻合良好,進一步驗證了所提載荷識別方法的有效性和可行性。

圖20 非線性位置識別的合力Fig.20 Identified resultant force at the nonlinear position

圖21 神經網絡模型預測的非線性恢復力Fig.21 Nonlinear restoring force predicted by neural network models

圖22 非線性位置識別的外載荷與實際測量的外載荷Fig.22 Identified and measured external force at the nonlinear position

5 結論

1)針對非線性描述函數選定困難,導致在不同激勵水平下時域非線性子空間方法辨識過程重復繁瑣的問題,采用多次重構的不同激勵水平下的非線性恢復力數據和非線性位置響應數據對神經網絡模型進行訓練,用以等效代替響應-非線性恢復力映射關系。該方法避免了辨識過程對系統模型的依賴性,使得系統不論受到的激勵水平大小如何,只需已知非線性位置的響應即可獲得非線性恢復力,提高了計算效率。

2)針對外載荷難以測量的情況,利用訓練后的響應-非線性恢復力神經網絡映射模型預測非線性恢復力,提出了基于神經網絡和子空間法的非線性系統載荷識別方法,并通過間隙非線性結構的數值與實驗研究驗證了所提載荷識別方法的有效性和可行性。

3)本研究僅基于間隙非線性系統對所提出的載荷識別方法進行了驗證,但理論上該方法也適用于其他類型的非線性系統。

猜你喜歡
模型系統
一半模型
Smartflower POP 一體式光伏系統
工業設計(2022年8期)2022-09-09 07:43:20
WJ-700無人機系統
ZC系列無人機遙感系統
北京測繪(2020年12期)2020-12-29 01:33:58
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
基于PowerPC+FPGA顯示系統
半沸制皂系統(下)
連通與提升系統的最后一塊拼圖 Audiolab 傲立 M-DAC mini
3D打印中的模型分割與打包
主站蜘蛛池模板: 欧美在线黄| 国产97公开成人免费视频| 免费高清毛片| 日韩av在线直播| 网友自拍视频精品区| 久久中文字幕不卡一二区| 日本亚洲欧美在线| 国产乱人伦偷精品视频AAA| 国产高清色视频免费看的网址| a级毛片网| 伊人久久青草青青综合| 国产女人在线视频| 麻豆国产原创视频在线播放| 亚洲一区二区三区麻豆| 国内毛片视频| 成人免费一级片| 天天综合色天天综合网| 在线日本国产成人免费的| 99精品在线看| 欧美一级夜夜爽| 一本久道久久综合多人| 手机永久AV在线播放| 色首页AV在线| 日韩精品久久久久久久电影蜜臀| av手机版在线播放| 国产欧美视频综合二区 | 女人一级毛片| 超清无码一区二区三区| 国产一区二区网站| 免费在线国产一区二区三区精品| 日本人妻丰满熟妇区| 日本www在线视频| 亚洲精品不卡午夜精品| 国产香蕉在线视频| 国产精品尤物在线| 国产成人精品18| 欧美激情视频一区二区三区免费| 亚洲欧美色中文字幕| 91美女视频在线| 久久久久国产精品免费免费不卡| 亚洲成在线观看| 久操中文在线| 午夜精品久久久久久久99热下载| 欧美综合激情| 综合人妻久久一区二区精品| 亚洲AⅤ波多系列中文字幕 | 国产激情在线视频| 天堂av高清一区二区三区| 久久国产高清视频| 最新国产网站| 日韩精品免费在线视频| 中文字幕有乳无码| 97视频免费在线观看| 在线免费无码视频| 97影院午夜在线观看视频| 亚洲国产一成久久精品国产成人综合| 亚洲人成电影在线播放| 国产丰满大乳无码免费播放 | 成人毛片免费观看| 欧美精品亚洲精品日韩专区| 5555国产在线观看| 97视频在线观看免费视频| 色欲国产一区二区日韩欧美| 青草视频在线观看国产| 亚洲国产精品无码AV| 黄色一及毛片| 一本色道久久88| 午夜福利无码一区二区| 午夜限制老子影院888| 国产香蕉一区二区在线网站| 一区二区在线视频免费观看| 亚洲国产精品日韩欧美一区| 国产熟睡乱子伦视频网站| 亚洲丝袜第一页| 国模私拍一区二区| 高清色本在线www| 久久亚洲国产最新网站| 欧美激情福利| 国产成人精品亚洲77美色| 欧美日韩91| 亚洲精品成人7777在线观看| 2024av在线无码中文最新|