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

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

0 引言

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

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

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

1 全波形反演基本理論

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

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

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

AU=S

(1)

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

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

(2)

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

將式(2)中的目標函數(shù)對m求導(dǎo)，可以導(dǎo)出目標函數(shù)在頻率域關(guān)于模型參數(shù)的梯度為

(3)

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

(4)

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

(5)

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

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

(6)

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

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

(7)

(8)

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

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

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

(9)

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

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

(10)

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

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

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

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

(11)

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

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

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

(12)

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

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

2 時域指數(shù)衰減波場

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

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

(13)

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

(14)

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

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

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

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

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

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

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

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

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

3 基于L2范數(shù)的時域加權(quán)目標函數(shù)

當(dāng)?shù)卣鹩涗浿幸霑r間阻尼因子后，遠炮檢距地震道的能量會被衰減得很小。Shin等[20]為了消除阻尼因子產(chǎn)生的這一負面影響，采用了對數(shù)目標函數(shù)，即

(15)

雖然對數(shù)目標函數(shù)可以消除指數(shù)衰減的負面影響，但對數(shù)目標函數(shù)本身存在固有的缺陷，即對小功率譜波場以及波場的振幅變化十分敏感，反演穩(wěn)定性較差。為此，Jun等[39]提出了加權(quán)的L2目標函數(shù)，史才旺[27]提出了基于L1范數(shù)的對數(shù)目標函數(shù)?；谇叭说难芯?，本文使用了基于時間域加權(quán)的L2目標函數(shù)。

本文的目標函數(shù)首先還是以傳統(tǒng)的L2差值目標函數(shù)為核心，但是與前人不同的是，本方法在時間域進行加權(quán)處理，因為相較于頻率域來說，時間域?qū)τ谟涗浀牟僮鞲臃奖愫椭庇^。時間域目標函數(shù)為

(16)

(17)

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

(18)

4 模型測試

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

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

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

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

表1 反演頻率分組參數(shù)

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

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

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

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

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

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

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

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

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

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

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

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

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

表2 復(fù)頻率域反演頻率分組及參數(shù)

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

圖16 密度模型

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

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

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

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

5 結(jié)束語

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