基于非线性高次拓频的时间域多尺度全波形反演方法

摘要

本发明涉及一种基于非线性高次拓频的时间域多尺度全波形反演方法，首先对地震数据取高次方；然后对高次方运算后的地震数据低通滤波，通过正传波场与反传波场的零延迟互相关运算得到梯度；通过共轭梯度法等更新模型，以此为初始模型，进行数据二次方反演，待迭代结束得到最终反演结果。本发明利用不同阶次重构的低频数据，实现了多尺度全波形反演，能有效克服全波形反演的跳周问题。避免了对数据原始频带的依赖，本发明提出的波形压缩和数据剥层方法实际上是基于能量的，因此无论地震记录的延时长短，进行高次方运算都能对每一个单炮记录进行低频重构和数据剥层，解决了由于时间延迟而导致的部分炮信息的损失问题。

说明书

技术领域：

本发明涉及一种利用非线性高次拓频方法进行缺失低频重构，并进行多尺度全波形反演的方法。非线性高次运算通过非线性方法拓展有限带宽地震数据的频带，提取和利用重构的低频信息能有效克服全波形反演的跳周问题。

背景技术：

全波形反演通过最小化观测数据与模拟数据的残差来驱动速度模型的更新。根据波恩近似，这个强非线性问题可以基于弱散射近似进行线性化，从而可以通过局部优化算法求解。这就要求反演有一个足够接近真实模型的初始模型。但是，实际处理中，常规的速度分析方法精度有限。另一个能降低反演对初始模型依赖性的办法是采用低频数据进行反演。低频数据的波恩近似误差较小，速度扰动和数据残差的线性关系更好。相应的，低频数据目标函数对应更少的局部极值，能有效防止反演的解陷入局部极小。因此低频数据对于全波形反演至关重要。

常规全波形反演已经有多尺度反演策略，即频带从低到高依次反演。用低频段反演的结果作为下一个高频段反演的初始模型，这样可以有效缓解反演的跳周问题。但是，常规多尺度反演的问题是，假如地震数据缺失低频信息，单纯靠低通滤波方法是无法获得反演需要的低频段信息的。因此，仍然无法缓解反演的跳周问题。

在实际采集中，炸药或气枪激发的震源都是有限带宽的，震源中往往不含有效的低频信息。常规震源中，一般认为5Hz以下的信息是不可靠的。低频震源的成本十分昂贵，目前还未见大规模应用。因此，目前来看，地震勘探采用的震源依然为低频缺失震源。如何在低频情况下反演介质模型的长波长信息是全波形反演必须解决的难题。进行缺失低频信息重构是一类重要的解决方法。在地震数据缺失低频情况下，前人已经提出一些方法重构缺失的低频信息，并用于全波形反演。Xie(2013.Recover certain low-frequencyinformation for full waveform inversion.83th SEG Annual InternationalMeeting,Expanded Abstracts,1053-1057.)在低频数据满足线性相位假设的情况下，研究了数据频率外推方法，成功重构出了部分缺失的低频信息，并能在一定程度上缓解全波形反演的跳周问题。在反射脉冲响应满足稀疏假设的情况下，Fei(2012.Full waveforminversion without low frequencies:a synthetic study.82th SEG AnnualInternational Meeting,Expanded Abstracts,641-645.)和Warner(Full-waveforminversion of cycle-skipped seismic data by frequency down-shifting.84th SEGAnnual International Meeting,Expanded Abstracts,903-907.)等人提出利用稀疏反演方法可以获得稀疏解，同时拓宽了地震记录的频带，基于此的全波形反演能有效避免局部极值解。Shin和Cha(2008.Waveform inversion in the Laplace domain.GeophysicalJournal International,173,922-931.)提出了Laplace域波形反演方法，根据地震数据的衰减特性提取出长波长信息，并用于更新介质的大尺度构造。Wu等人(2014.Seismicenvelope inversion and modulation signal model.Geophysics,79(3),WA13-WA24.)提出了地震数据包络反演，通过提取地震数据的包络，获得反演介质长波长分量的信息，并用于更新模型的宏观信息。

在以上方法中，虽然能在一定程度上缓解全波形反演的跳周问题，但是都存在不同的局限性。低频线性外推和稀疏反演方法分别需要相位线性假设和反射响应稀疏假设。Laplace波形反演计算量较大，而且衰减参数的提取复杂。包络反演伴随源频带较宽，会引入大量高频信息，影响长波长分量的构建。不同于前人的方法，本发明主要基于信号处理领域的相关知识进行低频重构。在时间域越窄的信号，对应越宽的频带。极端情况为脉冲函数，其在时间域对应无限窄的宽度，而其频带为无限宽。在频域无限窄的信号为单频信号，其对应的时间域信号为无限宽的正余弦信号。据此，我们可以采用压缩时间域波形的方法来拓宽信号的频带，然后利用拓展的低频信息来进行全波形反演。

发明内容：

本发明的目的是针对上述前人方法技术的不足，提出一种新的低频数据重构及全波形反演方法。

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

首先，对地震数据取高次方。考虑到波场传播时的几何扩散和地层吸收效应，地震记录随着时间的增加，波场能量越小。因此，从三次方开始，可以在保留有效信息的同时，较大化的压缩波形。然后，对高次方运算后的地震数据低通滤波，以提取重构的低频数据成分。构建新的伴随源，低通滤波后用以计算反传波场。计算震源的正传波场。正传波场与反传波场的零延迟互相关运算得到梯度。通过共轭梯度法等可以更新模型。当达到一定精度后，以此结果为初始模型，进行数据二次方的反演。反演流程与三次方相同。得到一定精度的结果后，进行数据一次方(数据本身)的全波形反演。待迭代结束后，得到最终的反演结果。本发明利用不同阶次重构的低频数据，实现了多尺度全波形反演，能有效克服全波形反演的跳周问题。

本发明所述的基于非线性高次拓频的时间域多尺度全波形反演方法是通过MATLAB平台实现的。

本发明所述的基于非线性高次拓频的时间域多尺度全波形反演方法，主要包括以下步骤：

a、安装MATLAB软件平台，要求采用MATLAB R2016a及以上版本。并且已配备并行工具包；

b、对数据进行静校正、去噪、去除多次波预处理，得到高质量的地震数据；

c、对地震数据进行子波估计，提取每一炮数据的震源子波；

d、通过速度分析，得到模型速度的最大和最小值，生成一个在速度范围内，随深度逐渐递增的线性梯度模型，作为反演的初始模型；

e、对于公式(1)所示的目标函数，设置n的值1≤n≤3，对观测数据进行n次方运算以拓频，通过频谱分析确定重构低频提取的低通滤波参数f_low，其中，σ_new表示目标函数值，v表示速度，n表示次数，d_cal与d_obs分别表示模拟数据与观测数据；

f、在初始模型上计算模拟数据，根据公式(1)计算目标函数值σ₀，设置最大迭代次数k_max；

g、进行第k(k≥1)次迭代。用震源子波在初始模型上模拟得到正传波场，用作为伴随源，进行参数为f_low的低通滤波后，进行反传得到伴随波场。计算正传波场与伴随波场的零延迟互相关得到本次更新的梯度g_k；

h、选择合适的步长α，更新模型m_k＝m_k-1+α·d_k。在新模型上，计算目标函数值σ_k；

i、如果σ_k＜σ_k-1且k＜k_max，则令k＝k+1，进行第g步；如果σ_k≥σ_k-1，则令α＝α/2，进行第h步；如果k＝k_max且n≥2，则令n＝n-1，进行第f步；如果k＝k_max且n＝1，则进行第j步；

j、输出最终结果m_k。

有益效果：

本发明利用高次方运算对时间域波形的压缩作用，成功重构出记录的缺失低频信息，并实现了新的多尺度全波形反演策略，能更有效的克服全波形反演的跳周问题，能在地震数据低频缺失情况下，依次反演得到速度模型的大尺度和小尺度构造信息。

本发明与现有技术相比具有的实质性特点：

1.本发明的缺失低频重构方法无需任何理论上的假设条件，可以应用于任何接收条件下的地震记录，也可以应用于任何低频缺失程度的地震数据。

2.根据地震记录的特点，可以采用不同的阶次进行低频重构，操控简便。

3.对地震记录进行高次方运算，在压缩波形的同时还兼具数据剥层的效果。高次方运算能突出数据间能量的差异。因此，越高次的运算，越能体现浅部记录的能量。阶次由高到低，对应能量的由浅到深，起到了时间域剥层法的作用。无论从低频角度，还是数据剥层角度，本发明对于克服全波形反演的跳周都是有帮助的。

4.通过调节数据的阶次，可以实现多尺度全波形反演。这样，避免了对数据原始频带的依赖，根据数据频带的不同，可以选择不同的运算阶次，以重构出足够的低频记录，用于全波形反演。

5.本发明可以用于混合采集地震记录的反演。混合采集的地震记录是混叠的，而本发明提出的波形压缩和数据剥层方法实际上是基于能量的，因此无论地震记录的延时长短，进行高次方运算都能对每一个单炮记录进行低频重构和数据剥层，避免了由于时间延迟而导致的部分炮信息的损失。

附图说明：

图1基于非线性高次拓频的时间域多尺度全波形反演方法的流程图。

图2缺失低频雷克子波(上)及其频谱(下)。

图3不同阶次低频重构效果。

(a)n＝1对应的地震子波波形；

(b)n＝2对应的地震子波波形；

(c)n＝3对应的地震子波波形；

(d)n＝11对应的地震子波波形；

(e)n＝1对应的地震子波频谱；

(f)n＝2对应的地震子波频谱；

(g)n＝3对应的地震子波频谱；

(h)n＝11对应的地震子波频谱。

图4目标函数性态曲线图。

(a)低通滤波前不同阶次目标函数性态曲线；

(b)30Hz低通滤波后不同阶次目标函数性态曲线；

(c)20Hz低通滤波后不同阶次目标函数性态曲线；

(d)10Hz低通滤波后不同阶次目标函数性态曲线。

图5速度模型。

(a)真实速度模型；

(b)初始速度模型。

图6不同阶次地震记录及其频谱。

(a)n＝1对应的地震记录；

(b)n＝2对应的地震记录；

(c)n＝3对应的地震记录；

(d)不同阶次记录的频谱。

图7不同阶次伴随源频谱。

图8反演结果。

(a)常规全波形反演结果；

(b)n＝3对应的反演结果；

(c)以(b)为初始，n＝2对应的结果；

(d)以(c)为初始，n＝1对应的结果。

图9反演结果抽道对比。

具体实施方式：

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

本发明所述的基于非线性高次拓频的时间域多尺度全波形反演方法，包括以下步骤：

a、在win7或Linux系统下安装MATLAB软件平台，要求采用MATLAB R2016a及以上版本，并且已配备并行工具包(Parallel Computing Toolbox)；

b、进行数据预处理。对数据进行静校正处理，校正起伏地表对反射同相轴的影响；对数据进行去噪处理，去除微震、低频和高频背景噪声及其他随机噪声；去除干扰波，包括声波、面波、工业电干扰、虚反射、多次反射、侧面波、底波、交混回响和鸣震。最终得到高质量的地震数据；

c、对地震数据进行子波估计，估计方法可以采用直达波估计法和自相关法等，提取出每一炮数据的震源子波；

d、通过速度分析，得到模型速度的最大和最小值，生成一个在速度范围内，随深度逐渐递增的线性梯度模型，作为反演的初始模型。初始模型的表达式可以写为(2)式

v₀(i)＝v_min+(i-1)*(v_max-v_min)/(n-1),>

其中，v₀为初始速度模型值，v_min和v_max分别为速度分析得到的最小和最大速度，i为模型纵向网格坐标，n为模型纵向最大网格点数；

e、对于公式(1)所示的目标函数，设置n的值(1≤n≤3)。对观测数据进行n次方运算以拓频，通过频谱分析确定重构低频提取的低通滤波参数f_low。

构建基于n次高阶数据的全波形反演目标函数为

其中，σ_new表示目标函数值，v表示速度，n表示次数，d_cal与d_obs分别表示模拟数据与观测数据。这里的和实际上分别代表对模拟数据的n次方和观测数据的n次方进行了参数为f_low的低通滤波处理。

对公式(1)所示的目标函数，计算其对速度的偏导数，得到

其中，是敏感核函数，与常规全波形反演的敏感核函数相同。不同于常规全波形反演的是，本发明中的伴随源不再是数据的残差，而是新的残差项

f、在初始模型上计算模拟数据，根据公式(1)计算目标函数值σ₀。设置最大迭代次数k_max；

g、进行第k(k≥1)次迭代。用震源子波在初始模型上模拟得到正传波场，用作为伴随源，进行参数为f_low的低通滤波后，进行反传得到伴随波场。计算正传波场与伴随波场的零延迟互相关得到本次更新的梯度g_k。

模型梯度的表达式可以写为

其中，g表示梯度，v表示速度，x_s表示震源索引，t表示时间，u表示震源正传波场，B表示伴随源的反传波场。由梯度表达式可知，反演的梯度是正传波场与反传波场的零延迟互相关，然后沿着时间积分，对所有炮的结果进行加和，最后作用上系数得到。

在得到梯度后，采用共轭梯度法求解更新方向，求解方法为

d(m_k)＝-g(m_k)+β·d(m_k-1),>

其中，

h、选择合适的步长α，更新模型m_k＝m_k-1+α·d_k。在新模型上，计算目标函数值σ_k；

i、如果σ_k＜σ_k-1且k＜k_max，则令k＝k+1，进行第g步；如果σ_k≥σ_k-1，则令α＝α/2，进行第h步；如果k＝k_max且n≥2，则令n＝n-1，进行第f步；如果k＝k_max且n＝1，则进行第j步；

j、输出最终结果m_k。

实施例1：

本发明的整体流程如图1所示。

采用如图2所示的雷克子波，从频谱分析结果可见，子波主频为20Hz，缺失8Hz以下的低频信息。采用本发明的高次低频重构方法。采用n＝1，2，3，11的重构结果分别如图3(a)-3(d)所示。可见，经过高次方运算后，地震波的波形变窄。图3(a)-3(d)的频谱分析结果分别如图3(e)-3(h)所示。可见，当n>2时，8Hz以下的能量均被不同程度的恢复。当n较大时，重构记录越来越接近于宽频脉冲函数。采用低通滤波可以提取重构出的低频成分。

下面分析目标函数的性态。真实模型采用100*100的均匀二维模型，网格间距为10m，速度为2500m/s，震源位于z＝900m，x＝500m位置。采样间隔为0.8ms，记录时间为1s。初始模型也为二维均匀速度模型，其他参数与真实模型一样，速度变化范围为2000m/s-3000m/s，速度间隔为20m/s。模拟观测记录和合成记录，取n＝1，2，3，计算目标函数曲线如图4a所示。可见，由于高次方运算后，在低频重构同时，也产生了高频信息，因此目标函数收敛区域更窄。采用低通滤波提取低频信息进行实验。对模拟记录和观测记录分别进行滤波参数为30Hz，20Hz和10Hz的低通滤波，计算目标函数曲线，分别如图4b，4c和4d所示。可见，重构出的低频信息可以帮助减少目标函数的极值，从而帮助克服反演的跳周问题。

实验采用的真实模型如图5a所示。初始模型为最大和最小速度的线性插值结果，如图5b所示。采用图2的低频缺失子波，进行正演模拟，得到的地震记录如图6a所示。对地震记录取高次方运算，n＝2和3的结果分别如图6b和6c所示。图6a-6c均采用了相同参数的低通滤波以压制高频成分。对图6a-6c进行频谱分析，结果如图6d所示。可见，n＝2和3时，均在低频缺失频段重构出了能量。根据公式(3)中伴随源的形式，计算三种情况的伴随源，并作相同参数的低通滤波，进行频谱，结果如图7所示。可见，n＝1，即常规全波形反演的伴随源仍然缺失低频成分，而n＝2和n＝3时，伴随源含有丰富的低频信息。

下面进行反演计算。首先采用常规全波形反演，在初始模型上迭代200次，得到结果如图8a所示。可见，反演结果在浅部出现严重假象，说明，由于低频缺失，反演过程跳周严重。然后，令n＝3，采用本发明的方法迭代200次，反演结果如图8b所示。可见，浅部的构造已经变得清晰。以图8b的结果为初始模型，进行n＝2的反演，迭代200次，结果如图8c所示。可见，宏观构造已经较清楚的构建出来。以图8c为初始模型，进行n＝1的反演，迭代200次的结果如图8d所示。可见，最终结果细节和深部构造清晰，反演效果明显好于常规反演结果。