• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    TTI介質(zhì)聲波方程分裂式PML吸收邊界條件研究

    2017-06-29 02:17:40丁仁偉李麗青李福元
    石油物探 2017年3期
    關(guān)鍵詞:波場(chǎng)橫波邊界條件

    張 衡,劉 洪,李 博,丁仁偉,李麗青,李福元

    (1.國土資源部海底礦產(chǎn)資源重點(diǎn)實(shí)驗(yàn)室,中國地質(zhì)調(diào)查局廣州海洋地質(zhì)調(diào)查局,廣東廣州510075;2.中國科學(xué)院地質(zhì)與地球物理研究所,中國科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室,北京100029;3.中國石油化工股份有限公司石油物探技術(shù)研究院,江蘇南京211103;4.山東科技大學(xué)地球科學(xué)與工程學(xué)院,山東青島266590)

    TTI介質(zhì)聲波方程分裂式PML吸收邊界條件研究

    張 衡1,劉 洪2,李 博3,丁仁偉4,李麗青1,李福元1

    (1.國土資源部海底礦產(chǎn)資源重點(diǎn)實(shí)驗(yàn)室,中國地質(zhì)調(diào)查局廣州海洋地質(zhì)調(diào)查局,廣東廣州510075;2.中國科學(xué)院地質(zhì)與地球物理研究所,中國科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室,北京100029;3.中國石油化工股份有限公司石油物探技術(shù)研究院,江蘇南京211103;4.山東科技大學(xué)地球科學(xué)與工程學(xué)院,山東青島266590)

    針對(duì)傾斜橫向各向同性(TTI)介質(zhì)聲波波動(dòng)方程的特點(diǎn),研究了TTI介質(zhì)分裂式完全匹配層(Split perfectly matched layer,SPML)吸收邊界條件。首先對(duì)常見的幾種TTI介質(zhì)聲波波動(dòng)方程進(jìn)行了歸納,并從TTI介質(zhì)波場(chǎng)傳播穩(wěn)定性的角度進(jìn)行對(duì)比分析,結(jié)果表明,引入橫波分量的TTI介質(zhì)縱橫波耦合方程適用于TTI介質(zhì)。然后從TTI介質(zhì)縱橫波耦合二階波動(dòng)方程出發(fā),推導(dǎo)得到其一階波動(dòng)方程的形式,進(jìn)而推導(dǎo)出一階波動(dòng)方程形式的SPML波動(dòng)方程,并給出了高階交錯(cuò)網(wǎng)格有限差分算法的具體實(shí)現(xiàn)過程。數(shù)值模擬結(jié)果表明,SPML吸收邊界條件能達(dá)到很好的人工邊界反射吸收效果,相比優(yōu)化海綿吸收邊界條件,其人工邊界反射吸收效果更好。

    TTI介質(zhì);分裂式完全匹配層;縱橫波耦合方程;一階波動(dòng)方程;高階交錯(cuò)網(wǎng)格有限差分

    Keywords:TTI media,split perfectly matched layer,P-wave and SV-wave coupled equation,first-order wave equation,high-order staggered-grid finite-difference

    各向異性正演數(shù)值模擬是各向異性逆時(shí)偏移和各向異性全波形反演的基礎(chǔ)[1]。1986年,THOMSEN[2]提出了在大多數(shù)情況下地球介質(zhì)為弱各向異性(各向異性參數(shù)為10%~20%)的觀點(diǎn),并在弱各向異性假設(shè)下提出弱各向異性橫向各向同性(TI)理論,即用3個(gè)各向異性參數(shù)——Thomsen各向異性參數(shù)ε(縱波各向異性參數(shù))、δ(變異系數(shù))和γ(橫波各向異性參數(shù))表述P波、SV波及SH波特征。此后30年,關(guān)于各向異性介質(zhì)的研究基本都是基于弱各向異性理論。

    1996年,TSVANKIN[3]基于THOMSEN的弱各向異性理論推導(dǎo)得到P-SV波VTI介質(zhì)精確相速度頻散方程。2006年,ZHOU等[4]從TSVANKIN相速度頻散方程出發(fā),推導(dǎo)得到TTI介質(zhì)qP波波動(dòng)方程,該方程雖然采用了令橫波速度為0的TTI聲波近似,但仍然是一個(gè)P-SV波耦合的波動(dòng)方程,因?yàn)楦飨虍愋圆▓?chǎng)傳播過程中偽SV波仍然存在[5],而且在角度劇變(即各向異性的對(duì)稱軸在空間變化較快)情況下波場(chǎng)傳播存在嚴(yán)重的不穩(wěn)定性問題[6]。為了解決TTI介質(zhì)聲波近似各向異性波場(chǎng)傳播存在的不穩(wěn)定性問題,F(xiàn)LETCHER等[7]引入適當(dāng)?shù)臋M波分量來減小角度劇變時(shí)產(chǎn)生的不穩(wěn)定性,但是同時(shí)也帶來了較強(qiáng)的偽SV波假象。DUVENECK等[8-9]從Hooke定律和運(yùn)動(dòng)學(xué)方程出發(fā)推導(dǎo)得到新的TTI介質(zhì)方程;ZHANG等[10-11]從VTI介質(zhì)彈性波方程出發(fā)采用自共軛旋轉(zhuǎn)方式推導(dǎo)出一個(gè)自共軛3DTTI波動(dòng)方程。但上述方程仍然存在偽SV波假象且在角度劇變時(shí)不穩(wěn)定。FOWLER等[12]從特征值矩陣分析的角度對(duì)各種TTI介質(zhì)耦合方程的精度和穩(wěn)定性進(jìn)行了綜合分析。BUBE等[13-14]討論了二階TI介質(zhì)波動(dòng)方程在TTI介質(zhì)聲波近似情況下不穩(wěn)定性產(chǎn)生的原因并推導(dǎo)了新的聲波和彈性波TI介質(zhì)一階波動(dòng)方程。

    TTI介質(zhì)正演時(shí),吸收邊界條件至關(guān)重要,因?yàn)槲者吔鐥l件極大地影響正演模擬效果。地下介質(zhì)被視為無限半空間,而計(jì)算區(qū)域有限,因此在人為截取的有限空間內(nèi)求解波動(dòng)方程會(huì)導(dǎo)致很強(qiáng)的人工邊界反射,必須采用吸收邊界條件對(duì)人工邊界反射進(jìn)行有效吸收。目前常用的吸收邊界條件包括單程波旁軸近似吸收邊界條件[15]和衰減吸收邊界條件[16-17]兩大類。CLAYTON等[15]提出的單程波旁軸近似吸收邊界條件應(yīng)用廣泛,在小角度反射時(shí)吸收效果較好,但對(duì)大角度入射的波吸收效果較差。CERJAN等[16]提出了添加阻尼層的海綿吸收邊界條件,該邊界條件不受波動(dòng)方程形式的限制,很容易實(shí)現(xiàn),但是邊界吸收的效果仍較差。BORDING[18]對(duì)CERJAN等[16]提出的阻尼系數(shù)進(jìn)行了優(yōu)化,計(jì)算得到優(yōu)化的海綿吸收系數(shù),但是計(jì)算效率相對(duì)CERJAN等的方法有所降低。BERENGER[17]最早在電磁學(xué)領(lǐng)域提出了完全匹配層(PML)吸收邊界條件,能達(dá)到完美的邊界吸收效果。CHEW等[19]將PML吸收邊界條件用于地震波數(shù)值模擬中,取得了很好的應(yīng)用效果。BéCACHE等[20]對(duì)各向異性介質(zhì)的PML吸收邊界條件進(jìn)行了深入研究,指出各向異性介質(zhì)應(yīng)用PML吸收邊界條件時(shí)存在不穩(wěn)定性。PML吸收邊界條件是目前應(yīng)用效果最好的吸收邊界條件,代表吸收邊界技術(shù)研究的前沿發(fā)展方向,但是其公式推導(dǎo)和編程實(shí)現(xiàn)復(fù)雜,需要對(duì)不同的波動(dòng)方程推導(dǎo)不同的PML方程形式。

    作為PML吸收邊界條件技術(shù)研究的一個(gè)重要方面,分裂式PML(SPML)吸收邊界條件已經(jīng)得到廣泛應(yīng)用。COLLINO等[21]提出二維彈性波介質(zhì)一階速度 應(yīng)力方程的SPML算法;裴正林[22-23]將該算法推廣到三維彈性波介質(zhì),實(shí)現(xiàn)了基于SPML吸收邊界條件的三維彈性波介質(zhì)波動(dòng)方程交錯(cuò)網(wǎng)格高階有限差分法數(shù)值模擬。

    本文研究的TTI介質(zhì)時(shí)間域SPML吸收邊界條件能很好地吸收邊界反射,實(shí)現(xiàn)TTI介質(zhì)高精度數(shù)值模擬。

    1 方法原理

    1.1 TTI介質(zhì)耦合qP波波動(dòng)方程

    P-SV波VTI介質(zhì)精確相速度頻散方程為[3]:

    式中:f=1-v2SZ/v2PZ,vSZ表示橫波速度,vPZ表示縱波速度;vP(θph)表示相速度,θph表示相速度傳播角度;ε和δ表示Thomsen各向異性參數(shù)[2]。

    在推導(dǎo)過程中直接引入TTI介質(zhì)聲波近似(TTI介質(zhì)聲波近似情況下直接令橫波速度vSZ=0)[24],此時(shí)推導(dǎo)得到的2DTTI介質(zhì)方程為[4]:

    式中:θ為對(duì)稱軸傾角;p為地震波場(chǎng);q為輔助波場(chǎng);H1和H2為輔助變量,其3DTTI方程形式為:

    式中:φ為三維方位角。當(dāng)θ=0時(shí),(3)式退化為3D VTI方程(垂直對(duì)稱軸的TI方程);當(dāng)θ=90°時(shí),(3)式退化為3DHTI方程(水平對(duì)稱軸的TI方程);當(dāng)三維方位角φ=0時(shí),(3)式簡化為2DTTI方程。VTI介質(zhì)和HTI介質(zhì)都可以看作是TTI介質(zhì)的某種特殊情形。

    本文將ZHOU等[4]、DUVENECK等[8-9]和ZHANG等[10-11]直接采用vSZ=0的TTI介質(zhì)聲波近似思想推導(dǎo)得到的TTI介質(zhì)方程歸納為基于TTI介質(zhì)聲波近似的TTI介質(zhì)方程。但是這種方程在實(shí)際計(jì)算中往往會(huì)產(chǎn)生數(shù)值不穩(wěn)定問題,尤其是在角度劇變的情況下[6]。FLETCHER等[7]完全放棄了采用TTI介質(zhì)聲波近似的思路,提出在方程推導(dǎo)中保留橫波分量。引入橫波分量的方式能有效消除TTI介質(zhì)聲波近似產(chǎn)生的不穩(wěn)定問題,F(xiàn)LETCHER等[7]引入橫波分量推導(dǎo)得到的2D TTI介質(zhì)縱橫波耦合qP波波動(dòng)方程為:

    TTI介質(zhì)波場(chǎng)模擬首先要滿足物理穩(wěn)定性條件,即P-SV波VTI介質(zhì)精確相速度頻散方程((1)式)中的相速度的平方不能小于0[26]:

    由(5)式推導(dǎo)得到TI介質(zhì)波場(chǎng)模擬中的物理穩(wěn)定性條件為:

    在TTI介質(zhì)聲波近似情況下,vSZ=0,則f=1,由(6)式推導(dǎo)得到TTI介質(zhì)聲波近似情況下的物理穩(wěn)定性條件為:

    在變速介質(zhì)情況下,由于對(duì)變系數(shù)求導(dǎo)會(huì)帶來微分算子的附加項(xiàng),影響算子的穩(wěn)定性,因此,我們借助對(duì)稱算子理論,將原始偏微分方程適當(dāng)對(duì)稱化,得到與常系數(shù)偏微分方程一樣穩(wěn)定的微分算子[11,27-28]。

    經(jīng)實(shí)際測(cè)試,TTI介質(zhì)波場(chǎng)傳播不穩(wěn)定的情況主要由角度劇變產(chǎn)生,包括對(duì)稱軸傾角θ或者三維方位角φ快速變化的區(qū)域[6]。本文設(shè)計(jì)了一個(gè)楔形2DTTI介質(zhì)模型用于算法穩(wěn)定性測(cè)試(圖1),模型的主要特點(diǎn)是角度劇變(圖1d),因此能很好地測(cè)試各向異性波動(dòng)方程算法的穩(wěn)定性。模型橫向網(wǎng)格數(shù)和縱向網(wǎng)格數(shù)都為300,橫向網(wǎng)格間距和縱向網(wǎng)格間距都為10m。

    圖1 楔形2DTTI介質(zhì)模型各向異性參數(shù)變化情況

    數(shù)值模擬時(shí),雷克震源子波主頻取25Hz,傳播時(shí)間為3s,時(shí)間采樣間隔為0.000 5s,吸收邊界條件采用優(yōu)化海綿吸收邊界條件[18],分別采用TTI介質(zhì)聲波近似方程(ZHOU等[4]方程)和引入橫波分量的TTI介質(zhì)縱橫波耦合方程(FLETCHER等[7]方程)來進(jìn)行數(shù)值模擬。

    從t=1s時(shí)刻的波場(chǎng)快照(圖2a,圖2b)可以看出,ZHOU等[4]方程基本上還能保持穩(wěn)定傳播,F(xiàn)LETCHER等[7]方程由于引入了橫波分量,產(chǎn)生了偽橫波(SV波)噪聲(圖2b中紅色箭頭所示)。而從t=1.5s時(shí)刻的波場(chǎng)快照(圖2c,圖2d)可以看出,ZHOU等[4]方程由于楔形界面處角度劇變已經(jīng)產(chǎn)生不穩(wěn)定現(xiàn)象,F(xiàn)LETCHER等[7]方程由于引入了橫波分量,能減小不穩(wěn)定情況的產(chǎn)生,仍然能保持穩(wěn)定傳播。因此FLETCHER等[7]各向異性波動(dòng)方程是一種穩(wěn)定的TTI介質(zhì)波動(dòng)方程,我們討論的二維和三維TTI介質(zhì)正演模擬將采用該方程。

    本文對(duì)常見的TTI介質(zhì)波動(dòng)方程進(jìn)行了歸納,認(rèn)為基于TTI介質(zhì)精確相速度頻散方程推導(dǎo)得到的TTI介質(zhì)波動(dòng)方程可以分為兩類:一類為TTI介質(zhì)聲波近似方程;另一類為TTI介質(zhì)縱橫波耦合方程。其中TTI介質(zhì)聲波近似方程可以視為TTI介質(zhì)縱橫波耦合方程在橫波速度為0時(shí)的某種特例。從TTI介質(zhì)波場(chǎng)傳播穩(wěn)定性的角度對(duì)這兩類方程進(jìn)行對(duì)比分析可知,引入橫波分量的TTI介質(zhì)縱橫 波耦合方程適用于TTI介質(zhì)。

    圖2 楔形2DTTI介質(zhì)模型波場(chǎng)快照

    1.2 2DTTI介質(zhì)時(shí)間域分裂式完全匹配層(SPML)波動(dòng)方程推導(dǎo)

    基于FLETCHER等[7]的二階2DTTI介質(zhì)縱橫波耦合波動(dòng)方程((4)式)來推導(dǎo)2DTTI介質(zhì)時(shí)間域SPML波動(dòng)方程。TTI介質(zhì)波動(dòng)方程的特點(diǎn)是含有比較復(fù)雜的交叉導(dǎo)數(shù)項(xiàng),如何有效處理TTI介質(zhì)方程中的交叉導(dǎo)數(shù)項(xiàng)是PML吸收邊界條件研究的難點(diǎn)。本文的思路是針對(duì)TTI介質(zhì)波動(dòng)方程的特點(diǎn)對(duì)二階方程中的交叉導(dǎo)數(shù)項(xiàng)進(jìn)行分裂,得到其一階方程的形式,進(jìn)而推導(dǎo)得到TTI介質(zhì)一階SPML波動(dòng)方程。

    引入4個(gè)輔助波場(chǎng)項(xiàng)ψx,ψz,Ωx,Ωz,從FLETCHER等[7]的二階2DTTI介質(zhì)波動(dòng)方程出發(fā),推導(dǎo)得到一階2DTTI介質(zhì)波動(dòng)方程的形式,對(duì)方程進(jìn)行傅里葉變換并引入X和Z方向的頻率域PML拉伸函數(shù),推導(dǎo)得到2DTTI介質(zhì)頻率域SPML波動(dòng)方程,將此方程反傅里葉變換到時(shí)間域,即可得到2DTTI介質(zhì)時(shí)間域SPML波動(dòng)方程:

    其中,輔助變量形式為:

    其中,引入了PML衰減因子dx,dz。本文采用的PML衰減因子表達(dá)式為[21]:

    式中:Lpml為PML邊界層的厚度;l為計(jì)算點(diǎn)距PML邊界的距離;Rcoeff為理論反射系數(shù),本文取Rcoeff=0.000 1。

    一系列公式推導(dǎo)研究可以發(fā)現(xiàn),所有2DTTI介質(zhì)的SPML方程都具有相同的形式,差分離散形式也相同,只是PML輔助變量表達(dá)形式不一致,其具體形式根據(jù)不同的TTI方程形式而定。本文只討論二維TTI介質(zhì)情形,但是這種方法也可以擴(kuò)展到三維TTI介質(zhì)(詳細(xì)推導(dǎo)見附錄A),在此不再展開討論。

    1.3 2DTTI介質(zhì)SPML波動(dòng)方程數(shù)值實(shí)現(xiàn)

    采用高階交錯(cuò)網(wǎng)格有限差分法來實(shí)現(xiàn)1.2節(jié)推導(dǎo)的2DTTI介質(zhì)SPML波動(dòng)方程。對(duì)(8)式進(jìn)行高階交錯(cuò)網(wǎng)格有限差分離散,具體實(shí)現(xiàn)公式如下:

    從二維SPML數(shù)值實(shí)現(xiàn)示意圖(圖3)可以看出,二維PML區(qū)域分為4個(gè)邊和4個(gè)角共8塊區(qū)域,具體實(shí)現(xiàn)時(shí),我們需要分別考慮各PML區(qū)域波場(chǎng)的衰減。理論上,無吸收邊界條件具有最強(qiáng)的邊界反射,而PML吸收層數(shù)越多,邊界反射吸收效果越好。但是增加的PML吸收層數(shù)也帶來了更多的計(jì)算量,從而在一定程度上降低了計(jì)算效率,因此在實(shí)際應(yīng)用中需要在PML吸收層數(shù)與計(jì)算效率之間尋求平衡。

    圖3 二維SPML數(shù)值實(shí)現(xiàn)示意圖解

    2 數(shù)值計(jì)算實(shí)例

    2.1 簡單均勻2DTTI介質(zhì)模型

    為了驗(yàn)證本文提出的TTI介質(zhì)SPML吸收邊界條件的有效性,首先采用簡單均勻2DTTI介質(zhì)模型進(jìn)行SPML點(diǎn)源響應(yīng)測(cè)試。模型橫向網(wǎng)格數(shù)和縱向網(wǎng)格數(shù)都為256,橫向網(wǎng)格間距和縱向網(wǎng)格間距都為10m,速度為3 000m/s。TTI介質(zhì)的各向異性參數(shù)ε=0.3,δ=0.1,對(duì)稱軸傾角為45°。

    均勻TTI介質(zhì)模型的單道記錄能很好地測(cè)試不同PML邊界層條件下的邊界反射吸收效果。數(shù)值模擬時(shí),雷克震源子波主頻取25Hz,震源位置設(shè)置在模型的中心,單檢波點(diǎn)取在震源位置左側(cè)50個(gè)網(wǎng)格點(diǎn)處,傳播時(shí)間為1.5s,時(shí)間采樣間隔為0.001 5s,采用12階高階交錯(cuò)網(wǎng)格有限差分格式進(jìn)行模擬。采用SPML吸收邊界條件,PML邊界吸收層數(shù)分別設(shè)為10層、20層和30層。取3種情形的單道記錄進(jìn)行對(duì)比(圖4)。由圖4可見,PML吸收層數(shù)為10層時(shí)邊界反射吸收效果不理想;而PML吸收層數(shù)為20層時(shí)已經(jīng)取得了比較好的邊界反射吸收效果(圖4a),但從局部放大圖可以看出,仍然存在著微弱的邊界反射;而PML吸收層數(shù)為30層時(shí)邊界反射導(dǎo)致的波場(chǎng)誤差降到僅為10-5,此時(shí)邊界反射基本被消除(圖4b)。從PML吸收層數(shù)為20層時(shí)SPML波動(dòng)方程在t=0.72s和t=0.90s的波場(chǎng)快照(圖5)也可以看到微弱的邊界反射。因此本例選定PML邊界吸收層數(shù)為30層。圖6a到圖6f分別給出了SPML波動(dòng)方程在PML吸收層數(shù)為30層時(shí)各個(gè)時(shí)刻的波場(chǎng)快照。由圖6可見,SPML對(duì)邊界反射吸收效果明顯。

    2.2 復(fù)雜BP2007 2DTTI介質(zhì)海洋標(biāo)準(zhǔn)模型

    以復(fù)雜BP2007 2DTTI介質(zhì)海洋標(biāo)準(zhǔn)模型為例,對(duì)比TTI介質(zhì)SPML吸收邊界條件與優(yōu)化海綿吸收邊界條件對(duì)邊界反射的吸收效果[18]。本文數(shù)值模擬時(shí)采用截取的部分BP2007 2DTTI介質(zhì)海洋標(biāo)準(zhǔn)模型(圖7),該模型主要模擬了近海岸的地質(zhì)構(gòu)造,特點(diǎn)是左部有一個(gè)被TTI各向異性地層包圍的高速各向同性鹽丘,模型參數(shù)包括縱波速度vPZ,Thomsen縱波各向異性參數(shù)ε,Thomsen變異系數(shù)δ,對(duì)稱軸角度θ,模型橫向網(wǎng)格數(shù)nx=601,縱向網(wǎng)格數(shù)nz=361,橫向網(wǎng)格間距和縱向網(wǎng)格間距都為6.25m,其中海水層為各向同性介質(zhì)。數(shù)值模擬時(shí),雷克震源子波主頻取25Hz,采用SPML吸收邊界條件,PML邊界設(shè)為30層,傳播時(shí)間為2s,時(shí)間采樣間隔為0.000 4s,采用12階高階交錯(cuò)網(wǎng)格有限差分格式進(jìn)行模擬。因?yàn)檎鹪粗糜诟飨蛲缘暮K畬又校琓TI數(shù)值模擬時(shí)不會(huì)出現(xiàn)影響qP波場(chǎng)傳播的偽橫波(SV波)噪聲,此種情況下不需要采取在震源處添加各向同性薄層消除偽SV波噪聲的策略[5,29]。從BP2007 2DTTI介質(zhì)海洋標(biāo)準(zhǔn)模型的正演單炮記錄(圖8)來看,SPML消除邊界反射效果良好(圖8a),而優(yōu)化海綿吸收邊界條件仍然有著比較強(qiáng)的邊界反射(圖8b中紅色箭頭標(biāo)示),SPML吸收邊界條件明顯優(yōu)于優(yōu)化海綿吸收邊界條件的邊界反射吸收效果。

    圖4 均勻2DTTI介質(zhì)模型不同PML吸收層數(shù)單道記錄

    圖5 均勻2DTTI介質(zhì)模型采用2DTTI介質(zhì)SPML波動(dòng)方程在PML吸收層數(shù)為20層時(shí)不同時(shí)刻的正演波場(chǎng)快照

    圖6 均勻2DTTI介質(zhì)模型采用2DTTI介質(zhì)SPML波動(dòng)方程在PML吸收層數(shù)為30層時(shí)不同時(shí)刻的正演波場(chǎng)快照

    圖7 截取的部分BP2007 2DTTI介質(zhì)海洋標(biāo)準(zhǔn)模型各向異性參數(shù)變化情況

    圖8 BP2007 2DTTI介質(zhì)海洋標(biāo)準(zhǔn)模型采用不同邊界條件正演的單炮記錄

    3 結(jié)束語

    本文對(duì)常見的幾種TTI介質(zhì)波動(dòng)方程進(jìn)行了歸納,并從波場(chǎng)傳播穩(wěn)定性的角度進(jìn)行分析,結(jié)果表明引入橫波分量的TTI介質(zhì)縱橫波耦合波動(dòng)方程適用于TTI介質(zhì)?;赥TI介質(zhì)縱橫波耦合二階波動(dòng)方程,推導(dǎo)得到一階方程形式的TTI介質(zhì)SPML波動(dòng)方程,研究發(fā)現(xiàn),所有2DTTI介質(zhì)的SPML方程都具有相同的形式,差分離散形式也相同,只是PML輔助變量表達(dá)形式不一致,其根據(jù)不同的TTI方程形式而定。數(shù)值模擬結(jié)果表明,本文研究的SPML吸收邊界條件能達(dá)到很好的人工邊界反射吸收效果。本文方法可以進(jìn)一步推廣到三維TTI介質(zhì)的高精度模擬中。至于2DVTI或3DVTI介質(zhì)情形,由于對(duì)稱軸角度為0,波場(chǎng)傳播較穩(wěn)定,此時(shí)適合采用噪聲相對(duì)較少的VTI介質(zhì)聲波近似方程,因?yàn)門TI介質(zhì)縱橫波耦合方程雖然由于引入了橫波分量解決了地震波在各向異性介質(zhì)中傳播的不穩(wěn)定性問題,但是也帶來了更多的橫波噪聲。VTI介質(zhì)波動(dòng)方程的特點(diǎn)是不含比較復(fù)雜的交叉導(dǎo)數(shù)項(xiàng),可以采用計(jì)算效率更高的非SPML吸收邊界條件來實(shí)現(xiàn)。

    [1] WARNER M,RATCLIFE A,NANGOO T,et al.Ani-sotropic 3Dfull-waveform inversion[J].Geophysics,2013,78(2):R59-R80

    [2] THOMSEN L.Weak elastic anisotropy[J].Geophysics,1986,51(10):1954-1966

    [3] TSVANKIN I.P-wave signatures and notation for transversely isotropic media:an overview[J].Geophysics,1996,61(2):467-483

    [4] ZHOU H B,ZHANG G Q,BLOOR R.An anisotropic acoustic wave equation for modeling and migration in 2DTTI media[J].Expanded Abstracts of 76thAnnual Internat SEG Mtg,2006:194-198

    [5] GRECHKA V,ZHANG L,RECTOR J W.Shear waves in acoustic anisotropic media[J].Geophysics,2004,69(2):576-582

    [6] BAKKER P M,DUVENECK E.Stability analysis for acoustic wave propagation in tilted TI media by finite differences[J].Geophysical Journal International,2011,185(2):911-921

    [7] FLETCHER R P,DU X,F(xiàn)OWLER P J.Reverse time migration in tilted transversely isotropic(TTI)media[J].Geophysics,2009,74(6):179-187

    [8] DUVENECK E,MILCIK P,BAKKER P M.Acoustic VTI wave equations and their application for anisotropic reverse-time migration[J].Expanded Abstracts of 78thAnnual Internat SEG Mtg,2008:2186-2190

    [9] DUVENECK E,BAKKER P M.Stable P-wave modeling for reverse-time migration in tilted TI media[J].Geophysics,2011,76(2):65-75

    [10] ZHANG Y,ZHANG H Z.A stable TTI reverse time migration and its implementation[J].Expanded Abstracts of 79thAnnual Internat SEG Mtg,2009:2794-2798

    [11] ZHANG Y,ZHANG H Z,ZHANG G Q.A stable TTI reverse time migration and its implementation[J].Geophysics,2011,76(3):WA3-WA11

    [12] FOWLER P J,DU X,F(xiàn)LETCHER R P.Coupled equations for reverse time migration in transversely isotropic media[J].Geophysics,2010,75(1):S11-S22

    [13] BUBE K P,NEMETH T,STEFANI J P,et al.Firstorder systems for elastic and acoustic variable-tilt TI media[J].Geophysics,2012,77(5):T157-T170

    [14] BUBE K P,NEMETH T,STEFANI J P,et al.On the instability in second-order systems for acoustic VTI and TTI media[J].Geophysics,2012,77(5):T171-T186

    [15] CLAYTON R,ENQUIST B.Absorbing boundary conditions for acoustic and elastic wave equations[J].Bulletin of the Seismological Society of America,1977,67(6):1529-1540

    [16] CERJAN C,KOSLOFF D,KOSLOFF R,et al.A nonreflecting boundary condition for discrete acoustic and elastic wave equations[J].Geophysics,1985,50(4):705-708

    [17] BERENGER J P.A perfectly matched layer for the absorption of electromagnetic waves[J].Journal of Computational Physics,1994,114(2):185-200

    [18] BORDING R P.Finite difference modeling-nearly optimal sponge boundary conditions[J].Expanded Abstracts of 74thAnnual Internat SEG Mtg,2004:1921-1924

    [19] CHEW W C,LIU Q H.Perfectly matched layers for elastodynamics:a new absorbing boundary condition[J].Journal of Computational Acoustics,1996,4(4):341-359

    [20] BéCACHE E,F(xiàn)AUQUEUX S,JOLY P.Stability of perfectly matched layers,group velocities and anisotropic waves[J].Journal of Computational Physics,2003,188(2):399-433

    [21] COLLINO F,TSOGKA C.Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[J].Geophysics,2001,66(1):294-307

    [22] 裴正林.三維各向異性介質(zhì)中彈性波方程交錯(cuò)網(wǎng)格高階有限差分法數(shù)值模擬[J].石油大學(xué)學(xué)報(bào)(自然科學(xué)版),2004,28(5):23-29

    PEI Z L.Three-dimensional numerical simulation of elastic wave propagation in 3-D anisotropic media withstaggered-grid high-order difference method[J].Journal of the University of Petroleum,2004,28(5):23-29

    [23] 裴正林.三維各向同性介質(zhì)彈性波方程交錯(cuò)網(wǎng)格高階有限差分法模擬[J].石油物探,2005,44(4):308-315

    PEI Z L.Numerical simulation of elastic wave equation in 3-D isotropic media with staggered-grid high-order difference method[J].Geophysical Prospecting for Petroleum,2005,44(4):308-315

    [24] ALKHALIFAH T.Acoustic approximations for processing in transversely isotropic media[J].Geophysics,1998,63(2):623-631

    [25] 李博,李敏,劉紅偉,等.TTI介質(zhì)有限差分逆時(shí)偏移的穩(wěn)定性探討[J].地球物理學(xué)報(bào),2012,55(4):1366-1375

    LI B,LI M,LIU H W,et al.Stability of reverse time migration in TTI media[J].Chinese Journal of Geophysics,2012,55(4):1366-1375

    [26] 李博.地震偏移成像和GPU超算技術(shù)[D].北京:中國科學(xué)院大學(xué),2011

    LI B.Seismic migration and GPU supercomputing technology[D].Beijing:University of Chinese Academy of Sciences,2011

    [27] CLAERBOUT J F.Imaging the earth’s interior[M].Palo Alto,California:Blackwell Scientific Publications,1985:1-398

    [28] 劉洪,王秀閩,曾銳,等.單程波算子積分解的象征表示[J].地球物理學(xué)進(jìn)展,2007,22(2):463-471

    LIU H,WANG X M,ZENG R,et al.Symbol description to integral solution of one-way wave operator[J].Progress in Geophysics,2007,22(2):463-471

    [29] 張衡,劉洪,唐祥德,等.基于平均導(dǎo)數(shù)優(yōu)化方法的VTI介質(zhì)頻率空間域正演[J].地球物理學(xué)報(bào),2015,58(9):3306-3316

    ZHANG H,LIU H,TANG X D,et al.Forward modeling of VTI media in the frequency-space domain based on an average-derivative optimal method[J].Chinese Journal of Geophysics,2015,58(9):3306-3316

    附錄A 3DTTI介質(zhì)時(shí)間域SPML波動(dòng)方程

    3DTTI介質(zhì)時(shí)間域SPML波動(dòng)方程的推導(dǎo)基于FLETCHER等[7]的二階3DTTI介質(zhì)縱橫波耦合波動(dòng)方程((4)式)。與二維TTI介質(zhì)相比,三維TTI介質(zhì)時(shí)間域SPML波動(dòng)方程的推導(dǎo)首先需要從FLETCHER等[7]的二階3DTTI介質(zhì)波動(dòng)方程出發(fā),推導(dǎo)得到一階3DTTI介質(zhì)波動(dòng)方程的形式,然后對(duì)方程進(jìn)行傅里葉變換并引入x,y和z 3個(gè)方向的頻率域PML拉伸函數(shù),將推導(dǎo)得到的3DTTI介質(zhì)頻率域SPML波動(dòng)方程反傅里葉變換,即可得到3DTTI介質(zhì)時(shí)間域SPML波動(dòng)方程:

    推導(dǎo)的輔助變量形式為:

    PML波動(dòng)方程推導(dǎo)過程中引入的頻率域PML拉伸函數(shù)ξx,ξy,ξz表示為:

    引入的PML衰減因子dx,dy,dz表達(dá)式為:

    式中:Lpml是PML邊界層的厚度;l是計(jì)算點(diǎn)距PML邊界的距離;Rcoeff是理論反射系數(shù)值,取0.000 1。

    對(duì)一系列公式推導(dǎo)研究可以發(fā)現(xiàn),所有3DTTI介質(zhì)方程的SPML方程都具有相同的形式,差分離散形式也相同,只是PML輔助變量表達(dá)形式不一致,其具體形式根據(jù)不同的TTI方程形式而定。

    (編輯:陳 杰)

    The research on split PML absorbing boundary conditions of acoustic equation for TTI media

    ZHANG Heng1,LIU Hong2,LI Bo3,DING Renwei4,LI Liqing1,LI Fuyuan1

    (1.MLR Key Laboratory of Marine Mineral Resources,Guangzhou Marine Geological Survey,China Geological Survey,Guangzhou510075,China;2.Key Laboratory of Petroleum Resources Research,Institute of Geology and Geophysics,Chinese Academy of Sciences,Beijing100029,China;3.Sinopec Geophysical Research Institute,Nanjing211103,China;4.College of Earth Science and Engineering,Shandong University of Science and Technology,Qingdao 266590,China)

    We study split perfectly matched layer(SPML)absorbing boundary conditions for TTI media based on the characteristics of TTI media acoustic wave equation.We firstly summarize several common TTI media acoustic wave equations and compare these wave equations in detail from the perspective of wavefield propagation stability for TTI media.We conclude that the TTI acoustic approximation equation introduced SV-wave component is applicable for TTI media.Then we derive the TTI media first-order acoustic wave equation based on the SPML boundary condition from the TTI P-wave and SV-wave coupled second-order wave equation.Afterwards we give the concrete implementation method for high-order staggered-grid finite-difference algorithm.The numerical modeling examples show that the SPML boundary condition can absorb the artificial boundary reflection very well.The boundary reflection absorbing effect of the proposed SPML algorithm is much better compared with the optimal sponge boundary condition.

    P631

    A

    1000-1441(2017)03-0349-13

    10.3969/j.issn.1000-1441.2017.03.005

    張衡,劉洪,李博,等.TTI介質(zhì)聲波方程分裂式PML吸收邊界條件研究[J].石油物探,2017,56(3):349-361

    ZHANG Heng,LIU Hong,LI Bo,et al.The research on split PML absorbing boundary conditions of acoustic equation for TTI media[J].Geophysical Prospecting for Petroleum,2017,56(3):349-361

    2016-02-14;改回日期:2016-04-29。

    張衡(1987—),男,博士,工程師,主要從事復(fù)雜各向異性介質(zhì)地震波場(chǎng)正演數(shù)值模擬和全波形反演研究工作。

    國家自然科學(xué)基金項(xiàng)目(41604110)、國土資源部海底礦產(chǎn)資源重點(diǎn)實(shí)驗(yàn)室開放基金項(xiàng)目(KLMMR-2015-A-14)、國家高技術(shù)研究發(fā)展計(jì)劃(863計(jì)劃)項(xiàng)目(2013AA092501)和國家科技重大專項(xiàng)(2011ZX05003-003)聯(lián)合資助。

    This research is financially supported by the National Natural Science Foundation of China(Grant No.41604110),the Open Fund of MLR Key Laboratory of Marine Mineral Resources(Grant No.KLMMR-2015-A-14),the National High-tech R &D Program of China(863Program)(Grant No.2013AA092501)and the National Science and Technology Major Project of China(Grant No.2011ZX05003-003).

    猜你喜歡
    波場(chǎng)橫波邊界條件
    橫波技術(shù)在工程物探中的應(yīng)用分析
    一類帶有Stieltjes積分邊界條件的分?jǐn)?shù)階微分方程邊值問題正解
    帶有積分邊界條件的奇異攝動(dòng)邊值問題的漸近解
    彈性波波場(chǎng)分離方法對(duì)比及其在逆時(shí)偏移成像中的應(yīng)用
    交錯(cuò)網(wǎng)格與旋轉(zhuǎn)交錯(cuò)網(wǎng)格對(duì)VTI介質(zhì)波場(chǎng)分離的影響分析
    基于Hilbert變換的全波場(chǎng)分離逆時(shí)偏移成像
    旋轉(zhuǎn)交錯(cuò)網(wǎng)格VTI介質(zhì)波場(chǎng)模擬與波場(chǎng)分解
    帶Robin邊界條件的2維隨機(jī)Ginzburg-Landau方程的吸引子
    揚(yáng)眉一顧,妖嬈橫波處
    橫波一顧,傲殺人間萬戶侯
    火花(2015年1期)2015-02-27 07:40:24
    777米奇影视久久| 免费av中文字幕在线| 久久久成人免费电影| 亚洲熟女精品中文字幕| 亚洲欧美成人综合另类久久久| 美女国产视频在线观看| 一级毛片 在线播放| 国产精品秋霞免费鲁丝片| 亚洲成人中文字幕在线播放| freevideosex欧美| 亚洲av中文av极速乱| 中文字幕av成人在线电影| 欧美精品一区二区免费开放| 国产极品天堂在线| freevideosex欧美| 七月丁香在线播放| 又粗又硬又长又爽又黄的视频| 有码 亚洲区| 国产亚洲最大av| 色视频在线一区二区三区| 亚洲,一卡二卡三卡| 欧美日韩国产mv在线观看视频 | 婷婷色av中文字幕| 国产精品久久久久久久久免| 国产中年淑女户外野战色| 一区二区三区乱码不卡18| 亚洲av免费高清在线观看| 亚洲精品日本国产第一区| 一本一本综合久久| 简卡轻食公司| 最近2019中文字幕mv第一页| 国产视频内射| 干丝袜人妻中文字幕| 久久精品国产亚洲av天美| 日韩,欧美,国产一区二区三区| 国产精品国产三级专区第一集| 麻豆国产97在线/欧美| 99久国产av精品国产电影| 新久久久久国产一级毛片| 女性生殖器流出的白浆| 免费不卡的大黄色大毛片视频在线观看| 欧美成人a在线观看| 晚上一个人看的免费电影| 亚洲,欧美,日韩| 2022亚洲国产成人精品| 激情五月婷婷亚洲| 欧美3d第一页| 黄色配什么色好看| 99久久精品热视频| 免费观看的影片在线观看| 久久精品国产a三级三级三级| 黄色配什么色好看| 久久久成人免费电影| 亚洲高清免费不卡视频| 观看免费一级毛片| av在线蜜桃| 国产亚洲欧美精品永久| 97精品久久久久久久久久精品| 男女无遮挡免费网站观看| 久久人妻熟女aⅴ| 少妇人妻久久综合中文| 国产精品熟女久久久久浪| 黄色欧美视频在线观看| 亚洲一区二区三区欧美精品| 久久影院123| 三级国产精品片| 欧美高清成人免费视频www| 男人舔奶头视频| 久久久精品94久久精品| 婷婷色av中文字幕| 色网站视频免费| 国产精品一及| 亚洲va在线va天堂va国产| 亚洲不卡免费看| h日本视频在线播放| 亚洲欧洲国产日韩| 久久国内精品自在自线图片| 人人妻人人看人人澡| 男女免费视频国产| 色哟哟·www| 插阴视频在线观看视频| 精品亚洲成a人片在线观看 | 欧美3d第一页| 老司机影院毛片| 亚洲av.av天堂| 久久午夜福利片| 免费黄网站久久成人精品| 黄色一级大片看看| 五月玫瑰六月丁香| 最近2019中文字幕mv第一页| 亚洲最大成人中文| 哪个播放器可以免费观看大片| 亚洲,一卡二卡三卡| 欧美高清性xxxxhd video| av天堂中文字幕网| 一边亲一边摸免费视频| 建设人人有责人人尽责人人享有的 | 乱码一卡2卡4卡精品| 色网站视频免费| 亚洲精品色激情综合| 五月天丁香电影| 精品一品国产午夜福利视频| 国产精品偷伦视频观看了| 日本av手机在线免费观看| 晚上一个人看的免费电影| 国产黄片美女视频| 大片电影免费在线观看免费| 欧美丝袜亚洲另类| 国产一区二区三区综合在线观看 | 老师上课跳d突然被开到最大视频| 欧美最新免费一区二区三区| 国产精品免费大片| 特大巨黑吊av在线直播| 亚洲国产高清在线一区二区三| 大码成人一级视频| 久久精品夜色国产| 久久久久久久亚洲中文字幕| 91久久精品电影网| 日本猛色少妇xxxxx猛交久久| 国产爱豆传媒在线观看| 777米奇影视久久| 卡戴珊不雅视频在线播放| 一本一本综合久久| 日韩大片免费观看网站| 国产黄片视频在线免费观看| 女的被弄到高潮叫床怎么办| 成人高潮视频无遮挡免费网站| 亚洲精品日韩在线中文字幕| 亚洲欧洲日产国产| 久久午夜福利片| 久热久热在线精品观看| 亚洲精品一二三| 成年美女黄网站色视频大全免费 | 伊人久久国产一区二区| 久久亚洲国产成人精品v| 亚洲av中文av极速乱| 伦理电影免费视频| 亚洲av电影在线观看一区二区三区| 国产亚洲精品久久久com| 亚洲精品一区蜜桃| 成人高潮视频无遮挡免费网站| 女人久久www免费人成看片| 中文字幕av成人在线电影| 成人国产麻豆网| 国产男人的电影天堂91| 日韩成人av中文字幕在线观看| 免费观看的影片在线观看| 三级国产精品片| 免费看日本二区| 直男gayav资源| 亚洲国产精品999| 大话2 男鬼变身卡| 国产精品秋霞免费鲁丝片| 偷拍熟女少妇极品色| 午夜老司机福利剧场| 国产精品久久久久久久久免| 亚洲人成网站在线观看播放| 国产精品久久久久久久电影| 一区二区三区精品91| 欧美性感艳星| 久久久a久久爽久久v久久| 女人久久www免费人成看片| av网站免费在线观看视频| 亚洲性久久影院| 国产精品一二三区在线看| 亚洲av成人精品一区久久| 国产深夜福利视频在线观看| 最后的刺客免费高清国语| 久久久久久久久久久丰满| 日韩在线高清观看一区二区三区| 在线观看av片永久免费下载| 免费高清在线观看视频在线观看| 成人国产麻豆网| 一级二级三级毛片免费看| 人妻少妇偷人精品九色| 欧美一级a爱片免费观看看| tube8黄色片| 国产男女超爽视频在线观看| av播播在线观看一区| 超碰av人人做人人爽久久| 亚洲av免费高清在线观看| 纯流量卡能插随身wifi吗| 狠狠精品人妻久久久久久综合| 欧美变态另类bdsm刘玥| 黑丝袜美女国产一区| 搡女人真爽免费视频火全软件| 最新中文字幕久久久久| 美女国产视频在线观看| 久久 成人 亚洲| 亚洲图色成人| 久久久久久久久久久免费av| 黄片无遮挡物在线观看| 视频区图区小说| av一本久久久久| 国产精品久久久久久av不卡| 熟女人妻精品中文字幕| 人人妻人人爽人人添夜夜欢视频 | av专区在线播放| 日本-黄色视频高清免费观看| 精品国产三级普通话版| 一级二级三级毛片免费看| 高清欧美精品videossex| 亚洲av综合色区一区| 亚洲性久久影院| 亚洲av国产av综合av卡| 男女边吃奶边做爰视频| 日韩亚洲欧美综合| 蜜桃亚洲精品一区二区三区| 人妻 亚洲 视频| 在线观看av片永久免费下载| 亚洲国产欧美在线一区| 毛片一级片免费看久久久久| 最近的中文字幕免费完整| 观看免费一级毛片| 美女脱内裤让男人舔精品视频| 99久久中文字幕三级久久日本| 国模一区二区三区四区视频| 18禁动态无遮挡网站| av.在线天堂| 亚洲av福利一区| 在线 av 中文字幕| 国产日韩欧美在线精品| 国产精品一区二区三区四区免费观看| 国产成人91sexporn| 日本爱情动作片www.在线观看| 欧美xxxx性猛交bbbb| 日本-黄色视频高清免费观看| 亚洲精品久久久久久婷婷小说| 国产成人91sexporn| 久久99热这里只频精品6学生| 99热6这里只有精品| 黄色一级大片看看| 久久久精品免费免费高清| 另类亚洲欧美激情| 欧美三级亚洲精品| 永久网站在线| 亚洲伊人久久精品综合| 嫩草影院入口| 日日啪夜夜爽| 人妻 亚洲 视频| 久久精品人妻少妇| 亚洲伊人久久精品综合| 久久精品久久精品一区二区三区| 中文字幕免费在线视频6| 国产在线视频一区二区| 国产精品无大码| 亚洲精品国产色婷婷电影| 久久6这里有精品| 亚洲成人一二三区av| 欧美激情国产日韩精品一区| 国产精品成人在线| 午夜老司机福利剧场| 国产精品人妻久久久久久| 亚洲av成人精品一二三区| 熟女人妻精品中文字幕| 亚洲中文av在线| 亚洲怡红院男人天堂| 久久97久久精品| 美女中出高潮动态图| 亚洲欧洲国产日韩| 欧美日韩亚洲高清精品| 一区二区三区精品91| 中文字幕精品免费在线观看视频 | 日本爱情动作片www.在线观看| h视频一区二区三区| 日日撸夜夜添| 久久精品国产a三级三级三级| 一区二区三区精品91| 老司机影院成人| 网址你懂的国产日韩在线| 大片电影免费在线观看免费| 99热这里只有是精品50| 日韩伦理黄色片| 久久久久国产精品人妻一区二区| 免费观看av网站的网址| 免费观看无遮挡的男女| 尾随美女入室| 一级毛片电影观看| 联通29元200g的流量卡| 一本—道久久a久久精品蜜桃钙片| 亚洲av综合色区一区| 久久久亚洲精品成人影院| 免费观看的影片在线观看| 国产成人精品福利久久| 日本av手机在线免费观看| 久久久精品免费免费高清| 一本色道久久久久久精品综合| 国产欧美日韩一区二区三区在线 | 久久久久久久久久成人| 97超视频在线观看视频| 又爽又黄a免费视频| 欧美日韩精品成人综合77777| 91久久精品国产一区二区成人| 日韩在线高清观看一区二区三区| 激情五月婷婷亚洲| 80岁老熟妇乱子伦牲交| 精品午夜福利在线看| 国产深夜福利视频在线观看| 精华霜和精华液先用哪个| 好男人视频免费观看在线| 日韩视频在线欧美| 蜜臀久久99精品久久宅男| 亚洲国产最新在线播放| 亚洲精品一二三| 国产成人免费观看mmmm| 97精品久久久久久久久久精品| 国产精品一区二区在线观看99| 精品久久久久久久久av| 晚上一个人看的免费电影| 精品酒店卫生间| 一级片'在线观看视频| 欧美精品亚洲一区二区| 欧美+日韩+精品| 日韩伦理黄色片| 最后的刺客免费高清国语| 国产大屁股一区二区在线视频| 久久久久久久久久人人人人人人| 久久99热这里只有精品18| 国产综合精华液| 韩国高清视频一区二区三区| 国产精品.久久久| 亚州av有码| 建设人人有责人人尽责人人享有的 | 国产成人精品一,二区| 深夜a级毛片| 免费观看的影片在线观看| 精品久久久久久久末码| 精品一品国产午夜福利视频| 成人国产麻豆网| 色婷婷久久久亚洲欧美| 久久久久精品性色| 亚洲久久久国产精品| 黑人高潮一二区| 亚洲欧美一区二区三区国产| 亚洲精华国产精华液的使用体验| 草草在线视频免费看| av国产久精品久网站免费入址| 成人二区视频| 麻豆成人av视频| 精品人妻偷拍中文字幕| 午夜免费鲁丝| 中文字幕亚洲精品专区| 国产在线免费精品| xxx大片免费视频| 国产精品一二三区在线看| 免费大片黄手机在线观看| 99久久精品热视频| 亚洲欧美中文字幕日韩二区| 久久精品国产亚洲av天美| 欧美日韩一区二区视频在线观看视频在线| 五月玫瑰六月丁香| videos熟女内射| 草草在线视频免费看| 99热这里只有精品一区| 久久99热这里只有精品18| 欧美成人午夜免费资源| 国产成人91sexporn| 国产在线一区二区三区精| 少妇人妻 视频| 国产精品蜜桃在线观看| 日韩精品有码人妻一区| 国产91av在线免费观看| 亚洲国产精品专区欧美| 内地一区二区视频在线| 亚洲av.av天堂| 国产日韩欧美亚洲二区| 99久久精品热视频| 老女人水多毛片| 高清日韩中文字幕在线| 中文在线观看免费www的网站| 亚洲国产精品专区欧美| 又黄又爽又刺激的免费视频.| 亚洲一级一片aⅴ在线观看| 久久久久久久久久久丰满| 亚洲伊人久久精品综合| 国产一区有黄有色的免费视频| 亚洲精品乱久久久久久| 蜜桃亚洲精品一区二区三区| 日韩欧美 国产精品| 高清不卡的av网站| 欧美一级a爱片免费观看看| 日韩精品有码人妻一区| 亚洲精品第二区| 人体艺术视频欧美日本| 亚洲av电影在线观看一区二区三区| av在线老鸭窝| 日韩制服骚丝袜av| 五月天丁香电影| 人人妻人人澡人人爽人人夜夜| 美女国产视频在线观看| 免费黄色在线免费观看| 国产高清不卡午夜福利| 亚洲av.av天堂| 国产精品久久久久久久久免| 91午夜精品亚洲一区二区三区| 最近手机中文字幕大全| 久久久欧美国产精品| 婷婷色综合大香蕉| 在现免费观看毛片| av福利片在线观看| 一本—道久久a久久精品蜜桃钙片| 国产精品蜜桃在线观看| 日产精品乱码卡一卡2卡三| 麻豆成人午夜福利视频| 插逼视频在线观看| 国产精品久久久久久av不卡| 内射极品少妇av片p| 亚洲国产欧美在线一区| 国产在线一区二区三区精| 精品久久久久久久久亚洲| 丝袜脚勾引网站| 久久久久久久久久久免费av| 久久99热这里只频精品6学生| 国产在线男女| 日韩三级伦理在线观看| 国产在线一区二区三区精| 哪个播放器可以免费观看大片| 国产视频内射| 男女边摸边吃奶| 18禁裸乳无遮挡免费网站照片| 亚洲久久久国产精品| 青春草亚洲视频在线观看| 久久这里有精品视频免费| 最近最新中文字幕免费大全7| 国产精品久久久久久久久免| 免费人妻精品一区二区三区视频| 99久久精品国产国产毛片| 免费看光身美女| 欧美精品国产亚洲| 国产乱人偷精品视频| 欧美bdsm另类| 亚洲精品456在线播放app| 久久99精品国语久久久| 亚洲欧美中文字幕日韩二区| 你懂的网址亚洲精品在线观看| 久久99热这里只频精品6学生| 午夜视频国产福利| 又粗又硬又长又爽又黄的视频| 如何舔出高潮| 深夜a级毛片| 一级黄片播放器| a级一级毛片免费在线观看| 欧美xxⅹ黑人| 亚洲国产最新在线播放| 黄色配什么色好看| 久久精品久久精品一区二区三区| 高清视频免费观看一区二区| 成人影院久久| 久久毛片免费看一区二区三区| 免费av中文字幕在线| 成人一区二区视频在线观看| www.色视频.com| 亚洲精品乱码久久久v下载方式| 国产黄频视频在线观看| 最近的中文字幕免费完整| 中文天堂在线官网| 国产老妇伦熟女老妇高清| 国产高清不卡午夜福利| 国产国拍精品亚洲av在线观看| 99热这里只有精品一区| 男女边摸边吃奶| 在线 av 中文字幕| 网址你懂的国产日韩在线| 久久青草综合色| 免费看光身美女| 中国国产av一级| 国产淫语在线视频| 在线观看三级黄色| 亚洲国产精品国产精品| 日韩伦理黄色片| 亚洲精品中文字幕在线视频 | 成人午夜精彩视频在线观看| 日韩av免费高清视频| 夜夜爽夜夜爽视频| 亚洲国产高清在线一区二区三| 精品国产三级普通话版| 日韩一区二区三区影片| 欧美日韩在线观看h| 国产精品久久久久久精品古装| 男女免费视频国产| 热99国产精品久久久久久7| 黄色日韩在线| 国语对白做爰xxxⅹ性视频网站| 99re6热这里在线精品视频| 亚洲国产精品一区三区| 久久久久久人妻| 一区二区av电影网| 国产深夜福利视频在线观看| 亚洲欧洲国产日韩| 新久久久久国产一级毛片| 爱豆传媒免费全集在线观看| 在线观看免费视频网站a站| 成人影院久久| 免费少妇av软件| 久久久成人免费电影| 99热这里只有精品一区| av免费观看日本| 国产精品久久久久久久久免| 成人黄色视频免费在线看| 妹子高潮喷水视频| 日韩av在线免费看完整版不卡| 熟女人妻精品中文字幕| 天天躁夜夜躁狠狠久久av| 乱系列少妇在线播放| 久久热精品热| 乱系列少妇在线播放| 国产爱豆传媒在线观看| 久久久久网色| 国产在线免费精品| 亚洲精品aⅴ在线观看| 日韩亚洲欧美综合| 在线天堂最新版资源| 午夜福利影视在线免费观看| 黄色日韩在线| 18禁在线无遮挡免费观看视频| 老司机影院毛片| 一级毛片我不卡| 美女脱内裤让男人舔精品视频| 最后的刺客免费高清国语| 精品国产露脸久久av麻豆| 一区二区三区免费毛片| 久久国产精品大桥未久av | 欧美日韩亚洲高清精品| 18禁裸乳无遮挡动漫免费视频| 在线播放无遮挡| 欧美日韩国产mv在线观看视频 | 国产爽快片一区二区三区| 国产一区亚洲一区在线观看| 免费看不卡的av| 久久午夜福利片| 亚洲精品亚洲一区二区| 亚洲人成网站在线播| 有码 亚洲区| 五月开心婷婷网| 色婷婷久久久亚洲欧美| 精品久久久久久久久亚洲| 男人添女人高潮全过程视频| 女人十人毛片免费观看3o分钟| 大话2 男鬼变身卡| 欧美日韩亚洲高清精品| 国产精品免费大片| 久久久久久久久久久免费av| 午夜激情福利司机影院| 日日撸夜夜添| 人人妻人人爽人人添夜夜欢视频 | 青春草国产在线视频| 久久久久久久久大av| 久久精品熟女亚洲av麻豆精品| 色视频在线一区二区三区| 99视频精品全部免费 在线| 国产精品一二三区在线看| 成人国产av品久久久| 超碰97精品在线观看| 国产一区亚洲一区在线观看| 国产一区二区三区综合在线观看 | 国产真实伦视频高清在线观看| 亚洲av成人精品一区久久| 国产黄片视频在线免费观看| 亚洲av欧美aⅴ国产| 中文字幕亚洲精品专区| 欧美少妇被猛烈插入视频| 三级经典国产精品| 亚洲国产色片| 亚洲成人中文字幕在线播放| 秋霞伦理黄片| 大片免费播放器 马上看| 午夜福利在线在线| 又黄又爽又刺激的免费视频.| 校园人妻丝袜中文字幕| 身体一侧抽搐| 国产一级毛片在线| 精品久久久久久久久av| av天堂中文字幕网| 欧美丝袜亚洲另类| 亚洲国产欧美在线一区| 女性生殖器流出的白浆| 少妇 在线观看| av免费在线看不卡| 尾随美女入室| 国产片特级美女逼逼视频| 成人国产麻豆网| 亚洲欧美成人精品一区二区| 欧美成人一区二区免费高清观看| 男女边摸边吃奶| 老女人水多毛片| 最近中文字幕高清免费大全6| 国内少妇人妻偷人精品xxx网站| 伊人久久国产一区二区| 色视频在线一区二区三区| 边亲边吃奶的免费视频| 国产在线一区二区三区精| 深爱激情五月婷婷| 夜夜看夜夜爽夜夜摸| 亚洲国产成人一精品久久久| 91久久精品国产一区二区三区| 欧美日韩一区二区视频在线观看视频在线| 嘟嘟电影网在线观看| 久久久国产一区二区| .国产精品久久| 超碰97精品在线观看| 麻豆国产97在线/欧美| 99精国产麻豆久久婷婷| 国产精品久久久久久久久免| 欧美少妇被猛烈插入视频| 全区人妻精品视频| 欧美极品一区二区三区四区| 97在线视频观看| 九九在线视频观看精品| 亚洲av综合色区一区| 亚洲欧洲日产国产| 精品酒店卫生间| 99热这里只有是精品在线观看| 大片电影免费在线观看免费| 久久久久久久久久人人人人人人| 蜜桃亚洲精品一区二区三区|