全波形反演技术(FWI) 于20世纪80年代初问世, LAILLY于1983年开展了理论和应用研究^[1]; TARANTOLA^[2]提出采用逆时传播记录残差计算速度下降方向的新方法, 即广义最小二乘反演理论的时间域全波形反演方法, 由于受到当时计算机硬件和技术发展水平的限制, 该方法一直未得到很好的推广和应用。近几十年来随着计算机硬件技术水平的迅速发展和提高, 该技术得到迅速发展, 并逐步走向实际应用。20世纪90年代PRATT^[3]提出了频率域全波形反演方法; 2008年SHIN等^[4]提出了利用阻尼波场零频分量反演低频模型作为频率域波形反演的初始模型(Laplace域全波形反演) 的方法。

全波形反演之所以成为行业的研究热点, 是因为该技术充分利用了波的动力学和运动学信息, 在低频速度模型建立(如层析成像) 和常规中高频地震成像处理(如偏移) 之间建立了桥梁。基于目前的发展水平, 全波形反演要求地震数据具备较低的频率成分和较大的偏移距。

全波形反演可以在时间域实现^{[2, 5]}, 也可以在频率域或拉普拉斯域实现^[3-4]。数值模拟方式可以是声波假设, 也可以是弹性波假设或粘弹性假设。不同实现域和正演模拟方式面临的困难有所不同。目前实际应用最多的是时间域声波方程假设。

国内也有许多学者开展了FWI研究, 并取得了一些进展。杨勤勇等^[6]总结归纳了全波形反演研究现状及发展趋势; 王毓玮等^[7]给出了弹性波FWI的实现策略; 王华忠等^[8]对FWI的实现给出自己的一套理念, 并且针对其中的核心问题进行了分析讨论; 张猛等^[9]研究了CPU与GPU协同加速FWI方法; 曲英铭等^[10]提出了一种多尺度双变网格的时间域全波形反演方法。

全波形反演技术在实际应用中遇到了很多难题, 比如初始模型的建立、子波估计、噪声压制、理论数据和实际数据的匹配以及局部寻优条件不满足等。除了上述问题, 在具体实现过程中, 还要面临计算时间和存储量问题。很多学者开始研究有效的计算和存储策略以达到降低计算时间和存储量的目的^[11-12]。FWI在海洋地震数据应用方面的主要难点除了上述常见的问题外, 还有多次波、鬼波等噪声压制问题。针对这些问题, 本文给出了海上地震数据处理流程, 并就FWI实际应用给出了一些解决问题的方法:鬼波、多次波预处理, 子波估计以及加速策略等, 海上实际地震数据的应用结果证明了所提方法的有效性和实用性。

1 数据准备及匹配

FWI反演框架中, 正演是基础, 是FWI成功的关键因素之一。必须采用先进的正演方法和足够精确的初始模型来模拟实际地震记录中出现的波组成分, 同时尽可能减少正演过程中产生的边界反射和频散等实际地震记录中没有的波组成分。目前的正演方法得到的地震记录很难和实际地震记录匹配, 所以还要对实际地震数据进行处理, 目的是消除正演算子不能模拟的波场成分。当然, 在这个过程中, 需要保护有效信号, 尤其是数据的低频成分。

做好上述工作有以下几个关键因素:建立初始模型、子波提取、原始地震数据的去噪滤波等处理及地震数据的振幅相位匹配。本文结合对FWI方法原理的研究和技术需求, 给出海上地震数据FWI处理流程, 见图 1。

FWI对初始模型要求较高, 在地震数据低频缺失的情况下, 很难获得满足FWI需要的低频(约0~7Hz) 速度模型。目前业内比较认可的做法是采用层析成像方法获得初始模型, 然后再进行FWI, 本文所用初始模型由层析成像方法获得。

FWI正演过程其实是模拟野外地震数据采集的过程, 所提取的子波要尽可能接近实际震源的子波, 因此近源点的直达波、浅层强反射(如海底反射) 都可以用来提取子波, 采用该子波模拟的地震数据才能较好地与实际地震记录相匹配。由于实际地震数据经过处理多为零相位, 所以子波也需要转化为零相位子波。

地震数据在进行FWI前, 必须经过去噪、去除多次波及滤波等处理。FWI抗干扰能力较差, 所以去噪工作一定要做好, 尤其是线性干扰噪声的去除。FWI从低频开始做起, 所以要对地震数据进行滤波处理。

只有理论地震记录和实际地震记录较好匹配后才能进行FWI。地震数据匹配包括振幅匹配和相位匹配。振幅匹配的原则是使得两者的能量在一个量级内, 具有可比性; 相位匹配是使两者主要地震同相轴的相位差小于半个子波长度。

2 加速策略

计算量巨大是FWI面临的另一个重要问题, 为此, 本文采用如下加速策略。

2.1 L_BFGS算法

常规全波形反演是接近病态的非线性反演, 解决这类问题一般有两种途径:梯度类算法和牛顿类算法。牛顿类算法是利用目标函数一阶偏导的梯度算子和二阶偏导的海森矩阵进行约束, 理论上比梯度类算法的收敛速度更快。但求解大规模问题时海森矩阵的计算量和存储量大。针对牛顿类算法中海森矩阵存在的问题, 进一步发展提出了拟牛顿类算法, 其中BFGS算法是求解无约束问题拟牛顿算法的最有效方法之一, 它不需要直接计算海森矩阵, 而是用另一个容易计算的近似矩阵来代替。而L_BFGS算法是BFGS算法的进一步优化。在迭代过程中, 采用保存部分梯度算子和模型修正量的信息来求取伪海森矩阵, 避免海森矩阵的直接存储和精确求解, 从而复杂度大大降低, 适合大规模的非线性迭代求解。

2.2 有效边界存储策略

FWI在实际应用中必须根据计算条件来平衡计算时间和内存开销, 计算时间过长和内存开销过大都不能高效完成FWI。时间域FWI的导数计算必须利用全部的正演波场值, 对大数据来说, 这些波场值不能存储在内存中。所以有学者提出有效边界存储策略, 即在进行FWI时, 只存储最后若干时刻的波场和各个时刻逆时重构所需的有效边界内的波场, 从而推出整个时间序列各个时刻的波场^[13]。该策略牺牲了一部分计算时间, 但可以使得内存占有量指数倍地减少。

2.3 超级炮集编码策略

时间域FWI的计算量和炮集数量成正比。常规梯度类或牛顿类FWI都需要按炮序逐个计算。超级炮集编码策略根据波的叠加性质, 将多个震源函数绑在一起, 然后求解波动方程, 所得到的超炮集记录等价于单个震源函数分别正演得到的记录的叠加。震源函数分别进行正演模拟, 生成对应的模拟合成超炮集记录。最后可用实际采集的超炮集记录和模拟合成的超炮集记录差的L2范数作为新目标函数。这样, 在用梯度法求解该目标函数时, 由于超炮集数量较小使得反演中的正演模拟次数减小, 从而提高全波形反演的效率。

实际地震数据FWI应用过程中, 应根据具体数据情况和硬件条件, 选择是否采用有效边界存储策略和超级炮集编码策略。如果内存满足FWI需求, 就没有必要采用有效边界存储策略, 因为该策略增加了一次正演计算。超级炮集编码的效率和一个超级炮集中的单炮数量成正比, 例如两炮组成一个超级炮集, 计算速度比单炮计算快一倍。超级炮集中单炮数量越多, 炮与炮之间的串扰噪声干扰就会越严重, 引入的噪声干扰最终影响波场模拟的精度。

3 实际数据应用

本文所用地震数据为海外某海上三维地震数据, 工区处于深水区, 地层中含有许多盐体, 速度建模是难点。该数据由拖缆采集, 最大炮检距约为8000m。针对海上地震数据多次波和鬼波发育的特征, 给出了适用于海上地震数据的FWI流程(图 1)。图 2a和图 2b分别为多次波压制前、后的炮集记录。对比图 2a和图 2b可见, 多次波能量被有效压制(图中箭头位置)。鬼波压制是海洋地震数据处理的难点及重点, 图 3a和图 3b分别为鬼波压制前、后的炮集记录。对比图 3a和图 3b可以明显看出, 自由界面的虚反射鬼波被有效剔除(箭头位置所示)。本文采用基于海底反射的方法来提取地震子波, 如图 4所示。图 4a为子波提取示意图。采用震源和检波点互换的方法提取地震子波的反传播, 求取圆形区域的波场, 经过加权平均得到所需子波。圆形区域的半径约等于炮点处海水深度的两倍。图 4b为采用该子波提取方法对试验工区实际资料提取的地震子波。图 5a和图 5b分别为由层析成像方法获得的初始速度模型和采用FWI得到的反演速度模型。对比图 5a和图 5b可以明显看出, 反演速度模型具有更高的精度, 尤其是盐丘内幕的速度呈现出更多的细节。图 6a和图 6b分别显示了采用初始速度模型和FWI反演速度模型得到的叠前深度偏移结果(采用Kirchhoff叠前深度偏移)。由于图 6b采用了精度更高的速度模型, 所以偏移成像效果更好。图 6c和图 6d分别为图 6a和图 6b中黄色框放大显示。由图 6可见, 经过FWI更新后的速度体再进行偏移, 其同相轴更为连续, 断层更为清晰。图 7a, 图 7b和图 7c分别展示了实际炮集记录、初始速度模型正演得到的炮集记录和FWI更新速度模型正演得到的炮集记录。由图 7可见, 随着逐步更新, 速度模型逐渐接近实际情况, 其正演炮集记录也更接近实际炮集记录。

4 结论

本文给出了海上实际FWI处理的流程以及实现策略。海上深水含盐区的速度建模是难点, 将层析成像与FWI相结合, 有效压制多次波等干扰是解决海上深水盐下速度建模和成像的有效手段。实际地震数据的FWI处理效果证明了本文处理流程和实现策略的有效性和实用性。

FWI在处理实际资料中存在很多问题, 但随着采集和处理手段的提升, 结合有效的FWI实现策略, FWI将在海上地震数据处理中发挥越来越重要的作用。