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

基于多域特征提取和深度學習的聲源被動測距?

2021-04-22 02:49:22肖旭王同王文博蘇林馬力任群言
應用聲學 2021年1期
關鍵詞:特征信號模型

肖旭 王同 王文博 蘇林 馬力 任群言

(1 中國科學院水聲環境特性重點實驗室 北京 100190)

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

(3 中國科學院大學 北京 100049)

0 引言

聲源被動測距作為聲吶系統的重要功能之一,一直是水聲工作者密切關注的問題[1]。由于水聲觀測信號受復雜的時、空、頻變及強多途、高噪聲和多普勒效應等因素影響,傳統的匹配場處理方法往往面臨計算量大和環境失配等問題。近年來,深度學習作為基于數據驅動方式的新興分支,以其強大的特征提取能力和高效處理復雜、高維、非線性系統問題的獨特優勢[2],為水聲被動測距提供了一種新思路。

深度神經網絡通過大量數據樣本建立高維參數之間的復雜非線性映射,適用于物理建模困難的問題,引發了水聲研究者的關注。Lefort等[3]通過水箱實驗研究了機器學習算法在水聲目標測距中的測距性能。Niu等[4]利用單水聽器和殘差卷積神經網絡對聲源進行定位。Wang等[5]提出了一種用于水下聲源測距的深度遷移學習方法,將從仿真環境獲得的預測能力遷移到實驗海域。

水聲接收信號包含聲源和聲信道的大量信息,其特征提取和構造是深度學習方法的關鍵環節。早期特征提取通常利用信號的自相關函數和功率譜估計,或采用時-頻分析方法來提取一些時頻聯合域特征,如短時傅里葉變換(Short-time Fourier transform,STFT)、Wingner-Ville分布等。然而,無論是功率譜分析還是時頻分析[6],它們包含大量與聲源位置不相關的信息或冗余信息,在形式上維數較大,一般無法直接應用于測距任務。而且,單一地采用某類特征通常會丟失掉部分特征,缺乏全面性。

針對以上問題,本文設計了一種基于多域特征提取和深度學習的方法來實現聲源被動測距。首先從聲信號中提取多域特征,包含時域波形結構特征、時域包絡特征、頻域譜特征和基于STFT的時頻聯合域特征。然后基于不同譜表達計算出一組聲學參數構成特征空間,在此基礎上采用最大相關-最小冗余準則(Maximum relevance and minimum redundancy,mRMR)選擇特征空間中重要度高的關鍵特征(與聲源位置相關性大)作為模型輸入,最后通過一種改進的深度神經網絡實現聲源距離的估計。神經網絡訓練時采用自適應矩估計(Adaptive moment estimation,Adam)優化算法和均方誤差(Mean squared error,MSE)代價函數進行更新模型參數,用L2和Dropout正則化策略實現網絡參數正則化。通過淺海環境仿真實例驗證了該方法的有效性,對比分析了波形參數對測距性能和模型收斂速度的影響。

1 理論模型

1.1 水聲信號多域特征提取

聲信號的多域特征可歸納為6類:時域波形結構特征、時域包絡特征、頻域譜特征、基于STFT的時頻聯合域特征、基于等效矩形帶寬的聽覺譜特征和基于正弦諧波模型的諧波譜特征[7]。每一類特征均包含了對應譜的多種聲音特性,這些特征屬性無法用單一尺度進行描述,只有在多維特征空間下才能表示。Peeters等[8]對這些聲學特征進行總結,并通過多維標度法提取了合適的聲學參數,使之與聲信號在每個維度上的坐標呈現較大的相關性。

這些聲學參數是基于不同譜表達計算的,其表征的物理含義各有不同,而水下聲源物理特征和聲場信息主要包含在信號的時域波形、頻域能量分布中。因此,本文綜合了以上更符合水聲信號特點的前4類特征,并將所對應的特征歸納為時域特征和時頻聯合域特征。

1.1.1 時域特征

水聲信號時域波形反映了信道對聲信號傳播的畸變作用,是獲取聲源總體特征最直接的來源。為了直接從時域提取特征,自相關系數是一種廣泛使用的分類特征。首先是對原始信號s(tn)求自相關系數,tn代表信號時刻,保留前N維的自相關系數(c∈{1,···,N})表示為[9]

其中,Ln是分幀時的窗長,c代表時間的滯后量。當聲源信號為瞬態聲信號時,其隨時間變化經歷ADSR過程,即激勵階段、衰減階段、穩態持續階段、釋音殘響階段,如圖1所示,其中激勵階段到衰減階段以信號振幅峰值處為分界點,后三階段統稱為下降階段。

圖1 聲信號ADSR過程Fig.1 ASDR process of sound signals

聲源時域信號的幅度強弱變化和接收距離的相關性較強,為提取描述聲源信號幅度變化的特征,對原始信號s(tn)進行Hilbert變換,然后使用截止頻率為5 Hz的三階Butterworth濾波器對振幅信號進行低通濾波,得到信號振幅包絡e(tn)。估計信號的起始時間tst、激勵階段和下降階段的終止時間tAend、tDend,在此基礎上定義時域的7個聲學參數,如表1所示。其中,為了估計激勵過程的時間長度,定義對數激勵時間為[10]

相對應地,將激勵階段和下降階段能量的平均時間斜率分別定義為激勵斜率和下降斜率。估計振幅包絡的最大值emax。時間質心給出了信號能量質心所在的時刻,其定義式為[11]

其中,n1和n2為聲信號的起始和終止時間對應的索引號,以濾除前后的空白時間。當信號聲能量大于一定閾值γemax時,其持續時間定義為信號的有效時間,根據許多經驗性的測試,γ取0.4性能較穩定,見文獻[8]。除以上參數外,其余參數的公式定義可在文獻[8]中一一得到。這些聲學參數反映了聲源的各種物理屬性,例如對數激勵時間反映了能量在上升過程中的時間長度,其與聲源距離呈正相關,聲源距離越遠,聲信號需達到穩態振動的時間越長。與波形包絡相關的特征(激勵時間、時域質心等)僅適用于瞬態信號,分析連續聲信號時,應提取自相關系數、載波信號調制和后文定義的頻域特征。

表1 時域特征Table 1 Temp oral features

1.1.2 時頻聯合域特征

對原始信號進行STFT,把每一幀進行快速傅里葉變換后的頻域信號在時間上堆疊起來得到時頻譜ak(tm),其中m代表幀數,k代表頻點數。對ak(tm)進行標準化處理:

令第k個頻點的頻率為fk,得到譜的前四階統計矩,分別定義為譜質心、譜延展、譜斜度和譜峰度[12],對于第m幀的頻譜,表示為

除了頻譜的統計特征,根據頻譜的斜率特征可以從STFT能量譜中提取譜斜率、譜衰減[8]、譜滾降[13]、譜通量[14]。最后根據線譜特征提取信號的平坦度及譜峰度[8],其定義和物理含義如表2所示。對于第m幀,上述物理量計算公式分別為

其中,?f為信號STFT后兩點的頻率差,fmax為奈奎斯特采樣決定的最高頻率,是令ak(tm)最大的k(tm)的值。以上參數被證實對聲信號的特征識別具有重要作用[15],被試信號對于它們的變化具有較高的敏感度。

表2 時頻聯合域特征Table 2 Temp oral-frequency features

1.2 基于mRMR準則的特征選擇

以上聲學參數作為水聲信號的輸入特征并不具有魯棒性,不同任務(如測距、識別)的訓練集擁有不同的最佳聲學參數。在1.1節中給出的特征中,有些特征可能是冗余的甚至是不相關的,導致機器學習算法的效率降低、性能損失。

最大相關-最小冗余準則(mRMR)是一種綜合考慮特征相關度和冗余度的特征重要性評價準則[16]。定義互信息I(A,B):

其中,變量A和B的概率密度分別是p(A)和p(B),其聯合概率密度是p(A,B)。設樣本數量為m,特征向量數量為n,特征向量fi=[f(i,1),f(i,2),···,f(i,m)]T,I(fi,fj)為樣本中第i個和第j個特征的相關性,其中i,j=1,2,3,···,n。設Om為類別標簽,I(fi,O)為特征與輸出類別O的相關性,其中向量O=[O1,O2,O3,···,Om]T。利用最大相關標準式選擇出與類別O相關性大的特征集合D:

式(8)中,|S|為集合S中所選特征的數量。利用最小冗余標準式剔除特征子集S中的冗余特征的集合R:

綜合以上條件,mRMR方法計算式為

給定具有N?1個特征的集合SN?1,總特征集合為F,計算集合{F?SN?1}中選擇第N個特征使得式(10)中的集合θ(D,R)最大:

利用mRMR準則對特征空間進行預處理,可以剔除冗余特征,降低計算代價,產生緊湊性和泛化能力更強的模型。算法流程如下:

(1)選擇令相關性最大的特征fn,即maxfnI(fn,O),將所選特性特征添加到空集合S中。

(2)在集合S的補集中找出具有非零相關性和零冗余的特征,如不包含,則轉步驟(4);否則,選出相關性最大的特征fk,即將選中的特征添加到集合S中。

(3)重復步驟(2),直到S的補集中所有特征的冗余不為零為止。

(4)選擇S的補集中互信息熵最大且具有非零相關性和非零冗余的特征fl,即,將選擇的特征加入集合S中。

(5)重復步驟(4),直到S的補集中所有特征的相關性為零。

(6)最后以隨機順序添加與S無關的特征。

1.3 改進的深度神經網絡模型

1.3.1 傳統前饋深度神經網絡

傳統的前饋深度神經網絡(Feedforward deep neural network,FF-DNN)根據內部的神經網絡層可以分為輸入層(輸入聲信號特征的層)、隱含層(所有中間層)和輸出層(輸出目標距離估計值的層)。單層網絡直接相互級聯,某一層的任意一個神經元與其上一層的每一個神經元相連。其局部模型可描述為是一個線性運算加上一個非線性轉移函數。

設深度神經網絡層數為L,l是每一層的索引號,x(l)、y(l)分別是第l層的輸入序列和輸出序列,網絡的輸入y(0)=x,w(l)和b(l)為第l層的權重矩陣和偏置向量,f(z)表示非線性轉移函數。標準的FF-DNN網絡描述如下(對于l∈{0,1,···,L?1層的第i個神經元})[17]:

訓練模型時,利用尋優算法對距離估計的代價函數進行迭代優化求極小值,找到合適的線性系數矩陣w和偏置向量b:

其中,J為代價函數,通常采用輸出層輸出的目標距離估計值與真實距離之間的均方誤差:

其中,N是聲信號訓練樣本數,j是其樣本索引號,zj是對應樣本的真實距離。關于上述求解優化問題,常使用梯度下降法、共軛梯度法、擬牛頓法等數值優化方法。由于求Hessian矩陣及其逆計算量十分巨大,最常用的優化算法仍然是梯度下降算法。

在聲源測距中,傳統DNN存在以下缺陷:

(1)測距誤差較大。聲源距離的代價函數可能高度非凸,迭代過程中容易陷入局部次優解或鞍點。

(2)算法收斂速度慢。梯度下降法的初始學習率和調整策略需人工調節,相同的學習率被應用于各個參數,效率低下。

(3)模型泛化性和魯棒性弱。全連接網絡的模型復雜度過高,參數稀疏度過低,易發生過擬合,以至于對環境變化和信號畸變過于敏感。

要解決以上問題,提高測距性能,重點在于如何改進尋優策略、加快收斂速度、防止過度學習。為此,本文引入1種自適應動態調整學習率的優化方法和2種網絡參數稀疏化技術來改進網絡模型。

1.3.2 Adam優化算法

由于水下噪聲、混響和水聲信道的多途干擾,目標距離的代價函數通常為非凸函數且局部次優解,從而需要較大的學習率來跳出局部最優。然而,當在全局最優值附近搜索時,學習率太大會導致過度學習,降低聲源測距精度。

Adam算法是一種動態調整參數學習率的自適應優化方法[18]。該方法通過梯度的一階和二階矩估計動態調整各網絡參數的學習率,在迭代過程中通過偏差糾正使學習率維持在一定范圍,從而獲得平穩的參數更新,這是解決聲源測距問題的理想方法。

設t為迭代次數,w為待估參數,J為代價函數,首先計算梯度的指數移動平均數mt。m0初值為0。綜合考慮之前時間步的梯度動量,設系數β1為指數衰減率,有

計算梯度平方的指數移動平均數,v0初始化為0。設系數β2為指數衰減率,有

m、v初始化為0會導致mt偏向于0,因此先進行偏差糾正再更新參數:

式(17)中,η為初始學習率。算法對更新的步長計算從梯度均值及梯度平方兩個角度進行自適應的調節[19],起到提高迭代效率和測距精度的作用。

1.3.3 網絡參數稀疏化

傳統的DNN往往受限于特定的水聲環境,對環境變化和信號畸變過于敏感,出現過擬合現象。具體表現在迭代過程中訓練誤差下降到一定程度時,測試誤差反而開始增大。為了生成泛用性強的模型,將數據映射到網絡特征后,網絡特征之間的重疊信息應盡可能少,相關性盡可能低,從而近似于標準正交基。其主要方法是使特征產生稀疏性:稀疏特征有更大可能線性可分,或者對非線性映射機制有更小的依賴[20]。

L2正則化是一種簡單且有效的網絡參數稀疏化方法。在式(14)加入懲罰項,通過懲罰因子λ控制網絡參數稀疏度:

Dropout正則化策略是另一種神經網絡稀疏化手段,其核心在于每次權重更新迭代中,對網絡的每一層,隨機將部分節點對應的權重值置零,使得線性系數矩陣和偏置向量達到稀疏化的效果,在一定程度上避免過擬合的問題。引入Dropout的神經網絡描述由式(12)變為[17]:

式(19)中,符號·表示向量點乘,r(l)是一個向量,其元素為服從伯努利隨機分布的隨機變量,分別以概率P和1?P取1和0為值,參數P是每個神經元的激活概率,通常P取[0.5,1.0]。使用該向量對上一層網絡的輸出y(l)進行采樣,產生一個約減的輸出?y(l)用于下一次網絡的輸入。這個操作依次進行,從而可以生成一個稀疏的網絡結構[17]。這種方法簡單易行、節省運算資源,且不會提升優化過程的復雜度。改進后的DNN能夠從有效的數據維度上學習到相對稀疏的特征,達到自動提取水聲信號關鍵特征的效果。

2 數值仿真

通過KRAKEN聲場計算工具,在聲速正梯度淺海環境參數下生成仿真數據。圖2描述了本文所使用的環境參數。

圖2 環境參數Fig.2 Enviromental parameters

仿真數據包括連續波(Continuous wave,CW)在50 Hz、150 Hz和300 Hz的信號,線性調頻(Linear frequency modulation,LFM)信號中心頻率為500 Hz、1000 Hz和2000 Hz,頻帶寬度范圍100~1000 Hz,信號長度0.2~1.0 s。將模擬接收信號拷貝100份并分成10組,每一組模擬接收信號分別添加信噪比(Signal-to-noise ratio,SNR)為1 dB、2 dB、3 dB、···、10 dB的高斯噪聲。接收點距離分布在1~10 km,深度是分布在5~145 m,網絡輸入訓練集占總樣本集80%,由16080個樣本組成,剩余20%的數據作為測試集,由4020個樣本組成。對生成樣本進行多域特征提取,對每一幀得到的特征序列求統計特征,得到均值和方差,對所有幀的自相關系數和頻域特征序列求統計特征,得到所有幀的時間均值和時間方差。最終提取到20100個樣本的36維特征,作為DNN網絡的輸入特征,特征空間如圖3所示。

圖3 特征空間Fig.3 Feature space

每個特征之間并不完全獨立,有些特征與聲源距離顯著相關,這里根據1.2節中的mRMR算法進行特征選擇。基于互信息的mRMR最高效和常用的[21],在步驟(4)中,輸出所選特征的互信息熵作為特征重要性的評價指標,對特征空間上所有特征進行重要度排序,結果如圖4所示。

圖4中的符號和表1、表2中一一對應,例如DS是下降斜率,SV_Mean是對信號每一幀的譜通量求的平均值;AC_Std是對每一幀自相關系數求的標準差。由排序結果可見與聲源距離相關性最強的前3項特征是激勵時間、下降斜率和譜通量均值,分別代表激勵階段的時間長度、衰減階段能量的平均時間斜率和頻譜包絡面積的均值,這些物理量恰是反映聲能量在傳播過程中衰減的基本物理量。通常,為了兼顧特征集合的多樣性和緊湊性,指標的閾值不宜過大或過小,經測試這里取0.03時,特征子集的維度為29,此時模型收斂性較好。最終得到與聲源距離相關性最高的10維時域特征與19維頻域特征。

圖4 聲學參數mRMR重要性排序Fig.4 The mRMR importance ordering of acoustic parameters

神經網絡引入Adam優化算法進行網絡訓練,初始學習率采用0.03,代價函數采用MSE函數,引入Dropout正則化處理,每次迭代神經元激活概率取85%,初始權重由截斷高斯分布模型產生,標準差為0.1。經測試雙曲正切函數、Sigmoid函數和ReLU函數作為隱含層激活函數均未出現梯度消失現象,其中雙曲正切函數在本問題中表現的收斂速度最快,且訓練過程中未出現死神經元。隱含層激活函數采用雙曲正切函數,輸出層采用200個Softmax節點,對應不同的距離的概率,輸出值最大的節點對應的距離為距離估計值。網絡迭代次數設置為20000次。不同波形參數的單頻信號和線性調頻信號作為發射信號訓練的DNN經20000次迭代后在測試集上的綜合測距精確率達到95%以上,最高達到98%以上。以其中一組波形參數(單頻信號,f0=150 Hz,zr=35 m,SNR=1~10 dB)的訓練和測試結果為例,圖5為該組信號的模型訓練和測試結果。其中圖5(a)給出了訓練完成后模型最終在測試集上的距離估計結果,紅線代表KRAKEN聲場模型中給定的聲源距離(即真實距離),藍圈代表網絡輸出的估計距離;圖5(b)給出了模型在訓練集和測試集上的測量精度隨迭代次數的變化曲線。

圖5 訓練和測試結果Fig.5 The ranging accuracy by iteration times on validation and training sets

以單頻信號作為發射信號,改變發射信號頻率(f=50 Hz,100 Hz,150 Hz)和聲源深度(zS=0~140 m,?z=20 m),經10000次迭代,對比不同發射條件下的估計精確率,如圖6所示,由圖6可見波形參數和聲源深度對模型性能的影響小,魯棒性較好。

圖6 不同發射頻率下各個深度的測距精確度Fig.6 Ranging accuracy at different depths with different transmission frequencies

以線性調頻信號作為發射信號,分析模型的收斂速度和測量精度,經10000次迭代,對不同中心頻率(fc=500 Hz,1000 Hz,2000 Hz)、不同頻帶寬度(fband=100 Hz,300 Hz,500 Hz,700 Hz,900 Hz)和不同時間長度(T=200 ms,400 ms,600 ms,800 ms,1000 ms)的信號源測距結果進行對比。圖7~圖9給出了不同波形參數的信號的網絡訓練損失曲線,即訓練過程中隨著迭代次數增大,測試集上代價函數值的變化曲線。這里的代價函數值為測量的均方誤差,下降越快,說明模型收斂速度越快。由此,圖7表明模型的收斂速度和測量精度與頻帶寬度呈負相關,說明頻帶越寬所含信息越豐富,模型訓練需要的時間越長;圖8表明模型的收斂速度和測量精度與信號持續時間呈正相關,可見信號持續時間越長,呈現出的特征越明顯;圖9表明模型的收斂速度與中心頻率和測量精度呈負相關,可見淺海中低頻的信號特征更加集中,高頻信號的特征更加分散,表明模型更適用于遠程探測聲吶。

圖7 信號時長為1.0 s時各頻帶寬度對應的網絡訓練損失曲線Fig.7 Network training loss curve with different frequency bandwidth when the signal duration is 1.0 s

圖8 信號頻帶寬度為500 Hz時不同時間長度對應的網絡訓練損失曲線Fig.8 Network training loss curves with different time lengths when the frequency bandwidth is 500 Hz

圖9 信號頻帶寬度為500 Hz、時長為1.0 s時不同中心頻率對應的網絡訓練損失曲線Fig.9 Network training loss curves with different center frequencies when the frequency bandwidth is 500 Hz and the duration is 1.0 s

3 討論和展望

相對已有的深度學習聲源測距方法,本文提出的方法可提取信號的多域特征及對特征空間進行篩選,有利于產生緊湊性和泛化能力更強的模型;其次,改進了神經網絡模型的尋優算法和參數稀疏化策略,可加快收斂速度并抑制模型過擬合。然而,機器學習是一種監督學習方法,其測距精度以建立功能全、信息豐富的數據庫為代價。對于未知的海洋環境、時變的環境參數,需建立大量拷貝數據進行訓練并加大網絡的復雜度,同時需要已知的聲源波形參數,對先驗知識和計算資源有較大依賴性。目前已有相關研究[5]通過在真實海洋測試數據的基礎上擴充仿真數據,以及采用遷移學習方法來減小深度學習方法對先驗知識的依賴性。本文下一步工作將針對海洋測試中遠距離測距的實際應用場景,通過建立數據庫、優化神經網絡結構、提高網絡復雜度、增加輸出節點數或改進輸出層標簽形式、增強自學習能力,以提升測距范圍和測距分辨率。此外,可通過對數據庫進一步擴充與聲源目標識別相關的特征和標簽,從而在聲源被動測距的基礎上同時執行目標識別相關的任務,如估計目標材料、形狀、運動姿態等,以上需在后續工作中進一步探究。

4 結論

采用一種基于多域特征提取的深度學習方法來實現聲源測距,通過淺海環境仿真實例驗證了該方法的有效性,并分析了波形參數對測距性能和模型收斂速度的影響。本方法構建了聲信號在時、頻域的多維感知特征量,采用最大相關-最小冗余準則(mRMR)提取了聲源和水下聲場的關鍵信息,在傳統深度神經網絡的基礎上引入自適應矩估計(Adam)優化、L2正則化和Dropout正則化處理,提升模型的收斂速度和泛用性。結果表明:此方法在模型訓練過程中收斂速度較快,預測性能較穩定,在所定條件下測試集上聲源的綜合測距精確率可達到95%以上,能夠實現對聲源距離的有效估計。此外,對不同發射信號訓練效率的對比表明,算法性能對波形參數和聲源深度具有良好的魯棒性,模型收斂速度和測距精度對于帶寬較小、持續時間長的瞬態發射信號較高。訓練后的模型在單次測距任務中僅需執行毫秒級運算,可實現數據的實時處理。

本文提出采用的聲源測距的深度學習方法,緊密結合了和海洋環境、傳輸距離相關的聲源信號時頻域多維感知特征,測試集上聲源測距精度優良、可信。未來,建立并不斷更新充實功能齊全、信息豐富的各種典型海洋環境參數、傳輸距離及聲源信號多域特征大數據庫,進一步優化算法和自學習能力后,本方法可望實現實時準確的水下目標被動定位、跟蹤和分類識別,是今后深入研究的目標和方向。

猜你喜歡
特征信號模型
一半模型
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
重要模型『一線三等角』
完形填空二則
重尾非線性自回歸模型自加權M-估計的漸近分布
如何表達“特征”
不忠誠的四個特征
當代陜西(2019年10期)2019-06-03 10:12:04
基于FPGA的多功能信號發生器的設計
電子制作(2018年11期)2018-08-04 03:25:42
抓住特征巧觀察
3D打印中的模型分割與打包
主站蜘蛛池模板: 午夜国产理论| 波多野结衣中文字幕久久| 国产玖玖视频| 高清视频一区| 亚洲精品视频免费| 美女毛片在线| 国产成人无码综合亚洲日韩不卡| 99国产在线视频| 五月综合色婷婷| 亚洲 成人国产| 狠狠色狠狠色综合久久第一次| 在线看AV天堂| 亚洲a级在线观看| 成人91在线| 国产精品对白刺激| 日韩人妻精品一区| 国产18页| 高清免费毛片| 久久毛片免费基地| 国产正在播放| 国产真实乱子伦精品视手机观看| www.亚洲国产| 中文字幕人妻av一区二区| 中文字幕伦视频| 欧美在线网| 久久一日本道色综合久久| 911亚洲精品| 婷婷亚洲综合五月天在线| 久久www视频| 色屁屁一区二区三区视频国产| 免费看a毛片| 美女内射视频WWW网站午夜| 怡红院美国分院一区二区| 亚洲av成人无码网站在线观看| 国产全黄a一级毛片| 欧洲亚洲欧美国产日本高清| 久久精品娱乐亚洲领先| 在线观看欧美国产| 91色在线视频| 在线中文字幕网| 一级毛片网| 日本成人在线不卡视频| 伊人久久青草青青综合| 麻豆国产原创视频在线播放| 欧美色丁香| 永久天堂网Av| 日本国产精品一区久久久| 五月天丁香婷婷综合久久| 亚洲品质国产精品无码| 天天婬欲婬香婬色婬视频播放| 91麻豆国产视频| 美女黄网十八禁免费看| 91麻豆久久久| 黄片在线永久| 中文字幕欧美日韩高清| 欧美精品一区在线看| 亚洲一区无码在线| 国内老司机精品视频在线播出| 精品無碼一區在線觀看 | jijzzizz老师出水喷水喷出| 综合色区亚洲熟妇在线| 丁香亚洲综合五月天婷婷| 91精品国产福利| 久久久久国产精品嫩草影院| 国产亚洲男人的天堂在线观看| 久久青青草原亚洲av无码| 成人精品在线观看| 国产嫩草在线观看| 欧美日韩中文国产va另类| 国产乱人伦偷精品视频AAA| 最新亚洲人成无码网站欣赏网| 中日韩一区二区三区中文免费视频| 浮力影院国产第一页| 激情午夜婷婷| 精品免费在线视频| 伊人91在线| 国产区网址| 好吊日免费视频| 久久青草免费91线频观看不卡| 九九九国产| 高清免费毛片| 精品成人免费自拍视频|