CN109655910B - 基于相位校正的探地雷达双参数全波形反演方法 - Google Patents

Abstract

本发明涉及一种基于相位校正的探地雷达双参数全波形反演方法。是在现有探地雷达双参数反演基础上，针对探地雷达数据对介电常数和电导率的敏感度差异及二者在反演过程中所存在的耦合效应，从探地雷达记录层面出发，依据损耗正切公式进行相位校正，有效解决了反演过程中双参数梯度不一致的问题。对TE模式下多参数的同步重建进行了针对性改进，显著提高了探地雷达双参数全波形反演结果的分辨率和精确度，为双参数反演过程中进行记录层面的有效解耦提供了可能。对于具有相同结构的介电常数模型和电导率模型，经过对应记录的相位校正后所计算得到的梯度与正确梯度之间基本不存在相位差。

Description

基于相位校正的探地雷达双参数全波形反演方法

技术领域：

本发明涉及一种电磁领域中用于同时反演介电常数和电导率参数的时间域全波形成像方法。针对探地雷达数据对介电常数和电导率的敏感度差异及二者在反演过程中所存在的折中效应，依据损耗正切公式对电磁波记录进行相位校正，有效降低全波形反演中由双参数之间的耦合所带来的影响，从而提高了反演的精确度。

背景技术：

全波形反演(Full Waveform Inversion,FWI)作为一种高分辨率的反演成像方法，是近几十年来勘探地球物理界的研究热点。上世纪80年代，Lailly(1983.The seismicinverse problem as a sequence of before stack migrations.CONFERENCE ONINVERSE SCATTERING,THEORY AND APPLICATION,SOCIETY FOR INDUSTRIAL AND APPLIEDMATHEMATICS,EXPANDED ABSTRACT.206-220)和Tarantola(1984.Inversion of seismicreflection data in the acoustic approximation.GEOPHYSICS 49(8).1259-1266)首次提出时间域全波形反演策略。他们将Claerbout的偏移成像理论发展成为一种最小二乘局部优化问题，利用正向传播波场与后向传播波场的互相关计算梯度方向，避免了Frechet导数的直接计算。随后全波形反演在大尺度的地震勘探和浅层地震工程中均得到了广泛的发展和应用。近年来，国内外的专家学者始终致力于将地震全波形反演的核心技术引入到探地雷达(Ground Penetrating Radar,GPR)成像研究领域。在跨孔雷达全波形反演方面，Ernst等人(2007.Full-waveform inversion of crosshole radar data based on 2-Dfinite-difference time-domain solutions of Maxwell's equations.IEEETRANSACTION ON GEOSCIENCE AND REMOTE SENSING 45(9).2807-2828)采用级联更新的方式交替反演介电常数和电导率；Meles等人(2010.A new vector waveform inversionalgorithm for simultaneous updating of conductivity and permittivityparameters from combination crosshole/borehole-to-surface GPR data.IEEETRANSACTION ON GEOSCIENCE AND REMOTE SENSING 48(9).3391-3407)提出了一种新颖的矢量全波形反演技术，在反演过程中实现了介电常数和电导率的同时迭代更新。在地面探地雷达全波形反演研究方面，Lavou′e等人(2013.2D full waveform inversion of GPRsurface data:permittivity and conductivity imaging.7^th INTERNATIONAL WORKSHOPON ADVANCED GROUND PENETRATING RADAR.)提出了一种基于地面多偏移距数据的频率域全波形反演方法，用于同时重建二维介电常数和电导率参数。冯晅等人(2017.Jointacoustic full-waveform inversion of crosshole seismic and ground-penetratingradar in the frequency domain.GEOPHYSICS 82(6).1-16)开展了基于探地雷达数据和地震数据的频率域交叉梯度联合全波形反演研究，实现了P波速度、介电常数和电导率的顺次更新。

与跨孔雷达相比，地面雷达在地下深部的覆盖范围有限，从而增加了反演问题的不适定性。此外，地面雷达测量以较小的反射角照明地下目标，可以提供更高分辨率的图像。然而，由于低波数的缺乏使得对精确的初始模型具有高度依赖性。降低的照明度(单边照明模式)增加了反演过程中介电常数和电导率这两种参数之间的权衡与折中(tradeoff)，使得多参数同步反演成像更具挑战性。此外，探地雷达数据表现出对介电常数和电导率具有不同灵敏度的特征。根据Lavou′e等人(2014.Two-dimensional permittivity andconductivity imaging by full waveform inversion of multioffset GPR data:afrequency-domain quasi-Newton approach.Geophysical Journal International197.248-268)在频率域的散射场分析结果，首先，在相同频率和相同扰动量的条件下，介电常数对散射场振幅变化的影响要远大于电导率的影响，因此，探地雷达数据本质上对介电常数比对电导率更敏感。一般而言，两参数振幅比值取决于损耗角的正切值及各自相对扰动量的比值。其次，数据对电导率的灵敏度随频率的增加而降低，因此，关于低电导率差异的高频信息可能隐藏在噪声水平之下。最后，另一个重要特征是介电常数的散射场和电导率的散射场之间存在90度的相位差。

针对开展全波形反演过程中，不同电参数之间所存在的量级和灵敏度差异，Meles等人(2010.A new vector waveform inversion algorithm for simultaneous updatingof conductivity and permittivity parameters from combination crosshole/borehole-to-surface GPR data.IEEE TRANSACTION ON GEOSCIENCE AND REMOTESENSING 48(9).3391-3407)通过在其算法中引入两个下降步长来处理这个问题。但是，这相当于将关于介电常数和电导率的优化分割为两个相对独立的问题，并且忽略了两种参数之间可能存在的折中。Lavou′e等人(2014.Two-dimensional permittivity andconductivity imaging by full waveform inversion of multioffset GPR data:afrequency-domain quasi-Newton approach.Geophysical Journal International197.248-268)考虑了两参数同时反演过程中存在的折中关系，并采用准牛顿优化策略，用L-BFGS来近似Hassen矩阵的逆。

然而，目前的探地雷达双参数反演研究中，在求取目标函数对模型参数的梯度时，都将观测记录与模拟记录的差值认为是受单一参数影响而得到的，这会导致得到的梯度与真实的梯度之间存在相位差。

发明内容：

本发明的目的就在于针对上述现有技术的缺陷，提供了一种基于相位校正的全波形反演方法，用于同步反演具有相同结构的介电常数与电导率参数。在常规探地雷达的双参数反演基础上，针对探地雷达数据对介电常数和电导率的敏感度差异及二者在反演过程中所存在的耦合效应，从探地雷达记录层面出发，依据损耗正切公式进行相位校正，有效解决了反演过程中双参数梯度不一致的问题，从而提高了反演结果的分辨率和精确度，为双参数反演过程中进行记录层面的有效解耦提供了可能。

本发明的目的是通过以下技术方案实现的：首先，设定双参数同步反演的反演频率以及每个频率的迭代次数，分别计算观测记录与只存在介电常数扰动的观测记录、以及与只存在电导率扰动的观测记录之间的相位差，并利用这两个相位差对观测记录进行两次相位校正，从而分别得到校正后的介电常数和电导率对应的观测记录；其次，利用当前的介电常数模型和电导率模型，进行电磁波TE模式下的正演模拟，以获得模拟记录，并利用上述两个相位差对模拟记录进行两次相位校正，从而分别得到校正后的介电常数和电导率对应的模拟记录；然后，计算目标函数值，并利用一阶伴随状态法求取目标函数对介电常数和电导率的梯度；最后，利用最速下降法与非精确线搜索确定模型更新的方向与步长，更新当前模型，对每一个频率的每一次迭代都重复上述操作，最终输出介电常数和电导率模型的反演结果

基于相位校正的探地雷达双参数全波形反演方法，包括以下步骤

a、输入在地面接收的多偏移距探地雷达观测记录、介电常数和电导率的初始模型以及子波；

b、设定双参数同步反演的M个反演频率[f₁,f₂,…,f_M]以及每个频率的迭代次数N；

c、分别计算观测记录与只存在介电常数扰动的观测记录之间的相位差δ_ε，以及观测记录与只存在电导率扰动的观测记录之间的相位差δ_σ，具体计算方法如下：

首先，利用先验信息分别估算电导率和介电常数的相对扰动量比值:

式中，δp代表两个参数的相对扰动量比值，δp^ζ为电导率模型的扰动场与背景场的比值，δp^ε为介电常数模型的扰动场与背景场的比值；

然后，计算损耗正切公式的值：

式中，tanδ代表损耗正切，ζ为介质的电导率参数，ε为介质的介电常数参数，ω＝2πf为角频率，f为反演频率；

接下来，根据计算得到的损耗正切值，计算出观测记录与只存在介电常数扰动的观测记录之间的相位差：

式中δ_ε即为观测记录与只存在介电常数扰动的观测记录之间的相位差；

最后，由于只存在介电常数扰动的记录与只存在电导率扰动的记录之间存在90度的相位差，因此，根据δ_ε可以进一步计算出观测记录与只存在电导率扰动的观测记录之间的相位差：

式中δ_σ即为观测记录与只存在电导率扰动的观测记录之间的相位差；

d、利用步骤c中计算的两个相位差δ_ε和δ_σ分别对观测记录各进行一次相位校正，即：利用步骤c中计算的相位差δ_ε，借助希尔伯特变换先对观测记录进行一次相位校正，得到与只存在介电常数扰动的观测记录相位一致的观测记录；利用步骤c中计算的相位差δ_σ，借助希尔伯特变换再对观测记录进行一次相位校正，得到与只存在电导率扰动的观测记录相位一致的观测记录；

e、利用当前的介电常数模型和电导率模型，进行电磁波TE模式的正演模拟，以获得模拟记录；

f、利用步骤c中计算的两个相位差δ_ε和δ_σ分别对模拟记录各进行一次相位校正，即：利用步骤c中计算的相位差δ_ε，借助希尔伯特变换先对模拟记录进行一次相位校正，得到与只存在介电常数扰动的模拟记录相位一致的模拟记录；利用步骤c中计算的相位差δ_σ，借助希尔伯特变换再对模拟记录进行一次相位校正，得到与只存在电导率扰动的模拟记录相位一致的模拟记录；

g、计算目标函数值，并将步骤d中经过相位校正后得到的“与只存在介电常数扰动的观测记录相位一致的观测记录”和“与只存在电导率扰动的观测记录相位一致的观测记录”，以及步骤f中经过相位校正后得到的“与只存在介电常数扰动的模拟记录相位一致的模拟记录”和“与只存在电导率扰动的模拟记录相位一致的模拟记录”共四个记录代入到伴随状态法中，求取目标函数值对介电常数和电导率的梯度；

h、利用最速下降法与非精确线搜索确定模型更新的方向与步长，并更新模型；

i、对每一个频率的每一次迭代都进行步骤c～h的处理，最终输出介电常数和电导率模型的反演结果。

有益效果：

本发明将时间域全波形反演技术引入到地面多偏移距探地雷达数据的成像当中，在雷达数据记录层面上，针对介电常数和电导率参数在反演过程中所存在的耦合效应进行解耦，对TE模式下多参数的同步重建进行了针对性改进，显著提高了利用基于相位校正的探地雷达双参数全波形反演方法进行成像在反演精度和分辨率两方面的效果，具体有以下优点：①探地雷达数据对介电常数和电导率的敏感度不同，且这种敏感度表现出较强的频率依赖性，由于两参数之间耦合折中效应的存在，使得在反演的不同阶段，两参数对记录的影响程度也不断变化，因此，从能否解耦的角度，在记录层面上进行校正大大提高了参数解耦的可能性，且解耦的计算难度也有所降低；②从解耦理论依据的角度来看，根据损耗正切公式对观测记录和模拟记录进行相位校正，较之根据两参数梯度公式所存在的π/2相位差进行简单的梯度相位校正更具合理性；③从解耦效果的角度来看，对于具有相同结构的介电常数模型和电导率模型，经过对应记录的相位校正后所计算得到的梯度与正确梯度之间基本不存在相位差；④从双参数同步反演效果的角度来看，利用该方法反演得到的背景场和扰动场的分辨率和精确度都有大幅度提升。

附图说明：

图1基于相位校正的探地雷达双参数全波形反演方法流程图。

图2线状扰动模型。

(a)介电常数真实模型；

(b)电导率真实模型。

图3线状扰动模型的初始模型。

(a)介电常数初始模型；

(b)电导率初始模型。

图4线状扰动模型在反演过程中某一次迭代的当前模型。

(a)介电常数当前模型；

(b)电导率当前模型。

图5对线状扰动模型的观测记录(实线)进行角度为δ_ε的相位校正(点线)，与只存在介电常数扰动的理论观测记录(虚线)进行对比，只展示了记录的某一道。

图6对线状扰动模型的观测记录(实线)进行角度为δ_σ的相位校正(点线)，与只存在电导率扰动的理论观测记录(虚线)进行对比，只展示了记录的某一道。

图7对当前线状扰动模型的模拟记录(实线)进行角度为δ_ε的相位校正(点线)，与只存在介电常数扰动的理论模拟记录(虚线)进行对比，只展示了记录的某一道。

图8对当前线状扰动模型的模拟记录(实线)进行角度为δ_σ的相位校正(点线)，与只存在电导率扰动的理论模拟记录(虚线)进行对比，只展示了记录的某一道。

图9不进行相位校正所得到的介电常数梯度。

图10不进行相位校正所得到的电导率梯度。

图11进行相位校正后所得到的介电常数梯度。

图12进行相位校正后所得到的电导率梯度。

图13介电常数的理论梯度。

图14电导率的理论梯度。

图15实例1所使用的真实模型。

(a)介电常数真实模型；

(b)电导率真实模型。

图16实例1所使用的初始模型。

(a)介电常数初始模型；

(b)电导率初始模型。

图17实例1得到的基于相位校正的反演结果。

(a)介电常数反演结果；

(b)电导率反演结果。

图18实例1得到的常规反演结果。

(a)介电常数反演结果；

(b)电导率反演结果。

图19实例1的真实模型、基于相位校正的反演结果、常规反演结果在不同位置处的单道对比。

(a)介电常数单道对比(x＝1.44m)；

(b)介电常数单道对比(x＝4.5m)；

(c)介电常数单道对比(x＝8.64m)；

(d)电导率单道对比(x＝1.44m)；

(e)电导率单道对比(x＝4.5m)；

(f)电导率单道对比(x＝8.64m)。

图20实例2所使用的真实模型。

(a)介电常数真实模型；

(b)电导率真实模型。

图21实例2所使用的初始模型。

(a)介电常数初始模型；

(b)电导率初始模型。

图22实例2得到的基于相位校正的反演结果。

(a)介电常数反演结果；

(b)电导率反演结果。

图23实例2得到的常规反演结果。

(a)介电常数反演结果；

(b)电导率反演结果。

图24实例2的真实模型、基于相位校正的反演结果、常规反演结果的单道对比。

(a)介电常数单道对比；

(b)电导率单道对比。

具体实施方式：

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

a、输入在地面接收的多偏移距探地雷达观测记录、介电常数和电导率的初始模型以及子波。在这些输入数据中，介电常数初始模型、电导率初始模型、子波都需要根据先验信息来确定。这三个参数中，介电常数和电导率是反演的目标参数而子波不是。更准确的初始模型和子波有助于获得更好的反演结果。

b、设定双参数同步反演的M个反演频率[f₁,f₂,…,f_M]以及每个频率的迭代次数N。M个频率的数值以及N的数值应该根据目标区域的大小、目标区域模型参数值大小、目标区域复杂程度等实际情况来综合确定。不同频率的观测记录是通过对原始观测记录进行维纳滤波来提取的。

c、分别计算观测记录与只存在介电常数扰动的观测记录之间的相位差δ_ε，以及观测记录与只存在电导率扰动的观测记录之间的相位差δ_σ，具体计算方法如下：

首先，利用先验信息分别估算电导率和介电常数的相对扰动量比值:

式中，δp代表两个参数的相对扰动量比值，δp^σ为电导率模型的扰动场与背景场的比值，δp^ε为介电常数模型的扰动场与背景场的比值。当真实电导率模型与真实介电常数模型具有相同的结构，并假设二者的初始模型与真实模型之间的偏差程度大致相同时，可以简单令该比值等于1。

然后，计算损耗正切公式的值：

式中，tanδ代表损耗正切，在探地雷达勘查中，介质一般呈现低损耗特征，即tanδ<<1，σ为介质的电导率参数，ε为介质的介电常数参数，ω＝2πf为角频率，f为反演频率。

接下来，根据计算得到的损耗正切，计算出观测记录与只存在介电常数扰动的观测记录之间的相位差：

式中δ_ε即为观测记录与只存在介电常数扰动的观测记录之间的相位差；

最后，由于只存在介电常数扰动的记录与只存在电导率扰动的记录之间存在90度相位差，因此，根据δ_ε可以计算出观测记录与只存在电导率扰动的观测记录之间的相位差：

式中δ_σ即为观测记录与只存在电导率扰动的观测记录之间的相位差。

d、利用步骤c中计算的两个相位差δ_ε和δ_σ分别对观测记录各进行一次相位校正，例如：通过如图2所示的介电常数和电导率模型生成的观测记录，根据δ_ε，借助希尔伯特变换先对观测记录进行角度为δ_ε的一次相位校正，校正后的观测记录与只存在介电常数扰动的理论观测记录相位接近相同(如图5所示)；根据δ_σ，借助希尔伯特变换再对观测记录进行角度为δ_σ的一次相位校正，校正后的观测记录与只存在电导率扰动的理论观测记录相位接近相同(如图6所示)。

e、利用如图4所示的当前介电常数模型和电导率模型，并结合子波信息，进行电磁波TE模式的正演模拟，以获得模拟记录。

f、利用步骤c中计算的两个相位差δ_ε和δ_σ分别对模拟记录各进行一次相位校正，例如：根据δ_ε，借助希尔伯特变换先对步骤e中的模拟记录进行角度为δ_ε的一次相位校正，校正后的模拟记录与只存在介电常数扰动的理论模拟记录相位接近相同(如图7所示)；根据δ_σ，借助希尔伯特变换再对步骤e中的模拟记录进行角度为δ_σ的一次相位校正，校正后的模拟记录与只存在电导率扰动的理论模拟记录相位接近相同(如图8所示)。

g、计算目标函数值，并求取目标函数分别对介电常数和电导率的梯度。

设定目标函数为观测记录与模拟记录之差的二范数的平方，则目标函数值的计算公式为：

式中,m表示模型参数,d_obs为观测记录，d_sim为模拟记录,表示计算二范数。

梯度计算采用伴随状态法，首先需要计算伴随源：将步骤d中经δ_ε校正后的观测记录与步骤f中经δ_ε校正后的模拟记录相减，作为介电常数的伴随源；将步骤d中经δ_σ校正后的观测记录与步骤f中经δ_σ校正后的模拟记录相减，作为电导率的伴随源。

在接收天线位置处分别反向激发介电常数的伴随源和电导率的伴随源，分别得到介电常数和电导率对应的伴随波场。采用电磁波TE模式下介电常数和电导率的梯度公式，计算目标函数对介电常数以及电导率的梯度值：

式中，E_y代表正传波场，E_y*代表反向传播的伴随波场，s为发射天线，r为接收天线，T为记录时长。

将两个梯度合并，即为总梯度：

式中，表示梯度算子。

例如：以图2为真实模型，图3为初始模型，在如图4所示的当前模型的基础上，进行下一次迭代；若不对观测记录和模拟记录分别进行步骤d和f的校正，则得到的介电常数梯度和电导率的梯度分别如图9和图10所示；对观测记录和模拟记录分别进行步骤d和f的校正，得到的介电常数梯度和电导率的梯度分别如图11和图12所示；分别利用只存在介电常数扰动和电导率扰动的模型计算得到的介电常数和电导率的理论梯度分别如图13和图14所示。可以观察到：经过步骤d和f的相位校正后所获得的梯度与理论梯度的变化形态和异常分布位置都高度吻合。

h、利用最速下降法与非精确线搜索确定模型更新的方向与步长，并更新模型。

最速下降法的更新方向为目标函数梯度的反方向：

式中，Δm^SD代表最速下降方向。

采用非精确线搜索方法计算更新步长α，是指选取步长α使目标函数C(m)得到可接受的下降量。即确定α使得C(m+αΔm)＜C(m)。

最后，同时更新介电常数和电导率模型：

m←m+αΔm^SD

该式表示将m+αΔm^SD作为更新后的模型赋值给m。

i、对每一个频率的每一次迭代都进行步骤c-h的处理，顺序为先在同一频率下进行N次迭代，再从低到高遍历所有频率，得到最终的介电常数和电导率模型输出。

实施例1：

a、输入在地面接收的多偏移距探地雷达观测记录、介电常数和电导率的初始模型以及子波。

采用简单的块体模型作为真实模型(如图15所示)，该模型内部含有三个尺度不同的块体，块体异常的相对介电常数为27，电导率为0.009；背景介质的相对介电常数为9，电导率为0.003。模型大小为9.54m*4.95m，在模型的顶部有h_a＝0.36m的空气层；有限差分网格间距h＝9cm，故该模型共计N_M＝107×60＝6420个网格点；地面观测系统由22个发射天线，间距5个网格点，54个接收天线，间距2个网格点组成，且每个天线发射的信号可被全部接收天线记录；将中心频率为120MHz的高斯子波作为合成记录和反演的子波，即认为子波已知；输入由真实模型计算的合成数据作为的观测数据；输入均匀背景模型作为初始模型(如图16所示)。

b、选取6个反演频率，分别为：10.0，19.3，31.4，45.8，62.0，80.0，单位为MHz，每个频率的迭代次数为10，用于完成双参数模型同步反演。

c、设定δp＝δp^σ/δp^ε＝1，计算观测记录与只存在介电常数扰动的观测记录之间的相位差δ_ε，以及观测记录与只存在电导率扰动的观测记录之间的相位差δ_σ：

d、利用步骤c中计算的两个相位差δ_ε和δ_σ分别对观测记录各进行一次相位校正，即：根据δ_ε，借助希尔伯特变换先对观测记录进行一次相位校正，得到与只存在介电常数扰动的观测记录相位一致的观测记录；根据δ_σ，借助希尔伯特变换再对观测记录进行一次相位校正，得到与只存在电导率扰动的观测记录相位一致的观测记录；

e、利用当前的介电常数模型和电导率模型，并结合子波信息，进行电磁波TE模式的正演模拟，以获得模拟记录。

f、利用步骤c中计算的两个相位差δ_ε和δ_σ分别对模拟记录各进行一次相位校正，即：根据δ_ε，借助希尔伯特变换先对模拟记录进行一次相位校正，得到与只存在介电常数扰动的模拟记录相位一致的模拟记录；根据δ_σ，借助希尔伯特变换再对模拟记录进行一次相位校正，得到与只存在电导率扰动的模拟记录相位一致的模拟记录；

g、计算目标函数值，并求取目标函数值对介电常数和电导率的梯度。

目标函数值的计算公式为：

式中,m表示模型参数,d_obs为观测记录，d_sim为模拟记录,表示计算二范数。

采用电磁波TE模式下介电常数和电导率的梯度公式，计算目标函数对介电常数以及电导率的梯度值：

式中，E_y代表正传波场，E_y*代表反向传播的伴随波场，s为发射天线，r为接收天线，T为记录时长。

h、利用最速下降法与非精确线搜索确定模型更新的方向与步长，并更新模型。

最速下降法的更新方向为目标函数梯度的反方向：

式中，Δm^SD代表最速下降方向。

采用非精确线搜索方法，选取步长α，使得C(m+αΔm)＜C(m)

最后，同时更新介电常数和电导率模型：

m←m+αΔm^SD

该式表示将m+αΔm^SD作为更新后的模型赋值给m。

i、对每一个频率的每一次迭代都进行步骤c-h的处理，顺序为先在同一频率下进行10次迭代，再从低到高遍历所有频率，得到最终的介电常数和电导率模型基于相位校正的反演结果(如图17所示)。

在完全相同的模型结构和反演参数设置下，若不进行步骤d和f的相位校正，则得到的介电常数和电导率模型的常规反演结果如图18所示。

可以看出，基于相位校正的双参数反演结果，三个块状目标体的基本形态和分布位置在两个参数的反演结果中均能被准确重建，反演结果与真实模型吻合度较高。相比之下，常规反演的结果(如图18所示)，从整体上看，介电常数反演的结果优于电导率，具体表现在异常体的介电常数反演结果在形状和数值上都与真实模型更接近，尽管如此，常规介电常数反演结果在数值上与真实模型仍然存在一定差异。此外，从常规反演的电导率结果可以看出，三个块体只有部分结构得到了恢复，且被恢复的异常位置也偏离了其各自的真实位置。整体的反演准确性较差，置信度较低。

分别对介电常数和电导率的真实模型(如图15所示)、基于相位校正的反演结果(如图17所示)、常规反演结果(如图18所示)抽取不同块体异常位置处的单道数据进行对比，如图19所示。从介电常数单道抽取对比中(如图19(a)-(c)所示)，不难发现，在三个不同的异常体位置处，基于相位校正的方法较常规反演方法更接近于真实模型。在电导率单道抽取对比中(如图19(d)-(f)所示)，常规反演结果与真实模型的位置偏离较远，且异常数值也与真实模型相差较大，相比之下，基于相位校正的反演结果，与真实模型在数值和位置方面均具有较好的拟合效果。

实施例2：

a、输入在地面接收的多偏移距探地雷达观测记录、介电常数和电导率的初始模型以及子波。

利用多参数耦合的随机土壤介质模型作为真实模型(如图20所示)，模型大小为9.54m*9m，在模型顶部具有h_a＝0.36m的空气层；有限差分网格间距h＝9cm，故该模型共计N_M＝107×105＝11235个网格点；地面观测系统由36个发射天线，间距3个网格点，107个接收天线，间距1个网格点组成，且每个天线发射的信号可被全部接收天线记录；将中心频率为120MHz的高斯子波作为合成记录的子波，即认为子波已知；输入由真实模型计算的合成数据作为的观测数据；输入真实模型的大尺度平滑结果作为初始模型(如图21所示)。

b、选取15个反演频率，分别为：10.0，12.9，16.6，21.0，26.0，31.6，37.6，44.1，51.0，58.3，66.0，74.0，82.3，91.0，100，单位为MHz，每个频率的迭代次数为40，用于完成双参数模型同步反演。

c、设定δp＝δp^σ/δp^ε＝1，计算观测记录与只存在介电常数扰动的观测记录之间的相位差δ_ε，以及观测记录与只存在电导率扰动的观测记录之间的相位差δ_σ：

d、利用步骤c中计算的两个相位差δ_ε和δ_σ分别对观测记录各进行一次相位校正，即：根据δ_ε，借助希尔伯特变换先对观测记录进行一次相位校正，得到与只存在介电常数扰动的观测记录相位一致的观测记录；根据δ_σ，借助希尔伯特变换再对观测记录进行一次相位校正，得到与只存在电导率扰动的观测记录相位一致的观测记录；

e、利用当前的介电常数模型和电导率模型，并结合子波信息，进行电磁波TE模式的正演模拟，以获得模拟记录。

f、利用步骤c中计算的两个相位差δ_ε和δ_σ分别对模拟记录各进行一次相位校正，即：根据δ_ε，借助希尔伯特变换先对模拟记录进行一次相位校正，得到与只存在介电常数扰动的模拟记录相位一致的模拟记录；根据δ_σ，借助希尔伯特变换再对模拟记录进行一次相位校正，得到与只存在电导率扰动的模拟记录相位一致的模拟记录；

g、计算目标函数值，并求取目标函数值对介电常数和电导率的梯度。

目标函数值的计算公式为：

式中,m表示模型参数,d_obs为观测记录，d_sim为模拟记录,表示计算二范数。

采用电磁波TE模式下介电常数和电导率的梯度公式，计算目标函数对介电常数以及电导率的梯度值：

式中，E_y代表正传波场，E_y*代表反向传播的伴随波场，s为发射天线，r为接收天线，T为记录时长。

h、利用最速下降法与非精确线搜索确定模型更新的方向与步长，并更新模型。

最速下降法的更新方向为目标函数梯度的反方向：

式中，Δm^SD代表最速下降方向。

采用非精确线搜索方法，选取步长α，使得C(m+αΔm)＜C(m)

最后，同时更新介电常数和电导率模型：

m←m+αΔm^SD

该式表示将m+αΔm^SD作为更新后的模型赋值给m。

i、对每一个频率的每一次迭代都进行步骤c-h的处理，顺序为先在同一频率下进行40次迭代，再从低到高遍历所有频率，得到最终的介电常数和电导率模型基于相位校正的反演结果(如图22所示)。

在完全相同的模型结构和反演参数设置下，若不进行步骤d和f的相位校正，得到最终的介电常数和电导率模型的常规反演结果如图23所示。

可以观察到常规反演结果中，介电常数异常的浅层分布特征，即0m-6m的区域可以得到较清晰地反映(如图23(a)所示)，然而模型中7m-9m范围内近水平分布的裂缝及局部不均匀随机扰动无法被准确识别与刻画。此外，未经相位校正的电导率的常规反演结果在反演稳定性和异常刻画的精确程度方面都较差(如图23(b)所示)。例如，对模型中间高值区域的局部小异常以及内部介质分布的非均匀性的识别方面非常模糊，数据连续性较差，异常过渡带无法准确识别。此外，基本无法重建模型底部渐变区的裂缝，无法提供高精度的清晰反演结果。相比之下，基于相位校正的双参数同步全波形反演结果不仅对模型整体异常分布带有清晰的刻画，对于异常带的形状和异常边界均有有效识别。模型中部区块异常的局部扰动和内部不均匀分布也被反演的清晰准确。最重要的是，与双参数的真实模型相比，反演结果中从浅到深的裂缝异常位置也得到了精确地表征。此外，对于这种具有相同结构的介电常数和电导率的复杂模型，经过相位校正后的同步反演结果中，两参数相同结构位置的对应完全一致。

分别对介电常数和电导率的真实模型(如图20所示)、基于相位校正的反演结果(如图22所示)、常规反演结果(如图23所示)进行单道抽取作为对比，如图24所示。可以看出模型从浅到深的各个位置，基于相位校正方法得到的介电常数值相较于常规反演更接近于真实模型(如图24(a)所示)；这一点在电导率的单道数据对比结果中尤为明显(如图24(b)所示)，不难发现，由常规反演方法得到的电导率结果与真实模型的差异较大，特别是在模型的中部和底部，这种差异更为突出，相比之下，基于相位校正方法得到的电导率值在变化趋势及幅值振荡方面都与真实模型拟合得很好。

Claims (1)

1.一种基于相位校正的探地雷达双参数全波形反演方法，其特征在于，包括以下步骤：

a、输入在地面接收的多偏移距探地雷达观测记录、介电常数和电导率的初始模型以及子波；

b、设定双参数同步反演的M个反演频率[f₁,f₂,…,f_M]以及每个频率的迭代次数N；

c、分别计算观测记录与只存在介电常数扰动的观测记录之间的相位差δ_ε，以及观测记录与只存在电导率扰动的观测记录之间的相位差δ_σ，具体计算方法如下：

首先，利用先验信息分别估算电导率和介电常数的相对扰动量比值:

式中，δp代表两个参数的相对扰动量比值，δp^σ为电导率模型的扰动场与背景场的比值，δp^ε为介电常数模型的扰动场与背景场的比值；

然后，计算损耗正切公式的值：

式中，tanδ代表损耗正切，σ为介质的电导率参数，ε为介质的介电常数参数，ω＝2πf为角频率，f为反演频率；

接下来，根据计算得到的损耗正切值，计算出观测记录与只存在介电常数扰动的观测记录之间的相位差：

式中δ_ε即为观测记录与只存在介电常数扰动的观测记录之间的相位差；

最后，由于只存在介电常数扰动的记录与只存在电导率扰动的记录之间存在90度的相位差，因此，根据δ_ε可以进一步计算出观测记录与只存在电导率扰动的观测记录之间的相位差：

式中δ_σ即为观测记录与只存在电导率扰动的观测记录之间的相位差；

d、利用步骤c中计算的两个相位差δ_ε和δ_σ分别对观测记录各进行一次相位校正，即：利用步骤c中计算的相位差δ_ε，借助希尔伯特变换先对观测记录进行一次相位校正，得到与只存在介电常数扰动的观测记录相位一致的观测记录；利用步骤c中计算的相位差δ_σ，借助希尔伯特变换再对观测记录进行一次相位校正，得到与只存在电导率扰动的观测记录相位一致的观测记录；

e、利用当前的介电常数模型和电导率模型，进行电磁波TE模式的正演模拟，以获得模拟记录；

f、利用步骤c中计算的两个相位差δ_ε和δ_σ分别对模拟记录各进行一次相位校正，即：利用步骤c中计算的相位差δ_ε，借助希尔伯特变换先对模拟记录进行一次相位校正，得到与只存在介电常数扰动的模拟记录相位一致的模拟记录；利用步骤c中计算的相位差δ_σ，借助希尔伯特变换再对模拟记录进行一次相位校正，得到与只存在电导率扰动的模拟记录相位一致的模拟记录；

g、计算目标函数值，并将步骤d中经过相位校正后得到的“与只存在介电常数扰动的观测记录相位一致的观测记录”和“与只存在电导率扰动的观测记录相位一致的观测记录”，以及步骤f中经过相位校正后得到的“与只存在介电常数扰动的模拟记录相位一致的模拟记录”和“与只存在电导率扰动的模拟记录相位一致的模拟记录”共四个记录代入到伴随状态法中，求取目标函数值对介电常数和电导率的梯度；

h、利用最速下降法与非精确线搜索确定模型更新的方向与步长，并更新模型；

i、对每一个频率的每一次迭代都进行步骤c～h的处理，最终输出介电常数和电导率模型的反演结果。