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

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

0 引言

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

為了進(jìn)一步提高逆時偏移的成像精度，發(fā)展了最小二乘逆時偏移方法(LSRTM)，其在反演的理論框架下，利用反偏移數(shù)據(jù)來不斷地與觀測數(shù)據(jù)相匹配，并通過優(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)行測試，很大程度上提高了LSRTM的計算效率.任浩然等(2013)將Hessian算子應(yīng)用于LSRTM中，獲得了相對保真的成像結(jié)果，改善了地震偏移的成像精度.Tan和Huang(2014)在原始LSRTM的基礎(chǔ)上，提出了波場分離成像條件，并通過更新震源波場實(shí)現(xiàn)陡傾角斷層成像.此外，為了克服LSRTM對偏移模型的依賴問題，劉玉金和李振春(2015)提出擴(kuò)展成像條件下的LSRTM方法，測試結(jié)果表明該方法可以獲得更準(zhǔn)確的振幅屬性信息.Huang等(2015)、周東紅等(2020)提出了帶地形的LSRTM方法，進(jìn)一步考慮地表起伏對波動方程成像結(jié)果的影響，指出無需對地震數(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ù)，在稀疏約束的理論框架下展開研究.Liu等(2020)在LSRTM中加入了Gabor反褶積濾波來提高鹽下構(gòu)造的成像精度.Li等(2020)通過引入權(quán)重因子來衰減地震數(shù)據(jù)中較強(qiáng)的反射波信息，進(jìn)而增強(qiáng)弱散射信號，實(shí)現(xiàn)改善深部構(gòu)造成像精度的目標(biāo).經(jīng)過多年的研究，LSRTM技術(shù)得到了快速地發(fā)展，逐漸成為當(dāng)前的研究熱點(diǎn)問題.但是地震波傳播是一個較為復(fù)雜的過程，使得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)透過上覆強(qiáng)散射地質(zhì)體來獲得深部精細(xì)構(gòu)造的高精度成像目標(biāo)，仍有待進(jìn)一步研究.

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

1 方法原理

1.1 線性化Born正演

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

(1)

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

(2)

u=u0+us.

(3)

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

(4)

+f，

(5)

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

(6)

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

1.2 最小二乘逆時偏移原理

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

(7)

(8)

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

(9)

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

(10)

(11)

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

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

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

(12)

(13)

(14)

式(14)可以簡化為(詳見附錄C)：

(15)

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

(16)

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

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

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

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

圖1 PA-LSRTM伴隨震源單道時頻域振幅譜(不同振幅權(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 時頻域振幅相位聯(lián)合的最小二乘逆時偏移(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速度模型來測試RTM、LSRTM和PA-LSRTM的成像效果.Marmousi速度模型如圖3a所示，經(jīng)過平滑以后的速度模型如圖3b所示，真實(shí)擾動模型如圖3c所示.Marmousi速度模型構(gòu)造復(fù)雜，同時在深部區(qū)域含有強(qiáng)散射構(gòu)造，增加了深部精細(xì)構(gòu)造的成像難度.Marmousi速度模型的橫向距離為15 km，縱向深度為3.5 km，網(wǎng)格間距為25 m.在地表以300 m間隔均勻分布50炮，每炮對應(yīng)600個檢波器全接收.雷克子波震源主頻為8 Hz，地震數(shù)據(jù)觀測最大時長為5 s，時間采樣間隔為2 ms.

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

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

圖3 (a) Marmousi速度模型； (b) 偏移速度模型； (c) 真實(shí)擾動模型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í)擾動數(shù)據(jù)進(jìn)行對比(圖7).將圖7a、b對比可以發(fā)現(xiàn)，RTM的成像結(jié)果對應(yīng)的模擬數(shù)據(jù)與真實(shí)擾動波場相差較大，不僅在深部區(qū)域地震擾動波場存在明顯的差異，在淺部區(qū)域RTM成像結(jié)果對應(yīng)的擾動波場同相軸也存在著明顯的噪聲干擾問題.將圖7a、c、d對比可以發(fā)現(xiàn)，LSRTM和PA-LSRTM成像結(jié)果對應(yīng)的模擬數(shù)據(jù)更接近真實(shí)的擾動波場.但是在LSRTM成像結(jié)果對應(yīng)的模擬數(shù)據(jù)中缺失一些弱同相軸信息，這主要是因?yàn)長SRTM的深部成像結(jié)果與真實(shí)擾動模型仍存在一定的差距.將圖7a、d對比可以發(fā)現(xiàn)，即使觀測數(shù)據(jù)中的一些微弱同相軸信息，在基于PA-LSRTM成像結(jié)果進(jìn)行Born正演的模擬數(shù)據(jù)中也有與之對應(yīng)的同相軸信息.同時提取圖4的單道成像剖面進(jìn)行對比分析，如圖8所示，可以看出RTM成像振幅與真實(shí)擾動模型振幅相差較大，LSRTM成像振幅有了明顯的改善.在圖8深部區(qū)域，可以明顯看到PA-LSRTM的振幅更加接近真實(shí)擾動模型振幅，也進(jìn)一步證明了PA-LSRTM方法在弱地震信號成像方面具有一定的優(yōu)勢.

圖5 Marmousi模型成像結(jié)果局部放大圖(圖4左下側(cè)紅色線框) (a) 真實(shí)擾動模型； (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í)擾動模型； (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 地震擾動波場數(shù)據(jù) (a) 真實(shí)擾動模型； (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)

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

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

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

圖11 鹽丘模型成像結(jié)果局部放大圖(圖10左下側(cè)紅色線框) (a) 真實(shí)擾動模型； (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í)擾動模型； (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í)擾動模型有一定的差距.在深度為3 km處(圖13)，PA-LSRTM成像振幅最接近真實(shí)擾動模型振幅.綜合鹽丘模型成像結(jié)果、局部放大對比圖及其單道成像剖面對比圖可以看出，PA-LSRTM方法在透過上覆強(qiáng)散射地質(zhì)體(鹽丘)對鹽下構(gòu)造成像方面具有一定優(yōu)勢.

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

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

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

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

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

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

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

附錄A

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

(A1)

(A2)

(A3)

(A4)

附錄B

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

Asus=A0u0,

(B1)

(B2)

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

(B3)

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

(B4)

(B5)

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

附錄C

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

(C1)

(C2)

(C3)

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

(C4)

(C5)

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

(C6)

(C7)

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

(C8)

則伴隨震源可以定義為

(C9)

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

附錄D

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

mk+1=mk-αkΔm，

(D1)

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

(D2)

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

(D3)

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

(D4)

(D5)

(D6)

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

(D7)

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

(D8)

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