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

無監督殘差網絡的地震數據重構方法

2022-08-02 07:36:02孟宏宇楊華臣張建中
石油地球物理勘探 2022年4期
關鍵詞:方法模型

孟宏宇 楊華臣 張建中*③

(①海底科學與探測技術教育部重點實驗室,266100;②中國海洋大學海洋地球科學學院,山東青島 266100;③青島海洋科學與技術試點國家實驗室海洋礦產資源評價與探測技術功能實驗室,山東青島 266100)

0 引言

在實際地震資料的采集過程中,受現場條件的限制,會遇到炮點和檢波點布設困難的情況,或者因設備損壞、故障等因素產生壞道,導致采集的地震數據在空間上不連續,出現地震道的缺失,給后續地震資料的進一步處理帶來困難。因此需使用一定的處理方法,對缺失的地震道數據進行重構。

地震數據的重構技術可分為傳統方法和基于機器學習的人工智能方法。傳統的地震數據重構方法可分為以下四種。

第一種是基于波動方程的地震數據重構方法[1],該方法使用偏移算子對地震數據進行重構,可用于處理隨機缺失和規則缺失的地震數據。但是,該方法需要一個近似的地下速度模型,并且計算量較大。

第二種是基于預測濾波的地震數據重構方法,該方法是通過設計預測誤差的濾波器對存在空間假頻的地震數據進行重構,通常用于存在規則道缺失的地震數據。Spitz[2]提出了經典的f-x域數據重構方法,能重構規則缺失的地震道記錄。Porsani[3]進一步提出半步f-x域地震道重構方法,提高了數據重構的速度和精度。在此基礎上,宜明理等[4]研究了f-k域抗假頻地震數據重構方法,張軍華等[5]采用并行算法在f-x域實現了三維地震數據重構。

第三種是基于多域變換的地震數據重構方法,是目前地震數據重構的主流方法。該類方法首先使用特定的數學變換,如傅里葉變換、拉東變換、小波變換等,將地震數據變換到其他域進行數據重構,然后對處理后的數據進行反變換,實現地震數據的重構。Duijndam等[5]首先提出了基于傅里葉變換的地震數據重構方法,使用最小二乘法估計傅里葉分量,對沿一個空間方向的不規則道缺失地震數據進行重構。Herrmann等[6]研發了一種基于Curvelet變換的地震數據重構方法,對存在多道缺失的地震數據進行重構。Fomel等[7]提出基于一種Seislet變換的反假頻迭代地震數據重構方法,利用壓縮感知Bregman迭代算法,對缺失的地震道重構。Ibrahim等[8]提出了一種基于拉東域的地震數據重構方法,利用全為零值的缺失地震道與未缺失地震道在拉東域中曲率范圍不同這一特點,通過去除缺失地震道對應的拉東系數,再利用拉東逆變換變回t-x域實現重構不規則缺失地震數據。Gan等[9]提出了在POCS算法框架下應用Seislet變換重構地震數據的方法,可重構非規則道缺失的地震數據。

第四種是基于矩陣降秩的地震數據重構方法[10],該類方法利用當數據存在缺失時,數據在頻率切片的Hankel矩陣的秩會升高這一原理,通過降低Hankel矩陣的秩或對數據降維重構數據。Zhang等[11]提出一種降秩與稀疏促進混合的數學模型,并將該模型用于三維地震數據的重構。Oropeza等[12]提出基于隨機奇異值分解方法重構地震數據,解決了原始截斷奇異值分解方法計算量較大的問題,提高了計算效率,低維地震數據的重構效果顯著。

近年來,隨著人工智能技術的發展,深度學習(Deep Learning,DL)方法越來越多地應用于地震數據重構。Wang等[13]利用存在道缺失的地震數據和完整的地震數據預先訓練殘差網絡,然后使用訓練后的網絡對道缺失的地震數據進行重構。該方法只適用于重構規則道缺失的地震數據,對不規則道缺失的地震數據的重構效果較差。鄭浩等[14]提出一種基于DL卷積神經網絡的智能化地震數據重構方法,該方法首先利用大量地震數據訓練搭建的卷積自編碼器模型,然后使用訓練后的卷積自編碼器模型重構缺失地震道的數據重構。Alwon[15]和Oli-veira等[16]提出使用生成對抗網絡重構道缺失的地震數據,但面對地質構造比較復雜的地震數據,會出現非地質偽影。Kaur等[17]進一步提出了使用生成對抗網絡重構三維地震數據的方法。傳統的有監督DL地震數據重構方法主要包括三個步驟:首先,構建訓練集,訓練集一般由道缺失的地震數據(輸入)和對應的完整地震數據(期望輸出,標簽)構成;然后,使用構建的訓練集訓練網絡模型,學習數據的特征;最后,使用訓練后的網絡模型重構道缺失的地震數據。一般利用DL解決地震數據重構問題的方法需使用大量合適的訓練集訓練網絡,并且訓練集的質量直接影響重構效果[18]。但是,從野外地震數據很難獲得訓練所用的標簽。為此,本文提出一種基于無監督DL地震數據重構方法。該方法可直接對存在道缺失的不完整地震數據進行重構,而無需構建訓練集預先訓練網絡,避免了有監督DL數據重構方法建立包含完整地震數據標簽的訓練集問題。

殘差網絡由He等[19]首先提出,其內部的殘差單元使用了跳躍連接,提升了信息傳遞路徑的數量、簡化了網絡的學習目標及難度、提升了網絡的準確度。同時,殘差網絡通過訓練使網絡中多余層的殘差為0,使得這些層變為恒等映射,解決了層數加深導致的網絡退化問題,提升了網絡性能。

與傳統的DL方法不同的是,本文提出的基于無監督殘差網絡DL的地震數據重構方法無需使用完整的地震數據作為標簽預先訓練殘差網絡,而是用隨機數據作為殘差網絡的輸入,以缺失道的原始地震記錄作為殘差網絡的期望輸出,通過反向傳播殘差網絡的預測輸出與期望輸出之間的誤差,迭代優化殘差網絡參數,使網絡輸出與原始地震記錄的差別最小,從而獲得參數最優的殘差網絡,再用該網絡計算缺失的地震道記錄。在網絡參數優化過程中,利用卷積的局部和平移不變性質,用卷積濾波器學習多尺度下地震數據鄰域之間的相似特征,并在網絡輸出中呈現。因無需使用訓練集,因此獲得的殘差網絡具有很強的泛化能力,適于不同特征地震數據重構。使用Marmousi模型模擬地震數據和實測海洋拖纜數據進行測試,結果驗證了本文方法的有效性。

1 深度殘差網絡地震數據重構原理

1.1 深度殘差網絡

廣泛應用的卷積神經網絡(Convolutional Neural Networks,CNN)主要由卷積層、池化層和全連接層組成。其中卷積層用于局部感知、特征提取;池化層用于合并語義相似的特征。通常,卷積層與池化層相互連接配合,組成多個卷積—池化組逐層提取特征,最后通過全連接層完成數據的特征提取、綜合與數據壓縮。在訓練階段,數據通過自動特征提取層和分類層向前傳播,得到CNN的預測輸出,然后通過損失函數計算總體誤差,通過反向傳播更新網絡參數,使總體誤差最小化。

提高CNN模型的深度可增加網絡中訓練參數的種類和數量,從而提高網絡模型的表達能力和特征提取能力。但是,當網絡層數超過一定數目后,會導致模型出現梯度爆炸和梯度消失等現象,網絡的訓練誤差升高、性能下降。為了克服這一缺陷,在網絡模型中引入跳躍連接結構,通過訓練使多余層的殘差為0,從而使這些層變為恒等映射,解決了網絡層數過大時引起的梯度消失和梯度爆炸等問題。

圖1為普通CNN的卷積單元與殘差網絡單元對比。如圖1a所示,假設某個神經網絡單元的輸入是x,實際輸出為F(x),期望輸出為H(x)。希望網絡單元學習到一個恒等映射函數H(x)=x,但隨著網絡層數的增加,這個恒等映射會變得越來越難以擬合,增大了直接學習的難度。為解決這一問題,殘差網絡利用跳躍連接將殘差單元的輸出變為F(x)+x。殘差單元不直接學習目標映射,而是轉換為學習一個殘差函數F(x)=H(x)-x,即網絡的學習內容變為殘差[20]。殘差網絡不僅解決了由于網絡層數增加造成的退化問題,還簡化了網絡的學習目標、降低了深度神經網絡的訓練難度[21]。

圖1 不同單元結構對比

深度殘差網絡(圖1b)由許多基本的殘差單元(Residual Block)組成,經典的兩層殘差單元包含兩個卷積層,每層卷積操作后使用ReLU激活函數激活后輸入到下一層。每個殘差單元的一般形式為

yt=xt+F(xt;Wt)

(1)

xt+1=h(yt)

(2)

式中:xt和xt+1分別是第t個殘差單元的輸入與輸出;Wt={Wt,k|1≤k≤K}是與第t個殘差單元相關聯的一組權重和偏差,其中K是單個殘差單元中的總層數(圖1b中K=2);k是當前層數;h(·)是激活函數。圖1b中常規基本殘差單元使用的是ReLU激活函數。

如果h(·)是一個恒等映射,即xt+1=yt,則可得到

xt+1=xt+F(xt;Wt)

(3)

將殘差網絡的輸入與輸出代入式(3),可得到殘差函數的一般形式

(4)

式中T為殘差單元的總數。式(4)具有良好的反向傳播特性,即使權重很小,層的梯度也不會消失。

1.2 地震數據重構方法

在地震數據重構過程中,設X∈Si×j為理想的完整地震數據(S代表單炮地震數據,i、j分別代表采樣點數和道數);X0∈Si×j為含缺失地震道的已知地震數據。是利用已知地震數據X0對缺失的地震道進行重構,得到完整的數據X。X0與X的關系可表示為

X0=X°s

(5)

式中:“°”表示哈達瑪積(Hadamard Product);s∈(0,1)i×j是采樣算子,是由0和1組成的矩陣,存在道缺失的位置為0,其余為1。式(5)的目的是區分出原始地震數據中已有地震道與需要重構的地震道。由X0求X是一個不適定的反演問題,為求得一個合適的解,必須加入某種先驗信息對反演進行約束。因此,將該反演問題轉換為求下列目標函數的最優解x*

(6)

式中:Q(·)是數據保真項,這里取Q(·)為二階范數;R(·)是一個先驗正則化函數;θ表示網絡參數;N為輸入網絡的原始隨機數據;fθ(·)代表網絡參數為θ的CNN模型,表示輸入與輸出間的非線性映射函數。本文中去掉正則化函數R(·),改用神經網絡參數來捕獲數據中的內在屬性。此時,上述反演問題便轉換為求解神經網絡的最優參數θ*

(7)

將式(7)改寫為如下最小化問題

(8)

網絡參數θ由神經網絡各層的深度、大小、核的個數、每個核中的元素以及每個核對應的偏差構成,其初始值隨機確定,使用優化算法迭代更新θ,通過令網絡輸出fθ(N)與s相乘,使X′與X0之間的誤差達到極小,獲得網絡最優參數θ*,從而得到完整的數據X

X=fθ*(N)

(9)

這樣構建的殘差網絡,建立了從N到X0的之間的復雜映射關系。在網絡參數優化過程中,通過卷積運算將地震數據相鄰道之間的相似特征施加到輸出的數據上。用該網絡預測的X′與X0最匹配,從而能夠重構缺失的地震道。

圖2、圖3對本文無監督學習方法與有監督學習方法進行了對比。一般的有監督學習方法地震數據重構流程如圖2所示。將完整的地震數據去掉一定的地震道后作為神經網絡的輸入,相應的完整地震數據作為標簽,即期望輸出。使用建立的大量訓練數據集,對神經網絡預先訓練,學習訓練數據集的特征,然后使用訓練后的網絡模型對存在道缺失的地震數據進行重構。

圖2 有監督學習方法數據重構流程

本文無監督學習方法地震數據重構流程如圖3所示。利用隱藏在已知地震道中的先驗特征進行重構缺失道數據。其中,網絡的輸入為隨機數據N,預測輸出為X′,期望輸出為X0。隨機初始化網絡參數θ,網絡從輸入N得到殘差網絡的輸出fθ(N)。將fθ(N)與s作用后得到X′,對X0與X′之間的均方誤差(MSE)進行反傳,優化網絡參數,使該均方誤差最小化,得到網絡的最優參數θ*。

圖3 無監督學習方法數據重構流程

1.3 深度殘差網絡架構

本文使用的殘差網絡整體框架如圖4所示,網絡的第二層到最后一層之間由多個殘差單元構成,并且使用了LeakyReLU激活函數,加入了批量歸一化操作。殘差單元的具體結構如圖5所示,每個殘差單元由兩個卷積層(Conv)、兩個批量歸一化層(BN)和一個LeakyReLU激活函數組成。殘差網絡的主體結構中,第一個卷積層的尺寸為3×3×32×64(包含64個卷積核為3×3的卷積濾波器,通道數為32),最后一個卷積層的尺寸為3×3×64×1(包含1個卷積核為3×3的卷積濾波器,通道數為64)。本文使用的殘差網絡一共包含8個殘差單元,每個殘差單元中的卷積層尺寸均為3×3×1×64(包含64個卷積核為3×3的卷積濾波器,通道數為1)。網絡輸入尺寸為512×512的隨機數據,輸出為重構的地震數據。

圖4 深度殘差網絡架構示意圖

圖5 殘差單元架構(Residual Block)示意圖

該模型在64位工作站上進行訓練,配置是NVDIA RTX 3090 GPU和Intel Core i9-10900k CPU,內存為64G,使用的是Python v3.6.1、Cuda11.1.0和Pytorch框架。

在網絡模型的訓練過程中,梯度消失和梯度爆炸會嚴重影響網絡模型的收斂。為了解決該問題,在原始殘差單元中加入了BN層。通過執行BN處理規范網絡中各層的數據,使網絡數據的均值為0、方差為1。上述處理可有效加快神經網絡的學習速度,在解決梯度彌散問題的同時增加網絡整體的穩定性和泛化能力[21]。

在殘差網絡中使用了性能更優異的LeakyReLU激活函數(圖6)。該函數的輸出對于輸入的負值具有很小的坡度,因此導數不會取到零值。這樣可減少靜默神經元的出現,避免了使用ReLU函數進入負區間后神經元不再繼續學習的問題。使用LeakyReLU函數對上一層節點的輸出做非線性變換,有助于殘差網絡更好地學習輸入與輸出間的復雜映射關系[22-23]。

圖6 LeakyReLU函數

LeakyReLU函數的具體公式如下

(10)

式中:z是LeakyReLU函數的輸入;L(z)表示LeakyReLU函數的輸出。

本文殘差網絡使用自適應性矩估計(Adam)優化算法。與常規優化算法相比,Adam優化算法更高效、更準確,對內存需求更少,且網絡參數的更新不受梯度變換影響,更適用于數據與參數較多的場景[24]。

網絡的收斂性與設置的學習率密切相關,學習率不宜過大。為了使殘差網絡在較快收斂同時使網絡的損失達到最小,本文采用0.001的學習率重構地震數據[25]。

2 模型數據測試

為驗證本文方法的地震道重構效果,利用理論模擬數據進行了測試。使用長10km、厚5km 的Marmousi模型。檢波器均勻分布于地表,檢波器間距為20m;炮點位于地表距離首個檢波器5km處;用20m×20m網格離散該模型,模型被離散成1000×500個矩形單元;震源采用主頻為8Hz雷克子波;時間采樣間隔為0.002s;記錄時長為4s。用波動方程有限差分法模擬地震波[26]。

首先從模擬地震記錄中按一定比例去除部分地震道,就得到存在缺失地震道的地震數據(圖7);然后用本文地震道重構方法重構整個地震數據。重構地震數據與原始地震數據有一定誤差,相當于對原始地震數據引入了噪聲。借用“信噪比(SNR)”定量說明對地震數據的重構效果。SNR計算式為

(11)

SNR越大,表示重構的地震數據與原始地震數據之間的差異越小,即重構效果越好。

本文設計的損失函數如下

(12)

2.1 網絡優化算法比較

本文的深度殘差網絡選用了Adam網絡優化算法,這里與Rprop、RMSprop算法進行對比測試,以說明選用Adam算法的優勢。對模擬地震數據(圖7a)隨機去除其中60%的地震道(圖7c),然后使用上述三個網絡優化算法對缺失的地震數據進行重構。設置學習率為0.001,迭代次數k為30000。使用不同優化算法的網絡重構數據的信噪比和MSE損失函數隨迭代次數的變化如圖8a、圖8b所示。由圖8a可見,Adam優化算法重構結果信噪比曲線達到穩定時,得到的重構數據信噪比最高,而且相比于其他算法,Adam優化算法重構地震數據的信噪比隨迭代次數上升得更快。由圖8b可見:使用Rprop優化算法的MSE損失曲線下降很慢;RMSprop的MSE損失曲線波動劇烈;使用Adam優化算法的MSE損失曲線更穩定、更快速收斂。

圖7 模擬地震數據

圖8 不同優化算法重構數據的SNR(a)與MSE(b)對比

因此,本文使用的Adam優化算法性能更好。

2.2 網絡結構比較

對本文使用的殘差網絡與不含跳躍連接結構的普通網絡模型進行對比試驗。分別使用殘差網絡和普通網絡對圖7c進行重構。兩種網絡模型的SNR和MSE損失函數隨迭代次數的變化如圖9、圖10所示。由圖可見,隨著迭代次數的增加,兩種網絡模型的MSE損失函數逐漸降低,重構地震數據信噪比逐漸增加,但殘差網絡相比于普通網絡模型的MSE損失更低,重構數據SNR更高。

圖9 不同網絡重構SNR對比

圖10 不同網絡MSE對比

2.3 重構地震數據隨迭代次數變化

為了直觀地展示使用本文方法對道缺失的地震數據重構的迭代過程,將圖7c作為原始地震數據進行完整數據重構。圖11所示是不同迭代次數的重構結果。由圖可見:因網絡輸入是隨機信號,迭代初期的輸出表現為隨機噪聲(圖11a);隨著迭代次數增加,能量強的地震波先逐漸重構出來(圖11b),然后,能量弱的地震波也被逐漸重構。當迭代次數達到170次時,地震資料的總體面貌開始被恢復,經過2000次迭代基本上恢復了地震資料。可見,以隨機信號作為輸入,通過對網絡輸出與原始含缺失道的地震數據之間的誤差進行反饋,不斷優化網絡模型,達到重構完整地震數據的目的,避免了常用監督學習網絡建立訓練集的繁重環節。

圖11 本文方法不同迭代次數((a)~(d)對應10、90、170、2000次迭代)的重構結果

2.4 模擬資料重構效果分析

對圖7a數據分別規則和隨機去除其中20%、40%和60%的地震道,使用本文方法對不同比例道缺失的地震數據重構。設置學習率為0.001、k為50000。重構地震數據的信噪比隨迭代次數的變化如圖12所示。由圖可見:道缺失數據的原始信噪比均低于10dB;經過1000次迭代后重構數據的SNR迅速提高,在20000次迭代后達到較高的重構質量并趨于穩定,說明本文方法在迭代次數較低時重構數據的SNR增長迅速,達到一定迭代次數后重構數據的SNR到達高位并保持穩定;數據缺失比例越大,重構數據的SNR越低,缺失的地震道越少重構效果越好;在數據缺失相同比例的情況下,對規則缺失數據的重構SNR要大于隨機缺失數據的重構信噪比,這是因為隨機缺失數據容易出現連續多道缺失的情況,使重構誤差較大。

圖12 本文方法對不同比例的規則缺失(a)、隨機缺失(b)數據重構SNR的對比

快速凸集投影軟閾值(FPOCS-Soft)方法是一種較新的地震道插值方法,可直接重構道缺失地震數據,具有較高計算精度和穩定性[27]。為了說明本文方法的有效性和精度,分別使用本文方法和FPOCS-Soft方法對隨機缺失40%地震道的數據進行重構,并對兩種重構結果進行對比。

圖13是隨機去除40%地震道的地震資料(圖7b)的重構結果。圖13a和圖13b分別為本文方法重構的完整地震記錄及誤差;圖13c和圖13d分別為使用FPOCS-Soft方法重構的完整地震記錄及誤差。由圖可見:本文方法在連續缺失道數較少時重構數據的誤差幾乎為0,在連續缺失道數較多的情況時重構數據的誤差會有所上升,但重構效果仍然遠好于FPOCS-Soft方法。FPOCS-Soft方法對連續缺失道的重構誤差更明顯,并且在無地震數據的區域還會出現一些假的地震數據。使用本文方法和FPOCS-Soft方法重構的地震數據的SNR分別為36.41dB和18.6dB,均方誤差分別為2.29×10-6、1.01×10-4。

圖13 隨機缺失40%地震數據兩種方法重構結果及誤差

圖14為對圖7b中第162道數據的重構結果對比。由圖可見,本文方法重構精度明顯高于FPOCS-Soft方法。另外,由于隨機缺失數據較規則缺失數據存在連續多道數據缺失,當連續缺失的地震道較多時,本文方法的重構精度會有所下降,但重構效果仍然明顯優于FPOCS-Soft方法。

圖14 兩種方法第162道模擬地震數據重構結果對比

3 實際資料應用

為了檢驗本文方法對實際數據的應用效果,將本文方法應用于南黃海海洋拖纜實測資料的單炮數據。對該地震數據分別進行規則和隨機去除40%的地震道(圖15),分別用本文方法和FPOCS-Soft方法對缺失地震道重構。設置學習率為0.001,迭代次數為20000。

圖16是對圖15c分別用本文方法和FPOCS-Soft方法重構結果及誤差分布,重構數據的信噪比分別為16.83dB和11.23dB,均方誤差分別為4.32×10-4和1.30×10-3。可以看出,本文方法可以很好地重構缺失的拖纜數據,同相軸變得更連續,重構的結果與原始拖纜數據(圖15a)基本一致,在沒有數據的位置不會出現假數據。而FPOCS-Soft方法對連續同相軸的重構結果存在較明顯的誤差,并且在重構數據中引入了新噪聲。

圖17是對圖15b和圖15c的第43道數據的重構結果及誤差對比。可見不論對規則缺失還是隨機缺失地震道的數據,與FPOCS-Soft方法結果相比,本文方法結果與原始地震數據(圖15a)基本一致,對缺失道數據的重構精度更高。對比規則缺失與隨機缺失地震數據的重構結果,可見因隨機缺失資料出現多道連續缺失的情況,本文方法的重構精度低于規則缺失數據,但仍然高于FPOCS-Soft方法對規則缺失數據的重構精度,說明了本文方法的良好應用效果和對不同缺失數據更強的適應性。

4 結論

本文提出了一種基于無監督殘差網絡的缺失地震數據重構方法。該方法以隨機數據作為網絡輸入,通過對網絡輸出與含缺失道的地震數據之間的誤差的反向傳播,通過迭代優化網絡參數;在迭代過程中,利用卷積的局部和平移不變性質,使用帶有跳躍連接的殘差網絡結構,學習多尺度地震數據鄰域之間的相似特征,并把這些先驗特征反映在網絡輸出中;使用性能優良的Adam優化算法,進一步提高了構建的深度殘差網絡模型的性能。本文方法只需使用地震道缺失的數據和優化的殘差網絡便可重構完整的地震數據,無需由訓練集訓練殘差網絡模型,避免了常規基于深度學習方法的構建大型訓練數據集和網絡模型泛化能力的問題,實用性更強。

利用模擬地震數據和實測海洋拖纜地震數據對本文方法進行了測試,并與效果較好的FPOCS-Soft地震道插值方法效果進行了比較。結果表明,不論對規則缺失還是隨機缺失地震道的地震數據,本文方法的重構精度都高于FPOCS-Soft方法,說明了本文方法的可行性和良好效果。對于存在連續缺失多道的地震數據,本文方法重構地震道的精度會降低,但仍然高于FPOCS-Soft方法。

猜你喜歡
方法模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
學習方法
3D打印中的模型分割與打包
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
主站蜘蛛池模板: 欧美国产菊爆免费观看| yy6080理论大片一级久久| 国产精品99r8在线观看 | 国产色伊人| 亚洲天堂视频在线观看| 亚洲国产天堂在线观看| 无码综合天天久久综合网| 一级毛片视频免费| 日韩a级片视频| 欧美成a人片在线观看| 日韩国产另类| 亚洲久悠悠色悠在线播放| 天天色综合4| 激情六月丁香婷婷| 亚洲成人在线网| 日韩高清一区 | 亚洲成人一区二区三区| 久久精品中文字幕免费| 亚洲中文字幕久久精品无码一区 | 国产网站免费观看| 久草视频中文| 在线看免费无码av天堂的| 毛片在线播放a| 美女无遮挡免费网站| 免费xxxxx在线观看网站| 无码AV高清毛片中国一级毛片| 亚洲一区二区无码视频| 精品久久久久久中文字幕女| 国产高清在线观看91精品| 99在线小视频| 乱人伦99久久| 狼友av永久网站免费观看| 在线人成精品免费视频| 日本91在线| 免费人欧美成又黄又爽的视频| 青草视频在线观看国产| 国产精品美乳| 婷婷五月在线视频| 熟妇丰满人妻av无码区| 国产99久久亚洲综合精品西瓜tv| 日韩在线第三页| 大学生久久香蕉国产线观看| 伊人无码视屏| 久久精品只有这里有| 中文字幕人成人乱码亚洲电影| 日韩成人免费网站| 日韩av无码精品专区| 亚洲最大在线观看| 日韩小视频网站hq| 天堂在线视频精品| 国产免费羞羞视频| 欧美劲爆第一页| 国产在线观看一区精品| 小13箩利洗澡无码视频免费网站| 欧美另类视频一区二区三区| 爆操波多野结衣| 亚洲品质国产精品无码| 亚洲人免费视频| 国产精品极品美女自在线| 综合成人国产| 国产亚洲精久久久久久久91| 青青青国产免费线在| 国产一级无码不卡视频| 久草性视频| 成人伊人色一区二区三区| 亚洲天堂高清| 久久精品嫩草研究院| 欧美不卡在线视频| 18黑白丝水手服自慰喷水网站| 老司机久久99久久精品播放| 在线欧美日韩国产| 国产精品密蕾丝视频| 无码区日韩专区免费系列 | 一本综合久久| 浮力影院国产第一页| 最新国语自产精品视频在| 四虎永久在线精品国产免费| 日韩精品无码一级毛片免费| 伊人丁香五月天久久综合| 国产综合另类小说色区色噜噜 | 中文字幕天无码久久精品视频免费 | 亚洲天堂首页|