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

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

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

0 引言

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

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

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

1 方法原理

1.1 線性化Born正演

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

(1)

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

(2)

u=u0+us.

(3)

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

(4)

+f，

(5)

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

(6)

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

1.2 最小二乘逆時偏移原理

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

(7)

(8)

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

(9)

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

(10)

(11)

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

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

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

(12)

(13)

(14)

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

(15)

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

(16)

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

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

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

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

圖1 PA-LSRTM伴隨震源單道時頻域振幅譜(不同振幅權重因子) (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ù)值試驗

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

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

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

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

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

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

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

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

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

圖10分別展示了RTM、LSRTM和PA-LSRTM的鹽丘模型成像結果.從圖10a中可以看出，RTM成像結果淺部能量較強.這主要是因為基于波動方程的互相關成像條件受到震源的影響較為嚴重.與此同時，RTM成像結果在深部區(qū)域存在明顯的成像不足等問題，這主要是因為地震波能量在傳播過程中快速衰減，深部反射的地震信號較弱，很難獲得高質(zhì)量的成像結果.將圖10a、b對比發(fā)現(xiàn)，LSRTM成像結果很好地呈現(xiàn)出了鹽丘速度模型的深部構造信息，使得整個鹽丘模型從淺部到深部成像振幅更加均衡，鹽丘構造的成像分辨率得到了較大的改善.將圖10a、b、c對比發(fā)現(xiàn)，PA-LSRTM方法在鹽丘下部區(qū)域精細構造成像分辨率上有了較大提高.此外，PA-LSRTM方法通過增強相位權重法，減弱振幅影響，有效緩解了震源處的強能量對淺部區(qū)域成像的干擾，使得PA-LSRTM成像結果從淺層到深層振幅更加均衡.

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

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

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

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

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

3 結論與認識

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

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

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

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

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

附錄A

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

(A1)

(A2)

(A3)

(A4)

附錄B

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

Asus=A0u0,

(B1)

(B2)

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

(B3)

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

(B4)

(B5)

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

附錄C

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

(C1)

(C2)

(C3)

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

(C4)

(C5)

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

(C6)

(C7)

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

(C8)

則伴隨震源可以定義為

(C9)

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

附錄D

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

mk+1=mk-αkΔm，

(D1)

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

(D2)

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

(D3)

當模型擾動量較小的情況下(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é)約了計算機的內(nèi)存.