陳 豪，李春香，魏艷敏

（貴州工程應(yīng)用技術(shù)學(xué)院礦業(yè)工程學(xué)院，貴州 畢節(jié) 551700）

地球物理問(wèn)題分為正演問(wèn)題和反演問(wèn)題。前者是已知地下形體幾何參數(shù)、密度、速度分布等，通過(guò)理論計(jì)算、模擬計(jì)算或者模型試驗(yàn)等方法，得到地下形體的異?；蛘呃碚搱?chǎng)值，通常稱之為“由源求場(chǎng)”；后者是已知地下的異?；驁?chǎng)值，求地下形態(tài)或地下形體的幾何參數(shù)等，稱為“由場(chǎng)求源”。

隨著礦產(chǎn)、油氣的日益開(kāi)采，淺部的資源幾近消耗殆盡，中深部礦產(chǎn)資源勘探已成為當(dāng)下及未來(lái)的發(fā)展趨勢(shì)。數(shù)值模擬研究能很好地探索各方法的地球物理響應(yīng)。具體到地震勘探，三維數(shù)值模擬較之二維能更好更全面地探索地震波場(chǎng)特征和傳播規(guī)律，從而有效指導(dǎo)地震數(shù)據(jù)采集、處理以及解釋。

井間地震是一種較新的地震勘探方法，其地震波場(chǎng)非常豐富，也非常復(fù)雜，這是井間地震的優(yōu)勢(shì)，同時(shí)也是井間地震的難點(diǎn)。井下致密地層對(duì)人工地震信號(hào)衰減較小，因此井間地震能夠接收到較高頻率的地震波，從而具有了高分辨率的優(yōu)勢(shì)，成為連接地面地震、測(cè)井和油藏地質(zhì)的橋梁。

常用的數(shù)值模擬方法有有限差分法、偽譜法、有限元法、傅里葉變換法等。由于有限差分法是用離散值表示連續(xù)值，用該方法進(jìn)行數(shù)值模擬，不可避免地出現(xiàn)高頻震蕩的問(wèn)題，編程中應(yīng)控制好時(shí)間步長(zhǎng)、網(wǎng)格間距等條件。相較于其它數(shù)值模擬方法，有限差分法具有推導(dǎo)過(guò)程簡(jiǎn)單、計(jì)算速度快、占用計(jì)算機(jī)內(nèi)存小等優(yōu)勢(shì)。

傳統(tǒng)的地震勘探是把地層當(dāng)作完全彈性介質(zhì)，但實(shí)際處理地震資料時(shí)發(fā)現(xiàn)，隨著傳播距離的增加，地震波頻率和振幅都會(huì)衰減，說(shuō)明實(shí)際地層并非完全弾性的。研究表明，實(shí)際地層更接近于粘彈性介質(zhì)。[1]關(guān)于地震波能量的衰減機(jī)理，當(dāng)前普遍的說(shuō)法是固體顆粒內(nèi)摩擦的熱效應(yīng)和博伊托（Boit）噴射流效應(yīng)，前者導(dǎo)致地震波的能量在傳播過(guò)程中以熱量的形式發(fā)生衰減和耗散，后者導(dǎo)致地震波的振幅產(chǎn)生衰減和速度的頻散。[2]

斯托克斯（Stockes）最早在1845年對(duì)粘彈性介質(zhì)中的地震波展開(kāi)了研究。國(guó)內(nèi)學(xué)者在粘彈性介質(zhì)方面也做了大量研究，為構(gòu)建相應(yīng)的模型來(lái)近似匹配地層介質(zhì)，有學(xué)者總結(jié)了開(kāi)爾芬（Kelvin）、達(dá)朗貝爾（D’Alembert）、玻爾茲曼（Boltzmann）和麥克斯韋爾（Maxwell）等幾種粘彈性介質(zhì)模型，給出了各種粘彈性介質(zhì)模型與完全彈性介質(zhì)模型之間彈性參數(shù)的對(duì)應(yīng)關(guān)系。[2]單啟銅等用偽譜法模擬了二維粘彈性介質(zhì)的地震波場(chǎng)，并討論分析了粘滯系數(shù)與頻率、振幅間的關(guān)系[3]；邵志剛等用有限差分法模擬了二維線性粘彈介質(zhì)中地震波場(chǎng)，并討論分析了相速度與粘滯系數(shù)的關(guān)系[4]；金勝文等用有限差分法在x-t域數(shù)值模擬了三維波動(dòng)方程[5]；王德利等用一種基于信息傳遞接口（MPI）的并行有限差分算法數(shù)值模擬了三維粘彈介質(zhì)的地震波場(chǎng)，該算法大大提高了三維數(shù)值模擬的效率[6]。本文用交錯(cuò)網(wǎng)格有限差分法在時(shí)間2階、空間8階有限差分法數(shù)值模擬開(kāi)爾芬粘彈性介質(zhì)，給定不同的品質(zhì)因子通過(guò)均勻模型試算，對(duì)記錄的Z分量進(jìn)行振幅和頻譜分析，得到了若干發(fā)現(xiàn)。

1 方法原理

在交錯(cuò)網(wǎng)格高階有限差分法中，波動(dòng)方程中位移速度、應(yīng)力變量分別定義在不同的半網(wǎng)格節(jié)點(diǎn)或整網(wǎng)格節(jié)點(diǎn)上。如圖1所示。

圖1 三維交錯(cuò)網(wǎng)格示意圖

2 數(shù)值模擬中的幾個(gè)問(wèn)題

2.1 吸收邊界條件

為消除數(shù)值模擬中人為引入邊界反射的影響，仍然采用完全匹配層（PML）處理二維邊界的思想來(lái)解決本文三維邊界反射的問(wèn)題。[8]即方程中任何一個(gè)變量都可以分為垂直于該面的一個(gè)分量和平行于該面的另外兩個(gè)分量，那么對(duì)垂直于該面的分量進(jìn)行衰減處理，平行于該面的另外兩個(gè)分量不做任何處理，而公共的8個(gè)小立方體邊界都要做衰減處理。用二維完全匹配層的思想及方程中不同變量定義在網(wǎng)格中的節(jié)點(diǎn)，就可導(dǎo)出吸收邊界的離散表達(dá)式。

2.2 震源加載

數(shù)值模擬時(shí)使用的震源激發(fā)方式不同，得到的波場(chǎng)也不相同。本文采用的激發(fā)方式為縱波震源，即將震源子波按波場(chǎng)模擬時(shí)間步長(zhǎng)進(jìn)行離散，當(dāng)震源波場(chǎng)向前計(jì)算一個(gè)時(shí)間步長(zhǎng)時(shí)，就在震源位置加入震源子波在該時(shí)刻的離散值的震源加載方式。采用零相位雷克（Ricker）子波作為震源的子波，其函數(shù)表達(dá)式為：

2.3 數(shù)值頻散

交錯(cuò)網(wǎng)格高階差分算法用有限的差分離散來(lái)代替微分方程，難免存在誤差項(xiàng)，因而模擬得到的結(jié)果不可避免地出現(xiàn)數(shù)值頻散。為了解決數(shù)值模擬所引起的頻散問(wèn)題，前人做了大量研究工作，主要有盡量減少空間網(wǎng)格間距、降低子波主頻、減少時(shí)間步長(zhǎng)、提高差分階數(shù)等。一般情況下，通過(guò)減小網(wǎng)格的空間間隔來(lái)降低數(shù)值頻散，每個(gè)最短波長(zhǎng)上離散網(wǎng)格點(diǎn)數(shù)是5～7個(gè)，就可以滿足計(jì)算的要求。[7]由于所研究的是井間地震波場(chǎng)數(shù)值模擬，主頻降得太低則與實(shí)際不相符，從而失去模擬的意義。因此，在保證計(jì)算量的前提下減少網(wǎng)格間的距離，用時(shí)間2階、空間8階的高差分格式來(lái)盡量降低數(shù)值頻散問(wèn)題。

3 模型試算

為了驗(yàn)證方法的正確性，設(shè)計(jì)了一個(gè)200m*200m*200m的三維均勻模型，直井接收，炮點(diǎn)坐標(biāo)S（50m,100m,100m），接收井井口坐標(biāo)（150m,100m,0m），網(wǎng)格間距、道間距皆選擇2.5m，共設(shè)置81道檢波器接收，邊界厚度為50m，采樣率為2.5ms，震源子波主頻為250Hz，取縱波速度為3000m/s，常密度為2.3g/cm3，取品質(zhì)因子QP=120、QS=100。得到如圖2所示的50ms時(shí)刻波場(chǎng)快照和如圖3所示的地震記錄。

圖2 50毫秒時(shí)刻均勻模型地震波場(chǎng)快照

圖3 均勻模型地震波場(chǎng)記錄

從圖2的波場(chǎng)快照中可以看到邊界得到了很好的吸收，縱波和橫波信息也能清晰顯示。由于XOZ面垂直于Vy分量，所以在XOZ面切片的快照中Vy分量只能顯示出橫波信息，YOZ面切片中的Vx分量、XOY面切片中的Vz分量亦是如此。

為探究粘彈性介質(zhì)中品質(zhì)因子與振幅、頻率的關(guān)系，依然使用上述模型，在其它參數(shù)不變的條件下，通過(guò)給定不同的品質(zhì)因子進(jìn)行模擬。抽取Z分量的第60道200ms-1000ms的記錄作振幅對(duì)比，結(jié)果如圖4所示。

圖4 不同品質(zhì)因子第60道地震記錄振幅對(duì)比

從圖中可以看到，隨著品質(zhì)因子的減小，地震記錄的振幅隨之降低，曲線變得更加光滑，頻散也隨之減弱，該結(jié)果與在完全匹配層吸收邊界條件下二維粘彈性介質(zhì)波場(chǎng)模擬的已有研究結(jié)論一致[3]，說(shuō)明粘彈性介質(zhì)對(duì)地震波振幅和能量存在吸收衰減的作用。可見(jiàn)，當(dāng)品質(zhì)因子取值較大時(shí)，相當(dāng)于沒(méi)有考慮介質(zhì)對(duì)地震波能量的吸收衰減，粘彈性介質(zhì)趨近于完全彈性介質(zhì)。

對(duì)上述模型地震記錄的Z分量進(jìn)行頻譜分析，結(jié)果如圖5所示。

圖5 不同品質(zhì)因子Z分量地震振幅譜對(duì)比

圖5中，較高曲線表示Z分量記錄整個(gè)81道的平均振幅譜，較低曲線表示第40道的振幅譜。從圖中可以看到，當(dāng)QP=500、QS=400（粘滯性很弱）時(shí)，主頻約為80Hz，頻帶范圍在40Hz到120Hz之間，頻帶寬約80，最大振幅值約73；當(dāng)QP=120、QS=100（粘滯性較強(qiáng)）時(shí)，主頻約為66Hz，頻帶范圍在40Hz到100Hz之間，頻帶寬60，最大振幅值約65；當(dāng)QP=80、QS=60（粘滯性很強(qiáng)）時(shí)，主頻約為57Hz，頻帶范圍在35Hz到90Hz之間，頻帶寬約55，最大振幅值約49。所以經(jīng)上述比較可以發(fā)現(xiàn)，隨著粘彈性介質(zhì)品質(zhì)因子的減小，主頻逐漸降低，振幅逐漸降低，頻帶范圍也逐漸變窄。對(duì)地震記錄的X分量、Y分量做頻譜分析，也能得到同樣的發(fā)現(xiàn)。

4 結(jié)論和認(rèn)識(shí)

本文使用時(shí)間2階、空間8階交錯(cuò)網(wǎng)格有限差分法模擬了三維開(kāi)爾芬粘彈性介質(zhì)的地震波場(chǎng)，結(jié)果表明在合適的網(wǎng)格間距和時(shí)間步長(zhǎng)等條件下，該方法可以滿足三維地震波場(chǎng)的數(shù)值模擬。但必須指出，由于三維地震波場(chǎng)數(shù)值模擬數(shù)據(jù)量龐大，為壓制頻散，本文所建立模型還是稍小，為壓制數(shù)值頻散，數(shù)值模擬時(shí)間還是稍短，總體計(jì)算速度也較緩慢。因此加快計(jì)算效率，并行算法或是更好的選擇。

粘彈性介質(zhì)模型有開(kāi)爾芬、達(dá)朗貝爾、玻爾茲曼和麥克斯韋爾等幾種介質(zhì)模型，本文僅對(duì)三維開(kāi)爾芬粘彈性介質(zhì)進(jìn)行了初步定性的模擬研究。為進(jìn)一步探索地震波在粘彈性介質(zhì)中的傳播規(guī)律，應(yīng)加深以上幾種粘彈性介質(zhì)模型的對(duì)比研究，以探究補(bǔ)償實(shí)際地震資料處理中粘滯性的影響因素，從而達(dá)到提高地震分辨率的目的。