彭煒颋,黃建平,2

(1.中國(guó)石油大學(xué)(華東),山東青島266580;2.海洋國(guó)家試驗(yàn)室海洋礦產(chǎn)資源評(píng)價(jià)與探測(cè)技術(shù)功能實(shí)驗(yàn)室,山東青島266071)

有限差分法[1-2]是一種經(jīng)典的數(shù)值模擬方法,能夠兼顧計(jì)算效率與模擬精度,目前已經(jīng)被廣泛應(yīng)用于勘探地震波場(chǎng)的數(shù)值模擬中。有限差分方法是用差分算子逼近空間偏導(dǎo)和時(shí)間偏導(dǎo)的方法,該方法可以實(shí)現(xiàn)離散模型的數(shù)值模擬,但是如果在對(duì)模型離散的過(guò)程中使用粗網(wǎng)格,模擬結(jié)果就會(huì)產(chǎn)生頻散誤差[3-4]。在有限差分方法中,常規(guī)的有限差分系數(shù)由泰勒級(jí)數(shù)展開推導(dǎo)而得[5],該差分系數(shù)在低波數(shù)段能有效且準(zhǔn)確地模擬地震波場(chǎng),但是在高波數(shù)段,模擬結(jié)果會(huì)出現(xiàn)嚴(yán)重的頻散誤差[6]。隨著地震技術(shù)的不斷發(fā)展,對(duì)地震波場(chǎng)正演模擬精度的要求不斷提高,優(yōu)化有限差分系數(shù)可以在不增加計(jì)算量,不改變差分格式的前提下,提高模擬的精度和計(jì)算效率。所以如何通過(guò)優(yōu)化差分系數(shù)來(lái)減小頻散誤差是有限差分方法研究領(lǐng)域一個(gè)重要的研究方向。

在優(yōu)化差分系數(shù)的研究中,用優(yōu)化算法最小化時(shí)空域頻散關(guān)系或者波數(shù)域頻散關(guān)系的誤差來(lái)獲得優(yōu)化差分系數(shù)是一種常用策略。雍鵬等[7]和鄒強(qiáng)等[8]在時(shí)空域中,通過(guò)最小化給定波數(shù)范圍內(nèi)的頻散誤差來(lái)計(jì)算優(yōu)化差分算子。雍鵬等[9]在時(shí)空域中實(shí)現(xiàn)了時(shí)空差分算子的同步優(yōu)化,并且采用共軛梯度法增加求解過(guò)程的穩(wěn)定性。LIU[10-11]提出用最小二乘法來(lái)計(jì)算優(yōu)化差分系數(shù),該算法的目標(biāo)函數(shù)是波數(shù)域頻散關(guān)系誤差或時(shí)空域頻散關(guān)系誤差。該方法通過(guò)求解一定波數(shù)范圍內(nèi)的非迭代解,實(shí)現(xiàn)了對(duì)頻散誤差的全局優(yōu)化。最小二乘方法雖然在低波數(shù)段產(chǎn)生了頻散誤差,卻有效拓寬了有限差分算子的有效帶寬,使得數(shù)值模擬中的高波數(shù)段頻散顯著降低。MIAO等[12]用交替方向乘子法解決了基于一范數(shù)目標(biāo)函數(shù)的差分系數(shù)優(yōu)化問題,且通過(guò)數(shù)值實(shí)驗(yàn)和理論分析驗(yàn)證了相較于二范數(shù)和無(wú)窮范數(shù)函數(shù),一范數(shù)目標(biāo)函數(shù)會(huì)獲得更小的低波數(shù)段頻散誤差。ZHANG等[13-14]采用模擬退火算法最小化波數(shù)域頻散關(guān)系絕對(duì)誤差來(lái)計(jì)算優(yōu)化差分系數(shù)。用模擬退火算法得到的優(yōu)化差分系數(shù)比用最小二乘方法得到的優(yōu)化差分系數(shù)具有更寬的有效帶寬,但是由于該方法在迭代過(guò)程中的不穩(wěn)定性,很難用該方法計(jì)算高階的優(yōu)化差分系數(shù)。YANG等[15]應(yīng)用抽樣方法來(lái)近似頻散關(guān)系,該方法的優(yōu)化差分系數(shù)在數(shù)值模擬中有著良好的表現(xiàn)。該抽樣方法類似無(wú)迭代的雷米茲交換算法。AN[16]用類雷米茲交換算法來(lái)計(jì)算真實(shí)頻散關(guān)系的統(tǒng)一近似,但是沒有根據(jù)誤差限來(lái)調(diào)整有效帶寬,導(dǎo)致該方法在數(shù)值模擬過(guò)程中產(chǎn)生較大的數(shù)值誤差。ERIK等[17]提出了利用振幅等紋波特性和最小二乘法利用計(jì)算任意樣本位置的任意階導(dǎo)數(shù)的優(yōu)化差分系數(shù)方法。該方法基于復(fù)值雷米茲交換算法,將復(fù)值雷米茲交換算法應(yīng)用于3個(gè)目標(biāo)函數(shù)(總誤差、相對(duì)誤差和群速度誤差)來(lái)獲得優(yōu)化差分系數(shù)。YANG等[18]利用雷米茲交換算法約束交錯(cuò)網(wǎng)格中的頻散關(guān)系,得到了一種全新的差分系數(shù),該差分系數(shù)能顯著拓寬模擬帶寬。但是,由于YANG等[18]對(duì)零波數(shù)處的嚴(yán)格約束,導(dǎo)致頻散曲線不具有振幅誤差等波紋的特性。

HE等[19]通過(guò)修改近零波數(shù)條件,用雷米茲交換算法約束波數(shù)域頻散關(guān)系的絕對(duì)誤差。該優(yōu)化差分系數(shù)的頻散曲線具有振幅誤差等波紋的特性,使得該優(yōu)化差分系數(shù)在數(shù)值模擬過(guò)程中具有最寬的有效帶寬和最大的低波數(shù)段誤差。為了兼顧雷米茲交換算法較寬的帶寬特性并且減少低波數(shù)段誤差,本文將拉格朗日乘數(shù)法[20-21]引入到差分系數(shù)優(yōu)化過(guò)程中,用拉格朗日乘數(shù)法求解優(yōu)化問題,得到新的優(yōu)化差分系數(shù),其中優(yōu)化問題由二范數(shù)目標(biāo)函數(shù)和由雷米茲交換算法導(dǎo)出的約束條件組成;再通過(guò)頻散測(cè)試分析該方法保留較寬的帶寬特性和降低低波數(shù)段頻散誤差的效果;最后用均勻模型和改進(jìn)的Marmousi模型的模型數(shù)據(jù)驗(yàn)證本文方法的有效性。

1 雷米茲交換法

首先考慮波場(chǎng)空間二階導(dǎo)數(shù)的差分形式,即:

(1)

式中:u為標(biāo)量波場(chǎng);h代表空間網(wǎng)格大小;am是有限差分方法的差分系數(shù);M代表差分階數(shù)的一半。

再根據(jù)平面波理論,可將波場(chǎng)表示為:

(2)

式中:k為波數(shù);τ為時(shí)間采樣間隔;ω為角頻率。

將公式(2)代入公式(1),并用歐拉公式化簡(jiǎn),可得到空間二階導(dǎo)數(shù)波數(shù)域頻散關(guān)系的絕對(duì)誤差近似式:

(3)

式中:β=kh。在奈奎斯特采樣定理下,β∈[0,π]。公式(3)代表了頻散誤差與波數(shù)之間的關(guān)系,可以由此計(jì)算差分系數(shù)的波數(shù)域頻散關(guān)系曲線。

應(yīng)用雷米茲交換算法之前,需要確定一個(gè)準(zhǔn)確的近零波數(shù)條件。該近零波數(shù)條件[19]為:

(4)

雷米茲交換法是一種尋找連續(xù)函數(shù)的最佳逼近多項(xiàng)式的方法,而最佳逼近多項(xiàng)式和連續(xù)函數(shù)之間的差值會(huì)在n+1個(gè)點(diǎn)上正負(fù)交替變化且絕對(duì)值相等[19],這種關(guān)系可以表示為:

(f-p)x1=-(f-p)x2=…=

(-1)n+1(f-p)xn+1

(5)

根據(jù)公式(3)可得到波數(shù)域頻散關(guān)系的絕對(duì)誤差公式,即:

(6)

結(jié)合公式(5)和公式(6),雷米茲交換算法計(jì)算優(yōu)化差分系數(shù)的線性方程組[19]為:

E(βi)=(-1)iAi=1,2,…,M+1

(7)

式中:A是變量;βi(i=1,2,…,M+1)是在[0,khmax_R]中的采樣點(diǎn),khmax_R是有效波數(shù)的最大值。

整理公式(7),可得計(jì)算優(yōu)化差分系數(shù)的線性方程組:

(8)

采用雷米茲交換算法求解方程(8)的過(guò)程分為以下兩步。

1) 將采樣點(diǎn)βi代入方程(8)中,并用高斯消元法求解差分系數(shù)。求解差分系數(shù)am時(shí),同時(shí)求出A。當(dāng)求得的A比誤差限小時(shí),khmax_R增加(誤差限是人為設(shè)置的,即在有效帶寬內(nèi)頻散誤差曲線的最大值,khmax_R的值就是有效帶寬的大小);當(dāng)求得的A比誤差限大時(shí),khmax_R減小。

2) 將步驟1)得到的差分系數(shù)代入公式(6),計(jì)算頻散關(guān)系的絕對(duì)誤差曲線。在頻散關(guān)系的絕對(duì)誤差曲線上取每個(gè)區(qū)間內(nèi)(由零點(diǎn)劃分的區(qū)間)最大絕對(duì)值點(diǎn)的波數(shù)值作為下一次迭代的采樣點(diǎn)。在頻散關(guān)系的絕對(duì)誤差函數(shù)(公式6)中,方程E(β)=0有M+1個(gè)零點(diǎn)(包括β=0)。這些零點(diǎn)將絕對(duì)誤差曲線分成M個(gè)區(qū)間,每個(gè)區(qū)間內(nèi)最大絕對(duì)值點(diǎn)的波數(shù)值作為下一次迭代的采樣點(diǎn),但是這樣只有M個(gè)新采樣點(diǎn),而下一次迭代需要M+1個(gè)新采樣點(diǎn),為了滿足下一次迭代的條件,將第M+1個(gè)采樣點(diǎn)設(shè)置為β=khmax_R。在第一次迭代中,采樣點(diǎn)通過(guò)在波數(shù)范圍[0,khmax_R]上均勻采樣得到。

重復(fù)步驟1)和步驟2)直到A和誤差限之間的差值小于一個(gè)足夠小的值。當(dāng)循環(huán)終止時(shí),通過(guò)對(duì)優(yōu)化差分系數(shù)計(jì)算絕對(duì)誤差曲線(E(β))得到在該曲線上的M+1個(gè)零點(diǎn)。

2 拉格朗日乘數(shù)法

空間二階導(dǎo)數(shù)波數(shù)域頻散關(guān)系的絕對(duì)誤差可從公式(3)變換為如下形式:

(9)

將公式(9)改寫為:

(10)

并將公式(10)在波數(shù)區(qū)間[0,khmax_L]上積分,得到目標(biāo)函數(shù)如下:

(11)

其中,khmax_L與雷米茲交換算法迭代終止時(shí)的khmax_R相同。

如前文所述,當(dāng)?shù)h(huán)終止時(shí),曲線E(β)具有M+1個(gè)零點(diǎn),用前M個(gè)零點(diǎn)構(gòu)建約束條件。先將前M個(gè)零點(diǎn)代入公式(3)中并化簡(jiǎn)可得:

(12)

且將公式(12)累加組成約束條件。約束條件形式如下:

(13)

目標(biāo)函數(shù)和約束條件都是基于波數(shù)域頻散關(guān)系的絕對(duì)誤差,兩者必須基于相同的頻散誤差形式,這樣可以保持約束條件和目標(biāo)函數(shù)的一致性。

之后,用拉格朗日乘數(shù)法將目標(biāo)函數(shù)和約束條件組合:

ψ(am)=φ(am)+λφ(am)

(14)

其中,λ是拉格朗日算子。將公式(14)對(duì)am(m=1,2,…,M)和λ求導(dǎo),可得如下表達(dá)式:

(15)

將公式(15)展開:

(16)

在波數(shù)區(qū)間[0,khmax_L]上,頻散關(guān)系絕對(duì)誤差E的最大值為:

(17)

表1 基于本文方法的二階導(dǎo)數(shù)的優(yōu)化差分系數(shù)(誤差限為0.0001)

3 理論誤差分析

為對(duì)比本文方法的優(yōu)缺點(diǎn),定義波數(shù)域頻散誤差為:

(18)

有效帶寬是優(yōu)化的波數(shù)范圍的大小,本文方法的有效帶寬是公式(11)中的khmax_L。有效帶寬也和頻散曲線上最后一個(gè)零點(diǎn)的波數(shù)大小成正比,頻散曲線上最后一個(gè)零點(diǎn)的波數(shù)值越大,有效帶寬越寬。

如圖1所示,在誤差限(0.001)相同的條件下,由本文方法的優(yōu)化差分系數(shù)有效帶寬小于由雷米茲交換的優(yōu)化差分系數(shù)有效帶寬,但是明顯寬于最小二乘法的優(yōu)化差分系數(shù)有效帶寬。圖2是在相同帶寬下本文方法與最小二乘方法在不同差分階數(shù)時(shí)的對(duì)比結(jié)果,可以看出,在相同的有效帶寬下,本文方法在低波數(shù)段的頻散誤差比最小二乘法更大。

圖1 3種優(yōu)化方法16階(a)和12階(b)差分系數(shù)的波數(shù)域頻散曲線對(duì)比結(jié)果

圖2 在相同帶寬下本文方法和最小二乘法的頻散曲線對(duì)比結(jié)果a 24階優(yōu)化差分系數(shù)頻散曲線; b 20階優(yōu)化差分系數(shù)頻散曲線; c 16階優(yōu)化差分系數(shù)頻散曲線; d 12階優(yōu)化差分系數(shù)頻散曲線

圖3為在相同誤差限下,模擬退火法、雷米茲交換算法和本文方法在12階(M=6)和16階(M=8)時(shí)的頻散曲線對(duì)比結(jié)果。當(dāng)M=6和M=8時(shí),黑線(本文方法)和紅線(模擬退火方法[13])最后一個(gè)零點(diǎn)的波數(shù)值大小幾乎一致,代表這兩種方法得到的差分系數(shù)有效帶寬近乎一致,藍(lán)線(雷米茲交換算法[19])的最后一個(gè)零點(diǎn)的波數(shù)值大于其它兩種方法,這是因?yàn)槔酌灼澖粨Q算法的誤差振幅符合等波紋條件,根據(jù)切比雪夫準(zhǔn)則,該算法具有最大的有效帶寬。模擬退火法的優(yōu)化差分系數(shù)具有較寬的有效帶寬[13],在該試驗(yàn)中(如圖3),兩種方法有效帶寬幾乎相同的現(xiàn)象表明了本文方法具有寬帶寬的特性。模擬退火算法由于其不穩(wěn)定性導(dǎo)致其無(wú)法求解高階差分系數(shù),而本文方法對(duì)于高階差分系數(shù)的穩(wěn)定求解使其在發(fā)展和應(yīng)用前景方面優(yōu)于模擬退火法。

圖3 相同誤差限下模擬退火法、雷米茲交換算法和本文方法的頻散曲線對(duì)比結(jié)果a 12階優(yōu)化差分系數(shù)頻散曲線; b 16階優(yōu)化差分系數(shù)頻散曲線

4 數(shù)值模擬

4.1 均勻模型

如圖4所示,在均勻模型數(shù)值模擬試驗(yàn)中,子波的主要能量集中分布在某一波數(shù)段內(nèi),而子波能量的分布又受主頻和網(wǎng)格間距的影響。一般而言,網(wǎng)格間距越小,子波能量越趨向于低波數(shù)段。另一方面,通過(guò)優(yōu)化方法得到的優(yōu)化差分系數(shù)在低波數(shù)段都會(huì)產(chǎn)生頻散,而用泰勒展開方法得到的常規(guī)差分系數(shù)在低波數(shù)段的頻散趨近于0。所以,在均勻介質(zhì)模擬且網(wǎng)格間距很小時(shí),常規(guī)差分系數(shù)的模擬結(jié)果優(yōu)于優(yōu)化差分系數(shù)的模擬結(jié)果。常規(guī)差分系數(shù)的計(jì)算公式為:

(19)

為了檢驗(yàn)本文方法,本文采用不同網(wǎng)格間距下的子波能量分布圖和頻散曲線對(duì)比圖對(duì)均勻模型模擬結(jié)果進(jìn)行預(yù)測(cè)分析。圖4a展示了3種優(yōu)化方法在誤差限為0.00001時(shí)的8階優(yōu)化差分系數(shù)頻散曲線。圖4b是不同網(wǎng)格間距下雷克子波能量歸一化結(jié)果。由圖4b可以看出,當(dāng)網(wǎng)格間距為5m時(shí),雷克子波的能量主要集中于kh=0.4波數(shù)段,并且在該波數(shù)段內(nèi),最小二乘法的頻散誤差比其它算法的頻散誤差更小,所以最小二乘法模擬結(jié)果的誤差會(huì)小于其它兩種方法的誤差。對(duì)于雷米茲交換法來(lái)說(shuō),在kh=0.4波數(shù)段上,雷米茲交換法比本文方法的頻散誤差更大,因此,當(dāng)網(wǎng)格間距為5m時(shí),雷米茲交換算法模擬結(jié)果的誤差大于本文方法的誤差。當(dāng)網(wǎng)格間距為10m時(shí),雷克子波的能量主要集中于kh=0.9波數(shù)段,根據(jù)圖4a所示,當(dāng)kh=0.9時(shí),3種方法模擬結(jié)果的誤差比較接近。當(dāng)網(wǎng)格間距為15m時(shí),雷克子波的能量主要集中于kh=1.2波數(shù)段,根據(jù)圖4a所示,在該波數(shù)段上,最小二乘法的頻散誤差也最大,最小二乘法模擬結(jié)果的誤差也最大。在圖4a中,雷米茲交換法的有效帶寬略寬于本文方法的有效帶寬,都在1.2左右,且兩種方法在kh=1.2附近的頻散曲線相似,所以這兩種方法模擬結(jié)果的誤差非常相近。

圖4 3種方法的8階優(yōu)化差分系數(shù)頻散曲線(a)以及不同網(wǎng)格間距下的雷克子波能量歸一化結(jié)果(b)

為了驗(yàn)證前文的分析結(jié)果,本文在200×200的網(wǎng)格上定義一個(gè)二維均勻模型,網(wǎng)格間距分別為5,10,15m。P波速度為1500m/s。震源位于模型中心。雷克子波的主頻為20Hz。圖5a是均勻模型,網(wǎng)格間距為5m,傳播時(shí)間為0.4s,時(shí)間步長(zhǎng)為0.0005s的波場(chǎng)快照對(duì)比結(jié)果。圖5b和圖5c是均勻模型,網(wǎng)格間距分別為10m和15m,傳播時(shí)間為2.0s,時(shí)間步長(zhǎng)為0.0005s的波場(chǎng)快照對(duì)比結(jié)果。圖5d是網(wǎng)格間距為10m時(shí)的波場(chǎng)快照殘差對(duì)比結(jié)果(波場(chǎng)快照殘差等于參考波場(chǎng)快照減去對(duì)應(yīng)方法的波場(chǎng)快照)。均勻模型測(cè)試中將30階常規(guī)差分算子的模擬結(jié)果作為參考波場(chǎng)。

圖5 4種方法的波場(chǎng)快照和波場(chǎng)殘差對(duì)比結(jié)果a 網(wǎng)格間距為5m時(shí)的波場(chǎng)快照對(duì)比結(jié)果; b 網(wǎng)格間距為10m時(shí)的波場(chǎng)快照對(duì)比結(jié)果; c 網(wǎng)格間距為15m時(shí)的波場(chǎng)快照對(duì)比結(jié)果; d 網(wǎng)格間距為10m時(shí)的波場(chǎng)殘差對(duì)比結(jié)果

為了直觀地對(duì)比不同優(yōu)化方法的頻散壓制能力,本文計(jì)算了這些波場(chǎng)快照頻散誤差總和與波場(chǎng)快照總誤差(圖6)。頻散誤差總和是波場(chǎng)快照殘差沿著深度方向的總和?？傉`差就是將頻散誤差總和相加。圖6a、圖6b、圖6c分別是圖5a、圖5b、圖5c中不同方法的波場(chǎng)快照頻散誤差總和對(duì)比結(jié)果。圖6d是圖6a、圖6b、圖6c的總誤差對(duì)比結(jié)果。由圖5和圖6可知,均勻模型數(shù)值模擬結(jié)果基本符合本文圖4的預(yù)測(cè)結(jié)果。

在網(wǎng)格間距為5m時(shí)(圖6a),由于子波主要能量集中于低波數(shù)段,所以模擬結(jié)果誤差與低波數(shù)段頻散呈正相關(guān),相較于雷米茲交換算法,本文方法的模擬結(jié)果誤差較小,說(shuō)明本文方法具有比雷米茲交換算法更小的低波數(shù)段誤差。圖6c中,由于子波能量主要集中于高波數(shù)段,所以模擬結(jié)果誤差與帶寬大小呈正相關(guān),在此前提下,本文方法和雷米茲交換算法具有相近的模擬結(jié)果誤差,說(shuō)明本文方法與雷米茲交換算法具有相似的有效帶寬。均勻模型的測(cè)試結(jié)果驗(yàn)證了本文方法相較于雷米茲交換算法具有較低的低波數(shù)段頻散和相似的有效帶寬。

圖6 不同網(wǎng)格間距的頻散誤差總和和總誤差對(duì)比結(jié)果a 網(wǎng)格間距為5m時(shí)的頻散誤差總和對(duì)比結(jié)果; b 網(wǎng)格間距為10m時(shí)的頻散誤差總和對(duì)比結(jié)果; c 網(wǎng)格間距為15m時(shí)的頻散誤差總和對(duì)比結(jié)果; d 不同網(wǎng)格間距下的總誤差對(duì)比結(jié)果

4.2 改進(jìn)的Marmousi模型

為了在炮記錄中更加直接和準(zhǔn)確地觀察直達(dá)波的頻散,本文在標(biāo)準(zhǔn)的Marmousi模型上加了一層速度為1500m/s的地層。整個(gè)模型被離散為1360×700個(gè)網(wǎng)格點(diǎn)。改進(jìn)的Marmousi模型如圖7所示。圖8、圖9和圖10分別為網(wǎng)格間距為10,20,30m時(shí),不同傳播時(shí)間的波場(chǎng)快照殘差以及炮記錄快照殘差對(duì)比結(jié)果。數(shù)值模擬試驗(yàn)采用的時(shí)間步長(zhǎng)為0.0005s,主頻為20Hz,炮點(diǎn)位置為模型表面的中點(diǎn)。數(shù)值模擬采用吸收邊界條件處理邊界反射。圖8、圖9和圖10 展示的炮記錄殘差由參考炮記錄減去優(yōu)化差分算子的模擬炮記錄得到,波場(chǎng)快照殘差由參考波場(chǎng)快照減去優(yōu)化差分算子的波場(chǎng)快照所得到。其中,優(yōu)化差分算子均是在誤差限為0.0001時(shí)計(jì)算得到的20階優(yōu)化差分算子。參考波場(chǎng)快照和參考炮記錄是40階常規(guī)差分系數(shù)的模擬結(jié)果。

圖7 改進(jìn)的Marmousi速度模型

圖8 網(wǎng)格間距10m時(shí)波場(chǎng)快照殘差(a)和炮記錄快照殘差(b)對(duì)比結(jié)果

圖9 網(wǎng)格間距20m時(shí)波場(chǎng)快照殘差(a)和炮記錄快照殘差(b)對(duì)比結(jié)果

圖10 網(wǎng)格間距30m時(shí)波場(chǎng)快照殘差(a)和炮記錄快照殘差(b)對(duì)比結(jié)果

當(dāng)網(wǎng)格間距為10m時(shí),本文方法的波場(chǎng)快照殘差和炮記錄殘差小于雷米茲交換法和最小二乘法的波場(chǎng)快照殘差和炮記錄殘差,說(shuō)明當(dāng)網(wǎng)格間距為10m時(shí),本文方法的優(yōu)化效果優(yōu)于其它兩種方法的優(yōu)化結(jié)果(圖8)。但當(dāng)網(wǎng)格間距繼續(xù)增大至20m時(shí),本文方法的波場(chǎng)快照殘差和炮記錄殘差與雷米茲交換法的波場(chǎng)快照殘差和炮記錄殘差相似,但是這兩種方法的波場(chǎng)快照殘差和炮記錄殘差小于最小二乘法方法的結(jié)果,說(shuō)明當(dāng)網(wǎng)格間距為20m時(shí),本文方法的優(yōu)化效果與雷米茲交換算法的優(yōu)化效果相似,且兩種方法的優(yōu)化效果優(yōu)于最小二乘法的優(yōu)化效果(圖9)。當(dāng)網(wǎng)格間距為30m時(shí),3種方法的波場(chǎng)快照殘差和炮記錄殘差相似,說(shuō)明網(wǎng)格間距為30m時(shí)3種優(yōu)化差分算子的優(yōu)化效果趨近于一致(圖10)。在改進(jìn)的Marmousi模型的模擬過(guò)程中,本文方法在網(wǎng)格間距較小的情況下,對(duì)于復(fù)雜模型中頻散誤差的壓制效果優(yōu)于最小二乘法和雷米茲交換算法的壓制效果。

5 結(jié)論

為了兼容雷米茲算法的優(yōu)點(diǎn)并減少該算法在低波數(shù)段產(chǎn)生的數(shù)值頻散,本文提出了一種拉格朗日乘數(shù)法與雷米茲交換算法耦合的差分系數(shù)優(yōu)化方法,引入雷米茲交換算法和拉格朗日乘數(shù)法來(lái)優(yōu)化顯式有限差分系數(shù)。頻散分析和數(shù)值模擬試驗(yàn)結(jié)果表明,本文方法具有較寬的有效帶寬和更低的低波數(shù)段數(shù)值頻散,這些特性拓寬了雷米茲交換算法的適用范圍。