復波場的實部和虛部對同時源波形反演的影響

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

1 Geology and Geophysics Program，Missouri University of Science and Technology，Rolla，MO 65409，USA 2 深部巖土力學與地下工程國家重點實驗室，中國礦業大學力學與土木工程學院，徐州 221116 3 東北石油大學地球科學學院，大慶 163318

0 引言

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

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

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

1 方法原理

1.1 無串擾的同時源全波形反演原理

對于無串擾同時源的FWI，采用最小化位于檢波點Xr的觀測數據dobs(Xr,Ω,t)和合成數據u(Xr,Ω,t;m)的殘差，其形式具體表示為：

(1)

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

(2)

(3)

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

(4)

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

(5)

由于震源編碼方法對應的梯度表達式和常規的梯度表達式一致，是正傳波場與反傳波場的特定形式的互相關.為了分析串擾噪聲的影響，將梯度表達式簡寫為：

(6)

其中，梯度簡寫形式G(x,m)中的ψs(x,Ω)和ψr(x,Ω)分別表示正傳波場與反傳波場.在式(6)的第二個等式中，第一項表示每一炮數據的自相關運算，第二項表示當前炮與其他炮的互相關運算，即串擾噪聲.由(6)式可知，只有從混疊的波場中解耦出所有炮的波場，然后只計算當前炮對應的正反傳波場的自相關，才能從根本上消除串擾噪聲.

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

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

+Enscos(ωnst+θns),

(7)

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

(8)

式(8)中，Ts表示波場達到穩態的時刻.

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

(9)

所以與參考信號相同頻率的已知信號的振幅和相位可以解耦出來，其計算公式為：

(10)

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

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

(11)

反傳的復波場也具有類似的表達形式：

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

(12)

然后通過梯度計算公式計算出目標函數對模型參數的梯度，其中速度模型的梯度計算式為：

(13)

式(13)中，Re表示取向量的實數部分.有了目標函數對應模型參數的梯度，可以通過優化算法求來迭代更新模型，具體的模型更新可以表示為:

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

(14)

式(14)中的k表示k次循環，tk表示選取的更新步長，Δmk+1表示利用收斂算法根據梯度信息構建的下降方向.本文采用共軛梯度算法來更新模型參數，其具體表達式為：

(15)

1.2 實波場、虛波場及復波場反演

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

基于式(11)、(12)和(13)的結果，可以自然地給出實波場、虛波場和復波場信息的反演梯度表達式.實波場反演是與波場相位的余弦有關，其梯度表達式為

(16)

其中greal表示實波場反演的梯度.虛波場反演是與波場相位的正弦有關，其梯度表達式為：

(17)

其中gimag表示虛波場反演的梯度.復波場的全波形反演方法，其梯度表達式為：

(18)

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

1.3 無串擾的同時源反演的算法流程

將PSD方法用到已知頻率的混疊波場的提取中，可以實現無串擾的編碼方法.通過提取當前炮的正傳波場，以及對應的反傳波場，通過對其進行特定形式的求取當前炮的梯度，以此方法求取編碼炮內其他炮對應的梯度，疊加形成總的編碼炮的梯度，進而求取所有編碼炮對應的梯度.最后通過前面所述的共軛梯度算法更新模型并迭代反演.總結來說，基于PSD的聲波同時源全波形反演步驟大致可分為以下6步：

(a)波動方程外推.選取一系列已知頻率的單頻信號并形成編碼源，采用式(3)進行數值模擬.

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

(c)伴隨態方程外推.使用(b)中求取的虛源代入式(4)進行數值模擬.

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

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

(f)重復(a)—(e)直到達到設定的迭代次數或設定的目標函數的門檻值，終止反演過程并輸出反演結果.

按照上述提到的6個步驟，設計了如算法1所示的具體程序實現框架.根據該算法，完成同時源的實波場、虛波場及復波場波形反演.

算法1 無串擾的同時源波形反演算法輸入:觀測數據和觀測系統信息,初始速度以及反演參數輸出:反演結果和目標函數值 do頻帶數循環 do迭代次數循環 do編碼炮數迭代·同時源模擬和波場解耦·伴隨源計算·伴隨波場外推和波場解耦 enddo編碼炮數迭代 計算(實波場、虛波場或復波場反演的)梯度并采用共軛梯度算法更新模型 enddo迭代次數循環 enddo頻帶數循環

2 數值試驗

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

2.1 Marmousi模型測試

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

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

表1 Marmousi模型的反演參數

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

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

圖3 Marmousi模型的梯度

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

圖4 Marmousi模型的反演結果

(19)

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

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

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

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

2.2 Overthrust模型測試

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

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

表2 Overthrust模型的反演參數

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

圖8 Overthrust模型的反演結果

圖9 Overthrust模型的反演結果局部顯示

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

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

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

3 結論與討論

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

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

(3)本文目前只開展了基于聲波方程的復波場實部和虛部反演研究.由于彈性波場更加復雜，目標函數的局部極值更多，下一步工作可以考慮在彈性波反演中做進一步研究，來增強彈性波動方程波形反演建模能力.