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

基于時變曲線模型的合成孔徑聲納圖像自動均衡方法研究

2014-06-27 05:41:29劉維江澤林劉紀元黃海寧
兵工學報 2014年3期
關鍵詞:方法模型

劉維,江澤林,劉紀元,黃海寧

(中國科學院聲學研究所,北京 100190)

基于時變曲線模型的合成孔徑聲納圖像自動均衡方法研究

劉維,江澤林,劉紀元,黃海寧

(中國科學院聲學研究所,北京 100190)

針對合成孔徑聲納(SAS)圖像不均衡問題,提出一種基于時變曲線模型的SAS圖像自動均衡方法。以聲傳播模型、水底后向散射模型和SAS成像模型為基礎,推導時變曲線(TVC)的表達式;結合SAS圖像的統(tǒng)計特征,推導TVC觀測量的計算方法;用非線性最小二乘擬合方法完成TVC估計;基于TVC進行了SAS圖像的自動均衡。用湖試和海試數(shù)據(jù)對該方法進行了驗證,結果表明推導的TVC表達式與試驗數(shù)據(jù)具有較好的吻合度,提出的自動均衡方法可有效地消除SAS圖像的不均衡問題。

信息處理技術;合成孔徑聲納;圖像均衡;時變曲線;威布爾分布

0 引言

受聲傳播損失、吸收損失、水底后向散射強度、聲納發(fā)射和接收系統(tǒng)引起的聲源級波動等因素的影響,聲納圖像強度(亮度)往往會有比較大的起伏。由于不同距離對應的虛擬孔徑長度不同,合成孔徑聲納(SAS)圖像的強度變化范圍會更大。聲納圖像強度的不均衡問題對聲納圖像的判讀和處理造成很大的影響。一方面,由于圖像輸出設備(顯示屏、打印機等)的動態(tài)范圍有限,SAS圖像強度不均衡可能會造成重要細節(jié)或目標丟失。另一方面,基于SAS圖像的計算機輔助檢測或分類(CAD/CAC)方法、目標自動檢測和識別(ATD/ATR)方法大都依賴目標與背景之間的強度差異,不均衡的聲納圖像會造成CAD/CAC和ATD/ATR方法失效[1-3]。因此,圖像均衡是聲納圖像處理中的一個關鍵步驟。

關于聲納圖像均衡的方法主要包括直方圖均衡化(HE)、局部直方圖均衡化(LHE)、局部高斯均衡方法(LGM)等[4-6]。這些方法大都是根據(jù)圖像局部統(tǒng)計特性對圖像進行均衡,沒有考慮聲傳播、后向散射和成像模型,因此應用于SAS圖像時效果并不理想。為了解決這一問題,提出一種基于時變曲線模型的SAS圖像自動均衡方法。在此方法中,綜合考慮聲傳播模型、后向散射模型以及SAS圖像統(tǒng)計模型,完成時變曲線(TVC)表達式的推導、TVC觀測量的構造和參數(shù)估計以及SAS圖像的自動均衡。

1 TVC模型

側視SAS像系統(tǒng)中常采用航跡方向和垂直航跡方向構建成像坐標系。其中垂直航跡方向又稱為距離向,距離向?qū)晜鞑シ较颉_@里假定r方向為距離向,y方向為航跡方向,則SAS圖像采用I(r, y)來表示。無論是聲傳播引起的聲波能量衰減、后向散射引起的衰減以及合成孔徑長度的變化引起的圖像強度變化,均可以轉換為距離的函數(shù)f(r).假定理想的SAS圖像為I0(r,y),則實際獲取的SAS圖像為

在成像過程中,距離向一般轉換為聲傳播時間進行處理,因此圖像中與距離相關的函數(shù)f(r)常常稱為TVC.在聲信號采集過程中,一般會采用時變增益(TVG)來補償f(r)帶來的影響。由于f(r)實際上與所處水域、水底底質(zhì)等諸多因素相關,在聲信號采集過程中,f(r)很難準確估計,所以TVG的作用非常有限。

(1)式中f(r)以乘積的方式存在,而乘性噪聲處理不如加性噪聲方便,因此對(1)式取對數(shù)將乘積轉化為求和操作,如(2)式所示。

lg[I(r,y)]=lg[I0(r,y)]+lg[f(r)].(2)

為了便于描述,后續(xù)說明中I、I0和f(r)均視為對數(shù)化以后的值。

f(r)與聲傳播距離、后向散射、合成孔徑長度等因素相關,分別推導各因素的表達式,最終求和后即得到f(r)的表達式。

1.1 聲傳播模型

SAS成像中,聲傳播引起的衰減包括擴展損失和吸收損失[7],其表達式如(3)式所示。其中對數(shù)項表示擴展損失,線性項表示吸收損失,常數(shù)項表示測繪帶近端對f1(r)的影響。

1.2 后向散射模型

假定聲線與水底的夾角為θ,根據(jù)Lambert定律[8],水底的后向散射系數(shù)與θ的關系如下所示:

式中:c2為常數(shù);b2為待求解系數(shù)。無論是拖曳式平臺,還是自主航式平臺,可以認為SAS基陣距水底的深度h是緩變的,而

則(4)式可以化簡為對數(shù)項和常數(shù)項之和,如(6)式所示:

1.3 合成孔徑長度模型

假定發(fā)射陣水平向開角為β,則合成孔徑長度與距離呈正比,即

SAS系統(tǒng)無論是采用時域算法[9-11],還是采用頻域算法[12],單像素點的強度均與合成孔徑長度呈正比。如果考慮(2)式對數(shù)項的影響,則合成孔徑長度引起的f(r)分量為

綜合上述(3)式、(6)式和(8)式可以得到f(r)的表達式為

1.4 其他因素

在SAS成像中,除了1.1~1.3節(jié)的因素之外,還有兩個值得注意的因素:一是水體;二是大掠射角的情況。下面分別進行說明。

1.4.1 水體

如圖1所示,A為聲納基陣,B、C分別為聲納波束與水底交點的近端和遠端。SAS工作在側視模式下。在近距離處(小于AB之間距離),聲納波束內(nèi)只包含水體(主要是體積混響)及懸浮的小目標(如魚群等),此時回波強度較弱,一般遠遠小于水底的回波強度,因此需要對水體部分的圖像進行特殊處理。采用線性項或常數(shù)項對其近似(這一近似的有效性可以從處理結果得到驗證),即

圖1 垂直波束示意圖Fig.1 Schematic diagram of vertical beam

1.4.2 大掠射角的情況

當距離大于AB時,聲納回波中開始包含水底回波分量。在聲納波束近端附近,掠射角較大;此時,水底回波強度有一個急速上升的過程。這一變化過程并不能采用(9)式進行描述。由于變化速度非常快,而且在峰值后下降,因此可以采用二次多項式進行描述,即

1.5 TVC表達式

綜合上述分析,可以給出TVC的表達式:

如1.1~1.4節(jié)所述,不同的距離上影響聲納回波強度的因素不同,故式中采用fw、fa和flg分別表示TVC.如圖1所示,水體回波與水底回波的分界點為B.理想情況下,r1為線段AB的長度;實際情況下,影響r1和r2的因素很多,如聲納基陣姿態(tài)、水底地形、水深、聲納基陣垂直向開角、SAS基陣安裝角等。

2 TVC估計方法

r1處TVC應當滿足連續(xù)性條件:

為了保證聲納圖像的連續(xù)性,式中fa與flg應當是連續(xù)變化的,故在r2處應當滿足連續(xù)性和光滑性條件,如(14)式所示:

利用(13)式和(14)式,可以將f(r)的參數(shù)數(shù)量降低為7個,即a3、a2、a1、b1、c1、r1和r2.通過SAS圖像得到時變曲線f(r)的觀測量g(r),在此基礎上利用非線性最小二乘擬合的方法[13]估計參數(shù)向量p,待估計量采用p*表示,則

(16)式可以采用置信區(qū)間(TS)法求解,其基本步驟如下:

1)在坐標點p的鄰域s內(nèi)采用函數(shù)q近似E(p,r);

2)限定s≤Δ,根據(jù)局部最小準則計算最佳步長s*;

3)如果

則接受p+s*作為新的坐標點,同時增大鄰域范圍Δ.否則,縮小鄰域范圍Δ;

4)重復執(zhí)行上述步驟1~步驟3直至得到p*.

對于函數(shù)E(p,r),可以采用E(p,r)的二階泰勒展開近似,即

式中:H為函數(shù)E(p,r)的二階導數(shù)矩陣(Hessian矩陣);g為函數(shù)E在p處的梯度,且

計算步長s是TS方法的關鍵,根據(jù)步驟2的局部最小準則和(19)式,可得

(21)式可采用特征方程法求解(具體見文獻[14])。

3 構造TVC的觀測量

獲取TVC測量最簡單的方法即從SAS圖像中沿距離方向取值。圖2中給出一個沿距離方向直接獲取的TVC觀測量的實例,可以看出直接獲取的TVC觀測量噪聲幅度非常大。圖3給出了一幅SAS圖像不同距離處像素值對應的概率密度曲線,可以看出不同距離的SAS圖像像素值的分布形式差異較大。由圖2和圖3可知,沿距離方向直接獲取TVC觀測量的方法存在較大缺陷,會對參數(shù)向量p的魯棒估計造成一定的困難。

圖2 TVC觀測量Fig.2 Measurement of TVC

圖3 SAS圖像像素值概率密度曲線Fig.3 Probability density curves of SAS image pixel values

為了獲取最佳的TVC觀測量,可以利用SAS圖像的統(tǒng)計特征。根據(jù)(1)式、(2)式可知,同一距離上的SAS圖像像素受f(r)的影響相同,因此可以認為其具有相同的分布。定義圖像像素集合Ir和像素值vr(y)?Ir,

則根據(jù)Ir可以估計Weibull分布的參數(shù)λr和kr.

Weibull分布為非對稱分布(如圖3所示),因此眾數(shù)是TVC最佳的觀測量,即

可以根據(jù)λr和kr計算Weibull分布的眾數(shù),但估計這兩個參數(shù)的計算量較大。由文獻[15]可知, Weibull分布的眾數(shù)和中值分別為

根據(jù)(1)式和(2)式可知,影響SAS圖像像素值TVC的f(r)隨距離變化,與航跡方向y無關。對于同一水下區(qū)域的聲納圖像,形狀參數(shù)kr沿距離方向的變化不大,而尺度參數(shù)λr沿距離方向的變化較大。因此,在實際應用中可以采用中值濾波代替眾數(shù)估計。此時,

式中:由于形狀參數(shù)kr變化不大,因此可以近似認為γr為常數(shù)。這一點可以從圖4和圖5中看出。圖4和圖5中給出了某一水底區(qū)域SAS圖像像素值Weibull分布的形狀參數(shù)kr和尺度參數(shù)λr隨距離的變化。從圖4和圖5中可以看出,形狀參數(shù)隨距離變化不明顯,而尺度參數(shù)隨距離變化明顯。

圖4 形狀參數(shù)krFig.4 Shape parameter kr

圖5 尺度參數(shù)λrFig.5 Scale parameter λr

4 SAS圖像均衡

利用(16)式、(24)式和(27)式可以得到TVC的f(r,p*),簡記為f*(r).利用f*(r)通過(28)式完成SAS圖像均衡,得到均衡后的圖像I0(r,y).

由于不同底質(zhì)對應的回波強度有差異,這種差異對于底質(zhì)分類等應用有比較重要的參考作用。因此,為了保證回波強度信息不丟失,在應用(28)式時,增加了β參數(shù)。

由于人眼或計算機可表示的動態(tài)范圍有限,因此經(jīng)過(28)式處理后的圖像I0(r,y)一般還需要經(jīng)過對比度增強才能達到最佳的效果。調(diào)整的依據(jù)即Weibull分布的尺度參數(shù)λ.假定圖像采用[0 1]之間的浮點數(shù)進行表示,用戶期望目標相對背景的信噪比優(yōu)于S,則預期的背景強度為

如果采用vmode表示I0(r,y)的眾數(shù),用于度量I0(r,y)背景的強度,圖像增強變換為

根據(jù)(23)式和(25)式可知,(31)式完成的圖像增強變換相當于改變圖像分布的尺度參數(shù)λ.經(jīng)過均衡處理以后,圖像I0(r,y)在不同距離上的像素值對應的Weibull分布參數(shù)相近,即不同距離上圖像像素的形狀參數(shù)和尺度參數(shù)相近。因此,只需要調(diào)整分布的尺度參數(shù)即可以達到對比度增強的目的。

5 試驗數(shù)據(jù)處理及分析

提出的SAS圖像自動均衡方法在SAS圖像顯示系統(tǒng)中已經(jīng)取得應用,應用的水域包括千島湖、渤海、東海、南海和阿曼灣等,涉及多種類型的水底底質(zhì)。下面選擇3組典型的湖試和海試試驗數(shù)據(jù)對提出的SAS圖像自動均衡方法的有效性進行說明。其中第一組數(shù)據(jù)為湖上試驗數(shù)據(jù),對應湖底地貌成像和水體中魚群的成像結果;第二組數(shù)據(jù)為湖上試驗數(shù)據(jù),對應嚴重不均衡的湖底地貌成像結果;第三組數(shù)據(jù)對應海上試驗數(shù)據(jù),對應海底地貌和海底管線成像結果。通過這3組湖試和海試試驗數(shù)據(jù),可以對提出的SAS圖像自動均衡方法的有效性進行較全面的驗證。

5.1 試驗數(shù)據(jù)處理結果

在試驗數(shù)據(jù)處理的結果中,給出了3類數(shù)據(jù):一是均衡前和均衡后的SAS圖像對比;二是TVC觀測量與估計量曲線;三是均衡前和均衡后SAS圖像中不同距離上像素值統(tǒng)計分布的概率密度曲線。

1)第一組試驗數(shù)據(jù)處理結果,包括圖6~圖9.

2)第二組試驗數(shù)據(jù)處理結果,包括圖10~圖12.

圖6 SAS圖像(橫軸為距離向,縱軸為沿航跡方向;圖像中黑色部分為水體,圖像寬度147 m,長度281 m)Fig.6 SAS image(abscissa:range,ordinate:track.Black regions in SAS images are related with water reflection.The images width is 174 m,and the image length is 281 m)

圖7 TVC觀測量和估計量Fig.7 Measurement and estimation of TVC

3)第三組試驗數(shù)據(jù)的處理結果,包括圖13~圖15.

圖8 圖6中魚群圖像A放大Fig.8 Enlarged image A of fish school in Fig.6

圖9 沿距離向圖像像素值的分布Fig.9 Distribution of image pixel values in range direction

5.2 對試驗數(shù)據(jù)處理結果的分析

在試驗數(shù)據(jù)分析時,從三個方面對提出的SAS圖像自動均衡方法的有效性進行檢驗,即:1)TVC觀測量與估計量的吻合程度;2)不同距離圖像像素值的統(tǒng)計分布;3)對SAS圖像判讀的影響。

圖10 SAS圖像(橫軸為距離向,縱軸為沿航跡方向;圖像中黑色部分為水體,圖像寬度113 m,長度196 m)Fig.10 SAS image(abscissa:range,ordinate:track.Black regions in SAS images are related with water frelection.The image width is 113 m,and the image length is 196 m)

圖11 TVC觀測量和估計量Fig.11 Measurement and estimation of TVC

通過對照TVC估計量與觀測量的吻合程度可以檢驗TVC模型的有效性和適用性。對照圖7、圖11和圖14可以看出,TVC的估計量與TVC觀測量的吻合程度都非常好。這說明,文中構造的TVC模型和提出的TVC表達式與試驗結果吻合較好,且具有較強的適用性。在圖14中,盡管有海底管線強干擾的存在,仍然可以得到TVC的估計量,這說明文中給出的TVC參數(shù)的估計方法具有較強魯棒性。

圖12 沿距離向圖像像素值的分布Fig.12 Distribution of image pixel values in range direction

從物理上看,對于同一成像區(qū)域,如果底質(zhì)類型大致相同,其SAS圖像像素值的分布應當趨于一致。通過均衡處理,消除聲傳播、后向散射、發(fā)射和接收電子系統(tǒng)引起的不均衡現(xiàn)象,對聲圖判讀和后續(xù)處理是非常有益的。對照圖9、圖12和圖15可以看出,在均衡前SAS圖像像素值的分布隨距離的變化劇烈波動;而均衡后SAS圖像像素值的分布趨于一致。

均衡對于SAS圖像判讀的影響在圖6、圖10和圖13中表現(xiàn)形式有所差異。在圖6中,均衡后的圖像可以清晰看到水體中目標強度比較弱的懸浮目標(魚群等,圖8為圖6中部分SAS圖像的局部放大)。圖10中,受水下地勢的影響,SAS圖像中距離較遠的部分回波較弱;通過均衡處理,使原本較暗的區(qū)域變得清晰可讀;圖10中,由于SAS圖像中近距離和遠距離上底質(zhì)存在較大差異,所以均衡后近距離和遠距離圖像像素值的分布仍然存在一定的不同。圖13中,由于水下管線的目標強度較高,再加上SAS圖像固有的不均衡問題,使SAS圖像中除管線之外的大部分細節(jié)被壓制,非常不利于圖像判讀(見圖13(a));通過均衡處理,有效地消除了這一問題,使管線附近的作業(yè)痕跡清晰地顯現(xiàn)出來。

綜合上述分析,可以看出文中推導的TVC模型合理,TVC表達式與試驗結果有較好的吻合度,提出的自動均衡方法有效地消除了SAS圖像中存在的不均衡現(xiàn)象,為SAS圖像判讀和后續(xù)處理打下了良好的基礎。

圖13 SAS圖像(橫軸為距離向,縱軸為沿航跡方向;SAS圖像中黑色部分為水體,圖像寬度188 m,長度389 m)Fig.13 SAS image(abscissa:range,ordinate:track.Black regions in SAS images are related with water reflection.The image width is 188 m,and the image length is 389 m)

圖14 TVC觀測量和估計量Fig.14 Measurement and estimation of TVC

圖15 沿距離向圖像像素值的分布Fig.15 Distribution of image pixel values in range direction

6 結論

密切結合物理模型,提出一種SAS圖像自動均衡方法,并在實際應用中具有良好的效果。此方法的關鍵有兩點:一是合理的TVC模型。通過湖試和海試數(shù)據(jù)的處理和分析,證明了推導的TVC表達式的合理性;二是TVC的觀測量和估計量的計算。結合SAS圖像統(tǒng)計特征,通過理論分析,表明可以采用Weibull分布的眾數(shù)作為TVC觀測量。利用非線性最小二乘擬合方法實現(xiàn)了TVC的魯棒估計。

該方法可以應用于SAS圖像后處理和實時處理。在實時處理環(huán)境中,聲納圖像隨聲納基陣的移動連續(xù)輸出,在這種情況下,可以通過沿航跡方向加窗的方式實現(xiàn)聲納圖像的實時自動均衡。

References)

[1] Coiras E,Myers V,Evans B.Reliable seabed characterization for MCM operations[C]∥OCEANS 2007.Vancouver,BC:IEEE, 2007:1-5.

[2] Fandos R,Zoubir A M.Optimal feature set for automatic detection and classification of underwater objects in SAS images[J].Selected Topics in Signal Processing,2011,5(3):454-468.

[3] Maussang F,Chanussot J,Hetet A.Mean-standard deviation representation of sonar images for echo detection:application to SAS images[J].Oceanic Engineering,2007,32(4):956-970.

[4] Dobeck G J.Image normalization using the serpentine forwardbackward filter:application to high-resolution sonar imagery and its impact on mine detection and classification[C]∥Defense and Security.Orlando,Florida,US:International Society for Optics and Photonics,2005:381-391.

[5] Capus C G,Banks A C,Coiras E,et al.Data correction for visualisation and classification of sidescan SONAR imagery[J].IET Radar,Sonar&Navigation,2008,2(3):155-169.

[6] Isaacs J C,Tucker J D.Signal diffusion features for automatic target recognition in synthetic aperture sonar[C]∥Digital Signal Processing Workshop and IEEE Signal Processing Education Workshop(DSP/SPE).Sedona,AZ:IEEE,2011:461-465.

[7] 李啟虎.數(shù)字式聲納設計原理.[M]合肥:安徽教育出版社, 2002:164-167.

LI Qi-hu.Principle of digital sonar design[M].Hefei:Anhui Educational Publishing House,2002:164-167.(in Chinese)

[8] 北大西洋公約組織水下研究中心,美國PSI公司.海洋水聲環(huán)境和聲納設計手冊[M].北京:海潮出版社,2010:129.

NATO Underwater Research Center,US PSI Corporation.Ocean acoustic environment and sonar design manual[M].Beijing:Haichao Press,2010:129.(in Chinese)

[9] 劉維.多子陣SAS三維模型仿真與數(shù)據(jù)處理算法研究[D].北京:中國科學院聲學研究所,2008:45-60.

LIU Wei.Research on multiple-receiver SAS three dimensional model simulation and data processing algorithms[D].Beijing:Institute of Acoustics,Chinese Academy of Sciences,2008:45-60. (in Chinese)

[10] 劉維,張春華,劉紀元.FFBP算法在合成孔徑聲納成像中的應用[J].聲學技術,2009,28(5):572-576.

LIU Wei,ZHANG Chun-hua,LIU Ji-yuan.Application of FFBP algorithm to synthetic aperture sonar imaging[J].Acoustical Technology,2009,28(5):572-576.(in Chinese)

[11] Ferguson B G,Wyber R J.Generalized framework for real aperture,synthetic aperture,and tomographic sonar imaging[J].O-ceanic Engineering,2009,34(3):225-238.

[12] 劉維,劉紀元,張春華.多子陣合成孔徑聲吶波數(shù)域算法不均勻采樣問題研究[J].聲學學報,2009,34(3):203-210.

LIU Wei,LIU Ji-yuan,ZHANG Chun-hua.Research on non-uniform sampling problem when adapting wave number algorithm to multiple-receiver synthetic aperture sonar[J].Acta Acustica, 2009,34(3):203-210.(in Chinese)

[13] Madsen K,Nielsen H B,Tingleff O.Methods for non-linear least squares problems[M].2nd ed.Copenhagen:University of Denmark,2004:17-50.

[14] Com A R,Gould N I M,Toint P L.Trust region methods[M]. Philadelphia:SIAM,2010:115-168.

[15] Pham Hoang.Springer handbook of engineering statistics[M]. New York:Springer,2006:64-72.

A Time-variant Curve Model-based Automatic Equalization Method for Synthetic Aperture Sonar Images

LIU Wei,JIANG Ze-lin,LIU Ji-yuan,HUANG Hai-ning
(Institute of Acoustics,Chinese Academy of Sciences,Beijing 100190,China)

A time-variant curve(TVC)model-based automatic equalization method is proposed for the intensity variation problem of synthetic aperture sonar(SAS)images.A theoretical expression of TVC is deduced based on sound transmission model,underwater backscattering strength model and SAS image model.A method to calculate the TVC observations is established based on the statistical model of SAS images.An estimation method based on non-linear least square fitting model(NL-LSFM)is used to acquire the optimized parameter of TVC.At last,the SAS images are automatically equalized and enhanced based on the optimized TVC.The method proposed has been validated by lake and sea trials.The test results show that the result calculated by the TVC expression is consistent with the experimental data,and the automatic equalization method can remove the intensity variation of SAS images properly.

information processing;synthetic aperture sonar;image equalization;time-variant curve; Weibull distribution

TB566

:A

1000-1093(2014)03-0347-08

10.3969/j.issn.1000-1093.2014.03.009

2013-04-08

國家自然科學基金項目(11204343);哈爾濱工程大學水下機器人技術重點實驗室基金項目(9140C27020112022601)

劉維(1980—),男,副研究員,博士。E-mail:liuwei@mail.ioa.ac.cn;江澤林(1985—),男,博士研究生。E-mail:jiangzelin1985@126.com

猜你喜歡
方法模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
學習方法
3D打印中的模型分割與打包
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
主站蜘蛛池模板: 中文字幕在线观| 91啦中文字幕| 亚洲第一av网站| 一级片免费网站| 波多野结衣中文字幕一区二区| 成人欧美日韩| 亚洲高清国产拍精品26u| 国产精品久久久久久久久kt| 茄子视频毛片免费观看| 91精品专区国产盗摄| 再看日本中文字幕在线观看| 国产96在线 | 国产综合欧美| 亚洲国产成人久久精品软件| 91福利一区二区三区| 亚洲激情99| 久久熟女AV| 亚洲欧美不卡视频| 日韩小视频在线观看| 久久综合伊人77777| 国产精品香蕉在线| av天堂最新版在线| 色综合天天综合| 无码专区国产精品一区| 99久久精品免费视频| 一级毛片网| 亚洲精品自产拍在线观看APP| 欧美日韩精品一区二区视频| 三级国产在线观看| 亚洲最新网址| 一级毛片中文字幕| 久久77777| 亚洲高清中文字幕在线看不卡| 久久精品国产精品青草app| 日本人真淫视频一区二区三区| 青青草欧美| 亚洲国产一区在线观看| 无码内射中文字幕岛国片 | 久久天天躁狠狠躁夜夜2020一| 五月婷婷导航| 久久综合成人| 国产在线观看第二页| 青青草a国产免费观看| 午夜毛片福利| 色视频国产| 无码高潮喷水在线观看| 美女内射视频WWW网站午夜 | 性视频久久| 伊人久久婷婷| 亚洲精品自产拍在线观看APP| 多人乱p欧美在线观看| 国产精品片在线观看手机版 | 无码专区在线观看| 亚洲福利片无码最新在线播放| 色综合久久无码网| 国产高清无码麻豆精品| 激情六月丁香婷婷| 91小视频在线播放| 国产a v无码专区亚洲av| 精品一区二区三区视频免费观看| 最新国产精品鲁鲁免费视频| 国产 在线视频无码| 亚洲精品无码成人片在线观看| 韩日午夜在线资源一区二区| 久久免费精品琪琪| 中文字幕亚洲专区第19页| 国产精品所毛片视频| 欧洲亚洲欧美国产日本高清| 国产欧美成人不卡视频| 亚洲区视频在线观看| 久久无码av三级| 福利一区在线| 国产精品亚洲αv天堂无码| 国产又爽又黄无遮挡免费观看 | 亚洲视频无码| 视频在线观看一区二区| 欧美综合激情| 99re热精品视频中文字幕不卡| 国产一区成人| 国产成人a毛片在线| aa级毛片毛片免费观看久| 国产日本欧美在线观看|