CN107765308A - 基于褶积思想与精确震源的重构低频数据频域全波形反演方法 - Google Patents

Abstract

本发明涉及一种基于褶积思想与精确震源的重构低频数据频域全波形反演方法，首先应用观测数据的直达波成分反演出可靠的震源函数，然后基于褶积思想重构出自激自收点的低频数据，最后应用此重构数据进行多尺度全波形反演的低频起始部分，使反演能有一个好的开始，最终得到一个好的结果。有效避免了周期跳跃的发生，使反演结果得到了极大改善，各种错误极小收敛均有效恢复。在最后的反演部分使用了动态随机编码的超级炮思想，有效提高了全波形反演的效率。本发明重构的自激自收点低频信息不仅是数学层面上的低频信息，只为全波形反演提供大区域扰动，更是客观物理层次上的真实低频信息，也可用于其他数据处理方向。

Description

基于褶积思想与精确震源的重构低频数据频域全波形反演 方法

技术领域

本发明涉及一种地震勘探的地下成像，尤其是基于褶积思想和精确震源重构地震数据缺失的低频信息，并结合超级炮和多尺度思想进行频域全波形反演。

背景技术：

随着油气开发的不断深入，常规的成像技术得到的速度模型已经无法满足在寻找油气矿过程中和开采过程中对地下高精度速度模型的需求。在此背景下，结合计算机技术的飞速更新，全波形反演(FWI)得到了广泛关注和迅速发展，成为当今地球物理学届的一大研究热点。上世纪80年代Tarantola(1984)提出了时间域全波形反演，首先提出了使用伴随状态法求解目标函数，使FWI在计算上成为可能，但受限于当时的计算机技术，时域FWI没有得到很大的发展进步。到了20世纪90年代，Pratt提出了频域全波形反演，使用有限个频率的数据结合从低频到高频的多尺度思想就能得到令人满意的反演结果。在效率上得到了质的飞跃，使FWI正式走到了世界前沿，但FWI是一个强非线性问题，使用局部优化算法导致FWI对初始模型和低频信息非常依赖，在实际数据中如果没有良好的低频信息和初始模型构建，那么反演就极其容易陷入局部极小中。

全波形反演是一个使用地震波场全部数据的数据拟合过程，几乎使用了地震记录中的所有有效信息包括其他方法所摒弃的多次波等等，采用数学的局部优化算法与非线性搜索技术，最终找到一个速度模型，使由其产生的模拟数据与实际数据差异最小。多尺度全波形反演策略在理论上很好地解决了反演中的跳周问题，从低频段开始反演，逐渐增加反演频率至高频，最终使得反演结果越来越精确。但是实际地震数据中并没有好的低频信息，反演不能从低频开始，这导致反演的结果很容易出现跳周的现象。Wu等(2013，2014)在时间域提出了利用地震记录的包络中的丰富低频信息进行多尺度反演，Hu等(2017)在前者的基础上提出了利用解调包络代替简单的Hilbert包络，得到了更加精确的结果，但时间域反演计算量巨大，效率太低。Shin和Cha(2009)提出了在Laplace域全波形反演来得到一个可靠的初始模型后，再在该模型的基础上进行频域FWI，在一定程度上解决了缺失低频信息的问题。Hu(2014)提出了利用波的干涉原理，使用两个相对高频的数据干涉得到一个相对低频的模拟数据来进行全波形反演，但在具体的数值运算时很难体现出干涉的物理意义。

现有的频域全波形反演都是基于多尺度方法进行，即先从低频开始反演，在实际数据中一般缺失5Hz以下频率地震信息，而5-7Hz信息有的情况下也不一定可靠，故在实际应用中，频域全波形反演没办法拥有好的开始，也就难以获得好的结果。

发明内容：

本发明的目的针对上述现有技术的不足，提供一种基于褶积思想与精确震源的重构低频数据的全波形反演方法。

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

基于褶积思想与精确震源的重构低频数据全波形反演方法的核心是在MATLAB2013a平台上，首先应用观测数据的直达波成分反演出可靠的震源函数，利用该精确震源与褶积反褶积理论重构出自激自收点的低频数据；然后基于褶积思想重构出自激自收点的低频数据；最后应用此重构数据进行多尺度全波形反演的低频起始部分，得到的含大尺度构造的速度模型作为新的初始模型进行普通的全波形反演。

基于褶积思想与精确震源的重构低频数据全波形反演方法，包括以下步骤：

a、安装MATLAB基础安装环境；

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

c、在现有资料的基础上，计算出炮检点所在层位的速度模型；

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

其中f为震源，U(ω)为实际采集的观测数据，S(ω)在速度模型上通过正演得到的模拟数据，对目标函数做关于模型参数f求导，并经过一系列变换后，得到最终的梯度表达式：

J＝[A^-1]^T(S(ω)-U(ω))

其中[A^-1]^T(S(ω)-U(ω))为残差反传波场，T代表共轭转置；

e、对预处理后地震记录，截取直达波部分，进行FFT变换得到频域观测数据，查看观测数据的频谱记录，由于实际数据中缺失5HZ以下成分，所以数据尽量从可靠的频率开始选用，按照频率从低到高的顺序要求，依次挑选对应频率的地震信号；

根据要求设定全波形反演相关参数，包括模型大小nz×nx，网格距dx,dz，最大采样时间Nt，时间采样间隔dt，反演起始频率f₀，反演频率个数nf，每个频率最大迭代次数itermax，迭代步长的最小值；

f、设定初始震源，如没有可靠先验信息，则全频率均设为10^-6；

g、给定截去低频数据的频域地震数据，炮检点所在层位的速度模型以及相关参数，进行频域全波形反演，得到不含低频的精确震源Z；

h、将自激自收记录近似看作是由褶积模型得到的记录，U＝Z·H，将地震记录U与得到的精确震源Z相除，得到地层反射系数H，观察地层反射系数随频率的变化可发现，在中低频段内地层反射系数是接近的，故利用已知雷克子波R的低频分量与选择的地层反射系数H相乘，可得重构低频数据UC，

i、利用现有的数据，在合理的区域速度范围内给出光滑的初始速度模型V₀；

j、使用雷克子波R与重构的低频数据进行全空间的频域全波形反演，直到有可靠观测数据的频率为止，得到改善后的初始速度模型V_X；

k、使用改善后的初始模型V_X，不含低频的精确震源Z，给定截去低频数据的频域地震数据和各种参数，进行频域全波形反演；

l、重复k步骤不断的增大频率直到反演频率达到设定的截止频率，终止计算并输出反演结果。

有益效果：本发明成功地通过直达波信息反演出震源形态，提出用褶积思想重构出真实地震数据中缺失的低频成分，最后将此重构数据成功应用于频域全波形反演中，得到了令人满意的成果。

该方法是利用区域现有资料建立炮检层的速度模型，再利用直达波信息使用伴随法反演出精确震源，然后利用该精确震源和褶积思想重构出缺失的低频信息，再利用重构出的低频数据反演出精度较高的初始速度模型。与常规频域全波形反演相比，本发明反演出了震源信息，并重构出了可靠的低频数据，有效避免了周期跳跃的发生，使反演结果得到了极大改善，各种错误极小收敛均有效恢复。解决了以下问题：

1、利用直达波信息和伴随法成功得到了精确的震源函数。不仅可以用于此处的重构低频信息，还可以用于其他地震处理解释中。

2、利用褶积思想与得到的震源信息重构出地震数据中缺失的低频成分，虽然只有自激自收点处的数据，但是却做到了零的突破，该低频成分与之前的重构数据不同之处在于它符合物理条件，即它是真实的物理层面上的低频信息而不是只能用于全波形反演的数学层面上的低频信息。

3、使用重构的低频数据，在初始模型的基础上进行低频段的全波形反演得到一个精度更高的速度模型，将其作为全波形反演的初始模型，利用一个高精度初始模型在一定程度上缓解了全波形反演的跳周现象。为全波形反演获得地下高精度速度模型奠定坚实基础。

4、动态随机震源编码策略有效地压制了因使用超级炮技术而产生的串扰噪声，使全波形反演能正确反演出真实的地下速度。

基于褶积思想与精确震源的重构低频数据有效缓减了全波形反演在实际应用中遇到的最大困难——缺失低频信息，在全波形反演的工业化进程中迈出了一大步。

附图说明

图1基于褶积思想与精确震源的重构低频数据全波形反演方法流程图。

图2真实模型图。

图3初始模型图。

图4震源函数反演的效果对比图。

图5某自激自收点的地层反射系数谱。

图6所有自激自收点的5HZ数据重构效果对比图。

图7缺失7HZ以下数据的全波形反演结果图。

图8使用2.5～7HZ重构数据的全波形反演结果图。

图9反演结果单道速度对比图

(左)50道速度对比，(右)140道速度对比图。

图10反演结果单道速度对比图

(左)250道速度对比，(右)320道速度对比图。

图11缺失7HZ以下数据的全波形反演结果图(抗噪能力测试结果)。

图12使用2.5～7HZ重构数据的全波形反演结果图(抗噪能力测试结果)。

具体实施方式

下面结合附图和实例对本发明进一步的详细说明基于褶积思想与精确震源的重构低频数据频域全波形反演方法，包括以下步骤：

a、程序是在MATLAB2013a软件框架下编写完成，安装基础包便可运行。

b、实际采集的地震观测数据需进行预处理，预处理包括：B1、多次波衰减、面波切除、消除交混回响和压制虚反射等不能模拟的地震波。B2、低频保护去噪、缺失地震道补偿、保幅处理等。

c、如果是海洋数据，并且在海水中放炮和安排检波点的话，那么炮检层的速度可确定为1500m/s；如果是陆地数据，那么需要通过采集炮点附近的检波器上的初至波到达时间，通过公式v＝dx/t求取炮检层的速度。从而建立起炮检层的速度模型，以反演震源函数。

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

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

二维频域常密度声波方程为

其中u为声波波场，F是震源矩阵，是震源函数f对应模型空间的一个矩阵。可以简化为

Au＝F

在反演震源时一般都是单炮反演，故F＝f。将上式两端同时对震源f求导得到：

其中阻抗矩阵A与震源f无关，故有所以有代入梯度公式中，最终得到目标函数的梯度为：

其中S(ω)-U(ω)即为伴随震源函数，根据以上公式发现目标函数的梯度可以通过入射波场与反向残差波场做互相关运算得到，这样避免了求取Jacobi矩阵，很大程度上节约了计算成本。

e、时间域地震记录中挑选震源附近的道集并截取其直达波部分进行FFT变换得到频域地震数据，查看地震数据的频谱记录，由于实际数据中缺失低频成分，所以在选择频率的时候尽量从可靠的低频开始(如7HZ)，按照频率从低到高的顺序要求，依次挑选对应频率的地震信号。

e、根据技术指标及工区要求设定全波形反演的相关参数，包括模型大小nz×nx，网格距dx,dz，最大采样时间Nt，时间采样间隔dt，反演起始频率f₀，反演频率的最大值，每个频率最大迭代次数，步长的最小值λ。

f、设定初始震源，如果没有可靠先验信息，则全频率均设为10^-6。

g、给定频域截取的直达波地震数据，炮检层的速度模型以及相关参数包括初始震源z₀、反演初始频率f₀(如7HZ)等等，进行反演震源的全波形反演。以0.1HZ为间隔不断加大频率至反演频率的最大值。得到不含低频信息的精确震源Z(ω)，如图4。

h、对自激自收点的数据按频率与精确震源相除得到自激自收点的地层反射系数，H＝U/Z。观察H(ω)在中频段的形态，挑选适中的H近似作为低频段的地层反射系数，如图5。

i、选择一已知雷克子波的低频段(2.5～7HZ)，与选择的地层反射系数H相乘，重构出自激自收点低频段的信息，UC＝R·H。图6显示了低频段中5HZ时的所有自激自收点的重构效果对比。可以看到虽然有一定的误差存在，但仍在可接受范围之内。

j、输入初始模型、重构的自激自收点低频数据、以及使用的已知雷克子波进行重构低频段(2.5～7)的全波形反演，以0.1HZ为间隔，通过伴随法求取梯度，非线性搜索确定迭代步长，每个频率迭代3次，将输出的速度模型作为新的初始模型。

k、将新的速度模型代入常规的频域全波形反演中去，进行中高频段内的数据反演。得到的速度更新结果如图8所示，对比于缺失低频的普通全波形反演结果，图7，我们可以看到左上区域和左中部出现的明显局部极小现象得到了很大的改善，左右侧的高速区块处错误反演出的低速区成功复原，底部的背斜构造层位速度正确归位。

实施例1

根据勘探要求，将MATLAB Distributed Computing Server(R2013a)在Windows 7旗舰版系统下进行安装，进行MATLAB运行平台的搭建。

利用Marmousi复杂模型进行测试，由于原始Marmousi模型巨大，而电脑硬件设备有限，对原始模型进行抽稀处理，对抽稀之后的Marmousi模型进行频域的全波形反演测试。真实模型(附图2)和初始模型(附图3)。

模型参数如下：

表1频域全波形反演测试参数

模型大小

网格距

横向距离

纵向深度

速度范围

起始频率

最大频率

121*384

25m

9600m

3025m

1.5～4.5km/s

7HZ

25HZ

模型网格大小为121×384，网格距dz＝dx＝25m，横向距离为9.6km，纵向深度为3.025km，Marmousi模型中速度范围从1.5km/s到4.5km/s，地震检波器安置在水下50m深处，检波器之间的间隔为25m，震源位于水下100m深处，炮点间隔25m。模拟震源选用10Hz主频的雷克子波，采样间隔为0.01s，实际采样总长度为10s，选用频率范围从7Hz到25Hz。在低频段重构数据时，采用单炮记录，因为要重构出所有的自激自收点数据；在中高频段的普通全波形反演中，采用超级炮思想，每个频率重新定义随机振幅和随机相位编码，以减少串扰噪声的影响。

基于褶积思想与精确震源的重构低频数据频域全波形反演参数如下：

选择主频为10HZ的雷克子波作为震源，并且舍弃7HZ以下数据，以模拟真实地震数据，甚至如此设计比一般真实数据缺少了更多的低频数据。在进行全波形反演时，以7HZ起始，0.25HZ为间隔直至25HZ。每个频率最多迭代8次，期间如果步长<0.1(初始步长为1)那么就退出迭代，进入下个频率的反演。

表2是电脑测试环境。

图11与图12是抗噪性测试的结果，从图中可以看到对于含噪声数据，全波形反演对低频信息显得更加依赖，使用基于褶积思想与精确震源的重构低频数据进行频域全波形反演可以很好地克服传统全波形反演的缺点，即使观测数据含有噪声，最终也能反演出基本准确的速度模型，为其他技术提供支持。

图1是整个反演过程的流程图，从流程图可以看出先利用震源附近的直达波信息反演出精确震源，然后重构出自激自收点的低频信息，再使用重构的数据得到较好的初始速度模型，最后代入混采超级炮多尺度频域全波形反演中得到理想的最终速度模型结果。

Claims (1)

1.一种基于褶积思想与精确震源的重构低频数据频域全波形反演方法，其特征在于，包括以下步骤：

a、安装MATLAB基础安装环境；

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

c、在现有资料的基础上，计算出炮检点所在层位的速度模型；

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

其中f为震源，U(ω)为实际采集的观测数据，S(ω)在速度模型上通过正演得到的模拟数据，对目标函数做关于模型参数f求导，经变换后得到最终的梯度表达式：

J＝[A^-1]^T(S(ω)-U(ω))

其中[A^-1]^T(S(ω)-U(ω))为残差反传波场，T代表共轭转置；

e、对预处理后地震记录，截取直达波部分，进行FFT变换得到频域观测数据，查看观测数据的频谱记录，由于实际数据中缺失5HZ以下成分，因此数据尽量从可靠的频率开始选用，按照频率从低到高的顺序要求，依次挑选对应频率的地震信号；

根据要求设定全波形反演相关参数，包括模型大小nz×nx，网格距dx,dz，最大采样时间Nt，时间采样间隔dt，反演起始频率f₀，反演频率个数nf，每个频率最大迭代次数itermax，迭代步长的最小值；

f、设定初始震源，如没有可靠先验信息，则全频率均设为10^-6；

g、给定截去低频数据的频域地震数据，炮检点所在层位的速度模型及相关参数，进行频域全波形反演，得到不含低频的精确震源Z；

h、将自激自收记录近似看作是由褶积模型得到的记录，U＝Z·H，将地震记录U与得到的精确震源Z相除，得到地层反射系数H，观察地层反射系数随频率的变化，在中低频段内地层反射系数是接近的，故利用已知雷克子波R的低频分量与选择的地层反射系数H相乘，得重构低频数据UC，

i、利用现有的数据，在合理的区域速度范围内给出光滑的初始速度模型V₀；

j、使用雷克子波R与重构的低频数据进行全空间的频域全波形反演，直到有可靠观测数据的频率为止，得到改善后的初始速度模型V_X；

k、使用改善后的初始模型V_X，不含低频的精确震源Z，给定截去低频数据的频域地震数据和各种参数，进行频域全波形反演；

l、重复k步骤不断的增大频率直到反演频率达到设定的截止频率，终止计算并输出反演结果。