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

一種低截獲背景下的集中式MIMO雷達快速功率分配算法

2023-07-04 09:51:46李正杰謝軍偉張浩為
雷達學報 2023年3期
關鍵詞:分配模型

李正杰 謝軍偉 張浩為* 溫 泉 劉 斌

①(空軍工程大學防空反導學院 西安 710051)

②(國防科技大學信息通信學院 西安 710106)

③(國防大學聯合作戰學院 北京 100091)

1 引言

作為一種新體制雷達,集中式多輸入多輸出(Multiple-Input Multiple-Output,MIMO)雷達引起了研究人員的廣泛關注。經過近20年的研究,出現了許多關于參數估計、主動抗干擾、射頻隱身等方面的研究成果,并已進入工程實踐階段[1]。相比于傳統相控陣雷達全相干的信號發射機制,集中式MIMO雷達可同時發射多個正交信號,從而具備優越的波形分集增益[2]、低截獲性能[3]、虛擬孔徑[4]以及目標檢測性能[5]。為充分發揮多發多收的體制優勢,進一步挖掘集中式MIMO雷達的作戰潛能,需要對有限的雷達系統資源進行優化分配。

近年來,圍繞資源分配技術在目標跟蹤上的應用問題,國內外學者基于認知雷達思想[6],逐漸形成了關于雷達資源分配的兩大原則[7]:(1)有限資源約束下提高跟蹤性能原則[8–16];(2)滿足跟蹤精度約束下減少資源消耗原則[17–23]。關于提高跟蹤性能問題,文獻[8]研究了集中式MIMO雷達在同時多波束工作模式下的資源分配問題,并證明了通過優化分配發射波束和功率可有效提高目標跟蹤精度。文獻[9]通過推導后驗克拉默-拉奧下界(Posterior Cramer-Rao Lower Bound,PCRLB),提出了一種集中式MIMO雷達同時多波束模式下的功率-帶寬聯合分配算法。文獻[10]推導了理想檢測條件下的預測條件克拉默-拉奧下界(Predicted Conditional CRLB,PC-CRLB),并研究了分布式雷達網絡中的節點選擇-功率分配問題。在文獻[8,9]的基礎上,文獻[11]進一步研究了集中式MIMO雷達多目標跟蹤中的功率-帶寬-波束聯合分配問題。文獻[12,13]對雜波環境下的跟蹤性能下界進行推導,并研究了多雷達系統中的功率分配問題以及子陣選擇-功率聯合分配問題。為進一步提高多目標跟蹤場景下的跟蹤質量,文獻[14]基于服務質量(Quality of Service,QoS)準則,提出了一種波束-功率聯合分配算法,以使各目標的跟蹤誤差逐漸逼近期望值。為提升滿足跟蹤精度要求的目標數量,文獻[15]提出了一種基于目標容量的同時多波束功率分配算法。文獻[16]考慮使目標跟蹤精度和滿足精度要求的目標容量同時最大化,研究了壓制式干擾場景下的波束選擇-功率分配問題。在減少發射資源消耗方面,文獻[17]提出一種面向射頻隱身的發射節點-駐留時間-功率聯合分配算法,該算法在滿足跟蹤精度要求的前提下盡可能減少對組網雷達系統的射頻資源消耗。針對機載組網雷達單目標跟蹤場景,文獻[18]提出一種基于低截獲概率(Low Probability of Intercept,LPI)性能準則的資源分配算法。文獻[19]針對集中式MIMO雷達多目標跟蹤問題,提出一種基于QoS準則的穩健功率分配算法。文獻[20]建立了集中式MIMO雷達體制下的空-時資源分配和波形選擇優化模型,在保證目標有效檢測的前提下,對系統資源消耗量和跟蹤性能同時進行優化。針對機會約束規劃的機會陣雷達網絡系統,文獻[21]研究了單批機動目標跟蹤下的發射功率最小化問題。文獻[22]同時考慮目標跟蹤精度和低截獲概率(LPI),提出一種相控陣雷達網絡多目標跟蹤中的目標指派和資源聯合分配算法。文獻[23]進一步提出一種基于LPI的功率-駐留時間-帶寬-脈沖長度多參數聯合分配算法。

上述工作為MIMO雷達資源分配問題的研究打下了良好基礎,但仍存在一些不足之處。如當前對于機動目標跟蹤的資源分配研究成果較少,且針對不同運動特性和電磁特性目標的威脅度評估問題缺乏應有的關注。此外,對于多目標跟蹤中同時考慮跟蹤精度和低截獲性能的資源分配研究成果也十分有限。基于上述情況,本文針對集中式MIMO雷達多機動目標跟蹤場景,提出了一種低截獲背景下的快速功率分配算法。本算法基于同時多波束工作模式,如圖1所示,雷達同時發射多個寬波束對空間進行監視,并在接收機中形成窄波束以提取目標信息。此時,每個目標都被一個波束獨立進行跟蹤,通過合理地向各發射波束分配功率,可提高資源利用率和雷達跟蹤性能。

圖1 集中式MIMO雷達同時多波束工作模式Fig.1 Simultaneous multi-beam working mode of the collocated MIMO radar

具體而言,本文首先將目標運動模型構建為自適應當前統計(Adaptive Current Statistical,ACS)模型。然后推導了PC-CRLB,并將其作為跟蹤性能下界。針對各目標運動特性和電磁特性的差別,構建了目標威脅評估模型。其次,分別計算了基于PC-CRLB的目標跟蹤誤差評估指數和雷達截獲概率,并通過引入線性化機制構建了關于發射功率的優化模型。最后,針對優化模型的凸性質和單調性,采用了一種低復雜度的基于序列松弛的求解算法進行快速求解。仿真結果驗證了所提算法的可行性和時效性。本文構建的基于發射功率分配的認知跟蹤系統如圖2所示。

圖2 認知跟蹤系統示意圖Fig.2 Schematic diagram of cognitive tracking system

2 系統建模

假設某集中式MIMO雷達位于x-y平面內,雷達采取均勻線陣排布,陣元數為M×N。將第m個發射陣元和第n個接收陣元分別表示為Tm和Rn,其中?m∈{1,2,...,M},?n∈{1,2,...,N},且對應陣元間距分別為dt和dr,雷達陣列模型如圖3所示。

圖3 集中式MIMO雷達陣列模型Fig.3 Collocated MIMO radar array model

為簡化模型,給出如下假設:

(1) 系統采用同時多波束工作模式,各發射陣元發射窄帶正交信號;

(2) 已知目標空域內存在Q個運動目標,各目標初始位置為先驗信息;

(3) 考慮遠場條件,認為目標散射系數對各發射、接收陣元相同;

(4) 為穩定系統內部各器件工作負荷,將各發射波束的功率平均分配給對應發射陣元。

2.1 信號模型

假設在k時刻,該集中式MIMO雷達向第q個目標發射的波形為

其中,Pk,q為發射功率;Ek,q(t)為歸一化的發射信號復包絡;fc為載波頻率。該信號的有效帶寬βk,q和有效時寬Tk,q分別滿足

以及

經目標反射后,接收信號可表示為

其中,ak,q為目標q的反射復增益;γk,q為信號衰減系數,與徑向距離的4次方成反比;τk,q為信號時延;為信號多普勒頻移;nk,q(t)為零均值的復高斯白噪聲。

2.2 運動模型

當前用于描述目標機動特性的運動模型主要分為單模型和多模型兩類[24]。由于單模型算法無須建立復雜的模型集來匹配目標機動,因此更具時效性。當前統計(Current statistical,CS)模型是一種比較切合實際的單模型,通過構建非零均值的加速度模型能較全面地描述目標機動特性[25],但存在自適應性不足的問題,需要進行改進。通過將Jerk輸入估計引入CS模型中,結合改進的輸入估計方法[26],對CS模型中狀態方程和機動加速度方差的調整方法進行改進,從而建立起ACS模型。將目標反射的復增益建模為一階馬爾可夫過程,構建擴展目標狀態向量,對應的增廣矩陣狀態方程可表示為

2.3 觀測模型

利用ESPRIT算法[29]對接收信號進行處理,可從中提取k時刻目標相對雷達的徑向距離、多普勒頻率、方位角以及反射增益等估計信息。對應的等效非線性量測模型為

其中,Rk,q,fk,q和θk,q分別代表目標q相對雷達工作中心的徑向距離、多普勒頻率和方位角;(x0,y0)代表雷達工作中心位置;λ表示雷達的工作波長;arctan 2(·)表示4 象限反正切函數;代表一個1×i維的向量,其中除第j個元素為1外其余元素均為0;?k,q~N(0,Gk,q)為零均值高斯白噪聲,其協方差矩陣可表示為

其中,矩陣對角線元素分別對應于式(10)中所有量測值的量測誤差協方差,且滿足[8]

其中,Bw為半功率波束寬度。由式(12)可知,觀測信息中所有參數的誤差協方差均與發射功率Pk,q成反比。因此,量測誤差協方差矩陣可以重新表示為

由式(13)可知,通過提高發射功率,能夠提高雷達量測精度。當存在多批目標時,需要對有限的發射功率進行合理分配。特別地,當雷達系統在作戰環境中執行多目標跟蹤任務時,其功率分配策略的制定應充分考慮各目標威脅程度的差異。因此,需要建立合理的目標威脅評估機制來指導優化模型建立。

3 跟蹤性能量化

3.1 PC-CRLB推導

標準PCRLB[9]可為無偏估計提供一個下界,并且已經被證明在高信噪比情況下非常接近實際跟蹤誤差[14]。然而,由于標準PCRLB基于以往所有時刻的觀測信息對目標狀態進行估計,因此并不適用于觀測信息變化劇烈的機動目標跟蹤場景[13]。為了更好地利用觀測信息,并結合機動跟蹤的非線性特點,采用粒子濾波對目標狀態進行估計,并且推導基于粒子濾波的PC-CRLB。首先,PC-CRLB滿足如下不等式[10]:

其中,JP(xk,q|z1:k-1,q)和JZ(xk,q|z1:k-1,q)分別為先驗信息的PC-BIM矩陣和數據信息的PC-BIM矩陣。為避免式(16)中復雜的數學期望運算,結合粒子濾波機制,PC-BIM可近似計算為[31]

以往跟蹤性能驅動的資源分配文獻可具體分為:最大化最差目標的跟蹤精度[10]和最大化目標整體跟蹤精度兩種類型[12]。然而在實際應用,各目標的威脅度并不相同,因此簡單地對各目標的PC-CRLB進行求和或尋找其最大值,并將其作為目標函數不夠合理。3.2節將針對目標運動特性和電磁特性,構建威脅度評估模型。

3.2 目標威脅評估模型

隨著航空兵器的發展,空襲樣式逐漸多樣、目標運動特征變化逐漸加劇,導致地面防空雷達面臨的戰場環境復雜多變[34]。構建準確及時的目標威脅評估模型,可有效提升地面防空雷達系統的資源分配效能,增強抗機動突防能力。

針對目標威脅評估,主要考慮以下4個方面:徑向距離、徑向速度、目標航向、電磁散射特性。具體而言,目標威脅度的量化過程如下:

(1) 徑向距離

(2) 徑向速度

(3) 目標航向

(4) 電磁散射特性

至此,目標q在k時刻的威脅度可建模為

其中,ηj(j=1,2,3,4)為徑向距離、徑向速度、目標航向和隱身性能的權重系數,滿足η1+η2+η3+η4=1。對χk,q進行歸一化后,可得

4 低截獲背景下的功率分配算法

4.1 截獲概率模型

假設目標攜帶攔截接收器,當接收器檢測到雷達信號的概率大于其閾值時,可能對地面雷達生存造成威脅,因此需要減小目標的截獲概率。為簡化模型,假設所有目標攜帶攔截接收器相同且工作在同種模式。由文獻[23]可知,目標q在k時刻對雷達發射信號的截獲概率為

其中,pfa為攔截接收器的虛警概率;Gt為雷達發射天線增益;GI為攔截接收器的接收天線增益;GIP為攔截接收器的信號處理增益;k0代表玻爾茲曼常數;T0代表雷達接收機的噪聲溫度;BI為攔截接收器的帶寬;FI為攔截接收器的噪聲因子;erfc(x)函數展開為

4.2 功率優化模型

由上述分析可知,MIMO雷達對目標的功率分配與跟蹤精度和截獲概率均密切相關。一般而言,為了保證低截獲性能,要求系統降低發射功率,而為了實現高精度跟蹤,應提高發射功率。在防空作戰時,精確跟蹤能力和低截獲能力都是雷達系統的追求目標;并且,隨著戰場態勢的實時變化,對跟蹤能力和低截獲能力的要求也發生改變。為建立關于目標跟蹤精度和雷達低截獲概率的優化模型,我們引入了線性約束機制[36]。考慮到雷達發射功率越高,跟蹤誤差越低[15]、截獲概率越高[37],因此,為平衡二者的單調性差別,令

其中,F2(Pk,q)表示雷達未被目標q截獲的概率。

綜上,優化模型可描述如下:

4.3 求解方法

式(28)的約束條件均是線性的,因此其性質由目標函數決定。在式(28)中,目標威脅值主要由目標運動狀態和電磁特性決定,在固定時刻可以近似看作常數。類似地,任務重要性權值和在固定時刻也可以看作常數。而F1(Pk,q)和 F2(Pk,q)均是關Pk,q的凸函數[38],因此,式(28)是一個關于發射功率的凸優化問題,通過采用內點法[17]便可輕松求解,但其算法復雜度為O(Q3.5)[39]。為提高時效性,基于目標函數關于功率的單調性特征和minmax問題[40]的特點,本節給出一種算法復雜度僅為O(Q)的快速求解方法。算法描述如下:

引理1:Dk,opt≥Dk。

引理2:若Dk≤,?q=1,2,...,Q。則Dk,opt=Dk,此時與之對應的解Pk,q=是最優的。

引理3:若Dk>,?q=1,2,...,Q。則最優解Pk,opt中的第q個元素的值為Pk,q,opt=。

引理1—引理3的證明過程見附錄。結合上述3個引理,式(28)可以通過反復求解式(29)并且直接令Pk,q,opt=得到,流程見表1。

表1 功率快速求解算法Tab.1 Fast power solving algorithm

5 仿真結果及分析

5.1 基本參數設置

本節仿真設置如下:考慮單部集中式MIMO雷達同時跟蹤Q=3個目標的場景。假設在x-y平面中,雷達工作中心位于原點處。采用粒子濾波算法對目標狀態進行實時估計,粒子數量為Np=200。設置蒙特卡羅試驗次數為Nsim=100,共有30幀數據用于每次試驗。雷達和攔截接收器的相關參數設置情況見表2,各目標初始運動參數如表3所示,目標與雷達的幾何位置如圖4所示。假設用發射功率為 0.4Ptotal的雷達波束對距離50 km外RCS為1的目標進行照射時,相應的接收端觀測誤差為G0=diag(1002,102,0.12,22,22)。在ACS模型中,設置各目標的初始機動頻率為u0=0.06,各目標在x方向和y方向的加速度變化如圖5所示。

表2 仿真參數設置Tab.2 Simulation parameter setting

表3 初始時刻目標運動參數Tab.3 Initial target motion parameters

圖4 雷達與目標的空間位置關系Fig.4 Spatial position relationship between radar and target

圖5 各機動目標的加速度變化情況Fig.5 Acceleration variation of each maneuvering target

5.2 任務場景構建

為探究參數變化對功率分配結果的影響以及驗證模型的魯棒性,設置了目標RCS起伏模型,具體如圖6所示。

圖6 目標RCS起伏模型Fig.6 Target RCS fluctuation model

根據仿真條件,運用層次分析法來確定目標威脅度模型中各目標特性對應的權重值,最終計算得到整個跟蹤任務期中各目標威脅值如圖7所示。此外,為驗證優化模型的正確性,考慮兩種任務重要性權值模型,具體如下:

圖7 目標威脅權重模型Fig.7 Task threat weight model

(1) 常數任務重要性權值模型?1。在該模型下,雷達對所有目標的跟蹤性能和低截獲性能重要性權值恒定不變,且滿足=0.8,=0.2,?q=1,2,3,?k=1,2,...,30。在這種情況下,資源分配的主要方向是提高跟蹤精度。

(2) 時變任務重要性權值模型?2。在此模型中,雷達對各目標所賦跟蹤性能和低截獲性能重要性權值各不相同,且各目標對應權值隨時間發生變化,的所有取值如圖8所示。因此,雷達對各目標的任務執行要求實時改變。

圖8 ?2模型中的任務重要性權值Fig.8 Task importance weight in model ?2

5.3 算法性能檢驗

圖9和圖10分別給出了?1模型和?2模型下的目標的真實軌跡及其估計軌跡。結合圖8可知,對目標1而言,在?2模型中,從第16 s開始跟蹤精度重要性權值降到小于0.4,而低截獲能力重要性權值上升到大于0.6,因此,在?2模型中對目標1的跟蹤效果將弱于在?1模型中的跟蹤效果,這一點在圖9和圖10的跟蹤軌跡中可以看出。此外,對于目標2而言,隨著徑向距離越來越遠,跟蹤誤差也逐漸增大。在?1模型中,由于更多關注目標跟蹤精度,雷達對目標2的跟蹤效果較為理想。但在?2模型中,盡管在第15 s之后關于目標2的跟蹤精度重要性權值上升到了0.6,但仍然比在?1模型中的跟蹤效果更差。由于在第15 s之后,目標3的跟蹤精度重要性為所有目標中最高,且達到了?1模型中的權值,因此必然導致更多的功率資源將分配給目標3,進而使其跟蹤效果比在?1模型中更好。

圖9 ?1模型下的目標跟蹤軌跡Fig.9 Target tracking trajectory in model ?1

圖10 ?2模型下的目標跟蹤軌跡Fig.10 Target tracking trajectory in model ?2

值得注意的是,在兩種模型下,雷達均能較好地完成對所有目標的跟蹤任務。為進一步描述目標跟蹤精度,定義跟蹤效果最差目標的均方根誤差(Root Mean Square Error,RMSE)為

圖11和圖12分別給出了在本文所提算法(任務模型1和任務模型2)、min-max PC-CRLB功率分配算法和功率平均分配3種資源分配策略下計算得到的跟蹤效果最差目標對應的PC-CRLB和RMSE。其中,min-max PC-CRLB功率優化分配算法表示的優化模型如下:

圖11 各算法關于最差情況的PC-CRLB性能對比Fig.11 PC-CRLB performance comparison of each algorithm on the worst case

圖12 各算法關于最差情況的RMSE性能對比Fig.12 RMSE performance comparison of each algorithm on the worst case

在min-max PC-CRLB功率分配算法中,雷達僅需要盡可能地提高最差目標的跟蹤精度,并采用內點法對式(32)進行求解。由圖11和圖12可知,功率平均分配策略在所有算法中表現最差,并且不能滿足跟蹤誤差閾值要求。min-max PC-CRLB功率分配算法由于僅需考慮跟蹤精度,從跟蹤性能角度而言,在所有算法中表現最佳。所提算法同時考慮了目標跟蹤性能和雷達低截獲能力,因此,在跟蹤性能方面表現稍差于min-max PC-CRLB功率分配算法。此外,由于任務模型1更加關注跟蹤性能,因此,所提算法在任務模型1下的跟蹤性能優于任務模型2下的跟蹤性能。

圖13展示了不同資源分配策略下關于雷達信號的最大截獲概率性能對比。由結果可知,在功率平均分配機制下,隨著目標1距離雷達越來越近,其對雷達信號的截獲概率逐漸增加,最終達到1。由于min-max PC-CRLB功率分配算法僅考慮跟蹤性能,未對抗截獲能力進行優化,min-max PC-CRLB功率分配策略下的最大截獲概率在3種功率優化分配算法中表現最差,且在初始時刻概率值接近1。而相較任務模型1而言,任務模型2更加關注雷達的抗截獲能力,因此,任務模型2下所提算法的最大截獲概率在所有算法中最低。

圖13 各算法關于最大截獲概率的性能對比Fig.13 Performance comparison of each algorithm for maximum intercept probability

圖14和圖15分別給出了在?1模型和?2模型下采用本文所提算法得到的功率分配結果。圖中不同網格的顏色代表歸一化的功率分配比率,定義為=Pk,q/Ptotal。為更好地分析功率分配結果,圖16給出了各目標與雷達的實時距離。

圖14 ?1模型中的雷達功率分配結果Fig.14 Results of radar power allocation in model ?1

圖15 ?2模型中的雷達功率分配結果Fig.15 Results of radar power allocation in model ?2

圖16 各目標相對雷達的徑向距離Fig.16 Radial distance of each target relative to radar

首先,對圖14的功率分配結果進行分析。在?1模型下,雷達主要考慮各目標的跟蹤性能。根據式(12)可知,距離越遠、RCS越小的目標對應的量測誤差越大,其跟蹤誤差也會相應增大。結合圖7、圖14、圖16可知,由于目標1始終距離雷達最近,理論上能夠實現較好的跟蹤精度,在min-max優化框架下相應的分配資源應較少[11]。但由于其目標威脅度較高,因此目標1獲得功率資源較多。但隨著徑向距離的拉大,目標1所獲功率資源呈現階梯式的下降特征。此外,由于目標2的目標威脅度大于目標3,且目標2的徑向距離逐漸增加,直至成為3個目標中距離雷達最遠目標,使得跟蹤精度變差,因此,目標2所獲得的功率分配資源較多,并且隨時間增加呈現一定的上升特征。盡管目標3距離雷達較遠,但由于其威脅程度最小,并且RCS值最大,使其實際跟蹤精度較好地滿足預期跟蹤精度要求,因此,目標3所獲功率資源最少。

由圖8、圖14和圖15可知,相較于在?1模型下的功率分配結果,由于目標1在?2模型中第15 s之后跟蹤性能重要性權值下降到低于0.4,導致目標1在15 s之后功率分配值的階梯式下降特征更加明顯。結合圖8可知,由于目標2在跟蹤周期的前半段跟蹤精度重要性權值為0.4,為所有目標中最低,因此,目標2在第15 s之前分配得到的功率值最低。由于在?2模型下第15 s后目標3的跟蹤精度重要性權值上升到0.8,因此,在第15 s后目標3成為各目標中跟蹤精度重要性權值最高的目標。相應地,如圖15所示目標3從16 s開始分配得到最多的功率資源。

為驗證所提求解算法的時效性,在一臺搭載i7-10750h處理器、16 GB雙通道內存的計算機上,采用MATLAB R2020b軟件對算法求解時間進行計算。經過100次蒙特卡羅試驗取平均值后得到算法平均計算時間如圖17。通過和內點法進行比較,可知所提算法具有更高的時效性,可使平均計算時間降低近50%。

圖17 算法平均計算時間Fig.17 Average calculation time of algorithm

6 結語

本文基于集中式MIMO雷達系統,提出了一種低截獲背景下針對多機動目標跟蹤的快速功率分配算法。該算法構建了目標綜合威脅度模型,并在此基礎上建立了關于機動目標跟蹤誤差和雷達低截獲性能的加權目標函數。通過在給定發射功率預算條件下動態調整各發射波束的功率大小,以實現提高目標跟蹤精度的同時保證雷達具備低截獲能力。為提高模型求解時效性,采用了一種基于序列松弛的凸優化算法進行求解。

仿真結果表明:(1)相比于功率平均分配策略和以提高跟蹤精度為目標的功率分配策略相比,所提算法能夠在提高目標跟蹤精度的基礎上保證雷達具備低截獲能力;(2)在任務模型為常數模型和時變模型兩種情況下,所提算法均能實現較好的跟蹤性能和低截獲性能。因此所提算法可通過調節任務重要性權值的方式,實現在一定可調節范圍內提高目標跟蹤精度和低截獲性能,從而具有良好的任務設計自由度和魯棒性;(3)本文采用的基于序列松弛的快速求解方法比內點法求解速度提高近50%。由于本文僅考慮發射功率分配情形,接下來將考慮對信號帶寬和發射波形等參數進行聯合優化分配,以進一步提高多目標跟蹤性能。

附錄

(1) 引理1證明:

根據非線性優化理論[41],由于Dk是式(32)的最優解,而式(29)又是式(28)的松弛形式。因此引理1顯然成立。

(2) 引理2證明:

(3) 引理3證明:

猜你喜歡
分配模型
一半模型
基于可行方向法的水下機器人推力分配
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
應答器THR和TFFR分配及SIL等級探討
遺產的分配
一種分配十分不均的財富
績效考核分配的實踐與思考
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 亚洲人成网站观看在线观看| 国内精品久久九九国产精品 | 欧美精品二区| 亚洲综合片| 亚洲av片在线免费观看| 亚洲成人在线免费观看| 粉嫩国产白浆在线观看| 国产人前露出系列视频| 亚洲日韩精品综合在线一区二区| 思思热精品在线8| 秋霞午夜国产精品成人片| 88av在线播放| 精品少妇人妻av无码久久| 欧美精品v| 精品少妇人妻av无码久久| 成人国产免费| 女人18一级毛片免费观看| 亚洲首页在线观看| 国产成人精品亚洲日本对白优播| 91精品视频播放| 久久这里只有精品国产99| 成人午夜天| 久草视频中文| 亚洲乱强伦| 亚洲无码免费黄色网址| 亚洲一级毛片在线观播放| 日韩在线1| 国产精品不卡永久免费| 色视频久久| 天堂av高清一区二区三区| 在线欧美一区| 67194成是人免费无码| 亚洲男人的天堂久久香蕉| 这里只有精品在线| 人人爽人人爽人人片| 首页亚洲国产丝袜长腿综合| 免费AV在线播放观看18禁强制| 日韩福利在线观看| 国产亚洲精久久久久久无码AV | 亚洲精品视频网| 香蕉网久久| 日韩精品一区二区三区大桥未久 | 国产成人免费| 欧美成人手机在线视频| 动漫精品啪啪一区二区三区| 国产福利一区视频| 亚洲婷婷丁香| 国产综合亚洲欧洲区精品无码| 免费无码一区二区| 天堂va亚洲va欧美va国产| 国产毛片高清一级国语| 亚洲第一色网站| 亚洲性日韩精品一区二区| 99国产精品一区二区| 亚洲黄色网站视频| 成人午夜网址| 亚洲无码A视频在线| 久青草国产高清在线视频| 欧美区国产区| 九月婷婷亚洲综合在线| 人妻无码中文字幕第一区| 制服丝袜亚洲| 欧美成人h精品网站| 国产丰满大乳无码免费播放 | 日韩福利在线观看| 国产成人精品一区二区三区| 特黄日韩免费一区二区三区| 九九热精品视频在线| 国产精品青青| 亚洲αv毛片| 成人福利在线视频| 亚洲第一成年网| 国产精品思思热在线| 欧美日韩精品在线播放| 伊人成人在线| 国产不卡一级毛片视频| 日韩在线影院| 国产精品成人第一区| 亚洲精品卡2卡3卡4卡5卡区| 国产主播喷水| 91免费国产在线观看尤物| 99视频在线免费看|