粘彈性介質(zhì)地震波傳播的褶積微分算子法數(shù)值模擬研究①

賀同江，劉紅艷，李小凡

（1.天津市地震局，天津300201；2.中國科學(xué)院地質(zhì)與地球物理研究所，北京100029）

粘彈性介質(zhì)地震波傳播的褶積微分算子法數(shù)值模擬研究①

賀同江1，劉紅艷1，李小凡2

（1.天津市地震局，天津300201；2.中國科學(xué)院地質(zhì)與地球物理研究所，北京100029）

早期的褶積微分算子法都是基于正反傅立葉變換而實現(xiàn)的，其精度比四階有限差分稍高。本文將計算數(shù)學(xué)中的Forsyte廣義正交多項式微分算子與褶積算子相結(jié)合，構(gòu)建了一個新的快速、高精度褶積微分算子，其計算結(jié)果非常接近實驗函數(shù)微分的精確值，精度與l6階有限差分相當。粘彈性波動方程更真實地描述了實際地下介質(zhì)中彈性波的傳播規(guī)律及其波場特征。本文以二維粘彈性波動方程為例，推導(dǎo)了粘彈性介質(zhì)波動方程的離散格式，用迭積微分算子法實現(xiàn)了粘彈性介質(zhì)的地震波場正演模擬，并對其波傳播特征進行了分析。計算結(jié)果表明該算法能正確模擬粘彈性介質(zhì)中的地震波，正確地反映粘彈性介質(zhì)中波場的傳播規(guī)律。

褶積微分算子；粘彈性介質(zhì)；地震波；數(shù)值模擬

Abstract：Early convolutional differentiators are all based on the Fourier transformation，their precision is a little higher than that of four－order finite difference.For improving the precision and efficiency of seismic modeling，a new modeling approach（Convolutional Forsyte Polynomial Differentiator Method，CFPD）is developed in this paper by using optimized convolutional operators for spatial differentiation and staggered－grid finite－difference for time differentiation in wave equation computation.The solution of this new method is much close to the exact value，and the precision is nearly equal to that of 16－order finite difference.Viscous－elastic wave equation can better explain wave phenomena of the true earth media.This paper applies the CFPD method to model 2－D viscous－elastic seismic wave field for the first time，and further derives the first－order velocity－stress discreate equations for viscous－elastic media.The numerical results show that the algorithm can bring reliable outcomes with high precision and fast speed，and viscous－elastic modeling is much efficient.

Key words：Convolutional Forsyte Polynomial Differentiator；Viscous－elastic media；Seismic wave；Numerical simulation

0 引言

眾所周知，在討論地震波傳播時絕大部分情況是把地震波看作在理想彈性介質(zhì)中傳播的。然而實際工作中所獲得的似正弦狀地震記錄與經(jīng)典彈性理論所預(yù)言的脈沖狀地震記錄之間存在巨大差異，人們發(fā)現(xiàn)實際巖層對在其中傳播的地震波有吸收作用，吸收激發(fā)脈沖波的某些頻譜使其能量發(fā)生損耗，因此實際地球介質(zhì)性質(zhì)更接近于粘彈性體。這對于解決油氣勘探面臨的問題，尤其是對于提高勘探精度都具有重要的意義［1］。近些年來圍繞粘彈介質(zhì)地震波場數(shù)值模擬的研究也為數(shù)不少［2－10］。迄今為止，地震波場正演模擬的數(shù)值算法有很多，每一種方法都存在其自身的優(yōu)越性與局限性［11－13］。盡管褶積微分算子不是一個新的概念，但早期的相關(guān)方法［14－18］都是基于傅氏變換的褶積微分算子，其精度僅稍高于四階有限差分的精度。本文是以計算數(shù)學(xué)中的Forsyte廣義正交多項式插值函數(shù)［19］為基礎(chǔ)，構(gòu)建一個新的褶積微分算子，由于該褶積微分算子具有相當與16階有限差分的高精度和短算子低階有限差分算法的高速度，故該方法已被廣泛應(yīng)用于復(fù)雜介質(zhì)的地震波場數(shù)值模擬中［20－24］。本文將該算子引入到粘彈性波一階速度－應(yīng)力方程的空間微分運算中去，采用時間錯格有限差分算子替代普通的差分算子以匹配高精度的褶積微分算子，從而構(gòu)造一種全新的粘彈性介質(zhì)地震波場數(shù)值模擬方法。本文重點研究二維粘彈性介質(zhì)中一階雙曲型粘彈性波方程基于Forsyte廣義正交多項式迭積微分算子法的數(shù)值模擬，首次把褶積微分算子法引入到粘彈性介質(zhì)的模擬中來，通過構(gòu)造不同的模型表明該方法能成功地模擬粘彈性介質(zhì)。

1 Forsyte廣義正交多項式迭積微分算子求導(dǎo)原理

本文利用錯格有限差分褶積微分算子法來計算粘彈性介質(zhì)中的地震波場，其主導(dǎo)思想是：利用基于Forsyte廣義正交多項式褶積微分算子有效表示計算波場對空間的偏導(dǎo)數(shù)，采用錯格有限差分法計算對時間的偏導(dǎo)數(shù)，其計算思路類似于偽譜法［20］。

首先我們給出Forsyte多項式微分算子。Forsyte多項式是一個廣義正交多項式，其插值函數(shù)可寫為

其中

f（xi）為被插值函數(shù)f（x）在點xi處的值。式（1）中的P0（x），……Pj＋1（x）定義為Forsyte多項式系統(tǒng)。

對式（1）中的x求導(dǎo)，可得

其中

Forsyte多項式微分算子可寫為

其中

此處的Cj不同于式（2）中的Cj。

將上述微分算子（3）離散化可得

這里，i為采樣指標；Δx為沿著x軸的采樣間隔。就實際應(yīng)用而言，須將微分算子截成短算子，勢必引起Gibbs現(xiàn)象。另外，多項式的引入還將引起Runge現(xiàn)象。為了消除這些現(xiàn)象，必須采用窗函數(shù)以截斷長微分算子。本文采用的是下列Gaussian窗函數(shù)：

其中mx為單邊截斷長度的采樣數(shù)；c為常數(shù)；a（0.1≤a≤0.75）為衰減因子。將微分算子（4）用式（5）截斷并鋸齒化后，可得如下實用的一階迭積微分算子：

類似地，二階迭積微分算子可寫為

通過對算子長度的調(diào)節(jié)及算子系數(shù)的優(yōu)化，可同時兼顧波場解的全局信息與局部信息。本文運用算子長度為9點的一階迭積微分算子求解波動方程。9點迭積算子的最優(yōu)權(quán)系數(shù)為－9.135 778 624 087 487E－004，6.215 276 667 089 320E－003，－2.450 472 852 706 014E－002，8.399 311 520 683 385E－002，0.000 000 000 000 000E＋000，－8.399 311 520 683 385E－002，2.450 472 852 706 014E－002，－6.215 276 667 089 320E－003，9.135 778 624 087 487E－004?？梢钥闯?，算子的系數(shù)是反對稱的。

2 2－D標準線性粘彈性體波動方程及其離散格式

標準線性粘彈性體2－D波動方程形式如下［3］：

對模型區(qū)間離散后，設(shè)n，m，k分別是沿著空間x軸、z軸，時間t軸的采樣點數(shù)，Δx、Δz、Δt分別是沿x軸、z軸、t軸的采樣間隔，mx、mz是沿x軸、z軸采樣數(shù)的半算子長度，則式（8）二維非均勻各向同性粘彈性波方程，其離散化的時間錯格差分空間迭積微分算子法格式為（體力為零）

其中應(yīng)力σxz，σzz的離散格式同應(yīng)力σxx的離散格式；記憶變量rzz，rxz的離散格式同記憶變量rxx的離散格式。

3 穩(wěn)定性條件及邊界條件

通過分析可知，式（9）所示的基于Forsyte廣義正交多項式褶積微分算子法的穩(wěn)定性條件類似于偽譜法的穩(wěn)定性條件［25］

其中Vmax為介質(zhì)速度最大值；Δt為時間步長

同理，用于偽譜法的Cerjan吸收邊界條件也同樣適用于基于Forsyte廣義正交多項式褶積微分算子法

該吸收邊界方法由Cerjon［26］、Kosloff［27］提出，其主導(dǎo)思想就是給模型的有效區(qū)域的邊界上加粘滯層，使得波的能量在粘滯層內(nèi)最大程度地衰減，從而消除邊界反射。式（12）中，α為衰減系數(shù)，N為粘滯層的網(wǎng)格點數(shù)。選取α＝0.015，N＝20。

4 數(shù)值實驗

模型如圖1所示為一個斷層模型，模型的網(wǎng)格剖分為256×256，以網(wǎng)格點為模型坐標，網(wǎng)格間距Δx＝Δz＝10m，采樣間隔Δt＝1ms。震源為45°傾斜集中力源，震源采用25Hz的Ricker子波，坐標為（128，60）。垂直斷層構(gòu)造的高度為50×10m，z＝100×10m為斷層結(jié)構(gòu)的上邊界，z＝150×10m為斷層結(jié)構(gòu)的下邊界，A點的坐標為（128，100），B點的坐標為（128，150）。模型參數(shù)如表1所示，彈性波場模擬采用表1未加入粘滯系數(shù)的參數(shù)。

表1 斷層模型參數(shù)

圖1 斷層模型Fig.1 Fault model.

圖2為粘彈性介質(zhì)和彈性介質(zhì)在500ms時刻兩個分量的波場快照，圖3為合成記錄圖，接收道位于z＝50的水平位置?？梢钥闯?，波在非均勻介質(zhì)中的傳播較為復(fù)雜：震源在上層中激發(fā)，產(chǎn)生了直達波P1、S1，直達波傳到斷層界面時發(fā)生了反射透射現(xiàn)象，從而產(chǎn)生了同類反射波P1P1、S1S1、轉(zhuǎn)換反射波P1S1、S1P1、同類透射波S1S2、P1P2、以及轉(zhuǎn)換透射波P1S2，在拐點A，B處產(chǎn)生繞射波。粘彈性介質(zhì)的波場快照波的震相與彈性介質(zhì)的一樣，從波場快照看出第二層介質(zhì)的品質(zhì)因子比較大，相對彈性介質(zhì)來說第二層介質(zhì)產(chǎn)生的透射波在粘彈性介質(zhì)中的衰減比較明顯，這與以前研究得出的結(jié)論一致［3－6］。數(shù)值實驗表明，基于Forsyte廣義正交多項式褶積微分算子法能正確模擬粘彈性介質(zhì)中的地震波。從圖2﹑圖3上看到比較明顯的邊界反射波，因此對于較復(fù)雜的非均勻模型，Cerjan吸收邊界條件的反射波吸收效果較差。

圖2 500ms時刻在粘彈性介質(zhì)和彈性介質(zhì)的波場快照Fig.2 Wave－field snapshots in visco－elastic media and elastic media（t＝500ms）.

圖4為（60，150）點處記錄到的水平分量和垂直分量進行了彈性波與粘彈性波對比分析，可見粘彈性波波場記錄的振幅能量明顯減弱，反射波波形和相位畸變嚴重，記錄中與S波相關(guān)的波，吸收和頻散都比P波嚴重。

從以上的模擬結(jié)果可以看出，褶積微分算子法有效地模擬了標準線性粘彈介質(zhì)中的地震波的傳播。構(gòu)造不同模型，地震記錄反應(yīng)了品質(zhì)因子的粘彈性衰減作用，S波的影響比P波的影響明顯。

5 結(jié)論

本文首次把基于Forsyte廣義正交多項式褶積微分算子法引入到粘彈性介質(zhì)的地震波場數(shù)值模擬中。數(shù)值試驗表明，該褶積微分算子法的正演結(jié)果正確，是一種快速的地震波場模擬方法；并且由于該方法在空間域可同時兼顧波場解的全局信息與局部信息，因而對于復(fù)雜介質(zhì)模型該數(shù)值模擬方法在空間域中能獲得較高的分辨率，精度較高。從而也證明了褶積微分算子法能成功的模擬粘彈性介質(zhì)的地震波場，是褶積微分算子法應(yīng)用的推廣，為褶積微分算子法的推出及后續(xù)研究的成功開展，將為高精度地震波模擬、地震波偏移、地震反演、地震波成像及地震波在復(fù)雜非均勻介質(zhì)中的傳播等研究提供更為廣泛的選擇。粘彈性介質(zhì)中的波的傳播更加接近地下的真實情況，通過對粘彈性介質(zhì)波場數(shù)值模擬，可以清楚地了解介質(zhì)的粘彈性對地震波傳播的影響。但本方法的邊界問題有待解決。

圖3 斷層模型波場合成記錄Fig.3 Synthetic record of the fault model in visco－elastic and elastic media.

圖4 點（60，150）處彈性波記錄和粘彈性波記錄對比（sls－粘彈性介質(zhì)，elastic－彈性介質(zhì)）Fig.4 Comparison of the record of elastic media and visco－elastic media at（60，150）point.

［1］ ［美］N H瑞克，著.許云，譯.粘彈性介質(zhì)中的地震波［M］.北京：地質(zhì)出版社，1981.

［2］ Carcione J M，D Kosloff，R Kosloff.Wave propagation simulation in a linear viscoacoustic medium［J］.Geophys.J.Roy.Astr.Soc.，1988，193：393－407.

［3］ Carcione J M，D Kosloff，R Kosloff.Wave propagation simulation in a linear viscoelastic medium［J］.Geophys.J.Roy.Astr.Soc.，1988，95：597－611.

［4］ Carcione J M.Wave propagation in anisotropic linear viscoelastic media［J］.Geophys.J.Int.，1990，101：739－950.

［5］ Carcione J M.Seismic modeling in viscoelastic media［J］.Geophysics，1993，58（1）：110－120.

［6］ Carcione J M.Constitutive model and wave equations for linear，viscoelastic，anisotropic media［J］.Geophysics，1995，60（2）：537－548.

［7］ 杜啟振，楊慧珠.方位各向異性黏彈性介質(zhì)波場有限元模擬［J］.物理學(xué)報，2003，52（8）：2010－201.

［8］ 張智，劉財，邵志剛，等.偽譜法在常Q粘彈介質(zhì)地震波場模擬中的應(yīng)用效果［J］.地球物理學(xué)進展，2005，20（4）：945－949.

［9］ 王德利，雍運動，韓立國，等.三維粘彈介質(zhì)地震波場有限差分并行模擬［J］.西北地震學(xué)報，2007，29（1）：30－34.

［10］ 唐啟軍，韓立國，王恩利，等.基于隨機各向同性背景的粘彈性單斜介質(zhì)二維三分量正演模擬［J］.西北地震學(xué)報，2009，31（1）：35－39.

［11］ Jose M，Carcione，Gerard C，Herman，et al..Review Article：Seismic modeling［J］.Geophysics，2002，67（4）：1304－1325.

［12］ 鄭洪偉，李廷棟，高銳，等.數(shù)值模擬在地球動力學(xué)中的研究進展［J］.地球物理學(xué)進展，2006，21（2）：360－369.

［13］ 劉魯波，陳曉非，王彥賓.切比雪夫偽譜法模擬地震波場［J］.西北地震學(xué)報，2007，29（1）：18－25.

［14］ Zhou B，Greenhalgh S A.Seismic scalar wave equation modeling by a convolutional diffentiator［J］.Bulletin of Seismological Society of America，1992，82（1）：289－303.

［15］ Zhou B，Greenhalgh S A.Numerical seismogram computations for inhomogeneous media using a short，variable length convolutional differentiator［J］.Geophysical Prospecting，1993，41：751－766.

［16］ 張中杰，滕吉文，楊頂輝.聲波與彈性波場數(shù)值模擬中的迭積微分算子法［J］.地震學(xué)報，1996，18（1）：63－69.

［17］ 戴志陽，孫建國，查顯杰.地震波場模擬中的褶積微分算子法［J］.吉林大學(xué)學(xué)報（地球科學(xué)版），2005，35（4）：520－524.

［18］ 戴志陽，孫建國，查顯杰.地震波混合階迭積算法模擬［J］.物探化探計算技術(shù)，2005 27（2）：111－114.

［19］ 謝靖.物探數(shù)據(jù)處理的數(shù)學(xué)方法［M］.北京：地質(zhì)出版社，1981：97－104.

［20］ 程冰潔，李小凡，龍桂華.基于廣義正交多項式褶積微分算子的地震波場數(shù)值模擬方法［J］.地球物理學(xué)報，2008，51（2）：531－537.

［21］ 劉紅艷，李小凡，張美根.非均勻介質(zhì)中地震波傳播的數(shù)值模擬［J］.物探化探計算技術(shù)，2008，30（3）：173－177.

［22］ 劉紅艷.復(fù)雜介質(zhì)中地震波場數(shù)值模擬研究——基于廣義正交多項式迭積微分算子法［D］.北京：中國科學(xué)院地質(zhì)與地球物理研究所，2008.

［23］ 程冰潔，李小凡.2.5維地震波場褶積微分算子法數(shù)值模擬［J］.地球物理學(xué)進展，2008，23（4）：1099－1105.

［24］ 李信富，李小凡.復(fù)雜介質(zhì)地震波傳播的褶積微分算子數(shù)值模擬［J］.地震學(xué)報，2008，30（4）：377－382.

［25］ Gazdag J.Modelling of the acoustic wave equation with transform method［J］.Geophysics，1981，46：854.

［26］ Cerjan C，Kosloff D，et al..A non－reflection boundary condition for discrete acoustic and elastic wave equation［J］.Geophysics，1985，50：705－708.

［27］ Kosloff R.Absorbing boundaries for wave propagation problems［J］.J.Comput.Phys.，1986，63：363－376.

Numerical Simulation of Seismic Wave Propagation in Viscous－elastic Media by the Convolutional Differentiator Method

HE Tong－jiang1，LIU Hong－yan1，LI Xiao－fan2
（1.Earthquake Administration of Tianjin Municipality，Tianjin 300201，China；2.Institute of Geology and Geophysics，Chinese Academy of Science，Beijing 100029，China）

P631.4＋12

A

1000－0844（2010）04－0318－07

2009－06－25

國家自然科學(xué)基金項目（編號40874024）；天津市地震安全基礎(chǔ)工程海上地震觀測試驗系統(tǒng)分項目

賀同江（1974－），男（漢族），碩士，工程師，現(xiàn)在天津市地震局監(jiān)測預(yù)報中心工作.