時(shí)頻域振幅相位聯(lián)合的最小二乘逆時(shí)偏移

胡勇， 潘冬明*， 吳如山， 韓立國(guó)， 段超然

1 中國(guó)礦業(yè)大學(xué)， 資源與地球科學(xué)學(xué)院， 徐州 221008 2 Modeling and Imaging Laboratory, University of California, Santa Cruz 95060, USA 3 吉林大學(xué)， 地球探測(cè)科學(xué)與技術(shù)學(xué)院， 長(zhǎng)春 130026 4 常州工學(xué)院， 土木建筑工程學(xué)院， 常州 213000

0 引言

隨著油氣勘探技術(shù)的不斷發(fā)展，同時(shí)對(duì)地下構(gòu)造探測(cè)精度要求也逐漸提高，如今正從構(gòu)造勘探階段逐步轉(zhuǎn)向巖性勘探階段，其中地震數(shù)據(jù)偏移成像方法在油氣勘探領(lǐng)域占有重要的地位.近年來(lái)Kirchhoff 偏移、高斯束偏移、單程波偏移和逆時(shí)偏移等地震成像方法得到了快速發(fā)展(秦寧等，2015；You et al.，2019；豆輝和徐逸鶴，2019；Sun et al., 2020).在復(fù)雜介質(zhì)中，逆時(shí)偏移(RTM)方法具有成像精度高的優(yōu)點(diǎn)，能夠適用于橫向速度變化劇烈的區(qū)域，獲取地下構(gòu)造高精度反射系數(shù)(Sun and Zhang, 2009；馬方正等，2016；You et al.，2017).但是逆時(shí)偏移的伴隨算子是正演算子的共軛轉(zhuǎn)置，而并非其逆算子(Claerbout，1991)，因此無(wú)法獲得真振幅成像結(jié)果.

為了進(jìn)一步提高逆時(shí)偏移的成像精度，發(fā)展了最小二乘逆時(shí)偏移方法(LSRTM)，其在反演的理論框架下，利用反偏移數(shù)據(jù)來(lái)不斷地與觀測(cè)數(shù)據(jù)相匹配，并通過(guò)優(yōu)化算法進(jìn)行迭代，最終獲取地下構(gòu)造的高精度成像結(jié)果.關(guān)于LSRTM的發(fā)展，首先由Nemeth等(1999)在最小二乘目標(biāo)函數(shù)的框架下實(shí)現(xiàn)了Kirchhoff偏移.隨后，Dai等(2012)、Dai和Schuster(2013)利用多源地震數(shù)據(jù)進(jìn)行測(cè)試，很大程度上提高了LSRTM的計(jì)算效率.任浩然等(2013)將Hessian算子應(yīng)用于LSRTM中，獲得了相對(duì)保真的成像結(jié)果，改善了地震偏移的成像精度.Tan和Huang(2014)在原始LSRTM的基礎(chǔ)上，提出了波場(chǎng)分離成像條件，并通過(guò)更新震源波場(chǎng)實(shí)現(xiàn)陡傾角斷層成像.此外，為了克服LSRTM對(duì)偏移模型的依賴問(wèn)題，劉玉金和李振春(2015)提出擴(kuò)展成像條件下的LSRTM方法，測(cè)試結(jié)果表明該方法可以獲得更準(zhǔn)確的振幅屬性信息.Huang等(2015)、周東紅等(2020)提出了帶地形的LSRTM方法，進(jìn)一步考慮地表起伏對(duì)波動(dòng)方程成像結(jié)果的影響，指出無(wú)需對(duì)地震數(shù)據(jù)進(jìn)行校正便可以獲得起伏地表高精度的成像結(jié)果.地震數(shù)據(jù)中的多次波同樣攜帶地下構(gòu)造信息，劉學(xué)建和劉伊克(2016)建立了基于多次波的LSRTM方法，并指出該方法能提供比反射波偏移成像更大的照明角度.陳生昌和周華敏(2018)在反射波成像的基礎(chǔ)上，研究了入射波與反射波在傳播算子上的差異，提出了新的反射波LSRTM方法.方修政等(2018)指出基于常規(guī)互相關(guān)成像條件的LSRTM梯度含有很強(qiáng)的低頻噪聲，為此提出了基于逆散射成像條件的LSRTM方法.田坤等(2018)、鞏向博等(2019)利用多源地震數(shù)據(jù)，在稀疏約束的理論框架下展開(kāi)研究.Liu等(2020)在LSRTM中加入了Gabor反褶積濾波來(lái)提高鹽下構(gòu)造的成像精度.Li等(2020)通過(guò)引入權(quán)重因子來(lái)衰減地震數(shù)據(jù)中較強(qiáng)的反射波信息，進(jìn)而增強(qiáng)弱散射信號(hào)，實(shí)現(xiàn)改善深部構(gòu)造成像精度的目標(biāo).經(jīng)過(guò)多年的研究，LSRTM技術(shù)得到了快速地發(fā)展，逐漸成為當(dāng)前的研究熱點(diǎn)問(wèn)題.但是地震波傳播是一個(gè)較為復(fù)雜的過(guò)程，使得LSRTM成像結(jié)果仍然受地震數(shù)據(jù)振幅信息影響較為嚴(yán)重.為了減弱振幅信息的影響，在互相關(guān)LSRTM目標(biāo)函數(shù)的基礎(chǔ)上提出了很多的改進(jìn)的方法(Zhang et al.，2015；Liu et al.,2016, 2017；李慶洋等，2016；Yi et al.，2019).但是如何在LSRTM方法基礎(chǔ)上，實(shí)現(xiàn)透過(guò)上覆強(qiáng)散射地質(zhì)體來(lái)獲得深部精細(xì)構(gòu)造的高精度成像目標(biāo)，仍有待進(jìn)一步研究.

為此，本文在前人的研究基礎(chǔ)上，提出時(shí)頻域振幅相位聯(lián)合的最小二乘逆時(shí)偏移方法(PA-LSRTM)，來(lái)實(shí)現(xiàn)對(duì)地下深部精細(xì)構(gòu)造的高精度成像.本文首先介紹PA-LSRTM目標(biāo)函數(shù)的構(gòu)建方法；其次，為了減弱地震數(shù)據(jù)振幅信息對(duì)成像結(jié)果的影響，在目標(biāo)函數(shù)中引入振幅權(quán)重因子，提高弱地震信號(hào)的可成像精度；然后，詳細(xì)推導(dǎo)了PA-LSRTM目標(biāo)函數(shù)對(duì)模型參數(shù)的梯度算子，并給出擾動(dòng)模型迭代表達(dá)式.最后，利用Marmousi模型和鹽丘模型進(jìn)行數(shù)值測(cè)試，并對(duì)RTM、LSRTM和PA-LSRTM的成像結(jié)果進(jìn)行對(duì)比分析.

1 方法原理

1.1 線性化Born正演

在常密度介質(zhì)中，聲波方程可以表示為

(1)

其中s為慢度場(chǎng)，即速度的倒數(shù)；u為聲壓波場(chǎng)；f為震源子波.式(1)中可以將慢度場(chǎng)的平方分解為背景慢度場(chǎng)平方與擾動(dòng)慢度場(chǎng)平方的和(李慶洋等，2016)：

(2)

u=u0+us.

(3)

根據(jù)式(3)，背景波場(chǎng)和擾動(dòng)波場(chǎng)滿足以下關(guān)系：

(4)

+f，

(5)

當(dāng)擾動(dòng)場(chǎng)us?u0時(shí)，可以用背景場(chǎng)代替總場(chǎng)u0≈u0+us.因此，可以獲得擾動(dòng)波場(chǎng)的表達(dá)式，即：

(6)

欲求解擾動(dòng)波場(chǎng)us，只需要在偏移模型上將背景波場(chǎng)u0代入公式(6)中，再做一次波動(dòng)方程正演模擬即可.

1.2 最小二乘逆時(shí)偏移原理

LSRTM是利用模擬數(shù)據(jù)來(lái)不斷地與觀測(cè)數(shù)據(jù)相匹配，通過(guò)最小化目標(biāo)函數(shù)來(lái)獲取擾動(dòng)模型的高精度成像結(jié)果.因此，LSRTM的目標(biāo)函數(shù)可以定義為

(7)

(8)

目標(biāo)函數(shù)對(duì)模型參數(shù)的梯度可以表示為

(9)

其中T表示矩陣轉(zhuǎn)置；ns、nr分別表示震源數(shù)目和檢波器數(shù)目；Rs=(us-ds)表示伴隨震源.因此，目標(biāo)函數(shù)對(duì)模型參數(shù)的梯度可以表示為反傳波場(chǎng)與正傳波場(chǎng)的零延遲互相關(guān)：

(10)

(11)

本文設(shè)定初始模型參數(shù)為m0=0，對(duì)應(yīng)的模擬數(shù)據(jù)us=0，此時(shí)Rs=-ds.因此，LSRTM第一次的梯度等價(jià)于常規(guī)逆時(shí)偏移成像結(jié)果.隨著模型參數(shù)不斷地被更新，模擬數(shù)據(jù)us逐漸趨近于觀測(cè)數(shù)據(jù)ds，同時(shí)對(duì)應(yīng)的LSRTM的目標(biāo)函數(shù)也逐漸減小，最終獲得地下構(gòu)造高精度偏移成像結(jié)果.

1.3 時(shí)頻域振幅相位聯(lián)合的最小二乘逆時(shí)偏移

時(shí)頻域振幅相位聯(lián)合的最小二乘逆時(shí)偏移(PA-LSRTM)是利用模擬數(shù)據(jù)的時(shí)頻域振幅相位信息來(lái)不斷地與觀測(cè)數(shù)據(jù)的時(shí)頻域振幅相位信息相匹配，通過(guò)最小化目標(biāo)函數(shù)來(lái)獲取擾動(dòng)模型的高精度成像結(jié)果.與常規(guī)LSRTM相似，PA-LSRTM的目標(biāo)函數(shù)定義為：

(12)

(13)

(14)

式(14)可以簡(jiǎn)化為(詳見(jiàn)附錄C)：

(15)

根據(jù)式(15)，伴隨震源可以定義為

(16)

式(16)為了避免分母為0，可在分母處加一個(gè)微小的常數(shù).其中，PA-LSRTM的梯度計(jì)算過(guò)程與常規(guī)LSRTM梯度計(jì)算相似，只需要將伴隨震源Rs反傳至模型空間并與正傳波場(chǎng)做零延遲互相關(guān).

圖1給出了伴隨震源某一道數(shù)據(jù)對(duì)應(yīng)的時(shí)頻域分析結(jié)果.從圖1a中可以看出，地震信號(hào)的振幅譜在淺部能量較強(qiáng)，但隨著傳播時(shí)間增大能量快速衰減.深部反射的弱振幅地震信號(hào)在目標(biāo)函數(shù)中占有很小的權(quán)重，因此增加了深部精細(xì)構(gòu)造的成像難度.當(dāng)ε=0.6時(shí)，伴隨震源的深部信號(hào)得到了明顯地增強(qiáng)(圖1b)，同時(shí)也均衡了淺部強(qiáng)能量，有助于提高深部弱反射信號(hào)在目標(biāo)函數(shù)中的權(quán)重.當(dāng)進(jìn)一步減小振幅權(quán)重因子ε=0.3，相應(yīng)的相位譜信息在PA-LSRTM目標(biāo)函數(shù)中得到了進(jìn)一步地增強(qiáng)(圖1c).當(dāng)ε=0時(shí)，從淺部到深部的振幅譜基本具有相同的權(quán)重，即使在5 s時(shí)，深部弱反射地震信號(hào)依然清晰可辨.但是，此時(shí)PA-LSRTM目標(biāo)函數(shù)中完全忽略了振幅信息，則此時(shí)目標(biāo)函數(shù)易受到噪聲相位的影響.因此，本文PA-LSRTM方法首先選用振幅權(quán)重因子ε=0.6來(lái)更新模型參數(shù)；然后，逐漸減小振幅權(quán)重因子，提高深部區(qū)域擾動(dòng)模型的成像精度，直至振幅權(quán)重因子為ε=0.3.這樣既能減弱振幅信息對(duì)深部區(qū)域成像不足的影響，同時(shí)也避免純相位目標(biāo)函數(shù)產(chǎn)生的噪聲干擾，實(shí)現(xiàn)最大化提高深部弱反射信號(hào)成像精度的目標(biāo).

1.4 時(shí)頻域振幅相位聯(lián)合的最小二乘逆時(shí)偏移流程

地震波在地下傳播過(guò)程較為復(fù)雜，常受到強(qiáng)散射介質(zhì)的屏蔽作用，致使深部構(gòu)造對(duì)應(yīng)的地震波場(chǎng)響應(yīng)較弱，難以獲得高精度成像結(jié)果.為了解決深部精細(xì)構(gòu)造成像問(wèn)題，PA-LSRTM方法通過(guò)增強(qiáng)時(shí)頻域相位譜的權(quán)重，同時(shí)減弱振幅信息的影響，提高深部弱地震信號(hào)的可成像精度，進(jìn)而實(shí)現(xiàn)深部區(qū)域和強(qiáng)散射地質(zhì)體下部構(gòu)造的高精度成像.PA-LSRTM算法的具體實(shí)施流程如圖2所示：首先，將觀測(cè)數(shù)據(jù)和模擬數(shù)據(jù)變換至?xí)r頻域，通過(guò)振幅權(quán)重因子在目標(biāo)函數(shù)中調(diào)整振幅相位占比，增強(qiáng)相位信息在成像過(guò)程中所占權(quán)重；然后，利用Gabor逆變換將時(shí)頻域伴隨震源變換至?xí)r間域，再利用伴隨算子把伴隨震源反傳至模型空間，獲得反傳波場(chǎng)；最后，計(jì)算正傳波場(chǎng)與反傳波場(chǎng)的零延遲互相關(guān)，即為模型參數(shù)的更新梯度，并結(jié)合局部?jī)?yōu)化算法對(duì)成像結(jié)果不斷地迭代.

圖1 PA-LSRTM伴隨震源單道時(shí)頻域振幅譜(不同振幅權(quán)重因子) (a) ε=1； (b) ε=0.6； (c) ε=0.3； (d) ε=0.Fig.1 Time-frequency domain amplitude information of the single trace PA-LSRTM adjoint source with different amplitude weight factors

圖2 時(shí)頻域振幅相位聯(lián)合的最小二乘逆時(shí)偏移(PA-LSRTM)流程圖Fig.2 Flow chart of joint least square reverse time migration of phase and amplitude in the time-frequency domain (PA-LSRTM)

2 數(shù)值試驗(yàn)

本文首先利用Marmousi速度模型來(lái)測(cè)試RTM、LSRTM和PA-LSRTM的成像效果.Marmousi速度模型如圖3a所示，經(jīng)過(guò)平滑以后的速度模型如圖3b所示，真實(shí)擾動(dòng)模型如圖3c所示.Marmousi速度模型構(gòu)造復(fù)雜，同時(shí)在深部區(qū)域含有強(qiáng)散射構(gòu)造，增加了深部精細(xì)構(gòu)造的成像難度.Marmousi速度模型的橫向距離為15 km，縱向深度為3.5 km，網(wǎng)格間距為25 m.在地表以300 m間隔均勻分布50炮，每炮對(duì)應(yīng)600個(gè)檢波器全接收.雷克子波震源主頻為8 Hz，地震數(shù)據(jù)觀測(cè)最大時(shí)長(zhǎng)為5 s，時(shí)間采樣間隔為2 ms.

圖4為Marmousi模型成像結(jié)果，與圖3c所示的真實(shí)擾動(dòng)模型對(duì)比可以看出，RTM成像結(jié)果基本能恢復(fù)Marmousi模型的構(gòu)造信息，但是在深部區(qū)域存在明顯的成像不足問(wèn)題，同時(shí)在振幅上與真實(shí)擾動(dòng)模型相差較大.將圖4a、b對(duì)比發(fā)現(xiàn)，LSRTM成像結(jié)果很好地恢復(fù)了Marmousi速度模型的深部構(gòu)造信息，使得深淺區(qū)域成像振幅更加均衡，分辨率也得到了很大的改善.將圖4a、b、c對(duì)比發(fā)現(xiàn)，PA-LSRTM成像結(jié)果在薄細(xì)層位成像精度上有了明顯地提高，同時(shí)成像振幅也更加接近真實(shí)擾動(dòng)模型.此外，在Marmousi模型的強(qiáng)散射構(gòu)造下部區(qū)域(左下側(cè)紅色線框)和儲(chǔ)層構(gòu)造區(qū)域(右下側(cè)紅色線框)的成像分辨率與LSRTM成像結(jié)果相比較，有了較大改善.

為了清晰地看出PA-LSRTM方法在深部精細(xì)構(gòu)造成像優(yōu)勢(shì)，提取圖4所示成像結(jié)果的局部放大圖來(lái)進(jìn)行對(duì)比分析，如圖5所示.可以看出PA-LSRTM成像結(jié)果相比于LSRTM成像結(jié)果在不整合面上的成像有了明顯地改善(左下側(cè)紅色線框)，同時(shí)強(qiáng)散射介質(zhì)底部的傾斜構(gòu)造也得到了清晰地成像.提取了圖4中Marmousi模型右側(cè)含油氣儲(chǔ)層構(gòu)造區(qū)域的成像結(jié)果進(jìn)行對(duì)比分析(右下側(cè)紅色線框)，如圖6所示.可以看出PA-LSRTM成像結(jié)果在儲(chǔ)層位置處的成像分辨率上相比于RTM和LSRTM成像結(jié)果有了明顯的提高.根據(jù)Marmousi模型成像結(jié)果的局部放大圖可以看出，PA-LSRTM方法相比于RTM和LSRTM方法在深部區(qū)域的成像方面有著一定的優(yōu)勢(shì).

圖3 (a) Marmousi速度模型； (b) 偏移速度模型； (c) 真實(shí)擾動(dòng)模型Fig.3 (a) Marmousi velocity model; (b) Migration velocity model; (c) True model perturbations

圖4 Marmousi模型成像結(jié)果 (a) RTM成像結(jié)果； (b) LSRTM成像結(jié)果； (c) PA-LSRTM成像結(jié)果.Fig.4 Imaging results of the Marmousi model (a) RTM result; (b) LSRTM result; (c) PA-LSRTM result.

進(jìn)一步分析RTM、LSRTM和PA-LSRTM的成像結(jié)果的可靠性，利用圖4所示成像結(jié)果進(jìn)行Born正演模擬，震源位置設(shè)置在橫向距離7.5 km處，并將產(chǎn)生的地震數(shù)據(jù)與真實(shí)擾動(dòng)數(shù)據(jù)進(jìn)行對(duì)比(圖7).將圖7a、b對(duì)比可以發(fā)現(xiàn)，RTM的成像結(jié)果對(duì)應(yīng)的模擬數(shù)據(jù)與真實(shí)擾動(dòng)波場(chǎng)相差較大，不僅在深部區(qū)域地震擾動(dòng)波場(chǎng)存在明顯的差異，在淺部區(qū)域RTM成像結(jié)果對(duì)應(yīng)的擾動(dòng)波場(chǎng)同相軸也存在著明顯的噪聲干擾問(wèn)題.將圖7a、c、d對(duì)比可以發(fā)現(xiàn)，LSRTM和PA-LSRTM成像結(jié)果對(duì)應(yīng)的模擬數(shù)據(jù)更接近真實(shí)的擾動(dòng)波場(chǎng).但是在LSRTM成像結(jié)果對(duì)應(yīng)的模擬數(shù)據(jù)中缺失一些弱同相軸信息，這主要是因?yàn)長(zhǎng)SRTM的深部成像結(jié)果與真實(shí)擾動(dòng)模型仍存在一定的差距.將圖7a、d對(duì)比可以發(fā)現(xiàn)，即使觀測(cè)數(shù)據(jù)中的一些微弱同相軸信息，在基于PA-LSRTM成像結(jié)果進(jìn)行Born正演的模擬數(shù)據(jù)中也有與之對(duì)應(yīng)的同相軸信息.同時(shí)提取圖4的單道成像剖面進(jìn)行對(duì)比分析，如圖8所示，可以看出RTM成像振幅與真實(shí)擾動(dòng)模型振幅相差較大，LSRTM成像振幅有了明顯的改善.在圖8深部區(qū)域，可以明顯看到PA-LSRTM的振幅更加接近真實(shí)擾動(dòng)模型振幅，也進(jìn)一步證明了PA-LSRTM方法在弱地震信號(hào)成像方面具有一定的優(yōu)勢(shì).

圖5 Marmousi模型成像結(jié)果局部放大圖(圖4左下側(cè)紅色線框) (a) 真實(shí)擾動(dòng)模型； (b) RTM成像結(jié)果； (c) LSRTM成像結(jié)果； (d) PA-LSRTM成像結(jié)果.Fig.5 Magnified view of Marmousi model imaging results (red rectangular at lower left side in Fig.4) (a) True perturbation model; (b) RTM result; (c) LSRTM result; (d) PA-LSRTM result.

圖6 Marmousi模型成像結(jié)果局部放大圖(圖4右下側(cè)紅色線框) (a) 真實(shí)擾動(dòng)模型； (b) RTM成像結(jié)果； (c) LSRTM成像結(jié)果； (d) PA-LSRTM成像結(jié)果.Fig.6 Magnified view of Marmousi model imaging results (red rectangular at lower right side in Fig.4) (a) True perturbation model; (b) RTM result; (c) LSRTM result; (d) PA-LSRTM result.

圖7 地震擾動(dòng)波場(chǎng)數(shù)據(jù) (a) 真實(shí)擾動(dòng)模型； (b) RTM成像結(jié)果； (c) LSRTM成像結(jié)果； (d) PA-LSRTM成像結(jié)果.Fig.7 Seismic perturbation data (a) True perturbation model; (b) RTM result; (c) LSRTM result; (d) PA-LSRTM result.

圖8 單道成像剖面圖(Marmousi模型11 km處)Fig.8 Single trace imaging profile (Marmousi model at 11 km)

利用鹽丘模型對(duì)RTM、LSRTM和PA-LSRTM方法進(jìn)行數(shù)值測(cè)試，旨在分析鹽丘底界面及其鹽下構(gòu)造的成像精度.真實(shí)鹽丘速度模型如圖9a所示，經(jīng)過(guò)平滑以后的偏移速度模型如圖9b所示，真實(shí)鹽丘擾動(dòng)模型如圖9c所示.鹽丘模型橫向距離為15 km，縱向深度為3.5 km，網(wǎng)格間距為25 m.在地表以300 m間隔均勻分布50炮，每炮對(duì)應(yīng)600個(gè)檢波器全接收.在測(cè)試過(guò)程中，采用震源主頻為8 Hz的雷克子波，地震數(shù)據(jù)的最大記錄時(shí)間為5 s，采樣時(shí)間間隔為2 ms.圖9所示模型中存在強(qiáng)散射鹽丘構(gòu)造，致使大部分地震波能量很難透過(guò)鹽丘體，攜帶鹽下構(gòu)造信息返回地表.此外，鹽下構(gòu)造的反射波至少需要往返兩次穿過(guò)鹽體，致使檢波器接收到鹽下構(gòu)造的地震波響應(yīng)信號(hào)較弱，使得鹽丘底界面及其鹽下精細(xì)構(gòu)造成像精度不足.

圖9 鹽丘模型 (a) 真實(shí)鹽丘速度模型； (b) 偏移速度模型； (c) 真實(shí)擾動(dòng)模型.Fig.9 Salt model (a) True salt velocity model; (b) Migration velocity model; (c) True perturbations model.

圖10 鹽丘模型成像結(jié)果 (a) RTM成像結(jié)果； (b) LSRTM成像結(jié)果； (c) PA-LSRTM成像結(jié)果.Fig.10 Imaging results of the salt model (a) RTM result; (b) LSRTM result; (c) PA-LSRTM result.

圖10分別展示了RTM、LSRTM和PA-LSRTM的鹽丘模型成像結(jié)果.從圖10a中可以看出，RTM成像結(jié)果淺部能量較強(qiáng).這主要是因?yàn)榛诓▌?dòng)方程的互相關(guān)成像條件受到震源的影響較為嚴(yán)重.與此同時(shí)，RTM成像結(jié)果在深部區(qū)域存在明顯的成像不足等問(wèn)題，這主要是因?yàn)榈卣鸩芰吭趥鞑ミ^(guò)程中快速衰減，深部反射的地震信號(hào)較弱，很難獲得高質(zhì)量的成像結(jié)果.將圖10a、b對(duì)比發(fā)現(xiàn)，LSRTM成像結(jié)果很好地呈現(xiàn)出了鹽丘速度模型的深部構(gòu)造信息，使得整個(gè)鹽丘模型從淺部到深部成像振幅更加均衡，鹽丘構(gòu)造的成像分辨率得到了較大的改善.將圖10a、b、c對(duì)比發(fā)現(xiàn)，PA-LSRTM方法在鹽丘下部區(qū)域精細(xì)構(gòu)造成像分辨率上有了較大提高.此外，PA-LSRTM方法通過(guò)增強(qiáng)相位權(quán)重法，減弱振幅影響，有效緩解了震源處的強(qiáng)能量對(duì)淺部區(qū)域成像的干擾，使得PA-LSRTM成像結(jié)果從淺層到深層振幅更加均衡.

接著提取鹽丘模型成像結(jié)果的局部放大圖進(jìn)行分析，如圖11和圖12所示.圖11是鹽下構(gòu)造成像結(jié)果左下側(cè)的局部放大對(duì)比圖，該區(qū)域主要是分析鹽下精細(xì)構(gòu)造的成像精度.根據(jù)圖11的對(duì)比結(jié)果可以看出，PA-LSRTM成像結(jié)果相比于RTM和LSRTM成像結(jié)果在鹽丘的下界面及其鹽下構(gòu)造成像更接近真實(shí)擾動(dòng)模型.圖12是鹽丘模型右下側(cè)的局部放大對(duì)比圖，相比于RTM和LSRTM成像結(jié)果，可以看出PA-LSRTM成像結(jié)果在深部區(qū)域精細(xì)構(gòu)造成像分辨率上有了很大的提高.圖13是單道成像剖面圖，可以看出RTM成像振幅與真實(shí)擾動(dòng)模型振幅相差較大，尤其是在深部區(qū)域，RTM幾乎無(wú)法獲得鹽下深部構(gòu)造信息.與RTM成像結(jié)果相比，LSRTM

圖11 鹽丘模型成像結(jié)果局部放大圖(圖10左下側(cè)紅色線框) (a) 真實(shí)擾動(dòng)模型； (b) RTM成像結(jié)果； (c) LSRTM成像結(jié)果； (d) PA-LSRTM成像結(jié)果.Fig.11 Magnified view of salt model imaging results (red rectangular at lower left side in Fig.10) (a) True perturbation model; (b) RTM result; (c) LSRTM result; (d) PA-LSRTM result.

圖12 鹽丘模型成像結(jié)果局部放大圖(圖10右下側(cè)紅色線框) (a) 真實(shí)擾動(dòng)模型； (b) RTM成像結(jié)果； (c) LSRTM成像結(jié)果； (d) PA-LSRTM成像結(jié)果.Fig.12 Magnified view of salt model imaging results (red rectangular at lower right side in Fig.10) (a) True perturbation model; (b) RTM result; (c) LSRTM result; (d) PA-LSRTM result.

圖13 單道成像剖面圖(鹽丘模型6.25 km處)Fig.13 Single trace imaging profile (salt model at 6.25 km)

成像振幅有了較大的提高，但是在深部區(qū)域的成像振幅仍然與真實(shí)擾動(dòng)模型有一定的差距.在深度為3 km處(圖13)，PA-LSRTM成像振幅最接近真實(shí)擾動(dòng)模型振幅.綜合鹽丘模型成像結(jié)果、局部放大對(duì)比圖及其單道成像剖面對(duì)比圖可以看出，PA-LSRTM方法在透過(guò)上覆強(qiáng)散射地質(zhì)體(鹽丘)對(duì)鹽下構(gòu)造成像方面具有一定優(yōu)勢(shì).

圖14 含噪擾動(dòng)波場(chǎng)(信噪比SNR=-3.59)Fig.14 Perturbation wavefield with Gaussian noise (SNR=-3.59)

圖15 含噪PA-LSRTM成像結(jié)果Fig.15 PA-LSRTM result with Gaussian noise

地震數(shù)據(jù)中相位信息相比于振幅信息具有與地下介質(zhì)更好的線性對(duì)應(yīng)關(guān)系(Fichtner et al.，2008；Djebbi and Alkhalifah，2014；胡勇等，2018)，但是純相位目標(biāo)函數(shù)容易受到噪聲的干擾，為此本文提出振幅相位聯(lián)合的最小二乘逆時(shí)偏移成像方法(PA-LSRTM).為了測(cè)試該方法的抗噪性，在地震數(shù)據(jù)中加入了較強(qiáng)高斯噪聲，含噪單炮記錄如圖14所示，其中深部區(qū)域的弱散射信號(hào)幾乎被噪聲掩蓋，這給深部區(qū)域成像帶來(lái)了一定的困難.在負(fù)信噪比的情況下(SNR=-3.59)，鹽丘模型的PA-LSRTM成像結(jié)果如圖15所示.將圖15與圖10c對(duì)比可以發(fā)現(xiàn)，鹽丘內(nèi)部區(qū)域成像結(jié)果受到噪聲影響較為嚴(yán)重，但依然能夠很好地對(duì)鹽丘邊界及其鹽下精細(xì)構(gòu)造進(jìn)行成像.這主要是因?yàn)檎穹鶛?quán)重因子開(kāi)始設(shè)定為ε=0.6，然后逐漸減小至ε=0.3，目標(biāo)函數(shù)中的這一部分振幅信息很好地緩解了噪聲相位對(duì)成像結(jié)果的干擾.因此，在PA-LSRTM目標(biāo)函數(shù)中重新分配振幅和相位的占比，既能很好的利用相位信息對(duì)深部構(gòu)造進(jìn)行成像，同時(shí)還能在一定程度上緩解噪聲相位對(duì)成像結(jié)果的干擾.從含強(qiáng)高斯噪聲的測(cè)試結(jié)果中可以看出，PA-LSRTM適合針對(duì)深部區(qū)域、強(qiáng)散射構(gòu)造、及其鹽下構(gòu)造進(jìn)行精細(xì)成像，同時(shí)該方法還具有一定的抗噪性.

3 結(jié)論與認(rèn)識(shí)

針對(duì)常規(guī)最小二乘逆時(shí)偏移鹽下精細(xì)構(gòu)造成像困難的問(wèn)題，本文提出時(shí)頻域振幅相位聯(lián)合的最小二乘逆時(shí)偏移方法，結(jié)合理論分析和模型測(cè)試得出以下幾點(diǎn)結(jié)論與認(rèn)識(shí)：

(1)LSRTM方法很好地解決了RTM成像振幅不均衡的問(wèn)題，同時(shí)大大提高了成像分辨率.但是由于LSRTM方法受振幅影響較為嚴(yán)重，因此很難獲得鹽下構(gòu)造等深部區(qū)域的高精度成像結(jié)果.

(2)時(shí)頻域目標(biāo)函數(shù)能夠很好地分離了地震信號(hào)的振幅和相位信息，通過(guò)引入振幅權(quán)重因子，調(diào)節(jié)振幅和相位信息在目標(biāo)函數(shù)中的權(quán)重，弱化振幅對(duì)成像結(jié)果的影響.在此基礎(chǔ)上，PA-LSRTM方法很好地實(shí)現(xiàn)了深部弱地震信號(hào)的高精度成像.

(3)純相位目標(biāo)函數(shù)容易受到噪聲的干擾，本文通過(guò)在時(shí)頻域目標(biāo)函數(shù)中聯(lián)合使用振幅和相位信息，在提高地下深部區(qū)域精細(xì)構(gòu)造的成像質(zhì)量的情況下，同時(shí)也保證了PA-LSRTM算法的抗噪性和穩(wěn)定性.

最后，Marmousi模型和鹽丘模型的數(shù)值測(cè)試結(jié)果表明，PA-LSRTM方法能夠很好地利用弱地震信號(hào)的時(shí)頻域振幅相位信息，實(shí)現(xiàn)透過(guò)上覆強(qiáng)散射地層獲得深部構(gòu)造高精度成像的目標(biāo).

附錄A

時(shí)頻域目標(biāo)函數(shù)對(duì)應(yīng)的振幅譜和相位譜，可以通過(guò)Gabor時(shí)頻變換獲得，其中觀測(cè)數(shù)據(jù)和模擬數(shù)據(jù)的Gabor變換為

(A1)

(A2)

(A3)

(A4)

附錄B

根據(jù)波動(dòng)方程的表達(dá)式，在計(jì)算擾動(dòng)波場(chǎng)的過(guò)程中，Born正演模擬也可以用矩陣的形式來(lái)表示：

Asus=A0u0,

(B1)

(B2)

由于As和u0與模型參數(shù)無(wú)關(guān)，則式(B2)變?yōu)?/p>

(B3)

因此，LSRTM的目標(biāo)函數(shù)相對(duì)于模型參數(shù)的偏導(dǎo)數(shù)可以表示為

(B4)

(B5)

因此，LSRTM和PA-LSRTM目標(biāo)函數(shù)對(duì)模型參數(shù)的梯度可以表示為正傳波場(chǎng)與反傳波場(chǎng)的零延遲互相關(guān).

附錄C

PA-LSRTM目標(biāo)函數(shù)是利用模擬數(shù)據(jù)的時(shí)頻域振幅相位信息來(lái)不斷地與觀測(cè)數(shù)據(jù)的時(shí)頻域振幅相位信息相匹配，則PA-LSRTM的目標(biāo)函數(shù)可以定義為

(C1)

(C2)

(C3)

將式(C3)代入式(C2)中，則目標(biāo)函數(shù)對(duì)模型參數(shù)的偏導(dǎo)數(shù)可以表示為

(C4)

(C5)

將式(C5)代入式(C4)中，則目標(biāo)函數(shù)對(duì)模型參數(shù)的梯度算子可以表示為

(C6)

(C7)

將式(A3)代入式(C7)中，目標(biāo)函數(shù)對(duì)模型參數(shù)的梯度可以簡(jiǎn)化為

(C8)

則伴隨震源可以定義為

(C9)

因此，PA-LSRTM目標(biāo)函數(shù)對(duì)模型參數(shù)的梯度與LSRTM梯度計(jì)算方法相同，只需將伴隨震源反傳至模型空間獲得的反傳波場(chǎng)，并與正傳波場(chǎng)做零延遲互相關(guān).

附錄D

利用L-BFGS優(yōu)化算法計(jì)算模型參數(shù)的更新方向，其迭代公式為

mk+1=mk-αkΔm，

(D1)

其中mk為第k步模型參數(shù)；αk為步長(zhǎng)；Δm為模型參數(shù)的更新方向.關(guān)于LSRTM優(yōu)化過(guò)程，實(shí)際上就是尋找目標(biāo)函數(shù)最小值的過(guò)程.假設(shè)目標(biāo)函數(shù)的起始點(diǎn)與最小值點(diǎn)在一個(gè)微小的鄰域內(nèi)，則LSRTM目標(biāo)函數(shù)的最小值點(diǎn)存在如下關(guān)系：

(D2)

舍去式(D2)中的高階項(xiàng)：

(D3)

當(dāng)模型擾動(dòng)量較小的情況下(m+δm≈m)，則模型參數(shù)的更新方向可以表示為

(D4)

(D5)

(D6)

sk=mk+1-mk，yk=gk+1-gk，

(D7)

其中Hk+1是根據(jù)向量對(duì){sk,yk}和Hk計(jì)算得到；Hkgk的乘積可以通過(guò)梯度gk與向量對(duì){sk,yk}之間一系列向量的內(nèi)積與向量的和來(lái)獲得的.其中的近似Hessian矩陣的逆Hk需滿足以下更新公式：

(D8)

因此，在L-BFGS 優(yōu)化算法中只需要保存少數(shù)的向量對(duì)，即可獲得Hessian 矩陣逆的近似，大大提高了計(jì)算效率，同時(shí)很大程度上節(jié)約了計(jì)算機(jī)的內(nèi)存.