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

李 娜

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

·偏移成像·

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

李 娜*

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

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

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

最小二乘逆時偏移 相位編碼 Huber范數(shù) 多震源 隨機(jī)最優(yōu)化

1 引言

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

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

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

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

2 方法原理

2.1 線性Born正演

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

(1)

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

(2)

(3)

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

(4)

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

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

(5)

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

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

ps=Lm

(6)

2.2 基于Huber模的LSRTM

常規(guī)LSRTM基于觀測數(shù)據(jù)與反偏移模擬數(shù)據(jù)的L2模的擬合，其目標(biāo)泛函為

(7)

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

上述基于L2范數(shù)擬合的LSRTM算法對噪聲非常敏感，尤其是當(dāng)?shù)卣饠?shù)據(jù)中含有奇異值時，LSRTM結(jié)果被噪聲嚴(yán)重干擾。L1范數(shù)相比L2范數(shù)更加穩(wěn)健，但當(dāng)數(shù)據(jù)誤差接近于0時，應(yīng)用L1范數(shù)準(zhǔn)則求取的目標(biāo)函數(shù)梯度會不穩(wěn)定，因而常用Huber模或混合模(hybrid)代替。兩種目標(biāo)泛函的形式如下

(8)

(9)

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

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

(10)

rnorm=Wnorm(Lm-pobs)

(11)

其中，不同的范數(shù)對應(yīng)不同的加權(quán)因子Wnorm，分別表示為

WL2=1

(12)

(13)

(14)

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

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

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

mk+1=mk-bkgk

(15)

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

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

2.3 隨機(jī)最優(yōu)化相位編碼技術(shù)

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

編碼函數(shù)多種多樣，考慮到簡便易實現(xiàn)的優(yōu)點，本文選用震源極性編碼，編碼函數(shù)為

(16)

式中：Nnit為總迭代次數(shù)；s為當(dāng)前炮號；k為當(dāng)前迭代次數(shù)；γ為隨機(jī)的+1或-1，在每次迭代中隨機(jī)生成。

隨機(jī)最優(yōu)化思想通過加權(quán)平均前期得到的迭代梯度減小了梯度的隨機(jī)波動，本文同樣采用指數(shù)衰減加權(quán)平均，優(yōu)化后的梯度表達(dá)式為

(17)

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

3 模型試算

3.1 簡單層狀模型

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

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

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

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

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

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

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

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

圖4 四種方法圖2左數(shù)據(jù)切除直達(dá)波后第50次迭代成像結(jié)果對比

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

圖6 采用9個超道集的含噪聲隨機(jī)最優(yōu)化Huber范數(shù)LSRTM結(jié)果

3.2 Marmousi模型

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

圖7 Marmousi速度模型

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

圖8 四種方法逆時偏移結(jié)果對比

4 結(jié)論

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

(2)相位編碼技術(shù)可顯著提升LSRTM的計算效率，但收斂較慢。隨機(jī)最優(yōu)化的動態(tài)相位編碼技術(shù)相比傳統(tǒng)相位編碼算法收斂更快，且可有效壓制串?dāng)_噪聲，進(jìn)一步節(jié)省計算成本、提高計算效率，降低LSRTM的計算成本(同常規(guī)逆時偏移的計算成本相同)。

[1] 李振春.地震偏移成像技術(shù)研究現(xiàn)狀與發(fā)展趨勢.石油地球物理勘探,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] 黃建平，李振春，孔雪等.碳酸鹽巖裂縫型儲層最小二乘偏移成像方法研究.地球物理學(xué)報,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] 劉玉金,李振春,吳丹等.局部傾角約束最小二乘偏移方法研究.地球物理學(xué)報,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的最小二乘偏移方法.吉林大學(xué)學(xué)報(地球科學(xué)版),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] 楊其強(qiáng)，張叔倫.最小二乘傅里葉有限差分法偏移.地球物理學(xué)進(jìn)展,2008,23(2):433-437. Yang Qiqiang,Zhang Shulun.Least-squares Fourier finite-difference migration.Progress in Geophysics,2008,23(2):433-437.

[11] 周華敏，陳生昌，任浩然等.基于照明補(bǔ)償?shù)膯纬滩ㄗ钚《似?地球物理學(xué)報，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] 黃建平，曹曉莉，李振春等.最小二乘逆時偏移在近地表高精度成像中的應(yīng)用.石油地球物理勘探,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] 李振春，郭振波，田坤.黏聲介質(zhì)最小平方逆時偏移.地球物理學(xué)報，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] 黃建平,李闖,李慶洋等.一種基于平面波靜態(tài)編碼的最小二乘逆時偏移方法.地球物理學(xué)報,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] 李慶洋,黃建平,李振春等.偽深度域聲波數(shù)值模擬方法及應(yīng)用.石油地球物理勘探,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] 李慶洋,黃建平,李振春等.優(yōu)化的多震源最小二乘逆時偏移.石油地球物理勘探,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年畢業(yè)于山東理工大學(xué)數(shù)學(xué)統(tǒng)計學(xué)專業(yè)，獲理學(xué)學(xué)士學(xué)位；2010年畢業(yè)于蘇州大學(xué)概率論與數(shù)理統(tǒng)計專業(yè)，獲理學(xué)碩士學(xué)位；2014年獲中國石油大學(xué)(華東)地質(zhì)資源與地質(zhì)工程專業(yè)工學(xué)博士學(xué)位；現(xiàn)在中原油田物探研究院從事地震波傳播正演模擬、偏移方法以及裂縫儲層預(yù)測和描述方法研究。

1000-7210(2017)05-0941-07

P631

A

10.13810/j.cnki.issn.1000-7210.2017.05.006

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

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

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