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

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

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

0 引言

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

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

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

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

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

1 方法原理

1.1 波動方程的有限差分離散

首先，將二維時間域聲波方程進行傅里葉變換，得到頻率域波動方程

(1)

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

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

(2)

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

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

(3)

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

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

Ap=b

(4)

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

1.2 混合共軛梯度迭代算法

計算數(shù)學中對大型不定方程組的求解，一般采用Krylov子空間迭代法，這類算法實際上是共軛梯度算法的推廣。算法從初始值p0開始，迭代計算更新pk(其中k為迭代次數(shù))，直到剩余誤差|b-Apk|達到足夠小時求得近似解。若沒有預(yù)處理器，這種方法收斂緩慢或根本不收斂[24]。為了提高計算穩(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)化計算性能，本文給出一種混合算法，即MLSCG，有效結(jié)合了式(5)和式(7)各自的特性，新的共軛梯度參數(shù)為

(8)

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

(9)

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

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

(12)

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

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

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

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

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

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

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

LU分解法和MLSCG算法的實際內(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)勢較為明顯。在模型實例中會有詳細的兩種方法的實際計算時間和內(nèi)存消耗對比。

2 數(shù)值算例

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

2.1 層狀介質(zhì)模型

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

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

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

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

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

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

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

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

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

圖7 兩種迭代算法收斂曲線對比

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

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

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

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

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

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

3 結(jié)束語

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