基于Huber范數的多震源最小二乘逆時偏移

李 娜

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

·偏移成像·

基于Huber范數的多震源最小二乘逆時偏移

李 娜*

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

李娜.基于Huber范數的多震源最小二乘逆時偏移.石油地球物理勘探，2017,52(5):941-947.

將Huber范數推廣到最小二乘逆時偏移中，提升了算法的穩健性;同時考慮到最小二乘逆時偏移的計算量較大，將隨機最優化相位編碼技術引入到基于Huber范數的多震源最小二乘逆時偏移中，可有效提高計算效率,壓制串擾噪聲;最后通過常規有限差分正演模擬記錄試算，驗證了算法的有效性和對復雜模型的適用性。

最小二乘逆時偏移 相位編碼 Huber范數 多震源 隨機最優化

1 引言

近年來，隨著油氣勘探精度的逐漸提高，地震波成像方法也在不斷發展[1]。相比于其他偏移方法，逆時偏移應用雙程波波動方程進行波場延拓，避免了對波動方程的近似，無傾角限制，是目前較為精確的常規偏移成像方法。但逆時偏移仍屬于常規偏移的范疇，其偏移算子是線性Born正演算子的共軛轉置，而不是它的逆[2]。當地震數據采集不足或不規則、地下構造情況復雜以及波場帶寬有限時，常規偏移方法只能對地下構造模糊成像，無法滿足巖性油氣藏勘探開發的需求。

最小二乘偏移將成像看作是最小二乘意義下的反演問題，通過迭代不斷擬合觀測數據和反偏移數據，可以得到振幅保真性更好、分辨率更高、偏移噪聲更少的成像結果。與常規偏移類似，最小二乘偏移主要可以分為三類：射線類最小二乘偏移、單程波類最小二乘偏移和雙程波類最小二乘偏移。射線類最小二乘偏移主要指最小二乘Kirchhoff偏移[3-5]和最小二乘高斯束偏移[6]。單程波類最小二乘偏移方法又可具體分為最小二乘分步傅里葉偏移[7,8]、最小二乘傅里葉有限差分偏移[9,10]和最小二乘廣義屏偏移[11]等。雙程波類最小二乘偏移即最小二乘逆時偏移(Least-Squares Reverse Time Migration,LSRTM)。Dai等[12]首先提出了線性和擬線性的LSRTM，并通過混合震源相位編碼技術測試了方法的正確性和有效性；黃建平等[13]將LSRTM應用到近地表保幅成像；郭振波等[14]提出了LSRTM的真振幅成像；李振春等[15]將LSRTM推廣到更加符合實際地下情況的黏聲介質中。

上述LSRTM的目標泛函都是基于L2范數擬合，考慮到地震數據常含有較多噪聲，常用的基于L2范數擬合的LSRTM算法對噪聲非常敏感，尤其是當地震數據中含有奇異值時，常規LSRTM結果被噪聲嚴重干擾。L1模對大噪聲的容忍性比L2模更好，但由于L1模在0值處不可導，常用Huber模或者混合模替代[16]，目前Huber模已在全波形反演中取得較好效果。為此，本文將Huber模發展到LSRTM中，顯著提升了LSRTM算法的穩定性，與常規LSRTM相比，本文算法在地震數據含有脈沖噪聲時仍能得到較好的反演結果。

LSRTM的計算量過大，限制了其進一步推廣應用。多震源相位編碼技術可顯著提高LSRTM的計算效率，Romero等[17]首次引入相位編碼的概念；Krebs等[18]將相位編碼技術應用于全波形反演中；Dai等[12]實現了基于相位編碼的多震源LSRTM，有效降低了計算成本; 黃建平等[19]提出了基于平面波靜態編碼的LSRTM; 李慶洋等[20]提出偽深度域聲波數值模擬方法；Moghaddam等[21]首先提出了隨機最優化相位編碼的全波形反演算法；隨后，李慶洋等[22]將隨機最優化思想推廣到相位編碼的LSRTM中，相比傳統相位編碼，LSRTM方法收斂更快。為此，本文將隨機最優化相位編碼技術應用于基于Huber泛函的多震源LSRTM中，與常規基于Huber范數的相位編碼多震源LSRTM相比收斂更快，在適用于含脈沖噪聲數據的同時顯著提高了計算效率。

2 方法原理

2.1 線性Born正演

二維常密度聲波方程可表示為

(1)

式中:s為慢度場；f為震源項；p是聲壓場。

(2)

(3)

將式(2)與式(3)相減，并應用Born近似，用背景波場p0替代總波場p0+ps，可得到擾動波場ps的控制方程

(4)

式(4)即為線性化Born正演(反偏移)方程，從中可以看出，擾動波場是背景波場與慢度擾動的相互作用，作為二次震源在背景介質中傳播產生的場，與慢度擾動呈線性關系，具有明確的物理含義。即實際編程實現時，線性化正演(反偏移)過程需要求解兩次正演方程：通過式(2)得到當前時刻的背景波場，然后再利用式(4)求得當前時刻的擾動波場。

利用背景介質中的格林函數G，可將式(4)中的擾動波場表示為

(5)

式中:Ω為積分域；m=-Δs2， 將其定義為模型參數。

為方便后續的推導，式(5)可寫成算子的形式，即

ps=Lm

(6)

2.2 基于Huber模的LSRTM

常規LSRTM基于觀測數據與反偏移模擬數據的L2模的擬合，其目標泛函為

(7)

式中：pobs為觀測記錄；ps=Lm為反偏移模擬記錄； ‖·‖2代表向量的L2范數。

上述基于L2范數擬合的LSRTM算法對噪聲非常敏感，尤其是當地震數據中含有奇異值時，LSRTM結果被噪聲嚴重干擾。L1范數相比L2范數更加穩健，但當數據誤差接近于0時，應用L1范數準則求取的目標函數梯度會不穩定，因而常用Huber模或混合模(hybrid)代替。兩種目標泛函的形式如下

(8)

(9)

式中： Δd為模擬數據與觀測數據的殘差；ε=λ×mean(|pobs|)為控制L1范數與L2范數權重的閾值,λ為0～1之間的數，一般取0.5左右。

采用梯度導引類算法(最速下降或共軛梯度)求解，需要計算目標泛函關于模型參數的梯度，三種不同范數的梯度公式都可統一寫成

(10)

rnorm=Wnorm(Lm-pobs)

(11)

其中，不同的范數對應不同的加權因子Wnorm，分別表示為

WL2=1

(12)

(13)

(14)

式中a為穩定因子，為避免分母為0，一般取0.01×mean(|pobs|)左右。

從式(10)可以看出，幾種不同范數的梯度公式具有相同的形式，僅殘差波場的加權因子不同。梯度求取過程也就是一個逆時偏移(RTM)的過程，與常規RTM不同的是，求梯度時的反傳波場為殘差記錄。考慮到波場存儲是RTM的瓶頸，本文計算梯度時采用震源波場重建方法，即每個時刻只在邊界存儲正傳波場信息，反傳時將其作為邊界條件重構震源波場，大幅度降低了存儲量[20]。

通過式(10)計算得到梯度后，模型的更新過程為

mk+1=mk-bkgk

(15)

式中：gk為第k次迭代的梯度；bk為第k次迭代的更新步長，步長可通過線性搜索等方法計算得到。

從梯度公式還可以看出，三種范數的唯一不同之處在于，L2模直接利用波場殘差計算梯度，而Huber模和混合模則是利用加權后的殘差記錄計算梯度，因而三種范數的更新流程相同，常規L2模下的各種加速方法都可直接應用于Huber模和混合模中。

2.3 隨機最優化相位編碼技術

相位編碼技術是目前非常實用的一項提高地震采集和數據處理效率的技術，其基本思想是通過設計編碼函數將原本多次采集的炮集數據壓縮成一次采集的混疊數據，從而減少數據規模，提高采集和處理效率。Moghaddam等[21]在多震源全波形反演中研究發現，通過加權平均前幾次迭代的相位編碼梯度可加快收斂速度。本文借鑒李慶洋等[22]提出的優化的多震源LSRTM方法，將隨機最優化思想推廣到基于Huber模或混合模的多震源LSRTM中。

編碼函數多種多樣，考慮到簡便易實現的優點，本文選用震源極性編碼，編碼函數為

(16)

式中：Nnit為總迭代次數；s為當前炮號；k為當前迭代次數；γ為隨機的+1或-1，在每次迭代中隨機生成。

隨機最優化思想通過加權平均前期得到的迭代梯度減小了梯度的隨機波動，本文同樣采用指數衰減加權平均，優化后的梯度表達式為

(17)

式中：j為加權的前期迭代次數，綜合考慮效果和效率，一般取為10；c為衰減因子，取為0.5。相比于常規LSRTM，基于相位編碼的LSRTM只需要修改線性正演算子和觀測數據即可，具體的計算流程與常規LSRTM相同，在此不再贅述。

3 模型試算

3.1 簡單層狀模型

首先以一個簡單的層狀模型為例驗證本文算法的正確性及有效性。模型速度場如圖1左所示，模型中既有大尺度的背斜、向斜構造，也有90°高陡構造，非常適合檢驗成像方法的優劣，對應的擾動模型如圖1右所示。網格數為401×271，縱、橫向網格間距均為10m。合成觀測數據所用計算參數為：在地表以50m間隔均勻分布81炮，每炮都是401個檢波點全接收，震源為主頻30Hz的雷克子波，時間采樣間隔為1ms，最大記錄時間為2.5s。

圖1 速度模型(左)及反射率模型(右)

LSRTM的正算子為線性Born正演算子，而一般生成炮記錄的正算子為常規有限差分正演算子，因而基于波形匹配策略的LSRTM需要首先測試兩種算子是否匹配。圖2左為常規有限差分正演模擬第41炮單炮記錄，圖2右為不含噪聲的線性Born正演模擬單炮記錄。可以看出，線性Born正演模擬記錄與常規有限差分正演模擬記錄相比，除了缺少明顯的直達波外，基本相同，振幅相位都基本一致。為了更加清晰地對比二者的相似度，從圖2左和圖2右中分別取出5道單道波形進行對比，如圖3所示。可以看出本文線性Born正演模擬算子與常規有限差分正演模擬算子匹配得較好，為其線性近似。

圖2 有限差分單炮記錄(左)及線性Born正演單炮記錄(右)

將81炮利用震源極性編碼方式組合成一個超道集，使計算量相當于單炮情形，緩解了計算需求。圖4為圖2左數據不同方法第50次迭代成像結果對比。可以看出，圖4b的結果被噪聲嚴重污染，而圖4c結果消除了脈沖噪聲的影響。對比圖4c和圖4d可以看出，本文方法可以有效改善成像效果，不僅能更好地壓制串擾噪聲，且收斂速度更快，成像結果分辨率更高。

圖5左為四種算法歸一化的數據殘差曲線。可以看出，在含脈沖噪聲的情況下常規L2模LSRTM不收斂，而基于Huber泛函的LSRTM仍然能較好地收斂。隨機最優化編碼方案相比常規相位編碼的收斂速度更快，更加接近不含噪聲情況下的收斂速度。由圖5右常規相位編碼和隨機最優化編碼策略的模型殘差曲線同樣可以看出隨機最優化策略的優勢。

需要注意的是，圖4d雖然沒有脈沖噪聲的干擾，但仍然含有較多由相位編碼引入的高頻串擾噪聲。為了壓制串擾噪聲，提高超道集的數目，將81炮利用震源極性編碼方式組合成9個超道集，得到的第50次迭代成像結果如圖6所示。對比圖4d和圖6可以看出，超道集數目的增加可以明顯壓制高頻串擾噪聲，得到同相軸更加清晰的成像結果。由于炮檢距較大，成像結果中仍然含有一些低頻噪聲，這是由于線性Born正演模擬與常規有限差分模擬結果大炮檢距不如小炮檢距吻合得好，因而不能完全壓制低頻噪聲，可以借助高通濾波、拉普拉斯濾波等進一步消除。

圖3 圖2抽取5道單道波形對比

圖4 四種方法圖2左數據切除直達波后第50次迭代成像結果對比

圖5 數據殘差收斂曲線(左)及模型殘差收斂曲線(右)

圖6 采用9個超道集的含噪聲隨機最優化Huber范數LSRTM結果

3.2 Marmousi模型

下面以Marmousi模型為例驗證算法對復雜模型的適用性。Marmousi模型速度場如圖7所示。模型參數如下：橫向網格間距為12.5m，網格點數為481，縱向網格間距為8m，網格點數為375，速度范圍為1.5～5.5km/s。采用時間二階空間八階有限差分正演模擬得到觀測記錄，觀測系統模擬標準數據的觀測系統，共有120炮，每一炮有101個接收道，炮間距為50m，道間距為25m，采用單邊放炮，最小炮檢距為50m，最大炮檢距為2550m，震源為主頻是20Hz的雷克子波，時間采樣間隔為0.6ms，最大記錄時間為3s。

圖7 Marmousi速度模型

圖8為四種方法逆時偏移成像結果對比。可以看出，不含噪聲數據的常規RTM成像剖面中各個同相軸基本清晰可見，但振幅能量不均衡、分辨率較低。對比圖8a和圖8b可以看出，LSRTM的成像結果有明顯的改善，不僅淺層的斷層刻畫更加清晰，而且深部的背斜也得到較好的凸顯，整體能量更加均衡、分辨率明顯提高(圖8b)。圖8c被噪聲嚴重污染，而Huber泛函相比L2模具有更高的噪聲容忍度，因而圖8d結果顯著消除了脈沖噪聲的影響，驗證了本文算法對復雜模型的適用性。

圖8 四種方法逆時偏移結果對比

4 結論

(1)Huber模或者混合模相比L2范數更穩健，能較好地消除地震數據中的脈沖噪聲，且計算流程與常規LSRTM相同，簡便易實現。

(2)相位編碼技術可顯著提升LSRTM的計算效率，但收斂較慢。隨機最優化的動態相位編碼技術相比傳統相位編碼算法收斂更快，且可有效壓制串擾噪聲，進一步節省計算成本、提高計算效率，降低LSRTM的計算成本(同常規逆時偏移的計算成本相同)。

[1] 李振春.地震偏移成像技術研究現狀與發展趨勢.石油地球物理勘探,2014,49(1):1-21. Li Zhenchun.Research status and development trends for seismic migration technology.OGP,2014,49(1):1-21.

[2] Claerbout J F.Earth Soundings Analysis: Processing versus Inversion.Blackwell Scientific Publications,1992.

[3] Nemeth T,Wu C J,Schuster G T.Least-squares migration of incomplete reflection data.Geophysics,1999,64(1),208-221.

[4] 黃建平，李振春，孔雪等.碳酸鹽巖裂縫型儲層最小二乘偏移成像方法研究.地球物理學報,2013,56(5):1716-1725. Huang Jianping,Li Zhenchun,Kong Xue et al.A study of least-squares migration imaging method for fractured-type carbonate reservoir.Chinese Journal of Geo-physics,2013,56(5):1716-1725.

[5] 劉玉金,李振春,吳丹等.局部傾角約束最小二乘偏移方法研究.地球物理學報,2013,56(3):1003-1011. Liu Yujin,Li Zhenchun,Wu Dan et al.The research on local slope constrained least-squares migration.Chinese Journal of Geophysics,2013,56(3):1003-1011.

[6] Hu H,Liu Y,Zheng Y et al.Least-squares Gaussian beam migration.Geophysics,2016,81(3):S87-S100.

[7] Kuehl H,Sacchi M D.Least-squares wave-equationmigration for AVP/AVA inversion.Geophysics,2003,68(1):262-273.

[8] 黃建平,孫隕松,李振春等.一種基于分頻編碼的最小二乘裂步偏移方法.石油地球物理勘探,2014,49(4):702-707. Huang Jianping,Sun Yunsong,Li Zhenchun et al.Least-squares split-step migration based on frequencydivision encoding.OGP,2014,49(4):702-707.

[9] 黃建平，薛志廣，步長城等.基于裂步DSR的最小二乘偏移方法.吉林大學學報(地球科學版),2014,44(1):369-374. Huang Jianping,Xue Zhiguang,Bu Changcheng et al.The study of least-squares migration method based on split-step DSR.Journal of Jilin University(Earth Science Edition),2014,44(1):369-374.

[10] 楊其強，張叔倫.最小二乘傅里葉有限差分法偏移.地球物理學進展,2008,23(2):433-437. Yang Qiqiang,Zhang Shulun.Least-squares Fourier finite-difference migration.Progress in Geophysics,2008,23(2):433-437.

[11] 周華敏，陳生昌，任浩然等.基于照明補償的單程波最小二乘偏移.地球物理學報，2014,57(8)：2644-2655. Zhou Huamin,Chen Shengchang,Ren Haoran et al.One-way wave equation least-squares migration based on illumination compensation.Chinese Journal of Geophysics,2014,57(8):2644-2655.

[12] Dai W,Fowler P,Schuster G T.Multi-source least-squares reverse time migration.Geophysical Prospecting,2012,60(4):681-695.

[13] 黃建平，曹曉莉，李振春等.最小二乘逆時偏移在近地表高精度成像中的應用.石油地球物理勘探,2014,49(1):107-112. Huang Jianping,Cao Xiaoli,Li Zhenchun et al.Least square reverse time migration in high resolution imaging of near surface.OGP,2014,49(1):107-112.

[14] 郭振波，李振春.最小平方逆時偏移真振幅成像.石油地球物理勘探,2014,49(1):113-120. Guo Zhenbo,Li Zhenchun.True-amplitude imaging based on least-squares reverse time migration.OGP,2014,49(1):113-120.

[15] 李振春，郭振波，田坤.黏聲介質最小平方逆時偏移.地球物理學報，2014，57(1):214-228. Li Zhenchun,Guo Zhenbo,Tian Kun.Least-squares reverse time migration in visco-acoustic medium.Chinese Journal of Geophysics,2014，57(1):214-228.

[16] Brossier R,Operto S and Virieux J.Which data resi-dual norm for robust elastic frequency-domain full waveform inversion.Geophysics,2010,75(3):R37-R46.

[17] Romero L A,Ghiglia D C,Ober C C et al.Phase encoding of shot records in prestack migration.Geophy-sics,2000,65(2):426-436.

[18] Krebs J R,Anderson J E,Hinkley.Fast full-wavefield seismic inversion using encoded sources.Geophysics,2009,74(6):WCC177-WCC188.

[19] 黃建平,李闖,李慶洋等.一種基于平面波靜態編碼的最小二乘逆時偏移方法.地球物理學報,2015,58(6):2046-2056. Huang Jianping,Li Chuang,Li Qingyang et al.Least-square reverse time migration with static plane-wave encoding.Chinese Journal of Geophysics,2015,58(6):2046-2056.

[20] 李慶洋,黃建平,李振春等.偽深度域聲波數值模擬方法及應用.石油地球物理勘探,2015,50(2):283-289. Li Qingyang,Huang Jianping,Li Zhenchun et al.Acoustic wave numerical simulation in pseudo-depth domain.OGP,2015,50(2):283-289.

[21] Moghaddam P P,Keers H,Herrmann F J et al.A new optimization approach for source-encoding full-waveform inversion.Geophysics,2013,78(3):R125-R132.

[22] 李慶洋,黃建平,李振春等.優化的多震源最小二乘逆時偏移.石油地球物理勘探,2016,51(2):334-341. Li Qingyang,Huang Jianping,Li Zhenchun et al.Optimized multi-source least-squares reverse time migration.OGP,2016,51(2):334-341.

(本文編輯：金文昱)

李娜 在站博士后，1985年生；2007年畢業于山東理工大學數學統計學專業，獲理學學士學位；2010年畢業于蘇州大學概率論與數理統計專業，獲理學碩士學位；2014年獲中國石油大學(華東)地質資源與地質工程專業工學博士學位；現在中原油田物探研究院從事地震波傳播正演模擬、偏移方法以及裂縫儲層預測和描述方法研究。

1000-7210(2017)05-0941-07

P631

A

10.13810/j.cnki.issn.1000-7210.2017.05.006

*河南省濮陽市華龍區慶山街3號中原油田物探研究院，457001。Email：lina19202@163.com

本文于2016年8月29日收到，最終修改稿于2017年6月23日收到。

本項研究受國家油氣重大專項項目(2016ZX05006-004)資助。