带限格林函数滤波多尺度全波形反演方法与流程
本发明涉及一种地震勘探的地下成像方法，尤其是利用梯度预条件方法实现多尺度全波形反演的带限格林函数滤波全波形反演方法。
背景技术：
：
全波形反演方法，是一个不断用模拟数据去拟合实际数据的过程，最终找到一个速度模型使得数据残差平方最小。它具有成像精度高、反演模型参数准确、复杂构造成像效果好等优点。为了顺应地震勘探高精度成像的要求，全波形反演迅速发展起来，并成为当今地球物理界的研究热点。20世纪80年代，塔兰托拉(1984)提出了时间域FWI方法，并通过伴随状态法求取目标函数的梯度，避免了雅克比矩阵的求取，大大缩短了FWI梯度计算所需要的时间。20世纪90年代，普拉特将FWI推广到了频率域，提出只需要几个离散的频率就可以得到高精度的反演结果，且低频到高频的多尺度反演策略可以解决陷入局部极小值的问题。
全波形反演就是利用地震数据中所有波形的信息，来进行地震模型反演，利用所有的信息约束得到的反演结果更接近真实的地质模型。地震全波形反演在实际数据应用中仍然面临着诸多难题，它是一个强非线性问题，对低频成分、震源函数、噪声和初始模型非常敏感(维里厄,2009)，这给FWI应用于实际数据处理带来了巨大的挑战。当地震数据中缺失低频成分或者伴随震源主频过高时，地震全波形反演结果都很容易出现周波跳跃现象。当前地震全波形反演的研究热点主要是如何缓解反演过程中的周波跳跃问题。导致地震全波形反演产生周波跳跃的主要因素有：一、地震数据为振荡的复频信号，且具有较高的主频，波形变化剧烈；二、地震数据中缺失低频成分；三、震源函数估计不准确，四、用于FWI的初始速度模型不准确。
为了解决地震全波形反演的周波跳跃问题，邦克斯等(1995)提出多尺度全波形反演策略，将地震数据按照反映地质体尺度的特性分成大，中，小三个尺度。并提出利用低通滤波的方法实现了从低频到高频的多尺度全波形反演，有效地缓解了全波形反演存在的周波跳跃问题。但是实际采集到的地震数据中常常不含有5Hz以下的低频信息，即使利用低通滤波滤除地震数据中的高频成分得到相对平缓的地震波形，然而当初始模型较差的情况下，该多尺度策略仍然会出现周波跳跃问题。吴等(2013)利用地震数据包络重构出了地震数据中缺失的低频成分，利用地震数据包络能够得到地下地质体的宏观构造。在利用包络全波形反演得到的速度模型作为常规全波形反演的初始速度模型，最终可以得到一个高精度的反演结果。蔡等(2016)利用模拟数据与观测数据互相关后的结果与观测数据自相关的结果做差，建立了基于优化互相关目标函数的全波形反演，该方法有效缓解了全波形反演中的数据不匹配问题。
技术实现要素：
：
本发明的目的针对上述现有技术的不足，提供一种不需要超低频成分便就能反演得到较准确的地质构造，有效缓解了全波形反演的周波跳跃问题的带限格林函数滤波多尺度全波形反演方法。
本发明的目的是通过以下技术方案实现的：
带限格林函数滤波多尺度全波形反演方法的核心是建立对应的全波形反演目标函数。本发明是利用模拟数据自褶积的结果再与观测数据和模拟数据互褶积的结果做差，来建立格林函数滤波全波形反演的目标函数。根据所建立格林函数滤波的目标函数再对反演参数求偏导数，便能得到目标函数对应参数的梯度。从实际数据提取子波的过程可能会与真实子波存在延时效应，本发明利用观测数据和模拟数据互相关得到峰值偏离中心的差值来校正模拟数据与观测数据存在的延时效应。利用时差校正以后的模拟数据带入到带限格林函数滤波多尺度全波形反演的目标函数中，继而求得模型更新的梯度。
全波形反演的核心是利用模拟数据不断地去拟合观测数据的过程，最终找到一个速度模型使得数据残差平方最小。在反演的过程中需要不断的进行正演模拟得到新的模拟数据与观测数据进行匹配，然而全波形反演的正演模拟过程相当耗时，这在一定程度上制约了地震全波形反演在大数据中的应用。为了加快全波形反演的进程，本发明利用炮与炮之间并行计算，充分发挥计算机多线程计算的优势。分炮并行时间域全波形反演方法的核心是一个主处理器和若干从处理器，主处理器监控整个全波形反演计算阶段，由主处理器将适应度的计算分配到各个从处理器上进行，计算完成之后再由主处理器收集正演模拟数据和正传波场用于计算全波形反演的梯度。时间域波动方程正演模拟炮与炮之间的计算是独立进行的互不干扰，利用分炮可以充分发挥当前计算机多线程并行计算的优秀性能，显著缩短了全波形反演的计算时间。
带限格林函数滤波多尺度全波形反演方法，包括以下步骤：
a、输入地震观测数据(Dobs)和观测系统，本发明以真实模型模拟得到的数据作为实际采集的地震数据输入到计算机系统中。观测系统模拟陆地地震数据采集的观测系统，即速度模型表面每一个网格分布一个检波器，每隔10个网格放置一个震源；
b、利用巴特沃斯高通滤波器滤除地震数据7Hz以下的低频信息，该步骤主要的目的是模拟真实地震数据采集过程中采集不到低频信息的情况。同时能够验证本发明提出的方法能够在缺失低频情况下，有效缓解全波形反演的周波跳跃问题；
c、提取实际速度模型的最大值和最小值，建立线性递增初始模型。该方法建立的初始速度模型不含有真实模型任何构造信息，大大提高了全波形反演的难度，同时能够验证本发明提出的带限格林函数滤波多尺度全波形反演方法的鲁棒性。实际数据处理过程中，该步骤可以利用速度分析来得到一个相对较好的初始速度模型，能够有效提高真实数带限格林函数滤波多尺度全波形反演的准确性；
d、从观测数据中提取震源子波，用于全波形反演的正演模拟。在实际数据处理中提取震源子波的方法有多种，例如：结合测井数据提取地震子波，利用地震数据自相关得到地震子波，利用最小二乘反演得到地震子波等；
e、根据实际计算机的配置打开并行进程，其中包括一个主进程和N个从进程，按照每个程序的需求分配给每个进程相应的内存空间。将主进程接收的观测数据，地震子波和初始速度模型按照既定的观测系统分配到各个从处理器，进行时间域有限差分正演模拟计算；
f、在每个从服务器中完成正演模拟之后，主处理器收集各个进程的计算结果，其中包括模拟数据(Dcal)和正传波场(Pf)；
g、利用观测数据与模拟数据互相关来计算观测数据与模拟数据的延时差。当观测数据与模拟数据不存在延时差的时候，两个数据互相关的峰值正好在互相关序列的中间位置。当观测数据与模拟数据存在一定延时差的时候，两个数据的互相关峰值会偏离中心位置，偏离中心位置的时间就是观测数据与模拟数据对应的延迟时间(Δτ)。利用将求得的延时差，对模拟数据进行延时校正；
h、利用观测数据与模拟数据互褶积(Dcal*Dobs)，即可得到带限格林函数滤波以后的观测数据。利用模拟数据的自褶积(Dcal*Dcal)，即可得到带限格林函数滤波以后的模拟数据。
i、根据带限格林函数滤波多尺度全波形反演方法的目标函数：
其中v为模型的速度参数，在求目标函数的梯度过程中对目标函数两端关于模型参数v求导，得到梯度表达式：
其中为伴随震源，由于震源附近地震道的能量较强，再经过褶积与互相关之后的值会变得更大，为了避免这种现象的发生，需要将震源附近的2道伴随震源数据充零；
j、将伴随波场与时差校正以后的正传波场做零延迟互相关，得到带限格林函数滤波多尺度全波形反演的梯度。
k、根据设定的预条件梯度系数(Alpha)，来对所求得的梯度进行预处理；
l、用L-BFGS优化算法计算速度模型更新方向，并通过Wolfe收敛准则寻找合适的步长；
m、判断是否满足终止条件，若满足则输出带限格林函数滤波多尺度全波形反演结果；若不满足终止条件，将反演结果作为下一个循环的初始模型，返回(5)步骤；
n、得到好的初始速度模型后，再进行常规时间域全波形反演，输出最终反演结果。
有益效果：本发明成功地将带限格林函数滤波多尺度等技术应用到了地震全波形反演中，提出的分炮并行时间域全波形反演，节约了大量的时间。
该方法是利用初始速度模型的格林函数来观测数据和模拟数据进行滤波，目的是让观测数据和模拟数据波形更加接近，避免出现周波跳跃问题。但是复杂介质的格林函数没办法直接表示出来，本发明把初始模型的格林函数与地震子波褶积之后再与观测数据和模拟数据作用，这样问题就简化成了一个褶积问题。在全波形反演的开始阶段利用格林函数对地震数据进行滤波，使得观测数据与模拟数据在波形上更加接近，结合预条件梯度方法对梯度做适当的平滑处理，便可以反演出模型的宏观信息。通过调整预条件梯度的系数来实现多尺度的反演策略。本发明同时利用分炮并行的反演策略，目的是缩短全波形反演过程中正演模拟所需要的时间。与常规频域全波形反演相比，带限格林函数滤波多尺度全波形反演方法，解决了以下问题：
①、利用互相关来校正模拟数据与观测数据的延时差，可以有效避免观测数据与模拟数据出现不匹配的现象，同时避免周波跳跃现象。
②、利用格林函数对观测数据和模拟数据进行滤波，可以有效减弱全波形反演的非线性。利用格林函数滤波可以理解为先反演得到一个过渡的速度模型，然后再以这个过渡的速度模型作为初始速度模型进行常规全波形反演。带限格林函数滤波策略，有效的缓解了因缺失低频而带来的跳周现象。
③、利用分炮并行的正演模拟策略，极大程度的加快了全波形反演的计算速度。原先时间域单炮正演模拟，需要等上一炮正演结束以后再进行下一炮正演模拟，这种串行的方式无法充分发挥当前计算机多线程的优秀性能。分炮并行的正演策略，显著缩短了全波形反演的计算时间，使得在时间允许的情况下可以加密全波形反演的炮数，以及扩大全波形反演的勘探面积。
④、利用预条件梯度多尺度策略，可以实现先更新模型的宏观信息，通过调整预条件系数来不断刻画模型的细节信息。该多尺度反演策略结合带限格林函数滤波反演方法，有效缓解了全波形反演的周波跳跃问题。
附图说明
图1带限格林函数滤波多尺度全波形反演方法流程图。
图2雷克子波形波图和雷克子波频谱图
(左)雷克子波波形图，(右)雷克子波频谱图。
图3互相关时差校正示意图。
图4速度模型
(左)真实速度模型，(右)线性递增初始模型。
图5全波形反演结果
(左)缺失低频常规时间域全波形反演结果，
(右)缺失低频带限格林函数滤波多尺度全波形反演。
图6带限格林函数滤波多尺度全波形反演中间过程图
(左)预条件梯度系数α＝0.3，(右)预条件梯度系数α＝0.6。
图7在图5(左)速度模型上正演数据的格林函数滤波地震记录。
图8在图5(右)速度模型上正演数据的格林函数滤波地震记录。
具体实施方式
下面结合附图和实例对本发明进一步的详细说明
带限格林函数滤波多尺度全波形反演方法，包括以下步骤：
a、输入地震观测数据(Dobs)和实际数据采集过程中的观测系统，本发明以真实模型模拟得到的数据作为实际采集的地震数据输入到计算机系统中。观测系统模拟陆地地震数据采集的观测系统，即速度模型表面每一个网格分布一个检波器，每隔10个网格放置一个震源。
利用野外实测数据进行全波形反演时需要对地震数据进行一些的预处理，其中包括：多次波衰减、面波切除、消除交混回响和压制虚反射等不能模拟的地震波。低频保护去噪、缺失地震道补偿、保幅处理，根据野外实测数居提取地震子波等；
b、利用巴特沃斯高通滤波器滤除地震数据7Hz以下的低频信息，该步骤主要的目的是模拟真实地震数据采集过程中，采集不到低频信息的情况。同时能够验证本发明提出的方法能够在缺失低频情况下，有效缓解全波形反演的周波跳跃问题；
巴特沃斯高通滤波器：以巴特沃斯函数来近似滤波器的系统函数。巴特沃斯滤波器是根据幅频特性在通频带内具有最平滑特性定义的滤波器。巴特沃思滤波器的高通模平方函数表示：
巴特沃斯滤波器的主要特征如下：
对所有的对所有的N，Ha(jΩ)＝0.707即20lg|Ha(jΩ)＝3dB，|Ha(jΩ)|2是Ω的单调上升函数，|Ha(jΩ)|2随着阶次N的增大而更接近于理想高通滤波器。滤波器的幅频特性随着滤波器阶次N的增加而变得越来越好，在截止频率Ωc处的函数值始终为1/2的情况下，通带内有更多的频带区的值接近于1；在阻带内更迅速的趋近于零。
c、提取实际速度模型的最大值和最小值，建立线性递增初始模型。该方法建立的初始速度模型不含有真实模型任何构造信息，大大提高了全波形反演的难度，同时能够验证本发明提出的带限格林函数滤波多尺度全波形反演方法的鲁棒性。实际数据处理过程中，该步骤可以利用速度分析来得到一个相对较好的初始速度模型，能够有效提高真实数带限格林函数滤波多尺度全波形反演的准确性；
根据技术指标及工区要求设定多域分频并行多尺度全波形反演的相关参数，包括模型大小nz×nx，网格距dx,dz，最大采样时间Nt，时间采样间隔dt，最大迭代次数Itermax，梯度最小的数量级gtol，目标函数要求精度tol，以及该工区速度反演的最大值与最小值vmax,vmin。
其中v0代表初始速度，vmin表示速度模型中的最小速度，vmax表示速度模型中的最大速度，nz代表z方向网格数目，Nz代表z方向网格总数目。
d、从观测数据中提取震源子波，用于全波形反演的正演模拟。在实际数据处理中提取震源子波的方法有多种，例如：结合测井数据提取地震子波，利用地震数据自相关得到地震子波，利用最小二乘反演得到地震子波等；
例如最小二乘反演的震源子波反演的目标函数：
其中表示震源附近直达波模拟记录，表示震源附近直达波观测记录。震源子波反演的目标函数对震源子波f的导数可表示为：
时间域声波方程两边同时对震源子波f求导得到：
其中正演算子L与震源子波无关，有即波场对震源子波f的导数：
得到震源子波反演的更新梯度：
其中震源函数伴随震源定义为：本文利用伴随状态法计算雅克比矩阵大大缩短了反演震源子波的时间。
fn＝fn-1+αΔf (10)
其中fn表示第n次迭代的震源子波值，α代表更新步长，Δf表示震源子波更新梯度。
e、根据实际计算机的配置打开并行进程，其中包括一个主进程和N个从进程，按照每个程序的需求分配给每个进程相应的内存空间。将主进程接收的观测数据，地震子波和初始速度模型按照既定的观测系统分配到各个从处理器，进行时间域有限差分正演模拟计算；
f、在每个从服务器中完成正演模拟之后，主处理器收集各个进程的计算结果，其中包括模拟数据(Dcal)和正传波场(Pf)；
g、利用观测数据与模拟数据互相关来计算观测数据与模拟数据的延时差。当观测数据与模拟数据不存在延时差的时候，两个数据互相关的峰值正好在互相关序列的中间位置。当观测数据与模拟数据存在一定延时差的时候，两个数据的互相关峰值会偏离中心位置，偏离中心位置的时间就是观测数据与模拟数据对应的延迟时间(Δτ)。利用将求得的延时差，对模拟数据进行延时校正；
互相关校正时差：本发明利用观测数据与模拟数据进行互相关，与中心点的差值，即为震源函数的延后值，再对观测记录进行整体校正，避免因延时过大，使得目标函数无法下降。互相关函数可以定义成如下所示：
其中A(xr,xs)obs是D(xr,t+τ,xs)obs的最大振幅。τ表示观测数据与模拟数据之间的时移。取互相关最大值的地方对应的时间，再与互相关以后时间序列的中心位置时刻做差，即是观测数据与模拟数据的延时差：
当Δτ＝0表示两者不存在延时问题。
h、利用观测数据与模拟数据互褶积(Dcal*Dobs)，即可得到带限格林函数滤波以后的观测数据。利用模拟数据的自褶积(Dcal*Dcal)，即可得到带限格林函数滤波以后的模拟数据。
格林函数滤波：利用格林函数滤波可以使得全波形的非线性程度减弱，避免出现周波跳跃现象。则格林函数滤波全波形反演的目标函数可以定义为：
其中G0表示初始速度所对应的格林函数。但是地震勘探是针对复杂地质模型而言的，复杂模型的格林函数很难被求出来，为了便于计算将目标函数褶积上正演模拟的子波f0。则原先的格林函数滤波与地震子波褶积以后得到了一个带限带的地震记录，如下所示：
其中有：
Dcal＝G0*f0 (15)
基于格林函数滤波的目标函数可以表示为：
其中v表示模型的速度参数，ns表示震源的个数，nr表示检波器的个数，T表示最大接收时间。
i、求目标函数的梯度过程中对目标函数两端关于模型参数v求导，得到梯度表达式：
得到伴随震源：
由于震源附近地震道的能量较强，再经过褶积与互相关之后的值会变得更大，为了避免这种现象的发生，需要将震源附近的2道伴随震源数据充零。其中伴随波场可以表示为：
j、将伴随波场与时差校正以后的正传波场做零延迟互相关，得到带限格林函数滤波多尺度全波形反演的梯度。
其中Pf表示经过时差校正以后正传波场，Pb伴随震源的反传波场。
k、根据设定的预条件梯度系数(Alpha)，来对所求得的梯度进行预处理；
预条件梯度的处理方法有很多种，例如对梯度做校正，平滑，或者利用交叉梯度等。在这里详细叙述一下本发明利用高斯频域平滑预条件梯度，来处理带限格林函数滤波全波形反演的梯度。高斯频域平滑函数为：
其中nz代表z方向网格数目，Nz代表z方向网格总数目，nx代表x方向网格数目，Nx代表x方向网格总数目。α为预条件梯度系数。对应的预条件梯度为：
g＝S*g0 (24)
l、用L-BFGS优化算法计算速度模型更新方向，并通过Wolfe收敛准则寻找合适的步长；
L-BFGS算法是一种迭代型的算法，迭代更新公式如下：
mk+1＝mk-αkHkgk
其中mk为模型更新参数，αk为Wolfe线性搜索得到的步长，Hk为近似Hessian矩阵的逆矩阵，gk为模型更新梯度。L-BFGS优化更新只需要保存少数的向量对用于更新Hessian矩阵，其更新公式如下：
Hk+1＝VkTHkVk+ρkskskT (25)
sk＝mk+1-mk，yk＝gk+1-gk (27)
其中Hk+1是根据向量对{sk,yk}和Hk计算得到，只储存n个向量对来隐式表达Hessian矩阵的逆矩阵。的乘积可以通过梯度gk与向量对{sk,yk}之间一系列向量的内积与向量的和来获得的。在每一次更新完成后，上一步向量对将被当前新向量对{sk+1,yk+1}取代。因此，向量对集合中包含最近n次迭代的曲率信号。在时实际中，当n≥5时都能获得较满意的结果。L-BFGS优化算法的初始近似Hessian矩阵每一次迭代中都不同。公式(20)中的近似Hessian矩阵的逆矩阵Hk需满足以下更新公式：
模型的更新方向，通过如下简易伪代码实现：
通过上式伪代码求得模型的更新量Hkgk，然后再通过Wolfe线性搜索获得合适的步长αk进行更新迭代。L-BFGS优化算法有效的克服了显式求取Hessian矩阵及其逆的困难，其具有超线性收敛速度，计算效率高，占用内存小，精度高等优点，较适合求解大规模非线性优化问题。
m、判断是否满足终止条件，若满足则输出带限格林函数滤波多尺度全波形反演结果；若不满足终止条件，将反演结果作为下一个循环的初始模型，返回(5)步骤；
判断是否满足终止条件：该频率最大迭代次数itermax，梯度最小的数量级gtol，目标函数要求精度tol。
n、得到好的初始速度模型后，再进行常规时间域全波形反演，输出最终反演结果。
实施例1
测试模型参数如下：
基于二维时间域常密度声波方程，对发明对带限格林函数滤波多尺度全波形反演方法进行测试。对Marmousi模型进行抽稀处理，测试模型网格横纵大小为69×192，模型网格间距为12.5m，横向距离为2.4km，纵向距离为0.8625km。震源沉放深度为50m，共20炮，均匀分布在地表。检波器沉放深度为12.5m，共192个检波器，即每个网格一个检波器。震源选用22Hz主频的雷克子波，采样间隔为0.001s，时间采样总长度为2.5s。本文利用水平线性模型作为FWI的初始速度模型，速度变化范围从1.5km/s到4km/s
分炮并行全波形反演参数如下：
当前计算机集群具有多线程同时处理优势，本发明利用分炮并行策略加快全波形反演的计算速度。选择4线程台式计算机进行测试，20个炮点均匀分布在地表。计算机的每个线程中可以进行单炮有限差分正演模拟，也就是说4个线程的计算机同时可以计算4炮的地震数据，20个炮点分五次分发的4个线程中进行正演模拟。全波形反演每一次模型更新迭代需要正演模拟2次正演模拟，那么分炮并行反演策略，在只利用4核计算机运算时就可以提高8倍的计算效率。如果对于大型计算机集群，其计算效率可以得到大大的改善。
预条件梯度多尺度全波形反演参数如下：
为了克服全波形反演过程中陷入局部极小和跳周等问题，同时也为了提高成像质量，采用预条件梯度多尺度全波形反演，通过调整预条件系数来实现先更新模型的长波长成分，在不断增大预条件梯度系数来逐渐恢复模型的精细结构。预条件系数首先选择为α＝0.3，得到全波形反演结果如图6(左)，然后再增大预条件系数(α＝0.6)如图6(右)。最终的反演结果如图5(右)所示，可以看出带限格林函数滤波多尺度全波形反演能够很好的恢复速度模型的细节构造。利用相同的预条件系数对常规全波形反演的梯度进行处理，得到如图5(左)所示的结果，可以看出带限格林函数滤波多尺度全波形反演较常规全波形反演结果，有效避免了周波跳跃问题。
图1是整个反演过程的流程图，从流程图可以看出带限格林函数滤波全波形反演结合预条件梯度方法能够从一个较差的初始速度模型中得到一个较精确的反演结果(图5)。对比图7和图8可以看出，利用带限格林函数滤波多尺度全波形反演结果作为速度模型在进行正演模拟得到的观测数据的格林函数滤波后的数据更加相近。
- 还没有人留言评论。精彩留言会获得点赞！