CN108254731A - 探地雷达数据的多尺度阶梯式层剥离全波形反演方法 - Google Patents

Abstract

本发明涉及一种多尺度阶梯式层剥离全波形反演方法，首先，选取模型分层数目、反演频率及单频单层的反演次数，依据新的反演序列确定各层各频率的反演序号；其次,采用维纳滤波器将雷达观测数据与子波滤到不同频段，采用介电常数梯度公式，对计算出的梯度添加对应层的空间窗，获取当前层的梯度值；最后计算模型更新量，逐层逐频完成迭代运算。采用多尺度策略减少了对初始模型的依赖，有效解决了周波跳跃问题；有别于现有的偏移距‑时间窗，采用空间窗口获取特定层梯度，可快速实现收敛；新的阶梯式反演序列，有效地减小两个相邻层反演频率之间差距，从而获得更好的反演结果，有效改善模型深部的反演效果，提高全波形反演的稳定性。

Description

探地雷达数据的多尺度阶梯式层剥离全波形反演方法

技术领域：

本发明涉及一种探地雷达领域中用于重建介电常数的时间域全波形成像方法。结合频率多尺度策略和常规层剥离反演序列，给出了一个新的阶梯式反演序列来提高全波形反演的稳定性，并减少周波跳跃问题。

背景技术：

探地雷达(Ground Penetrating Radar,GPR)数据的全波形反演(Full WaveformInversion,FWI)是近地表参数成像的一种有前景的技术。基于广义最小二乘法的全波形反演方法最早出现在Tarantola(1984.Linearized inversion of seismic reflectiondata.GEOPHYSICAL PROSPECTING 32(6).998-1015)的开创性论文中，该方法为时域全波形反演技术奠定了坚实的理论基础。随后被广泛应用于地震勘探数据处理中，但在电磁数据处理应用方面仍比较少见。近年来，国内外的专家学者始终致力于将地震全波形反演的核心技术引入到探地雷达成像研究领域。Moghaddam等人(1991.Comparison of the borniterative method and Tarantola's method for an electromagnetic time-domaininverse problem.INTERNATIONAL JOURNAL OF IMAGING SYSTEMS&TECHNOLOGY 3(4).318-333)利用全波形对小目标介电常数成像做了一些尝试，但没有取得令人满意的效果。全波形反演的理论与技术以其高精度、多参数建模的能力，已成为当前探地雷达领域的研究热点。如对层状土壤或分层结构中电导率和介电常数反演，以及对土壤含水量的估算。20世纪90年代，Pratt将Tarantola的理论发展到了频率域，由此奠定了频率域全波形反演发展的基础。而后，采用准牛顿优化策略，学者们提出了用于同时重构介电常数和电导率的二维全波形反演。而后又提出了基于探地雷达和钻孔地震数据的联合全波形反演的创新方法，并继续针对地面雷达数据和浅层地震反射资料的联合全波形反演进行研究，应用于土木工程成像领域。

尽管全波形反演可以提供地下的高分辨率图像,但是存在严重的缺陷，即对初始模型的高度依赖性。此外,由于数据量和非线性程度均较大，常常采用基于梯度类优化算法来实现反演过程的收敛。然而,梯度类算法很容易陷入局部极小值。为了避免这种情况,必须使用合适的初始模型来确保收敛,从而生成与至少半个波长的观测数据重叠的合成数据。为了解决这些限制，Bunks等人(1995.Multiscale seismic waveforminversion.GEOPHYSICS 60(5).1457-1473)提出了一种时域全波形反演的多尺度方案，即利用低频数据,在全波形反演的初始阶段获得准确的平滑背景和大尺度结构，再利用高频数据进一步刻画精细结构，从而保证结果的稳定性和置信度，更重要的是，它减少了对初始模型的依赖。随后，学者们提出了一种高度稳定的策略，称为串行反演。对低频数据进行波形反演，并将其结果作为高频数据后续反演的初始模型。

一般而言，地下模型深部的反演精度取决于浅部的反演精度。浅部反演不准确会严重影响深部的反演效果，这个问题在多尺度全波形反演中无处不在。为了解决这一问题，主要采用一种从上至下进行分层反演的基于数据的层剥离方法，通常结合了偏移和时间窗。该方法的实质在于从上至下使每一层的目标函数最小化，以确保总体目标函数值最小。随后，利用阻抗地震反射数据，开发了一种全新的剥层全波形反演策略来控制有效更新深度(EUD)，避免了持续更新深度控制(CUDC)和射线追踪。Wang等人(2014.Reflectionseismic waveform tomography.JOURNAL OF GEOPHYSICAL RESEARCH SOLID EARTH 114(B3)B03304)又使用波形层析成像技术来获取地下模型的高分辨率图像。然而，选择适当的偏移-时间窗来提取某一固定深度以上的特定波形通常并不容易，并且可能导致分层全波形反演处理在强非均匀介质中的不可靠。此外，在这种策略下，当两个相邻层的反演频率存在较大差异，则两层之间边界附近可能会出现局部最小值。

发明内容：

本发明的目的就在于针对上述现有技术的不足，提供一种多尺度阶梯式层剥离全波形反演方法，用于重建介电常数。在常规的基于数据的层剥离反演的基础上，在由整个观察数据计算出的梯度上添加空间窗口，因此该方法也称为基于模型的层剥离。同时结合多尺度策略，给出了全新的阶梯式反演序列，以提高全波形反演的稳定性，并缓解周波跳跃问题。

本发明的目的是通过以下技术方案实现的：首先，确定并选取模型分层数目、反演频率以及单频单层的反演次数，依据基于模型的多尺度阶梯式层剥离全波形反演的新的反演次序确定各层各频率的反演序号；其次，采用维纳低通滤波器将雷达观测数据与子波滤到不同的频率段；然后，采用TM模式下所推导的介电常数梯度公式，计算整个模型空间的梯度值，并对计算出的梯度添加对应层的空间窗，获取当前层的梯度值；最后，利用最速下降法计算模型更新方向，并采用非精确线搜索方法计算更新步长，更新当前模型，将最后一层内最后一个频带的最后一次反演结果作为整个全波形反演的最终结果。

探地雷达数据的多尺度阶梯式层剥离全波形反演方法是通过MATLAB平台实现的。

探地雷达数据的多尺度阶梯式层剥离全波形反演方法，包括以下步骤：

a、输入作为观测数据的地面雷达记录、初始模型以及子波；

b、选取模型分层数目、反演频率以及单频单层的反演次数；

c、依据多尺度阶梯式层剥离全波形反演的基本原理，确定各层各频率的反演序号；

d、设定初始反演序号为1；

e、确定当前第k次反演所对应的层号与频率(k＝1,2,3……,N)；

f、构造维纳滤波器，用于高频观测数据与高频子波向低频的转换：

式中，f_Wiener代表维纳低通滤波器，W_target为低频目标震源子波，W_original为原始震源子波，ε为控制数值溢出因子；

g、利用步骤f中的滤波器计算低频观测数据与低频模拟数据；

h、参考一阶伴随状态法，推导并采用TM模式下的介电常数梯度公式，计算整个模型空间的梯度值；

式中E_Z代表前向波场，E_Z*代表伴随波场，σ为电导率(单位S/m)，S为震源项，nt则为时间域的采样点；

i、对步骤h中计算出的梯度添加对应层的空间窗，获取当前层的梯度值；

j、利用最速下降法计算模型更新方向；

k、利用步骤i的梯度与步骤j的更新方向，并采用非精确线搜索方法计算更新步长，更新当前模型；

l、反演序号加一，以当前模型作为下一次反演的初始模型，重复步骤e-k，直到反演序号结束，输出最终反演结果。

有益效果：本发明将全波形反演方法引入到探地雷达数据的成像当中，结合频率多尺度策略和常规层剥离的反演频率序列，给出了一个新的阶梯式反演序列，对在TM模式下重建介电常数的应用进行了大量的针对性改进。显著提高了利用多尺度阶梯式层剥离全波形反演方法进行成像在反演精度和分辨率两方面的效果，具体有以下优点：①采用时间域多尺度策略，利用低频数据在反演初始阶段获得准确的平滑背景和大尺度结构，再利用高频数据进一步刻画精细结构，从而减少了对初始模型的依赖，有效缓解了周波跳跃问题；②区别于传统的偏移距-时间窗，由整个观测数据计算出的梯度上添加了空间窗口来获取特定层的梯度，对于避免周波跳跃及实现反演的快速收敛至关重要；③本发明提出的新的阶梯式反演序列，可以有效地减小两个相邻层反演频率之间的差距，从而获得更好的反演结果；④频率多尺度方案与阶梯式反演序列相结合，有效改善模型深部的反演效果，提高全波形反演的稳定性。

附图说明：

图1基于探地雷达数据的多尺度阶梯式层剥离全波形反演方法流程图。

图2两种窗口的比较。

(a)应用于观测数据的组合偏移距-时间窗口；

(b)应用于由整个观测数据计算出的梯度上的基于模型的空间窗口。

图3两种层剥离全波形反演方法的反演频率次序。

(a)常规层剥离全波形反演的反演频率次序；

(b)多尺度阶梯式层剥离全波形反演的反演频率次序。

图4简单介电常数模型。

(a)输入的真实模型；

(b)初始模型，真实模型去除中间两个高值异常体的大尺度平滑结果。

图5简单介电常数模型的三种反演结果对比。

(a)仅使用频率多尺度全波形反演的结果；

(b)常规反演序列的层剥离全波形反演方法的结果；

(c)多尺度阶梯式层剥离全波形反演的结果。

图6简单介电常数模型的两种数据拟合结果。

(a)单道数据拟合，三种不同的反演策略中各在同一位置(x＝17.82m)选择一个单道数据，并与真实模型数据进行数据拟合；

(b)1维深度剖面，选择固定深度h＝5.76m的数据，得到三种不同反演策略的一维深度反演剖面，并与真实模型相同深度的数据进行对比。

图7随机等效介质土壤介电常数模型。

(a)输入的真实模型；

(b)初始模型，真实模型的大尺度平滑结果。

图8随机等效介质土壤介电常数模型的三种反演结果对比。

(a)仅使用频率多尺度全波形反演的结果；

(b)常规反演序列的层剥离全波形反演方法的结果；

(c)多尺度阶梯式层剥离全波形反演的结果。

图9随机等效介质土壤介电常数模型的两种数据拟合结果。

(a)单道数据拟合，三种不同的反演策略中各在同一位置(x＝7.29m)选择一个单道数据，并与真实模型数据进行数据拟合；

(b)1维深度剖面，选择固定深度h＝8.46m的数据，得到三种不同反演策略的一维深度反演剖面，并与真实模型相同深度的数据进行对比。

具体实施方式：

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

探地雷达数据的多尺度阶梯式层剥离全波形反演方法，包括以下步骤：

a、输入作为观测数据的地面雷达记录、初始模型以及子波，精确的初始模型有助于获得更好的反演结果；

b、选取模型分层数目、反演频率以及单频单层的反演次数。这三个参数需要在综合考虑模型大小、构造复杂程度、计算量等因素之后加以确定；

c、依据多尺度阶梯式层剥离全波形反演的基本原理，确定各层各频率的反演序号；基于模型的多尺度阶梯式层剥离全波形反演的基本原理概述如下：假设已选定M个反演频率，记作f₁,f₂,…,f_M，反演模型被等分为N层，记作L₁，L₂，...，L_N，且每层内每个频率的迭代次数均相同，记作i。对于第一层，首先使用f₁进行迭代更新i次，然后利用f₂重复上述反演。接下来我们用f₁在第二层进行第一次迭代i次，然后返回第一层用f₃开始反演，随后第二层用f₂开展反演，接下来第三层开始采用f₁实现第一次反演迭代更新i次。依此反演次序，我们将继续完成对每一层的迭代更新，直到最后一个频率f_M在最后一层L_N中完成最后一次反演运算(例如图3(b)所示)。与常规层剥离全波形反演的反演频率次序相比(例如图3(a)所示),本发明提出的这种反演次序确保了任意两个相邻层之间参与反演的最大频率序号差异为2，大大缩短了频率间的差异；

d、设定初始反演序号为1；

e、确定当前第k次反演所对应的层号与频率(k＝1,2,3……,N)；

f、构造维纳滤波器(Winner Filter)，用于高频观测数据与高频子波向低频的转换：

式中，f_Wiener代表维纳低通滤波器，W_target为低频目标震源子波，W_origina为原始震源子波，ε为控制数值溢出因子；

g、利用步骤f中的维纳滤波器计算低频观测数据与低频模拟数据，所采用的维纳低通滤波器对观测数据与所选雷达子波进行滤波，从而得到低频带信息。低频带反演结果将继续作为高频带反演的初始模型，然后按步骤c中的阶梯式层剥离全波形反演次序进行逐层逐频反演；

h、在反演算法框架中，关键点则是计算梯度算子。在地震勘探的全波形反演中，梯度是由一阶伴随状态法确定的，即梯度是由前向波场和伴随波场的零延迟互相关。其中伴随波场是在接收器位置所激励的时间反转数据残差计算获得的。因此，本专利参考地震勘探中的一阶伴随状态法，推导并采用电磁波TM模式下的介电常数梯度公式，计算整个模型空间的梯度值；

式中E_Z代表前向波场，E_Z*代表伴随波场，σ为电导率(单位S/m)，S为震源项，nt则为时间域的采样点；

i、对步骤h中计算出的梯度添加对应层的空间窗，获取当前层的梯度值；

常规的层剥离全波形反演是基于数据的，即在观测数据上添加组合的偏移距-时间窗口(例如图2(a)所示)来选择观测和合成数据的特定部分来计算相应的梯度值，然后执行层剥离反演。常规层剥离方法从根本上保证了目标函数由浅部到深部是最小的，从而确保整体的目标函数也是最小的。但是，确定这种组合的偏移距-时间窗口的范围通常并不容易；

本发明在利用全部观测数据计算的梯度值上增加空间窗口(例如图2(b)所示)，以获取特定反演层的梯度值，这种窗口避免了对同相轴的识别和提取，从而取代了常规方法在观测数据上添加的偏移距-时间窗口；

j、利用最速下降法计算模型更新方向：

式中，m代表模型参数，这里为介电常数，Δm_k ^SD为下降方向，则第k次迭代的下降方向为目标函数梯度的相反值；

k、利用步骤i的梯度与步骤j的更新方向，并采用非精确线搜索方法计算当前迭代更新步长因子α_k，所谓非精确线搜索，是指选取步长α_k使目标函数f得到可接受的下降量，即对于当前迭代k，存在Δf_k＝f(m_k)-f(m_k+α_kΔm_k)>0是可接受的，其中f(m_k)代表当前模型的目标函数值，最后计算模型更新量，更新当前第k次迭代的模型：

式中，m代表模型参数，这里为介电常数，常量α_k代表迭代更新步长因子，Δm_k ^SD为下降方向；

l、反演序号加一，以当前模型作为下一次反演的初始模型，重复步骤e-k，直到反演序号结束，输出最终反演结果。

实施例1：

a、输入作为观测数据的地面雷达记录、初始模型以及子波；本发明建立了反演的真实模型(图4(a))，来测试和评估阶梯式层剥离全波形反演的工作流程及其效果。真实的介电常数模型：在具有分层背景的模型中心处具有两个高值异常目标体。除两个较高的异常区域外，介电常数随深度增加，且模型介电常数最大值与最小值分别为ε_max＝29，ε_min＝6。模型大小为8.64m×34.02m，在模型顶部具有h_a＝0.45m的空气层，有限差分网格间距h＝9cm，故该模型共计N_M＝97×379＝36763个网格点；地面观测系统由76个发射天线，间距5个网格点，190个接收天线，间距1个网格点组成，且每个天线发射的信号可被全部接收天线记录；输入由真实模型计算的合成数据作为反演的观测数据；输入真实模型去除中间两个高值异常体的大尺度平滑结果作为初始模型(图4(b))；输入中心频率为120MHz的子波；

b、在综合考虑模型大小、构造复杂程度、计算量等因素之后，确定并选取步骤a中模型的平均分层数目为20层；从5-80MHz的范围中选择15个频率；单频单层的反演次数为5；

c、依据多尺度阶梯式层剥离全波形反演的基本原理，确定各层各频率的反演序号；基于模型的多尺度阶梯式层剥离全波形反演的基本原理：根据步骤b中所选定的参数，已选定的15个反演频率，记作f₁,f₂,…,f₁₅，平均分层数目20层，记作L₁，L₂，...，L₃₀，且单频单层的反演次数为i＝5。对于第一层，首先使用f₁进行迭代更新5次，然后利用f₂重复上述反演。接下来我们用f₁在第二层进行第一次迭代5次，然后返回第一层用f₃开始反演，随后第二层用f₂开展反演，接下来第三层开始采用f₁实现第一次反演迭代更新5次。依此反演次序，我们将继续完成对每一层的迭代更新，直到最后一个频率f₁₅在最后一层L₃₀中完成最后一次反演运算(图3(b))。与常规层剥离全波形反演的反演频率次序相比(图3(a)),本发明提出的这种反演次序确保了任意两个相邻层之间参与反演的最大频率序号差异为2，大大缩短了频率间的差异；

d、设定初始反演序号为1；由于本发明仅开展单参数介电常数ε的反演研究，因此反演过程中电导率σ始终保持不变；

e、确定当前第k次反演所对应的层号与频率(k＝1,2,3……,N，其中，N＝15*5*20＝1500)；

f、构造维纳滤波器(Winner Filter)，用于高频观测数据与高频子波向低频的转换：

式中，f_Wiener代表维纳低通滤波器，W_target为低频目标震源子波，W_origina为原始震源子波，ε为控制数值溢出因子；

g、利用步骤f中的维纳滤波器计算低频观测数据与低频模拟数据，所采用的维纳低通滤波器对观测数据与所选雷达子波进行滤波，从而得到低频带信息。低频带反演结果将继续作为高频带反演的初始模型，然后按步骤c中的多尺度阶梯式层剥离全波形反演次序进行逐层逐频反演；

h、在反演算法框架中，关键点则是计算梯度算子。在地震勘探的全波形反演中，梯度是由一阶伴随状态法确定的，即梯度是由前向波场和伴随波场的零延迟互相关。其中伴随波场是在接收器位置所激励的时间反转数据残差计算获得的。因此，本专利参考地震勘探中的一阶伴随状态法，推导并采用电磁波TM模式下的介电常数梯度公式，计算整个模型空间的梯度值；

式中E_Z代表前向波场，E_Z*代表伴随波场，σ为电导率(单位S/m)，S为震源项，nt则为时间域的采样点；

i、对步骤h中计算出的梯度添加对应层的空间窗，获取当前层的梯度值；

常规的层剥离全波形反演是基于数据的，即在观测数据上添加组合的偏移距-时间窗口(图2(a))来选择观测和合成数据的特定部分来计算相应的梯度值，然后执行层剥离反演。常规层剥离方法从根本上保证了目标函数由浅部到深部是最小的，从而确保整体的目标函数也是最小的。但是，确定这种组合的偏移距-时间窗口的范围通常并不容易；

本发明在利用全部观测数据计算的梯度值上增加空间窗口(图2(b))，以获取特定反演层的梯度值，这种窗口对于避免周波跳越问题及实现收敛至关重要，从而取代了常规方法在观测数据上添加的偏移距-时间窗口。因此，本发明的多尺度阶梯式层剥离全波形反演亦可称作基于模型的层剥离全波形反演；

j、利用最速下降法计算模型更新方向：

式中，m代表模型参数，这里为介电常数，Δm_k ^SD为下降方向，则第k次迭代的下降方向为目标函数梯度的相反值；

k、利用步骤i的梯度与步骤j的更新方向，并采用非精确线搜索方法计算当前迭代更新步长因子α_k，所谓非精确线搜索，是指选取步长α_k使目标函数f得到可接受的下降量，即对于当前迭代k，存在Δf_k＝f(m_k)-f(m_k+α_kΔm_k)>0是可接受的，其中f(m_k)代表当前模型的目标函数值，最后计算模型更新量，更新当前第k次迭代的模型：

式中，m代表模型参数，这里为介电常数，常量α_k代表迭代更新步长因子，Δm_k ^SD为下降方向；

l、反演序号加一；以当前模型作为下一次反演的初始模型，重复步骤e-k，直到反演序号结束，输出最终反演结果；

为了突出本发明：基于探地雷达数据的多尺度阶梯式层剥离全波形反演方法的效果，将本发明的反演结果(图5(c))与仅使用频率多尺度全波形反演方法以及采用传统反演序列的层剥离全波形反演方法的结果进行比较(图5(a)和图5(b))。与真实模型相比(图4(a))，可以观察到上述三种方法均能清楚地反映模型的整体异常分布特征，即模型中部两个楔形的高介电常数异常块体及成层分布的背景结构。由仅使用频率多尺度全波形反演方法的结果(图5(a))，可以看出模型中间两个高值楔形异常体的下边界反演得不够准确，在真实模型中，两个块体的下边界约位于深度5.8m，而在仅使用频率多尺度全波形反演方法的结果中，该边界已超过深度为6m的位置，出现明显的向下偏移现象，此外，从(图4(a))与(图5(a))的对比不难发现，仅使用频率多尺度全波形反演方法的反演重建不够准确，中部的两个高值楔形异常的反演值明显低于真实模型。由传统反演序列的层剥离全波形反演方法的结果(图5(b))，不难发现模型中间的高值异常目标体的边界轮廓很模糊，背景的分层效果较差。相比之下，本发明提出的多尺度阶梯式层剥离全波形反演方法的反演(图5(c))，其结果的质量明显优于上述两种方法获得的质量。按照本发明所阐述的此种新的反演次序，经连续反演之后，两个中心目标高值区域被精确重建。另外，模型深部的各个分层界面也都得到很好的表征；

为了验证本发明结果的可靠性，我们从上述三种不同的反演策略中各在同一位置选择一个单道数据，并与真实模型数据进行数据拟合(图6(a))。可以看出，多尺度阶梯式层剥离全波形反演的介电常数的振幅在三者中与真模型最接近。此外，选择固定深度h＝5.76m的数据，得到三种不同反演策略的一维深度反演剖面(图6(b))，不难发现多尺度阶梯层剥离全波形反演的结果也更接近于真实模型。

实施例2：

a、输入作为观测数据的地面雷达记录、初始模型以及子波；本发明基于多参数耦合方程建立了随机土壤介质模型作为反演的真实模型(图7(a))。模型大小为9.54m*9m，在模型顶部具有h_a＝0.36m的空气层；有限差分网格间距h＝9cm，故该模型共计N_M＝107×105＝11235个网格点；地面观测系统由36个发射天线，间距3个网格点，107个接收天线，间距1个网格点组成，且每个天线发射的信号可被全部接收天线记录；输入由真实模型计算的合成数据作为反演的观测数据；输入真实模型的大尺度平滑结果作为初始模型(图7(b))；输入中心频率为120MHz的子波；

b、在综合考虑模型大小、构造复杂程度、计算量等因素之后，确定并选取步骤a中模型的平均分层数目为30层；从5-80MHz的范围中选择15个频率；单频单层的反演次数为4；

c、依据多尺度阶梯式层剥离全波形反演的基本原理，确定各层各频率的反演序号；基于模型的多尺度阶梯式层剥离全波形反演的基本原理：根据步骤b中所选定的参数，已选定的15个反演频率，记作f₁,f₂,…,f₁₅，平均分层数目30层，记作L₁，L₂，...，L₃₀，且单频单层的反演次数为i＝4。对于第一层，首先使用f₁进行迭代更新4次，然后利用f₂重复上述反演。接下来我们用f₁在第二层进行第一次迭代4次，然后返回第一层用f₃开始反演，随后第二层用f₂开展反演，接下来第三层开始采用f₁实现第一次反演迭代更新4次。依此反演次序，我们将继续完成对每一层的迭代更新，直到最后一个频率f₁₅在最后一层L₃₀中完成最后一次反演运算(图3(b))。与常规层剥离全波形反演的反演频率次序相比(图3(a)),本发明提出的这种反演次序确保了任意两个相邻层之间参与反演的最大频率序号差异为2，大大缩短了频率间的差异；

d、设定初始反演序号为1；由于本发明仅开展单参数介电常数ε的反演研究，因此反演过程中电导率σ始终保持不变；

e、确定当前第k次反演所对应的层号与频率(k＝1,2,3……,N，其中，N＝15*4*30＝1800)；

f、构造维纳滤波器(Winner Filter)，用于高频观测数据与高频子波向低频的转换：

式中，f_Wiener代表维纳低通滤波器，W_target为低频目标震源子波，W_origina为原始震源子波，ε为控制数值溢出因子；

g、利用步骤f中的维纳滤波器计算低频观测数据与低频模拟数据，所采用的维纳低通滤波器对观测数据与所选雷达子波进行滤波，从而得到低频带信息。低频带反演结果将继续作为高频带反演的初始模型，然后按步骤c中的多尺度阶梯式层剥离全波形反演次序进行逐层逐频反演；

h、在反演算法框架中，关键点则是计算梯度算子。在地震勘探的全波形反演中，梯度是由一阶伴随状态法确定的，即梯度是由前向波场和伴随波场的零延迟互相关。其中伴随波场是在接收器位置所激励的时间反转数据残差计算获得的。因此，本专利参考地震勘探中的一阶伴随状态法，推导并采用电磁波TM模式下的介电常数梯度公式，计算整个模型空间的梯度值；

式中E_Z代表前向波场，E_Z*代表伴随波场，σ为电导率(单位S/m)，S为震源项，nt则为时间域的采样点；

i、对步骤h中计算出的梯度添加对应层的空间窗，获取当前层的梯度值；

常规的层剥离全波形反演是基于数据的，即在观测数据上添加组合的偏移距-时间窗口(图2(a))来选择观测和合成数据的特定部分来计算相应的梯度值，然后执行层剥离反演。常规层剥离方法从根本上保证了目标函数由浅部到深部是最小的，从而确保整体的目标函数也是最小的。但是，确定这种组合的偏移距-时间窗口的范围通常并不容易；

本发明在利用全部观测数据计算的梯度值上增加空间窗口(图2(b))，以获取特定反演层的梯度值，这种窗口对于避免周波跳越问题及实现收敛至关重要，从而取代了常规方法在观测数据上添加的偏移距-时间窗口。因此，本发明的多尺度阶梯式层剥离全波形反演亦可称作基于模型的层剥离全波形反演；

j、利用最速下降法计算模型更新方向：

式中，m代表模型参数，这里为介电常数，Δm_k ^SD为下降方向，则第k次迭代的下降方向为目标函数梯度的相反值；

k、利用步骤i的梯度与步骤j的更新方向，并采用非精确线搜索方法计算当前迭代更新步长因子α_k，所谓非精确线搜索，是指选取步长α_k使目标函数f得到可接受的下降量，即对于当前迭代k，存在Δf_k＝f(m_k)-f(m_k+α_kΔm_k)>0是可接受的，其中f(m_k)代表当前模型的目标函数值，最后计算模型更新量，更新当前第k次迭代的模型：

式中，m代表模型参数，这里为介电常数，常量α_k代表迭代更新步长因子，Δm_k ^SD为下降方向；

l、反演序号加一；以当前模型作为下一次反演的初始模型，重复步骤e-k，直到反演序号结束，输出最终反演结果；

为了突出本发明：基于探地雷达数据的多尺度阶梯式层剥离全波形反演方法的效果，将本发明的反演结果(图8(c))与仅使用频率多尺度全波形反演方法以及采用传统反演序列的层剥离全波形反演方法的结果进行比较(图8(a)和图8(b))。与真实模型相比(图7(a))，可以观察到上述三种方法均能清楚地反映模型的整体异常分布特征。但是，由仅使用频率多尺度全波形反演方法的结果(图8(a))，不难发现模型中间高值区域的下边界反演得不准确，有明显的向下偏移，模型底部渐变区的裂缝位置也发生偏移。此外，从采用传统反演次序的层剥离全波形反演方法的结果来看(图8(b))，模型整体反演结果不稳定，深部异常反演偏差较大，特别是深部裂缝位置反演得不清晰，且出现很明显的分层界面。相比之下，本发明提出的多尺度阶梯式层剥离全波形反演方法(图8(c))，其反演结果的整体异常分布均匀，数据连续性好，特别是异常过渡带可以准确的表征。此外，模型中部区块异常的局部扰动异常和内部不均匀分布也反演的清晰准确。最重要的是，异常带的各个边界以及裂缝从浅到深的位置的反演精度也很高；

为了验证本发明结果的可靠性，我们从上述三种不同的反演策略中各在同一位置选择一个单道数据，并与真实模型数据进行数据拟合(图9(a))。可以看出，多尺度阶梯式层剥离全波形反演的介电常数的振幅在三者中与真模型最接近。此外，选择固定深度h＝8.46m的数据，得到三种不同反演策略的一维深度反演剖面(图9(b))，不难发现多尺度阶梯层剥离FWI的反演结果也更接近于真实模型。

Claims (1)

1.一种探地雷达数据的多尺度阶梯式层剥离全波形反演方法，其特征在于，包括以下步骤：

a、输入作为观测数据的地面雷达记录、初始模型以及子波；

b、选取模型分层数目、反演频率以及单频单层的反演次数；

c、依据多尺度阶梯式层剥离全波形反演的基本原理，确定各层各频率的反演序号；

d、设定初始反演序号为1；

e、确定当前第k次反演所对应的层号与频率(k＝1,2,3……,N)；

f、构造维纳滤波器，用于高频观测数据与高频子波向低频的转换：

式中，f_Wiener代表维纳低通滤波器，W_target为低频目标震源子波，W_original为原始震源子波，ε为控制数值溢出因子；

g、利用步骤f中的滤波器计算低频观测数据与低频模拟数据；

h、参考一阶伴随状态法，推导并采用TM模式下的介电常数梯度公式，计算整个模型空间的梯度值；

式中E_Z代表前向波场，E_Z*代表伴随波场，σ为电导率(单位S/m)，S为震源项，nt则为时间域的采样点；

i、对步骤h中计算出的梯度添加对应层的空间窗，获取当前层的梯度值；

j、利用最速下降法计算模型更新方向；

k、利用步骤i的梯度与步骤j的更新方向，并采用非精确线搜索方法计算更新步长，更新当前模型；

l、反演序号加一，以当前模型作为下一次反演的初始模型，重复步骤e-k，直到反演序号结束，输出最终反演结果。