CN103439707A - 一种探地雷达加窗距离偏移成像方法 - Google Patents

Abstract

本发明公开了一种探地雷达加窗距离偏移成像方法，先对原始雷达记录剖面进行预处理，逐列进行加窗和一维傅里叶变换。窗函数的中心点根据各列绝对值最大值所在的时间采样点而定；再进行空间维的加窗和一维傅里叶变换。窗函数的中心点根据各列能量值最大值所在的空间采样点而定；最后进行谱域填充和二维逆傅里叶变换，得到二维成像结果。该探地雷达加窗距离偏移成像方法能提高成像质量。

Description

一种探地雷达加窗距离偏移成像方法

技术领域

本发明属于探地雷达成像技术领域，涉及一种探地雷达加窗距离偏移成像方法。

背景技术

探地雷达是一种有效的无损探测技术。它通过空域扫描向探测区域发射电磁波并接收散射回波，可实现对未知区域内部的成像处理，获得未知区域中的隐蔽目标参数，即目标分布信息和散射强度信息，有效应用于市政工程、考古、地雷探测、反恐等多种场合。探地雷达的空域扫描有沿线的一维扫描和在表面的二维扫描。一维扫描时，发射天线和接收天线分别以一定的间隔沿线移动。在每个位置处，发射天线向探测区域发射电磁波，接收天线接收探测区域的散射回波。然后移动发射天线和接收天线到下一个位置，又可以获得一道散射回波。通过在整个测线上移动发射天线和接收天线，便可以获得多道散射回波。发射天线和接收天线可以装配在一起同时移动，也可以分别移动。接收天线还可以选择为阵列天线的形式。这些配置方式分别对应于探地雷达应用中不同的扫描方式。本专利适用于各种探地雷达天线一维沿线等间距扫描方式下的成像方法。探地雷达成像的目的是从多道散射回波(即原始记录剖面)中恢复出地下区域的散射强度分布信息，即成像结果。一维扫描可以获得二维成像结果，其中一维为横向扫描维，另一维为纵向深度维。本专利针对一维扫描下的探地雷达数据进行加窗距离偏移成像。设一维测线沿地表布置，测线方向设为x方向，测线范围为[A，B]，该测线上共有L个测点，坐标分别为x_l=l·Δx，l=0，…，L-1，Δx表示空间采样间隔，即：空间维采样矢量为

发射天线和接收天线装配在一起同时移动。在测点x_l处，发射天线向探测区域发射电磁波，接收天线接收地下探测区域的散射回波，该点处的一维散射回波记为s_l(t)=[s_l(t₀)…s_l(t_k)…s_l(t_K-1)]^T，其中K表示时间维采样点数，上标T表示转置。时间维采样矢量为

其中t_k=k·Δt，k=0，…，K-1，Δt表示时间采样间隔，则探地雷达的采样时窗为W=t_K-1-t₀。则整个记录剖面数据E₀(x，t)可表示为E₀(x，t)=[s₁(t)…s_l(t)…s_L(t)]，即E₀(x，t)为一个二维矩阵，其尺寸为K×L，纵向维为时间维，横向维为空间维。成像的目的便是通过探地雷达记录剖面数据E₀(x，t)获得整个探测区域的散射强度值。

探地雷达的成像方法有多种，距离偏移成像算法适用于等间距采样下对均匀背景介质中点散射型目标的成像处理，广泛应用于探地雷达信号处理中【参考文献：雷文太.脉冲GPR高分辨成像算法研究.国防科学技术大学，博士学位论文，2006】。成像处理前需要对原始数据进行均衡、解振荡、去直达波、零点校正等预处理【参考文献：Jol H M主编.雷文太，童孝忠，周旸译.探地雷达理论与应用.北京：电子工业出版社，2011】，设预处理后的记录剖面数据为E₁(x，t)=[s′₁(t)…s′_l(t)…s′_L(t)]，仍为L列。此外，成像处理前，需设定成像区域并预知探测区域中电磁波传播的波速。对于下视探地雷达系统而言，成像区域的横向维矢量的取值区间一般取为原始扫描的测线范围。对于前视或斜前视的探地雷达系统而言，该取值区间需要根据具体的探测场景加以确定，此处统一记为[h_a，h_b]。纵向深度维矢量需根据探地雷达的探测深度进行选取，与时窗W有关，记为[g_a，g_b]。

传统的距离偏移成像方法中，首先对预处理后的记录剖面E₁(x，t)进行二维傅里叶变换，然后在空间谱域进行坐标变换和插值处理，获得目标成像区域的二维空间谱域值，再对该二维空间谱域值进行二维逆傅里叶变换，从而获得最终的二维成像结果。具体的成像步骤如下：

(1)对设定的成像区域进行网格划分

将该区域分别沿纵向深度维和横向维等间隔地划分为L_z和L_x个网格，则整个成像区域划分为L_z×L_x个网格。成像的目的便是获得该L_z×L_x个网格处的散射强度值，即O(z_m，x_n)m=1，…，L_z；n=1，…，L_x。在下视探地雷达成像处理中，一般设L_x=L，即：成像区域中的横向维网格划分数目与探地雷达记录剖面中的横向维采样点数相同，成像区域中的横向维矢量和探地雷达扫描时的横向维矢量相同。成像区域的纵向维矢量记为

其中M=L_z。

(2)对预处理后的记录剖面E₁(x，t)进行二维傅里叶变换

探地雷达的预处理剖面在空间维和时间维均为离散形式，其离散形式可表示为E₁(l·Δx，k·Δt)，1≤l≤L，0≤k≤K-1，其中l和k分别表示预处理剖面的横向(空间维)和纵向(时间维)的离散序号，Δx和Δt分别表示预处理剖面的横向(空间维)采样间隔和纵向(时间维)采样间隔。该记录剖面在空间——时间域(oxt平面)中为二维实矩阵。预处理剖面E₁(x，t)的二维离散傅里叶变换见下式：

$F_1(u,v)=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}\underset{k=0}{\overset{K-1}{\mathrm{\Σ}}}{E}_1(k,l)\·e^{-j(\frac{2\mathrm{\π}}{K}\mathrm{uk}+\frac{2\mathrm{\π}}{L}\mathrm{vl})}$

进行二维傅里叶变换后，空间——时间域(oxt平面)中的二维实矩阵转换为空间谱--频域(og_xf平面)中的二维复矩阵。记录剖面中的空间维矢量

经过二维傅里叶变换后对应于横向空间谱矢量

记录剖面中的时间维矢量

经过二维傅里叶变换后对应于频率矢量

矢量

和

的计算公式分别如下：

${\stackrel{\→}{g}}_x={\mathrm{\Δg}}_x\·[\mathrm{ceil}(-\frac{L}{2}):\mathrm{floor}\left(\frac{L-1}{2}\right)]$

$\stackrel{\→}{f}=\mathrm{\Δf}\·[\mathrm{ceil}(-\frac{K}{2}):\mathrm{floor}\left(\frac{K-1}{2}\right)]$

其中，

运算符ceil(α)表示大于等于α的最小整数，运算符floor(β)表示小于等于β的最大整数，矢量生成算符[λ：μ]表示以自然数1为间隔，产生从整数λ到整数μ的一维矢量。

(3)计算二维空间谱域(og_xg_z平面)中各网格点的复数值

该平面中成像区域的横向维(g_x维)和纵向维(g_z维)的矢量分别为

和 ${\stackrel{\→}{g}}_z,$ ${\stackrel{\→}{g}}_z=\mathrm{\Δ}{g}_z\·[0:M-1],$ 其中 $\mathrm{\Δ}{g}_z=\frac{1}{z_{M-1}-z_0}.$

对任一网格点(m·Δg_z，l·Δg_x)，根据下面的式子计算对应的频率值

$f_l\left(m\right)=\sqrt{{(m\·\mathrm{\Δ}{g}_z)}^2+{(l\·\mathrm{\Δ}{g}_x)}^2}$

然后，在空间谱——频域(og_xf平面)中，读取或计算第l列中频率值为f_l(m)所对应的复数值，作为该网格点(m·Δg_z，l·Δg_x)处的值，填充到该网格点。

遍历所有的网格点，从而获得二维空间谱域(og_xg_z平面)中各网格点的复数值，记为F₂(m，l)。

(4)计算二维空域(oxz平面)中的成像结果

对二维空间谱域(og_xg_z平面)中的填充结果F₂(m，l)进行二维逆傅里叶变换，并取其实部，公式如下：

$O(p,q)=\frac{1}{\mathrm{ML}}\·\mathrm{real}[\underset{l=0}{\overset{L-1}{\mathrm{\Σ}}}\underset{m=0}{\overset{M-1}{\mathrm{\Σ}}}{F}_2(m,l)\·e^{j(\frac{2\mathrm{\π}}{M}\mathrm{pm}+\frac{2\mathrm{\π}}{L}\mathrm{ql})}]$

其中，p和q分别表示二维空域深度维和横向维的采样点，运算符real[·]表示取实部操作。该二维矩阵O(p，q)即为距离偏移成像算法的成像结果。

传统的距离偏移成像算法中，对探地雷达记录剖面E₁(x，t)进行二维傅里叶变换时，并没有考虑到记录剖面在横向维和纵向维的局部特征，而是直接进行了二维傅里叶变换，这实际上相当于用与空间采样矢量和时间采样矢量等长的矩形窗分别进行横向维和纵向维的一维傅里叶变换。在探地雷达对地探测时，扫描区域中的目标分布通常是稀疏的。即在整个成像区域中，目标仅占据很少的空域部分。因此，目标的散射数据在探地雷达记录剖面中具有显著的局部特征。传统的距离偏移成像算法将雷达记录剖面中的所有数据视为同等重要的，实际上并不能获得高质量的成像结果。

因此，有必要设计一种新型的探地雷达距离偏移成像方法。

发明内容

本发明所要解决的技术问题是提供一种探地雷达加窗距离偏移成像方法，该探地雷达加窗距离偏移成像方法可提高成像质量。

发明的技术解决方案如下：

一种探地雷达加窗距离偏移成像方法，包括以下步骤：

步骤1：采用探地雷达进行一维沿线扫描，雷达记录的剖面数据为E₀(x，t)，对E₀(x，t)预处理后的探地雷达记录剖面数据为E₁(x，t)；其中x是检测点的横向位置，t表示时间；

步骤2：对E₁(x，t)逐列进行一维加窗傅里叶变换，得到二维复矩阵F₁ ^w(x_l，f)，0≤，≤L-1；

其中，L表示空间维的横向采样点数，l为空间维的横向采样点的序号；

步骤3：对F₁ ^w(x_l，f)，0≤，≤L-1进行空间维加窗傅里叶变换，得到H₁ ^w(g_x，f_k)；

步骤4：对H₁ ^w(g_x，f_k)进行谱域插值和二维逆傅里叶变换，从而获得最终的扫描区域的二维成像结果。【该步骤和传统距离偏移成像算法的(3)、(4)步相同，参见背景技术，为现有技术。】【原有技术是分别在时间维和空间维进行一维傅里叶变换，没有“加窗”二字。此处，本专利所提的方法，都是采用的“加窗”傅里叶变换，这也正是本专利的核心思想。】

所述的步骤2包括以下子步骤：

1)对任一列即第l列对应的一维矢量E₁(x_l，k)，先搜索其中绝对值最大值对应的时间采样点，记为k_s，即k_s满足：

E₁(x_l，k_s)=max[abs(E₁(x_l，k))]，k=0，…，K-1，K表示时间维的采样点数，k为时间维采样点的序号；

2)利用中心位置位于k_s采样点的窗函数对该一维矢量进行加窗处理，得到一维矢量E₁ ^w(x_l，k)=E₁(x_l，k)·w₁(k)，其中窗函数

$w_1\left(k\right)=0.5[1-\cos \left(2\mathrm{\π}\frac{k-k_s}{K}\right)],0\≤k\≤K-1;$

3)再对E₁ ^w(x_l，k)进行一维傅里叶变换，即

${F_1}^w(x_l,f)=\underset{k=0}{\overset{K-1}{\mathrm{\Σ}}}{E_1}^w(x_l,k)e^{-j\frac{2\mathrm{\π}}{K}\mathrm{kf}};$

遍历预处理后雷达记录剖面E₁(x，t)的所有列，逐列进行绝对值最大值位置的确定和加窗处理，再进行一维傅里叶变换，得到二维复矩阵F₁ ^w(x_l，f)，0≤l≤L-1。【此处的f无需定义，与权利要求书1中步骤2中的F₁ ^w(x_l，f)，0≤l≤L-1的f的含义相同。即：时间维采样序列t经过傅里叶变换后对应于频率维采样序列f】

所述的步骤3包括以下子步骤：

1)计算E₁(x，t)沿空间维的能量分布值

$A\left(l\right)=\underset{k=0}{\overset{K-1}{\mathrm{\Σ}}}{\left|E_1(l,t_k)\right|}^2;$

2)搜索该能量分布值中绝对值最大值对应的空间采样点，记为l_s，即l_s满足

A(l_s)=max[abs(A(l))]，l=0，…，L-1

3)对F₁ ^w(x_l，f)，0≤l≤L-1进行横向维的加窗处理，得到二维复矩阵G₁ ^w(l，f)=w₂(l)·F₁ ^w(l，f)，0≤l≤L-1，其中窗函数 $w_2\left(l\right)=0.5[1-cos\left(2\mathrm{\π}\frac{l-l_s}{L}\right)],0\≤l\≤L-1;$

4)对二维复矩阵G₁ ^w(l，f)的任一行G₁ ^w(l，f_k)，进行一维傅里叶变换，得到

g_x，是与横向维采样序列0≤l≤L-1相对应的空间谱序列。

此处的k表示是频率域的序号，参见以前的描述“步骤3：对F₁ ^w(x₁，f)，0≤l≤L-1进行空间维加窗傅里叶变换，得到H₁ ^w(g_x，f_k)；”，f_k实际上是频率变量f的离散化。

时间维采样序列t经过傅里叶变换后对应于频率维采样序列f，离散化后即为f_k。此处，横向维采样序列0≤l≤L-1，经过空间傅里叶变换后对应于空间谱序列g，离散化后即为g_x。

探地雷达进行一维沿线扫描，所获得的雷达记录剖面记为E₀(x，t)，E₀表示散射场的电场强度值，x是检测点的横向位置，t表示时间。该雷达记录剖面为一个二维矩阵，其尺寸为K×L，K表示时间维的采样点数，L表示空间维的横向采样点数。成像的目的便是通过探地雷达记录剖面数据E₀(x，t)获得整个探测区域的散射强度值；成像处理前需要对原始数据进行均衡、解振荡、去直达波、零点校正等预处理，设预处理后的记录剖面数据为E₁(x，t)=[s′₁(t)…s′_l(t)…s′_L(t)]，仍为L列；根据探地雷达空间扫描区域，确定二维成像区域，将该区域分别沿纵向深度维和横向维等间隔地划分为L_z和L_x个网格，则整个成像区域划分为L_z×L_x个网格。成像的目的便是获得该L_z×L_x个网格处的散射强度值，即O(z_m，x_n)，m=0，…，L_z-1；n=0，…，L_x-1。在下视探地雷达成像处理中，一般设L_x=L，即：成像区域中的横向维网格划分数目与探地雷达记录剖面中的横向维采样点数相同，成像区域中的横向维矢量和探地雷达扫描时的横向维矢量相同。成像区域的纵向维矢量记为

其中M＝L_z。

本发明的构思是：探地雷达对某个探测区域进行探测成像时，目标通常是稀疏分布的，目标的散射回波在整个雷达记录剖面中具有典型的局部特征。进行距离偏移成像时，分别对雷达记录剖面在时间维和空间维进行加窗的一维傅里叶变换，窗函数的中心位置根据目标散射回波所处的大致位置而定。从而较大程度上抑制了非目标散射回波的影响，提高了成像质量。

有益效果：

本发明的探地雷达加窗距离偏移成像方法，先对原始雷达记录剖面进行预处理，逐列进行加窗和一维傅里叶变换。窗函数的中心点根据各列绝对值最大值所在的时间采样点而定；再进行空间维的加窗和一维傅里叶变换。窗函数的中心点根据各列能量值最大值所在的空间采样点而定；最后进行谱域填充和二维逆傅里叶变换，得到二维成像结果。该探地雷达加窗距离偏移成像方法能提高成像质量。

本发明的独创性体现在：对某个探测区域进行距离偏移成像处理时，先进行时间维的一维加窗和傅里叶变换。窗函数的中心位置是根据雷达记录剖面中散射强度的绝对值最大值的位置而定的。然后再进行空间维的加窗处理和傅里叶变换。窗函数的中心位置是根据雷达记录剖面中各列散射回波能量值最大值的位置而定的。再进行空间谱域的填充和二维逆傅里叶变换。通过加窗处理，可以较大程度抑制非目标散射回波，提高成像质量。实验结果显示，本发明提出的方法与现有的距离偏移成像方法相比，成像质量有明显提高。

附图说明

图1示出了合成孔径天线对探测区域的扫描示意图；

图2示出了图1所示的探测区域中单个点目标的含噪空域散射回波；

图3示出了图2所示的原始散射回波经预处理后的雷达记录剖面；

图4示出了传统距离偏移成像算法的二维成像结果；

图5示出了加窗距离偏移成像算法的二维成像结果；

图6示出了图1所示的探测区域中多个点目标的分布图；

图7示出了图6所示的探测区域中多个点目标的含噪空域散射回波；

图8示出了图7所示的原始散射回波经预处理后的雷达记录剖面；

图9示出了传统距离偏移成像算法的二维成像结果；

图10示出了加窗距离偏移成像算法的二维成像结果；

标号说明：1-收发一体天线。

具体实施方式

以下将结合附图和具体实施例对本发明做进一步详细说明：

实施例1：

本实例是针对均匀背景媒质中目标的合成孔径探测成像，但本技术并不局限于合成孔径扫描，对实孔径扫描和多发多收扫描方式也适用。

首先运用探地雷达对探测区域进行一维扫描，如图1所示。收发一体天线沿x方向进行合成孔径扫描，扫描起止点分别为0.2m和1.8m，总的扫描长度为1.6m，扫描间隔为0.02m，一共在81个位置处进行探测记录。即：当收发一体天线移动到一个孔径点处后保持静止，发射天线向下辐射电磁波，接收天线接收地下探测区域的散射回波。完成一次扫描后，收发一体天线移动到下一个。探测区域如图1所示，在x=1m、z=0.4m处放置一个半径为0.02m的细长型铁棒，长度远大于其半径，导体棒的取向垂直于纸面。发射天线为线源形式，向地下发射电磁波。散射回波被阵列天线接收，共81道数据，如图2所示。图2中的原始散射数据经过去直达波和时窗截断预处理后的数据如图3所示。下面分别运用传统距离偏移成像算法和加窗距离偏移成像算法对预处理后数据进行成像处理。横向成像区域为[0.2，1.8]m，纵向成像区域为[0，1.3]m。传统距离偏移成像结果O₁(x，z)如图4所示。采用hann窗分别对预处理雷达记录剖面的纵向时间维和横向空间维进行一维加窗和傅里叶变换，然后再进行谱域填充和二维逆傅里叶变换，得到加窗距离偏移成像算法的处理结果O₂(x，z)，如图5所示。两个成像结果中，都给出了点目标的成像结果，但图4中的成像结果的旁瓣要高于图5中的旁瓣，即：传统距离偏移成像算法的聚焦效果要差于加窗距离偏移成像。对二维成像结果O(x_q，z_p)，0≤p≤M-1，0≤q≤L-1，定义聚焦度因子

$\mathrm{\Ω}=\frac{{\left[\underset{p=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{q=0}{\overset{L-1}{\mathrm{\Σ}}}{\left[O(x_q,z_p)\right]}^2\right]}^2}{\underset{p=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{q=0}{\overset{L-1}{\mathrm{\Σ}}}{\left[O(x_q,z_p)\right]}^4}$

该聚焦因子Ω可用来表示二维图像的聚焦程度。理想情况下，一个点目标的聚焦度因子为1。用上式分别对传统距离偏移成像结果和加窗距离偏移成像结果进行计算，得到二者的聚焦度因子分别为Ω₁=17.31，Ω₂=7.56。与传统距离偏移成像结果相比，加窗距离偏移成像结果的聚焦度改善了56.3％。

实施例2：

探地雷达的天线配置同实施例1。探测区域如图6所示，在三个位置处放置细长型铁棒，长度远大于其半径，导体棒的取向垂直于纸面。发射天线为线源形式，向地下发射电磁波。散射回波被阵列天线接收，共81道数据，原始记录剖面和预处理后的记录剖面分别如图7和图8所示。二维成像区域的设置同实施例1。传统距离偏移成像算法和加窗距离偏移成像算法的成像结果分别如图9和图10所示，同样采用hann窗函数。从成像结果同样可以看出，加窗距离偏移成像算法有更好的聚焦效果。计算二者的聚焦度因子，分别为38.59和13.28。加窗后，成像结果的聚焦度改善了65.6％。

Claims (3)

1.一种探地雷达加窗距离偏移成像方法，其特征在于，包括以下步骤：

步骤1：采用探地雷达进行一维沿线扫描，雷达记录的剖面数据为E₀(x，t)，对E₀(x，t)预处理后的探地雷达记录剖面数据为E₁(x，t)；其中x是检测点的横向位置，t表示时间；

步骤2：对E₁(x，t)逐列进行一维加窗傅里叶变换，得到二维复矩阵F₁ ^w(x_l，f)，0≤l≤L-1；

其中，L表示空间维的横向采样点数，l为空间维的横向采样点的序号；

步骤3：对F₁ ^w(x_l，f)，0≤l≤L-1进行空间维加窗傅里叶变换，得到H₁ ^w(g_x，f_k)；

步骤4：对H₁ ^w(g_x，f_k)进行谱域插值和二维逆傅里叶变换，从而获得最终的扫描区域的二维成像结果。

2.根据权利要求1所述的探地雷达加窗距离偏移成像方法，其特征在于，所述的步骤2包括以下子步骤：

1)对任一列即第l列对应的一维矢量E₁(x_l，k)，先搜索其中绝对值最大值对应的时间采样点，记为k_s，即k_s满足：

E₁(x_l，k_s)=max[abs(E₁(x_l，k))]，k=0，…，K-1，K表示时间维的采样点数，k为时间维采样点的序号；

2)利用中心位置位于k_s采样点的窗函数对该一维矢量进行加窗处理，得到一维矢量E₁ ^w(x_l，k)=E₁(x_l，k)·w₁(k)，其中窗函数

0≤k≤K-1：

3)再对E₁ ^w(x_l，k)进行一维傅里叶变换，即

${F_1}^w(x_l,f)=\underset{k=0}{\overset{K-1}{\mathrm{\Σ}}}{E_1}^w(x_l,k)e^{-j\frac{2\mathrm{\π}}{K}\mathrm{kf}};$

遍历预处理后雷达记录剖面E₁(x，t)的所有列，逐列进行绝对值最大值位置的确定和加窗处理，再进行一维傅里叶变换，得到二维复矩阵F₁ ^w(x_l，f)，0≤l≤L-1。

3.根据权利要求2所述的探地雷达加窗距离偏移成像方法，其特征在于，所述的步骤3包括以下子步骤：

1)计算E₁(x，t)沿空间维的能量分布值

$A\left(l\right)=\underset{k=0}{\overset{K-1}{\mathrm{\Σ}}}{\left|E_1(l,t_k)\right|}^2;$

2)搜索该能量分布值中绝对值最大值对应的空间采样点，记为l_s，即l_s满足

A(l_s)=max[abs(A(l))]，l=0，…，L-1

3)对F₁ ^w(x_l，f)，0≤l≤L-1进行横向维的加窗处理，得到二维复矩阵G₁ ^w(l，f)=w₂(l)·F₁ ^w(l，f)，0≤l≤L-1，其中窗函数 $w_2\left(l\right)=0.5[1-cos\left(2\mathrm{\π}\frac{l-l_s}{L}\right)],0\≤l\≤L-1;$

4)对二维复矩阵G₁ ^w(l，f)的任一行G₁ ^w(l，f_k)，进行一维傅里叶变换，得到g_x，是与横向维采样序列0≤l≤L-1相对应的空间谱序列。

Images