亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        基于Huber函數(shù)的頻率域全波形反演

        2013-12-01 09:07:08呂曉春顧漢明成景旺
        石油物探 2013年5期
        關(guān)鍵詞:范數(shù)高斯反演

        呂曉春,顧漢明,成景旺

        (1.中國地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院,湖北武漢430074;2.中國地質(zhì)大學(xué)(武漢)構(gòu)造與油氣資源教育部重點(diǎn)實(shí)驗(yàn)室,湖北武漢430074)

        利用人工地震獲取的地震波信息反演地下介質(zhì)物性參數(shù)是地球物理學(xué)的重要研究內(nèi)容。地震波場反演方法很多,如相位反演成像、振幅反演成像以及全波形反演成像等等。其中全波形反演不同于走時(shí)反演、振幅反演等傳統(tǒng)的成像方法,它利用地震記錄的整體信息來進(jìn)行反演,是一種高精度的成像方法。20世紀(jì)80年代,直接通過求取模型參數(shù)擾動(dòng)計(jì)算Frechet導(dǎo)數(shù)的做法工作量極為龐大,使全波形反演局限于一維模型或簡單模型。Tarantola[1]借鑒逆時(shí)偏移的思想,利用正向傳播波場與殘差波場的互相關(guān)求取梯度方向,進(jìn)行了時(shí)間域的全波形反演研究。Pratt等[2-4]將全波形反演引入到頻率域中,根據(jù)炮點(diǎn)和接收點(diǎn)可互易性原理,導(dǎo)出了求解目標(biāo)函數(shù)梯度的高效算法,并利用高斯牛頓法對非均勻介質(zhì)的速度進(jìn)行了全波形反演。至此,地球物理學(xué)者開始對全波形反演開展更深入的研究。國外近十年中頻率域全波形反演已被廣 泛 應(yīng) 用 于 建 立 高 分 辨 率 的 速 度 模 型[5-7]。Virieux[8]總結(jié)了當(dāng)前全波形反演所取得的一些成果和認(rèn)識,指出了該方法所面對的難點(diǎn)以及下一步研究方向和趨勢。國內(nèi)很多學(xué)者也做了大量的研究工作,吳國忱等[9-10]以及殷文等[11]對頻率域的聲波和彈性波有限差分正演模擬做了很詳細(xì)的研究,為頻率域反演提供了前提手段;許焜等[12-13]基于頻域率有限元方法進(jìn)行了反演研究;龍桂華等[14]對粘彈聲波介質(zhì)進(jìn)行了頻率域反演;卞愛飛等[15]提出了自適應(yīng)剝層法全波形反演。

        隨著計(jì)算機(jī)技術(shù)的發(fā)展,以及地震處理要求的提高,全波形反演方法的高分辨率特性引起了人們的注意。然而,實(shí)際地震信號的低頻缺失以及噪聲的存在制約著全波形反演的應(yīng)用,其中噪聲的存在是一個(gè)非常關(guān)鍵的問題,因此有必要找到一個(gè)合適的目標(biāo)函數(shù)。L1范數(shù)相對于L2范數(shù)抗噪性強(qiáng)[16-17],但數(shù)據(jù)誤差接近于零時(shí),L1范數(shù)準(zhǔn)則求取的目標(biāo)函數(shù)梯度會(huì)出現(xiàn)不穩(wěn)定性,而L2范數(shù)則比較穩(wěn)定,不會(huì)出現(xiàn)這個(gè)問題。Guitton[18]在使用擬牛頓法進(jìn)行地震反演時(shí),通過求取一個(gè)實(shí)數(shù)型的Huber目標(biāo)函數(shù)的梯度來實(shí)現(xiàn)。Ha等[19]基于Huber函數(shù)構(gòu)建反演目標(biāo)函數(shù)來提高算法的抗噪性和穩(wěn)定性。韓淼等[20]也進(jìn)行了反演目標(biāo)函數(shù)的幾種范數(shù)準(zhǔn)則的對比研究。

        為了既滿足穩(wěn)定性,又滿足抗噪性,我們將L1范數(shù)和L2范數(shù)結(jié)合來建立目標(biāo)函數(shù),提出一種基于Huber函數(shù)的頻率域全波形反演新方法。頻率域全波形反演是在復(fù)數(shù)域進(jìn)行的,因此引入了復(fù)數(shù)形式的Huber函數(shù),并推導(dǎo)了基于Huber準(zhǔn)則的目標(biāo)函數(shù)的梯度表達(dá)式。通過在幾個(gè)簡單的合成記錄中加入隨機(jī)脈沖噪聲、高斯噪聲和連續(xù)噪聲,驗(yàn)證了該方法的應(yīng)用效果。

        1 復(fù)數(shù)形式的Huber函數(shù)

        Huber[21]給出實(shí)數(shù)形式的Huber函數(shù)表達(dá)式:

        其中,r為計(jì)算值與觀測值的誤差,是一個(gè)實(shí)數(shù)變量;ε是一個(gè)閾值,用來判定不同大小的誤差所對應(yīng)的目標(biāo)函數(shù)表達(dá)式(圖1)。類似實(shí)數(shù)域的Huber函數(shù),對于復(fù)數(shù)r=a+bi,我們定義兩種不同的范數(shù),其中L1范數(shù)計(jì)算表達(dá)式為,L2范數(shù)計(jì)算表達(dá)式為,就可以得到復(fù)數(shù)域的Huber函數(shù)表達(dá)式:

        圖1 Huber函數(shù)示意

        從公式(2)中可以看出,在誤差大于閾值ε時(shí),目標(biāo)函數(shù)選取L1范數(shù)準(zhǔn)則。我們知道L1范數(shù)是所有范數(shù)中對數(shù)據(jù)加權(quán)最小的[22],如果數(shù)據(jù)非常精確,那么某一預(yù)測值遠(yuǎn)離其觀測值的地方便是重要的影響因素,這時(shí)應(yīng)選擇L2范數(shù),因?yàn)樗鼘^大的誤差加權(quán)也大。而對于包含噪聲的合成記錄來說,噪聲的存在使得數(shù)據(jù)不是很準(zhǔn)確,應(yīng)使用L1范數(shù)來減小這部分較大誤差的權(quán)值。由于L1范數(shù)總是大于L2范數(shù),故(2)式還可以改寫成本文反演中采用公式(3)來建立目標(biāo)函數(shù),其中的閾值ε的選取非常關(guān)鍵,我們采用Bube等[23]提出的方法,閾值ε的大小取為計(jì)算值與觀測值誤差數(shù)列的標(biāo)準(zhǔn)差的0.6倍。

        2 頻域率全波形反演基本理論

        2.1 頻率域波動(dòng)方程正演計(jì)算

        頻率域聲波波動(dòng)方程可以表示為

        式中:ρ(x)為密度;k(x)為體變模量;ω 為角頻率;u(x,ω)和s(x,ω)分別為波場壓力值與震源。將方程(4)加入PML邊界條件后采用有限差分算子將其離散,則(4)式可簡單表示為

        式中:A(ω)為與頻率和介質(zhì)速度有關(guān)的阻抗矩陣。假如模型空間網(wǎng)格數(shù)為(Nx,Nz),那么A(ω)是一個(gè)(Nx×Nz)×(Nx×Nz)的有限差分算子矩陣。該系數(shù)矩陣為大型稀疏矩陣,其非零元素對稱地分布在對角線上,但是它并非對稱矩陣,其值的大小是不一樣的。求解該大型稀疏矩陣方程組常用的方法有迭代法和直接分解法。對于多震源問題,我們采用直接分解(如LU分解)來求解,因?yàn)橄禂?shù)矩陣A(ω)只需分解一次,其LU分解結(jié)果就可以用于更多的新的震源求解,可以大量節(jié)省計(jì)算量,而且在全波形反演中的迭代反演上更為重要,因?yàn)槠溆写罅康恼鹪匆约疤撜鹪辞蠼?,每一次迭代過程中系數(shù)矩陣A(ω)不變,其LU分解可以多次重復(fù)使用。

        2.2 目標(biāo)函數(shù)的建立以及梯度的求解

        反問題的求解在數(shù)學(xué)上其實(shí)是一個(gè)最優(yōu)化求解問題,其中目標(biāo)函數(shù)的構(gòu)造是一個(gè)關(guān)鍵。根據(jù)上述復(fù)數(shù)形式的Huber表達(dá)式,我們構(gòu)造如下的目標(biāo)函數(shù):

        式中:p為模型變量;Ns為炮點(diǎn)數(shù);Nr為接收點(diǎn)數(shù);rs,r為計(jì)算波場與觀測波場的殘差:

        結(jié)合(8)式和(9)式,我們將基于 Huber函數(shù)得到的目標(biāo)函數(shù)梯度表達(dá)式統(tǒng)一表示為

        式中:r=(r1,r2,…,rNr),

        從(10)式中我們看到梯度求解的關(guān)鍵在于Jacobin矩陣的求解。對(5)式兩邊參數(shù)pk進(jìn)行求導(dǎo),有

        由于(13)式和(5)式的形式一樣,而且它的阻抗矩陣也一樣,故可以直接利用正演的LU分解結(jié)果求解偏導(dǎo)數(shù)波場。將(13)式代入(10)式可得

        其中A[]-1Tr*為殘差的反傳播場,我們只需將各個(gè)接收點(diǎn)處的殘差組成新的震源再正演一次便可得到。從(14)式可以看出目標(biāo)函數(shù)對速度變量的梯度可以通過計(jì)算虛震源與殘差反傳播場的零延遲卷積得到。

        2.3 預(yù)梯度法反演

        地球物理反演通過求取目標(biāo)函數(shù)的最小值,不斷對模型參數(shù)進(jìn)行修改,求取滿足條件的最佳模型。本文反演過程中采用近似Hessian矩陣的對角線來做梯度的預(yù)條件算子[14],其模型更新量可以表示為:

        式中:α為待定的步長;diag(Ha)表示近似 Hessian矩陣Ha的對角線矩陣,用來對梯度進(jìn)行修正,以加快收斂速度;λ為可調(diào)節(jié)的正則化因子;I表示單位矩陣為目標(biāo)函數(shù)的梯度。其中Ha可以表示為

        從(16)式不難發(fā)現(xiàn)Hessian矩陣可以表示成偏導(dǎo)數(shù)的互相關(guān),根據(jù)(13)式可求出Hessian矩陣,令i=j(luò)就可求出它的對角元素。

        3 數(shù)值模型試算

        為了驗(yàn)證基于Huber函數(shù)的頻率域全波形反演方法的實(shí)用性以及正確性,在原始合成記錄上分別加入了隨機(jī)脈沖噪聲、連續(xù)噪聲和高斯噪聲來進(jìn)行模型試算。

        3.1 含隨機(jī)脈沖噪聲的數(shù)值試算

        通過一個(gè)簡單的地質(zhì)模型來進(jìn)行隨機(jī)脈沖噪聲的數(shù)值試算。模型計(jì)算范圍為5km×2km,一共有5層,中間有一個(gè)高速體(圖2)。原始合成記錄通過頻率域有限差分正演模擬得到,正演過程中采用主頻為10Hz的雷克子波,頻率間隔為0.24Hz,最大頻率為30Hz。共有51炮,炮間距100m。每炮有251道接收,道間距為20m。圖3a為炮點(diǎn)位置位于0處的原始合成記錄,圖3b為加了噪聲的記錄,其中有9個(gè)充零道和4個(gè)脈沖噪聲。

        圖2 簡單模型

        圖3 震源位于0處的模擬炮集記錄

        圖4 含隨機(jī)脈沖噪聲數(shù)據(jù)反演結(jié)果

        分別采用基于L2范數(shù)的最小平方反演法和本文提出的Huber反演方法對不含噪聲和含噪聲的記錄進(jìn)行反演。反演過程中從2.4到20.0Hz共選取20個(gè)反演頻率,每個(gè)頻率反演40次。反演采用的初始模型為一個(gè)線性模型(圖4a),初始速度為1 500m/s,最深處的速度為3 000m/s,速度計(jì)算為v=1 500+kz,k為15m/s,z為網(wǎng)格點(diǎn)數(shù)。圖4b和4c分別為采用最小平方法對原始合成記錄(圖3a)和含隨機(jī)噪聲數(shù)據(jù)(圖3b)的反演結(jié)果。由圖4b可以看到,最小平方法對沒有噪聲的數(shù)據(jù)進(jìn)行反演可以得到較好的速度模型,根據(jù)Snell射線原理,高速體下方透射波較少,故該界面與理論有所差異,但模型的大部分層位都較好地反映出來。相對而言,對含有噪聲的數(shù)據(jù)最小平方法反演結(jié)果較差(圖4c),除了第一層較好地被反演出來,其他各層的分界面都反演較差。圖4d為采用本文方法對含隨機(jī)噪聲數(shù)據(jù)的反演結(jié)果,對比圖4c和圖4d可看出,基于Huber函數(shù)的全波形反演受噪聲影響小,反演結(jié)果較好。

        3.2 含連續(xù)噪聲的數(shù)值試算

        采用和3.1節(jié)同樣的模型(圖2)、觀測系統(tǒng)以及正、反演參數(shù)。如圖5所示,在原始合成單炮記錄(圖5a)上加入了2個(gè)不同斜率的連續(xù)噪聲(圖5b和圖5c)。為了和原始記錄的頻譜相一致,連續(xù)噪聲同樣采用主頻為10Hz的雷克子波,其中連續(xù)噪聲的振幅值隨著偏移距而改變[24]。

        從2個(gè)含連續(xù)噪聲數(shù)據(jù)的反演結(jié)果(圖6)來看,連續(xù)噪聲對反演結(jié)果的影響較大。第一種斜率噪聲存在的情況下,利用最小平方反演法(圖6a)沒有將高速體反演出來,只反演出上界面,而且出現(xiàn)了假的層位;采用本文方法(圖6c)可以大致反演出高速體,而且下面也沒有出現(xiàn)假的層位。第二種噪聲存在的情況下,最小反演法反演(圖6b)的第一層的深度與理論不吻合,高速體同樣沒有反演出來;而采用本文方法得到的反演結(jié)果(圖6d)第一層的埋深有所改善,且高速體也被反演出來,但與不含噪聲反演結(jié)果(圖4b)相比,還有所差別。

        圖5 震源位于0處的模擬炮集記錄

        圖6 含連續(xù)噪聲數(shù)據(jù)反演結(jié)果

        3.3 含高斯噪聲的數(shù)值試算

        我們采用Marmousi模型(圖7)的合成記錄加入高斯噪聲來進(jìn)行反演測試。Marmousi模型計(jì)算范圍取為9.00km×3.18km,網(wǎng)格大小為451×159。原始合成記錄也是通過頻率域有限差分正演模擬得到的,正演過程中采用主頻為10Hz的雷克子波,頻率間隔為0.24Hz,正演最大頻率為30Hz。共有91炮,炮間距為100m。每一炮451道接收,道間距為20m。

        圖8分別給出了沒有加噪聲(圖8a)和加了信噪比為5dB的高斯噪聲(圖8b)的合成記錄。高斯噪聲的產(chǎn)生采用Box-Muller方法[25]。

        圖7 Marmousi模型

        圖8 震源位于3 100m處的模擬炮集記錄

        分別采用最小平方反演法和本文提出的Huber反演方法對圖8所示不含噪聲和含有噪聲的原始合成記錄進(jìn)行反演。反演過程中從2.4到20.0Hz共選取20個(gè)反演頻率,每個(gè)頻率反演40次。反演采用的初始模型為一個(gè)線性模型(圖9a),速度計(jì)算為v=1 500+kz,k為17m/s,z為網(wǎng)格點(diǎn)數(shù)。圖9b和圖9c分別為采用最小平方法對原始合成記錄和含高斯噪聲數(shù)據(jù)的反演結(jié)果,圖9d為本文方法對含高斯噪聲數(shù)據(jù)的反演結(jié)果。從3種反演結(jié)果來看,最小平方法對高斯噪聲敏感度較低,合成記錄中加入高斯噪聲也可以得到好的反演結(jié)果。由此可以看出最小平方法和Huber函數(shù)反演法受高斯噪聲的影響都較小,除了速度模型光滑度較差外,各個(gè)層的分界面均可以較好地反演出來。

        圖9 含高斯噪聲數(shù)據(jù)反演結(jié)果

        4 結(jié)束語

        針對傳統(tǒng)的全波形反演采用L2范數(shù)準(zhǔn)則構(gòu)建目標(biāo)函數(shù)抗噪性差和L1范數(shù)則可能導(dǎo)致反演不穩(wěn)定的問題,我們引入復(fù)數(shù)形式的Huber函數(shù)準(zhǔn)則來建立目標(biāo)函數(shù),同時(shí)考慮L2范數(shù)和L1范數(shù)的優(yōu)勢,提出了一種頻率域全波形反演新方法。通過以上理論研究以及加入不同噪聲的數(shù)值模型數(shù)據(jù)試算,可以得到以下認(rèn)識:

        1)在頻率域全波形反演中,基于L2范數(shù)的最小平方法反演在無噪聲的情況下能夠得到精確的地下模型結(jié)構(gòu)。但是L2范數(shù)對噪聲的敏感度較高,在噪聲存在的情況下不能得到理想的反演結(jié)果。

        2)基于Huber函數(shù)準(zhǔn)則建立的目標(biāo)函數(shù),具有L1范數(shù)的抗噪性以及L2范數(shù)的穩(wěn)定性。在噪聲存在的情況下,依然能夠得到較好的反演結(jié)果。

        3)L2范數(shù)對隨機(jī)脈沖噪聲和連續(xù)噪聲敏感度較高,而對高斯噪聲的敏感度相對來說較低。因此,對加入高斯噪聲的合成記錄,基于L2范數(shù)和基于Huber函數(shù)的方法反演效果均不錯(cuò);而對加入隨機(jī)脈沖噪聲和連續(xù)噪聲的合成記錄,基于Huber函數(shù)反演要比基于L2范數(shù)反演效果好。

        [1]Tarantola A.Inversion of seismic reflection data in the acoustic approximation[J].Geophysics,1984,49(8):1259-1266

        [2]Pratt R G,Shin C S,Hicks G J.Guass-Newton and full Newton methods in frequency-space seismic waveform inversion[J].Geophysics,1998,63(2):341-362

        [3]Pratt R G.Seismic waveform inversion in the frequency domain,Part 1:theory and verification in a physical scale model[J].Geophysics,1999,64(3):888-901

        [4]Pratt R G.Seismic waveform inversion in the frequency domain,Part 2:fault delineation in sediments using crosshole data[J].Geophysics,1999,64(3):902-914

        [5]Hicks G J,Pratt R G.Reflection waveform inversion using local decent methods:Estimating attenuation and velocity over a gas-sand deposit[J].Geophysics,2001,66(2):598-612

        [6]Malinowski M,Operto S.Quantitative imaging of the Permo-Mesozoic complex and its basement by frequency domain waveform tomography of wide-aperture seismic data from the Polish Basin[J].Geophysical Prospecting,2008,56(6):805-825

        [7]Ben-Hadj-Ali H,Operto S,Virieux J.Velocity model building by 3Dfrequency-domain,full-waveform inversion of wide-aperture seismic data[J].Geophysics,2008,73(5):VE101-VE117

        [8]Virieux J,Operto S.An overview of full-waveform inversion in exploration geophysics[J].Geophysics,2009,74(6):WCC1-WCC26

        [9]吳國忱,梁鍇.VTI介質(zhì)頻率-空間域準(zhǔn)P波正演模擬[J].石油地球物理勘探,2005,40(5):535-545 Wu G C,Liang K.Quasi P-wave forward modeling in frequency-space domain in VTI media[J].Oil Geophysical Prospecting for Petroleum,2005,40(5):535-545

        [10]吳國忱,梁鍇.VTI介質(zhì)準(zhǔn)P波頻率空間域組合邊界條件研究[J].石油物探,2005,44(4):301-307 Wu G C,Liang K.Combined boundary conditions of quasi P-wave within frequency-space domain in VTI media,Oil Geophysical Prospecting for Petroleum,2005,44(4):301-307

        [11]殷文,印興耀,吳國忱,等.高精度頻率域彈性波方程有限差分方法及波場模擬[J].地球物理學(xué)報(bào),2006,49(2):561-568 Yin W,Yin X Y,Wu G C,et al.The method of finite difference of high precision elastic wave equation in the frequency domain and wavefiled simulation[J].Chinese Journal Geophysics(in Chinese),2006,49(2):561-568

        [12]許琨,王妙月.聲波方程頻率域有限元參數(shù)反演[J].地球物理學(xué)報(bào),2001,44(6):852-864 Xu K,Wang M Y.Finite element inversion of the coefficients of acoustic equation in frequency domain[J].Chinese Journal Geophysics(in Chinese),2001,44(6):852-864

        [13]許琨,王妙月.利用地質(zhì)規(guī)則塊體建模方法的頻率域有限元彈性波速度反演[J].地球物理學(xué)報(bào),2004,47(4):708-717 Xu K,Wang M Y.Frequency-domain finite-element inversion of elastic-wave velocity using the geological regular-blocky-model method[J].Chinese Journal Geophysics(in Chinese),2004,47(4):708-717

        [14]龍桂華,李小凡,張美根,等.頻率域粘彈性聲波透射波形速度反演[J].地震學(xué)報(bào),2009,31(1):32-41 Long G H,Li X F,Zhang M G,et al.Visco-acoustic transmission waveform inversion for velocity structure in space-frequency domain[J].Acta Seismologica Sinica,2009,31(1):32-41

        [15]Bian A F,Yu W H.Layer-Stripping full waveform inversion with damped seismic reflection data[J].Journal of Earth Science,2011,22(2):241-249

        [16]Claerbout J,Muir J.Robust modeling with erratic data[J].Geophysics,1973,38(5):326-344

        [17]Crase E,Pica A M,Noble M,et al.Robust elastic nonlinear waveform inversion:application to real data[J].Geophysics,1990,55(5):527-538

        [18]Guitton A,Symes W.Robust inversion of seismic data using the Huber norm[J].Geophysics,2003,68(4):1310-1318

        [19]Ha T,Chung W,Shin C.Waveform inversion using a back-propagation algorithm and a Huber function norm[J].Geophysics,2009,74(3):R15-R24

        [20]韓淼,韓立國,崔杰,等.頻率域全波形反演目標(biāo)函數(shù)的幾種范數(shù)準(zhǔn)則的對比研究[J].中國地球物理學(xué)會(huì)二十七屆年會(huì)論文集,2011,730 Han M,Han L G,Cui J,et al.The contrastive research of several norm-criterions of the objective function in FWI[J].Chinese Geophysical Society 27thannual meeting proceedings,2011,730

        [21]Huber J.Robust regression:Asymptotics,conjectures and Monte Carlo:annals of statistics[J].Geophysics,1973,38(10):799-821

        [22]姚姚.地球物理反演基本原理與理論研究[M].湖北武漢:中國地質(zhì)大學(xué)出版社,2002:17-18 Yao Y.Geophysical inversion of the basic theory and method of application[M].Wuhan:China University of Geosciences Press,2002:17-18

        [23]Bube K,Nemeth T.Fast line searches for the robust solution of linear systems in the hybrid L1/L2and Huber norms[J].Geophysics,2007,72(2):A13-A17

        [24]Lu W,Zhang W,Liu D.Local linear coherent noise attenuation based on local polynomial approximation[J].Geophysics,2006,71(6):V163-V169

        [25]Box G,Muller M.A note on the generation of random normal deviates[J].Annals of Mathematical Statistics,1958,29(2):610-611

        猜你喜歡
        范數(shù)高斯反演
        小高斯的大發(fā)現(xiàn)
        反演對稱變換在解決平面幾何問題中的應(yīng)用
        天才數(shù)學(xué)家——高斯
        基于低頻軟約束的疊前AVA稀疏層反演
        基于自適應(yīng)遺傳算法的CSAMT一維反演
        基于加權(quán)核范數(shù)與范數(shù)的魯棒主成分分析
        矩陣酉不變范數(shù)H?lder不等式及其應(yīng)用
        有限域上高斯正規(guī)基的一個(gè)注記
        一類具有準(zhǔn)齊次核的Hilbert型奇異重積分算子的范數(shù)及應(yīng)用
        疊前同步反演在港中油田的應(yīng)用
        日韩精品一区二区三区中文| 少妇隔壁人妻中文字幕| 青青草在线免费观看视频| 国产女人18毛片水真多18精品| 尤物视频在线观看| 日本亚洲欧美在线观看| 国内精品极品久久免费看| 免费亚洲老熟熟女熟女熟女| 一本色道久久88综合日韩精品| 国产99r视频精品免费观看| 亚洲国产精品成人久久av| 亚洲一区二区三区偷拍女| 又爽又黄又无遮挡的视频| 国产v视频| 国产女人高潮的av毛片| 精品激情成人影院在线播放| 国产婷婷色综合av蜜臀av| 亚洲国产精品久久久久秋霞1| 亚洲最大的av在线观看| 日韩性爱视频| 337人体做爰大胆视频| 亚洲中文字幕久爱亚洲伊人| 蜜桃视频一区视频二区| 国产激情视频在线观看的| 黄色成人网站免费无码av| 色婷婷一区二区三区四区| 日本不卡的一区二区三区中文字幕| 一区二区三区人妻无码| 精品亚洲国产探花在线播放| 国产人妖直男在线视频| 人人妻人人澡人人爽国产一区| 国产av无码专区亚洲av手机麻豆 | 男女做那个视频网站国产| 免费人成激情视频在线观看冫| 在线观看av中文字幕不卡| 少妇裸淫交视频免费看| 国产成人av无码精品| 人与嘼交av免费| 午夜av内射一区二区三区红桃视| 亚洲av日韩一卡二卡| 国产综合久久久久|