基于組合震源編碼的多尺度全波形反演方法

國 運(yùn) 東

(中國石化集團(tuán)中原油田分公司 物探研究院，河南 濮陽 457001)

0 引言

全波形反演方法(FWI)綜合利用地震數(shù)據(jù)的走時(shí)、相位、振幅等全波場(chǎng)信息，通過地震數(shù)據(jù)與模擬數(shù)據(jù)的差值建立目標(biāo)泛函，不斷更新地下參數(shù)模型，使得模擬地震數(shù)據(jù)越來接近觀測(cè)數(shù)據(jù)，從而達(dá)到建立不同參數(shù)場(chǎng)的目標(biāo)，是一種高精度、高分辨率的物性參數(shù)建模方法。20世紀(jì)80年代，Tarantola[1-2]首先建立了FWI基于L2范數(shù)的觀測(cè)數(shù)據(jù)和模擬數(shù)據(jù)差值泛函，奠定了波形反演的理論框架，F(xiàn)WI逐步成為地球物理工作者的研究熱點(diǎn)。但是由于波形反演是一種迭代的計(jì)算過程，計(jì)算量巨大并受當(dāng)時(shí)計(jì)算硬件的限制，導(dǎo)致FWI方法發(fā)展甚是緩慢。20世紀(jì)90年代，Song等和Pratt等將全波形反演方法從時(shí)間域推廣到到頻率域，建立了一套完整的頻域FWI理論[3-4]。近年來隨著計(jì)算機(jī)的迅速發(fā)展，特別是大型CPU集群以及GPU顯卡的發(fā)展，全波形反演方法得到更多學(xué)者的廣泛關(guān)注與發(fā)展，介質(zhì)從簡單到復(fù)雜，從聲介質(zhì)到彈性介質(zhì)[5-6]，再到各向異性[7]以及各向異性彈性介質(zhì)[8-10]；反演參數(shù)逐步增多，從單一速度參數(shù)反演到速度、密度、各向異性、黏滯系數(shù)等多參數(shù)聯(lián)合反演[11]；反演過程在不同反演域中進(jìn)行，從時(shí)間域發(fā)展到頻率域，再到Laplace域[12]以及Laplace-Fourier域[13]等；近年來針對(duì)FWI不同問題的解決方法也相繼被提出[14-17]。

全波形反演中參數(shù)模型與地震觀測(cè)數(shù)據(jù)存在強(qiáng)的非線性關(guān)系，當(dāng)?shù)卣饠?shù)據(jù)主頻較高時(shí)，其對(duì)初始模型的依賴性高，當(dāng)初始模型與真實(shí)模型差異過大時(shí)，則會(huì)產(chǎn)生周波跳躍的現(xiàn)象。20世紀(jì)80年代Claerbout等指出地震數(shù)據(jù)中的低頻成分對(duì)于背景速度場(chǎng)中的反演較為敏感[18]，Jannane等分析得到數(shù)據(jù)中的高頻成分對(duì)于速度擾動(dòng)(反射系數(shù))較為敏感[19]。Bunks等提出分頻多尺度反演方法與網(wǎng)格多尺度反演方法，指出在大網(wǎng)格、低頻尺度下反演的背景速度場(chǎng)，可以有效地解決全波形反演的周波跳躍問題[20]；Sirgue等給出了頻域多尺度FWI的頻率選取策略[21]；Boonyasiriwat等構(gòu)建了地震數(shù)據(jù)的維納濾波器[22]，在獲取準(zhǔn)確地震子波的基礎(chǔ)上，實(shí)現(xiàn)了時(shí)域FWI的多尺度方法。

波形反演過程中，需要多次迭代進(jìn)行參數(shù)更新，每次除了進(jìn)行求取梯度的波場(chǎng)正演與殘差反傳外，還需要多次正演模擬求取步長等額外的計(jì)算量，因此計(jì)算量巨大的問題也一直是限制FWI實(shí)用化的重要因素。早在1992年Rietveld等人指出應(yīng)用控制照明的方法將單炮地震記錄合成平面波編碼的地震波記錄，然后應(yīng)用傳統(tǒng)偏移方法處理，面炮成像在目標(biāo)區(qū)域取得較好的效果，同時(shí)比傳統(tǒng)逐炮偏移方法計(jì)算效率大大提高[23-25]。Romero等通過相位編碼的方法來壓制不同炮間的串?dāng)_噪聲，通過應(yīng)用幾種不同的編碼方法，獲得了較好的偏移剖面[26-27]。Vigh和Starr將多震源技術(shù)引入到三維全波形反演中來提高反演效率，但會(huì)引入很多串?dāng)_噪聲[28]。Krebs 等在全波形反演過程中使用了動(dòng)態(tài)極性編碼方法[29]，在得到與傳統(tǒng)方法相近反演結(jié)果的同時(shí)，減小了一個(gè)數(shù)量級(jí)的計(jì)算量。隨后，Tang等通過相位編碼的多震源最小二乘方法，大大提高了反演成像的效率[30]。Schuster等通過詳細(xì)分析多震源編碼的各種影響因素，根據(jù)編碼函數(shù)和相干項(xiàng)強(qiáng)度的關(guān)系，建立了多震源反演中編碼函數(shù)與信噪比的定量關(guān)系[31-32]。Huang等和黃建平等實(shí)現(xiàn)了分頻編碼的多震源LSM，加快了反演成像的收斂速度[33-35]，并且其適用于移動(dòng)觀測(cè)系統(tǒng)。

通過分頻方法可以有效地實(shí)現(xiàn)速度場(chǎng)的多尺度反演方法，克服反演過程中存在周波跳躍現(xiàn)象，提高全波形反演計(jì)算的穩(wěn)定性，并且基于維納濾波方法反演效果較好[34]，但是維納濾波需要準(zhǔn)確的地震子波。為避免在數(shù)據(jù)域?yàn)V波，胡春輝和曲英銘等學(xué)者提出了時(shí)移多尺度反演方法[36-38]，其避免了對(duì)數(shù)據(jù)進(jìn)行低通濾波的預(yù)處理步驟，通過計(jì)算梯度時(shí)使用多次時(shí)移成像，并應(yīng)用加權(quán)求和獲取低波數(shù)的更新梯度，這種時(shí)間域反演方法更加自然和巧妙，但是這種方法再求取更新梯度時(shí)需要多次時(shí)移成像。本節(jié)進(jìn)一步發(fā)展了一種組合震源的方法，首先對(duì)子波和地震數(shù)據(jù)進(jìn)行時(shí)移組合疊加，再進(jìn)行互相關(guān)梯度求取，只需要一次逆時(shí)偏移的計(jì)算量，就可以完成梯度的求取，實(shí)現(xiàn)多尺度反演。

1 方法原理

基于L2范數(shù)全波形反演理論的目標(biāo)函數(shù)為：

(1)

其中E(v)為目標(biāo)函數(shù)，p(v,sn)為用震源sn對(duì)速度模型v正演模擬的炮記錄，dn為野外觀測(cè)的地震記錄，Ns為炮記錄的總數(shù)量。

為提高反演的計(jì)算效率將編碼技術(shù)引入到反演過程，并對(duì)震源和炮記錄進(jìn)行編碼，并且本文將所有炮都編成一個(gè)超級(jí)炮，得到基于編碼的多震源波形反演的目標(biāo)函數(shù)[29]：

(2)

式中:en為編碼矩陣序列,?表示時(shí)間域褶積。需要指出的是，對(duì)于正交極性編碼序列en來源于正交矩陣e滿足或近似滿足以下條件：eeT=I，其中T表示矩陣轉(zhuǎn)置運(yùn)算，I為對(duì)角單位陣。

對(duì)子波和地震數(shù)據(jù)進(jìn)行時(shí)移組合疊加，得到新的目標(biāo)函數(shù)：

(3)

其中，τi為第i次的時(shí)移量，Cmax為組合疊加的個(gè)數(shù)，本文發(fā)現(xiàn)一般選取3就可以達(dá)到理想的反演效果。

梯度求取公式不變，只是換成組合震源的伴隨源即可：

(4)

本文對(duì)地震子波的組合震源進(jìn)行了時(shí)域和頻域的分析，說明組合震源方法在多尺度反演的有效性。假設(shè)震源子波是主頻為15 Hz的雷克子波，其波形如圖1a所示，對(duì)震源子波進(jìn)行時(shí)移組合，其中時(shí)移參數(shù)的選取至關(guān)重要，本文后面會(huì)進(jìn)行討論，發(fā)現(xiàn)半波長的時(shí)移組合可以達(dá)到理想的效果。在這先做出半波長的時(shí)移波形圖(圖1b)和頻譜圖(圖2b)，一般記波峰到波谷選取半波長(在組合時(shí)移中一般選取半波長減去幾個(gè)點(diǎn)為最佳選取組合)。通過組合時(shí)移的波形圖可以看出通過組合時(shí)移，子波波形明顯變寬。本文中僅需兩次時(shí)移組合，權(quán)值矩陣為[0.25,0.5,0.25],通過頻譜分析，發(fā)現(xiàn)組合震源子波的頻帶明顯往低頻端移動(dòng)，從15 Hz的主頻位置移動(dòng)到10 Hz主頻位置。通過分析一維目標(biāo)泛函的形態(tài)曲線(圖3)，得出組合震源子波的目標(biāo)泛函互相關(guān)函數(shù)，出現(xiàn)局部極值的寬度明顯增加。

a—主頻15 Hz的雷克子波波形;b—半波長的組合震源波形a—Rake wavelet waveform with dominant frequency of 15 Hz;b—half wavelength combined source waveform圖1 不同子波時(shí)間域波形示意Fig.1 Waveform diagram in time domain of different wavelet

a—主頻15 Hz的雷克子波頻譜;b—時(shí)移半波長的組合震源頻譜a—Rake wavelet spectrum with dominant frequency of 15 Hz;b—half wavelength combined source spectrum圖2 不同子波頻譜示意Fig.2 Spectrogram of different wavelet spectrum

圖3 不同子波互相關(guān)示意Fig.3 Cross correlation of different wavelets

從組合震源編碼的全波形反演方法流程圖4可知：相比于傳統(tǒng)的全波形反演方法，其每次反演都需要選取不同的組合參數(shù)，對(duì)炮記錄先進(jìn)行編碼形成超級(jí)炮，然后對(duì)超級(jí)炮進(jìn)行時(shí)移組合，可以有效地減少計(jì)算量；對(duì)震源子波先進(jìn)行組合得到組合震源，再用組合震源進(jìn)行編碼，形成超級(jí)炮正演，并且在不同的迭代過程中選擇動(dòng)態(tài)編碼方法，壓制部分串?dāng)_。并且相對(duì)于傳統(tǒng)的時(shí)移全波形反演，再進(jìn)行互相關(guān)梯度求取，只需要一次逆時(shí)偏移的計(jì)算量移的計(jì)算量，就可以完成梯度的求取，大大減少了計(jì)算量，實(shí)現(xiàn)多尺度反演。

圖4 組合震源編碼的全波形反演流程Fig.4 The flow of combine coding FWI

2 模型測(cè)試

使用Marmousi模型(圖5a)進(jìn)行測(cè)試，來驗(yàn)證本文方法的有效性。速度模型大小為4.6 km×1.5 km，橫縱向網(wǎng)格間距均為10 m，速度從1.5～5.5 km/s變化。放炮數(shù)量為46，炮間隔100 m，為避免反演假像和檢波點(diǎn)分布在表面每個(gè)網(wǎng)格點(diǎn)上。采用一個(gè)線性梯度速度模型，從1.5 km/s增加到4 km/s(平均速度略低于真實(shí)模型)，作為反演的初始模型(圖5b)。反演過程中水層速度為1.5km/s，并設(shè)為恒定。采用交錯(cuò)網(wǎng)格一階有限差分方程進(jìn)行正演模擬和伴隨波場(chǎng)的計(jì)算，時(shí)間精度為2階，空間精度為8階。

圖5 反演用的真實(shí)速度模型(a)和初始速度模型(b)Fig.5 The waveform inversions velocity model with the real model(a) and the initial model(b)

a—Marmousi模型雷克子單炮地震記錄;b—組合震源單炮記錄;c—極性編碼超級(jí)炮記錄;d—組合震源超級(jí)炮記錄a—single shot gather generated by the Marmousi model using the Rake wavelet;b—single shot gather with the combined source;c—super shot gather with polarity coding;d—super shot gather with the combined source圖6 不同觀測(cè)地震記錄以及對(duì)應(yīng)的組合震源炮記錄Fig.6 Different shot gathers generated by the Marmousi model using the Rake wavelet and the combine gather of different shot gathers

為說明所提出的組合震源FWI方法的特征，將組合震源的震源項(xiàng)與傳統(tǒng)方法產(chǎn)生的震源炮記錄進(jìn)行比較，其中權(quán)值矩陣為[0.25,0.50,0.25]，組合時(shí)移參數(shù)選取為半波長。圖6a展示了單炮地震觀測(cè)記錄，圖6b顯示了組合震源單炮地震記錄；圖6c顯示了基于極性編碼的記錄超級(jí)炮，圖6d顯示了組合震源FWI的觀測(cè)記錄超級(jí)炮。通過傳統(tǒng)炮記錄和組合震源炮記錄可以看出不管是單炮記錄或者編碼超級(jí)炮記錄組，組合震源炮記錄的地震同相軸明顯變寬，說明其數(shù)據(jù)更加凸顯低波數(shù)成分，對(duì)大尺度構(gòu)造的反演能力相對(duì)增強(qiáng)。

使用主頻15 Hz的Ricker子波作為震源項(xiàng)的常規(guī)FWI的反演結(jié)果如圖7a所示。由于地震數(shù)據(jù)的主頻較高，在反演過程中存在周波跳躍問題，在較淺的區(qū)域中含有一些較強(qiáng)的反演噪聲。對(duì)于較深的區(qū)域，反演的速度模型中僅看出少部分結(jié)構(gòu)背景以及含有較強(qiáng)的背景噪聲，這是因?yàn)榘殡S源中高頻組分能量要強(qiáng)得多。為驗(yàn)證本文方法的多尺度反演效果，利用上面提出的反演策略進(jìn)行速度場(chǎng)構(gòu)建，組合參數(shù)從半波長到0變化，采用相同的迭代次數(shù)，反演結(jié)果如圖7b所示。由于初始速度場(chǎng)采用的是線性梯度模型，其與真實(shí)模型存在顯著差異，因此組合震源FWI反演結(jié)果其深部結(jié)構(gòu)的反演效果也不夠好。但是，與常規(guī)FWI的反演結(jié)果相比，可以看出在很大程度上克服了周波跳躍問題，較淺部分(深度小于0.5 km，水平位置：2～4 km)的反演結(jié)果與真實(shí)的速度模型非常接近。通過收斂的誤差曲線(圖8)可以看出本文方法可以收斂到更低的水平。

圖7 常規(guī)FWI反演(a)和使用組合震源的FWI(b)反演結(jié)果對(duì)比Fig.7 The inversion results for with conventional FWI(a) and FWI using combine source inversion strategy(b)

圖8 不同反演方法的數(shù)據(jù)誤差曲線Fig.8 Convergence curves of different multi-source FWIs

在測(cè)試中發(fā)現(xiàn)，組合震源反演對(duì)于組合時(shí)移參數(shù)的選取十分關(guān)鍵，發(fā)現(xiàn)時(shí)移參數(shù)選取過小，不能改善的反演效果有限，當(dāng)時(shí)移參數(shù)τi=0時(shí)，組合震源的反演與常規(guī)方法一致，但是當(dāng)時(shí)移參數(shù)選取過大時(shí)，反演不僅不能取得好的效果，而且有可能會(huì)變差。本文建議時(shí)移參數(shù)選取半波長，并且越接近反演效果越好。為對(duì)比反演參數(shù)對(duì)結(jié)果的影響，選取與2中維納濾波結(jié)合組合震源反演方法一致，權(quán)值系數(shù)矩陣[c1,c2,c3]=[0.25,0.50,0.25]，測(cè)試結(jié)果如圖9所示。通過對(duì)比發(fā)現(xiàn)當(dāng)時(shí)移組合參數(shù)選取0.25半波長時(shí)，反演的結(jié)果(圖9a)比常規(guī)的維納濾波方法具有改善，但是不如半波長的反演效果(圖7b)；當(dāng)組合時(shí)移參數(shù)選取一個(gè)波長的時(shí)間長度，反演結(jié)果甚至不如常規(guī)的維納濾波方法，出現(xiàn)了周波跳躍問題，復(fù)雜構(gòu)造的反演變得不清晰，斷層的反演也不準(zhǔn)確(圖9b)。

通過分頻方法可以有效地實(shí)現(xiàn)速度場(chǎng)的多尺度反演，提高全波形反演計(jì)算的穩(wěn)定性，且基于維納濾波方法反演效果較好[34]，但是維納濾波需要準(zhǔn)確的地震子波。而組合震源方法可以與不依賴子波的方法相結(jié)合達(dá)到多尺度反演的目的。首先，假設(shè)原始子波是主頻率為15 Hz的Ricker子波，采用組合震源方法對(duì)數(shù)據(jù)進(jìn)行處理，處理完的數(shù)據(jù)直接采用不依賴子波的反演方法，不依賴子波的方法選取基于卷積的方法，可以通過將收集的數(shù)據(jù)與來自正演波場(chǎng)的參考道和來自觀測(cè)波場(chǎng)的參考道卷積從而達(dá)到消除震源子波的影響[39]，其中參考道的選取采用極性編碼加權(quán)的參考道的方法[40]。圖10展示了不依賴子波的組合震源FWI的反演結(jié)果，可以看出反演的速度模型可以平穩(wěn)地收斂到正確的位置。利用組合震源的多尺度反演方法，盡管分辨率相對(duì)于準(zhǔn)確子波有所降低，但是成功避免了因?yàn)橹懿ㄌS問題產(chǎn)生的反演假象。多尺度反演策略通過恢復(fù)速度模型中更多的波數(shù)分量，可以更好地應(yīng)用所有數(shù)據(jù)成分，反演更加復(fù)雜的構(gòu)造。

a—組合時(shí)移參數(shù)τ=0.25λ;b—組合時(shí)移參數(shù)τ=1.0λa—combined time shift parameters τ=0.25λ;b—combined time shift parameters τ=1.0λ圖9 不同組合時(shí)移參數(shù)反演Fig.9 The inversion results with different combined time shift parameter

圖10 不依賴子波的組合震源反演結(jié)果Fig.10 The inversion results using combine-source multi-scale source independent FWI

3 結(jié)論及討論

本文針對(duì)FWI在中低波數(shù)反演中，由于強(qiáng)非線性導(dǎo)致反演過程中存在周波跳躍問題，提出了新的多尺度反演策略，通過利用對(duì)子波和地震數(shù)據(jù)進(jìn)行時(shí)移組合疊加，發(fā)展了一種組合震源的方法，通過模型試算對(duì)比分析，得出如下幾點(diǎn)認(rèn)識(shí)：

1)通過組合震源反演方法，可以實(shí)現(xiàn)從中低波數(shù)到高波數(shù)多尺度的反演策略，其可以有效地克服周波跳躍現(xiàn)象。基于組合震源的FWI方法與傳統(tǒng)的維納濾波相結(jié)合，突出低波數(shù)的成分，使得反演結(jié)果更加穩(wěn)定；與不依賴子波的方法相結(jié)合實(shí)現(xiàn)多尺度的反演，反演結(jié)果相對(duì)準(zhǔn)確。

2)基于組合震源的方法與多震源相結(jié)合，在進(jìn)行互相關(guān)梯度求取時(shí)，只需要一次逆時(shí)偏移的計(jì)算量，就可以完成梯度的求取，與在梯度場(chǎng)的互相關(guān)疊加的方法相比，可以有效減少計(jì)算量。

3)組合震源的方法對(duì)于組合時(shí)移參數(shù)的選取比較關(guān)鍵，一般要小于波長并接近半波長，當(dāng)?shù)卣饠?shù)據(jù)的頻率成分較高時(shí)，可以采用多次組合迭代策略，其可以進(jìn)一步突出低頻的成分，反演的結(jié)果更穩(wěn)定。

雖然本文的組合震源方法可以特定地突出中低波數(shù)分量，但是如果地震數(shù)據(jù)中不存在低頻成分也無能為力。因此如果地震數(shù)據(jù)缺少低頻成分，可以采用包絡(luò)、反射波反演等其他方法先進(jìn)行背景速度更新。此外,本文采用極性編碼多震源方法，在減少計(jì)算量的同時(shí)也存在串?dāng)_噪聲影響，如果受觀測(cè)系統(tǒng)的限制，比較難合成多震源數(shù)據(jù)，下一步可以發(fā)展不同觀測(cè)系統(tǒng)下與組合震源相結(jié)合的方法，提高反演的穩(wěn)定性。