T分布起伏地表最小二乘逆時偏移

李 慶 洋

(中國石化中原油田物探研究院，河南濮陽 457001)

0 引言

隨著油氣勘探、開發的深入，勘探難度逐漸加大，在復雜構造、復雜巖性、復雜地表條件等地區進行油氣地震勘探面臨諸多難題，其中劇烈起伏的地表已成為制約上述地區地震勘探發展的瓶頸之一[1]。傳統的偏移成像方法大多基于水平地表假設，對起伏地形的適應性不強，容易導致復雜地表條件下的地下成像質量不高[2]。

相比于其他偏移方法，逆時偏移(RTM)應用雙程波動方程延拓波場，避免了對波動方程的近似，無傾角限制，是目前較為精確的偏移成像方法，但仍屬于常規偏移的范疇，其偏移算子是線性正演算子的共軛轉置，而不是它的逆[3]。當地震數據采集不足或不規則、地下構造復雜、波場帶寬有限時，常規偏移方法只能對地下構造模糊成像，無法滿足巖性油氣藏勘探、開發的需求。最小二乘偏移(Least-Squares Migration,LSM)將成像看作最小二乘意義下的反演問題，通過不斷擬合誤差泛函，得到振幅保真性更好、分辨率更高、偏移噪聲更少的成像結果[4-6]。早期的LSM主要基于射線理論[7-8]和單程波理論[9-13]，近年來基于雙程波理論的LSRTM得到廣泛關注[14-19]，但目前的LSRTM算法大多基于水平地表假設，沒有考慮起伏地形的影響，在一定程度上限制了其推廣應用。

正演模擬是LSRTM的組成單元，有限差分法計算效率高、占用內存小且易簡單實現，是當前應用最廣泛的求解雙程波動方程方法，但在模擬復雜地表及地下復雜界面時，在應用笛卡爾坐標系下的規則網格剖分時，易出現階梯狀邊界而形成虛假繞射波。貼體網格能很好地擬合任意起伏地形，宜于處理起伏地表問題[20]。但貼體網格需要坐標變換，而求解變換后的曲坐標系波動方程時，標準交錯網格會引入插值誤差而降低模擬精度。為此，本文將全交錯網格[21]引入貼體網格中，實現了起伏地表有限差分正演模擬，進一步推導了貼體網格線性Born正演方程，在此基礎上提出了基于貼體全交錯網格的起伏地表LSRTM算法，較好地克服了起伏地形的影響。

此外，起伏地形條件的地震數據常含有較多噪聲，常用的基于L2模擬合的LSRTM算法對噪聲非常敏感，尤其是當地震數據中含有奇異值時，常規LSRTM結果被噪聲嚴重干擾。L1模對強噪聲的容忍性較L2模更好，但由于L1模在零值處不可導，因而常用Huber范數或L1、L2混合模替代[22]。T分布是另一種較L2模更穩健的方法，已被成功用于噪聲數據全波形反演[23]。Aravkin等[24]認為，T分布相比Huber范數和混合模，在缺失數據的條件下更穩健，且沒有多余參數，因而簡單實用，而Huber范數和混合模的結果嚴重依賴參數選取，需要大量的嘗試。本文將T分布推廣到起伏地表LSRTM，通過處理含噪數據，證實了方法的有效性。同時考慮到LSRTM的計算量，引入動態相位編碼技術，將多炮數據組合成一個超道集，明顯節省了計算量、提高了計算效率。

1 貼體網格正演模擬

眾所周知，RTM及LSRTM算法的核心為波場延拓，波場模擬方法的優劣直接決定LSRTM的成敗，在起伏地形條件下更是如此。因此，在起伏地形條件下開展LSRTM，首先要研究起伏地形條件的正演模擬。傳統有限差分法在面對起伏地形時存在困難，而貼體網格則是一種很好的解決方案，在山前帶具有廣闊的應用前景。

貼體網格是一種適合復雜地表介質的網格離散方法，網格生成的原則是使離散后的網格邊界與地表形態吻合，以避免人為產生的階梯邊界引起的虛假散射，貼體網格可以由計算空間到物理空間的坐標變換獲得(圖1)。

圖1 貼體網格映射示意圖

利用鏈式法則和常規笛卡爾坐標系波動方程，可導出曲線坐標系二維聲波波動方程[25]

式中:vx和vz為質點振動速度；p為聲壓；ρ為介質密度；v為聲波速度。

由式(1)可知，曲坐標系下每個變量都要在同一網格點的兩個方向計算空間導數，常規的標準交錯網格已不滿足需求(圖2a)，計算vx的垂向導數與vz的橫向導數均需復雜的波場插值，降低了模擬精度。由于同位網格的一階中心差分缺少耗散特性，引起高頻振蕩現象，因此需要人為濾波，增加了實現的復雜性、降低了模擬精度。全交錯網格機制的主要思想是速度和聲壓的不同分量交錯定義在同一網格點(圖2b)。目前全交錯網格主要用于水平地表和矩形規則網格，本文將全交錯網格引入曲線坐標系[25]。由全交錯網格的網格定義機制可知，其滿足式(1)的交錯分布特性，避免了插值和濾波運算，提高了模擬精度、降低了實現復雜性。由于同一變量定義在同一網格的兩個不同位置，所以需要分別更新、計算，具體差分格式與常規標準交錯網格相同，因而簡單、易實現。

圖2 網格剖分示意圖

邊界條件是正演模擬的關鍵，受計算機內存和計算速度等因素限制，數值模擬的范圍必須是有限的，除上邊界外的其他三個邊界都會產生人工邊界問題。為此，本文采用完全匹配層吸收邊界[26]，使波能夠自由地穿過邊界，而不產生反射；上邊界為自由邊界，本文采用牽引力鏡像法實施[20]。

2 起伏地表LSRTM基本原理

2.1 起伏地表線性化波動方程

為了推導的簡便，仿照一階速度—應力方程與二階彈性波方程的相互轉換方式，可將式(1)重新整理為

(2)

(3)

由場的疊加原理可知，總波場p可理解為由背景介質產生的背景波場p0和由擾動介質產生的擾動波場ps疊加而成，即

p=p0+ps

(4)

p0與p都滿足波動方程，分別將其代入式(2)并相減，然后應用Born近似，由p0替代p0+ps，可得ps的控制方程

(5)

式(5)即為曲線坐標系線性化Born正演(反偏移)方程。可見，ps是由p0與Δs2的相互作用作為二次震源在背景介質中傳播的波場，與Δs2呈線性關系，具有明確的物理含義。

定義背景介質中的格林函數G滿足

(6)

利用格林函數可將式(5)中的ps表示為

(7)

式中：m=-Δs2，將其定義為模型參數； Ω為積分空間范圍。為方便后續的推導，式(7)可寫成算子的形式

ps=Lm

(8)

2.2 基于T分布的起伏地表LSRTM

基于反演的成像方法尋求最優的地下介質模型，以使正演波場與觀測波場殘差的模最小，是一個最小范數問題。常規基于L2模的LSRTM目標泛函為

(9)

式中pobs為觀測記錄。

L1模相比L2模更穩健，但當數據誤差接近于零時，應用L1范數準則求取目標函數梯度會出現不穩定，因而常用Huber范數或混合模代替，但由于最后的反演結果嚴重依賴于參數的選取，一般需要多次嘗試，大大增加了計算量[27]。

與混合模不同的是，T分布對參數的依賴性較弱，因而更簡單、高效。基于T分布的LSRTM目標泛函為

(10)

式中：α為自由度參數；σ為尺度參數。特別地，當α=1時，式(10)為柯西分布。大量試算表明，α=2、σ=1是一組較好的參數組合，適用于絕大多數情況。

采用梯度導引類算法(最速下降或共軛梯度)求解，需要計算目標泛函關于模型參數的梯度，L2模和T分布目標泛函的梯度公式可分別寫為

(11)

(12)

式中上角“*”表示共軛轉置。由式(11)、式(12)可以看出，兩種范數的不同之處在于：L2模直接利用波場殘差計算梯度，T分布則利用加權后的殘差記錄計算梯度，其他更新流程完全相同。因此，基于常規L2模的各種加速方法都可直接應用于T分布，如相位編碼技術。

LSRTM的計算量過于龐大，從而限制了其推廣應用。考慮到LSRTM的計算量與炮數成線性關系，因而通過相位編碼技術[28]將多個炮集組合成一個超道集，可有效減小計算量，在此基礎上應用共享存儲并行編程(OpenMP)技術，可進一步提高計算效率，本文選用動態震源極性編碼[29]。

3 模型試算

本文給出三個模型算例，以驗證本文算法的有效性和適用性。

3.1 傾斜界面模型

以傾斜界面的網格化離散為例，說明貼體網格在刻畫不規則界面時的優勢。圖3為傾斜界面模型剖分示意圖。由圖可見：常規矩形網格不能很好地離散傾斜界面，在地表界面處產生很多階梯狀的毛刺(圖3a)；貼體網格則可很好地逼近復雜界面(圖3b)，不僅可以完全擬合真實界面，且在界面處滿足很好的正交性，有利于實施邊界條件。

圖4為0.3s時刻波場快照。由圖可見：常規有限差分算法對傾斜界面的離散作用較差，產生較強的虛假繞射和散射，嚴重干擾了有效波場(圖4a)；貼體全交錯網格算法很好地擬合了傾斜界面(圖4b)。

圖3 傾斜界面模型剖分示意圖

圖4 0.3s時刻波場快照

3.2 起伏地表洼陷模型

圖5為起伏地表洼陷模型，圖6為貼體網格剖分及其山峰處的局部放大圖。可見，貼體網格對復雜界面的適應性較強。圖7為本文算法得到的RTM結果及其Laplace濾波結果。由圖可見：首先，在本文算法得到的RTM結果中地下構造基本被低頻噪聲掩蓋(圖7a)。其次，RTM結果的Laplace濾波結果雖然可正確成像地下構造，但仍然存在如下問題(圖7b)：①偏移噪聲大。Laplace濾波不能徹底去除低頻噪聲，且還引入了高頻噪聲；②反射同相軸中間能量強、兩側能量較弱，即振幅均衡性不佳；③由于地下照明強度隨深度的增大而減弱，因而RTM結果深部能量較弱，振幅保真性差。

采用本文提出的起伏地表LSRTM算法可以有效地解決常規RTM存在的問題。然而LSRTM的計算量過于龐大，目前的計算機資源無法滿足運算需要。多震源技術可有效緩解計算量問題，但會引入串擾噪聲，采用動態相位編碼技術可很好地壓制串擾噪聲，在大幅降低計算量的同時，可得到與常規LSRTM算法相當的結果。將101炮地震數據利用震源極性編碼方式組合成一個超道集，使計算量相當于單炮情形，從而大大緩解了計算需求。圖8為基于相位編碼的起伏地表LSRTM成像結果。由迭代 30次的LSRTM成像結果可見，地下構造清晰，在振幅保真性、均衡性、壓制低頻噪聲等方面較常規RTM結果明顯改善，但存在由編碼引入的較強高頻串擾噪聲(圖8a)；迭代80次的LSRTM成像結果與理論反射率模型非常接近，有效地壓制了串擾噪聲(圖8b)。

圖5 起伏地表洼陷模型

圖6 貼體網格剖分(a)及其山峰處的局部放大圖(b)

圖7 本文算法得到的RTM結果(a)及其Laplace濾波(b)

圖8 基于相位編碼的起伏地表LSRTM成像結果

圖9為殘差曲線與起伏地表洼陷模型上界面中點處振幅譜。由圖可見：①隨著迭代次數增加，數據殘差和模型殘差都逐漸減小，開始下降趨勢明顯，隨著迭代次數的進一步增大，下降趨勢減慢；由于目標泛函是數據空間的擬合，因此數據殘差能收斂到更低值(圖9a)。②隨著迭代次數增大，高頻成分逐漸得到恢復，振幅趨近于真值，表明LSRTM可以反演高頻成分(圖9b)。

測試發現，三種范數對隨機噪聲的容忍度基本相同，噪聲較弱時效果較好，噪聲較強時效果較差。當數據中含有異常值時，混合模及T分布較L2模具有較大優勢。圖10為含噪單炮記錄，圖11為含脈沖噪聲的LSRTM結果。由圖可見： L2模LSRTM結果含有較多干擾(圖11a)； 混合模(圖11b)和T分布(圖11c)LSRTM結果顯著消除了脈沖噪聲的影響，但前者需要多次嘗試以尋求最優參數，因此T分布LSRTM更簡單、實用。

3.3 起伏地表中原模型

圖12為速度模型及其反射率模型。由圖可見，起伏地表中原模型不僅地表起伏劇烈，且地下構造復雜，可檢驗偏移算法的成像效果。采用動態相位編碼技術將218炮數據組合成一個超道集，然后利用本文算法迭代計算。圖13為起伏地表中原模型LSRTM結果。由圖可見：在1次迭代的LSRTM結果中存在非常強的近地表低頻噪聲，不能識別地下構造，且由編碼引入的高頻串擾也很強(圖13a)；20次迭代的LSRTM結果(圖13b)已明顯壓制了低頻噪聲，基本可識別地下構造，但仍存在較強的高頻串擾，且橫向振幅均衡性較差；隨迭代次數進一步增大，LSRTM結果質量變好，基本消除了高頻串擾和低頻噪聲，顯著提高了振幅均衡性和保真性(圖13c、圖13d)。

圖9 殘差曲線(a)與起伏地表洼陷模型上界面中點處振幅譜(b)

圖10 含噪單炮記錄

圖11 含脈沖噪聲的LSRTM結果

圖12 速度模型(a)及其反射率模型(b)

圖13 起伏地表中原模型LSRTM結果

4 結束語

針對起伏地表地震偏移成像存在的問題，本文發展了基于貼體全交錯網格的起伏地表LSRTM算法。通過理論分析及模型試算，得到以下認識：

(1)貼體網格能很好地擬合起伏界面，避免了矩形網格中由于階梯離散產生的虛假繞射波，且可適應任意復雜地表情形。在曲線坐標系下采用全交錯網格，避免了標準交錯網格的波場插值和同位網格的高頻振蕩問題，提高了模擬精度，減小了實現復雜度。

(2)與常規RTM相比，LSRTM可實現真振幅成像，提高了分辨率、減小了偏移噪聲，顯著提高了振幅均衡性和保真性。本文測試的三種范數在隨機噪聲下的表現基本相同，但相比常規L2模，混合模與T分布能較好地壓制脈沖噪聲的影響，且T分布不依賴于參數選取，因此更簡單、實用。

(3)在起伏地表LSRTM算法中，加入動態相位編碼技術不僅能極大地降低計算成本，且可有效壓制串擾噪聲，在此基礎上應用OpenMP技術，可顯著提高計算效率，從而將LSRTM的計算成本降低到同常規RTM相同的水平。

尚需指出，本文算法僅進行了理論模型測試，還未用于實際資料處理，且采用的模型近地表變化較為平緩，在橫向變速劇烈的近地表條件下算法的適應性還需進一步探討。此外，如何將其推廣到GPU等快速計算設備上進一步提高計算效率，如何通過預條件和正則化算子等增加算法的穩定性并加快收斂速度等是下一步的研究方向。