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

        ?

        滿秩分解最小二乘法船舶航向模型辨識(shí)

        2022-02-18 08:12:30包政凱朱齊丹劉永超
        智能系統(tǒng)學(xué)報(bào) 2022年1期
        關(guān)鍵詞:船舶模型

        包政凱,朱齊丹,劉永超

        (哈爾濱工程大學(xué) 智能科學(xué)與工程學(xué)院, 黑龍江 哈爾濱 150001)

        隨著數(shù)字化和智能化技術(shù)的發(fā)展,智能船舶發(fā)展呈現(xiàn)出日益加快的趨勢(shì)[1-3]。我國(guó)目前處于《智能船舶發(fā)展行動(dòng)計(jì)劃2019-2021年》提到的三年目標(biāo)階段,正在實(shí)現(xiàn)智能船1.0向智能船2.0的過渡。自動(dòng)靠離泊、遠(yuǎn)程遙控、自主航行等智能船典型應(yīng)用場(chǎng)景需要精確的船舶航向自適應(yīng)控制算法[4-7]。國(guó)內(nèi)已產(chǎn)品化的航向控制器與國(guó)外先進(jìn)產(chǎn)品相比,航向精度、能耗、轉(zhuǎn)向快速性和轉(zhuǎn)向平穩(wěn)性等性能指標(biāo)仍有較大差距,主要原因有:1) 航向模型參數(shù)在線辨識(shí)的實(shí)時(shí)性、收斂性和辨識(shí)精度難以同時(shí)保證,在線估計(jì)時(shí)變參數(shù)和擾動(dòng)的能力不夠;2) 海洋高頻擾動(dòng)、測(cè)量噪聲和測(cè)量病態(tài)引起的航向角速度信號(hào)波動(dòng)處理困難;3) 舵角控制環(huán)應(yīng)對(duì)航速和海況變化的自適應(yīng)能力不足;4) 操舵執(zhí)行機(jī)構(gòu)性能受限,與國(guó)外接近實(shí)現(xiàn)連續(xù)控制效果有較大差距。

        針對(duì)上述第一個(gè)問題,目前常用的船舶航向模型參數(shù)辨識(shí)算法主要分為2類:自校正自適應(yīng)和模型參考自適應(yīng)[8]。自校正自適應(yīng)常用的辨識(shí)算法包括最小二乘法[9-14]、支持向量機(jī)[15-21]、神經(jīng)網(wǎng)絡(luò)[22-24]、極大似然估計(jì)[25]、遺傳算法[26]等。文獻(xiàn)[11]采用最小二乘算法辨識(shí)船舶操縱模型的參數(shù),并引入新型的學(xué)習(xí)率,提升了辨識(shí)精度。文獻(xiàn)[12]針對(duì)船舶響應(yīng)模型參數(shù)辨識(shí)問題,在最小二乘算法遞推迭代過程中將單新息拓展為多新息,大大改進(jìn)了辨識(shí)算法的準(zhǔn)確性。文獻(xiàn)[13]為了提升辨識(shí)過程的快速性,通過構(gòu)造辨識(shí)誤差的評(píng)價(jià)函數(shù),提出了一種動(dòng)態(tài)遺忘因子最小二乘辨識(shí)算法。文獻(xiàn)[14]在船舶模型辨識(shí)最小二乘算法中引入一個(gè)非線性的反正切函數(shù)來(lái)處理新息,有效地提升了算法的辨識(shí)精度。文獻(xiàn)[18-20]采用支持向量機(jī)對(duì)船舶航向模型參數(shù)進(jìn)行辨識(shí)研究。支持向量機(jī)具有全局最優(yōu)搜索能力,得到的模型泛化能力強(qiáng),結(jié)合優(yōu)化算法可以提升辨識(shí)精度并簡(jiǎn)化調(diào)參。文獻(xiàn)[21]調(diào)查了船舶操縱運(yùn)動(dòng)的水動(dòng)力模型,將神經(jīng)網(wǎng)絡(luò)應(yīng)用到船舶操縱模型和三自由度模型辨識(shí)中。文獻(xiàn)[22]針對(duì)大型油船的參數(shù)辨識(shí)問題,提出了基于人工神經(jīng)網(wǎng)絡(luò)的參數(shù)辨識(shí)方法。文獻(xiàn)[23]提出了油輪轉(zhuǎn)向數(shù)學(xué)模型,通過廣義橢球基函數(shù)神經(jīng)網(wǎng)絡(luò)算法辨識(shí)轉(zhuǎn)向模型參數(shù)。神經(jīng)網(wǎng)絡(luò)辨識(shí)模型參數(shù)對(duì)初值沒有依賴性,并且可以保證辨識(shí)參數(shù)的收斂性。文獻(xiàn)[25]針對(duì)水面無(wú)人艇操縱性預(yù)報(bào)問題,利用極大似然法辨識(shí)無(wú)人艇參數(shù),采用辨識(shí)的結(jié)果預(yù)報(bào)了無(wú)人艇的操縱運(yùn)動(dòng)。文獻(xiàn)[26]利用船舶操縱模型的航向角數(shù)據(jù)構(gòu)造目標(biāo)函數(shù),然后選擇遺傳算法進(jìn)行尋優(yōu),提升了辨識(shí)效率。模型參考自適應(yīng)[27]設(shè)計(jì)方法中,采用構(gòu)造李雅普諾夫函數(shù),并與航向控制器聯(lián)合設(shè)計(jì),可以保證模型參數(shù)辨識(shí)李雅普諾夫意義下穩(wěn)定,同時(shí)實(shí)現(xiàn)航向漸進(jìn)穩(wěn)定。

        上述方法對(duì)于提高模型參數(shù)辨識(shí)精度、收斂速度和實(shí)時(shí)性有很多貢獻(xiàn),但缺乏對(duì)辨識(shí)算法收斂性分析。工程中通常使用滿足數(shù)據(jù)激勵(lì)條件的回轉(zhuǎn)、Z型、自由航行實(shí)驗(yàn)數(shù)據(jù)來(lái)保證模型辨識(shí)算法的收斂性,但在自適應(yīng)控制過程中,持續(xù)激勵(lì)條件往往無(wú)法滿足,使得參數(shù)辨識(shí)瞬時(shí)收斂性問題成為自校正自適應(yīng)控制的難點(diǎn)。對(duì)于數(shù)據(jù)長(zhǎng)時(shí)間欠激勵(lì)、受擾條件下如何保證或提升在線辨識(shí)算法的收斂性問題還沒有被解決。

        本文針對(duì)數(shù)據(jù)欠激勵(lì)受擾系統(tǒng),設(shè)計(jì)了一種快速、準(zhǔn)確、穩(wěn)定的線性模型在線辨識(shí)通用方法,并將該算法應(yīng)用于船舶航向模型參數(shù)在線辨識(shí)過程中。將動(dòng)態(tài)降維的思想引入遺忘因子最小二乘法 (forgetting factor least square algorithm, FFLS),提出了一種基于滿秩分解的最小二乘法(full rank decomposition least square algorithm, FRDLS),使得欠激勵(lì)數(shù)據(jù)對(duì)降維后的辨識(shí)模型重新滿足激勵(lì)條件。所提供的算法能夠根據(jù)死區(qū)閾值動(dòng)態(tài)調(diào)整滿秩分解矩陣,使最小二乘法選擇性更新模型參數(shù),限制了擾動(dòng)導(dǎo)致的參數(shù)辨識(shí)誤差累計(jì)。在樣本數(shù)據(jù)長(zhǎng)時(shí)間欠激勵(lì), 海洋環(huán)境干擾工況下,將實(shí)船數(shù)據(jù)用于航向模型參數(shù)辨識(shí),對(duì)比FFLS和FRDLS算法、多新息最小二乘法[12](multi-innovation least square algorithm, MILS)、最小二乘支持向量機(jī)[19](least square support vector machines algorithm, LSSVM)辨識(shí)船舶操縱模型參數(shù)的辨識(shí)結(jié)果曲線和航向角速度預(yù)報(bào)誤差曲線,可以發(fā)現(xiàn)滿秩分解最小二乘法能有效抑制數(shù)據(jù)欠激勵(lì)和存在海洋環(huán)境擾動(dòng)時(shí)參數(shù)辨識(shí)過程中的漂移和發(fā)散問題,并提高參數(shù)辨識(shí)的精度,同時(shí)具有計(jì)算量小、實(shí)時(shí)性強(qiáng)、簡(jiǎn)單易實(shí)現(xiàn)的優(yōu)點(diǎn)。

        1 船舶航向模型和基本辨識(shí)算法

        遠(yuǎn)洋航行的船舶多工作在航向保持模式下,該工況下,船舶自動(dòng)舵通過給定小舵角抵抗海洋環(huán)境干擾,使實(shí)際航向保持在期望航向誤差帶內(nèi),辨識(shí)準(zhǔn)確的航向模型參數(shù)是實(shí)現(xiàn)高精度航向控制的前提。

        1.1 船舶航向模型

        船舶操縱運(yùn)動(dòng)模型可以采用一階線性Nomoto 模型[13]:

        式中:r為船體艏搖角速度;r˙為艏搖角加速度;φ為船舶艏向角;δ為舵角;K、T為船舶回轉(zhuǎn)能力和回轉(zhuǎn)慣性參數(shù),受到船舶載荷條件和航行速度的影響。最小二乘改進(jìn)模型辨識(shí)算法非常依賴于經(jīng)典狀態(tài)空間模型參數(shù),實(shí)際工程中考慮到外界擾動(dòng)和模型攝動(dòng)對(duì)參數(shù)辨識(shí)的不利影響,此處引入擾動(dòng)項(xiàng)對(duì)Nomoto模型進(jìn)行改進(jìn):

        式(2)中 δd為干擾力矩等效舵角,包含了舵的零位誤差、船舶航向模型非線性、測(cè)量噪聲中的非零均值部分、海洋環(huán)境擾動(dòng)力矩等。為辨識(shí)做準(zhǔn)備,將改進(jìn)Nomoto模型整理為標(biāo)準(zhǔn)最小二乘形式。

        為了加快參數(shù)辨識(shí)收斂速度,引入歸一化因子 αr˙和 αr,辨識(shí)前通常使 用歸一化因子使數(shù)據(jù) 向量各維度具有相同幅值范圍,令:

        選取輸出量zk、參數(shù)向量 θk、數(shù)據(jù)向量hk,下標(biāo)k代表第k時(shí)刻:

        整理得最小二乘測(cè)量方程為

        式中wk為k時(shí)刻測(cè)量噪聲。

        1.2 遺忘因子最小二乘法

        遞推最小二乘算法因其簡(jiǎn)單、實(shí)時(shí)性強(qiáng)的優(yōu)點(diǎn),廣泛應(yīng)用于模型參數(shù)在線辨識(shí)。設(shè)給定數(shù)據(jù)集,其中,hk=[hk1hk2···hkn]T為一次測(cè)量數(shù)據(jù) 向 量 , 令 數(shù) 據(jù) 矩 陣 為H=[h1h2···hN]T∈ RN×n,輸出矩陣為Z=[z1z2···zN]T,則線性回歸方程可以表示為

        式中: θ =[θ1θ2···θn]T為待估計(jì)系數(shù)向量;w=[w1w2···wN]T為N次測(cè)量的模型擬合誤差向量;θ=[θ1θ2···θn]T為待辨識(shí)參數(shù)數(shù)量。最小二乘法(least square algorithm, LS)采用誤差向量二范數(shù)為損失函數(shù):

        矩陣H的秩 R ank(H)=n時(shí),HTH可逆,由式(9)可得 θ的最小二乘估計(jì):

        為了將最小二乘轉(zhuǎn)化為遞推形式,可令k時(shí)刻數(shù)據(jù)矩陣、測(cè)量矩陣、方差矩陣和辨識(shí)結(jié)果矩陣分別為

        為避免隨著數(shù)據(jù)增多,新數(shù)據(jù)作用逐漸消失,引入遺忘因子 β1(0.9< β1<1) 得:

        2 基于滿秩分解的最小二乘算法設(shè)計(jì)

        當(dāng)式(8)中矩陣H的秩r=Rank(H)

        式中:

        正交投影變換矩陣V的選取不唯一,本文結(jié)合遞推最小二乘法給出V的一種求解方法。設(shè)第k次測(cè)量數(shù)據(jù) {hk,zk}對(duì)應(yīng)的變換陣為Vk=[Vk1Vk2]。

        式中:ei為第i行為1,其余行為0的n維列向量;hs為死區(qū)閾值,用于削弱干擾的影響;hki為第k次測(cè)量數(shù)據(jù)向量的第i行數(shù)據(jù),當(dāng)i=1,2,···,n時(shí),有p=1,2,···,r或q=1,2,···,n?r,最終可得:

        滿秩分解最小二乘法算法實(shí)現(xiàn)步驟為

        1) 初始化式(13)中方差陣Pk、參數(shù)向量 θk;初始化上一時(shí)刻正交變換陣Vk?1為n階單位陣;初始化滿秩分解后方差陣等于Pk,滿秩分解后參數(shù)向量等于 θk;

        2) 通過測(cè)量值向量hk根據(jù)式(19)、(20)求取正交變換矩陣Vk;

        3) 若Vk不等于Vk?1,根據(jù)式(18)求和,否則跳到4);

        4) 根據(jù)式(13)求Pk+1;根據(jù)式(17)求和;根據(jù)式(16)由求,令等于

        6) 循環(huán)2)~5)可實(shí)現(xiàn)參數(shù)在線辨識(shí)。

        3 實(shí)驗(yàn)驗(yàn)證

        用于驗(yàn)證辨識(shí)算法的實(shí)驗(yàn)數(shù)據(jù)為三級(jí)海況下,一艘貨船的自動(dòng)舵數(shù)據(jù)記錄儀從大連到新加坡航行過程中記錄獲得,系統(tǒng)采樣頻率為50 Hz。艏搖角速度r、舵角 δ和航向 ψ 數(shù)據(jù)如圖1所示,角速度信號(hào)測(cè)量獲得后,經(jīng)過IIR低通濾波器濾除高頻測(cè)量噪聲后用于模型辨識(shí)。

        圖1 實(shí)驗(yàn)數(shù)據(jù)曲線Fig.1 Experiment data curves

        3.1 船舶航向模型參數(shù)辨識(shí)

        分別采用LS、FFLS、MILS、LSSVM和FRDLS算法,基于圖1中的實(shí)驗(yàn)數(shù)據(jù)對(duì)式(5)中的模型參數(shù)進(jìn)行辨識(shí),對(duì)比辨識(shí)結(jié)果,分析在辨識(shí)數(shù)據(jù)充分激勵(lì)和欠激勵(lì)工況下,辨識(shí)算法的辨識(shí)精度、收斂性和實(shí)時(shí)性。

        MILS算法新息數(shù)選為10,F(xiàn)RDLS算法與FFLS辨識(shí)算法遺忘因子均為 β =0.99994,F(xiàn)RDLS算法閾值參數(shù)hs=1,LSSVM算法選擇懲罰因子C為20,數(shù)據(jù)滑動(dòng)窗口長(zhǎng)度為500 s,K、T的初值選為0.2和250,擾動(dòng) δd初值為0,式(4)中歸一化系數(shù)分別為 αr˙=25和 αr=50。對(duì)比5種辨識(shí)算法,參數(shù)辨識(shí)結(jié)果曲線如圖2、3所示。

        圖2 K辨識(shí)結(jié)果曲線Fig.2 K-coefficient curves

        圖3 T辨識(shí)結(jié)果曲線Fig.3 T-coefficient curves

        從圖2、3可以看出,F(xiàn)FLS、MILS算法辨識(shí)K、T參數(shù)過程中,在200 s和3 000 s附近均出現(xiàn)參數(shù)發(fā)散,LSSVM算法數(shù)據(jù)滑動(dòng)窗口進(jìn)入1 200秒前后的船舶兩次轉(zhuǎn)向數(shù)據(jù)段時(shí),K、T參數(shù)辨識(shí)收斂,滑動(dòng)窗口進(jìn)入1 800s之后航向保持?jǐn)?shù)據(jù)段,辨識(shí)參數(shù)開始發(fā)散;FRDLS算法辨識(shí)K、T參數(shù)和擾動(dòng)過程相對(duì)平穩(wěn),在船舶轉(zhuǎn)向數(shù)據(jù)段和長(zhǎng)時(shí)間航向保持?jǐn)?shù)據(jù)段均未出現(xiàn)參數(shù)發(fā)散。辨識(shí)算法引入滿秩分解后提高了遺忘因子最小二乘法參數(shù)辨識(shí)過程的魯棒性。參數(shù)辨識(shí)結(jié)果見表1。

        表1 參數(shù)辨識(shí)結(jié)果Table 1 Parameters identification results

        3.2 船舶運(yùn)動(dòng)預(yù)報(bào)

        為了對(duì)比上述算法辨識(shí)船舶航向模型參數(shù)的精度,可以將圖2~圖4中參數(shù)辨識(shí)結(jié)果代入式(3)中,以舵角數(shù)據(jù)為輸入,對(duì)式(3)微分方程進(jìn)行求解,從而預(yù)報(bào)船舶艏搖角速度。將預(yù)報(bào)結(jié)果與實(shí)驗(yàn)數(shù)據(jù)作差,可得到預(yù)報(bào)誤差:

        圖4 δd辨識(shí)結(jié)果曲線Fig.4 δd-coefficient curves

        式中ri和rip分別為第i時(shí)刻航向角速度測(cè)量值和預(yù)報(bào)值。預(yù)報(bào)誤差越小,則辨識(shí)精度越高。為了保證航向角速度預(yù)報(bào)精度,此處采用四階龍格庫(kù)塔法求解航向模型微分方程。

        式中:x=[rψ]T;f(x,t)=[(?r+Kδ )/Tr]T;Ts為系統(tǒng)采樣周期,此處取Ts=0.02s。角速度預(yù)報(bào)曲線和預(yù)報(bào)誤差曲線如圖5、6。

        圖5 角速度預(yù)報(bào)曲線Fig.5 Yaw prediction curve

        圖5為航向模型參數(shù)辨識(shí)結(jié)果預(yù)報(bào)航向角速度曲線和實(shí)驗(yàn)測(cè)量角速度曲線。圖6中給出了參數(shù)辨識(shí)結(jié)果角速度預(yù)報(bào)誤差。從0~600 s預(yù)報(bào)誤差曲線對(duì)比可知,F(xiàn)RDLS、MILS、FFLS相比于LS、LSSVM算法預(yù)報(bào)誤差下降更快,有更快的辨識(shí)收斂速度。為了定量描述參數(shù)辨識(shí)精度,采用狀態(tài)預(yù)報(bào)均方根誤差(root mean square error,RMSE)作為航向模型辨識(shí)結(jié)果的精度指標(biāo),可根據(jù)式(23)求航向角速度預(yù)報(bào)均方根誤差。

        圖6 角速度預(yù)報(bào)誤差曲線Fig.6 Yaw prediction error curve

        式中ei為式(21)中第i組數(shù)據(jù)航向角速度預(yù)報(bào)誤差,角速度預(yù)報(bào)誤差結(jié)果和辨識(shí)過程收斂性見表2。

        表2 角速度預(yù)報(bào)誤差Table 2 Yaw prediction error

        根據(jù)表2,F(xiàn)FLS相比于LS算法,引入遺忘因子后參數(shù)辨識(shí)精度有所提高,同時(shí)導(dǎo)致辨識(shí)過程產(chǎn)生了發(fā)散現(xiàn)象;FRDLS算法引入滿秩分解后,提升了FFLS算法辨識(shí)精度的同時(shí),消除了辨識(shí)過程中的發(fā)散,說(shuō)明遺忘因子和滿秩分解兩種新息處理方式相結(jié)合具有合理性。LSSVM算法在1 800 s后發(fā)散部分較多,對(duì)應(yīng)的預(yù)報(bào)均方根誤差結(jié)果為 94 943°·s?1,F(xiàn)RDLS 算法在 1 800 s后未出現(xiàn)參數(shù)發(fā)散現(xiàn)象,說(shuō)明了FRDLS算法在數(shù)據(jù)欠激勵(lì)、環(huán)境擾動(dòng)工況下具有更好的收斂性。

        3.3 算法實(shí)時(shí)性分析

        辨識(shí)船舶航向模型參數(shù)的過程中,使用了179 996組測(cè)量數(shù)據(jù)對(duì)參數(shù)估計(jì)結(jié)果進(jìn)行了更新,采用Qt軟件C++編程語(yǔ)言,在window10系統(tǒng)下記錄了FFLS、MILS、LSSVM和 FRDLS算法每次使用測(cè)量數(shù)據(jù)從計(jì)算數(shù)據(jù)向量到參數(shù)辨識(shí)值更新的時(shí)間,得到算法單次循環(huán)耗時(shí)曲線如圖7所示。

        圖7 算法單次循環(huán)耗時(shí)曲線Fig.7 Algorithm time consumption for a single loop

        算法單次循環(huán)最大耗時(shí)和平均耗時(shí)見表3。從算法耗時(shí)統(tǒng)計(jì)結(jié)果可知,F(xiàn)RDLS算法單次循環(huán)平均耗時(shí)為27.39 μ s,遠(yuǎn)小于LSSVM算法,耗時(shí)介于FFLS算法和MILS算法之間。由于滿秩分解過程會(huì)額外消耗算法單次循環(huán)時(shí)間,犧牲了FRDLS算法部分實(shí)時(shí)性,以提升收斂性,耗時(shí)為FFLS算法的3倍左右,實(shí)時(shí)性在可接受范圍以內(nèi)。

        表3 算法單次循環(huán)耗時(shí)Table 3 Algorithm time consumption for a single loop

        4 結(jié)束語(yǔ)

        本文針對(duì)數(shù)據(jù)欠激勵(lì)受擾動(dòng)系統(tǒng)參數(shù)在線辨識(shí)算法實(shí)時(shí)性、參數(shù)收斂性和辨識(shí)精度難以同時(shí)保證的問題,提出了一種滿秩分解最小二乘算法,并使用該算法對(duì)長(zhǎng)時(shí)間數(shù)據(jù)欠激勵(lì)和海洋環(huán)境擾動(dòng)工況下的船舶航向運(yùn)動(dòng)模型參數(shù)進(jìn)行辨識(shí)。辨識(shí)結(jié)果表明,存在數(shù)據(jù)欠激勵(lì)時(shí),該算法能抑參數(shù)發(fā)散,并獲得更高的辨識(shí)精度。且改進(jìn)后的算法仍然計(jì)算量小、實(shí)現(xiàn)簡(jiǎn)單、實(shí)時(shí)性強(qiáng)、通用性好,具有較高實(shí)用價(jià)值,為實(shí)現(xiàn)準(zhǔn)確穩(wěn)定的航向自適應(yīng)控制和欠驅(qū)動(dòng)船自動(dòng)靠泊奠定了基礎(chǔ)。

        猜你喜歡
        船舶模型
        一半模型
        計(jì)算流體力學(xué)在船舶操縱運(yùn)動(dòng)仿真中的應(yīng)用
        基于改進(jìn)譜分析法的船舶疲勞強(qiáng)度直接計(jì)算
        重要模型『一線三等角』
        重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
        船舶!請(qǐng)加速
        BOG壓縮機(jī)在小型LNG船舶上的應(yīng)用
        船舶壓載水管理系統(tǒng)
        3D打印中的模型分割與打包
        FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
        亚洲天堂在线视频播放| 国产一级一厂片内射视频播放| 日本国产一区二区在线| 中文字幕有码在线亚洲| 亚洲97成人在线视频| 天天做天天爱夜夜夜爽毛片 | 精品久久中文字幕系列| 国产亚洲成人精品久久| 天天做天天爱夜夜爽女人爽| 国产精选污视频在线观看| 久久久久亚洲av无码观看| 欧美日韩亚洲成人| 丰满熟妇人妻av无码区| 国产女人av一级一区二区三区| 麻豆最新国产av原创| 国产精品国产三级国产av品爱网 | 91短视频在线观看免费| 亚洲国产高清在线视频| 中文字幕中文字幕人妻黑丝| 一本色道久久88加勒比综合| 西川结衣中文字幕在线| 天堂а在线中文在线新版| 毛片内射久久久一区| 亚洲性啪啪无码AV天堂| 亚洲欧美成人久久综合中文网| 中国亚洲av第一精品| 影音先锋色小姐| 人妻系列无码专区久久五月天| 99re在线视频播放| 国产最新一区二区三区| 亚洲成人一区二区av| 亚洲精品国产精品乱码在线观看| 精品免费久久久久久久| 午夜丰满少妇性开放视频| 亚洲亚洲网站三级片在线| 国产不卡一区二区三区视频| 亚洲一区二区三区精品久久av| 少妇被黑人整得嗷嗷叫视频| 99久久超碰中文字幕伊人| 国产内射在线激情一区| 毛片一级精油按摩无码|