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

基于增廣卡爾曼濾波器的時域傳遞路徑分析方法

2024-12-31 00:00:00朱雨何智成趙亮
中國機械工程 2024年9期

摘要:時域傳遞路徑分析方法用于解決瞬態工況下復雜系統振動噪聲問題不僅未能解決自然頻率附近頻響函數矩陣的病態問題,而且利用現有頻域信息轉換提取得到的所需時域信息精度較低,因此提出一種基于增廣卡爾曼濾波器的時域傳遞路徑分析方法。該方法采用增廣卡爾曼濾波器輔以遺傳算法估計時域工況載荷,通過最小二乘算法辨識單位脈沖響應函數,將時域工況載荷和對應的單位脈沖響應函數進行線性卷積以計算各傳遞路徑的時域貢獻量。算例表明,所提方法采用的增廣卡爾曼濾波器載荷識別誤差小于傳統方法的去卷積濾波器所識別載荷的誤差,最小二乘算法辨識的單位脈沖響應函數誤差小于對頻響函數直接進行快速逆傅里葉變換或者構造有限單位脈沖響應濾波器的誤差,且所提方法在復雜結構上也同樣具有較小的誤差。

關鍵詞:時域傳遞路徑分析方法;增廣卡爾曼濾波器;遺傳算法;最小二乘算法

中圖分類號:TB535

DOI:10.3969/j.issn.1004-132X.2024.09.020

Time-domain TPA Method Based on AKF

ZHU Yu1 HE Zhicheng1 ZHAO Liang2

1.School of Mechanical and Vehicle Engineering,Hunan University,Changsha,410082

2.NIO Technology(Anhui) Co.,Ltd.,Hefei,230061

Abstract: When time-domain TPA methods used for tackling vibrations and noise problems in complex systems during transient conditions, which failed to overcome the ill-conditioned problems of the frequency response functions matrix near natural frequencies, and showed low accuracy in extracting the required time-domain information from the existing frequency-domain data. Therefore, a new time-domain TPA method was developed based on AKF. This method started with estimating time-domain operational loads using AKF supplemented by GA, and then identified the impulse response functions using LS algorithm. Finally, time-domain contributions for each transfer path was computed by linearly convolving time-domain operational loads with respective impulse response functions. Case study demonstrates that the load-identification errors of AKF applied in the proposed method are smaller than that of traditional deconvolution filters. Additionally, the errors of impulse response functions identified by LS algorithm are smaller than those obtained by direct inverse fast Fourier transform or creating finite impulse response filters from frequency response functions. Furthermore, the proposed method achieves small errors even in complex structures.

Key words: time-domain transfer path analysis(TPA) method; augmented Kalman filter(AKF); genetic algorithm(GA); least square(LS) algorithm

0 引言

機械系統在運作時難免會產生振動噪聲(noise,vibration,harshness,NVH)現象。過度的振動噪聲不僅對產品的穩定性、可靠性和使用壽命造成威脅,還會導致較差的用戶體驗。在對復雜機械系統進行NVH問題分析時,目標測點的響應往往是由各激勵通過不同傳遞路徑產生的影響線性疊加而來。為了深入解析并優化這些問題,提升產品的整體性能,需對激勵源及各傳遞路徑進行全面的分析。傳遞路徑分析(transfer path analysis,TPA)作為一種關鍵的試驗性分析方法,能夠準確識別出對目標測點振動噪聲貢獻最大的主要傳遞路徑,為機械系統的NVH問題提供一種優化策略和解決方案。

近年來TPA方法不斷發展,目前已經形成了傳統TPA方法、基于組件的TPA方法和基于傳遞率的TPA方法三類[1]。根據處理數據類型的不同,又可分為頻域TPA方法和時域TPA方法兩類。前者適用于穩態或準穩態工況,能夠較好地反映機械系統整體的頻域特性,由于其方法成熟且計算效率高,故在NVH研究中應用廣泛[2-5]。后者既適用于穩態工況又適用于瞬態工況,甚至可以進行路徑回放及聲品質分析,是對頻域TPA方法的補充,特別是在處理更為復雜的系統問題時具有較明顯的優越性。

最經典的時域TPA方法通過構造有限單位脈沖響應(finite impulse response,FIR)濾波器獲取單位脈沖響應函數,采用去卷積技術求解時域工況載荷,并且已被集成到LMS Test.Lab軟件,成為應用最廣泛的時域TPA方法。郝鵬等[6]提出一種新的運動聲源時域TPA方法,較頻率分析法具有較高的分析精度和更高的計算效率。褚志剛等[7]建立某汽車發動機對車內副駕駛位置噪聲的時域TPA模型,識別出了主要傳遞路徑。謝小平等[8]創新性地提出了一種基于時域TPA的瞬態聲品質分析方法和流程,成功地找到了對聲品質貢獻較大的路徑,并且通過虛擬修改各路徑傳遞函數值,優化了車內聲品質。呂來等[9]將頻域工況傳遞路徑分析(operational transfer path analysis,OTPA)方法推廣到時域,利用指示點和目標點的時域工況數據實現路徑分解,提出了時域OTPA方法。隨著消費者對產品性能和舒適度期望的增長,時域TPA方法有了一些發展與應用,然而,該方法在實踐中仍然面臨諸多挑戰。一是在構建去卷積濾波器的過程中需要對頻響函數(frequency response function,FRF)矩陣求逆,但在自然頻率附近,FRF矩陣往往呈現病態性。二是此方法需要進行大量的頻時轉換操作,不僅增加了計算負擔,同時難以保證轉換過程的精確度。因此,目前時域TPA方法仍然存在著計算復雜以及分析精度較低等問題。

增廣卡爾曼濾波器(augmented Kalman filter,AKF)[10]是一種將未知輸入融入狀態變量,進而與狀態一起進行最優估計的方法。該方法具有實時分析且未知輸入估計精度較高的優點[11-12],因此本文基于AKF提出一種新的時域TPA方法,并采用軟件仿真分析了懸臂薄板模型和某車輛模型在時域載荷激勵下目標點的振動傳遞情況;通過與傳統方法比較,說明了新方法的優越性。

1 傳統時域TPA方法

TPA理論由“源-路徑-接受體”模型演變而來,對于線性時不變系統,目標點處的響應可視為多個激勵源通過不同傳遞路徑產生影響的線性疊加。對于時域TPA方法,其分析模型如圖1所示,目標點處的響應可以表達為如下形式:

式中,t為時間變量;ye(t)為第e個目標點的總響應;yypi(t)為第i個結構載荷對目標點的貢獻,上標“yp”表示結構載荷至目標點區域;n為結構傳遞路徑的數量;yyqj(t)為第j個聲學載荷對目標點的貢獻,上標“yq”表示聲學載荷至目標點區域;w為聲學傳遞路徑的數量。

由于線性時不變系統目標點響應是各激勵信號與對應單位脈沖響應函數(impulse response function,IRF)卷積的線性疊加,故式(1)又可寫為

式中,τ為時間變量;pi(τ)、qj(τ)分別為第i條、第j條路徑上的結構載荷和聲學載荷;hypei(t-τ)、hyqej(t-τ)分別為結構載荷、聲學載荷至目標點e的非耦合IRF。

顯然,時域工況載荷的識別和IRF的獲取是進行目標點振動或噪聲時域分析的基礎,時域TPA方法主要包括以下步驟。

(1)工況數據采集。在實驗對象上安裝傳感器,采集運行工況下各感興趣位置處的響應信號,包括指示點響應信號d(t)、a(t)以及目標點響應信號y(t)。

(2)濾波器構造。使用力錘(激振器)或體積聲源獲得感興趣位置的FRF并構造矩陣,包括結構載荷、空氣載荷到指示點的FRF矩陣HDP(ω)、HAQ(ω)以及到目標點的FRF矩陣HYP(ω)、HYQ(ω),對HDP(ω)、HAQ(ω)進行奇異值分解或正則化處理并求其廣義逆以構造去卷積濾波器,對HYP(ω)、HYQ(ω)進行相應處理以構造相應的卷積濾波器,其中上標“DP”和“AQ”分別表示結構載荷和聲學載荷至指示點區域,“YP”和“YQ”同“yp”和“yq”,大小寫不同是為了在時頻域做出區分。

(3)時域工況載荷識別。得到步驟(1)的指示點響應信號d(t)和a(t)以及步驟(2)的指示點去卷積濾波器后,可以采用去卷積的方法進行時域工況載荷識別:

式中,上標“+”表示求廣義逆;IFFT(·)表示構造去卷積濾波器;p(t)和q(t)表示各路徑結構載荷和聲學載荷構成的向量。

(4)時域傳遞路徑分析。利用步驟(2)得到的目標點卷積濾波器和步驟(3)得到的時域工況載荷p(t)、q(t),通過式(2)便可進行目標點的時域傳遞路徑分析。

2 基于AKF的時域TPA方法

本研究提出一種新的時域TPA方法,該方法首先通過AKF估計時域下的工況載荷,隨后采用最小二乘算法(least squares algorithm,LS)來識別各路徑下的IRF,最后將載荷與相應的IRF進行線性卷積計算,得到各傳遞路徑下的時域貢獻量。

2.1 基于AKF的時域工況載荷識別

對于一個已經在空間上離散的線性系統,僅考慮結構載荷作用下的影響,它的連續時間動力學方程可以用二階微分方程表示:

式中,u(t)為系統位移向量;K、C和M分別為系統剛度、阻尼和質量矩陣;f(t)為系統載荷向量;Sp為指定各結構載荷空間分布的矩陣。

其中,wk和νk分別具有相關的協方差矩陣Q(w)=E(wkwTl)(E(·)表示求期望值,上標“w”代表過程誤差wk,下標“k”和“l”表示不同時刻)和R=E(vkvTl),分別用來表征模型和測量的不確定性,這里假設Q(w)和R均不隨時間變化。

AKF使用包含輸入力的增廣狀態向量,將結構響應和輸入力一起估計為狀態的一部分。為建立AKF的完整系統方程,需要未知輸入力模型,隨機游走模型是輸入力建模的一般方法,連續域的輸入力隨機游走模型狀態方程定義為

其中,ξ*0為初始狀態(可以通過測量或者先驗知識獲取);ξ^*0為初始狀態估計;P^0為初始估計協方差(一般ξ^*0和P^0預先設置,通常ξ^*0估計得越準,P^0設置得越小);ξ^*+k-1為k-1時刻的后驗狀態估計;ξ^*-k為k時刻的先驗狀態估計;P^+k-1為k-1時刻的后驗估計協方差;P^-k為k時刻的先驗估計協方差;Q*定義為Q(w)和Q(η)組成的塊對角矩陣,分別是模型和輸入不確定性的協方差矩陣,假設Q*不隨時間變化(則Q*k-1為常值);Kk為k時刻的卡爾曼增益;Rk為k時刻的測量噪聲協方差矩陣(因為R不隨時間變化,因此Rk是常值)。

由于ξ^*+k=ξ^+kp^+k包含未知輸入的后驗狀態估計,所以AKF可對時域載荷進行最優估計。

2.2 基于遺傳算法的指示點響應參數選擇

由于AKF及諸多載荷識別方法均受到指示點響應參數的影響,因此擬采用遺傳算法(GA)確定指示點響應參數,進一步提高載荷識別精度。指示點響應參數一般包含指示點個數、類型和位置,其中指示點個數和指示點類型在分析前均可提前確定。

一般來說,指示點個數越多,捕獲的結構響應信息越豐富,載荷識別的結果往往越準確,但是過多的指示點常常會導致數據處理上的復雜性增加,一般選擇指示點個數為激勵數的2倍(或者更多)。指示點類型有加速度、速度和位移(應變),現有的文獻資料表明[13-15],當使用加速度和位移(應變)響應信息時,AKF估計結果最穩定,也最準確,前期的仿真分析結果表明該組合同樣適用于傳統方法。指示點位置通常利用計算得到,這里提出應用GA來確定最佳指示測點位置,流程如圖2所示。其中第i(i=1,2,…,n,n為結構載荷個數)個結構載荷的均方根誤差公式定義為

這里通常通過計算所有結構載荷的均方根誤差和來搜索最佳指示測點位置,該公式也可以用來計算IRF和各路徑時域貢獻量的均方根誤差。

2.3 基于LS的IRF辨識及傳遞路徑貢獻量分析

對于時域TPA方法,其傳遞路徑一般用IRF來表征。在信號處理領域,對于單輸入單輸出的線性因果系統,如果輸入是因果序列,則系統e位置的響應與i位置激勵的關系在時域上可用式(2)的卷積形式來表示(這里假設圖1僅有第i個結構載荷作用于系統):

其中,m為輸出信號的時間索引,表示輸出信號的不同時刻;l為卷積和的中間索引,用于遍歷輸入信號或者IRF的各個值;pi(m)、hypei(m)、ye(m)分別為離散時間激勵、離散單位脈沖響應函數和離散系統響應,可以將式(26)展開成矩陣形式,令pi(m)=pim,hypei(m)=heim,ye(m)=yem,其余類似,并假設hypei(l)=0(l≥m),則式(26)的部分數據矩陣形式可寫為

3 數值算例

為驗證基于AKF的時域TPA方法的準確性與適用性,分別對懸臂薄板模型施加隨機激勵、沖擊激勵和正弦激勵,明確新方法的應用流程,并與傳統時域TPA方法進行對比分析。隨后,將此新方法運用于某車輛模型,以此證明該方法在處理復雜結構問題中的有效性。

3.1 懸臂薄板模型

如圖3所示,懸臂薄板模型的長、寬和厚度分別為768mm×368mm×8mm,其材料參數設置如下:彈性模量E=71 GPa,泊松比ν=0.31,材料密度ρ=2700 kg/m3,模態阻尼比ξj為0.1(1~20 Hz),0.04(gt;20 Hz)。采用4416個四邊形網格對模型進行劃分,網格尺寸為8 mm×8 mm。令模型的一側完全約束,讓其他側處于自由狀態。

在懸臂薄板模型上選取3個激勵點,沿Z向加載隨機激勵,方案1:任取6個指示點(6925、7269、7495、7555、8580、9431),方案2:基于GA選取6個指示點(5054、6228、8555、8648、9290、9254)。如圖4所示,任取指示點可能會導致AKF識別載荷出現波動,

而基于GA選取指示點可以改善這一情況,提高載荷識別精度。

為進一步驗證所提出的基于AKF的時域TPA方法的優越性,在均使用GA選取指示點的前提下,利用Moore-Penrose偽逆方法、Tikhonov正則化方法和奇異值截斷方法[16]對指示點頻響函數矩陣求逆,并分別采用快速逆傅里葉變換法和頻率采樣法獲得單位脈沖響應函數,構造了6種傳統時域TPA方法對目標點進行時域傳遞路徑分析并和所提方法進行比較。

首先,沿激勵點Z向分別加載隨機激勵、正弦激勵和沖擊激勵,將不同方法識別的時域工況載荷與真實載荷進行了比較,如圖5、圖6和圖7所示。圖中,True LOAD表示真實載荷,AKFLOAD表示增廣卡爾曼濾波器法,Pinv表示Moore-Penrose偽逆法,Tikh表示Tikhonov正則化法,Tsvd表示奇異值截斷法,IFFT表示快速逆傅里葉變換法,FSM表示頻率采樣法。則Pinv IFFT的含義為:首先采用Moore-Penrose偽逆法對指示點頻響函數矩陣求逆,然后利用快速逆傅里葉變換法構造去卷積濾波器,其他表示同上。結果表明,傳統時域TPA方法的去卷積濾波器僅在激勵點2沖擊載荷的識別上一度超過了AKF,而在其余激勵點處各時域工況載荷的識別上普遍存在精度較低的問題;AKF減小了載荷識別誤差,提高了幾乎所有激勵點處各時域工況載荷的識別精度。由于采用GA確定了最佳指示點位置,并且選擇加速度和位移(應變)作為指示點響應信息,因此確保了動態力和系統狀態的可識別性、穩定性和唯一性。

利用頻率采樣法、快速逆傅里葉變換法和提出的最小二乘算法分別求解激勵點到目標點的單位脈沖響應函數并與實際單位脈沖響應函數進行比較,結果在圖8中給出。圖中,TRUE表示真

實單位脈沖響應函數,LS表示最小二乘算法,FSM表示頻率采樣法,IFFT表示快速逆傅里葉變換法。可知,在不考慮力錘激勵誤差和周圍環境噪聲的前提下,雖然最小二乘算法的計算量更大,但其求解精度高于傳統方法。

在得到各激勵點的時域工況載荷和對應的單位脈沖響應函數后,分別采用所提方法和傳統方法計算各傳遞路徑的時域貢獻量并與真實值比較,如圖9、圖10和圖11所示。結果表明,所提方法除了在初始時刻估計精度比較低外,其他時刻估計精度要高于傳統方法。這主要是因為AKF的初始參數是隨機設置的,隨著時間更新會逐漸調整至最佳,從而實現最優估計。

3.2 車輛模型

在實際工程環境中,車輛承受各種隨機動載荷,其中路面激勵是影響汽車NVH的主要因素之一,由于輪胎的復雜特性,在實踐中往往提取輪心載荷代替輪胎進行計算分析。本節利用所提出的時域TPA方法獲取輪心瞬態激勵并分析車內振動以驗證所提方法的有效性。

某車輛模型如圖12所示,該模型主要由發動機、底盤和車身等部分組成,其材料密度、彈性模量和泊松比等參數參照有關標準定義。

需要強調的是,由于4個輪心共計有12個載荷(沿X、Y、Z方向),為避免占用過多篇幅,

這里只分析了輪心的Z向載荷及其振動傳遞情況。將AKF識別的時域載荷與傳統方法進行比較,如圖13所示。結果表明,6種傳統方法幾乎無法識別輪心Z向的隨機載荷,而對于AKF,即使是整車這種復雜結構,其估計結果除了在部分時刻有些許波動外,大部分時刻能較準確地反映隨機載荷的真實值。

將基于最小二乘算法辨識的輪心Z向至座椅Z向振動單位脈沖響應函數和傳統方法進行比較,如圖14所示。結果表明,頻率采樣法、快速逆傅里葉變換法幾乎無法轉換提取得到所需的IRF,雖然在復雜結構最小二乘算法不可避免地引入了些許誤差,但其曲線擬合度仍可以滿足要求。

將AKF估計的時域工況載荷和LS辨識的單位脈沖響應函數進行線性卷積計算各傳遞路徑的時域貢獻量并與傳統方法進行比較,如圖15所示。結果表明,傳統方法僅能勉強反映輪心Z向載荷對目標點振動貢獻量的變化趨勢,而基于AKF的時域傳遞路徑分析方法雖然在一些時刻效果不理想,但整體上可以較準確地反映各路徑對目標點的時域貢獻量。

4 結論

針對傳統時域傳遞路徑分析(TPA)方法分析精度較低的問題,本文提出了一種新的時域TPA方法,從理論上對該方法及其求解思路進行了闡明,并用仿真驗證了該方法的有效性,得出了以下結論。

(1)該方法采用遺傳算法確定指示點位置,增廣卡爾曼濾波器估計時域工況載荷,識別精度較高。算例表明,即使均采用遺傳算法,該方法對各種類型載荷的識別精度也要高于傳統時域TPA方法的去卷積濾波器。

(2)該方法采用最小二乘法辨識單位脈沖響應函數,算例表明,在不考慮錘擊誤差和環境噪聲的前提下,精度要高于快速逆傅里葉變換法和頻率采樣法。各路徑時域貢獻量僅在一些時刻分析精度較低,其余時刻和真實值較為吻合。

參考文獻:

[1]van der SEIJS M V, de KLERK D, RIXEN D J. General Framework for Transfer Path Analysis:History, Theory and Classification of Techniques[J]. Mechanical Systems and Signal Processing, 2016, 68:217-244.

[2]YE Shaogan, HOU Liang, ZHANG Pandeng, et al. Transfer Path Analysis and Its Application in Low-frequency Vibration Reduction of Steering Wheel of a Passenger Vehicle[J]. Applied Acoustics, 2020, 157:107021.

[3]王少華, 譚博歡, 張邦基, 等. 純電動客車振動試驗分析及動力總成隔振優化[J]. 振動與沖擊, 2021, 40(1):226-232.

WANG Shaohua, TAN Bohuan, ZHANG Bangji, et al. Vibration Test Analysis and Powertrain Vibration Isolation Optimization of Pure Electric Bus[J]. Journal of Vibration and Shock, 2021, 40(1):226-232.

[4]de SITTER G, DEVRIENDT C, GUILLAUME P, et al. Operational Transfer Path Analysis[J]. Mechanical Systems and Signal Processing, 2010, 24(2):416-431.

[5]SHIN T, KIM Y S, AN K, et al. Transfer Path Analysis of Rumbling Noise in a Passenger Car Based on In-situ Blocked Force Measurement[J]. Applied Acoustics, 2019, 149:1-14.

[6]郝鵬, 鄭四發, 連小珉. 運動噪聲源的時域傳遞路徑模型及貢獻率分析[J]. 機械工程學報, 2012, 48(8):104-109.

HAO Peng, ZHENG Sifa, LIAN Xiaomin. Time-domain Transfer Path Modal and Contribution Analysis of Moving Noise Sources[J]. Journal of Mechanical Engineering, 2012, 48(8):104-109.

[7]褚志剛, 熊敏, 楊洋, 等. 車內噪聲時域傳遞路徑分析[J]. 振動與沖擊, 2015, 34(17):161-166.

CHU Zhigang, XIONG Min, YANG Yang, et al. Time-domain Transfer Path Analysis of Automobile Interior Noise[J]. Journal of Vibration and Shock, 2015, 34(17):161-166.

[8]謝小平, 王洪波, 梁煬烊. 瞬態聲品質時域傳遞路徑方法與應用研究[J]. 振動與沖擊, 2020, 39(7):253-259.

XIE Xiaoping, WANG Hongbo, LIANG Yangyang. Analysis of Transient Sound Quality in Time Domain Transfer Path Method and Its Application[J]. Journal of Vibration and Shock, 2020, 39(7):253-259.

[9]呂來, 唐中華, 張志飛, 等. 時域工況傳遞路徑分析方法[J]. 振動與沖擊, 2023, 42(15):92-100.

LYU Lai, TANG Zhonghua, ZHANG Zhifei, et al. Analysis Method of Time Domain Operating Conditions Transfer Path[J]. Journal of Vibration and Shock, 2023, 42(15):92-100.

[10]LOURENS E, REYNDERS E, de ROECK G, et al. An Augmented Kalman Filter for Force Identification in Structural Dynamics[J]. Mechanical Systems and Signal Processing, 2012, 27:446-460.

[11]CUMBO R, TAMAROZZI T, JANSSENS K, et al. Kalman-based Load Identification and Full-field Estimation Analysis on Industrial Test Case[J]. Mechanical Systems and Signal Processing, 2019, 117:771-785.

[12]NAETS F, CUADRADO J, DESMET W. Stable Force Identification in Structural Dynamics Using Kalman Filtering and Dummy-measurements[J]. Mechanical Systems and Signal Processing, 2015, 50:235-248.

[13]CUMBO R, MAZZANTI L, TAMAROZZI T, et al. Advanced Optimal Sensor Placement for Kalman-based Multiple-input Estimation[J]. Mechanical Systems and Signal Processing, 2021, 160:107830.

[14]ZOU Donglin, ZHAO Han, LIU Gaoyu, et al. Application of Augmented Kalman Filter to Identify Unbalance Load of Rotor-bearing System:Theory and Experiment[J]. Journal of Sound and Vibration, 2019, 463:114972.

[15]MAES K, LOURENS E, van NIMMEN K, et al. Design of Sensor Networks for Instantaneous Inversion of Modally Reduced Order Models in Structural Dynamics[J]. Mechanical Systems and Signal Processing, 2015, 52:628-644.

[16]彭博, 鄭四發, 廖祥凝, 等. 室內聲場重現中的時域去卷積方法[J]. 振動與沖擊, 2016, 35(4):115-120.

PENG Bo, ZHENG Sifa, LIAO Xiangning, et al. Time Domain Deconvolution for Reproducing Sound Field in an Ordinary Room[J]. Journal of Vibration and Shock, 2016, 35(4):115-120.

(編輯 王旻玥)

作者簡介:

朱 雨,男,1998年生,碩士研究生。研究方向為汽車NVH。E-mail:ncu_zhuyu@163.com。

何智成(通信作者),男,1983年生,教授、博士研究生導師。研究方向為先進結構與智能設計,智能汽車與智能控制。E-mail:hezhicheng815@163.com。

收稿日期:2023-11-12

基金項目:湖南省杰出青年基金(2021JJ10016);廣西科技重大專項(桂科AA22068108);柳州市科技計劃(2022AAA0101,2021AAA0111)

主站蜘蛛池模板: 亚洲手机在线| 国产精品成人免费综合| 漂亮人妻被中出中文字幕久久| 亚洲人成网站观看在线观看| 亚洲日韩精品无码专区97| 国产欧美性爱网| 澳门av无码| 看你懂的巨臀中文字幕一区二区| 伊人福利视频| 91啦中文字幕| 一级毛片a女人刺激视频免费| 国产伦片中文免费观看| 国产理论最新国产精品视频| 亚洲av片在线免费观看| 亚洲国产综合自在线另类| 成人精品亚洲| 午夜精品国产自在| 国产91色在线| 亚洲系列无码专区偷窥无码| 99成人在线观看| 内射人妻无套中出无码| 亚洲日韩久久综合中文字幕| 午夜福利无码一区二区| 国产精品爽爽va在线无码观看| 国产精品美女自慰喷水| 精品無碼一區在線觀看 | 色国产视频| 国产精品无码影视久久久久久久| 91精品免费久久久| 精品视频一区在线观看| 日韩a在线观看免费观看| 中文无码伦av中文字幕| 国产黑丝一区| 国产精品亚洲αv天堂无码| 欧美午夜久久| 国产精品无码制服丝袜| 无码福利视频| 国产精品久久久免费视频| 国产精品开放后亚洲| 欧美色图久久| 国产精品丝袜在线| 国产精品亚欧美一区二区| 久久精品国产91久久综合麻豆自制| 华人在线亚洲欧美精品| 熟女日韩精品2区| 9久久伊人精品综合| 国产麻豆精品在线观看| 日韩东京热无码人妻| 最新日韩AV网址在线观看| 国产精品亚洲а∨天堂免下载| 国产成人精彩在线视频50| 99偷拍视频精品一区二区| 午夜电影在线观看国产1区| 亚洲天堂精品在线观看| 日韩精品资源| 特级精品毛片免费观看| 久久精品人人做人人爽电影蜜月| 蜜臀av性久久久久蜜臀aⅴ麻豆| 午夜福利无码一区二区| 国产精品白浆无码流出在线看| 国产成人亚洲毛片| 中文字幕无码制服中字| 国产超碰一区二区三区| 国产亚洲精品自在久久不卡| 亚洲小视频网站| 久久黄色影院| 天堂av综合网| 欧美福利在线观看| 日韩欧美国产精品| 视频二区国产精品职场同事| 无码内射在线| 色老二精品视频在线观看| 中国成人在线视频| 国产成人毛片| 在线观看亚洲国产| 欧美日韩国产在线人成app| 日韩国产 在线| 国产AV无码专区亚洲精品网站| 成人国产精品2021| 国产精欧美一区二区三区| 天天综合网站| 亚洲天堂在线视频|