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

運載火箭自適應制導及在線軌跡重構方法研究

2023-03-15 02:04:12汪軼俊梁艷遷馬振偉史會濤
上海航天 2023年1期
關鍵詞:故障

汪軼俊,梁艷遷,周 鼎,馬振偉,史會濤

(1.上海航天技術研究院,上海 201109;2.上海宇航系統工程研究所,上海 201109;3.中山大學 系統科學與工程學院,廣東 廣州 510006)

0 引言

作為運載火箭的關鍵核心子系統,制導系統承擔著將有效載荷準確送入預定軌道的任務。在運載火箭“高可靠、高任務適應性和高經濟性”的演進目標指引下[1],伴隨行業態勢、基礎理論、硬件技術、系統架構的快速發展[2-4],制導方法正經歷著計算效率升級與創新突破。

本文回顧運載火箭傳統制導方法的發展歷程,并立足運載火箭主動適應和自主決策能力的增強需求[5],討論分析基于在線軌跡規劃的制導方法研究進展。在此基礎上,針對運載火箭剩余能量損失導致原目標軌道不可達的問題,提出一種任務降級在線軌跡規劃方法,并基于在線軌跡重構開展自適應制導方法的研究。

1 運載火箭制導方法發展歷程

1.1 傳統制導方法發展歷程

運載火箭制導方法的發展一直與硬件條件和基礎理論密切相關,導航設備、計算裝置的升級,以及控制理論與計算方法的創新持續支撐著制導方法的更新。

在導航設備和計算能力相對不足的條件下,我國早期運載火箭如長征一號配置陀螺和加速度計的捷聯導航系統,在模擬計算裝置下利用姿態角偏差與視加速度積分反饋的方式控制橫法向位置,對關機時間、起飛質量、偏差以及隨機干擾形成補償[6-7],這種制導方法是基于經典控制理論對特征量偏差的反饋干擾補償。

為滿足遠程火箭的發展對導航精度提升的需求,長征二號系列、風暴一號等運載火箭開始采用平臺+計算機的制導方案,平臺將測量的三個方向視速度增量轉換成電脈沖送入計算機依次積分,利用獲取的視速度、視位置、視位置積分項進行橫法向導引方程計算,通過與裝訂值的偏差反饋將火箭控制到標準彈道附近,并實現對關機特征量的控制[8-10],這就是所謂的隱式攝動制導,其在長征三號和長征四號系列運載火箭的早期任務中也得以沿用[11-13]。該方法在當時中等速度、小容量數字增量式箭載計算機條件下能夠完成相對高精度的制導[7]。

隨著箭載計算機技術的發展,運載火箭的速度、位置等信息可顯式地直接計算獲取,制導方法也變為顯式攝動制導[8]。長征四號和長征二號丁運載火箭先后開始采用全量型數字計算機[14-16],同時慣性導航開始采用動力調諧陀螺平臺[14-16]。在全新的慣性平臺和計算設備支撐下,運載火箭制導進入全數字化階段。

在基于經典控制理論偏差控制的攝動制導日趨完善的前提下,研究的重點開始轉向系統冗余和性能優化。硬件方面,長征系列運載火箭逐步開展了平臺+激光慣性主從冗余診斷[17]、箭載計算機三CPU 冗余診斷[18]、雙八表 捷聯慣性主從冗余診斷[19]、單十表捷聯慣性診斷[20]等增長可靠性的工作;而理論方法也隨著箭載計算機能力的發展開始向最優控制過渡。

由于新任務對運載火箭入軌精度、較大偏差適應性、可靠性等多方面的高要求,基于最優控制理論實時生成最優程序角的迭代制導方法成為當時一種較為合適的工程設計選擇[21-22]。迭代制導利用當前實時導航狀態與終端入軌約束,計算最優性能彈道,并給出程序角指令,其不需要依賴標稱軌跡且實時動態更新,具有較高的制導精度并能適應一定程度的推力下降等偏差。經過10 多年的工程型號應用,迭代制導方法已日臻成熟,并且根據不同的發射任務需求發展出考慮入軌姿態約束、滿足多終端約束、跨滑行段、預測修正等多種改進型[23-26],相關工程應用也正逐步推進。

1.2 基于在線軌跡規劃的制導方法概述

近年來,高性能箭載計算機、數值計算理論與方法的快速發展,為運載火箭的制導發展開辟了新的方向,即基于在線軌跡規劃進行計算制導與控制[3,27]。盡管迭代制導也是一種最優軌跡控制,但基于線性動力學邊值問題解析解的方式使得其無法考慮復雜狀態/控制約束,并且需要預先裝訂終端目標。航天發射任務的高可靠性要求對運載火箭的自主化、智能化和故障適應性提出了更高的要求,典型的場景就是運載火箭發生非致命故障導致剩余能力不足、載荷無法進入原目標軌道時,飛行控制能否實現重構以完成任務救援或減損,而對飛行軌跡和目標軌道進行在線規劃重構成為主要的備選方案之一,也是運載火箭制導自適應的主要途徑[28]。

2012 年10 月,“德爾塔4”在發射GPS-2F 衛星的過程中發生推力下降,火箭通過在線生成新的飛行軌跡,充分利用剩余推進劑成功完成衛星入軌任務;同月,SpaceX 公司的“獵鷹9”火箭在執行空間站貨運補給任務時,一臺發動機發生故障而關機,通過任務重構,另外8 臺發動機多工作近30 s 以彌補推力損失,成功完成飛船發射[1]。2016 年國際宇航大 會(IAC),俄羅斯專家指出2010 年12 月GLONASS 衛星發射失利中,若及時在線規劃基礎級與上面級新的交接班條件,上面級利用自身能力仍可將衛星送入預定軌道[5]。進一步表明在線任務重構對提高任務完成率、降低損失的重要性。

在線軌跡規劃的主要思路是考慮系統動力學及其他路徑約束、任務或階段終端約束,根據任務需求以某種性能指標最優作為條件,應用先進的問題處理策略及規劃算法,實現滿足計算速度和計算精度前提下的最優軌跡生成與制導控制指令輸出[29]。基于在線軌跡規劃制導的關鍵點為:①規定計算條件下的求解效率;② 算法求解精度;③偏差條件下算法的收斂性,即魯棒性[30]。當運載火箭發生動力系統故障需要在線軌跡重構(即在線軌跡優化)時,因不存在較為精確的優化初值,且須盡快完成飛行軌跡和制導指令生成,傳統制導方法不再適用,故亟須研究魯棒性強、計算效率高的在線軌跡重構方法以提升任務冗余度與可靠性。據此,國內外從模型約束處理、離散方法、求解算法開發、制導魯棒性增強等方面開展了相關研究[31-36]。

結合相關研究進展及前期工作[37],本文最終考慮基于偽譜離散與凸優化理論,提出一種運載火箭故障后剩余能力不足的彈道軌道聯合在線重構算法框架,并進行自適應制導方法的研究,為新一代運載火箭的自主化及可靠性增長提供基礎支撐和技術途徑參考。

2 運載火箭任務降級在線軌跡重構方法

本文主要討論運載火箭真空飛行段動力系統非致命故障致使有效載荷無法進入原目標軌道的情況,即發生任務降級。這種情況下,基于預先裝訂標稱軌跡或目標軌道的傳統制導策略與方法將不再適用,需要火箭能夠快速實現飛行彈道及目標軌道的在線規劃,進而重構飛行任務及制導以減小非致命故障帶來的損失。針對這類任務降級問題,本文基于偽譜離散理論與凸優化方法,提出一種彈道和目標軌道聯合在線規劃方法,期望開發一套泛化的算法框架適用于某型長征系列運載火箭多種軌道類型任務,實現簡潔高效的飛行軌跡重構。

2.1 任務降級規劃問題過程約束

任務降級軌跡規劃問題可描述為一個動力系統剩余能力、狀態及控制路徑、終端軌道根數關系約束下3 自由度動力學系統的最優控制問題。對于真空飛行段,系統動力學可通過變量代換構建成一個控制-仿射系統[37],即

式中:Vx、Vy、Vz和x、y、z分別為原入軌點軌道慣性坐標系(G系,原點位于地心O、y軸指向原入軌點、z軸垂直于軌道面、x軸形成的右手坐標系)下火箭的速度和位置分量;gx、gy、gz為引力加速度分量;T為故障后的發動機推力幅值;m為火箭質量(故障后可根據診斷和參數辨識確定秒耗量,進而確定質量變化規律,無需作為狀態變量);ux、uy、uz為發動機推力方向等效分量,其定義為

式中:?、ψ分別為俯仰角和偏航角。

該等效代換引入一個關于推力方向的附加約束:

相應地,δj(j=x,y,z)為推力方向分量的變化率,其滿足:

2.2 終端約束及性能指標

火箭狀態正常情況下,載荷入軌的終端約束可表示為狀態變量與標稱軌道根數之間的強非線性函數關系。然而,在動力系統非致命故障后剩余能力不足時,終端約束難以直接施加,需要構造新的終端約束和性能指標進行任務重構。

任務降級情況下,期望盡可能保證星箭留軌生存,即對軌道近地點提出了下限要求;保證生存的前提下,再考慮重構的新目標軌道與原目標軌道盡可能相似。在該需求下,一種自然的策略是“狀態觸發”多階段優化[38]。然而,復雜的優化邏輯分支及強非線性的根數加權指標不利于保證在線計算的收斂性和實時性??紤]對不同任務目標軌道的適應性,并兼顧在線快速求解的可靠性與簡潔性,在圓形目標軌道基礎上,本文提出一組松弛形式的終端約束為

式中:σh、σq∈R 為松弛變量;wh、wq∈R+為松弛項的懲罰因子;GM 為地球引力常數;xf、yf為終端位置分量;Vxf、Vyf為終端速度分量。

構建上述終端約束和性能指標組合式(5)和式(6),可從物理問題和數值求解兩個方面給出解釋:物理方面,如果火箭剩余能量足以形成圓形軌道,則松弛變量趨向于零,性能指標退化為圓軌道對應的入軌高度最高;如果剩余能量不足以形成圓軌道,終端約束與性能指標組合即等價于加權優化偏心率和半長軸,通過懲罰項夾逼可形成近圓的橢圓軌道,即形成了一套“能則入圓,否則成橢”的算法框架。數值求解方面,上述松弛-懲罰的約束和指標組合建模形式與序列凸化迭代方法中常用的“虛擬控制”措施相似,在早期迭代中虛擬控制可有效避免由凸化或線性化近似誤差引起的“偽不可行”現象,并且隨序列迭代下降不致影響最終的收斂解。

綜上所述,任務降級軌跡規劃的狀態變量為x=[Vx,Vy,Vz,x,y,z,ux,uy,uz]T∈R9,控制變量為u=[δx,δy,δz]T∈R3,相應的最優控制問題為

2.3 約束處理與分析

針對式(7)構建的任務降級彈道-軌道聯合規劃最優控制問題,本節分別對系統動力學、過程約束、終端約束和性能指標函數進行凸化處理,并進行相應的等價性分析。

首先,式(1)描述系統動力學通過變量代換轉化成了形如=f(x,t)+B·u的典型控制-仿射系統,序列迭代初始參考軌跡無需控制變量,實現了狀態與控制的解耦,利于算法收斂。故式(1)中僅剩引力加速度為非凸項,為保證較長飛行時間過程中的模型精度并兼顧算法的收斂性能,本文采用序列近似方法對引力加速度分量進行凸化處理,即

其次,變量代換引入的推力方向附加約束式(3)可通過無損凸化方法進行松弛,即

這樣處理相當于將“球殼”狀的非凸可行域凸化為“實心球”狀可行域,即二階錐約束,其等價性可由龐特里亞金極大值原理證明[39]。

最后,關于終端約束式(5)和性能指標式(6)中的非凸項,終端約束只需在終端時刻單點上得到滿足,即離散后每個終端約束只有1 維;性能指標無需嚴格滿足,其近似形式與原函數僅需在趨向性或指向性上等價與一致。

所以,本文對上述終端約束和性能指標進行序列線性化近似處理,終端約束式(5)可凸化為

對于性能指標式(6),通過線性化及引入非負輔助變量的方式,將二次指標函數轉換為便于凸優化算法求解的線性函數與二階錐約束的組合,即

式中:σJ為非負輔助變量,通過加權系數wJ進行懲罰。

在序列迭代過程中,當σh、σq和σJ全部收斂至零時,則凸化的性能指標與原指標式(6)同解,可保證降級任務軌道為圓軌道;若最終σh、σq和σJ不為零,說明當前狀態和故障條件下的剩余能力無法形成圓軌道,則獲得針對原問題的次優解,對應軌道高度與偏心率折中優化盡可能接近圓形的橢圓軌道。即可通過一套算法框架支撐“能則入圓,否則成橢”的在線重構策略。

此外,對上述終端約束和性能指標函數線性化近似處理后,為保證近似的精度和算法的收斂,引入相關變量二階錐形式的信賴域約束,即

式中:εx、εy、εVx和εVy為信賴域半徑。

綜上所述,通過無損凸化、序列近似、變量增廣及序列線性化措施,任務降級彈道-軌道聯合規劃問題可以被凸化成用于序列迭代的如下子問題:

進一步,利用Radau 偽譜法對問題式(15)進行離散化處理并排布成如下標準的二階錐規劃(SOCP)形式的子問題:

2.4 在線重構算法設計及數值試驗

結合上述討論,以某型長征系列運載火箭末級推力下降故障為例,本節給出基于偽譜離散和凸優化的任務降級彈道軌道聯合在線重構序列迭代算法,見表1。

表1 在線重構算法偽代碼Tab.1 Pseudo-code of the online reconstruction algorithm

針對本文涉及的真空段飛行工況,由于動力學系統非線性較弱,序列迭代算法對初始參考軌跡的要求相對較低。為提高算法在線執行的自主性和簡便性,所有狀態變量在離散點上的初始值設置為標稱軌跡在對應秒點上的狀態值,由于系統動力學為控制-仿射形式,無需給控制變量賦初值。

為驗證算法1 的正確性和快速收斂性能,同時兼顧算法的嵌入式平臺移植能力,本節給出任務降級在線重構的數值實驗與嵌入式平臺性能測試。

數值實驗采用標準C 語言,基于Visual Studio 2015 集成開發環境在桌面計算機(Intel-i7-10510U@1.8 GHz,內存16 GB)上進行。對象選取某型運載火箭末子級[37],原目標軌道為圓形標稱軌道,參考軌跡選取末子級飛行標稱彈道,偽譜離散配點數為30 個,序列迭代終止的無量綱門限為5×10-4?;诋斍八憷x取的運載火箭末子級彈道余量,考慮推力故障剩余比例為0.86~1.00,從推力比例為1 開始以0.02 為臺階逐步下調推力確定工況1~工況8,時間零點為三級起飛時刻,考核驗證算法對不同故障程度的適應性,以及收斂性能的一致性。不同故障條件下降級重構彈道優化結果在地心發慣系下的分量如圖1~圖6 所示。

圖1 X 方向不同剩余推力比例飛行距離結果對比Fig.1 Comparison of the flight distance results in the Xdirection under different residual thrust ratios

圖2 Y 方向不同剩余推力比例飛行距離結果對比Fig.2 Comparison of the flight distance results in the Ydirection under different residual thrust ratios

圖3 X 方向不同剩余推力比例飛行速度結果對比Fig.3 Comparison of the flight velocity results in the Xdirection under different residual thrust ratios

圖4 Y 方向不同剩余推力比例飛行速度結果對比Fig.4 Comparison of the flight velocity results in the Ydirection under different residual thrust ratios

圖5 不同剩余推力比例下俯仰角結果對比Fig.5 Comparison of the pitch angle results under different residual thrust ratios

圖6 不同剩余推力比例下偏航角結果對比Fig.6 Comparison of the yaw angle results under different residual thrust ratios

相應地,聯合重構獲得的不同故障條件下所形成的目標軌道如圖7 所示,具體軌道根數見表2。由圖7 和表1 可知,推力故障下降幅度越大,形成的降級軌道半長軸越低;當推力比例下降到0.86 時,降級軌道近地點已低于地表,物理上已不可行;當前狀態下火箭剩余能力足以形成圓軌道時,算法可以得到相應的最高圓軌道,在無法進入圓軌道的情況下,算法也能夠給出橢圓重構軌道,且從松弛變量的大小可以發現近地點高度接近最優。這表明,本文的算法對于不同的故障推力工況具有良好的適應性。

表2 不同剩余推力比例的重構軌道參數Tab.2 Reconfigured orbital parameters of different residual thrust ratios

此外,利用優化獲取的控制變量解驅動原始動力學進行積分,并與優化獲取的狀態變量解進行對比,可量化評估算法的求解精度,位置和速度的最終的積分精度以及序列迭代求解次數與計算耗時,見表3。可以看出,優化的入軌位置精度在0.6 km 以下,速度精度小于2 m/s,算法序列迭代6 次以內可收斂,求解耗時為百毫秒量級且相對穩定。這表明,本文的任務降級在線重構算法具備較高的精度、良好的求解穩定性與快速性,具備在線應用基礎。

表3 在線重構算法求解精度與性能Tab.3 Accuracy and performance of the online reconfiguration algorithm

為考核本文的任務降級在線軌跡重構算法對不同硬件平臺的適應性以支撐未來工程應用,進一步地在x86 架構和ARM 架構的嵌入式計算平臺上對算法進行計算耗時與內存占用等性能測試,相關測試結果見表4??梢钥闯?,隨著離散配點數的增加,即問題規模擴張,算法的內存占用基本呈線性增加趨勢,而計算耗時的增加幅度逐步變大,為后續工程研制提供了改進方向與參考;結合表3,針對同一工況,在精度滿足要求的同時,本文的算法在多個嵌入式平臺上運行具備1 s 以內收斂的能力,且迭代收斂次數對離散點數選取不敏感,具備良好的跨平臺求解穩定性,可以支撐基于在線軌跡重構的自適應制導。

表4 不同硬件平臺下性能對比測試Tab.4 Performance tests with different hardware platforms

3 基于在線軌跡重構的自適應制導方法

在線軌跡重構對于制導來說屬于基于模型的“一次性”前饋指令,即利用當前觸發時刻的系統狀態構造最優降級彈道和目標軌道,并生成程序角剖面用于跟蹤,或者直接以重構軌道為目標進行迭代制導計算直至到達關機條件。然而,考慮后續飛行過程中可能的故障狀態變化,如果模型信息不進行更新,動力學未建模及近似誤差、控制系統跟蹤誤差、故障診斷誤差等因素可能會引起制導精度不足乃至優化獲取的軌跡最終不可達,這對充分利用實時導航和故障診斷信息的自適應制導提出了需求,即需在線滾動執行軌跡重構并進行誤差修正。

3.1 滾動時域自適應制導算法設計

本節以任務降級彈道軌道聯合重構的偽譜-序列凸化算法為核心組件,將其嵌入滾動時域框架,基于在線軌跡重構提出一種解決火箭推力故障后任務降級的自適應制導方法?;谠诰€軌跡重構的自適應制導算法如下。

初始化在線軌跡優化算法所需參數,停止更新終端死區時間[Tbland];置制導周期索引[i=0]。

步驟1故障注入判斷。

如有故障注入,轉至步驟2;否則,使用標稱軌跡制導指令,轉至步驟5。

步驟2故障注入開始最優制導。

在故障注入時刻t0,采集火箭導航信息xR(t0)和故障信息;執行算法1 得到u*(t,xR(t0)),t∈[t0,tf];在制導周期[t0,t1]內,使用標稱軌跡制導指令量應用于系統;置i=1。

步驟3停止滾動時域更新判斷。

如果tf-ti≤Tbland,停止滾動時域更新,將u*(t,xR(ti-1))作用于系統直至關機,轉至步驟5;否則,轉至步驟4。

步驟4滾動時域軌跡優化。

到達制導周期觸發時刻ti時,判斷第(i-1)次軌跡優化結果是否可用;如可用,將u*(t,xR(ti-1))作用于系統,并采集火箭導航信息xR(ti)和故障信息,執行算法1 得到u*(t,xR(ti)),t∈[ti,tf];否則,將(i-2)次軌跡優化結果按照時序作用于系統,并采集火箭導航信息xR(ti)和故障信息,執行算法1 得到u*(t,xR(ti)),t∈[ti,tf]。

步驟5關機判斷。

如果未達到關機條件,轉至步驟1;如果達到關機條件,結束。

注意到,在上述算法的步驟3 中,進行了停止滾動時域軌跡優化更新的判斷,即在最后的飛行時間Tbland內,不再進行優化更新計算,在此將Tbland稱為“停止更新終端死區時間”。如此設置的目的在于防止在數值優化時域過短的情況下出現計算不穩定、不收斂的情況,防止臨近關機時間附近的復雜邏輯判斷。同時,在較短的死區Tbland之內,軌跡優化重構帶來的性能指標利益已經極小,系統誤差和故障檢測參數對制導軌跡可行性的影響規模也可忽略。

3.2 數值試驗與分析

本節選取與2.4 節相同的對象參數和數值實驗環境,設置滾動時域周期為1 s,指令程序角插值步長為0.01 s,仿真動力學模型積分步長為1 ms,停止更新終端死區時間為5 s,通過數值實驗對3.1 節算法的基本正確性和有效性進行驗證。

首先,在不施加偏差的條件下,通過不同故障時刻-剩余推力比例工況的閉環仿真實驗,對比分析單次軌跡重構前饋與滾動時域自適應制導的效果,相關結果見表5。可以發現,在無偏差情況下,單次重構前饋與自適應制導的入軌結果相近。一方面再次驗證了單次軌跡重構算法的較高計算精度;另一方面驗證了滾動時域自適應制導算法的正確性與有效性。

表5 滾動時域自適應制導與單次重構結果對比Tab.5 Comparison of the single reconfiguration and self-adaptive guidance results in the rolling time domain

最后,在閉環仿真的火箭動力學模型中加入推力偏差,而滾動時域制導算法中配置一種推力下降工況,考核基于在線軌跡重構的制導算法對故障診斷結果偏差的適應能力。分別考慮故障發生較早(第3 s)和較晚(第150 s)兩種情況,設置推力剩余比例為0.95,分別對動力學模型施加±2%、±3%和±5%的偏差,仿真的入軌結果見表6。結果表明:基于不準確的推力模型,滾動時域最優制導算法仍可以根據火箭實時狀態得到最優推力指令,且保持100%的遞歸可行性。在實際推力施加不同偏差的從小到大變化過程中,優化模型中的推力參數始終不變,入軌近地點高度卻隨之增大,這一結果充分驗證了本文的自適應制導算法利用系統實時狀態修正誤差的能力,反映了算法對模型失配和外部干擾較好的魯棒性,也從另一個角度表明動力系統故障診斷精度對制導結果的顯著影響。

表6 含有推力偏差的滾動時域自適應制導仿真結果Tab.6 Simulation results of the self-adaptive guidance in the rolling time domain with thrust bias

4 結束語

本文從運載火箭的可靠性和適應性提升需求出發,結合我國運載火箭發展歷程梳理了工程實踐制導方法的發展歷程,并基于運載火箭故障冗余任務重構問題,分析了基于在線軌跡規劃的制導方法研究進展。立足工程實踐,基于偽譜離散和凸優化理論,提出了一種運載火箭任務降級彈道軌道聯合在線重構的統一算法框架,據此給出了一種基于滾動時域軌跡重構的自適應制導技術途徑,并通過數值實驗和嵌入式測試,進行了正確性和有效性驗證。本文的研究可為制導技術的更新升級和運載火箭的自主化、智能化發展提供參考依據及基礎支撐。

猜你喜歡
故障
故障一點通
奔馳R320車ABS、ESP故障燈異常點亮
WKT型可控停車器及其故障處理
基于OpenMP的電力系統并行故障計算實現
電測與儀表(2016年5期)2016-04-22 01:13:50
故障一點通
故障一點通
故障一點通
故障一點通
故障一點通
江淮車故障3例
主站蜘蛛池模板: 国产不卡国语在线| 国产精品久久久久无码网站| 国产呦视频免费视频在线观看| 777午夜精品电影免费看| 亚洲欧洲一区二区三区| www精品久久| 蜜臀av性久久久久蜜臀aⅴ麻豆| 成人综合久久综合| 国产97区一区二区三区无码| 国产成人精品高清在线| 亚洲专区一区二区在线观看| 欧美曰批视频免费播放免费| 日本一本在线视频| 免费人成视频在线观看网站| 免费在线观看av| 成年免费在线观看| 欧美一区国产| 国产剧情一区二区| 99国产在线视频| 国产99精品久久| 免费看久久精品99| 91色综合综合热五月激情| aⅴ免费在线观看| 亚洲视频二| 在线观看国产精品日本不卡网| 国产毛片一区| 日韩成人在线网站| 91热爆在线| AV片亚洲国产男人的天堂| 98超碰在线观看| 日本爱爱精品一区二区| 精品久久久久久中文字幕女| 国产成人欧美| 久青草国产高清在线视频| 亚洲精品免费网站| 精品综合久久久久久97超人| 精品国产成人国产在线| 青青青草国产| 色网站在线免费观看| 久久99这里精品8国产| 国产免费怡红院视频| 久久黄色视频影| 国产日韩精品欧美一区喷| 久久香蕉国产线看观看式| 97免费在线观看视频| 热99精品视频| 亚洲av日韩av制服丝袜| 精品乱码久久久久久久| 亚洲天堂精品视频| 国产精品页| 亚洲综合色婷婷| 精品午夜国产福利观看| 午夜福利在线观看入口| 99热这里只有免费国产精品| a国产精品| 成人在线综合| 重口调教一区二区视频| 色悠久久久久久久综合网伊人| 97国产一区二区精品久久呦| 久久国产V一级毛多内射| 久久国产精品嫖妓| 日韩欧美国产三级| 国内精品自在自线视频香蕉| 国产精品久久久精品三级| 国产成人高清在线精品| 毛片免费试看| 国产欧美网站| 国产综合另类小说色区色噜噜| 免费女人18毛片a级毛片视频| 亚洲成a人片| 国产成人精品高清不卡在线 | 日本成人不卡视频| 好吊色国产欧美日韩免费观看| 欧美日韩在线国产| 一区二区理伦视频| 国产国模一区二区三区四区| 欧美福利在线播放| 青青青视频91在线 | 免费一级成人毛片| 成人夜夜嗨| 中文字幕调教一区二区视频| 又爽又大又黄a级毛片在线视频|