聲波最小二乘偏移不同精確伴隨算子對(duì)的定量關(guān)系分析

王建森 任玉曉 陳 磊 嚴(yán) 冬 楊傳根 許新驥

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

0 引言

在反射波地震勘探中，偏移是利用地震記錄對(duì)地下結(jié)構(gòu)進(jìn)行成像的重要手段。近幾十年來(lái)，地球內(nèi)部及近地表構(gòu)造研究的需求推動(dòng)了偏移技術(shù)的迅速發(fā)展。然而，當(dāng)面對(duì)地震觀(guān)測(cè)數(shù)據(jù)不足或不規(guī)則、地下結(jié)構(gòu)復(fù)雜和波場(chǎng)帶寬有限時(shí)，傳統(tǒng)偏移技術(shù)如Kirchhoff偏移和逆時(shí)偏移(RTM)很難將數(shù)據(jù)偏移到準(zhǔn)確的位置，難以滿(mǎn)足實(shí)際地震勘探需求。

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

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

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

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

1 方法原理

1.1 聲波波場(chǎng)延拓及其伴隨算子矩陣形式

1.1.1 波場(chǎng)延拓算子

一維時(shí)空域的二階聲波方程為

(1)

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

(2)

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

(3)

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

(4)

Bp=s

(5)

(6)

地震記錄可以通過(guò)觀(guān)測(cè)矩陣D提取各時(shí)刻對(duì)應(yīng)的波場(chǎng)。因此，觀(guān)測(cè)數(shù)據(jù)可以表示為

d=Dp=DB-1s

(7)

二維波場(chǎng)延拓算子與一維類(lèi)似，詳細(xì)推導(dǎo)見(jiàn)附錄A。

1.1.2 波場(chǎng)延拓伴隨算子

如上所述，推導(dǎo)一對(duì)伴隨算子最直接的方法之一是通過(guò)矩陣轉(zhuǎn)置。因此，了解波場(chǎng)延拓算子B-1及其伴隨算子B-T對(duì)于后續(xù)推導(dǎo)正演—偏移算子對(duì)有較大幫助。為了簡(jiǎn)單起見(jiàn)，只給出一維矩陣分析，對(duì)于高維可得出類(lèi)似的結(jié)論。

應(yīng)用線(xiàn)性代數(shù)方法，B-1可以顯式表示為

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

(8)

定義一系列矩陣算子{Pn}，n∈{1,2,…,nt}，將零時(shí)刻波場(chǎng)映射到n時(shí)刻波場(chǎng)?？紤]式(6)中的定義，令P0=I，可得

(9)

因此，波場(chǎng)延拓p=B-1s可以寫(xiě)成子矩陣求和的形式

(10)

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

(11)

(12)

通過(guò)引入速度加權(quán)波場(chǎng)p(z,t)=c2(z)q(z,t)，可以將式(12)轉(zhuǎn)化為常規(guī)波動(dòng)方程

(13)

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

(14)

1.2 三組精確伴隨算子對(duì)

1.2.1 基于Born近似的精確伴隨算子對(duì)

根據(jù)Born近似理論，在小速度擾動(dòng)假設(shè)下，Born正演方程為

(15)

(16)

式中⊙表示Hadamard乘積。因此Born正演獲得的觀(guān)測(cè)數(shù)據(jù)可以表示為

(17)

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

(18)

1.2.2 基于RTM的精確伴隨算子對(duì)

經(jīng)典的RTM過(guò)程包括正向波場(chǎng)延拓、逆向波場(chǎng)延拓和成像。利用矩陣形式表示RTM全過(guò)程是推導(dǎo)基于RTM精確伴隨算子對(duì)的關(guān)鍵。正向和逆向波場(chǎng)延拓可表示為

(19)

(20)

(21)

通過(guò)零延遲互相關(guān)成像條件，偏移過(guò)程可表示為

(22)

式中FRTM為逆時(shí)偏移算子。相應(yīng)的反偏移過(guò)程可表示為

(23)

1.2.3 基于混合匹配的精確伴隨算子對(duì)

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

綜上所述，得到了LSM的三組精確伴隨算子對(duì)，分別是利用常規(guī)離散波動(dòng)方程的Born-AdjBorn算子對(duì)、利用常規(guī)離散波動(dòng)方程的DeRTM-RTM算子對(duì)和利用自伴隨離散波動(dòng)方程的Born-RTM算子對(duì)。對(duì)應(yīng)的LSM過(guò)程分別稱(chēng)為L(zhǎng)SBM、LSRTM和自伴隨LSBRTM。此外，考慮自伴隨波場(chǎng)u、常規(guī)離散化波場(chǎng)p及其伴隨波場(chǎng)q之間的關(guān)系，即u=c-1p=cq，可以得到三組精確伴隨算子對(duì)成像結(jié)果之間的定量關(guān)系。

1.3 三組精確伴隨算子對(duì)的定量關(guān)系

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

其次，根據(jù)附錄C中的推導(dǎo)，可以得出以下定量結(jié)論。

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

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

1.4 基于精確伴隨算子對(duì)的LSM算法

LSM的目標(biāo)是速度擾動(dòng)模型，該模型描述了觀(guān)測(cè)數(shù)據(jù)所隱含的地下構(gòu)造信息。因此求解過(guò)程通常包括正演算子、偏移算子與最小二乘算法。由于最小二乘算法在大多數(shù)LSRTM文獻(xiàn)中都有講述[27，40-41]，因此本文將重點(diǎn)關(guān)注精確伴隨算子對(duì)的正確實(shí)現(xiàn)。

通過(guò)上述推導(dǎo)，得到了LSBM、LSRTM和自伴隨LSBRTM伴隨算子對(duì)的矩陣形式。然而，它們的矩陣運(yùn)算是非常昂貴的，矩陣維度為nx×nz×nt。因此，需要將矩陣運(yùn)算轉(zhuǎn)化可實(shí)現(xiàn)的程序語(yǔ)言。表1總結(jié)了各算子、矩陣的實(shí)際意義。

LSM算法的目標(biāo)是通過(guò)最小化目標(biāo)函數(shù)求解速度擾動(dòng)模型m

(24)

式中F表示前面推導(dǎo)出的三個(gè)正演算子之一。目標(biāo)函數(shù)的第一項(xiàng)表示數(shù)據(jù)殘差，第二項(xiàng)是由阻尼參數(shù)λ平衡的模型正則化項(xiàng)，并采用了預(yù)條件共軛梯度算法獲得LSBM、LSRTM和自伴隨LSBRTM成像結(jié)果。

表1 伴隨算子對(duì)中各矩陣或算子的實(shí)際意義

2 數(shù)值算例

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

2.1 層狀模型

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

與一次迭代偏移結(jié)果相比，多次迭代能夠以更高的分辨率、更均衡的振幅和更少的偽影刻畫(huà)出貼近真實(shí)速度擾動(dòng)模型的成像結(jié)果。在波速分布由低速變化至高速的界面位置處，其速度擾動(dòng)模型由負(fù)值變?yōu)檎?，在一次迭代偏移結(jié)果中界面刻畫(huà)均為正值—負(fù)值—正值，在LSM結(jié)果中界面均收斂為接近真實(shí)速度擾動(dòng)模型的負(fù)值—正值。

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

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

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

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

此外，圖4c中的差值剖面是一個(gè)沒(méi)有明顯偽影的層狀模型圖像。然而，在圖2a和圖2c中，自伴隨RTM和傳統(tǒng)RTM結(jié)果在地表附近都存在較強(qiáng)的低頻偽影，其與速度變化和反射結(jié)構(gòu)無(wú)關(guān)，而成像結(jié)果的差異只與速度變化有關(guān)，即目標(biāo)反射體。因此，這可能是一個(gè)有效減少RTM成像結(jié)果中低頻偽影的策略。

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

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

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

應(yīng)用平滑半徑為100個(gè)網(wǎng)格距的偏移速度模型(圖6a)，LSBM、LSRTM和自伴隨LSBRTM成像結(jié)果如圖7所示，可見(jiàn)，在速度誤差的影響下，LSBM、LSRTM和自伴隨LSBRTM均可對(duì)地下層狀構(gòu)造較準(zhǔn)確成像。

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

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

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

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

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

2.2 Marmousi模型

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

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

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

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

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

圖10 Marmousi模型(a)真實(shí)速度； (b)平滑速度； (c)速度擾動(dòng)

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

色矩形中的地質(zhì)構(gòu)造更加清晰； 并且速度加權(quán)項(xiàng)LSRTM對(duì)深層的照明更強(qiáng)，紅色箭頭處的地層信息與速度擾動(dòng)模型基本一致。說(shuō)明在Marmousi模型算例中，顯式加入速度加權(quán)項(xiàng)后會(huì)增強(qiáng)深層的成像效果。

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

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

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

圖14 Mairmousi模型三種LSM方法的數(shù)據(jù)殘差對(duì)比(a)第15號(hào)炮； (b)第35號(hào)炮； (c)所有50炮

3 討論

一般認(rèn)為Born正演和RTM分別是LSM的正演和偏移算子。然而，在數(shù)值離散化后，它們可能無(wú)法準(zhǔn)確地通過(guò)點(diǎn)積測(cè)試。當(dāng)采用非精確伴隨算子對(duì)進(jìn)行LSM時(shí)，會(huì)導(dǎo)致收斂性差或累積誤差較大的偏移成像結(jié)果。因此，在聲波方程矩陣表達(dá)式和矩陣運(yùn)算的基礎(chǔ)上，本文從理論上推導(dǎo)了基于Born近似理論和RTM過(guò)程的三組精確伴隨算子對(duì)，并找出了它們之間的等價(jià)條件和定量關(guān)系。同時(shí)，開(kāi)發(fā)了一套無(wú)矩陣代碼，通過(guò)數(shù)值算例驗(yàn)證該推論。

正如在精確伴隨算子對(duì)推導(dǎo)中所證明的，LSBM和LSRTM成像結(jié)果之間的主要區(qū)別源于相應(yīng)算子對(duì)中的波場(chǎng)延拓矩陣B-1和B-T，即式(1)及其伴隨波動(dòng)方程(式(14))。后者可以進(jìn)一步擴(kuò)展為

(25)

式(25)意味著算子B-T實(shí)際上是傳統(tǒng)的聲波方程B-1加上與速度空間擾動(dòng)有關(guān)的兩項(xiàng)。如果偏移速度模型為均勻速度模型，即?c2/?z=0、?2c2/?z2=0，則LSBM和LSRTM會(huì)得到相同的結(jié)果。因此，可以從不同的角度對(duì)Born-AdjBorn和DeRTM-RTM伴隨算子對(duì)進(jìn)行對(duì)比，并得到相同的結(jié)論。并且，在二階聲波方程中(以一維情況為例)，若將速度項(xiàng)移至方程左側(cè)，即變?yōu)槁软?xiàng)，即

(26)

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

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

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

4 結(jié)論

最小二乘偏移需要正演和偏移兩個(gè)過(guò)程實(shí)現(xiàn)成像結(jié)果的迭代更新，并且正演和偏移過(guò)程需滿(mǎn)足精確伴隨關(guān)系?；贐orn近似理論和RTM過(guò)程，本文首先在常規(guī)離散化聲波方程的基礎(chǔ)上推導(dǎo)出了兩對(duì)理論伴隨算子。采用自伴隨方法進(jìn)行波動(dòng)方程離散時(shí)，這兩組算子對(duì)會(huì)產(chǎn)生另一組精確伴隨算子對(duì)。這三組算子對(duì)均可通過(guò)無(wú)矩陣算法編程實(shí)現(xiàn)，并能準(zhǔn)確通過(guò)點(diǎn)積測(cè)試。此外，還推導(dǎo)了這三組伴隨算子對(duì)之間的定量關(guān)系，并應(yīng)用二維模型進(jìn)行了驗(yàn)證。定量關(guān)系為：采用同一源在均勻速度模型中進(jìn)行波場(chǎng)延拓時(shí)，Born-AdjBorn算子對(duì)和DeRTM-RTM算子對(duì)是相同的； 速度加權(quán)的LSRTM結(jié)果等于LSBM和自伴隨LSBRTM的結(jié)果。

附錄A 二維聲波波場(chǎng)延拓算子

對(duì)于時(shí)空域二維聲波方程

(A-1)

假設(shè)Δz=Δx，其有限差分格式為

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

(A-2)

將波場(chǎng)表示為矩陣形式

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方向的樣點(diǎn)數(shù)。在二維情況下，矩陣T變成了大小為(nx×nz)×(nx×nz)的三對(duì)角塊矩陣

(A-4)

(A-5)

附錄B AdjBorn算子和RTM算子的 等價(jià)條件

將偏移過(guò)程分為以下三個(gè)步驟分析它們的等價(jià)條件。

(B-1)

(B-2)

(B-3)

(B-4)

通過(guò)上述分析可知，當(dāng)且僅當(dāng)采用同一震源子波，偏移速度模型均勻的情況下，AdjBorn算子和RTM算子等價(jià)，即

FAdjBorn=FRTM

(B-5)

(B-6)

附錄C 三組伴隨算子對(duì)之間的關(guān)系

(C-1)

(C-2)

=FAdjBorn-sd=FRTM-sd

(C-3)

(C-4)

(C-5)

FRTM-sd，可知FRTM=CFRTM-s。因此，C可以看作是一個(gè)正則化因子，并給出相應(yīng)的最小二乘解

=C-1mLSRTM-s

(C-6)

因此，三種LSM結(jié)果之間的定量關(guān)系為

(C-7)