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

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

?

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

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

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

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

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

1引言

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

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

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

2窗函數(shù)逼近法

2.1差分算子及截?cái)嗾`差

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

(1)

(2)

(3)

截?cái)嗪?，可以表示?/p>

(4)

(5)

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

(6)

(8)

(9)

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

(10)

(11)

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

(12)

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

(13)

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

如今信號(hào)領(lǐng)域廣泛使用的窗函數(shù)包括：Hanning窗、Kaiser窗、Chebyshev窗等，而在利用窗函數(shù)截?cái)啾平邢薏罘炙阕臃矫妫灿胁煌瑢W(xué)者使用了不同的窗函數(shù)，如：廣義加權(quán)的Hanning窗(Zhou and Greenhalgh，1992)、高斯窗(Igel，1995)、基于自褶積的Chebyshev窗(王之洋，2015)和本文將重點(diǎn)對(duì)比的改進(jìn)二項(xiàng)式窗(Chu and Stoffa，2012)，二項(xiàng)式窗的窗函數(shù)形如公式(14).

(14)

Chu和Stoffa認(rèn)為，當(dāng)N=N+M時(shí)，窗函數(shù)的形式會(huì)發(fā)生改變，其截?cái)嘈Ч?，可以更好地控制?shù)值頻散，而當(dāng)M=0時(shí)，截?cái)鄥?shù)為傳統(tǒng)的有限差分算子系數(shù)，從這個(gè)角度來說，改進(jìn)二項(xiàng)式窗是有限差分法和偽譜法連接的橋梁.

在空間域，將偽譜法的空間褶積序列截?cái)嗑偷玫搅擞邢薏罘址?，如果從信?hào)處理層面來講，截?cái)嘞喈?dāng)于加了一個(gè)窗函數(shù)(王之洋等，2015).這就是說，我們可以將時(shí)間域的窗函數(shù)看成一個(gè)有限長的濾波器，濾波器主瓣寬度和旁瓣衰減也就決定了其性質(zhì).Diniz等(2012)的分析結(jié)論表明，主瓣寬度與截?cái)嗟淖V范圍呈負(fù)相關(guān)關(guān)系，即主瓣越寬，則截?cái)嗟淖V范圍越低.而旁瓣衰減與截?cái)嗾`差的穩(wěn)定性呈正相關(guān)關(guān)系，即旁瓣衰減越大，誤差穩(wěn)定性越高.

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

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

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

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

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

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

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

(15)

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

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

通過實(shí)驗(yàn)分析，在N+1=25時(shí)確定p=3、e=3有著較優(yōu)的截?cái)嘈Ч?而N+1=9時(shí)確定p=3.5、e=1有較優(yōu)的截?cái)嘈Ч?圖5為不同調(diào)制參數(shù)下一階導(dǎo)數(shù)的截?cái)嗾`差.

圖5　不同調(diào)制參數(shù)下的截?cái)嗾`差Fig.5　Truncation errors for different modulation parameters

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

圖6　調(diào)制窗函數(shù)截?cái)嗾`差效果(N=24) (a) 截?cái)嗾`差對(duì)比； (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ù)截?cái)嗾`差效果(N=8) (a) 截?cái)嗾`差對(duì)比； (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ù)截?cái)嗾`差對(duì)比Fig.8　Comparison of truncation errors between modulated Chebyshev window and auto-convolution Chebyshev window

3模型測試

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

(16)

(17)

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

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

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

(18)

(19)

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

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

(20)

3.1簡單模型

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

圖9　均勻ISO介質(zhì)測試 (a) 8階改進(jìn)二項(xiàng)式窗波場快照X分量; (b) 8階改進(jìn)二項(xiàng)式窗波場快照Z分量; (c) 8階余弦調(diào)制窗波場快照X分量; (d) 8階余弦調(diào)制窗波場快照Z分量; (e) 12階改進(jìn)二項(xiàng)式窗波場快照X分量; (f) 12階改進(jìn)二項(xiàng)式窗波場快照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階改進(jìn)二項(xiàng)式窗波場快照X分量; (b) 12階改進(jìn)二項(xiàng)式窗波場快照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模型

為了進(jìn)一步測試基于余弦調(diào)制有限差分算子的穩(wěn)定性，我們在復(fù)雜的BP模型上進(jìn)行了彈性波數(shù)值模擬.對(duì)比了常規(guī)有限差分算子、基于最小二乘優(yōu)化方法(Yang，2014)和本文余弦調(diào)制方法在正演時(shí)的穩(wěn)定性.模型規(guī)模為3792×1896，空間采樣間隔為6.25 m，時(shí)間采樣間隔為0.5 ms，時(shí)間采樣點(diǎn)數(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型號(hào)為GTX960，該圖形處理器擁有4 GB顯存，192個(gè)計(jì)算核心，每個(gè)流處理器可最多容納2048個(gè)線程.CPU為Inter(R)i7-4790k@4GHZ，計(jì)算機(jī)內(nèi)存為8 G.表4顯示了相同炮點(diǎn)不同階數(shù)的常規(guī)有限差分算子和優(yōu)化有限差分算子在BP模型上的計(jì)算耗時(shí)，可以看出，優(yōu)化后的有限差分算子并沒有增加過多的計(jì)算花費(fèi).

表4　各階算子計(jì)算效率

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

圖13　BP模型炮記錄 (a) 常規(guī)8階算子； (b) 最小二乘8階算子； (c) 調(diào)制窗8階算子； (d)(e)(f)分別為(a)(b)(c)左側(cè)紅框相同位置的細(xì)節(jié)放大； (g)(h)(i)分別為(a)(b)(c)右側(cè)紅框相同位置的細(xì)節(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ù)截?cái)鄠巫V法空間褶積序列的基礎(chǔ)上，分析了現(xiàn)有的多種三角類窗函數(shù)和改進(jìn)二項(xiàng)式窗函數(shù)，對(duì)比分析不同窗函數(shù)的截?cái)嘈Ч?，提出了一種基于余弦調(diào)制的Chebyshev窗.從時(shí)間域來看，窗函數(shù)可以看成一個(gè)有限長的濾波器，濾波器主瓣寬度和旁瓣衰減決定了其性質(zhì).具體來講，主瓣寬度控制了截?cái)嗟念l譜覆蓋范圍，主瓣越窄，截?cái)嗟淖V范圍越寬，即負(fù)相關(guān)關(guān)系.而旁瓣衰減則控制了誤差的穩(wěn)定性，即偏差程度，旁瓣衰減越大，誤差穩(wěn)定性越好.通過對(duì)幾種傳統(tǒng)窗函數(shù)的截?cái)嗾`差分析我們發(fā)現(xiàn)，無論從截?cái)嗟淖V范圍還是誤差穩(wěn)定性來看，Chebyshev窗都是一個(gè)較優(yōu)的選擇，但是實(shí)際上，其誤差穩(wěn)定性仍較改進(jìn)二項(xiàng)式窗差.

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

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

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.

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

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

董良國, 馬在田, 曹景忠等. 2000. 一階彈性波方程交錯(cuò)網(wǎng)格高階差分解法. 地球物理學(xué)報(bào), 43(3): 411-419.

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

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

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

(本文編輯何燕)

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

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

doi:10.6038/cjg20160728 中圖分類號(hào)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é)報(bào)，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.