一种利用小波变换进行探地雷达属性特征融合方法与流程

本发明属于近地表地球物理领域，涉及一种基于小波变换的探地雷达属性特征融合方法。

背景技术：

探地雷达属性分析通过对目标的几何和物理性质进行定量描述，有助于提取雷达记录中隐藏的信息，因此通过提取不同的探地雷达属性，可以更好地了解地下结构和性质。然而，单一雷达属性只能部分显示地下目标的细节，多个属性可以获取更多关于地下性质的信息量。不同数据成分的综合解释是一项艰巨的任务，需要采用数据融合的方法来优化有用信息的提取。

技术实现要素：

本发明的目的是针对现有技术的不足，提供一种基于小波变换的探地雷达属性特征融合方法。该方法是利用小波变换将来自同一探地雷达数据的不同属性剖面融合为一个复合的输出显示。

本发明具体采用的技术方案如下：

一种利用小波变换进行探地雷达属性特征融合方法，该方法包括如下步骤：

s1：对探地雷达的采集数据进行预处理后，得到探地雷达的目标剖面图像；

s2：从探地雷达的目标剖面图像中提取出待融合的属性特征图像；

s3：基于小波变换对两个属性特征图像a和b进行融合，融合方法如下：

s31：分别对每个属性特征的图像进行小波变换，小波变换的公式如下：

其中：下标j表示分解尺度，cj+1为图像信号cj的低频分量，分别为图像信号cj沿水平、垂直和对角线方向上的高频分量；当j＝0时，c0是未经过分解的原始图像信号；h和g分别表示与尺度函数和小波函数相对应的滤波系数矩阵；h’和g’分别为h和g的共轭转置矩阵；

s32：针对两个属性特征图像经过小波变换后得到的高频分量，基于极大值融合规则进行特征融合，公式如下：

式中，dm(i，j)为融合图像在(i，j)位置处的细节系数值；da(i，j)、db(i，j)分别为属性特征图像a和b在(i，j)位置处的细节系数值，即小波变换得到的图像高频分量。

s33：针对两个属性特征图像经过小波变换后得到的低频分量，基于边缘提取融合规则进行特征融合，公式如下：

式中：cm(i，j)为融合图像在(i，j)位置处的近似系数值，ca(i，j)和cb(i，j)分别为属性特征a和属性特征b的图像在(i，j)位置处的近似系数值，即小波变换得到的图像低频分量；ga和gb分别为属性特征a和属性特征b的近似系数的梯度矩阵，ga(i，j)和gb(i，j)表示(i，j)位置处的梯度矩阵ga和gb；

属性特征a和属性特征b的梯度矩阵gz计算通式为：

其中，z∈[a，b]，“*”’表示卷积，和分别是sobel边缘检测算子在水平和垂直方向上的模板；

s4：对融合图像的近似系数和细节系数进行小波逆变换，得到融合后的图像。

作为优选，所述的s1中，对探地雷达采集数据的预处理包括数据编辑、静校正、直流成分去除、增益、带通滤波和一致性噪声去除。

作为优选，所述的s2中，从探地雷达的目标剖面图像中提取出的属性特征为振幅一阶导数和图像纹理同质性特征。

进一步的，所述振幅一阶导数的提取方法如下：

s211：首先给出复道数信号：

c(t)＝f(t)+ig(t)

及其希尔伯特变换为：

其中，f(t)为输入信号，虚部部分g(t)为f(t)经线性卷积变换而来，g(t)和f(t)的相位差为90°；t为信号采样时间；

s212：计算信号的瞬时振幅：

a(t)＝[f(t)²+g(t)²]^1/2

s213：对瞬时振幅按时间进行求导，得到振幅一阶导数，即包络线的时间变化率。

进一步的，所述图像纹理同质性特征通过灰度共生矩阵(glmc)方法提取，其计算公式为：

式中，pi，j表示所述目标剖面图像信号在(i，j)位置处的灰度共生矩阵值，n为图像的灰度级。

本发明可对地下潜在目标体的物理和几何性质进行多种属性特征的提取和优化选择，并通过小波变换将这些不同的属性剖面综合成一个复合的输出显示，从而同时利用振幅、相干度、反射特征、相位和光谱内容等信息来改善地下成像，提高数据解释的精度和效率。

因此，相对于现有技术而言，本发明通过小波变换将来自同一探地雷达数据的不同属性进行融合，优化有效信息的提取，帮助我们更好地描述特征和定位潜在目标，进一步打开自动或计算机辅助解释的道路。

附图说明

图1为利用小波变换进行探地雷达属性特征融合方法流程图；

图2为某处地下冰川层位探测的500mhz振幅图。

图3为该处500mhz振幅图中提取的振幅一阶导数属性。

图4为该处500mhz振幅图中提取的同质性属性。

图5为振幅一阶导数和同质性的属性融合图像。

具体实施方式

下面结合附图和具体实施例本发明做进一步阐述和说明。

如图1所示，为本实施例中的利用小波变换进行探地雷达属性特征融合方法流程，该方法包括如下步骤：

第一步、对采集的数据进行常规处理：

探地雷达的采集数据需要预先进行数据处理，得到探地雷达的目标剖面图像。具体的数据处理过程可根据即使的测量对象进行调整，以能够得到准确、完整的信号为准。在本发明中，数据处理包括数据编辑、静校正、直流成分去除、增益、带通滤波和一致性噪声去除等；此外如果地下条件需要，还包括偏移。

第二步、属性特征的提取：

当获得探地雷达的目标剖面图像后，需要根据实际的分析和应用，对需提取的属性特征进行优化选择。根据选择的属性特征，从目标剖面图像中提取出相应的属性特征图像，记为属性a、属性b……属性n，用于进行后续融合。已知的探地雷达属性多达上百种，具体提取的属性特征需要根据实际需要进行确定。

在本实施例中，从探地雷达的目标剖面图像中提取出的属性特征有两个，分别为振幅一阶导数和图像纹理同质性特征，记为属性特征a、属性特征b。下面分别对两个属性特征的具体提取方法进行叙述。

1)振幅一阶导数的提取方法

本发明通过计算复道数信号，提取出探地雷达信号的振幅一阶导数属性。其具体过程如下：

首先给出复道数信号：

c(t)＝f(t)+ig(t)

及其希尔伯特变换为：

其中，f(t)为输入信号，虚部部分g(t)为f(t)经线性卷积变换而来，g(t)和f(t)的相位差为90°；t为信号采样时间；

则信号的瞬时振幅的计算式如下：

a(t)＝[f(t)²+g(t)²]^1/2

而振幅一阶导数为瞬时振幅的时间导数，即包络线的时间变化率，它显示了反射事件能量的变化。因此，对瞬时振幅按时间进行求导，得到振幅一阶导数。

2)图像纹理特征的提取方法

图像纹理的纹理特征可以通过灰度共生矩阵(glmc)方法来提取。而同质性(homogeneity)属于纹理属性的一种，它是图像整体平滑度的度量，其值越大表示图像局部纹理越均匀，其计算如下：

式中，pi，j表示所述目标剖面图像信号在(i，j)位置处的灰度共生矩阵值，n为图像的灰度级。

第三步、基于小波变换对两个属性特征图像a和b进行融合，融合步骤如下：

1)分别对每个属性特征图像(a、b)进行小波变换：

假设h和g分别表示与尺度函数和小波函数相对应的滤波系数矩阵，则二维的小波变换可以用以下公式表示：

其中：下标j表示分解尺度，cj+1为图像信号cj的低频分量，分别为图像信号cj沿水平、垂直和对角线方向上的高频分量；当j＝0时，c0是未经过分解的原始图像信号；h′和g′分别为h和g的共轭转置矩阵。

信号经过小波变换后，其高频分量实际上代表了细节系数，而低频分量则代表了近似系数。后续可分别在高频域和低频域上采用不同的规则进行融合。

2)高频域的细节系数融合规则：

细节系数反映了图像中较大对比度变化的特征，如图像边缘。因此，针对两个属性特征图像信号经过小波变换后得到的高频分量，基于极大值融合规则进行特征融合，公式如下：

式中，dm(i，j)为融合图像信号在(i，j)位置处的细节系数值；da(i，j)、db(i，j)分别为属性特征a和属性特征b的图像信号在(i，j)位置处的细节系数值，即小波变换得到的高频分量。

3)低频域的近似系数融合规则：

使用sobel边缘检测算子，可以按下式计算属性特征图像a的近似系数的梯度矩阵ga：

其中，ga(i，j)表示(i，j)位置处的梯度矩阵ga，ca(i，j)为属性特征图像a在(i，j)位置处的近似系数值，即小波变换得到的低频分量；“*”’表示卷积，和分别是sobel边缘检测算子在水平和垂直方向上的模板。

同样，可以得到属性特征图像b的近似系数的梯度矩阵gb：

其中，gb(i，j)表示(i，j)位置处的梯度矩阵gb，cb(i，j)为属性特征图像b在(i，j)位置处的近似系数值，即小波变换得到的低频分量。

由于梯度矩阵g能够突出探地雷达图像中的边缘信息，所以我们选择在每个图像中具有较大g值的近似系数作为在同一空间位置上的待融合近似系数。具体而言，针对两个属性特征图像经过小波变换后得到的低频分量，基于边缘提取融合规则进行特征融合，公式如下：

式中：cm(i，j)为融合图像信号在(i，j)位置处的近似系数值，ca(i，j)和cb(i，j)分别为属性特征图像a和b在(i，j)位置处的近似系数值，即小波变换得到的低频分量；ga和gb分别为属性特征a和属性特征b的梯度矩阵；

第四步、最后对融合图像的近似系数和细节系数进行小波逆变换，得到融合后的图像。小波逆变换可以根据前述的小波变换公式逆向进行，不再赘述。

下面将上述方法应用于某一具体的实例中，以展示其技术效果。

图2是该处地下冰川层位探测的500mhz振幅图。白色箭头指示层位及其内部的反射。

图3是按照前述方法从该500mhz图中提取的振幅一阶导数属性，它可以用于检测尖锐的界面和不连续性。图中白色曲线为对地下分界面的追踪，以确定地下的层位分布；而白色虚线指示了潜在的层内连续/不连续的反射。

图4是按照前述方法从该500mhz图中提取的同质性属性，它在局部可以突出显示探地雷达图像的均匀性。从图中可以看到白色箭头指示的位置呈现出图像中强振幅反射和低同质性的一致趋势的特征；同时，它可以帮助我们进行层位的识别，对于连续/不连续的层位反射，根据低同质性的聚类分布可以进行层位追踪。

图5是按照前述方法将振幅一阶导数和同质性这两个属性特征进行融合后的融合图像。通过一个多属性的复合显示图，可以快速获得地下综合的视图。白色曲线指示了层位的分布情况，而白色箭头指示了潜在的层理结构。此外，强振幅反射和低同质性表现出高度的一致性，这强化了对地下目标的形态和展布的刻画，增强了我们的视觉效果，有利于数据的解释。

以上所述的实施例只是本发明的一种较佳的方案，然其并非用以限制本发明。有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型。因此凡采取等同替换或等效变换的方式所获得的技术方案，均落在本发明的保护范围内。