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

基于多模型網絡的激光末制導炮彈諸元解算方法

2023-10-07 01:49:28劉暢雷紅波林時堯范世鵬王江
兵工學報 2023年9期
關鍵詞:模型

劉暢, 雷紅波, 林時堯, 范世鵬*, 王江

(1.北京理工大學 宇航學院, 北京 100081; 2.北京理工大學 中國-阿聯酋智能無人系統“一帶一路”聯合實驗室, 北京 100081;3.中國兵器科學研究院, 北京 100089)

0 引言

半主動激光末制導炮彈是一種由火炮發射、在末端利用激光導引頭捕獲目標并對低速目標進行精確打擊的制導彈藥[1],由于集常規火炮的低成本與導彈的精確打擊能力等優點于一身,得到現代陸軍炮兵的廣泛應用[2]。為在末制導作用時捕獲目標,目前使用射表來確定射擊諸元,其中表尺(射角)、程裝(慣性陀螺解鎖時刻)、延時(激光照射器開始照射的時間)為基本諸元,修正諸元是對氣象偏差(包括氣壓、氣溫、風向、風速)等進行修正[3]。諸元解算是否合理,關系到分系統能否正常工作(能否成功捕獲目標、有限修正能力能否滿足需求以及慣性陀螺是否發生碰框),直接影響彈藥的命中概率。

針對傳統射表應用插值法求解射擊諸元精度不高、解算速度慢、無法準確處理精準的分層實時氣象信息等問題,學者們在射表編制理論、射表數據采集、射表射擊試驗等方面做出了諸多貢獻,提高了命中概率。趙明陽等[4]為精確計算氣象對諸元的影響,應用貪心優化算法對氣象修正表自動分區并計算擬合系數,逼近函數系數精確度較高,實現了自動分區擬合的功能。褚進等[5]對高原條件下某型末制導炮彈彈道進行仿真分析,對比了不同諸元對名義彈道及捕獲域的特點,為高原射表編擬提供了參考。丁天寶等[6]為改進以實彈射擊為主的射表編制方法,提出了高原與平原交接階梯抑制方法,減小了平原與高原射表在高度交接段的原理誤差,適用于寬海拔作戰任務。陳瑞軍等[7]針對基本諸元解算速度慢、不能滿足有效攻擊區的戰術指標、命中率低的問題,提出了簡化的彈道模型并應用攻擊區中心為表載射程的解算方法,提高了制導彈藥射表的編制精度。賈波等[8]對氣動參數進行靈敏度分析并進行辨識,給出了基于彈道精確測量的射表編擬方法。

然而,目前基于傳統射表的諸元解算存在以下問題:

1)射擊諸元與彈道參數存在強非線性的映射關系,本質上是一個復雜的多因素耦合問題。目前通過線性插值法求取,會帶來較大的截斷誤差。

2)若將實時獲取的分層氣象信息作為輸入,則其維度大大增加,射表規模隨之顯著增大,且插值計算量較大,難以滿足快速計算需求。

3)目前雖有一些基于傳統射表的射擊諸元計算軟件,但忽略諸元相關的因素較多,求解精度差。

隨著人工智能技術的不斷發展,諸多學者將其應用在彈道解算的問題中。文獻[9]應用深度神經網絡,將升力、阻力系數等作為輸入,將預期落點作為輸出,解決了故障條件下高超飛行器容錯制導問題。文獻[10]為解決空空導彈無法同時解算多種攻擊區的規劃及解算問題,提出了深度擬合網絡模型及訓練策略,相對于傳統方法提高了解算的準確性及實時性。文獻[11]針對彈道軌跡預測精度不足的問題,構建混合輸出單元,應用增強上下文信息長短期記憶模型進行預測,結果表明該方法相對于傳統方法在實時性及泛化能力上均有較大優勢。文獻[12]應用卷積及長短期記憶混合神經網絡,解決了彈丸非線性軌跡預測問題;應用龍格-庫塔法解算6自由度彈道模型得到大量樣本,應用混合網絡模型對軌跡進行預測,結果表明在預測精度上優于單一的網絡模型,為軌跡預測提供了一定的參考。

本文針對快速精確求取射擊諸元的實際需求,利用神經網絡在表達映射關系及求解非線性擬合問題上的優勢,提出一種基于多模型網絡的諸元解算方法,使諸元解算求解精度更高。該諸元解算方法的主要創新點為:1)可對實時氣象進行分析,并將其作為多模型網絡的輸入,相對于僅考慮地面氣象或標準氣象,精度更高;2)將諸元解算這類多因素耦合問題進行功能解耦,使多模型網絡中的個體僅需實現單一功能,改善了各網絡的可解釋性與求解精度,大大降低了網絡優化與設計的復雜度。

1 彈道模型及氣象分析

1.1 激光末制導炮彈彈道模型

激光末制導炮彈由無控段、程控段及末制導段組成,程裝及延時是各段之間轉換的重要標志。若炮彈僅以無控及程控段飛行,則此時的落點稱為名義落點。激光末制導炮彈彈道如圖1所示。

圖1 飛行軌跡圖

對于彈道軌跡,通常采用非線性6自由度模型描述彈道[13],氣動特性由風洞實驗獲取。非線性6自由度彈道模型通常由質心動力學方程、繞質心轉動動力學方程、質心運動學方程、繞質心轉動運動學方程及幾何關系等構成[14]。建立的微分方程組[15]如下:

(1)

(2)

(3)

(4)

式中:m為彈丸實時質量;vm、θ、ψv分別表示彈丸速度、彈道傾角及彈道偏角;wx、wy、wz表示速度在彈體坐標系下的轉動角速度分量;x、y、z表示彈丸在地面坐標系下的3個位置分量;?、ψ、γ表示俯仰角、偏航角及滾轉角;p為彈丸所受的推力;α、β分別為攻角和側滑角;Fx、Fy、Fz分別為彈丸受到的阻力、升力及側向力;g為重力加速度;γv為速度傾斜角;Jx、Jy、Jz分別為彈丸在彈體坐標系下沿3個坐標軸的轉動慣量;Mx、My、Mz分別為彈丸在彈體坐標系下沿3個坐標軸的力矩。相關坐標系、力和力矩的定義見文獻[16]。其中,角度幾何關系如式(5)所示:

(5)

本文暫不考慮末制導段,后續建立樣本時僅以名義落點作為期望射程,其原因為:1)當名義落點與實際期望落點相差較小時,依靠彈丸的過載能力可對較小的偏差進行修正,精確命中目標;2)若考慮末制導段,則由于彈丸具有一定的過載能力,會出現基本諸元對應多個射程的情況,此時會出現奇異;3)射表編擬時,同樣采用名義落點作為期望射程。因此,合適的基本諸元(表尺、程裝)對命中目標有著關鍵性的作用。

當給定程裝、表尺、初速等初始條件時,應用4階龍格-庫塔法對式(1)~式(4)的6自由度微分方程組進行求解,可得到對應的彈道軌跡。其余條件不變,表尺及程裝的變化對彈道軌跡的影響如圖2 所示。

圖2 彈道軌跡

由圖2可見:軌跡1和軌跡2的程裝相同,軌跡2的表尺大于軌跡1,可以看出軌跡2的射程大于軌跡1;軌跡1和軌跡3的表尺相同,軌跡1的程裝大于軌跡3,軌跡3的射程大于軌跡1;表尺和程裝對射程的影響較大。表尺即為微分方程組的初始彈道傾角θ0;程裝即慣性制導段的初始時刻Tcz,此時彈上慣性系統開始工作,通過重力補償使彈丸基本按照直線飛行,其目的為增大射程及以較小的速度傾角進入末制導。由此可知升力Y變化為

(6)

由式(6)可以看出,程裝對微分方程中的升力影響較大,可見程裝變化對射程影響較大。

圖2中軌跡1與軌跡4射程相近,但二者表尺和程裝大小均不相同,因此無法單一從某個參數判斷諸元對射程的影響,表明諸元解算是一個多維度耦合問題。

1.2 氣象條件分析

氣象條件的精確分析對彈丸的射擊精度有著決定性影響,因此對氣象條件進行充分的分析是非常必要的。其中,氣壓、氣溫、風向及風速對彈道的影響較大,下面對其進行逐一分析。

1.2.1 氣溫分析

氣溫及氣壓主要對馬赫數有很大影響,式(1)~式(4)中的升力系數等均與馬赫數有關。目前,都是通過簡易或精密分析方式處理氣溫。簡易法是假定當前環境符合標準氣象環境,即氣溫按線性變化,變化趨勢如式(7)所示:

τ=τ0-Gy

(7)

式中:τ0為地面氣溫,τ0=299.8 K;G為氣溫變化系數,G=6.328×10-3。從式(7)中可以看出,在采用上述方法進行氣溫簡易計算時,需要假設氣溫變化趨勢是恒定的,但實際中氣象情況復雜,標準氣象的線性變化假設與實際情況有較大差距。精密氣象是測量不同彈道層下的氣溫,雖然此方法精度較高,但計算較為復雜,且每層間的溫度按線性插值求得。若要對氣溫進行精確分析,則需使用實際氣溫。以我國某地區近4年6~7月份為例,其溫度變化曲線如圖3所示。

圖3 氣溫變化趨勢

由圖3可以看出,相同地區相近月份的氣溫變化趨勢是類似的,因此可以將相近月份的氣溫進行非線性擬合。應用此方法不僅可以對實際氣溫進行分析,得到比線性方法更高的擬合精度,而且計算量相比于傳統高精度計算方法計算量更小。

1.2.2 氣壓分析

與氣溫分析類似,以相同地區近4年6~7月份為例,其溫度變化曲線如圖4所示。

圖4 氣壓變化趨勢

由圖4可以看出,相同地區相近月份的變化趨勢一致性較好,同樣可以對氣壓進行擬合處理。

1.2.3 風向及風速分析

風對彈道的影響主要是引起附加攻角Δα及側滑角Δβ,使得在求解彈道軌跡時,在攻角和側滑角的基礎上還需要考慮附加角度。傳統射表不考慮風的影響,然而實際射擊中,風對方向及射程的影響較大。為便于分析修正,通常采用彈道平均值的方法,即求得一個恒定的平均風代替實際風。由于風在不同彈道層的影響不同,各層風所賦的權重(層權)不同。根據風向和風速可將風分解為橫風(垂直于彈道方向)及縱風(平行于彈道方向)。兩者求取層權的方式相似,因此以縱風為例進行求取層權的說明。

設彈道高分為n層,首先假設全彈道均有相同的平均縱風,計算出相對于無風狀態下的射程偏差xw;其次計算出第1層無風狀態下的射程偏差xw1,此時風偏在第1層所起作用x1大小為

x1=xw-xw1

(8)

依次類推,可求出每層所起作用大小xi,進而可得每層的層權為

Qi=xi/xw

(9)

式中:Qi為每層的層權。于是可得彈道縱風為

(10)

式中:xwind為恒定的平均縱風;xi,wind為每層的實際縱風。若要求取準確的層權,則不僅需要考慮射程,還需要考慮初始速度、初始彈道傾角、氣動參數等,而傳統方法中的多維插值計算量較大,難以滿足快速計算需求;若計算時僅考慮射程影響,則無法求取精確的層權。因此,如何在多維參數影響下求取層權,對分析風的影響是至關重要的。

2 基于多網絡的諸元解算模型

2.1 多模型網絡分析

由第1節分析可知,所需解決的問題為:

1)如何對氣溫、氣壓進行擬合。擬合完畢后,僅將參數作為網絡模型的輸入即可,在保證精度的同時減小計算量。

2)不考慮風的影響,如何求解基本諸元。將氣溫、氣壓的擬合參數及期望射程作為網絡模型的輸入,可輸出表尺及程裝。

3)當有風干擾時,在考慮多方面因素下,如何快速準確地求出彈道的平均風(包括橫風及縱風)。

4)求出平均風后,如何對基本諸元修正,獲取實際諸元。

通過第1節分析可知,諸元解算作為一個耦合問題,其輸入輸出之間存在一定的映射關系,神經網絡在解決映射及非線性擬合問題上有明顯的優勢,因此可應用神經網絡解決此問題。另一方面,通過上述分析可知,各個問題之間存在一定的耦合關系,無法將所有問題通過單一的網絡求解,因此本文提出一種基于多模型網絡的諸元解算方法。本文求解框架如圖5所示。

圖5 模型整體框架

由圖5可以看出,多模型網絡的流程為:

步驟1將實際氣象中的氣壓、氣溫進行擬合,將擬合參數、期望射程等作為諸元解算神經網絡的輸入,輸出為對應的表尺、程裝、彈道高度。建立樣本后,對該網絡進行訓練。

步驟2將諸元解算神經網絡輸出的無風狀態下的表尺及彈道高度作為縱風層權神經網絡的輸入,輸出為對應的層權,再將實際的風向、風速與層權按式(10)進行計算,求得平均彈道縱風。

步驟3與步驟2類似,可求得平均彈道橫風,再根據橫風及射程大小得到方向修正量。

步驟4將得到的平均彈道縱風作為諸元修正網絡的輸入,將表尺的變化量作為輸出,將變化量與步驟1得到的表尺相加,得到實際表尺。由于表尺、程裝之間存在耦合,僅考慮變化表尺。

根據上述步驟,得到的表尺、程裝及方向修正量即為諸元解算的最終結果。

2.2 樣本建立

各網絡模型的樣本建立類似,因此以諸元解算為例進行說明。由于實際飛行實驗數據較少,無法構建大量樣本,可根據實際飛行數據對仿真中的氣動參數進行擬合。應用擬合后的參數進行仿真,得到大量樣本,此時獲得的樣本與實際飛行試驗數據相近,諸元計算精度高。

2.2.1最小二乘法及樣本建立

由1.2.1節及1.2.2節分析可知,相同地區相近月份氣溫、氣壓的一致性較好,因此可將相近月份如每間隔兩個月進行參數擬合。以氣壓為例,將測得的實際氣壓作為若干離散點,坐標分別為(y1,p1), (y2,p2),…,(yi,pi),…,(yn,pn),通過最小二乘法中的線性回歸法[17]擬合。線性回歸法主要是通過最小化誤差的平方和,尋找最佳的函數匹配。對實際氣壓的誤差εi及殘差i的求取如式(11)、式(12)所示:

(11)

i=pi-i

(12)

(13)

擬合函數形式如式(14)所示:

g=0.989 8e-0.123y

(14)

根據式(14)進行擬合后,J2即誤差平方和為0.003 026,R-square為0.998 9,由此可看出應用式(14)的擬合效果較好,因此僅需將[0.989 8,-0.123]T作為氣壓輸入即可,其余月份方法類似。氣溫的擬合方法相同,不再贅述。

本文以某型激光末制導炮彈為例,其理論初速與裝藥號有關,為相應定值。因此以初速為585 m/s為例,初始參數表尺(單位為密位,6 000 mil為360°)和程裝(單位為分劃,1分劃=0.2 s)的采樣范圍如表1所示。

表1 參數的采樣范圍

其余初始參數不變,為實際發射條件,此時環境為無風狀態。根據不同初始條件,對式(1)~式(4)的6自由度彈道方程進行求解,即可得到對應的彈道軌跡。氣動參數已經根據實際飛行數據進行了擬合,可保證仿真結果與實際飛行彈道相近。通過表1的采樣間隔共生成1 050條彈道,實現對實際飛行軌跡的高泛化度覆蓋。

2.2.2 網絡訓練

神經網絡是一種具有很強非線性擬合能力的數學方法,可以解決一些傳統推理模型無法解決的問題[18]。通過提取、學習和訓練數據中的隱藏模式,并調整每個節點的權重,神經網絡的輸出值趨于預期輸出。本文采用深度神經網絡[19](DNN)對模型進行訓練。由于本文采用多模型網絡進行擬合,每個網絡僅需要解決單一簡單問題,大大降低了網絡設計復雜度,因此本文設計的DNN僅使用2個全連接層作為隱藏層,每層有10個神經元。以諸元解算神經網絡為例,其結構如圖7所示。

圖7 DNN網絡結構

由圖7可知,本文所設計的DNN包括輸入層、隱藏層和輸出層,其中,輸出層為表尺、程裝、最大彈道高度,輸入層數據為期望射程、氣溫及氣壓擬合參數。兩個相鄰層之間的數據傳輸以權重的形式進行。

神經元的數學模型如式(15)所示:

(15)

式中:σ為神經元的輸出;f(·)為激活函數;di為第i個神經元的輸入;ωi、δ為對應神經元的權重及閾值。本文采用ReLU作為輸入層和隱層的激活函數,具有簡化計算過程、防止梯度消失及訓練網絡快速等優勢,輸出層不設置激活函數,其表達式為

f(d)=max(0,d)

(16)

式中:d為當前神經元的輸入。

本文選擇表尺、程裝、最大彈道高度等期望輸出值與實際輸出值的最小均方差作為損失函數,如式(17)所示:

(17)

(18)

對于式(18),根據鏈式法,有

(19)

(20)

DNN的訓練過程實際上就是通過上述反向傳播和梯度下降算法不斷優化網絡參數,使誤差平方和最小化。當誤差達到預期時停止對式(20)更新,完成訓練。采用ADAM優化器[20]對E(σ,)進行優化,優化后的參數為

(21)

式(18)和式(21)中的學習率為訓練預設初值,在DNN訓練過程中ADAM優化器會自適應調節該參數。

本文訓練過程基于TensorFlow-1.13.0-GPU版構建DNN模型,顯卡型號為RTX3060 12G,學習率選取0.001。以諸元解算神經網絡為例,根據2.2.1節中得到的1 050條彈道,按照92%的數據作為訓練集、4%的數據作為測試集、4%的數據作為驗證集的比例,對網絡進行離線訓練、測試及驗證。當網絡精度達到預設要求時訓練結束。

3 基于多網絡的諸元解算模型驗證

為驗證多模型網絡的精度,本節給出諸元解算的數值仿真實驗。首先根據驗證集的效果論證多模型網絡的精度;其次選取3組不同的期望射程及實際氣象,將其輸入模型中進行驗證;再次進行蒙特卡洛仿真實驗,進而論證模型的精度及魯棒性;然后與射表解算出的諸元進行對比,論證模型求解精度更高;最后,通過實際打靶論證本文提出的模型在實際飛行試驗中同樣適用。

3.1 多模型網絡精度

根據第2節1 050條彈道中4%、4%的數據作為驗證集、測試集,即分別有42條彈道作為驗證及測試,訓練中損失函數的變化情況如圖8所示。

圖8 損失函數變化

設定訓練1 000輪或損失函數值小于0.000 01為終止條件,由圖8可以看出,當訓練輪數到453時,已經達到期望的損失函數值。通過上述變化情況可以看出,DNN損失值下降迅速,在第60輪前下降較為迅速,后續變化較為平穩,因此可以說明網絡已收斂。使用測試集測試訓練好的神經網絡,以諸元解算模型中的程裝為例,其效果如圖9所示。

圖9 測試集精度

由圖9可以看出,模型預測的程裝值與實際程裝值相近,絕對誤差為0.031 9,均方誤差為0.002 9,由此可見訓練完畢的諸元解算模型精度較高。

3.2 多模型網絡驗證

為方便分析驗證,初速選取585 m/s,其余仿真條件如表2所示。

表2 仿真條件

3.2.1 驗證諸元解算網絡模型

在暫時不考慮風的情況下,將3組期望射程及氣溫、氣壓擬合的參數作為模型輸入,對模型輸出得到表尺及程裝進行驗證,結果如圖10所示。

圖10 諸元解算模型輸出結果

由圖10可見,模型輸出射程與期望射程相近,具體誤差及精度如表3所示。

表3 誤差與精度

由表3可知,3組期望射程下誤差值在20 m附近,誤差率在0.1%附近,由此可進一步論證諸元解算網絡的精度較高。

3.2.2 驗證縱風及橫風修正網絡

以縱風修正網絡為例,將諸元修正網絡輸出的彈道高度及期望射程作為輸入,輸出為縱風層權。如期望射程16.3 km時,將彈道高度按400 m分層,分為10層,縱風層權神經網絡輸出的層權為[0.035 294,0.070 588,0.070 588,0.082 352,0.129 411,0.141 176,0.105 882,0.2,0.164 705,0],將縱風層權及實際風速代入式(10)中,可得平均縱風為7.855 6 m,平均橫風獲取方式相同。

3.2.3 輸出修正后的諸元

將平均橫風、縱風與期望射程輸入諸元修正網絡,經過修正后的諸元結果如圖11及表4、表5所示。

表4 期望射程的誤差與精度

表5 橫偏的誤差與精度

圖11 多模型網絡輸出結果

由圖11可以看出,多模型網絡的輸出射程、橫偏與期望射程、橫偏相近。

由表4、表5可知,修正后的射程誤差值仍然在20 m附近,誤差率約為0.15%,修正后的橫偏誤差值小于1 m,因此可證明本文設計的多模型網絡精度較高。

3.3 蒙特卡洛實驗

為進一步驗證多模型網絡的精度,進行蒙特卡洛仿真實驗,初始條件與3.2節相同,選擇第2組的實際氣象,期望射程從[13 km,18 km]中隨機選取。分別重復100次蒙特卡洛實驗,繪制各次實驗的射程誤差,結果如圖12所示。

圖12 蒙特卡洛仿真實驗結果

由圖12可知,誤差最大值為29.6 m,方差為9 m,最大相對誤差不到0.1%,本文設計的多模型網絡誤差較小,精度較高,可滿足發射要求。100次 蒙特卡洛實驗射表解算的平均時間為 46 ms,本文方法解算時間為34 ms。由此可見,本文方法不僅在精度上有一定的提升,同時也可實現快速解算。

3.4 仿真對比驗證

為對比驗證基于多模型網絡與射表解算的精度,初始條件與3.2節相同,期望射程為15.3 km及16.3 km,二者關于諸元的求解結果如圖13及表6所示。

表6 誤差對比

圖13 仿真結果對比

由圖13可知,相對于射表解算,多模型網絡更接近期望射程,具體誤差及精度如表6所示。

由表6可知,兩組期望射程下多模型網絡相對于射表解算精度更高,可提高100 m。在[13 km,18 km]射程范圍內每間隔500 m做上述實驗,多模型網絡的誤差方差為21.7 m,射表解算的誤差方差為117.6 m。由此可見,基于多模型網絡求取的諸元相對于射表解算精度更高。

3.5 實際飛行驗證

為驗證多模型網絡在實際飛行實驗中的精度,2022年7月,對多模型網絡進行飛行試驗論證。此次飛行試驗期望射程為16.00 km,將實際氣象輸入到多模型網絡中,網絡得到的實際射擊諸元如表7所示。

表7 實際射擊諸元

根據上述諸元與射表解算進行理論仿真驗證,誤差如表8所示。

通過上述仿真驗證可以看出,多模型網絡相對于射表解算與期望射程更為接近。采用上述諸元進行的兩次飛行試驗均中靶,如圖14所示。

圖14 飛行試驗效果

由圖14可知,兩次飛行均命中,驗證了基于多模型網絡的諸元解算方法的有效性。

4 結論

本文針對基于射表解算射擊諸元精度差的問題,提出了一種基于多模型神經網絡的諸元解算模型。首先建立了末制導炮彈6自由度彈道模型,利用數學仿真手段構建訓練、測試及驗證樣本。其次分析了氣象條件,對氣溫、氣壓進行擬合,對風引入層權求解。再次,對多模型網絡中每個DNN的訓練樣本進行訓練。得到主要結論如下:

1)訓練完畢的多模型網絡在不同期望射程下都表現出較高精度,在100次蒙特卡洛試驗中誤差較小。

2)通過仿真對比可知,多模型網絡相對于射表解算出的諸元精度高,進一步證明了本文方法的優勢。

3)在實際飛行試驗中兩次均中靶,驗證了多模型網絡在實際飛行試驗中的有效性。本文提出方法通用性良好,可推廣到各類精確制導武器。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 日本高清在线看免费观看| 成年人视频一区二区| 精品视频一区在线观看| 久久久久免费看成人影片 | 国产在线拍偷自揄观看视频网站| 手机精品福利在线观看| 中国精品久久| 91亚洲视频下载| 亚洲成aⅴ人在线观看| 久久夜夜视频| 亚洲午夜天堂| 精品无码人妻一区二区| 欧美激情第一区| 成人亚洲视频| 国产免费黄| 国产va免费精品观看| 久久久久久久久18禁秘| 九九热精品视频在线| 亚洲浓毛av| 国产香蕉在线视频| 五月六月伊人狠狠丁香网| 日本不卡免费高清视频| 国产日韩精品欧美一区灰| 亚洲综合精品第一页| 国产美女在线免费观看| 国产美女91呻吟求| 成人福利在线观看| 免费看黄片一区二区三区| 午夜综合网| 狠狠躁天天躁夜夜躁婷婷| 毛片久久网站小视频| 久久a毛片| 国产第一页第二页| 国产精品欧美激情| www欧美在线观看| 中国国产一级毛片| 99精品在线看| 国产精品久久久久久搜索| 日韩无码视频网站| 国产成人精品在线| a级毛片视频免费观看| 五月婷婷综合色| 欧美影院久久| 色呦呦手机在线精品| 国产免费高清无需播放器| 91网址在线播放| 日韩欧美中文字幕一本| 日韩精品一区二区深田咏美| 99久久精品国产自免费| 亚洲欧美日韩另类在线一| 亚洲精品无码不卡在线播放| 国产草草影院18成年视频| 999国产精品永久免费视频精品久久 | 久久综合九九亚洲一区| 欧美怡红院视频一区二区三区| 欧美福利在线播放| 18禁黄无遮挡网站| 男人天堂伊人网| 无码国内精品人妻少妇蜜桃视频| 中文字幕亚洲另类天堂| 日本人真淫视频一区二区三区| 亚洲三级视频在线观看| a级毛片视频免费观看| 国产第一页亚洲| 中文字幕乱妇无码AV在线| 日本一区二区不卡视频| 免费视频在线2021入口| 成人免费视频一区二区三区| 久久 午夜福利 张柏芝| 日韩色图区| 国产全黄a一级毛片| 久热这里只有精品6| 成人av手机在线观看| 欧美日本视频在线观看| 免费在线色| 亚洲精品无码AⅤ片青青在线观看| 欧美国产成人在线| 最新国语自产精品视频在| 色噜噜狠狠狠综合曰曰曰| 久久成人国产精品免费软件| 日韩第九页| 国产一区二区免费播放|