孫昌瀟， 毛偉建， 張慶臣， 石星辰

1 中國科學院精密測量科學與技術創(chuàng)新研究院計算與勘探地球物理研究中心； 大地測量與地球動力學國家重點實驗室， 武漢 430077 2 中國科學院大學， 北京 100049

0 引言

地震疊前深度偏移成像是揭示地下油氣藏構造信息的關鍵技術，傳統的偏移方法通過求解正演算子的伴隨算子，將地表接收點采集到的時間域地震波場反傳到地下來獲得散射點的位置.由于偏移算子不是正演算子的逆，并且受到采集系統局限、深部地層構造復雜、實際地震數據通常伴有噪聲等因素的影響，傳統的偏移方法通常只能得到一個模糊的成像，分辨率較低.反演是獲取地下物性參數的重要手段，并且隨著計算機技術的飛速發(fā)展，反演正在越來越多地應用到地震數據處理中，能夠為地震資料解釋提供強有力的證據.最小二乘偏移聯合偏移和反演，依據線性化反演的思路(Tarantola，1984)，將成像作為一個反問題進行求解，在最小二乘目標函數約束下，修正偏移成像的振幅，反演地層真實的反射信息.

Nemeth等(1999)提出最小二乘偏移成像方法，用于減少成像中采集腳印的影響和迭代恢復缺失地震數據.隨后，最小二乘偏移得到了國內外地球物理學家廣泛的研究和應用.通過在同樣觀測系統下模擬數據與實際觀測數據的迭代更新擬合，最小二乘偏移可以解決傳統構造成像分辨率低、照明不均衡的問題(Chavent and Plessix，1999；Duquet et al.，2000；Aoki and Schuster，2009；Dai and Schuster，2013；馬方正等，2016)，同時可以壓制成像串擾(Dai et al.，2011，2012；劉玉金等，2013；Dutta，2017；Zhang et al.，2019；張攀和毛偉建，2018).但是最小二乘偏移往往需要幾十次迭代，每次迭代都包含一次從震源到地下散射點的正向傳播，和一次來自接收點位置的反向傳播，計算時間是傳統偏移方法的2n(n為迭代次數)倍，因此大量的計算耗時是需要解決的問題.

高斯束偏移基于高斯束模擬非奇異、振幅處處正則的地震波場，是一種靈活、準確、高效的深度域偏移成像方法.它兼具射線類和波動方程類偏移成像的優(yōu)勢，能夠解決兩點射線追蹤存在的多路徑的問題(Hill，1990，2001；Gray，2005；Gray and Bleistein，2009)，適用于三維復雜地質構造的成像.20世紀90年代初，Hill(1990)提出高斯束偏移成像的思路，給出了高斯束的計算公式，以及基于高斯束疊加的波場延拓公式，用于疊后偏移.在此之后，高斯束偏移經歷了從疊后到疊前、2D到3D、聲波均勻各向同性介質到復雜各向異性介質中的高斯束偏移(Popov et al.，2010 ；Protasov，2015；Li et al.，2018；Yang et al.，2018b).最小二乘高斯束偏移結合最小二乘和高斯束偏移成像的優(yōu)勢，能夠提高成像分辨率，均衡成像振幅(Hu et al.，2016；Yuan et al.，2017).Yang等(2018a)提出時間域最小二乘高斯束偏移成像的方法，推導了時間域高斯束Born正演公式和偏移公式，并且運用正則化約束方法增強反演的穩(wěn)定性.Yue等(2019)將最小二乘高斯束偏移方法應用到二維彈性介質成像，提高了PP和PS波成像的分辨率.目前最小二乘高斯束偏移的相關研究與討論較少，并且還未涉及三維介質中最小二乘高斯束偏移的應用.最小二乘高斯束偏移具有靈活度高、計算效率快等優(yōu)勢，因此該方法值得我們研究和開發(fā).

本文提出了一種最小二乘彈性高斯束偏移成像方法，并把它應用于三維彈性波成像.首先，在炮點和檢波點附近稀疏位置分別向下進行基于高斯束疊加表示的多波型波場延拓，三維空間下表達為P波、S1波和S2波三種波型的波場延拓形式.然后根據高斯束的走時和振幅信息，將三維地下介質的擾動點能量映射到地表稀疏位置并得到多波型局部平面波.利用不同的應力邊界條件下，多波和多分量地震記錄之間的關系，構建波型波矢量轉換矩陣(栗學磊和毛偉建，2016)，合成多分量局部平面波.然后在地表劃分一系列重疊排列的高斯窗，應用局部逆傾斜疊加公式(岳玉波等，2019a，2019b)，在接收點合成多分量地震記錄.偏移算子取Born正演(反偏移)算子的伴隨算子，采用共軛梯度算法更新成像結果來擬合反偏移數據與實際觀測數據，使成像振幅逐漸趨近于地下真實的反射系數，從而提高偏移成像的分辨率，增強地下照明.本文通過構建的波型波矢量轉換矩陣來進行每次迭代過程中多波型、多分量局部平面波之間的轉換，從而降低最小二乘偏移PP和PS成像串擾.為了更好地平衡PP和PS偏移反演的結果，我們在反演過程中分別控制修正PP和PS成像的步長，使迭代過程穩(wěn)定收斂.通過模型測試表明，最小二乘彈性高斯束偏移方法能夠有效提高成像分辨率.復雜Marmousi2模型的測試結果表明對于陡傾角構造和深層射線覆蓋率低、成像振幅弱的區(qū)域，最小二乘彈性高斯束偏移可以提高成像保幅能力.從三維彈性波成像的縱向剖面和橫向截面上看，最小二乘彈性高斯束偏移能夠得到更為清晰的成像結果.

1 工作原理

1.1 彈性高斯束Born正演

根據一階Born近似模擬散射波場的理論(Beylkin and Burridege，1990；Bleistein et al.，2001)，在彈性介質中，由震源s激發(fā)，經地下介質傳播到r處接收的vr波型彈性波場uvr(r,s,ω)表示為

(1)

定義激發(fā)震源為爆炸源，因此vs代表P波型，vr代表P、S1或S2波型(栗學磊和毛偉建，2016)；S(ω)為頻率域震源函數；M(x)表示x點的反射率，它與介質的彈性參數有關；格林函數Gvsvr(r,s,ω)表示在地表r位置接收到的由s位置激發(fā)vs波型單位源的響應.將公式(1)中Gvsvr(r,s,ω)表示為Gvs(x,s,ω)Gvr(r,x,ω)的形式，并且根據互易性定理Gvr(r,x,ω)=Gvr(x,r,ω)，公式(1)彈性波場uvr(r,s,ω)的表達式變?yōu)?/p>

(2)

因此，一階Born正演彈性波場表達為震源到地下散射點的格林函數、檢波點位置處接收的格林函數、頻率域震源函數，以及反射率在頻率域的積分.

彈性高斯束Born正演的格林函數表示為彈性高斯束疊加的形式(Hill，2001)：

(3)

高斯束Born正演過程是將地下散射點能量映射到地表合成局部平面波，再在地表劃分一系列重疊排列的高斯窗，將(p,ω)域數據變換到(r,ω)域，該過程可以看作局部傾斜疊加的逆過程.Hill(2001)給出了高斯窗劃分的歸一化表達式：

(4)

其中，常數a表示相鄰兩個高斯束的間隔，ωl和wl分別表示參考頻率和初始束寬度.將公式(4)代入彈性高斯束疊加表示的檢波點格林函數Gvr(x,r,ω)中，并且引入相移校正因子exp[-iωp·(r-L)]，則

(5)

(6)

(7)

將公式(5)、(6)代入一階Born近似公式(2)中，并且將震源格林函數也表達為彈性高斯束疊加的形式，則震源s點激發(fā)地震波場，到r處接收的彈性波矢量場un(r,s,ω)的表達式為

其中Dvr(L,s,p,ω)表示地表稀疏高斯束中心位置合成的vr波型局部平面波，具體表達式如下：

(9)

1.2 最小二乘彈性高斯束偏移

最小二乘偏移的本質在于運用偏移的手段，通過接收點采集到的地震記錄反演背景地球物理模型下未知散射點或擾動點的值，其中重要的環(huán)節(jié)是計算高斯束Born正演(反偏移)算子和它的伴隨算子，即高斯束偏移算子.對一階Born近似正演過程求共軛轉置，偏移過程表示為

(10)

Gvs*(x,s,ω)、Gvr*(x,r,ω)和S*(ω)分別表示震源格林函數、檢波點格林函數和震源函數的復數共軛.格林函數表示為彈性高斯束疊加的形式，則

(11)

(12)

假定l表示彈性高斯束Born正演(反偏移)算子，彈性高斯束偏移算子用lT表示，則彈性高斯束反偏移過程可以用算子l作用于地下介質擾動模型m的形式表示.依照線性化反演的思路，構建最小二乘框架下的目標函數：

φ(m)=‖d-lm‖2，

(13)

式中，d表示地表接收點采集到的實際地震數據.因此，最小二乘偏移的目的在于通過正演模擬地震數據來擬合實際的觀測地震數據.上式中，令?φ(m)/?m=0，整理得到

m=(lTl)-1lTd，

(14)

其中，lTl稱為Hessian矩陣.由于嚴格意義上求解Hessian矩陣計算量大、成本高，通常研究Hessian矩陣逆的近似作為預條件，來幫助反演地下介質的反射率(Guitton，2004；Tang，2009；Ayeni and Biondi，2010；任浩然等，2013).

1.3 共軛梯度反演算法

用梯度導向的迭代算法求解目標函數：

mk+1=mk+αgk+1，

(15)

其中，k表示迭代次數；α表示步長；gk+1表示第k+1次迭代的目標梯度；mk和mk+1分別表示第k次和第k+1次反演的地下介質擾動，在最小二乘彈性高斯束偏移成像中，mk和mk+1表示第k次和第k+1次迭代的PP和PS成像結果.從表達式中可以得出，反演的準確性是依據模擬數據與實際觀測數據的擬合程度來判定.

共軛梯度算法通過加入一個對梯度的修正，使反演過程更加穩(wěn)定收斂.假設初始模型m0=0，第一次迭代的成像結果m1=αlTd，相當于對觀測數據做一次偏移成像.第k+1(k=1，2，3…)次迭代得到偏移成像mk+1的計算流程如下：

(16)

其中sk+1和β分別表示第k+1次迭代的共軛梯度和它的步長.

最小二乘彈性高斯束偏移需要同時反演PP和PS成像，在數據域表現為多分量地震數據的擬合.為了更好地控制PP和PS成像反演的速率和下降方向，本文采用了一種分別擬合PP成像和PS成像反偏移多分量模擬數據與實際多分量觀測數據的方法.在每次迭代過程中，分別控制PP和PS成像修正的步長，記作αPP和αPS，該方法使PP和PS成像反演更穩(wěn)定.

2 模型試算

為了驗證最小二乘彈性高斯束偏移成像方法的有效性，分別采用二維凹陷、Marmousi2模型以及三維SEG/EAGE Salt模型進行測試，實際地震數據通過有限差分方法正演得到.分別對比了傳統彈性高斯束偏移方法和最小二乘彈性高斯束偏移方法的成像結果，并且分析了成像頻譜變化和目標函數的收斂情況.

2.1 2D凹陷模型

圖1 凹陷模型(a) P波速度； (b) PP反射率.Fig.1 The sunken model(a) P wave velocity； (b) PP reflectivity.

圖2 不同方法正演的單炮地震記錄有限差分正演z分量(a)、x分量(c)數據；彈性高斯束Born正演z分量(b)、x分量(d)數據.Fig.2 The synthetic data using different modeling methodsz-component (a) and x-component (c) data using finite difference forward modeling method；z-component (b) and x-component (d) data using elastic Gaussian beam Born modeling method.

圖3表示傳統的彈性高斯束偏移成像與最小二乘彈性高斯束偏移成像的對比結果.如該圖顯示，傳統彈性高斯束偏移PP和PS成像(圖3a、3c)的分辨率較低，照明不均衡，并且由于觀測系統的局限性，凹陷兩側傾斜斷層構造的射線覆蓋率低，成像振幅比較弱. 10次迭代后最小二乘彈性高斯束偏移PP、PS成像如圖3b、3d所示，結果表明，最小二乘偏移成像的分辨率有了明顯的提高，地下照明不足或者照明不到的地方也得到了一定的補償，斷層構造地區(qū)的成像振幅得到增強.

圖3 (a)和(c)是傳統彈性高斯束PP和PS波偏移結果，(b)和(d)是最小二乘彈性高斯束PP和PS波偏移結果Fig.3 (a) and (c) are conventional elastic Gaussian beam migration results for PP and PS waves;(b) and (d) least-squares elastic Gaussian beam migration results for PP and PS waves

圖4表示同一個尺度范圍下顯示的觀測數據與反偏移模擬數據殘差，1次迭代和10次迭代后的z分量數據殘差如圖4a、4b所示，x分量數據殘差如圖4c、4d所示.從圖中可以看出，迭代10次后x分量和z分量的數據殘差都明顯減小.圖4b、4d中顯示仍然有很小一部分數據殘差未隨著迭代次數的增加而減小，這是因為反偏移無法合成超過射線覆蓋范圍以外的地震記錄.為了進一步比較目標函數隨迭代次數的收斂情況，我們給出隨迭代次數增加的歸一化目標函數收斂曲線，如圖5.經過10次迭代后，目標函數收斂為第1次迭代的20%左右.

圖4 觀測和模擬單炮地震記錄的殘差1次迭代后z分量(a)、x分量(c)的數據殘差；10次迭代后z分量(b)、x分量(d)的數據殘差.Fig.4 The data residuals between simulated and observed dataz-component (a) and x-component (c) residuals at the 1st iteration；z-component (b) and x-component (d) residuals at the 10th iteration.

圖5 歸一化的目標函數收斂曲線Fig.5 The convergence curve of the normalized objective function

2.2 2D Marmousi2模型

圖6 Marmousi2模型(a) P波速度； (b) PP反射率.Fig.6 The Marmousi2 model(a) P wave velocity； (b) PP reflectivity.

圖7中分別給出了傳統彈性高斯束偏移PP(圖7a)和PS(圖7c)成像剖面以及10次迭代后的最小二乘彈性高斯束偏移PP(圖7b)和PS(圖7d)成像剖面.從圖中可以觀察到，傳統偏移方法得到的PP成像剖面照明不均衡，特別是在地層深部和地層結構變化劇烈的區(qū)域，照明度明顯地削弱，并且成像分辨率低，即使PS成像的分辨率比PP成像的分辨率要高，但也存在照明不均衡的問題.相比傳統偏移PP和PS成像，最小二乘彈性高斯束偏移方法提高了PP和PS成像分辨率，并且能夠一定程度上均衡地層照明，補償深部弱信號，增強成像的保幅性.然而隨著迭代次數的增加，最小二乘彈性高斯束偏移成像受到部分假象和噪聲影響信噪比降低，這是因為目標函數解的非光滑性使得迭代后的成像結果常伴有一定的低頻噪聲.并且由于最小二乘偏移是一種線性化反演方法，復雜的地層構造增加了地震數據擬合的難度，進而一定程度上影響成像反演的穩(wěn)定性.

圖7 (a)和(c)是傳統彈性高斯束PP和PS波偏移結果; (b)和(d)是最小二乘彈性高斯束PP和PS波偏移結果Fig.7 (a) and (c) are conventional elastic Gaussian beam migration results for PP and PS waves;(b) and (d) least-squares elastic Gaussian beam migration results for PP and PS waves

我們選取模型中結構變化劇烈的部分(圖6b黑色方框所示)放大進行更清晰地對比分析.從圖8虛線方框部分的成像效果來看，與傳統彈性高斯束偏移相比，最小二乘彈性高斯束偏移方法的成像分辨率更高，地層結構更清晰.另外，在圖中黑色箭頭所指部分，由于地層結構復雜，傳統的高斯束偏移不足以對該區(qū)域進行精確成像，而最小二乘彈性高斯束偏移能夠補償該區(qū)域的照明，提高成像的保幅性.

圖8 圖6b中黑色方框部分的偏移結果比較(a)、(b) PP、PS反射率； (c)、(d) 傳統的彈性高斯束PP、PS偏移； (e)、(f) 最小二乘彈性高斯束PP、PS偏移.Fig.8 The comparison among the migration results in the black box of Fig.6b(a) and (b) PP and PS reflectivity； (c) and (d) Conventional elastic Gaussian beam migration for PP and PS images； (e) and (f) Least-squares elastic Gaussian beam migration for PP and PS images.

分別抽取傳統彈性高斯束偏移和最小二乘彈性高斯束偏移成像剖面其中的一道，繪制單道頻譜曲線，如圖9a、9b所示.由于震源為Ricker子波，傳統偏移成像的頻譜范圍有限，通過多次迭代，彈性高斯束最小二乘偏移PP和PS成像的頻譜范圍都得到一定的拓寬，進而說明迭代后的成像分辨率有所提高.

圖9 1次迭代和10次迭代后PP、PS成像的頻譜Fig.9 The spectrums of PP and PS images at the 1st and 10th iteration

2.3 3D彈性SEG/EAGE Salt模型

圖10 截斷SEG/EAGE Salt模型(a) P波速度； (b) PP反射率.Fig.10 The truncated SEG/EAGE Salt model(a) P wave velocity； (b) PP reflectivity.

圖11 不同方法正演的三分量地震記錄有限差分正演z分量(a)、x分量(b)、y分量(c)數據；彈性高斯束Born正演z分量(d)、x分量(e)、y分量(f)數據.Fig.11 The three-component data using different modeling methodsz-component (a)，x-component (b) and y-component (c) data using finite difference forward modeling method；z-component (d);x-component (e) and y-component (f) data using elastic Gaussian beam Born modeling method.

為了比較三維最小二乘彈性高斯束偏移PP、PS成像與傳統彈性高斯束偏移PP、PS成像，圖12給出三維偏移結果x=1.2 km、y=2 km、z=2.4 km上的成像切片.通過圖中黑色箭頭所指位置可以看出，相比傳統高斯束偏移，迭代后的PP和PS成像的分辨率都有了明顯地提高，深層照明也得到了補償，并且更清晰地刻畫出地層的特征.深度切片能夠幫助我們了解地下介質的分布特征，圖13展示了z=1.3 km處的深度切片.圖13a、13c表示傳統彈性高斯束偏移PP和PS的深度切片，圖13b、13d表示最小二乘彈性高斯束偏移PP和PS的深度切片，如圖所示，最小二乘彈性高斯束偏移深度切片的分辨率更高，照明范圍更廣.

圖12 (a)和(c)是傳統彈性高斯束PP和PS波偏移結果，(b)和(d)是最小二乘彈性高斯束PP和PS波偏移結果Fig.12 (a) and (c) are conventional elastic Gaussian beam migration results for PP and PS waves; (b) and (d) least-squares elastic Gaussian beam migration results for PP and PS waves

圖13 z=1.3 km切片上傳統彈性高斯束PP和PS偏移(a、c)，最小二乘彈性高斯束PP和PS偏移(b、d)的成像結果Fig.13 Depth slice at z=1.3 km.(a) and (c) are conventional elastic Gaussian beam migration results for PP and PS waves; (b) and (d) least-squares elastic Gaussian beam migration results for PP and PS waves

分別作出y=2 km處傳統彈性高斯束偏移和最小二乘彈性高斯束偏移PP、PS成像剖面的二維頻譜，如圖14所示.對比圖14a、14c，PS成像的分辨率比PP成像的分辨率要高，頻譜范圍更廣.分別比較圖14a、14b和圖14c、14d中方框所示部分，最小二乘彈性高斯束偏移方法一定程度上拓寬了PP和PS成像的頻譜范圍，因此可以證明，最小二乘彈性高斯束偏移方法可以提高成像分辨率.

圖14 1次迭代(a、c)和10次迭代(b、d)后PP(a、b)和PS(c，d)波成像的二維頻譜Fig.14 The 2D spectrums of images for PP (a，b) and PS (c，d) waves at the 1st (a，c) and 10th (b，d) iteration

對比圖15給出的1次迭代和10次迭代后三分量反偏移模擬數據與觀測數據殘差可知，地震記錄殘差中的絕大部分一次反射波經過迭代后收斂.我們在最小二乘彈性高斯束偏移之前對數據進行了去直達波處理，但仍存在部分殘余的直達波，這部分直達波殘差在迭代過程中未收斂.由于地層構造復雜，存在高斯束無法完全覆蓋的區(qū)域，因此反偏移也無法模擬該部分正演的地震波.

圖15 觀測和模擬單炮地震記錄的殘差1次迭代后z分量(a)、x分量(b)、y分量(c)的數據殘差；10次迭代后z分量(d)、x分量(e)、y分量(f)的數據殘差.Fig.15 The data residuals between simulated and observed data. z-component (a)，x-component (b) and y-component (c) residuals at the 1st iteration；z-component (d); x-component (e) and y-component (f) residuals at the 10th iteration

3 結論

本文推導了彈性高斯束Born正演(反偏移)模擬多分量地震記錄的公式，通過在數據域進行反偏移模擬數據與實際觀測數據的擬合實現最小二乘意義下的三維彈性高斯束偏移成像，并采用共軛梯度法迭代修正偏移成像振幅.本文提出在每次偏移和反偏移迭代過程中通過構建波型波矢量轉換矩陣完成多波型、多分量局部平面波之間的轉換，從而降低最小二乘偏移PP和PS成像串擾，并且提高反演效率.為了更好地平衡PP和PS成像反演的結果，我們在迭代過程中分別控制修正PP和PS成像的步長，使迭代過程穩(wěn)定收斂.模型測試結果表明，與傳統彈性高斯束偏移方法相比，最小二乘彈性高斯束偏移方法可以顯著地提高成像結果的分辨率，提高成像剖面的保幅能力，以及增強復雜構造的成像照明度.