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

自適應時窗多道相關的疊后地震波阻抗反演方法

2024-08-22 00:00:00彭真許輝群張偉
石油地球物理勘探 2024年4期
關鍵詞:方法

摘要:目前提高反演橫向連續性的主要手段仍是基于模型參數先驗信息假設的正則化方法,當先驗信息與實際情況矛盾時,正則化方法會出現問題。為此,基于多道反演理論和方法,提出一種自適應時窗多道相關的疊后地震波阻抗反演方法。首先,利用主動學習的思想,通過查詢原始地震數據篩選地震數據關鍵特征的位置;其次,根據上述位置將地震數據劃分為不同尺度的時窗,有助于斷層信息保護以及異常體邊界識別;再次,根據點乘計算自適應時窗內地震道之間的相關性;最后,通過“加權求和”地震道得到重構地震道——反演目標泛函的約束條件。所提方法借助主動學習以及地震數據局部相關的思想,較好地反映了自適應時窗內的多道間的相關性,并且可并行化處理。模型測試和實際數據測試表明,所提方法反演結果的橫向連續性較好,且可以兼顧特殊構造信息,但計算效率較低,在淺、深層速度差異較大的區域存在一定局限性,適用于具有特殊構造的狹小地區。

關鍵詞:自適應時窗,多道相關,橫向連續性,相關性,地震波阻抗反演

中圖分類號:P631文獻標識碼:A DOI:10.13810/j.cnki.issn.1000-7210.2024.04.020

Post?stack seismic impedance inversion method with multi?channel correlation under adaptive window

PENG Zhen1,XU Huiqun1,ZHANG Wei2

(1.College of Geophysics and Petroleum Resources,Yangtze University,Wuhan,Hubei 430100,China;2.New Resources Geophysical Exploration Division,BGP Inc.,CNPC,Zhuozhou,Hebei 072751,China)

Abstract:At present,the main way to improve the lateral continuity of inversion is still the regularization method based on the assumption of prior information of model parameters.When the prior information contra?dicts the actual situation,the regularization method will have problems.Therefore,based on the theory and method of multi?channel inversion,a post?stack seismic impedance inversion method with multi?channel correla?tion under an adaptive window is proposed.Firstly,using the idea of active learning,the position of key fea?tures of seismic data is screened by querying the original seismic data.Secondly,the seismic data are divided into time windows of different scales according to the above locations,which is helpful for fault information pro?tection and abnormal body boundary identification.Thirdly,the correlation between seismic traces in the adap?tive window is calculated according to the dot multiplication.Finally,the“weighted sum”seismic trace is used to obtain the reconstructed seismic trace,which is the constraint condition of the inversion objective function.With the help of the idea of active learning and local correlation of seismic data,the proposed method can better reflect the correlation between multiple channels in the adaptive window and can be processed in pa?rallel.The model test and the field data test show that the proposed method has good lateral continuity of inversion results,and can take into account the special structural information.However,the computational efficiency is low,and there are some limitations in the shallow and deep areas with large speed differences,making the proposed method suitable for narrow areas with special structures.

Keywords:adaptive window,multi?channel correlation,lateral continuity,correlation,seismic impedance inversion

彭真,許輝群,張偉.自適應時窗多道相關的疊后地震波阻抗反演方法[J].石油地球物理勘探,2024,59(4):828-836.

PENG Zhen,XU Huiqun,ZHANG Wei.Post-stack seismic impedance inversion methodwith multi-channel correlation under adaptive window[J].Oil Geophysical Prospecting,2024,59(4):828-836.

0引言

地震波阻抗反演是獲取地下介質信息的有效途徑,反演結果可以反映地下介質的巖性變化,是儲層預測的一項關鍵技術[1-2]。地震波阻抗具有高分辨率、地質信息豐富的特點,疊后地震波阻抗反演在儲層預測中具有重要意義[3-5]。然而,由于地震數據的頻帶限制以及數據噪聲等影響,地震波阻抗反演問題是一個不適定問題。

為了獲取較穩定和精確的反演結果,可以使用正則化方法約束反演過程,從而獲取高精度反演結果[6]。目前有平滑約束和塊約束兩種典型的正則化地震反演方法[7-9]。平滑約束有Tikhonov[10]提出的L2范數約束反演參數的方法,通過低頻初始模型得到平滑反演結果。然而平滑約束會模糊邊界和斷層信息,并不能有效提高反演結果的分辨率。為此,王寧等[11]提出數據驅動的塊排列正則化方法,利用疊后地震數據高信噪比的特點,通過提取一個塊排列的矩陣記錄地下構造信息,利用地震數據局部相似性記錄相鄰采樣點的位置,通過構造地震道的塊排列正則項有效提高了反演結果的分辨率。然而,上述方法是逐道進行的,忽略了地層構造信息的空間連續性,使單道反演結果的橫向連續性較差[12-14]。

現階段地震波阻抗反演問題可以轉化為常規的單道反演,單道反演的計算效率高、穩定性強[15-16]。由于單道反演只考慮當前道的地震響應,缺乏相鄰道的橫向約束,從而不易控制反演結果的橫向連續性[17]。因此,有人嘗試將單道地震波阻抗反演轉變為多道反演。Gholami[18]利用全變分約束多道非線性波阻抗反演,反演結果具有塊狀特征,增強了空間連續性。Hamid等[5]將所有地震道連成一道數據,并在水平方向平滑約束相鄰道反演結果;但是當真實情況與假設前提不一致時,不能得到準確的反演結果。隨后,人們發展了數據驅動的多道地震反演方法,如Hamid等[19]利用觀測數據計算地震同相軸的局部傾角,并將傾角信息和測井數據引入地震反演,提高了反演結果的橫向連續性。Huang等[20]根據地震數據估算斜率屬性,并引入全變分正則化反演,獲得了高分辨率反演結果。印興耀等[21]提出了地震數據互相關驅動的多道反演方法,通過計算地震數據的相關性描述地震結構特征,并作為地震反演的一個約束項,該方法極大地提高了反演結果的橫向連續性。

目前提高反演橫向連續性的主要手段仍是基于模型參數先驗信息假設的正則化方法,當先驗信息與實際情況矛盾時,正則化方法會出現問題[22]。為此,文中基于多道反演理論和方法,提出一種自適應時窗多道相關的疊后地震波阻抗反演方法。首先,利用主動學習的思想[23],通過查詢原始地震數據篩選地震數據關鍵特征的位置;其次,根據上述位置將地震數據劃分為不同尺度的時窗,有助于斷層信息保護以及異常體邊界識別;再次,根據點乘計算自適應時窗內地震道之間的相關性;最后,通過“加權求和”地震道得到重構地震道——反演目標泛函的約束條件。所提方法借助主動學習以及地震數據局部相關的思想,較好地反映了自適應時窗內的多道間的相關性,并且可并行化處理。文中利用模型和實際數據測試了所提方法的效果。

1方法原理

所提方法在計算地震道相關性時根據數據特征變化確定參與計算的地震道,而不是人為選取地震道。鑒于此,采用主動學習(Active Learning)[23]的整體思想自動選取地震道,主要利用機器學習算法查詢地震數據并篩選地震數據關鍵特征的位置,通過上述位置劃分地震數據,得到一系列不同道數的地震數據,由此計算地震數據的相關性更符合地質變化規律,從而獲得多道時窗自適應效果。所提方法避免了固定時窗導致的斷層以及異常體邊界信息的破壞。

由于地震數據可以視為一種序列信息,并且地震道之間具有一定相關性,因此可以計算時窗內的每一道與所有道的相關性,從而得到每一道占所有道的權重,通過不同權重的“加權求和”得到重構地震道。重構地震道包含了其他道的信息,即權重越小的道包含的信息越少,權重越大的道包含的信息越多。

1.1自適應時窗的創建

為更好地保護斷層、異常體等邊界信息,采用查詢算法得到地震數據的關鍵特征位置,并利用該位置劃分地震數據,從而識別關鍵地質構造并實現自適應時窗多道相關。

首先,建立機器學習模型。利用已知地震數據和測井波阻抗數據預訓練徑向基函數(RBF)神經網絡得到一個初始模型。其次,加入查詢算法。采用機器學習的最小置信度算法查詢地震數據,即

式中:y(?)為對于給定的地震數據樣本x,通過模型預測得到概率最大的波阻抗數據;pθ(?)為在模型參數為θ時計算的概率值;arg max(?)為最大概率值的位置;arg min(?)為最小概率值的位置;xL(?)C為x中最小置信度值位置的地震數據。

通過上述算法可以找到關鍵地震特征的位置,利用該位置即可將原始地震數據劃分為不同尺度的時窗(圖1)。

1.2地震數據相關性計算

計算任一時窗地震數據的相關性的步驟分述如下。

(1)將時窗內的地震數據排序,得到seis1,seis2,…,seisj,…,seisn,j=1,2,…,n,n為總CDP數。

(2)設置單位矩陣WQ、W K、W V作為求取權重的基礎矩陣,然后將時窗內的每道數據seisj分別與WQ、W K、W V相乘得到qj、kj、vj,即通過

qj=WQ seisj(2)

kj=WK seisj(3)

vj=W Vseisj(4)

初始化矩陣。WQ、W K、W V是單位矩陣,主要利用WQ、W K求取時窗內每一道與其他道的權重系數,W V主要用于“加權求和”。設置三個單位矩陣可以確保并行化計算。

(3)求解每一道與其他道的相關性,即

qj(T)?ki

式中d為qj(T)?ki的維度,目的是使qj(T)?ki的結果滿足期望為0、方差為1的正態分布,類似歸一化操作。通過qj(T)查詢j道與第i(i≠j)道的相關性,即求取每一道與其他道的相關性并Softmax歸一化,可以得到權重矩陣

A=A n1(A 21)A n2(A 22)A nn(A 2n)」|(6)

該權重矩陣元素的數值范圍為0~1,以進行加權求和。

(4)將步驟(2)得到的vj轉換為

V=[v1,v2,…,v n](7)

將A與V相乘重構一個具有相鄰地震道信息的多道地震數據,即

A seis=AVT(8)

通過多道信息約束可保證地震波阻抗反演結果的橫向連續性。

通過上述步驟即可實現自適應時窗的多道相關性計算。

1.3自適應時窗多道相關的疊后地震反演目標函數構建

觀測地震數據d是由反射系數r與震源產生的帶限子波w褶積而成[24],即

d=wr+n(9)

式中n為隨機噪聲。

常規單道地震反演的目標函數為

f(r)=||d-wr||2(2)(10)

基于式(10)的解,可以得到r,進一步通過遞推公式

Z(u)=Z(0)exp 2 r(u)(11)

估計波阻抗Z(u)。式中:Z(0)為第一層波阻抗;r(u)為r的離散形式,u=1,2,…,U,U為地層層數。

為了充分利用各種先驗地質信息提高反演精度,在寬帶約束反演(BCI,Broadband Constrained Inversion)方法的基礎上,筆者提出自適應時窗多道相關的疊后地震波阻抗反演方法。該方法通過查詢將地震數據劃分為不同尺度的時窗,分別計算每個時窗內地震道之間的相關性,并重構地震道作為反演約束條件,構建了基于多道相關的目標函數

f(Z)=min{||A seis-GZ||2(2)+λ||D-GZ||2(2)

+μ||Zpior-Z||2(2)}

式中:Z為Z的矩陣形式;G為地震子波時移矩陣與一階差分矩陣的乘積;D為原始地震數據;λ為正則化參數;Zpior為通過測井資料沿層內插、外推建立的先驗波阻抗模型;μ為模型正則化參數,可控制Zpior對反演的貢獻率。

2模型測試

2.1無噪模型測試

采用Overthrust模型(圖2)[25]測試所提方法的可行性。利用主頻為30 Hz的Ricker子波與圖2褶積得到正演記錄(圖3)。由于地震數據頻帶限制,反演缺少低頻成分,因此將圖2進行0~5 Hz的低通濾波得到低頻初始阻抗模型,以補充反演結果的低頻成分。設置主動學習模型的初始樣本集為圖2中CDP 50、CDP 75、CDP 150、CDP 200、CDP 225、CDP 275、CDP 300、CDP 350、CDP 400處的地震道及其波阻抗。選取初始樣本集的原則為:對于模型數據而言,選取地震數據中差異大的區域中的一個代表性樣本;對于實際數據,選取井位置的地震數據和阻抗(井數據的差異明顯)。根據初始樣本集查詢的自適應時窗如圖4所示。

圖5、圖6分別為模型不加噪反演結果、殘差。由圖5可見:相對于單道寬帶約束反演結果(圖5a),固定時窗多道相關(圖5b)和自適應時窗多道相關(圖5c)反演結果在一定程度上提高了橫向連續性;圖5b斷層邊界模糊,也未能較好地識別異常體,圖5c斷層邊界清晰,且較好地識別了異常體。因此,自適應時窗多道相關反演不僅可以提高反演結果的橫向連續性,還能刻畫斷層邊界和異常地質體邊界。由圖6可見,自適應時窗多道相關反演的殘差(圖6c)明顯小于固定時窗多道相關(圖6b)以及單道寬帶約束反演(圖6a),反演結果精度較高,但也存在一定局限性,如在深層高速體位置(150 ms)存在一定誤差。文中采用皮爾遜相關系數(Pearson?Correlation Coefficient,PCC)、決定系數(Coefficient of Determination,R2)和均方誤差(Mean Squared Er?ror,MSE)定量分析反演結果(表1)。PCC反映了兩個量之間的相關程度,越接近1相關程度越高;R2表示回歸分析中自變量對因變量的解釋,越接近1解釋效果越好;MSE反映了預測值與真實值之間的差異程度,越接近0預測精度越高。可見,自適應時窗多道相關反演的精度高于單道寬帶約束反演和固定時窗多道相關反演(表1),但計算效率較低(表2),主要原因是計算過程涉及大量矩陣存儲和計算。

2.2加噪模型測試

對圖3添加10%、20%和30%的高斯噪聲(表3),結果表明所提方法具有一定抗噪性,但隨著噪聲增加,抗噪性能下降。為了更直觀地了解自適應時窗多道相關反演方法的抗噪性,分析加入30%高斯噪聲正演記錄(圖7)的反演結果。圖8為圖7的重構合成記錄,圖9為圖8的反演結果。由圖可見,自適應時窗多道相關反演結果(圖9c)的橫向連續性強于單道寬帶約束反演結果(圖9a),并且對復雜構造的識別效果好于固定時窗多道相關反演結果(圖9b)。圖10為CDP 150處的單道反演結果。由圖可見,自適應時窗多道相關反演結果更接近真實阻抗,具有一定的抗噪性。上述模型測試結果表明:單道寬帶約束反演運行效率高,但反演結果的橫向連續性差,僅適用于地層結構較平滑地區;固定時窗多道相關反演結果的橫向連續性較好,但難以兼顧特殊構造(斷層、異常體等)信息,適用于無特殊構造、橫向連續性較差的地區;自適應時窗多道相關反演結果的橫向連續性較好,可以兼顧特殊構造信息,但計算效率較低,且在淺、深層速度差異較大的區域存在一定局限性(圖6c),適用于具有特殊構造的狹小地區。

3實際資料處理

A區現有30口井的聲波、密度測井資料。圖11為A區地震剖面,通過主動學習得到自適應時窗(圖12)。圖13為不同方法的反演結果。由圖可見,自適應時窗多道相關反演結果(圖13d)的橫向連續性優于單道寬帶約束(圖13a)、固定時窗多道相關(圖13b)及稀疏脈沖(圖13c)反演結果。井震標定結果表明,圖13d與測井資料吻合程度更好,并且與測井波阻抗的相關系數為0.85,且能提高反演的縱向分辨率。圖14為單道反演結果。由圖可見,自適應時窗多道相關反演結果與測井曲線基本吻合,驗證了方法的可行性和有效性。

4結束語

本文提出了一種自適應時窗多道相關的疊后地震波阻抗反演方法。該方法不僅考慮地震道之間的相關性,還同時兼顧特殊構造信息,使反演結果更合理。模型測試和實際數據測試表明,該方法反演結果的橫向連續性較好,且可以兼顧特殊構造信息,但計算效率較低,且在淺、深層速度差異較大的區域存在一定局限性,適用于具有特殊構造的狹小地區。

參考文獻

[1]李慶忠.論地震約束反演的策略[J].石油地球物理勘探,1998,33(4):423?438.

LI Qingzhong.On strategy of seismic restricted inver?sion[J].Oil Geophysical Prospecting,1998,33(4):423?438.

[2]LI G,LI H,MA Y,et al.Analysis of the ambiguity of log?constrained seismic impedance inversion[J].Petro?leum Science,2011,8(2):151?156.

[3]黃捍東,趙迪,任敦占,等.基于貝葉斯理論的薄層反演方法[J].石油地球物理勘探,2011,46(6):919?924.

HUANG Handong,ZHAO Di,REN Dunzhan,et al.A thin bed inversion method based on Bayes theory[J].Oil Geophysical Prospecting,2011,46(6):919?924.

[4]宋磊,印興耀,宗兆云,等.基于先驗約束的深度學習地震波阻抗反演方法[J].石油地球物理勘探,2021,56(4):716?727.

SONG Lei,YIN Xingyao,ZONG Zhaoyun,et al.Deep learning seismic impedance inversion based on prior constraints[J].Oil Geophysical Prospecting,2021,56(4):716?727.

[5]HAMID H,PIDLISECKY A.Multitrace impedanceinversion with lateral constraints[J].Geophysics,2015,80(6):M101?M111.

[6]FOMEL S.Shaping regularization in geophysical?esti?mation problems[J].Geophysics,2007,72(2):R29?R36.

[7]耿偉恒,陳小宏,李景葉,等.基于L1?2正則化的地震波阻抗“塊”反演[J].石油地球物理勘探,2022,57(6):1409?1417.

GENG Weiheng,CHEN Xiaohong,LI Jingye,et al.Seismic“blocky”acoustic impedance inversion based on L1?2 regularization[J].Oil Geophysical Prospecting,2022,57(6):1409?1417.

[8]ULRYCH T J,SACCHI M D.Information?Based In?version and Processing with Applications[M].Elsevier Science,New York,America,2005.

[9]楊俊,尹成,代榮獲,等.疊后地震數據自適應正則化參數稀疏約束反演方法[J].地球物理學進展,2020,35(6):2259?2264.

YANG Jun,YIN Cheng,DAI Ronghuo,et al.Sparse constrained post?stack seismic data inversion with adap?tive selection of regularization parameters[J].Progress in Geophysics,2020,35(6):2259?2264.

[10]TIKHONOV A N.Solution of incorrectly formulatedproblems and the regularization method[J].Soviet Math,1963,5(4):1035?1038.

[11]王寧,周輝,王玲謙,等.數據驅動的塊排列正則化多道疊前地震反演[J].地球物理學報,2022,65(7):2681?2692.

WANG Ning,ZHOU Hui,WANG Lingqian,et al.Data?driven multichannel pre?stack seismic inversion via patch?ordering regularization[J].Chinese Journal of Geophysics,2022,65(7):2681?2692.

[12]ZHOU D,YIN X,ZONG Z.Multi?trace basis?pursuitseismic inversion for resolution enhancement[J].Geo?physical Prospecting,2019,67(3):519?531.

[13]YUAN S,WANG S,TIAN N,et al.Stable inversion?based multitrace deabsorption method for spatial conti?nuity preservation and weak signal compensation[J].Geophysics,2016,81(3):V199?V212.

[14]MA M,WANG S,YUAN S,et al.Multichannel spa?tially correlated reflectivity inversion using block sparse Bayesian learning[J].Geophysics,2017,82(4):V191?V199.

[15]段友祥,崔樂樂,孫歧峰,等.波動方程正演引導的深度學習地震波形反演[J].石油地球物理勘探,2023,58(3):485?494.

DUAN Youxiang,CUI Lele,SUN Qifeng,et al.Deep learning seismic waveform inversion based on the for?ward modeling guidance of wave equation[J].Oil Geo?physical Prospecting,2023,58(3):485?494.

[16]宗兆云,印興耀,吳國忱.基于疊前地震縱橫波模量直接反演的流體檢測方法[J].地球物理學報,2012,55(1):284?292.

ZONG Zhaoyun,YIN Xingyao,WU Guochen.Fluid identification method based on compressional and shear modulus direct inversion[J].Chinese Journal of Geophysics,2012,55(1):284?292.

[17]霍國棟,杜啟振,王秀玲,等.縱向和橫向同時約束AVO反演[J].地球物理學報,2017,60(1):271?282.

HUO Guodong,DU Qizhen,WANG Xiuling,et al.AVO inversion constrained simultaneously in vertical and lateraldirections[J].Chinese Journal of Geo?physics,2017,60(1):271?282.

[18]GHOLAMI A.A fast automatic multichannel blindseismic inversion for high?resolution impedance reco?very[J].Geophysics,2016,81(5):V357?V364.

[19]HAMID H,Pidlisecky A.Structurally constrained im?pedance inversion[J].Interpretation,2016,4(4):T577?T589.

[20]HUANG G,CHEN X,QU S,et al.Directional totalvariation regularized high?resolution prestack AVA in?version[J].IEEE Transactions on Geoscience and Re?mote Sensing,2022,DOI:10.1109/TGRS.2021.3078293.

[21]印興耀,楊亞明,李坤,等.地震數據互相關驅動的多道反演方法[J].地球物理學報,2020,63(10):3827?3837.

YIN Xingyao,YANG Yaming,LI Kun,et al.Multi?trace inversion driven by cross?correlation of seismic data[J].Chinese Journal of Geophysics,2020,63(10):3827?3837.

[22]HUANG G,CHEN X,LI J,et al.The slope?attribute?regularized high?resolution prestack seismic inversion[J].Surveys in Geophysics,2021,42(3):625?671.

[23]GEIFMAN Y,El?Yaniv R.Deep active learning witha neural architecture search[C].Proceedings of the 33rd International Conference on Neural Information Processing Systems,2019,5976?5986.

[24]ROBINSON E A.Predictive decomposition of timeseries with application to seismic exploration[J].Geo?physics,1967,32(3):418?484.

[25]WEIMER P,DAVIS T L.Applications of 3?D Seis?mic Data to Exploration and Production[M].Ameri?can Association of Petroleum Geologists,1996.

(本文編輯:劉勇)

作者簡介

彭真碩士研究生,1999年生;2021年獲長江大學計算機科學學院軟件工程專業學士學位;目前在該校地球物理與石油資源學院攻讀地質工程專業碩士學位,主要研究方向為智能地震反演、地震相控建模以及儲層預測。

猜你喜歡
方法
中醫特有的急救方法
中老年保健(2021年9期)2021-08-24 03:52:04
高中數學教學改革的方法
河北畫報(2021年2期)2021-05-25 02:07:46
化學反應多變幻 “虛擬”方法幫大忙
變快的方法
兒童繪本(2020年5期)2020-04-07 17:46:30
學習方法
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
最有效的簡單方法
山東青年(2016年1期)2016-02-28 14:25:23
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
主站蜘蛛池模板: 人妻少妇久久久久久97人妻| 毛片网站在线看| 青青草原国产精品啪啪视频| 欧美日韩午夜| AV不卡国产在线观看| 国产欧美专区在线观看| 欧美精品另类| 在线视频一区二区三区不卡| 亚洲国产在一区二区三区| 精品视频福利| 亚洲精品国偷自产在线91正片| 日韩毛片基地| 自拍偷拍欧美日韩| 自慰网址在线观看| 无码综合天天久久综合网| 成年A级毛片| 亚洲综合第一页| 婷婷色丁香综合激情| 福利小视频在线播放| 99精品热视频这里只有精品7| 国产精品久久久久久久久| 在线观看91香蕉国产免费| 色丁丁毛片在线观看| 欧美日韩免费在线视频| 四虎国产成人免费观看| 她的性爱视频| 制服丝袜在线视频香蕉| 国产成人一区在线播放| 國產尤物AV尤物在線觀看| 国产亚洲高清视频| 久久亚洲黄色视频| 色AV色 综合网站| 99ri精品视频在线观看播放| 久久免费视频6| 国产三级韩国三级理| 99视频在线免费| 国产一区二区三区视频| 国产在线一区视频| 亚洲天堂网在线观看视频| 美女国产在线| 免费女人18毛片a级毛片视频| 国产一级二级三级毛片| 成人福利视频网| 国产成人91精品| 成人国产精品一级毛片天堂| 三级国产在线观看| 精品三级网站| 最新日韩AV网址在线观看| 国产成人高清精品免费| 亚洲va视频| 色综合中文字幕| 国产精品视频猛进猛出| 国产在线八区| 亚洲天堂首页| 亚洲系列无码专区偷窥无码| 美女被躁出白浆视频播放| 国产一级在线播放| 中文无码影院| 亚洲视频一区在线| 人妻少妇久久久久久97人妻| 五月六月伊人狠狠丁香网| 毛片免费在线视频| 99久久精品免费看国产电影| 中文字幕人成人乱码亚洲电影| 国产又粗又爽视频| 青青青视频蜜桃一区二区| 亚洲欧美在线综合图区| 毛片在线播放a| 亚洲精品无码AⅤ片青青在线观看| 久久综合色88| 亚洲人成日本在线观看| 特级aaaaaaaaa毛片免费视频 | 三上悠亚一区二区| 亚洲精品第一在线观看视频| 午夜a级毛片| 欧美精品在线看| 国产理论最新国产精品视频| 伊人色天堂| 54pao国产成人免费视频| 扒开粉嫩的小缝隙喷白浆视频| 亚洲欧美日韩中文字幕在线| 九九久久99精品|