劉立民,劉定進(jìn),李 博

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

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

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

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

1 方法原理

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

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

(1)

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

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

(2)

初始條件:

(3)

(4)

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

(5)

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

(6)

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

(7)

(8)

式中:Me,Ce和Ke分別為單元質(zhì)量矩陣、單元阻尼矩陣和單元?jiǎng)偠染仃嚒Ｍㄟ^單元?jiǎng)偠染仃嚨暮铣?整體平衡方程為:

(9)

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

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

(10)

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

(11)

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

(12)

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

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

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

(13)

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

(14)

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

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

圖2 譜元法波場(chǎng)模擬求解域分塊方案示意

(15)

(16)

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

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

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

(17)

參照?qǐng)D3,分塊內(nèi)的節(jié)點(diǎn)總數(shù)為:

(18)

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

(19)

(20)

其中有:

(21)

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

(22)

(23)

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

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

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

2 性能分析

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

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

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

表1 算法測(cè)試使用集群條件

2.1 收斂性分析

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

2.2 可擴(kuò)展性及并行效率分析

在擴(kuò)展性方面,測(cè)試以下兩類情況。

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

圖6 圖5算例中相對(duì)殘差變化曲線

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

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

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

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

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

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

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

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

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

有限差分法逆時(shí)偏移及譜元法逆時(shí)偏移的成像結(jié)果對(duì)比如圖11所示。有限差分法逆時(shí)偏移由于采用了雙程波方程描述波場(chǎng)傳播,沒有傾角限制,可以實(shí)現(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)造成像時(shí),達(dá)不到理想的成像效果。

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

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

圖11 起伏地表模型偏移成像對(duì)比

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

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

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

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

圖12 SEAM三維起伏地表模型

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

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

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

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

4 結(jié)論

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

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

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