CN102830401A - 一种探地雷达加窗加权后向投影成像方法 - Google Patents

Abstract

本发明公开了一种探地雷达加窗加权后向投影成像方法，对成像区域中的每个成像单元，根据扫描场景计算该成像单元时延曲线，然后提取原始记录剖面经预处理后的数据剖面在该时延曲线上的散射数据，生成一维散射数据矢量。根据探地雷达天线的波瓣宽度和成像单元的深度，计算截断长度W_a，然后结合成像单元的横向位置截取该一维散射数据矢量中长度为W_a的一段数据。根据所截取的数据，设计加权因子进行成像运算。遍历所有的成像单元，从而完成成像运算。该探地雷达加窗加权后向投影成像方法能提高成像质量。

Description

一种探地雷达加窗加权后向投影成像方法

技术领域

本发明属于探地雷达成像技术领域，涉及一种探地雷达加窗加权后向投影成像方法。

背景技术

探地雷达是一种有效的无损探测技术。它通过空域扫描向探测区域发射电磁波并接收散射回波，可实现对未知区域内部的成像处理，获得未知区域中的隐蔽目标参数，即目标分布信息和散射强度信息，有效应用于市政工程、考古、地雷探测、反恐等多种场合。探地雷达的空域扫描有沿线的一维扫描和在表面的二维扫描。一维扫描时，发射天线和接收天线分别以一定的间隔沿线移动。在每个位置处，发射天线向探测区域发射电磁波，接收天线接收探测区域的散射回波。然后移动发射天线和接收天线到下一个位置，又可以获得一道散射回波。通过在整个测线上移动发射天线和接收天线，便可以获得多道散射回波。发射天线和接收天线可以装配在一起同时移动，也可以分别移动。接收天线还可以选择为阵列天线的形式。这些配置方式分别对应于探地雷达应用中不同的扫描方式。本专利适用于各种探地雷达天线一维沿线扫描方式下的成像方法。探地雷达成像的目的是从多道散射回波（即原始记录剖面）中恢复出地下区域的散射强度分布信息，即成像结果。一维扫描可以获得二维成像结果，其中一维为横向扫描维，另一维为纵向深度维。

本专利针对一维扫描下的探地雷达数据进行多尺度后向投影成像。设一维测线沿地表布置，测线方向设为x方向，测线范围为[A,B]，该测线上共有L个测点，坐标分别为x_i，i＝1,…,L。发射天线和接收天线装配在一起同时移动。在测点x_i处，发射天线向探测区域发射电磁波，接收天线接收地下探测区域的散射回波，该点处的一维散射回波记为s_i(t)=[s_i(t₀)…s_i(t_k)…s_i(t_K-1)]^T，其中K表示时间维采样点数，上标T表示转置。其采样时窗为W=t_K-1-t₀。则整个记录剖面数据E₀(x,t)可表示为E₀(x,t)=[s₁(t)…s_i(t)…s_L(t)]，即E₀(x,t)为一个二维矩阵，其尺寸为K×L。成像的目的便是通过探地雷达记录剖面数据E₀(x,t)获得整个探测区域的散射强度值。

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

传统的加权后向投影成像方法中，是将成像区域分别沿纵向深度维和横向维等间隔地划分为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。具体的成像步骤如下：

（1）对成像区域中的每一网格位置(z_m,x_n)，根据探测扫描场景计算各测点x_i处对应的电磁波双程传播时延τ_m，n，i，提取该测点处的一维回波信号s_i′(t)在时刻τ_m，n,i处的值q_m，n，i，生成一维信号w_m,n=[q_m，n，1，…,q_m，n，L]。

（2）将该一维信号的加权求和作为网格(z_m,x_n)处的成像结果，即(z_m,x_n)处的散射强度为

其中加权因子α(z_m,x_n)定义为 $\mathrm{\α}(z_m,x_n)=\left\{\begin{array}{cc}1,& s=0\\ \frac{u}{s},& s\≠0\end{array}\right..$ 其中s和u分别为该成像点对应的时延曲线处的散射回波w_m,n的标准差和均值。

（3）遍历成像区域中所有的网格，分别计算各点的散射强度O(z_m,x_n)即可获得整个成像区域的成像结果O(z_m,x_n)m=1,…,L_z;n=1,…,L_x。

传统的加权后向投影成像方法进行成像处理时，对成像区域中的每个单元，提取该成像单元对应双曲线上的所有散射数据，并进行加权累加成像处理。而在GPR（地质雷达）应用中，扫描区域中通常存在多个目标，各目标的绕射曲线通常存在部分重叠的现象。进行成像处理时，如果对每个成像单元采用整个时延曲线上散射回波数据进行加权成像，将会使得相邻目标的散射数据交叉重叠的现象，造成多目标成像质量的下降。

因此，有必要设计一种探地雷达加窗加权后向投影成像方法。

发明内容

本发明所要解决的技术问题是提供一种探地雷达加窗加权后向投影成像方法，该探地雷达加窗加权后向投影成像方法能显著提高多目标区域的成像质量。

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

一种探地雷达加窗加权后向投影成像方法，将整个成像区域划分为L_z×L_x个网格，一维测线沿地表布置，一维测线方向设为x方向，该一维测线上共有L个测点，坐标分别为x_i，i＝1,…,L；通过合成孔径天线对探测区域的扫描；加窗加权后向投影成像的方法包括以下步骤：

步骤1：提取各成像单元对应的时延曲线上散射回波数据；

对成像区域中的每一网格位置(z_m,x_n)，m=1,…,L_z;n=1,…,L_x，根据探测扫描场景计算第i个孔径测点x_i处对应的电磁波双程传播时延τ_m，n，i，提取该测点x_i处的一维回波信号s_i′(t)在时刻τ_m，n,i【解释：下标i表示第i个孔径的测点编号，该测点的横向坐标为x_i。对每个孔径测点i，遍历成像区域中的每个网格位置(z_m,x_n)，都可以计算获得传播时延值，因此，此处的传播时延值用τ_m，n,i表示。】处的幅度值q_m，n，i【解释：q_m，n，i即为一维信号s_i′(t)在t=τ_m，n，i时的幅度值，即：q_m，n，i=s_i′(t=τ_m，n，i)，此处用q进行了替代】，生成一维信号w_m,n=[q_m，n，1,…,q_m，n，L]；【传播时延τ_m，n，i是通过深度、空间位置以及介电常数计算得来的，为现有成熟计算方法】

步骤2：计算有效扫描宽度W_a并截取数据：

W_a的计算公式为

其中，z_m为成像单元的深度，θ₁表示探地雷达天线的零功率波瓣宽度；【零功率波瓣宽度可根据实际的天线形式加以计算获得【参考文献：Jol H M主编.雷文太,童孝忠,周旸译.探地雷达理论与应用.北京:电子工业出版社,2011】，round(·)表示四舍五入取整；

截取数据是指截取一维信号w_m,n中长度为W_a的一段数据

作为成像单元(z_m,x_n)对应的时延曲线加窗后的散射回波数据；

所截取的这一段散射数据

为：

式中，

二者满足W_l+W_r+1=W_a，其中fix(·)表示截尾取整；

步骤3：对成像单元(z_m,x_n)对应的时延曲线加窗范围中的散射数据进行加权求和作为该成像单元的散射强度值，即

$O(z_m,x_n)=\mathrm{\α}(z_m,x_n)\mathrm{\Σ}{w}_{m,n}^c=\mathrm{\α}(z_m,x_n)\underset{l=l_1}{\overset{l_2}{\mathrm{\Σ}}}{q}_{m,n,l};$

式中l₁和l₂分别表示对一维信号w_m,n=[q_m，n，1，…,q_m，n，L]进行加窗截取时，窗口的左端和右端的序号值；加权因子α(z_m,x_n)的计算方法如下：

其中，e_lf、e_mid和e_rg分别表示将加窗后的散射回波矢量平均分为三段后各段散射数据的能量值，计算方法为 $e_{\mathrm{lf}}=\underset{l=1}{\overset{a}{\mathrm{\Σ}}}{\left[w_{m,n}^c\left(l\right)\right]}^2,e_{\mathrm{mid}}=\underset{l=a+1}{\overset{b}{\mathrm{\Σ}}}{\left[w_{m,n}^c\left(l\right)\right]}^2,$ $e_{\mathrm{rg}}=\underset{l=b+1}{\overset{W_a}{\mathrm{\Σ}}}{\left[w_{m,n}^c\left(l\right)\right]}^2,$ 其中 $a=\mathrm{fix}\left(\frac{W_a}{3}\right),b=\mathrm{fix}\left(\frac{{2W}_a}{3}\right)$ 分别表示将长度为W_a的加窗后散射回波矢量

平均分为三段时各段的分界点；

${\mathrm{\α}}^{\′}=\left\{\begin{array}{cc}1,& s=0\\ \frac{u}{s},& s\≠0\end{array}\right.,$ 其中s和u分别为该成像点对应的时延曲线加窗后的散射回波矢量

的标准差和均值；

按照前述的步骤1-步骤3遍历成像区域中所有的网格，分别计算各点的散射强度O(z_m,x_n)即可获得整个成像区域的成像结果O(z_m,x_n)m=1,…,L_z;n=1,…,L_x。

L取值范围为[3,10000]。

本发明的构思是：探地雷达对某个探测区域进行探测成像时，探测区域中目标的绕射曲线通常是有限宽度的，而不是占据整个扫描长度。当探测区域中存在多个目标时，各个目标的绕射曲线存在部分重叠。因此，传统的基于整个扫描长度的时延曲线上散射回波的加权成像存在误差。为提高成像质量，本发明提出了加窗加权的后向投影成像方法。通过计算获得不同深度处成像单元的时延曲线截断长度。基于该长度，对整个扫描长度范围内的时延曲线进行截断处理。然后对该截断区间内的散射回波数据进行加权因子的计算和累加成像运算。

有益效果：

本发明的探地雷达加窗加权后向投影成像方法，对每个成像单元进行成像处理时，首先提取该成像单元对应的时延曲线上的散射回波数据；然后根据探地雷达天线的波瓣宽度和该成像单元所处的空间位置仅截取一维散射回波数据中的一段数据；再基于该段数据的统计特征设计加权因子并进行加权处理。遍历所有的成像单元从而完成成像运算。该方法对某个成像单元进行成像处理时，仅对该成像单元对应的有效范围内的散射数据进行加权相加处理，避免了传统方法对整个扫描范围内散射数据的加权处理，从而抑制了该成像单元有效范围之外的其他目标的散射数据的干扰，提高了成像质量。实验结果显示，本发明提出的方法与现有的成像方法相比，成像质量显著提高。

附图说明

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

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

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

图4示出了传统的加权后向投影成像算法中某些成像单元对应的时延曲线；

图5示出了传统的加权后向投影成像算法的成像结果；

图6示出了加窗加权后向投影成像算法中某些成像单元对应的时延曲线；

图7示出了加窗加权后向投影成像算法的成像结果；

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

具体实施方式

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

实施例1：

实施例1：

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

首先运用探地雷达对探测区域进行一维扫描，如图1所示。区域z>0为空气，区域z<0为介质层，z=0为空气—介质的交界面。在介质层中的7个位置上，分别放置了半径为0.02m的细长型铁棒，长度远大于其半径，导体棒的取向垂直于纸面。收发一体天线沿着z=0.1m的直线在地表上方对地下区域进行合成孔径扫描，整个扫描孔径从x=0.2m到x=1.98m，各孔径间距为Δx=0.02m，一共有90个孔径点。在每个孔径点处，发射天线向下发射电磁波，接收天线接收从地下散射的回波信号。随着收发一体天线从左到右移动扫描，一共可以获得90道散射数据，如图2所示。图2中的原始散射数据经过去直达波的预处理后的数据如图3所示。下面分别运用传统的加权后向投影算法和加窗加权后向投影算法进行成像运算。图4给出了传统的加权后向投影算法中的时延曲线处散射回波的提取示意图。图中分别画出了成像区域中三个成像单元所对应的三条时延曲线，占据了整个扫描孔径。从图中可以清楚地看出，真实目标的时延曲线仅占据了某个孔径区间，而成像时是采用了整个时延曲线上的散射回波。当探测区域中存在多目标时，各目标的绕射曲线就会产生部分重叠。如果采用整个孔径上时延曲线的散射回波进行成像处理的话，将造成成像质量的劣化。图5给出了传统的加权后向投影成像算法的成像结果。从图中可见，除了对7个目标的正确成像外，引入了较强的旁瓣干扰。而本专利提出的加窗加权后向投影成像算法仅截取整个时延曲线上的部分长度进行成像运算。成像区域中三个成像单元对应的截断时延曲线如图6所示。本实例中，根据探地雷达天线的波瓣宽度和成像单元的深度，截取窗口宽度W_a统一取为31。运用加窗加权后向投影成像算法对图3所示的预处理后数据进行成像处理，结果如图7所示。从图中可见，成像质量明显提高，旁瓣干扰得以有效抑制。

Claims (2)

1.一种探地雷达加窗加权后向投影成像方法，其特征在于，将整个成像区域划分为L_z×L_x个网格，一维测线沿地表布置，一维测线方向设为x方向，该一维测线上共有L个测点，坐标分别为x_i，i＝1,…,L；通过合成孔径天线对探测区域的扫描；加窗加权后向投影成像的方法包括以下步骤：

步骤1：提取各成像单元对应的时延曲线上散射回波数据；

对成像区域中的每一网格位置(z_m,x_n)，m=1,…,L_z;n=1,…,L_x，根据探测扫描场景计算第i个孔径测点x_i处对应的电磁波双程传播时延τ_m，n，i，提取该测点x_i处的一维回波信号s_i′(t)在时刻τ_m，n，i，生成一维信号w_m,n=[q_m，n，1,…,q_m,n，L]；

步骤2：计算有效扫描宽度W_a并截取数据：

W_a的计算公式为

其中，z_m为成像单元的深度，θ₁表示探地雷达天线的零功率波瓣宽度；round(·)表示四舍五入取整；

截取数据是指截取一维信号w_m,n中长度为W_a的一段数据

作为成像单元(z_m,x_n)对应的时延曲线加窗后的散射回波数据；

所截取的这一段散射数据

为：

式中，

二者满足W_l+W_r+1=W_a，其中fix(·)表示截尾取整；

步骤3：对成像单元(z_m,x_n)对应的时延曲线加窗范围中的散射数据

进行加权求和作为该成像单元的散射强度值，即

$O(z_m,x_n)=\mathrm{\&alpha;}(z_m,x_n)\mathrm{\&Sigma;}{w}_{m,n}^c=\mathrm{\&alpha;}(z_m,x_n)\underset{l=l_1}{\overset{l_2}{\mathrm{\&Sigma;}}}{q}_{m,n,l};$

式中l₁和l₂分别表示对一维信号w_m,n=[q_m，n，1,…,q_m，n，L]进行加窗截取时，窗口的左端和右端的序号值；加权因子α(z_m,x_n)的计算方法如下：

其中，e_lf、e_mid和e_rg分别表示将加窗后的散射回波矢量平均分为三段后各段散射数据的能量值，计算方法为

$e_{\mathrm{rg}}=\underset{l=b+1}{\overset{W_a}{\mathrm{\&Sigma;}}}{\left[w_{m,n}^c\left(l\right)\right]}^2,$ 其中 $a=\mathrm{fix}\left(\frac{W_a}{3}\right),$ $b=\mathrm{fix}\left(\frac{{2W}_a}{3}\right)$ 分别表示将长度为W_a的加窗后散射回波矢量

平均分为三段时各段的分界点；

${\mathrm{\&alpha;}}^{\&prime;}=\left\{\begin{array}{cc}1,& s=0\\ \frac{u}{s},& s\&NotEqual;0\end{array}\right.,$ 其中s和u分别为该成像点对应的时延曲线加窗后的散射回波矢量

的标准差和均值；

按照前述的步骤1-步骤3遍历成像区域中所有的网格，分别计算各点的散射强度O(z_m,x_n)即可获得整个成像区域的成像结果O(z_m,x_n)m=1,…,L_z;n=1,…,L_x。

2.根据权利要求1所述的探地雷达加窗加权后向投影成像方法，其特征在于，L取值范围为[3,10000]。

Images