一种基于衍射和成像与最小熵技术的探地雷达波速估计方法.pdf

本发明公开了一种基于衍射和成像与最小熵技术的探地雷达波速估计方法，首先利用均值法去除图像直达波分量；接下来计算图像二维方向的归一化能量，选择合适的阈值确定目标感兴趣区域；然后根据介质的介电常数确定波速范围，对于每个波速值，在目标感兴趣区域内根据目标双曲线特征利用衍射和方法对图像进行合成孔径成像处理，计算处理后图像的熵值；最后选取最小熵值对应的图像为最佳成像，其对应的波速为最佳波速。本发明通过选择目标感兴趣区域和设置累加点水平位置范围减小衍射和处理的复杂度，利用图像熵值自动确定最佳成像和最佳速度，在较完整保留目标信息的情况下有效降低计算量，计算简单、鲁棒性强、估计精度高，适合于工程应用。

1.一种基于衍射和成像与最小熵技术的探地雷达波速估计方法，其特征在于，包括以
下步骤：
(1)输入探地雷达二维B扫描图像e(xi,tj)，1≤i≤M，1≤j≤N，其中横坐标xi＝i·Δx
为水平位置，纵坐标tj＝j·Δt为信号往返时间，M为道数，N为每道的数据样点数，Δx为采
样的道间距，Δt为采样时间间隔；
(2)利用均值法对B扫描图像进行处理，去除地表直达波，如下所示：
$\begin{array}{cc}{e}_1(x_i,t_j)=e(x_i,t_j)-\frac{1}{M}\underset{i=1}{\overset{M}{\Σ}}e(x_i,t_j)& 1\≤i\≤M,1\≤j\≤N\end{array}---\left(1\right)$
其中e1(xi,tj)为去除地表直达波后的图像；
(3)对去除直达波后的图像，计算水平位置和时间两个方向的归一化能量，如下所示：
$\begin{array}{cc}{P}_x\left(i\right)=\underset{j=1}{\overset{N}{\Σ}}{e}_1(x_i,t_j)\·e_1(x_i,t_j)& 1\≤i\≤M,1\≤j\≤N\end{array}---\left(2\right)$
$\begin{array}{cc}\stackrel{\‾}{P_x\left(i\right)}=\frac{P_x\left(i\right)}{max\left(P_x\right(i))}& 1\≤i\≤M\end{array}---\left(3\right)$
$\begin{array}{cc}{P}_t\left(j\right)=\underset{i=1}{\overset{M}{\Σ}}{e}_1(x_i,t_j)\·e_1(x_i,t_j)& 1\≤i\≤M,1\≤j\≤N\end{array}---\left(4\right)$
$\begin{array}{cc}\stackrel{\‾}{P_t\left(j\right)}=\frac{P_t\left(j\right)}{max\left(P_t\right(j))}& 1\≤j\≤N\end{array}---\left(5\right)$
其中max(·)表示取最大值，为水平位置方向的归一化能量，时间方向的归一
化能量，选择阈值对和进行处理，得到矩形目标感兴趣区域R为[xi min≤xi≤
xi max，tj min≤tj≤tj max]，其中，imin为水平位置的起始点坐标序号，imax为水平位置的终
点坐标序号，jmin为时间的起始点坐标序号，jmax为时间的终点坐标序号；
(4)根据介质的特性，确定介电常数er取值范围er min≤εr≤εr max，由波速值v与介电常
数εr的关系确定波速范围为：
$\frac{c}{\sqrt{{\mathrm{\ϵ}}_{rmax}}}\≤v\≤\frac{c}{\sqrt{{\mathrm{\ϵ}}_{rmin}}}---\left(6\right)$
(5)选择波速步进为Δv，对于波速范围内每一个波速值v，在目标感兴趣区域R内利用
衍射和方法对图像进行处理，如下所示：
$\begin{array}{cc}{e}_2(x_i,t_j,v)=\underset{k}{\Σ}{e}_1(x_k,t_{i,j,k})& imin\≤i\≤imax,jmin\≤j\≤jmax\end{array}---\left(7\right)$
其中e2(xi,tj,v)为衍射和方法处理后图像，k为测量点(xk,0)的水平位置坐标序号，
ti,j,k是电磁波从测量点(xk,0)到目标点(xi,tj)的来回时间，定义如下：
$t_{i,j,k}=\frac{2R_{i,j,k}}{v}---\left(8\right)$
其中Ri,j,k为测量点(xk,0)到(xi,tj)的距离；
(6)对于波速范围内每一个波速值v，计算衍射和成像处理后的图像熵值Q(v)，如下所
示：
$Q\left(v\right)=\frac{{\[\underset{j=jmin}{\overset{jmax}{\Σ}}\underset{i=imin}{\overset{imax}{\Σ}}{e_2}^2(x_i,t_j,v)\]}^2}{\underset{j=jmin}{\overset{j\max }{\Σ}}\underset{i=imin}{\overset{imax}{\Σ}}{e_2}^4(x_i,t_j,v)}---\left(9\right)$
(7)选取最小熵值对应的衍射和成像为最佳成像e2(xi,tj,v)opt：
$e_2{(x_i,t_j,v)}_{opt}=\arg \underset{e_2}{min}\[Q\left(v\right)\]---\left(10\right)$
其中表示当图像熵Q(v)取得最小值时，对应的衍射和成像结果e2(xi,tj,
v)。
最佳成像e2(xi,tj,v)opt对应的波速为最佳波速：
$v_{opt}=\arg \underset{v}{min}\[Q\left(v\right)\]---\left(11\right)$
其中表示当图像熵Q(v)取得最小值时，对应的波速v。
2.根据权利要求1所述的基于衍射和成像与最小熵技术的探地雷达波速估计方法，其
特征在于，所述步骤(3)中，目标感兴趣区域R的确定方法如下：
i)分别确定水平位置方向归一化能量的阈值T1和时间方向归一化能量的阈值T2，如下
所示：
$T1=p1\·max\left(\stackrel{\‾}{P_x\left(i\right)}\right)=p1---\left(12\right)$
$T2=p2\·max\left(\stackrel{\‾}{P_t\left(i\right)}\right)=p2---\left(13\right)$
其中max(·)表示取最大值，p1和p2分别为水平位置和时间方向归一化能量的阈值系
数，一般取0<p1<0.1，0<p2<0.1；
ii)对于水平位置方向的归一化能量从位置零点开始，向后递推，将各位置点能
量数据与阈值T1比较，取第一个大于阈值T1的数据坐标序号为起始点坐标序号imin，其对
应的坐标为：
$\begin{array}{cc}{x}_{imin}=x_i|\stackrel{\&OverBar;}{P_x\left(i\right)}\&gt;T1& i=1,2,...,M\end{array}---\left(14\right)$
iii)对于水平位置方向归一化能量从位置终点开始，向前递推，将各位置点能量
数据与阈值T1比较，取第一个大于阈值T1的数据坐标序号为终点坐标序号i max，其对应的
坐标为：
$\begin{array}{cc}{x}_{imax}=x_i|\stackrel{\&OverBar;}{P_x\left(i\right)}\&gt;T1& i=M,M-1,...,1\end{array}---\left(15\right)$
iv)对于时间方向归一化能量从时间零点开始，向后递推，将各时间点能量数据
与阈值T2比较，取第一个大于阈值T2的数据坐标序号为起始点坐标序号j min，其对应的坐
标为：
$\begin{array}{cc}{t}_{jmin}=t_j|\stackrel{\&OverBar;}{P_t\left(j\right)}\&gt;T2& j=1,2,...,N\end{array}---\left(16\right)$
v)对于时间方向归一化能量从时间终点开始，向前递推，将各时间点能量数据与
阈值T2比较，取第一个大于阈值T2的数据坐标序号为终点坐标序号jmax，其对应的坐标为：
$\begin{array}{cc}{t}_{jmax}=t_j|\stackrel{\&OverBar;}{P_t\left(j\right)}\&gt;T2& j=N,N-1,...,1\end{array}---\left(17\right)$
vi)根据水平位置和时间两个方向的起始点和终点坐标，得到目标感兴趣区域R为矩形
区域：
[xi min≤xi≤xi max，tj min≤tj≤tj max] (18)。
3.根据权利要求1所述的基于衍射和成像与最小熵技术的探地雷达波速估计方法，其
特征在于，所述步骤(5)中测量点(xk,0)的水平位置坐标序号k的取值范围为：
$k=\left\{\begin{array}{cc}\&lsqb;i-id,i+id\&rsqb;,& id+1\&lt;i\&lt;M-id\\ \&lsqb;1,i+id\&rsqb;,& 1\&le;i\&le;id+1\\ \&lsqb;i-id,M\&rsqb;,& M-id\&le;i\&le;M\end{array}\right.---\left(19\right)$
其中i为目标点(xi,tj)的水平位置坐标序号，取值范围为[i min,i max]，id为衍射和
算法水平位置范围坐标序号的阈值。
4.根据权利要求1所述的基于衍射和成像与最小熵技术的探地雷达波速估计方法，其
特征在于，所述步骤(5)中Ri,j,k计算如下：
$R_{i,j,k}=\sqrt{{(x_k-x_i)}^2+{\left(\frac{{\mathrm{vt}}_j}{2}\right)}^2}---\left(20\right)$
由此得到ti,j,k的计算式为：
$t_{i,j,k}=\frac{2R_{i,j,k}}{v}=\sqrt{\frac{4{(x_k-x_i)}^2}{v^2}+{t_j}^2}---\left(21\right)$
将Δx和Δt代入式(21)，得到目标点(xi,tj)分布的双曲线中，水平位置坐标序号k对应
的时间坐标序号ji,j,k为：
$j_{i,j,k}=\sqrt{\frac{4{\mathrm{\&Delta;x}}^2{(k-i)}^2}{{\mathrm{\&Delta;t}}^2{v}^2}+j^2}---\left(22\right).$

一种基于衍射和成像与最小熵技术的探地雷达波速估计方法

技术领域

本发明属于数字信号处理领域，特别涉及到探地雷达的B扫描图像处理，具体涉及
一种基于衍射和成像与最小熵技术的探地雷达波速估计方法。

背景技术

探地雷达基于电磁波传播与散射原理，通过向地下发射电磁波信号并接收地下介
质不连续处散射的回波实现对地下目标的探测。与电阻率法、低频电磁感应法和地震法等
地下探测方法相比，探地雷达具有探测速度快、探测过程连续、分辨率高、操作方便灵活、探
测费用低、探测范围广(能探测金属和非金属)等优点，在地质、资源、环境、工程和军事等领
域得到广泛的应用。

探地雷达对地下目标进行探测时，电磁波在不同介质中的传播速度不同，波速是
合成孔径成像和目标定位的关键参数之一，是解释探测数据的基础。波速估计不准确将导
致合成孔径成像分辨率下降，影响对目标的探测，因此波速估计是探地雷达的一个重要研
究领域。

目前探地雷达的波速估计方法主要有三种：介电常数法、反射系数法和目标信息
法：

1)介电常数法是直接利用波速和介质介电常数的关系进行计算，其核心是测量介
质的介电常数，由于相对介电常数的测量复杂，而且随着环境水分含量的不同变化较大。该
方法存在估计误差大、操作复杂的缺点；

2)反射系数法是通过在地面放置金属板，首先通过测量无金属板时地面回波强度
和放置金属板后地面回波强度得到反射系数，然后利用反射系数和相对介电常数之间的关
系计算相对介电常数，最后通过波速和相对介电常数的关系计算波速值，该方法只能测量
地下浅层介质的波速，适用场合较窄；

3)目标信息法主要分两类：延时估计法和双曲线估计法。延时估计法是已知目标
深度的情况下，通过测量目标回波延时估计波速。这种方法虽然简单，但是在浅地层应用中
目标回波延时估计误差较大，影响波速估计的精度。双曲线估计法是利用B扫描图像中目标
分布的双曲线特征，在一定波速范围内对利用双曲线约束关系，对图像进行参数空间变换
或合成孔径成像处理，根据处理效果确定波速估计值；常见的参数空间变换方法有Hough变
换和Radon变换等；常见的合成孔径成像方法有频率－波数域偏移成像法、衍射和成像法和
微波全息成像法等。

在上述几种方法中，介电常数法、反射系数法和延时估计法估计误差较大，可操作
性不强；参数空间变换方法和合成孔径成像方法利用目标的双曲线特征进行波速估计，具
有实现方便的优点，是目前研究的热点。现有的合成孔径成像方法一般先利用合成孔径处
理实现对目标的聚焦，然后通过人工比较不同的波速下的成像效果得到最佳波速，存在计
算量大和误差大的缺陷。因此，如何以较小的计算量实现对波速的高精度自动估计，对于改
善合成孔径成像效果，提高探地雷达的探测性能具有重要意义。

发明内容

本发明要解决的技术问题是，为了克服现有的合成孔径成像波速估计方法的缺
陷，提供一种基于衍射和成像与最小熵技术的探地雷达波速估计方法，该方法采用计算复
杂度低的衍射和合成孔径成像方法，通过选择目标感兴趣区域和设置水平位置范围减少衍
射和处理的计算量，利用图像最小熵技术自动选择最佳成像和最佳波速，解决现有方法计
算量大、精度低和鲁棒性差的问题。

为了解决上述技术问题，本发明采用的技术方案如下：

一种基于衍射和成像与最小熵技术的探地雷达波速估计方法，包括以下步骤：

(1)输入探地雷达二维B扫描图像e(xi,tj)，1≤i≤M，1≤j≤N，其中横坐标xi＝i·
Δx为水平位置，纵坐标tj＝j·Δt为信号往返时间，M为道数(总列数)，N为每道的数据样
点数(总行数)，Δx为采样的道间距(水平间隔距离)，Δt为采样时间间隔；

(2)利用均值法对B扫描图像进行处理，去除地表直达波，如下所示：

其中e1(xi,tj)为去除地表直达波后的图像；

(3)对去除直达波后的图像，计算水平位置和时间两个方向的归一化能量，如下所
示：

其中max(·)表示取最大值，为水平位置方向的归一化能量，时间方向的
归一化能量，选择阈值对和进行处理，得到矩形目标感兴趣区域R为[ximin≤xi≤
ximax，tjmin≤tj≤tjmax]，其中，imin为水平位置的起始点坐标序号，imax为水平位置的终点
坐标序号，jmin为时间的起始点坐标序号，jmax为时间的终点坐标序号；

(4)根据介质的特性，确定介电常数εr取值范围εrmin≤εr≤εrmax，由波速值v与介电
常数εr的关系确定波速范围为：

(5)选择波速步进为Δv，对于波速范围内每一个波速值v，在目标感兴趣区域R内
利用衍射和方法对图像进行处理，如下所示：

其中e2(xi,tj,v)为衍射和方法处理后图像，k为测量点(xk,0)的水平位置坐标序
号，ti,j,k是电磁波从测量点(xk,0)到目标点(xi,tj)的来回时间，定义如下：

其中Ri,j,k为测量点(xk,0)到(xi,tj)的距离；

(6)对于波速范围内每一个波速值v，计算衍射和成像处理后的图像熵值Q(v)，如
下所示：

(7)选取最小熵值对应的衍射和成像为最佳成像e2(xi,tj,v)opt：

其中表示当图像熵Q(v)取得最小值时，对应的衍射和成像结果e2(xi,
tj,v)。

最佳成像e2(xi,tj,v)opt对应的波速为最佳波速：

其中表示当图像熵Q(v)取得最小值时，对应的波速v。

进一步，所述步骤(3)中，目标感兴趣区域R的确定方法如下：

i)分别确定水平位置方向归一化能量的阈值T1和时间方向归一化能量的阈值T2，
如下所示：

其中max(·)表示取最大值，p1和p2分别为水平位置和时间方向归一化能量的阈
值系数，一般取0<p1<0.1，0<p2<0.1；

ii)对于水平位置方向的归一化能量从位置零点开始，向后递推，将各位置
点能量数据与阈值T1比较，取第一个大于阈值T1的数据坐标序号为起始点坐标序号imin，
其对应的坐标为：

iii)对于水平位置方向归一化能量从位置终点开始，向前递推，将各位置点
能量数据与阈值T1比较，取第一个大于阈值T1的数据坐标序号为起始点坐标序号imax，其
对应的坐标为：

iv)对于时间方向归一化能量从时间零点开始，向后递推，将各时间点能量
数据与阈值T2比较，取第一个大于阈值T2的数据坐标序号为起始点坐标序号jmin，其对应
的坐标为：

v)对于时间方向归一化能量从时间终点开始，向前递推，将各时间点能量数
据与阈值T2比较，取第一个大于阈值T2的数据坐标序号为起始点坐标序号jmax，其对应的
坐标为：

vi)根据水平位置和时间两个方向的起始点和终点坐标，得到目标感兴趣区域R为
矩形区域：

[ximin≤xi≤ximax，tjmin≤tj≤tjmax] (18)。

进一步，所述步骤(5)中测量点(xk,0)的水平位置坐标序号k的取值范围为：

其中i为目标点(xi,tj)的水平位置坐标序号，取值范围为[imin,imax]，id为衍射
和算法水平位置范围坐标序号的阈值。

进一步，所述步骤(5)中Ri,j,k计算如下：

由此得到ti,j,k的计算式为：

将Δx和Δt代入式(21)，得到目标点(xi,tj)分布的双曲线中，水平位置坐标序号k
对应的时间坐标序号ji,j,k为：

本发明的有益效果：

1、本发明中目标感兴趣区域选择为矩形区域，在确定区域范围时，分别从二维坐
标(水平位置和时间)的零点和终点数据开始与阈值比较，向中间部分数据递推，可降低目
标感兴趣区域的搜索计算量；

2、本发明在衍射和合成孔径成像计算中，通过设置水平位置范围选择用于累加的
双曲线目标像素点，可进一步降低算法的累加计算量；

3、本发明通过计算不同波速下衍射和成像的熵值，利用最小熵技术自动选择最佳
成像和最佳波速，可在较完整保留目标信息的情况下有效降低计算量，具有鲁棒性强、估计
精度高的优点，适合于工程应用。

附图说明

图1为本发明的实现流程图；

图2为探地雷达原始B扫描图像；

图3为去除直达波后的图像；

图4为水平位置方向的归一化能量；

图5为时间方向的归一化能量；

图6为目标感兴趣区域内图像熵值随波速变化曲线；

图7最小熵值对应的衍射和成像结果。

具体实施方式

下面结合附图和实施例，对本发明进行进一步详细说明。如图1所示，本发明基于
衍射和成像与最小熵技术的探地雷达波速估计方法，具体步骤如下：

1)采用FDTD方法来仿真生成探地雷达的B扫描图像，如图2所示。仿真模型参数如
下：

(a)地下介质为干砂，其相对介电常数为εr＝4，电磁波中心频率为900MHz，真实波
速为1.50×108m/s；

(b)仿真区域宽度为3m,深度为2m，两个目标均为理想圆柱形导体，半径为0.1m，水
平位置分别为1.0m和1.8m，深度分别为0.2m和0.3m；

(c)道间距Δx为0.01m，采样时间间隔Δt为0.01179ns，每道数据有1018个采样
点，总采样时间为12ns；

2)利用均值法对原始图像进行处理，得到去除地表直达波后的图像，如图3所示；

3)对去除直达波后的图像，搜索目标感兴趣区域，步骤如下：

a)计算水平位置和时间两个方向的归一化能量；

b)设置水平位置方向归一化能量阈值p1＝0.05，时间方向归一化能量阈值p2＝
0.03；

c)对于两个方向的归一化能量，分别从坐标零点和终点的数据开始与阈值比较，
向中间部分数据递推，得到目标感兴趣区域R为[0.6m≤xi≤2.2m，2.0756ns≤tj≤
6.1679ns]，结果分别如图4和图5所示；

4)根据干砂的相对介电常数范围3≤εr≤5，得到波速v的范围为：

1.34×108m/s≤v≤1.73×108m/s；

5)选择波速步进为Δv＝0.01×108m/s，对于波速范围内每一个波速值v，设置衍
射和算法水平位置范围坐标序号阈值id＝25，在目标感兴趣区域R内利用衍射和方法对图
像进行合成孔径成像处理；

6)对于波速范围内每一个波速值v，计算衍射和成像处理后的图像熵值，得到目标
感兴趣区域内不同波速下图像熵值变化的曲线，如图6所示，最小熵值为396.3，此时对应的
波速为1.60×108m/s，估计误差为6.67％；

7)选取最小熵值对应的衍射和成像为最佳成像，如图7所示。

显然，上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的
实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其
它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。而这些属于本发
明的精神所引伸出的显而易见的变化或变动仍处于本发明的保护范围之中。