孫林潔,劉春成,張世鑫

(中海油研究總院,北京100028)

PML邊界條件下孔隙介質(zhì)彈性波可變網(wǎng)格正演模擬方法研究

孫林潔,劉春成,張世鑫

(中海油研究總院,北京100028)

孔隙介質(zhì)彈性波正演模擬對于儲層預(yù)測和油氣檢測具有重要的指導(dǎo)意義。為提高孔隙介質(zhì)波場正演模擬的精度,將可變網(wǎng)格方法應(yīng)用于基于Biot理論的雙相介質(zhì)彈性波高階有限差分正演模擬方法中,有效解決了二階混合偏導(dǎo)數(shù)差分近似的可變網(wǎng)格過渡問題,推導(dǎo)了基于PML邊界條件的可變網(wǎng)格差分格式,實(shí)現(xiàn)了空間與時間步長任意整數(shù)倍變化的雙相介質(zhì)彈性波可變網(wǎng)格正演模擬。對于細(xì)微地質(zhì)結(jié)構(gòu)或橫向變化劇烈的介質(zhì)區(qū)域,利用可變網(wǎng)格對區(qū)域精細(xì)重采樣可有效提高正演的模擬精度;對于介質(zhì)相對穩(wěn)定的區(qū)域,采用大步長采樣可保證正演計算效率。數(shù)值試驗(yàn)結(jié)果表明,與常規(guī)正演模擬方法相比,雙相介質(zhì)可變網(wǎng)格正演模擬方法的精度、效率有明顯提高,對精細(xì)研究孔隙介質(zhì)地震響應(yīng)有著重要意義。

Biot理論;孔隙介質(zhì);PML邊界條件;可變網(wǎng)格;有限差分

近年來,隨著石油勘探開發(fā)技術(shù)的發(fā)展,對孔隙介質(zhì)的關(guān)注程度日益提高。當(dāng)孔隙中含有流體時,不同流體的存在會對介質(zhì)參數(shù)產(chǎn)生不同的影響,從而影響到地震反射波的走時、振幅和頻率等特征[1]。因此,根據(jù)油氣儲層在不同地層條件和流體相態(tài)下的地震波反射特征進(jìn)行正演模擬,分析不同孔隙流體的地球物理特性在各種地震信息上的表現(xiàn),總結(jié)規(guī)律和認(rèn)識,對指導(dǎo)儲層預(yù)測和油氣檢測等工作具有重要意義。

波動方程正演模擬方法是探索地震波在介質(zhì)中傳播規(guī)律的重要手段。常規(guī)波動方程正演模擬方法基于等效理論,不單獨(dú)考慮介質(zhì)孔隙中的流體,將介質(zhì)中固、液兩相對地震波的作用簡化為等效介質(zhì)參數(shù),所得到的模擬結(jié)果無法體現(xiàn)流體響應(yīng)特征。Biot提出的雙相介質(zhì)理論[2-3]已經(jīng)得到廣泛應(yīng)用,Biot模型假設(shè)巖石由骨架和孔隙中的流體構(gòu)成,這種模型假設(shè)考慮了孔隙流體的作用[4-5],比等效單相介質(zhì)模型更為符合實(shí)際情況[6],適合討論、分析含流體儲層的地震響應(yīng)特征。

有限差分方法是一種常用的波動方程數(shù)值解法,具有精度高、計算速度快的優(yōu)點(diǎn),但是該方法也具有一定局限性,如存在數(shù)值頻散,且必須滿足方法穩(wěn)定性條件等[7]。常規(guī)差分方法使用規(guī)則網(wǎng)格對模型進(jìn)行剖分,一般在同一模型內(nèi)部使用統(tǒng)一的采樣步長。而孔隙介質(zhì)多為非均勻介質(zhì),模型參數(shù)在不同區(qū)域變化并不均勻,使用單一尺度的網(wǎng)格步長對模型進(jìn)行剖分,必然會導(dǎo)致局部過采樣影響計算效率或局部采樣不足影響模擬精度。Moczo[8]首次提出可變網(wǎng)格的思路,數(shù)值模擬時根據(jù)模型參數(shù)變化調(diào)整計算網(wǎng)格步長,針對尺度較小的地質(zhì)體和模型參數(shù)橫向變化劇烈的區(qū)域采用精細(xì)網(wǎng)格剖分;針對模型參數(shù)變化平緩的區(qū)域采用大網(wǎng)格剖分。這種方法既可以保證正演模擬的精度,同時在一定程度上提高了計算效率。

可變網(wǎng)格方法提出至今已被成功應(yīng)用于常規(guī)單相介質(zhì)聲波[9-11]、彈性波[12-14]的有限差分正演模擬當(dāng)中。本文將可變網(wǎng)格方法推廣至基于Biot理論的雙相介質(zhì)彈性波二階位移方程的正演模擬中,成功地解決了二階混合偏導(dǎo)數(shù)常規(guī)差分近似的可變網(wǎng)格過渡問題。文中正演模擬使用孔隙介質(zhì)的完全匹配層(PML)邊界條件,推導(dǎo)了可變網(wǎng)格區(qū)域靠近模型邊界時的差分格式,從而得到任意倍數(shù)、任意可變網(wǎng)格位置的孔隙介質(zhì)可變網(wǎng)格正演模擬算法,實(shí)現(xiàn)了對孔隙介質(zhì)的高效、精確數(shù)值模擬。

1 方法原理

1.1 雙相介質(zhì)彈性波方程PML差分格式

Biot給出的彈性波方程形式為[2]:

(1)

式中:u和U分別為固相與液相位移;e和ξ分別為固相與液相介質(zhì)的體應(yīng)變;A,N,Q,R為雙相介質(zhì)彈性參數(shù),表征固相、液相介質(zhì)的彈性性質(zhì)與相互的耦合作用;ρ11,ρ22,ρ12為雙相介質(zhì)質(zhì)量參數(shù),表征固相、液相介質(zhì)單位體積的等效質(zhì)量與相互作用對質(zhì)量的影響;B為耗散參數(shù),表征液相粘滯性對地震波的衰減效果。將方程(1)展開成為二階位移形式,在二維情況下,以x方向?yàn)槔?有:

(2)

一般而言,求解彈性波方程PML差分格式需要將原方程加入衰減項(xiàng)后根據(jù)空間導(dǎo)數(shù)方向拆分,并通過坐標(biāo)變換的方式得到與原方程同解的穩(wěn)定邊界格式[15]。與彈性波方程中僅包含二階時間導(dǎo)數(shù)的情況不同,考慮介質(zhì)耗散的雙相彈性波方程中同時存在一階與二階時間導(dǎo)數(shù),分解方程會形成欠定方程組。為簡化這個問題,從PML邊界條件吸收衰減方式入手,在保持原波動方程形式的同時加入衰減項(xiàng)。加入PML衰減項(xiàng)的波動方程形式需要進(jìn)行以下變換[15-17]:

(3)

式中:βi為i方向的衰減項(xiàng)。則加入衰減項(xiàng)的PML方程為:

(4)

該方程解的離散差分形式求解方式詳見附錄A。

至此我們得到了雙相介質(zhì)彈性波方程的PML差分格式。在正演模擬過程中,通過多次數(shù)值試驗(yàn),合理選取模型邊界和角點(diǎn)位置的PML吸收參數(shù),而令模型區(qū)域的吸收參數(shù)為0[15],即可消除模型邊界反射,得到精確的模擬結(jié)果。

1.2 雙相介質(zhì)可變網(wǎng)格PML差分格式

可變網(wǎng)格有限差分方法的主要特點(diǎn)為:在需要精細(xì)模擬的區(qū)域使用小網(wǎng)格對模型進(jìn)行重采樣計算,其它區(qū)域使用大網(wǎng)格計算[11]?？勺兙W(wǎng)格模型剖分的示意圖如圖1所示。

常規(guī)差分格式具有對稱性,針對求導(dǎo)點(diǎn)選取距離對稱的點(diǎn)參與導(dǎo)數(shù)計算。而可變網(wǎng)格算法使用不同步長剖分模型,導(dǎo)致采用常規(guī)差分方法計算可變網(wǎng)格邊界區(qū)域的導(dǎo)數(shù)時會產(chǎn)生比較大的誤差。同時空間步長影響算法的數(shù)值頻散程度,也使得不同步長網(wǎng)格交界處由于算法問題而產(chǎn)生弱反射[18]。

如圖2所示處理方式,以x方向過渡區(qū)第n個點(diǎn)為例,當(dāng)大小網(wǎng)格步長比為α?xí)r,空間M階可變網(wǎng)格過渡差分格式為[18]:

(5)

同理可以得到z方向的過渡差分格式。

為實(shí)現(xiàn)二維可變網(wǎng)格計算,則需要解決x,z方向的混合偏導(dǎo)數(shù)過渡問題。M階精度差分格式計算混合偏導(dǎo)需要使用其周圍M2個點(diǎn)的函數(shù)值,對于高階差分而言,維持這個系統(tǒng)的對稱性求解過渡差分格式顯然非常復(fù)雜。考慮到求解二階混合偏導(dǎo)問題可以轉(zhuǎn)化為求解兩個嵌套的一階導(dǎo)數(shù)問題。具體推導(dǎo)過程參見附錄B。

圖1 可變網(wǎng)格剖分示意圖解

圖2 可變空間步長網(wǎng)格交界區(qū)域處理方式示意圖解(星形標(biāo)記代表導(dǎo)數(shù)計算點(diǎn))

僅在某一方向上改變網(wǎng)格步長,上述過渡處理方式不涉及插值計算。而在兩個空間維度改變網(wǎng)格步長時,由于計算過渡區(qū)導(dǎo)數(shù)值所需要的部分波場值不存在,需要插值得到。圖3以二維可變網(wǎng)格算法為例,標(biāo)出了二維可變網(wǎng)格計算時需要進(jìn)行插值的點(diǎn)(圖3中黑色圓點(diǎn)為需要進(jìn)行插值的點(diǎn))。

時間域顯示差分格式需要空間步長與時間采樣率滿足一定關(guān)系來保證差分算法的穩(wěn)定性,因此在調(diào)整空間網(wǎng)格步長時,時間采樣間隔必須滿足小網(wǎng)格的穩(wěn)定性條件。若模擬計算采用統(tǒng)一的時間采樣間隔,必然會造成大網(wǎng)格區(qū)的時間過采樣。因此需要在時間域改變采樣步長,小網(wǎng)格區(qū)加密時間片循環(huán),再對計算結(jié)果采用大時間步長重采樣,從而保證計算效率[18-19]。

圖3 可變網(wǎng)格計算中需要進(jìn)行插值的點(diǎn)

由于加密時間片循環(huán)會造成部分參與計算的時間片缺失,因此在改變時間步長時需要設(shè)置過渡帶,求出大時間步長循環(huán)下當(dāng)前時刻和下一時刻的波場值,插值得到小時間步長計算所需的各時刻波場值,如圖4所示。由于時間差分算子精度遠(yuǎn)小于空間差分算子精度,因此時間維度上的插值計算引入的誤差不會造成明顯的干擾現(xiàn)象。

圖4 可變時間步長過渡區(qū)處理方式示意圖解(“十”字形記號表示需要插值得到的點(diǎn))

2 模型試算

2.1 均勻介質(zhì)模型

首先建立均勻介質(zhì)模型,測試PML邊界的吸收效果與可變網(wǎng)格方法的適用性。

均勻雙相介質(zhì)模型的彈性參數(shù)與質(zhì)量參數(shù)如表1所示。模型大小為1km×1km;大網(wǎng)格空間步長為10m,時間步長為1.0ms;小網(wǎng)格空間步長為1m,時間步長為0.1ms。震源為等能量源,位于模型中心,使用主頻為35Hz的Ricker子波作為震源子波。將小網(wǎng)格區(qū)置于模型邊緣測試PML邊界的吸收效果與可變網(wǎng)格過渡區(qū)的處理效果。

圖5和圖6分別給出了該均勻雙相介質(zhì)模型180ms的波場快照和炮記錄。由圖5和圖6可以看到,雙相介質(zhì)中存在快縱波與慢縱波。固相介質(zhì)中的慢縱波能量遠(yuǎn)小于同條件下快縱波的能量,而液相介質(zhì)中的慢縱波能量與快縱波能量基本處于一個量級,但由于介質(zhì)的耗散作用,慢縱波能量衰減迅速,因此實(shí)際情況下很難觀測到慢縱波的存在。正演模擬中應(yīng)用的PML邊界條件可以很好地吸收傳播到邊界的快縱波與慢縱波,幾乎不會產(chǎn)生任何邊界反射,而將可變網(wǎng)格區(qū)域放置于模型角點(diǎn)處,使得部分小網(wǎng)格點(diǎn)參與計算PML邊界的波場,計算結(jié)果中,波場快照和炮記錄中都沒有明顯的波場異?，F(xiàn)象,說明該過渡方法不會引入明顯人為干擾,該模型驗(yàn)證了可變網(wǎng)格方法的可行性。

表1 均勻介質(zhì)模型雙相介質(zhì)參數(shù)

圖5 均勻雙相介質(zhì)模型180ms波場快照(黑色虛線框代表變網(wǎng)格區(qū)域)

圖6 均勻雙相介質(zhì)模型炮記錄

2.2 溶洞模型

使用一個溶洞模型來測試雙相介質(zhì)可變網(wǎng)格方法在精度與計算效率方面的優(yōu)勢。

溶洞模型如圖7所示,模型中涉及到的介質(zhì)參數(shù)如表2所示。均勻單相介質(zhì)背景下分布著兩個溶洞,其中較大溶洞半徑為50m,較小溶洞半徑為4m。溶洞內(nèi)部介質(zhì)為雙相介質(zhì)。模型大小為1km×1km;大網(wǎng)格空間步長為10m,時間采樣率為1.0ms;小網(wǎng)格空間步長為1m,時間采樣率為0.1ms?？勺兙W(wǎng)格區(qū)域如圖7中虛線框所示,圖中星形記號表示炮點(diǎn)位置。震源設(shè)置為組合縱波源,采用主頻為50Hz的Ricker子波作為震源子波。對比同一模型的大網(wǎng)格采樣模擬結(jié)果(圖8a)和小網(wǎng)格采樣模擬結(jié)果(圖8b)可以看到,小網(wǎng)格采樣不僅不會因?yàn)椴介L間距過大而錯過較小溶洞的樣點(diǎn),同時更小的采樣間隔也精細(xì)地刻畫了較大溶洞的邊緣輪廓特征,使得模擬結(jié)果更精確。

圖7 溶洞模型示意圖解

表2 溶洞模型雙相介質(zhì)參數(shù)

模型參數(shù)A/GPaN/GPaQ/GPaR/GPaρ11/(kg·m-3)ρ22/(kg·m-3)ρ12/(kg·m-3)B/(kg·m-3·s-1)溶洞內(nèi)介質(zhì)33．7821．113．1242．41301572．3512．5-307．50．15×105背景介質(zhì)194．80162．40002078．4000

圖8 溶洞模型采樣模擬結(jié)果

在其它模擬參數(shù)相同的條件下,針對同一模型,計算空間步長為10m,時間步長為1.0ms的常規(guī)大網(wǎng)格與空間步長為1m,時間步長為0.1ms的常規(guī)小網(wǎng)格模擬結(jié)果,并與可變網(wǎng)格方法模擬結(jié)果進(jìn)行比較與分析。

圖9,圖10和圖11分別為采用大網(wǎng)格采樣、可變網(wǎng)格采樣與小網(wǎng)格采樣正演模擬得到的波場快照與炮記錄。由于激發(fā)與接收都在單相介質(zhì)中進(jìn)行,因此這里只給出固相介質(zhì)中波場快照與單相介質(zhì)炮記錄結(jié)果。大網(wǎng)格采樣模擬由于采樣不足,較大溶洞邊緣呈階梯狀,造成模擬結(jié)果中存在很多角點(diǎn)的繞射波。同時溶洞內(nèi)部介質(zhì)波速較低,流體造成能量耗散,導(dǎo)致溶洞內(nèi)部數(shù)值頻散嚴(yán)重。由于大網(wǎng)格步長大于較小溶洞半徑,沒有采集到較小溶洞的樣點(diǎn),因此模擬結(jié)果中也沒有體現(xiàn)較小溶洞的地震響應(yīng)。

圖9 溶洞模型大網(wǎng)格采樣模擬結(jié)果

圖10 溶洞模型可變網(wǎng)格采樣模擬結(jié)果

可變網(wǎng)格采樣模擬方法對溶洞分布區(qū)域使用小步長模擬,由于針對溶洞的采樣步長較為精細(xì),溶洞邊緣較為平滑,有效降低了角點(diǎn)繞射波的能量,同時小步長計算也有效壓制了大溶洞內(nèi)部低速介質(zhì)的數(shù)值頻散。由于減小了采樣步長,較小溶洞的響應(yīng)也在模擬結(jié)果中體現(xiàn)出來?？勺兙W(wǎng)格方法計算的結(jié)果幾乎與整個模型均使用小網(wǎng)格模擬的結(jié)果(圖11)完全一致,從而驗(yàn)證了可變網(wǎng)格方法模擬的精確性。

表3給出了對溶洞模型采用3種不同模擬方法消耗的機(jī)時與占用的內(nèi)存情況。正演模擬使用的機(jī)型為HP Compaq dx2318MT,CPU為Intel CPU E8400,主頻為3.0GHz,內(nèi)存為4G,在Window XP系統(tǒng)下運(yùn)行。由統(tǒng)計結(jié)果可以看到,整個模型使用小網(wǎng)格采樣計算,計算機(jī)時是大網(wǎng)格采樣計算的570倍,內(nèi)存消耗是大網(wǎng)格采樣模擬的88倍,而可變網(wǎng)格模擬的計算機(jī)時約為大網(wǎng)格采樣計算機(jī)時的5倍,內(nèi)存消耗約為大網(wǎng)格采樣模擬的2倍,相對于小網(wǎng)格采樣模擬大大提高了計算效率。

圖11 溶洞模型小網(wǎng)格采樣模擬結(jié)果

表3 溶洞模型3種模擬方法計算成本統(tǒng)計

模擬方法計算機(jī)時/s占用內(nèi)存/KB空間步長為10m，時間步長為1．0ms的大網(wǎng)格常規(guī)有限差分方法3811732空間步長為1m，時間步長為0．1ms的小網(wǎng)格常規(guī)有限差分方法216431035240可變網(wǎng)格有限差分方法20225868

3 結(jié)束語

本文將可變網(wǎng)格方法應(yīng)用于基于Biot理論的雙相介質(zhì)彈性波正演模擬,解決了空間二階混合偏導(dǎo)數(shù)的可變網(wǎng)格過渡問題,并得到了可變網(wǎng)格過渡區(qū)PML差分格式,從而實(shí)現(xiàn)了任意倍數(shù)、任意位置變步長的可變網(wǎng)格有限差分正演模擬方法。該方法可以針對模型中重點(diǎn)區(qū)域進(jìn)行精細(xì)刻畫,有效降低矩形網(wǎng)格采樣率低引起的人為散射問題,壓制由于大網(wǎng)格離散引起的數(shù)值誤差。文中多個模型驗(yàn)證了基于PML邊界條件的孔隙介質(zhì)彈性波可變網(wǎng)格正演模擬方法是一種精確、高效的雙相介質(zhì)彈性波正演模擬方法。將該方法應(yīng)用于復(fù)雜精細(xì)雙相介質(zhì)模型中,能夠大大提高正演模擬效率,對研究儲層流體正演響應(yīng)具有積極的作用。同時該方法能夠推廣至其它復(fù)雜波動理論,值得進(jìn)一步深入研究。

[1] 楊頂輝,張中杰,滕吉文,等.雙相各向異性研究、問題與應(yīng)用前景[J].地球物理學(xué)進(jìn)展,2000,15(2):7-21 Yang D H,Zhang Z J,Teng J W.The study of two-phase anisotropy questions and applied prospects[J].Progress in Geophysics,2000,15(2):7-21

[2] Biot M A.Theory of propagation of elasticwaves in a fluid-saturated porous solid,I:low-frequency range[J].The Journal of the Acoustical Society of America,1956,28(2):16-178

[3] Biot M A.Theory of propagation of elastic waves in a fluid-saturated porous solid,II:higher-frequency range[J].The Journal of the Acoustical Society of America,1956,28(2):179-191

[4] Biot M A.Mechanics of deformations and acoustic propagation in porous media[J].Journal of Applied Physics,1962,33(4):1482-1498

[5] Biot M A.Generalized theory of acoustic propagation in porous dissipative media[J].The Journal of the Acoustical Society of America,1962,34(9):1254-1264

[6] 趙海波,王秀明,王東,等.完全匹配層吸收邊界在孔隙介質(zhì)彈性波模擬中的應(yīng)用[J].地球物理學(xué)報,2007,50(2):581-591 Zhao H B,Wang X M,Wang D,et al.Applications of the boundary absorption using a perfectly matched layer for elastic wave simulation inporoelastic media[J].Chinese Journal of Geophysics,2007,50(2):581-591

[7] Alford R M,Kelly K R,Boore D M.Accuracy of finite difference modeling of the acoustic wave equation[J].Geophysics,1974,39(6):834-842

[8] Moczo P.Finite-difference technique for SH waves in 2D media,using irregular girds:application to the seismic response problem[J].Geophysical Journal International,1989,99(2):321-329

[9] 黃超,董良國.可變網(wǎng)格與局部時間步長的高階差分地震波數(shù)值模擬[J].地球物理學(xué)報,2009,52(1):176-186 Huang C,Dong L G.High-order finite-difference method in seismic wave simulation with variable grids and local time-steps[J].Chinese Journal of Geophysics,2009,52(1):176-186

[10] 孫林潔,印興耀.基于PML邊界條件的高倍可變網(wǎng)格有限差分?jǐn)?shù)值模擬方法[J].地球物理學(xué)報,2011,54(6):1614-1623 Sun L J,Yin X Y.A finite difference scheme based on PML boundary condition with high power grid step variation[J].Chinese Journal of Geophysics,2011,54(6):1614-1623

[11] 孫成禹,李勝軍,倪長寬,等.波動方程變網(wǎng)格步長有限差分?jǐn)?shù)值模擬[J].石油物探,2008,47(2):123-128 Sun C Y,Li S J,Ni C K,et al.Elastic wave modeling on a variable gird finite-difference method[J].Geophysical Prospecting for Petroleum,2008,47(2):123-128

[12] 黃超,董良國.可變網(wǎng)格與局部時間步長的交錯網(wǎng)格高階差分彈性波模擬[J].地球物理學(xué)報,2009,52(11):2870-2878 Huang C,Dong L G.Staggered-grid high-order finite-difference method in elastic wave simulation with variable grids and local time-steps[J].Chinese Journal of Geophysics,2009,52(11):2870-2878

[13] 朱生旺,曲壽利,魏修成,等.變網(wǎng)格有限差分彈性波動方程數(shù)值模擬方法[J].石油地球物理勘探,2007,42(6):634-639 Zhu S W,Qu S L,Wei X C,et al.Numeric simulation by grid-various finite-difference elastic wave equation[J].Oil Geophysical Prospecting,2007,42(6):634-639

[14] 李振春,李慶洋,黃建平,等.一種穩(wěn)定的高精度雙變網(wǎng)格正演模擬與逆時偏移方法[J].石油物探,2014,53(2):127-136 Li Z C,Li Q Y,Huang J P.A stable and high-precision dual-variable grid forward modeling and reverse time migration method[J].Geophysical Prospecting for Petroleum,2014,53(2):127-136

[15] 王守東.聲波方程完全匹配層吸收邊界[J].石油地球物理勘探,2003,38(1):31-34 Wang S D.Absorbing boundary condition for acoustic wave equation by perfectly matched layer[J].Oil Geophysical Prospecting,2003,38(1):31-34

[16] 陳可洋.聲波完全匹配層吸收邊界條件的改進(jìn)算法[J].石油物探,2009,48(1):76-79 Chen K Y.The improvement of perfect matched layer absorbs condition[J].Geophysical Prospecting for Petroleum,2009,48(1):76-79

[17] 陳可洋.完全匹配層吸收邊界條件研究[J].石油物探,2010,49(5):472-477 Chen K Y.The research of perfect matched layer condition[J].Geophysical Prospecting for Petroleum,2010,49(5):472-477

[18] 吳國忱,王華忠.波場模擬中的數(shù)值頻散分析與校正策略[J].地球物理學(xué)進(jìn)展,2005,20(1):58-65 Wu G C,Wang H Z.Analysis of numerical dispersion in wave-field simulation[J].Progress in Geophysics,2005,20(1):58-65

[19] 劉洋,李承楚,牟永光.任意偶數(shù)階精度有限差分法數(shù)值模擬[J].石油地球物理勘探,1998,33(1):1-10 Liu Y,Li C C,Mou Y G.Finite-difference numerical modeling of any even-order accuracy[J].Oil Geophysical Prospecting,1998,33(1):1-10

附錄A

加入衰減項(xiàng)的PML方程為:

(A1)

將時間導(dǎo)數(shù)的二階中心差分格式代入公式(A1):

(A2)

整理并求解這個方程組,得到雙相介質(zhì)彈性波PML起始方程:

(A3)

其中,

(A4)

將方程(A4)依照求導(dǎo)方向拆分,加入不同方向的衰減系數(shù),利用二階空間偏導(dǎo)的M階差分格式進(jìn)行離散近似即可得到雙相介質(zhì)彈性波方程的PML差分格式。

以x方向?yàn)槔?空間二階偏導(dǎo)數(shù)的M階精度常規(guī)差分格式為[15]:

(A5)

式中:cm為差分系數(shù)。

空間二階混合偏導(dǎo)數(shù)的M階精度常規(guī)差分格式為[17]:

(A6)

對于x方向,有:

(A7)

其中,

(A8)

對z方向也有相似的結(jié)論。

對于x,z方向:

(A9)

其中,

(A10)

附錄B

以x方向過渡區(qū)的第n個點(diǎn)為例,當(dāng)大小網(wǎng)格步長比為α?xí)r,空間M階可變網(wǎng)格差分格式為[9]:

(B1)

同理得到z方向的過渡差分格式。

x,z方向的混合偏導(dǎo)數(shù)的過渡問題可以轉(zhuǎn)化為求解兩個嵌套的一階導(dǎo)數(shù)問題,即:

(B2)

同理,過渡區(qū)的二階混合偏導(dǎo)數(shù)也可以通過兩次嵌套的一階導(dǎo)數(shù)過渡差分格式求解。即:

(B3)

則基于PML邊界過渡區(qū)的時間二階、空間M階精度聲波方程可變網(wǎng)格差分格式為:

(B4)

其中,

(B5a)

(B5b)

(B5c)

(B5d)

(B5e)

(B5f)

(B5g)

(B5h)

(B5i)

(B5j)

(B5k)

(B5l)

(B5m)

(編輯：陳 杰)

A variable grid finite-difference scheme with PML boundary condition in elastic wave of porous media

Sun Linjie,Liu Chuncheng,Zhang Shixin

(CNOOCResearchInstitute,Beijing100028,China)

In order to improve the accuracy of forward modeling in porous media,a grid variation finite-difference simulation method is applied to Biot’s two-phase elastic wave equation.The difference approximation schemes in grid-variation zone for second-order mixed partial derivative are solved and the variable-grid finite difference scheme with PML boundary is derived,which realizes the two-phase media elastic wave variable-grid forward modeling with integral-multiple variable spatial and time grid steps.Numerical tests show that the grid step variation method would improve the simulation accuracy by using smaller gird step to resample the subtle geological structure and the area with rapid change of geological conditions.Moreover,the algorithm can also ensure the computational efficiency by suppressing the sampling rate to the region with relatively stable model parameters.Compared with conventional forward modeling method,the numerical testing result indicates that the two-phase variable-grid forward modeling method obviously improves the simulation accuracy and efficiency.

Biot theory,porous media,PML boundary condition,variable grid,finite difference

2015-03-30;改回日期：2015-07-08。

孫林潔(1987—),女,工程師,主要從事地震儲層預(yù)測研究工作。

國家科技重大專項(xiàng)“海洋深水區(qū)油氣勘探關(guān)鍵技術(shù)”項(xiàng)目(2011ZX05025-001)資助。

P631

A

1000-1441(2015)06-0652-13

10.3969/j.issn.1000-1441.2015.06.003