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

基于經驗模態分解的沖擊載荷高效識別方法

2025-08-28 00:00:00劉玲楊曉明張力
機械強度 2025年8期

中圖分類號:TB123;O327 DOI: 10.16579/j.issn.1001. 9669.2025.08.010

0 引言

在運行或維護階段,飛機復合機翼、衛星太陽能蜂窩面板以及風力發電機葉片等結構易受外部碎片撞擊或工具跌落沖擊,導致結構的功能性和完整性被破壞,嚴重威脅運行安全[。以往沖擊損傷檢測大多依靠人工目視以及超聲檢測等地面無損檢測手段完成。然而,由于結構功能和復雜程度日益增加,沖擊事件引發的結構基體損傷通常難以用肉眼直接觀察。超聲和紅外熱成像檢測雖能探測肉眼無法觀察的內部損傷,但在面對大型結構時,實施過程耗時耗力,且成本高昂[2]。因此,為了提高沖擊事件的監測效率及降低監測成本,在結構上安裝傳感器并通過結構振動響應(如應變、加速度、位移等)間接識別沖擊事件已成為當前主流的監測方法。

由結構振動響應間接識別沖擊事件屬于結構動力學第2類反問題,其識別結果與使用的方法密切相關。趙剛等4提出了一種使用電阻應變片和改進三角測量法的低能沖擊定位新方法,能夠使定位誤差小于10cm 。趙發剛等5基于光纖光柵傳感網絡和小波包分解獲取了蜂窩板潛在沖擊位置的能量譜特征,并基于此特征實現了蜂窩板的全場沖擊定位。BOUKRIA等利用Tikhonov正則化方法識別了作用于圓形鋼板的沖擊載荷,其中L曲線法用于選擇最優正則化參數。

QIU等首先通過模式識別方法定位沖擊事件,然后使用Tikhonov正則化方法重構了沖擊載荷的時間歷程。喬百杰等[8-9]利用沖擊載荷的時空稀疏性,提出了沖擊載荷稀疏識別方法。

盡管對沖擊載荷識別方法的研究已經取得很大的進展,但仍有不足。例如,三角測量法要求傳感網絡捕獲彈性波的首達時間,采樣頻率極高,顯著增加了數據采集和處理設備的運行負擔[1];利用正則化方法識別沖擊載荷通常需提前建立所有潛在沖擊位置與測點的傳遞函數矩陣,導致求解模型維度巨大,計算耗時。LI等提出了一種兩步迭代算法,以較小的模型維度實現任意位置載荷的定位和辨識,但正則化方法通常需在最優正則化參數的設定下才能獲取最優解,當沖擊響應信號變化時,需反復選擇最優正則化參數。此外,為了改善反問題的不適定性,采用正則化方法也需布置多個測點[12]。

針對上述不足,提出了一種基于經驗模態分解的沖擊載荷識別新方法。首先,利用經驗模態分解技術將未知沖擊加速度響應信號分解為多個模態加速度響應;然后,通過度量低階模態加速度響應中未校正振型向量與振型矩陣列向量的共線性即可實現沖擊定位;最后,使用高斯基函數擬合沖擊載荷的時間歷程,并通過二維梯度下降法快速優化目標函數并獲取最優基函數參數。在懸臂板上進行試驗驗證,結果表明,所提出的方法僅需單個加速度計的測量數據即可精準、高效識別沖擊載荷。

1基本理論

1.1沖擊響應前向傳遞模型

考慮一個線彈性結構,假設其除邊界外的任意位置 q 受橫向沖擊載荷 f(t)=fqg(t) 作用,在零初始條件下,測點 p 的位移響應可由模態疊加原理[13]表示為

式中, n 為模態截斷數 ;fq 為沖擊載荷的幅值; ?r(p) 和?r(q) 分別為測點 p 和沖擊載荷作用點 q 的第 r 階質量歸一化振型; g(t) 為單位沖擊載荷的時間歷程(其最大值被設置為1); ? 表示卷積運算; hr(t) 為第 r 階位移脈沖響應函數,其表達式為

式中, 分別為第 r 階模態阻尼比和第 r 階固有頻率; ωd,r 為結構的第 r 階有阻尼固有頻率,且有 ωd,r=

對式(1)關于時間 t 進行兩次求導,可得結構在沖擊載荷作用下的加速度響應,為

式中, 為第 r 階加速度脈沖響應函數,其表達式為

1. 2 經驗模態分解

經驗模態分解是根據自適應篩選過程,將信號分解為有限數量的具有多尺度振蕩模式的固有模態函數和具有單調趨勢的殘差[14]。利用經驗模態分解,測量的完整沖擊加速度響應 ap,q(t) 可以表示為 m 個本征模態函數(模態加速度)和殘差項 r(t) 的總和,即

式中, ap,q,r(t) 為第 r 階本征模態函數或第 r 階模態加速度響應,其表達式為

為降低模型的復雜程度并提高沖擊定位的效率,沖擊載荷的時間歷程在低頻段可由幅值為 fq 的理想沖擊函數近似表示,即 f(t)=fqδ(t-t0) ,其中 t0 表示沖擊事件發生的時間瞬間。因此,式(5)中的低頻分量可表示為

將式(6代人式(4),分離模態加速度響應中的低頻和高頻分量,式(4)可進一步表示為

ap,qlow(t)+ap,qhigh(t)+r(t)

式中, ?m 為分解的本征模態函數(IntrinsicModeFunction,IMF)的總數; N 為低階模態響應的個數。

由于實測的加速度響應信號中含有噪聲,從原始加速度響應信號提取的各階模態加速度響應中可能包含多個瀕率成分,因此無法精確地表征結構真實的模態加速度響應。為了獲得準確的模態加速度響應,可使用帶通濾波器對模態加速度響應信號進行濾波,從而確保每階模態加速度響應信號的主頻率為對應階次的結構固有頻率[15]]

2沖擊載荷識別方法

2.1 沖擊定位

在式(6)中,結構固有頻率、模態阻尼比以及測點的振型可通過試驗模態分析或有限元分析獲取。根據閾值法,沖擊事件發生的時間瞬間 t0 可設置為加速度響應絕對值第一次達到最大絕對值 5% 的瞬間,即 t0= 因此,未知沖擊加速度響應的低頻分量中僅含有沖擊載荷幅值和沖擊位置的振型兩個未知參數。將模態加速度響應按照固有頻率由低到高排列,其低頻分量可進一步描述為

式中, 為一個與測點振型和系統模態參數有關的時間函數; r(q) 為未校正振型向量。

由式(8)可知,向量 r(q) 由沖擊載荷幅值與沖擊載荷作用位置振型向量的乘積表示,由于結構各個潛在沖擊位置存在差異性,當選擇合適的參數 N 時,向量r(q) 的方向唯一。

考慮到加速度測量數據為離散數據,令采樣點的個數為 ns ,采樣間隔為 Δt ,則 r(q) 的表達式可根據式(6)和式(8)推導為

式中, θr 為向量的夾角; 為符號判定函數,其表達式為

令結構上所有潛在沖擊位置的前 N 階振型構成的振型矩陣為 ? (其可由有限元分析或試驗模態分析獲取)

式中, M 為潛在沖擊點的數量。

對比式(9)和式(11)可知,理想情況下,當未知沖擊載荷作用在位置 qm 時,未校正振型向量 r(q) 與振型矩陣 ? 的第 m 列共線。因此,沖擊載荷定位問題可以被轉換成搜索未校正振型向量與振型矩陣列向量共線性的問題。由于建模誤差和測量噪聲的存在, r(q) 與 ? 的第 m 列不可能完全共線,因此可認為共線性最大的列所在的位置為未知沖擊載荷的作用位置。此外,由于結構振型存在對稱性,為確保沖擊定位不定位到結構對稱點,振型矩陣 ? 應滿足列滿秩條件,這也是選擇參數 N 的基本原則。

2.2沖擊載荷時間歷程重構

根據沖擊載荷的時間分布特征,其時間歷程可由如下高斯函數近似擬合為

式中, μ,σ 分別為高斯函數最大值對應的時間坐標和高斯函數的方差。沖擊載荷的幅值可在沖擊載荷位置確定后估算為

當沖擊載荷幅值被確定,沖擊載荷時間歷程重構的問題便轉化成了以下求解高斯基函數 g(t) 最優參數μ 和 σ 的優化問題。

式中, T(μ,σ) 為損失函數,表示實測加速度響應與計算得到的加速度響應之間的累計最小二乘誤差;xp,q 分別為實測的加速度響應序列和由式(3)計算得到的加速度響應序列; xp,q(ti) 和 分別為實測的加速度響應序列和計算的加速度響應序列的第 i 個采樣點。

根據沖擊載荷時間歷程的分布特征,本文將 σ 的取值范圍設置為 0~1×10-3 ,初始值設定為 5×10-4 ,取值間隔設置為 1×10-5 。變量 τ 的取值范圍被限制為 [t0,tc] ,初始值設定為 (t0+tc)/2 ,并以 Δt 為間隔進行離散化。其中, tc 表示 xp,q 絕對值最大值對應的時間。為了減少計算成本,采用二維梯度下降法來求解式(14)所表述的優化問題。由于高斯基函數的可微性,損失函數 T(μ σ 迭代步 下的梯度可以根據以下等式求解

式中, w=[τ,σ]T 。下一次迭代更新的參數可表示為

式中, γ 為梯度下降的步長,其取值設置為0.01。此外,停止迭代的閾值設為 1×10-6 O

2.3 實施流程

上述沖擊載荷識別方法的總體流程如圖1所示。由圖1可知,所述方法包含3個具體實施步驟:結構模態參數獲取、沖擊載荷定位以及沖擊載荷時間歷程重構。其中,步驟1在離線準備階段實施,用于獲取結構的模態參數和振型矩陣。步驟2和步驟3屬于在線實施步驟。當未知沖擊載荷作用在結構上時,數據采集系統記錄下加速度計測量的沖擊響應信號;隨后,利用經驗模態分解技術將加速度響應信號分解成多個模態加速度響應,并利用式(6)和式(10)計算得到向量 r(q) ,通過度量 r(q) 與振型矩陣 ? 列向量的余弦相似度即可獲取未知沖擊載荷的作用位置信息;最后,根據定位結果計算未知沖擊載荷的幅值并將其代入式(14)建立優化目標函數,利用梯度下降法迭代獲取最優基函數參數,將最優參數代入式(12)即可獲取沖擊載荷的時間歷程。

3 試驗驗證

懸臂板是工程中廣泛應用的結構,飛機復合機翼、衛星太陽能蜂窩面板以及風力發電機葉片等結構在一定條件下都可簡化為懸臂板結構。因此,為了驗證所提沖擊載荷識別方法在實際結構上的有效性,將一個尺寸為 600mm×300mm×3mm 的鋁合金懸臂板作為測試對象,并進行了一系列試驗分析。懸臂板物理參數如表1所示,試驗裝置如圖2所示。

由圖2可知,懸臂板表面被劃分為12個正方形區域,并且每個區域的中心被假設為潛在的沖擊載荷作用點,為了清楚區分各個潛在沖擊點,所有點被按照1~12的順序進行了編號[圖2(a)」。在本試驗中,沖擊載荷通過型號為INV9311的力錘對板表面實施敲擊產生。力錘前端配備了靈敏度為 10mV/N 的力傳感器,用于精確測量沖擊載荷的時間歷程數據。懸臂板的橫向振動加速度信號由型號為PCB365A03的壓電加速度計進行采集,加速度計的靈敏度為 10mV/g 。沖擊載荷信號與加速度信號通過NICompact-DAQ數據采集系統進行同步采集,采樣頻率設定為 5kHz 。所有數據的后續處理及未知沖擊載荷的識別分析在一臺配備i7-8550U處理器、8GB內存的筆記本電腦上完成。

圖1利用模態加速度響應的沖擊載荷識別方法整體流程圖

Fig.1Overallflow chart of impact load identification method using modalacceleration response

表1懸臂板的物理參數

Tab.1 Physicalparametersofcantileverplate

3.1懸臂板沖擊載荷定位

為確定懸臂板的模態參數,首先對懸臂板進行了試驗模態分析,振型、固有頻率、模態阻尼比如圖3所示。根據模態分析結果,結合振型矩陣 Φ 的列滿秩條件,確定了N的取值范圍為 N?4 ,即利用經驗模態分解獲取板的前4階模態加速度響應即可實現該板的沖擊載荷定位。

使用尼龍錘頭敲擊了板上12個潛在沖擊點,得到了12組加速度響應。圖4所示為在懸臂板上6#點受到垂直于板表面的沖擊載荷作用時,測點的加速度響應。

根據圖1所示的沖擊載荷定位步驟,采用經驗模態分解法對加速度信號進行處理,成功提取了板的前4階模態加速度響應信號。對每個信號進行帶通濾波后,其頻率主成分均與板的固有頻率相對應。由經驗模態分解提取的前4階模態加速度如圖5所示。由圖5可知,分解得到的各階模態加速度響應在時間坐標上的分布規律與式(6所描述的關系一致。

圖2 懸臂板試驗裝置Fig.2Testsetup of thecantileverplate

圖46#點受沖擊載荷作用時的加速度響應

Fig.4Acceleration response of point6# under impact load

根據模態分析結果,計算得到了模態加速度脈沖響應函數,通過將模態加速度脈沖響應函數以及測點的振型代入式(8),得到了向量 進而,根據式(9)計算得到未校正振型向量 r(q) 。最后,通過計算向量 r(q) 與振型矩陣 ? 列向量的共線性,得到了沖擊載荷定位結果,如圖6所示。由圖6可知,未校正振型向量 r(q) 與振型矩陣 ? 的第6列余弦相似度最大,表明沖擊位置被精準定位到。此外,觀察圖6可得,余弦相似度的值在坐標軸上呈近似對稱分布,這是因為板的振型存在對稱性。在實際工程應用中,為了進一步提高沖擊載荷定位結果的準確性,可適當增加模態加速度響應的階數。

圖5板的模態加速度響應

Fig.5Modalaccelerationresponsesoftheplate

圖7所示為板上剩余11個潛在沖擊點的沖擊定位結果。由圖7可知,真實沖擊載荷作用位置均與余弦相似度最大值對應的橫坐標取值相符,表明所有的沖擊載荷均被準確定位到。此外,懸臂板上12個潛在沖擊點沖擊定位的平均耗時約為0.8s,表明所提沖擊載荷定位方法具有很高的計算效率。

為了進一步驗證所提方法在不同錘頭激勵下定位沖擊源的有效性,還使用橡膠錘頭和鋁合金錘頭敲擊了懸臂板,并進行了沖擊定位測試,沖擊定位結果如表2所示。由表2可知,當錘頭材質更堅硬時,沖擊定位成功率更高。這是因為堅硬錘頭產生沖擊力的持續時間更短,其與理想脈沖更為接近。

表2不同錘頭類型下的沖擊定位結果

3.2懸臂板沖擊載荷時間歷程重構

當未知沖擊載荷定位完成后,還需根據定位結果將沖擊載荷的時間歷程進行重構。為了衡量重構效果,選用了峰值相對誤差 ePRE 以及相對誤差 eRE 兩種評價指標[1],可表示為

由式(17)可知,指標 ePRE 則側重于評估重構的沖擊載荷峰值與真實沖擊載荷峰值的差異,而指標 eRE 側重于評估重構的沖擊載荷與真實沖擊載荷整體的差異,

為了降低模態截斷誤差,本研究在重構過程中將模態截斷數設置為6。根據沖擊載荷定位結果,利用式(13)計算得到了沖擊載荷的幅值。通過將幅值以及板的前6階加速度脈沖響應函數代人式(14),得到了沖擊載荷時間歷程重構的優化目標函數。隨后,使用二維梯度下降法求解最優擬合參數,重構了作用在12個潛在沖擊點的沖擊載荷。表3所示為沖擊載荷時間歷程重構的結果,其中包括沖擊載荷的幅值、最優參數 、ePRE 指標以及 eRE 指標。由表3可知,所有重構的沖擊載荷與真實沖擊載荷的 ePRE 指標均小于 10% ,并且 eRE 指標也小于 40% ,表明所提沖擊載荷重構方法具有良好的精度。圖8所示為作用在2、5、9和12號點的沖擊載荷的時間歷程及對應的重構結果。由圖8可知,重構的沖擊載荷與真實沖擊載荷的曲線幾乎重合,表明了重構結果的準確性。

表3懸臂板沖擊載荷時間歷程重構結果

Tab.3 Reconstructionresults of impactload time history for

4結論

基于以上的理論分析和試驗驗證,得到如下主要結論:

1)為了降低傳感器布線的復雜度和數據采集系統的采樣成本,提出了一種基于經驗模態分解技術的沖擊載荷識別新方法。該方法僅需單個加速度計即可實施。

2)在總體尺寸為 600mm×200mm×3mm 的鋁合金懸臂板結構上驗證了所提出方法的有效性,36次沖擊定位測試的成功率為91. 67% 。

3)采用高斯函數近似擬合沖擊載荷的時間歷程,重構結果的 ePRE 指標和 eRE 指標分別小于 10% 和 40% ,表明重構結果十分精準。

參考文獻(References)

[1]XIAO D,SHARIF-KHODAEI Z,ALIABADI M H. Impact force identification for composite structures using adaptive waveletregularised deconvolution[J]. Mechanical System and Signal Processing,2024,220:111608.

[2]劉以龍,劉杰,劉江南.基于子結構分析的動態載荷和模型參數 復合反演研究[J].機械強度,2013,35(5):553-558. LIU Yilong,LIU Jie,LIU Jiangnan. Research on composite inversion of dynamic loads and structural parameters based on substructure analysis[J]. Journal of Mechanical Strength,2013,35 (5):553-558.(In Chinese)

[3]襲著有,閆云聚,常曉通.基于遺傳算法的動態載荷識別優化方 法[J].機械強度,2015,37(4):593-597. XI Zhuyou,YAN Yunju,CHANG Xiaotong. Optimization method of dynamic load identification based on genetic algorithm[J]. Journal ofMechanical Strength,2015,37(4):593-597.(In Chinese)

[4]趙剛,李書欣,劉立勝,等.應變片在復合材料低能量沖擊定位中 的應用[J].振動、測試與診斷,2018,38(3):526-530. ZHAO Gang,LI Shuxin,LIU Lisheng,etal. Low energy impact localization of composite materials based on resistance strain gauge [J].Journal of Vibration,Measurement amp; Diagnosis,2018,38(3): 526-530.(In Chinese)

[5]趙發剛,周春華,梁大開,等.衛星典型復合材料蜂窩結構板的沖 擊定位方法[J].振動、測試與診斷,2016,36(6):1204-1209. ZHAO Fagang,ZHOU Chunhua,LIANG Dakai,et al. Impact and locating method research on satelite’s typical composite honeycomb structure panel[J].Journal of Vibration,Measurement amp; Diagnosis,2016,36(6):1204-1209.(In Chinese)

[6]BOUKRIA Z, PERROTIN P, BENNANI A,et al. Experimental impact force location and identification using inverse problems:: application fora circular plate[J].International Journal of Mechanics,2011,5(1):48-55.

[7]QIU BB,ZHANG M,XIEY G,et al. Localisation of unknown impact loads on a steel plate using a pattern recognition method combined with the similarity metric via structural stress responses in the time domain[J].Mechanical Systems and Signal Processing,2019,128:429-445.

[8]喬百杰,陳雪峰,劉金鑫,等.機械結構沖擊載荷稀疏識別方法研 究[J].機械工程學報,2019,55(3):81-89. QIAO Baijie,CHEN Xuefeng,LIU Jinxin,et al. Sparse identification of impact force acting on mechanical structures[J]. Journal of Mechanical Engineering,2019,55(3):81-89.(InChinese)

[9]LIU JJ,QIAO BJ,WANG Y N,et al. Group sparsity extension of Non-convex sparse regularization via convex optimization for impact force identification[J].Mechanical Systems and Signal Processing,2023,201:110661.

[10]GOUTAUDIER D,GENDRE D,KEHR-CANDILLEV,et al. Long-range impact localizationwitha frequencydomaintriangulation technique:application to a large aircraft composite panel[J]. Composite Structures,2020,238:111973.

[11]LI QF,LUQH. Force localization and reconstruction usinga twostep iterativeapproach[J].Journal ofVibrationand Control,2018, 24(17):3830-3841.

[12]LIURX,DOBRIBANE,HOUZC,et al.Dynamic load identification formechanical systems:areview[J].ArchivesofComputationalMethods inEngineering,2022,29(2):831-863.

[13]蔡芳盛,嚴剛,曾捷,等.基于非負貝葉斯正則化的沖擊載荷光纖 識別法[J].振動、測試與診斷,2023,43(3):427-434. CAIFangsheng,YANGang,ZENGJie,etal.Anoptical fiber techniqueof impact load identification method based on non-negative Bayesianregularization[J].Journal ofVibration,Measurementamp;Diagnosis,2023,43(3):427-434.(InChinese)

[14]LI C,WANG XL,TAO ZY,et al. Extraction of time varying information from noisy signals:anapproach based ontheempirical modedecomposition[J].Mechanical Systemsand Signal Processing,2011,25(3):812-820.

[15]YANGJS,FU ZY,ZOUYF,etal.A response reconstruction method based on empirical mode decompositionand modal synthesismethod[J].Mechanical Systems and Signal Processing,2023, 184:109716.

[16]邱雨晴,王磊,王曉宇,等.基于改進函數擬合法的沖擊載荷識別 研究[J].機械工程學報,2022,58(3):157-166. QIUYuqing,WANG Lei,WANG Xiaoyu,et al.Researchon impact force reconstruction based on improved function fitting method[J]. Journal ofMechanical Engineering,2022,58(3):157- 166.(In Chinese)

Abstract: Aiming atthe problemsof traditional impact load identification methods,suchas the requirement foralarge numberofsensors,highamplingfrequencyandlowidentificationaccuracy,anewimpactloadidentificationmethodbased onempirical mode decomposition(EMD)technology was proposed.TheEMDtechnologywasused todecompose the complete impactresponsetoobtain the modalaccelerationresponse.Theimpactlocation wasquicklyrealized by measuring thecolinearity betweentheuncorrectedmodeshape vectorandthecolumn vectorof themodeshapematrixinthemodal accelerationresponse.According tothepositioningresults,anoptimization objective function wasconstructed.The time historyof the impact load wasfitted byusing the Gausian basis function,and the optimal fiting parameters werequickly solved by using the two-dimensional gradient descent method.Testsconducted onacantilever platewith dimensionsof (2號 600mm×200mm×3mn n show that with only one accelerometer,the success rate of 36 impact positioning tests is 91.67% .The peak relative error and relative error index of the reconstruction results are less than 10% and 40% ,respectively.

Keywords:Impactloadidentification;Empiricalmodedecomposition;Modalacelerationresponse;Uncorrctedoe vector; Basis function fiting Correspondingauthor:ZHANGLi,E-mail:zhangli@jcut.edu.cn

Fund:Natural Science Foundation ofHubei Province (025AFCoo5); Jingmen Major Scienceand Technology Innovation PlanProject (2024ZDYF004);Jingmen ScienceandTechnologyPlan Project (2024YDKY23);Jingchu Universityof

Technology Doctoral Startup Fund Project(YY202444) Received:2024-12-19 Revised:2025-02-24

主站蜘蛛池模板: 在线日韩一区二区| 日本黄色不卡视频| 91 九色视频丝袜| 麻豆精品视频在线原创| 精品在线免费播放| 日本欧美精品| 97青青青国产在线播放| 国产情侣一区二区三区| 国产白浆视频| 国产香蕉在线视频| 国产女人在线观看| 国产在线拍偷自揄观看视频网站| 国产乱子伦精品视频| 不卡午夜视频| 99精品高清在线播放| 亚洲国产精品久久久久秋霞影院 | 亚洲va在线观看| 先锋资源久久| 色偷偷男人的天堂亚洲av| 伊人久久大香线蕉影院| 日本高清视频在线www色| 欧美全免费aaaaaa特黄在线| 欧美翘臀一区二区三区| 538国产视频| 亚洲国产91人成在线| 色色中文字幕| 亚洲人成网站在线播放2019| 日韩午夜伦| 日本成人一区| 思思热在线视频精品| 91在线中文| 成人免费视频一区| 人妻少妇乱子伦精品无码专区毛片| 国产精品色婷婷在线观看| 国产福利大秀91| 色妞www精品视频一级下载| 在线欧美a| 国产精品v欧美| 国产玖玖视频| 伊人天堂网| 有专无码视频| 亚洲自拍另类| 91在线无码精品秘九色APP | 日韩无码黄色| 噜噜噜久久| 国产无遮挡猛进猛出免费软件| 午夜啪啪福利| 欧美一区精品| 伊人AV天堂| 亚洲欧洲日韩综合| 久久香蕉国产线| 成人免费一级片| 日本精品视频一区二区| 无码国产伊人| 亚洲三级a| 99999久久久久久亚洲| 欧美日韩动态图| 国产在线无码av完整版在线观看| 操操操综合网| 国产丰满大乳无码免费播放 | 国产免费久久精品99re丫丫一| 欧美亚洲激情| 国产成人成人一区二区| 99精品这里只有精品高清视频| 国产在线一二三区| 欧美日韩中文字幕二区三区| 伊人成人在线| 精品乱码久久久久久久| 欧美成人手机在线视频| 日韩精品中文字幕一区三区| 又污又黄又无遮挡网站| 国产高清在线丝袜精品一区| 91在线无码精品秘九色APP| 国产亚洲精| 一级毛片无毒不卡直接观看| 91色在线视频| 精品福利视频导航| 久久青草精品一区二区三区| 99无码中文字幕视频| 国产精品一区在线观看你懂的| 久久久国产精品无码专区| 国内精品九九久久久精品|