基于小波域最小平方滤波的多尺度自适应全波形反演

白 璐， 韩立国， 张 盼， 胡 勇

(吉林大学 地球探测科学与技术学院，长春 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