采用混合共軛梯度迭代的頻域地震數(shù)值模擬

劉文革 周覓路 彭浩天 劉福烈 牟其松

(①西南石油大學(xué)地球科學(xué)與技術(shù)學(xué)院,四川成都 610500； ②中國石油西南油氣田分公司勘探開發(fā)研究院，四川成都 610000)

0 引言

地震全波形反演利用完整的波場信息估算地下介質(zhì)的屬性模型，與傳統(tǒng)方法相比，精度更高。全波形反演可在時(shí)間域或頻率域進(jìn)行，但無論在哪個(gè)域，計(jì)算成本都很高，這限制了它在實(shí)際生產(chǎn)中的應(yīng)用。地震數(shù)值模擬作為全波形反演的重要組成部分，在很大程度上影響著反演精度和效率。當(dāng)前波形反演中所使用的數(shù)值模擬方法按計(jì)算域可分為時(shí)間域和頻率域。頻率域的優(yōu)勢在于頻率分量選取靈活、計(jì)算成本低。在具體的正演方法實(shí)現(xiàn)上，有限差分方法具有實(shí)現(xiàn)簡單、占用內(nèi)存小等特點(diǎn)被廣泛使用[1]。

1972年，Lysmer等[2]提出了頻率—空間域正演模擬方法； 隨后Marfurt[3]定量對(duì)比了時(shí)間域與頻率域的聲波和彈性波數(shù)值模擬方法，指出頻率—空間域正演模擬比較適合多炮激發(fā)，且容易解決計(jì)算的不穩(wěn)定性問題。由于頻率域方法相較于時(shí)間域方法，頻散現(xiàn)象更嚴(yán)重，為了降低頻散誤差，研究人員做了許多嘗試。Jo等[4]提出了頻率域波動(dòng)方程的最優(yōu)化9點(diǎn)差分格式，頻散現(xiàn)象得到了較好壓制，即每一波長內(nèi)分布有4 個(gè)網(wǎng)格點(diǎn)就能使相速度誤差限制在1% 以內(nèi)； Shin等[5]在此基礎(chǔ)上提出了25點(diǎn)優(yōu)化差分算子，使每一波長內(nèi)的網(wǎng)格點(diǎn)數(shù)可進(jìn)一步降低到2.5； 吳國忱等[6]在此基礎(chǔ)上建立了25點(diǎn)優(yōu)化差分算子，實(shí)現(xiàn)了VTI介質(zhì)準(zhǔn)P波的數(shù)值模擬； 劉璐等[7]提出了優(yōu)化15點(diǎn)差分格式，有較高的計(jì)算效率和較好的頻散壓制效果； 范娜等[8]提出了一種三維聲波方程有限差分系數(shù)優(yōu)化算法，在給定有限差分模板情況下即可求解出優(yōu)化系數(shù)，較傳統(tǒng)的差分系數(shù)具有更低的數(shù)值頻散。唐超等[9]基于L1范數(shù)建立目標(biāo)函數(shù)，采用交替方向乘子法(ADMM)求解交錯(cuò)網(wǎng)格有限差分系數(shù)，提出了一種新的聲波方程交錯(cuò)網(wǎng)格優(yōu)化有限差分正演模擬方法，數(shù)值試驗(yàn)表明該方法能將頻散控制得更低。

以上方法大都基于方形網(wǎng)格，而實(shí)際模型通常是矩形網(wǎng)格。為此，Tang等[10]構(gòu)建了一種二維自適應(yīng)17點(diǎn)有限差分格式，該方法不僅適用于正方形網(wǎng)格，也適應(yīng)于矩形網(wǎng)格，且將每一波長內(nèi)的網(wǎng)格點(diǎn)數(shù)降低到2.4； Xu等[11]提出了一種自適應(yīng)的9點(diǎn)有限差分格式，不僅適用于矩形網(wǎng)格，而且還具有一般25點(diǎn)格式的精度和最優(yōu)化9點(diǎn)差分格式的效率。

另一方面，為提高頻率域正演的計(jì)算效率，也出現(xiàn)了一些較實(shí)用的技術(shù)策略，如：殷文[12]針對(duì)頻率—空間域有限差分彈性波的數(shù)值模擬，采用了流水線技術(shù)和分治策略進(jìn)行并行計(jì)算，提高了計(jì)算效率； 有研究人員實(shí)現(xiàn)了基于多重網(wǎng)格預(yù)條件的迭代法正演，并在標(biāo)準(zhǔn)模型上得到了較精確的結(jié)果[13-14]； Du等[15]基于雙共軛梯度穩(wěn)定(BiCGstab)算法進(jìn)行二維聲波方程的數(shù)值模擬，在保證一定精度的條件下，相較于LU直接分解法大大減少了內(nèi)存消耗和計(jì)算時(shí)間； Belonosov等[16]同樣基于BiCGstab算法，提出了適用于三維各向同性非均質(zhì)彈性波模擬的迭代求解器，結(jié)合混合并行技術(shù)(MPI與OpenMP相結(jié)合)，加快了收斂速度、提高了效率； 熊治濤等[17]在大地電磁響應(yīng)計(jì)算中采用不完全LU直接分解法預(yù)條件因子的BiCGstab算法求解有限元方程，也得到了可靠的結(jié)果； Jiang等[18]開發(fā)了一種頻率域多尺度有限差分方法，能夠以較小的內(nèi)存和較短的時(shí)間求解Helmholtz方程，可在低頻情況下得到精確的解； Jaysaval等[19]使用Schur補(bǔ)迭代方法(即廣義極小殘量法)，降低了求解Helmholtz方程的復(fù)雜度，提高了頻率域數(shù)值模擬的計(jì)算效率； 顏紅芹[20]在利用聲波散射理論進(jìn)行VSP波場模擬的過程中，引入預(yù)條件算子和最小二乘方法，在計(jì)算優(yōu)化的同時(shí)使過程更加穩(wěn)定。上述策略均是基于固定網(wǎng)格進(jìn)行的優(yōu)化，有學(xué)者受到時(shí)間域數(shù)值模擬變網(wǎng)格技術(shù)的啟發(fā)，還提出了頻率域的變網(wǎng)格模擬方法，為提高頻率域地震正演效率提供了新的思路：韓利等[21]根據(jù)頻率域正演中不同頻率間相互獨(dú)立的特點(diǎn)，在低頻段采用大網(wǎng)格、高頻段采用小網(wǎng)格計(jì)算，實(shí)現(xiàn)了變網(wǎng)格步長模擬，在保證模擬精度的同時(shí)，減少了計(jì)算量和內(nèi)存消耗； 郭振波等[22]進(jìn)一步發(fā)展了變網(wǎng)格技術(shù)，構(gòu)建了起伏地表?xiàng)l件下的頻域正演算法，計(jì)算效率可提高5倍以上，同時(shí)內(nèi)存占用大幅減少； 吳悠等[23]結(jié)合交錯(cuò)網(wǎng)格和混合網(wǎng)格完成了頻率—空間域非均質(zhì)聲波方程有限差分模擬，有效提高了計(jì)算效率。

在頻率域正演迭代法中，基于Krylov子空間的方法應(yīng)用較為廣泛。基于Krylov子空間法的BiCGstab算法因內(nèi)存消耗小、收斂速度快，成了許多學(xué)者的選擇，但在實(shí)現(xiàn)時(shí)往往需要設(shè)計(jì)較為復(fù)雜的預(yù)處理器，且在處理復(fù)雜模型時(shí)容易出現(xiàn)迭代不收斂的情況。針對(duì)現(xiàn)有頻率域地震正演方法在求解時(shí)諧波動(dòng)方程所遇到的問題，同時(shí)考慮計(jì)算的效率和穩(wěn)定性，本文在傳統(tǒng)的共軛梯度(CG)算法基礎(chǔ)上，采用合理、優(yōu)化的最小二乘共軛梯度參數(shù)，發(fā)展了一種混合共軛梯度迭代 (MLSCG) 算法。試驗(yàn)結(jié)果表明，本文算法在保證計(jì)算精度的同時(shí)可高效地實(shí)現(xiàn)數(shù)值模擬，并且克服了傳統(tǒng)迭代方法計(jì)算不穩(wěn)定的缺點(diǎn)。

1 方法原理

1.1 波動(dòng)方程的有限差分離散

首先，將二維時(shí)間域聲波方程進(jìn)行傅里葉變換，得到頻率域波動(dòng)方程

(1)

式中：ω為角頻率；P(x,z,ω)為壓力場；v為介質(zhì)速度；S(ω)表示頻率域震源項(xiàng)；xs、zs分別為震源橫、縱坐標(biāo)。

本文使用最優(yōu)化9點(diǎn)有限差分格式[4]，其最終的差分方程為

(2)

式中：h是空間網(wǎng)格間隔；Pm,n表示點(diǎn)(m,n)處的壓力場；a=0.5461、b=0.6248、c=0.09381是最優(yōu)化差分系數(shù)。

時(shí)諧、常密度聲波方程也可用Helmholtz方程表示為

(3)

式中：L為波算子；K=ω/v為波數(shù)； Δ是拉普拉斯算子。

使用有限差分將式(3)離散，得到線性系統(tǒng)

Ap=b

(4)

式中：A是大型稀疏復(fù)值矩陣； 待求單頻波場向量p由P離散得到； 向量b由S離散得到。設(shè)計(jì)算網(wǎng)格的離散點(diǎn)數(shù)為i，空間維數(shù)是D，則A的尺度是iD×iD。頻率域聲波方程轉(zhuǎn)化為線性方程組之后，如何求解該方程系統(tǒng)十分關(guān)鍵。一般的直接法是使用LU分解法，該方法計(jì)算精度高但占用內(nèi)存大，耗費(fèi)時(shí)間長，且難以模擬高維度、高密度采集； Krylov子空間迭代法中使用較多的是BiCGstab算法，該方法計(jì)算效率高，收斂速度快，但對(duì)計(jì)算穩(wěn)定性條件要求較高，且容易出現(xiàn)迭代不收斂。為此，本文引入最小二乘混合共軛梯度迭代求解方法，在保持計(jì)算效率的條件下有效提高了計(jì)算的穩(wěn)定性。

1.2 混合共軛梯度迭代算法

計(jì)算數(shù)學(xué)中對(duì)大型不定方程組的求解，一般采用Krylov子空間迭代法，這類算法實(shí)際上是共軛梯度算法的推廣。算法從初始值p0開始，迭代計(jì)算更新pk(其中k為迭代次數(shù))，直到剩余誤差|b-Apk|達(dá)到足夠小時(shí)求得近似解。若沒有預(yù)處理器，這種方法收斂緩慢或根本不收斂[24]。為了提高計(jì)算穩(wěn)定性和收斂速度，往往需要構(gòu)造合理的預(yù)處理器。

CG算法是由Hestenes等[25]提出，其共軛梯度參數(shù)(CG(HS))定義為

(5)

(6)

式中M為終止迭代次數(shù)。

由于CG算法的收斂速度較慢， Dai等[26]提出了修正的共軛梯度參數(shù)(CG(DY))

(7)

該算法在一定程度上可以提高收斂速度[27]。

為了優(yōu)化計(jì)算性能，本文給出一種混合算法，即MLSCG，有效結(jié)合了式(5)和式(7)各自的特性，新的共軛梯度參數(shù)為

(8)

式中θk為最小二乘系數(shù)，可以通過最小二乘問題求出

(9)

表1 最小二乘系數(shù)計(jì)算時(shí)間對(duì)比統(tǒng)計(jì)

使用迭代算法求解線性系統(tǒng)Ap=b，需要設(shè)置迭代誤差ε用于控制迭代過程。迭代誤差定義為

(12)

由于矩陣A具有對(duì)角占優(yōu)且高度稀疏的特點(diǎn)，所以本文使用雅各比預(yù)處理算子進(jìn)行優(yōu)化，以加速迭代收斂。

為了說明MLSCG算法的特性，設(shè)計(jì)了1000×1000的稀疏矩陣求解問題。準(zhǔn)確解為單位向量，迭代初值為零向量，最大相對(duì)誤差是1.0×10-7。使用混合算法、BiCGstab算法等四種方法進(jìn)行計(jì)算，所得到的迭代收斂曲線如圖1所示。由圖可見：CG(HS)算法收斂速度較慢，但迭代過程相對(duì)穩(wěn)定； CG(DY)收斂速度快，但穩(wěn)定性條件較為苛刻； BiCGstab算法從曲線上看并不是很穩(wěn)定，有一定的鋸齒現(xiàn)象(在后續(xù)的模型試驗(yàn)中會(huì)有說明)； MLSCG算法在提高收斂速度的同時(shí)，保持了算法的穩(wěn)定性，并且在處理復(fù)雜線性系統(tǒng)時(shí)能夠正確地得到近似解。

圖1 不同算法的迭代收斂曲線

1.2.1 時(shí)間復(fù)雜度分析

時(shí)間復(fù)雜度是指執(zhí)行算法所需要的計(jì)算工作量。設(shè)N是模型離散化后系數(shù)矩陣的階數(shù)，M是總迭代次數(shù)。如果僅考慮單一炮記錄的數(shù)值模擬，LU分解法和MLSCG算法的時(shí)間復(fù)雜度對(duì)比分析如表2所示。可以看出，LU分解法的時(shí)間復(fù)雜度不管是在2D還是3D情況下都遠(yuǎn)大于MLSCG算法，這是因?yàn)镸LSCG算法并不需要對(duì)線性系統(tǒng)進(jìn)行分解，同時(shí)值得注意的是，在3D情況下，LU分解法的時(shí)間復(fù)雜度過高，導(dǎo)致其難以推廣應(yīng)用。

表2 兩種方法時(shí)間復(fù)雜度、內(nèi)存需求對(duì)比

1.2.2 內(nèi)存需求分析

LU分解法和MLSCG算法的實(shí)際內(nèi)存需求如表2所示。兩種方法的內(nèi)存需求與N階矩陣規(guī)模緊密相關(guān)，而MLSCG算法的內(nèi)存需求和迭代次數(shù)無直接關(guān)系。在2D情況下，MLSCG算法的內(nèi)存需求稍小于LU算法； 但在3D情況下，MLSCG的內(nèi)存優(yōu)勢較為明顯。在模型實(shí)例中會(huì)有詳細(xì)的兩種方法的實(shí)際計(jì)算時(shí)間和內(nèi)存消耗對(duì)比。

2 數(shù)值算例

本文所有數(shù)值模擬的硬件參數(shù)為Intel(R)Xeon (R) Silver 4110 CPU @ 2.1GHz(8核16線程)，RAM 32GB，模擬過程在單線程下進(jìn)行。

2.1 層狀介質(zhì)模型

為了驗(yàn)證頻域數(shù)值方法的正確性和適用性，本文設(shè)計(jì)了層狀聲學(xué)介質(zhì)模型(圖2)進(jìn)行數(shù)值模擬。模型網(wǎng)格數(shù)為300×300； 網(wǎng)格尺寸為5m×5m； 震源為Ricker子波，主頻為30Hz，位于(750m，750m)； 模型四周采用完全匹配層(Perfectly Matched Layer，PML)吸收邊界條件，寬度為50個(gè)網(wǎng)格。通過數(shù)值模擬，可以得到不同方法的單頻波場數(shù)據(jù)片(圖3)和合成波場快照(圖4)。

圖2 層狀介質(zhì)模型

由圖3可見：①M(fèi)LSCG迭代法在ε=0.01%條件下，得到的10、30、50Hz的波場切片與LU直接法的結(jié)果基本一致； ②從波場殘差來看，其振幅的相對(duì)變化范圍在2%～4%，由此可證明本文算法具有較好的計(jì)算精度； ③10、30、50Hz的波場幅值符合震源傅里葉變換后對(duì)應(yīng)頻率的振幅值； ④在模型的速度分界面處波形出現(xiàn)相應(yīng)的變化，上層波場受分界面處反射波的影響出現(xiàn)了局部的干涉現(xiàn)象； ⑤30和50Hz單頻波場殘差變化具有一定的空間分布特征，以震源為中心在局部呈輻射狀。

圖3 不同頻率MLSCG法(左)、LU法(中)單頻實(shí)部波場切片及其對(duì)應(yīng)的殘差對(duì)比(右)

將不同的單頻波場加權(quán)組合后可以得到時(shí)間域的合成波場快照(圖4)。圖4比較了3種迭代誤差條件下MLSCG算法與LU分解法波場快照的差異性。需要說明的是，此處的波場記錄是由60個(gè)單頻波場分量(1～60Hz)疊加得到，波場快照時(shí)刻為200ms。由圖4可見：①當(dāng)ε=0.100%時(shí)，與LU法結(jié)果相比，MLSCG迭代法波場快照出現(xiàn)了較為明顯的頻散現(xiàn)象； ②當(dāng)ε=0.010%時(shí)，兩種方法的波場快照基本一致，相對(duì)偏差為2%； ③當(dāng)ε=0.001%時(shí)，兩種方法的相對(duì)偏差僅為0.2%。

圖4 不同迭代誤差條件下MLSCG法(左)、LU法(中)波場快照及其對(duì)應(yīng)的殘差對(duì)比(右)

同時(shí)，為了進(jìn)一步考察算法的計(jì)算效率，以及為MLSCG算法最大相對(duì)誤差參數(shù)的設(shè)置提供參考，將MLSCG迭代法最大相對(duì)誤差分別設(shè)置為0.1000%、0.0100%、0.0010%、0.0001%，與 LU算法(相對(duì)誤差為0)的計(jì)算時(shí)間對(duì)比如圖5所示。可以看出：ε=0.0100%時(shí)的計(jì)算時(shí)間是ε=0.1000%時(shí)的 6倍；ε=0.0010%時(shí)的計(jì)算時(shí)間是ε=0.0100%時(shí)的 1.6倍；ε=0.0001%計(jì)算時(shí)間僅是ε=0.0010%時(shí)的 1.4倍； 即使ε=0.0001%時(shí)所用時(shí)間也只是LU方法的48.72%。

應(yīng)用層狀模型對(duì)算法的穩(wěn)定性進(jìn)行分析。當(dāng)相對(duì)誤差分別達(dá)到10%、1%、0.1%、0.01%、0.001%、0.0001%時(shí)，數(shù)值模擬所需迭代次數(shù)如圖6所示。由圖可見，在經(jīng)過530次迭代后，相對(duì)誤差迅速地從10.0% 降到0.1%，之后收斂速度放緩，當(dāng)相對(duì)誤差達(dá)到0.0001%時(shí)，則需要近7500次迭代。圖7為BiCGstab法和本文算法相對(duì)誤差隨迭代次數(shù)變化曲線的對(duì)比，其中兩種迭代方法初值均設(shè)為零，并采用雅各比預(yù)處理。可以看出：本文算法在前20次迭代收斂迅速，隨后收斂速度逐漸變慢，在500次迭代過后相對(duì)誤差接近0.1%，而且整體的曲線變化較為平穩(wěn)，迭代誤差基本隨著迭代次數(shù)增加而減少； BiCGstab算法的收斂曲線波動(dòng)很大，尤其是前300次迭代，基本不收斂，直到400～500次迭代相對(duì)誤差才逐漸收斂到10%左右。值得注意的是，若將誤差設(shè)置為0.0100%，本文算法需要進(jìn)行3200次迭代，而BiCGstab算法只需2000次，因此BiCGstab算法的后續(xù)收斂速度較快，但穩(wěn)定性難以保證，如在層狀模型試驗(yàn)中有相當(dāng)一部分的頻率分量使用BiCGstab算法時(shí)并不能收斂。

圖6 相對(duì)誤差隨迭代次數(shù)的變化曲線

圖7 兩種迭代算法收斂曲線對(duì)比

2.2 復(fù)雜介質(zhì)模型

為進(jìn)一步驗(yàn)證本文頻域模擬方法的適用性，選取部分Marmousi-1模型(圖8)進(jìn)行試算。模型網(wǎng)格數(shù)為300×300，網(wǎng)格尺寸為5m×5m； 震源是主頻為30Hz的Ricker子波位于(750m，5m)； 檢波器均勻分布在模型表面，間隔為5m； 四周同樣采用PML邊界，厚度為50層； 記錄長度為1.5s。圖9為本文算法和LU算法模擬的單炮記錄，可以看出，本文算法的單炮記錄與LU方法基本一致，說明本文算法對(duì)于復(fù)雜模型仍然有較高的精度。

圖8 Marmousi-1速度模型內(nèi)部黑色矩形框代表實(shí)際數(shù)值計(jì)算區(qū)域，下三角表示激發(fā)點(diǎn)位置，上三角表示檢波點(diǎn)位置

圖9 部分Marmousi-1模型兩種算法模擬的單炮記錄對(duì)比

為了分析模型的網(wǎng)格數(shù)對(duì)頻域數(shù)值模擬的影響，對(duì)Marmousi-1模型進(jìn)行重采樣。原始模型網(wǎng)格數(shù)為737×751，記為模型1； 對(duì)其橫、縱向都進(jìn)行1/4重采樣，網(wǎng)格數(shù)變?yōu)?84×187，記為模型2； 橫、縱向都進(jìn)行1/2重采樣，網(wǎng)格數(shù)變?yōu)?68×375，記為模型3； 網(wǎng)格尺寸都設(shè)置為10m×10m，PML邊界厚度設(shè)置為20層，除正演時(shí)長根據(jù)模型大小等比例延長外，其余參數(shù)設(shè)置相同。設(shè)置ε=0.01%，針對(duì)30Hz單頻波場，三個(gè)模型的模擬所占內(nèi)存及計(jì)算時(shí)間如表3所示。由表可見：①隨著模型網(wǎng)格數(shù)的增大，LU算法和MLSCG迭代算法所占用內(nèi)存均呈線性增長，但MLSCG法占用內(nèi)存相較LU減少了15%～20%，且模型越大，內(nèi)存減少越多； ②模型網(wǎng)格數(shù)增大4倍，LU方法計(jì)算時(shí)間會(huì)增加接近15倍，而MLSCG迭代算法由于其良好的收斂性，增加的計(jì)算時(shí)間并不多，相比LU方法，對(duì)于模型3計(jì)算時(shí)間減少了70%以上，對(duì)于模型1甚至減少了87.98%。

表3 三個(gè)模型、兩種方法30Hz波場模擬內(nèi)存占用和用時(shí)統(tǒng)計(jì)

3 結(jié)束語

本文提出了一種基于混合共軛梯度迭代(MLSCG)的頻率域地震數(shù)值模擬方法。該方法使用最小二乘方法優(yōu)化共軛梯度參數(shù)，提高了頻率域正演的計(jì)算穩(wěn)定性。相比傳統(tǒng)直接求解算法，該方法內(nèi)存占用更小、計(jì)算速度更快； 相比常規(guī)迭代算法，只需要簡單的預(yù)處理，不僅能保持較快的收斂速度，還有更高的計(jì)算穩(wěn)定性。