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

基于最小一乘的虛擬陣元波束形成仿真研究

2020-05-09 08:57:58張偉民郭海濤金其余田原嫄焦圣喜
聲學技術 2020年2期
關鍵詞:信號方法

張偉民,郭海濤,金其余,田原嫄,焦圣喜

基于最小一乘的虛擬陣元波束形成仿真研究

張偉民1,郭海濤2,金其余3,田原嫄4,焦圣喜5

(1. 內蒙古大學電子信息工程學院,內蒙古呼和浩特 010021;2. 海南熱帶海洋學院海洋信息工程學院,海南三亞 572022;3. 內蒙古大學數學科學學院,內蒙古呼和浩特 010021;4. 東北電力大學機械工程學院,吉林吉林 132012;5. 東北電力大學自動化工程學院,吉林吉林 132012)

水下多波束成像系統成像質量的關鍵主要在于波束形成的質量。在信號頻率一定時,常規相移波束形成只有通過增加基陣的孔徑來得到更窄的波束,但是在有限的應用條件下,這種方法又受到實際工程的限制。因此,提出一種基于最小一乘估計的虛擬陣元波束形成方法。在低信噪比的情況下,利用最小一乘估計的穩健性可以得到更加準確的估計信號,從而在不增加基陣尺寸的情況下獲得更窄和抗干擾性更強的波束,提高陣增益和聲吶圖像的角度分辨率。通過Matlab仿真并與其他方法比較驗證,證明了上述方法的優越性。

多波束成像系統;最小一乘估計;穩健性;虛擬陣元;波束形成

0 引言

多波束成像聲吶通常在輻射角度區域形成大量的窄波束,其特點是高分辨率[1]、高數據率、寬覆蓋[2-4]且能夠實時成像[5]。而這些波束的性能直接影響著成像聲吶的方位分辨率[6-7]。常規波束形成只有通過提高信號的發射頻率和接收陣列的孔徑得到較窄的波束來提高圖像的分辨率。由于實際情況中,高頻信號衰減嚴重,陣列孔徑不能做得很大。因此,在保證陣列增益、波束寬度等指標優良的條件下,減少陣元使用數成為陣列發展的趨勢之一[8],虛擬陣元波束形成就是在這種情況下產生的。此類算法能夠使用較少的陣元數目得到較好的波束性能。人們在數理統計中用“穩健性”這個概念來描述異常值對某種方法的影響程度。虛擬陣元接收信號估計算法的準確性和穩健性決定了虛擬陣元波束形成的性能和應用的可靠性。

虛擬陣元波束形成是虛擬陣列波束形成方法中的一類,主要有線性預測法[9]、最小二乘法[9]和基于遺傳算法的聲矢量陣虛擬陣元波束形成法[10]等。線性預測法是將虛擬陣元的前個(或后個)陣元數據作為參考信號來實現虛擬陣元的構建。隨著虛擬陣元數目的增加,估計結果會產生明顯的累積誤差。最小二乘法通過處理真實陣元接收信號來估計虛擬陣元接收信號。最小二乘法相比線性預測法估計結果更加準確且不產生累積誤差。但實際海洋環境中的噪聲較為復雜,接收信號中必然會存在異常值,而最小二乘法的穩健性較差,使估計誤差增大。基于遺傳算法的聲矢量陣虛擬陣元波束形成是通過遺傳算法優化得到估計系數,從而構建出虛擬陣元接收信號。但在回波信號信噪比較低時,實陣元接收到的信號復雜且奇異值較多,通過遺傳算法優化得到的解容易陷入局部極值,使估計誤差增大。

本文利用最小一乘(Least Absolute Deviation, LAD)在估計回歸系數時穩健的優點[11-12],給出一種基于最小一乘穩健估計的虛擬陣元波束形成方法。由于最小一乘計算困難,所以利用文獻[13]線性規劃方法快速準確地估計出最小一乘線性回歸系數,從而解決了最小一乘求解難的問題。

1 最小一乘估計

1.1 線性預測模型

用向量及矩陣形式表示為

其中:為元素全為1的維列向量;

1.2 最小一乘線性回歸系數估計算法

最小一乘的線性回歸系數需要通過求解一個無約束不可微最優化問題來獲得,即

2 虛擬陣元波束形成的構建

2.1 信號模型與噪聲模型

對于單頻信號,基陣的指向性函數與陣元間隔有關,反映了信號在陣元之間存在一定的空間相關性[15]。這是虛擬陣元可以實現的理論依據。通過處理實陣元接收信號,估計構建出虛擬陣元接收信號,從而在理論上擴展了陣列的孔徑。因為在應用最小一乘進行最優系數的估計時,需要在回波角度附近進行觀測來獲得多組觀測方程,并根據這些方程求出最優解系數,所以當接收基陣形成單個波束時,構建虛擬陣元首先要對回波進行方位估計,這在很大程度上決定了信號估計的準確度。而多波束聲吶具有很多的波束形成器,使各個波束順序排列充滿整個搜索空間,并進行空間分割,相當于許多部單波束聲吶同時工作,每個波束形成器只負責指向一個角度,所以在進行虛擬陣元波束形成時無需進行方位估計。

現假設陣列為元均勻直線陣,處于目標回波聲場的遠場區,目標回波(,)可近似為平面波,回波信號頻率為0,背景噪聲為加性高斯限帶噪聲。不考慮各個陣元的指向性,將接收陣的第一個陣元作為基準陣元,回波與基陣垂直方向的夾角為0。如圖1所示,斜線表示回波信號,圖中黑色圓點為真實陣元,白色圓點為虛擬陣元。

在確定了信號頻率與陣元數時,由線陣的指向性函數可知,波束受到入射角度的影響,所以忽略時間因素。回波信號可表示為

式中:s(θ)是無噪聲的回波信號,即 ; ,為陣列相鄰陣元接收信號之間的時延差;為加性高斯限帶噪聲。

用矩陣形式可表示為

其中:

同理,可構建多個虛擬陣元的接收信號,即:

其中,值越小,說明估計的虛擬陣元信號越準確可靠。

2.2 陣增益

如圖1所示,0表示回波信號入射方向與基陣法線的夾角,各陣元靈敏度一致,信號與噪聲互不相關,且不同陣元接收的噪聲也是不相關的。兩個陣元信號之間的互相關系數為

分子為兩個陣元信號的時間平均,分母為歸一化因子。利用虛擬陣元接收信號估計算法,使陣元數為的實陣列陣元數虛擬擴展至+個。由于信號環境復雜,觀測數據不可避免地存在噪聲,在構建虛擬陣元信號時必然會產生估計誤差,誤差大小可用值反映。各個實陣元信號中的噪聲是非相關的,所以各虛擬陣元信號中的噪聲是由各實陣元噪聲的線性組合構成的,虛擬擴展后的陣列噪聲包括實陣列中的非相關噪聲和個虛擬陣元中的相關噪聲。根據波束形成原理,實陣元和虛擬陣元的接收信號線性疊加形成指向性,虛擬陣元信號中的相關噪聲和實陣元的非相關噪聲疊加,最終陣列噪聲仍為個不相關的噪聲。回波信號為式(12)中所示,信號的幅值為1,若噪聲為平穩信號,則基陣輸出端的平均功率為

式中前兩項為基陣輸出信號平均功率,最后一項為基陣輸出噪聲平均功率。則基陣的輸出信噪比為

單個無方向陣元的輸入信噪比為

陣增益定義為

由陣增益可知,當虛擬陣元接收信號估計準確時,式(25)得到最大值。當信號完全相關時,陣增益為

從式(26)可以看出,虛擬陣列輸出信噪比的大小主要由的值決定。值表示估計誤差,這說明估計算法的準確性與虛擬陣元波束形成的性能緊密相關。從式(25)陣增益表達式上看出,陣增益的大小不僅與估計算法的準確度有關,還與實陣元和虛擬陣元接收信號之間的相關性有關。相關性越強,分子前一部分的值越大,估計信號也越準確,的值也會隨之減小,陣增益得到提高。反之,則陣增益逐漸下降。這表明虛擬陣元接收信號的數目存在約束,不能無限制地擴展。

3 計算機仿真

下面通過計算機仿真來驗證算法的性能,并與其他相關算法,如常規相移法和文獻[10]中的線性預測法、最小二乘法進行比較。仿真中,陣列采用元等間隔均勻線陣,所有陣元各向同性,置于目標回波聲場的遠場區,陣元間隔=0.01 m,虛擬陣元波束形成方法皆前向構建個虛擬陣元接收信號。回波信號中心頻率為100 kHz,入射角度為0,背景為加性高斯限帶噪聲。

3.1 波束性能及陣增益分析

通過比較最小一乘(LAD)虛擬陣元波束形成與常規相移波束形成的波束圖和陣列增益隨著虛擬陣元信號數目變化曲線,來驗證算法的性能和有效性。圖2是在無噪聲和有噪聲時,16個實陣元采用最小一乘法構建8個虛擬陣元和常規相移法的波束比較圖。圖2(a)中顯示,最小一乘法在理想情況下得到的波束寬度相比相移法明顯變窄4o,兩種方法形成的旁瓣高度相差不大,第一旁瓣均在-13 dB左右,其他旁瓣高度基本保持在-18 dB附近。圖2(b)和圖2(c)是將兩種方法在回波信號信噪比為-4 dB和4 dB時的比較結果。從比較的波束圖中可以看出,在低信噪比的情況下,最小一乘法得到的波束主瓣更窄。由于最小一乘估計的穩健性,所以在信號環境較差的時候依然能估計出較為可靠的值,構建的虛擬陣元接收信號與實際陣元接收信號的誤差更小。從圖2(b)和圖2(c)中可以看出,最小一乘法相比常規相移法在不同信噪比時都能獲得較窄的波束,而旁瓣高度兩者相差不大。通過上述實驗可以看出最小一乘法在低信噪比的環境下仍可得到較窄的波束。

圖2 無噪聲和有噪聲時常規相移和LAD虛擬陣的波束形成圖

為了能直觀地體現出陣增益的變化規律,在陣增益的實驗中將陣元間距增大至0.1 m,信號頻率更改為1.5 kHz,實陣列增益曲線的陣元間距設為0.1 m,并與不同陣元間距的虛擬陣列增益的變化曲線進行對比。圖3是實陣列增益和虛擬陣列的陣增益隨著虛擬陣元信號數目變化的曲線圖。虛擬陣列分別設置了3種不同陣元間距:=0.1 m,=0.2 m,=0.3 m。仿真中回波信噪比為0 dB,回波角度0為30o,實陣陣元數為8個,信號脈沖寬度為10 ms。從圖3的曲線圖中可以看出。選定恰當的虛擬陣元數時,最小一乘法能夠得到比較高的陣增益。由式(25)和仿真曲線可以驗證,實陣列增益值在上述參數確定后為一固定值。由圖3中的陣增益曲線圖可以看出,通過構建虛擬陣元可以有效地增加陣增益,但隨著構建的虛擬陣元信號數目增加到某一臨界值時,陣增益曲線逐漸下降。如陣元間距為0.1 m的陣增益曲線,當虛擬陣元信號數目為45時陣增益達到最高值,隨后陣增益開始減少。觀察陣元間距和虛擬陣元數目臨界值的變化可以看出,增大陣元間距后,虛擬陣元數的臨界數減小,這反映了實陣元接收信號與虛擬陣元接收信號之間的相關性與陣元的間距有關,且呈反比關系。

圖3 陣增益隨著虛擬陣元信號數的變化曲線

3.2 不同方法對比分析

在信噪比為-4 dB、0 dB和4 dB的回波環境下,對比最小一乘(LAD)、最小二乘法和線性預測法3種虛擬陣元波束形成方法的波束圖,結果如圖4所示。圖4(a)~4(c)是三種方法隨著信噪比的增加的波束圖,隨著信噪比的增加,最小一乘法(LAD)和最小二乘法的波束主瓣寬度和旁瓣高度逐漸接近。這是因為最小一乘與最小二乘算法在理論上的估計準則導致的,當回波信噪比較高時,估計結果基本相似。在不同信噪比的環境下,線性預測法的旁瓣高度相比其他兩種算法總能保持最低。在低信噪比時,相比其他算法,最小一乘法總能得到較窄的波束,波束指向性和抗干擾能力更好,但旁瓣的高度不如其他兩種虛擬陣元波束形成方法的旁瓣低,結果如圖4(a)所示。

圖4 信噪比為4 dB、0 dB和4 dB時LAD法、線性預測法和最小二乘法三3種虛擬陣的波束形成圖

通過上述仿真可知,虛擬陣元接收信號估計算法的準確度和實現方式直接影響波束性能。若算法估計誤差較大或者在低信噪比的環境中下算法的穩健性較差,則陣列輸出波束的形狀也會產生明顯畸變,波束的性能將受到極大影響。

3.3 多波束性能分析

3.1和3.2節中的波束仿真是在單個波束的情況下,研究了最小一乘虛擬陣元波束形成算法的優缺點。圖5是常規相移法與最小一乘法進行多波束仿真的結果。為了研究最小一乘法角度分辨的性能,以波束主瓣在坐標軸上交點之間的夾角為間隔形成三個順序排列的波束,并與常規相移法進行比較,仿真中采用16個實陣元的線陣并構建8個虛擬陣元接收信號。從圖5(a)中可以看出,常規相移法的三個波束主瓣的極大值分別在5o、15o和25o方向上,角度分辨率為10o。圖5(b)最小一乘虛擬陣元波束形成方法畫出三個波束圖,三個波束主瓣的極大值分別出現在5o、12o和19o方向上,角度分辨率為7o。相比常規相移波束形成,最小一乘法的角度分辨率提高了3o。角度分辨率的提升將有利于對目標的識別和成像。如圖5(b)所示最小一乘虛擬陣元多波束形成的旁瓣相比常規相移多波束有所升高,主瓣角度附近的旁瓣高度相比常規相移波束形成提高了約3 dB。

圖5 常規相移多波束形成和LAD虛擬陣元多波束形成的波束圖

3.4 虛擬陣元數目對波束性能的影響

下面在信噪比為-4 dB、實陣元和觀測方程個數確定后,分析波束性能受虛擬陣元數目的影響情況。為了使現象明顯,本次仿真實際陣元數為8個,觀測方程以1o為間隔在0=0o兩側得到31組方程,其他仿真條件同上,仿真結果如圖6所示。

圖6 實陣元為8個時含0、8、32和128個虛擬陣元的線陣波束圖

從圖6中可以看出,增加虛擬陣元接收信號的數目可以明顯地提升波束的指向性。仿真中虛擬陣元接收信號個數由0增加到32時,形成的波束主瓣形狀逐漸尖銳且寬度逐漸變窄,旁瓣高度逐步降低。但是當虛擬陣元接收信號個數增加至128個時,主瓣兩側的旁瓣高度明顯增高,在進行目標方位的判斷時會造成模糊。旁瓣高度增加是因為利用最小一乘構建的虛擬陣元接收信號是完全根據真實陣元接收信號構建的,虛擬陣元數目的不斷增加,虛擬陣元接收信號與實際接收信號之間的相關性開始減弱,估計的誤差逐漸增大,導致波束的形狀發生畸變。以上結果也與圖3中曲線的變化相統一,當構建的虛擬陣元接收信號數目達到臨界值時,陣增益開始下降,波束的性能也開始變差。若虛擬信號數目大于臨界值時,波束性能將不能有效提升,因此虛擬陣元接收信號數目不能無限制地增加。

4 結論

對虛擬陣元接收信號的估計構建是虛擬陣元波束形成算法的關鍵,估計算法直接影響波束形成的性能。本文研究了現有虛擬陣元波束形成方法構建虛擬陣元接收信號的區別,提出了一種基于最小一乘穩健估計的波束形成方法,該方法在低信噪比的環境下構建的虛擬陣元接收信號更加準確、穩健。

相比其他虛擬陣元波束形成方法,本文方法形成的波束抗干擾能力更強且波束寬度更窄。該方法形成多波束時,可在不增加基陣孔徑的同時得到的較高的角度分辨率,這將有利于成像聲吶的目標識別。從理論分析和實驗仿真都可看出,虛擬陣列比實際陣列可以獲得更高的陣增益。在實際中,噪聲并非理想加性高斯限帶噪聲且信噪比往往較低,波束形成方法抗干擾的能力尤為重要,因此本文提出的最小一乘虛擬陣元波束形成在-4 dB 這樣較低信噪比的情況下更具優勢。

[1] XIA W J, JIN X, DOU F W. Thinned array design with minimum number of transducers for multi-beam imaging sonar[J]. IEEE Transactions on Oceanic Engineering, 2017, 42(4): 892-900.

[2] BROWN C J, BLONDEL P. Developments in the application of multibeam sonar backscatter for seafloor habitat mapping[J]. Applied Acoustics, 2009, 70(10): 1242-1247.

[3] FONSECA L, MAYER L. Remote estimation of surficial seafloor properties through the application angular range analysis to multibeam sonar data[J]. Marine Geophysical Researches, 2007, 28(2): 119-126.

[4] de MOUSTIER C, LONSDALE P F, SHOR A N. Simultaneous operation of the sea beam multibeam echo-sounder and the seaMARC II bathymetric sidescan sonar system[J]. IEEE Transactions on Oceanic Engineering, 1990, 15(2): 84-94.

[5] CHO H, KIM B, YU S C. AUV-based underwater 3-D point cloud generation using acoustic lens-based multibeam sonar[J]. IEEE Transactions on Oceanic Engineering, 2017, 43(4): 856-872.

[6] XU C, LI H, CHEN B, WANG X. Angular response classification of multibeam sonar based on multi-angle interval division[C]//IEEE/OES China Ocean Acoustics (COA), Harbin, China: IEEE, 2016: 1-4.

[7] XU W, SHI H,ZHANG H. Sparse-reconstruction-basedhigh resolution beamforming and its application to multi-beam systems[C]//2012 Oceans–Yeosu, Yeosu, Korea (south): IEEE, 2012: 1-4.

[8] 倪淑燕, 程乃平, 倪正中. 共軛虛擬陣列波束形成方法[J]. 電子學報, 2011, 39(9): 2120-2124.

NI Shuyan, CHENG Naiping, NI Zhengzhong. Conjugate virtual array beamforing methed[J]. Acta Electronica Sinica, 2011, 39(9): 2120-2124.

[9] 胡鵬. 虛擬陣元波束形成方法研究[D]. 西安: 西北工業大學, 2006.

HU Peng. Study on beamforing algorithm of array with virtual elements[D]. Xi’an: Northwestern Polytechnical University, 2006.

[10] 陳歡, 楊德森, 張攬月, 等. 基于遺傳算法的聲矢量陣虛擬陣元波束形成[J]. 信號處理, 2009, 25(10): 1498-1501.

CHEN Huan, YANG Desen, ZHANG Lanyue, et al. Virtual array beam forming algorithm based on genetic algorithm of vector hydrophones[J]. Signal Processing, 2009, 25(10): 1498-1501.

[11] 陳希孺. 最小一乘線性回歸(上)[J]. 數理統計與管理, 1989(5): 48-55.

CHEN Xiru. Least absolute deviation regression(vol. 1 of 2)[J] Journal of applied statistics and management, 1989(5): 48-55.

[12] LI Y B, ARCE G R. A maximum likelihood approach to least absolute deviation regression[J]. EURASIP Journal on Advances in Signal Processing, 2004(12): 1762-1769.

[13] 王福昌, 胡順田, 張艷芳. 最小一乘回歸系數估計及其MATLAB實現[J]. 防災科技學院學報, 2007(4): 85-89.

WANG Fuchang, HU Shuntian, ZHANG Yanfang. The coefficient estimation of least absolute deviation regression and implementation in MATLAB[J]. Journal of Institute of Disaster-Prevention Science and Technology, 2007(4): 85-89.

[14] 陳希孺. 最小一乘線性回歸(下)[J]. 數理統計與管理, 1989(6): 48-56.

CHEN Xiru. Least absolute deviation regression(2 of 2)[J]. Journal of Applied Statistics and Management, 1989(6): 48-56.

[15] 李貴斌. 聲吶基陣設計原理[M]. 北京: 海洋出版社, 1993.

LI Guibin. The principle of sonar array design[M]. Beijing: China Ocean Press, 1993.

Simulation of virtual array beamforming based on least absolute deviation

ZHANG Weimin1, GUO Haitao2, JIN Qiyu3, TIAN Yuanyuan4, JIAO Shengxi5

(1. College of Electronic Information Engineering, Inner Mongolia University, Hohhot 010021, Inner Mongolia, China; 2. School of Marine Information Engineering, Hainan Tropical Ocean University, Sanya 572022, Hainan, China; 3. School of Mathematical Sciences, Inner Mongolia University, Hohhot 010021, Inner Mongolia, China; 4. School of Mechanical Engineering, Northeast Dianli University, Jilin 132012, Jilin, China; 5. College of Automation Engineering, Northeast Dianli University, Jilin 132012, Jilin, China)

The quality of beamforming is the key to the imaging quality of underwater multi-beam imaging system. When the signal frequency is constant, the conventional phase-shifting beamforming algorithm can only get a narrower beam by increasing the aperture of the array, but this method is restricted by practical conditions in engineering. Therefore, a virtual array beamforming algorithm based on least absolute deviation (LAD) estimation is presented, by which a narrower and more robust beam can be obtained without increasing the array size, moreover the array gain and the angle resolution of sonar image can be increased under low signal to noise ratio (SNR). The superiority of the above method is verified by using Matlab simulation and comparing with other methods.

multi-beam imaging system; least absolute deviation estimation; robustness; virtual array; beamforming

TN911.7

A

1000-3630(2020)-02-0134-07

10.16300/j.cnki.1000-3630.2020.02.002

2018-10-04;

2018-12-17

國家自然科學基金(61661038,41076060)、內蒙古自然科學基金(2014MS0601)資助項目

張偉民(1993-), 男, 內蒙古包頭人, 碩士研究生, 研究方向為信號處理。

郭海濤,E-mail: ghtpaper@126.com

猜你喜歡
信號方法
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
學習方法
孩子停止長個的信號
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
基于LabVIEW的力加載信號采集與PID控制
一種基于極大似然估計的信號盲抽取算法
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
主站蜘蛛池模板: 午夜啪啪网| 日韩无码黄色网站| 毛片大全免费观看| 麻豆精品在线| 日本免费a视频| 国产麻豆福利av在线播放| 国产精品视频系列专区| 国产麻豆福利av在线播放| 亚洲系列无码专区偷窥无码| 中字无码av在线电影| 色国产视频| 国产精品污污在线观看网站| 色妞永久免费视频| 久久黄色免费电影| 在线不卡免费视频| 中文字幕在线永久在线视频2020| 波多野结衣一区二区三视频| 无码免费视频| 91精品专区| 欧美有码在线| 国产成人无码综合亚洲日韩不卡| 999国内精品久久免费视频| 国产精品第三页在线看| 国产精品人人做人人爽人人添| 久久精品亚洲中文字幕乱码| 欧美日韩精品一区二区视频| 91色在线观看| 国产自在线播放| 久久人体视频| 亚洲AV无码乱码在线观看裸奔 | 国产成人91精品| 日韩精品久久久久久久电影蜜臀| 国产精品久久国产精麻豆99网站| 久久人人97超碰人人澡爱香蕉| 免费观看精品视频999| 99人体免费视频| 欧美视频免费一区二区三区| 久久鸭综合久久国产| 精品久久777| 精品人妻无码中字系列| 日韩精品专区免费无码aⅴ| 最新亚洲av女人的天堂| 四虎影视8848永久精品| 国产成人亚洲日韩欧美电影| 久久久久久久久久国产精品| 亚洲日本在线免费观看| 伊人久热这里只有精品视频99| 免费无遮挡AV| 国产毛片高清一级国语| 国产最爽的乱婬视频国语对白| 国产欧美精品专区一区二区| 亚洲精品成人片在线观看| 在线精品亚洲国产| 亚洲AV成人一区二区三区AV| 精品国产中文一级毛片在线看| 成人字幕网视频在线观看| 精品人妻一区无码视频| 国产欧美日韩一区二区视频在线| 国产一级二级三级毛片| 视频国产精品丝袜第一页| 欧美在线天堂| 亚洲成aⅴ人在线观看| 国产精品污视频| 爱色欧美亚洲综合图区| 香蕉在线视频网站| 亚洲男人天堂2020| 天堂成人在线视频| 国产成人高精品免费视频| 国产成人高清精品免费| 毛片网站观看| 久久精品一卡日本电影 | 亚洲国产精品无码久久一线| 国产成人禁片在线观看| 久久久受www免费人成| 国产精品无码一区二区桃花视频| 亚洲精品不卡午夜精品| 欧美色视频日本| 中文字幕在线欧美| 日韩无码真实干出血视频| 亚欧美国产综合| 国产97视频在线| 成人国产一区二区三区|