基于小波域最小平方濾波的多尺度自適應全波形反演

白 璐， 韓立國， 張 盼， 胡 勇

(吉林大學 地球探測科學與技術學院，長春 130011)

?

基于小波域最小平方濾波的多尺度自適應全波形反演

白 璐， 韓立國， 張 盼， 胡 勇

(吉林大學 地球探測科學與技術學院，長春 130011)

常規的全波形反演對初始模型和低頻數據的依賴性較強，反演精確度經常會受到“跳周”現象的嚴重影響。將小波域最小平方濾波引入到全波形反演中，并利用小波變換的多尺度特性，有效地提高了反演過程的穩定性，減小反演受到“跳周”現象影響的可能性。小波域最小平方濾波具有更高的精度和更好的性能，利用該算法調整模擬波場的相位，從而改善模擬波場和觀測波場之間的相位差異，并由此構造一個新的目標函數。同時利用小波變換的多尺度特性將地震數據分解為不同頻帶數據，實現多尺度反演。數值模擬實驗結果表明，基于小波域最小平方濾波的多尺度自適應全波形反演，對初始模型和低頻數據的依賴程度降低，較常規全波形反演具有更好的穩定性，受到“跳周”現象影響的可能性大大降低。

FWI； 最小平方濾波； 小波變換； 多尺度； 目標函數

0 前言

全波形反演是一種利用疊前地震記錄的波形信息來重建地下介質的物性參數的高精度建模方法。Gauthier等[1]和Mora[2]于上世界80年代實現了二維地震資料的全波形反演。實踐證明了全波形反演具有精細地刻畫地下地質構造的能力，但同時也存在著一些重要問題(采集方式等對偏移距的限制；實際數據的頻率成分不能提供全波形反演所需的低頻信息；該方法對初始模型的依賴性較強等)。近年來，全波形反演作為勘探地球物理領域的研究熱點，發展迅速，針對傳統全波形反演存在的問題提出了很多改進方式和新的方法，方法理論和實際應用方面都取得了很大的進展，Shin等[3]提出了拉普拉斯傅立葉全波形反演，是解決反演所需低頻信息缺失問題的重要突破；Warner[4]提出了自適應全波形反演(AWI)，使反演更為穩定,有效地克服了“跳周”現象的影響；應用上，Christian等[5]對西非地區海相碳酸鹽數據進行了拉普拉斯傅立葉全波形反演，得到了較好的結果。

小波變換是一種有效的多尺度時頻分析方法。它對確定的信號有一種“集中”的能力[6-7]，并且信號和噪聲在小波域具有不同的表現特征。由于這些時頻域局部特性和多分辨分析特點，以及小波變換的低熵性、去相關性等，小波變換被廣泛地用于圖像處理和信號去噪方面[8-12]?？紤]到單一方法的不足，一些學者提出將小波變換與各種濾波算法相結合的新算法[13-17]；趙艷明等[13]提出了一種小波-維納濾波算法，證明了該組合算法比單一的小波閾值去噪或維納濾波去噪的效果更好；汪魯才等[14]提出了一種有效的基于小波變換與中值濾波相結合的干涉圖濾波算法；蔡國林等[16]提出了一種小波-維納濾波器，證明了利用小波變換的優勢和維納濾波的統計特性進行干涉圖濾波可以取得很好的效果。

為了解決全波形反演存在的問題，將小波域最小平方濾波引入到全波形反演過程中，提出了時間域多尺度自適應全波形反演。利用小波變換的多尺度特性對波場數據進行分頻，從低頻帶數據開始進行反演。同時在每個頻帶中利用小波域最小平方濾波對模擬波場數據進行濾波處理，達到調整模擬波場相位的目的。最終使得基于小波域最小平方濾波的多尺度自適應全波形反演對初始模型和低頻數據的依賴性降低，反演過程更為穩定。

1 時間域全波形反演

二維全波形反演的目標函數可以寫為式(1)。

E=‖d-dobs‖2

(1)

其中：dobs為觀測到的地震數據；d為正演模擬得到的地震數據。d可以表示參數m的函數即：

d=F(m)

(2)

其中：m為描述各種地球物理參數的矢量，如速度、密度、彈性參數等；F(·)表示依賴于m的地震波傳播過程。則目標函數對m的一階導數為：

(3)

考慮時間域離散，式(1)目標函數可改寫為：

(4)

式中：Ns、Nr、Nt分別代表震源數目、檢波點數目、時間采樣點數。

伴隨狀態法求取的梯度可表示為：

(5)

其中：v表示速度；d為正傳波場；R為拾取檢波點處波場的算子。

得到梯度后，我們可以選取適當的優化方法，確定合適的步長a，利用式(6)迭代更新初始模型(Rn為模型更新量)，即可得到全波形反演的結果。

mn+1=mn+a Rn

(6)

2 小波變換及小波域最小平方濾波

假設存在不同尺度的細節空間Wj和逼近空間Vj，

(7)

(8)

其中，任意子空間都是正交的，即Wj⊥Vj，由多分辨率分析可知:

V0= V1⊕W1=V2⊕W2⊕W1=

V3⊕W3⊕W2⊕W1=Λ

(9)

則對于一個信號f(t)，可以把它分解為細節部分W1和逼近部分V1，然后再將逼近部分V1進一步分解到下一個尺度空間，如此循環，就可以得到任意尺度上的細節部分和逼近部分。將信號f(t)向不同的尺度空間Wj、Vj投影，可得到該信號在不同尺度上的細節信號和逼近信號，即：

(10)

(11)

dj,k=[f(t),ψj,k(t)]

(12)

cj,k=[f(t),φj,k(t)]

(13)

其中：ψj,k(t)為小波函數；φj,k(t)為尺度函數；dj,k、cj,k分為細節系數和近似系數。若分解尺度為J，則有

(14)

式(14)即為離散小波變換的重構公式，式(12)、式(13)就是小波變換的分解公式。

小波域中的最小平方濾波就是先對信號進行小波變換，分解為近似系數和各個尺度上的細節系數，然后對各尺度系數分別進行濾波的過程。若將小波變換后的近似系數和細節系數cJ,k、dj,k作為濾波器的輸入信號，對應的濾波算子分別為u、fj，則實際輸出分別為：

AJ,k=cJ,k*u, Dj,k=dj,k*fj

(15)

若對應的期望輸出信號分別為yJ,k、sj,k，則相應的誤差能量為：

(16)

根據最小平方原理，當誤差能量取最小值時，即可得到最佳濾波算子u、fj。

圖1 濾波前后記錄對比圖Fig.1 Comparison of the record before and after filtering(a)濾波前；(b)濾波后

3 多尺度自適應全波形反演

考慮對不同尺度地震信號地反演，對模擬波場和觀測波場分別進行小波變換，抽取近似系數和不同尺度的細節系數，即抽取不同頻帶的地震數據。對不同尺度的模擬波場數據進行最小平方濾波來調整它的相位，對于J級多尺度小波變換，共需要進行J+1次最小平方濾波：

AJ,k=cJ,ku,Dj,k=dj,kfj,j=1,2,…,J

(17)

其中:AJ,k為最小平方濾波后的近似系數；Dj,k為濾波后的細節系數；J代表分解尺度。

圖1為濾波前后記錄的對比圖，可以看到濾波后記錄間的相位差明顯減小。

則式(1)目標函數可改寫為如式(18)、式(19)所示。

(18)

(19)

(20)

其中：v表示速度；d為正傳波場；R為拾取檢波點處波場的算子。得到梯度后，使用拋物線法確定最佳步長(公式(21))，目標函數φ在α0、α0+h、α0+2h三點處滿足φ0>φ1<φ2，確定步長后依然按照公式(6)迭代更新初始模型。

(21)

4 數值模擬實驗

4.1 Marmousi模型實驗

從Marmousi模型中抽取大小為127×384個網格的模型，在模型表面增加11個網格的水層作為真實模型(圖2(a))，網格距為25 m，真實大小為3 450 m×9 600 m。為了使常規全波形反演方法易“跳周”，反演結果對比明顯，給出一個較差的初始模型(圖2 (b))。實驗中使用混合震源方式，每個震源均為主頻相同的雷克子波，具有隨機設定的不同的激發時間，震源位置隨機分布于模型表面。

首先，用一個比較低的反演主頻，提供足夠的低頻信息進行全波形反演。將震源主頻設置為6 Hz,采樣間隔為2 ms，接收時間為8 s。檢波器排列位于模型表面，各檢波點間距為25 m，共接收384道。采用有限差分算法正演，并應用矩陣快速運算，反演采用Fletcher-Reeves形式的共軛梯度法迭代更新，每次迭代至少需要兩次正演，反演過程中應用2級小波變換。

反演結果如圖3、圖4所示。可以看到，在低頻信息充足的情況下，兩種算法都可以很好地重建速度模型，反演精確度較高，反演結果相差無幾。圖5所示為模型橫向6.5 km處，反演結果與真實模型及初始模型的縱向速度對比。從圖4(a)中可以看到，模型淺中部的反演速度與真實速度基本一致(<2 km)，但對于深部，只能重建其構造，而不能準確地恢復速度信息。這是因為模型中地震波的能量分布不均勻，淺部的能量比深部要強，且全波形反演對振幅變化很敏感，加之實驗中所用的初始模型較差，初始模型速度與真實速度值相差較大，尤其是深部高阻帶部分，這對恢復真實速度值也有很大的影響。

圖2 模擬實驗所用速度模型Fig.2 Model used in numerical simulation experiment(a)真實模型；(b)初始模型

圖3 常規FWI反演結果Fig.3 Model recovered using conventional FWI

圖4 縱向速度對比Fig.4 Longitudinal velocity comparison of the result of multi-scales adaptive FWI based on wavelet transform(a)最終結果；(b)多尺度反演不同階段

圖5(a)～圖5(c)中，展示了多尺度自適應反演的不同頻帶的結果。由圖5(a)可見，低頻帶反演結果已經大致得出模型的基本輪廓，而經過中低頻反演和全頻帶反演之后，細節信息得到改善，界面和構造越來越清晰。三個階段的反演中，低頻帶信息主要能夠恢復模型輪廓，而中低頻帶反演和全頻帶反演則對速度的正確更新有更重要的貢獻(圖4(b))。

4.2 缺少低頻情況下的對比實驗

為了得到初始模型差且缺失低頻信息下的對比結果，將反演主頻提高到10 Hz,其他參數設置均保持不變。這樣提供反演的數據最低頻率大約為5 Hz左右，不能滿足常規全波形反演對低頻信息的需求。

常規全波形反演的結果如圖6所示，由于初始模型較差，造成模擬波場與觀測波場間的差異較大，加之此時低頻信息不足，導致常規l2范數陷入局部極小值，反演結果精確度較低，“跳周”現象嚴重。相比之下，可以看到圖7(c)所示的基于小波變換的多尺度自適應全波形反演結果較好，并沒有出現嚴重的“跳周”現象，雖然相比6 Hz主頻時反演結果精確度有所降低，但依然能成功地反演出模型的主要構造。實驗過程中的反演結果同真實模型的誤差函數隨迭代次數的變換如圖8所示，常規FWI算法在大概200次迭代后，受“跳周”影響誤差函數很難再下降，而基于小波域最小平方濾波的全波形反演算法的誤差函數下降更多，反演結果模型精確度更高。

為了說明不同反演頻帶的范圍，對小波變換后得到的不同頻帶的記錄進行頻譜分析。抽取不同頻帶記錄的第181道(4.5 km處)數據分別進行頻譜分析，結果如圖9所示。低頻帶范圍在5 Hz～7 Hz,中低頻帶范圍在5 Hz～9 Hz,而全頻帶范圍大致為5 Hz～13 Hz。

此外，全部測試都是在PC上完成的，其配置為Intel(R)Core(TM)i7-4790 CUP @3.60GHz和32 GB RAM。常規FWI算法迭代一次至少需要30 s，多尺度自使用全波形迭代一次至少需要50 s。雖然小波變換以及濾波過程增加了計算負擔，效率有所降低，但考慮到算法的實用性，在反演穩定性大幅增加的情況下,這些計算量是可以接受的。

圖6 多尺度自適應全波形反演結果Fig.6 Model recovered by multi-scales adaptive FWI(a)低頻段；(b)低頻+中頻段；(c)全頻帶

圖7 常規FWI反演結果Fig.7 Model recovered using conventional FWI

圖8 目標函數曲線Fig.8 The error function curve(a)多尺度自適應全波形反演；(b)常規FWI

圖9 不同頻帶數據頻譜分析Fig.9 The spectrum of the trace extracted from data in different frequency bands(a)低頻；(b)中低頻；(c)全頻帶

5 結束語

對時間域FWI和小波變換以及小波域最小平方濾波的基礎理論進行了研究，將小波域最小平方濾波應用于時間域全波形反演中，充分利用小波變換的多尺度特性進行數據分頻處理，實現了基于小波域最小平方濾波的多尺度自適應全波形反演，并得到了較好的實驗結果。通過上述Marmousi模型數值模擬實驗結果的對比，可以得出以下結論：

1)相比于常規全波形反演，基于小波變換的多尺度自適應全波形反演對初始模型和低頻信息的依賴程度更低。當初始模型較差時，常規FWI受“跳周”現象影響嚴重，反演結果精確度較低，而多尺度自適應全波形反演依然能夠很好地重建真實模型。

2)雖然增加的小波變換和濾波過程增加了計算量，對計算效率有所影響，但與常規FWI算法相比，基于小波變換的多尺度自適應全波形反演的穩定性更強，考慮到處理實際數據存在的初始模型差、低頻缺失等問題，這里算法具有更好的實用性，效率上也是可以接受的。

綜上所述，小波域最小平方濾波對調整模擬波場的相位信息具有很好的效果，基于小波域最小平方濾波的多尺度自適應全波形反演從數據處理和多尺度策略兩方面提高了反演穩定性，有效地克服了“跳周”現象的影響，具有更好的實用性。

[1] GAUTHIER O, VIRIEUX J, TARANTOLA A.Two-dimensional nonlinear inversion of seismic waveform:mumeical results[J]. Geophysics, 1986, 51(7):1387-1403.

[2] MORA P. Nonlinear two-dimensional elastic inversion of multi-offset seismic data[J]. Geophysics, 1987,52(9):1211-1228.

[3] SHIN C,Y.HO CHA.Waveform inversion in the Laplace-Fourier domain[J].Geophysical Journal International,2009,177(3),1067-1079.

[4] WARNER,M.GUASCH,L.Adaptive waveform inversion - FWI without cycle skipping:Theory[R].EAGE Extended Abstracts,Amsterdam,2014.

[5] CHRISTIAN A.,RIVERA*,BERTRAND DUQUET.Laplace-fourier FWI as an alternative model building tool for depth imaging studies:Application to marine carbonates field[C].SEG,New Orleans,2015:1054-1058.

[6] I DAUBECHIES. The wavelet transform, time-frequency localization and signal analysis[J]. IEEE Trans. on IT.1990,36(5):960-1006.

[7] BURRUS C S, GOPINATH R A, GUO H T. Introduction to wavelets and wavelet transform[M]. Upper Saddle River, Prentic Hall, 1998.

[8] MALLAT S.Multi-frequency channel decompositions of images and wavelet models[J].IEEE Trans. Speech Signal Processing,1989,37:2091-2110.

[9] MALLAT S. A theory for multiresolution signal decomposition: The wavelet representation[J]. IEEE Trans on PAMI,1989,11(7):764-693.

[10]張旭東,詹毅,馬永琴.不同信號的小波變換去噪方法[J].石油地球物理勘探,2007,42(增刊):118-123.

ZHANG X D, ZHAN Y, MA Y Q. Approaches of denoise by wavelet transform of different signals[J]. Oil Geophysical Prospecting, 2007, 42(Supplement):118-123. (In Chinese)

[11]柳建新,韓世禮,馬捷.小波分析在地震資料去噪中的應用[J].地球物理學進展,2006,21(2):541-545.

LIU J X, HAN S L, MA J. Application of wavelet analysis in seismic data denoising[J]. Progress in Geophysics, 2006,21(2):541-545. (In Chinese)

[12]張華,潘冬明,張興巖.二維小波變換在去除面波干擾中的應用[J].石油物探,2007,46(2):147-150.

ZHANG H, PAN D M, ZHANG X Y. Application of 2-D wavelet transformation in eliminating surface wave interference[J]. Geophysical Prospecting for Petroleum, 2007,46(2):147-150. (In Chinese)

[13]趙艷明,全子一.一種有效的小波- Wiener濾波去噪算法[J].北京郵電大學學報, 2004,27(4):41-45.

ZHAO Y M, QUAN Z Y. An efficient wavelet-wiener denoising Algorithm[J]. Journal of Beijing University of Posts and Telecommunications, 2004,27(4):41-45. (In Chinese)

[14]汪魯才,王耀南,毛六平.基于小波變換和中值濾波的InSAR干涉圖像濾波方法[J].測繪學報,2005,34(2):108-112.

WANG L C, WANG Y N, MAO L P, et al. An algorithm of interferometric phase filter of InSAR based on wavelet analysis and median filter algorithm[J]. Acta Geodaetica et Cartographica Sinica,2005,34(2):108-112. (In Chinese)

[15]田沛,李慶周,馬平,等.一種基于小波變換的圖像去噪新方法[J].中國圖形圖像學報, 2008, 13(3): 394-399.

TIAN P,LI Y Z,MA P,et al.A new method based on wavelet transform for image denoising[J].Journal of Image and Graphics,2008,13(3):394-399.(In Chinese)

[16]蔡國林,李永樹,劉國祥.小波-維納組合濾波算法及其在InSAR干涉圖去噪中的應用[J]. 遙感學報，2009,13(1):129-136.

CAI G L, LI S G, LIU G X. Wavelet-wiener combined filter and its application on InSAR interferogram[J]. Journal of Remote Sensing,2009, 13(1):129-136. (In Chinese)

[17]張涇周,張光磊,戴冠中.自適應算法與小波變換在心電信號濾波中的應用[J].生物醫學工程學雜志,2006, 23(5): 977-980.

ZHANG J Z, ZHANG G L, DAI G Z. The application of adaptive algorithm and wavelet transform in the filtering of ECG signal[J]. Journal of Biomedical Engineering, 2006,23(5): 977-980. (In Chinese)

[18]PLESSIX R E.A review of the adjoint-state method for computing the gradient of a functional with geophysical applications[J]. Geophysical Journal International, 2006,167(2):495-503.

Multi-scales adaptive full waveform inversion based on the wavelet domain least square filter

BAI Lu, HAN Li-guo, ZHANG Pan, HU Yong

(Jilin University, Changchun 130011, China)

Conventional full waveform inversion (FWI) has a strong dependence on the initial model and low frequency data, it is often suffers from cycle skipping when this two conditions are not met. In order to solve the problem, we introduced the wavelet domain least square filter to FWI, and took advantage of the multiscale characteristic of the wavelet transform. It can effectively avoid the influence of cycle skipping in the inversion procedure, and improves the stability of FWI. Wavelet domain least square filter has the higher accuracy than the time domain, with this can narrow the phase difference between the predicted data and observed data, and construct a new objective function, to make the inversion procedure steadily converge to the global minimum. Meanwhile, with the multi-scales characteristic of wavelet transform, the data can be divided into different frequency bands, and implement the multi-scales inversion. The result of numerical simulation experiment demonstrates that multi-scales adaptive FWI based on the wavelet transform is much less dependent on initial model and low-frequency data. The method can be immune to cycle skipping, and more robust than conventional FWI.

FWI; least square filter; wavelet transform; multi-scales; objective function

2016-03-23 改回日期：2016-04-19

國家“863”計劃重大項目課題(2014AA06A605)

白璐(1991-)，女，碩士，研究方向為地震波場模擬與波形反演，E-mail：1124739332@qq.com。

1001-1749(2016)05-0618-08

P 631.4

A

10.3969/j.issn.1001-1749.2016.05.07