黃建平,婁璐烽,彭煒颋,李慶洋

(1.中國(guó)石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院,山東青島 266580;2.海洋國(guó)家實(shí)驗(yàn)室海洋礦產(chǎn)資源評(píng)價(jià)與探測(cè)技術(shù)功能實(shí)驗(yàn)室,山東青島 266071;3.中國(guó)石化中原油田分公司物探研究院,河南濮陽(yáng) 457001)

由于差分離散,有限差分法在數(shù)值模擬過(guò)程中會(huì)出現(xiàn)明顯數(shù)值頻散[1-2],影響地震波場(chǎng)模擬精度[3-4],差分系數(shù)優(yōu)化是壓制地震波場(chǎng)數(shù)值頻散的有效策略。差分系數(shù)優(yōu)化通過(guò)調(diào)整網(wǎng)格點(diǎn)上的差分系數(shù)重新分配“權(quán)重”,能更為準(zhǔn)確地從當(dāng)前時(shí)刻波場(chǎng)值計(jì)算下一時(shí)刻地震波場(chǎng),在不增加計(jì)算量的情況下降低數(shù)值頻散[5],提高地震波場(chǎng)數(shù)值模擬精度[6-8]。根據(jù)實(shí)現(xiàn)方式的差異,差分系數(shù)優(yōu)化方法分為窗函數(shù)和最優(yōu)化方法,前者通過(guò)尋找合適的窗函數(shù)來(lái)快速求取差分系數(shù),后者則是引入最優(yōu)化策略計(jì)算差分系數(shù)。用于差分系數(shù)優(yōu)化的窗函數(shù)法包括海寧窗、高斯窗、可縮放二項(xiàng)式窗、Chebyshev窗等,它將偽譜法的卷積級(jí)數(shù)截?cái)嗟玫酱八阕觼?lái)代替差分系數(shù)。早期,Zhou[9]和Igel[10]等分別使用海寧窗和高斯窗截?cái)嗟玫酱八阕?以獲得更高精度的模擬波場(chǎng);隨后,Chu和Stoffa[11]提出縮放二項(xiàng)式窗,得到旁瓣衰減迅速的差分系數(shù),提高了計(jì)算穩(wěn)定性。Chebyshev窗亦可單獨(dú)用于窗算子設(shè)計(jì),王之洋等[12]提出一種Chebyshev自褶積組合窗,通過(guò)調(diào)節(jié)紋波率、自褶積次數(shù)、加權(quán)系數(shù)3個(gè)參數(shù),能夠直觀調(diào)節(jié)窗形態(tài),改變差分系數(shù)精度。鑒于不同窗算子的特性差異,有學(xué)者設(shè)計(jì)組合窗算子[13-15],Wang等[16]在縮放二項(xiàng)式窗基礎(chǔ)上提出余弦調(diào)制二項(xiàng)式窗,將海寧窗、Chebyshev窗、縮放二項(xiàng)式窗組合,通過(guò)改變調(diào)制時(shí)間和調(diào)制范圍,最終獲得不同的主瓣旁瓣振幅響應(yīng)。也有學(xué)者通過(guò)優(yōu)化方法來(lái)提升窗算子精度,Ren等[17]采用最小二乘優(yōu)化方法迭代求取更高精度的差分系數(shù)。最優(yōu)化策略的差分系數(shù)優(yōu)化方法主要包括模擬退火法、雷米茲交換算法、最小二乘法等。模擬退火法優(yōu)化的差分系數(shù)擁有更高的頻譜覆蓋范圍和較小的誤差限,朱遂偉等[18]采用模擬退火法將傅里葉有限差分系數(shù)中的兩個(gè)系數(shù)擴(kuò)展到4個(gè),能顯著提升成像效果。相較于模擬退火法,雷米茲交換算法計(jì)算效率更高,只需幾次迭代就能求得高精度差分系數(shù)[19-20]。He等[21]應(yīng)用此方法求解差分系數(shù),在獲得更寬頻帶的同時(shí)降低了計(jì)算成本,特別是對(duì)大尺度模型的誤差累積減少有積極作用。最小二乘法是一種最小二范數(shù)全局優(yōu)化方法,它依據(jù)頻散關(guān)系的誤差形式構(gòu)建目標(biāo)函數(shù)來(lái)求解優(yōu)化的差分系數(shù)[22],與之類(lèi)似的還有最大范數(shù)優(yōu)化[23]、最小范數(shù)優(yōu)化[24]等。其中Liu首先提出基于最小二乘法的差分系數(shù)優(yōu)化方法,該方法能明顯壓制高波數(shù)域頻散,其低階差分系數(shù)可以取代高階傳統(tǒng)差分系數(shù),基于此,Liu[25]和任志明[26]進(jìn)一步發(fā)展了交錯(cuò)網(wǎng)格差分系數(shù),Yong等[27]則發(fā)展了彈性波差分系數(shù)優(yōu)化。雖然最小二乘法是差分系數(shù)優(yōu)化研究的一個(gè)熱點(diǎn),但其低波數(shù)段誤差仍然較大,存在頻散誤差以類(lèi)似于正弦曲線的形態(tài)小范圍偏離理想值的問(wèn)題,無(wú)法準(zhǔn)確逼近0誤差。因此如何在保持最小二乘法寬頻帶的優(yōu)勢(shì)下,減小其低波數(shù)段頻散誤差,成為改進(jìn)最小二乘差分系數(shù)優(yōu)化方法的一個(gè)研究方向。筆者提出一種基于拉格朗日乘子的優(yōu)化新方法,為保證寬頻帶,采用絕對(duì)誤差二范數(shù)構(gòu)造目標(biāo)函數(shù),同時(shí)加入絕對(duì)誤差和作為約束條件,以增加目標(biāo)函數(shù)穩(wěn)定性并減小低波數(shù)段頻散誤差。

1 方法原理

1.1 最小二乘優(yōu)化原理

最小二乘法差分系數(shù)優(yōu)化原理從空間偏導(dǎo)數(shù)的差分格式入手,有限差分法的二階空間偏導(dǎo)數(shù)的2M階精度差分格式可以表示為

(1)

式中,x代表水平方向;h為空間采樣間隔;p為在m網(wǎng)格點(diǎn)的波場(chǎng)值;cm為差分系數(shù);c0為中心差分系數(shù);M為差分系數(shù)長(zhǎng)度的一半或差分階數(shù)的一半。差分系數(shù)可以認(rèn)為是不同網(wǎng)格點(diǎn)的“權(quán)重”,即如何利用當(dāng)前時(shí)刻中心點(diǎn)及其周?chē)W(wǎng)格點(diǎn)的波場(chǎng)值來(lái)計(jì)算下一時(shí)刻中心點(diǎn)波場(chǎng)值。需要指出的是,當(dāng)前時(shí)刻中心點(diǎn)處的波場(chǎng)值對(duì)下一時(shí)刻該中心點(diǎn)的波場(chǎng)值影響最大,中心差分系數(shù)c0也最大,隨著網(wǎng)格點(diǎn)距中心點(diǎn)的距離增加,這些網(wǎng)格點(diǎn)波場(chǎng)值的“權(quán)重”也減小,差分系數(shù)值隨之變小。此外,考慮的范圍趨于無(wú)窮遠(yuǎn),將有無(wú)限個(gè)網(wǎng)格點(diǎn)波場(chǎng)值對(duì)下一時(shí)刻波場(chǎng)產(chǎn)生影響,即cm無(wú)窮多,此時(shí)波場(chǎng)計(jì)算較為準(zhǔn)確,但計(jì)算量巨大。為此,需要尋找新的差分系數(shù),實(shí)現(xiàn)有限個(gè)點(diǎn)對(duì)下一時(shí)刻波場(chǎng)相對(duì)精確估算。

將波場(chǎng)p表示為簡(jiǎn)諧波形式,將波數(shù)引入推導(dǎo)中,波場(chǎng)形式如下:

p(x+mh)=p0exp(ik(x+mh)).

(2)

式中,k為波數(shù);p0為一常數(shù)。將式(2)代入式(1)可得

(3)

為得到僅包含差分系數(shù)的等式,取k=0,式(3)可以表示為

(4)

式(4)反映了不同網(wǎng)格點(diǎn)“權(quán)重”之和存在平衡狀態(tài),差分系數(shù)在一維線網(wǎng)格、二維面網(wǎng)格中皆圍繞中心對(duì)稱(chēng)分布,網(wǎng)格點(diǎn)處波場(chǎng)值的累加對(duì)下一時(shí)刻中心波場(chǎng)值的影響需達(dá)到平衡,否則當(dāng)前時(shí)刻的波場(chǎng)值將引起下一時(shí)刻波場(chǎng)值的異常突變。

將式(4)代入式(3)推導(dǎo)可得

(5)

令式(5)中

β=kh,

(6)

f(β)=β2,

(7)

φm(β)=2(1-cos(mβ)).

(8)

聯(lián)立式(5)～(8)可得

(9)

為使kh在Nyquist頻率上等價(jià)于π,β的取值范圍設(shè)置為0到π,取式(9)的絕對(duì)誤差二范數(shù),即最小二乘法,構(gòu)造出對(duì)變量β進(jìn)行積分的誤差函數(shù)E:

(10)

不同于傳統(tǒng)泰勒展開(kāi)獲取差分系數(shù)的策略,最小二乘法利用誤差最小化,使差分格式逼近偏導(dǎo)數(shù)。積分上限b代表優(yōu)化波數(shù)區(qū)間,b越大,優(yōu)化波數(shù)區(qū)間越大,但低波數(shù)段頻散誤差也相應(yīng)變大,反之優(yōu)化波數(shù)區(qū)間變小,能滿足更小頻散誤差需求。需要注意的是,當(dāng)b取值變小時(shí),計(jì)算結(jié)果可能不穩(wěn)定,對(duì)差分系數(shù)求解有一定影響。

對(duì)式(10)求極小,即E對(duì)cm求導(dǎo)數(shù)值為0,可得

(11)

式(11)展開(kāi)可得

(12)

求解式(12)即可得到最小二乘法優(yōu)化的差分系數(shù)。從式(12)矩陣方程可知,矩陣中的元素皆為積分形式極為復(fù)雜,前人研究證實(shí)僅由該式求取的差分系數(shù)可能出現(xiàn)頻散誤差不穩(wěn)定問(wèn)題。本文中新方法主要通過(guò)增加絕對(duì)誤差和,引入拉格朗日乘子穩(wěn)定性條件,減弱最小二乘法頻散誤差不穩(wěn)定問(wèn)題。

1.2 拉格朗日法

基于拉格朗日法的差分系數(shù)優(yōu)化從空間偏導(dǎo)數(shù)的差分格式入手,與最小二乘法類(lèi)似得到包含差分系數(shù)的等式。在誤差函數(shù)的構(gòu)建上,首先同時(shí)推導(dǎo)出絕對(duì)誤差二范數(shù)及其誤差和的表達(dá)式,隨后構(gòu)造拉格朗日函數(shù),以絕對(duì)誤差二范數(shù)為目標(biāo)函數(shù),保證差分系數(shù)寬頻帶,以絕對(duì)誤差和為約束條件,增加目標(biāo)函數(shù)穩(wěn)定性,減小差分系數(shù)低波數(shù)段頻散誤差,最終得到包含差分系數(shù)和拉格朗日乘子的矩陣方程。該誤差函數(shù)可以在保證絕對(duì)誤差二范數(shù)最小的基礎(chǔ)上通過(guò)可變參數(shù)拉格朗日乘子增加函數(shù)穩(wěn)定性,約束由二范數(shù)引起的低波數(shù)段頻散誤差不穩(wěn)定,兼顧寬頻帶和低頻散誤差的效果。

在相同階數(shù)下,拉格朗日法與最小二乘法得到的矩陣方程矩陣元素的表達(dá)式和矩陣大小不同,拉格朗日法得到的矩陣方程的行數(shù)和列數(shù)總比最小二乘法多一行和多一列,多出的是為了求解拉格朗日乘子。

與式(3)相同推導(dǎo),新方法將平面波的簡(jiǎn)諧波解代入空間導(dǎo)數(shù)的差分格式中,在二維情況下引入波數(shù),得到與式(3)類(lèi)似的公式:

(13)

(14)

式中,kx和kz分別為x、z方向的波數(shù)。在二維波動(dòng)方程中,地震波在xz平面中傳播,其振幅、傳播速度等物理量都可以分解到x、z方向上,同理波數(shù)k也可以分解為kx、kz,分別為總波數(shù)在x、z方向上的投影值,θ為地震波傳播方向與x軸的夾角,可得

kx=kcosθ,

(15)

kz=ksinθ.

(16)

根據(jù)最小二乘法優(yōu)化差分系數(shù)的思想,首先要得到一個(gè)代表誤差的量,將式(6)、(15)和(16)代入式(13)與(14),定義誤差函數(shù)如下:

(17)

(18)

式中,誤差函數(shù)Am(β,θ)、Bm(β,θ)分別為式(13)、(14)的絕對(duì)誤差,當(dāng)二者同時(shí)取極小值時(shí),此時(shí)的差分系數(shù)cm滿足在二維情況下同時(shí)減小x、z方向上的頻散。若絕對(duì)誤差無(wú)限逼近零,則優(yōu)化的差分系數(shù)無(wú)限逼近真實(shí)空間偏導(dǎo)數(shù)。為簡(jiǎn)化形式,引入中間變量:

ψm(β,θ)=2(1-cos(mβcosθ)),

(19)

φm(β,θ)=2(1-cos(mβsinθ)).

(20)

為達(dá)到既保持最小二乘法頻帶寬的優(yōu)勢(shì),又達(dá)到減小低波數(shù)段頻散誤差的目的,采用拉格朗日乘數(shù)法構(gòu)造誤差函數(shù)。誤差函數(shù)φ(β,θ,λ)表達(dá)式如下:

φ(β,θ,λ)=

(21)

式中,λ為拉格朗日乘子,代表絕對(duì)誤差和在誤差函數(shù)中所占權(quán)重,在計(jì)算差分系數(shù)的同時(shí)λ也會(huì)計(jì)算出來(lái)。該誤差函數(shù)包含傳播角度θ,將波數(shù)劃分為x方向和z方向,能夠更好模擬聲波在二維模型中的傳播過(guò)程,并且充分利用兩個(gè)方向的波數(shù)。其中前半部分為x方向的絕對(duì)誤差二范數(shù),繼承了最小二乘法頻帶寬的優(yōu)勢(shì),后半部分為z方向的絕對(duì)誤差和,體現(xiàn)了新方法為約束最小二乘法在低波數(shù)段產(chǎn)生較大頻散誤差所采取的策略。前后部分都采取二重積分的形式,優(yōu)化角度θ取0～2π以保證差分系數(shù)能滿足所有傳播方向,積分上限b取0～π,b取值越大,優(yōu)化波數(shù)區(qū)間越廣。

為使誤差函數(shù)取極小,令誤差函數(shù)φ(β,θ,λ)對(duì)差分系數(shù)cm和拉格朗日乘子λ的偏導(dǎo)數(shù)都為0,同時(shí)滿足最小誤差條件,條件表示如下:

(22)

(23)

化簡(jiǎn)可得

(24)

(25)

式(24)、(25)分別為滿足式(22)、(23)的矩陣,式(24)包含M個(gè)等式,且同時(shí)含有M個(gè)差分系數(shù)以及拉格朗日乘子λ,式(25)僅包含M個(gè)差分系數(shù)。

將式(24)、(25)聯(lián)立展開(kāi)可得

(26)

從矩陣方程中可以看到,需要求解的未知量為差分系數(shù)cm以及拉格朗日乘子λ共M+1個(gè)變量,方程最左端為(M+1)×(M+1)的方陣,右端為(M+1)×1的長(zhǎng)方陣。求解過(guò)程發(fā)現(xiàn),矩陣方程中用二重積分表達(dá)的元素值對(duì)矩陣方程求解的穩(wěn)定性有重要影響。計(jì)算的差分階數(shù)越高,矩陣方程越大,對(duì)積分精確計(jì)算的要求也越高。

表1為最小二乘法和拉格朗日法利用最小化空間域頻散關(guān)系絕對(duì)誤差優(yōu)化差分系數(shù)的b取值，b代表優(yōu)化波數(shù)區(qū)間，誤差限η為1×10-4。通過(guò)限制誤差限，本文中對(duì)比了兩種方法在相同誤差限下的b取值，表1中最小二乘法記為L(zhǎng)SM，拉格朗日法記為L(zhǎng)agrange?？梢钥吹?，除M=3，其他階數(shù)下拉格朗日法的b取值普遍比最小二乘法小，其原因可能為拉格朗日法相對(duì)于最小二乘法在低波數(shù)段的頻散誤差較小，但同時(shí)其有效頻帶寬度變小，頻散誤差最大值變大，為保證符合誤差限要求，必須降低b值。在表1的b取值下，利用拉格朗日法計(jì)算得到不同階數(shù)的優(yōu)化差分系數(shù)，如表2所示。當(dāng)M從2增加至10時(shí)，相應(yīng)的拉格朗日乘子λ分別為-3.694×10-5、2.912×10-5、-1.221×10-5、1.131×10-5、-9.017×10-6、1.765×10-6、-6.347×10-6、5.175×10-6、-1.563×10-6。

表1 最小二乘法和拉格朗日法的b取值Table 1 Value of b in the least square method and Lagrange method

表2 拉格朗日法的優(yōu)化差分系數(shù)Table 2 Optimal finite-difference coefficient of Lagrange method

2 頻散測(cè)試

在正演模擬中不同地震子波具有不同主頻,實(shí)際地震波更加復(fù)雜,往往覆蓋多種頻率,從幾赫茲到幾十赫茲不等[28],且炸藥震源每次激發(fā)時(shí)產(chǎn)生的波形不盡相同,因此差分系數(shù)的有效頻帶寬度需較寬,即使是固定波形的激發(fā)震源,也期望差分系數(shù)在某一頻段內(nèi)頻散誤差小。由此,差分系數(shù)優(yōu)化須同時(shí)追求寬頻帶和低頻散誤差,本文的目標(biāo)就是在保持最小二乘法寬頻帶的基礎(chǔ)上,進(jìn)一步減小其低波數(shù)段頻散誤差。

采用絕對(duì)誤差ε進(jìn)行頻散分析,在頻散曲線中,ε等于0時(shí),數(shù)值模擬無(wú)空間頻散,ε偏離0越大,頻散誤差越大,在正演模擬中空間頻散越明顯。

為對(duì)比最小二乘法和新方法的頻散誤差大小,本文中以kh為橫軸,絕對(duì)誤差ε為縱軸,作圖1所示最小二乘法與新方法優(yōu)化的差分系數(shù)在不同階數(shù)下的絕對(duì)誤差頻散分析圖。誤差限皆為1×10-4,重點(diǎn)關(guān)注低波數(shù)段。圖1中,虛線表示誤差限,以限制最大頻散誤差值,并調(diào)整參數(shù)b使兩條頻散曲線的最大頻散誤差一致,便于比較二者在相同誤差限下的頻散誤差大小。通過(guò)圖1可以看到,在低波數(shù)段,粉色曲線明顯更靠近零值,曲線更加平緩,說(shuō)明新方法頻散誤差比最小二乘法小,隨著波數(shù)變大,兩條曲線漸趨一致,達(dá)到誤差限之后快速下降,頻散誤差急速增大,新方法比最小二乘法更快超出誤差限,說(shuō)明新方法的有效頻帶寬度較小。

綜上可以看出:①在相同誤差限下,新方法優(yōu)化的差分系數(shù)在有效頻帶寬度內(nèi)的頻散誤差小于最小二乘法;②新方法對(duì)低階和高階差分系數(shù)都有優(yōu)化作用,都能減小頻散誤差,在階數(shù)較低時(shí),如M=3、5,優(yōu)化效果好,階數(shù)較高時(shí),如M=8,優(yōu)化效果減弱;③在低波數(shù)段,新方法的優(yōu)化效果更加明顯,減小頻散誤差效果最好,隨著波數(shù)變大,效果減弱;④新方法的有效頻帶寬度較小,但與最小二乘法差距很小。頻散分析證明,本文方法不但繼承了最小二乘法寬頻帶的優(yōu)點(diǎn),同時(shí)減小了低波數(shù)段頻散誤差,在低階、低波數(shù)段的頻散壓制效果好。

圖1 最小二乘法與新方法的頻散分析圖Fig.1 Dispersion analysis of least square method and new method

由于最小二乘法缺少限制條件約束,也沒(méi)有使用一范數(shù)模型作為目標(biāo)函數(shù),導(dǎo)致其在有效頻帶寬度展寬的同時(shí)低波數(shù)段頻散誤差增加。新方法能夠改善這一問(wèn)題,其原因在于新方法采用拉格朗日乘數(shù)法構(gòu)造誤差函數(shù),該誤差函數(shù)包含絕對(duì)誤差二范數(shù),能夠擴(kuò)寬有效頻帶,且增加了絕對(duì)誤差和,能增強(qiáng)誤差函數(shù)穩(wěn)定性,在頻散曲線中體現(xiàn)為低波數(shù)段頻散誤差的減小。如Miao等[24]提出的一種基于一范數(shù)構(gòu)造目標(biāo)函數(shù)的差分系數(shù)優(yōu)化方法,能減小誤差累積,增加方法穩(wěn)定性。由于本文中將兩種方法的誤差限調(diào)整一致,因此新方法在低波數(shù)段頻散誤差減小明顯,隨后兩種方法漸趨一致。差分階數(shù)較低時(shí),新方法減小低波數(shù)段頻散誤差的效果明顯(圖1(a)、(b))。隨著階數(shù)升高,最小二乘法精度提升,因此新方法在高階時(shí)減小頻散誤差的效果減弱(圖1(d))。新方法也為減小頻散誤差犧牲了部分有效頻帶寬度。

圖2為傳統(tǒng)方法與新方法優(yōu)化的差分系數(shù)在不同差分階數(shù)下的絕對(duì)誤差頻散分析。在差分階數(shù)選取上,分別選取10階新方法,10、16、24階傳統(tǒng)方法進(jìn)行對(duì)比。首先觀察4條曲線在全波數(shù)域的頻散誤差,可以發(fā)現(xiàn)kh在1.5～3之間,傳統(tǒng)方法的絕對(duì)誤差頻散曲線隨波數(shù)增大而增大,圖中黑色曲線整體高于紅色曲線,新方法優(yōu)于同階數(shù)傳統(tǒng)方法,且與16階傳統(tǒng)方法相當(dāng)。為觀察頻散曲線前部的微小變化,重點(diǎn)放大kh在0～2之間的區(qū)域,且只顯示頻散誤差在±2×10-4之間的曲線?？梢钥吹?新方法有效頻帶更寬,具體表現(xiàn)為新方法有效頻帶寬度比傳統(tǒng)方法10、16階廣,可與24階比擬,同時(shí),新方法的低波數(shù)段頻散誤差大于傳統(tǒng)方法,頻散曲線在給定的誤差范圍內(nèi)振蕩較快,以類(lèi)似于正弦函數(shù)的形態(tài)小范圍偏離理想值,無(wú)法準(zhǔn)確達(dá)到0誤差,傳統(tǒng)方法則表現(xiàn)為一段0直線。

分析說(shuō)明:①新方法具有寬頻帶優(yōu)勢(shì),能極大擴(kuò)寬有效頻帶,低階新方法可以達(dá)到高階傳統(tǒng)方法的有效頻帶寬度;②在低波數(shù)段,新方法頻散誤差大于傳統(tǒng)方法,可能是因?yàn)閮?yōu)化的差分系數(shù)具有更寬的有效頻帶寬度,但代價(jià)是在有效頻帶寬度內(nèi)分布有更大的誤差,通過(guò)大量模型和參數(shù)測(cè)試,這些誤差都在可以接受的范圍。該結(jié)果同樣由新方法誤差函數(shù)的構(gòu)造導(dǎo)致,新方法的目標(biāo)函數(shù)為絕對(duì)誤差二范數(shù),其獲得的最小二乘解不是零誤差,而是在“二范數(shù)最小”意義下達(dá)到的殘差最小。因此與最小二乘法類(lèi)似,能夠極大展寬有效頻帶,同時(shí)低波數(shù)段頻散誤差較大。

在測(cè)試過(guò)程中,隨著階數(shù)升高,本文方法也存在一定不穩(wěn)定性,隨著參數(shù)b減小,新方法的最大頻散誤差、低波數(shù)段頻散誤差、有效頻帶寬度同步減小,但發(fā)現(xiàn)部分異?，F(xiàn)象:①在某些較小的b值上,新方法符合誤差限,同時(shí)低波數(shù)段頻散誤差減小非常明顯,如M=5,b=1.56時(shí),M=6,b=1.87時(shí);②b小于一定值時(shí)頻散曲線不穩(wěn)定,b的微小變化可能引起曲線劇烈變化,如M=5,b<1.55時(shí),M=6,b<1.70時(shí),M=7,b<2.02時(shí),M=8,b,2.25時(shí);③在某些b值上,新方法頻散誤差突然變大,如M=6,b=2.04時(shí),M=7,b=2.15時(shí),M=8,b=2.265時(shí),M=9,b=2.39時(shí);④差分階數(shù)較大時(shí),隨著b減小,最大頻散誤差不減小,如M=7,b>2.10時(shí)。總體而言,新方法在低階時(shí)很穩(wěn)定,隨著階數(shù)升高,逐漸不穩(wěn)定,頻散曲線更容易隨b減小發(fā)生突變,如M>6時(shí)。究其原因,可能是隨著矩陣階數(shù)升高,矩陣穩(wěn)定性變差,導(dǎo)致求解變得困難,且矩陣元素皆為積分形式,積分上限b的微小變化能夠影響整個(gè)矩陣,特別是b值較小時(shí)高階系數(shù)的求解。當(dāng)優(yōu)化波數(shù)范圍較小時(shí)可能存在不穩(wěn)定,可將正則化算子[29]引入系數(shù)求取過(guò)程,雍鵬等[30]采用牛頓法快速求解優(yōu)化差分系數(shù)。

圖2 傳統(tǒng)方法與新方法的頻散分析Fig.2 Dispersion analysis of traditional method and new method

3 模型試算

3.1 均勻模型

考慮到研究重點(diǎn)為低波數(shù)段,取模型網(wǎng)格點(diǎn)數(shù)為200×200,模型速度為1 500 m/s,震源位于模型中心,計(jì)算主頻為30 Hz的雷克子波,空間采樣間隔5 m,時(shí)間采樣間隔0.1 ms,以保證數(shù)值頻散以空間頻散為主,時(shí)間頻散幾乎不出現(xiàn)。

為便于對(duì)比3種方法數(shù)值模擬的差異,采用四宮格形式展示波場(chǎng)快照計(jì)算結(jié)果,每個(gè)宮格對(duì)應(yīng)一種情況,同時(shí)對(duì)應(yīng)給出波場(chǎng)殘差圖以比較不同方法的模擬精度和穩(wěn)定性。圖3為傳統(tǒng)方法、最小二乘法、新方法分別在傳播時(shí)間t=0.3,1.6 s時(shí)的均勻模型波場(chǎng)快照以及波場(chǎng)快照殘差對(duì)比,其中圖3(a)、(c)為波場(chǎng)快照對(duì)比,圖3(b)、(d)為波場(chǎng)快照殘差對(duì)比,參考波場(chǎng)采用50階傳統(tǒng)方法(可視為精確解)。沿逆時(shí)針?lè)较?左上角編號(hào)1采用8階傳統(tǒng)方法,左下角編號(hào)2采用8階最小二乘法,右下角編號(hào)3采用8階新方法,右上角編號(hào)4采用12階傳統(tǒng)方法,由此,3種方法可依次循環(huán)對(duì)比。

對(duì)比圖3(b)發(fā)現(xiàn),低階傳統(tǒng)方法的殘差表現(xiàn)為黑白交替的深色條帶狀,殘差值大,最小二乘法的殘差表現(xiàn)為邊緣模糊的帶狀,中部殘差值最大,且在部分區(qū)域殘差覆蓋范圍擴(kuò)大,新方法則能改進(jìn)這一問(wèn)題,殘差值明顯減小,殘差覆蓋范圍縮小,最接近于高階傳統(tǒng)方法模擬精度,在圖3(d)箭頭與紅框所示處同樣能觀察到頻散依次改善的情況。計(jì)算結(jié)果表明在相同階數(shù)下,新方法的空間頻散比傳統(tǒng)方法和最小二乘法小,波場(chǎng)快照殘差更小,模擬精度更高,低階新方法可以替代高階傳統(tǒng)方法的模擬精度。這與頻散分析中新方法寬頻帶且低波數(shù)段頻散誤差小的結(jié)論相符,新方法有效頻帶寬度比傳統(tǒng)方法寬,因此殘差小于傳統(tǒng)方法,低波數(shù)段頻散誤差比最小二乘法小,因此殘差小于最小二乘法,能量更加集中于反射波主波形上,減小能量分散。

為精細(xì)對(duì)比波場(chǎng)快照殘差圖,抽取單道進(jìn)行分析。圖4為3種方法不同階數(shù)在不同傳播時(shí)刻的單道殘差對(duì)比分析,圖4(a)、(b)分別來(lái)自圖3(b)、(d)500 m處的單道記錄。由下到上依次繪制8階傳統(tǒng)方法、最小二乘法、新方法和12階傳統(tǒng)方法,并繪制參考道單道殘差(值為零),用以量化3種方法的殘差值。在圖中相同深度處選點(diǎn)測(cè)量殘差值,選段測(cè)量殘差覆蓋范圍,進(jìn)行數(shù)值標(biāo)注。

圖4(a)中,低階傳統(tǒng)方法的殘差曲線變化幅度最劇烈,呈明顯的鋸齒形,而同階最小二乘法和新方法的殘差曲線變化幅度小,較平滑,且新方法殘差值明顯比最小二乘法小。在圖4(b)中,同階數(shù)時(shí)新方法殘差曲線最平緩、殘差值最小,特別是在中部深度520 m處,優(yōu)化效果顯著,而深度750 m處新方法殘差值最小,即頻散最小。從殘差覆蓋范圍的測(cè)量中也可以證實(shí)新方法能夠減小殘差覆蓋范圍。

簡(jiǎn)單模型測(cè)試結(jié)果表明,新方法地震波場(chǎng)快照殘差小于最小二乘法,同時(shí)約束了殘差覆蓋范圍。在同階數(shù)下,新方法的頻散小于傳統(tǒng)方法和最小二乘法,證明了新方法能夠在明顯減小傳統(tǒng)方法頻散的優(yōu)勢(shì)上進(jìn)一步降低最小二乘法頻散的結(jié)論,能夠?qū)崿F(xiàn)更高的數(shù)值模擬精度。

圖3 不同方法的波場(chǎng)快照和殘差對(duì)比Fig.3 Comparison of snapshots and residuals of snapshots of different methods

3.2 Marmousi模型

為驗(yàn)證新方法對(duì)復(fù)雜模型的適應(yīng)性,采用SEG提供的國(guó)際Marmousi模型進(jìn)行正演模擬,模型網(wǎng)格點(diǎn)數(shù)為1 360×700,震源位于表面中心,計(jì)算主頻為30 Hz的雷克子波,空間采樣間隔5 m,時(shí)間采樣間隔0.1 ms,地震波場(chǎng)模擬時(shí)間2.0 s。

圖5為最小二乘法、新方法在傳播時(shí)間t=2.0 s時(shí)的Marmousi模型波場(chǎng)快照和波場(chǎng)快照殘差對(duì)比。其中圖5(a)為波場(chǎng)快照對(duì)比,圖5(b)為波場(chǎng)快照殘差對(duì)比,參考波場(chǎng)采用50階傳統(tǒng)方法。為使對(duì)比直觀,采取左右對(duì)比的方式,左、右上方分別為8階最小二乘法、新方法,左、右下方分別為16階最小二乘法、新方法。

圖5(a)中模擬波場(chǎng)非常復(fù)雜,外圍波前面能量最強(qiáng)但無(wú)規(guī)則彎曲,內(nèi)部存在眾多反射波形,波形混疊難以分辨,左下角存在明顯淺色區(qū)域,可能是模型構(gòu)造復(fù)雜、存在速度異常體導(dǎo)致。通過(guò)圖5(a)無(wú)法分辨出兩種方法的頻散大小,因此重點(diǎn)對(duì)比波場(chǎng)快照殘差,可以看出,在相同階數(shù)下新方法波場(chǎng)快照殘差小于最小二乘法,在箭頭所指波前面處和紅框所示內(nèi)部區(qū)域處殘差明顯減小。說(shuō)明低階新方法優(yōu)化效果好,8階新方法的模擬精度高于16階最小二乘法,但高階新方法優(yōu)化效果減弱,應(yīng)用效果不顯著。Marmousi模型正演表明,新方法在復(fù)雜模型中同樣能夠減小頻散,將地震波能量集中于波前面上,提升模擬精度。

為更加精細(xì)對(duì)比兩種方法的波場(chǎng)快照情況,抽取單道進(jìn)行分析。圖6為兩種方法不同階數(shù)的單道殘差對(duì)比分析,其中,圖6(a)、(b)分別為圖5(b)中3 000、5 000 m處的單道殘差。將4條殘差曲線按階數(shù)分為兩組進(jìn)行對(duì)比,下方藍(lán)色組、上方紅色組分別為8階和16階最小二乘法和新方法。按組對(duì)比可見(jiàn),新方法殘差曲線總體上較最小二乘法平緩,尤其是藍(lán)色低階方法,具體表現(xiàn)為波峰波谷殘差值明顯減小,曲線整體更靠近零值,從標(biāo)記的數(shù)值點(diǎn)上可見(jiàn),藍(lán)色低階時(shí)殘差值減小60%以上,紅色高階時(shí)殘差值減小40%以上,新方法減小數(shù)值頻散的效果優(yōu)于最小二乘法。綜上所述,新方法對(duì)于復(fù)雜模型也具有較好適應(yīng)性,在低階、高階時(shí)都能提高數(shù)值模擬精度,低階時(shí)改善效果明顯,隨著階數(shù)升高改善效果減弱。Marmousi模型單道殘差分析同樣證明新方法壓制數(shù)值頻散的效果優(yōu)于最小二乘法。

圖4 不同方法的單道殘差對(duì)比Fig.4 Comparison of single-channel residuals of different methods

圖5 Marmousi模型的波場(chǎng)快照和殘差對(duì)比Fig.5 Comparison of snapshots and residuals of snapshots of Marmousi model

圖6 兩種方法的單道殘差對(duì)比Fig.6 Comparison of single-channel residuals of two methods

4 結(jié)束語(yǔ)

針對(duì)最小二乘差分系數(shù)優(yōu)化法在低波數(shù)段頻散誤差大的問(wèn)題,在采用絕對(duì)誤差二范數(shù)的同時(shí),有效利用誤差和約束頻散誤差,構(gòu)造拉格朗日函數(shù)求解差分系數(shù),提出了一種基于拉格朗日乘子的空間域差分系數(shù)優(yōu)化新方法。均勻模型和Marmousi模型測(cè)試結(jié)果表明:新方法的低波數(shù)段頻散誤差小于最小二乘法,低階、低波數(shù)時(shí)優(yōu)化效果明顯,階數(shù)升高效果減弱,且相比傳統(tǒng)方法,新方法的有效頻帶更寬;新方法能夠減小波場(chǎng)殘差同時(shí)縮小殘差覆蓋范圍,同時(shí)其數(shù)值頻散也優(yōu)于最小二乘法,能在一定程度上提高地震波場(chǎng)模擬精度;新方法對(duì)復(fù)雜模型適應(yīng)性較好,能在壓制數(shù)值頻散的同時(shí)集中波前面能量,減少地震波能量分散從而提升模擬精度。