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

1中國科學(xué)院地質(zhì)與地球物理研究所 中國科學(xué)院油氣資源研究重點實驗室，北京 100029

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

1 引言

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

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

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

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

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

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

頻率域二維標量波動方程可寫為

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

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

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

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

其中

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

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

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

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

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

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

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

其中

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

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

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

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

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

將平面波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時即可得到常規(guī)9點差分的相速度頻散關(guān)系式：

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

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

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

為驗證ADM 25點格式系數(shù)分別在Δx≥Δz和Δz＞Δx兩種情況下的幾何對稱性，以r=Δz／Δx=2.5為例，將Δx和Δz對應(yī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點格式在r=Δz／Δx=2.5情況下的相速度頻散曲線對比圖，分析得出ADM 25點格式只需要每個波長2.78個網(wǎng)格點就能將數(shù)值頻散誤差限制到1%以內(nèi)，因此可以根據(jù)ADM 25點格式系數(shù)的幾何對稱性從Δx≥Δz情況直接得到Δz＞Δx情況的優(yōu)化系數(shù).

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

5.1 構(gòu)建ADM 25點格式PML波動方程

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

圖2 當Δx≥Δz時對于不同r=Δx／Δz常規(guī)四階9點格式和四階ADM 25點格式相速度頻散曲線圖：左圖為常規(guī)四階9點格式相速度頻散曲線圖；右圖為四階ADM 25點格式相速度頻散曲線圖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

二維頻率域標量PML聲波波動方程可寫為（任浩然等，2009）

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

式中

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

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

5.2 均勻模型測試

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

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

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

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

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

5.3 Marmousi模型測試

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

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

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

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

圖3 當r=Δz／Δx=2.5時常規(guī)四階9點格式和四階ADM 25點格式相速度頻散曲線圖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點格式、ADM 25點格式與解析方法單道記錄對比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點格式計算的15Hz單頻波場；（c）常規(guī)9點格式計算得到的時間域地震記錄；（d）ADM 9點格式計算得到的時間域地震記錄；（e）ADM 25點格式計算得到的時間域地震記錄；（f）時間域12階高階有限差分方法計算得到的時間域地震記錄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)坐標系的混合網(wǎng)格算法雖然提高了精度，但是它的局限性在于只適用于相同的方向采樣間隔，從而限制了其實際應(yīng)用的靈活性.而基于平均導(dǎo)數(shù)方法的頻率域正演方法綜合了常規(guī)算法和混合網(wǎng)格算法的優(yōu)點，能在保持混合網(wǎng)格精度的同時適用于不同的方向采樣間隔.Chen（2012）的二階平均導(dǎo)數(shù)9點算法能將每個波長所需的網(wǎng)格點數(shù)降低到4個網(wǎng)格點數(shù)，但是二階的精度仍然較差.本文將二階的平均導(dǎo)數(shù)頻率域正演算法推廣到四階，發(fā)展了四階平均導(dǎo)數(shù)頻率域正演的算法，構(gòu)造得到一個優(yōu)化的四階精度25點格式，通過系數(shù)優(yōu)化將每個波長所需的網(wǎng)格點數(shù)降低到2.78個，相比二階平均導(dǎo)數(shù)格式和四階常規(guī)格式明顯提高了計算精度.數(shù)值模擬結(jié)果表明本文的基于平均導(dǎo)數(shù)方法的頻率域高階正演方法不僅具有很好的精度和效率，而且具有很好的適用性和靈活性.

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

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

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

對式（A1）右端項的 Pm+1，n+1、Pm-1，n+1、Pm+1，n-1和Pm-1，n-1分別進行泰勒展開：

則有

對式（A2）右端項的Pm，n+1和Pm，n-1再分別進行泰勒展開：

則有

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

即可得到

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

其他角度旋轉(zhuǎn)坐標系二階空間導(dǎo)數(shù)差分離散公式推導(dǎo)依此類推.顯然上述公式只在Δx=Δz的情況下才成立，因此證明旋轉(zhuǎn)坐標系的二階空間導(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.

附中文參考文獻

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

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

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

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

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