低頻缺失數(shù)據(jù)的彈性波全波形反演策略?

郭啟民，何兵壽??，史才旺

(1.中國海洋大學(xué)海底科學(xué)與探測技術(shù)教育部重點實驗室,山東 青島 266100；2.青島海洋科學(xué)與技術(shù)試點國家實驗室 海洋礦產(chǎn)資源評價與探測技術(shù)功能實驗室,山東 青島 266237；3.南方科技大學(xué) 地球與空間科學(xué)學(xué)院,廣東 深圳 518000)

地震全波形反演是利用非線性最優(yōu)化方法反演地震波形記錄，獲取地下彈性參數(shù)分布的一種方法，它具有揭示復(fù)雜地質(zhì)背景下構(gòu)造細(xì)節(jié)和巖性的潛力[1]。因此，若全波形反演技術(shù)在實際生產(chǎn)中能得到成功應(yīng)用，必然能為波動方程偏移等處理環(huán)節(jié)提供更為準(zhǔn)確的彈性參數(shù)模型，提高地震數(shù)據(jù)的成像精度，其還可為后續(xù)的層位劃分、儲層預(yù)測、井位標(biāo)定等工作提供依據(jù)，降低勘探風(fēng)險。

地震全波反演的研究最早起源于Tarantola等[2]提出的基于最小二乘理論的二維聲波方程時間域全波形反演方法，之后學(xué)術(shù)界對這種方法進行了補充與完善，先后研究了聲波方程的時間域多尺度反演策略[3]，基于相位編碼技術(shù)的時間域多炮同時反演方法[4]，全波形反演的邊界存儲問題[5]及并行加速策略[6-9]等，這些成果豐富了時間域全波形反演的理論體系，推動了該技術(shù)的實用化進程。但時間域多尺度反演存在巨大的計算量問題，且當(dāng)時很難為波場延拓提供一個滿意的單頻波場。為解決上述問題，Pratt 等[10]將全波形反演理論拓展到頻率域，他應(yīng)用LU 分解算法求解阻抗矩陣，提高了多炮反演的計算效率，且由于在頻率域更易于實現(xiàn)從低頻到高頻的多尺度反演，因此頻率域全波形反演得到了廣泛關(guān)注，相關(guān)理論與算法不斷完善[11-14]，并在二維反射縱波數(shù)據(jù)的全波形反演中取得了應(yīng)用。但時間域反演和頻率域反演都存在嚴(yán)重的低頻依賴問題，當(dāng)野外記錄中缺少有效的低頻信息時，多尺度反演往往會失效，同時，三維情況下頻率域全波形反演還存在存儲消耗大和矩陣迭代計算效率低等缺陷。為解決此問題，Shin等[15-16]又基于拉普拉斯域?qū)︻l率不敏感的特性發(fā)展了拉普拉斯域全波形反演方法，降低了全波形反演對于低頻數(shù)據(jù)的依賴程度，但這種方法仍存在計算效率過低的問題，不利于工業(yè)應(yīng)用。由此，學(xué)術(shù)界又提出了基于時間域正演和拉普拉斯域反演的混合域反演算法[17-19]，這種算法將波場模擬放在時間域進行，避開了波阻抗矩陣的求逆運算，再將模擬結(jié)果變換到拉普拉斯域進行梯度尋優(yōu)計算。混合域反演算法聯(lián)合了拉普拉斯域反演和時間域正演的優(yōu)點，在波場模擬過程中利用傅里葉變換直接求解頻率域的解，不增加額外的計算量，可同時獲得多個頻率的波場信息，不需要超大內(nèi)存空間，方便實現(xiàn)多尺度算法且便于并行處理。

彈性波全波形反演技術(shù)的實現(xiàn)思路與聲波全波形反演類似，其發(fā)展歷程也與聲波全波反演技術(shù)相似。Tarantola[20]和Mora[21]首先將二維時間域全波形反演從聲波擴展到彈性波領(lǐng)域，實現(xiàn)了基于多分量數(shù)據(jù)的縱、橫波速度和密度反演。Pratt[22]將頻率域聲波波形反演理論推廣到了彈性波反演領(lǐng)域，Pyun等[23]發(fā)展了Laplace域的全波反演技術(shù)，Jun H等[24]研究了彈性波方程混合域全波反演方法。由于彈性波全波形反演和聲波全波形反演具有相似的實現(xiàn)思路與方程，因此聲波反演中的大部分技術(shù)(如多尺度算法、相位編碼技術(shù)、邊界存儲技術(shù)和并行加速技術(shù)等)都能擴展到彈性波形反演中，在實現(xiàn)過程中只需要將彈性波方程與對應(yīng)技術(shù)結(jié)合即可。同時，時間域、頻率域、拉普拉斯域以及混合域的彈性波全波反演技術(shù)的優(yōu)缺點與相同條件下的聲波全波形反演類似。

彈性波全波形反演的特殊性主要表現(xiàn)在兩方面：(1)多參數(shù)反演。彈性波全波形反演需要同時反演縱波速度和橫波速度，多種參數(shù)的耦合效應(yīng)增加了反演難度，其不適定性和計算量過大的問題也比聲波全波形反演更嚴(yán)重，同時，橫波速度的反演需要輸入數(shù)據(jù)中包含更豐富的低頻信息；(2)反演隱含了輸入的多分量地震記錄具有相同頻譜的假設(shè)，由于波場模擬時是縱波源激發(fā)，波在模型空間中以彈性波的形式傳播，此時炮點波場中各分量在時間域必然具有相同的頻譜，這就要求實際接收到的三分量地震記錄也必須具有相同的頻譜。但實際上，由于地層對縱、橫波的吸收機理不同，橫波高頻成分的衰減遠(yuǎn)較縱波嚴(yán)重，導(dǎo)致縱橫波的頻譜出現(xiàn)明顯差異，這對彈性波全波形反演中的震源設(shè)置問題造成了困擾，且目前沒有合適的解決方案。

本文主要研究彈性波全波形反演的低頻依賴問題，利用Laplace-Fourier域反演與時間域反演相結(jié)合的方法，降低低頻能量缺失對彈性波全波形反演的影響。具體步驟為：首先通過Laplace-Fourier域反演對低頻缺失的地震數(shù)據(jù)指數(shù)衰減獲得低頻分量，由此重建縱、橫波速度的中、長波長分量模型；然后將Laplace-Fourier域的反演結(jié)果作為混合域反演[24]的初始模型采用多尺度反演策略逐步恢復(fù)出縱、橫波速度模型的短波長分量，實現(xiàn)了低頻缺失時的彈性波全波形反演。

1 Laplace域波場的獲取及其特征分析

1.1 指數(shù)衰減對時間域數(shù)據(jù)低頻成分的恢復(fù)作用

現(xiàn)有一信號d(t)，對其施加一個隨時間的指數(shù)衰減項，得到一個衰減后的信號：

(1)

現(xiàn)假設(shè)d(t)為一雷克子波信號，中心頻率為5 Hz。利用衰減因子α=2，4，8對此信號進行衰減，圖1顯示了衰減前和衰減后的信號。圖中可以看出，衰減因子的作用不僅改變了信號的振幅，也改變了信號的波形。隨著衰減因子的增大，信號的第二個波谷逐漸消失。

圖1 原信號及衰減信號Fig.1 Original signal and attenuated signal

對圖1中的4個信號進行傅里葉變換得到其振幅譜(見圖2)，圖中可以看出，隨著衰減因子的增大，信號低頻成分產(chǎn)生了明顯的變化，并且衰減因子越大，低頻率成分的相對振幅越強。同時，當(dāng)衰減因子α較大時，從衰減后波場的振幅譜中能夠明顯看到零頻率處有值，而原信號的零頻率振幅幾乎為零。由此可見，時間域指數(shù)衰減使得原信號中產(chǎn)生了原信號不存在的零頻率成分和低頻成分，Ha和Shin的研究也證明了這一點[25]。Shin等[25]的研究表明：這些由于對地震波場指數(shù)衰減而新產(chǎn)生的低頻成分對于地震波全波反演來說是有效并且可用的。因此，在對缺失低頻的地震數(shù)據(jù)進行全波反演時，如果在反演前對地震數(shù)據(jù)進行指數(shù)衰減使其產(chǎn)生新的低頻成分，就可以依靠這些低頻信息反演出彈性參數(shù)的中、長波長成分，并由此建立一個比較可靠的全波反演初始模型。

圖2 圖1的振幅譜Fig.2 Amplitude spectrum of Figure 1

1.2 Laplace域波場的獲取

時間域Laplace變換由下式定義[15-16]：

(2)

由式(2)可知，Laplace變換的本質(zhì)是對指數(shù)衰減后的信號進行傅里葉變換，它是一種線性變換，可直接作用于波動方程。因此，對時間域的二階彈性波動方程[16]沿時間維度作Laplace變換可得Laplace域波動方程：

(3)

式中：x、z分別為直角坐標(biāo)系的兩個方向；U和W為Laplace域彈性波場的x分量和z分量；ρ為密度；λ和μ為拉梅常數(shù)；Fx和Fz為Laplace域震源項。

本文求解式(3)的方法如下：首先采用有限差分技術(shù)在時間域直接求解時間域彈性波方程，在求解的時間循環(huán)中加入離散Laplace變換,以最小成本提取任意頻率的Laplace域波場。根據(jù)式(2)，離散Laplace變換及其反變換表示為[26]：

(4)

式中：Δt為采樣間隔；nt為時間序列的長度；nω為頻率序列的長度；Sn=σn+iωn。

式(4)和離散傅里葉變換的區(qū)別僅在于一個時間衰減因子，時間域波動方程正演過程中在每一個時間步均利用(4)式進行變換即可得到Laplace域的波場值。Shin等[15-16]的研究表明，對于聲波方程，利用Laplace域的零頻率波場進行反演能夠獲得模型的長波長分量，同時使用零頻和低頻分量進行反演能夠獲得模型的長、中波長分量。Shin等[15-16]將這種方法稱為Laplace-Fourier域反演。由于聲波和彈性波全波形反演具有相同的數(shù)學(xué)基礎(chǔ)，這種方法很容易推廣到彈性波全波反演領(lǐng)域。

1.3 Laplace域波場特征分析

Laplace域波場的形態(tài)較常規(guī)頻率域波場抽象，本文以一個簡單模型為例來說明Laplace域波場特征。模型網(wǎng)格規(guī)模為150×200，網(wǎng)格大小10 m×10 m。炮點位于模型表面的中間位置，采用主頻8 Hz的Ricker子波，衰減常數(shù)為2。

圖3 零頻率的Laplace域波場形態(tài)Fig.3 Wave field pattern in laplacedomain with zero frequency

圖3為衰減常數(shù)2對應(yīng)的零頻率Laplace域波場，可以看出，零頻率Laplace域波場的主要特征為：波場峰值位于炮點處，隨著與炮點間距也來越遠(yuǎn)，波場振幅逐步降低，呈指數(shù)規(guī)律衰減。

接下來討論Laplace域波場與介質(zhì)速度的關(guān)系，用一個兩層速度模型來計算Laplace域波場，上層速度為1 500 m/s，下層速度分別為1 600和2 500 m/s,界面深度為100 m,縱橫波速度比為1.73。接收點位于模型表面，計算得到零頻率Laplace域地震記錄見圖4。圖中不同下層速度對應(yīng)的Laplace域零頻率波場振幅不一致，說明衰減波場的零頻率分量對速度值高低是有響應(yīng)的，它包含了介質(zhì)速度信息。

圖4 零頻率的Laplace域地震記錄(水平分量)Fig.4 Seismic records in laplacedomain with zero frequency(horizontal component)

Shin等[15-16]的研究結(jié)果表明對信號進行指數(shù)衰減會增加原本不存在的低頻成分，下面通過一個地震記錄實例來分析這些低頻成分的有效性。圖5(a)為某一單炮記錄，圖5(b)是圖5(a)高通濾波(截斷頻率3 Hz)后的結(jié)果。圖6為圖5的歸一化振幅譜和相位譜，其中歸一化振幅譜是對傅里葉變換計算出的振幅譜的每一個頻率對應(yīng)的多道振幅值作歸一化得到，由圖6可見，對于3 Hz以上的頻率段，濾波后的地震記錄保持了原地震記錄的頻譜特征，但3 Hz以下的頻譜產(chǎn)生了很大變化，此頻率段并不是原信號的有效成分，對反演不可用。

圖5某一單炮地震記錄以及高通濾波后的結(jié)果Fig.5 A single shot seismic record andits results after high pass filtering

利用衰減常數(shù)2對缺失了3 Hz以下信號的炮記錄(見圖5(b))進行衰減并求取振幅譜和相位譜(見圖7)，圖7中振幅譜是通過對圖6中歸一化振幅譜取對數(shù)得到，由圖7可知，對于濾波后的地震記錄，其衰減后的頻譜與原始記錄存在明顯的相似性，因此可以認(rèn)為衰減后波場包含了有效的低頻成分，合理利用這些人工操作產(chǎn)生的低頻成分則能夠恢復(fù)模型的大尺度信息。

2 彈性波Laplace-Fourier域反演方法

在Laplace-Fourier域反演過程中會出現(xiàn)遠(yuǎn)道信息缺失、探測深度變淺的問題。這一現(xiàn)象的產(chǎn)生原因是遠(yuǎn)道初至?xí)r間大，致使Laplace域能量衰減劇烈，波場振幅變?nèi)?。此時如果采用L2范數(shù)作為目標(biāo)函數(shù)進行反演會失去大部分遠(yuǎn)道信息，使探測深度變淺。本文中采用L1范數(shù)作為Laplace-Fourier域的反演目標(biāo)函數(shù)。

基于L1范數(shù)的目標(biāo)函數(shù)為：

(5)

式中:E為目標(biāo)函數(shù);U、D分別為正演模擬記錄和觀測記錄的Laplace變換;Ns為炮數(shù)。

此目標(biāo)函數(shù)對應(yīng)的伴隨震源表示為：

(6)

圖6 圖5的振幅譜和相位譜Fig.6 Amplitude spectrum and phase spectrum of Fig.5

圖7 對圖5進行指數(shù)衰減后的振幅譜和相位譜Fig.7 Amplitude spectrum and phase spectrum after exponential attenuation of Fig.5

上式中Δdi為第i炮的波場殘差，ri可視為歸一化后的殘差波場。為避免因Δdi數(shù)值過小而產(chǎn)生發(fā)散，通常利用一個較小的實數(shù)對ri作如下處理：

(7)

由式(7)可知，由于歸一化處理消除了殘差波場的振幅差異，因此它不會對伴隨震源的振幅產(chǎn)生影響，這實際上起到了增強Laplace域遠(yuǎn)道波場振幅的作用，增大了反演深度。

解決了上述理論問題后，采用如下步驟實現(xiàn)彈性波的Laplace-Fourier域全波形反演，獲取縱、橫波的長波長模型：①讀入初始模型，地震記錄，子波等；②進行多炮的時間域正演，過程中插入Laplace變換以獲得Laplace域波場；③對觀測地震記錄進行Laplace變換，計算Laplace域的殘差波場和伴隨震源；④將伴隨震源從Laplace域反變換到時間域；⑤將伴隨震源在時間域進行反傳，過程中插入Laplace變換以獲得Laplace域的反傳波場；⑥計算梯度；⑦計算步長，并完成模型更新；⑧檢驗是否達(dá)到迭代終止條件，若沒有則返回步驟②。

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

Laplace-Fourier域反演通過對低頻缺失的地震數(shù)據(jù)指數(shù)衰減獲得了地震記錄的低頻分量，并可由此反演出模型的中、長波長成分。以Laplace-Fourier域反演的反演結(jié)果作為初始模型即可進行彈性波混合域全波形反演。

混合域反演的特征是在時間域進行正演和反傳，在頻率域進行梯度的計算和處理。其中炮點波場的正向延拓和殘差波場的反傳采用高階有限差分法實現(xiàn)，時間域正演中的炮點波場重構(gòu)采用有效邊界存儲策略實現(xiàn)[5]。頻率域的梯度計算需要頻率域波場，因此混合域反演實現(xiàn)的關(guān)鍵在于如何在時間域正演中獲得頻率域的波場。本文在正演的時間循環(huán)中用離散傅里葉變換來提取所需的頻率域波場：

(8)

炮點波場和殘差波場的時間域延拓過程中，應(yīng)用(8)式進行變換即可獲得梯度計算所需的頻率域波場。

在混合域中，當(dāng)使用幾個頻率同步反演時，傅里葉變換可以用很小的運算代價在一次時間域正演中得到多個頻率的波場，不需要像頻率域正演那樣單獨計算每一個頻率波場。

獲得頻率域波場后，即可采用史才旺[27]的方法進行梯度的計算與模型的更新，直至得到最終反演結(jié)果。

4 數(shù)值算例

用圖8所示模型的正演記錄驗證本文算法的有效性，該模型中縱波速度模型是休斯頓大學(xué)應(yīng)用地球物理實驗室的Marmousi2模型的一部分，橫波速度模型是將縱波速度除以1.73得到的。該模型大小為8 980 m×3 480 m，網(wǎng)格大小20 m×20 m。正演所用的觀測系統(tǒng)如下：地表縱波源激發(fā)，震源為主頻5 Hz的Ricker子波，中間放炮，炮點均勻分布于地表，第一炮位于0 m處，炮間距100 m，每炮450道接收，道間距20 m，在模型兩邊排列固定不動，炮點滾動，共得到90炮兩分量記錄，各炮記錄長度均為8 s，采樣間隔2 ms。

圖8 縱橫波速度模型Fig.8 P-and S-wave velocity model

首先對兩分量正演記錄進行高通濾波，濾除3 Hz以下的頻率，利用濾波后的記錄進行常規(guī)頻率域多尺度反演。反演所用初始模型中縱、橫波速度均在深度方向上線性變換，在水平方向上無變化(見圖9)。

由于輸入數(shù)據(jù)的頻帶限制，反演過程中所采用的頻率組信息見表1。

圖9 反演初始縱、橫波速度模型Fig.9 Inversion of initial P-wave and S-wave velocity models

表1 用于反演的頻率組信息Table 1 frequency group information for inversion

每個頻帶的迭代次數(shù)設(shè)為20次，最終得到常規(guī)混合域多尺度反演結(jié)果見圖10?？梢钥吹?，反演結(jié)果與真實模型在形態(tài)上偏差較大，這說明對于缺乏低頻成分的多分量數(shù)據(jù)，常規(guī)混合域多尺度反演方法是失敗的。

同樣以圖9作為初始模型，利用本文給出的Laplace-Fourier域反演算法來對濾除了低頻成分的多分量數(shù)據(jù)進行反演。Laplace-Fourier域反演需要事先確定使用的頻率和衰減常數(shù)，本文使用的頻率為0、0.5、0.8、1.4、2 Hz，衰減常數(shù)為2。圖11顯示了Laplace-Fourier域的反演結(jié)果。容易看出，Laplace-Fourier域反演結(jié)果正確反映了原模型的大尺度信息。

圖10 混合域頻率多尺度反演結(jié)果Fig.10 Multiscale frequency inversion in mixed domain

圖11 Laplace-Fourier域反演結(jié)果Fig.11 Laplace-fourier domain inversion results

把圖11所示的Laplace-Fourier域反演結(jié)果作為初始模型，利用混合域多尺度反演方法得到的最終反演結(jié)果見圖12。Laplace-Fourier域反演和混合域多尺度聯(lián)合反演的最終結(jié)果明顯優(yōu)于單純的混合域多尺度反演方法(見圖10)，證實了本文的反演方法針對低頻缺失數(shù)據(jù)的有效性。

圖12 以圖11作為初始模型的混合域頻率多尺度反演結(jié)果Fig.12 Multiscale frequency inversion in mixed domain using Figure 11 as the initial model

5 結(jié)論

(1)彈性波全波形反演中，當(dāng)?shù)卣饠?shù)據(jù)缺失低頻成分時，低精度的初始模型極易使反演陷入局部極值，得到錯誤的反演結(jié)果。對時間域波場施加指數(shù)衰減項能夠產(chǎn)生人工的低頻成分，基于此的Laplace-Fourier域反演能夠在地震數(shù)據(jù)缺失低頻信息時恢復(fù)模型的中、長波長分量。

(2)當(dāng)多分量數(shù)據(jù)缺失低頻信息時，可以通過Laplace-Fourier域反演獲得一個相對可靠的初始模型，在此基礎(chǔ)上使用混合域頻率多尺度反演能夠更好地重建地下介質(zhì)的彈性參數(shù)，得到更為準(zhǔn)確的縱、橫波速度模型。