CN110954945B - 一种基于动态随机震源编码的全波形反演方法 - Google Patents

Abstract

本发明公开了一种基于动态随机震源编码的全波形反演方法，属于地球物理领域，通过对震源的振幅、相位、极性、炮点位置4个方面进行随机选取，提出了一种动态随机震源编码新策略，包括以下主要步骤：随机选取实际炮与时间编码序列卷积，合成超级炮剖面；计算超级炮残差、梯度方向及最优步长；进行反演迭代，更新目标参数；判断是否达到反演终止条件。它避免了炮与炮之间的串扰噪声，在不影响反演精度的前提下大幅提高反演效率。将该新动态随机震源编码策略引入到变密度声波方程全波形反演的L‑BFGS算法求解中，有效减小了同步反演时介质密度与声波速度的串扰影响，提高了全波形反演速度。

Description

一种基于动态随机震源编码的全波形反演方法

技术领域

本发明属于地球物理领域，具体涉及一种基于动态随机震源编码的全波形反演方法。

背景技术

在常规的声波波动方程全波形反演中，通常假设介质密度为已知的常数值，仅反演声波速度。但其实速度和密度在地震数据的定量解释和油气解释中都扮演着重要角色，其中密度信息可预测储层流体饱和度。因此全波形反演中获取可靠的密度参数已成为迫切需求。随着大规模科学计算能力的提升，多参数FWI逐渐成为可能，并迅速成为FWI的一个重要研究方向。地震数据成像研究表明，多参数全波形反演存在参数间的串扰和收敛性问题的难点，另一方面，全波形反演计算量巨大，在反演过程中需要平衡计算精度和计算量之间的关系，因此在保证反演精度的基础上有效提高全波形反演的计算效率至关重要。

与传统方法对比，随机时移震源编码可有效提高叠前偏移效率，取得了良好的反演效果。Vigh等^[1]使用平面波震源编码实现了地震数据的全波形反演，Feng等^[2]提出了基于CPU/GPU和震源随机编码技术的混合域全波形反演技术，相比传统的相位编码技术有了进一步的改进，Zhang等^[3]提出了一种基于时域有限差分和相位敏感检测方法的时频混合域同步源全波形反演算法，避免了由于串扰噪声降低反演结果质量的问题，Matharu与Sacchi^[4]通过检查多参数震源编码的Hessian矩阵来研究源编码和参数权衡之间的关系，实现了多参数全波形反演中的源编码，并提供了弹性各向同性反演的应用。

变密度声波波动方程的全波形双参数反演过程中，由于声波速度与介质密度两参数的摄动对反演结果的影响敏感程度不一致，以及两参数之间存在的串扰现象，增加了多参数全波形反演的非线性程度，使反演问题的性态变差，严重影响介质密度的反演精度。

发明内容

针对现有技术的不足，本发明的目的在于提供一种拟合精度高、计算快速的变密度声波方程多参数全波形反演方法，能够在保证反演精度的情况下，有效提高反演效率的基于动态随机震源编码的全波形反演方法。

为了达到上述目的，本发明提供以下技术方案：

一种基于动态随机震源编码的全波形反演方法，包括以下步骤：

S1.将观测数据的地震时间域共剖点记录作为输入，建立初始模型并设定迭代精度；

S2.随机选取实际炮与时间编码序列卷积，合成超级炮剖面；

S3.构建反演目标函数；

S4.对目标函数进行求解；

S5.引入伴随场；

S6.计算更新步长，并更新当前模型；

S7.以步骤S6得到的模型作为新的初始模型，重复步骤S2～S6直至反演结束，得到最终的反演结果。

优选的方案，所述步骤S2中合成超级炮剖面，采用以下4种动态随机策略：

(1)横向炮点位置随机；

(2)纵向炮点相位随机，对各个震源加上一个随机长度的延时；

(3)纵向炮点极性随机，分别对各个震源乘以均匀随机的1或-1；

(4)震源振幅随机，实现动态随机震源编码策略；

其序列e_n可以表示为：

式(1)中P_random为改变编码序列极性而随机给出的整数序列；A_random为改变编码序列振幅而随机给出的实值序列；t_random为改变编码序列相位而随机给出的实值序列，以通过添加随机时延改变编码序列振幅；f为主频，t为时间。

对3者的不同设置对应了不同的编码方法，当A_random、t_random为常数，P_random为随机序列时，对应随机极性编码策略；当P_random、t_random为常数，A_random为随机序列时，对应随机振幅编码策略；当A_random、P_random为常数，t_random为随机序列时，对应随机相位编码策略；若同时对3者随机取值，对应随机极性—振幅—相位编码策略，所得到的时域编码序列具有极性随机、振幅随机、相位随机的特点。

优选的方案，所述步骤S3中构建反演目标函数，具体为采用如下步骤构建反演目标函数：

构建目标函数E(ρ,v)：

式(2)中H代表观测区域，p(x,z,t)代表对于猜测模型正演模拟得到的数据，p^obs(x,z,t)代表对于真实模型正演模拟得到的数据(观测数据)，ρ为介质密度，v为声波速度，t为时间序列。

优选的方案，所述步骤S4中对目标函数进行求解，具体为采用拉格朗日算子法及局部求导类的有限内存BFGS法求解：

应用拉格朗日算子法，目标函数变为：

式(3)中，φ、

分别为ρ、v_x、v_z的逆时参数，拉格朗日算子

满足以下方程：

采用局部求导类的有限内存BFGS法(L-BFGS)求解，对目标函数J(ρ,v)分步积分：

式(5)中：

优选的方案，所述步骤S5中引入伴随场，具体为：

令

伴随方程定义如下：

相应的梯度公式为：

优选的方案，所述步骤S6中，采用非精确线搜索方法Wolfe准则计算更新步长。

优选的方案，所述步骤S7中重复步骤S2～S6直至反演结束，采用以下准则判断是否达到反演终止条件：

若全波形反演迭代次数小于5次，则重新计算超级炮残差，返回步骤S3；每5次全波形反演迭代后，采用最大迭代次数作为迭代停止条件，若没有达到迭代停止条件，返回步骤S2。

相对于现有技术，本发明具有以下有益技术效果：

本发明提供的这种基于动态随机震源编码的全波形反演方法，提出了一种声波方程全波形双参数反演算法，并将该算法成功用于变密度声波方程的介质密度与声波速度双参数同步反演中，有效减小了同步反演时介质密度与声波速度的串扰影响；在多尺度全波形双参数反演的基础上，采用由混合震源编码策略和动态随机策略组成的动态随机震源编码策略，将不同位置的多个单一震源组合成一个混合震源，同时对时域编码序列的振幅、相位和极性随机取值，所得到的时域编码序列具有振幅随机、相位随机、极性随机的特点；从而减少了地震波场正演模拟的次数，提高了全波形反演速度，同时，能更好地压制串扰噪声，具有更高的加速比。

附图说明

图1为本发明方法的方法流程示意图。

图2为本发明方法的模型的参数分布示意图。

图3为本发明方法的模型的反演初始模型。

图4为本发明方法的模型的单炮反演结果。

图5为本发明方法的模型的超级炮反演结果。

图6为本发明方法的随机极性编码策略。

图7为本发明方法的模型的随机极性震源编码策略反演结果。

图8为本发明方法的随机时延编码策略。

图9为本发明方法的模型的随机时延震源编码策略反演结果。

图10为本发明方法的随机振幅编码策略。

图11为本发明方法的模型的随机振幅震源编码策略反演结果。

图12为本发明方法的动态随机震源编码策略。

图13为本发明方法的模型的动态随机震源编码策略反演结果。

图14为本发明方法的6种反演策略下的单道剖面对比图。

图15为本发明方法的6种反演策略下的重构误差及总耗时。

具体实施方式

下面结合具体实施例和附图对本发明进行进一步说明：

如图1所示，本发明一种基于动态随机震源编码的全波形反演方法，包括以下步骤：

步骤1：将观测数据的地震时间域共剖点记录作为输入，建立初始模型并设定迭代精度；

步骤2：随机选取实际炮与时间编码序列卷积，合成超级炮剖面；

本发明在4种随机编码策略的基础上，对随机编码策略进行了改进，即每5次全波形反演迭代后，便重新从参与反演的所有震源中随机抽出若干炮进行随机编码，使参与编码的独立震源也是随机挑选的，具有炮点位置动态随机、炮点组合动态随机的特点，4种随机策略综合动态策略便形成了动态随机震源编码策略。该动态随机震源编码策略从多个方面综合压制炮间的串扰噪声，削弱了串扰噪声对反演结果的影响，其反演目标函数可改写为：

式中e_n为时域编码序列，

表示时域卷积。由于计算第1项时可以将所有源编码入超级炮同时激发，因此只需要计算一次p，而不需要进行多次声波模拟，T_s代表原始剖面的时间长度，T_c代表序列e_n的持续时间，编码后的源具有T_s+T_c-1个时间样本，为了计算正演模拟必须更多地运行T_c个时间步长，计算效率η由下式给出：

式中，η_conv为卷积的计算效率，M为源的数目；

步骤3：构建反演目标函数；

声波全波形反演实质是利用已知的实测数据来重构地下介质中的模型参数：密度ρ和声速v的空间分布，根据正演模拟数据与实测数据之间的拟合最优，目标函数可以定义为：

式中H代表观测区域，p(x,z,t)代表对于猜测模型正演模拟得到的数据，p^obs(x,z,t)代表对于真实模型正演模拟得到的数据(观测数据)，ρ为介质密度，v为声波速度，t为时间序列。

步骤4：对目标函数进行求解；

应用拉格朗日算子法，目标函数变为：

式(3)中，φ、

分别为ρ、v_x、v_z的逆时参数，拉格朗日算子

满足以下方程：

采用局部求导类的有限内存BFGS法(L-BFGS)求解，对目标函数J(ρ,v)分步积分：

式(5)中：

步骤5：引入伴随场；

令

伴随方程定义如下：

相应的梯度公式为：

步骤6：计算更新步长，并更新当前模型；具体为采用非精确线搜索方法Wolfe准则计算更新步长；

步骤7：以步骤S6得到的模型作为新的初始模型，重复步骤S2～S6直至反演结束，得到最终的反演结果。

应用例

以下结合一个实施例对本发明方法进行进一步说明：

以1994BP国际标准模型为例，考虑不同策略情形下变密度声波方程全波形双参数反演，均选取如图3所示反演初始模型。

表1震源编码策略代号及对应策略内容

首先设置该加拿大推覆体模型网格大小为nx＝400，nz＝300，网格间隔为dx＝dz＝5m，震源为30Hz的零相位Ricker子波，在地表激发40炮，第一炮的激发点位置为0m，炮间距为50m，地表布置400个检波点，全孔径接收，记录时间总长度为1.25s，模型采样间隔为0.5ms，采用10Hz低频与30Hz高频相结合的多尺度反演。

如图4所示为策略(S1)常规单炮反演得到的结果。图中可见，模型反演结果与真实模型十分接近，总体反演效果达到预期。仅由于仅地面位置激发和接收，深部区域反演效果相对浅层稍差，而且由于模型两侧源点较模型中部少，推覆体两端的反演效果较中部差，但与真实模型都保持了一致的趋势。

图5为策略(S2)将剖面数据直接叠加为超级炮的反演结果。图中可见，由于超级炮在直接叠加时，各炮之间的信息相互影响产生的串扰噪声，导致反演剖面中褶皱界面杂乱无章，断层界面也不清晰、信噪比较差，反演效果非常不理想。

图6(a)为随机极性编码策略(S3)的时域编码序列，图6(b)为与原始剖面做卷积运算并叠加得到的超级炮剖面，图中可见，时域编码序列和超级炮剖面具有随机极性和炮点位置随机的特点，但是由于极性反转的存在，在叠加得到的超级炮剖面中存在零值，因而在特定位置会丢失部分信息，引入了新的串扰噪声，这将导致反演结果的不准确。

图7为随机极性编码策略(S3)得到的双参数的反演结果。图中可见，由于随机极性震源编码策略对炮间的串扰噪声具有一定的压制作用，相对直接叠加超级炮剖面反演精度大幅提高，噪声得到了较好的压制，异常体大体轮廓能够清晰反映，但反演剖面的中深层部位依旧存在较大的串扰噪声，反演结果亦有明显的震荡现象。

图8(a)-(b)分别为随机时移编码策略(S4)的时域编码序列及对应的超级炮剖面，图中可见，时域编码序列和超级炮剖面具有随机时移和炮点位置随机的特点，但是由于加入了随机时延，导致深层信息不够完善，且炮间的串扰噪声依然不能很好地压制，因此将导致深层反演结果较差的问题。

图9为随机时移编码策略(S4)得到的双参数的反演结果。图中可见，由于采用了随机时移震源编码策略，对炮间的串扰噪声具有了一定的压制作用，异常体大体轮廓能够清晰显示，其反演结果与随机极性编码策略类似，在深层部位仍然存在一定的串扰噪声，反演结果亦有震荡现象，且由于深层部位的信息缺失将导致相比S3策略更差的深层反演结果。

图10(a)-(b)分别为随机振幅编码策略(S5)的时域编码序列及对应的超级炮剖面，图中可见，时域编码序列和超级炮剖面具有随机振幅和炮点位置随机的特点，由于加入了随机振幅，因此炮间的串扰噪声将得到较好的改善，地下信息不存在缺失的情况，但是由于振幅随机将导致参与反演的各炮间信息的不对称，因此反演结果不能得到较大的提升。

图11为随机振幅编码策略(S5)得到的双参数的反演结果。图中可见，由于采用了随机振幅震源编码策略，对炮间的串扰噪声具有了一定的压制作用，异常体大体轮廓能够清晰显示，其反演结果与随机极性编码策略类似，在深层部位仍然存在一定的串扰噪声，反演结果亦有震荡现象，异常体的某些细节没有得到很好的突出。

图12(a)-(b)分别为动态随机编码策略(S6)的时域编码序列及对应的超级炮剖面，图中可见，时域编码序列和超级炮剖面具有随机极性、随机时移、随机振幅和炮点位置随机的特点，由于综合考虑了上述各随机方法的优点，且通过进行每5次全波形反演迭代后，动态的重新从参与反演的所有震源中随机抽出若干炮进行随机编码，使参与编码的独立震源也是随机挑选的，具有炮点位置动态随机、炮点组合动态随机特点的策略，避免了相同的独立震源导致的相同的串扰噪声作用在每次反演迭代中，从而达到进一步压制串扰噪声的目的，提高反演精度。

图13为动态随机编码策略(S6)得到的双参数的反演结果。图中可见，由于采用了动态随机震源编码策略，对炮间的串扰噪声具有了很好的压制作用，异常体轮廓能够清晰显示，深层反演结果的串扰噪声基本消除，震荡现象也基本被压制，反演精度与单炮反演结果接近。

综上所述，震源编码策略有效提高了全波形反演的计算效率，但也牺牲了部分反演精度。如何在反演过程中平衡计算精度和计算量之间的关系，在提高反演效率的基础上最大程度保证反演精度至关重要。

BP1994地形偏移国际标准模型测试表明：动态随机震源编码反演策略在给定较少的反演迭代次数时，可以达到与单炮反演策略几乎相同的精度，大幅提高了反演效率，而且该策略相比随机极性、随机时延、随机振幅震源编码反演策略可以很好地压制串扰噪声，获得了更高的反演精度，且总耗时加速比较其他震源编码反演策略有了大幅提升，综合考量，动态随机震源编码反演策略为最优策略。

本发明方法的优势如下：

为了更清晰地说明6种不同策略的反演精度，选取了模型图2中位置所示的纵向截面所得到的单道曲线，图14为6种策略反演结果与真实模型对比结果如所示。图中可见，单炮反演策略(S1)总体趋势与真实模型吻合较好，浅层偏差较少，由于观测方式为地面观测方式，得到的深层信息不，导致深层出现了细微偏差，因此单炮反演策略能够达到反演的要求，但耗时巨大，反演效率低下；超级炮反演策略(S2)与真实模型虽然具有大体相同的曲线形态，但是数值上具有很大的偏差，且反演结果具有剧烈的震荡和很大的串扰噪声，反演精度很差；随机极性震源编码反演策略(S3)与随机时延震源编码反演策略(S4)反演结果类似，与真实模型大体趋势较为吻合，相对S2策略精度有了很大的提高，但是在细节方面依然存在着较大的震荡现象，说明串扰噪声依然存在，并没有得到彻底的压制；动态随机震源编码反演策略(S5)反演结果与S1策略反演基本重合，在突变界面处比S1策略反演更为准确，曲线基本不存在剧烈的震荡现象，说明串扰噪声得到了彻底的压制，曲线形态更加平滑，且在界面分界处突变更加明显。

为了定量说明各种反演策略的效率和精度，6种不同反演策略的BP1994模型反演参数如表1所示，各种策略正演方法均采用炮并行的FDTD正演方法，最优化方法均采用L-BFGS最优化算法，反演使用的计算机配置为Intel Xeon(R)Platimum 8168CPU 2.76GHz×96，RAM Memory 503.6GiB，系统环境为Ubuntu 18.04.1TLS，Gnome3.28.2，编程环境为Anaconda3 Spyder4 Python3.6.2，反演过程均采用CPU并行策略，使用CPU核心数为50。

表2不同反演策略的BP1994模型多参数反演参数表

^*反演迭代次数为人为设置，其他迭代次数均由程序达到终止条件自行终止

重构误差的计算公式为：

式中p_True代表真实模型参数，p_Initial代表初始模型参数，p代表反演模型参数。

从表2可以得到，S1表示的单炮反演策略具有较好的重构误差，但单次迭代耗时巨大，且需要较高的迭代次数才能够达到终止条件，反演效率十分低下；S2表示的直接叠加超级炮反演策略具有较多的低频反演迭代次数，证明在反演模型大体轮廓的过程中由于串扰噪声的影响降低了反演的收敛性，反演结果不准确，具有很高的重构误差；S3表示的随机极性震源编码反演策略和S4表示的随机时延震源编码反演策略都降低了低频反演迭代次数，且S4策略在低频反演时更为高效，但都同时提高了高频反演迭代次数，证明在对低频模型进行细节反演时由于串扰噪声的影响降低了反演的收敛性，单次反演耗时加速比大约为2.14，但由于增加了迭代次数，因此总耗时加速比仅为1.32和1.46，最终反演结果重构误差略高于S1策略，且S3策略重构误差略低于S4策略；S5表示的动态随机时延震源编码反演策略给定了较少的反演迭代次数，加速比为1.99，总耗时加速比为2.44，最终反演结果重构误差显示S5策略十分接近S1策略。

为了更好地选取最优策略，定义反演策略评价函数：

式中φ_Evaluation为反演策略评价函数值，σ为介质密度重构误差与声波速度重构误差的标准差，Error_ρ、Error_v分别为介质密度和声波速度的重构误差，T_i为反演策略总耗时，T₀为参考反演策略总耗时，即传统单炮反演策略总耗时，当重构误差为零时，评价函数值为零，说明反演结果与真实模型完全一致，具有最好的反演效果；当T_i很小时，表明反演策略具有很高的加速比，此时log(T_i)/log(T₀)很小，说明反演效率很高；σ的存在保证了介质密度重构误差与声波速度重构误差相差不大。该反演评价函数既反映了密度、速度同时反演的准确性，又反映了反演效率，其中反演结果的准确性占主导因素，经计算，不考虑时间项时，当Error_ρ＝1，Error_v＝1时函数取得最大值，当Error_ρ＝0，Error_v＝0时函数取得最小值，总体而言，评价函数越小证明反演策略越优。

图15为6种反演策略的重构误差及总耗时。从图15中可以明显看出，S1表示的单炮反演策略具有很好的重构误差，但反演总耗时最长；S2表示的直接叠加超级炮反演策略具有相当大的重构误差，但反演总耗时明显减少；S3表示的随机极性震源编码反演策略、S4表示的随机时延震源编码反演策略与S5表示的随机振幅震源编码反演策略具有较低、总体而言都大于S1策略的重构误差，反演总耗时也有明显的减少，但存在着两个参数重构误差相差较大的现象；S6表示的动态随机时延震源编码反演策略具有与S1策略接近的重构误差，两个参数重构误差几乎相等，相比于其他策略具有大幅的降低。因此，针对反演精度和反演效率两个方面的综合考量，S6策略为最优的反演策略。从反演评价函数值曲线也可以看出，S1策略由于具有较好的参数重构结果，因此具有较小的反演评价函数值；S2策略虽然具有较少的反演总耗时，但重构结果很差，因此反演评价函数值很高；S3策略、S4策略与S5策略反演总耗时相差不大，但是由于重构结果的差异，因此反演评价函数值由小到大分别为S5策略、S3策略与S4策略，说明三种策略中S5策略更为有效；S6策略的反演重构结果与S1结果相近，但具有更高的反演效率，因此获得了最小的反演评价函数值，因此S6策略为最优的反演策略。

参考文献

[1]Vigh D,Starr E W,2007.3D prestack plane-wave,full-waveforminversion.Geophysics,73(5),VE135-VE144.

[2]Feng H,Liu H,Sun J.et al,2017.Hybrid domain full waveforminversion based on GPU/CPU and source random coding technique.GeophysicalProspecting for Petroleum(in Chinese),56(1):107-115.

[3]Zhang Q C,Mao W J,Zhou H,et al.,2018.Hybrid-domain simultaneous-source full waveform inversion without crosstalk noise.Geophysical JournalInternational,215,1659-1681.

[4]Matharu G,Sacchi M D,2018.Source encoding in multiparameter fullwaveform inversion.Geophysical Journal International,214,792–810.

Claims (6)

1.一种基于动态随机震源编码的全波形反演方法，其特征在于，包括以下步骤：

S1.将观测数据的地震时间域共剖点记录作为输入，建立初始模型并设定迭代精度；

S2.随机选取实际炮与时间编码序列卷积，合成超级炮剖面；

S3.构建反演目标函数；

S4.对目标函数进行求解；

S5.引入伴随场；

S6.计算更新步长，并更新当前模型；

S7.以步骤S6得到的模型作为新的初始模型，重复步骤S2～S6直至反演结束，得到最终的反演结果；

所述步骤S2中合成超级炮剖面，采用以下4种动态随机策略：

(1)横向炮点位置随机；

(2)纵向炮点相位随机，对各个震源加上一个随机长度的延时；

(3)纵向炮点极性随机，分别对各个震源乘以均匀随机的1或-1；

(4)震源振幅随机，实现动态随机震源编码策略；

其序列e_n可以表示为：

式(1)中P_random为改变编码序列极性而随机给出的整数序列；A_random为改变编码序列振幅而随机给出的实值序列；t_random为改变编码序列相位而随机给出的实值序列，以通过添加随机时延改变编码序列振幅；f为主频，t为时间；

对3者的不同设置对应了不同的编码方法，当A_random、t_random为常数，P_random为随机序列时，对应随机极性编码策略；当P_random、t_random为常数，A_random为随机序列时，对应随机振幅编码策略；当A_random、P_random为常数，t_random为随机序列时，对应随机相位编码策略；若同时对3者随机取值，对应随机极性—振幅—相位编码策略，所得到的时域编码序列具有极性随机、振幅随机、相位随机的特点。

2.根据权利要求1所述的基于动态随机震源编码的全波形反演方法，其特征在于，所述步骤S3中构建反演目标函数，具体为采用如下步骤构建反演目标函数：

构建目标函数E(ρ,v)：

式(2)中H代表观测区域，p(x,z,t)代表对于猜测模型正演模拟得到的数据，p^obs(x,z,t)代表对于真实模型正演模拟得到的数据，ρ为介质密度，v为声波速度，t为时间序列。

3.根据权利要求1所述的基于动态随机震源编码的全波形反演方法，其特征在于，所述步骤S4中对目标函数进行求解，采用拉格朗日算子法及局部求导类的有限内存BFGS法求解：

应用拉格朗日算子法，目标函数变为：

式(3)中，φ、

分别为ρ、v_x、v_z的逆时参数，拉格朗日算子

满足以下方程：

采用局部求导类的有限内存BFGS法求解，对目标函数J(ρ,v)分步积分：

式(5)中：

4.根据权利要求1所述的基于动态随机震源编码的全波形反演方法，其特征在于，所述步骤S5中引入伴随场，具体为：

令

伴随方程定义如下：

相应的梯度公式为：

5.根据权利要求1所述的基于动态随机震源编码的全波形反演方法，其特征在于，所述步骤S6中，采用非精确线搜索方法Wolfe准则计算更新步长。

6.根据权利要求1所述的基于动态随机震源编码的全波形反演方法，其特征在于，所述步骤S7中重复步骤S2～S6直至反演结束，采用以下准则判断是否达到反演终止条件：

若全波形反演迭代次数小于5次，则重新计算超级炮残差，返回步骤S4；每5次全波形反演迭代后，采用最大迭代次数作为迭代停止条件，若没有达到迭代停止条件，返回步骤S2。

Images