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

一種基于自適應(yīng)濾波的海雜波背景下多目標(biāo)檢測方法

2022-07-07 01:58:21馬紅光郭金庫姜勤波劉志強(qiáng)
現(xiàn)代信息科技 2022年4期

馬紅光 郭金庫 姜勤波 劉志強(qiáng)

摘? 要:文章提出一種基于自適應(yīng)濾波的多目標(biāo)檢測方法。將雷達(dá)回波等分成回波矩陣Xi,計算Xi的協(xié)方差矩陣并進(jìn)行特征值分解;利用特征值矩陣D計算奇異譜,估計主分量個數(shù)Nev,以Nev>3作為門限判斷回波矩陣Xi是否包含目標(biāo);通過特征矢量矩陣V構(gòu)成的自適應(yīng)濾波器對Xi濾波,估算濾波后回波脈沖的Pareto模型參數(shù),生成Pareto隨機(jī)序列;采用K-L散度識別目標(biāo)回波,用峰值檢測法確定各個目標(biāo)位置。通過實(shí)測海雜波數(shù)據(jù)實(shí)驗,驗證了所提方法的有效性。

關(guān)鍵詞:海雜波;多目標(biāo)檢測;特征值分解;自適應(yīng)濾波;K-L散度

中圖分類號:TN957.51? ? ? ? 文獻(xiàn)標(biāo)識碼:A文章編號:2096-4706(2022)04-0072-05

A Multi-target Detection Method Based on Adaptive Filtering in

Sea Clutter Background

MA Hongguang, GUO Jinku, JIANG Qinbo, LIU Zhiqiang

(Xian Daheng Tiancheng IT Co. Ltd. Xian? 710026, China)

Abstract: A multi-target detection method is proposed based on adaptive filtering. The radar echo is firstly evenly partitioned into echo matrices Xi. The covariance matrix of Xi is calculated and then the eigenvalue decomposition is performed. The singular spectrum is calculated via the eigenvalue matrix D, and the number of principal components Nev is determined. The threshold Nev>3 judges if targets are contained in Xi. The adaptive filtering is applied to Xi with the eigenvector matrix V. The model parameters of Pareto distribution are estimated for each filtered echo. The random series of Pareto distribution are generated. The target echoesare identified via K-L Divergence. The positions of targets are determined by peak finding technique. Trials have been conducted on the measured sea clutter datasets and the effectiveness of the proposed method is validated.

Keywords: sea clutter; multi-target detection; eigenvalue decomposition; adaptive filtering; K-L divergence

0? 引? 言

海雜波是雷達(dá)波束照射海面時波束覆蓋區(qū)域產(chǎn)生的回波。海浪宏觀運(yùn)動及海面不規(guī)則的毛細(xì)管狀運(yùn)動,使得海雜波具有較強(qiáng)的“非線性、非平穩(wěn)和非高斯”等特性,致使傳統(tǒng)的恒虛警檢測器(CFAR)不能適用于海雜波背景下的多目標(biāo)檢測[1]。學(xué)者們對海雜波特性的研究由來已久,早期工作主要集中在對海雜波統(tǒng)計學(xué)模型的研究上,典型的海雜波統(tǒng)計學(xué)模型有Weibull[2]、log-normal和K分布[3,4]。近年來,對實(shí)測海雜波信號“尖峰”特性的研究表明,Pareto分布能更加準(zhǔn)確地描述這一特性[5]。衡量模型有效性的一個重要前提條件是:在雷達(dá)相干處理時間內(nèi)海雜波必須是平穩(wěn)的隨機(jī)過程,以確??梢杂煤愣ǖ哪P蛥?shù)實(shí)現(xiàn)CFAR目標(biāo)檢測。雷達(dá)探測海面上弱小目標(biāo)時,需要增加相干處理時間以增強(qiáng)目標(biāo)回波能量,與此同時海雜波的非平穩(wěn)特性也顯著增強(qiáng),以至于不再滿足上述前提條件。因此,基于統(tǒng)計學(xué)模型的目標(biāo)檢測方法僅適用于相干處理時間短的場合。

為解決上述問題,已有大量研究成果見諸報端:Wu提出了一種基于奇異譜分析的海面小目標(biāo)檢測方法[6];Chen提出了一種利用SVD-FRFT濾波的海雜波壓制方法[7];Li提出了一種在強(qiáng)海雜波中檢測動目標(biāo)狀態(tài)的自適應(yīng)探測方法[8];翟東奇等提出了基于非線性自適應(yīng)濾波器的海雜波抑制技術(shù)[9];Zhang提出了一種基于改進(jìn)的快速聚類分段的海雜波中弱小目標(biāo)檢測方法[10];Lang針對高分辨SAR圖像目標(biāo)識別問題,提出了一種基于像素聚類的艦船目標(biāo)檢測方法[11];Su提出了基于混沌理論、徑向基神經(jīng)網(wǎng)絡(luò)及k-means聚類的海雜波預(yù)測方法[12];針對海雜波的非平穩(wěn)混沌特性,Ma提出了一種非平穩(wěn)混沌時間序列相空間重構(gòu)方法[13],在此基礎(chǔ)上提出一種基于多尺度有向Lyapunov指數(shù)的海雜波中弱小目標(biāo)檢測方法[14]。

上述方法的共同特點(diǎn)是不再依賴海雜波的統(tǒng)計學(xué)模型,有效克服了傳統(tǒng)方法的不足,但存在計算復(fù)雜度較高、部分方法需要運(yùn)用大量人為標(biāo)定的訓(xùn)練數(shù)據(jù)集,這些問題仍是其工程應(yīng)用的桎梏。為此,本文提出一種基于自適應(yīng)濾波的海雜波背景下目標(biāo)檢測方法。9C24E143-2BAE-460D-8C58-FABBD0E48C9A

1? 方法與步驟

本文所提的基于自適應(yīng)濾波的海雜波背景下目標(biāo)檢測方法,具有算法簡單快捷、不依賴任何目標(biāo)先驗知識,可自動區(qū)分海雜波和目標(biāo)回波等特點(diǎn)。所提方法的流程圖如圖1所示。

所提方法的計算步驟為:

Step 1:以雷達(dá)掃描整個監(jiān)測海面的回波為對象,根據(jù)對海觀測雷達(dá)天線轉(zhuǎn)動速度ΩE、雷達(dá)的波束方位寬度θA和脈沖重復(fù)頻率PRF,計算每個方位向回波脈沖數(shù)量:

(1)

Step 2:按照np將雷達(dá)掃描整個監(jiān)測海面的回波脈沖N等分為Nbin個方位向單元回波,構(gòu)成對應(yīng)各個方位向的回波矩陣Xi,Xi為np×nt的矩陣(i=1,2,…Nbin),Nbin=floor(N/np),floor(·)為取整函數(shù),nt為回波脈沖采樣點(diǎn)數(shù)。

Step 3:計算回波矩陣Xi的協(xié)方差矩陣:

(2)

其中,H表示回波矩陣的共軛轉(zhuǎn)置,C為np×np的正定Hermitian矩陣。

Step 4:對協(xié)方差矩陣進(jìn)行特征值分解:

[V? D]=eig(C)(3)

其中,D為對角矩陣,其對角線元為C的特征值,V各列為各個特征值對應(yīng)的特征矢量。

Step 5:提取D的對角線元,按降序重排C的特征值,并相應(yīng)調(diào)整V各列的位置,計算奇異譜:

(4)

其中,dj為重新排序后C的特征值。

Step 6:對回波矩陣進(jìn)行自適應(yīng)濾波:

Xi=VTXi? (5)

其中,T為矩陣轉(zhuǎn)置。

選擇一個門限Thr,將滿足σj≥Thr的特征值個數(shù)Nev作為主分量(Principal component, PC),將其余特征值作為次分量(Minor component, MC),將濾波后的回波Xi分為Sp和Sm。

Step 7:若Nev≤3,判定回波Xi為海雜波,反之,對Sp和Sm做進(jìn)一步處理(Sm可用于弱小目標(biāo)檢測)。

Step 8:采用最大似然估計依次估算Sp(或Sm)各行的Pareto分布模型參數(shù)(a,b):

(6)

其中,x為回波的瞬時幅度,nt為回波脈沖采樣點(diǎn)數(shù),Γ(·)為伽馬分布函數(shù),a為形狀參數(shù),b為尺度參數(shù)。

利用估計的模型參數(shù)(a,b),通過Pareto隨機(jī)數(shù)發(fā)生器生成與回波信號等長的隨機(jī)序列y,采用Kullback-Leibler(K-L)散度識別目標(biāo)回波脈沖:

(7)

其中,Px為回波信號幅度的概率密度函數(shù),Qy為隨機(jī)序列y的Pareto分布函數(shù)。

將K-L散度最大值對應(yīng)的回波信號識別為目標(biāo)回波。

Step 9:采用Matlab的findpeaks函數(shù)確定目標(biāo)所在距離:

(8)

其中,x=p或m,Sx(i,:)中,i為回波矩陣Sp或Sm被識別為目標(biāo)回波所在的行號,輸入選項‘MinPeakProminence、minP為最小峰值顯著性;函數(shù)的輸出pks、locs、w、P分別為被確認(rèn)有效的譜峰高度、位置、寬度和顯著性。

最小峰值顯著性的估計方法是,首先用沒有最小峰值顯著性約束的findpeaks測量回波信號的全部峰值顯著性P,計算P的均值μ和標(biāo)準(zhǔn)離差σ,令:

minP=μ+nσ,n?[1,3]? (9)

可得自適應(yīng)minP,確保目標(biāo)檢測門限值的合理性。

Step 10:重復(fù)S2~S9,即可獲得整個監(jiān)測海面上目標(biāo)的位置。

在步驟S3回波矩陣Xi中,每一行對應(yīng)一個回波基帶信號的采樣值vi=xi+jyi,回波的瞬時相位θi=tan-1(yi/xi),所構(gòu)造的協(xié)方差矩陣C的對角線元為每個回波脈沖瞬時功率,非對角線元,保留了回波的瞬時幅度、相位特征與脈沖間的相干特性,經(jīng)步驟S4~S6后,所得到是根據(jù)海雜波與目標(biāo)回波在幅度和相位上的不同而進(jìn)行的信號正交分離的結(jié)果,與僅從回波幅度特征區(qū)分雜波與目標(biāo)的方法相比具有明顯優(yōu)勢。

在步驟S6~S7中,選擇門限Thr和確定區(qū)分海雜波與目標(biāo)回波的原則是所提方法的核心問題。對大量海雜波數(shù)據(jù)的研究表明,海雜波主要由Bragg、白冠(Whitecap)和突發(fā)(Burst)散射形成[4],其中,Bragg散射為入射電波在海面產(chǎn)生的諧振現(xiàn)象,在海表面毛細(xì)管狀波和海浪重力波這兩個不同尺度的波長上都會發(fā)生。在一定的海況下,海浪會按其運(yùn)動周期發(fā)生后向散射,來自多個海面波的回波就會進(jìn)行相干疊加形成雷達(dá)接收的回波,是海雜波的主要成分。白冠是由波浪上端變窄下端變寬后形成類似劈狀的結(jié)構(gòu)而產(chǎn)生的散射。但并非所有劈狀結(jié)構(gòu)都能形成白冠散射,只有當(dāng)劈狀結(jié)構(gòu)破碎后才形成白冠,因此白冠在海雜波中所占比例較低。突發(fā)散射是由即將破碎的波浪面產(chǎn)生的后向散射,波面在破碎前形狀呈陡峭的平面狀,由此產(chǎn)生較強(qiáng)的鏡面反射。突發(fā)散射在HH極化下具有很強(qiáng)的功率;但在VV極化下,海面動態(tài)結(jié)構(gòu)形成的多徑效應(yīng)使得突發(fā)散射變得極其微弱,有時候甚至可以忽略。圖2為煙臺養(yǎng)馬島實(shí)驗場實(shí)測海雜波的奇異譜[15]。

由圖2可知,Bragg散射占比80%以上,白冠散射約占10%~20%,而突發(fā)散射的占比在0.01%左右。因此,門限Thr=0.001是我們在實(shí)驗中選取主分量的門限。

主分量個數(shù)Nev≤3時,回波矩陣Xi中不包含目標(biāo),反之,回波矩陣Xi中包含目標(biāo)。圖3為煙臺養(yǎng)馬島海雜波數(shù)據(jù)20200722150408_798_scanning.mat的主分量數(shù)量Nev在雷達(dá)掃描方位向的變化情況,雷達(dá)主要參數(shù)及實(shí)驗場地詳見[16]。如圖3所示,T1脈沖為雷達(dá)工作于“模式2”發(fā)射的第1個單載頻脈沖(載頻fc=9.3 GHz,脈寬τ=40 ns),T2脈沖為雷達(dá)工作于“模式2”發(fā)射的第2個LFM脈沖(載頻fc=9.3 GHz,脈寬τ=3 μs,調(diào)頻斜率K=75),2個脈沖瞬時帶寬相同,均為Bw=25 MHz,脈沖重復(fù)頻率PRF=1.6 kHz,天線水平波束寬度θA=1.2°,天線轉(zhuǎn)速ΩE=24 r/min,因此,各方位向的脈沖數(shù)np=14。由圖3的實(shí)驗結(jié)果并對照[15]中的目標(biāo)場景,選擇Nev>3作為目標(biāo)檢測的門限是合理的。9C24E143-2BAE-460D-8C58-FABBD0E48C9A

在確定回波中包含目標(biāo)之后,從自適應(yīng)濾波后的回波中識別目標(biāo)回波是所提方法的創(chuàng)新點(diǎn)之一。傳統(tǒng)方法主要是在歐氏空間通過計算幅度之間的差異來區(qū)分海雜波與目標(biāo)回波,虛警和漏報概率較高。相較于傳統(tǒng)方法,采用信息幾何方法處理這一問題有明顯的優(yōu)勢,海雜波與目標(biāo)回波被視為黎曼空間中隨機(jī)流形上的點(diǎn)集[17],利用K-L散度可正確估算隨機(jī)流形上點(diǎn)之間的差異,具有穩(wěn)定可靠的識別率。

這里需要說明步驟S8中為何要分別估計每個濾波后回波的Pareto分布模型參數(shù)(a,b)的原因,圖4為濾波后回波Pareto分布模型參數(shù)(a,b)隨方位向變化的曲線,為準(zhǔn)確判斷回波屬性,需要分別估計其模型參數(shù),生成對應(yīng)的Pareto隨機(jī)序列,降低誤判概率。

圖4? Pareto分布模型參數(shù)(a,b)在方位向的變化圖

在步驟S9中,采用回波信號中最小峰值顯著性minP作為檢測目標(biāo)的門限,峰值顯著性P=pks/w,即峰值高度與寬度之比,這一門限綜合考查峰值高度和寬度,并根據(jù)回波的動態(tài)特性自適應(yīng)地調(diào)整門限,在確保對目標(biāo)可靠檢測的同時還可以有效降低“海尖峰”造成的虛警概率,使整個目標(biāo)檢測過程具有恒虛警檢測特性。

2? 實(shí)驗結(jié)果及有效性分析

為驗證所提方法的有效性,我們對2020年7月8日15時04分08秒在煙臺養(yǎng)馬島實(shí)驗場實(shí)測海雜波數(shù)據(jù)集798~808進(jìn)行了實(shí)驗研究,每個數(shù)據(jù)集的目標(biāo)檢測結(jié)果基本一致,這里僅給出20200722150408_798_scanning.mat的結(jié)果。

表1、表2僅給出回波矩陣主分量對應(yīng)的目標(biāo)檢測結(jié)果,其中方位角為回波矩陣對應(yīng)方位角的均值。表1為T1脈沖回波全部的目標(biāo)檢測結(jié)果,觀察濾波后的回波可知,T1回波存在較強(qiáng)的近程雜波,為此,將距離小于2 km的目標(biāo)判斷為近程雜波,如圖5所示。

表1? T1脈沖主分量回波目標(biāo)檢測結(jié)果

由于T2脈沖功率遠(yuǎn)大于T1脈沖功率,所檢測到的目標(biāo)數(shù)量遠(yuǎn)多于表1中的目標(biāo)數(shù)量,為便于比較,表2僅列出了與表1目標(biāo)相近的方位向目標(biāo)檢測結(jié)果。圖5為雷達(dá)波束方位角為5.65°時的PC、MC回波,經(jīng)計算K-L散度,圖5(a)的第1行被判定為目標(biāo)回波,目標(biāo)檢測結(jié)果如表1所示,除近程雜波外,僅探測到6.73 km處的目標(biāo);對比表2的目標(biāo)檢測結(jié)果可知,除6.365 km的目標(biāo)外,11.252~13.886 km這一范圍內(nèi)存在連續(xù)12個峰值點(diǎn),對照[15]給出的目標(biāo)場景,該方位向存在一座航道浮標(biāo)與浮標(biāo)相距一段距離的芝罘島,因為T1脈沖功率小,該島嶼為弱目標(biāo),不能在PC回波中被檢測到,但在MC回波中,該島嶼的回波清晰可見。

由上述實(shí)驗結(jié)果可知,本文所提方法能夠準(zhǔn)確區(qū)分海雜波與目標(biāo)回波,具有探測海面弱小目標(biāo)的能力。

3? 結(jié)? 論

海雜波背景下的目標(biāo)檢測、識別與跟蹤是雷達(dá)信號處理領(lǐng)域中的熱點(diǎn)問題,經(jīng)過長期的研究與探索,已取得豐碩成果。近十年來,隨著非線性科學(xué)的不斷進(jìn)步,尤其是人工智能技術(shù)的快速發(fā)展,出現(xiàn)了許多新穎的信號處理方法,但受限于雷達(dá)信號處理對實(shí)時性等方面的要求,這些方法多數(shù)仍處于理論探索階段,尚不能用于工程實(shí)踐。西安大成科技的研究團(tuán)隊長期致力于該項研究工作,在關(guān)注新理論、新方法的同時,更加注重解決工程中的實(shí)際問題,本文所提方法正是基于此目的而取得的階段性成果,我們對該成果的工程應(yīng)用充滿信心。

參考文獻(xiàn):

[1] 逯旺旺,楊勇,張斌,基于譜峭度特征識別的海雜波弱小目標(biāo)檢測方法 [J].現(xiàn)代信息科技,2020,4(4):31-35.

[2] FAY F A,Clarke J,Peters R S. Weibull distribution applied to sea-clutter [C]. London:Proc. IEE Conf. Radar77,1977.

[3] SAYAMAS,SHUJI,SekineM,et al. Log-normal, log-Weibull and K-distributed sea clutter [J].IEICE Transactions on communications,2002,85(7):1375-1381.

[4] WARD K D,TOUGH R J A,WATTS S. Sea Clutter: Scattering, the K Distribution and Radar Performance [M].London:The Institution of Engineering and Technology,2013.

[5] GRAHAM,VICTOR,WEINBERG. Noncoherent Radar Detection in Correlated Pareto Distributed Clutter [J].IEEE Transactions on Aerospace and Electronic Systems,53(5):2628-2636.

[6] WU X J,DING H,LIU N B,et al. A Method for Detecting Small Targets in Sea Surface based on Singular Spectrum Analysis [J].IEEE Transactions on Geoscience and Remote Sensing,2021,60:1-17.

[7] CHEN Z,HE C,ZHAO C,et al. Using SVD-FRFT Filtering to Suppress First-Order Sea Clutter in HFSWR [J].IEEE Geoscience and Remote Sensing Letters,2017,14(7):1076-1080.9C24E143-2BAE-460D-8C58-FABBD0E48C9A

[8] LI Y,SIRA S P,MORAN B,et al. Adaptive Sensing of Dynamic Target State in Heavy Sea Clutter [C]//2007 2nd IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing. St. Thomas:IEEE,2007:9-12.

[9] 翟東奇,江朝抒,鄧曉波,等.基于非線性自適應(yīng)濾波器的海雜波抑制技術(shù) [J].航空科學(xué)技術(shù),2018,29(6):73-78.

[10] ZHANGH H Y,ZHAO Z,XIAO F X. Robust Detection Method of Small Targets in Sea-Clutter via Improved Fast clustering Segmentation [C]//2016 8th International Conference on Intelligent Human-Machine Systems and Cybernetics (IHMSC). Hangzhou:2016:123-126.

[11] LANG H T,XI Y Y,ZHANG X. Ship Detection in High-Resolution SAR Images by clustering Spatially Enhanced Pixel Descriptor [J].IEEE Transactions on Geoscience and Remote Sensing,2019,57(8):5407-5423.

[12] SU X H,SUO J D. Prediction of Sea Clutter Based on Chaos Theory with RBF and K-mean Clustering [C]//2006 CIE International Conference on Radar. Shanghai:IEEE,2006:1-4.

[13] MA H G,ZHANG C L,LI F. State space reconstruction for nonstationary time-series [J].Journal of Computational and Nonlinear Dynamics,2017,12(3):031009.

[14] WANG R,LI X Y,MA H G. Detection of small target in sea clutter via multiscale directional Lyapunov exponents [J].Sensor Review,2019,39(6):752-762.

[15] 劉寧波,丁昊,黃勇,等.X波段雷達(dá)對海探測試驗與數(shù)據(jù)獲取年度進(jìn)展 [J].雷達(dá)學(xué)報,2021,10(1):173-182.

[16] 劉寧波,董云龍,王國慶,等.X波段雷達(dá)對海探測實(shí)驗與數(shù)據(jù)獲取 [J].雷達(dá)學(xué)報,2019,8(5):656-667.

[17] 趙興剛,王首勇.基于K-L散度和散度均值的改進(jìn)矩陣CFAR檢測器 [J].中國科學(xué):信息科學(xué),2017,47(2):247-259.

作者簡介:馬紅光(1959—),男,漢族,河南鄭州人,教授(總師),博士,主要研究方向:非線性信息處理、目標(biāo)探測與識別;郭金庫(1980—),男,漢族,山東菏澤人,副教授(CEO),博士,主要研究方向:復(fù)雜電磁環(huán)境下目標(biāo)檢測與識別;姜勤波(1976—),男,漢族,江蘇丹陽人,副教授,博士,主要研究方向:復(fù)雜系統(tǒng)建模與仿真、電子對抗;劉志強(qiáng)(1979—),男,漢族,四川崇州人,副教授,博士,主要研究方向:復(fù)雜系統(tǒng)建模與仿真、電子對抗。9C24E143-2BAE-460D-8C58-FABBD0E48C9A

主站蜘蛛池模板: 波多野结衣视频一区二区| 国产欧美日韩一区二区视频在线| 女同久久精品国产99国| 久久青草视频| 午夜a视频| 蝴蝶伊人久久中文娱乐网| 亚洲最新在线| 男人天堂亚洲天堂| 视频二区中文无码| 四虎永久免费在线| 亚洲日韩高清在线亚洲专区| 国产精品一区二区在线播放| 国产精品一线天| 亚洲国产成人麻豆精品| 国产丝袜啪啪| 久久中文电影| 国产视频a| 92精品国产自产在线观看| 永久毛片在线播| 亚洲国产日韩在线成人蜜芽| 欧美日韩亚洲国产主播第一区| 亚洲精品无码AⅤ片青青在线观看| 亚洲三级成人| 男女性午夜福利网站| 青青青视频蜜桃一区二区| 日本亚洲最大的色成网站www| 欧美一区二区三区国产精品| 精品人妻AV区| 国产日本欧美在线观看| 国产啪在线91| 久久精品波多野结衣| 欧美亚洲国产精品第一页| 亚洲视频四区| 一本二本三本不卡无码| 欧美激情综合| 日韩美一区二区| 国产色婷婷视频在线观看| 欧美特黄一级大黄录像| 欧美中文字幕在线二区| 国产在线观看第二页| 永久成人无码激情视频免费| 日韩欧美中文字幕在线韩免费| 波多野结衣的av一区二区三区| 91精品久久久久久无码人妻| 久久中文字幕不卡一二区| 666精品国产精品亚洲| 3D动漫精品啪啪一区二区下载| 日本国产精品| 美女视频黄频a免费高清不卡| 免费国产好深啊好涨好硬视频| 日韩第九页| 欧美有码在线观看| 国产区网址| 亚洲国产精品一区二区第一页免| 久久久久久尹人网香蕉| 国产成人91精品免费网址在线| 成人在线亚洲| 国产在线91在线电影| 亚洲无卡视频| 亚洲欧美极品| 狠狠做深爱婷婷综合一区| 婷婷丁香在线观看| 欧美日韩亚洲国产| 国产a v无码专区亚洲av| 色婷婷啪啪| 丰满少妇αⅴ无码区| 在线观看91香蕉国产免费| 国产免费久久精品99re丫丫一| 中文字幕乱妇无码AV在线| 免费不卡在线观看av| 亚洲国产中文欧美在线人成大黄瓜| 亚洲av无码人妻| 99热免费在线| 欧美国产日韩在线| 免费日韩在线视频| 中文字幕人成人乱码亚洲电影| 午夜性爽视频男人的天堂| 欧美一区二区福利视频| 在线a网站| 精品少妇人妻一区二区| 国产99免费视频| 亚洲欧美综合精品久久成人网|