一种基于反射信号重构的薄层介质参数反演方法

1.本发明涉及探地雷达信号处理技术领域，具体涉及对含有薄层的分层介质场景的参数反演。


背景技术：

2.探地雷达是一种无损检测工具，探地雷达使用高频电磁波来检测地下结构的特征。该技术已广泛应用于浅层地下地球物理勘探，如确定地下目标、含水量、层厚和物质密度等。层状介质厚度反演是探地雷达在实际工程中的一个重要应用。它根据不同材料之间的介电不连续性产生的反射信号来确定测量数据的厚度和介电参数。然而，受限于雷达本身的纵向分辨率，当层厚很薄时，反射信号会发生重叠，因此很难在时域上区分出薄层，尤其当厚度小于四分之一波长。
3.基于子空间分析的多信号分类(music)和独立成分分析(ica)的方法被用于实现超分辨。然而，由于多次迭代，这两种方法都很耗时。反卷积是时域或频域的常用方法，但由于计算复杂，不容易应用于实时处理。机器学习，如支持向量机回归(svr)，在预测薄层的时间延迟和介电常数方面有很好的结果，但需要大量数据来进行模型训练。利用梯度下降、牛顿迭代等优化算法进行反演的方法可以在实时探测中被应用，但是会因反演参数过多而使结果陷入局部最优解。
4.因此可以看出探地雷达薄层介质参数的反演中，由于雷达反射信号在薄层中产生混叠，传统的方法不能很好解决薄层的时延估计问题，并且一些复杂的算法会由于运算量过大而无法做到实时探测。


技术实现要素：

5.有鉴于此，本发明提供了一种基于反射信号重构的薄层介质参数的方法，适用于含有薄层的分层介质模型，能够解决薄层的时延估计问题，且该方法精度和收敛速度均较高，能够获得更好的薄层参数反演性能，是一种高效准确的反演方法。
6.为达到上述目的，本发明的技术方案为：一种基于反射信号重构的薄层介质参数反演方法，该方法用于在未知先验信息的情况下，反演地下薄层介质参数，先验信息包括探测场景材料和实际厚度，方法包括如下步骤：
7.首先通过广义反射系数的定义以及地下分层的褶积模型重构出反射信号模型，该模型包含了层状介质的物理参数。
8.然后构建优化问题，优化问题的代价函数为模型和实际反射信号之间的误差平方，通过遗传算法求解优化问题，使误差平方最小化；通过遗传算法求解优化问题在初值范围内进行优化。
9.进一步地，该方法适用于公路路面、桥梁内部分层结构以及其他包含薄层的分层介质探测场景，同时也适用于非薄层的介质探测。
10.进一步地，通过广义反射系数的定义以及地下分层的褶积模型重构出反射信号模
型，具体为：
11.地下各介质层都是线性、均匀、各向同性的，且为非磁性无损材料。
12.则第i层和第i 1层之间界面的反射系数r
i,i 1
和透射系数分别为τ
i,i 1
：
[0013][0014]
电磁波的衰减常数为：
[0015][0016]
其中：ω表示角频率；ε
ri
表示第i层的相对介电常数，ε
ri 1
表示相对介电常数；σi表示第i层的电导率；μi表示第i层的磁导率。
[0017]
定义广义反射系数为接收到的每层分界面的反射波振幅与发射波的比值，记为r
i,i 1
；由于第一层顶部上方为空气介质，因此有r
0,1
＝r
0,1
，r
0,1
表示第0层和第1层之间界面的反射系数；得第i个广义反射系数r
i,i 1
为：
[0018][0019]
式中dk为第k层的层厚；k取值为1～i；αk为第k层的电磁波的衰减常数；r
k,k 1
表示第k层和第k 1层之间界面的反射系数。
[0020]
层状介质的反射系统简化为线性时不变lti系统，则接收信号yr(t)为：
[0021][0022]
式中，x(t)为t时刻的雷达发射信号，x(t-ti)为t-ti时刻的雷达发射信号；n为总分层个数；n(t)表示高斯白噪声。
[0023]
进一步地，构建优化问题，优化问题的代价函数为模型和实际反射信号之间的误差平方，通过遗传算法求解优化问题，使误差平方最小化，具体为：
[0024]
最优化问题为：
[0025]
其中，o(l)为包含参数的代价函数，l为待优化的参数矩阵；o(l)表达式为：
[0026][0027]
其中：yr(t)为模型信号；y(t)为实际反射信号；||||2为l2范数。
[0028]
优化问题的物理意义为：求取最优参数，使得重构的反射信号与实际接收的反射信号拟合度最高。
[0029]
进一步地，通过遗传算法求解优化问题在初值范围内进行优化，具体为：
[0030]
将介质层双程发射时间ti与电性参数εi和σi作为优化对象，在给定初值范围的条件下得到参数最优解；利用广义反射系数频谱估计参数，得到反演初值：通过对广义反射系数序列的频谱进行分析，初步估计出待反演参数的值，并以此作为遗传算法的输入参数初值，在采用遗传算法求解优化问题在初值范围内进行优化。
[0031]
其中分层介质第一层底面和第二层底面的反射系数序列为r(t)：
[0032]
r(t)＝r
12
δ(t-t1) r
23
δ(t-t2)
[0033]
其中r
12
和r
23
为广义反射系数，δ(t-t1)和δ(t-t2)分别为t1，t2时刻的冲激响应。
[0034]
对反射系数系列作傅里叶变换即得到反射系数序列谱r(f)：
[0035][0036]
其中r(f)为时间零点在中间层的中点位置处的广义反射系数频谱；j为复数符号。
[0037]
将反射系数序列谱r(f)分为幅度谱|r(f)|和相位谱θ(f)：
[0038][0039]
θ(f)＝arctan[tan(πfδt)ro/re]-2πft
1-πfδt
[0040]
其中：t1为第一层介质的电磁波双程反射时延，δt＝t
1-t0为第二层介质的电磁波双程反射时延，re和ro为广义反射系数对进行奇偶分解的结果：
[0041][0042]
对幅度谱|r(f)|进行求导并令其等于0：
[0043][0044]
解得频率为n为正整数，n＝1,2,...。
[0045]
f为幅度谱极大值点和极小值点所对应的频率，且为周期性分布，即相邻的两个极大/小值的间隔为幅度谱的凹陷周期δf：
[0046][0047]
利用时间厚度与幅度谱凹陷周期之间的反比关系，利用反射系数频谱的周期来估计出时间厚度。
[0048]
θ(f)分为两个部分，第一部分包含参数δt和反射系数奇偶分量，影响相位谱的振幅大小；第二部分包含参数δt和t1，影响相位谱的变化速度；反射系数奇偶分量包含介电常数和与电导率相关的衰减常数，则通过相位谱估计出三层介质的介电常数和中间层的电导率。
[0049]
得到的反演初值包括上述估计的时间厚度、三层介质的介电常数和中间层的电导率。
[0050]
有益效果：
[0051]
1、本发明提出了一种基于反射信号重构的薄层介质参数的方法。本发明针对薄层
反射波在时域无法区分的问题，首先利用包含各层的介电常数和双向传播时间的广义反射系数构建了探地雷达反射信号的正演模型，然后将参数反演转化为优化问题，定义上述信号模型和实际反射信号之间的误差平方为代价函数，通过遗传算法使误差最小化。同时，针对解决优化算法在多参数反演情况下容易产生局部最优解的问题，通过分析广义反射系数谱来初步求解待反演的参数，作为算法的初始值，从而提高算法的精度和收敛速度。
附图说明
[0052]
图1是本发明实施方式的信号处理流程图；
[0053]
图2是广义反射系数模型示意图；
[0054]
图3是幅度谱凹陷周期与时间厚度t的关系示意图；
[0055]
图4是仿真模型示意图。
具体实施方式
[0056]
下面结合附图并举实施例，对本发明进行详细描述。
[0057]
本发明目的在于克服探地雷达在薄层介质参数反演方法的不足，提供了一种基于反射信重构的参数反演方法。本发明适用于含有薄层的分层介质模型，为工程实践提供了高效精确的方法来检测路面、桥梁等工程质量问题。
[0058]
图1是本发明实施方式的信号处理流程图。本发明通过以下步骤实现：
[0059]
步骤1：建立探地雷达接收反射信号模型；
[0060]
假设地下各层都是线性、均匀、各向同性的，且为非磁性无损材料。则第i层和第i 1层之间界面的反射系数r
i,i 1
和透射系数分别为τ
i,i 1
：
[0061][0062]
电磁波的衰减常数可写为：
[0063][0064]
其中：ω表示角频率；ε
ri
表示第i层的相对介电常数，ε
ri 1
表示相对介电常数；σi表示第i层的电导率；μi表示第i层的磁导率。
[0065]
定义广义反射系数为接收到的每层分界面的反射波振幅与发射波的比值，如图2所示，记为r
i,i 1
。由于第一层顶部上方为空气介质，因此有r
01
＝r
01
。经推导可得第i个广义反射系数r
i,i 1
为：
[0066][0067]
式中dk为第k层的层厚；k取值为1～i；αk为第k层的电磁波的衰减常数；rk,
k 1
表示第k层和第k 1层之间界面的反射系数。
[0068]
将层状介质的反射系统通常简化为线性时不变(lti)系统，那么接收信号yr(t)为：
[0069][0070]
式中，式中，x(t)为t时刻的雷达发射信号，x(t-ti)为t-ti时刻的雷达发射信号；n为总分层个数；n(t)表示高斯白噪声。
[0071]
步骤2：构造代价函数，优化求解参数；
[0072]
将参数反演问题转化为以下最优化问题：
[0073][0074]
其中，o(l)为包含参数的代价函数，l为待优化的参数矩阵，参数矩阵中包含[ε
r1
，ε
r2
ε
r3
，σ1，σ2，d1,d2]，即包含三层介质的介电常数，两层的电导率和两层的层厚；o(l)其表达式为：
[0075][0076]
||||2为l2范数。该优化问题的物理意义为：求取最优参数，使得重构的反射信号能够与实际接收的反射信号拟合度最高。
[0077]
利用遗传算法对上述优化问题进行求解，将介质层双程发射时间ti与电性参数εi和σi作为优化对象，在给定初值范围的条件下得到参数最优解。
[0078]
步骤3：利用广义反射系数频谱估计参数，得到反演初值；
[0079]
将介质层双程发射时间ti与电性参数εi和σi作为优化对象，在给定初值范围的条件下得到参数最优解；利用广义反射系数频谱估计参数，得到反演初值：通过对广义反射系数序列的频谱进行分析，初步估计出待反演参数的值，并以此作为遗传算法的输入参数初值，在采用遗传算法求解优化问题在初值范围内进行优化。
[0080]
通过对广义反射系数序列的频谱进行分析，可以初步估计出待反演参数的值，并以此作为优化算法的输入参数初值。以三层介质模型为例。
[0081]
分层介质第一层底面和第二层底面的反射系数序列可表示为：
[0082]
r(t)＝r
12
δ(t-t1) r
23
δ(t-t2)
[0083]
其中r
12
和r
23
为广义反射系数，δ(t-t1)和δ(t-t2)分别为t1，t2时刻的冲激响应。
[0084]
对反射系数系列作傅里叶变换即可得到反射系数序列谱：
[0085][0086]
其中r(f)为时间零点在中间层的中点位置处的广义反射系数频谱；j为复数符号。
[0087]
将反射系数序列谱r(f)分为幅度谱|r(f)|和相位谱θ(f)：
[0088][0089]
θ(f)＝arctan[tan(πfδt)ro/re]-2πft
1-πfδt
[0090]
其中：t1为第一层介质的电磁波双程反射时延，δt＝t
1-t0为第二层介质的电磁波双程反射时延，re和ro为广义反射系数对进行奇偶分解的结果：
[0091][0092]
进行求导并令其等于0：
[0093][0094]
解得：
[0095][0096]
根据导数的定义可知，f为幅度谱极大值点和极小值点所对应的频率，且为周期性分布，即相邻的两个极大/小值的间隔为幅度谱的凹陷周期：
[0097][0098]
图3展示了广义反射系数的幅度谱凹陷周期与时间厚度t的关系，其中图3中的(a)对应时间厚度为1.33ns，图3中的(b)对应时间厚度为4ns，利用时间厚度与幅度谱凹陷周期之间的反比关系，即可利用反射系数频谱的周期来估算出时间厚度。
[0099]
θ(f)分为两个部分，第一部分包含参数δt和反射系数奇偶分量，主要影响了相位谱的振幅大小。第二部分包含参数δt和t1，主要影响了相位谱的变化速度。反射系数奇偶分量包含介电常数和与电导率相关的衰减常数，则通过相位谱就可以估计出三层介质的介电常数和中间层的电导率。
[0100]
得到的反演初值包括上述估计的两层时间厚度、三层介质的介电常数和两层的电导率；将估计出的反演初值作为优化算法的初始输入值，实现快速准确的薄层参数反演。
[0101]
自此，就完成了一种基于反射信号重构的薄层介质参数反演方法。
[0102]
为了验证本发明提出的一种基于反射信号重构的薄层介质参数反演方法，设计仿真实验进行分析。如图4所示，仿真采用gprmax软件，构建三层介质模型，仿真参数如表格1所示。探地雷达被简化为彼此靠近的发射器和接收器，由金属板全反射波来获得雷达发射波，且需要确保天线到金属板的距离与到模型的距离相同。
[0103]
表格1仿真参数设置
[0104][0105]
反演结果如表2所示，可见本发明所提方法取得了较高的反演精度。经试验验证本
发明具有较快的收敛速度。
[0106]
表2仿真反演结果
[0107][0108]
同时，也分析了在不同信噪比条件下的反演结果，证明本发明方法的抗噪性能较强。
[0109]
本发明实施例中采用2ghz天线的ids c-thrue探地雷达，进行测试，场景底部铺有吸波材料，实验采用四层介质模型：空气-pvc材料-蒸压加气混凝土-蒸压灰砂砖。表4展示了实测模型的层厚反演结果，其误差均在10％以内，相对介电常数的结果也在合理的范围内。
[0110]
表3实测参数设置
[0111][0112]
表4实测反演结果
[0113][0114]
综上，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护
范围之内。