劉立民,劉定進,李 博

(中國石油化工股份有限公司石油物探技術(shù)研究院,江蘇南京211103)

目前,國內(nèi)外常用的逆時偏移成像算法包括有限差分法和有限元法。有限差分法因?qū)崿F(xiàn)過程簡單且計算高效而被廣泛應(yīng)用,但是,當(dāng)山地地表起伏較大或地質(zhì)構(gòu)造復(fù)雜時,難以確保成像精度,且有限差分逆時偏移的適應(yīng)性和靈活性較差。有限元法具有良好的邊界適應(yīng)性及模擬精度,可以適應(yīng)山地起伏地表等復(fù)雜邊界問題,但面對大規(guī)模復(fù)雜三維問題時,難以在波場模擬精度和運算效率方面取得平衡。國內(nèi)外學(xué)者從地下介質(zhì)特征、波場傳播算子、吸收邊界條件以及并行計算效率等方面對逆時偏移方法進行了大量研究,逐步提高成像精度及運行效率。冀國強等[1]研究基于粘聲介質(zhì)的逆時偏移方法,同時考慮因地層吸收而導(dǎo)致振幅衰減與相位頻散,提出一種正則化形式的穩(wěn)定傳播粘聲逆時偏移方法,提升深層構(gòu)造的分辨率與可信度。曲英銘等[2]提出面向高陡構(gòu)造的粘聲棱柱波逆時偏移方法,改善了高陡構(gòu)造的照明和成像效果。鞏向博等[3]研究基于稀疏約束的最小二乘逆時偏移方法,提高小尺度散射體的成像分辨率。李青陽等[4]提出一種新的卷積完全匹配層(NCPML)吸收邊界條件,用于逆時偏移后,計算效率和內(nèi)存占用上較常規(guī)的完全匹配層(SPML)吸收邊界條件更優(yōu)。段心標(biāo)[5]提出基于GPU的傅里葉有限差分逆時偏移方法,保護了低頻信息的成像且提高了成像效率。

PATERA[6]在研究Navier-Stokes方程的數(shù)值解法時提出譜元法(spectral element method,SEM),該方法不僅具備有限元法處理復(fù)雜構(gòu)造邊界的適應(yīng)性,還有偽譜法高階插值和快速收斂的優(yōu)勢。譜元法的基本思路是:首先,將求解域剖分成相互連接互不重疊的有限個單元;然后,在每個單元上應(yīng)用偽譜法,將單元的近似解表示成不同格式的截斷正交多項式展開式;最后,采用伽遼金法求解得到整體問題的近似解,實現(xiàn)求解域內(nèi)波場傳播的模擬[7-8]。譜元法的核心是采用的插值基函數(shù)形式,目前常見的用于譜元法單元插值的基函數(shù)有切比雪夫(Chebyshev)正交多項式和勒讓德(Legendre)正交多項式。PRIOLO等[9]在研究聲波波場傳播特性時,采用切比雪夫正交多項式進行單元插值,形成了切比雪夫譜元法,之后該方法逐步被應(yīng)用于不同情況下的波場模擬中。MADAY等[10]在譜元法單元分析時采用勒讓德正交多項式作為單元插值基函數(shù),同時結(jié)合GLL(gauss lobatto legendre)積分,將積分節(jié)點設(shè)定為插值節(jié)點,形成了勒讓德譜元法。李孝波等[11]采用勒讓德譜元法進行局部工程地震的波場模擬研究,認為譜元法能夠較好地反映地形條件和近地表特征及斷裂構(gòu)造樣式,具有較高的精度。林偉軍等[12]對復(fù)雜非均勻介質(zhì)下譜元法波場模擬效率進行研究,通過在單元內(nèi)引入獨立的輔助網(wǎng)格,在較稀疏的主網(wǎng)格上進行波場求解,獲得較高的精度。

與傳統(tǒng)方法相比,譜元法雖然在理論上有待進一步完善,但由于結(jié)合了有限元法處理復(fù)雜構(gòu)造的幾何適應(yīng)性和偽譜法高精度快速收斂的特性,目前已成為復(fù)雜地表、復(fù)雜構(gòu)造地震模擬與成像的重要工具[13]。隨著油氣地震勘探的不斷深入,逐步向復(fù)雜山地、復(fù)雜構(gòu)造及復(fù)雜儲層探區(qū)發(fā)展,這些探區(qū)地表起伏大,橫向速度變化劇烈,多處地表高速體出露,地下高速火成巖及低速膏泥巖不規(guī)則穿插,常規(guī)的有限差分逆時偏移成像方法不能完全適應(yīng)起伏地表復(fù)雜斷裂精確成像,無法滿足勘探開發(fā)及井位部署需求,急需研究適應(yīng)起伏地表的復(fù)雜構(gòu)造高精度成像方法。本文以起伏地表復(fù)雜構(gòu)造高精度地震成像為目標(biāo),以起伏地表譜元逆時偏移成像方法為研究內(nèi)容,推導(dǎo)了波動方程切比雪夫譜元逆時偏移成像算法,同時為了提高該方法的并行計算效率,基于自由度凝聚和局部松弛的思想,推導(dǎo)了切比雪夫譜元逆時偏移方法的多級并行算法(HEP-CSE-RTM),運行效率大幅度提高,達到與常規(guī)有限差分法逆時偏移相當(dāng)?shù)挠嬎阈省Ｓ眉幽么笃鸱乇砟嫜跀鄬佣S模型數(shù)據(jù)進行測試,驗證了本文方法對起伏地表探區(qū)復(fù)雜構(gòu)造成像的適應(yīng)性;同時利用SEAM三維起伏地表復(fù)雜構(gòu)造模型數(shù)據(jù)進行了測試,成像結(jié)果表明,該方法在起伏地表復(fù)雜構(gòu)造區(qū)成像精度明顯高于常規(guī)的有限差分法逆時偏移,驗證了該方法的有效性和科學(xué)性。

1 方法原理

1.1 切比雪夫譜元逆時偏移方法(CSE-RTM)

采用切比雪夫譜元方法并結(jié)合隱式Newmark時間積分方法求解波動方程,得到波場傳播算子,進行逆時偏移成像,該解法在空間上具有譜精度,在時間上達到二階精度。均勻介質(zhì)空間Ω中,時間T內(nèi)二維波動方程如下:

(1)

式中:Ω為分析空間R2中有界區(qū)域,Ω≡Γ為分析域Ω的邊界;u(x,z)為聲波位移場;c0(x,z)為聲波速度場;x,z和t分別為空間和時間坐標(biāo);f為震源項。

時間T內(nèi),Ω空間邊界Γ上Clayton-Engquist吸收邊界條件(CE-ABC)[14]:

(2)

初始條件:

(3)

(4)

式中:s為求解域Ω邊界Γ的切線方向;v為H1(Ω)空間中的任意函數(shù)。

(5)

在一個標(biāo)準(zhǔn)單元內(nèi),選取切比雪夫多項式的極值點(Chebyshev-Gauss-Lobatto點)為插值點,使用高階切比雪夫正交多項式對試函數(shù)進行逼近,標(biāo)準(zhǔn)單元內(nèi)的試函數(shù)為:

(6)

插值函數(shù)表示為:

(7)

(8)

式中:Me,Ce和Ke分別為單元質(zhì)量矩陣、單元阻尼矩陣和單元剛度矩陣。通過單元剛度矩陣的合成,整體平衡方程為:

(9)

圖1 局部坐標(biāo)與標(biāo)準(zhǔn)坐標(biāo)的相互轉(zhuǎn)換

考慮到算法的穩(wěn)定性,在時間域處理時,采用隱式Newmark時間積分方案,在一個時間步長t～t+Δt內(nèi)采用下列假設(shè):

(10)

式中:α和δ是由數(shù)值精度和穩(wěn)定性要求決定的參數(shù),δ稱為Newmark因子,α和δ取不同數(shù)值,代表不同差分方案;當(dāng)α≥0.25(0.5+δ)2且δ≥0.5時是絕對穩(wěn)定的[15-16]。因此,常微分方程組(9)可以離散為如下形式:

(11)

其中初始值可以通過如下方程求得:

(12)

關(guān)于Newmark積分方法有兩點值得注意:一是該積分方法有多種表示形式,限于篇幅,本文不做贅述;二是盡管該積分方法是絕對穩(wěn)定的,但是時間步長也要足夠小,因為時間步長是波場模擬及成像精度的決定性因素之一。

1.2 譜元逆時偏移多級并行算法(HEP-CSE-RTM)

在譜元法波場模擬及成像中,通過結(jié)構(gòu)網(wǎng)格空間離散及經(jīng)典Newmark方法時間離散將偏微分方程化成如下代數(shù)方程組:

(13)

對一個超大規(guī)模的譜元法波場模擬及成像問題,假設(shè)其中有nd個節(jié)點,每個節(jié)點的自由度數(shù)(NDOF)是s,這樣整體的自由度數(shù)是s·nd,公式(13)中的kij是整體剛度矩陣,si和dj分別是載荷向量和待求的位移向量。在譜元法波場模擬及成像問題中矩陣kij通常是對稱正定的,這樣線性方程組的求解等價于尋找對應(yīng)泛函極小時的dj,泛函形式如下:

(14)

根據(jù)最小勢能原理,若能設(shè)計一種迭代格式使得泛函Π的值不斷降低,那么必然能夠得到方程組(13)的解。由此思路出發(fā)設(shè)計譜元法逆時偏移多級并行迭代算法。

為了進行大規(guī)模并行計算,首先需將求解域劃分成塊(Set),把nd個節(jié)點按照初始位置分成B個分塊(對應(yīng)B個處理器),各分塊內(nèi)的自由度約為s·nd/B。擬采取的分割方法與域分解思路不同處是,將分割面放在單元內(nèi),而不是節(jié)點上(圖2)。分割面穿過的單元定義為界面單元,相鄰分塊共享界面單元,以此保證位移連續(xù)條件,該方法與域分解方法的重要不同點是不需要再引入罰函數(shù)或拉格朗日乘子,但方程組迭代求解過程中,分塊需要與相鄰分塊通信獲取界面單元上外部節(jié)點的解。

圖2 譜元法波場模擬求解域分塊方案示意

(15)

(16)

以二維起伏地表模型譜元法波場模擬及成像問題為例,對于起伏地表而言,三角形網(wǎng)格具有更好的適應(yīng)性(圖3),在形狀較為規(guī)則的區(qū)域,為保持一致,我們同樣使用三角形網(wǎng)格,相當(dāng)于在矩形網(wǎng)格圖形中心位置增加一個離散點,形成1-3-2、2-3-5、5-3-4和4-3-1(平面內(nèi)按逆時針排列)4個三角形單元,構(gòu)建一個滿足坐標(biāo)二次函數(shù)的物理量分布,實現(xiàn)較高的計算精度。

于是,每個分塊內(nèi)的網(wǎng)格數(shù)量為:

圖3 網(wǎng)格與三角形單元對應(yīng)示意

(17)

參照圖3,分塊內(nèi)的節(jié)點總數(shù)為:

(18)

在求解域的每個分塊內(nèi),使用位移增量模式的線性組合來表示節(jié)點位移值,如第I分塊Set(I)內(nèi)節(jié)點的位移值可寫作:

(19)

(20)

其中有:

(21)

分別是與高階自由度對應(yīng)的剛度矩陣和非平衡力向量。得到凝聚自由度的整體平衡方程為:

(22)

(23)

低于要求的精度ε便得到方程組(13)的解。其中,‖·‖表示對向量取模運算。

由公式(19)至公式(22)可知,本算法不需要組合整體剛度陣。但考慮到僅使用很少的位移增量模式表示精細自由度,相當(dāng)于引入了附加限制,因而迭代中需要通過分塊內(nèi)松弛計算自適應(yīng)地調(diào)整位移增量模式。譜元法中節(jié)點只與周圍數(shù)目非常少的節(jié)點有聯(lián)系,得到剛度陣通常是帶狀分布和稀疏的,因此,圖2 中的每一個分塊僅需和鄰近分塊交換界面信息,產(chǎn)生的通信絕大部分都是局部的,避免出現(xiàn)大規(guī)模的整體信息交換情況,通信負擔(dān)低。一個典型二級結(jié)構(gòu)迭代方案的求解流程如圖4所示,編程中可以使用MPI非阻塞通信實現(xiàn)通信與計算的重疊,減少處理器等待時間,圖中虛線箭頭表示通信數(shù)據(jù)流向。上述方法還可通過層層自由度凝聚,建立多級平衡求解格式,實現(xiàn)更大規(guī)模問題模擬。目前使用的位移增量模式包含常規(guī)位移增量模式和自適應(yīng)位移增量模式兩類,其中常規(guī)位移增量模式用于捕捉分塊的整體運動變形趨勢;自適應(yīng)松弛位移增量模式用于捕捉分塊內(nèi)不均勻變形,放在各分塊所用的從處理器進行計算。

圖4 譜元法二級自由度凝聚迭代方案流程

2 性能分析

選擇加拿大起伏地表逆掩斷層速度模型(圖5)作為算例進行性能分析,該模型橫向24930m,縱向10000m,空間采樣為15m×15m,包含落差達1500m的起伏地表、高速體及低速體同時出露地表、眾多高陡傾角斷層以及逆掩推覆構(gòu)造。在其地表高點進行切比雪夫譜元法單炮正演模擬來考察本文算法的收斂性、可擴展性及并行效率,時間步長為1ms,模擬長度為10s。

圖5 加拿大起伏地表逆掩斷層速度模型

測試所用集群環(huán)境見表1,最多使用了2001個處理器,求解有10×108個未知數(shù)的波場傳播問題。

表1 算法測試使用集群條件

2.1 收斂性分析

首先設(shè)計簡單算例I分析算法的收斂性。采用的處理器數(shù)目是2400個核。將圖5模型剖分成具有10×108個節(jié)點規(guī)模的單元,系統(tǒng)共有10×108個未知數(shù),分成2000塊,每塊有50×104個自由度,共使用了2001個處理器,進行單炮的正演模擬,得到該問題相對殘差隨迭代變化曲線(圖6),可知每兩個殘差下降一個數(shù)量級,收斂速度保持穩(wěn)定。

2.2 可擴展性及并行效率分析

在擴展性方面,測試以下兩類情況。

首先,固定求解域分塊數(shù)目,增加分塊內(nèi)自由度數(shù)。將圖5所示模型分為8×8共64個分塊,不斷增加分塊內(nèi)的未知數(shù)數(shù)量,結(jié)果如圖7所示,分塊內(nèi)的未知數(shù)從400個增長到了980000個,可以看出,無論是每個分塊的未知數(shù)數(shù)量還是整個問題的規(guī)模都發(fā)生了巨大的變化,但收斂所需迭代次數(shù)達到28次,就不再增長。

圖6 圖5算例中相對殘差變化曲線

其次,固定每個分塊內(nèi)的自由度數(shù),在每個分塊內(nèi)放置50×104個自由度,像砌墻一樣不斷增加分塊的數(shù)目,由2個變化到2000個(10×108個未知數(shù)),從圖8中可以看出,當(dāng)問題規(guī)模足夠大時,收斂所需迭代步數(shù)也不再變化,同時,每步迭代所需要的平均計算時間的增長速度也遠遠低于分塊增長速度。

圖7 收斂所需迭代次數(shù)與分塊內(nèi)未知數(shù)關(guān)系

綜合迭代次數(shù)和迭代時間說明,本文提出的多級并行切比雪夫譜元逆時偏移方法有很好的擴展性,借助集群系統(tǒng)能夠有效提高成像問題大規(guī)模并行計算效率。

最后,在并行效率方面,使用800×104個未知數(shù)和3200×104未知數(shù)算例測試了譜元逆時偏移多級并行算法的并行加速比狀況(圖9),本文提出的算法具有很好的并行效率,800×104的算例只有使用超過625個進程后,加速比才出現(xiàn)下降的情況,且對比兩個算例可以看出,規(guī)模越大,加速比下降的位置出現(xiàn)越晚。

圖8 收斂所需迭代次數(shù)(a)、迭代步平均計算時間(b)與處理器數(shù)目關(guān)系

圖9 超大規(guī)模稀疏線性方程組并行算法加速比曲線

3 模型數(shù)據(jù)測試

3.1 加拿大起伏地表逆掩斷層二維模型

為了驗證譜元逆時偏移多級并行算法在解決起伏地表復(fù)雜構(gòu)造精確成像上的能力及效率,選取圖5所示模型及數(shù)據(jù)包,該數(shù)據(jù)包采用480道中間放炮對稱觀測系統(tǒng)進行正演模擬,共采集278炮單炮數(shù)據(jù),道距為15m,炮點距為90m,最小炮檢距為15m,最大炮檢距為3600m,每道4ms采樣,道長為5s(圖10)。單元內(nèi)采用4階切比雪夫多項式插值,模型按照4×4分塊方案調(diào)用17個計算節(jié)點并行計算。

有限差分法逆時偏移及譜元法逆時偏移的成像結(jié)果對比如圖11所示。有限差分法逆時偏移由于采用了雙程波方程描述波場傳播,沒有傾角限制,可以實現(xiàn)回轉(zhuǎn)波成像,在近地表復(fù)雜構(gòu)造、中深層復(fù)雜構(gòu)造及逆掩推覆構(gòu)造區(qū)都能成像,特別是在深部逆掩推覆構(gòu)造帶,構(gòu)造邊界刻畫清晰干脆。然而,由于有限差分法網(wǎng)格剖分通常采用矩形規(guī)則網(wǎng)格,其在處理近地表速度突變區(qū)復(fù)雜構(gòu)造成像時,達不到理想的成像效果。

圖10 加拿大起伏地表逆掩斷層正演單炮記錄

譜元逆時偏移方法由于采用譜元法求解波場傳播問題,求解域單元離散方式多樣,可以處理諸如起伏地表及速度突變區(qū)復(fù)雜構(gòu)造成像,其成像結(jié)果相對于有限差分發(fā)逆時偏移精度明顯提高(圖11b)。兩種方法的成像對比結(jié)果表明:譜元逆時偏移方法具有較高的成像精度,適應(yīng)起伏地表及復(fù)雜構(gòu)造,可以精確刻畫斷裂系統(tǒng),具有明顯成像優(yōu)勢。

圖11 起伏地表模型偏移成像對比

從單炮成像及全數(shù)據(jù)體成像兩方面分析有限差分逆時偏移方法(FD-RTM)和譜元逆時偏移多級并行算法的成像耗時,單炮成像有限差分逆時偏移方法耗時7.28min,譜元逆時偏移多級并行算法耗時6.62min;有限差分逆時偏移方法全部278炮整體耗時33.88h,而譜元逆時偏移多級并行算法全部728炮整體耗時30.86h?？梢钥闯?由于譜元法采用切比雪夫高階多項式插值,不論是單炮偏移,還是全數(shù)據(jù)體偏移,譜元法逆時偏移相對于有限差分法逆時偏移耗時均有所減少。

3.2 SEAM三維起伏地表復(fù)雜構(gòu)造模型

為了進一步驗證譜元逆時偏移多級并行算法在解決三維起伏地表復(fù)雜構(gòu)造精確成像上的能力及效率,選取如圖12所示的SEAM三維起伏地表速度模型,該速度模型x方向為12510m,y方向為14460m,z方向為10000m,網(wǎng)格大小為12.5m×12.5m,包含落差達1550m的起伏地表、出露地表的高速山體、三角洲低速沉積層及地下高陡斷層推覆體構(gòu)造。

正演單炮數(shù)據(jù)共25355炮,采用正交對稱觀測系統(tǒng)采集,全排列接收,道距為12.5m,線距為125m,炮點距為25m,炮線距為250m,最小炮檢距為12.5m,橫向最大炮檢距為12500m,縱向最大跑檢距為14500m,每道采樣間隔為4ms,道長為5s。選取圖12極高點附近部分炮集數(shù)據(jù)共1005炮,進行過極高點剖面(圖13)局部三維成像,單元內(nèi)采用4階切比雪夫多項式插值,模型按照8×8分塊方案調(diào)用65個計算節(jié)點并行計算。

圖12 SEAM三維起伏地表模型

有限差分法逆時偏移及譜元法逆時偏移的成像結(jié)果如圖14所示。同樣地,由于有限差分逆時偏移方法三維網(wǎng)格剖分通常采用正六面體等規(guī)則網(wǎng)格,在處理三維近地表速度突變區(qū)復(fù)雜構(gòu)造成像時,采用該方法不能達到理想的成像效果。譜元逆時偏移方法求解域單元離散方式多樣,可以處理諸如三維起伏地表及速度突變區(qū)復(fù)雜構(gòu)造成像,其近地表淺層成像結(jié)果相對于有限差分法逆時偏移精度明顯提高(綠框),同時,深層目標(biāo)區(qū)(藍框)背斜構(gòu)造及其內(nèi)幕成像也有明顯的優(yōu)勢(圖14b)。兩種方法的成像對比結(jié)果表明:譜元逆時偏移方法適用于三維起伏地表及復(fù)雜構(gòu)造成像且具有較高的成像精度。

同樣地,從單炮成像及全數(shù)據(jù)體成像兩方面分析有限差分逆時偏移方法和譜元逆時偏移多級并行算法的成像耗時,單炮成像有限差分逆時偏移方法耗時12.5min,譜元逆時偏移多級并行算法耗時13.8min;有限差分逆時偏移方法整體耗時8.8d,而譜元逆時偏移多級并行算法整體耗時9.7d。可以看出,由于譜元法采用稀疏網(wǎng)格高階多項式插值,以及高擴展性多級并行計算方案,同樣地,不論是單炮偏移,還是全數(shù)據(jù)體偏移,譜元法逆時偏移相對于有限差分法逆時耗時沒有較大的增長,并且通過進一步研究基于GPU集群系統(tǒng)的譜元逆時偏移,該方法的運算效率有望進一步提高。

圖13 SEAM三維起伏地表模型過最高點速度剖面(a)及正演單炮記錄(b)

圖14 SEAM起伏地表模型偏移成像結(jié)果對比

4 結(jié)論

本文結(jié)合絕對穩(wěn)定的隱式Newmark時間積分方法,將切比雪夫譜元法應(yīng)用于起伏地表復(fù)雜構(gòu)造區(qū)地震波場模擬及成像,推導(dǎo)了采用切比雪夫譜元方法模擬波場傳播的逆時偏移算法。

為了提升該方法的并行效率,基于自由度凝聚的思想,推導(dǎo)出針對譜元波場模擬及逆時偏移成像處理的多級并行迭代求解算法,該算法在控制收斂所需迭代次數(shù)的同時,具有并行效率不隨處理器數(shù)目增多而降低的優(yōu)點。由于譜元法既有有限元方法的邊界適應(yīng)性,又具有偽譜法的高精度快速收斂特性,再結(jié)合多級并行計算方案,譜元逆時偏移多級并行算法對起伏地表復(fù)雜構(gòu)造的成像適應(yīng)性及成像精度優(yōu)于常規(guī)的有限差分逆時偏移方法,同時其運算效率也與常規(guī)有限差分法逆時偏移相當(dāng)。

雖然譜元逆時偏移方法具有起伏地表適應(yīng)性強、精度高和收斂快等優(yōu)點,但也存在較低的時間精度與較高的空間譜精度不匹配問題,這將是今后需要重點解決的問題。