聲波最小二乘偏移不同精確伴隨算子對的定量關系分析

2022-06-11 01:24:28王建森任玉曉楊傳根許新驥

石油地球物理勘探 2022年3期

關鍵詞：模型

王建森任玉曉陳磊嚴冬楊傳根許新驥

(①山東大學齊魯交通學院，山東濟南 250061； ②山東大學巖土與結構工程研究中心，山東濟南 250061； ③華能西藏水電安全工程技術研究中心，西藏林芝 860000)

0 引言

在反射波地震勘探中，偏移是利用地震記錄對地下結構進行成像的重要手段。近幾十年來，地球內部及近地表構造研究的需求推動了偏移技術的迅速發展。然而，當面對地震觀測數據不足或不規則、地下結構復雜和波場帶寬有限時，傳統偏移技術如Kirchhoff偏移和逆時偏移(RTM)很難將數據偏移到準確的位置，難以滿足實際地震勘探需求。

隨著計算能力的飛速發展，提出了基于最小二乘優化框架的地震偏移方法，在真實速度擾動模型成像方面取得了良好的應用效果。在最小二乘偏移(LSM)中，地震偏移成像結果通常通過最小化數據殘差的L2范數獲得，該殘差表征了正演模擬數據與實際觀測數據之間的差異[1]。LSM的思想首先應用于Kirchhoff偏移[2-3]，然后推廣到單程波動方程偏移[4-5]。目前，該思想已應用于雙程波動方程，形成了最小二乘逆時偏移(LSRTM)[6-11]。與傳統的RTM 成像[12-14]相比，LSRTM有助于生成高分辨率、高保真度和較少低頻噪聲的偏移成像結果[15-17]。

一般來說，常規的偏移方法利用地震正演算子的共軛轉置而非真正意義的逆算子成像，導致由成像結果重建的地震數據與原始地震數據之間存在誤差。LSM通過基于線性反演迭代的方法逼近逆算子，最終獲得同預處理地震數據相匹配的地下構造模型，這種線性反演方法需要通過正演算子(或反偏移算子)和偏移算子建立地下構造模型與地震觀測數據之間的線性關系[18-19]，正演—偏移算子對的伴隨性是影響迭代收斂和最終偏移成像結果的重要因素。通常，正演算子通過Born近似理論導出[20-21]，偏移算子最常用的是RTM。然而，Xu等[22]證明了使用一階聲波方程時該算子對(Born正演和RTM)經過數值離散計算后，很難滿足精確伴隨關系[18]。因此，構建滿足精確伴隨關系的正演—偏移算子對是最小二乘逆時偏移中的一個關鍵步驟[23-25]，也是提高偏移質量的有效手段。

針對伴隨算子對存在的問題，國內外學者主要采用兩種途徑進行研究：一是基于伴隨狀態法推導連續伴隨波動方程，再謹慎進行離散獲得波場延拓伴隨算子，保證正演—偏移算子對的精確伴隨性[26-27]；二是基于矩陣形式先進行離散后直接推導伴隨矩陣的策略。通常采用點積測試[19]作為一種便捷的數值方法檢驗兩個算子之間的伴隨性，而點積測試最簡單、直接的方法之一是利用算子矩陣的轉置。不同于傳統的疊后RTM算子，Ji[28]提出了一種適用于時域聲波方程有限差分模擬的矩陣計算方法并推導了伴隨矩陣； Yao等[29]基于矩陣形式對頻率域LSRTM精確伴隨算子對進行了研究； Xu等[22]詳細描述了基于聲波Born近似理論的精確伴隨算子對的推導過程，并應用于預條件LSM。許多學者將這種基于Born正演的精確伴隨算子對應用于聲波[22，30-31]、彈性波[15，32]，甚至各向異性[33-34]、黏滯介質[7，9，35]的LSM。由于RTM是最常用的偏移算子，因此，另一種構造精確伴隨算子對的方法是先通過RTM過程的矩陣表示，再將其轉置，得到RTM反偏移算子(DeRTM算子)，DeRTM-RTM算子對同樣有較好的成像效果[36]。此外，一些學者利用自伴隨波動方程將這兩種算子對混在一起，并展示了Born-RTM算子對在各種數值算例下的良好應用效果[30，37-39]。然而，這三組伴隨算子對之間的區別與聯系尚不清楚，需要進一步研究和探索。

為此，本文首先在Ji[28]的研究基礎上，將時域二階聲波方程離散成矩陣形式，得到Born正演和RTM過程的矩陣表示。將矩陣轉置即產生兩組精確伴隨算子對，分別表示為Born-AdjBorn算子對和RTM-DeRTM算子對。然后，基于聲波方程的自伴隨離散化，對這兩個算子進行混合匹配，得到了第三組精確伴隨算子對，即自伴隨Born-RTM算子對。最后進一步推導了這三組精確伴隨算子對之間的定量關系，并通過二維數值算例進行了驗證。

1 方法原理

1.1 聲波波場延拓及其伴隨算子矩陣形式

1.1.1 波場延拓算子

一維時空域的二階聲波方程為

(1)

式中：p為壓力場；c(z)為聲波速度；t和z分別表示時間和深度。根據中心二階有限差分格式，式(1)的離散形式可以寫為

(2)

式中：n=1，2，…，nt，nt為時間方向樣點數；i=1，2，…，nz，nz為深度方向樣點數；α=Δt2/Δz2，其中Δt和Δz分別為時間間隔和深度網格間隔。將式(2)表示為矩陣形式

(3)

式中：pn=[p1,n,p2,n,…,pnz,n]T為n時刻波場的向量形式；I是單位矩陣；T是一個帶狀矩陣，具體為

(4)

Bp=s

(5)

(6)

地震記錄可以通過觀測矩陣D提取各時刻對應的波場。因此，觀測數據可以表示為

d=Dp=DB-1s

(7)

二維波場延拓算子與一維類似，詳細推導見附錄A。

1.1.2 波場延拓伴隨算子

如上所述，推導一對伴隨算子最直接的方法之一是通過矩陣轉置。因此，了解波場延拓算子B-1及其伴隨算子B-T對于后續推導正演—偏移算子對有較大幫助。為了簡單起見，只給出一維矩陣分析，對于高維可得出類似的結論。

應用線性代數方法，B-1可以顯式表示為

B-1=I+A+A2+…+Ant

(8)

定義一系列矩陣算子{Pn}，n∈{1,2,…,nt}，將零時刻波場映射到n時刻波場。考慮式(6)中的定義，令P0=I，可得

(9)

因此，波場延拓p=B-1s可以寫成子矩陣求和的形式

(10)

上式可以理解為，在n時刻的波場pn是由此前源項si通過Pn-i映射的所有波場的總和，這符合惠更斯原理和波場疊加原理。其伴隨過程p=B-Ts，也可以表示為子陣求和形式

(11)

(12)

通過引入速度加權波場p(z,t)=c2(z)q(z,t)，可以將式(12)轉化為常規波動方程

(13)

式(13)表明伴隨波場q可以通過在常規波動方程中加載速度加權源項c2s獲得。在速度模型不均勻的情況下，忽略速度加權因子將導致伴隨波場的計算不準確，但這與二階聲波方程自伴隨的性質不一致。為了保持空間離散中的自伴隨性，引入另一個中間波場u(z,t)=c-1(z)p(z,t)，滿足

(14)

1.2 三組精確伴隨算子對

1.2.1 基于Born近似的精確伴隨算子對

根據Born近似理論，在小速度擾動假設下，Born正演方程為

(15)

(16)

式中⊙表示Hadamard乘積。因此Born正演獲得的觀測數據可以表示為

(17)

Born正演的伴隨算子是式(17)中正演算子FBorn的轉置，即

(18)

1.2.2 基于RTM的精確伴隨算子對

經典的RTM過程包括正向波場延拓、逆向波場延拓和成像。利用矩陣形式表示RTM全過程是推導基于RTM精確伴隨算子對的關鍵。正向和逆向波場延拓可表示為

(19)

(20)

(21)

通過零延遲互相關成像條件，偏移過程可表示為

(22)

式中FRTM為逆時偏移算子。相應的反偏移過程可表示為

(23)

1.2.3 基于混合匹配的精確伴隨算子對

為了保證混合匹配后算子對的伴隨性，簡單直接的方法是使用式(14)，而不是傳統的聲波方程，可以保證FBorn-s=FDeRTM-s和FAdjBorn-s=FRTM-s，即利用自伴隨聲波方程得到的兩組伴隨算子對Born-RTM和DeRTM-AdjBorn實際上是相同的。因此，在一般情況下，可以使用自伴隨Born-RTM算子對作為LSM中的另一組精確伴隨算子對。

綜上所述，得到了LSM的三組精確伴隨算子對，分別是利用常規離散波動方程的Born-AdjBorn算子對、利用常規離散波動方程的DeRTM-RTM算子對和利用自伴隨離散波動方程的Born-RTM算子對。對應的LSM過程分別稱為LSBM、LSRTM和自伴隨LSBRTM。此外，考慮自伴隨波場u、常規離散化波場p及其伴隨波場q之間的關系，即u=c-1p=cq，可以得到三組精確伴隨算子對成像結果之間的定量關系。

1.3 三組精確伴隨算子對的定量關系

因此，在大多數采用常規離散波動方程的情況下，當且僅當同一源在均勻速度模型中延拓時，Born-AdjBorn算子對和DeRTM-RTM算子對是相同的。在此條件下，Born近似正演和RTM可以作為一組精確的伴隨算子對。

其次，根據附錄C中的推導，可以得出以下定量結論。

(1)FBorn=FBorn-s和FAdjBorn=FRTM-s。常規離散波動方程的Born-AdjBorn算子對與自伴隨離散波動方程的Born-RTM算子對相同。因此，AdjBorn結果等于自伴隨RTM結果，LSBM和自伴隨LSBRTM的結果是相同的。

由于式(C-6)中的Hessian矩陣通常是高度不適定的，其逆矩陣很難通過數值計算。因此，在式(C-6)中的最小二乘解通常采用迭代優化算法計算，例如本文中使用的預條件共軛梯度法。然而，由于迭代優化算法難以收斂至極小值，導致數值解可能與最小二乘解并不完全一致，式(C-7)中LSM結果之間的定量關系較難滿足。

1.4 基于精確伴隨算子對的LSM算法

LSM的目標是速度擾動模型，該模型描述了觀測數據所隱含的地下構造信息。因此求解過程通常包括正演算子、偏移算子與最小二乘算法。由于最小二乘算法在大多數LSRTM文獻中都有講述[27，40-41]，因此本文將重點關注精確伴隨算子對的正確實現。

通過上述推導，得到了LSBM、LSRTM和自伴隨LSBRTM伴隨算子對的矩陣形式。然而，它們的矩陣運算是非常昂貴的，矩陣維度為nx×nz×nt。因此，需要將矩陣運算轉化可實現的程序語言。表1總結了各算子、矩陣的實際意義。

LSM算法的目標是通過最小化目標函數求解速度擾動模型m

(24)

式中F表示前面推導出的三個正演算子之一。目標函數的第一項表示數據殘差，第二項是由阻尼參數λ平衡的模型正則化項，并采用了預條件共軛梯度算法獲得LSBM、LSRTM和自伴隨LSBRTM成像結果。

表1 伴隨算子對中各矩陣或算子的實際意義

2 數值算例

在數值算例中，應用Ricker子波、基于這些算子對的有限差分代碼計算合成數據。由于所推導的算子對不考慮任何吸收邊界條件，本文通過復制模型的外邊界層擴展模型，并在數值邊界反射波到達之前停止波場延拓。在預條件共軛梯度算法中，阻尼參數λ設置為0.001。

2.1 層狀模型

自伴隨LSBRTM、LSBM、LSRTM的一次迭代偏移結果分別如圖2a～圖2c所示，分別等價于自伴隨RTM、傳統Born偏移、傳統RTM結果。LSM的30次迭代偏移結果如圖3所示。

與一次迭代偏移結果相比，多次迭代能夠以更高的分辨率、更均衡的振幅和更少的偽影刻畫出貼近真實速度擾動模型的成像結果。在波速分布由低速變化至高速的界面位置處，其速度擾動模型由負值變為正值，在一次迭代偏移結果中界面刻畫均為正值—負值—正值，在LSM結果中界面均收斂為接近真實速度擾動模型的負值—正值。

圖1 層狀模型(a)原始速度； (b)平滑速度； (c)速度擾動

圖2 層狀模型不同方法的一次迭代偏移結果(a)自伴隨RTM； (b)AdjBorn； (c)RTM

圖3 層狀模型四種LSM方法的30次迭代結果(a)自伴隨LSBRTM； (b)LSBM；(c)LSRTM； (d)速度加權LSRTM

為了驗證三組伴隨算子對結果之間的定量關系，計算了不同方法偏移成像結果的差值，如圖4所示。圖4a、圖4b和圖4d中極小的振幅證明了傳統Born偏移結果(AdjBorn)與自伴隨RTM的等價性、速度加權自伴隨RTM與傳統RTM之間的等價性以及LSBM與自伴隨LSBRTM的等價性。對于自伴隨LSBRTM與速度加權LSRTM結果的差異，由于存在誤差，無法收斂至零，圖4e難以證明其對應的定量關系。

此外，圖4c中的差值剖面是一個沒有明顯偽影的層狀模型圖像。然而，在圖2a和圖2c中，自伴隨RTM和傳統RTM結果在地表附近都存在較強的低頻偽影，其與速度變化和反射結構無關，而成像結果的差異只與速度變化有關，即目標反射體。因此，這可能是一個有效減少RTM成像結果中低頻偽影的策略。

圖4 層狀模型不同方法偏移結果的差值剖面(a)傳統Born偏移與自伴隨RTM； (b)速度加權自伴隨RTM與傳統RTM； (c)傳統RTM與自伴隨RTM；(d)自伴隨LSBRTM與LSBM； (e)速度加權LSRTM與自伴隨LSBRTM

在均勻偏移速度模型(2500m/s)下測試了LSBM與LSRTM，偏移成像結果如圖5所示。由于存在速度誤差，第二個界面位置存在偏差，但兩個偏移成像結果之間的差值為零，驗證了本文的推論，即當偏移速度模型是均勻的時候，LSBM和LSRTM會呈現相同的成像結果。

在速度存在誤差、數據含有噪聲、數據有缺失的情況下，對三組精確伴隨算子對進行敏感性分析。

應用平滑半徑為100個網格距的偏移速度模型(圖6a)，LSBM、LSRTM和自伴隨LSBRTM成像結果如圖7所示，可見，在速度誤差的影響下，LSBM、LSRTM和自伴隨LSBRTM均可對地下層狀構造較準確成像。

在原始地震數據中添加白噪聲，數據信噪比為-20dB(圖6b)。三組精確伴隨算子對的LSM成像結果如圖8所示，受噪聲的影響，均產生了較多的偽影。

圖5 均勻偏移速度模型下兩種方法的LSM結果(a)LSBM； (b)LSRTM

圖6 低精度的偏移速度速度模型(a)、震源位于1000m的含噪聲數據(b)及隨機缺失87.5%的道集數據(c)

每炮地震數據隨機缺失175道，即缺失87.5%(圖6c)，三組精確伴隨算子對的LSM成像結果如圖9所示。數據缺失對LSBM、LSRTM和自伴隨LSBRTM偏移結果影響一致，均會在淺層產生些許低頻噪聲。

在速度存在誤差、數據含有噪聲、數據有缺失的情況下，LSBM與自伴隨LSBRTM之間差值極小，在10-4數量級； LSBM與LSRTM、自伴隨LSBRTM與LSRTM之間的差值較大，與成像結果的幅值在同一數量級。因此，速度誤差、數據噪聲以及數據缺失對三組精確伴隨算子對之間的定量關系基本沒有影響。

2.2 Marmousi模型

應用更復雜的Marmousi模型(圖10)測試三種的精確伴隨算子對。模型網格數為1000×300，網格間距為10m。50個震源均勻分布在模型表面，震源主頻為10Hz。檢波器以10m間距均勻分布在模型表面。一次迭代偏移結果如圖11所示，經過10次迭代，最終的LSM結果如圖12所示。

圖7 低精度速度模型的三種方法的LSM結果(a)LSBM； (b)LSRTM； (c)自伴隨LSBRTM

圖8 數據含有噪聲時三種方法的LSM結果(a)LSBM； (b)LSRTM； (c)自伴隨LSBRTM

圖9 數據有缺失時三種方法的LSM結果(a)LSBM； (b)LSRTM； (c)自伴隨LSBRTM

與層狀模型結果類似，圖11和圖12中所有偏移結果符合預期，LSM結果展現了更少的低頻噪聲和更清晰的界面信息。圖12中紅色矩形和箭頭突出了成像結果之間的差異，通過與速度擾動模型進行對比可以看出，LSBM與自伴隨LSBRTM更關注淺層的精細成像； LSRTM相較于LSBM和自伴隨LSBRTM對深層(2km以下)成像效果更好，紅

圖10 Marmousi模型(a)真實速度； (b)平滑速度； (c)速度擾動

圖11 Marmousi模型不同方法的一次迭代偏移結果(a)自伴隨RTM； (b)AdjBorn； (c)RTM

色矩形中的地質構造更加清晰；并且速度加權項LSRTM對深層的照明更強，紅色箭頭處的地層信息與速度擾動模型基本一致。說明在Marmousi模型算例中，顯式加入速度加權項后會增強深層的成像效果。

在定量關系上，圖13中極小的振幅分別驗證了一次迭代偏移結果和LSM結果相應的定量關系，與層狀模型算例類似，由于存在迭代誤差，難以證明速度加權LSRTM與LSBM(或自伴隨LSBRTM)之間的定量關系。圖14的數據收斂曲線表明三組伴隨算子對的收斂速度快且相似，經過10次迭代后，三種LSM方法的歸一化數據殘差都小于0.001。

圖12 Marmousi模型四種LSM方法的10次迭代結果(a)自伴隨LSBRTM； (b)LSBM； (c)LSRTM； (d)速度加權LSRTM

圖13 Marmousi模型不同方法偏移結果的差值剖面(a)傳統Born偏移和自伴隨RTM； (b)速度加權自伴隨RTM和傳統RTM； (c)自伴隨LSBRTM與LSBM

圖14 Mairmousi模型三種LSM方法的數據殘差對比(a)第15號炮； (b)第35號炮； (c)所有50炮

3 討論

一般認為Born正演和RTM分別是LSM的正演和偏移算子。然而，在數值離散化后，它們可能無法準確地通過點積測試。當采用非精確伴隨算子對進行LSM時，會導致收斂性差或累積誤差較大的偏移成像結果。因此，在聲波方程矩陣表達式和矩陣運算的基礎上，本文從理論上推導了基于Born近似理論和RTM過程的三組精確伴隨算子對，并找出了它們之間的等價條件和定量關系。同時，開發了一套無矩陣代碼，通過數值算例驗證該推論。

正如在精確伴隨算子對推導中所證明的，LSBM和LSRTM成像結果之間的主要區別源于相應算子對中的波場延拓矩陣B-1和B-T，即式(1)及其伴隨波動方程(式(14))。后者可以進一步擴展為

(25)

式(25)意味著算子B-T實際上是傳統的聲波方程B-1加上與速度空間擾動有關的兩項。如果偏移速度模型為均勻速度模型，即?c2/?z=0、?2c2/?z2=0，則LSBM和LSRTM會得到相同的結果。因此，可以從不同的角度對Born-AdjBorn和DeRTM-RTM伴隨算子對進行對比，并得到相同的結論。并且，在二階聲波方程中(以一維情況為例)，若將速度項移至方程左側，即變為慢度項，即

(26)

將該聲波方程進行離散，并加載震源項s，其結果與加載速度加權源項c2s的式(13)相似。這意味著，在LSM過程中，將慢度擾動模型作為收斂目標即可利用該聲波方程的自伴隨性質進行波場延拓模擬。因此，在最小二乘框架下的LSM方法中，波動方程的選擇及模型參數化對于迭代收斂性能和最終的偏移成像結果影響較大。

此外，不論對于連續波動方程還是波動方程的離散形式，由于Born-AdjBorn和DeRTM-RTM算子對都是與一次反射波對應的線性化算子，實際觀測數據中可能包含更復雜的波場分量，如多次波或繞射波等。因此，通過正演模擬(Born正演和DeRTM過程)很難得到與實際觀測數據準確匹配的合成數據，所以不應追求獲得非常小的數據殘差，特別是對于復雜的Marmousi模型。這使得驗證LSRTM與LSBM(或自伴隨LSBRTM)結果之間的定量關系(式(C-7))變得更加困難。

在實際地震勘探過程中，地震觀測數據存在較嚴重的噪聲，且實際子波與波場延拓模擬子波之間存在差異，正演地震數據與實際觀測數據本身具有一定誤差。若采用非精確伴隨算子對進行最小二乘偏移，會引入新的誤差，且在迭代收斂過程中容易累計誤差，難以保證迭代收斂，導致實際資料偏移結果存在偽影、假象。因此，在偏移速度模型較為準確的前提下，采用精確伴隨算子對進行最小二乘偏移成像可有效提升實際地震資料的成像質量。

4 結論

最小二乘偏移需要正演和偏移兩個過程實現成像結果的迭代更新，并且正演和偏移過程需滿足精確伴隨關系。基于Born近似理論和RTM過程，本文首先在常規離散化聲波方程的基礎上推導出了兩對理論伴隨算子。采用自伴隨方法進行波動方程離散時，這兩組算子對會產生另一組精確伴隨算子對。這三組算子對均可通過無矩陣算法編程實現，并能準確通過點積測試。此外，還推導了這三組伴隨算子對之間的定量關系，并應用二維模型進行了驗證。定量關系為：采用同一源在均勻速度模型中進行波場延拓時，Born-AdjBorn算子對和DeRTM-RTM算子對是相同的；速度加權的LSRTM結果等于LSBM和自伴隨LSBRTM的結果。

附錄A 二維聲波波場延拓算子

對于時空域二維聲波方程

(A-1)

假設Δz=Δx，其有限差分格式為

pi+1,j,n+pi,j-1,n+pi,j+1,n)-pi,j,n-1

(A-2)

將波場表示為矩陣形式

pn=[p1,1,n,p2,1,n,…,pnz,1,n,p1,2,n,p2,2,n,…,

pnz,2,n,…,p1,nx,n,p2,nx,n,…,pnz,nx,n]T

(A-3)

式中nx、nz分別為x和z方向的樣點數。在二維情況下，矩陣T變成了大小為(nx×nz)×(nx×nz)的三對角塊矩陣

(A-4)