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

基于EMD能量占比的海面漂浮小目標特征檢測

2021-01-26 10:46:10時艷玲劉子鵬張學良顧為亮
系統工程與電子技術 2021年2期
關鍵詞:特征檢測

時艷玲, 劉子鵬, 張學良, 顧為亮

(南京郵電大學通信與信息工程學院, 江蘇 南京 210003)

0 引 言

有效探測海面小目標是雷達檢測領域的熱點問題之一。針對非線性非平穩的海雜波,有效的統計模型有很多,比如K分布、對數正態分布、韋伯爾分布、逆高斯分布等[1-8],但由于海雜波的非線性非平穩性,統計模型存在一定的誤差。

同時,除了建模這一有效的檢測方法外,眾多學者還提出了基于特征的海面目標檢測算法。海內外的研究者關注較多的特征研究理論主要包括混沌理論和分形理論[9-11]。在特征的提取研究中,混沌理論相關的特征是否有效依然存在爭議[12]。一般認為,在海雜波的短時序列中并不能使用單分量的混沌模型來提取特征,但海雜波可能是多個混沌模型復合的。另一方面,分形理論由于具有較好的可解釋性[13],在海雜波特征研究方面廣受歡迎。

另外,提取特征的方法還包括時頻分析法[14-15]。時頻分析法主要包括:快速傅里葉變換、短時傅里葉變換、分數階傅里葉變換、小波變換等。于曉涵等人提出了基于短時稀疏分數階變換和短時稀疏分數階模糊函數的雷達機動目標檢測和估計方法[16]。但是,傳統的傅里葉變換分析方法缺乏對信號某一特定頻率的研究,短時傅里葉變換與小波分析僅能對信號某段特定頻率進行研究,不能同時準確地表達時間和頻率上的信號變化,因此這幾類時頻處理方法難以處理非線性非平穩的海雜波信號。經驗模態分解(empirical mode decomposition, EMD)[17]克服了上述時頻處理方法的缺陷,是真正意義上的自適應時頻分析方法,尤為適合解決現實生活中遇到的非線性非平穩的信號[18-19]。

近些年的研究和發展,使得EMD方法取得了非常廣泛的應用。張建等人提出了一種基于EMD的低頻固有模態函數(intrinsic mode function, IMF)分量的能量比目標檢測方法[20],通過采用高階IMF分量與低階IMF分量的能量比作為特征來實現目標檢測,效果顯著。關鍵等人提出了一種基于固有模態能量熵的微弱目標檢測算法[21],利用能量熵特征描述各階IMF分量的變化,有效增強了雷達對海雜波中微弱目標的檢測能力。張林等人提出了一種基于分形特性改進的EMD目標檢測算法[22],并采用快速傅里葉變換進行去噪,實現了性能優良的檢測算法。

海雜波的頻率一般位于低頻位置,而小目標信雜比較低,速度較慢,多普勒頻率窄,且目標多普勒頻譜易被海雜波的頻譜掩蓋。EMD方法能將接收的回波數據從高頻分解到低頻,實現從高頻到低頻由粗到細的過程。通過EMD的方法,可以細化海雜波和目標的頻率并分解成不同的IMF分量,利用不同的IMF分量觀察海雜波和目標的差異。但是,大多數采用EMD的研究方法主要存在兩處缺陷:第一,對于IMF分量的選取沒有固定的選取標準,導致很多IMF分量提取的特征并不能區分目標和雜波。第二,大多數研究者提取的特征性能有限,無法更進一步地提升EMD帶來的優勢。針對這兩點,本文的創新點主要是:第一,針對IMF分量的選取,提出了一系列的篩選標準;第二,根據提出的篩選標準,利用篩選出的IMF分量,提出了一種基于能量占比特征的檢測方法。經過大量的實驗對比分析,發現提出的檢測方法具有優良的檢測性能。

本文首先采用EMD將接收回波數據分解為若干個IMF分量,實現對接收回波的頻率從高頻到低頻的分解;在采用EMD后,純海雜波的能量主要集中在前幾個IMF分量中,而目標的前幾個IMF分量的能量相對較少。同時,由于海雜波的非線性非平穩特性,回波中分解出的IMF分量并不都適合用來做檢測。為此,分別建立IMF分量與原始數據的線性相關性,求出各分量的相關系數,并利用平均均值與標準差之比作為篩選IMF分量的準則,實現對目標所在的IMF分量的自動篩選。最后,確定目標所在的IMF分量在原始信號中的能量占比情況,并將能量占比作為特征來實現海面漂浮目標的異常檢測。最終,通過實測數據進行實驗,驗證了所提的能量占比算法較對比算法具有更好的目標檢測性能。

1 EMD算法的特征分析與提取

1.1 EMD算法簡介

EMD算法可以根據信號本身的局部時間尺度特征來進行平穩化分解,是希爾伯特黃變換(Hilbert-Huang transform, HHT)的核心,具有良好的自適應性,適合處理非線性非平穩的復雜信號,如海雜波信號。EMD可以將信號分為多個IMF分量和殘留信號,每個IMF包含同等數目的極值點和零交叉數,連續兩個極值點之間定義了信號局部波動特征,反映了信號在不同尺度上的特性。每個IMF都關于時間軸局部對稱。

設雷達接收回波數據的幅度為y(n),n=1,2,…,N,經EMD后[23],可表示成:

(1)

式中,xi(n)是第i個IMF分量;R(n)為殘差分量;I為IMF分量的個數。

據此,接收回波信號y(n)被分解為具有不同頻譜特性的IMF分量xi(n),低階IMF分量對應信號的高頻部分,高階IMF分量對應信號的低頻部分。

1.2 能量占比特征提取

相關系數是表征兩序列之間的線性相關程度的量,常作為統計特征值來進行數據線性相關程度的判斷,可應用于故障診斷、目標探測等。假設已經接收了K段長度為N在第p個距離單元的回波數據y(n,k,p),n=1,2,…,N,k=1,2,…,K,p=1,2,…,P,則第i個IMF分量xi(n,k,p)與接收回波y(n,k,p)的相關系數Corrxi,y(k,p)可表示為

Corrxi,y(k,p)=

(2)

(3)

(4)

由于IMF分量反映了數據在不同頻率范圍的分布情況,為了獲得目標所在的IMF分量,需要對IMF分量進行篩選,首先定義了均值-標準差之比(mean-standard deviation ratio, MSR),即

(5)

然后,當對所有P個距離單元的MSRi(p)進行平均時,定義平均MSR(mean of MSR, MMSR)作為篩選準則,即

(6)

式中,ε為一個較大的正數,ε的取值將依據具體的數據確定??紤]到適合的IMF分量不僅需要特征值上的差異明顯,且在整體上的分布要比較均衡,以減少海雜波和目標特征值上的混疊,因此希望海雜波和目標的標準差盡可能小。

通過式(6),得到了第p個距離單元的兩個IMF分量,假設為第j個和第l個IMF分量xj(n,k)和xl(n,k),j,l∈{1,2,…,I}。為了避免公式復雜,在后續的公式中,默認數據的第3維是指第p個距離單元。目標的存在不僅改變了原始海雜波的IMF分量的頻率結構信息,而且也改變了原始海雜波IMF分量的能量結構。通過大量的實驗分析,雜波的頻譜寬而目標的頻譜窄,且都在低階的IMF分量上,目標的頻譜被淹沒在雜波的頻譜上。同時,目標的IMF分量的能量占比特征比同階數雜波的IMF分量的能量占比要低。那么,xj(n,k)在接收回波數據中的能量占比(energy ratio, ER)為

(7)

考慮到目標的運動性,目標會影響多個IMF分量,由上述篩選過程可以求出其中受影響最嚴重的兩個IMF分量,經過大量實驗,這兩個IMF分量的能量占比特征均有較好的檢測性能。于是,采用均值能量占比(average energy ratio, AER)作為檢測統計量來實現海面漂浮小目標檢測,AER檢測統計量的計算公式為

(8)

式中,j和l為被篩選出的IMF分量的階數。利用AER作為檢測統計量,如果待檢測單元的檢測統計量小于判決門限則認為有目標;否則認為當前單元為純雜波單元。由于式(8)中的檢測統計量涉及許多復雜操作,實驗中采用蒙特卡羅方法來確定給定虛警下的檢測門限。

依據上面的特征提取過程,本文提出的基于EMD的AER特征檢測的流程圖,如圖1所示。圖1中,pf表示給定的虛警概率,Hi表示有目標,H0表示無目標。

圖1 基于EMD的能量占比特征檢測流程圖Fig.1 Flow chart of energy proportion feature detection based on EMD

依據圖1,本文提出的基于EMD的AER特征檢測算法的具體步驟如下。

步驟 1接收K段長度為N在第p個距離單元的回波數據的幅度為y(n,k,p)。

步驟 2將第k段接收回波數據進行EMD,計算I個IMF分量xi(n,k,p)與接收回波幅度y(n,k,p)的相關系數Corrxi,y(k,p)。

步驟 3根據式(6),篩選出第j和第l個IMF分量xj(n,k)和xl(n,k)。

步驟 4根據式(7)和式(9),計算篩選出的xj(n,k)和xl(n,k)在接收回波y(n,k)中的能量占比,并求AER(k)作為特征檢測量。

步驟 5采用蒙特卡羅方法來確定給定虛警下的檢測門限,并利用AER特征實現海面漂浮小目標的檢測。

2 實驗結果分析

本節首先給出相關系數分析,然后介紹ER特性分析,隨后給出利用MMSR篩選目標所在的IMF分量的合理性分析,以說明本文提出的AER檢測器的合理性,最后將AER與快速傅里葉變換(fast Fourier transform, FFT)、頻域香農熵(shannon entropy, SE)檢測器[24]以及低頻IMF能量比檢測器(proportion of low frequency IMFs energy, LIMF)[18]進行對比實驗。

采用全相參X波段(ice multiparameter imaging X-band, IPIX)雷達數據庫駐留模式下的10組實測雷達數據[25]以及南非科學與工業研究院(the Council for Scientific and Industrial Research, CSIR)雷達數據TFA10_001,其數據說明如表1所示。其中,前10組數據每組由14個相鄰距離單元組成,每個距離單元采樣點數為131 072(即131.072 s),目標是用金屬絲網包裹、直徑約1 m的聚苯乙烯泡沫塑料球。目標所在的距離單元稱為目標單元,目標單元周圍的2個或3個單元受漂浮目標的影響,稱為受影響單元,其他距離單元稱為雜波單元。這10組數據的距離分辨率為30 m,其中第一組數據為3~4級海況,屬于高海況,第7組數據為中海況,其他數據屬于低海況。第11組數據由64個相鄰距離單元組成,中心頻率為6.9 GHz,距離分辨率為15 m,海況等級是4級海況,屬于高海況。

表1 IPIX雷達數據和CSIR數據說明

2.1 相關系數分析

針對相關系數分析,采用第4組數據(#31)的HH極化以及第11組(TFA10_001)數據。圖2顯示了所有距離單元的前5個IMF分量與接收回波的平均相關系數。

圖2 前5個IMF分量的相關系數均值對比Fig.2 Comparison of correlation coefficient mean values of the first five IMF components

由表1可知,第4組數據(#31)的目標單元為第7個距離單元,受影響單元為第6、8、9個距離單元,剩余的單元為雜波單元。由圖2可知,首先目標單元和受影響單元的前幾個分量與接收回波數據相關系數比雜波單元的相關系數要小很多,說明采用相關系數可以粗略地分出目標單元和雜波單元。其次,由圖2(a)可知,在目標單元中,IMF1、IMF2、IMF3這3個分量的相關系數最小,說明目標嚴重影響了低階IMF分量。最后,隨著IMF分量的階數增加,目標單元的相關系數也增加,導致與海雜波的區分性小。在圖2(b)中,目標前3個IMF分量的相關系數小于雜波單元,而在IMF4和IMF5中這種規律并不存在。出現該現象的原因是海雜波主要集中在低頻部分,IMF分量是由高頻到低頻變化排列的,即低階分量代表高頻部分,而高階分量代表低頻部分。階數越高,海雜波的主導性越明顯,則相關系數會變大,導致海雜波與目標的可分性變差。而在低階分量上,海雜波的主導性弱,目標的主導性強,導致分解前后的相關系數小,海雜波與目標的可分性變好。所以,通過這個實驗可以得知,對IMF分量的選取不是隨機任意的,需要對IMF分量進行篩選,通過一定的篩選標準自適應地選取兩組IMF分量來進行后續的目標檢測,該兩組IMF分量含有表征海雜波和目標之間差異的特性。本文給出的其余組數據也顯示出相類似的實驗結果,不再一一列出。

2.2 ER特征分析

ER特征的實驗參數與相關系數的實驗參數相同。圖3顯示了所有14個距離單元的前5個IMF分量在原始信號中的ER值。從圖3中可以看出,首先,目標單元的IMF1在原始信號中的ER值非常小;其次,隨著IMF分量的階數增加,目標所在單元的ER值也隨之提高;最后,海雜波所在單元的各階IMF分量(除IMF1分量外)的ER值分布較均勻,且雜波所在單元前幾個IMF的ER值要高于目標單元的IMF的ER值。出現上述3個現象的主要是:首先,在目標單元,接收回波的高頻部分受到了目標小球的影響,由于目標的頻譜比較窄,使得能量比較集中在后幾個IMF分量中,因此前幾個IMF分量的ER比較小。其次,隨著IMF分量階數的增加,分解到高階IMF分量的能量增加,使得ER值增加。最后,雜波所在單元前幾個IMF的ER值要高于目標單元的IMF的ER值,主要是因為雜波的頻譜相對于目標的頻譜較寬,使得分解的能量較為分散地落到各階IMF分量上,與此同時,由于海尖峰的存在,也會有較多的能量分解到較低階的IMF分量上,使得雜波的較低階IMF分量的ER值高于目標的ER值。隨著IMF階數的增加,這種差距越來越小,甚至有些雜波的高階IMF分量的ER值會小于目標的ER值。鑒于此,可選用ER值作為特征區分海雜波和目標。

圖3 #31數據HH極化的前5個IMF分量能量占比的均值對比Fig.3 Comparison of the mean energy proportion of the first five IMF components for data #31 under HH polarization

經過式(6)對IMF分量的篩選,確定出IMF2和IMF3能體現目標與海雜波的區別。圖4顯示了IMF2和IMF3能量占比的檢測性能,縱軸表示檢測概率pd。利用蒙特卡羅方法來確定給定虛警范圍內的檢測門限,待檢測單元的檢測統計量小于判決門限則認為有目標;否則認為當前單元為純雜波單元。

圖4 #31數據HH極化下IMF2和IMF3能量占比檢測性能對比Fig.4 Detection performance comparison of IMF2 and IMF3 energy proportion for data #31 under HH polarization

兩個IMF分量已顯示出較好的檢測性能??紤]到目標的運動性,目標可能被分到兩個IMF分量中,所以單一選擇一個IMF分量是不合理的。于是,給出平均能量占比的方案,即AER,更加符合實際情況。不選擇IMF1的原因是,在圖2(a)和圖3中,海雜波所在單元IMF1與目標所在單元IMF1所體現出的相關性相當,ER值也不相上下,不能體現出海雜波與目標的明顯區別,故在實際應用中不宜選擇IMF1。

2.3 篩選合理性分析

本節主要對利用MMSR來選擇IMF分量的合理性進行分析。采用第1組數據(#17)和第4組數據(#31)進行對比分析,其中第1組數據是3~4級海況,屬于高海況,第4組數據是2~3級海況,屬于低海況。通過對高低海況的對比分析,來說明選擇IMF分量的合理性。圖5顯示了兩組數據HH極化下的MSR和各IMF分量的MMSR。

圖5 #17數據在HH極化下的MSR和MMSR對比Fig.5 Comparison of MSR and MMSR for data #17 under HH polarization

由表1可知,第1組數據(#17)的目標單元為第9個距離單元,受影響單元為第8、10和11個距離單元,剩余的單元為雜波單元。第4組數據(#31)的目標單元為第7個距離單元,受影響單元為第6、8和9個距離單元,剩余的單元為雜波單元。從圖5(a)和圖6(a)可以看出,目標所在單元的MSR較低,而雜波單元的MSR較高。這是因為在同一標準差水平下,目標的相關系數均值較低。圖5(b)和圖6(b)是對各IMF分量的MSR在所有距離單元上取平均,即MMSR,橫坐標表示分量的階數。從圖5(b)和圖6(b)可以看出,不管是高海況還是低海況,IMF2和IMF3均體現出較大的MMSR,這說明綜合使用IMF2和IMF3是合理的,而單一選擇某個IMF分量會導致無法達到最優的檢測性能。

圖6 #31數據在HH極化下的MSR和MMSR對比Fig.6 Comparison of MSR and MMSR for data #31 under HH polarization

表2給出了兩組數據HH極化下單個IMF分量能量占比的檢測概率。由表2可以看出,在各種虛警概率下,#17數據的IMF2和IMF3分量擁有較高的檢測概率,#31數據的IMF2和IMF3分量擁有相對較高的檢測概率。這主要是因為在不同海況下,采用EMD后,同一數據的不同IMF分量描述不同的頻率,而不同數據的同階IMF分量描述的頻率也不同。又因為不同海況下,目標所在單元受影響的頻率段不同,主要集中在高頻部分。因此,需要根據不同的海況來找出目標回波中受影響最嚴重的頻率段,也就是找出能夠凸顯海雜波和目標回波之間差異的IMF分量,然后提取ER特征做檢測。

表2 第1組數據(#17)和第4組數據(#31)在HH極化下不同IMF分量的能量占比檢測概率

針對第11組數據,做了同樣的分析處理,IMF分量的篩選結果如圖7所示。

圖7 第11組數據的篩選結果Fig.7 Filtering results of the eleventh set of data

由圖7可知,通過MMSR準則可以篩選出兩組分量分別是IMF1和IMF2,這說明提出的篩選準則能夠根據數據來做出合理選擇。

本文提出的MMSR準則在平均尺度和離散程度之間衡量,是針對具有不同均值的各距離單元在離散程度上的描述,即希望目標數據和海雜波數據在各自的均值水平上都比較聚集,防止目標和海雜波數據過多的混疊,從而能夠得出平穩且合理的檢測效果。通過上述實驗證明,所提出的MMSR準則能夠非常準確地自動篩選出最適合做檢測的IMF分量,也進一步驗證了綜合使用兩個篩選出的IMF分量是合理可行的。

2.4 檢測性能分析

本節采用表1給出的11組數據的4種極化方式,對比本文提出的AER與對比算法FFT、SE和LIMF 4種檢測器的性能。實驗參數為:K=10 000,I=8,N分別取512,1 024,2 048和4 096。虛警概率pf=0.001。其中,為了保證獲得足夠的數據,對相鄰數據段采取每次滑動10個點的處理。每種極化方式的觀測時長均為0.512 s,1.024 s,2.048 s和4.096 s。

2.4.1 AER檢測性能分析

圖8顯示了11組數據在4種極化方式下的檢測概率圖。為了便于分析,圖9顯示了11組數據的平均信雜比(average signal to clutter ratio, ASCR)圖。如圖8所示,在不同極化方式下有不同的檢測概率。其中HV和VH極化下的檢測概率相似,而HH和VV極化方式的檢測概率差異明顯。造成這種現象的主要原因是ASCR的不同,一般來說,ASCR越高,檢測概率越高,大部分的實驗結果與此結論是相符的。對于第11組數據,ASCR很高,所以具有非常好的檢測概率。在IPIX的10組數據中,除第1組數據(#17)和第7組數據(#280)外,其余8組均在低海況下采集。在整體上,由于提出的AER實質上反映的就是ASCR,所以在實測數據中,ASCR越高,檢測概率就越高。但這并不是嚴格成立的,檢測概率還要受到風速、浪高和雷達方向角的影響。對于第1組數據,HH極化下的ASCR最高,卻沒有達到非常高的檢測概率。同樣,HV和VH極化下的檢測概率也沒有達到與ASCR相應的水平,這主要是因為第1組數據是在高海況下采集的,掠射角度較小,浪高2.2 m,產生了遮擋效應,目標經常會被海浪掩蓋,從而影響了檢測性能。

圖8 11組數據在4種極化方式下AER特征檢測器檢測性能Fig.11 Detection performance diagram of AER feature detector with four polarization modes of 11 sets of data

圖9 11組數據的ASCRFig.9 ASCR of 11 sets of data

表3是11組數據在4種極化方式的ASCR與檢測概率的關系表,其中,N=4 096,虛警概率pf=0.001。ASCR的計算公式如下:

(9)

表3 檢測概率關于ASCR分布表

2.4.2 4種檢測算法的檢測性能對比

將分析本文提出的EMD下的AER特征檢測器與FFT、頻域SE檢測器以及LIMF檢測器的實驗性能。4種檢測算法針對4種極化方式的數據觀測時長為0.512 s,1.024 s,2.048 s和4.096 s,即采用同樣的方法對數據進行分段和選取。圖10~圖13顯示了4種檢測算法的檢測性能。從圖10~圖13中可以直觀地看出,在4種極化方式下,本文所提出的AER方法均有較好的檢測概率。FFT的檢測方法不能有效地提取特征,特別是在ASCR很低的數據上完全喪失性能。頻域SE檢測算法在采用FFT下,不能完全挖掘出海雜波的非平穩和非均勻特性。LIMF檢測算法同樣用EMD后的IMF分量提取特征,但在提取的特征上直接將前3組IMF分量混合在一起,忽略了分解時各IMF分量的差異性,雖然在某些特定數據上表現效果較好,但整體上遜色于提出的AER。本文所提方法挖掘出海雜波的非平穩和非均勻特性,又克服了FFT的缺點,所以在各組數據上都有較好的表現。這同時也說明了較之于傳統的傅里葉變換下的特征檢測,采用EMD以及所選取的AER特征在海雜波目標檢測上有明顯的優勢。圖14為第11組數據采用4種不同算法下的檢測概率圖,觀測時長分別為0.064 s,0.128 s,0.256 s和0.512 s。其中,需要指出的是,在觀測時長超過1.024 s后,4種檢測方法的檢測概率均已達到1,因而此處未將觀測時長超過1.024 s的對比放入文中。由圖14可以看出,本文采用的方法與其他方法相比略有優勢,不過微乎其微,但是相當高的檢測概率也足以說明本文采用的算法是行之有效的。

圖10 N=512時不同極化方式下4種極化方法的檢測結果對比(pf=0.001)Fig.10 Comparison of detection results of four methods under different polarization modes when N=512 (pf=0.001)

圖11 N=1 024時不同極化方式下4種極化方法的檢測結果對比(pf=0.001)Fig.11 Comparison of detection results of four methods under different polarization modes when N=1 024 (pf=0.001)

圖12 N=2 048時不同極化方式下4種極化方法的檢測結果對比(pf=0.001)Fig.12 Comparison of detection results of four methods under different polarization modes when N=2 048 (pf=0.001)

圖13 N=4 096時不同極化方式下4種極化方法的檢測結果對比(pf=0.001)Fig.13 Comparison of detection results of four methods under different polarization modes when N=4 096 (pf=0.001)

圖14 第11組數據在不同觀測時長下4種方法的檢測結果對比(pf=0.001)Fig.14 Comparison of detection results of four methods under different observation duration of the data for the 11th set of data (pf=0.001)

3 結 論

本文主要對接收回波數據的幅度進行了EMD,并對得到的IMF分量進行了研究分析。首先,從分解后的IMF分量與原數據的相關性角度切入,分析了各層IMF分量與原數據的相關性,同時對比雜波數據和目標數據之間的相關性差異,采用MMSR準則篩選出能充分體現出目標與雜波差異的兩組IMF分量。然后,利用IMF分量在原數據中的ER值作為提取的特征,對篩選出的IMF分量進行ER特征的提取與分析,發現目標單元的ER值比較低,而海雜波單元的ER值比較高。最后,對選出的兩組IMF分量的ER值進行平均,得到AER作為特征進行檢測,并與傳統的FFT檢測、FFT下的頻域SE檢測以及LIMF檢測做了對比,實驗結果顯示本文所提出的方法具有比較顯著的效果,說明所提出的方法是可行的。

猜你喜歡
特征檢測
抓住特征巧觀察
“不等式”檢測題
“一元一次不等式”檢測題
“一元一次不等式組”檢測題
“幾何圖形”檢測題
“角”檢測題
新型冠狀病毒及其流行病學特征認識
如何表達“特征”
不忠誠的四個特征
當代陜西(2019年10期)2019-06-03 10:12:04
抓住特征巧觀察
主站蜘蛛池模板: 欧美19综合中文字幕| 久久精品视频亚洲| 婷婷亚洲综合五月天在线| 亚洲国产精品美女| 蜜桃臀无码内射一区二区三区 | 天天躁夜夜躁狠狠躁图片| 国产午夜人做人免费视频中文 | 国产精品午夜福利麻豆| 美女内射视频WWW网站午夜| 日韩欧美国产中文| 亚洲国内精品自在自线官| 精品国产香蕉伊思人在线| 精品欧美视频| 国产极品美女在线| 成人av手机在线观看| 丝袜高跟美脚国产1区| 午夜国产精品视频| 四虎免费视频网站| 国产午夜精品鲁丝片| 青青操国产| 国产91视频观看| 亚洲国产清纯| 最近最新中文字幕免费的一页| 无码一区二区三区视频在线播放| a天堂视频| 久久亚洲高清国产| 亚洲永久色| 青青国产视频| 亚洲bt欧美bt精品| 国产丝袜第一页| 精品自窥自偷在线看| 欧美亚洲欧美| 2021精品国产自在现线看| 拍国产真实乱人偷精品| 99人体免费视频| 亚洲一级毛片免费看| 8090成人午夜精品| 男女男精品视频| 日韩一二三区视频精品| 丝袜无码一区二区三区| 亚洲色偷偷偷鲁综合| 亚洲中文字幕97久久精品少妇| 精品91视频| 亚洲综合香蕉| 亚洲男人的天堂久久精品| 天天摸夜夜操| 国产精品亚洲va在线观看| 国产视频资源在线观看| 国产成人亚洲无吗淙合青草| 国产一级妓女av网站| 91精品伊人久久大香线蕉| 一级高清毛片免费a级高清毛片| 欧美激情视频一区| 色综合激情网| 尤物国产在线| 99精品一区二区免费视频| a欧美在线| 国产精品天干天干在线观看| 久草视频福利在线观看| 国产精品久久久久久久久久98| 国产成人凹凸视频在线| 国产靠逼视频| 国产精品成| 亚洲美女一级毛片| 国产一线在线| 伊人丁香五月天久久综合 | 亚洲AV无码一区二区三区牲色| 国产在线91在线电影| 手机永久AV在线播放| 国产免费久久精品99re不卡| 五月婷婷丁香综合| 无码aaa视频| 凹凸国产熟女精品视频| 欧美亚洲日韩中文| 在线欧美a| 在线一级毛片| 高清视频一区| 天天综合亚洲| 伊人久综合| 久久国产黑丝袜视频| 婷婷激情亚洲| 免费在线a视频|