BFO-PSO算法下的彈性波數(shù)值模擬

趙平起， 何書(shū)梅， 倪天祿， 趙明， 張家良， 吳吉忠,2，魏朋朋， 李聞達(dá)， 白文磊*

1 中國(guó)石油大港油田分公司， 天津 300280 2 東北石油大學(xué)， 大慶 163318 3 中國(guó)科學(xué)院地質(zhì)與地球物理研究所， 北京 100029 4 中國(guó)科學(xué)院地球科學(xué)研究院， 北京 100029 5 中國(guó)科學(xué)院油氣資源研究院重點(diǎn)實(shí)驗(yàn)室， 北京 100029 6 中國(guó)科學(xué)院大學(xué)， 北京 100049

0 引言

地震波正演模擬是地震勘探反演和成像的基礎(chǔ)和關(guān)鍵，提高彈性波數(shù)值模擬的精度和效率具有非常重要的理論和工程意義.目前廣泛運(yùn)用于地震波數(shù)值模擬的方法主要是有限差分法(Chu and Stoffa，2012；Di Bartolo et al.，2012；Yan et al.，2016；Yang et al.，2017；Liang et al.，2018；He et al.，2019；杜澤源等，2019；Miao et al.，2020)，其適用于GPU并行計(jì)算，計(jì)算速度快并且易于編程實(shí)現(xiàn).

應(yīng)用有限差分算法進(jìn)行彈性波數(shù)值模擬時(shí)，使用差分算子來(lái)代替微分算子，因而必然會(huì)產(chǎn)生數(shù)值誤差，如果忽略這一誤差，將會(huì)嚴(yán)重影響數(shù)值模擬的精度，導(dǎo)致彈性波數(shù)值模擬過(guò)程中出現(xiàn)數(shù)值頻散問(wèn)題(Chu and Stoffa，2012；王之洋等，2015).為了解決差分算子頻散問(wèn)題，可以在數(shù)值模擬時(shí)選擇主頻較低的子波，較小的離散網(wǎng)格間距，或者優(yōu)化差分算子.然而，隨著子波主頻的降低，高頻成分缺失，由此也將導(dǎo)致地震反演和成像的分辨率大大降低；而減小離散網(wǎng)格的間距則會(huì)增加計(jì)算量，對(duì)于大模型來(lái)說(shuō)，由此帶來(lái)的計(jì)算量、存儲(chǔ)量問(wèn)題進(jìn)一步制約勘探技術(shù)的發(fā)展(王之洋等，2015；He et al.，2019；劉立彬等，2020).優(yōu)化差分算子，提高低階數(shù)差分算子的精度，可以在不增加計(jì)算量，保證計(jì)算效率的同時(shí)，減小差分算子對(duì)微分算子的逼近誤差，從而避免數(shù)值模擬過(guò)程中的數(shù)值計(jì)算誤差，因此，優(yōu)化有限差分算子在彈性波動(dòng)方程的數(shù)值求解中具有非常重要的意義.

Chu和Stoffa(2012)指出，根據(jù)Taylor級(jí)數(shù)展開(kāi)，可以使用二項(xiàng)式窗函數(shù)截?cái)鄠巫V法的空間褶積序列推導(dǎo)得出有限差分算子，然而由于截?cái)嘈?yīng)存在，頻譜泄露不可避免，為了減弱頻譜泄露，可以選擇優(yōu)化的窗函數(shù).應(yīng)用窗函數(shù)法優(yōu)化有限差分算子，其主要目標(biāo)是設(shè)計(jì)一種通帶窄，阻帶衰減大的窗函數(shù)，以盡可能的減弱截?cái)嗨鶎?dǎo)致的頻譜泄露(王之洋等，2015；Wang et al.，2017；Ren et al.，2018)，然而，不可忽略的問(wèn)題是，窗函數(shù)的通帶寬度與阻帶衰減是相互矛盾的，很難同時(shí)滿(mǎn)足通帶窄、阻帶衰減大的要求.因此，另一種更常用的優(yōu)化有限差分算子的方法是最優(yōu)化方法(Liu，2013；Zhang and Yao，2013；Yang et al.，2017；He et al.，2019；Miao and Zhang，2020)，通過(guò)構(gòu)造包含有限差分系數(shù)的目標(biāo)函數(shù)，將有限差分算子的優(yōu)化問(wèn)題轉(zhuǎn)化為多參數(shù)優(yōu)化問(wèn)題，進(jìn)而應(yīng)用最優(yōu)化算法進(jìn)行優(yōu)化求解.Zhang和Yao(2013)首次使用模擬退火算法優(yōu)化有限差分算子.He等(2019)將Remez交換算法應(yīng)用于常規(guī)網(wǎng)格和交錯(cuò)網(wǎng)格的有限差分算子優(yōu)化.Miao和Zhang(2020)基于L1范數(shù)構(gòu)造了包含有限差分系數(shù)的目標(biāo)函數(shù)，并應(yīng)用交替方向乘子算法(ADMM)獲得了優(yōu)化的有限差分算子.

應(yīng)用最優(yōu)化方法優(yōu)化有限差分算子，其目標(biāo)是在保證精度誤差的同時(shí)，獲得最大的譜覆蓋范圍(Liu，2013；Zhang and Yao，2013；He et al.，2019)，常規(guī)有限差分優(yōu)化方法往往難以兼顧精度和譜覆蓋范圍兩個(gè)方面，尤其是在高波數(shù)情況下，難以滿(mǎn)足數(shù)值模擬精度的要求.因此，需要尋求一種更適用于高維多參數(shù)問(wèn)題的快速收斂且易于跳出局部極值的高效優(yōu)化算法.

啟發(fā)式優(yōu)化算法是目前廣泛應(yīng)用的一種優(yōu)化算法(雷秀娟等，2012；Beheshti and Shamsuddin，2013；Pillay，2016；Qu et al.，2015；Li et al.，2016；Hussain et al.，2018)，不同于傳統(tǒng)優(yōu)化算法，啟發(fā)式算法通過(guò)在解空間范圍內(nèi)隨機(jī)搜索以獲得最優(yōu)結(jié)果，是一種基于概率計(jì)算的領(lǐng)域隨機(jī)搜索算法，該類(lèi)算法主要包括：模擬退火算法(SA)，粒子群算法(PSO)，遺傳算法(GA)，等等.相比于其他算法，PSO算法在迭代更新時(shí)利用了當(dāng)前最優(yōu)解的信息，從而可以更快收斂到最優(yōu)解(Kennedy and Eberhart，1995；Hazra and Sinha，2011；Jadoun et al.，2015)，然而，PSO算法的全局搜索能力對(duì)初始參數(shù)的要求較高，當(dāng)初始參數(shù)隨機(jī)性較大時(shí)，算法很容易陷入局部極值，尤其是對(duì)于一些參數(shù)維度較高的問(wèn)題、多峰問(wèn)題和病態(tài)問(wèn)題，會(huì)出現(xiàn)過(guò)早收斂甚至難以收斂的問(wèn)題.細(xì)菌覓食算法(BFO)是一種仿生類(lèi)算法(Passino，2002；李珺等，2013；Daryabeigi and Dehkordi，2014)，其模擬大腸桿菌覓食行為，算法主體結(jié)構(gòu)包括三層循環(huán)，分別對(duì)應(yīng)趨化、復(fù)制、驅(qū)散三個(gè)步驟，具有較強(qiáng)的并行搜索能力，同時(shí)，驅(qū)散操作加強(qiáng)了搜索的隨機(jī)性，從而有助于提高算法的全局搜索能力，然而，BFO算法包含三層循環(huán)嵌套，使得算法結(jié)構(gòu)相對(duì)復(fù)雜，引入了較多參數(shù)，導(dǎo)致算法的收斂速度較慢.因此，結(jié)合PSO算法和BFO算法的思想，使其優(yōu)勢(shì)互補(bǔ)，提升性能，改進(jìn)算法可以在快速收斂的同時(shí)搜索到全局最優(yōu)解，以同時(shí)滿(mǎn)足有限差分算子對(duì)精度和頻譜覆蓋范圍的要求.

本文提出一種BFO-PSO算法下的有限差分優(yōu)化方法，并采用優(yōu)化的有限差分算子進(jìn)行彈性波數(shù)值模擬.針對(duì)PSO算法全局搜索能力較弱的問(wèn)題，引入BFO算法中的趨化、復(fù)制、驅(qū)散三個(gè)步驟，形成BFO-PSO混合優(yōu)化算法；構(gòu)造包含有限差分系數(shù)的目標(biāo)函數(shù)，并運(yùn)用BFO-PSO混合優(yōu)化算法求取全局最優(yōu)解，獲得優(yōu)化的有限差分算子，通過(guò)理論頻散分析，比較優(yōu)化性能；進(jìn)而應(yīng)用此優(yōu)化的差分方法分別在層狀介質(zhì)模型和復(fù)雜模型上進(jìn)行彈性波數(shù)值模擬，對(duì)比分析合成地震記錄.

1 理論分析

1.1 BFO-PSO算法下的有限差分算子優(yōu)化

根據(jù)采樣定理，連續(xù)信號(hào)f(x)的一階空間導(dǎo)數(shù)可以表示為(Chu and Stoffa，2012)：

亞當(dāng)·斯密認(rèn)為，人類(lèi)的“愛(ài)”是“有層次的”，而自由市場(chǎng)通過(guò)平等自愿的交換，讓人們能夠“以自私為目的，達(dá)成利他的結(jié)果，最終使每個(gè)人都互惠互利”。從制度經(jīng)濟(jì)學(xué)的角度看，破除以地方政府為代表的政府機(jī)會(huì)主義，重塑政治生態(tài)，有耐于持之以恒、與時(shí)俱進(jìn)地優(yōu)化“規(guī)則”。正如布坎南所指出的那樣：“要改變一種游戲或競(jìng)賽的結(jié)果，改變參加競(jìng)賽的人并不重要，而改變競(jìng)賽規(guī)則最為重要。”

(1)

其中，Δx為空間采樣間隔，Δx/π為Nyquist波數(shù)，fn=f(nΔx).

將x=0代入公式(1)并截?cái)?，可推?dǎo)得出一階空間導(dǎo)數(shù)的有限差分算子：

(2)

對(duì)公式(2)左右兩端同時(shí)應(yīng)用Fourier變換，可推導(dǎo)出：

(3)

其中，kx為波數(shù).

根據(jù)公式(3)，可得出有限差分算子的頻散關(guān)系：

(4)

式中，考慮波數(shù)趨于0的情況，滿(mǎn)足：

進(jìn)一步，應(yīng)用最優(yōu)化算法優(yōu)化有限差分系數(shù)，可以構(gòu)造目標(biāo)函數(shù)為：

(5)

因此，本文結(jié)合BFO算法的全局尋優(yōu)能力，將BFO算法的趨化、復(fù)制、驅(qū)散三個(gè)步驟引入PSO算法中的粒子速度和位置更新策略，以形成BFO-PSO混合優(yōu)化算法，進(jìn)一步應(yīng)用BFO-PSO算法對(duì)包含有限差分系數(shù)的目標(biāo)函數(shù)(公式(5))進(jìn)行求解，從而得到精度誤差較小且頻譜覆蓋范圍較大的有限差分算子.

PSO算法是一種基于概率計(jì)算的領(lǐng)域隨機(jī)搜索算法，種群中的每一個(gè)粒子都代表目標(biāo)函數(shù)的一組解，在該算法中，每個(gè)粒子都以隨機(jī)的方式向當(dāng)前最優(yōu)粒子的方向移動(dòng)，從而使得整個(gè)種群向最優(yōu)解的方向運(yùn)動(dòng)，進(jìn)而實(shí)現(xiàn)在目標(biāo)函數(shù)解空間范圍內(nèi)以隨機(jī)搜索的方式求解最優(yōu)解.常規(guī)PSO算法的粒子速度和位置更新策略如下：

(6)

為了克服PSO算法對(duì)初始參數(shù)的高度依賴(lài)性以及難以跳出局部極值的缺陷，本文在PSO算法的粒子速度和位置更新方式中增加BFO算法的趨化、復(fù)制、驅(qū)散三種步驟，以提高算法在鄰域搜索時(shí)的隨機(jī)性，形成BFO-PSO混合優(yōu)化算法，從而提高PSO算法的全局搜索能力，降低算法對(duì)初始參數(shù)的依賴(lài)性，同時(shí)進(jìn)一步提高算法的收斂速度.

首先引入趨化步驟，在粒子速度和位置更新過(guò)程中，粒子種群向著全局最優(yōu)解的方向運(yùn)動(dòng)，但是，對(duì)于一些參數(shù)維度較高的問(wèn)題、多峰問(wèn)題和病態(tài)問(wèn)題，粒子種群很容易趨向局部最優(yōu)，從而導(dǎo)致算法全局搜索能力較差.為了避免這一問(wèn)題，在粒子位置更新時(shí)增加隨機(jī)方向性，使粒子隨機(jī)地向一個(gè)方向運(yùn)動(dòng)或者翻轉(zhuǎn)，從而增加粒子種群在全局范圍內(nèi)的搜索隨機(jī)性，使其具有跳出局部最優(yōu)的能力，增加算法收斂到全局最優(yōu)解的概率.引入趨化步驟的粒子位置更新公式為：

(7)

其次是復(fù)制步驟，將BFO算法中細(xì)菌優(yōu)勝劣汰的繁殖過(guò)程引入到PSO算法粒子速度和位置更新策略中，對(duì)種群中的每個(gè)粒子計(jì)算適應(yīng)度值，并根據(jù)目標(biāo)函數(shù)要求按照從優(yōu)到劣的順序?qū)ζ渑判?，并將排序在末尾的一半種群粒子淘汰，而將剩余的種群粒子進(jìn)行完全復(fù)制，使得種群大小不變.通過(guò)這一步驟，可以保證粒子種群在每一次迭代更新時(shí)都向著部分較優(yōu)解的方向運(yùn)動(dòng)，從而有效提高粒子種群向最優(yōu)解方向移動(dòng)的速度，提高算法的收斂速度.

最后是驅(qū)散步驟，為了進(jìn)一步放置粒子種群過(guò)早的收斂到局部最優(yōu)位置，避免種群粒子聚集在局部極值位置處，按照一定的概率將種群粒子隨機(jī)移動(dòng)到解空間中的新區(qū)域，以此提高算法在鄰域內(nèi)搜索的隨機(jī)性，從而增加粒子種群的隨機(jī)性，避免種群過(guò)早的收斂到局部極值，提高全局搜索能力.

根據(jù)上述分析，將BFO算法的趨化、復(fù)制、驅(qū)散三種步驟引入PSO算法，改進(jìn)PSO算法的粒子速度和位置更新策略，形成BFO-PSO混合優(yōu)化算法，算法流程圖如圖1所示.

圖1 BFO-PSO算法流程圖Fig.1 Flowchart of the BFO-PSO algorithm

1.2 理論頻散分析

根據(jù)公式(3)所示有限差分算子的頻散關(guān)系式，分析比較BFO-PSO算法，PSO算法以及Remez交換算法(He et al.，2019)下的優(yōu)化有限差分算子的數(shù)值頻散曲線(xiàn)，比較其精度誤差.

分別應(yīng)用BFO-PSO算法和PSO算法對(duì)包含有限差分系數(shù)的目標(biāo)函數(shù)(公式(5))進(jìn)行優(yōu)化求解，以獲得BFO-PSO算法下的優(yōu)化有限差分算子.本文中，BFO-PSO算法和PSO算法的初始參數(shù)設(shè)定如下：粒子位置表示有限差分系數(shù)，并進(jìn)行隨機(jī)初始化；粒子速度表示有限差分系數(shù)的調(diào)整方向；種群大小為N=300，最大迭代次數(shù)為T(mén)=1500，權(quán)重系數(shù)為ω=0.7，粒子驅(qū)散概率為0.25.

如圖2、3所示，分別為BFO-PSO算法和PSO算法優(yōu)化得到的一階空間導(dǎo)數(shù)有限差分算子的數(shù)值頻散曲線(xiàn)，其中，圖2應(yīng)用BFO-PSO算法，圖3應(yīng)用PSO算法，圖中黑色實(shí)線(xiàn)表示不同階數(shù)下常規(guī)有限差分算子的數(shù)值頻散曲線(xiàn).對(duì)比分析，可以明顯看出，采用BFO-PSO算法得到的優(yōu)化有限差分算子在保證精度誤差的同時(shí)，其頻譜覆蓋范圍遠(yuǎn)大于常規(guī)有限差分算子.8階BFO-PSO算法優(yōu)化的有限差分算子的頻散曲線(xiàn)甚至接近16階常規(guī)有限差分算子的頻散曲線(xiàn)，兩者的頻譜覆蓋范圍幾乎一致，而12階BFO-PSO算法優(yōu)化的有限差分算子的頻散曲線(xiàn)則遠(yuǎn)遠(yuǎn)優(yōu)于20階常規(guī)有限差分算子的頻散曲線(xiàn)，在12階時(shí)就可以滿(mǎn)足高波數(shù)的精度要求.

圖2 不同階數(shù)的BFO-PSO算法下的優(yōu)化有限差分算子的頻散曲線(xiàn)Fig.2 Dispersion curves of the optimized FD operator based on the BFO-PSO algorithm for different values of N (different orders)

圖3 不同階數(shù)的PSO算法下的優(yōu)化有限差分算子的頻散曲線(xiàn)Fig.3 Dispersion curves of the optimized FD operator based on the PSO algorithm for different values of N (different orders)

進(jìn)一步分析比較BFO-PSO混合優(yōu)化算法的全局搜索能力和收斂速度，圖4給出了BFO-PSO算法與PSO算法優(yōu)化獲得的有限差分算子的數(shù)值頻散曲線(xiàn)對(duì)比.可以看出，在8階、12階、16階、20階的情況下，BFO-PSO算法和PSO算法優(yōu)化得到的有限差分算子，其數(shù)值頻散曲線(xiàn)幾乎完全一致，即兩種優(yōu)化算法得到的有限差分算子是一致的.然而，當(dāng)差分算子的階數(shù)繼續(xù)增加到24階時(shí)，應(yīng)用BFO-PSO算法得到的優(yōu)化有限差分算子在保持精度誤差不變的同時(shí)，具有更大的頻譜覆蓋范圍.基于這一分析，可以說(shuō)明，引入趨化、復(fù)制、驅(qū)散三種步驟后形成的BFO-PSO混合優(yōu)化算法很好的結(jié)合了PSO算法的快速收斂性和BFO算法的全局搜索能力，對(duì)于多參數(shù)優(yōu)化問(wèn)題，BFO-PSO混合優(yōu)化算法可以快速收斂到較好的全局最優(yōu)解，就本文而言，應(yīng)用BFO-PSO算法可以很快的獲得精度誤差小且適應(yīng)高波數(shù)情況的優(yōu)化有限差分算子.

圖4 BFO-PSO算法與PSO算法優(yōu)化的有限差分算子的頻散曲線(xiàn)對(duì)比Fig.4 Comparison of numerical dispersion curves between optimized FD operator based on the BFO-PSO algorithm and the PSO algorithm

圖5為分別應(yīng)用BFO-PSO算法和Remez交換算法得到的優(yōu)化有限差分算子的數(shù)值頻散曲線(xiàn)對(duì)比，其中，彩色實(shí)線(xiàn)分別對(duì)應(yīng)不同階數(shù)的BFO-PSO算法優(yōu)化的有限差分算子的頻散曲線(xiàn)，彩色虛線(xiàn)分別對(duì)應(yīng)不同階數(shù)的基于Remez交換算法的有限差分算子的頻散曲線(xiàn).對(duì)比可得，在保證精度誤差控制在一定范圍內(nèi)的同時(shí)，BFO-PSO算法優(yōu)化后的有限差分算子具有更大的頻譜覆蓋范圍，更適用于高波數(shù)情況下的數(shù)值模擬.表1列出了應(yīng)用BFO-PSO算法得到的優(yōu)化有限差分算子的系數(shù).

圖5 BFO-PSO算法與Remez交換算法(He et al.，2019)優(yōu)化的有限差分算子的頻散曲線(xiàn)對(duì)比Fig.5 Comparison of numerical dispersion curves between the optimized FD operator based on the BFO-PSO algorithm and the Remez algorithm (He et al., 2019)

表1 BFO-PSO算法下的優(yōu)化有限差分算子系數(shù)Table 1 Optimized FD coefficients based on the BFO-PSO algorithm

應(yīng)用有限差分法進(jìn)行數(shù)值模擬時(shí)，差分格式的穩(wěn)定性是一個(gè)非常重要的問(wèn)題.這里假定空間步長(zhǎng)均勻，即Δx=Δz，介質(zhì)最大波速為v，時(shí)間步長(zhǎng)為Δt，給出常規(guī)網(wǎng)格一階有限差分算子的簡(jiǎn)單穩(wěn)定性條件(Lines et al.，1999)：

(8)

2 數(shù)值模擬與分析

2.1 層狀介質(zhì)模型

首先在應(yīng)用BFO-PSO算法優(yōu)化的有限差分算子在各向同性雙層介質(zhì)模型上進(jìn)行彈性波數(shù)值模擬測(cè)試，模型示意圖如圖6所示，模型上下兩層的縱波速度、橫波速度、密度分別在圖中標(biāo)出，模型寬和高均為4000 m，空間步長(zhǎng)設(shè)置為Δx=Δz=8 m，時(shí)間采樣間隔為Δt=0.0005 s.震源坐標(biāo)為(2000 m,1000 m)，采用主頻為25 Hz的Ricker子波.

圖6 雙層介質(zhì)模型Fig.6 Double-layer velocity model

圖7和圖8為分別應(yīng)用BFO-PSO算法優(yōu)化得到的有限差分算子和常規(guī)有限差分算子進(jìn)行彈性波數(shù)值模擬所得到的合成波場(chǎng)快照的X分量和Z分

量，其中，圖7a和圖8a應(yīng)用8階常規(guī)有限差分算子，圖7b和圖8b應(yīng)用BFO-PSO算法優(yōu)化得到的8階有限差分算子，圖7c和圖8c應(yīng)用16階常規(guī)有限差分算子.對(duì)比分析發(fā)現(xiàn)，應(yīng)用8階常規(guī)有限差分算子獲得的合成波場(chǎng)快照，無(wú)論是X分量還是Z分量，均存在明顯的數(shù)值頻散現(xiàn)象，如圖中黑色箭頭指示.而應(yīng)用BFO-PSO算法優(yōu)化得到的有限差分算子進(jìn)行彈性波數(shù)值模擬，所獲得合成波場(chǎng)快照中幾乎沒(méi)有數(shù)值頻散，數(shù)值模擬的精度甚至達(dá)到了16階常規(guī)差分算子的效果.因此，采用BFO-PSO算法可以有效的提高有限差分算子的精度和效率，優(yōu)化獲得的低階有限差分算子在逼近微分算子時(shí)就具有較小的數(shù)值誤差，進(jìn)而可以在不增加計(jì)算量，保證計(jì)算效率的同時(shí)，達(dá)到常規(guī)高階有限差分算子的精度，滿(mǎn)足高精度數(shù)值模擬的要求，有效提高彈性波數(shù)值模擬的精度和效率.

圖7 應(yīng)用不同有限差分算子得到的合成波場(chǎng)快照(X分量)(a) 8階常規(guī)有限差分算子； (b) 8階BFO-PSO算法優(yōu)化的有限差分算子； (c) 16階常規(guī)有限差分算子.Fig.7 Synthetic wavefield snapshots (X component) using different FD operators(a) Eighth-order conventional FD operator; (b) Eighth-order FD operator optimized based on the BFO-PSO algorithm; (c) Sixteenth-order conventional FD operator.

圖8 應(yīng)用不同有限差分算子得到的合成波場(chǎng)快照(Z分量)(a) 8階常規(guī)有限差分算子； (b) 8階BFO-PSO算法優(yōu)化的有限差分算子； (c) 16階常規(guī)有限差分算子.Fig.8 Synthetic wavefield snapshots (Z component) using different FD operators(a) Eighth-order conventional FD operator; (b) Eighth-order FD operator optimized based on the BFO-PSO algorithm; (c) Sixteenth-order conventional FD operator.

2.2 Marmousi模型

之后，我們?cè)谳^為復(fù)雜的Marmousi模型上進(jìn)行彈性波數(shù)值模擬，以進(jìn)一步對(duì)比分析BFO-PSO算法優(yōu)化得到的有限差分算子在壓制數(shù)值頻散、提升彈性波數(shù)值模擬效率方面的性能.如圖9所示，為Marmousi縱波速度模型，在本文數(shù)值模擬實(shí)驗(yàn)中，模型橫波速度和縱波速度的比值為vP/vS=1.7，密度設(shè)置為ρ=2000 kg·m-3.模型空間離散網(wǎng)格大小設(shè)置767×376，空間步長(zhǎng)設(shè)置為Δx=12 m,Δz=8 m，時(shí)間采樣間隔為Δt=0.001 s，記錄時(shí)長(zhǎng)為t=5 s.震源坐標(biāo)為(4600 m,1 m)，采用主頻為25 Hz的Ricker子波.

圖9 Marmousi模型Fig.9 Marmousi model

圖10為分別應(yīng)用BFO-PSO算法優(yōu)化得到的有限差分算子和常規(guī)有限差分算子在Marmousi模型上進(jìn)行彈性波數(shù)值模擬所得到的合成地震記錄.其中，圖10a、b、c分別為應(yīng)用8階常規(guī)有限差分算子、8階BFO-PSO算法優(yōu)化得到的有限差分算子、16階常規(guī)有限差分算子得到的合成地震記錄，圖10d、e分別為區(qū)塊1和2的放大記錄.根據(jù)合成地震記錄的對(duì)比分析，可以發(fā)現(xiàn)，應(yīng)用BFO-PSO算法優(yōu)化的有限差分算子得到的地震記錄具有更好的精度和效果，幾乎與16階常規(guī)有限差分算子的精度一致，而8階常規(guī)算子得到的地震記錄上則可以清晰看到數(shù)值頻散存在.因此，通過(guò)在層狀介質(zhì)模型和復(fù)雜的Marmousi模型上進(jìn)行的彈性波數(shù)值模擬測(cè)試，可以充分說(shuō)明，BFO-PSO算法優(yōu)化得到的有限差分算子具有較高的精度，可以有效壓制用差分算子替代微分算子時(shí)導(dǎo)致的數(shù)值頻散，優(yōu)化后的低階差分算子就可以達(dá)到常規(guī)高階算子的精度，從而可以有效提高對(duì)彈性波動(dòng)方程求取數(shù)值解時(shí)的計(jì)算效率，降低對(duì)于計(jì)算資源的要求.

圖10 在Marmousi模型應(yīng)用不同有限差分算子進(jìn)行彈性波數(shù)值模擬得到的合成地震記錄(Z分量)(a) 8階常規(guī)有限差分算子； (b) 8階BFO-PSO算法優(yōu)化的有限差分算子； (c) 16階常規(guī)有限差分算子； (d) 區(qū)塊1放大記錄； (e) 區(qū)塊2放大記錄.圖(d)、(e)中，左邊為8階常規(guī)有限差分算子，中間為8階BFO-PSO算法優(yōu)化的有限差分算子，右邊為16階常規(guī)有限差分算子.Fig.10 Synthetic shot records (Z component) obtained by numerical modelling of elastic waves using different FD operators on Marmousi model(a) Eighth-order conventional FD operator; (b) Eighth-order FD operator optimized based on the BFO-PSO algorithm; (c) Sixteenth-order conventional FD operator; (d) Zoomed view of block 1; (e) Zoomed view of block 2. In (d) and (e), the left is eighth-order conventional FD operator; middle is eighth-order FD operator based on the BFO-PSO algorithm; and right is sixteenth-order conventional FD operator.

3 結(jié)論

本文提出一種BFO-PSO算法下的有限差分算子優(yōu)化方法，并應(yīng)用此優(yōu)化的有限差分算子分別在層狀介質(zhì)模型和復(fù)雜模型上進(jìn)行彈性波數(shù)值模擬，對(duì)比分析合成地震記錄.首先，針對(duì)PSO算法對(duì)初始參數(shù)依賴(lài)性大及算法易于陷入局部極值的缺點(diǎn)，引入BFO算法的趨化、復(fù)制、驅(qū)散三種步驟，以改進(jìn)PSO算法的粒子速度和位置更新策略，形成BFO-PSO混合優(yōu)化算法.其次，構(gòu)造包含有限差分系數(shù)的目標(biāo)函數(shù)，并應(yīng)用BFO-PSO混合優(yōu)化算法求取最優(yōu)解，獲得優(yōu)化的有限差分算子，根據(jù)理論頻散曲線(xiàn)分析，BFO-PSO算法可以在快速收斂的同時(shí)獲得較好的全局最優(yōu)解，對(duì)于低階差分算子，在保證精度誤差的同時(shí)，有效擴(kuò)大了頻譜覆蓋范圍，滿(mǎn)足彈性波數(shù)值模擬對(duì)高波數(shù)的要求；對(duì)于高階差分算子，BFO-PSO算法優(yōu)化得到的差分系數(shù)精度更好，即BFO-PSO算法更適用于復(fù)雜的多參數(shù)優(yōu)化問(wèn)題.最后，應(yīng)用BFO-PSO算法優(yōu)化后的有限差分算子進(jìn)行彈性波數(shù)值模擬，對(duì)比分析合成波場(chǎng)快照及合成地震記錄，進(jìn)一步驗(yàn)證了BFO-PSO混合優(yōu)化算法獲得的優(yōu)化有限差分算子在壓制數(shù)值頻散方面的有效性，降低了彈性波動(dòng)方程數(shù)值求解的誤差，提高數(shù)值模擬的精度和效率.