一种地质雷达B-scan数据处理方法与流程

本发明属于地质雷达信号的信号处理领域，特别涉及一种地质雷达b-scan数据处理方法。

背景技术：

隧道施工是公路和铁路建设中的一个重要问题，其中隧道衬砌施工是目前隧道施工中的一个重要环节。后期的隧道衬砌检测是对施工质量的一个重要评判。另一方面，铁路路基病害诱发路基失稳，是影响铁路安全运营的重要因素，开展有效的铁路路基检测方法具有重要意义。布置于检测目标表面的地质雷达(ground-penetratingradar，gpr)通过发射天线发射宽频带短脉冲的高频电磁波，电磁波在衬砌、路基及更深层介质的传播过程中经过存在电性差异的介质体或者界面时产生的反射信号被布置于表层的接收天线接收。地质雷达具有超浅层勘探的独特优势，是一种快速的高分辨率无损检测方法，在隧道以及路基检测中具有广泛应用。

目前基于地质雷达b-scan数据的信号处理方法主要有

现有技术1：平均值减去法

(1)在b-scan剖面上选取连续的若干道，将这些道对应时间点的值累加并除以道数，从而得到平均道；

(2)对b-scan上的每一道减去(1)中的平均道，得到最终的数据处理结果；

现有技术1的缺点：

由于实际地质雷达数据较为复杂，用于计算平均道的道数较难选择，选择不当时易对有用信号造成损伤，从而影响后续的衬砌异常检测；该方法假设杂波及界面反射信号在b-scan剖面上呈现为严格的水平结构，而实际的接收数据通常不满足该假设，因此该方法的处理效果会受到一定的限制；

现有技术2：滤波方法

在b-scan剖面上，杂波的主要成分及界面反射信号通常表现为水平的带状结构，且杂波的主要成分通常在剖面的最上方，而钢筋反射信号为抛物线形态，利用它们在时间-空间域的特征区别或者在频率-波数域的不同分布区域，可以设计时间-空间域滤波器或者频率-波数域滤波器，从而将两种成分进行分离；

现有技术2的缺点：

这些滤波器通常是基于数值模拟设定滤波参数，而实际的地质雷达b-scan剖面和模拟剖面会有一定的偏差，因为将设定的滤波器参数用于实际的b-scan剖面处理时，会导致滤波参数难以设置合适的问题，从而影响后续的分离效果以及检测效果。

技术实现要素：

基于此，本公开揭示了一种地质雷达b-scan数据处理方法，其特征在于，该方法包括以下步骤：

s100、对b-scan剖面上的每一道数据(a-scan)，采用谱白化拓频方法进行处理，得到高时间分辨率的b-scan剖面；

s200、基于相位移偏移成像方法，对高时间分辨率b-scan数据进行二维偏移成像，得到b-scan成像剖面；

s300、基于形态成分分析方法，对步骤s200得到的b-scan成像剖面进行迭代稀疏分离；

s400、对经过迭代稀疏分离出的结果，分析其点状特性和层状结构。

本发明具有如下有益效果：

1)本公开的谱白化频谱拓展方法可以有效的提高b-scan剖面的时间分辨率以及增强深层弱反射的强度；

2)本公开采用相位移偏移可以对b-scan剖面高效地进行成像，一方面使得钢筋产生的双曲线型的散射信号聚焦回钢筋的地下位置，使得钢筋成像为点状特征，而混凝土层和黏土层或围岩的界面产生的反射波成像以后仍然为层状结构，从而实现了对两者特征的有效区分，这表明相位移偏移可以提高b-scan剖面的横向分辨率，有益于目标的横向精细检测；

3)本公开的分离方法利用了不同结构稀疏表示特征的不同，对不同形态成分进行稀疏分离，具体采用无下采样的二维离散小波变换对点状结构进行稀疏表示，二维曲波变换对层状结构进行稀疏表示。这实现了变换基与所分析特征的一致性，因此可以有效的分离点状结构和层状结构。

附图说明

图1为本公开一个实施例中的地质雷达b-scan数据处理流程；

图2为本公开一个实施例中的基于谱白化方法的傅里叶域b-scan数据频谱拓展流程；

图3(a)为本公开一个实施例中的形态成分分析的稀疏分离流程；

图3(b)为本公开一个实施例中的无下采样小波原子示意；

图3(c)为本公开一个实施例中的曲波原子示意；

图4(a)为本公开一个实施例几何模型示意；

图4(b)为本公开一个实施例中利用有限差分方法生成的b-scan剖面；

图4(c)为本公开一个实施例偏移后的b-scan剖面；

图4(d)为本公开一个实施例中分离的点状结构分量；

图4(e)为本公开一个实施例中分离的层状结构分量；

图5(a)为本公开一个实施例中实测的b-scan剖面；

图5(b)为本公开一个实施例中谱白化拓频处理后的高分辨率b-scan剖面；

图5(c)为本公开一个实施例偏移后的b-scan剖面；

图5(d)为本公开一个实施例中分离的点状结构分量；

图5(e)为本公开一个实施例中分离的层状结构分量；

图6(a)为本公开一个实施例中实测的b-scan剖面；

图6(b)为本公开一个实施例中分离的点状结构分量；

图6(c)为本公开一个实施例中分离的层状结构分量。

具体实施方式

下面结合附图1至附图6(c)和具体实施方式对本发明做进一步详细的说明。

在一个实施例中，本公开揭示了一种地质雷达b-scan数据处理方法，该方法包括以下步骤：

s100、对b-scan剖面上的每一道数据(a-scan)，采用谱白化拓频方法进行处理，得到高时间分辨率的b-scan剖面；

s200、基于相位移偏移成像方法，对高时间分辨率b-scan数据进行二维偏移成像，得到b-scan成像剖面；

s300、基于形态成分分析方法，对步骤s200得到的b-scan成像剖面进行迭代稀疏分离；

s400、对经过迭代稀疏分离出的两个成分，分析其点状特性和层状结构。

进一步的，如果b-scan剖面的某一道数据的时间分辨率较高则步骤s100省略。

在本实施例中，为了解决上述现有技术所存在的问题，根据需要检测的目标信号特征，本实施例提供了一种地质雷达b-scan剖面信号处理方法，以服务于后期的衬砌或路基检测。首先采用基于谱白化的信号频谱拓展(拓频)技术提高二维地质雷达b-scan数据信号时间方向的分辨率，若接收到的二维地质雷达数据时间分辨率较高，即所述某一道数据的带宽大于5个倍频程，该步骤可以省略。然后将拓频结果进行相位移偏移成像，对混凝土层内的钢筋进行聚焦成像，最后基于钢筋等点目标与反射界面(混凝土与黏土层或围岩界面)等层特征的区别，采用基于形态成分分析的稀疏分离方法分别得到钢筋网以及反射界面等目标。

在一个实施例中，所述步骤s100中的拓频方法具体包括以下步骤：

s101、计算b-scan剖面上的每一道数据的一维傅里叶变换，得到所述每一道数据的振幅谱和相位谱；

s102、在所述每一道数据的有效频宽范围内，设定若干个带通滤波器，将振幅谱分为若干个子带；

s103、计算每个子带内的数据振幅谱加权系数；

s104、将每个子带内的振幅谱乘以对应的加权系数，相位谱不变，得到拓频后的频谱；

s105、基于拓频后的频谱，计算b-scan剖面上的每一道数据的一维傅里叶反变换，得到拓频后的b-scan剖面上的每一道数据。

在本实施例中，所述步骤s102中带通滤波器的数量为5-8个。

在一个实施例中，所述步骤s200中的相位移偏移成像方法包括以下步骤：

s201、对高时间分辨率处理后的b-scan数据进行二维傅里叶变换，对应的关系式为：

式中，f(x，z，t)|z＝0表示在检测目标表层接收到的b-scan数据，x表示空间坐标，z表示深度，t表示时间，为虚数单位；f(kx，z，ω)|z＝0为b-scan信号的二维傅里叶变换，kx表示空间频率(空间波数)，ω为角频率；

s202、将相移算子作用于深度为z的波场f(kx，z，ω)，得到深度为z+δz处的波场；

s203、基于二维傅里叶反变换，并令时间t＝0，得到b-scan的相位移偏移成像结果。

在一个实施例中，所述步骤s300中的迭代稀疏分离由以下关系式决定：

式中，f为偏移成像后的b-scan剖面，f1表示其中的钢筋及其他点结构分量，f2表示其中的水平或近似水平层状分量，为二维无下采样离散小波变换基，为f1小波变换系数的l1范数，为二维曲波变换基，为f2曲波变换系数的l1范数，为对应误差的误差项，λ为正则化参数，满足λ＞0。

在一个实施例中，所述步骤s400中的稀疏分离出的结果包括两个成分：一个是包括钢筋或其他点状特征的剖面；另一个为平层结构，其中最浅的层为发射天线到接收天线的直达波以及发射天线在检测目标表层反射波成像得到的层，深部的层为衬砌与围岩或黏土及其他介质层的分界面。

在一个实施例中，所述步骤s103具体为：对每个子带，计算其振幅均方根值，用一固定常数除以该均方根值，即得到对应子带内振幅谱的加权系数。

在本实施例中，所述固定常数为整体比例因子，选择2000，3000均可，没有固定要求。

在一个实施例中，步骤s202中的计算公式为：

式中，为相移算子，δz为延拓步长，υ为电磁波在隧道衬砌内的传播速度的一半。

在一个实施例中，步骤s202中b-scan的成像结果计算公式如下

然后根据深度z和υ的换算关系，将成像剖面换算为时间-空间剖面

在一个实施例中，本发明提供一种用于隧道衬砌异常以及路基病害检测的地质雷达数据信号处理方法，包括对地质雷达数据的高分辨拓频处理、偏移成像以及基于形态成分分析的稀疏分离等步骤，以得到用于异常检测的点状结构分量和层状结构分量。

如图1所示，本发明所提供的地质雷达数据处理方法，具体通过如下步骤实施：

准备地质雷达b-scan剖面，进行如下处理：

步骤1：对b-scan剖面上的每一道，即每一个a-scan，f0(t)进行傅里叶域频谱拓频处理，以提高其分辨率；

参照图2，傅里叶域频谱拓展处理包括如下步骤：

(1)对a-scanf0(t)进行一维傅里叶变换，公式如下：

式中，t表示时间，i为虚数单位，f0(ω)为a-scan信号的傅里叶变换，ω为角频率。由于采集数据为离散数据，因此一维傅里叶变换的实施可采用其离散傅里叶变换的快速算法实现，以提高计算效率。

(2)在信号的有效频带范围内，设置k个窄带带通滤波器hk，每个带通滤波器的频率范围为(ωk-1，ωk)，k＝1，2，l，k，对信号进行频率滤波，得到k个子带谱为hk{f0(ω)}，k＝1，2，l，k；

(3)设定常数因子s，在每个带通滤波器hk(k＝0，1，l，k)的频率范围(ωk-1，ωk)内，首先计算信号频谱的平均能量密度pk为

其中，|hk{f0(ω)}|表示hk{f0(ω)}的模值；

然后，计算该频率范围内频谱的加权系数wk为

(4)对每个带通滤波器hk(k＝0，1，l，k)频率范围内的信号频谱乘以相应的加权系数wk，得到拓展后的频谱，公式如下：

(5)将拓频后的频谱合并为做一维傅里叶反变换，得到拓频后的信号，公式如下：

步骤2：基于相位移偏移成像方法，对高分辨拓频处理后的b-scan数据进行二维偏移成像；

首先将b-scan信号f(x，z，t)|z＝0作为输入信号，这里z表示深度，z＝0表示在表层接收到的b-scan数据，x表示空间坐标，相位移偏移处理主要包括如下步骤：

(1)对f(x，z，t)|z＝0做二维傅里叶变换，对应的关系式为：

式中，f(kx，z，ω)|z＝0为b-scan信号的二维傅里叶变换，kx表示空间频率(空间波数)。由于采集数据为离散数据，因此二维傅里叶变换的实施采用其离散傅里叶变换的快速算法来实现，通常情况下需要将数据在空间和时间方向上补零，使得空间采样点数nx和时间上的采样点数nt均为2的整次幂，从而以提高计算效率。dt和dx为时间和空间方向上的采样间隔，相应的kx和ω均为离散的序列值，其采样间隔dkx＝1/(dxgnx)，dω＝2π/(dtgnt)，kx的取值为0，dkx，2dkx，l，(nx-1)dkx，ω的取值为0，dω，2dω，l，(nt-1)dω。

(2)将相移算子作用于深度z处的波场f(kx，kδz，ω)，可以延拓得到深度为(k+1)δz处的波场为：

式中，δz为延拓步长，通常取值为0.1厘米或者0.2厘米等，深度离散化为z＝0，δz，2δz，l，kδz，l，(k-1)δz，延拓步数为k，其取值取决于检测深度。υ为电磁波在隧道衬砌内传播速度的一半，跟钢筋混凝土的相对介电常数有关，设置为其中c0＝299792458米/秒为光在真空中的传播速度，εr为介质的相对介电常数。εr通常取为已有的经验参数，但是受到具体施工条件的约束以及混凝土介质的变化，在具体的计算中，通常需要微调εr，标准为使b-scan中的钢筋反射双曲线信号聚焦能量最强。为相移算子，具体计算中，当时，令

(c)利用二维傅里叶反变换，并令时间t＝0，我们可以重构b-scan信号为：

式中空间频率kx上的傅里叶反变换采用快速傅里叶变换算法实现，在频率ω方向上则直接进行叠加运算。对每个深度处的波场均进行上述操作。

步骤3：基于形态成分分析，对b-scan进行迭代稀疏分离，基于形态成分分析的迭代稀疏分离由以下关系式决定：

式中，f以及已知，f1，f2为待求的分量。f为偏移后的b-scan剖面，f1为待分离的分量1，表示其中的钢筋等点结构分量，f2为分量2，表示其中的水平层状分量，这些均为向量化表示。为二维的无下采样离散小波变换基，小波函数选取为6阶的symmlet小波，也可以选取为其他阶次的symmlet小波或者daubechies小波等，为二维曲波变换基，尺度数选取为8，图3(b)和图3(c)分别给出了二维无下采样离散小波变换的基原子示意和二维曲波原子示意。为f1小波变换系数向量的l1范数，为f2曲波变换系数向量的l1范数，定义为系数向量所有元素的绝对值之和，为对应的误差项，定义为f-f1-f2所有元素的平方和。λ为正则化参数，满足λ＞0。如图3(a)所示，f1的初始值定义为f，f2为0初始化。每次迭代中，执行以下步骤：

(1)固定f2，采用迭代软阈值(iterativesoft-thresholding)的方法求解如下关系式定义的反问题：

得到f1；

(2)固定f1，采用迭代软阈值(iterativesoft-thresholding)的方法求解如下关系式定义的反问题：

得到f2。

如图3(a)所示，迭代终止条件用最大迭代次数限定，选取为300次，即可以满足计算精度的要求，又能显著提高计算效率。

步骤4：对分离出的两个成分，进行分析；

步骤3分离出的两个分量f1和f2分别表示偏移后b-scan剖面中的点状结构分量和层状结构分量，根据隧道衬砌的结构设置，点状分量剖面主要包括钢筋等点状特征，层状结构分量主要有浅层的层状杂波，主要由发射天线到接收天线的直达波层以及衬砌表面的反射波构成，深部的层状结构为衬砌与回填层或者黏土层、围岩的分界面。

下面利用本公开提供的数据处理方法对(1)模拟的隧道衬砌地质雷达b-scan数据和(2)实测的b-scan数据进行了处理；对于有限差分合成的模型数据，由于介质较为均匀，反射波衰减少，时间分辨率较高，因此省略了步骤1的高分辨拓频处理，仅进行了步骤2-步骤4的处理。对于实测的b-scan数据进行了步骤1-步骤4的处理分析。

图4(a)-图4(e)为模型数据及其处理结果：

如图4(a)所示，其为模型的几何结构示意，混凝土层厚度为0.45米，混凝土下方为黏土层。钢筋距离表层0.05米，每根钢筋的半径为0.01米，相邻钢筋的横向中心距离为0.1米。发射源信号为脉冲源，发射天线的中心频率为1ghz，发射天线和接收天线的横向距离为0.01米，天线的横向移动间隔为0.01米。图4(b)为生成的b-scan剖面，横坐标为水平方向，纵坐标表示接收时间，可以看到早期接收的水平强杂波以及钢筋反射的双曲线波形。钢筋反射的双曲线波形相互干扰，降低了剖面的横向分辨率，不利于钢筋的检测。在7ns附近的混凝土和黏土层界面的水平反射层能量较弱；并且有部分反射信号被钢筋反射的双曲线波形所覆盖，如图中黑框所示，因此无法进行有效识别。图4(c)为相位移偏移成像后的b-scan剖面，钢筋反射的双曲线波形被聚焦为增强的点状结构，不同于杂波和深层反射波的水平结构。基于这两种不同的结构成分，利用步骤5的形态成分分析稀疏分离得到图4(d)和图4(e)的两个分量。图4(d)为点状结构分量，图4(e)为水平层状结构分量。可以看出，钢筋等点状结构主要分布在图4(d)上，而水平杂波以及深层反射界面主要分布在图4(e)上，而且在图4(b)上被覆盖的部分水平反射层在图4(e)上也得到了清晰显示。

图5(a)-图5(e)为实测隧道衬砌地质雷达数据处理结果：

图5(a)为实测的用于衬砌异常检测的b-scan数据，发射天线的主频为900mhz，隧道管片厚度为35厘米，其中分布两层钢筋，主筋半径为0.01米，配筋半径为0.005米。从图5(a)上可以较为明显地看到浅层的强能量层状杂波以及管片中钢筋反射的双曲线波形。而由于实际介质的频散作用，管片与后方围岩的界面反射在图5(a)上较弱，很难用于对隧道衬砌异常进行有效检测。图5(b)为高分辨拓频处理结果，可以看出剖面的时间方向分辨率得到了有效提高，而且深层的反射波能量得到了有效增强。图5(c)为偏移成像结果，钢筋反射的双曲线波形得到了有效聚焦，深层反射界面表现更为平滑，这使得两者呈现出不同的形态结构。图5(d)和图5(e)为利用步骤3的分离结果。图5(d)为点状结构分量，可以清晰的看到上下两组钢筋的分布。图5(e)为水平层状结构分量，深层的管片与围岩的反射界面得到了清晰刻画。

图6(a)-图6(c)为实测铁路路基地质雷达数据处理结果：

图6(a)为实测的用于铁路路基检测的b-scan数据，其中发射天线的主频为900mhz，地下分布有网状的钢筋结构。图中可以看到浅层的强能量杂波和钢筋网反射的双曲线结构。而且，地下界面反射的层状结构由于和钢筋反射的双曲线结构相互干扰，而无法清晰识别。利用本发明专利的方法，首先对图6(a)的b-scan数据进行基于谱白化的高分辨处理，然后利用步骤3进行稀疏分离。图6(b)和图6(c)为对应的分离结果。图6(b)为分离的点状结构，可以清晰分辨出钢筋的分布。图6(c)为分离的层状结构，除浅层的强杂波外，地下反射界面也得到了显著增强。这为后面的路基病害检测提供了有效的处理结果。

本公开在实际应用中，随着深度的增加，电磁波的传播能量快速衰减，且频率变低，因此从衬砌界面反射回的信号能量较弱，主频较低，直接在b-scan剖面上利用现有技术进行识别难度大，本发明基于谱白化的频率拓展方法可以有效的提高成像剖面的时间分辨率以及弱反射的幅度，有利于后期的界面检测。本发明对拓频后的b-scan进行高效的相位移偏移成像，使得钢筋反射的抛物线型的散射信号聚焦回钢筋的地下位置，且表现为点状特征，空间分辨率得到了有效提高，衬砌界面产生的反射波成像以后仍然为线状结构。这说明相比较于现有技术，相位移偏移算法可以使得它们表现为区别更大的结构特征，为后续的分离提供良好前提。本发明利用形态成分分析稀疏分离法将钢筋等点状结构和界面的层状结构进行分离，其中采用无下采样的二维离散小波变换对点状结构进行稀疏表示，二维曲波变换对层状结构进行稀疏表示。本发明的分离方法利用了不同结构的稀疏表示特点，因此可以有效的分离点状结构和层状结构。相比较于现有技术1，本发明的数据处理方法能否适应中等非均匀介质的情况；且相比较于现有技术2的固定参数滤波法，本发明的分离方法能够自适应于分离数据，自动调整参数，不需人为设置参数。

最后需要说明的是，以上模型和实际资料算例对本发明的目的，技术方案以及有益效果提供了进一步的验证，这仅属于本发明的具体实施算例，并不用于限定本发明的保护范围，在本发明的精神和原则之内，所做的任何修改，改进或等同替换等，均应在本发明的保护范围内。