CN105891888B - 多域分频并行多尺度全波形反演方法 - Google Patents

Abstract

本发明涉及一种多域分频并行多尺度全波形反演方法，是用时间域经验模态分解全波形反演建立平滑初始模型，然后用分频并行Laplace‑Fourier域全波形反演来建立精度较高初始模型，最后用分频并行频域多尺度全波形反演重建地下高精度速度模型。整个过程采用动态随机编码策略，在加快反演速率的同时压制超级炮内部的串扰噪声。本发明在反演地下速度过程中用了多域联合反演、分频并行和动态随机震源编码策略。在低频缺失，初始模型不好的情况下，用时间域经验模态分解全波形反演建立平滑初始模型，缓解全波形反演跳周现象。在反演过程中采用主从式并行计算方法，极大地提高了全波形反演的计算效率，高了反演精度。

Description

多域分频并行多尺度全波形反演方法

技术领域

本发明涉及一种地震勘探的地下成像，尤其是分频并行和多域多尺度反演方法，并利用动态随机震源编码策略进行全波形反演。

背景技术：

随着石油工业的发展和勘探开发的不断深入，从构造勘探阶段逐渐走向岩性勘探阶段，其勘探难度逐渐加大。为了顺应这种要求，全波形反演迅速发展起来，并成为当今地球物理界的研究热点。20世纪80年代，塔兰托拉提出了时间域全波形反演理方法，并将目标函数对模型参数求导通过残差反传波场和正传波场互相关运算得到，从而避免了求取雅克比矩阵。20世纪90年代，普拉特将全波形反演推广到了频率域，提出只需要几个离散的频率就可以得到高精度的反演结果，而且低频到高频的反演策略可以解决陷入局部极小值的问题。在频率域全波形反演的基础上申等人发展了拉普拉斯-傅里叶域全波形反演，该方法能够在缺失低频的情况下可以获得一个高精度初始模型，然后再利用分频并行全波形反演得到最终反演结果。

全波形反演是一个基于地震全波场模拟的数据拟合过程，几乎使用了地震记录中所有有效信息，采用优化算法不断进行搜索最终找到模拟数据与实际数据拟合差最小的速度模型。多尺度全波形反演一定程度上解决了反演中的跳周问题，从低频开始逐渐增加反演的频率，可使得反演结果由原来较差的初始模型不断向真实模型靠近。但是常规频域多尺度全波形反演是一种串行反演的过程，即上一频率的反演结果作为下一频率反演的初始模型。当前大型集群可以实现多核同时工作，设想在当前的集群上进行串行全波形反演必然不能完全发挥当前计算机性能的优势。

当实际地震数据中缺失低频时，反演不能从低频开始，这会导致反演的结果很容易出现跳周的现象。在时间域吴如山等人提出利用包络全波形反演可以得到一个好的初数模型，然后再进行常规全波形反演，但是在时间域波动方程正演耗时较长，计算量较大。在频率域申等人提出利用拉普拉斯-傅里叶域全波形反演，该方法在一定程度上避免因缺失低频而导致的跳周现象。拉普拉斯-傅里叶域全波形反演其实相当于剥层法，开始将拉普拉斯衰减系数设定为一个较大的值，这样先反演浅层部分，然后不断减小拉普拉斯衰减系数，逐渐向深层进行反演。这样的反演策略可以有效避免同时反演较多的参数，由浅入深的反演策略明显减少了每一次反演的参数，即减弱了全波形反演的非线性。在实际数据中缺失5Hz以下频率地震信息，但5-7Hz信息有的情况下也不一定可靠，如果直接用5Hz信息来进行拉普拉斯-傅里叶域全波形反演，这样可能得到的反演结果并不一定是地下真实的结果，不一定符合当前地质背景。当发现当前频率信息不可靠，然后再换稍高的频率进行拉普拉斯-傅里叶域全波形反演，这样在试算不同频率的过程中，寻找最可靠频率信息可能会浪费大量的时间。

现有的频域全波形反演都是基于串行方法进行，即先从5Hz开始反演，且必须要等5Hz反演完了之后才能给下一频率6Hz进行反演，同样要等6Hz反演之后才能进行7Hz的反演。

发明内容：

本发明的目的针对上述现有技术的不足，提供一种提高全波形反演效率的多域分频并行多尺度全波形反演方法。

本发明的目是通过以下技术方案实现的：

多域分频并行多尺度全波形反演方法的核心是在计算软件2013a平台上将系统分为一个主处理器和若干从处理器，主处理器监控整个全波形反演计算阶段，由主处理器将适应度的计算分配到各个从处理器上进行，计算完成之后再由主处理器收集结果，在拉普拉斯-傅里叶域全波形反演之后由技术人员按照工区地质背景以及技术参数进行选择相对可靠的结果作为下一步分频并行全波形反演的初始模型，从而完成拉普拉斯-傅里叶域全波形反演建立初始模型的步骤。利用分频并行拉普拉斯-傅里叶域全波形反演可以同时得到多个初始模型的反演结果，有效地避免了重复操作而带来的时间浪费。同样的原理，分频并行频域全波形反演，在串行反演的基础上增加多个频率进行并行全波形反演，这样既能实现由低频到高频的多尺度反演策略，又可以充分发挥当前集群的优秀性能，显著缩短了全波形反演的计算时间。

多域分频并行多尺度全波形反演方法，包括以下步骤：

a、搭建计算软件并行工作库的安装环境，并安装计算软件并行计算工具箱；

b、对实际采集的地震数据进行预处理；

c、首先在预估速度范围建立线性递增初始模型，利用经验模态分解重构地震记录中的低频信息，其次利用时间域经验模态分解全波形反演得到平滑初始模型，再利用拉普拉斯-傅里叶域全波形反演获得高精度初始模型；

d、根据最小二乘原理构造目标函数：

其中m为模型参数，d_obs为实际采集的观测数据，d_cal在速度模型通过正演得到的计算数据，在求目标函数的梯度过程中对目标函数两端关于模型参数m求导，得到梯度表达式：

其中(A^-1)^T(d_obs-d_cal)为残差反传波场，u表示正传波场；

e、时间域地震记录进行FFT变换得到频域地震数据，查看地震数据的频谱记录，由于实际数据中缺失低频成分，所以在选择频率的时候尽量从低频开始，按照频率从低到高的顺序要求，依次挑选对应频率的地震信号。

根据要求设定多域分频并行多尺度全波形反演相关参数，包括模型大小nz×nx，网格距dx,dz，最大采样时间Nt，时间采样间隔dt，拉普拉斯衰减系数σ，反演起始频率f₀，反演频率个数nf，每个频率最大迭代次数itermax，梯度最小的数量级gtol，目标函数要求精度tol，速度反演的最大值vmax与最小值vmin；

f、利用计算软件开启系统Num个并行进程，包括一个主进程和Num个从进程，按照每个程序的需求分配给每个进程相应的内存空间；

g、给定频域地震数据，时间域经验模态分解全波形反演得到的平滑初始模型以及相关参数，进行拉普拉斯-傅里叶域全波形反演，将平滑初始模型v₀、不同反演初始频率nf₀在拉普拉斯-傅里叶域利用分频并行全波形反演得到多个不同低频反演结果，选择符合地质背景的反演结果用于下一步分频并行频域全波形反演；

h、在每个从服务器中完成拉普拉斯-傅里叶域全波形反演之后，主处理器收集各个进程反演的结果，按照技术指标及工区地质背景选择初始模型；

i、主处理器将拉普拉斯-傅里叶域全波形反演得到的高精度初始模型lapv₀重新分配到各个从处理器中，并将第一组反演频率1f₁,1f₂…1f_Num以及这几个频率所对应的地震记录信息分配到各个从处理器中，还需要分配一些相应的参数如：该频率最大迭代次数，梯度最小的数量级gtol，目标函数要求精度tol，以及速度反演的最大值vmax与最小值vmin；

j、在完成第一组频率的全波形反演之后，需要对混合震源的编码方式进行从新编码，在编码的过程中满足随机要求，主处理器将新的震源编码方式分配到各个从处理器中；

k、主处理器收集各个进程上不同频率的反演结果并取各个进程结果的平均值meanv₀，然后主处理器再将meanv₀重新分配到各个从处理器中，并将第二组反演频率2f₁,2f₂…2f_Num以及这几个频率所对应的地震记录信息分配到各个从处理器中，在进行分频并行频域全波形反演的过程中由MATLAB(科学计算软件)的主处理器按照对应的集中调度动态任务分配策略，分配以上信息到各个从处理器上，按照d步骤构造的目标函数和梯度利用超记忆梯度法在每个进程中互不干扰的进行更新迭代运算，从处理器完成第二组频率的分频并行全波形反演；

l、重复k步骤不断的增大频率直到反演频率达到设定的要求nf₁,nf₂…nf_Num，或满足精度要求则终止计算并输出反演结果；

有益效果：本发明成功地将多域分频并行多尺度等技术应用到了全波形反演中，提出的分频并行拉普拉斯-傅里叶全波形反演不仅能够减小对初始模型的依赖，缓解因低频缺失带来的跳周问题，节约了大量的时间。利用分频并行和动态随机震源编码策略显著提高了全波形反演的计算效率。

该方法是利用时间域经验模态分解全波形反演建立平滑初始模型，然后利用分频并行拉普拉斯-傅里叶域全波形反演来建立精度较高初始模型，最后利用分频并行频域全波形反演重建地下高精度速度模型。在整个计算过程中震源都采用动态随机编码策略，目的是在加快反演速率的同时压制超级炮内部的串扰噪声。与常规频域全波形反演相比，本发明反演地下速度的过程中用了多域反演和分频并行策略，通过将多个频率分频并行的进行反演，解决了以下问题：

1、利用时间域经验模态分解全波形反演方法获得平滑初始模型，该平滑初始模型优于线性递增初始模型，基本恢复出了速度模型的宏观信息，有效的缓解了因缺失低频而带来的跳周现象。

2、利用拉普拉斯-傅里叶域全波形反演可以得到一个高精度初始模型，并作为分频并行全波形反演的初始模型，一个高精度初始模型在一定程度上缓解了全波形反演的跳周现象。为全波形反演获得地下高精度速度模型奠定坚实基础。

3、在拉普拉斯-傅里叶域进行全波形反演时，起始频率对应的一些地震信息可能不可靠，需要在拉普拉斯-傅里叶域对多个不同频率的地震信息进行反演。利用分频并行策略之后，可以同时得到多个拉普拉斯-傅里叶域反演结果，然后按照相应的技术指标以及工区预先了解的地质背景选择符合实际情况的高精度初始模型。该方法有效避免了因重复操作而带来的时间上的浪费。

4、在频域全波形反演串行的基础上增加多个频率进行并行全波形反演，这样既能实现由低频到高频的多尺度反演策略，又可以充分利用当前计算集群的优秀性能，而且还加密了反演的频率，在提高效率的同时反演精度得到有效的提高。该方法显著缩短了全波形反演的计算时间。

5、动态随机震源编码策略能有效压制超级炮内部的串扰噪声，该方法在保证反演精度的同时一定程度上减少了正演的次数，有效地缩短了全波形反演的计算时间。

多域分频并行多尺度全波形反演算法比较直观，适合大规模工业开发，主从模式可以得到接近线性的加速比，能够很好的解决全波形反演计算量大的问题。

附图说明

图1多域分频并行多尺度全波形反演方法流程图。

图2串并行反演图

a串行反演，b并行反演。

v0代表初始模型，v1代表第一个频率组反演结果

v11代表第一个频率组中的第一个频率反演结果

v12代表第一个频率组中的第二个频率反演结果

v13代表第一个频率组中的第三个频率反演结果

v14代表第一个频率组中的第四个频率反演结果。

图3真实模型图。

图4初始模型图

(左)线性递增初始模型图，(右)时间域经验模态分解FWI平滑初始模型图。

图5动态随机震源编码产生的超级炮波场图

(左)随机震源组合成超级炮波场图，

(右)随机震源组合成超级炮波场加高斯噪声。

图6单炮波场快照，频域和拉普拉斯域对照图

(左)单炮频域波场快照，(右)单炮拉普拉斯-傅里叶域波场快照图。

图7分频并行拉普拉斯域初始模型建立和分频并行频域全波形反演结果图(无噪测试结果)

(左)分频并行拉普拉斯-傅里叶域初始模型建立图，

(右)分频并行频域全波形反演结果图。

图8分频并行拉普拉斯域反演模型归一化误差和分频并行频域反演模型归一化误差图(无噪测试结果)

(左)分频并行拉普拉斯-傅里叶域反演模型归一化误差图，

(右)分频并行频域反演模型归一化误差图。

图9分频并行频域反演结果单道速度对比图(无噪测试结果)

(左)120道速度对比，(右)270道速度对比图。

图10分频并行拉普拉斯域初始模型建立和分频并行频域全波形反演结果图(初始模型依赖测试结果)

(左)分频并行拉普拉斯-傅里叶域初始模型建立图，

(右)分频并行频域全波形反演结果图。

图11分频并行拉普拉斯域反演模型归一化误差和分频并行频域反演模型归一化误差图(初始模型依赖测试结果)

(左)分频并行拉普拉斯傅里叶域反演模型归一化误差图，

(右)分频并行频域反演模型归一化误差图。

图12分频并行拉普拉斯域初始模型建立和分频并行频域全波形反演结果图(抗噪能力测试结果)

(左)分频并行拉普拉斯傅里叶域初始模型建立图，

(右)分频并行频域全波形反演结果图。

图13分频并行拉普拉斯域反演模型归一化误差和分频并行频域反演模型归一化误差图(抗噪能力测试结果)

(左)分频并行拉普拉斯-傅里叶域反演模型归一化误差图，

(右)分频并行频域反演模型归一化误差图。

图14分频并行拉普拉斯域初始模型建立和分频并行频域全波形反演结果图(缺失低频测试结果)

(左)分频并行拉普拉斯-傅里叶域初始模型建立图，

(右)分频并行频域全波形反演结果图。

图15分频并行拉普拉斯域反演模型归一化误差和分频并行频域反演模型归一化误差图(缺失低频测试结果)

(左)分频并行拉普拉斯-傅里叶域反演模型归一化误差图，

(右)分频并行频域反演模型归一化误差图。

具体实施方式

下面结合附图和实例对本发明进一步的详细说明

多域分频并行多尺度全波形反演方法，包括以下步骤：

a、程序是在MATLAB(科学计算软件)2013a软件框架下编写完成，根据相应的并行计算要求，搭建MATLAB(科学计算软件)并行工作库的安装环境，并安装MATLAB(科学计算软件)并行计算工具箱。

b、实际采集的地震数据首先需进行预处理，其中预处理包括：

B1、多次波衰减、面波切除、消除交混回响和压制虚反射等不能模拟的地震波。

B2、低频保护去噪、缺失地震道补偿、保幅处理等。

B3、对地震记录进行子波估计。

c、预估速度范围并建立线性递增初始模型，利用经验模态分解重构地震记录中的低频信息，经过经验模态分解之后的地震记录中含有丰富的低频信息，而且整个地震记录的主频较低，有效地避免了因低频缺失而导致的跳周现象。然后利用时间域经验模态分解全波形反演得到一个平滑初始模型(如附图3(右))。该平滑初始模型再用拉普拉斯傅里叶域全波形反演获得一个更精确的初始模型。

经验模态分解：将单道地震波形信号中所有局部极大值点和局部极小值点识别出来，然后用三次样条曲线将所有局部极大值点连接起来构成原始波形的上包络线upenv(u(t)，同样再用三次样条曲线将所有局部极小值点连接起来构成原始波形的下包络线lowenv(u(t)，上下包络线应将原始波形包在中间。求出上下包络的平均值：

即第一次经验模态分解得到的地震波形。本专利只用第一次分解的波形即可，因为其中已经含有我们所需要的低频成分。

d、根据最小二乘原理构造目标函数：

其中m代表模型参数这里主要指地质体的速度参数，d_obs实际采集的观测数据，d_cal在速度模型通过正演模拟得到的计算数据。地震数据中的全波形反演现在仍然是一个局部优化的过程，在求目标函数的梯度过程中需要对目标函数两端关于模型参数m求导，得到梯度表达式：

波动方程Au＝s两边同时对m求导得到Jacobi(雅克比)矩阵：

其中A为阻抗矩阵，m为模型参数在本专利中代表地下速度参数，最终得到目标函数的梯度为：

其中(A^-1)^T(d_obs-d_cal)为残差反传波场(虚震源波场)，u代表正传波场，根据以上公式发现目标函数的梯度可以通过入射波场与反向残差波场做互相关运算得到，这样避免了求取Jacobi(雅克比)矩阵，很大程度上节约了计算成本。

e、时间域地震记录进行FFT变换得到频域地震数据，查看地震数据的频谱记录，由于实际数据中缺失低频成分，所以在选择频率的时候尽量从低频开始，按照频率从低到高的顺序要求，依次挑选对应频率的地震信号。

e、根据技术指标及工区要求设定多域分频并行多尺度全波形反演的相关参数，包括模型大小nz×nx，网格距dx,dz，最大采样时间Nt，时间采样间隔dt，拉普拉斯衰减系数σ，反演起始频率f₀，反演频率个数Nf，每个频率最大迭代次数，梯度最小的数量级gtol，目标函数要求精度tol，以及该工区速度反演的最大值与最小值vmax,vmin。

f、利用MATLAB(科学计算软件)开启系统Num个并行进程，包括一个主进程和Num个从进程，按照每个程序的需求分配给每个进程相应的内存空间。

g、给定频域地震数据，时间域经验模态分解全波形反演得到的初始模型以及相关参数，首先进行拉普拉斯傅里叶域全波形反演。将初始模型v₀、不同反演初始频率nf₀(选择Num较低频率，因为地震数据有些低频不可靠，在拉普拉斯傅里叶域利用分频并行反演得到多个不同低频反演结果，最后选择符合地质背景的反演结果用于下一步常规频域全波形反演)。在进行拉普拉斯傅里叶域全波形反演的过程中选择多组拉普拉斯衰减系数σ，从σ＝5每隔0.2逐渐减小一直到σ＝0为止，在每次改变σ的同时对混合震源进行从新编码，在编码的过程中满足随机要求。由MATLAB(科学计算软件)的主处理器按照相应的集中调度动态任务分配策略，分配不同的初始频率到各个从处理器上，按照d步骤构造的目标函数和梯度利用超记忆梯度优化算法在每个进程中互不干扰的进行更新迭代运算。

拉普拉斯傅里叶域全波形反演：时间域经过变换得到拉普拉斯傅里叶域，然后在拉普拉斯傅里叶进行正演模拟，求取梯度进行反演具体变换公式如下所示：

按照上式的变换公式可以得到拉普拉斯傅里叶域声波方程：

其中s＝iw+σ，σ为拉普拉斯域衰减系数，fs表示震源，v表示地下地震波是速度。

h、在每个从服务器中完成分频并行拉普拉斯傅里叶域全波形反演之后，主处理器收集各个进程反演的结果，按照相应的技术指标以及工区预先了解的地质背景选择符合实际情况的高精度初始模型。

i、主处理器将拉普拉斯傅里叶域全波形反演得到的平滑初始模型lapv₀重新分配到各个从处理器中，并将第一组反演频率1f₁,1f₂…1f_Num以及这几个频率所对应的地震记录信息分配到各个从处理器中。还需要分配一些相应的参数如：该频率最大迭代次数，梯度最小的数量级gtol，目标函数要求精度tol，以及该工区速度反演的最值与最小值vmax,vmin。在进行分频并行频域全波形反演的过程中由MATLAB(科学计算软件)的主处理器按照相应的集中调度动态任务分配策略，分配如上信息到各个从处理器上，按照d步骤构造的目标函数和梯度利用超记忆梯度法在每个进程中互不干扰的进行更新迭代运算，从处理器完成第一组频率的分频并行全波形反演。

超记忆梯度法：利用当前梯度及之前多步梯度的信息来对当前梯度的信息进行校正，以加快目标函数的收敛速度。其迭代形式为：

m_k+1＝m_k+α_kd_k

超记忆梯度法下降方向可以通过如下形式来表示：

当k＞mem，β_k-i+1∈[0,s_k ⁱ],(i＝1,…mem)

这里的s_k ¹＝1,i＝2,…mem；

其中，0＜ρ＜1，mem＞0是一个正整数，mem我们称之为记忆度。不同记忆度目标函数收敛速度及模型反演精度不同，记忆梯度数量过多则会导致梯度方向校正过度，目标函数下降变慢，记忆梯度数量少则退化为共轭梯度法。经过实验验证得到，在地震全波形反演中最好的记忆度大致在mem＝4,5。超记忆梯度法有着一定的优势，但当梯度记忆过多，计算速度变慢，内存需求增大。所以选择合适的记忆度才能发挥超记忆梯度法最大的优势。

j、在完成第一组频率的全波形反演之后，需要对混合震源编码方式进行从新编码，在编码的过程中满足随机要求。主处理器将新的震源编码方式重新分配到各个从处理器中，这样不断更新混合震源的编码信息我们称之为动态编码。

震源编码：随机编码就是地表超级炮激发出的震源相位和振幅都是随机分配，设随机相位编码信息为：

其中s'∪s＝{1,2…N_s}，s中元素的个数一定，表示一个超级炮中所含单炮数目一定，可根据实际需求调整混合比例可适当增减超级炮中单炮个数，s中元素大小随机分配。randn表示在[-1,1]之间随机分布的数字。那么所得到的随机编码超级炮可以表示为：

spfs＝Γ·fs

其中fs为震源。动态编码的意思就是在超记忆梯度法每迭代几次之后重新更换随机编码信息，最终可实现地表全覆盖。动态随机编码策略能够有效地增加地表炮的覆盖密度，减弱串扰噪声对全波形反演的影响，同时减少波动方程正演模拟次数。

k、主处理器收集各个进程上不同频率的反演结果并取各个进程结果的平均值meanv₀，然后主处理器再将meanv₀重新分配到各个从处理器中，并将第二组反演频率2f₁,2f₂…2f_Num以及这几个频率所对应的地震记录信息分配到各个从处理器中。还需要分配一些相应的参数如：该频率最大迭代次数，梯度最小的数量级gtol，目标函数要求精度tol，以及该工区速度反演的最值与最小值vmax,vmin。在进行分频并行频域全波形反演的过程中由MATLAB(科学计算软件)的主处理器按照相应的集中调度动态任务分配策略，分配如上信息到各个从处理器上，按照d步骤构造的目标函数和梯度利用超记忆梯度法在每个进程中互不干扰的进行更新迭代运算，从处理器完成第二组频率的分频并行全波形反演。

l、重复k步骤不断的增大频率直到反演频率达到设定的要求nf₁,nf₂…nf_Num，或满足精度要求则终止计算并输出结果。

实施例1

根据勘探要求，将行工具箱和计算软件分散式运算引擎(R2013a)在Windows 7旗舰版系统下进行安装，进行计算软件并行平台的搭建。

利用速度模型进行测试，因原始速度模型较大考虑到硬件设备有限，对原始模型进行抽稀处理，利用抽稀之后的速度模型进行全波形反演测试。真实模型(附图2)和时间域经验模态分解全波形反演结果(图3(右))。

模型参数如下：

表1 多域分频并行多尺度全波形反演测试参数

模型大小

网格距

横向距离

纵向深度

速度范围

起始频率

频带范围

128*384

50m

19.2km

6.4km

1.5-4km/s

5Hz

5-12Hz

模型网格大小为128×384，网格距dz＝dx＝50m，横向距离为19.2km，纵向深度为6.4km，模型中地震波速度范围从1.5km/s到4km/s，地震检波器安置在模型表面，每一个网格点都是一个检波器，且检波器之间的间隔为50m，模型表面震源位置随机分布，一个超级炮中含有38个震源。震源选用12Hz主频的雷克子波，采样间隔为0.001s，实际采样总长度为4.5s，频率范围从5Hz到12Hz。

分频并行拉普拉斯傅里叶域全波形反演参数如下：

选择四个初始频率，分别为4.9Hz，5Hz，5.1Hz和5.2Hz，这么选择的目的是为了防止实际数据中因缺失低频成分或者是低频成分不可靠而导致的反演出现错误现象发生，同时选择多个频率进行反演，最后选择符合地质背景的结果，这样能明显节约寻找可靠起始频率的时间。由于本次进行的只是模型测试，以及硬件设备有限，只能开四核并行计算，所以只选择了4个初始频率，当然如果能开更多的核，可以再多加一些初始频率这样能很大程度上节约选择起始频率的时间。设定拉普拉斯衰减项从σ＝5开始，每隔0.2计算一次，一直到σ＝0为止(σ＝5:-0.2:0)，每个衰减系数σ最大迭代次数为30次。

分频并行频域全波形反演参数如下：

为了克服全波形反演过程中陷入局部极小和跳周等问题，同时也为了提高成像质量，采用多尺度频率全波形反演，将反演频率分为20个频率组，每个频率组包含有4个频率，每个频率迭代30次，相邻频率组之间重叠度为50％，即每个频带中包含的频率分别为f₁→f₄，f₂→f₆，f₄→f₈，f₆→f₁₀……f₄₀→f₄₄(其中脚标表示44个频率的索引位置)。

设定频域全波形的起始频率为5Hz，频率每隔0.2Hz计算一次，一直到频率f＝12Hz停止(f＝5:0.2:12Hz)，每个频率最大迭代次数为30次。

基于以上参数，在表2所示的测试环境下进行全波形反演。

表2 性能测试环境

表3 多域分频并行多尺度全波形反演对比结果

从图4中可以看出时间域经验模态分解全波形反演得到的平滑初始模型相对于线性初始模型构造清晰了很多。利用时间域经验模态分解全波形反演的平滑初始模型可以简化全波形反演的非线性，让反演结果更接近真实的模型。结合表3可以得出，常规全波形反演与本专利提出分频并行随机震源编码全波形反演在模型反演在误差上差别不大，证明了本发明对于处理全波形精度问题的有效性，但是如果从计算时间的角度上可以看出常规的频域全波形反演的计算时间接近6小时，而经过本发明的方法处理之后时间可以缩短30分钟，可见从计算时间的角度来分析，该方法显著提高了计算效率，为实际数据大模型全波形反演打下良好计算基础。

基于以上参数，在表2所示的测试环境下进行多种影响因素测试：

表4 不同影响因子多域分频并行多尺度全波形反演测试结果

注：表3和表4中模型反演归一化误差计算公式为：

从表4中可以看出本发明利用超记忆梯度类优化算法能有效避免噪声的影响，从含噪声数据模型反演误差可以看出，噪声对全波形反演有一定的影响在可允许的范围类全波形反演，其反演结果不会受到很大的影响。当信噪比低于5时，整个全波形反演将不能进行下去，因为所得到的反演结果远远偏离真实模型，反演的结果甚至比初始模型还差，这就表明地震数据在进行全波形反演之前一定要把相应的噪声从地震数据中清除掉，以免它对全波形反演造成严重影响。

从表4中可以看出初始模型的不同对反演结果有着一定的影响，初始模型给的越好最终反演的精度越高，当然初始模型很差则全波形反演这种局部优化的迭代方式将不能满足当前反演的要求。所以本文给出两种合乎地下情况的初始模型，一种是时间域经验模态分解全波形反演的结果，从图4中可以看出时间域经验模态分解全波形反演建立的平滑初始模型结果比线性初始模型好得多，好的初始模型全波形反演很有利。另一种是线性递增初始模型，线性递增初始模型与真实模型误差为0.0089，从图9、图10中都可以看出，线性递增初始模型的反演结果明显不如时间域经验模态分解全波形反演建立的初始模型好，但也在基本能把速度模型反演出来，证明了本发明即使在初始模型较差的情况下也能得到很好的反演结果，证明了该方法的有效性。

从表4中可以看出缺失低频的条件下对反演结果有着一定的影响，从图14中拉普拉斯-傅里叶域全波形反演建立的初始模型，可以明显看出图14中的高精度初始模型与图12所示高精度初始模型存在明显的却别，这是由于当低频缺失的情况下对整个全波形反演都有着严重的影响。本发明利用分频并行拉普拉斯-傅里叶反演来建立初始模型，该反演方法的好处在于能够由浅入深逐渐进行反演，这样能够减少每一步的反演参数，即可有效减小全波形反演的非线性，从而得到一个高精度初始模型，然后再利用分频并行频域多尺度全波形反演可有效避免跳周现象。

图1是整个反演过程的流程图，从流程图可以看出先利用时间域经验模态分解全波形反演可以得到一个较好的平滑初始模型，然后再利用分频并行拉普拉斯-傅里叶域全波形反演获得一个高精度初始模型，最后利用分频并行频域多尺度全波形反演获得最终的反演结果。

Claims (1)

1.一种多域分频并行多尺度全波形反演方法，其特征在于，包括以下步骤：

a、搭建计算软件并行工作库的安装环境，并安装计算软件并行计算工具箱；

b、对实际采集的地震数据进行预处理；

c、首先在预估速度范围建立线性递增初始模型，利用经验模态分解重构地震记录中的低频信息，其次利用时间域经验模态分解全波形反演得到平滑初始模型，再利用拉普拉斯-傅里叶域全波形反演获得高精度初始模型；

d、根据最小二乘原理构造目标函数：

其中m为模型参数，d_obs为实际采集的观测数据，d_cal在速度模型通过正演得到的计算数据，在求目标函数的梯度过程中对目标函数两端关于模型参数m求导，得到梯度表达式：

$\▿C\left(m\right)=-\mathrm{Re}\[\frac{\∂A}{\∂m}{u}^T{\left(A^{-1}\right)}^T(d_{obs}-d_{cal})\]$

其中(A^-1)^T(d_obs-d_cal)为残差反传波场，u表示正传波场；

e、时间域地震记录进行FFT变换得到频域地震数据，查看地震数据的频谱记录，由于实际数据中缺失低频成分，所以在选择频率的时候尽量从低频开始，按照频率从低到高的顺序要求，依次挑选对应频率的地震信号；

根据要求设定多域分频并行多尺度全波形反演相关参数，包括模型大小nz×nx，网格距dx,dz，最大采样时间Nt，时间采样间隔dt，拉普拉斯衰减系数σ，反演起始频率f₀，反演频率个数nf，每个频率最大迭代次数itermax，梯度最小的数量级gtol，目标函数要求精度tol，速度反演的最大值vmax与最小值vmin；

f、利用开放的并行计算软件开启系统Num个并行进程，包括一个主进程和Num个从进程，按照每个程序的需求分配给每个进程相应的内存空间；

g、给定频域地震数据，时间域经验模态分解全波形反演得到的平滑初始模型以及相关参数，进行拉普拉斯-傅里叶域全波形反演，将平滑初始模型v₀、不同反演初始频率nf₀在拉普拉斯-傅里叶域利用分频并行全波形反演得到多个不同低频反演结果，选择符合地质背景的反演结果用于下一步分频并行频域全波形反演；

h、在每个从服务器中完成拉普拉斯-傅里叶域全波形反演之后，主处理器收集各个进程反演的结果，按照技术指标及工区地质背景选择初始模型；

i、主处理器将拉普拉斯-傅里叶域全波形反演得到的高精度初始模型lapv₀重新分配到各个从处理器中，并将第一组反演频率1f₁,1f₂…1f_Num以及这几个频率所对应的地震记录信息分配到各个从处理器中，还需要分配一些相应的参数如：该频率最大迭代次数，梯度最小的数量级gtol，目标函数要求精度tol，以及速度反演的最大值vmax与最小值vmin；

j、在完成第一组频率的全波形反演之后，需要对混合震源的编码方式进行从新编码，在编码的过程中满足随机要求，主处理器将新的震源编码方式分配到各个从处理器中；

k、主处理器收集各个进程上不同频率的反演结果并取各个进程结果的平均值meanv₀，然后主处理器再将meanv₀重新分配到各个从处理器中，并将第二组反演频率2f₁,2f₂…2f_Num以及这几个频率所对应的地震记录信息分配到各个从处理器中，在进行分频并行频域全波形反演的过程中由计算软件的主处理器按照对应的集中调度动态任务分配策略，分配以上信息到各个从处理器上，按照d步骤构造的目标函数和梯度利用超记忆梯度法在每个进程中互不干扰的进行更新迭代运算，从处理器完成第二组频率的分频并行全波形反演；

l、重复k步骤不断的增大频率直到反演频率达到设定的要求nf₁,nf₂…nf_Num，或满足精度要求则终止计算并输出反演结果。