起伏地表直接成像技術研究進展

曲英銘

(1.中國石油大學(華東)地球科學與技術學院地球物理系,山東青島266580;2.中國石油化工股份有限公司地球物理重點實驗室,江蘇南京211103)

隨著我國油氣資源的日益短缺以及東部油氣勘探程度的不斷提高,地震勘探的重點逐步轉向了勘探難度更大的西部與南方。這些地區起伏劇烈的地表及復雜的近地表結構給地震勘探工作帶來了極大的困難與挑戰[1]:由于地表劇烈起伏、近地表橫向變速劇烈,地表一致性假設不再滿足,高程靜校正后的成像結果會發生畸變,成像精度無法滿足地質解釋的需求。

起伏地表的地震資料處理經歷了以下發展過程:常規靜校正、改進的高精度靜校正、基準面校正[2-4]以及起伏地表直接成像方法。其中基于起伏地表構造直接進行成像是解決復雜近地表地區成像問題的根本方法。

起伏地表直接成像方法總體上包括以下3種:射線類方法[5-6]、單程波方法[7-9]及雙程波方法[10-11]。射線類成像方法對復雜近地表條件具有較好的適應性,可直接在起伏地表上進行波場延拓與偏移成像[12-14]。該方法主要分Kirchhoff偏移法和束偏移法兩種。Kirchhoff偏移法由繞射掃描疊加方法衍生而來,基于地震記錄的加權繞射疊加,應用波動方程的積分解進行地震波場的反向傳播與成像[12]。為了壓制Kirchhoff偏移方法產生的偏移噪聲,GRAY[13]將束偏移引入到Kirchhoff方法中,既保留了Kirchhoff方法計算速度快、對復雜地表適應性強的優點,又提高了偏移剖面的信噪比。岳玉波等[15]基于局部平面波分解,提出了一種起伏地表保幅高斯束成像方法,不僅能夠克服起伏地表的影響,還能得到較為保幅的地下構造成像結果。楊繼東等[16]將改進的高斯束-菲涅爾束應用到起伏地表束偏移方法中,提高了起伏地表平面的分解精度,均衡了淺、中、深層的成像精度。黃建平等[17]將岳玉波等[15]提出的起伏地表保幅高斯束成像方法發展到起伏地表彈性波。起伏地表單程波成像方法主要包括零速層法[18]、逐步累積法[19]、直接下延法[20]和波場上延法[21]等。直接下延法和波場上延法是由零速層法與逐步累積法相結合發展而來的。學者們對波場延拓算子做了許多改進:程玖兵等[22]提出了基于優化系數傍軸近似方程的頻率-空間域偏移算子,該算子能較好地適用于橫向變速情況下的復雜近地表構造成像;王成祥等[23]提出了一種混合法起伏地表單程波偏移方法,提高了起伏地表條件下的成像精度;葉月明等[24-25]在頻率-空間域實現了帶誤差補償的起伏地表有限差分單程波深度偏移方法,壓制了常規頻率-空間域有限差分單程波偏移方法產生的偏移噪聲,并引入保幅算子,提出了雙復雜構造單程波保幅偏移方法。基于雙程波的逆時偏移方法[26]對于高陡構造等復雜地質體成像精度高,近年來隨著計算機技術的發展已成為學者們研究的熱點。關于逆時偏移的研究主要針對提高計算效率[27-28]、優化成像條件[29]、壓制成像噪聲[30-33]以及起伏地表條件下的改善成像質量[34-36]等。

本文對起伏地表直接成像技術的研究進展進行綜述,主要針對射線類、單程波及雙程波起伏地表疊前深度偏移方法。

1 基于射線理論的起伏地表偏移成像

1984年,WIGGINS[12]提出了一種基于起伏地表地震數據的Kirchhoff積分公式及偏移方法,GRAY[13]驗證了Kirchhoff起伏地表直接成像方法的成像優勢,JAGER等[14]優化并發展了起伏地表Kirchhoff成像方法,提出了保幅Kirchhoff起伏地表直接成像方法。Kirchhoff成像方法計算效率高,能夠較好地處理起伏地表,因此也是目前工業界常用的偏移方法。但該方法對復雜構造成像精度較低,且存在嚴重的偏移噪聲。

束偏移方法實際上是一種改進的Kirchhoff成像方法。該方法通過傾斜疊加將局部區域地震數據分解為不同方向的波數并進行成像[37-41]。2005年,GRAY[13]提出了一種基于局部靜校正的起伏地表高斯束成像方法(圖1),但當地表起伏劇烈時,靜校正會導致地震波場產生畸變,影響地下構造成像的精度。針對這個問題,岳玉波等[15]提出了針對起伏地表的保幅高斯束成像方法,將高斯束格林函數(圖2a)與起伏地表條件下利用瑞雷積分計算的逆時傳播的地震波場(圖2b)結合起來,推導并實現了起伏地表條件下的高斯束逆時延拓算子,利用反褶積成像條件實現了起伏地表保幅高斯束成像。該方法的優點是在復雜的近地表處直接進行局部平面波分解,消除了幾何擴散對地震波振幅的影響,從而消除了劇烈起伏地表的影響,且對成像結果進行了保幅。圖3對比了采用岳玉波等[15]提出的起伏地表保幅高斯束成像方法和GRAY[13]提出的基于局部靜校正的起伏地表高斯束成像方法對Canadian Foothills模型成像的結果。Canadian Foothills模型存在起伏劇烈的地表和復雜的近地表速度變化,地震成像極為困難。由圖3可見,相對基于局部靜校正的起伏地表高斯束成像結果(圖3c),基于保幅延拓法的起伏地表高斯束成像結果(圖3b)成像精度和信噪比更高。

采用起伏地表保幅高斯束成像方法對A探區起伏地表實際資料(圖4)進行了偏移處理。圖4a是該探區的速度場,圖4b是某一炮的單炮記錄,圖4c是起伏地表的高程變化,該地區最大高程為58m,局部變化劇烈。實際資料共計329炮,最大道數為320道,采樣點數為1501,采樣率為4ms,道間距為50m,炮間距為100m,CDP間距為25m;速度場橫向1297個CDP,縱向深度7560m,深度采樣間隔是4m。圖5為A探區起伏地表實際資料保幅高斯束偏移成像結果,可以看出起伏地表的影響得到了很好的消除。

圖1 局部靜校正法[13]

圖2 高斯束格林函數(a)和復雜地表條件下的反向延拓波場(b)[15]

圖3 Canadian Foothills模型及疊前深度偏移結果對比a Canadian Foothills模型; b 基于保幅延拓法的起伏地表高斯束成像; c 基于局部靜校正的起伏地表高斯束成像[13]

圖4 A探區起伏地表實際資料a 速度場; b 單炮記錄; c 起伏地表高程

圖5 A探區起伏地表實際資料保幅高斯束偏移成像結果

與波動方程類偏移方法相比,基于射線理論的偏移方法本身就能夠很好地適應起伏地表構造,因此,各種新的射線類偏移成像方法都可以發展為起伏地表直接成像方法。隨著彈性波多分量起伏地表高斯束成像方法的發展,射線類起伏地表直接成像方法已經較為成熟。此外,射線類成像方法具有計算速度快的優勢,因此,射線類起伏地表直接成像方法在工業界得到了廣泛應用。但是,該方法基于射線理論,主要反映地下構造的運動學特征,難以準確地對地下復雜構造的動力學特征進行刻畫,且該方法無法對高陡構造等復雜構造進行成像。

2 基于單程波方程的起伏地表保幅疊前深度偏移

隨著波動方程疊前深度偏移技術的逐漸完善以及在復雜構造地區的廣泛應用,波動方程偏移已經被證實是復雜構造成像較為理想的途徑。相對Kirchhoff偏移方法,波動方程疊前深度偏移方法在構造成像方面有著明顯的優勢。基于波動方程的保幅偏移成像分兩步:第一步是波場傳播,第二步是建立共成像點道集的成像原則。SAVA等[42]基于波動方程偏移建立了角度域共成像點道集(CIGs),并建立了一種保幅成像原則。朱緒峰[43]對保幅單程波傳播理論進行了研究。

基準面校正是地震數據處理必不可少的一步,在近地表橫向速度變化劇烈以及地表起伏地區,基準面校正更加重要。鑒于常規高程基準面校正無法適應起伏地表和近地表速度橫向變化劇烈地區,BERRYHILL[44]提出了針對疊后數據的波動方程基準面校正方法,并將該方法發展到疊前數據[45]。波動方程波場外推技術可以將野外地震數據從地面延拓到任一平面,從而消除了地形劇烈變化對構造成像的影響。BEASLEY等[18]提出“零速層”方法,將水平基準面挪到了地表最高點或最高點以上,在該水平基準面和實際地表之間插入虛擬層,虛擬層的速度用常速度(零速度或非常小的速度)填充,然后從水平基準面開始進行偏移,這期間只考慮地震波的垂向傳播,遇到實際地層再恢復正常。該方法的優勢在于無需做高程靜校正,只對速度模型做微小改變,巧妙地化解了起伏地表的影響。1991年,RESHEF[19]提出了“逐步-累加”的波場延拓思想。近年來,學者們基于BEASLEY等[18]和RESHEF[19]的方法提出了“波場上延”法和“直接下延”法,有效地壓制了地表劇烈變化的影響。如2006年,田文輝等[2]提出了復雜地表和地下地質條件下的直接下延波動方程疊前深度偏移方法;葉月明等[24-25]應用直接下延法研究雙復雜介質條件下的疊前深度偏移,在保持較高計算效率的同時,有效地改善了雙復雜介質的成像效果。

常規傅里葉有限差分法(FFD)在速度變化很大的情況下不穩定。為了克服這一缺點,MILLINAZZO等[3]對基于Pade逼近得到的平方根算子做復Pade逼近,推導出基于復Pade逼近的傅里葉有限差分算子,通過旋轉平方根的分支截斷得到較好的成像效果。該方法的優點在于減少了偏移噪聲,同時對延拓步長的要求不高,提高了計算效率。圖6對比了該方法和常規FFD方法用于2D-SEG起伏地表模型(圖6a)的成像結果[20],相對常規FFD(圖6b),基于復Pade逼近的FFD方法(圖6c)不僅壓制了偏移噪聲,而且很好地改善了深層成像效果。

采用基于復Pade逼近的FFD方法對A探區起伏地表實際資料(圖4)進行了偏移處理測試。圖7a至圖7c分別是基于常規FFD算子、旋轉角度α=5°的復Pade逼近FFD算子和復Pade逼近的保幅FFD算子,采用直接下延法進行疊前深度偏移的成像結果。可以看出:基于復Pade逼近的FFD算子的成像效果要好于常規FFD算子的成像效果,偏移噪聲得到了明顯壓制;基于復Pade逼近的保幅FFD算子的成像效果要好于常規FFD算子的成像效果和基于復Pade逼近的FFD算子的成像效果,在壓制偏移噪聲的同時,有效保持了深層反射能量,淺中層和深層的構造形態都得到了較好的成像。

圖6 2D-SEG起伏地表模型疊前深度偏移a 2D-SEG速度模型; b 常規FFD偏移方法; c 基于復Pade逼近的FFD偏移方法 [20]

圖7 A探區起伏地表實際資料疊前深度偏移a 常規FFD偏移方法; b 基于復Pade逼近的FFD方法; c 基于復Pade逼近的保幅FFD方法

波動方程成像方法相對射線理論成像方法能夠更好地對復雜地下構造進行成像。在計算機技術無法承受逆時偏移巨大的計算量和內存壓力之前,單程波成像方法得到了廣泛應用。但基于單程波方程的成像方法對起伏地表的適應性差,需要對起伏地表進行特殊處理。當學者們認識到起伏地表直接成像的必要性之后,基于單程波方程的起伏地表直接成像方法發展迅速,出現了很多用于消除劇烈起伏地表影響的單程波成像方法。隨著計算機技術的快速發展及逆時偏移技術的推廣應用,近些年基于單程波方程的起伏地表直接成像方法發展較為緩慢,這是因為單程波成像方法在成像精度上不如逆時偏移方法,在計算效率和對起伏地表的適應性上不如射線類方法。

3 基于雙程波方程的起伏地表疊前逆時偏移

3.1 雙程波波場延拓算子

逆時偏移跟正演密不可分,正演技術的發展一定程度上帶動了逆時偏移的發展。關于起伏地表對地震波場的影響研究[46]可以追溯到20世紀40年代,只是因為當時技術不發達,這項工作一直沒有得到展開。到了20世紀70年代,隨著計算機技術的迅速發展,地震模擬技術也隨之迅速發展,20世紀八、九十年代掀起了起伏地表正演模擬研究的熱潮。正演模擬總體上分為兩類:射線類和波動方程類。射線類方法最早由WIGGINS[12]提出,證明了起伏地表情況下Kirchhoff積分偏移的適應性,但該方法終因基于高頻近似、在速度橫向變化劇烈地區適應性差等缺陷發展緩慢。波動方程類方法是對建立的地質模型進行網格劃分,劃分后的模型由有限個離散點組成。該方法對介質沒有橫向變化的限制,如果網格足夠小,得到的解會非常精確,而且綜合考慮了地震波的運動學和動力學特征,因此得到了廣泛應用。

幾十年來,解決起伏地表問題的波動方程方法基本可以歸為四大類:有限差分法、有限元法、譜元法和邊界元法。其中有限差分法因其良好的模擬精度和計算效率應用最為廣泛,該方法用差分算子代替微分算子,將波動方程和計算區域離散化,通過不斷更新迭代獲得各個離散點上的波場值。20世紀60年代,ALTERMAN等[47]開創性提出了層狀介質中的彈性波有限差分方法,自此有限差分方法在地震勘探領域不斷發展。1986年,VIRIEUX[48]提出了一階速度-應力波動方程交錯網格有限差分法,隨后又實現了SH波、P-SV波的波場模擬,極大地提高了波場模擬精度和計算效率,且沒有增加存儲空間。針對起伏地表頻散嚴重的問題,JASTRAM等[49]采用可變網格的方法進行網格剖分,TESSMER等[50]將坐標變換的思想應用于起伏地表。董良國[51]用縱向坐標變換實現了起伏地表條件下的彈性波數值模擬,之后又基于泰勒展開發展了非規則網格差分方法。褚春雷[52]在前人研究的基礎上采用non-Sibson插值方法進行聲波方程正演模擬以及逆時偏移,取得了不錯的效果。有限元法是變分和剖分插值,該方法的核心部分是有限元網格的生成,目前有各向異性網格、貼體網格、并行網格和自適應網格等。1974年,THOMPSON等[53]系統研究了貼體網格,提出了偏微分數值求解方法,證明了Laplace方程和Poisson方程可以實現保角變換。孫建國等[54]采用自適應網格技術生成網格,并將其應用于復雜地表情況下的地震波數值模擬中。總體來說,有限元法的優點是能夠很好地處理起伏地表問題,缺點是內存要求高,計算量大。譜元法是結合有限元和偽譜法發展起來的一種數值模擬方法。KOMATITSCH等[55]用譜元法實現了起伏地表條件下的二維和三維彈性波數值模擬,并且分析了起伏地表引起的各種干擾波,但是未研究地表劇烈變化條件下的有效反射問題。CHE等[56]通過修改Chebyshev譜元算法實現了起伏地表彈性波數值模擬,證明了該方法計算效率高、收斂速度快,并對產生的干擾波進行了討論,但未研究如何消除干擾波。總體來說,譜元法具有精度高、穩定性好、耗時少等優點,但因其計算量仍然較大未得到廣泛應用。邊界元法的計算量較大,而且其本身的半解析性質難以適應地表劇烈變化的情況,因此其應用受到了限制。

3.2 起伏地表逆時偏移

3.2.1 逆時偏移的基本原理

逆時偏移的研究歷史可以追溯到20世紀80年代。BAYSAL等[57]和WHITMORE[26]率先提出了RTM的思想,證明了RTM的優越性,DAN[58]和MUFTI[59]將這一思想應用到地震勘探領域,然而早期的RTM研究都是針對疊后資料的。后期針對疊前數據的RTM研究可以分為聲波方程RTM和彈性波方程RTM兩類。LESAGE等[60]和JONES等[61]將聲波RTM應用于復雜構造的成像,證明了聲波RTM對于復雜構造成像的優勢。隨著多分量地震采集技術和計算機技術的發展,彈性波RTM也得到了廣泛的研究[62-65]。彈性波RTM主要包括兩種類型:一是標量RTM,該方法對縱橫波地震記錄分別進行RTM成像[66],由于準確分離縱橫波記錄還存在一定的難度,因此該方法的應用受到限制。二是矢量彈性波RTM,將縱橫波波場的分離放在應用成像條件之后[67]。在各向同性介質中,對波場應用亥姆霍茲分解定理可以得到縱橫波分離的波場[68]。另外,在縱橫波成像(PS成像和SP成像)中,需要校正橫波的極性[69-70]。雖然RTM對于復雜構造成像具有很大的優勢,但關于如何提高RTM的成像精度和計算效率問題,還需要進一步研究。近年來,學者們通過改造波動方程或進行波場分解[71-74]、修改成像條件[75-76]和成像后濾波[34]等提高了RTM的成像精度和計算效率。

對于二維各向同性介質,聲波方程可以通過降階表示為:

式中:u表示應力;vx、vz分別表示x方向和z方向質點的速度;ρ表示密度;vP表示縱波速度。其時間二階、空間十階的交錯網格有限差分格式可以表示為:

式中:Δx和Δz分別代表水平方向x、垂直方向z的網格間距;k表示時間上的離散值;U、P、Q分別代表u、vx、vz;Δt表示時間步長;i、j分別代表水平方向x和垂直方向z的離散值;C為有限差分系數。由此可以進行聲波方程交錯網格正演,同時對檢波點處波場進行反傳,根據成像條件完成成像。疊前RTM的成像條件包括三類:一是激發時間成像條件;二是互相關成像條件;三是震源波場與檢波點波場振幅比成像條件。其中互相關成像條件可表示為:

式中:x、y是水平坐標;z為深度坐標;t為時間。彈性波動方程逆時偏移的基本原理與聲波相似,這里不再贅述。

3.2.2 變網格逆時偏移

在起伏地表情況下,采用常規有限差分地震波模擬方法必須取很小的網格間距才能保證計算精度,而這樣會帶來巨大的計算量以及局部過采樣問題。為了解決這一問題,JASTRAM等[77]基于二維聲波方程提出了變網格步長的算法。其后,馬光克等[78]實現了可變網格高階有限差分算法,并將其應用于逆時偏移成像。當空間網格變化時,時間采樣需要滿足精細網格對應的穩定性條件,這會增加偏移計算的時間。為此,TAE-SEOB[79]、HUANG等[80]、李振春等[81]先后實現了基于交錯網格的時空雙變算法(圖8),并進行了虛假反射誤差分析。

在進行變網格RTM時,需要注意的是波場外推過程中局部精細網格和粗糙網格之間的轉移。圖9為變網格RTM流程圖[81]。利用近地表起伏模型進行了聲波RTM測試,分別采用粗網格和雙變網格進行了逆時偏移,得到的成像結果如圖10所示[81]。可見雙變網格RTM的成像效果(圖10c)比粗網格RTM的成像效果(圖10b)好。

圖8 雙變網格原理[81]

圖9 變網格逆時偏移流程[81]

3.2.3 坐標變換的有限差分逆時偏移

相對于變網格方法,曲網格對于起伏地表更具優勢。曲網格地震正演模擬最初只適用于內部界面[82],隨后有學者根據HESTHOLM等[83]推導的自由邊界條件公式,推導出適用于起伏地表的曲網格顯式精確條件,提高了計算精度。如魯雁翔等[84]推導了彈性波正向傳播和逆時傳播曲線網格差分格式,實現了曲線交錯網格有限差分法逆時偏移成像。

近年來,貼體網格因其對起伏地表良好的適應性而受到關注,學者們將貼體網格應用于起伏地表正演模擬,取得了不錯的效果。貼體網格生成的原則是使離散后的網格邊界與地表形態吻合,以避免人為產生的階梯邊界引起的虛假散射[53]。ZHANG等[85]采用DRP/opt MacCormack格式計算空間導數,實現了貼體同位網格上的一階速度-應力方程有限差分正演模擬,但該方法計算量較大。隨后,APPELO等[86]、蘭海強等[87]推導并實現了貼體網格上的二階波動方程正演模擬算法。然而二階位移方程在泊松比較大的介質條件下容易出現不穩定,且該算法不易推廣到高階。CHENG等[88]采用貼體旋轉交錯網格實現了起伏地表條件下的地震波正演模擬,李慶洋等[89]發展了曲線坐標系下基于仿真型有限差分的貼體全交錯網格正演模擬算法。如圖11所示,貼體網格可以通過計算空間到物理空間的坐標變換來獲得[90]。

利用含有高斯山峰、山谷的洼陷模型進行了逆時偏移成像測試(圖12)。圖12a為含有高斯山峰、山谷的洼陷模型,該模型共4層,速度分別為2500、3500、4000、4500m/s,模型大小為301×201個網格點,網格間距為10m,最大高程差可達400m。圖12b 為貼體網格剖分圖,可以看出,貼體網格對復雜界面具有強適應性及正交性。觀測系統:301個檢波器均勻分布在地表以下10m處,道間距10m,第一炮位于(0,20m),炮間距30m,共101炮。計算參數:時間間隔0.6ms,最大記錄時間1.8s,震源為主頻 20Hz的雷克子波。圖12c為聲波逆時偏移成像結果,可以看出,起伏地表的影響得到了很好的消除,洼陷構造得到了較準確的成像。

圖10 逆時偏移成像結果a 近地表起伏模型; b 粗網格; c 雙變網格[81]

3.2.4 三角網格

非結構性網格適用于復雜區域,可以控制網格的疏密,較好地模擬曲界面。目前,非結構性網格有很多種形式,其中不規則三角網格作為一種很典型的非結構性網格形式(圖13)[36]受到了國內外學者的廣泛關注。如GUO等[92]使用三角網格法處理起伏邊界,其本質就是對任意形狀的自由表面邊界進行坐標變換,通過引入不規則網格差分算子,使應力-速度方程能夠適用于任意形狀的邊界,將波動方程的差分計算轉換到新坐標系下進行。孫小東等[36]利用三角網格進行差分離散,并將其應用于復雜地表疊前逆時偏移。

對起伏地表模型(圖14a)進行三角網格剖分,然后做RTM,圖14b和圖14c分別是模擬記錄和逆時偏移的結果。可以看出,模型深部同相軸得到了很好的成像,成像剖面的信噪比也比較高。

圖13 起伏地表三角網格剖分a 凹形地表剖分; b 凸形地表剖分[36]

3.2.5 起伏地表粘聲逆時偏移

實際地下介質存在明顯的粘彈性,對地震波有強烈的吸收衰減作用。消除粘彈性對地震波傳播的影響一直是學者們研究的重要課題。到目前為止,主要發展了兩大類方法:反Q濾波和反Q偏移。其中反Q濾波技術在對地震波的吸收衰減進行補償的同時改善了相位特征,早期應用較為廣泛[93-94]。隨著偏移技術的發展,反Q偏移技術逐漸成為主流。根據所依據的理論不同,反Q偏移可分為基于射線理論、基于單程波理論和基于雙程波理論3種。基于射線理論的反Q偏移主要有帶衰減補償的Kirchhoff偏移方法[95]、變換域衰減補償方法[96]以及疊前深度偏移衰減補償方法[97]。1994年,DAI等[98]率先提出了基于單程波方程的反Q偏移方法,在進行偏移的同時補償吸收衰減的影響。隨后,學者們使用單程波方法進行了深度偏移補償。陳雪等[99]依據單程波方程對粘彈介質進行了地震波場數值模擬,并對正演及照明結果進行了對比分析。隨著逆時偏移的發展,學者們又提出了基于雙程波方程的逆時偏移Q補償算法,提高了成像精度。如LIU等[100]推導了一種適用于起伏地表的粘彈介質Q補償逆時偏移方法。簡單說來,粘聲RTM的基本原理就是在常規聲波RTM正、反向延拓過程中同時進行吸收衰減補償。對于起伏地表粘聲RTM,需要同時校正起伏地表和粘彈性對成像的影響。QU等[91]提出了一種粘聲起伏地表逆時偏移方法,推導了適應起伏地表貼體網格的輔助坐標系粘聲擬微分方程,該方法既能夠校正地下粘彈性引起的振幅衰減及相位畸變,又能夠克服劇烈起伏地表對成像的影響。通過實際資料試處理對起伏地表粘聲逆時偏移方法進行了測試。該資料共1180炮,每炮906個檢波點。偏移速度模型和Q模型如圖15所示,速度范圍為1910～6480m/s,Q的范圍為36～110。震源函數為25Hz的雷克子波,時間采樣間隔為2ms,記錄時間為8s。圖16a 和圖16b分別為粘聲和聲波起伏地表RTM成像結果,可見聲波RTM的深部構造成像能量弱,且信噪比較低,而起伏地表粘聲RTM成像結果經過Q補償,斷層、不整合面和其它構造都得到了很好的成像。圖17為從圖16中水平位置4000m處抽取的波數譜曲線[91],證明了粘聲RTM經過衰減補償,成像結果的能量和分辨率都得到了提高。

圖15 實際資料的速度場(a)和Q模型(b)[91]

圖16 起伏地表逆時偏移成像結果對比a 粘聲; b 聲波[91]

圖17 從圖16水平位置4000m處抽取的波數譜曲線[91]

3.2.6 起伏地表各向異性逆時偏移

不同于粘聲逆時偏移在波場補償過程中需要引入特殊處理壓制指數級增長的高頻噪聲,各向異性逆時偏移的核心是各向異性介質的正演模擬,因此這里只介紹各向異性介質數值模擬技術。各向異性介質也是當前學者們研究的熱點,其中最受關注的主要有水平橫向各向同性介質(HTI)、垂直橫向各向同性介質(VTI)和傾斜橫向各向同性介質(TTI)三類,統稱為TI介質。在各向同性介質中,RTM、LSRTM以及全波形反演(FWI)使用的都是標量縱波方程,但是在TI介質中,由于存在縱橫波耦合問題,簡單的標量縱波方程難以表述地震波的運動學和動力學特征。針對這一問題,國內外學者提出了很多解決方法:TSVANKIN[101]推導了VTI介質qP-qSV波精確相速度公式,ALKHALIFAH[102]對該公式進行了簡化,提出了著名的“聲學近似”公式。簡單說來,各向異性逆時偏移的發展歷程,就是qP波方程的簡化歷程。根據前人的研究成果,其大致可分為兩類:第一類是qP-qSV波耦合波動方程,該方程在傾角急劇變化時極不穩定。為了解決這一難題,FLETCHER等[103]引入非零橫波速度來改進qP-qSV波耦合波動方程,從而保證成像的穩定性,DUVENECK等[104]提出一種新的穩定波動方程。第二類是完全不含偽橫波干擾的TI介質純qP波控制方程。DU等[105]在平方根近似和弱各向異性近似條件下,推導出一種TTI介質時間-波數域純qP波解耦方程,該方程在運動學上較為準確,但是在時空域為擬微分形式,導致數值計算存在困難。為此,SONG等[106]、 LIU等[107]提出了時空域最佳分離近似(OSA)方法,該方法雖然計算精度高,但成本也非常高;PESTANA等[108]提出了快速展開法(REM),該方法在VTI介質逆時偏移成像中得到了廣泛應用;ZHAN等[109]提出了基于REM的TTI介質混合有限差分-偽譜法,該方法需要進行多次傅里葉變換和插值,計算效率很低;CHU等[110]推導了一種不含分數階算子的時空域純qP波近似方程。針對起伏地表條件下的TI介質RTM成像,需要在校正各向異性的同時應用起伏地表成像方法。蘭海強等[87]實現了VTI介質起伏地表地震波數值模擬,李慶洋[111]基于貼體網格實現了起伏地表TTI介質地震波正演模擬。

利用如圖18所示的高斯山峰、山谷模型進行了各向異性介質地震波數值模擬測試。高斯山峰、山谷模型的高程表達式為:

3.2.7 起伏地表最小二乘逆時偏移

常規起伏地表逆時偏移雖然可以得到較正確的地下成像結果,滿足構造成像的需求,但仍然存在如下問題:①偏移噪聲大。Laplacian濾波雖然能有效去除低頻噪聲,但去除不徹底,且引入了高頻噪聲。②偏移剖面上反射同相軸中間能量強、兩側能量較弱,即振幅均衡性不佳。③由于地下照明強度隨深度的增加而減弱,因而RTM結果深部能量較弱,振幅保真性差。采用起伏地表LSRTM算法可以有效解決常規RTM存在的問題。

近些年來,最小二乘偏移因其具有信噪比高、保真性好、分辨率高等優點受到國內外學者的廣泛關注,其中最小二乘逆時偏移的研究最為廣泛,多震源最小二乘逆時偏移可以有效提高計算效率。為了壓制多震源串擾噪聲,國內外學者提出了多種編碼技術,主要包括極性編碼[112-113]、振幅編碼[113]、分頻編碼[114]、平面波編碼[115-116]等。2016年,HOU等[117]提出了基于擴展模型空間的最小二乘逆時偏移,降低了最小二乘偏移對偏移速度的依賴性。隨后,徐凱等[118]和ZHU等[119]分別提出了校正粘彈性和各向異性的最小二乘逆時偏移方法。起伏地表最小二乘逆時偏移需要將聲波方程轉化到曲坐標系下,根據Born近似理論得到迭代公式,通過不斷迭代得到成像結果。

利用3.2.3節圖12a所示含高斯山峰、山谷的洼陷模型進行了起伏地表LSRTM測試。觀測系統和計算參數與3.2.3節相同。采用震源極性編碼方式將101炮地震數據組合成一個超道集,使其計算量相當于單炮情形,從而大大緩解了計算的需求。圖20a為反射系數模型,圖20b為基于相位編碼的起伏地表LSRTM得到的第30次迭代結果。可以看出,地下構造清晰可見,相對起伏地表RTM結果(圖12c)在振幅保真性、均衡性、壓制低頻噪聲等方面都有了較大改善,但由編碼引入的高頻串擾噪聲也清晰可見。圖20c 為起伏地表LSRTM得到的第80次迭代結果,該結果與理論反射系數模型(圖20a)非常接近,相對第30次迭代結果(圖20b) 有效壓制了串擾噪聲。

圖18 高斯山峰、山谷模型速度場(a)及其網格剖分(b)

圖19 高斯山峰、山谷模型0.45s時刻的水平分量波場快照a VTI介質; b TTI介質

圖20 起伏地表LSRTM測試a 反射系數模型; b LSRTM第30次迭代結果; c LSRTM第80次迭代結果

近年來,彈性波最小二乘逆時偏移也是一個熱點方向。為了實現多分量數據成像,JIN等[120]建立了能夠對縱橫波進行同時反演的彈性反演理論。QU等[121]提出了基于波場分離技術的彈性波最小二乘逆時偏移,對縱橫波速度及密度分量進行分別成像,并提出了基于波場分離的起伏地表彈性波最小二乘逆時偏移[122]。

4 起伏地表速度反演方法

速度分析的準確性對于地震成像非常重要,疊前深度偏移對速度分析的要求更高,常規的速度分析方法對于起伏地表條件下的地震成像顯得力不從心,為此,國內外學者研究了針對疊前深度偏移的速度反演方法。主要包括兩類:第一類是起伏地表層析速度反演方法,第二類是起伏地表全波形反演方法。兩種方法均可以獨立地獲取疊前偏移速度場,也可以進行聯合,實現多級優化反演、層析反演用以獲取較高質量的低波數速度信息,波形反演用以恢復速度場中的高波數成分,從而得到優化的疊前偏移速度場。

4.1 起伏地表層析反演

針對起伏地表條件,基于角道集的成像域走時層析是行之有效的速度反演方法。ZHANG等[123]和秦寧等[124-125]研究了基于角道集的層析速度反演方法,給出了起伏地表成像域層析速度反演流程(圖21)[124]。首先利用高斯束偏移方法提取角度域共成像點道集,由SAVA等[42]提出的成像條件將不同偏移張角對應的成像值累加到共成像點道集所對應的角度范圍內即可;然后根據實際的層析目標地質體,選擇最合適的參數化形式,以最低的計算成本得到精確的層析速度。走時殘差的拾取是成像域走時層析速度反演的關鍵步驟,基于角道集能夠獲得精度較高的走時殘差。走時殘差和深度殘差轉換關系式為:

Δt=2sΔzcosαcosβ

(5)

式中:Δt為走時殘差向量;Δz為深度殘差向量;s為成像點處的局部慢度值;α為反射層傾角;β為射線入射角,對應角度域共成像點道集的角度。圖22為走時殘差與深度殘差轉化關系示意圖[124]。另外,在反演過程中,可以通過加入正則化和井約束來提高層析反演的穩定性和精度。圖23a和圖23b為起伏地表模型及其炮集記錄,圖23c和圖23d為層析速度場及深度偏移剖面[126],可以看出,起伏地表模型得到了較為清晰的成像。

圖21 起伏地表成像域層析速度反演流程[124]

圖22 走時殘差與深度殘差轉化關系[124]

4.2 起伏地表全波形反演

在全波形反演過程中,通常采用從低頻到高頻逐步反演的多尺度反演策略。Wiener濾波器能夠將源信號轉化為十分接近目標信號的形式,BOONYASIRIWAT等[126]提出了基于Wiener濾波器的高效時間域全波形反演方法,可以將地震數據的頻率轉化為較低頻率,以實現多尺度反演。考慮到計算量和空間采樣問題,采用變網格可以很好地適應起伏地表的情況,因此,時間域多尺度雙變網格全波形反演可以較好地解決起伏地表的速度反演問題。采用如圖24a 所示的起伏地表模型對時空雙變網格全波形反演方法進行了測試,圖24b和圖24c分別為時空雙變網格FWI和常規粗網格FWI經過50次迭代得到的反演速度[127]。可以看出,常規粗網格FWI的結果在起伏地表處出現了較大誤差,而時空雙變網格FWI的結果較為準確,深部速度界面和起伏地表附近的速度值都得到了正確的反演。與時間域多尺度雙變網格全波形反演類似,在頻率域同樣可以實現多尺度雙變網格全波形反演。

圖23 起伏地表模型層析反演a 起伏地表速度模型; b 炮集記錄; c 層析速度場; d 層析后的深度偏移剖面[126]

圖24 50次迭代的反演速度a 起伏地表模型; b 雙變網格FWI; c 常規粗網格FWI[127]

常規的矩形網格模擬起伏地表時因其階梯狀離散會引起虛假散射和繞射,因此難以準確模擬起伏地表的波場。針對這一問題,我們采用曲網格來進行離散化。垂向曲網格-矩形網格耦合方法是將起伏地表附近的區域網格化為曲網格,而深層區域則網格化為矩形網格,可以較好地模擬起伏地表的波場。然而,垂向曲網格坐標變換技術只能在垂向上進行網格映射[49],使該方法受到嚴格的穩定性條件限制,因此無法穩定地模擬地震波在劇烈起伏地表中的波場。QU等[128]根據地表附近正交曲網格的坐標變換方程[82],提出了曲線系統下的正交曲網格-矩形網格耦合方法,用以反演起伏地表條件下的速度。采用SEG逆掩斷層模型(圖25)對正交曲網絡-矩形網格方法進行了測試,圖26為最終速度反演結果[128]。該反演結果表明,正交曲網格-矩形網格耦合方法可以在復雜近地表情況下較好地反演出速度值。

圖25 SEG逆掩斷層模型a vP; b vS[128]

圖26 SEG逆掩斷層模型反演結果a vP; b vS[128]

5 應用中的瓶頸問題及技術方案

起伏地表射線類方法能夠很好地適應起伏地表,但主要反映地下構造的運動學信息,難以對地下復雜構造的動力學信息進行準確的刻畫,且該方法難以對高陡構造等復雜構造進行準確成像。因此,在研究起伏地表射線類方法時,要盡可能考慮地下構造的復雜性。①圍繞高斯束格林函數構建方式,考慮模型介質水平及縱向速度變化對束形態的影響,發展基于速度驅動的起伏地表自適應高斯束偏移方法。同時考慮淺部及深部構造成像分辨率,引入聲學中的“菲涅爾帶”構建思想,發展起伏地表投影菲涅爾束偏移方法。②考慮復雜地下構造的粘彈性和各向異性,發展起伏地表粘彈性及各向異性射線偏移方法。③引入反演思想,發展起伏地表最小二乘Kirchhoff偏移或高斯束偏移。

起伏地表單程波成像方法在成像精度上不如逆時偏移方法,在計算效率和對起伏地表的適應性上不如射線類方法,因此后續的研究重點是突出兩者的優點,弱化兩者的缺點。

基于雙程波的逆時偏移在近十年來得到了飛速發展。該方法從理論上可以對任意復雜構造進行精確成像,近年來大量的應用實例也證明該方法比射線類和單程波類成像方法得到更好的成像結果,但該方法對劇烈起伏地表的適應性較差,需要進行特殊處理。起伏地表逆時偏移的關鍵是起伏地表波場延拓算子的構建方法,從求解波動方程的角度可將其分為有限差分法、有限元法、偽譜法、譜元法等。其中,有限差分法的突出優勢是計算速度快、占用內存少,但其對劇烈起伏地表的適應性最差。制約起伏地表逆時偏移方法的瓶頸問題主要有兩個:①對速度的依賴性強;②計算量龐大,特別是近些年基于反演的最小二乘逆時偏移,對地下構造粘彈性和各向異性的引入更是給該方法帶來了巨大的計算負擔。起伏地表逆時偏移方法對速度的依賴性強,但在山前帶等劇烈起伏地表、復雜構造區域,精確的速度建模非常困難,為此,地球物理學家們研究了多種對起伏地表進行準確速度反演的方法。然而,任意一種速度建模方法都不能獲得準確的速度,需要將多種速度建模方法配合使用,充分發揮速度分析、層析反演、波形反演等多種速度建模方法的優勢,并采用多信息融合的淺—中—深層、低—中—高頻逐級優化的反演策略,提高起伏地表速度反演的精度和效率。

在計算效率提升方面,可以采取以下策略:①利用各種近似,在保證計算精度的前提下最大限度地降低計算量。比如校正各向異性時,使用弱各向異性假設,校正粘彈性時,只補償對成像影響相對大的振幅衰減的影響,忽略相位頻散的影響。②對起伏地表波場模擬算子進行優化,提高波場模擬的效率。③充分利用GPU/CPU協同加速、MPI/Openmp并行加速及波場區域分解等技術。④在起伏地表最小二乘逆時偏移中,利用編碼技術壓縮地震數據,并發展提高最小二乘逆時偏移收斂速度的反演算法等。

6 結論與展望

地震成像是一門不斷發展、不斷深入的學科。成像技術不斷發展的最終目的是使地下復雜構造精確成像。起伏地表直接偏移技術不僅具有處理靈活、計算效率高的特點,而且能對復雜構造進行精確成像,消除常規偏移方法中的噪聲影響,對于我國西南地區復雜山地、丘陵構造的成像有著重要意義。然而,目前起伏地表直接偏移技術尚處于實驗探索階段,存在很多問題需要解決。地下巖性及儲層精細描述的需求必然會帶動起伏地表直接偏移技術的發展,筆者認為起伏地表直接偏移技術的發展方向有以下兩個方面。

在提高偏移成像精度方面,此前的成像方法都是基于各種假設使模型簡單化,但地下介質受構造運動的影響其性質多種多樣,需要對復雜介質進行不斷探索才能實現起伏地表的精確成像。偏移速度的準確性也會對成像質量有很大的影響,如何通過反演得到準確的速度也是目前需要深入探索的方向。

在提高成像精度的同時,計算量的增加降低了工作效率也是一個需要考慮的實際問題。為了消除起伏地表的影響,對起伏地表進行網格剖分的方式會影響成像質量。雙變網格計算速度快、效率高,但無法消除起伏地表噪聲的影響。曲網格成像質量高,但進行坐標變換無疑會增加計算量,降低工作效率。

綜上所述,如何在提高復雜構造成像效果的同時,保證算法的靈活簡便,降低計算成本,是廣大地球物理學家需要繼續探索的方向。