基于平均導(dǎo)數(shù)方法的聲波方程頻率域高階正演

張衡，劉洪，劉璐，金維浚，史小東

1中國(guó)科學(xué)院地質(zhì)與地球物理研究所 中國(guó)科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室，北京 100029

2中國(guó)科學(xué)院大學(xué)，北京 100049

1 引言

頻率域正演是頻率域全波形反演的基礎(chǔ)（Virieux and Operto，2009；劉璐等，2013b），它最早是用有限元的方法來(lái)實(shí)現(xiàn)的，由Lysmer和Drake（1972）提出用于模擬地震波傳播并由 Marfurt（1984）和Shin（1988）進(jìn)一步發(fā)展.Pratt（1990）將有限差分算法用于頻率域正演中來(lái)進(jìn)行跨孔層析成像和地震波形反演的研究.

頻率域正演的優(yōu)勢(shì)在于多個(gè)頻率可以單獨(dú)模擬，適于多炮同時(shí)模擬，且不存在累積誤差.在頻率域中正演問題演變成對(duì)于每個(gè)給定的頻率求解大型稀疏線性方程組的問題，一般需要較大的內(nèi)存存儲(chǔ)空間，尤其是在三維頻率域大尺度模擬中大型矩陣往往不易求解，從而限制了三維頻率域全波形反演的應(yīng)用（Operto et al.，2007）.

以相速度誤差1%為標(biāo)準(zhǔn)，Pratt和 Worthington（1990）提出的5點(diǎn)聲波頻率域正演算法需要每個(gè)波長(zhǎng)13個(gè)網(wǎng)格點(diǎn)才能正確模擬，計(jì)算精度差.為了解決這個(gè)問題，Jo等（1996）提出了旋轉(zhuǎn)坐標(biāo)系的9點(diǎn)格式來(lái)求解標(biāo)量波動(dòng)方程，引入45°旋轉(zhuǎn)坐標(biāo)系和原水平直角坐標(biāo)系組成的混合網(wǎng)格來(lái)近似拉普拉斯算子，并對(duì)加速度項(xiàng)進(jìn)行加權(quán)平均，采用優(yōu)化方法求取優(yōu)化系數(shù).這種優(yōu)化9點(diǎn)混合網(wǎng)格差分格式能使得每個(gè)波長(zhǎng)所需要的網(wǎng)格點(diǎn)數(shù)在1%的誤差范圍內(nèi)降低到4個(gè)網(wǎng)格點(diǎn)數(shù)，明顯提高了計(jì)算精度.Shin和Sohn（1998）隨后將這種方法擴(kuò)展到4階25點(diǎn)格式，能使得每個(gè)波長(zhǎng)所需要的網(wǎng)格點(diǎn)數(shù)在1%的誤差范圍內(nèi)降低到2.5個(gè)網(wǎng)格點(diǎn)數(shù).曹書紅和陳景波（2012）推導(dǎo)的17點(diǎn)格式在常規(guī)四階9點(diǎn)格式的基礎(chǔ)上引入了45°方向的旋轉(zhuǎn)坐標(biāo)系，可將一個(gè)波長(zhǎng)內(nèi)所需的網(wǎng)格點(diǎn)數(shù)降低到2.56個(gè).劉璐等（2013a）則從減小系數(shù)矩陣帶寬的角度出發(fā)推導(dǎo)得到15點(diǎn)格式，在減小帶寬的同時(shí)保持了相當(dāng)?shù)挠?jì)算精度，從而提高了計(jì)算效率.

但是，Jo等（1996）提出的混合二階9點(diǎn)格式只適用于橫縱向等間隔情況，這在實(shí)際波場(chǎng)模擬中受限很大.針對(duì)旋轉(zhuǎn)坐標(biāo)系的局限性，Chen（2012）利用平均導(dǎo)數(shù)方法（ADM）提出了二階九點(diǎn)格式頻率域正演方法.該方法將聲波頻率域方程中空間導(dǎo)數(shù)項(xiàng)的差分近似表示為正交方向上3個(gè)網(wǎng)格點(diǎn)的加權(quán)平均形式，能適用于不同的方向的空間采樣間隔，因此能夠作為二階精度頻率域有限差分正演的一般格式，增加了方法的靈活性并拓展了應(yīng)用范圍.Chen（2012）的平均導(dǎo)數(shù)二階9點(diǎn)算法能將每個(gè)波長(zhǎng)所需的網(wǎng)格點(diǎn)數(shù)降低到4個(gè)網(wǎng)格點(diǎn)數(shù)，但是二階的精度仍然較差，因此需要發(fā)展高階的算法.為了降低逆矩陣的帶寬和存儲(chǔ)量，在保持精度較高的情況下導(dǎo)數(shù)的差分逼近需要階數(shù)盡可能的低，因此本文選擇四階算法開展研究.

本文首先介紹了四階常規(guī)頻率域正演算法和旋轉(zhuǎn)坐標(biāo)頻率域正演算法，討論了四階情況下旋轉(zhuǎn)坐標(biāo)系的局限性并做了證明.針對(duì)傳統(tǒng)四階格式的局限性，在 Chen（2012）二階平均導(dǎo)數(shù)方法（ADM）的基礎(chǔ)上提出四階頻率域ADM 25點(diǎn)有限差分正演差分格式，并證明這種方法能作為高階頻率域正演的一般格式，四階常規(guī)格式和經(jīng)典的四階25點(diǎn)混合格式都可以作為這種方法的一種特例.隨后通過最小二乘方法求取優(yōu)化系數(shù)并做了頻散分析，表明ADM 25點(diǎn)算法能夠?qū)⒚總€(gè)波長(zhǎng)所需的網(wǎng)格點(diǎn)數(shù)降低到2.78個(gè)，相比二階平均導(dǎo)數(shù)格式和四階常規(guī)格式明顯提高了計(jì)算精度.最后，結(jié)合PML吸收邊界條件進(jìn)行數(shù)值模擬，模型測(cè)試表明ADM 25點(diǎn)算法能對(duì)不同橫縱向間距的模型進(jìn)行高精度模擬.

2 傳統(tǒng)四階差分格式及其局限性

頻率域二維標(biāo)量波動(dòng)方程可寫為

其中，P是位移分量，ω是角頻率，v為速度.

二階聲波頻率域正演精度低，數(shù)值頻散嚴(yán)重，需要較小的網(wǎng)格間距.而四階正演能提高正演的精度并降低數(shù)值頻散.傳統(tǒng)四階差分格式有常規(guī)四階9點(diǎn)差分格式和混合四階25點(diǎn)差分格式.

將式（1）中的二階空間導(dǎo)數(shù)項(xiàng)寫成常規(guī)四階差分的形式即可得到如下的四階常規(guī)9點(diǎn)差分格式，圖1a為這種差分格式的示意圖.

圖1 （a）四階常規(guī)9點(diǎn)差分格式示意圖；（b）四階混合網(wǎng)格25點(diǎn)差分格式示意圖；（c）四階ADM 25點(diǎn)差分格式示意圖Fig.1 （a）Fourth-order conventional nine-point scheme；（b）Fourth-order mixed grid 25-point scheme；（c）Fourth-order ADM 25-point scheme

引入旋轉(zhuǎn)坐標(biāo)系構(gòu)建混合網(wǎng)格的思想是由Jo等（1996）首先提出來(lái)的，Jo引入旋轉(zhuǎn)坐標(biāo)系構(gòu)建了9點(diǎn)的二階混合網(wǎng)格差分格式，Shin和Sohn（1998）將之推廣到四階，構(gòu)建了25點(diǎn)格式.將45°、26.6°、63.4°旋轉(zhuǎn)坐標(biāo)系引入四階9點(diǎn)常規(guī)網(wǎng)格差分格式構(gòu)造，即可得到25點(diǎn)混合網(wǎng)格差分格式，圖1b為這種差分格式的示意圖.

其中

式（3）前兩項(xiàng)即為常規(guī)9點(diǎn)差分離散是拉普拉斯項(xiàng)的形式，后四項(xiàng)即為由不同旋轉(zhuǎn)坐標(biāo)系產(chǎn)生的高階差分離散項(xiàng)，后四項(xiàng)通過泰勒公式展開并化簡(jiǎn)得到：

從式（4）—（7）可以看出，只有當(dāng)Δx=Δz時(shí)上述四項(xiàng)才能成為拉普拉斯項(xiàng)的近似，此時(shí)：

而當(dāng)Δx≠Δz時(shí)混合25點(diǎn)格式并不成立，這就是混合25點(diǎn)格式的局限性，這種局限性正是由其為了提高精度所引入的旋轉(zhuǎn)坐標(biāo)系所導(dǎo)致的（證明見附錄A）.推及開去，Shin和Sohn（1998）的混合25點(diǎn)格式引入了三個(gè)角度的旋轉(zhuǎn)坐標(biāo)系，而曹書紅和陳景波（2012）引入的聲波四階17點(diǎn)格式做了簡(jiǎn)化，只將45°的坐標(biāo)系引入到常規(guī)9點(diǎn)格式，但是由于其仍然是基于旋轉(zhuǎn)坐標(biāo)系構(gòu)建的差分格式，因此也只能適用于Δx=Δz等間隔的情形.

Shin和Sohn（1998）指出常規(guī)9點(diǎn)差分格式要求每個(gè)波長(zhǎng)內(nèi)所需的網(wǎng)格點(diǎn)數(shù)為5個(gè)，精度較低，而混合25點(diǎn)網(wǎng)格則能將每個(gè)波長(zhǎng)內(nèi)所需的網(wǎng)格點(diǎn)數(shù)降低到2.5個(gè)，相比常規(guī)9點(diǎn)差分格式而言提高了數(shù)值模擬的精度.但是從以上分析可知，混合25點(diǎn)格式的局限性在于只適用于Δx=Δz等間隔的情形，而并不適用于Δx≠Δz不等間隔的情形，而常規(guī)9點(diǎn)格式雖然適用于Δx≠Δz不等間隔的情形，但是精度較低.

3 基于平均導(dǎo)數(shù)方法的聲波方程頻率域四階25點(diǎn)格式

由于實(shí)際應(yīng)用中經(jīng)常會(huì)出現(xiàn)橫向和縱向不等間隔的情況，此時(shí)基于旋轉(zhuǎn)坐標(biāo)系的混合網(wǎng)格雖然精度較高，但是不再適用，因此需要提出一種新的既能不受間隔限制從而具有廣泛適用性又精度較高的頻率域正演算法.基于Chen（2012）的平均導(dǎo)數(shù)法構(gòu)建差分格式的思路，本文構(gòu)建得到四階25點(diǎn)聲波ADM有限差分格式（圖1c）.

對(duì)式（1）頻率域二維標(biāo)量波動(dòng)方程中的空間導(dǎo)數(shù)項(xiàng)和加速度項(xiàng)分別引入加權(quán)優(yōu)化系數(shù)，空間導(dǎo)數(shù)項(xiàng)的四階差分近似的波場(chǎng)值表示為正交方向上5個(gè)網(wǎng)格點(diǎn)波場(chǎng)值的加權(quán)平均形式，而加速度項(xiàng)則表示為ADM 25點(diǎn)差分格式中所有25點(diǎn)的加權(quán)平均.據(jù)此構(gòu)建得到的四階25點(diǎn)ADM有限差分格式為

其中

式 （8）中 α1、α2、α3、β1、β2、β3、b1、b2、b3、b4、b5、b6、b7、b8、b9是加權(quán)系數(shù)，且有如下關(guān)系：

當(dāng)α1=β1=1，α2=β2=0且b1=1，b2=b3=b4=b5=b6=b7=b8=0時(shí)ADM差分格式（8）可寫成公式（2）的形式，即常規(guī)9點(diǎn)差分是ADM 25點(diǎn)格式的特例.而當(dāng)Δx=Δz=Δ，α1=β1，α2=β2，α3=β3時(shí)ADM差分格式（8）可寫成公式（3）的形式，即混合25點(diǎn)差分也是ADM 25點(diǎn)格式的特例.

因此ADM方法不僅適用于Δx=Δz等間距的情況，也適用于Δx≠Δz的情況，而很多實(shí)際復(fù)雜模型都往往不是等間距的，因此ADM方法是一種相比常規(guī)網(wǎng)格和混合網(wǎng)格更一般性的方法，能適用于不同的網(wǎng)格間距，具有廣泛的適用性.將ADM從二階推到四階，改善了數(shù)值頻散，提高了計(jì)算精度，從而達(dá)到了聲波頻率域的高精度模擬.

4 系數(shù)優(yōu)化與頻散分析

采用經(jīng)典的頻散分析方法，引入一個(gè)平面波表示P（x，z，ω）=P0e-i（kxx+kzz）來(lái)進(jìn)行頻散分析的研究.

將平面波P（x，z，ω）=P0e-i（kxx+kzz）代入 ADM差分格式即式（8），求得相速度頻散關(guān)系式如下：

其中

代人式（9）即可得到Δx≥Δz情況下的相速度頻散關(guān)系式：

其中

令α1=β1=1，α2=β2=0，且b1=1，b2=b3=b4=b5=b6=b7=b8=0時(shí)即可得到常規(guī)9點(diǎn)差分的相速度頻散關(guān)系式：

本文通過最小二乘方法來(lái)求取ADM 25點(diǎn)格式的優(yōu)化系數(shù).求取系數(shù)時(shí)令1／G的取值范圍為 ［0，0.4］，間隔取為0.0001，傳播角度θ的傳播范圍取為角度間隔取為1°，采用最小二乘的方法求取的不同情況下的優(yōu)化系數(shù)如表1所示（以， ， ，， ，，r=11.21.522.533.125為例）.

表1 當(dāng)Δx≥Δz時(shí)對(duì)于不同r=求得的優(yōu)化系數(shù)Table 1 Optimization coefficients for different r=whenΔx≥Δz

表1 當(dāng)Δx≥Δz時(shí)對(duì)于不同r=求得的優(yōu)化系數(shù)Table 1 Optimization coefficients for different r=whenΔx≥Δz

優(yōu)化系數(shù) r=1 r=1.2 r=1.5 r=2 r=2.5 r=3 r=3.125 α1 0.991997726 0.805233989 0.619957247 0.413925303 0.373903717 0.394354979 0.407408592 α2 0.009593831 0.112399090 0.205383107 0.381893427 0.476554836 0.596182744 0.633290989 α3 -0.005592694-0.013647785-0.013944500-0.090884052-0.163637558-0.293022734-0.336605396 β1 0.991997726 1.060121932 1.118246442 0.876284789 0.786932676 0.697703728 0.675165649 β2 0.009593831-0.023780585-0.054965284 0.087325937 0.134458446 0.181993685 0.194014042 β3 -0.005592694-0.007364847-0.004827885-0.025008643-0.027916654-0.030889070-0.031642021 b1 0.889849126 0.864910920 0.865809648 0.751500449 0.717526005 0.674541008 0.662508511 b2 0.023342617 0.039022204 0.043656325 0.101934600 0.126218557 0.156146763 0.164567318 b3 0.023342617 0.037842560 0.041959903 0.088981227 0.102441283 0.121283998 0.126658830 b4 -0.014507490-0.015400714-0.015821289-0.020612454-0.025805466-0.033875380-0.036452751 b5 -0.014507490-0.015905813-0.018837732-0.010689362-0.007258385-0.004748004-0.004141499 b6 0.022740682 0.011553916 0.004760254-0.017542515-0.029133236-0.043937436-0.048151536 b7 0.000889472-0.000391572-0.001168792-0.000068165-0.000872688-0.002058724-0.002438015 b8 -0.003009849-0.001596517 0.000121029 0.002365050 0.005485087 0.010580676 0.012228197 b9 -0.003009849 0.000401680 0.003963411-0.002820826-0.002843110-0.002736470-0.002684763

圖2分別給出了ADM 25點(diǎn)格式和常規(guī)9點(diǎn)格式在不同r=Δx／Δz情況下根據(jù)表1求取的優(yōu)化系數(shù)計(jì)算得到的頻散曲線.通過頻散曲線分析得出，常規(guī)9點(diǎn)格式當(dāng)G大于5的時(shí)候相速度誤差才能控制在±1%以內(nèi)，而達(dá)到同樣的誤差精度ADM 25點(diǎn)格式只需要G=2.78，ADM 25點(diǎn)格式的頻散誤差小于常規(guī)9點(diǎn)格式的頻散誤差，即ADM 25點(diǎn)格式精度更高.

為驗(yàn)證ADM 25點(diǎn)格式系數(shù)分別在Δx≥Δz和Δz＞Δx兩種情況下的幾何對(duì)稱性，以r=Δz／Δx=2.5為例，將Δx和Δz對(duì)應(yīng)項(xiàng)的系數(shù)互換，從r=Δx／Δz=2.5得到r=Δz／Δx=2.5的優(yōu)化系數(shù)：

α1=0.786932676， α2=0.134458446，

α3=-0.027916654， β1=0.373903717，

β2=0.476554836， β3=-0.163637558，

b1=0.717526005， b2=0.102441283，

b3=0.126218557， b4=-0.007258385，

b5=-0.025805466， b6=-0.029133236，

b7=-0.000872688， b8=-0.002843110，

b9=0.005485087.

圖3為根據(jù)這些優(yōu)化系數(shù)求得的頻散曲線及與常規(guī)9點(diǎn)格式在r=Δz／Δx=2.5情況下的相速度頻散曲線對(duì)比圖，分析得出ADM 25點(diǎn)格式只需要每個(gè)波長(zhǎng)2.78個(gè)網(wǎng)格點(diǎn)就能將數(shù)值頻散誤差限制到1%以內(nèi)，因此可以根據(jù)ADM 25點(diǎn)格式系數(shù)的幾何對(duì)稱性從Δx≥Δz情況直接得到Δz＞Δx情況的優(yōu)化系數(shù).

5 數(shù)值模擬與分析

5.1 構(gòu)建ADM 25點(diǎn)格式PML波動(dòng)方程

吸收邊界條件的效果極大地影響正演模擬求解的效果，因此采用一個(gè)好的吸收邊界條件非常重要.目前應(yīng)用的吸收邊界條件包括單程波旁軸吸收邊界條件和衰減吸收邊界條件兩大類，Pratt（1990）和Chen（2012）在頻率域模擬中采用了Clayton和Enquist（1977）的45°單程波旁軸吸收邊界條件，Shin（1995）提出了添加阻尼層的頻率域海綿吸收邊界條件，但是這些吸收邊界條件在實(shí)際應(yīng)用時(shí)邊界吸收效果較差.為了更好地吸收邊界反射，本文數(shù)值模擬邊界處理時(shí)采用目前應(yīng)用效果最好的PML吸收邊界條件來(lái)消除人工邊界反射.PML是一種衰減型的吸收邊界條件，最早由Berenger（1994）提出.本文采用二維頻率域標(biāo)量PML聲波波動(dòng)方程，并構(gòu)建得到適于ADM 25點(diǎn)格式的PML波動(dòng)方程.

圖2 當(dāng)Δx≥Δz時(shí)對(duì)于不同r=Δx／Δz常規(guī)四階9點(diǎn)格式和四階ADM 25點(diǎn)格式相速度頻散曲線圖：左圖為常規(guī)四階9點(diǎn)格式相速度頻散曲線圖；右圖為四階ADM 25點(diǎn)格式相速度頻散曲線圖Fig.2 Phase velocity curves of the conventional Nine-point scheme and ADM 25-point scheme for different r= Δx／Δz whenΔx≥Δz：Left figure is the phase velocity curves of the conventional Nine-point scheme and The Right figure is the ADM 25-point scheme

二維頻率域標(biāo)量PML聲波波動(dòng)方程可寫為（任浩然等，2009）

式（12）結(jié)合ADM 25點(diǎn)差分公式即式（8）得到ADM 25點(diǎn)格式的PML波動(dòng)方程如下：

式中

從而得到系數(shù)矩陣M的形式：

根據(jù)速度模型不同的縱橫向采樣間隔求得的優(yōu)化系數(shù)，結(jié)合構(gòu)建的ADM 25點(diǎn)格式PML波動(dòng)方程，即可實(shí)現(xiàn)ADM 25點(diǎn)格式頻率域數(shù)值模擬.下面通過一個(gè)簡(jiǎn)單均勻模型和復(fù)雜Marmousi模型來(lái)分析ADM 25點(diǎn)格式在不同縱橫向采樣間隔情況下的正演精度和正演效率，常規(guī)四階9點(diǎn)算法和ADM 9點(diǎn)算法因?yàn)橐材苓m用于不同的網(wǎng)格間隔，作為對(duì)比方法進(jìn)行比較.

5.2 均勻模型測(cè)試

均勻模型的單道記錄能更好地測(cè)試不同格式數(shù)值頻散的效果及計(jì)算精度.取模型為200×200的網(wǎng)格數(shù)，取波速為2000m·s-1的一個(gè)均勻模型（圖4）.數(shù)值模擬時(shí)，雷克震源子波主頻取20Hz，震源位置設(shè)置在模型的中心，單檢波點(diǎn)取在震源位置左側(cè)50個(gè)網(wǎng)格點(diǎn)處.

均勻模型測(cè)試采用常規(guī)9點(diǎn)格式和ADM 25點(diǎn)格式進(jìn)行頻率域正演模擬并與解析方法進(jìn)行對(duì)比，解析解采用下面公式（Alford et al.，1974）：

其中，F(xiàn)（ω）為頻率域震源，H0（2）為零階二類Hankel函數(shù)，r為觀測(cè)點(diǎn)到震源點(diǎn)的距離.

圖4 均勻模型示意圖Fig.4 Sketch map of homogeneous velocity model

第一種情況設(shè)橫向網(wǎng)格間距為9.6m，縱向網(wǎng)格間距為8m，此時(shí)橫向縱向網(wǎng)格間距比為1.2，傳播時(shí)間0.6s，時(shí)間采樣間隔為2ms.分別采用常規(guī)9點(diǎn)格式和ADM 25點(diǎn)格式進(jìn)行頻率域正演模擬，橫向縱向網(wǎng)格間距比為1.2情況下的ADM 25點(diǎn)差分格式優(yōu)化系數(shù)如表1所示.取常規(guī)9點(diǎn)格式和ADM 25點(diǎn)格式兩種格式的單道記錄并與解析方法進(jìn)行對(duì)比（圖5a），如圖在這種情況下兩種格式都能精確模擬，都沒有出現(xiàn)數(shù)值頻散誤差，且與解析解對(duì)應(yīng)良好.第二種情況如果將橫向網(wǎng)格間距和縱向網(wǎng)格間距分別增大到13.2m和11m，橫向縱向網(wǎng)格間距比保持1.2不變，除傳播時(shí)間增加到0.8s外，其他參數(shù)不變.分別取常規(guī)9點(diǎn)格式和ADM 25點(diǎn)格式兩種格式的單道記錄并與解析方法進(jìn)行對(duì)比（圖5b），此時(shí)ADM 25點(diǎn)的單道記錄幾乎沒有頻散與解析解對(duì)應(yīng)良好，而常規(guī)9點(diǎn)的單道記錄出現(xiàn)了很嚴(yán)重的頻散，與解析解相比出現(xiàn)較大誤差.這是由于兩種方法數(shù)值模擬精度差異較大而導(dǎo)致常規(guī)9點(diǎn)格式在大網(wǎng)格間距下數(shù)值頻散誤差較大，而ADM 25點(diǎn)仍能保持高精度.具體而言由于常規(guī)9點(diǎn)格式達(dá)到無(wú)頻散模擬所需的每個(gè)波長(zhǎng)網(wǎng)格點(diǎn)數(shù)為5個(gè)，則最大網(wǎng)格間距大概為2000／40／5=10m，而ADM 25點(diǎn)格式達(dá)到無(wú)頻散模擬所需的每個(gè)波長(zhǎng)網(wǎng)格點(diǎn)數(shù)僅為2.78個(gè)，則最大網(wǎng)格間距可達(dá)到2000／40／2.78=18m，因此對(duì)于第一種情況，橫縱向最大網(wǎng)格間距為9.6m，此時(shí)兩種格式都能精確模擬，而對(duì)于第二種情況，橫縱向最大網(wǎng)格間距為13.2m，已經(jīng)超過了常規(guī)9點(diǎn)格式精確模擬所允許的10m間距，因此產(chǎn)生了比較嚴(yán)重的數(shù)值頻散，而ADM 25點(diǎn)格式仍能精確模擬.

5.3 Marmousi模型測(cè)試

下面采用較復(fù)雜的Marmousi模型來(lái)進(jìn)行模擬.標(biāo)準(zhǔn)Marmousi模型橫向網(wǎng)格間距為12.5m，而縱向網(wǎng)格間距為4m，橫縱向網(wǎng)格間距比達(dá)到了3.125，而且 Marmousi模型是一個(gè)非均勻復(fù)雜模型，能進(jìn)一步測(cè)試ADM 25點(diǎn)算法的精度和效率.本文截取了Marmousi速度模型網(wǎng)格大小為301×301的一個(gè)區(qū)域（圖6a）.數(shù)值模擬時(shí)，雷克震源子波主頻取15Hz，震源位置設(shè)置在縱向第11層的橫向40個(gè)網(wǎng)格點(diǎn)處，坐標(biāo)為（500m，40m），檢波點(diǎn)排列置于地表，每個(gè)網(wǎng)格點(diǎn)設(shè)置一個(gè)檢波器，全排列共301個(gè)檢波器，PML層數(shù)設(shè)為40層，傳播時(shí)間為2 s，時(shí)間采樣間隔為0.004s.

分別采用常規(guī)9點(diǎn)格式、ADM 9點(diǎn)格式和ADM 25點(diǎn)格式進(jìn)行頻率域正演模擬，橫向縱向網(wǎng)格間距比為3.125情況下的ADM 25點(diǎn)差分格式優(yōu)化系數(shù)如表1所示.從15Hz的單頻波場(chǎng)切片（圖6b）可以看出，在左側(cè)速度漸變帶，頻率域波場(chǎng)變化不大，而在中間和右側(cè)速度變化復(fù)雜，入射波和反射波的相互干涉導(dǎo)致了反射界面上頻率域波場(chǎng)的劇烈擾動(dòng).將頻率域波場(chǎng)通過反傅里葉變換轉(zhuǎn)換到時(shí)間域得到時(shí)間域的波場(chǎng)快照和地表檢波器采集的地震記錄.從常規(guī)9點(diǎn)格式的地震記錄（圖6c）、ADM 9點(diǎn)格式的地震記錄（圖6d）以及ADM 25點(diǎn)格式的地震記錄（圖6e）對(duì)比可以清晰看出，ADM 25點(diǎn)能精確模擬，而常規(guī)9點(diǎn)和ADM 9點(diǎn)都出現(xiàn)了較強(qiáng)的數(shù)值頻散（如圖6（c，d，e）黑色箭頭標(biāo)示），因此常規(guī)9點(diǎn)格式和ADM 9點(diǎn)格式在橫縱向網(wǎng)格間距不一致的情況下波場(chǎng)模擬精度差，相較本文方法更容易出現(xiàn)數(shù)值頻散誤差，本文方法體現(xiàn)出在克服數(shù)值頻散誤差方面的優(yōu)越性.另外從地震記錄上來(lái)看PML吸收邊界條件的引入很好地消除了人工邊界反射（圖6（c，d，e））.通過與圖6f中時(shí)間域12階高階有限差分方法（Yan and Liu，2013）的結(jié)果進(jìn)行對(duì)比，驗(yàn)證了本文方法與時(shí)間域高階方法精度相當(dāng).

對(duì)Marmousi模型相同計(jì)算區(qū)域同等計(jì)算精度下做了一個(gè)計(jì)算效率測(cè)試，計(jì)算平臺(tái)為聯(lián)想Think Centre M8300t臺(tái)式機(jī)電腦（酷睿i7八核，3.4GHz）.單炮測(cè)試結(jié)果為常規(guī)9點(diǎn)格式頻率域單炮模擬耗時(shí)8427s，ADM 9點(diǎn)格式頻率域單炮模擬耗時(shí)3650s，而ADM 25點(diǎn)格式頻率域單炮模擬耗時(shí)2772s.因此在相同的精度要求下，ADM 25點(diǎn)的計(jì)算效率要高于ADM 9點(diǎn)，并明顯高于常規(guī)9點(diǎn)，這是由于ADM 25點(diǎn)格式精度相比常規(guī)9點(diǎn)格式和ADM 9點(diǎn)方法提升明顯，因此可以通過采用更大的網(wǎng)格間距來(lái)提高計(jì)算效率.需要特別說(shuō)明的是，更高階的頻率域正演方法帶來(lái)的精度提升空間有限，并不能顯著提升計(jì)算精度，反而因?yàn)橄禂?shù)矩陣帶寬增加會(huì)明顯增加計(jì)算成本，降低計(jì)算效率.

就內(nèi)存存儲(chǔ)空間而言，如果常規(guī)9點(diǎn)格式需要nx×nz（nx和nz分別為橫向和縱向網(wǎng)格點(diǎn)數(shù)）個(gè)網(wǎng)格點(diǎn)來(lái)存儲(chǔ)系數(shù)矩陣，則ADM 9點(diǎn)格式需要個(gè)網(wǎng)格點(diǎn)，本文的ADM 25點(diǎn)格式只需要個(gè)網(wǎng)格點(diǎn).因此ADM 25點(diǎn)格式耗費(fèi)的內(nèi)存存儲(chǔ)空間要少于常規(guī)9點(diǎn)和ADM 9點(diǎn)，從而減少了存儲(chǔ)要求.

圖3 當(dāng)r=Δz／Δx=2.5時(shí)常規(guī)四階9點(diǎn)格式和四階ADM 25點(diǎn)格式相速度頻散曲線圖Fig.3 Phase velocity curves of the conventional nine-point scheme and ADM 25-point scheme for r=Δz／Δx=2.5

圖5 均勻模型常規(guī)9點(diǎn)格式、ADM 25點(diǎn)格式與解析方法單道記錄對(duì)比Fig.5 Comparison of single traces obtained by conventional nine-point scheme、ADM 25-point scheme and analytical method on homogeneous model

圖6 （a）截取的Marmousi速度模型；（b）ADM 25點(diǎn)格式計(jì)算的15Hz單頻波場(chǎng)；（c）常規(guī)9點(diǎn)格式計(jì)算得到的時(shí)間域地震記錄；（d）ADM 9點(diǎn)格式計(jì)算得到的時(shí)間域地震記錄；（e）ADM 25點(diǎn)格式計(jì)算得到的時(shí)間域地震記錄；（f）時(shí)間域12階高階有限差分方法計(jì)算得到的時(shí)間域地震記錄Fig.6 （a）Truncated Marmousi model；（b）40Hz monochromatic wavefield computed by ADM 25-point scheme；（c）Time-domain seismograms computed with the conventional nine-point scheme；（d）Time-domain seismograms computed with the ADM 9-point scheme；（e）Time-domain seismograms computed with the ADM 25-point scheme；（f）Time-domain seismograms computed with the time domain 12-order high order finite difference scheme

6 結(jié)論

頻率域常規(guī)算法能適用于不同的采樣間隔，但是精度較低，而基于旋轉(zhuǎn)坐標(biāo)系的混合網(wǎng)格算法雖然提高了精度，但是它的局限性在于只適用于相同的方向采樣間隔，從而限制了其實(shí)際應(yīng)用的靈活性.而基于平均導(dǎo)數(shù)方法的頻率域正演方法綜合了常規(guī)算法和混合網(wǎng)格算法的優(yōu)點(diǎn)，能在保持混合網(wǎng)格精度的同時(shí)適用于不同的方向采樣間隔.Chen（2012）的二階平均導(dǎo)數(shù)9點(diǎn)算法能將每個(gè)波長(zhǎng)所需的網(wǎng)格點(diǎn)數(shù)降低到4個(gè)網(wǎng)格點(diǎn)數(shù)，但是二階的精度仍然較差.本文將二階的平均導(dǎo)數(shù)頻率域正演算法推廣到四階，發(fā)展了四階平均導(dǎo)數(shù)頻率域正演的算法，構(gòu)造得到一個(gè)優(yōu)化的四階精度25點(diǎn)格式，通過系數(shù)優(yōu)化將每個(gè)波長(zhǎng)所需的網(wǎng)格點(diǎn)數(shù)降低到2.78個(gè)，相比二階平均導(dǎo)數(shù)格式和四階常規(guī)格式明顯提高了計(jì)算精度.數(shù)值模擬結(jié)果表明本文的基于平均導(dǎo)數(shù)方法的頻率域高階正演方法不僅具有很好的精度和效率，而且具有很好的適用性和靈活性.

附錄A：旋轉(zhuǎn)坐標(biāo)系二階空間導(dǎo)數(shù)差分離散

下面證明旋轉(zhuǎn)坐標(biāo)系的二階空間導(dǎo)數(shù)差分離散只適用于等間隔情形，以45°旋轉(zhuǎn)坐標(biāo)系為例.

45°旋轉(zhuǎn)坐標(biāo)系二階空間導(dǎo)數(shù)差分離散公式為

對(duì)式（A1）右端項(xiàng)的 Pm+1，n+1、Pm-1，n+1、Pm+1，n-1和Pm-1，n-1分別進(jìn)行泰勒展開：

則有

對(duì)式（A2）右端項(xiàng)的Pm，n+1和Pm，n-1再分別進(jìn)行泰勒展開：

則有

式（A3）式代人式（A2）式得：

即可得到

結(jié)合式（A1）得到45°旋轉(zhuǎn)坐標(biāo)系的二階空間導(dǎo)數(shù)差分離散公式為

其他角度旋轉(zhuǎn)坐標(biāo)系二階空間導(dǎo)數(shù)差分離散公式推導(dǎo)依此類推.顯然上述公式只在Δx=Δz的情況下才成立，因此證明旋轉(zhuǎn)坐標(biāo)系的二階空間導(dǎo)數(shù)差分離散只適用于等間隔情形.

Alford R M，Kelly K R，Boore D M.1974.Accuracy of finitedifference modeling of the acoustic wave equation.Geophysics，39（6）：834-842，doi：10.1190／1.1440470.

Berenger J P.1994.A perfectly matched layer for the absorption of electromagnetic waves.Journal of Computational Physics，114（2）：185-200，doi：10.1006／jcph.1994.1159.

Cao S H，Chen J B.2012.A 17-point scheme and its numerical implementation for high-accuracy modeling of frequency-domain acoustic equation.Chinese J.Geophys.（in Chinese），55（10）：3440-3449，doi：10.6038／j.issn.0001-5733.2012.10.027.

Chen J B.2012.An average-derivative optimal scheme for frequency-domain scalar wave equation.Geophysics，77（6）：T201-210，doi：10.1190／GEO2011-0389.1.

Clayton R，Enquist B.1977.Absorbing boundary conditions for acoustic and elastic wave equations.Bulletin of the Seismological Society of America，67（6）：1529-1540.

Jo C H，Shin C S，Suh J H.1996.An optimal 9-point，finitedifference，frequency-space，2-D scalar wave extrapolator.Geophysics，61（2）：529-537，doi：10.1190／1.1443979.

Liu L，Liu H，Liu H W.2013a.Optimal 15-point finite difference forward modeling in frequency-space domain.Chinese J.Geophys.（in Chinese），56（2）：644-652，doi：10.6038／cjg20130228.

Liu L，Liu H，Zhang H，et al.2013b.Full waveform inversion based on modified quasi-Newton equation.Chinese J.Geophys.（in Chinese），56（7）：2447-2451，doi：10.6038／cjg20130730.

Lysmer J，Drake L A.1972.A finite-element method for seismology.∥Bolt B A.Methods in computational physics，vol.11：Seismology：Surface waves and earth oscillations.New York：Academic Press，11：181-216.

Marfurt K J.1984.Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations.Geophysics，49（5）：533-549，doi：10.1190／1.1441689.

Operto S，Virieux J，Amestoy P，et al.2007.3Dfinite-difference frequency-domain modeling of visco-acoustic wave propagation using a massively parallel direct solver：A feasibility study.Geophysics，72（5）：SM195-SM211，doi：10.1190／1.2759835.

Pratt R G.1990.Frequency-domain elastic wave modeling by finite differences：A tool for cross-hole seismic imaging.Geophysics，55（5）：626-632，doi：10.1190／1.1442874.

Pratt R G，Worthington M H.1990.Inverse theory applied to multi-source cross-hole tomography，part 1：Acoustic wave equation method.Geophysical Prospecting，38（3）：287-310，doi：10.1111／j.1365-2478.1990.tb01846.x.

Ren H R，Wang H Z，Gong T.2009.Seismic modeling of scalar seismic wave propagation with finite-difference scheme in frequency-space domain.Geophysical Prospecting for Petroleum（in Chinese），48（1）：20-26.

Shin C.1988.Nonlinear elastic wave inversion by blocky parameterization［Ph.D.thesis］.Tulsa，OK：University of Tulsa.

Shin C.1995.Sponge boundary condition for frequency-domain modeling.Geophysics，60（6）：1870-1874，doi：10.1190／1.1443918.

Shin C S，Sohn H J.1998.A frequency-space 2-D scalar wave extrapolator using extended 25-point finite-difference operator.Geophysics，63（1）：289-296，doi：10.1190／1.1444323.

Virieux J，Operto S.2009.An overview of full-waveform inversion in exploration geophysics.Geophysics，74（6）：WCC1-WCC26，doi：10.1190／1.3238367.

Wu G C，Luo C M，Liang K.2007.Frequency-space domain finite difference numerical simulation of elastic wave in TTI media.Journal of Jilin University （Earth Science Edition）.（in Chinese），37（5）：1023-1033.

Yan H Y，Liu Y.2013.Visco-acoustic prestack reverse-time migration based on the time-space domain adaptive high-order finite-difference method.Geophysical Prospecting，61（5）：941-954，doi：10.1111／1365-2478.12046.

附中文參考文獻(xiàn)

曹書紅，陳景波.2012.聲波方程頻率域高精度正演的17點(diǎn)格式及數(shù)值實(shí)現(xiàn).地球物理學(xué)報(bào)，55（10）：3440-3449，doi：10.6038／j.issn.0001-5733.2012.10.027.

劉璐，劉洪，劉紅偉.2013a.優(yōu)化15點(diǎn)頻率-空間域有限差分正演模擬.地球物理學(xué)報(bào)，56（2）：644-652，doi：10.6038／cjg20130228.

劉璐，劉洪，張衡等.2013b.基于修正擬牛頓公式的全波形反演.地 球 物 理 學(xué) 報(bào)，56（7）：2447-2451，doi：10.6038／cjg20130730.

任浩然，王華忠，龔婷.2009.標(biāo)量地震波頻率-空間域有限差分法數(shù)值模擬.石油物探，48（1）：20-26.

吳國(guó)忱，羅彩明，梁楷.2007.TTI介質(zhì)彈性波頻率-空間域有限差分?jǐn)?shù)值模擬.吉林大學(xué)學(xué)報(bào)（地球科學(xué)版），37（5）：1023-1033.