地球物理学报  2017, Vol. 60 Issue (12): 4790-4800   PDF    
基于Student's t分布的不依赖子波最小二乘逆时偏移
李庆洋, 黄建平 , 李振春     
中国石油大学(华东)地球科学与技术学院, 青岛 266580
摘要:最小二乘逆时偏移(Least-Squares Reverse Time Migration,LSRTM)与常规偏移相比具有更高的成像分辨率、振幅保真性及均衡性等优势,是当前研究的热点之一.震源子波的估计直接影响LSRTM结果的好坏,在实际情况下考虑到震源子波的空变特性,其估计十分困难.为了消除子波对LSRTM结果的影响,本文发展了基于卷积目标泛函的不依赖子波LSRTM算法.目标泛函由观测记录卷积模拟记录的参考道以及模拟记录卷积观测记录的参考道组成,由于观测子波和模拟子波在目标泛函的两项中同时存在,从而消除了子波的影响.此外,常用的基于L2范数拟合的LSRTM算法对噪声非常敏感,尤其是当地震数据中含有异常值时,常规LSRTM无法得到满意的结果.Student's t分布相比L2范数具有更好的稳健性,本文将其推广到不依赖子波LSRTM中,提升了算法的稳健性,最后通过理论模型及实际资料试算验证了算法的有效性和对复杂模型的适应性.
关键词: 最小二乘逆时偏移      不依赖子波      Student's t分布      目标泛函      卷积     
Source-independent least-squares reverse time migration using student's t distribution
LI Qing-Yang, HUANG Jian-Ping, LI Zhen-Chun     
School of Geosciences of China University of Petroleum(East China), Qingdao 266580, China
Abstract: Compared to the conventional migration method, the least-squares reverse time migration (LSRTM) has a lot of advantages, such as higher imaging resolution, amplitude preservation and amplitude balance. It is the focus of current research. However, the source wavelet estimation for LSRTM is a very difficult task. The challenge of determining source strength, which can vary from source to source, is even greater. In this paper, we developed a source-independent LSRTM using convolved wavefields. The misfit function consists of the convolution of the observed wavefields with a reference trace from the modeled wavefields, plus the convolution of the modeled wavefields with a reference trace from the observed wavefield. In this case, the source wavelet of the observed and the modeled wavefields are equally convolved with both terms in the misfit function, and thus, the effects of the source wavelets are eliminated. In addition, the field data often contain a lot of noise. The L2 norm based LSRTM algorithm is very sensitive to noise, especially when the data contains outliers. In this case, the conventional LSRTM result is seriously contaminated by noise. Compared to L2 norm, Student's t distribution has better robustness. We extend the Student's t distribution to the SILSRTM algorithm. Theoretical models and field data processing verify the effectiveness of the algorithm and suitability for complex models.
Key words: Least-squares reverse time migration    Source-independent    Student's t distribution    Object function    Convolution    
1 引言

近年来,随着油气勘探精度的逐渐提高,地震波成像方法也在不断发展.从早期的叠后偏移发展到叠前偏移,由时间偏移发展到深度偏移,由单程波方法发展到双程波算法,且适用介质也越来越复杂,由声波到弹性波、各向异性、黏滞性等.相比于其他偏移方法,逆时偏移(Reverse Time Migration, RTM)应用双程波动方程进行波场延拓,避免了对波动方程的近似,无倾角限制,是目前较为精确的常规偏移成像方法.然而RTM仍属于常规偏移的范畴,其偏移算子是线性正演算子的共轭转置,而不是它的逆(Claerbout,1992).当地震数据采集不足或不规则、地下构造情况复杂以及波场带宽有限时,常规偏移方法只能对地下构造模糊成像,无法满足岩性油气藏勘探开发的需求.

最小二乘偏移(Least-Squares Migration, LSM)将成像看作是最小二乘意义下的反演问题,通过迭代不断拟合观测数据和模拟(反偏移)数据,可以得到振幅保真性更好、分辨率更高、偏移噪声更少的成像结果.Tarantola(1984)首先给出了最小二乘反演问题的理论框架.随后,Nemeth等(1999)提出了最小二乘Kirchhoff偏移算法,并证实了LSM在成像缺失道和不规则数据方面的优势.Dai等(2011)将正则化项引入到LSM中,加快了收敛速度、提高了计算效率.黄建平等(2013)研究了最小二乘Kirchhoff偏移在碳酸盐岩储层中的应用.刘玉金等(2013)实现了局部倾角约束最小二乘偏移方法.Kühl和Sacchi (2003)Clapp等(2005)杨其强和张叔伦(2008)Tang(2009)沈雄君和刘能超(2012)周华敏等(2014)陈云峰等(2014)分别在最小二乘单程波方面做了深入的探讨.近年来基于双程波动方程的最小二乘逆时偏移算法逐渐成为研究热点(Dai et al., 2012郭振波和李振春,2014黄建平等,2014李振春等,2014李庆洋等,2016a).

虽然LSRTM相比常规偏移具有较大的优势,但其在实际资料处理中仍然面临着诸多困难.LSRTM的影响因素主要包括:偏移速度场的精度,震源子波估计的好坏以及地震数据噪声的强弱.实际震源子波的特征(能量、振幅及时延等)难以估计,且炸药震源的可重复性差,不同炮点位置的震源特征也不尽相同,因而子波估计一直是LSRTM及全波形反演的难点.在全波形反演中,常用的解决子波估计问题的方法主要分为两种:卷积方法与反卷积方法,且这类算法大多是在频率域实现(敖瑞德等,2015).Choi和Alkhalifah(2011)Zhang等(2015a, 2016)将卷积算法推广到时间域全波形反演中,取得了较好的效果.然而,LSRTM中子波问题研究较少(Zhang et al., 2015b),Tu和Herrmann(2015)给出了一种子波估计方法,并将其应用到最小二乘多次波成像中,但实现起来较为复杂.为此,本文将卷积算法发展到LSRTM中,提出了不依赖子波的LSRTM算法(Source-independent LSRTM, SILSRTM).SILSRTM的目标泛函由观测波场卷积模拟波场的参考道以及模拟波场卷积观测波场的参考道组成,由于观测子波和模拟子波在目标泛函的两项中同时存在,因而消除了子波的影响.

同时,考虑到地震数据常含有较多噪声,常用的基于L2范数拟合的LSRTM算法对噪声非常敏感,尤其是当地震数据中含有异常值时,常规LSRTM结果被噪声严重污染.L1模对大噪声的容忍性比L2模更好,但由于L1模在零值处不可导,因而常用Huber模或者混合模替代(Brossier et al., 2010).Student′s t分布是另一种相比L2范数具有更好稳健性的方法,已被成功用于噪声数据的全波形反演中(Jeong et al., 2015).Aravkin等(2011)的研究表明Student′s t分布相比Huber和混合模,在缺失数据下具有更好的表现,且Huber和混合模的结果严重依赖参数的选取,需要大量的尝试,而Student′s t分布没有太多参数选取的困扰,因而简单实用.本文将Student′s t分布发展到不依赖子波LSRTM中,提出了基于Student′s t分布的不依赖子波LSRTM算法(TSILSRTM).本算法在地震数据含有异常值且子波未知的情况下,仍能得到较好的结果.同时考虑到LSRTM的计算量,以上算法皆被推广到多震源相位编码下,在保证效果的同时显著提高了效率.

2 方法原理 2.1 线性化波动方程

二维常密度声波方程可表示为如下形式:

(1)

其中,s为慢度场,即速度的倒数;f为震源项;p是声压场.

将慢度场的平方s2分解为背景慢度场平方s02与扰动慢度场平方Δs2的叠加,同时总波场p可理解为由背景介质产生的背景波场p0和由扰动介质产生的扰动波场ps叠加而成.考虑到背景波场与总波场都满足波动方程,将其分别代入公式(1)并相减,然后借助线性Born近似,用背景波场p0替代总波场p0+ps,可得到扰动波场ps的控制方程:

(2)

公式(2)即为线性化正演(反偏移)方程,从中可以看出,扰动波场是背景波场与慢度扰动的相互作用作为二次震源在背景介质中传播产生的场,与慢度扰动呈线性关系,具有明确的物理含义.

利用背景介质中的格林函数G,可将公式(2)中的扰动波场表示为

(3)

式中Ω为空间范围,m=-Δs2,将其定义为模型参数.

为方便后续的推导,公式(3)可写成算子的形式:

(4)

其中,为扰动波场的线性正演算子.

2.2 不依赖子波LSRTM

常规LSRTM基于观测数据与反偏移模拟数据的L2模拟合,其目标泛函为

(5)

其中,pobs为观测记录,pcal=Lm为利用2.1节的反偏移方程计算得到的模拟记录,‖ ‖2代表向量的L2范数.

地震波场可以表示为格林函数与子波的卷积,因而公式(5)可写成如下形式:

(6)

其中,G为格林函数,s为子波,*代表卷积运算,下标cal和obs分别代表模拟数据与观测数据.实际地震资料偏移成像时,偏移速度场是已知的,对应于公式(6)中的Gcal已知,因而只有当子波估计正确时才能有效拟合逼近观测数据pobs,说明常规LSRTM结果严重依赖于子波的估计好坏.

Choi和Alkhalifah(2011)指出通过交叉卷积波场可以消除子波的影响,从而无需子波估计.具体的方法为:每炮模拟数据与同一炮观测数据的某一参考道进行卷积,获得模拟卷积波场;观测数据与同一炮模拟数据的相应参考道进行卷积,获得观测卷积波场.该思想目前主要应用于全波形反演中,本文将其发展到LSRTM中,研究不依赖子波的LSRTM(SILSRTM)算法,相应的目标泛函为

(7)

其中,pobsref为从观测数据中选出的参考道,pcalref为从模拟数据中选出的参考道.通过交叉卷积使得目标泛函的两项中均包含观测子波和模拟子波,从而消除了子波对最后结果的影响.需要指出的是,参考道的选取对结果有一定的影响,敖瑞德等(2015)以及Zhang等(2016)详细讨论了参考道的选取办法,一般选取近偏移距道或者加权平均道即可.

采用梯度导引类算法(最速下降或共轭梯度)求解公式(7),需要计算目标泛函关于模型参数的梯度,详细的推导过程参见Choi和Alkhalifah(2011),在此直接给出梯度公式的形式:

(8)

(9)

其中LTL的伴随算子,rSI为互相关后的残差记录,代表互相关运算.从rSI可以看出,虽然卷积运算会引起相位的改变,但rSI中的互相关运算可以消除由卷积引起的时间延迟,从而保证算法的正确性.

公式(8)与常规LSRTM的梯度公式具有相同的形式,仅残差波场不同.梯度求取过程也就是一个逆时偏移(RTM)的过程,与常规RTM不同的是,求梯度时的反传波场为残差记录.考虑到波场存储是RTM的一大瓶颈,本文计算梯度时采用震源波场重建方法,即只存储最后一个时刻的波场和每个时刻的边界信息,反传时将其分别作为初始条件和边界条件重构震源波场,大幅度降低了存储量.

通过公式(8)计算得到梯度后,模型的更新过程为

(10)

其中,gk为第k次迭代的梯度,αk为第k次迭代的更新步长,步长可通过线性搜索等方法计算得到.

2.3 基于Student′s t分布的LSRTM

上述基于L2范数拟合的LSRTM算法对噪声非常敏感,尤其是当地震数据中含有异常值时,LSRTM结果被噪声严重干扰.L1范数相比L2范数更加稳健,但当数据误差接近于零时,应用L1范数准则求取的目标函数梯度会出现不稳定,因而常用Huber模或混合模代替.这两类目标泛函的形式如下:

(11)

(12)

其中,Δd为模拟数据与观测数据的残差,ε=α×mean(|pobs|)为控制L1范数与L2范数权重的阈值,α为0~1之间的常数.由于最后的反演结果严重依赖于α的选取,因而一般需要多次尝试,在一定程度上增加了计算量和算法的复杂性.

Student′s t分布是另一种相比L2范数具有更好稳健性的方法,已被成功用于噪声数据的全波形反演中(Jeong et al., 2015).Aravkin等(2011)的研究表明,t分布相比Huber模和混合模,在缺失数据下具有更好的表现,且Student′s t分布对参数的依赖性较弱,因而简单实用.本文将Student′s t分布发展到不依赖子波LSRTM中,提出了基于Student′s t分布的不依赖子波LSRTM算法(TSILSRTM),相应的目标泛函为

(13)

其中,为交叉卷积后的残差,s为自由度参数,σ为尺度参数.特别的当s=1时,公式(13)变成柯西分布;当s趋于无穷时,公式(13)则退化为常规的标准正态分布.大量试算表明,针对含噪地震数据,s=2,σ=1是一组较好的参数组合,适用于绝大多数情况,因而Student′s t分布方法不需要纠结于参数的选取.

公式(13)对应的梯度为

(14)

其中,rSI=pobsrefr为互相关后的残差记录.同样,公式(14)与常规LSRTM的梯度公式具有相同的形式,因而常规L2范数下的各种加速方法都可直接应用到Student′s t分布中.

2.4 随机最优化相位编码技术

LSRTM的计算量过于庞大,从而限制了其推广应用.考虑到LSRTM的计算量与炮数成线性关系,因而通过相位编码技术(Romero et al., 2000)将多个炮集组合成一个超道集,可有效地减小计算量.相比于静态编码技术,动态编码在迭代过程中可以更好地压制由相位编码引入的串扰噪声(Schuster et al., 2011).但相位编码LSRTM的收敛速度仍然太慢,Moghaddam等(2013)在多震源全波形反演中研究发现:相位编码算法的目标泛函是真实目标泛函的随机无偏估计,其梯度也是如此.通过加权平均之前的梯度,可减小梯度的随机波动,加快收敛速度(李庆洋等,2016b).

为此,本文将随机最优化思想应用到基于Student′s t分布的不依赖子波LSRTM算法中,用以进一步提高算法的计算效率.鉴于当前迭代次数下计算的梯度相比之前迭代计算的梯度更加准确,因而本文采用指数衰减加权平均,优化后的梯度表达式为

(15)

其中,k为当前迭代次数;m为加权的前期迭代次数,综合考虑效果和效率,令其等于10;α为衰减因子,取为0.5.

3 模型试算 3.1 简单平层模型

首先以一个简单的平层模型为例来验证本文算法的正确性及有效性.模型速度场如图 1a所示,对应的扰动模型如图 1b所示,网格数为201×201,纵横向网格间距均为10 m.合成观测数据所用计算参数为:在地表以100 m间隔均匀分布21炮,每炮都是201个检波点全接收,震源为主频20 Hz的雷克子波(图 2a,高斯函数的二阶导数),时间采样间隔0.5 ms,最大记录时间1 s.

图 1 速度模型(a)和反射率模型(b) Fig. 1 (a) Velocity model; (b) Slowness perturbation
图 2 真实子波(a)和错误子波(b) Fig. 2 (a) The correct source wavelet; (b) The incorrect source wavelet

为了测试不依赖子波算法的有效性,所用错误子波如图 2b所示,形态为高斯函数的一阶导数,主频为30 Hz.可以看出,错误子波与真实子波不仅振幅和相位差异巨大,且时间延迟也有较大差别.图 3a3b分别为采用真实子波(图 2a)和错误子波(图 2b)时常规LSRTM(公式(5))计算得到的第30次迭代结果,图 3c为采用错误子波时SILSRTM(公式(7))第30次迭代结果,图 3所用色标相同.对比图 3可知,常规LSRTM在子波错误时,不仅损坏了振幅的保真性,且同相轴的相位和位置都产生了偏差,对后续的解释反演会造成严重的错误指导,而SILSRTM算法即使在子波严重错误时仍能得到正确的结果,说明了本文算法的有效性和实用性.

图 3 采用正确子波时常规LSRTM结果(a); 采用错误子波时常规LSRTM结果(b); 采用错误子波时SI-LSRTM结果(c) Fig. 3 (a) Conventional LSRTM result with correct wavelet; (b) Conventional LSRTM result with incorrect wavelet; (c) SILSRTM result with incorrect wavelet

图 4为两种算法的初始残差记录,其中4a对应常规LSRTM,4b为SILSRTM.常规LSRTM的残差是原始记录道集,不论子波对错形态都不发生任何变化,而SILSRTM的残差经过交叉卷积和互相关后,考虑了所用子波的信息,因而图 4b的同相轴相比4a明显向下移动,这也正是SILSRTM在子波错误时仍能得到较好成像结果的原因所在.

图 4 初始残差记录 (a)常规LSRTM; (b) SILSRTM. Fig. 4 Initial residual shot (a) Conventional LSRTM; (b) SILSRTM.

仍然以平层模型为例,测试本文算法对脉冲噪声的适应性.图 5a为某一不含噪声单炮记录,5b为含有少量脉冲噪声的单炮记录.图 6ab分别为常规LSRTM与基于Student′s t分布的LSRTM在含脉冲噪声数据下的成像结果,图 6a的结果被噪声严重污染,而Student′s t分布相比L2模具有更高的噪声容忍度,因而6b结果显著消除了脉冲噪声的影响,且Student′s t分布相比Huber模或混合模不受参数选取的困扰,因而简单实用.

图 5 单炮记录 (a)不含噪声; (b)脉冲噪声. Fig. 5 Single shot record without noise (a) and with spike noise (b)
图 6 含噪数据成像结果 (a)常规LSRTM; (b)基于t分布的LSRTM. Fig. 6 Image result at 30th iteration with noisy data (a) Conventional LSRTM, (b) Student′s t based LSRTM.
3.2 盐丘模型

下面以SEG/EAGE Salt模型为例,来验证算法对复杂模型的适用性.盐丘模型速度场如图 7a所示,对应的扰动模型如图 7b所示.从图 7可以看出,盐丘模型中存在高速异常体,通过对盐下弱照明区域的成像分析,可用于检验成像算法的优劣.模型参数如下:横向网格点数为645,纵向网格点数为150,纵横向网格间距皆为20 m.正演模拟计算参数为:在地表以80 m间隔均匀分布160炮,每炮都是645个检波点全接收,震源为主频16 Hz的雷克子波(图 8a),时间采样间隔1 ms,最大记录时间4 s.

图 7 速度模型(a)和反射率模型(b) Fig. 7 (a) Velocity model; (b) Slowness perturbation
图 8 真实子波(a)和错误子波(b) Fig. 8 (a) The correct source wavelet; (b) The incorrect source wavelet

由于LSRTM的计算量过于庞大,即使目前的计算机技术日新月异仍然无法满足计算需求.多震源技术可有效缓解计算量问题,但会引入串扰噪声,采用随机最优化的动态相位编码技术可很好地压制串扰噪声,在大幅度降低计算量的同时可得到与常规LSRTM算法相当的结果.本文将160炮利用震源极性编码方式组合成一个超道集,使计算量相当于单炮情形,从而大大缓解了计算需求.

图 9a为采用真实子波的常规多炮RTM叠加剖面,其成像结果被严重的低频噪声覆盖,特别是在盐丘顶部.图 9b为采用真实子波的常规相位编码LSRTM计算得到的第50次的迭代结果,与RTM相比LSRTM显著压制了成像结果中的低频和高频噪声,基本消除了震源效应,使整个剖面的振幅更加均衡,且成像分辨率也得到了极大的提高,盐下成像质量得到很大的提升.为了测试在多震源相位编码下不依赖子波算法的有效性,所用错误子波如图 8b所示,主频为10 Hz的高斯函数一阶导数,同样错误子波与真实子波振幅和相位都有较大差异.图 9c为采用错误子波时常规相位编码LSRTM计算得到的第50次迭代结果,图 9d为采用错误子波时相位编码SILSRTM算法第50次迭代结果.对比图 9c9d可知,常规相位编码LSRTM在子波错误时,不仅损坏了振幅的保真性、同相轴的相位和位置,也引入了许多偏移噪声,而相位编码SILSRTM算法即使在子波严重错误时仍能得到正确的结果.在子波错误的基础上加入少量脉冲噪声,此时基于L2模的相位编码SILSRTM算法得到的第50次迭代结果如图 9e所示,而基于Student′s t分布的相位编码SILSRTM(TSILSRTM)相应结果如图 9f所示.图 9e的有效成像剖面被脉冲噪声污染严重,而图 9f的成像结果仍然可以媲美采用真实子波且不含噪声情况下常规LSRTM的结果,从而说明了本文算法能够很好地解决震源未知且地震数据含有异常噪声的情形.

图 9 采用正确子波时常规多炮RTM结果(a); 采用正确子波时常规LSRTM结果(b); 采用错误子波时常规LSRTM结果(c); 采用错误子波时SILSRTM结果(d); 含脉冲噪声且采用错误子波时SILSRTM结果(e); 含脉冲噪声且采用错误子波时TSILSRTM结果(f) Fig. 9 (a) Conventional RTM result with correct wavelet, (b) Conventional LSRTM result with correct wavelet, (c) Conventional LSRTM result with incorrect wavelet, (d) SILSRTM result with incorrect wavelet. (e) SILSRTM result with incorrect wavelet and noisy data. (f) TSILSRTM result with incorrect wavelet and noisy data

图 10为归一化后的数据残差随迭代次数变化曲线,其中黑色星线代表采用真实子波的常规LSRTM残差收敛曲线,以此作为参考;当子波错误时常规LSRTM误差一直很大,基本不收敛、且不断震荡(粉色箭头线),因而无法得到正确的成像结果;绿色圆圈线为SILSRTM在错误子波下的残差曲线,可以看出收敛很好,基本和真实子波LSRTM收敛曲线(黑色星线)重合,间接说明了SILSRTM的有效性;蓝线代表含有异常噪声时采用错误子波的SILSRTM残差曲线,由于基于L2模的目标泛函无法处理异常噪声,因而蓝线震荡异常剧烈且不收敛;TSILSRTM在含有异常噪声且子波错误时的残差曲线为红色曲线,相比蓝线收敛得非常好,无法收敛到很小值是由于观测数据中含有异常值,而TSILSRTM可以很好地消除异常噪声的影响,因而成像结果不受噪声干扰.

图 10 残差曲线 Fig. 10 The normalized data residual curve
3.3 实际资料试算

最后将本文算法应用到西部某探区实际资料中,图 11为经走时层析建模得到的偏移速度场,左侧含有高陡构造,离散模型横向10 km,最大深度6 km,纵横向网格间距都是20 m.观测道集共97炮,且炮点分布极不规则,炮记录长度为5 s,采样点数为1251,采样间隔4 ms,道间距40 m,图 13a为采集到的某一单炮记录,黑色箭头所指倾斜同相轴为高陡构造的反射波.

图 11 速度模型 Fig. 11 Velocity model
图 12 常规RTM剖面(a)和TSILSRTM第20次迭代结果(b) Fig. 12 (a) RTM result, (b) Image result obtined by TSILSRTM at 20th iteration
图 13 真实炮记录(a)和反偏移记录(b) Fig. 13 (a) True shot record, (b) Shot record obtained by demigration

采用本文算法进行成像试算,正演模拟的时间采样间隔为1 ms,时间采样点5001.由于本文算法无需估计子波,因而只需考虑波场延拓的频散及稳定性条件,给定的模拟子波为主频15 Hz的雷克子波.图 12a为不依赖子波逆时偏移成像结果,图 12b为本文TSILSRTM迭代20次成像结果.对比图 12可以清晰发现RTM能基本反映出地下构造,但振幅均衡性较差,上部能量较弱基本被深部能量掩盖,而TSILSRTM成像质量相比于RTM有了明显改善,深部和浅部的振幅更加均衡(图中椭圆),且分辨率得到明显提高,弱信号得到一定的加强.尤其是对高陡构造的刻画,如图 12中箭头所指位置,本文算法相比RTM结果有很大提升,同相轴更加连续、分辨率更高,且RTM结果中深部的高陡能量较弱,在图 12b中也得到了明显的改善.

为了进一步检验本文算法的正确性,利用TSILSRTM第20次迭代的结果做线性正演模拟(反偏移),得到的炮记录如图 13b所示,该炮对应的真实记录如图 13a所示.对比图 13可以看出,本文算法反偏移得到的炮记录与野外记录到的炮记录吻合得非常好,主要同相轴的走时相位等都基本相同,例如图中箭头所指的高陡反射信息等,从而可在一定程度上说明本算法的正确性.

图 14为归一化后的数据残差随迭代次数变化曲线,可以看出在子波未知的情况下本文算法在实际资料下也收敛得较好,进一步验证了本算法的有效性.

图 14 残差曲线 Fig. 14 The normalized data residual curve
4 结论与认识

为了更好地适用于地震子波未知及含有异常噪声情况,本文发展了基于Student′s t分布的不依赖子波最小二乘逆时偏移算法.通过理论分析及模型试算,得到了如下几点认识:

(1) 不依赖子波LSRTM算法的目标泛函由观测记录卷积模拟记录的参考道以及模拟记录卷积观测记录的参考道组成,将真实子波与模拟子波同时引入到目标泛函的两项中,消除了子波的影响,即使模拟子波相比真实子波有较大的相位振幅及时间延迟差异,SILSRTM算法仍能得到正确的结果.

(2) Student′s t分布相比L2范数具有更好的稳健性,且与Huber模或混合模相比,Student′s t分布对参数的选取不敏感,因而简单实用.将Student′s t分布推广到不依赖子波LSRTM中,提升了算法的稳健性,能够同时解决子波未知且地震数据含有脉冲噪声的情形.

(3) 随机最优化的动态相位编码技术不仅能极大地降低计算成本,且可有效压制串扰噪声,从而将TSILSRTM的计算成本降低到与常规逆时偏移相同的水平.

目前,本文算法仍然比较初步,如下几个方向同样值得深入探索:研究如何将其推广到GPU等快速计算设备上,进一步提高计算效率;研究预条件和正则化算子等,增加算法的稳定性和加快收敛速度.

致谢

感谢中国石油大学(华东)SWPI课题组的帮助,感谢中国石油大学(北京)张庆臣博士的讨论和建议.特别感谢两位评审专家为本文完善提出的宝贵意见,感谢编辑对本文修改给予的帮助.

参考文献
Ao R D, Dong L G, Chi B X. 2015. Source-independent envelope-based FWI to build an initial model. Chinese Journal of Geophysics (in Chinese), 58(6): 1998-2010. DOI:10.6038/cjg20150615
Aravkin A, Leeuwen T V, Herrmann F J. 2011. Robust full-waveform inversion using the student's t-distribution.//81st Annual International Meeting, SEG, Expanded Abstracts, 2669-2673.
Brossier R, Operto S, Virieux J. 2010. Which data residual norm for robust elastic frequency-domain full waveform inversion. Geophysics, 75(3): R37-R46. DOI:10.1190/1.3379323
Chen Y F, Wang H Z, Ren H R, et al. 2014. Matrix representation of least square prestack migration based on linearization. Oil Geophysical Prospecting (in Chinese), 49(6): 1091-1096.
Choi Y, Alkhalifah T. 2011. Source-independent time-domain waveform inversion using convolved wavefields:Application to the encoded multisource waveform inversion. Geophysics, 76(5): R125-R134. DOI:10.1190/geo2010-0210.1
Claerbout J F. 1992. Earth Soundings Analysis:Processing Versus Inversion. Boston, MA:Blackwell Scientific Publications.
Clapp M L, Clapp R G, Biondi B L. 2005. Regularized least-squares inversion for 3-D subsalt imaging.//75th Annual International Meeting, SEG, Expanded Abstracts, 1814-1817.
Dai W, Wang X, Schuster G T. 2011. Least-squares migration of multisource data with a deblurring filter. Geophysics, 76(5): R135-R146. DOI:10.1190/geo2010-0159.1
Dai W, Fowler P, Schuster G T. 2012. Multi-source least-squares reverse time migration. Geophysical Prospecting, 60(4): 681-695. DOI:10.1111/gpr.2012.60.issue-4
Guo Z B, Li Z C. 2014. True-amplitude imaging based on least-squares reverse time migration. Oil Geophysical Prospecting (in Chinese), 49(1): 113-120.
Huang J P, Li Z C, Kong X, et al. 2013. A study of least-squares migration imaging method for fractured-type carbonate reservoir. Chinese Journal of Geophysics (in Chinese), 56(5): 1716-1725. DOI:10.6038/cjg20130529
Huang J P, Cao X L, Li Z C, et al. 2014. Least square reverse time migration in high resolution imaging of near surface. Oil Geophysical Prospecting (in Chinese), 49(1): 107-112.
Jeong W, Kang M, Kim S, et al. 2015. Full waveform inversion using student's t distribution:a numerical study for elastic waveform inversion and simultaneous-source method. Pure and Applied Geophysics, 172(6): 1491-1509. DOI:10.1007/s00024-014-1020-7
Kühl H, Sacchi M D. 2003. Least-squares wave-equation migration for AVP/AVA inversion. Geophysics, 68(1): 262-273. DOI:10.1190/1.1543212
Li Q Y, Huang J P, Li Z C, et al. 2016a. Mean-residual normalized cross-correlation least-squares reverse time migration and its application. Chinese Journal of Geophysics (in Chinese), 59(8): 3006-3015. DOI:10.6038/cjg20160823
Li Q Y, Huang J P, Li Z C, et al. 2016b. Optimized multi-source least-squares reverse time migration. Oil Geophysical Prospecting (in Chinese), 51(2): 334-341.
Li Z C, Guo Z B, Tian K. 2014. Least-squares reverse time migration in visco-acoustic medium. Chinese Journal of Geophysics (in Chinese), 57(1): 214-228. DOI:10.6038/cjg20140118
Liu Y J, Li Z C, Wu D, et al. 2013. The research on local slope constrained least-squares migration. Chinese Journal of Geophysics (in Chinese), 56(3): 1003-1011. DOI:10.6038/cjg20130328
Moghaddam P P, Keers H, Herrmann F J, et al. 2013. A new optimization approach for source-encoding full-waveform inversion. Geophysics, 78(3): R125-R132. DOI:10.1190/geo2012-0090.1
Nemeth T, Wu C J, Schuster G T. 1999. Least-squares migration of incomplete reflection data. Geophysics, 64(1): 208-221. DOI:10.1190/1.1444517
Romero L A, Ghiglia D C, Ober C C, et al. 2000. Phase encoding of shot records in prestack migration. Geophysics, 65(2): 426-436. DOI:10.1190/1.1444737
Schuster G T, Wang X, Huang Y, et al. 2011. Theory of multisource crosstalk reduction by phase-encoded statics. Geophysical Journal International, 184(3): 1289-1303. DOI:10.1111/gji.2011.184.issue-3
Shen X J, Liu N C. 2012. Split-step least-squares migration. Progress in Geophysics (in Chinese), 27(2): 761-770. DOI:10.6038/j.issn.1004-2903.2012.02.044
Tang Y X. 2009. Target-oriented wave-equation least-squares migration/inversion with phase-encoded hessian. Geophysics, 74(6): WCA95-WCA107. DOI:10.1190/1.3204768
Tarantola A. 1984. Linearized inversion of seismic reflection data. Geophysical Prospecting, 32(6): 998-1015. DOI:10.1111/gpr.1984.32.issue-6
Tu N, Herrmann F J. 2015. Fast least-squares imaging with surface-related multiples:Application to a North Sea data set. The Leading Edge, 34(7): 788-794. DOI:10.1190/tle34070788.1
Yang Q Q, Zhang S L. 2008. Least-squares fourier finite-difference migration. Progress in Geophysics (in Chinese), 23(2): 433-437.
Zhang Q C, Zhou H, Chen H M, et al. 2015a. Robust source-independent elastic full waveform inversion in the time domain.//77th EAGE Conference and Exhibition. EAGE.
Zhang Q C, Zhou H, Chen H M, et al. 2015b. Source-independent least-squares reverse time migration.//77th EAGE Conference and Exhibition. EAGE.
Zhang Q C, Zhou H, Li Q Q, et al. 2016. Robust source-independent elastic full-waveform inversion in the time domain. Geophysics, 81(2): R29-R44. DOI:10.1190/geo2015-0073.1
Zhou H M, Chen S C, Ren H R, et al. 2014. One-way wave equation least-squares migration based on illumination compensation. Chinese Journal of Geophysics (in Chinese), 57(8): 2644-2655. DOI:10.6038/cjg20140823
敖瑞德, 董良国, 迟本鑫. 2015. 不依赖子波、基于包络的FWI初始模型建立方法研究. 地球物理学报, 58(6): 1998–2010. DOI:10.6038/cjg20150615
陈云峰, 王华忠, 任浩然, 等. 2014. 线性反演最小二乘叠前偏移的矩阵形式解析. 石油地球物理勘探, 49(6): 1091–1096.
郭振波, 李振春. 2014. 最小平方逆时偏移真振幅成像. 石油地球物理勘探, 49(1): 113–120.
黄建平, 李振春, 孔雪, 等. 2013. 碳酸盐岩裂缝型储层最小二乘偏移成像方法研究. 地球物理学报, 56(5): 1716–1725. DOI:10.6038/cjg20130529
黄建平, 曹晓莉, 李振春, 等. 2014. 最小二乘逆时偏移在近地表高精度成像中的应用. 石油地球物理勘探, 49(1): 107–112.
李庆洋, 黄建平, 李振春, 等. 2016a. 去均值归一化互相关最小二乘逆时偏移及其应用. 地球物理学报, 59(8): 3006–3015. DOI:10.6038/cjg20160823
李庆洋, 黄建平, 李振春, 等. 2016b. 优化的多震源最小二乘逆时偏移. 石油地球物理勘探, 51(2): 334–341.
李振春, 郭振波, 田坤. 2014. 黏声介质最小平方逆时偏移. 地球物理学报, 57(1): 214–228. DOI:10.6038/cjg20140118
刘玉金, 李振春, 吴丹, 等. 2013. 局部倾角约束最小二乘偏移方法研究. 地球物理学报, 56(3): 1003–1011. DOI:10.6038/cjg20130328
沈雄君, 刘能超. 2012. 裂步法最小二乘偏移. 地球物理学进展, 27(2): 761–770. DOI:10.6038/j.issn.1004-2903.2012.02.044
杨其强, 张叔伦. 2008. 最小二乘傅立叶有限差分偏移. 地球物理学进展, 23(2): 433–437.
周华敏, 陈生昌, 任浩然, 等. 2014. 基于照明补偿的单程波最小二乘偏移. 地球物理学报, 57(8): 2644–2655. DOI:10.6038/cjg20140823