• <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
    日日夜夜操网爽| 操出白浆在线播放| 青草久久国产| 亚洲精品久久国产高清桃花| av超薄肉色丝袜交足视频| 18禁裸乳无遮挡免费网站照片 | 亚洲成人久久爱视频| 狂野欧美激情性xxxx| 国产精品久久电影中文字幕| 最新美女视频免费是黄的| 亚洲 欧美一区二区三区| 女人爽到高潮嗷嗷叫在线视频| 777久久人妻少妇嫩草av网站| 午夜视频精品福利| 美女午夜性视频免费| 黄色视频不卡| 99久久无色码亚洲精品果冻| 制服丝袜大香蕉在线| 露出奶头的视频| 搡老妇女老女人老熟妇| а√天堂www在线а√下载| 国产欧美日韩一区二区精品| 黄色丝袜av网址大全| 丁香六月欧美| 免费在线观看视频国产中文字幕亚洲| 免费无遮挡裸体视频| 91九色精品人成在线观看| 欧美亚洲日本最大视频资源| 波多野结衣av一区二区av| 国产欧美日韩一区二区三| 黄色视频,在线免费观看| 桃红色精品国产亚洲av| 欧美丝袜亚洲另类 | 日韩欧美在线二视频| 亚洲中文日韩欧美视频| 精品高清国产在线一区| 一边摸一边做爽爽视频免费| 成人一区二区视频在线观看| 国产精品乱码一区二三区的特点| 中文亚洲av片在线观看爽| 91字幕亚洲| 欧美黄色片欧美黄色片| 亚洲aⅴ乱码一区二区在线播放 | 欧美乱色亚洲激情| 国产99白浆流出| 久久人人精品亚洲av| 亚洲精品国产精品久久久不卡| 在线观看午夜福利视频| 久久久久久亚洲精品国产蜜桃av| 脱女人内裤的视频| 国产视频一区二区在线看| 一区福利在线观看| 国产一区二区三区在线臀色熟女| 好男人在线观看高清免费视频 | 久久99热这里只有精品18| 叶爱在线成人免费视频播放| 国产欧美日韩一区二区精品| 中文字幕最新亚洲高清| 欧美国产精品va在线观看不卡| 亚洲国产欧美日韩在线播放| 国产又爽黄色视频| 亚洲人成77777在线视频| 国产色视频综合| 特大巨黑吊av在线直播 | 精品久久久久久久毛片微露脸| 国产亚洲av嫩草精品影院| 91麻豆av在线| 国产精品亚洲美女久久久| 男人舔女人的私密视频| 亚洲av成人av| 久久精品成人免费网站| 中出人妻视频一区二区| 热re99久久国产66热| 亚洲熟女毛片儿| 男人舔女人下体高潮全视频| 久久香蕉国产精品| 黄频高清免费视频| 1024视频免费在线观看| 国产成人欧美在线观看| 看黄色毛片网站| 国产成人精品久久二区二区免费| 久久久久免费精品人妻一区二区 | 欧美日本亚洲视频在线播放| 99热只有精品国产| 十八禁人妻一区二区| 老司机福利观看| 精品卡一卡二卡四卡免费| 亚洲人成网站在线播放欧美日韩| 村上凉子中文字幕在线| 亚洲欧美日韩无卡精品| 国产又爽黄色视频| 国产91精品成人一区二区三区| 日韩三级视频一区二区三区| 国产精品 国内视频| 日韩有码中文字幕| 两个人看的免费小视频| 老司机在亚洲福利影院| 国产av在哪里看| 久久精品影院6| 亚洲精品国产区一区二| 男女下面进入的视频免费午夜 | 人人妻人人澡人人看| 一本精品99久久精品77| 欧美 亚洲 国产 日韩一| 亚洲国产欧美日韩在线播放| 此物有八面人人有两片| 亚洲自拍偷在线| 男人舔女人下体高潮全视频| 久久中文看片网| 天天添夜夜摸| 欧美性猛交╳xxx乱大交人| 一个人免费在线观看的高清视频| 国产精品精品国产色婷婷| 视频区欧美日本亚洲| videosex国产| 国产黄色小视频在线观看| 亚洲无线在线观看| 亚洲人成77777在线视频| 最近最新中文字幕大全免费视频| 国产成人精品久久二区二区91| 香蕉丝袜av| 一个人免费在线观看的高清视频| 亚洲九九香蕉| 亚洲av电影在线进入| 欧美日韩中文字幕国产精品一区二区三区| 美女高潮喷水抽搐中文字幕| 亚洲中文av在线| 久久人妻福利社区极品人妻图片| 欧美性猛交黑人性爽| 亚洲熟妇熟女久久| 国产一区二区在线av高清观看| 他把我摸到了高潮在线观看| 黄色女人牲交| 国产精品免费一区二区三区在线| 99热只有精品国产| 露出奶头的视频| 少妇裸体淫交视频免费看高清 | 亚洲最大成人中文| 中文字幕精品免费在线观看视频| 欧美乱色亚洲激情| 麻豆成人av在线观看| 在线观看66精品国产| 日本熟妇午夜| 亚洲五月色婷婷综合| 亚洲国产欧美日韩在线播放| 一区二区三区高清视频在线| www.自偷自拍.com| 人人妻人人澡欧美一区二区| 国产又色又爽无遮挡免费看| АⅤ资源中文在线天堂| 在线天堂中文资源库| 视频区欧美日本亚洲| 国产精品精品国产色婷婷| av有码第一页| 国产亚洲精品久久久久5区| 国产欧美日韩精品亚洲av| 99国产综合亚洲精品| 1024视频免费在线观看| bbb黄色大片| 美女扒开内裤让男人捅视频| 一级毛片女人18水好多| 亚洲 欧美一区二区三区| 搡老熟女国产l中国老女人| 国产精品国产高清国产av| 久久国产精品影院| 亚洲精品国产区一区二| 国产爱豆传媒在线观看 | 亚洲av片天天在线观看| 真人做人爱边吃奶动态| netflix在线观看网站| 成人国产综合亚洲| 亚洲一区中文字幕在线| 黄网站色视频无遮挡免费观看| 欧美亚洲日本最大视频资源| av中文乱码字幕在线| 久久久国产欧美日韩av| 亚洲中文字幕一区二区三区有码在线看 | 久久 成人 亚洲| 国产精品国产高清国产av| 国产精品98久久久久久宅男小说| 女同久久另类99精品国产91| 亚洲国产欧美日韩在线播放| 18禁黄网站禁片免费观看直播| av在线播放免费不卡| 久久中文看片网| a级毛片a级免费在线| 国产精品日韩av在线免费观看| 亚洲av日韩精品久久久久久密| 免费女性裸体啪啪无遮挡网站| 91成年电影在线观看| 日韩大码丰满熟妇| 免费在线观看成人毛片| 无人区码免费观看不卡| 精品午夜福利视频在线观看一区| 夜夜躁狠狠躁天天躁| 色播在线永久视频| 免费一级毛片在线播放高清视频| 久久久久久久精品吃奶| 欧美zozozo另类| 十分钟在线观看高清视频www| 免费看日本二区| 成人18禁在线播放| 真人做人爱边吃奶动态| 一二三四在线观看免费中文在| 国产午夜福利久久久久久| 国产真人三级小视频在线观看| 琪琪午夜伦伦电影理论片6080| videosex国产| 亚洲五月天丁香| 性欧美人与动物交配| 男男h啪啪无遮挡| 亚洲国产高清在线一区二区三 | 亚洲美女黄片视频| 午夜福利欧美成人| 两个人免费观看高清视频| 黄色视频不卡| 香蕉丝袜av| 两性夫妻黄色片| 桃色一区二区三区在线观看| 一卡2卡三卡四卡精品乱码亚洲| 人人澡人人妻人| 在线十欧美十亚洲十日本专区| 国产午夜福利久久久久久| 亚洲国产看品久久| 色综合站精品国产| 老鸭窝网址在线观看| 日本一区二区免费在线视频| 成熟少妇高潮喷水视频| 免费一级毛片在线播放高清视频| 手机成人av网站| 日韩免费av在线播放| 国产一卡二卡三卡精品| 免费女性裸体啪啪无遮挡网站| 人人妻人人澡人人看| √禁漫天堂资源中文www| 日日爽夜夜爽网站| 少妇熟女aⅴ在线视频| 日本成人三级电影网站| 亚洲成人久久性| 国产精品久久久av美女十八| 在线永久观看黄色视频| 女生性感内裤真人,穿戴方法视频| 一夜夜www| 在线观看免费午夜福利视频| www.熟女人妻精品国产| 啦啦啦韩国在线观看视频| 欧美 亚洲 国产 日韩一| 亚洲三区欧美一区| 国产成人欧美| 色综合亚洲欧美另类图片| 嫁个100分男人电影在线观看| 久久久久久久久中文| 国产av在哪里看| 国产精品av久久久久免费| 国产成人一区二区三区免费视频网站| 精品国产国语对白av| 久久久久亚洲av毛片大全| 嫩草影院精品99| 法律面前人人平等表现在哪些方面| 一本大道久久a久久精品| 亚洲成人精品中文字幕电影| 人成视频在线观看免费观看| 国产av一区在线观看免费| 亚洲国产高清在线一区二区三 | 女性生殖器流出的白浆| 久久久久久亚洲精品国产蜜桃av| 成人国产综合亚洲| 亚洲激情在线av| 国产成人影院久久av| 免费观看精品视频网站| 亚洲电影在线观看av| av有码第一页| 亚洲人成电影免费在线| 男男h啪啪无遮挡| 久久人人精品亚洲av| 国产午夜精品久久久久久| 亚洲国产欧美一区二区综合| 最新在线观看一区二区三区| 韩国av一区二区三区四区| 757午夜福利合集在线观看| 午夜免费成人在线视频| 久久中文看片网| 国产熟女午夜一区二区三区| 久久久久免费精品人妻一区二区 | 美女高潮喷水抽搐中文字幕| 麻豆国产av国片精品| 精品久久久久久久久久久久久 | 日本黄色视频三级网站网址| 岛国视频午夜一区免费看| 欧美日韩乱码在线| 在线观看日韩欧美| 99精品在免费线老司机午夜| 一级黄色大片毛片| 女人高潮潮喷娇喘18禁视频| 黄片小视频在线播放| 亚洲黑人精品在线| 成年女人毛片免费观看观看9| 观看免费一级毛片| 欧美日韩精品网址| 99国产极品粉嫩在线观看| 啦啦啦 在线观看视频| 在线观看免费日韩欧美大片| 午夜久久久在线观看| 在线十欧美十亚洲十日本专区| 精华霜和精华液先用哪个| a级毛片在线看网站| 国产在线精品亚洲第一网站| 一区二区三区激情视频| 国产精品,欧美在线| 成年人黄色毛片网站| 精品一区二区三区av网在线观看| 女同久久另类99精品国产91| 一进一出好大好爽视频| 男人的好看免费观看在线视频 | 两性午夜刺激爽爽歪歪视频在线观看 | 日韩高清综合在线| 午夜久久久在线观看| 中文亚洲av片在线观看爽| 欧美黄色片欧美黄色片| 欧美日韩福利视频一区二区| 免费在线观看日本一区| 国产一区二区三区视频了| 波多野结衣高清作品| 91国产中文字幕| √禁漫天堂资源中文www| 一级a爱视频在线免费观看| 深夜精品福利| 一本一本综合久久| 国产伦在线观看视频一区| 国产精品,欧美在线| 夜夜夜夜夜久久久久| 国产亚洲av高清不卡| 亚洲国产欧洲综合997久久, | 侵犯人妻中文字幕一二三四区| 亚洲av日韩精品久久久久久密| 日本在线视频免费播放| 国产av在哪里看| 男人的好看免费观看在线视频 | 久热爱精品视频在线9| 天堂动漫精品| 日韩一卡2卡3卡4卡2021年| 香蕉久久夜色| www日本黄色视频网| 国产高清有码在线观看视频 | 天天添夜夜摸| 无人区码免费观看不卡| 亚洲av成人一区二区三| 一级a爱片免费观看的视频| 亚洲五月天丁香| 级片在线观看| 悠悠久久av| 国产三级黄色录像| 免费在线观看成人毛片| 一本综合久久免费| 男女午夜视频在线观看| 日本免费a在线| 丰满人妻熟妇乱又伦精品不卡| 精品国产超薄肉色丝袜足j| 91成年电影在线观看| 在线观看免费午夜福利视频| 中文字幕av电影在线播放| 午夜影院日韩av| 性色av乱码一区二区三区2| 国产单亲对白刺激| 一级黄色大片毛片| 亚洲中文av在线| 国产麻豆成人av免费视频| 欧美日本视频| 男人舔女人下体高潮全视频| 亚洲精品粉嫩美女一区| 亚洲欧美激情综合另类| 少妇熟女aⅴ在线视频| 黄片小视频在线播放| 黄色片一级片一级黄色片| 中文字幕高清在线视频| 国产又爽黄色视频| a级毛片在线看网站| xxx96com| 悠悠久久av| 国产三级黄色录像| 91麻豆av在线| 亚洲 欧美 日韩 在线 免费| 亚洲av日韩精品久久久久久密| 女人被狂操c到高潮| 欧美成人午夜精品| 久久久久久国产a免费观看| 久久精品国产99精品国产亚洲性色| 色婷婷久久久亚洲欧美| 欧美激情 高清一区二区三区| 91国产中文字幕| 搡老熟女国产l中国老女人| 久久婷婷人人爽人人干人人爱| 国产成人影院久久av| 黄色丝袜av网址大全| 久久精品91蜜桃| 一个人免费在线观看的高清视频| 亚洲一区高清亚洲精品| 777久久人妻少妇嫩草av网站| 波多野结衣高清作品| 女人被狂操c到高潮| 欧美国产精品va在线观看不卡| 国产主播在线观看一区二区| www.www免费av| 午夜免费鲁丝| 日韩欧美国产一区二区入口| 一边摸一边抽搐一进一小说| 可以在线观看毛片的网站| 亚洲精品在线观看二区| 亚洲五月天丁香| 日日夜夜操网爽| 久久久久亚洲av毛片大全| 亚洲av熟女| 高潮久久久久久久久久久不卡| 成人国产一区最新在线观看| 亚洲五月天丁香| 可以在线观看毛片的网站| 国产在线观看jvid| 丝袜人妻中文字幕| 免费一级毛片在线播放高清视频| 欧美又色又爽又黄视频| а√天堂www在线а√下载| 观看免费一级毛片| 久久国产乱子伦精品免费另类| 色播亚洲综合网| 日日摸夜夜添夜夜添小说| 无遮挡黄片免费观看| 午夜福利视频1000在线观看| 热re99久久国产66热| 久久狼人影院| av福利片在线| 国产免费男女视频| 国产主播在线观看一区二区| av福利片在线| 黄色成人免费大全| 亚洲av日韩精品久久久久久密| 国产一级毛片七仙女欲春2 | 亚洲av美国av| 国产av又大| 91九色精品人成在线观看| 国产不卡一卡二| 成年免费大片在线观看| netflix在线观看网站| 久久精品aⅴ一区二区三区四区| 国产成人av教育| 女同久久另类99精品国产91| 欧美黑人精品巨大| 18美女黄网站色大片免费观看| 黑人欧美特级aaaaaa片| 国内少妇人妻偷人精品xxx网站 | 欧美一级a爱片免费观看看 | 2021天堂中文幕一二区在线观 | 欧美乱码精品一区二区三区| 久久婷婷人人爽人人干人人爱| 美女 人体艺术 gogo| 曰老女人黄片| 黄色视频不卡| 国产av一区在线观看免费| 最好的美女福利视频网| 99riav亚洲国产免费| 亚洲av片天天在线观看| 美女免费视频网站| 午夜福利视频1000在线观看| 亚洲av成人一区二区三| 午夜免费观看网址| 桃色一区二区三区在线观看| 国产99白浆流出| 久久精品国产亚洲av高清一级| 精品一区二区三区av网在线观看| 美女免费视频网站| 最近在线观看免费完整版| 亚洲色图av天堂| 嫩草影院精品99| 黄色a级毛片大全视频| 99re在线观看精品视频| 久久精品国产99精品国产亚洲性色| 成人av一区二区三区在线看| 这个男人来自地球电影免费观看| 国产伦一二天堂av在线观看| 免费看a级黄色片| 日本免费a在线| 91麻豆av在线| 亚洲国产欧美日韩在线播放| 日日干狠狠操夜夜爽| 精品久久久久久久人妻蜜臀av| 草草在线视频免费看| av片东京热男人的天堂| 老汉色∧v一级毛片| 亚洲国产欧美一区二区综合| 美国免费a级毛片| 国产一区二区激情短视频| 亚洲人成网站高清观看| 日韩 欧美 亚洲 中文字幕| 国内精品久久久久久久电影| 亚洲国产高清在线一区二区三 | 视频在线观看一区二区三区| 亚洲av成人一区二区三| 99热只有精品国产| 午夜福利成人在线免费观看| 天天躁夜夜躁狠狠躁躁| 久久天躁狠狠躁夜夜2o2o| 中文字幕精品免费在线观看视频| 变态另类丝袜制服| xxx96com| 久久国产乱子伦精品免费另类| 国产成人一区二区三区免费视频网站| 精品一区二区三区视频在线观看免费| 欧美在线一区亚洲| 熟妇人妻久久中文字幕3abv| 女生性感内裤真人,穿戴方法视频| 制服人妻中文乱码| 久热这里只有精品99| 亚洲成国产人片在线观看| ponron亚洲| 真人一进一出gif抽搐免费| 黄色 视频免费看| 99riav亚洲国产免费| 又黄又粗又硬又大视频| 国产欧美日韩精品亚洲av| 色尼玛亚洲综合影院| 日韩成人在线观看一区二区三区| 99久久精品国产亚洲精品| 国产精品综合久久久久久久免费| 1024视频免费在线观看| av电影中文网址| 国产精品爽爽va在线观看网站 | 成年版毛片免费区| 黄频高清免费视频| 国产亚洲精品综合一区在线观看 | 日本 av在线| 成年免费大片在线观看| 18禁观看日本| 热99re8久久精品国产| 日本免费一区二区三区高清不卡| 中文字幕人妻熟女乱码| 麻豆久久精品国产亚洲av| 18美女黄网站色大片免费观看| 亚洲av五月六月丁香网| 亚洲av熟女| 美女国产高潮福利片在线看| 亚洲七黄色美女视频| 在线十欧美十亚洲十日本专区| 免费观看人在逋| 精品高清国产在线一区| 99在线人妻在线中文字幕| 黄色丝袜av网址大全| 黄片小视频在线播放| 国产av一区在线观看免费| 日本 av在线| 特大巨黑吊av在线直播 | 久久午夜亚洲精品久久| 亚洲人成网站在线播放欧美日韩| 午夜福利在线观看吧| 亚洲成人免费电影在线观看| 国产国语露脸激情在线看| 欧美av亚洲av综合av国产av| √禁漫天堂资源中文www| 男女下面进入的视频免费午夜 | 一级a爱视频在线免费观看| 午夜精品在线福利| 国产精品影院久久| 中出人妻视频一区二区| 在线永久观看黄色视频| 窝窝影院91人妻| 欧美激情极品国产一区二区三区| 精品卡一卡二卡四卡免费| av在线天堂中文字幕| 久久精品国产综合久久久| 少妇裸体淫交视频免费看高清 | 美女 人体艺术 gogo| 男人舔奶头视频| 非洲黑人性xxxx精品又粗又长| 伊人久久大香线蕉亚洲五| bbb黄色大片| 国产亚洲精品第一综合不卡| 法律面前人人平等表现在哪些方面| 99热只有精品国产| 精品久久久久久久人妻蜜臀av| 一个人免费在线观看的高清视频| 亚洲国产日韩欧美精品在线观看 | 三级毛片av免费| 日韩精品青青久久久久久| 免费看日本二区| 一区二区日韩欧美中文字幕| 欧美丝袜亚洲另类 | 少妇粗大呻吟视频| 人人妻人人澡人人看| 日韩成人在线观看一区二区三区| 欧美一级毛片孕妇| 在线永久观看黄色视频| 精品久久久久久久久久免费视频| 啦啦啦 在线观看视频| 久久精品国产综合久久久| 欧美日韩乱码在线| 亚洲国产欧美网| а√天堂www在线а√下载| 久久国产精品影院| 欧美成狂野欧美在线观看| 欧美又色又爽又黄视频| 久久久久国产一级毛片高清牌| 在线观看一区二区三区| 制服人妻中文乱码| 亚洲av成人av| 亚洲国产欧美网| 免费女性裸体啪啪无遮挡网站| 黑人巨大精品欧美一区二区mp4| 欧美色视频一区免费| 亚洲专区中文字幕在线| 久久久国产成人精品二区| 欧美午夜高清在线| 免费一级毛片在线播放高清视频| 少妇裸体淫交视频免费看高清 | 性欧美人与动物交配| 日韩精品免费视频一区二区三区| 日韩欧美国产在线观看| 久久精品影院6|