基于時域加權的拉普拉斯—頻率域彈性波全波形反演

劉張聚 童思友 方云峰 賈君蓮

(①中國海洋大學海底科學與探測技術教育部重點實驗室，山東青島 266100；②青島海洋科學與技術試點國家實驗室海洋礦產資源評價與探測技術功能實驗室，山東青島 266061；③浙江大學地球科學學院，浙江杭州 310027；④自然資源部第二海洋研究所海底科學重點實驗室，浙江杭州310012；⑤東方地球物理公司物探技術研究中心，河北涿州 072751；⑥東方地球物理公司研究院，河北涿州 072751)

0 引言

全波形反演(Full Waveform Inversion，FWI)是利用地震記錄的全波形信息，通過匹配模擬記錄與實際記錄反演地下介質參數的一種高分辨率建模方法。Tarantola[1-2]最先提出了時間域全波形反演方法，該方法基于廣義最小二乘理論，使用理論記錄與實際記錄的殘差來構建目標函數。此后，全波形反演得到了全面的發展。Bunks等[3]提出了時間域多尺度反演，降低了反演的非線性影響；Pratt[4]提出了在頻率域進行全波形反演，不需要對地震記錄進行預處理，且使用有限個頻率信息即可得到高精度的反演結果；Sirgue等[5]綜合時間域反演和頻率域反演的優勢，提出了時間—頻率混合域反演方法，避免了因時間域求梯度所需要的震源波場重構以及頻率域正演所需計算機內存過大等問題。

另外，還有一種研究思路是降低全波形反演的非線性，降低陷入局部極小值的可能性，從而在缺失低頻或者初始模型不準確的情況下，反演依然能逼近正確結果。Choi等[28]提出的基于全局互相關全波形反演可有效降低反演的非線性，其強調的是記錄之間整體的相關程度，并非是記錄中點與點的匹配程度；Liu等[29]詳細對比了全局互相關目標函數與L2范數目標函數之間的差異，說明了互相關目標函數的優勢；Tao等[30]將全局互相關應用到了全波形層析，其反演結果同樣優于基于旅行時的層析成像；李青陽等[31]將互相關與反射波波形反演相結合降低了反演的非線性。Luo等[32]建立了反褶積目標函數，與互相關目標函數類似，降低了陷入局部極小值的可能性，但相對而言對帶限或非脈沖源的靈敏度降低。Warner等[33-35]建立了自適應波形反演理論，通過將維納濾波器引入目標函數，尋找將兩種記錄相匹配的濾波器，再建立濾波器的目標函數，使其逐漸逼近零延遲的單位脈沖函數；Guasch等[36]將該方法推廣到三維，并且在缺失低頻的情況下得到了較好的反演結果。van Leeuwen等[37]提出了波形重構反演，將波動方程本身作為懲罰項引入到目標函數中，通過重構波場更新模型參數，將解空間擴展到了數據空間和模型空間，有效減小了局部極小值的影響；Li等[38]將該方法引入時間域，提出了時間域波形重構法，相比于頻率域可有效減少計算成本。

針對如何在缺失低頻信息的情況下建立較高精度的初始模型這一難題，本文建立了一種基于時域加權的拉普拉斯—頻率域彈性波全波形反演方法(Time-Weighted Laplace-Frequency Domain Elastic Full Waveform Inversion，TWLF-EFWI)。該方法在時間—頻率域彈性波全波形反演的基礎上，通過引入拉普拉斯阻尼因子，降低了全波形反演對于初始模型的依賴程度，形成了時間—拉普拉斯—頻率域的混合域全波形反演方法。TWLF-EFWI方法仍使用傳統的L2范數作為目標函數，為消除因引入拉普拉斯阻尼因子帶來的地震道能量不均衡、遠炮檢距地震道能量衰減嚴重的影響，在目標函數中引入了時域加權因子，消除了拉普拉斯阻尼因子的負面影響。Marmousi 2部分模型測試結果表明，在缺失低頻信息以及初始模型精度較低的情況下，本文方法可得到含有模型宏觀構造的反演結果，足以滿足后續反演對初始模型的精度要求。

1 全波形反演基本理論

1.1 頻率域彈性波全波形反演

全波形反演可以描述為一個基于地震全波場模擬的數據擬合過程，目的是實現正演記錄和實際觀測記錄的完全擬合，也可以將全波形反演重新轉換為線性化的最小二乘問題。

在頻率域，波動方程矩陣形式可以表示為

AU=S

(1)

式中：A為阻抗矩陣；U為全部網格點的波場值所組成的列向量；S為震源項。

基于廣義最小二乘的全波形反演目標函數定義為

(2)

式中：m為模型參數；Ns為總炮數；Δd=dcal-dobs，為模擬地震記錄dcal與實際地震記錄dobs的差；上標“*”表示矩陣的伴隨形式。全波形反演就是從目標函數出發，求解模型的更新量，然后通過對模型的迭代更新，使目標函數達到最小，或者使前后兩次的目標函數的相對誤差達到限定值，進而得到與其對應的模型參數的最優解。

將式(2)中的目標函數對m求導，可以導出目標函數在頻率域關于模型參數的梯度為

(3)

式中J為Fréchet導數矩陣，可由波動方程兩邊同時對參數m求導得到

(4)

式中R表示波場U到接收記錄d的映射關系。則頻率域梯度可表示為

(5)

上式表明目標函數的梯度由正演波場U、阻抗矩陣對模型參數的導數?A/?m和以地震記錄殘差共軛為震源的正演波場A-1RTΔd*三部分組成。

對于二維彈性波方程，在頻率域表達式為

(6)

式中：Ufx和Ufz分別為頻率域速度波場的水平和垂直分量；Fx和Fz為對應的頻率域震源項；ω為角頻率；ρ為密度；λ和μ為拉梅常數。

聯合式(5)和式(6)，可得目標函數關于拉梅常數的梯度在頻率域的表達式

(7)

(8)

式中vP和vS分別為縱波速度和橫波速度。

1.2 拉普拉斯—頻率域彈性波全波形反演

拉普拉斯—頻率域的二維彈性波方程與頻率域的二維彈性波方程極其相似，為

(9)

式中：ULx和ULz分別為拉普拉斯—頻率域波場的水平和垂直分量；Lx和Lz為對應的拉普拉斯—頻率域震源；s=σ+iω，其中σ為拉普拉斯衰減常數。

對比式(9)與式(6)可以發現，二者的唯一區別是ω與s，所以拉普拉斯—頻率域也可以看作是“復頻率域”。兩個域的目標函數梯度也具有相同的形式，在拉普拉斯—頻率域具體形式為

(10)

同樣，根據式(8)，可得目標函數關于縱波速度和橫波速度的梯度。

1.3 混合域彈性波全波形反演

全波形反演在時間域和頻率域里的基本理論是相通的，二者并無差異。時間—頻率混合域全波形反演就是在時間域模擬正、反傳波場，然后轉換到頻率域計算梯度。兩種域的結合只需使用離散傅里葉變換(DFT)即可實現，即在時間域求取波場時，在時間循環中加入DFT，當循環結束時也就得到了對應的頻率域的單頻波場，之后再用頻率域單頻波場計算梯度。混合域反演結合了時間域模擬波場的高效性和頻率域計算梯度的靈活性，其流程如圖1所示。

DFT和逆離散傅里葉變換(IDFT)的表達式為

(11)

式中：u(kΔt)和UF(ωn)分別為時間域波場和頻率域波場；Δt為時間采樣間隔；NT為時間域序列長度；NF為頻率域序列長度。

圖1 時間—頻率域全波形反演流程

在上述混合域的基礎上，使用拉普拉斯—頻率域代替頻率域，形成了時域—拉普拉斯—頻率域的全波形反演(圖2)。與圖1流程不同的是，需要將DFT改為離散拉普拉斯變換(DLT)進行時域與拉普拉斯—頻率域間的轉換，并且在求取伴隨震源時需要做一次DLT和逆離散拉普拉斯變換(IDLT)。另外本文將伴隨震源在時域中進行了加權處理，后文會對此進行詳細描述。DLT和IDLT的表達式為

(12)

式中：UL(sn)為拉普拉斯—頻率域波場；NL為拉普拉斯—頻率域序列長度。

圖2 時間—拉普拉斯—頻率域全波形反演流程

2 時域指數衰減波場

相對于傅里葉變換，拉普拉斯變換只是增加了一個時間指數衰減項。正是這一指數衰減項，在拉普拉斯—頻率域反演降低了對低頻信息的依賴。

時間域信號u(t)的拉普拉斯變換為

(13)

將拉普拉斯變換中的衰減項單獨取出，在時間域分析該項對信號頻率的影響。時間域信號衰減可表示為

(14)

在時間域的乘積等價于頻率域的卷積，衰減后信號的離散傅里葉變換與原始信號和阻尼因子在頻域的循環卷積相同，所以衰減后信號的零頻、低頻分量是原始信號頻率分量的加權和。Shin等[20]曾描述衰減波場中包含 “mirage-like”(海市蜃樓式)低頻成分；Ha等[23]詳細證明了這種衰減波場中確實存在零頻和低頻信息，并證明了這些信息的有效性，指出這些信息可以看作是原始信號的模糊處理結果。

為展示阻尼因子對低頻成分的作用和效果，分別用地震子波及其振幅譜、地震記錄及其相位譜加以說明。圖3對比了子波經不同阻尼因子衰減后的結果。原始子波(圖3a)是雷克子波經過切比雪夫濾波，去除了0～3Hz低頻信息。阻尼因子為2、8和20的衰減子波如圖3b～圖3d所示。原始子波的最高波峰在中間，兩側分別有兩個較小的波峰(圖3a)；當阻尼因子為20時，子波的第一個波峰已經超過了中間的波峰值，第三個波峰幾乎消失(圖3d)。子波衰減是隨著時間變化的，到達時間越早振幅增益作用越強，到達時間越晚振幅抑制作用越強，而且阻尼因子越大，作用強度越大。

圖3 不同阻尼因子衰減后的子波歸一化對比(a)原始； (b)σ=2； (c)σ=8； (d)σ=20

圖4為圖3四個子波的振幅譜，可以看出，隨著阻尼因子的增大，子波中不存在的0～3Hz低頻信息慢慢出現，當阻尼因子為20時，頻譜中0～3Hz頻率成分已十分明顯。所以，對信號進行阻尼衰減可以重構出低頻信息。

圖5a為單炮地震記錄，圖5b為濾除3Hz以下頻率成分后的記錄，圖6為兩個記錄的相位譜。相比于原始記錄，濾波后記錄的相位譜(圖6b)中，0～3Hz的部分已經扭曲畸變。然后將兩個記錄分別進行阻尼因子為2、8和20的衰減，其相位譜如圖7所示。濾波記錄經過衰減后的相位譜(圖7右)，0～3Hz部分不再扭曲，與原始記錄衰減后的相位譜(圖7左)明顯相似。可以認為，經過衰減后的低頻成分包含有地下介質的有效信息，可以使用該低頻成分進行反演。

圖4 不同阻尼因子衰減后的子波歸一化振幅譜對比(a)原始； (b)σ=2； (c)σ=8； (d)σ=20

圖5 原始地震記錄(a)及其濾除0～3Hz低頻成分后的記錄(b)

圖6 原始地震記錄(a)及其濾波后記錄(b)的歸一化相位譜

圖7 原始地震記錄(左)及其濾波后記錄(右)經不同阻尼因子衰減后的歸一化相位譜(a)σ=2； (b)σ=8； (c)σ=20

3 基于L2范數的時域加權目標函數

當地震記錄中引入時間阻尼因子后，遠炮檢距地震道的能量會被衰減得很小。Shin等[20]為了消除阻尼因子產生的這一負面影響，采用了對數目標函數，即

(15)

雖然對數目標函數可以消除指數衰減的負面影響，但對數目標函數本身存在固有的缺陷，即對小功率譜波場以及波場的振幅變化十分敏感，反演穩定性較差。為此，Jun等[39]提出了加權的L2目標函數，史才旺[27]提出了基于L1范數的對數目標函數。基于前人的研究，本文使用了基于時間域加權的L2目標函數。

本文的目標函數首先還是以傳統的L2差值目標函數為核心，但是與前人不同的是，本方法在時間域進行加權處理，因為相較于頻率域來說，時間域對于記錄的操作更加方便和直觀。時間域目標函數為

(16)

(17)

不論是用對數目標函數還是加權系數目標函數，其目的都是消除阻尼因子的負面影響，增強遠炮檢距地震道的能量，使其在反演中發揮更大的作用。本文應用地震道差值模的倒數作為加權系數，即

(18)

4 模型測試

使用Marmousi 2部分模型測試本文的反演方法。為節約反演運算時間，將Marmousi 2縱波速度模型進行了截取和抽稀，得到了網格節點數為401×160、網格尺寸為15m×15m的模型作為真實縱波速度模型(圖8a)；縱橫波速度比設為1.7，獲得其橫波速度模型(圖8b)；密度設為常數2000kg/m3。為提高反演效率，默認陸地淺部水平層的速度已知(vP=1500m/s)。

正演采用主頻為6Hz的雷克子波作為激發震源，記錄長度為6s，采樣間隔為1.5ms。炮點和檢波點均勻布置在地面，炮點范圍為0～6km，共41炮；檢波點范圍為0.15～5.85km，每炮381道。圖9為模擬的第20炮兩分量記錄。

圖8 Marmousi 2部分模型(a)縱波速度； (b)橫波速度

使用了如圖10所示的一維線性梯度模型作為初始模型，模型中不包含任何構造信息。采用多尺度反演策略，參照激發子波的主頻，選取了4個頻率組(表1)，前兩組的頻率低于3Hz，后兩組的頻率高于3Hz，每組分別迭代40次。

表1 反演頻率分組參數

圖9 第20炮觀測記錄(a)水平分量； (b)垂直分量

全頻帶數據常規時間—頻率域反演最終結果如圖11所示。與圖8對比可以看出，反演結果十分接近真實模型，兩側層狀介質、中部的異常體、深部的斷裂和背斜構造都得到了精確重構。為了更為直觀、定量地反映反演效果，抽取x=3.9km處(中部隆起處)的縱、橫波速度曲線(圖12)，可以發現，反演結果十分接近真實值。

濾除模擬記錄3Hz以下的成分，以圖10的一維線性梯度模型作為初始模型進行常規時域—頻率域反演。因記錄中的可用頻率成分全部在3Hz以上，所以選取表1中的后兩組頻率進行反演。

圖10 初始速度模型(a)縱波速度； (b)橫波速度

圖11 全頻帶常規時間—頻率域反演結果(a)縱波速度； (b)橫波速度

缺失低頻成分的常規時間—頻率域反演結果如圖13所示。因沒有精確的低波數更新，所以很難避免周波跳躍現象，以至于陷入局部極值。與真實模型(圖8)對比可見，模型的大尺度構造幾乎沒有得到重構，只得到了一些小尺度的細節更新，并且產生了一些構造假象；斷層、斷塊和模型深部的不整合面都不可識別。對比圖13與圖11可見，低頻信息對常規反演的重要性。在初始模型精度不高的情況下，就需要使用記錄中的低頻信息重構出模型的大尺度構造，進而得到最終的高精度反演結果。若記錄中缺失低頻信息，反演就會極易陷入局部極小，導致反演失敗。

應用圖10的初始模型對缺失0～3Hz低頻的模擬記錄進行時間—拉普拉斯—頻率域全波形反演。因低頻已經通過拉普拉斯變換重構，所以可以使用0～3Hz的頻率成分，經測試后選擇了5個尺度，每個尺度進行10次迭代，具體參數如表2所示。

圖14為時間—拉普拉斯—頻率域反演的結果，可以看到，模型中大尺度的構造已經得到恢復，其形狀與真實模型的輪廓十分相近。然后將該反演結果作為初始模型，進行常規時間—頻率域反演，反演尺度參數與表1中的后兩組相同。由圖15的反演結果可以看出，與真實模型十分接近，各個層位都得到了精確的反演，斷層、斷塊等構造也十分清晰。與圖11相比，兩種結果相差無幾，說明該方法可以在缺失低頻信息的情況下，達到含低頻信息時的反演精度。

圖12 x=3.9km處全頻帶常規時間—頻率域反演速度曲線與真實及初始速度曲線的對比

為了測試本文方法對密度反演的適用性，將常密度修改為空間變化(圖16)，與縱、橫波速度同時參與反演。其他正、反演參數與以上相同。

圖13 缺失低頻成分時的常規時間—頻率域反演結果(a)縱波速度； (b)橫波速度

圖14 缺失低頻成分時的時間—拉普拉斯—頻率域反演結果(a)縱波速度； (b)橫波速度

對缺失0～3Hz低頻成分的記錄進行時間—拉普拉斯—頻率域全波形反演，結果如圖17所示，可以看出，縱、橫波速度和密度模型的大尺度構造同樣

表2 復頻率域反演頻率分組及參數

圖15 以時間—拉普拉斯—頻率域反演結果為初始模型的常規時間—頻率域反演結果

圖16 密度模型

得到了比較好的恢復。將該反演結果作為初始模型，進行時間—頻率域反演，結果如圖18所示。

從圖18可以看出，縱、橫波速度的反演結果同樣十分接近真實模型。根據Ben-Hadj-Ali等[40]的模型誤差公式，計算圖15的縱、橫波速度相對誤差分別為0.2437、0.2355，圖18的縱、橫波速度相對誤差分別為0.2497、0.2380，說明密度是否參與反演對縱、橫波速度的反演結果幾無影響，同時也說明本文方法同樣適用于密度反演。對于密度項，由于與其他參數的耦合作用，其最終反演結果比速度反演略差，但是可以看出本文方法的密度反演結果仍有較高的可信度，可以體現出真實密度模型的構造形態與數值大小。上述測試結果表明，在缺失低頻信息的情況下，本文方法仍然可以重構宏觀構造模型，然后以此反演結果作為初始模型，進行常規時間—頻率域反演，就可避免陷入局部極小，得到精確的反演結果。

圖17 缺失低頻成分時的時間—拉普拉斯—頻率域縱波速度(a)、橫波速度(b)和密度(c)反演結果

圖18 以時間—拉普拉斯—頻率域反演結果為初始模型的縱波速度(a)、橫波速度(b)和密度(c)常規時間—頻率域反演結果

5 結束語

本文建立了一種基于時域加權的拉普拉斯—頻率域的彈性波全波形反演方法。該方法在結合時域和頻率域各自優勢的基礎上，通過引入拉普拉斯阻尼因子，降低了全波形反演對于低頻信息的依賴性。改進的目標函數，通過時域加權的形式消除了拉普拉斯阻尼因子對地震記錄能量的衰減影響，構建形式方便、靈活。另外在變密度情況下，該方法仍能得到較高精度的縱、橫波速度反演結果，密度反演結果由于受到其他參數的串擾，使得其精度略低于速度反演結果，因此如何改善該方法的密度反演精度應是繼續研究的方向。