王之洋，劉洪，唐祥德，王洋

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

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

1 引言

隨著我國對產(chǎn)油區(qū)的勘探深度不斷加深，對更高精度的成像和反演的需求也越來越迫切，而作為成像和反演方法基礎(chǔ)的地震波場數(shù)值模擬技術(shù)也變得尤為重要.常用的數(shù)值模擬方法分為三大類：幾何射線法、積分方程法和波動方程法，波動方程法的模擬結(jié)果因為包含了波場的運動學(xué)與動力學(xué)特征，是最常用的（程冰潔等，2008）.而使用較多的波動方程數(shù)值模擬方法則是有限差分方法（Alterman and Karal，1968；Kelly et al.，1976；Dablain，1986；Igel et al.，1995）和偽譜法（Gazdag，1981；Kosloffand Baysal，1982；Carcione et al.，1992）.有限差分數(shù)值模擬技術(shù)起源于20世紀60年代末，Alterman和Karal（1968）首次使用顯式有限差分技術(shù)來模擬層狀介質(zhì)中二階彈性波波場.Kelly等（1976）提出和發(fā)展了適合非均勻介質(zhì)的二階彈性波方程有限差分數(shù)值模擬方法.Dablain（1986）實現(xiàn)了聲波方程高階有限差分的數(shù)值模擬.劉洋等（1998）從Taylor級數(shù)展開出發(fā)，推導(dǎo)出了任意偶數(shù)階導(dǎo)數(shù)的任意偶數(shù)精度的有限差分算子.董良國等 （2000）實現(xiàn)了彈性波方程有限差分在時間和空間上任意高階精度的差分解法.裴正林和牟永光 （2003）推導(dǎo)了一階空間導(dǎo)數(shù)交錯網(wǎng)格的任意偶數(shù)階精度差分系數(shù).Liu和Sen（2009a，2010）提出了一種基于時空域頻散關(guān)系的有限差分法，該差分法能有效地改善波動方程數(shù)值解精度.而后，Yan和 Liu（2013a，2013b）將該時空域有限差分法推廣發(fā)展到黏滯聲波和各向異性聲波方程的數(shù)值模擬與逆時偏移中.

數(shù)值頻散問題是有限差分數(shù)值模擬不可避免的一個問題，會直接影響模擬的精度.數(shù)值頻散的產(chǎn)生是由于使用差分算子來逼近微分算子，截斷之后導(dǎo)致了誤差；為了降低數(shù)值頻散，可以采用更細的網(wǎng)格或者降低子波主頻.但是更細的網(wǎng)格意味著海量的計算量，降低了計算效率，且對存儲和計算性能都提出了挑戰(zhàn)；降低子波主頻則會損失高頻成分，影響分辨率.前人在研究如何降低數(shù)值頻散方面做了許多工作，其一是利用通量校正技術(shù)（FCT）來壓制數(shù)值頻散，Boris和Book（1973）提出利用FCT技術(shù)來壓制數(shù)值頻散.楊頂輝和騰吉文（1997）將FCT技術(shù)應(yīng)用到各向異性介質(zhì)中的三分量地震數(shù)值模擬.其二是通過最優(yōu)化方法來壓制數(shù)值頻散，Holberg（1987）和 Robertsson等（1994）利用最優(yōu)化方法，最小化頻散誤差，計算優(yōu)化的有限差分系數(shù).Zhang和Yao（2013）通過優(yōu)選設(shè)定誤差的容許范圍，利用模擬退火法得到優(yōu)化算子，其擁有更高的頻譜覆蓋范圍和較小的誤差限.Liu（2013）基于最小二乘方法，得到一種全局優(yōu)化的有限差分算子.Yang等（2014）基于數(shù)值頻散關(guān)系和最小二乘理論推導(dǎo)了適合求解空間一階導(dǎo)數(shù)的交錯網(wǎng)格優(yōu)化差分系數(shù)，并實現(xiàn)了彈性波數(shù)值模擬.其三是優(yōu)選窗函數(shù)壓制數(shù)值頻散，F(xiàn)ornberg（1987）證明了有限差分方法隨著階數(shù)的增加，逼近偽譜法.偽譜法因為采用了所有的點，解決了數(shù)值頻散的問題，換句話說，在空間域，可以采用不同的截斷窗函數(shù)去截斷偽譜法的空間褶積序列推導(dǎo)出有限差分算子.Zhou和Greenhalgh（1992）使用了廣義加權(quán)的Hanning窗截斷得到了有限差分算子.Igel等（1995）使用了高斯窗截斷得到了有限差分算子.Chu和Stoffa（2012）使用了二項式窗統(tǒng)一了有限差分方法和偽譜法，使用二項式窗截斷不僅可以推導(dǎo)出常規(guī)的有限差分算子系數(shù)，而且稍加改進，就可以設(shè)計出一種優(yōu)化的有限差分算子，Chu和Stoffa（2012）認為，這種改進的二項式窗截斷方案與Liu和Sen（2009b）提出的一種截斷的有限差分方法本質(zhì)上是相同的.

在空間域，將偽譜法的空間褶積序列截斷就得到了有限差分法，如果從信號處理層面來講，截斷相當(dāng)于加了一個窗函數(shù).截斷會產(chǎn)生邊緣效應(yīng)，造成頻譜泄露.為了壓制數(shù)值頻散，選擇合適形狀的有限長的窗函數(shù)使數(shù)據(jù)逐漸截斷；不同的窗函數(shù)會產(chǎn)生不同的結(jié)果，具體來講，窗函數(shù)相當(dāng)于有限長的濾波器，不同的窗函數(shù)，其幅值響應(yīng)擁有不同的主瓣和旁瓣形狀，主瓣的形狀控制著過渡帶的范圍，也就是頻譜覆蓋范圍；旁瓣的形狀決定了有限差分算子逼近微分算子的偏差程度.如果有種窗函數(shù)的幅值響應(yīng)主瓣窄，旁瓣衰減大，由該種窗函數(shù)截斷的有限差分算子的精度誤差譜覆蓋范圍大，誤差穩(wěn)定性高.譜覆蓋范圍大，意味著采用低階窗函數(shù)截斷得到的有限差分算子就可以達到高階常規(guī)有限差分算子的精度；精度誤差穩(wěn)定性高意味著逼近精度不會出現(xiàn)大幅波動；滿足這兩者，表明有限差分算子逼近微分算子的程度就越好，數(shù)值頻散越小.本文就是從窗函數(shù)的性質(zhì)出發(fā)，設(shè)計一種窗函數(shù)，截斷偽譜法的空間褶積序列得到優(yōu)化的有限差分算子，抑制數(shù)值頻散.

Zhou和Greenhalgh（1992）提出的廣義加權(quán)的Hanning窗，其本質(zhì)是三角函數(shù)類窗，三角函數(shù)類窗函數(shù)的幅值響應(yīng)有相對較窄的主瓣，但是旁瓣衰減程度不夠高，使用該三角函數(shù)類窗函數(shù)截斷逼近的有限差分算子，其截斷逼近的精度誤差波動相對較大，即使有較大的譜覆蓋范圍，但是精度誤差的波動對逼近的效果的影響還是很大（Zhang and Yao，2013）.Chu和Stoffa（2012）的改進的二項式窗是一個可調(diào)節(jié)的窗函數(shù)，但是其主瓣還是過寬，導(dǎo)致過渡帶過長，使用該改進的二項式窗函數(shù)截斷逼近的有限差分算子，其精度誤差譜覆蓋范圍較小，且對于低階有限差分算子，其精度誤差存在較大波動，特別是12階以下的差分算子.本文分析了窗函數(shù)幅值響應(yīng)的主旁瓣屬性對有限差分算子逼近微分算子的精度的影響，研究了自褶積組合的效應(yīng)，設(shè)計了一種基于Chebyshev自褶積組合的窗函數(shù)，該窗函數(shù)的幅值響應(yīng)在維持較窄主瓣的前提下，可以獲得更好的旁瓣衰減，從而進一步提高了窗函數(shù)截斷逼近的有限差分算子的精度誤差穩(wěn)定性，并且可以通過只改變?nèi)齻€參數(shù)，更直觀和可視化地調(diào)節(jié)主瓣和旁瓣的形狀，進而控制有限差分算子逼近微分算子的精度，保持了較大的靈活性.

2 窗函數(shù)截斷的逼近誤差分析

2.1 常規(guī)有限差分算子

我們使用sinc函數(shù)插值理論推導(dǎo)出有限差分算子，根據(jù)離散信號的采樣理論（Diniz et al.，2012），一個帶限的連續(xù)信號f（x）可以被以一個均勻采樣的信號fn通過sinc函數(shù)插值重建：

其中，Δx為采樣間隔，為截止波數(shù).

如果對公式（1）左右兩邊分別求一階導(dǎo)數(shù)和二階導(dǎo)數(shù)，并取x＝0處的導(dǎo)數(shù)值，可以得到：

Chu和Stoffa（2012）認為存在一個長度為N＋1點的窗函數(shù)，N為偶數(shù)，去截斷公式（2）和公式（3），得到常規(guī)有限差分算子.

2.2 窗函數(shù)對逼近效果的影響

前人做了很多的工作去選擇不同的截斷窗函數(shù)，Zhou和Greenhalgh（1992）提出了廣義加權(quán)的Hanning窗，Igel等（1995）使用了高斯窗，Chu和Stoffa（2012）提出了改進的二項式窗.Zhou和Greenhalgh（1992）提出的廣義加權(quán)的 Hanning窗本質(zhì)上是一類三角函數(shù)窗.本文比較四種窗函數(shù)，分別是Hanning窗，改進的二項式窗（Chu and Stoffa，2012），Kaiser和Chebyshev窗，分析窗函數(shù)幅值響應(yīng)的主旁瓣屬性對有限差分算子逼近微分算子的精度的影響.

圖1是窗函數(shù)長度為N＋1＝25時的系數(shù)分布，圖2展示了對應(yīng)窗函數(shù)幅值響應(yīng).從圖中可以看出，圖1中窗函數(shù)系數(shù)分布越細長，在圖2中，其幅值響應(yīng)，主瓣越寬.Diniz等（2012）指出較寬的主瓣意味著加窗處理后會出現(xiàn)較寬的過渡帶，隨著過渡帶的增寬，頻譜覆蓋范圍會變小，對高波數(shù)部分的逼近精度就越低，出現(xiàn)較為嚴重的頻散效應(yīng)；相比于其他窗函數(shù)，改進的二項式窗擁有最寬的主瓣，這表明該改進的二項式窗頻譜范圍沒有其余三個窗大，對高波數(shù)部分逼近精度不夠；另外改進的二項式窗擁有衰減最大的旁瓣，旁瓣衰減越大，證明其逼近的穩(wěn)定性越高，逼近精度不會出現(xiàn)較大的波動.Hanning窗雖然擁有較小的主瓣，旁瓣衰減也是最快的，但是其前面幾個頻率點的旁瓣衰減幅度較小，逼近精度會出現(xiàn)較大的波動.相對而言，Kasier窗和Chebyshev窗是比較自由的窗函數(shù)，通過調(diào)節(jié)β和r，可以擁有不同的幅值響應(yīng).這里選擇r＝60，Chebyshev窗的旁瓣衰減固定在－60dB，擁有較窄的主瓣.而Kasier窗和Hanning類似，選擇Kaiser窗的參數(shù)β＝5，由圖2可知，在前面幾個頻率點旁瓣衰減幅度較小，低波數(shù)時會出現(xiàn)逼近精度的波動.

圖1 四種窗函數(shù)對比，N＝24Fig.1 Comparison of window functions，N＝24

為了進一步確認可調(diào)參數(shù)的Kasier窗和Chebyshev窗的幅值響應(yīng)，圖3分別示意β＝2，5，9的Kasier窗的幅值響應(yīng)和r＝35，60，95的Chebyshev窗的幅值響應(yīng).窗函數(shù)長度N＋1＝25.

圖3表明，當(dāng)固定窗函數(shù)長度時，對于Kaiser窗，隨著β增大，其幅值響應(yīng)的變化表現(xiàn)為：主瓣寬度增加，旁瓣衰減增大.而對于Chebyshev窗，隨著r的增大，其幅值響應(yīng)的主旁瓣也有相同的變化；對于這兩種窗函數(shù)來說，Chebyshev窗擁有衰減更大的旁瓣，逼近精度相對來說更加穩(wěn)定.

本文以二階導(dǎo)數(shù)為例，說明不同窗函數(shù)截斷逼近的精度誤差.

因為n＝0是公式（2）和公式（3）的一個奇異點，根據(jù)Lee和Seo（2002），公式（2）和公式（3）可以表示為：

圖2 四種窗函數(shù)的幅值響應(yīng)，頻率點M＝512Fig.2 Amplitude response of four window functions，frequency points M＝512

圖3 改變β和r，兩種窗函數(shù)幅值響應(yīng)，頻率點M＝512（a）Kaiser窗；（b）Chebyshev窗.Fig.3 Amplitude response of two window functions，differentβand r，frequency points M＝512（a）Kaiser window；（b）Chebyshev window.

以二階導(dǎo)數(shù)為例，引入誤差函數(shù)：

圖4a表明，Hanning窗和Kaise窗截斷逼近的有限差分算子擁有更大的譜范圍；改進的二項式窗截斷逼近的有限差分算子譜覆蓋范圍較小，Chebyshev窗截斷逼近的有限差分算子譜覆蓋范圍介于它們之間，但更接近Hanning窗和Kaise窗截斷逼近的有限差分算子的譜覆蓋范圍.圖4b是將圖4a的精度誤差放大1000倍的示意圖，從圖4b中可以發(fā)現(xiàn)Hanning窗和Kaiser窗截斷的有限差分算子的逼近精度誤差曲線圍繞零點波動得比較大，而Chebyshev窗和改進的二項式窗截斷逼近的有限差分算子的精度誤差曲線波動最小.這也與前面的分析相符，旁瓣衰減越大，逼近精度誤差出現(xiàn)波動越小，結(jié)果越準確.綜上，在上述提到的四種窗函數(shù)里，折中主瓣寬度和旁瓣衰減，Chebyshev窗是一個較優(yōu)的選擇.

2.3 自褶積組合窗截斷的誤差分析

N是偶數(shù)，將一個長度為N＋1的Chebyshev窗函數(shù)作自褶積運算，可以得到長度為2 N＋1的Chebyshev自褶積窗函數(shù)，同理，將L個相同長度為N＋1的Chebyshev窗函數(shù)作L－1次自褶積，即可得到L階Chebyshev自褶積窗函數(shù).

對Chebyshev窗函數(shù)應(yīng)用傅里葉變換，并令ω＝kxΔx，即可得其幅值響應(yīng)函數(shù)為：

因為空間域的褶積計算相當(dāng)于波數(shù)域的乘積，所以，空間域自褶積L次的Chebyshev窗函數(shù)，在波數(shù)域表現(xiàn)L個Chebyshev窗譜函數(shù)相乘.自褶積之后的窗函數(shù)的幅值響應(yīng)為L個W（ω）相乘.

圖5是不同褶積階數(shù)的Chebyshev自褶積窗幅值響應(yīng)，頻率點M＝512.展示了經(jīng)過自褶積之后，窗函數(shù)的幅值響應(yīng)的變化.進行自褶積的Chebyshev窗長度N＋1＝5，r＝60，將此Chebyshev窗分別作1－5次的自褶積運算，構(gòu)建長度N＋1分別為9，13，17，21，25的Chebyshev自褶積窗；在對Chebyshev窗進行自褶積之后，調(diào)整Chebyshev自褶積窗的系數(shù)為這樣可以使1，且對于不同的n，wc（n）之間的距離也隨之增加.定義1次褶積運算L＝2，依次類推，圖中L最大為6.

隨著褶積次數(shù)的增大，自褶積之后Chebyshev窗函數(shù)的旁瓣衰減得很快，主瓣越來越窄.圖5中綠色線為Chebyshev窗函數(shù)的幅值響應(yīng)，其旁瓣衰減位于－60dB，經(jīng)過一次自褶積后，觀察圖5中青色線，其旁瓣衰減位于－120dB，可見自褶積擁有很好的壓制旁瓣的效果.

如果固定窗函數(shù)的長度N＋1為一個常數(shù)，對于不同的褶積階數(shù)，窗函數(shù)的幅值響應(yīng)見圖6，固定窗函數(shù)長度為N＋1＝25，綠色線為未褶積前長度為25點Chebyshev窗的幅值響應(yīng)曲線，紅色線和青色線分別為自褶積1次的L＝2的Chebyshev自褶積窗和自褶積3次的L＝4的Chebyshev自褶積窗.由褶積原理可知，L＝2的Chebyshev自褶積窗是長度為13的Chebyshev窗自褶積1次生成；L＝4的Chebyshev自褶積窗是長度為7的Chebyshev窗自褶積3次生成.從圖6可以看到，固定窗函數(shù)長度，主瓣寬度取決于自褶積階數(shù)，成一個正比的關(guān)系，階數(shù)越高，主瓣越寬.而隨著自褶積階數(shù)的提高，旁瓣衰減迅速增大，自褶積階數(shù)越高，旁瓣衰減越大，更有效抑制頻譜泄露.

仍以二階導(dǎo)數(shù)為例，使用公式（12）計算圖6所示三種窗函數(shù)截斷逼近的有限差分算子的精度誤差，這時N＝24.

圖4 四種窗函數(shù)截斷逼近的有限差分算子的精度誤差曲線（二階導(dǎo)數(shù)）（N＝24）（a）精度誤差；（b）精度誤差放大1000倍.Fig.4 Accuracy errorfor second derivative of four window truncated finite－difference operators（N＝24）（a）Accuracy error；（b）Magnified 1000times.

圖5 Chebyshev自褶積窗函數(shù)幅值響應(yīng)（N＝4），L＝1，2，3，4，5，6Fig.5 Amplitude response ofChebyshev auto－convolution window function（N＝4），L＝1，2，3，4，5，6

圖6 Chebyshev自褶積窗函數(shù)幅值響應(yīng)（固定N＝24），L＝2，4Fig.6 Amplitude response of Chebyshev auto－convolution window function（Fix N＝24），L＝2，4

圖7所示不同自褶積次數(shù)的Chebyshev自褶積窗截斷逼近的有限差分算子精度誤差曲線，圖7a是三種窗函數(shù)截斷逼近的精度誤差曲線，圖7b是將圖7a放大1000倍之后的精度誤差曲線.其中綠色線是Chebyshev窗函數(shù)截斷逼近的有限差分算子的精度誤差曲線，其擁有更高的譜覆蓋范圍，青色線是自褶積4次的Chebyshev窗函數(shù)截斷逼近的有限差分算子的精度誤差曲線，其譜覆蓋范圍很低；紅色線代表自褶積2次的Chebyshev窗函數(shù)截斷逼近的有限差分算子的精度誤差曲線，其譜覆蓋范圍介于前兩者之間.圖7a中無法看出逼近精度誤差曲線的穩(wěn)定性；圖7b展示了三種窗函數(shù)截斷逼近的有限差分算子精度誤差曲線的穩(wěn)定性，可以觀察到綠色線Chebyshev窗截斷逼近的有限差分算子的精度誤差曲線波動相對較大，紅色線和青色線波動較小.由此，對窗函數(shù)自褶積之后再去截斷偽譜法的空間褶積序列得到的有限差分算子的逼近精度誤差穩(wěn)定性更高，對壓制頻譜泄露有更好的效果.綜合截斷窗函數(shù)的幅值響應(yīng)的主瓣和旁瓣的性能及不同窗函數(shù)截斷逼近的有限差分算子的精度誤差的比較，折中譜覆蓋范圍和逼近精度誤差穩(wěn)定性，窗函數(shù)自褶積2次，可以得到幅值響應(yīng)屬性更優(yōu)的窗函數(shù).

但是，由圖7可知，自褶積2次的Chebyshev窗函數(shù)截斷逼近的有限差分算子的精度誤差雖然有較高的穩(wěn)定性，但是其譜覆蓋范圍小于Chebyshev窗函數(shù)截斷逼近的有限差分算子的譜覆蓋范圍.基于前面的分析，本文的目的是為了設(shè)計一種窗函數(shù)，使截斷得到的有限差分算子在合適的譜覆蓋范圍的前提下，獲得更好的精度誤差穩(wěn)定性.因此，考慮使用加權(quán)函數(shù)，對兩種窗函數(shù)做一個加權(quán)組合.本文選用一個簡單的線性加權(quán)函數(shù).

wc（n）為Chebyshev窗函數(shù)，L為自褶積之后的Chebyshev窗函數(shù)，這里L(fēng)＝2，λ為加權(quán)系數(shù)，λ∈0，［1］.因為根據(jù)前面的分析，窗函數(shù)截斷逼近的有限差分算子的精度主要由窗函數(shù)的幅值響應(yīng)的主瓣和旁瓣特性決定，Chebyshev窗函數(shù)的幅值響應(yīng)有較窄的主瓣，自褶積之后的Chebyshev窗函數(shù)的幅值響應(yīng)則有更大的旁瓣衰減，只要確定λ，就可以確定組合的窗函數(shù).

本文提出的基于Chebyshev窗函數(shù)自褶積組合窗，僅需要三個參數(shù)，就可以直觀地優(yōu)選窗函數(shù)去截斷逼近偽譜法的空間褶積序列得到優(yōu)化的有限差分算子；第一個參數(shù)是Chebyshev窗函數(shù)的紋波率r；第二個參數(shù)是窗函數(shù)自褶積的次數(shù)L；第三個參數(shù)是線性組合的加權(quán)系數(shù)λ.本文選擇λ＝0.89，r＝60，L＝2，圖8展示了不同階的自褶積組合窗函數(shù)截斷逼近的有限差分算子的精度誤差曲線（二階導(dǎo)數(shù)）.觀察圖8a，很明顯，8階的自褶積組合窗函數(shù)截斷逼近的有限差分算子相比常規(guī)8階的有限算子，擁有更高的精度，達到了常規(guī)12階的精度，其精度誤差曲線基本重合.而12階的自褶積組合窗函數(shù)截斷逼近的有限差分算子的準確度遠遠高于常規(guī)12階有限差分算子的精度，甚至超過了常規(guī)24階的精度.圖8b是將圖8a放大1000倍之后的精度誤差曲線，相比常規(guī)有限差分算子的精度誤差曲線，自褶積組合窗函數(shù)截斷逼近的有限差分算子的精度誤差都有波動，但是，其精度誤差控制在0.0004之內(nèi)，擁有較高的精度誤差穩(wěn)定性.

將該自褶積組合窗函數(shù)截斷逼近的有限差分算子與改進的二項式窗函數(shù)截斷逼近的有限差分算子（Chu and Stoffa，2012）的精度誤差曲線相比較，對比結(jié)果見圖9.圖9a說明8階的自褶積組合窗函數(shù)截斷逼近的有限差分算子精度誤差曲線的譜覆蓋范圍不如8階的改進的二項式窗函數(shù)截斷逼近的有限差分算子精度誤差曲線，但是其精度誤差穩(wěn)定性大大提高，8階的改進的二項式窗函數(shù)截斷逼近的有限差分算子精度誤差曲線圍繞零點出現(xiàn)巨大的波動.對于16階和24階的情況，自褶積組合窗函數(shù)截斷逼近的有限差分算子精度誤差曲線譜覆蓋范圍要大于改進的二項式窗函數(shù)截斷逼近的有限差分算子精度誤差曲線；16階的自褶積組合窗函數(shù)截斷逼近的有限差分算子精度高于24階的改進的二項式窗函數(shù)截斷逼近的有限差分算子精度；圖9b是將圖9a放大1000倍之后的精度誤差曲線，自褶積組合窗函數(shù)截斷逼近的有限差分算子在低波數(shù)時的精度誤差曲線波動要小于同條件的改進的二項式窗函數(shù)截斷逼近的有限差分算子的精度誤差曲線.

表1列出了該自褶積組合窗函數(shù)截斷逼近的有限差分算子系數(shù)（二階導(dǎo)數(shù)）.

對于一階導(dǎo)數(shù)的情況，引入誤差函數(shù)：

圖10展示了該自褶積組合窗函數(shù)截斷逼近的有限差分算子在不同階數(shù)時的精度誤差曲線（一階導(dǎo)數(shù)）.對于一階導(dǎo)數(shù)，自褶積組合窗函數(shù)截斷逼近的有限差分算子的精度誤差曲線在頻譜覆蓋范圍和誤差穩(wěn)定性方面相比于常規(guī)有限差分算子也有較好的結(jié)果.

圖7 自褶積窗函數(shù)截斷逼近的有限差分算子精度誤差曲線（二階導(dǎo)數(shù)）（N＝24，L＝2，4）（a）精度誤差曲線；（b）精度誤差曲線放大1000倍.Fig.7 Accuracy errorfor second derivative of auto－convolution window truncated finite－difference operators（N＝24，L＝2，4）（a）Accuracy error；（b）Magnified 1000times.

圖8 不同N的自褶積組合窗函數(shù)截斷逼近的有限差分算子精度誤差曲線（二階導(dǎo)數(shù)）（a）精度誤差曲線；（b）精度誤差曲線放大1000倍.Fig.8 Accuracy errorfor second derivative of auto－convolution window truncated finite－difference operators with different N （different order）（a）Accuracy error；（b）Magnified 1000times.

圖9 不同N的自褶積組合窗函數(shù)截斷逼近的有限差分算子與改進的二項式窗函數(shù)截斷逼近的有限差分算子的精度誤差比較（二階導(dǎo)數(shù)）（a）精度誤差曲線；（b）精度誤差曲線放大1000倍.Fig.9 Comparison of accuracy error for second derivative between auto－convolution window truncated finite－difference operators and improved Binomial window truncated finite－difference operators with different N （different order）（a）Accuracy error；（b）Magnified 1000times.

圖10 不同N的自褶積組合窗函數(shù)截斷逼近的有限差分算子精度誤差曲線（一階導(dǎo)數(shù)）（a）精度誤差曲線；（b）精度誤差曲線放大1000倍.Fig.10 Accuracy error for first derivative of auto－convolution window truncated finite－difference operators with different N （different order）（a）Accuracy error；（b）Magnified 1000times

表1 自褶積組合窗函數(shù)截斷逼近的有限差分算子系數(shù)（二階導(dǎo)數(shù)）Table 1 Auto－convolution window truncated finite－difference operates coefficients for the second derivative

表2列出了該自褶積組合窗函數(shù)截斷逼近的有限差分算子系數(shù)（一階導(dǎo)數(shù)）.

表2 自褶積組合窗函數(shù)截斷逼近的有限差分算子系數(shù)（一階導(dǎo)數(shù)）Table 2 Auto－convolution window truncated finite－difference operates coefficients for the first derivative

3 模型測試

線性彈性理論的假設(shè)下，可以建立應(yīng)變張量和位移之間的聯(lián)系，可得：

εkl是應(yīng)變張量，ul和uk代表位移分量.線性彈性體內(nèi)，應(yīng)變張量和應(yīng)力張量有如下的本構(gòu)關(guān)系：

σij為應(yīng)力張量，cijkl為彈性系數(shù)張量.

從應(yīng)力表示的動力學(xué)平衡方程出發(fā)，不考慮體積力的影響，得到彈性動力學(xué)方程：

在公式（16），（17），（18）上，分別應(yīng)用Chebyshev自褶積組合窗截斷逼近的優(yōu)化有限差分算子和常規(guī)的有限差分算子，進行彈性波數(shù)值模擬.

3.1 各向同性均勻介質(zhì)模型

在這個部分，我們做一個脈沖響應(yīng)的數(shù)值模擬，比較8階Chebyshev自褶積組合窗截斷逼近的優(yōu)化有限差分算子和常規(guī)8階，12階有限差分算子數(shù)值模擬的效果.定義二維各向同性介質(zhì)，網(wǎng)格大小為255×300，網(wǎng)格間距為5m，縱波速度為2000m·s－1，橫波速度為1500m·s－1，ρ＝1000kg·m－3.點源在中間激發(fā)，采用主頻為50Hz的Ricker子波，Δt＝0.5ms，nt＝3000.

圖11是250ms的波場快照.從波場快照上可以清晰地看出，圖11a的X和Z分量都存在較明顯的數(shù)值頻散，圖11b采用我們的優(yōu)化有限差分算子，能夠有效地壓制數(shù)值頻散，相比圖11c采用常規(guī)12階有限差分算子的數(shù)值模擬結(jié)果，我們的優(yōu)化8階有限差分算子有更好的效果.

脈沖響應(yīng)的數(shù)值模擬證明，Chebyshev自褶積組合窗截斷逼近的優(yōu)化有限差分算子在壓制頻散上較常規(guī)方法更有優(yōu)勢.

3.2 Marmousi模型

為了進一步檢測我們的優(yōu)化的有限差分算子，對復(fù)雜的Marmousi模型進行數(shù)值模擬測試.速度模型大小為737×750，橫向采樣間隔為12.5m，縱向采樣間隔為4m采用主頻為30Hz的Ricker子波，震源位置在表面，Δt＝0.5ms，nt＝10000.圖12是 Marmousi速度模型.

采用的GPU型號是GTX 550TI，有192個計算核心和1GB的顯存，CUDA版本是4.2，CPU型號是Intel（R）Core i5 2300＠2.8GHz，配有8GB內(nèi)存.分別應(yīng)用8階，12階，16階，24階的常規(guī)和優(yōu)化的有限差分算子進行計算效率對比測試，并且選擇5個不同的炮點位置，炮點1對應(yīng)x＝0.0m，炮點2對應(yīng)x＝2303.1m，炮點3對應(yīng)x＝4606.3m，炮點4對應(yīng)x＝6909.4m，炮點5對應(yīng)x＝9212.5m.

表3和表4分別列出了不同階數(shù)的優(yōu)化和常規(guī)的有限差分算子在Marmousi模型上的測試時間.

圖13顯示了常規(guī)8階，12階有限差分算子和8階Chebyshev自褶積組合窗截斷逼近的有限差分算子的炮記錄的Z分量，炮點在中間位置.

圖11 Chebyshev自褶積組合窗優(yōu)化有限差分算子和常規(guī)有限分算子脈沖響應(yīng)波場快照對比（250ms）（a）采用常規(guī)8階算子；（b）采用優(yōu)化8階算子；（c）采用常規(guī)12階算子.左圖為X分量，右圖為Z分量.Fig.11 Comparison of snapshot of impulse responses before and after using our optimized operators（250ms）（a）Using the 8th conventional operator；（b）Using the 8th optimized operator；（c）Using the 12th conventional operator.left is Xcomponent，right is Zcomponent.

如果將圖13a整體放大，可以觀察到圖13a左側(cè)圖有明顯的數(shù)值頻散，中間圖采用我們的優(yōu)化算子，數(shù)值頻散得到較大的壓制，相比右側(cè)圖采用常規(guī)12階算子有更好的數(shù)值模擬結(jié)果.現(xiàn)將圖13a中頻散較明顯的三個區(qū)塊放大，繪制圖13b，13c，13d.圖13b是圖13a中區(qū)塊1的放大，可以清晰地觀察到，左側(cè)圖存在比較明顯的頻散現(xiàn)象；而中間圖數(shù)值頻散得到了較好的壓制，且要略優(yōu)于右側(cè)采用常規(guī)12階算子的數(shù)值模擬結(jié)果.圖13c是圖13a中區(qū)塊2的放大，比較圖13c的左中右三圖，可以發(fā)現(xiàn)采用8階Chebyshev自褶積組合窗截斷逼近的優(yōu)化有限差分算子壓制數(shù)值頻散效果更明顯.圖13d是圖13a中區(qū)塊3的放大，圖13d左側(cè)圖采用常規(guī)8階有限差分算子模擬，可以發(fā)現(xiàn)同相軸存在較為明顯的頻散現(xiàn)象，即出現(xiàn)較長的“拖尾”，中間圖采用8階Chebyshev自褶積組合窗截斷逼近的優(yōu)化有限差分算子模擬，很大程度地抑制了數(shù)值頻散，壓制效果是圖13d三圖中最好的.

圖12 Marmousi速度模型Fig.12 Marmousi velocity model

表3 不同階數(shù)的優(yōu)化有限差分算子在Marmousi模型上的計算時間對比Table 3 Comparison of computation cost using different orders auto－convolution window truncated finite－difference operators

Marmousi模型的數(shù)值模擬也證明，Chebyshev自褶積組合窗截斷逼近的優(yōu)化有限差分算子在壓制頻散上較常規(guī)方法更有優(yōu)勢.

4 結(jié)論

本文分析了改進二項式窗函數(shù)和三角函數(shù)類窗函數(shù)的優(yōu)缺點，提出了一種基于Chebyshev自褶積組合窗函數(shù)截斷逼近的有限差分算子優(yōu)化方法.在空間域，有限差分法可看作是偽譜法空間褶積序列的截斷.截斷窗函數(shù)對有限差分算子逼近微分算子精度的影響，從本質(zhì)上來講，是由其幅值響應(yīng)的主瓣和旁瓣屬性所決定的；其主瓣越窄，旁瓣衰減越大，逼近的效果越好，數(shù)值頻散壓制得更徹底.基于Chebyshev自褶積組合窗函數(shù)擁有較好的主旁瓣屬性，與改進的二項式窗相比較，該窗函數(shù)截斷逼近的有限差分算子的精度誤差曲線擁有更大的譜范圍，另外，有比三角函數(shù)類窗函數(shù)截斷逼近的有限差分算子更大的精度誤差穩(wěn)定性，且可只調(diào)節(jié)三個參數(shù)，直觀可視化地控制主瓣和旁瓣的形狀，進而調(diào)整有限差分算子逼近微分算子的精度，這也是窗函數(shù)截斷逼近優(yōu)化的一個優(yōu)點.

表4 不同階數(shù)的常規(guī)有限差分算子在Marmousi模型上的計算時間對比Table 4 Comparison of computation cost using different orders conventional finite－difference operators

本文提出的優(yōu)化有限差分算子比常規(guī)有限差分算子有更高的精度.本文定義了自褶積窗函數(shù)的三個參數(shù)為λ＝0.89，r＝60，L＝2，得到8階的Chebyshev自褶積組合窗截斷逼近的有限差分算子的精度超過了12階常規(guī)有限差分算子；12階的Chebyshev自褶積組合窗截斷逼近的有限差分算子的精度超過了24階常規(guī)有限差分算子；另外，通過該自褶積組合窗截斷逼近的有限差分算子的精度誤差在0.0004之內(nèi)；在脈沖響應(yīng)和Marmousi模型上的測試，都證實了該方案在壓制數(shù)值頻散的有效性.也可以通過選擇不同的窗函數(shù)參數(shù)，調(diào)節(jié)加權(quán)組合系數(shù)和自褶積次數(shù)，達到譜覆蓋范圍和穩(wěn)定性的一個優(yōu)選.

圖13 Chebyshev自褶積組合窗優(yōu)化的有限差分算子和常規(guī)有限差分算子模擬炮記錄對比（Z分量）（a）單炮記錄（左側(cè)圖采用常規(guī)8階算子，中間圖采用優(yōu)化8階算子，右側(cè)圖采用常規(guī)12階算子）；（b）區(qū)塊1放大炮記錄；（c）區(qū)塊2放大炮記錄；（d）區(qū)塊3放大炮記錄.Fig.13 Comparison of a shot record for the Marmousi model computed by the conventional finite－difference operators and auto－convolution window truncated finite－difference operators（Zcomponent）（a）One shot record（left using the 8th conventional operator，middle using the 8th optimized operator，right using the 12th conventional operator）；（b）Zoomed in view of a shot record（Block 1）；（c）Zoomed in view of a shot record（Block 2）；（d）Zoomed in view of a shot record（Block 3）.

Alterman Z，Karal F C Jr.1968.Propagation of elastic waves in layered media by finite difference methods.Bulletinof SeismologicalSocietyofAmerica，58（1）：367－398.

Boris J P，Book D L.1973.Flux－corrected transport.I.SHASTA，a fluid transport algorithm that works.JournalofComputational Physics，11（1）：38－69.

Carcione J M，Kosloff D，Behle A，et al.1992.A spectral scheme for wave propagation simulation in 3－D elastic－anisotropic media.Geophysics，57（12）：1593－1607，doi：10.1190／1.1443227.

Cheng B J，Li X F，Long G H.2008.Seismic waves modeling by convolutional Forsyte polynomial differentiator method.Chinese J.Geophys.（in Chinese），51（2）：531－537.

Chu C L，Stoffa P L.2012.Determination of finite－difference weights using scaled binomial windows.Geophysics，77（3）：W17－W26，doi：10.1190／GEO2011－0336.1.

Dablain M A.1986.The application of high－order differencing to the scalar wave equation.Geophysics，51（1）：54－66，doi：10.1190／1.1442040.

Diniz P S R，da Silva E A B，Netto S L.2012.Digital Signal Processing System Analysis and Design.Beijing：China Machine Press.

Dong L G，Ma Z T，Cao J Z，et al.2000.A staggered－grid highorder difference method of one－order elastic wave equation.Chinese J.Geophys.（in Chinese），43（3）：411－419.

Fornberg B.1987.The pseudospectral method：Comparisons with finite differences for the elastic wave equation.Geophysics，52（4）：483－501，doi：10.1190／1.1442319.

Gazdag J.1981.Modeling of the acoustic wave equation with transform methods.Geophysics，46（6）：854－859，doi：10.1190／1.1441223.

Holberg O.1987.Computational aspects of the choice of operator and sampling interval for numerical differentiation in large－scale simulation of wave phenomena.Geophysical Prospecting，35（6）：629－655，doi：10.1111／j.1365－2478.1987.tb00841.x.

Igel H，Mora P，Riollet B.1995.Anisotropic wave propagation through finite－difference grids.Geophysics，60（4）：1203－1216，doi：10.1190／1.1443849.

Kelly K R，Ward R W，Treitel S，et al.1976.Synthetic seismograms：A finite－difference approach.Geophysics，41（1）：2－27，doi：10.1190／1.1440605.

Kosloff D D，Baysal E.1982.Forward modeling by a Fourier method.Geophysics，47（10）：1402－1412，doi：10.1190／1.1441288.

Lee C，Seo Y.2002.A new compact spectral scheme for turbulence simulations.Journal of Computational Physics，183（2）：438－469，doi：10.1006／jcph.2002.7201.

Liu Y，Li C C，Mou Y G.1998.Finite－difference numerical modeling of any even－order accuracy.OGP （in Chinese），33（1）：1－10.

Liu Y，Sen M K.2009a.A new time－space domain high－order finitedifferent method for the acoustic wave equation.Journal of Computational Physics，228（23）：8779－8806，doi：10.1016／j.jcp.2009.08.027.

Liu Y，Sen M K.2009b.Numerical modeling of wave equation by a truncated high－order finite－difference method. Earthquake Science，22（2）：205－213，doi：10.1007／s11589－009－0205－0.

Liu Y，Sen M K .2010.Acoustic VTI modeling with a time－space domain dispersion－relation－based finite－difference scheme.Geophysics，75（3）：A11－A17，doi：10.1190／1.3374477.

Liu Y.2013.Globally optimal finite－difference schemes based on least squares.Geophysics，78（4）：T113－T132，doi：10.1190／geo2012－0480.1.

Pei Z L，Mou L G.2003.A staggered－grid high－order difference method for modeling seismic wave propagation in inhomogeneous media.Journal of China University of Petroleum （Edition of Natural Science）（in Chinese），27（6）：17－27.

Robertsson J O A，Blanch J O，Symes WW，et al.1994.Galerkinwavelet modeling of wave propagation：Optimal finite－difference stencil design.Mathematical and Computer Modelling，19（1）：31－38，doi：10.1016／0895－7177（94）90113－9.

Yan H Y，Liu Y.2013a.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.

Yan H Y，Liu Y.2013b.Pre－stack reverse－time migration based on the time－space domain adaptive high－order finite－difference method in acoustic VTI medium.Journal of Geophysics and Engineering，10（1）：015010，doi：10.1088／1742－2132／10／1／015010.

Yang D H，Teng J W.1997.FCT finite difference modeling of three－component seismic records in anisotropic medium.OGP（in Chinese），32（2）：181－190.

Yang L，Yan H Y，Liu H.2014.Least squares staggered－grid finite－difference for elastic wave modelling.Exploration Geophysics，45（4）：255－260，doi：10.1071／EG13087.

Zhang J H，Yao Z X.2013.Optimized finite－difference operator for broadband seismic wave modeling.Geophysics，78（1）：A13－A18，doi：10.1190／GEO2012－0277.1.

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

附中文參考文獻

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

董良國，馬在田，曹景忠等.2000.一階彈性波方程交錯網(wǎng)格高階差分解法.地球物理學(xué)報，43（3）：411－419.

劉洋，李承楚，牟永光.1998.任意偶數(shù)階精度有限差分法數(shù)值模擬.石油地球物理勘探，33（1）：1－10.

裴正林，牟永光.2003.非均勻介質(zhì)地震波傳播交錯網(wǎng)格高階有限差分法模擬.石油大學(xué)學(xué)報（自然科學(xué)版），27（6）：17－27.

楊頂輝，騰吉文.1997.各向異性介質(zhì)中三分量地震記錄的FCT有限差分模擬.石油地球物理勘探，32（2）：181－190.