李 慶 洋

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

0 引言

隨著油氣勘探、開發(fā)的深入，勘探難度逐漸加大，在復(fù)雜構(gòu)造、復(fù)雜巖性、復(fù)雜地表條件等地區(qū)進(jìn)行油氣地震勘探面臨諸多難題，其中劇烈起伏的地表已成為制約上述地區(qū)地震勘探發(fā)展的瓶頸之一[1]。傳統(tǒng)的偏移成像方法大多基于水平地表假設(shè)，對起伏地形的適應(yīng)性不強(qiáng)，容易導(dǎo)致復(fù)雜地表條件下的地下成像質(zhì)量不高[2]。

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

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

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

1 貼體網(wǎng)格正演模擬

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

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

圖1 貼體網(wǎng)格映射示意圖

利用鏈?zhǔn)椒▌t和常規(guī)笛卡爾坐標(biāo)系波動方程，可導(dǎo)出曲線坐標(biāo)系二維聲波波動方程[25]

式中:vx和vz為質(zhì)點(diǎn)振動速度；p為聲壓；ρ為介質(zhì)密度；v為聲波速度。

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

圖2 網(wǎng)格剖分示意圖

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

2 起伏地表LSRTM基本原理

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

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

(2)

(3)

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

p=p0+ps

(4)

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

(5)

式(5)即為曲線坐標(biāo)系線性化Born正演(反偏移)方程?？梢姡琾s是由p0與Δs2的相互作用作為二次震源在背景介質(zhì)中傳播的波場，與Δs2呈線性關(guān)系，具有明確的物理含義。

定義背景介質(zhì)中的格林函數(shù)G滿足

(6)

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

(7)

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

ps=Lm

(8)

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

基于反演的成像方法尋求最優(yōu)的地下介質(zhì)模型，以使正演波場與觀測波場殘差的模最小，是一個最小范數(shù)問題。常規(guī)基于L2模的LSRTM目標(biāo)泛函為

(9)

式中pobs為觀測記錄。

L1模相比L2模更穩(wěn)健，但當(dāng)數(shù)據(jù)誤差接近于零時，應(yīng)用L1范數(shù)準(zhǔn)則求取目標(biāo)函數(shù)梯度會出現(xiàn)不穩(wěn)定，因而常用Huber范數(shù)或混合模代替，但由于最后的反演結(jié)果嚴(yán)重依賴于參數(shù)的選取，一般需要多次嘗試，大大增加了計算量[27]。

與混合模不同的是，T分布對參數(shù)的依賴性較弱，因而更簡單、高效?；赥分布的LSRTM目標(biāo)泛函為

(10)

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

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

(11)

(12)

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

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

3 模型試算

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

3.1 傾斜界面模型

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

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

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

圖4 0.3s時刻波場快照

3.2 起伏地表洼陷模型

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

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

圖5 起伏地表洼陷模型

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

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

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

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

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

3.3 起伏地表中原模型

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

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

圖10 含噪單炮記錄

圖11 含脈沖噪聲的LSRTM結(jié)果

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

圖13 起伏地表中原模型LSRTM結(jié)果

4 結(jié)束語

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

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

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

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

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