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

        ?

        改進(jìn)SO優(yōu)化KELM的液體火箭發(fā)動機(jī)故障檢測*

        2023-09-04 02:33:14王聞浩
        航天控制 2023年4期
        關(guān)鍵詞:學(xué)習(xí)機(jī)權(quán)值適應(yīng)度

        王聞浩,許 亮

        1.天津理工大學(xué)電氣工程與自動化學(xué)院,天津 300384 2.天津市復(fù)雜系統(tǒng)控制理論與應(yīng)用重點(diǎn)實(shí)驗(yàn)室,天津 300384

        0 引言

        液體火箭發(fā)動機(jī)作為液體火箭的核心,其運(yùn)行工況的正常與否直接決定了火箭是否可以成功發(fā)射。液體火箭發(fā)動機(jī)屬于多學(xué)科、多領(lǐng)域的復(fù)雜耦合系統(tǒng),當(dāng)其運(yùn)行于高溫高壓的極端環(huán)境下,極易出現(xiàn)難以預(yù)料的故障[1]。因此,對液體火箭發(fā)動機(jī)進(jìn)行故障檢測顯得尤為重要。

        分類算法是液體火箭發(fā)動機(jī)故障檢測的主要手段和方法。早期的故障檢測方法以神經(jīng)網(wǎng)絡(luò)為主。王克昌[2]將BP神經(jīng)網(wǎng)絡(luò)應(yīng)用于火箭發(fā)動機(jī)故障診斷,指出在大量訓(xùn)練樣本的支撐下,診斷結(jié)果的正確性可以得到保證。但其局限性在于,對于訓(xùn)練樣本中不存在的故障情況,診斷結(jié)果的可信度大大降低。同時(shí),當(dāng)存在參數(shù)相近,而故障模式不同的樣本時(shí),訓(xùn)練過程較慢,且難以獲得理想的結(jié)果。王平等[3]采用神經(jīng)網(wǎng)絡(luò)對滾動軸承故障包絡(luò)信號進(jìn)行自動識別。從包絡(luò)信號的時(shí)域和頻域提取故障特征信息,作為神經(jīng)網(wǎng)絡(luò)的輸入。結(jié)果表明,該神經(jīng)網(wǎng)絡(luò)可以有效識別滾動軸承的各類運(yùn)行故障。隨著分類算法的快速發(fā)展,火箭發(fā)動機(jī)故障檢測方法也更加多元。韓泉東等[4]將支持向量機(jī)應(yīng)用于火箭發(fā)動機(jī)故障檢測。采用支持向量機(jī)對所有的故障類型設(shè)計(jì)一個(gè)分類器。根據(jù)所有分類器的輸出結(jié)果,對樣本故障類型進(jìn)行判斷,仿真結(jié)果表明該方法可以有效識別故障。俞剛等[5]將支持向量機(jī)和卡爾曼濾波器結(jié)合。首先采用最小二乘支持向量機(jī),對部件故障進(jìn)行分類。最后由卡爾曼濾波器估計(jì)故障部件參數(shù),從而進(jìn)行故障檢測。

        以上算法均可實(shí)現(xiàn)火箭發(fā)動機(jī)的故障檢測,但也存在一定的問題。BP神經(jīng)網(wǎng)絡(luò)采用反向傳播學(xué)習(xí)算法,其網(wǎng)絡(luò)結(jié)構(gòu)、初始權(quán)值和閾值對網(wǎng)絡(luò)訓(xùn)練的影響很大。支持向量機(jī)分類的重要環(huán)節(jié)是懲罰因子和核函數(shù)參數(shù)的選取。對于這些問題,人們只能依靠已有的經(jīng)驗(yàn),或?qū)㈦S機(jī)參數(shù)帶入模型,通過觀察輸出結(jié)果選取最優(yōu)值。為了解決這些問題,多數(shù)學(xué)者開始使用群體算法來優(yōu)化分類模型。楊晉朝[6]采用遺傳算法優(yōu)化BP神經(jīng)網(wǎng)絡(luò)的權(quán)值和閾值,結(jié)果表明該模型可以有效避免網(wǎng)絡(luò)陷入局部最優(yōu)。Xu等[7]采用量子遺傳算法優(yōu)化BP神經(jīng)網(wǎng)絡(luò)。許亮等[8]采用粒子群算法對小波神經(jīng)網(wǎng)絡(luò)的權(quán)值和基函數(shù)參數(shù)進(jìn)行了尋優(yōu),仿真結(jié)果表明所提出的模型訓(xùn)練時(shí)間最少,且精度最高。陳占國等[9]建立了松鼠搜索算法優(yōu)化支持向量機(jī)的故障檢測模型,實(shí)驗(yàn)表明該模型具有較高的準(zhǔn)確率。馬碩[10]采用改進(jìn)后的哈里斯鷹算法優(yōu)化概率神經(jīng)網(wǎng)絡(luò)(Probabilistic Neural Networks,PNN),并與PNN和BP算法的分類結(jié)果進(jìn)行了對比。結(jié)果表明所提出的故障檢測模型要優(yōu)于另外兩種算法。

        核極限學(xué)習(xí)機(jī)(Kernel Extreme Learning Machine,KELM)自提出以來,已被應(yīng)用于多個(gè)領(lǐng)域。馬超等[11]將KELM用于液壓泵在線預(yù)測,為液壓泵的預(yù)測控制和故障預(yù)測提供了參考。陳紹煒等[12]建立了基于KELM的模擬電路故障診斷模型,并且與支持向量機(jī)和極限學(xué)習(xí)機(jī)進(jìn)行對比,結(jié)果顯示準(zhǔn)確率優(yōu)于另外兩個(gè)模型。劉鷺航等[13]采用半監(jiān)督極限學(xué)習(xí)機(jī)進(jìn)行故障檢測,解決了由于控制力矩陀螺故障數(shù)據(jù)不足導(dǎo)致的結(jié)果不準(zhǔn)確的問題。馮磊華等[14]用KELM優(yōu)化電站鍋爐,有效減少了有害氣體的排放。

        本文依據(jù)KELM在分類和擬合方面具有強(qiáng)大適應(yīng)性的優(yōu)點(diǎn),將其應(yīng)用于液體火箭發(fā)動機(jī)故障檢測。依據(jù)蛇優(yōu)化算法(Snake Optimizer,SO)較強(qiáng)的尋優(yōu)能力和存在的一些不足,對SO進(jìn)行了改進(jìn),并將改進(jìn)后的算法用于優(yōu)化核函數(shù)參數(shù)和懲罰因子。最后基于歷史樣本數(shù)據(jù),對建立的故障檢測模型進(jìn)行了實(shí)驗(yàn)仿真。

        1 核極限學(xué)習(xí)機(jī)

        單隱層前饋神經(jīng)網(wǎng)絡(luò)在進(jìn)行訓(xùn)練時(shí),多采用梯度下降法。傳統(tǒng)的學(xué)習(xí)算法存在訓(xùn)練速度慢、易陷入局部最優(yōu)、學(xué)習(xí)速率難以選擇等缺點(diǎn)。鑒于此,Huang教授提出了極限學(xué)習(xí)機(jī)(Extreme Learning Machine,ELM)算法[15]。該算法僅需設(shè)置隱含層的個(gè)數(shù)即可,輸入層與隱含層之間的權(quán)值和隱含層的閾值由隨機(jī)數(shù)生成。該方法具有訓(xùn)練速度快、泛化性能好等優(yōu)點(diǎn)。算法結(jié)構(gòu)如圖1所示。

        圖1 極限學(xué)習(xí)機(jī)結(jié)構(gòu)

        通過圖1可以看出,網(wǎng)絡(luò)的輸出受隱含層的輸出和隱含層與輸出層之間連接權(quán)值β影響。隱含層的輸出由輸入層與隱含層之間的連接權(quán)值a、隱含層的閾值b以及隱含層的激活函數(shù)G(·)決定。而輸入層與隱含層之間的權(quán)值a和隱含層的閾值b由隨機(jī)數(shù)生成。在確定激活函數(shù)后,可以求得隱含層的輸出。因此,我們的目標(biāo)是求得隱含層與輸出層之間的權(quán)值β。

        以單個(gè)訓(xùn)練樣本為例。假設(shè)網(wǎng)絡(luò)有n個(gè)輸入層,m個(gè)輸出層,l個(gè)隱含層。輸入層與隱含層之間的權(quán)值矩陣為al×n,隱含層閾值矩陣為bl×1,隱含層與輸出層之間的權(quán)值矩陣為βl×m。隱含層激活函數(shù)為G(·)。則網(wǎng)絡(luò)的輸出可以表示為:

        (1)

        式(1)可簡化為:

        Hβ=YT

        (2)

        式中:H1×l為隱含層輸出矩陣,HT為輸出矩陣Y的轉(zhuǎn)置。為使得預(yù)測輸出和實(shí)際輸出誤差盡可能小。對以下方程組進(jìn)行求解:

        (3)

        得:

        (4)

        其中:H-1為輸出矩陣H的逆矩陣。

        在2012年,Huang等將支持向量機(jī)中核函數(shù)的思想引入ELM,提出了核極限學(xué)習(xí)機(jī)[16]。在保證網(wǎng)絡(luò)結(jié)構(gòu)泛化能力的前提下,提高了系統(tǒng)的魯棒性和運(yùn)算速度,且算法的實(shí)現(xiàn)也較為簡單。

        首先引入懲罰因子C,將隱含層與輸出層之間的權(quán)值矩陣β重新定義為:

        (5)

        將核函數(shù)矩陣ΩKELM代替HHT,可得:

        (6)

        那么,KELM的輸出函數(shù)為:

        (7)

        本文核函數(shù)采用徑向基(RBF)核函數(shù),其數(shù)學(xué)表達(dá)形式為:

        (8)

        其中:σ為核函數(shù)參數(shù)。

        在實(shí)際應(yīng)用中,核極限學(xué)習(xí)機(jī)(KELM)的核函數(shù)參數(shù)σ和懲罰因子C對算法的輸出結(jié)果影響較大,但其具體數(shù)值的確定又較為困難。學(xué)者們一般通過選取不同的數(shù)值進(jìn)行多次嘗試,選擇輸出結(jié)果較好的一組數(shù)值作為參數(shù)。但該方法需要消耗大量的時(shí)間,且存在隨機(jī)性。鑒于此,本文采用一種新的群體智能算法對KELM模型參數(shù)進(jìn)行尋優(yōu),以確定模型的最優(yōu)參數(shù)。

        2 ISO-KELM故障檢測模型

        1.1 蛇優(yōu)化算法

        蛇優(yōu)化算法(SO)是由Hashim和Hussien兩位學(xué)者根據(jù)蛇的交配行為所提出的一種新的元啟發(fā)式算法[17]。該算法根據(jù)蛇在低溫和食物充足的條件下進(jìn)行交配這一特性,將算法流程分為探索和開發(fā)兩個(gè)階段。首先將所有個(gè)體劃分為雄性種群和雌性種群,在整個(gè)算法流程中,兩個(gè)種群的位置更新公式完全相同。因此下文僅對雄性種群更新公式進(jìn)行介紹。

        食物數(shù)量Q和溫度T的更新公式為:

        (9)

        (10)

        式中:N為最大迭代次數(shù),n為當(dāng)前迭代次數(shù),c1為常數(shù),設(shè)置為0.5??梢钥闯?食物數(shù)量隨著迭代次數(shù)的增加而增加,溫度隨著迭代次數(shù)的增加而降低。

        當(dāng)Q<0.25時(shí),算法處于探索階段,種群個(gè)體通過移動自身位置來尋找食物,更新公式為:

        Xi(t+1)=Xrand(t)±c2×A×

        ((Xmax-Xmin)×r+Xmin)

        (11)

        (12)

        式中:Xi為更新后的雄性個(gè)體位置,fi為該個(gè)體的適應(yīng)度。Xrand為隨機(jī)雄性個(gè)體,frand為隨機(jī)雄性個(gè)體的適應(yīng)度。Xmax和Xmin分別為求解問題的上邊界和下邊界。r為0~1之間的隨機(jī)數(shù)。c2為常數(shù),設(shè)置為0.05。

        當(dāng)Q>0.25時(shí),算法進(jìn)入開發(fā)階段。在該階段,當(dāng)T>0.6時(shí),雄性個(gè)體在食物周圍進(jìn)行開發(fā),其更新公式為:

        Xi(t+1)=Xfood±c3×T×r×(Xfood-Xi(t))

        (13)

        式中:Xi為雄性個(gè)體位置,Xfood為所有種群中最優(yōu)個(gè)體的位置。c3為常數(shù),設(shè)置為2。

        當(dāng)T≤0.6時(shí),雄性個(gè)體將進(jìn)入戰(zhàn)斗模式或交配模式,由隨機(jī)數(shù)決定。戰(zhàn)斗模式更新公式為:

        Xi(t+1)=Xi(t)+c3×Fm×r×

        (Q×Xbest,f-Xi(t))

        (14)

        (15)

        式中:Xi為雄性個(gè)體位置,fi為其適應(yīng)度。Xbest,f為雌性群體中的最優(yōu)個(gè)體位置,fbest,f為其適應(yīng)度。

        交配模式更新公式為:

        Xi,m(t+1)=Xi,m(t)+c3×Mm×r×

        (Q×Xi,f(t)-Xi,m(t))

        (16)

        (17)

        式中:Xi,m為雄性個(gè)體位置,fi,m為其適應(yīng)度。Xi,f為雌性個(gè)體位置,fi,f為其適應(yīng)度。

        交配之后有50%的概率會產(chǎn)生后代,用以代替雄性種群中的最差個(gè)體。

        2.2 蛇優(yōu)化算法改進(jìn)

        蛇優(yōu)化算法(SO)的初始種群由軟件的隨機(jī)數(shù)函數(shù)生成,可能會使得初始個(gè)體離最優(yōu)點(diǎn)較遠(yuǎn),導(dǎo)致算法尋優(yōu)能力不足;SO探索和開發(fā)階段的轉(zhuǎn)換由食物數(shù)量(Q)來決定。由于Q的大小僅受迭代次數(shù)影響,使得SO進(jìn)入開發(fā)階段較晚,導(dǎo)致算法后期尋優(yōu)能力不足;同時(shí),SO后期種群多樣性減少,難以跳出局部最優(yōu)點(diǎn)。

        針對以上問題,分別提出以下3點(diǎn)算法改進(jìn)策略。

        2.2.1 Tent映射

        針對初始化種群離散度不高,導(dǎo)致算法尋優(yōu)能力下降的問題。采用Tent混沌映射對種群進(jìn)行初始化操作。Tent映射表達(dá)式為:

        (18)

        本文將α的數(shù)值設(shè)置為0.499。因?yàn)楫?dāng)α等于0.5時(shí),在經(jīng)過50次左右的迭代后,種群會收斂至0,失去遍歷性。

        2.2.2 動態(tài)策略

        針對算法較晚進(jìn)入開發(fā)階段,導(dǎo)致尋優(yōu)能力較弱的特點(diǎn),采用動態(tài)策略轉(zhuǎn)換探索和開發(fā)階段。為了使SO提前結(jié)束探索階段,增加開發(fā)階段的迭代次數(shù),采用以下動態(tài)策略控制階段轉(zhuǎn)換:在探索階段,當(dāng)種群最優(yōu)適應(yīng)度連續(xù)6代不發(fā)生變化的時(shí)候,SO進(jìn)入開發(fā)階段。動態(tài)策略公式為:

        D(t+1)=c4×D(t)+|fbest(t)-fbest(t-1)|

        (19)

        式中:D為動態(tài)參數(shù)。fbest為種群最優(yōu)適應(yīng)度。c4為常數(shù),設(shè)置為0.1。當(dāng)D的值小于1×10-6時(shí),SO進(jìn)入開發(fā)階段。

        此外,將常數(shù)c4引入探索階段雌性種群的位置更新公式,加快種群的收斂。由于SO的特殊性,在雌性種群快速收斂的同時(shí),雄性種群依舊可以在全局進(jìn)行遍歷。

        2.2.3 柯西變異

        針對SO后期易陷入局部最優(yōu)的特點(diǎn),采取柯西變異策略對陷入局部最優(yōu)的種群進(jìn)行處理。

        柯西變異策略來源于柯西分布??挛鞣植嫉母怕拭芏群瘮?shù)為:

        (20)

        當(dāng)a=1時(shí),稱為標(biāo)準(zhǔn)柯西分布。在SO中,柯西分布隨機(jī)變量生成函數(shù)為:

        μ=tan(π×(r-0.5))

        (21)

        r為0~1之間的隨機(jī)數(shù)。個(gè)體位置擾動公式為:

        Xworst(t+1)=Xworst(t)×(1+μ)

        (22)

        式中:Xworst為最優(yōu)適應(yīng)度較差的種群。

        進(jìn)入開發(fā)階段后,當(dāng)全局最優(yōu)值連續(xù)6代沒有發(fā)生變化的時(shí)候,對兩個(gè)種群中最好適應(yīng)度較差的群體進(jìn)行柯西變異。

        改進(jìn)后的算法流程如圖2所示。

        圖2 改進(jìn)蛇優(yōu)化算法流程圖

        2.3 算法優(yōu)化流程

        針對KELM參數(shù)難以確定的問題,采用上文改進(jìn)后的蛇優(yōu)化算法對參數(shù)進(jìn)行尋優(yōu)。

        整個(gè)算法流程如圖3所示。

        圖3 ISO-KELM模型流程圖

        具體尋優(yōu)過程為:

        1)導(dǎo)入數(shù)據(jù),包括數(shù)據(jù)樣本和標(biāo)簽;

        2)將樣本劃分為訓(xùn)練集和測試集;

        3)對訓(xùn)練集數(shù)據(jù)進(jìn)行歸一化處理,再使用訓(xùn)練集歸一化時(shí)的映射對測試集進(jìn)行歸一化;

        4)初始化ISO算法參數(shù),包括種群數(shù)量,種群維度,種群上下限和測試函數(shù);

        5)建立ISO-KELM模型,對KELM參數(shù)進(jìn)行尋優(yōu)。具體過程如圖2所示;

        6)將所求得的最優(yōu)懲罰因子C和核函數(shù)參數(shù)σ賦予KELM;

        7)使用優(yōu)化后的KELM進(jìn)行故障檢測;

        8)對尋優(yōu)結(jié)果進(jìn)行評價(jià)對比。

        3 實(shí)驗(yàn)仿真

        3.1 故障模式

        本次實(shí)驗(yàn)數(shù)據(jù)選取某型號液體火箭發(fā)動機(jī)歷史試車數(shù)據(jù),共420組。正常數(shù)據(jù)和故障數(shù)據(jù)分別為210組。故障數(shù)據(jù)包括渦輪葉片損壞、燃燒室喉部燒蝕、氧化劑泵汽蝕、氧化劑泵管路堵塞和氧增壓系統(tǒng)減壓閥故障5種典型故障模式。同時(shí),為了使所提出算法具有較高的準(zhǔn)確性和適用性,選取14個(gè)代表性強(qiáng)且相關(guān)性弱的發(fā)動機(jī)參數(shù)作為輸入。

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

        本次實(shí)驗(yàn)采用仿真軟件進(jìn)行仿真。選取前42組數(shù)據(jù)作為測試集,剩余數(shù)據(jù)作為訓(xùn)練集。訓(xùn)練集中的前42組數(shù)據(jù)作為驗(yàn)證集,并將驗(yàn)證集的準(zhǔn)確率作為ISO的適應(yīng)度。

        ISO參數(shù)設(shè)置如下:種群數(shù)量設(shè)置為30,種群維度為2,第一維代表懲罰因子C,第二維代表核函數(shù)參數(shù)σ。種群上限和下限分別設(shè)置為100和0.001。

        3.3 結(jié)果對比

        我們將核極限學(xué)習(xí)機(jī)(KELM)、蛇優(yōu)化算法優(yōu)化核極限學(xué)習(xí)機(jī)(SO-KELM)、改進(jìn)后的蛇優(yōu)化算法優(yōu)化核極限學(xué)習(xí)機(jī)(ISO-KELM)3種模型的故障檢測結(jié)果和傳統(tǒng)BP神經(jīng)網(wǎng)絡(luò)模型故障檢測結(jié)果進(jìn)行對比。部分模型輸出結(jié)果如圖4所示。

        圖4 部分模型輸出結(jié)果

        為便于直觀了解各個(gè)故障檢測模型的性能,采用誤報(bào)率、漏報(bào)率和準(zhǔn)確率對4種模型進(jìn)行評估。具體結(jié)果如表1所示。

        表1 模型檢測結(jié)果對比

        通過表1可以看出,4種模型的誤報(bào)率均為0。漏報(bào)率的高低直接影響了模型的準(zhǔn)確率。其中,KELM模型的漏報(bào)率相比于傳統(tǒng)BP神經(jīng)網(wǎng)絡(luò)較低。同時(shí),分別采用SO和ISO優(yōu)化KELM后的模型均能提高檢測結(jié)果的準(zhǔn)確率,且ISO-KELM模型的漏報(bào)率更低,檢測結(jié)果也更加準(zhǔn)確。

        為了驗(yàn)證改進(jìn)蛇優(yōu)化算法具有更好的性能,采用適應(yīng)度曲線作為算法的評估條件。

        圖5分別為ISO和SO的適應(yīng)度變化曲線??梢钥闯?在算法開始運(yùn)行時(shí),ISO的適應(yīng)度低于SO的適應(yīng)度,原因是ISO的初始種群遍歷性較強(qiáng),在前期有更強(qiáng)的尋優(yōu)能力。由于ISO在探索階段尋優(yōu)能力有限,因此在前幾次迭代過程中適應(yīng)度沒有發(fā)生變化。此時(shí)算法提前結(jié)束探索階段,進(jìn)入開發(fā)階段。且在第18代尋得最優(yōu)值。而SO整體收斂速度較慢,在第25代完成收斂,且尋得的最優(yōu)值大于ISO。

        圖5 算法適應(yīng)度曲線

        4 結(jié)論

        使用改進(jìn)蛇優(yōu)化算法優(yōu)化核極限學(xué)習(xí)機(jī)(ISO-KELM),并將其應(yīng)用于液體火箭發(fā)動機(jī)故障檢測??梢缘贸鲆韵陆Y(jié)論:

        (1)改進(jìn)后的蛇優(yōu)化算法可以提高種群前期的尋優(yōu)能力,使種群提前進(jìn)入開發(fā)階段,加快后期收斂,且收斂能力更強(qiáng);

        (2)ISO可以和KELM模型相結(jié)合,通過ISO獲得KELM模型的最優(yōu)參數(shù),從而實(shí)現(xiàn)對液體火箭發(fā)動機(jī)故障的有效檢測;

        (3)仿真實(shí)驗(yàn)表明,所提出ISO-KELM故障檢測模型的準(zhǔn)確率優(yōu)于傳統(tǒng)BP神經(jīng)網(wǎng)絡(luò)故障檢測算法,實(shí)驗(yàn)結(jié)果更加理想。

        猜你喜歡
        學(xué)習(xí)機(jī)權(quán)值適應(yīng)度
        改進(jìn)的自適應(yīng)復(fù)制、交叉和突變遺傳算法
        一種融合時(shí)間權(quán)值和用戶行為序列的電影推薦模型
        CONTENTS
        極限學(xué)習(xí)機(jī)綜述
        基于極限學(xué)習(xí)機(jī)參數(shù)遷移的域適應(yīng)算法
        分層極限學(xué)習(xí)機(jī)在滾動軸承故障診斷中的應(yīng)用
        基于權(quán)值動量的RBM加速學(xué)習(xí)算法研究
        基于空調(diào)導(dǎo)風(fēng)板成型工藝的Kriging模型適應(yīng)度研究
        中國塑料(2016年11期)2016-04-16 05:26:02
        一種基于AdaBoost的極限學(xué)習(xí)機(jī)分類方法
        少數(shù)民族大學(xué)生文化適應(yīng)度調(diào)查
        久久午夜av一区二区三区| 亚洲国产字幕| 亚洲国产人成自精在线尤物| 自拍偷自拍亚洲精品第按摩| 色妞色视频一区二区三区四区| 日韩欧美亚洲综合久久影院d3| 日韩偷拍视频一区二区三区| 97色综合| 亚洲综合精品一区二区| 亚洲av无码偷拍在线观看| 四虎国产精品免费久久| 国产精品乱一区二区三区| 久久国产精品色av免费看| 日本丰满少妇裸体自慰| 亚洲av无码不卡| 久久夜色精品国产亚洲噜噜| 久亚洲一线产区二线产区三线麻豆 | 国产精品自拍网站在线| 亚洲av无码成h在线观看| 久久老子午夜精品无码怎么打| 老熟女熟妇嗷嗷叫91| 日本一区二区三区综合视频| 欧美综合天天夜夜久久| 国产精品-区区久久久狼| 亚洲一区二区国产精品视频| 丰满少妇人妻久久精品| 人妻夜夜爽天天爽一区| 中文字幕亚洲好看有码| 亚洲国产最新免费av| 亚洲国产精品久久电影欧美| 尤物yw无码网站进入| 久久色悠悠亚洲综合网| 亚洲精品宾馆在线精品酒店| 免费无码毛片一区二区三区a片| 久久精品国产只有精品96| 国产一区二区三区成人| 久久99精品久久水蜜桃| 性导航app精品视频| 高清成人在线视频播放| 色哟哟最新在线观看入口| 国农村精品国产自线拍|