起伏地表譜元逆時偏移方法

劉立民,劉定進,李 博

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

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

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

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

1 方法原理

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

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

(1)

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

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

(2)

初始條件:

(3)

(4)

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

(5)

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

(6)

插值函數表示為:

(7)

(8)

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

(9)

圖1 局部坐標與標準坐標的相互轉換

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

(10)

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

(11)

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

(12)

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

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

在譜元法波場模擬及成像中,通過結構網格空間離散及經典Newmark方法時間離散將偏微分方程化成如下代數方程組:

(13)

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

(14)

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

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

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

(15)

(16)

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

于是,每個分塊內的網格數量為:

圖3 網格與三角形單元對應示意

(17)

參照圖3,分塊內的節點總數為:

(18)

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

(19)

(20)

其中有:

(21)

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

(22)

(23)

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

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

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

2 性能分析

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

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

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

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

2.1 收斂性分析

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

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

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

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

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

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

圖7 收斂所需迭代次數與分塊內未知數關系

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

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

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

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

3 模型數據測試

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

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

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

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

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

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

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

3.2 SEAM三維起伏地表復雜構造模型

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

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

圖12 SEAM三維起伏地表模型

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

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

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

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

4 結論

本文結合絕對穩定的隱式Newmark時間積分方法,將切比雪夫譜元法應用于起伏地表復雜構造區地震波場模擬及成像,推導了采用切比雪夫譜元方法模擬波場傳播的逆時偏移算法。

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

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