一種雷米茲與拉格朗日耦合的有限差分系數(shù)優(yōu)化方法

彭煒颋,黃建平,2

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

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

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

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

1 雷米茲交換法

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

(1)

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

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

(2)

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

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

(3)

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

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

(4)

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

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

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

(5)

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

(6)

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

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

(7)

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

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

(8)

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

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

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

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

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

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

(9)

將公式(9)改寫為:

(10)

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

(11)

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

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

(12)

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

(13)

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

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

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

(14)

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

(15)

將公式(15)展開:

(16)

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

(17)

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

3 理論誤差分析

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

(18)

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

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

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

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

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

4 數(shù)值模擬

4.1 均勻模型

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

(19)

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

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

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

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

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

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

4.2 改進的Marmousi模型

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

圖7 改進的Marmousi速度模型

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

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

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

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

5 結(jié)論

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