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

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

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

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

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

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

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

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

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

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

(1)

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

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

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

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

1.2 Laplace域波場的獲取

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

(2)

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

(3)

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

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

(4)

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

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

1.3 Laplace域波場特征分析

Laplace域波場的形態(tài)較常規(guī)頻率域波場抽象，本文以一個(gè)簡單模型為例來說明Laplace域波場特征。模型網(wǎng)格規(guī)模為150×200，網(wǎng)格大小10 m×10 m。炮點(diǎn)位于模型表面的中間位置，采用主頻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域波場的主要特征為：波場峰值位于炮點(diǎn)處，隨著與炮點(diǎn)間距也來越遠(yuǎn)，波場振幅逐步降低，呈指數(shù)規(guī)律衰減。

接下來討論Laplace域波場與介質(zhì)速度的關(guān)系，用一個(gè)兩層速度模型來計(jì)算Laplace域波場，上層速度為1 500 m/s，下層速度分別為1 600和2 500 m/s,界面深度為100 m,縱橫波速度比為1.73。接收點(diǎn)位于模型表面，計(jì)算得到零頻率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é)果表明對信號進(jìn)行指數(shù)衰減會增加原本不存在的低頻成分，下面通過一個(gè)地震記錄實(shí)例來分析這些低頻成分的有效性。圖5(a)為某一單炮記錄，圖5(b)是圖5(a)高通濾波(截?cái)囝l率3 Hz)后的結(jié)果。圖6為圖5的歸一化振幅譜和相位譜，其中歸一化振幅譜是對傅里葉變換計(jì)算出的振幅譜的每一個(gè)頻率對應(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))進(jìn)行衰減并求取振幅譜和相位譜(見圖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)酢４藭r(shí)如果采用L2范數(shù)作為目標(biāo)函數(shù)進(jìn)行反演會失去大部分遠(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進(jìn)行指數(shù)衰減后的振幅譜和相位譜Fig.7 Amplitude spectrum and phase spectrum after exponential attenuation of Fig.5

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

(7)

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

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

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

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

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

(8)

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

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

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

4 數(shù)值算例

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

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

首先對兩分量正演記錄進(jìn)行高通濾波，濾除3 Hz以下的頻率，利用濾波后的記錄進(jìn)行常規(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

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

同樣以圖9作為初始模型，利用本文給出的Laplace-Fourier域反演算法來對濾除了低頻成分的多分量數(shù)據(jù)進(jìn)行反演。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í)了本文的反演方法針對低頻缺失數(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ù)缺失低頻成分時(shí)，低精度的初始模型極易使反演陷入局部極值，得到錯誤的反演結(jié)果。對時(shí)間域波場施加指數(shù)衰減項(xiàng)能夠產(chǎn)生人工的低頻成分，基于此的Laplace-Fourier域反演能夠在地震數(shù)據(jù)缺失低頻信息時(shí)恢復(fù)模型的中、長波長分量。

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