國 運 東

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

0 引言

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

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

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

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

1 方法原理

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

(1)

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

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

(2)

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

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

(3)

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

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

(4)

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

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

a—主頻15 Hz的雷克子波頻譜;b—時移半波長的組合震源頻譜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 不同子波互相關示意Fig.3 Cross correlation of different wavelets

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

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

2 模型測試

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

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

a—Marmousi模型雷克子單炮地震記錄;b—組合震源單炮記錄;c—極性編碼超級炮記錄;d—組合震源超級炮記錄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 不同觀測地震記錄以及對應的組合震源炮記錄Fig.6 Different shot gathers generated by the Marmousi model using the Rake wavelet and the combine gather of different shot gathers

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

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

圖7 常規(guī)FWI反演(a)和使用組合震源的FWI(b)反演結果對比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

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

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

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

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

3 結論及討論

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

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

2)基于組合震源的方法與多震源相結合，在進行互相關梯度求取時，只需要一次逆時偏移的計算量，就可以完成梯度的求取，與在梯度場的互相關疊加的方法相比，可以有效減少計算量。

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

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