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

基于MIC特征提取與BO-CatBoost的航空發(fā)動機(jī)RUL預(yù)測

2024-02-23 09:16:04李東君李東文朱貴富
空軍工程大學(xué)學(xué)報 2024年1期
關(guān)鍵詞:發(fā)動機(jī)方法模型

李東君, 李 亞, 李東文, 朱貴富

(1.昆明理工大學(xué)信息工程與自動化學(xué)院,昆明,650504;2.昆明理工大學(xué)信息化建設(shè)管理中心,昆明,650504)

航空發(fā)動機(jī)作為飛機(jī)的核心部件之一,內(nèi)部結(jié)構(gòu)復(fù)雜,在運行期間,其健康狀況易受環(huán)境及內(nèi)部構(gòu)造等多種因素的影響。相關(guān)研究表明,在飛機(jī)因機(jī)械故障引起的飛行事故中,發(fā)動機(jī)故障導(dǎo)致的飛行事故數(shù)量占比最高[1]。飛機(jī)發(fā)動機(jī)一旦發(fā)生故障或失效,將會嚴(yán)重威脅人們的生命健康及財產(chǎn)安全。為提升發(fā)動機(jī)運行的安全性、可靠性及穩(wěn)定性,并有效保障人們的生命及財產(chǎn)安全,對發(fā)動機(jī)開展視情維修與健康監(jiān)測具有重要意義。

故障預(yù)測和健康管理(prognostics and health management,PHM)技術(shù)主要利用機(jī)械設(shè)備在運行期間產(chǎn)生的各類數(shù)據(jù),運用數(shù)據(jù)處理與數(shù)據(jù)分析等手段,實現(xiàn)對復(fù)雜設(shè)備的健康狀態(tài)檢測、預(yù)測管理及維修策略支持,有效降低故障發(fā)生率和維修成本[2]。航空發(fā)動機(jī)剩余使用壽命(remaining useful life,RUL)維護(hù)預(yù)測是PHM的關(guān)鍵任務(wù)之一。RUL預(yù)測的主要目標(biāo)是根據(jù)機(jī)械設(shè)備健康狀態(tài)的有效信息,對設(shè)備到達(dá)安全操作限制的剩余壽命時間進(jìn)行估計[3]。由于航空發(fā)動機(jī)長期工作在條件惡劣的環(huán)境中,為保證發(fā)動機(jī)運行的安全性,需要定期對發(fā)動機(jī)進(jìn)行維護(hù),如何有效利用發(fā)動機(jī)的歷史退化信息對發(fā)動機(jī)的健康狀態(tài)進(jìn)行RUL預(yù)測是當(dāng)下研究的重點。

目前,因預(yù)測原理的差異,國內(nèi)外關(guān)于航空發(fā)動機(jī)的RUL預(yù)測方法主要分為2類:基于物理失效模型的RUL預(yù)測方法和基于數(shù)據(jù)驅(qū)動的RUL預(yù)測方法[4]。雖然基于物理失效模型的RUL預(yù)測方法具有良好的預(yù)測精度,但由于航空發(fā)動機(jī)內(nèi)部系統(tǒng)構(gòu)造極為復(fù)雜,在對部件進(jìn)行退化過程建模時需要掌握部件原理的專業(yè)知識,使得預(yù)測難度增大且普適性較差。基于數(shù)據(jù)驅(qū)動的RUL預(yù)測方法無需建立復(fù)雜的數(shù)學(xué)模型,它借助機(jī)器學(xué)習(xí)手段挖掘和分析退化數(shù)據(jù)隱藏的內(nèi)部信息來獲得數(shù)據(jù)間內(nèi)在的關(guān)聯(lián)特征,且預(yù)測的性能和準(zhǔn)確率也較高[5-6]。

隨著人工智能技術(shù)及智能算法的飛速發(fā)展,基于數(shù)據(jù)驅(qū)動的RUL預(yù)測方法逐漸成為國內(nèi)外學(xué)者研究的熱點。文獻(xiàn)[7]對歷史監(jiān)測數(shù)據(jù)進(jìn)行融合,將貝葉斯理論(bayesian)應(yīng)用到航空發(fā)動機(jī)RUL預(yù)測中。文獻(xiàn)[8]利用最大信息系數(shù)(maximal information coefficient,MIC)對短期負(fù)荷數(shù)據(jù)進(jìn)行特征選擇,將篩選的特征輸入構(gòu)建的模型中進(jìn)行預(yù)測。文獻(xiàn)[9]提出了一種利用皮爾遜相關(guān)系數(shù)方法對發(fā)動機(jī)數(shù)據(jù)子集進(jìn)行協(xié)變量篩選獲得最優(yōu)協(xié)變量表達(dá)式。文獻(xiàn)[10]采用類別特征梯度提升算法(categorical boosting,CatBoost)構(gòu)建孔隙壓力預(yù)測模型,有效提升了預(yù)測精度。文獻(xiàn)[11]提出了一種在聯(lián)邦學(xué)習(xí)框架中使用全局健康退化表示(a global health degradation representation ,GHDR)的新的RUL預(yù)測方法。

對于復(fù)雜機(jī)械設(shè)備的RUL預(yù)測問題,雖然現(xiàn)有文獻(xiàn)所提方法在航空發(fā)動機(jī)剩余壽命方面已經(jīng)取的一定的成績,但綜合考慮發(fā)動機(jī)自身復(fù)雜的物理結(jié)構(gòu)及運行環(huán)境的惡劣條件等因素,對于發(fā)動機(jī)的RUL預(yù)測還存在一些待解決的:發(fā)動機(jī)傳感器由于受外部環(huán)境因素的影響,使監(jiān)測的數(shù)據(jù)受到噪聲干擾,無法有效去除環(huán)境噪聲;針對航空發(fā)動機(jī)監(jiān)測參數(shù)非線性特點,現(xiàn)有方法無法充分提取對發(fā)動機(jī)運行周期影響較大的關(guān)鍵特征,在RUL預(yù)測準(zhǔn)確度上仍有提升空間。

針對以上問題,本文提出了一種基于MIC特征提取與貝葉斯優(yōu)化類別特征梯度提升(bayesian optimization categorical boosting,BO CatBoost)相結(jié)合的發(fā)動機(jī)RUL預(yù)測方法。考慮到不同監(jiān)測特征對發(fā)動機(jī)壽命影響的差異,引入MIC算法計算各個特征對發(fā)動機(jī)RUL的相關(guān)性,篩選出相關(guān)性較大的特征作為CatBoost預(yù)測模型的輸入,引入貝葉斯優(yōu)化(bayesian optimization,BO)算法對CatBoost預(yù)測模型中的超參數(shù)進(jìn)行訓(xùn)練和調(diào)優(yōu),得到RUL預(yù)測結(jié)果的最優(yōu)值。利用評價指標(biāo)均方根誤差(root mean square error,RMSE)、判定系數(shù)(coefficient of determination,R2)和平均絕對誤差(mean absolute error,MAE)對預(yù)測模型進(jìn)行性能分析和評估,有效驗證了本文所提方法的可行性。

1 RUL預(yù)測方法

1.1 最大信息系數(shù)

最大信息系數(shù)MIC以互信息(mutual information,MI)為基礎(chǔ),是一種用于衡量2個變量之間線性或非線性相關(guān)性強(qiáng)弱的算法[12]。設(shè)X={x1,x2,…,xn}與Y={y1,y2,…,yn}分別為數(shù)據(jù)集中的隨機(jī)變量,n為樣本數(shù)量,則X與Y之間的MI為:

(1)

式中:p(x,y)為X與Y之間的聯(lián)合概率密度;p(x)與p(y)分別為X與Y之間的邊緣概率密度。

MIC克服了MI在計算連續(xù)變量的聯(lián)合概率密度函數(shù)困難的缺陷,能最大程度地找到兩變量之間的相關(guān)性[13]。MIC的計算公式為:

式中:B為樣本數(shù)量;N為樣本變量;I(x;y)為x與y之間的MI。2個變量間的MIC值越接近1,則其相關(guān)性越強(qiáng),MIC∈[0,1]。

1.2 CatBoost算法

類別特征梯度提升算法(CatBoost)是在梯度提升決策樹(gradient boosting decision tree,GBDT)算法的基礎(chǔ)上改進(jìn)的算法[14]。CatBoost算法是由類別特征和梯度提升組成的一種高準(zhǔn)確性梯度提升框架。該算法以對稱決策樹作為基學(xué)習(xí)器,解決了梯度偏差及預(yù)測偏移問題,有效防止模型的過度擬合。該算法將數(shù)據(jù)的類別特征進(jìn)行編碼,每層分裂時,設(shè)置分裂閾值,并將所有的類別特征與指定的特征進(jìn)行組合,參與下一層分裂。分裂結(jié)束后,使用梯度無偏來估計預(yù)測偏移,尋找最優(yōu)目標(biāo)。

建立分類特征樹的過程中,因需要考慮監(jiān)測數(shù)據(jù)之間的相關(guān)性,本文引入MIC算法來計算監(jiān)測參數(shù)之間的非線性關(guān)聯(lián)性。MIC能有效避免互信息在計算連續(xù)變量的聯(lián)合概率密度函數(shù)困難的問題,且最大程度地挖掘2個變量之間的非線性相關(guān)性。

本文使用CatBoost回歸算法來解決發(fā)動機(jī)RUL的預(yù)測問題。設(shè)原始數(shù)據(jù)集為|D|={(x1,y1),(x2,y2),…,(xn,yn)},則σ=(σ1,σ2,…,σn)為|D|經(jīng)過重新排序后的序列狀態(tài)。

(3)

式中:xσj,k為數(shù)據(jù)集σj的第k個特征;p為先驗概率,用來減少噪聲數(shù)據(jù)的干擾;?為大于0時的權(quán)重系數(shù)值,用于調(diào)節(jié)p的影響程度。

1.3 貝葉斯優(yōu)化算法

貝葉斯優(yōu)化(BO)算法是一種基于搜索函數(shù)和高斯過程(gaussian processes,GP)的參數(shù)更新優(yōu)化算法,根據(jù)定義的目標(biāo)函數(shù)迭代評估新的參數(shù)[15]。利用高斯過程進(jìn)行調(diào)參,設(shè)置最優(yōu)參數(shù)后,不斷迭代更新先驗值和模型參數(shù),直到找到最優(yōu)的超參數(shù)組合。本文使用貝葉斯算法中的GP及采集函數(shù)(acquisition function,AF)優(yōu)化CatBoost中的超參數(shù)。GP回歸表達(dá)式:

f(x)~GP(m(x),K(x,x))

(4)

式中:K(x,x)為協(xié)方差矩陣;m(x)為均值向量函數(shù)。

采集函數(shù)對樣本的候選值進(jìn)行評估后得到最優(yōu)解。其計算公式為:

式中:σ(x)為GP的方差;u(x)為樣本均值。

φ(*)為標(biāo)準(zhǔn)正態(tài)分布的累積分布函數(shù),α為超參數(shù)。貝葉斯優(yōu)化CatBoost算法超參數(shù)流程如圖1所示。

圖1 貝葉斯優(yōu)化CatBoost算法超參數(shù)流程圖

1.4 評價指標(biāo)

為更好地驗證本文所提RUL預(yù)測模型的準(zhǔn)確性,選取實驗中2個常用的模型性能評價指標(biāo),即判定系數(shù)R2[16]、均方根誤差RMSE[17-19]與平均絕對誤差MAE[20]來衡量模型的預(yù)測性能。

RMSE通常用于度量模型的預(yù)測值和真實值之間的總體偏差,其值越小,則預(yù)測性能精度越高。RMSE的計算式為:

(6)

R2用于反映由樣本回歸線做出解釋的離差平方和中的比重,其值越接近1,說明模型擬合度越好,其表達(dá)式如下:

(7)

MAE表示絕對誤差的平均值,其計算式為:

(8)

2 基于MIC與BO-CatBoost的預(yù)測模型框架

本文所提出的航空發(fā)動機(jī)RUL預(yù)測模型框架總體介紹:利用MIC算法篩選出能夠表征發(fā)動機(jī)退化性能的傳感器監(jiān)測參數(shù),為降低監(jiān)測參數(shù)量綱對預(yù)測模型的影響,將監(jiān)測參數(shù)歸一化在[0,1]之間。構(gòu)建基于貝葉斯超參數(shù)尋優(yōu)的CatBoost預(yù)測模型進(jìn)行訓(xùn)練,將測試集輸入訓(xùn)練結(jié)束的預(yù)測模型中預(yù)測發(fā)動機(jī)的RUL。并通過評價指標(biāo)RMSE、R2、MAE來評價預(yù)測模型的性能。

1)數(shù)據(jù)預(yù)處理

在數(shù)據(jù)預(yù)處理階段,利用MIC算法分析各個監(jiān)測參數(shù)與發(fā)動機(jī)運行壽命間的相關(guān)性強(qiáng)弱,篩選出相關(guān)性較強(qiáng)的監(jiān)測參數(shù)用于實驗驗證;為降低監(jiān)測參數(shù)量綱對預(yù)測模型的影響,對監(jiān)測參數(shù)進(jìn)行歸一化處理。

2)CatBoost模型訓(xùn)練

構(gòu)建基于CatBoost算法的發(fā)動機(jī)RUL預(yù)測訓(xùn)練模型,設(shè)置超參數(shù)及其尋優(yōu)范圍;利用BO算法對CatBoost訓(xùn)練模型中的超參數(shù)進(jìn)行優(yōu)化,尋優(yōu)過程中返回RMSE最小值及尋優(yōu)參數(shù)的取值。利用篩選出的最優(yōu)超參數(shù)作為預(yù)測模型的最終超參數(shù)組合。

3)剩余壽命預(yù)測

將經(jīng)過數(shù)據(jù)預(yù)處理后的監(jiān)測參數(shù)輸入到構(gòu)建的BO-CatBoost預(yù)測模型中進(jìn)行訓(xùn)練,返回預(yù)測評價結(jié)果RMSE,R2及MAE的值,實現(xiàn)航空發(fā)動機(jī)的RUL預(yù)測。

3 實驗驗證

3.1 數(shù)據(jù)集介紹

本文用于實驗的數(shù)據(jù)是來自美國國家航空航天局(national aeronautics and space administration,NASA)發(fā)布的渦扇發(fā)動機(jī)商用模塊化航空推進(jìn)系統(tǒng)仿真數(shù)據(jù)集(commercial modular aero-propulsion system simulation,C-MAPSS)[21]。該數(shù)據(jù)集中包含4個子數(shù)據(jù)集FD001~FD004,每個子數(shù)據(jù)集都由訓(xùn)練集與測試集組成。訓(xùn)練集中包含了發(fā)動機(jī)從初始運行狀態(tài)到磨損失效后的全壽命周期數(shù)據(jù),測試集僅包含發(fā)動機(jī)在故障發(fā)生前的部分運行周期數(shù)據(jù)[22]。每個子數(shù)據(jù)集中包含飛行高度、馬赫數(shù)與油門桿解算器角度3種操作條件,運行周期及21種航空發(fā)動機(jī)傳感器監(jiān)測參數(shù)。本文選取FD001~FD004共4個子數(shù)據(jù)集進(jìn)行發(fā)動機(jī)RUL實驗。C-MAPSS數(shù)據(jù)集見表1。

表1 C-MAPSS數(shù)據(jù)集

圖2展示了FD001數(shù)據(jù)集中100臺發(fā)動機(jī)的最大運行周期分布情況。從圖2可知,在FD001數(shù)據(jù)集100臺發(fā)動機(jī)中最小的運行周期為128,最大的運行周期為362,其他發(fā)動機(jī)的運行周期大部分分布在[145,250]范圍內(nèi)。

圖2 FD001數(shù)據(jù)集100臺發(fā)動機(jī)的最大運行周期分布情況

3.2 數(shù)據(jù)預(yù)處理

3.2.1 特征提取

不同的監(jiān)測參數(shù)對構(gòu)建的訓(xùn)練模型具有不同程度的影響,為有效提升發(fā)動機(jī)RUL預(yù)測的精度,采用基于MIC的特征選擇方法對影響航空發(fā)動機(jī)壽命周期的監(jiān)測參數(shù)進(jìn)行篩選。表2為FD001不同監(jiān)測參數(shù)與發(fā)動機(jī)RUL的MIC計算結(jié)果。由表2可知,編號為1、5、10、16、18與19的監(jiān)測參數(shù)MIC值為0,且編號6的值接近于0,說明這7個監(jiān)測參數(shù)對航空發(fā)動機(jī)的壽命運行周期相關(guān)性極弱,將這7種監(jiān)測參數(shù)從21種監(jiān)測參數(shù)中剔除,利用剩余的14種監(jiān)測進(jìn)行航空發(fā)動機(jī)的RUL預(yù)測實驗。

3.2.2 數(shù)據(jù)預(yù)處理

因航空發(fā)動機(jī)不同的監(jiān)測參數(shù)具有不同的量綱,為了縮小監(jiān)測參數(shù)數(shù)值之間的差異,提高預(yù)測的效率及準(zhǔn)確率,本文選取最小最大歸一化公式對發(fā)動機(jī)監(jiān)測參數(shù)進(jìn)行歸一化處理。其計算公式如下[23]:

3.3 基于貝葉斯優(yōu)化的CatBoost模型

因CatBoost訓(xùn)練模型中的超參數(shù)取值范圍不同,會對發(fā)動機(jī)的RUL預(yù)測值帶來不同程度的影響。本文使用BO算法對CatBoost回歸預(yù)測模型中的超參數(shù)進(jìn)行尋優(yōu)。搭建BO-CatBoost訓(xùn)練模型并設(shè)置超參數(shù)取值范圍。表3為BO-CatBoost超參數(shù)尋優(yōu)結(jié)果。

表3 BO-CatBoost超參數(shù)尋優(yōu)結(jié)果

為獲得最佳的超參數(shù)值,本文分別訓(xùn)練預(yù)測模型中的超參數(shù)iterations與learning_rate在不同取值組合下的評估指標(biāo)RMSE值。圖3為預(yù)測模型不同超參數(shù)值的RMSE。由圖3可知,當(dāng)iterations為999,learning_rate為0.07時,模型訓(xùn)練效果最好,RMSE值最小。

圖3 預(yù)測模型不同超參數(shù)值的RMSE

3.4 實驗預(yù)測結(jié)果

圖4展示了FD001~FD004數(shù)據(jù)集下的航空發(fā)動機(jī)RUL預(yù)測結(jié)果。圖4中,紅色曲線為發(fā)動機(jī)的RUL真實值,藍(lán)色曲線為RUL預(yù)測值。由于FD001與FD003為在單一操作條件和故障模式下采集的數(shù)據(jù),發(fā)動機(jī)數(shù)量少,預(yù)測結(jié)果較為稀疏。FD002與FD004為在多工況環(huán)境中采集的數(shù)據(jù),發(fā)動機(jī)數(shù)量較多,與其他2個數(shù)據(jù)集的RUL預(yù)測相比,預(yù)測難度具有一定挑戰(zhàn)性。由對比結(jié)果可知,發(fā)動機(jī)的RUL真實值與RUL預(yù)測值比較貼合,預(yù)測誤差較小,說明本文所提模型的預(yù)測效果較好。

(a)FD001測試集

(b)FD002測試集

(c)FD003測試集

(d)FD004測試集

3.5 比較分析

為全面評估本文所提算法的性能,實驗過程中選擇RMSE、MAE及R23個性能評價指標(biāo)來衡量各個預(yù)測模型的性能。表4展示了不同預(yù)測模型在FD001~FD004數(shù)據(jù)集中的RUL預(yù)測結(jié)果。由表4可知,與其他預(yù)測方法相比,采用MIC-BO-CatBoost方法進(jìn)行預(yù)測時,評價指標(biāo)RMSE與MAE的值最小,R2也最接近于1。

圖5展示了MIC-BO-CatBoost預(yù)測模型在FD001~FD004數(shù)據(jù)集中的預(yù)測誤差箱線圖。箱子中間的實線表示RUL預(yù)測的期望,箱子的規(guī)模越小,表示模型在預(yù)測過程中的RUL預(yù)測結(jié)果的不確定性越低,準(zhǔn)確率越高。由圖5可知,本文所提預(yù)測模型MIC-BO-CatBoost在FD001與FD003數(shù)據(jù)集上的箱線圖規(guī)模較小,預(yù)測精確度較高。

表4 不同預(yù)測模型在FD001~FD004數(shù)據(jù)集中的RUL預(yù)測結(jié)果

圖5 MIC-BO-CatBoost預(yù)測模型在FD001~FD004數(shù)據(jù)集中的預(yù)測誤差箱線圖

為了對比不同模型在FD001~FD004數(shù)據(jù)集中的預(yù)測效果,本文構(gòu)建嶺(Ridge)回歸[24]、K近鄰回歸(KNN)算法[25]、極端梯度提升(XGBoost)算法[26]與本文所提方法進(jìn)行比較分析。分別選取FD001~FD004數(shù)據(jù)集中的第32號、3號、76號及56號發(fā)動機(jī)進(jìn)行RUL預(yù)測。圖6展示了FD001~FD004在不同預(yù)測模型下的RUL預(yù)測結(jié)果。灰色區(qū)間代表RUL預(yù)測分布的95%置信區(qū)間。由圖6可知,其他預(yù)測模型的RUL預(yù)測結(jié)果大部分位于RUL真實值之上或之下,說明預(yù)測結(jié)果表現(xiàn)出滯后或超前預(yù)測,而本文所提方法的RUL預(yù)測值緊密圍繞RUL真實值波動,擬合程度較高。與其他3種預(yù)測模型的對比結(jié)果可知,本文所提方法的RUL預(yù)測值與真實值的偏差最小,且預(yù)測結(jié)果基本都被95%置信區(qū)間覆蓋,有效說明本文所提方法能夠更好地用于航空發(fā)動機(jī)的剩余使用壽命預(yù)測。

(a)測試集FD001第32號發(fā)動機(jī)

(b)測試集FD002第3號發(fā)動機(jī)

(c)測試集FD003第76號發(fā)動機(jī)

(d)測試集FD004第56號發(fā)動機(jī)

4 結(jié)論

本文針對航空發(fā)動機(jī)非線性特征提取困難及預(yù)測精度不高等問題,提出了一種基于MIC-BO-CatBoost的航空發(fā)動機(jī)預(yù)測模型,使用C-MPASS發(fā)動機(jī)退化數(shù)據(jù)集進(jìn)行驗證和分析,得到以下結(jié)論:

1)通過MIC方法分析各個監(jiān)測參數(shù)對發(fā)動機(jī)壽命運行周期的相關(guān)性,提取能表征發(fā)動機(jī)退化過程的非線性特征。

2)利用CatBoost算法對RUL預(yù)測模型進(jìn)行訓(xùn)練,有效解決了預(yù)測中的梯度偏差及預(yù)測偏移問題。

3)使用BO-CatBoost算法建立預(yù)測模型,通過在FD001~FD004數(shù)據(jù)集上進(jìn)行實驗評估后,結(jié)果表明與其他預(yù)測模型相比,本文所提方法的預(yù)測精度更加貼近發(fā)動機(jī)的RUL真實值,證明了MIC-BO-CatBoost預(yù)測模型的有效性。

猜你喜歡
發(fā)動機(jī)方法模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
發(fā)動機(jī)空中起動包線擴(kuò)展試飛組織與實施
3D打印中的模型分割與打包
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
捕魚
新一代MTU2000發(fā)動機(jī)系列
新型1.5L-Eco-Boost發(fā)動機(jī)
主站蜘蛛池模板: 成年人视频一区二区| 黄色网页在线播放| 午夜免费小视频| 日韩精品成人在线| 99无码中文字幕视频| 欧美不卡视频一区发布| 欧美日韩精品一区二区在线线 | 国产另类视频| 免费人成在线观看视频色| 性网站在线观看| 婷婷色丁香综合激情| 粗大猛烈进出高潮视频无码| 午夜日b视频| 玖玖精品视频在线观看| 日本a级免费| 特级精品毛片免费观看| 国产精品嫩草影院视频| 狼友视频一区二区三区| 色悠久久综合| 男女男免费视频网站国产| 久久精品无码一区二区日韩免费| 狠狠干综合| 91精品综合| 亚洲精品无码高潮喷水A| 在线国产资源| 亚洲国产清纯| 国产美女在线免费观看| 亚洲欧美不卡视频| 国产理论一区| 亚洲综合第一区| 亚洲中文无码av永久伊人| 亚洲男人的天堂视频| 久久综合色视频| 国产成人亚洲精品无码电影| 久久无码免费束人妻| 激情爆乳一区二区| a亚洲天堂| 欧美区一区| 久久无码免费束人妻| 午夜福利免费视频| 色哟哟国产精品| 国产精品久久久久久影院| 日韩毛片免费视频| 高清码无在线看| 四虎国产成人免费观看| 日韩av高清无码一区二区三区| 欧美亚洲国产一区| 欧美精品v欧洲精品| 欧美精品v| 精品少妇人妻av无码久久| 久热re国产手机在线观看| 中文字幕人成人乱码亚洲电影| 人妻精品久久无码区| 亚洲首页在线观看| 国产H片无码不卡在线视频| 免费无遮挡AV| 亚洲欧美在线综合图区| 日本在线视频免费| 国产精品美女网站| 免费国产小视频在线观看| 麻豆国产原创视频在线播放| 91成人免费观看| 中国黄色一级视频| 亚洲一区精品视频在线| 欧美日韩一区二区在线播放| AV在线天堂进入| 精品天海翼一区二区| 日韩最新中文字幕| 亚洲侵犯无码网址在线观看| 亚洲午夜福利在线| 国产成人精品一区二区| 久久久精品无码一二三区| 四虎亚洲精品| 青青草原国产| 免费视频在线2021入口| 国产97公开成人免费视频| 嫩草影院在线观看精品视频| 欧美亚洲国产精品久久蜜芽 | 国产午夜福利亚洲第一| 日本成人在线不卡视频| 欧美在线三级| 亚洲AV人人澡人人双人|