基于余弦調(diào)制Chebyshev窗的彈性波高精度正演

鄭婉秋， 孟小紅*， 劉建紅， 王建

1 中國地質(zhì)大學(xué)(北京)地球物理與信息技術(shù)學(xué)院， 北京　100083　2 中國石油集團東方地球物理勘探有限責(zé)任公司物探技術(shù)研究中心, 河北涿州　072751

?

基于余弦調(diào)制Chebyshev窗的彈性波高精度正演

鄭婉秋1， 孟小紅1*， 劉建紅2， 王建1

1 中國地質(zhì)大學(xué)(北京)地球物理與信息技術(shù)學(xué)院， 北京1000832 中國石油集團東方地球物理勘探有限責(zé)任公司物探技術(shù)研究中心, 河北涿州072751

摘要有限差分時間域正演是彈性波逆時偏移和全波形反演的基礎(chǔ)，正演的計算精度也控制著偏移結(jié)果的準(zhǔn)確性，若精度不高，則在偏移、反演后會帶來假象.為了有效提高正演精度，本文結(jié)合窗函數(shù)優(yōu)化方法，在窗函數(shù)截斷偽譜法空間褶積序列以逼近有限差分算子的基礎(chǔ)上，提出了一種基于Chebyshev窗的余弦調(diào)制模型，在原始Chebyshev窗的基礎(chǔ)上引入了調(diào)制次數(shù)和調(diào)制范圍，通過調(diào)節(jié)這兩個參數(shù)可以人工可視化的調(diào)節(jié)截斷誤差，新的窗函數(shù)繼承了Chebyshev窗的特點，在不明顯降低截斷譜范圍的基礎(chǔ)上明顯降低了截斷誤差.本文針對不同正演階數(shù)N，給出了一組經(jīng)驗調(diào)制系數(shù)，并通過數(shù)值模擬方法，對比了新方法、改進二項式窗和基于最小二乘優(yōu)化方法的正演效果.結(jié)果表明，基于余弦調(diào)制的Chebyshev窗控制數(shù)值頻散的能力更強，在大網(wǎng)格下可以得到更精確的正演結(jié)果.從經(jīng)濟角度分析，該方法減小了計算花費，提高了計算效率.

關(guān)鍵詞有限差分； 數(shù)值頻散； 窗函數(shù)； 彈性波； 余弦調(diào)制

1引言

彈性波逆時偏移是依據(jù)彈性波動理論，以彈性波方程為基礎(chǔ)的逆時偏移方法，彈性矢量波場逆時偏移技術(shù)采用全波波動方程，可獲得更加豐富的細節(jié)信息，不論在各向同性介質(zhì)還是各向異性介質(zhì)，彈性波對地下介質(zhì)的成像效果要優(yōu)于單純聲波的成像效果(Chang and McMechan,1987; Sun and McMechan,2001).而彈性波正演是彈性波逆時偏移的基礎(chǔ)，準(zhǔn)確的正演可以避免成像假象，提高成像精度.常用的數(shù)值模擬分為三類：幾何射線法、積分方程法和波動方程法，波動方程法的模擬結(jié)果包含了地震波的運動學(xué)和動力學(xué)特征，是最常用的(程冰潔等，2008).

較多使用的時域模擬方法有偽譜法(Gazdag,1981;Kosloff and Baysal,1982;Carcione et al.,1992)和有限差分方法(Alterman and Karal,1968;Kelly et al.,1976;Dablain,1986;Igel et al.,1995).

偽譜法的計算精度較高，在時間軸上其采用有限差分方法進行近似，而在空間方向采用傅里葉變換求解，理論上其利用了計算區(qū)域的所有點，這樣的算法避免了網(wǎng)格頻散的現(xiàn)象，但是大大增加了正演耗時(Gazdag，1981；Kosloff and Baysal，1982；Carcione et al.,1992).有限差分方法即在時間方向和空間方向皆采用差分近似微分的方法，自20世紀(jì)60年代起就被地球物理學(xué)家廣泛使用，逐步發(fā)展出常規(guī)網(wǎng)格有限差分、普通交錯網(wǎng)格有限差分、旋轉(zhuǎn)交錯網(wǎng)格有限差分.Alterman和Karal (1968)首次利用有限差分法進行彈性波正演.Madariaga(1976)首次實現(xiàn)了基于交錯網(wǎng)格的彈性波數(shù)值模擬.Saenger等(2002，2004)提出了適應(yīng)于各向異性介質(zhì)的旋轉(zhuǎn)交錯網(wǎng)格有限差分方法.劉洋等(1998)實現(xiàn)了任意偶數(shù)階彈性波的正演模擬.董良國等(2000)實現(xiàn)了一階彈性波的交錯網(wǎng)格有限差分.但是對于有限差分?jǐn)?shù)值模擬來說，數(shù)值頻散是一個不可回避的問題.隨著震源子波主頻的升高，空間采樣間隔增大，數(shù)值頻散會愈發(fā)嚴(yán)重.根本上來講，利用有限點的差分近似代替連續(xù)的微分算子，在截斷后就產(chǎn)生了數(shù)值頻散.降低數(shù)值頻散的簡單辦法就是采用更細致的網(wǎng)格并降低子波主頻，但是這會帶來海量的計算.前人在如何降低數(shù)值頻散方面做了諸多工作，其一是通量校正技術(shù)(FCT)，自Boris 和Book (1973)首先提出通量校正技術(shù)(FCT)后，F(xiàn)ei和Larner(1995)將該技術(shù)推廣到彈性波正演，隨后楊頂輝和騰吉文(1997)將FCT技術(shù)運用到各向異性介質(zhì)的模擬中.其二是通過最優(yōu)化方法，如模擬退火法(Zhang and Yao，2013)、最小二乘法(Yang，2014)等.再有就是本文所利用的窗函數(shù)法，據(jù)Fornberg (1987)的研究，偽譜法可以表示為有限差分法的高階近似，換而言之，即利用時間域窗函數(shù)可以截斷偽譜法序列以得到有限差分系數(shù).前人利用了不同的窗函數(shù)以得到優(yōu)先差分系數(shù)，如加權(quán)Hanning窗(Zhou and Greenhalgh，1992)、高斯窗(Igel，1995)、改進二項式窗(Chu and Stoffa，2012)、基于自褶積的Chebyshev窗(王之洋等，2015).

在前人的研究基礎(chǔ)上，本文提出了一種基于余弦調(diào)制的Chebyshev窗，即將Chebyshev窗與一定范圍的余弦函數(shù)相乘數(shù)次.這里的范圍即調(diào)制范圍，次數(shù)即調(diào)制次數(shù).研究表明，經(jīng)過余弦調(diào)制后的Chebyshev窗提高了原有窗函數(shù)的截斷效果，在不明顯降低截斷譜范圍的基礎(chǔ)上，大大提高了截斷誤差的穩(wěn)定性.通過實驗分析，本文給出了針對不同階數(shù)的調(diào)制參數(shù)，并通過數(shù)值模擬證明，該方法可適用于較大網(wǎng)格的正演模擬，且在壓制數(shù)值頻散方面較其他窗函數(shù)和優(yōu)化方法有較大的優(yōu)勢，從經(jīng)濟上講，我們可以利用改進后的低階模擬代替原始高階計算，大大降低了計算花費.

2窗函數(shù)逼近法

2.1差分算子及截斷誤差

據(jù)前人研究表明，有限差分算子可以通過sinc函數(shù)插值理論推導(dǎo)得到，sinc函數(shù)可以重構(gòu)均勻采樣后的帶限信號(Diniz et al.,2012),sinc函數(shù)可以表示為

(1)

(2)

(3)

截斷后，可以表示為

(4)

(5)

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

(6)

(8)

(9)

將公式(8)和公式(9)變換到頻率域后為

(10)

(11)

對于一階導(dǎo)數(shù)，可以將誤差函數(shù)表示成

(12)

對于二階導(dǎo)數(shù)，可以將誤差函數(shù)表示成

(13)

2.2傳統(tǒng)窗函數(shù)下的逼近

如今信號領(lǐng)域廣泛使用的窗函數(shù)包括：Hanning窗、Kaiser窗、Chebyshev窗等，而在利用窗函數(shù)截斷逼近有限差分算子方面，也有不同學(xué)者使用了不同的窗函數(shù)，如：廣義加權(quán)的Hanning窗(Zhou and Greenhalgh，1992)、高斯窗(Igel，1995)、基于自褶積的Chebyshev窗(王之洋，2015)和本文將重點對比的改進二項式窗(Chu and Stoffa，2012)，二項式窗的窗函數(shù)形如公式(14).

(14)

Chu和Stoffa認為，當(dāng)N=N+M時，窗函數(shù)的形式會發(fā)生改變，其截斷效果更佳，可以更好地控制數(shù)值頻散，而當(dāng)M=0時，截斷參數(shù)為傳統(tǒng)的有限差分算子系數(shù)，從這個角度來說，改進二項式窗是有限差分法和偽譜法連接的橋梁.

在空間域，將偽譜法的空間褶積序列截斷就得到了有限差分法，如果從信號處理層面來講，截斷相當(dāng)于加了一個窗函數(shù)(王之洋等，2015).這就是說，我們可以將時間域的窗函數(shù)看成一個有限長的濾波器，濾波器主瓣寬度和旁瓣衰減也就決定了其性質(zhì).Diniz等(2012)的分析結(jié)論表明，主瓣寬度與截斷的譜范圍呈負相關(guān)關(guān)系，即主瓣越寬，則截斷的譜范圍越低.而旁瓣衰減與截斷誤差的穩(wěn)定性呈正相關(guān)關(guān)系，即旁瓣衰減越大，誤差穩(wěn)定性越高.

圖1是不同窗函數(shù)在N=24時的時間域展布，圖2為各窗函數(shù)的幅值響應(yīng).從圖中可以發(fā)現(xiàn)，當(dāng)時間域窗函數(shù)系數(shù)變化率較高時，其幅值響應(yīng)有著更寬的主瓣，而窗函數(shù)尾端曲率越大，則旁瓣衰減越明顯.在圖1所示各種窗函數(shù)中，改進二項式窗系數(shù)有著最快的下降速率，也對應(yīng)了最寬的主瓣范圍，但也擁有最大的旁瓣衰減.對于Chebyshev窗，當(dāng)參數(shù)γ增大時，其主瓣將變寬緩，旁瓣衰減變大.Kaiser窗與Chebyshev窗有著類似的性質(zhì).

圖1　窗函數(shù)對比Fig.1　Comparison of window functions

圖2　各窗函數(shù)的幅值響應(yīng)Fig.2　Amplitude responses of window functions

綜合以上分析可以得出結(jié)論：若想有效地控制數(shù)值頻散，應(yīng)該尋求主瓣較窄、旁瓣衰減較大的窗函數(shù)截斷褶積序列以得到有限差分系數(shù).圖3顯示了圖1中所示各種傳統(tǒng)窗函數(shù)在二階導(dǎo)數(shù)下的截斷誤

圖3　傳統(tǒng)窗函數(shù)二階導(dǎo)數(shù)下的截斷誤差 (a) 截斷誤差對比； (b) 將(a)縱軸放大1000倍結(jié)果.Fig.3　Truncation errors for second derivative of traditional window function (a) Comparison of truncation errors； (b) The longitudinal axis of (a) amplified 1000 times.

差，可以看出，Kaiser窗有著最大的截斷譜范圍，但從圖3b中看出，其截斷誤差的穩(wěn)定性是幾種窗函數(shù)中最差的，這一點也印證了我們之前對主瓣范圍和旁瓣衰減在誤差控制方面的分析.在另外三種窗函數(shù)中，改進二項式窗的誤差穩(wěn)定性最好，但是譜范圍也是最小的，這就意味著其對高波數(shù)部分的逼近效果較差，會產(chǎn)生較強的數(shù)值頻散.Hanning窗與Chebyshev窗的譜范圍大致相同，其在低波數(shù)范圍的誤差較小，但是其誤差隨波數(shù)的增加而迅速增加，在高波數(shù)部分表現(xiàn)出了強烈的動蕩，從這個角度來講，雖然Chebyshev窗在整個區(qū)域內(nèi)都存在著誤差，但其較為穩(wěn)定，不存在誤差的急劇變化，在這幾種傳統(tǒng)窗函數(shù)中為較優(yōu)的選擇.

2.3余弦調(diào)制Chebyshev窗下的逼近

所謂余弦調(diào)制，即將Chebyshev窗與一定范圍的余弦函數(shù)在時間域相乘數(shù)次，可以表示為

(15)

其中，wdcs代表余弦調(diào)制后的窗函數(shù)，wdche表示Chebyshev窗函數(shù)，e表示調(diào)制次數(shù)，p表示調(diào)制范圍，N+1表示窗函數(shù)長度.可以看出，我們在余弦調(diào)制的過程中定義了兩個參數(shù)，即調(diào)制范圍和調(diào)制次數(shù).

圖4顯示了不同調(diào)制參數(shù)下窗函數(shù)的幅值響應(yīng)，可以看出，當(dāng)p值增大，e值不變時，窗函數(shù)的主瓣寬度逐漸減小，旁瓣衰減變化不大.而當(dāng)e值增大，p值不變時，主瓣寬度有所增加，但旁瓣衰減也逐步增大.也就是說，可以通過對這兩個參數(shù)的動態(tài)選擇，以控制所調(diào)制窗函數(shù)的主瓣寬度和旁瓣衰減，最終達到優(yōu)化原有Chebyshev窗、優(yōu)化有限差分系數(shù)的目的.以24階有限差分算子為例，計算了不同調(diào)制范圍和調(diào)制次數(shù)下的二階導(dǎo)數(shù)截斷誤差，如圖5.可以看出，調(diào)制范圍相同，調(diào)制次數(shù)增加時，截斷譜范圍降低、誤差穩(wěn)定性增強.而當(dāng)調(diào)制次數(shù)不變，調(diào)制范圍增大時，誤差穩(wěn)定性降低、譜范圍增大.所以調(diào)制次數(shù)和調(diào)制范圍為兩個相互制約且起反作用的參數(shù).為平衡兩者，本文采用的做法是對不同階數(shù)進行多次誤差分析和正演模擬試驗，以選擇較為合適的經(jīng)驗調(diào)制參數(shù)，表1為通過實驗分析得到的不同階數(shù)下截斷效果較好的調(diào)制范圍參數(shù)和調(diào)制次數(shù)參數(shù).

圖4　不同調(diào)制參數(shù)下的窗函數(shù)幅值響應(yīng)Fig.4　Amplitude responses for different modulation parameters of window function

通過實驗分析，在N+1=25時確定p=3、e=3有著較優(yōu)的截斷效果,而N+1=9時確定p=3.5、e=1有較優(yōu)的截斷效果.圖5為不同調(diào)制參數(shù)下一階導(dǎo)數(shù)的截斷誤差.

圖5　不同調(diào)制參數(shù)下的截斷誤差Fig.5　Truncation errors for different modulation parameters

從圖6中可以看出，對于高階數(shù)N=24時，雖然調(diào)制后的窗函數(shù)犧牲了一定的譜范圍，但是余弦調(diào)制Chebyshev窗在控制誤差穩(wěn)定性上有了很大的提升，通過和改進二項式窗(M=16)的對比可見，新改造的窗函數(shù)對誤差穩(wěn)定性的控制已不亞于改進二項式窗，而且有更寬的譜范圍.圖7顯示了對于低階N=8時，幾種窗函數(shù)的截斷效果，由圖中可以看出，改進二項式窗(M=16)和Chebyshev窗出現(xiàn)了不穩(wěn)定的情況，而此時的余弦調(diào)制窗函數(shù)則依舊保持了誤差的穩(wěn)定性，其截斷的譜范圍也沒有明顯降低.綜合圖3、圖5及圖6分析并結(jié)合正演階數(shù)、截斷譜范圍和誤差穩(wěn)定性三方面綜合考慮，基于余弦調(diào)制Chebyshev窗是最優(yōu)的選擇.表2及表3分別是一階導(dǎo)數(shù)、二階導(dǎo)數(shù)下的有限差分系數(shù)，圖8中整合了不同階數(shù)下優(yōu)化有限差分算子和基于自褶積Chebyshev窗(王之洋等，2015)截斷的有限差分算子的二階導(dǎo)數(shù)截斷誤差，從圖中不難發(fā)現(xiàn)，余弦調(diào)制后的Chebyshev窗有更好的誤差穩(wěn)定性.

圖6　調(diào)制窗函數(shù)截斷誤差效果(N=24) (a) 截斷誤差對比； (b) 將(a)縱軸放大1000倍結(jié)果.Fig.6　Truncation errors of modulating window function (N=24) (a) Comparison of truncation errors； (b) The longitudinal axis of (a) amplified 1000 times.

圖7　調(diào)制窗函數(shù)截斷誤差效果(N=8) (a) 截斷誤差對比； (b) 將(a)縱軸放大1000倍結(jié)果.Fig.7　Truncation errors of modulating window function (N=8) (a) Comparison of truncation errors； (b) The longitudinal axis of (a) amplified 1000 times.

4階8階12階16階p5.53.544e0.4123

圖8　余弦調(diào)制Chebyshev窗和自褶積Chebyshev窗二階導(dǎo)數(shù)截斷誤差對比Fig.8　Comparison of truncation errors between modulated Chebyshev window and auto-convolution Chebyshev window

3模型測試

對于有限差分法，正演時要對時間軸和空間軸分別進行差分表示，再按照時間軸方向遞推，據(jù)此本文進行了基于應(yīng)變-位移方程的各向同性和各向異性彈性波模擬，其中時間軸二階導(dǎo)數(shù)可以表示為

(16)

(17)

對于空間軸，本文采用常用的2N階離散格式：

表2　一階導(dǎo)數(shù)下的有限差分系數(shù)

表3　二階導(dǎo)數(shù)下的有限差分系數(shù)

(18)

(19)

在公式(16)—(19)中，ux和uz代表離散的時間域波場，am為有限差分系數(shù)，Δx和Δz代表空間步長.

可以說，彈性波動理論的基礎(chǔ)是本構(gòu)方程、運動微分方程和幾何方程，想得到最終的離散正演公式，只需將公式(16)—(19)帶入這三個方程即可，以x方向為例，最終的離散公式可以表示為

(20)

3.1簡單模型

對于均勻模型的彈性波數(shù)值模擬是驗證窗函數(shù)截斷所得到有限差分系數(shù)質(zhì)量的有力手段，我們以40 Hz Ricker子波為震源，將震源置于模型中間，模型縱波速度設(shè)定為2000 m·s-1，橫波速度為1300 m·s-1，模型的空間采樣間隔設(shè)定為10 m，網(wǎng)格大小為600×600，并采用海綿吸收邊界改善模型邊界反射.本文測試了改進二項式窗8階、余弦調(diào)制窗8階、改進二項式窗12階算子在各向同性(ISO)和各向異性(VTI)介質(zhì)模型下的正演效果，如圖9，數(shù)值驗證表明：余弦調(diào)制窗8階的ISO介質(zhì)正演效果遠遠超過了改進二項式窗8階的正演效果，其控制數(shù)值頻散的能力甚至超過了改進二項式窗12階算子，而新方法8階算子的VTI正演效果明顯優(yōu)于改進二項式窗12階算子，如圖10.數(shù)值模擬的結(jié)果與誤差譜范圍的覆蓋及誤差穩(wěn)定性分析結(jié)果有著統(tǒng)一的結(jié)論，在誤差允許范圍內(nèi)，譜范圍越大，誤差穩(wěn)定性越強，控制數(shù)值頻散的能力越好.這就意味著較低階的有限差分算子可以代替高階有限差分算子，在控制誤差的同時，大大節(jié)約了計算成本.

圖9　均勻ISO介質(zhì)測試 (a) 8階改進二項式窗波場快照X分量; (b) 8階改進二項式窗波場快照Z分量; (c) 8階余弦調(diào)制窗波場快照X分量; (d) 8階余弦調(diào)制窗波場快照Z分量; (e) 12階改進二項式窗波場快照X分量; (f) 12階改進二項式窗波場快照Z分量.Fig.9　Tests of homogeneous isotropic media (a) Using the 8th improved binomial window function (X component); (b) Using the 8th improved binomial window function (Z component); (c) Using the 8th modulated window function (X component); (d) Using the 8th modulated window function (Z component); (e) Using the 12th improved binomial window function (X component); (f) Using the 12th improved binomial window function (Z component).

圖10　均勻VTI介質(zhì)測試 (a) 12階改進二項式窗波場快照X分量; (b) 12階改進二項式窗波場快照Z分量; (c) 8階余弦調(diào)制窗波場快照X分量; (d) 8階余弦調(diào)制窗波場快照Z分量.Fig.10　Tests of homogeneous anisotropic media (a) Using the 12th improved binomial window function (X component); (b) Using the 12th improved binomial window function (Z component); (c) Using the 8th modulated window function (X component); (d) Using the 8th modulated window function (Z component).

3.2BP模型

為了進一步測試基于余弦調(diào)制有限差分算子的穩(wěn)定性，我們在復(fù)雜的BP模型上進行了彈性波數(shù)值模擬.對比了常規(guī)有限差分算子、基于最小二乘優(yōu)化方法(Yang，2014)和本文余弦調(diào)制方法在正演時的穩(wěn)定性.模型規(guī)模為3792×1896，空間采樣間隔為6.25 m，時間采樣間隔為0.5 ms，時間采樣點數(shù)為14000.震源位于表層中間位置，Ricker子波主頻為30 Hz，圖11和圖12分別為BP縱波速度模型和密度模型.

圖11　BP縱波速度模型Fig.11　BP longitudinal wave velocity model

圖12　BP密度模型Fig.12　BP density model

(21)

其中VS為橫波速度，VP為縱波速度.

我們將正演算法并行化，采用的GPU型號為GTX960，該圖形處理器擁有4 GB顯存，192個計算核心，每個流處理器可最多容納2048個線程.CPU為Inter(R)i7-4790k@4GHZ，計算機內(nèi)存為8 G.表4顯示了相同炮點不同階數(shù)的常規(guī)有限差分算子和優(yōu)化有限差分算子在BP模型上的計算耗時，可以看出，優(yōu)化后的有限差分算子并沒有增加過多的計算花費.

表4　各階算子計算效率

從正演道記錄的局部細節(jié)可以清晰地看出，基于余弦調(diào)制的Chebyshev窗在抑制數(shù)值頻散方面有突出的優(yōu)勢.由圖13d，13e，13f可見，在直達波下，傳統(tǒng)方法出現(xiàn)了較明顯的“拖尾”現(xiàn)象，而基于余弦調(diào)制的Chebyshev窗很好地控制了直達波下的數(shù)值頻散，且控制頻散的能力強于最小二乘優(yōu)化方法.對比圖13g，13h，13i中有明顯反射的局部細節(jié)，可以發(fā)現(xiàn)余弦調(diào)制方法可以使有效反射更加清晰，其抑制頻散的能力也是這三種方法中最好的.

圖13　BP模型炮記錄 (a) 常規(guī)8階算子； (b) 最小二乘8階算子； (c) 調(diào)制窗8階算子； (d)(e)(f)分別為(a)(b)(c)左側(cè)紅框相同位置的細節(jié)放大； (g)(h)(i)分別為(a)(b)(c)右側(cè)紅框相同位置的細節(jié)放大.Fig.13　Shot records for the BP model (a) Conventional 8th operator; (b) Least squares 8th operator; (c) Modulated 8th operator; (d)(e)(f) Zoomed in left red box of (a)(b)(c); (g)(h)(i) Zoomed in right red box of (a)(b)(c).

4結(jié)論

本文在窗函數(shù)截斷偽譜法空間褶積序列的基礎(chǔ)上，分析了現(xiàn)有的多種三角類窗函數(shù)和改進二項式窗函數(shù)，對比分析不同窗函數(shù)的截斷效果后，提出了一種基于余弦調(diào)制的Chebyshev窗.從時間域來看，窗函數(shù)可以看成一個有限長的濾波器，濾波器主瓣寬度和旁瓣衰減決定了其性質(zhì).具體來講，主瓣寬度控制了截斷的頻譜覆蓋范圍，主瓣越窄，截斷的譜范圍越寬，即負相關(guān)關(guān)系.而旁瓣衰減則控制了誤差的穩(wěn)定性，即偏差程度，旁瓣衰減越大，誤差穩(wěn)定性越好.通過對幾種傳統(tǒng)窗函數(shù)的截斷誤差分析我們發(fā)現(xiàn)，無論從截斷的譜范圍還是誤差穩(wěn)定性來看，Chebyshev窗都是一個較優(yōu)的選擇，但是實際上，其誤差穩(wěn)定性仍較改進二項式窗差.

據(jù)此，本文提出的基于余弦調(diào)制的Chebyshev窗在保證原始截斷譜范圍的基礎(chǔ)上，大大增強了Chebyshev窗截斷誤差的穩(wěn)定性，使其提升了兩個數(shù)量級.在余弦調(diào)制的過程中，我們引入了兩個參數(shù)，分別為調(diào)制范圍和調(diào)制次數(shù)，調(diào)制范圍主要控制了窗函數(shù)的主瓣范圍，而調(diào)制次數(shù)影響了旁瓣衰減.經(jīng)過對不同調(diào)制參數(shù)進行誤差分析和數(shù)值模擬，本文總結(jié)了一套經(jīng)驗系數(shù).

通過對均勻模型的試算，驗證了在理論誤差分析中得出的結(jié)論：余弦調(diào)制窗8階的均勻介質(zhì)正演效果遠遠超過了改進二項式窗8階的正演效果，其控制數(shù)值頻散的能力甚至超過了改進二項式窗12階算子，而新方法8階算子的VTI介質(zhì)正演效果明顯優(yōu)于改進二項式窗12階算子，對于BP模型的模擬也驗證了本文方法的有效性.在保證誤差穩(wěn)定的前提下，可以以低階數(shù)的運算代替較高階數(shù)的運算，這減小了計算耗時，提高了經(jīng)濟效率.

References

Alterman Z, Karal F C Jr. 1968. Propagation of elastic waves in layered media by finite difference methods. Bull. Seism. Soc. Am., 58(1): 367-398.

Boris J P, Book D L. 1973. Flux-corrected transport. I. SHASTA, a fluid transport algorithm that works. Journal of Computational 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.

Chang W F, McMechan G A. 1987. Elastic reverse time migration. Geophysics,52:1367-1375.

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 high-order difference method of one-order elastic wave equation. Chinese J. Geophys. (in Chinese), 43(3): 411-419.

Fei T, Larner K. 1995. Elimination of numerical dispersion in finite-difference modeling and migration by flux-corrected transport. Geophysics, 60(6):1830-1842.

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.

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. Oil Geophysical Prospecting (in Chinese), 33(1): 1-10.

Madariaga R. 1976. Dynamics of an expanding circular fault. Bull. Seism. Soc. Am., 65: 163-188.

Saenger E H, Shapiro S. 2002. An effective velocities in fractured media: A numerical study using the rotated staggered finite-difference grid. Geophysical Prospecting, 50: 183-194.

Saenger E H，Thomas B. 2004. Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rots-ted staggered grid. Geophysics, 69(2): 583-591.

Sun R, McMechan G A. 2001. Scalar reverse-time depth migration of prestack elastic seismic data. Geophysics, 66 (5): 1519-1527.

Wang Z Y, Liu H, Tang X D, et al. 2015. Optimized finite-difference operators based on Chebyshev auto-convolution combined window function. Chinese J. Geophys. (in Chinese), 58(2): 628-642, doi: 10.6038/cjg20150224.

Yang D H, Teng J W. 1997. FCT finite difference modeling of three-component seismic records in anisotropic medium. Oil Geophysical Prospecting (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: 255-260.

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. Bull. Seism. Soc.Am., 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.

王之洋, 劉洪, 唐祥德等. 2015. 基于Chebyshev自褶積組合窗的有限差分算子優(yōu)化方法. 地球物理學(xué)報, 58(2): 628-642, doi: 10.6038/cjg20150224.

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

(本文編輯何燕)

基金項目國家重大科研裝備研制項目(ZDYZ2012-1-02-04)和國家自然科學(xué)基金(41474106)聯(lián)合資助.

作者簡介鄭婉秋，在讀碩士，主要從事復(fù)雜介質(zhì)中地震波模擬及偏移成像方面的研究.E-mail:453294543@qq.com *通訊作者孟小紅，教授，主要從事計算與綜合地球物理研究.E-mail:mxh@cugb.edu.cn

doi:10.6038/cjg20160728 中圖分類號P631

收稿日期2015-11-30，2016-06-11收修定稿

High precision elastic wave equation forward modeling based on cosine modulated Chebyshev window function

ZHENG Wan-Qiu1, MENG Xiao-Hong1*， LIU Jian-Hong2, WANG Jian1

1SchoolofGeophysicsandInformationTechnology,ChinaUniversityofGeosciences,Beijing100083,China2BGPResearchandDevelopmentCenter,CNPC,HebeiZhuozhou072751,China

AbstractThe finite difference forward modeling is the basis of elastic wave reverse-time migration and full waveform inversion in the time domain. The accuracy of forward modeling also controls the accuracy of seismic imaging and inversion. The migration or inversion will bring illusion if the accuracy is not high. We can get optimized explicit finite difference operators by using the window function to truncate spatial convolution counterpart of the pseudo-spectral method. Based on this, a cosine modulated Chebyshev window is designed. On the basis of the original Chebyshev window, the modulation times and modulation domain are introduced, and we can adjust truncation error visually by controlling these two parameters. As the new window function inherits the character of Chebyshev window, we observe that the spectral range using the modulated window function for truncation is significantly broader than using the conventional window function with stable error. For different forward modeling orders N, we give a set of empirical modulation factors and compare the forward modeling effect of the new method and improved binomial window by the numerical simulation method. The results demonstrate that the operators based on the cosine modulated Chebyshev window can efficiently suppress the numerical dispersion and get more accurate forward modeling results on the large grid. From economic perspective, this method reduces the computational cost and improves efficiency.

KeywordsFinite difference; Numerical dispersion; Window function； Elastic wave； Cosine modulation

鄭婉秋， 孟小紅， 劉建紅等. 2016. 基于余弦調(diào)制Chebyshev窗的彈性波高精度正演.地球物理學(xué)報，59(7):2650-2662,doi:10.6038/cjg20160728.

Zheng W Q, Meng X H, Liu J H, et al. 2016. High precision elastic wave equation forward modeling based on cosine modulated Chebyshev window function. Chinese J. Geophys. (in Chinese),59(7):2650-2662,doi:10.6038/cjg20160728.