楊育臣，方金偉，王寧，李松齡

1 Geology and Geophysics Program，Missouri University of Science and Technology，Rolla，MO 65409，USA 2 深部巖土力學(xué)與地下工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，中國(guó)礦業(yè)大學(xué)力學(xué)與土木工程學(xué)院，徐州 221116 3 東北石油大學(xué)地球科學(xué)學(xué)院，大慶 163318

0 引言

近年來(lái)，發(fā)展了寬頻帶、寬方位、高密度的地震數(shù)據(jù)采集技術(shù)，通過(guò)地震成像技術(shù)實(shí)現(xiàn)地下參數(shù)的高信噪比、高分辨率和高保真度的建模.其中，具有最高成像精度的反演成像以全波形反演(FWI)為代表，實(shí)現(xiàn)構(gòu)造信息和巖性油氣藏的高精度參數(shù)建模.全波形反演通過(guò)建立觀測(cè)和合成數(shù)據(jù)的某種距離度量的目標(biāo)泛函，利用優(yōu)化理論計(jì)算出參數(shù)下降方向和更新步長(zhǎng)，不斷修改模型直至合成數(shù)據(jù)接近觀測(cè)數(shù)據(jù)，從而獲得準(zhǔn)確的模型參數(shù).全波形反演能夠充分利用地震數(shù)據(jù)中的振幅和相位信息，獲得高分辨率高精度的速度模型.

計(jì)算效率低下是制約全波形反演進(jìn)行大規(guī)模計(jì)算以及實(shí)際資料處理的主要因素，因此提升全波形反演的效率就顯得十分重要(邢貞貞等，2019).除了使用高性能計(jì)算平臺(tái)加速之外(Yang et al.，2015；Gokhberg and Fichtner，2016；張猛等，2014；桂生等，2017)，通常使用算法類的效率優(yōu)化方案來(lái)加速反演.基于算法類的加速方法通常使用震源編碼策略提升全波形反演效.震源編碼的思想是通過(guò)某種震源組合方式將很多單炮形成一個(gè)超級(jí)震源(編碼炮)，實(shí)現(xiàn)幾十炮甚至上百炮的同時(shí)模擬，從而減少全波形反演中波場(chǎng)模擬的次數(shù)，極大地減少全波形反演的計(jì)算量.然而，相比于傳統(tǒng)的同時(shí)源模擬(震源編碼)方法，比如隨機(jī)源編碼(Krebs et al.，2009)、隨機(jī)時(shí)移(Dai and Schuster，2009)、平面波編碼(Zhang et al.，2003；Vigh and Starr，2008)等算法，一種基于頻率選擇的同時(shí)源方法不僅能夠保證計(jì)算效率，而且可以獲得不受編碼串?dāng)_噪聲影響的反演結(jié)果.Huang和Schuster(2012)第一次在頻率域全波形反演處理海洋觀測(cè)系統(tǒng)數(shù)據(jù)時(shí)使用了該方法.而時(shí)間-頻率域(混合域)的同時(shí)源反演方法最早由Dai等(2013)提出并用來(lái)加速最小二乘逆時(shí)偏移成像.在混合域同時(shí)源模擬方法中，將一系列不同頻率的單頻源組合起來(lái)形成了一個(gè)編碼源，并在時(shí)間域完成波場(chǎng)的模擬.編碼炮的梯度計(jì)算則是在頻率域完成的，通過(guò)計(jì)算每一個(gè)單頻子波對(duì)應(yīng)的梯度然后疊加形成編碼炮對(duì)應(yīng)的梯度.這種方法的關(guān)鍵技術(shù)是如何準(zhǔn)確地解耦混疊波場(chǎng).Dai等(2013)采用離散傅里葉變換的方法實(shí)現(xiàn)了波場(chǎng)解耦，Huang和Schuster(2018)也是借鑒于這種濾波方法，將同時(shí)源方法運(yùn)用到聲波全波形反演來(lái)處理海洋拖攬數(shù)據(jù).同時(shí)，Zhang等(2018)在聲波同時(shí)源全波形反演中引入了一種全新的波場(chǎng)解耦方法.在該方法中，相敏檢測(cè)(PSD)方法(Nihei and Li，2007)用來(lái)實(shí)現(xiàn)單頻波場(chǎng)的完全解耦.同時(shí)，在他們的工作中，也給出了一些關(guān)鍵的實(shí)現(xiàn)步驟，比如采用隨機(jī)的震源頻率，可變的編碼炮數(shù)等.基于這種PSD波場(chǎng)解耦的同時(shí)源編碼策略，Tromp和Bachmann(2009)實(shí)現(xiàn)了同時(shí)源層析成像.為了進(jìn)一步消除子波對(duì)反演的影響，Zhang等(2020b)將卷積型目標(biāo)函數(shù)引入同時(shí)源反演中.同時(shí)，一種地震數(shù)據(jù)正則化的策略用來(lái)增強(qiáng)同時(shí)源反演算法的魯棒性(Zhang et al.，2020a).

全波形反演就是同時(shí)使用波場(chǎng)的運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)特征.眾所周知，如果初始模型不是很理想，多信息匹配很容易使得全波形反演陷入局部極小值，所以開(kāi)展數(shù)據(jù)子集或者波場(chǎng)屬性反演就顯得十分必要(Sun and Schuster，1993；Fu et al.，2018；孟鴻鷹和劉貴忠，1999；馬堅(jiān)偉等，2000；董良國(guó)等，2015；梁展源等，2019).本文以同時(shí)源模擬與解耦的復(fù)波場(chǎng)為基礎(chǔ)，拓展出實(shí)波場(chǎng)反演、虛波場(chǎng)反演及復(fù)波場(chǎng)反演策略，并通過(guò)數(shù)值結(jié)果來(lái)對(duì)比三種反演的優(yōu)缺點(diǎn).具體文章結(jié)構(gòu)如下：在引言之后，回顧了無(wú)串?dāng)_的同時(shí)源反演方法；緊接著給出了實(shí)波場(chǎng)反演、虛波場(chǎng)反演以及復(fù)波場(chǎng)全波形反演的梯度表達(dá)式；接下來(lái)通過(guò)數(shù)值試驗(yàn)的方式說(shuō)明復(fù)波場(chǎng)中的實(shí)部和虛部信息對(duì)波形反演的影響；最后一部分給出結(jié)論與建議.

1 方法原理

1.1 無(wú)串?dāng)_的同時(shí)源全波形反演原理

對(duì)于無(wú)串?dāng)_同時(shí)源的FWI，采用最小化位于檢波點(diǎn)Xr的觀測(cè)數(shù)據(jù)dobs(Xr,Ω,t)和合成數(shù)據(jù)u(Xr,Ω,t;m)的殘差，其形式具體表示為：

(1)

其中，Ω表示角頻率，x表示計(jì)算區(qū)域，并且Xr?x.給定一系列位于Xs的單頻隨機(jī)信號(hào)s(Xs,Ω,t)，則編碼炮表示為：

(2)

(3)

其中Xs?x，L[·]是聲波正演模擬算子，相應(yīng)的伴隨方程表示為:

(4)

其中?表示伴隨算子，并且：

(5)

由于震源編碼方法對(duì)應(yīng)的梯度表達(dá)式和常規(guī)的梯度表達(dá)式一致，是正傳波場(chǎng)與反傳波場(chǎng)的特定形式的互相關(guān).為了分析串?dāng)_噪聲的影響，將梯度表達(dá)式簡(jiǎn)寫為：

(6)

其中，梯度簡(jiǎn)寫形式G(x,m)中的ψs(x,Ω)和ψr(x,Ω)分別表示正傳波場(chǎng)與反傳波場(chǎng).在式(6)的第二個(gè)等式中，第一項(xiàng)表示每一炮數(shù)據(jù)的自相關(guān)運(yùn)算，第二項(xiàng)表示當(dāng)前炮與其他炮的互相關(guān)運(yùn)算，即串?dāng)_噪聲.由(6)式可知，只有從混疊的波場(chǎng)中解耦出所有炮的波場(chǎng)，然后只計(jì)算當(dāng)前炮對(duì)應(yīng)的正反傳波場(chǎng)的自相關(guān)，才能從根本上消除串?dāng)_噪聲.

有效地解耦混疊波場(chǎng)是解決編碼噪聲的關(guān)鍵，因此引入相敏檢測(cè)方法來(lái)解決這個(gè)問(wèn)題(Zhang et al.，2018；Tromp and Bachmann，2019).相敏檢測(cè)方法的主要思想是利用兩個(gè)與未知信號(hào)具有相同頻率但彼此相位相差90°的參考信號(hào)(aref0°和aref90°)沿著時(shí)間積分來(lái)提取振幅為E相位為θ的目標(biāo)信號(hào)as(t).為了簡(jiǎn)潔表示，取參考信號(hào)為sin(ωit)和其90°相移信號(hào)為cos(ωit).空間某一點(diǎn)的混疊波場(chǎng)可以表示為:

as(t)=E1cos(ω1t+θ1)+E2cos(ω2t+θ2)+…

+Enscos(ωnst+θns),

(7)

使|ωi-ωj|≥pΔω(i,j∈ns,p∈Z)，且Δω=2π/T.這里的T是地震記錄的長(zhǎng)度，p是最小的整數(shù)，pΔω表示編碼信號(hào)之間的最小頻率間隔.當(dāng)定義最小的積分區(qū)間為Tp=T/p，對(duì)于任何整數(shù)q(p∈[1,p])倍的Tp區(qū)間qTp都可以作為合理的積分區(qū)間.通常，需要一定的時(shí)間等待波場(chǎng)達(dá)到穩(wěn)態(tài)狀況，所以q比p要小.在積分區(qū)間qTp上，參考信號(hào)sin(ωit)與混疊信號(hào)相乘并積分：

(8)

式(8)中，Ts表示波場(chǎng)達(dá)到穩(wěn)態(tài)的時(shí)刻.

相似地，參考信號(hào)的90°相移信號(hào)cos(ωit)與混疊信號(hào)相乘并積分：

(9)

所以與參考信號(hào)相同頻率的已知信號(hào)的振幅和相位可以解耦出來(lái)，其計(jì)算公式為：

(10)

式(8)和(9)積分主要利用了三角函數(shù)的正交性，只有與參考信號(hào)相同頻率的信號(hào)的積分不為零，與其他頻率的信號(hào)的積分為零，從而實(shí)現(xiàn)目標(biāo)信號(hào)的提取.依次對(duì)空間中的所有點(diǎn)做上述積分計(jì)算，就可以求解出整個(gè)空間某一頻率波場(chǎng)的振幅Ei和相位θi.同理，可以求取空間上其它不同頻率的空間波場(chǎng)所對(duì)應(yīng)的振幅和相位.則對(duì)應(yīng)于每一個(gè)頻率的空間正傳復(fù)波場(chǎng)的共軛波場(chǎng)為：

u(x,ωi,t;m)=Eicos(θi)-iEisin(θi),

(11)

反傳的復(fù)波場(chǎng)也具有類似的表達(dá)形式：

u′(x,ωi,t;m)=E′icos(θ′i)+iE′isin(θ′i).

(12)

然后通過(guò)梯度計(jì)算公式計(jì)算出目標(biāo)函數(shù)對(duì)模型參數(shù)的梯度，其中速度模型的梯度計(jì)算式為：

(13)

式(13)中，Re表示取向量的實(shí)數(shù)部分.有了目標(biāo)函數(shù)對(duì)應(yīng)模型參數(shù)的梯度，可以通過(guò)優(yōu)化算法求來(lái)迭代更新模型，具體的模型更新可以表示為:

mk+1=mk+tkΔmk+1,

(14)

式(14)中的k表示k次循環(huán)，tk表示選取的更新步長(zhǎng)，Δmk+1表示利用收斂算法根據(jù)梯度信息構(gòu)建的下降方向.本文采用共軛梯度算法來(lái)更新模型參數(shù)，其具體表達(dá)式為：

(15)

1.2 實(shí)波場(chǎng)、虛波場(chǎng)及復(fù)波場(chǎng)反演

無(wú)串?dāng)_噪聲的同時(shí)源反演的前提條件是混疊波場(chǎng)的完全解耦，而實(shí)現(xiàn)混疊波場(chǎng)(正傳波場(chǎng)和反傳波場(chǎng))的有效解耦方式是PSD方法.PSD方法只需要在波場(chǎng)外推的過(guò)程中沿著時(shí)間方向積分，如式(8)和(9)所示.波場(chǎng)信息的恢復(fù)如式(10)、(11)和(12)所示，從式子可以看到，PSD方法可以有效地得到空間波場(chǎng)的振幅和相位信息，并構(gòu)建出顯式的復(fù)波場(chǎng)，這為實(shí)波場(chǎng)、虛波場(chǎng)及復(fù)波場(chǎng)信息的反演奠定了基礎(chǔ).

基于式(11)、(12)和(13)的結(jié)果，可以自然地給出實(shí)波場(chǎng)、虛波場(chǎng)和復(fù)波場(chǎng)信息的反演梯度表達(dá)式.實(shí)波場(chǎng)反演是與波場(chǎng)相位的余弦有關(guān)，其梯度表達(dá)式為

(16)

其中g(shù)real表示實(shí)波場(chǎng)反演的梯度.虛波場(chǎng)反演是與波場(chǎng)相位的正弦有關(guān)，其梯度表達(dá)式為：

(17)

其中g(shù)imag表示虛波場(chǎng)反演的梯度.復(fù)波場(chǎng)的全波形反演方法，其梯度表達(dá)式為：

(18)

從式(18)可以看出，在全波形反演中，復(fù)波場(chǎng)的實(shí)部和虛部信息均參與了反演.

1.3 無(wú)串?dāng)_的同時(shí)源反演的算法流程

將PSD方法用到已知頻率的混疊波場(chǎng)的提取中，可以實(shí)現(xiàn)無(wú)串?dāng)_的編碼方法.通過(guò)提取當(dāng)前炮的正傳波場(chǎng)，以及對(duì)應(yīng)的反傳波場(chǎng)，通過(guò)對(duì)其進(jìn)行特定形式的求取當(dāng)前炮的梯度，以此方法求取編碼炮內(nèi)其他炮對(duì)應(yīng)的梯度，疊加形成總的編碼炮的梯度，進(jìn)而求取所有編碼炮對(duì)應(yīng)的梯度.最后通過(guò)前面所述的共軛梯度算法更新模型并迭代反演.總結(jié)來(lái)說(shuō)，基于PSD的聲波同時(shí)源全波形反演步驟大致可分為以下6步：

(a)波動(dòng)方程外推.選取一系列已知頻率的單頻信號(hào)并形成編碼源，采用式(3)進(jìn)行數(shù)值模擬.

(b)數(shù)據(jù)殘差(虛源)計(jì)算.按照式(5)所示，分別從觀測(cè)數(shù)據(jù)和合成數(shù)據(jù)中抽取某一炮的數(shù)據(jù)對(duì)應(yīng)的單頻地震記錄，在觀測(cè)系統(tǒng)的約束下求取該炮對(duì)應(yīng)的數(shù)據(jù)殘差，同理求取其他炮對(duì)應(yīng)的數(shù)據(jù)殘差并疊加形成編碼炮的數(shù)據(jù)殘差.

(c)伴隨態(tài)方程外推.使用(b)中求取的虛源代入式(4)進(jìn)行數(shù)值模擬.

(d)使用PSD方法進(jìn)行波場(chǎng)解耦.從正傳和反傳的混疊波場(chǎng)中分別提取編碼炮內(nèi)的所有炮對(duì)應(yīng)的波場(chǎng)響應(yīng)，然后通過(guò)每一炮的正傳波場(chǎng)的共軛波場(chǎng)與反傳波場(chǎng)相關(guān)求得每一炮的梯度，通過(guò)疊加所有炮的梯度并形成總梯度.其中三種反演對(duì)應(yīng)的梯度表達(dá)式如式(16)—(18)所示.

(e)通過(guò)共軛梯度算法更新模型.

(f)重復(fù)(a)—(e)直到達(dá)到設(shè)定的迭代次數(shù)或設(shè)定的目標(biāo)函數(shù)的門檻值，終止反演過(guò)程并輸出反演結(jié)果.

按照上述提到的6個(gè)步驟，設(shè)計(jì)了如算法1所示的具體程序?qū)崿F(xiàn)框架.根據(jù)該算法，完成同時(shí)源的實(shí)波場(chǎng)、虛波場(chǎng)及復(fù)波場(chǎng)波形反演.

算法1 無(wú)串?dāng)_的同時(shí)源波形反演算法輸入:觀測(cè)數(shù)據(jù)和觀測(cè)系統(tǒng)信息,初始速度以及反演參數(shù)輸出:反演結(jié)果和目標(biāo)函數(shù)值 do頻帶數(shù)循環(huán) do迭代次數(shù)循環(huán) do編碼炮數(shù)迭代·同時(shí)源模擬和波場(chǎng)解耦·伴隨源計(jì)算·伴隨波場(chǎng)外推和波場(chǎng)解耦 enddo編碼炮數(shù)迭代 計(jì)算(實(shí)波場(chǎng)、虛波場(chǎng)或復(fù)波場(chǎng)反演的)梯度并采用共軛梯度算法更新模型 enddo迭代次數(shù)循環(huán) enddo頻帶數(shù)循環(huán)

2 數(shù)值試驗(yàn)

在數(shù)值試驗(yàn)部分中，使用Marmousi模型和Overthurust模型進(jìn)行反演測(cè)試.程序設(shè)計(jì)為：CPU完成控制流和數(shù)據(jù)的輸入輸出，GPU負(fù)責(zé)核心的差分模擬計(jì)算.本文涉及的測(cè)試硬件為英特爾的i5-4460的CPU和GeForce RTX 3090的GPU.

2.1 Marmousi模型測(cè)試

首先采用Marmousi模型進(jìn)行反演試算，模型如圖1a所示.模型大小為6.63 km × 2.37 km，模型離散的空間步長(zhǎng)為10 m.其中采用交錯(cuò)網(wǎng)格有限差分方法實(shí)現(xiàn)波動(dòng)方程的離散，并在邊界處引入吸收邊界條件.有限差分的時(shí)間精度為2階，空間差分精度為4階.觀測(cè)數(shù)據(jù)是采用15 Hz的雷克子波生成均勻分布的80炮數(shù)據(jù)，炮間距為80 m，每一炮均有663個(gè)檢波點(diǎn)接收數(shù)據(jù).地震記錄的采樣間隔為1 ms，模擬時(shí)間為8 s.反演所使用的初始模型如圖1b所示，該初始模型是采用空間 2 km × 2 km對(duì)真實(shí)模型進(jìn)行空間平均得到的.采用多尺度反演策略，總共劃分5個(gè)頻帶進(jìn)行反演，總共迭代100次.其中Marmousi模型對(duì)應(yīng)的詳細(xì)反演參數(shù)在表1中列出.

圖1 (a)真實(shí)Marmousi模型；(b)初始模型

表1 Marmousi模型的反演參數(shù)

在反演之前，首先展示混疊波場(chǎng)的解耦過(guò)程.圖2a是編碼炮對(duì)應(yīng)的最大時(shí)刻的混疊波場(chǎng)，圖2b、c是解耦出的某一炮的振幅和相位信息；圖2d、e是解耦出的另外一炮的振幅和相位信息.從圖2可以看出，相敏檢測(cè)方法可以在波場(chǎng)模擬中，通過(guò)積分的方式完全解耦出不同頻率子波對(duì)應(yīng)的波場(chǎng)信息，實(shí)現(xiàn)混疊波場(chǎng)的完全解耦.對(duì)于反傳波場(chǎng)的波場(chǎng)解耦也是類似的過(guò)程.接下來(lái)比較三種反演策略計(jì)算的梯度.圖3a—c分別展示了實(shí)波場(chǎng)反演、虛波場(chǎng)反演以及復(fù)波場(chǎng)反演在相同迭代次數(shù)的梯度.對(duì)比圖3中的3個(gè)子圖，可以看到實(shí)波場(chǎng)反演的梯度有著明顯的振幅能量變化，虛波場(chǎng)反演的梯度更多是層的位置信息，整體能量均衡性比較好，復(fù)波場(chǎng)反演的梯度則同時(shí)具備實(shí)波場(chǎng)反演和虛波場(chǎng)反演的特征，梯度信息更加豐富，構(gòu)造層次感比較強(qiáng).

圖2 (a)混疊波場(chǎng)，相敏檢測(cè)解耦出的某一炮的(b)振幅信息和(c)相位信息，以及解耦出的另外一炮的(d)振幅信息和(e)相位信息

圖3 Marmousi模型的梯度

在分析完梯度特征的基礎(chǔ)上，開(kāi)展反演試算.對(duì)于三種反演方法來(lái)說(shuō)，除了選用的梯度計(jì)算不同之外，其余的反演參數(shù)均相同.此處展示了三種反演策略在第二個(gè)頻帶、第五個(gè)頻帶的反演結(jié)果，如圖4所示.從最終的反演結(jié)果可以看到，三種方法均取得了比較好的反演結(jié)果，模型中的構(gòu)造信息都得到了有效的恢復(fù).通過(guò)仔細(xì)對(duì)比，相比于復(fù)波場(chǎng)全波形反演，無(wú)論是實(shí)波場(chǎng)反演還是虛波場(chǎng)反演，其反演結(jié)果的背景速度均有一定的模糊作用，分辨率有所降低.為了定量的評(píng)估反演結(jié)果的可靠性，本文采用式(19)來(lái)衡量計(jì)算結(jié)果的誤差：

圖4 Marmousi模型的反演結(jié)果

(19)

其中erra表示絕對(duì)誤差，errr表示相對(duì)誤差，nx和nz表示模型在二維空間的離散點(diǎn)數(shù)，|·|表示取絕對(duì)值運(yùn)算.此處計(jì)算了反演結(jié)果的絕對(duì)誤差，初始模型對(duì)應(yīng)的絕對(duì)誤差為263.68 m·s-1，實(shí)波場(chǎng)反演結(jié)果的絕對(duì)誤差為203.82 m·s-1，虛波場(chǎng)反演結(jié)果的誤差198.18 m·s-1，復(fù)波場(chǎng)反演的誤差為193.90 m·s-1.對(duì)比最終反演結(jié)果的絕對(duì)誤差，可以發(fā)現(xiàn)復(fù)波場(chǎng)反演的精度最高，其次是虛波場(chǎng)反演.為了更清楚地展示反演結(jié)果，抽取位于x=2.2 km和x=4.2 km的反演結(jié)果并分別展示在圖5a、b中.從反演剖面上可以看出，復(fù)波場(chǎng)反演在某些局部位置可以得到更加準(zhǔn)確的速度值.

圖5 Marnousi模型位于(a)x=2.2 km和(b)x=4.2 km處的反演剖面對(duì)比

圖6a、b分別展示了三種方法對(duì)應(yīng)的目標(biāo)函數(shù)值及模型誤差值隨著迭代次數(shù)變化的曲線，此處模型誤差值是通過(guò)式(19)中的相對(duì)誤差計(jì)算得到的.從圖6中可以看出，隨著迭代次數(shù)的增加，目標(biāo)函數(shù)和模型的相對(duì)誤差都呈現(xiàn)明顯的下降趨勢(shì).就目標(biāo)函數(shù)的下降曲線而言，盡管虛波場(chǎng)反演的目標(biāo)函數(shù)下降曲線與復(fù)波場(chǎng)反演的曲線吻合度較高一些，但三種方法的差異不大；從模型參數(shù)的下降曲線可以看出，反演的初始階段，虛波場(chǎng)反演與復(fù)波場(chǎng)反演的誤差有著很高的一致性，且優(yōu)于實(shí)波場(chǎng)反演方法；隨著迭代次數(shù)的進(jìn)一步增加，相位反演的模型收斂性變差，這說(shuō)明復(fù)波場(chǎng)反演中的實(shí)波場(chǎng)信息在高波數(shù)反演中貢獻(xiàn)突出.

圖6 目標(biāo)函數(shù)值(a)和模型誤差值(b)隨著迭代次數(shù)變化曲線

2.2 Overthrust模型測(cè)試

接下來(lái)，在Overthrust模型上進(jìn)行反演測(cè)試.模型大小170×800，網(wǎng)格間距10 m.采用主頻為15 Hz的雷克子波生成100炮數(shù)據(jù)，數(shù)據(jù)的時(shí)間采樣點(diǎn)為1 ms，地震記錄長(zhǎng)8 s.模型的真實(shí)速度如圖7a所示，反演所使用的初始模型如圖7b所示.初始速度模型是一個(gè)隨著深度遞增的線性速度模型.相比于Marmousi模型反演所用的初始模型，此處使用的初始模型比較粗糙，因此會(huì)使得反演目標(biāo)函數(shù)的非線性增強(qiáng)，所以依舊使用多尺度反演策略來(lái)降低反演的多解性.反演中使用了5個(gè)反演尺度，每個(gè)頻帶的頻率范圍及編碼炮個(gè)數(shù)不相同，具體的反演參數(shù)如表2所示.

圖7 (a)真實(shí)Overthrust模型；(b)初始模型

表2 Overthrust模型的反演參數(shù)

依舊設(shè)計(jì)了三種反演方法試算，在這三種反演算法中除了梯度計(jì)算不同外，其余的反演參數(shù)均相同.最終的反演結(jié)果如圖8a—c所示，其中圖8a—c分別為實(shí)波場(chǎng)反演、虛波場(chǎng)反演以及復(fù)波場(chǎng)反演結(jié)果.從反演結(jié)果可以看到，在初始模型比較差的情況下，實(shí)波場(chǎng)反演不能得到有效的模型參數(shù)，尤其是在黑框所標(biāo)記的逆沖斷層處成像效果很差；對(duì)于虛波場(chǎng)反演而言，即使是初始速度較差的情況下，依舊能取得比較好的反演結(jié)果，并且其成像分辨率基本上與復(fù)波場(chǎng)反演的分辨率相當(dāng)，這就說(shuō)明虛波場(chǎng)反演對(duì)速度模型的敏感度較低.由于全波形反演是實(shí)波場(chǎng)和虛波場(chǎng)信息的綜合使用，這說(shuō)明全波形反演的低波數(shù)成分主要由虛波場(chǎng)信息來(lái)恢復(fù)，從而避免陷入局部極值的問(wèn)題；隨著模型的低波數(shù)成分被有效地恢復(fù)，全波形反演可以綜合使用實(shí)波場(chǎng)和虛波場(chǎng)信息實(shí)現(xiàn)高分辨率的反演結(jié)果.于是，本文開(kāi)展了一種組合的反演試算，即在低波數(shù)反演中，采用虛波場(chǎng)反演，在高波數(shù)反演中，使用復(fù)波場(chǎng)反演，這樣可以降低實(shí)波場(chǎng)反演對(duì)低波數(shù)的依賴性.組合反演的最終結(jié)果如圖8d所示，從反演結(jié)果上可以看到，相比于虛波場(chǎng)反演，反演精度有了明顯的提升，并達(dá)到了與復(fù)波場(chǎng)反演接近的精度.為了更好地對(duì)比編碼方法的反演結(jié)果，圖8e展示了常規(guī)序列源的反演結(jié)果.對(duì)比發(fā)現(xiàn)，采用序列源帶寬子波反演方法，反演結(jié)果在層內(nèi)具有更好的光滑性，但是層界面不是很清晰.相似地，初始結(jié)果和五種反演結(jié)果的絕對(duì)誤差也被計(jì)算：初始模型對(duì)應(yīng)的絕對(duì)誤差為439.80 m·s-1，實(shí)波場(chǎng)反演結(jié)果的絕對(duì)誤差為349.11 m·s-1，虛波場(chǎng)反演結(jié)果的誤差216.80 m·s-1，復(fù)波場(chǎng)反演的誤差為166.45 m·s-1，組合反演結(jié)果的誤差為191.28 m·s-1，以及常規(guī)序列源反演的誤差為207.10 m·s-1.對(duì)比最終反演結(jié)果的絕對(duì)誤差，可以發(fā)現(xiàn)虛波場(chǎng)反演的精度可以通過(guò)組合反演策略得到提升；同時(shí)，也可以從圖9中的局部反演結(jié)果對(duì)比得到相同的結(jié)論.

圖8 Overthrust模型的反演結(jié)果

圖9 Overthrust模型的反演結(jié)果局部顯示

相似地，五種反演方法對(duì)應(yīng)的目標(biāo)函數(shù)值及模型誤差值隨著迭代次數(shù)的變化曲線被繪制在圖10中.在圖10a，可以看到在五個(gè)反演尺度中，目標(biāo)函數(shù)均呈現(xiàn)出比較好的下降趨勢(shì)；除了常規(guī)序列源反演包含著豐富的數(shù)據(jù)量，所以目標(biāo)函數(shù)的殘差整體比編碼方法的殘差要大；在編碼方法中，相比于其他幾種反演方法，實(shí)波場(chǎng)反演的目標(biāo)函數(shù)在低頻下降緩慢，同時(shí)對(duì)比圖10b，實(shí)波場(chǎng)反演的模型不收斂，這說(shuō)明實(shí)波場(chǎng)反演已經(jīng)陷入局部極值；雖然復(fù)波場(chǎng)反演在低頻有著最大的下降率，但是隨著迭代次數(shù)的增加，其目標(biāo)函數(shù)的下降曲線逐漸與虛波場(chǎng)反演相一致；對(duì)于組合的反演方案，在開(kāi)始階段與虛波場(chǎng)反演吻合度很高，隨著迭代次數(shù)的增加，其目標(biāo)函數(shù)值逐漸低于虛波場(chǎng)反演，直至最后的頻帶，組合反演的目標(biāo)函數(shù)值是四種反演算法中最低的，這說(shuō)明這種策略可以降低目標(biāo)函數(shù)的復(fù)雜度，使得目標(biāo)函數(shù)下降的更多.

從模型下降曲線(圖10b)可以看出，實(shí)波場(chǎng)反演的模型收斂性最差，虛波場(chǎng)反演、組合反演，復(fù)波場(chǎng)反演，和序列源反演均呈現(xiàn)出比較好的模型收斂性，且在模型的收斂性方面復(fù)波場(chǎng)反演優(yōu)于組合反演，組合反演優(yōu)于虛波場(chǎng)反演和序列源反演.值得注意的是，第四個(gè)頻帶處，復(fù)波場(chǎng)反演出現(xiàn)了小幅度的模型誤差增加，這充分說(shuō)明了復(fù)波場(chǎng)反演的高度非線性性.

圖10 目標(biāo)函數(shù)值(a)和模型誤差值(b)隨著迭代次數(shù)變化曲線

3 結(jié)論與討論

(1)在地震建模中，全波形反演具有高精度建模的潛在優(yōu)勢(shì).主要有兩方面挑戰(zhàn)限制了該方法的大規(guī)模應(yīng)用.第一個(gè)方面是計(jì)算效率的問(wèn)題，本文引入的無(wú)串?dāng)_的同時(shí)源反演方法，可以有效地提高計(jì)算效率且不影響成像質(zhì)量；第二個(gè)就是該反演方法是高度非線性的，因此對(duì)初始模型有著很強(qiáng)的依賴性.本文從復(fù)波場(chǎng)信息解耦的角度出發(fā)，給出了實(shí)波場(chǎng)、虛波場(chǎng)以及復(fù)波場(chǎng)反演的方法，在反演中分別或者共同使用復(fù)波場(chǎng)的實(shí)部和虛部信息，這對(duì)復(fù)雜介質(zhì)速度建模有著重要的意義.

(2)從實(shí)波場(chǎng)反演、虛波場(chǎng)反演、到復(fù)波場(chǎng)反演的試算可以看出，當(dāng)初始速度比較準(zhǔn)確時(shí)，三種方法都能得到比較好的反演結(jié)果；當(dāng)初始模型準(zhǔn)確性較低時(shí)，實(shí)波場(chǎng)反演很快就陷入局部極值；而虛波場(chǎng)反演依舊可以得到比較好的反演結(jié)果；這說(shuō)明實(shí)波場(chǎng)反演對(duì)初始模型的依賴性較強(qiáng)，綜合使用實(shí)波場(chǎng)和虛波場(chǎng)反演的全波形反演，在低頻反演階段主要是虛波場(chǎng)反演的作用保證反演不陷入局部極值.因此，一種組合的反演策略，低頻采用虛波場(chǎng)反演，高頻采用復(fù)波場(chǎng)反演可以在復(fù)雜模型建模中有一定的應(yīng)用潛力.

(3)本文目前只開(kāi)展了基于聲波方程的復(fù)波場(chǎng)實(shí)部和虛部反演研究.由于彈性波場(chǎng)更加復(fù)雜，目標(biāo)函數(shù)的局部極值更多，下一步工作可以考慮在彈性波反演中做進(jìn)一步研究，來(lái)增強(qiáng)彈性波動(dòng)方程波形反演建模能力.