CN102621548A - 一种探地雷达多尺度后向投影成像方法 - Google Patents

Abstract

本发明公开了一种探地雷达多尺度后向投影成像方法，步骤1：设置成像区域，设定第一尺度的成像网格，获得第一尺度的二维成像结果；尺度号i＝1；步骤2：对第i尺度下的二维成像结果进行局部极大值的提取，确定处于局部极大值的成像区域；步骤3：设定第i+1尺度的成像网格，对处于局部极大值区域的细成像网格，再次运用后向投影成像方法进行成像处理；步骤4：对非局部极大值区域的细成像网格进行成像结果的填充；步骤5：判断细成像网格是否满足设置的条件：如不满足，则使i递增1，并返回步骤2；如满足，则成像结束。该探地雷达多尺度后向投影成像方法在成像质量和精度相当的情况下，能显著提高成像效率。

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】。成像处理前需要对原始数据进行均衡、解振荡、去直达波、零点校正等预处理【参考文献：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)处的散射强度为 $O(z_m,x_n)=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{q}_{m,n,l}.$

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

传统的后向投影成像方法将成像区域进行单次均匀划分，再分别计算各成像网格的散射强度值。为提高成像质量，通常将成像网格划分的很密。L_z和L_x的取值一般要满足如下关系式： $\left\{\begin{array}{c}{L}_z\≥K\\ L_x\≥L\end{array}\right..$

对探地雷达扫描探测来说，扫描区域中的目标分布通常是稀疏的。即在整个成像区域中，目标仅占据很少的空域部分。经过对成像区域的划分，仅有很少的部分网格代表了目标。传统的后向投影成像方法将所有的网格依次进行计算获得成像结果。在非目标区域的大部分网格上，也进行了复杂的“延迟——累加”运算。实际上，在这些区域中，无需进行精细的网格划分和成像。因此，传统的后向投影成像方法的成像效率是很低的，不能有效应用于探地雷达实时成像场合。

发明内容

本发明所要解决的技术问题是提供一种探地雷达多尺度后向投影成像方法，该探地雷达多尺度后向投影成像方法在成像质量和精度相当的情况下，能显著提高成像效率。

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

一种探地雷达多尺度后向投影成像方法，包括以下步骤：

步骤1：设置成像区域，设定第一尺度的成像网格，运用后向投影成像方法进行成像处理，获得第一尺度的二维成像结果；尺度号i＝1；

步骤2：对第i尺度下的二维成像结果进行局部极大值的提取，确定处于局部极大值的成像区域；

步骤3：保持整个成像区域不变，设定第i+1尺度的成像网格，第i+1尺度的成像网格的尺寸小于第i尺度的成像网格，记第i尺度的成像网格和第i+1尺度的成像网格分别为粗成像网格和细成像网格；对处于局部极大值区域的细成像网格，再次运用后向投影成像方法进行成像处理；

步骤4：对非局部极大值区域的细成像网格进行成像结果的填充：这些网格处的成像值直接由与之对应的粗成像网格处的成像结果给出；从而获得第i+1尺度的二维成像结果；

步骤5：判断细成像网格是否满足设置的条件：

如不满足，则使i递增1，即i＝i+1，并返回步骤2；

如满足，则成像结束。

探地雷达进行一维沿线扫描，所获得的雷达记录剖面记为E₀(x，t)；E₀表示散射场的电场强度值，x是检测点的横向位置，t表示时间；该雷达记录剖面为一个二维矩阵，其尺寸为K×L，K表示时间维的采样点数，L表示空间维的横向采样点数；成像的目的便是通过探地雷达记录剖面数据E₀(x，t)获得整个探测区域的散射强度值；根据探地雷达空间扫描区域，确定横向维成像范围为[h_a，h_b]，纵向维成像范围为[g_a，g_b]；

第一尺度的成像网格的尺寸按照以下方法确定：

将该探地雷达空间扫描区域分别沿纵向深度维和横向维等间隔地划分为

和

个网格；

和

按照如下方式加以确定： $\left\{\begin{array}{c}{L}_{z_1}=\mathrm{round}\left(\frac{1}{M}K\right)\\ L_{x_1}=\mathrm{round}\left(\frac{1}{N}L\right)\end{array}\right.;$ 其中round(·)函数表示取整运算，M和N取大于2的整数，同时保证 $\left\{\begin{array}{c}{L}_{z_1}\≥3\\ L_{x_1}\≥3\end{array}\right..$

步骤2的具体过程为：

提取第i尺度下的网格成像结果O₁(z_m，x_n)的最大值和最小值，分别记为

和分以下三种情况提取局部极大值对应的成像区域：【这里是对整个成像区域的处理，即在整个成像区域中，提取满足一定条件的某些成像区域，作为下一步进行更细网格成像的区域。】

(1)如果 $O_{1_{\max }}>0,O_{1_{\min }}<0,$

则设置阈值系数S₁和S₂，分别提取满足式

和式

的所有成像区域为局部极大值区域；

(2)如果 $O_{1_{\max }}>0,O_{1_{\min }}\&GreaterEqual;0,$

则设置阈值系数S₁，提取满足式

的所有成像区域为局部极大值区域；

(3)如果 $O_{1_{\max }}\&le;0,O_{1_{\min }}<0,$

则设置阈值系数S₂，提取满足式

的所有成像区域为局部极大值区域；

阈值系数S₁和S₂应满足 $\left\{\begin{array}{c}{S}_1\&Element;\left(\mathrm{0.3,0.9}\right)\\ S_2\&Element;\left(\mathrm{0.3,0.9}\right)\end{array}\right.;$ 局部极大值区域占据了部分成像区域，记为{(m，n)|z_m∈P；x_n∈Q}；其中，P和Q分别为纵向维成像范围[g_a，g_b]和横向维成像范围[h_a，h_b]的子集，分别表示了第i尺度下的成像结果局部极大值所对应的纵向和横向范围。【P和Q可分别表示为多个纵向和横向区间的并集，即在整个成像区域中，所提取的局部极大值区域，可能在横向上和纵向上体现为多段不连续的区域。因此，一般表示方式下，可以表示为多个纵向和横向区间的并集。】

步骤5中所述的设置的条件为L_xf≥L，其中L_xf表示当前小尺度细网格成像结果的横向维网格划分数目。

步骤3中，第i+1尺度的网格数量为第i尺度的网格数量的G₁×G₁倍，G₁为细化系数，G₁的取值为G₁∈[2，7]；在细化的过程中，整个成像区域的空间范围不变。

细成像网格按照成像坐标的位置可以分为两部分，分别为{(m，n)|z_m∈P；x_n∈Q}和

。下面分别获得这两部分范围内各成像网格的散射强度值，从而完成细网格的成像处理。

(1){(m，n)|z_m∈P；x_n∈Q}网格处的成像

运用后向投影成像方法，遍历所有满足{(m，n)|z_m∈P；x_n∈Q}的细网格，获得整个区域{(m，n)|z_m∈P；x_n∈Q}的成像结果，记为O₂(z_m，x_n)，m∈P；n∈Q。(2)

网格处的成像

在粗成像网格成像结果的非局部极大值区域中，粗网格同样被细化为了细网格。粗成像网格的某个网格单元，细化为了G₁×G₁个网格单元。这G₁×G₁个网格单元的成像结果(散射强度值)可以直接用粗网格中某个网格单元的成像结果来填充。在粗成像网格中，遍历所有的非局部极大值区域，可以实现对细成像网格中所有非局部极大值区域的成像填充，该区域的成像结果表示为O₂(z_m，x_n)，

。

对上述两部分的成像结果进行空域组合，即可获得细成像网格的成像结果，记为O₂(z_m，x_n)，

即：对局部极大值区域的细成像网格，进行后向投影成像获得散射强度值；对其他区域的细成像网格，直接将粗成像网格的成像结果进行填充。

本发明的构思是：探地雷达对某个探测区域进行探测成像时，目标通常仅位于若干个有限区域范围内。因此，可以仅对这些有限区域范围进行精细成像，而对目标范围外的其他区域进行粗略成像，无需对整个成像区域进行精细成像。因此，可以通过多尺度成像方法大幅降低运算量。首先对整个成像区域进行大尺度粗网格成像，获得目标的大致位置。然后对目标的大致位置范围进行小尺度精细网格划分并成像，其他区域的小尺度成像结果直接用相应位置处的大尺度成像结果填充。

有益效果：

本发明提出了一种用于探地雷达的多尺度后向投影成像方法。其特点在于：对某个探测区域进行成像处理时，先进行大尺度(粗网格)的成像，确定目标的大致分布位置；然后仅对这些大致分布位置进行小尺度(细网格)精细成像，对其他的非目标区域不再进行精细成像。该成像方法一方面保证了目标区域范围内的精细成像，另一方面大大降低了整个成像过程的运算量。实验结果显示，本发明提出的方法与现有的成像方法相比，成像效率显著提高。

附图说明

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

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

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

图4示出了大尺度粗网格(10×10)的后向投影成像结果；

图5示出了大尺度成像结果提取的局部极大值区域(S₁＝0.4)；

图6示出了小尺度细网格(40×40)的成像结果(G₁＝4)；

图7示出了小尺度成像结果提取的局部极大值区域(S₁′＝0.5)；

图8示出了更小尺度(120×120)的成像结果(G₂＝3)；

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

图10示出了多个点目标的空间分布位置；

图11示出了原始雷达记录剖面经过预处理后的结果；

图12示出了大尺度粗网格(15×15)的后向投影成像结果；

图13示出了大尺度成像结果提取的局部极大值区域(S₁＝0.3)；

图14示出了小尺度细网格(90×90)的成像结果(G₁＝6)；

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

标号说明：1-发射天线，2-接收阵列天线。

具体实施方式

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

实施例1：

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

首先运用探地雷达对探测区域进行一维扫描，如图1所示。发射天线固定在x＝1m、z＝0m处，接收阵列放置在z＝0m处，其横向扩展为[0.2 1.8]m，两相邻阵元间距0.02m，一共有81个阵元。探测区域如图1所示，在x＝1m、z＝0.6m处放置一个半径为0.02m的细长型铁棒，长度远大于其半径，导体棒的取向垂直于纸面。发射天线为线源形式，向地下发射电磁波。散射回波被阵列天线接收，共81道数据，如图2所示。图2中的原始散射数据经过去直达波的预处理后的数据如图3所示。下面运用多尺度后向投影成像方法获得该探测区域的二维成像结果。确定横向成像区域为[0.2，1.8]m，纵向成像区域为[0，1.3]m。设置大尺度(对于第i尺度和第i+1尺度，将前者称为大尺度，后者称为小尺度，后同。)粗网格，取运用后向投影成像方法对该大尺度粗网格进行成像，成像结果如图4所示。按照步骤3中的第二种情况，设置阈值系数S₁＝0.4，提取目标局部极大值区域，如图5所示。然后，设置细化系数G₁＝4，对大尺度粗网格进行细化处理。再次运用后向投影成像算法对局部极大值处的细网格进行成像运算，其他区域的网格直接进行成像结果的填充，获得小尺度成像结果，如图6所示。此时，小尺度成像结果不满足退出条件，因此继续进行成像运算。按照步骤3中的第二种情况，设置阈值系数S₁′＝0.5，提取目标局部极大值区域，如图7所示。然后，设置细化系数G₂＝3，对小尺度细网格再次进行细化处理。同样地，对局部极大值处的细网格进行后向投影运算，对其他区域的细网格直接进行成像结果的填充，获得更小尺度成像结果，如图8所示。该成像结果的空间网格数目满足退出条件，因此此处的成像结果即为多尺度后向投影算法的最终结果。在CPU为Pentium(R)Dual-Core 2.00GHz，内存为2.00GB的PC上运行，该成像算法共耗时7.99s。

此处同时给出运用传统的后向投影成像算法对原始散射数据的成像处理。运用传统的后向投影算法进行成像处理时，首先确定成像区域，区域设置同上，横向成像区域为[0.2，1.8]m，纵向成像区域为[0，1.3]m。然后进行成像网格划分。为了与上述的多尺度后向投影成像方法进行对比，此处将成像区域的横向维和纵向维分别离散化为和

个网格，即将整个成像区域划分为120×120个网格。然后运用后向投影成像方法对各网格逐一进行成像运算，所获得的成像结果如图9所示。在相同的计算机配置条件下，该算法耗时504.7s。对比图8和图9可见，在目标位置附近，成像结果都获得了目标的散射强度值和目标的精确位置。在非目标区域，多尺度后向投影算法仅进行了大尺度粗网格的成像。而传统的后向投影成像算法对所有的成像网格都进行了精细成像，在非目标区域范围内浪费了大量的处理时间。而这些区域实际上并不需要进行精细成像。运用多尺度后向投影成像方法，运算时间是传统后向投影成像算法的1.58％，成像结果完全能够满足对目标定位和精细成像的要求。

实施例2：

探地雷达的天线配置同实施例1。探测区域如图10所示，在三个位置处放置细长型铁棒，长度远大于其半径，导体棒的取向垂直于纸面。发射天线为线源形式，向地下发射电磁波。散射回波被阵列天线接收，共81道数据，原始记录剖面经过预处理后的数据如图11所示。成像区域的粗网格设置为 $\left\{\begin{array}{c}{L}_{z_1}=15\\ L_{x_1}=15\end{array}\right.,$ 大尺度成像结果如图12所示。根据图12的大尺度成像结果，进行局部极大值提取，对应于前述的步骤2的第二种情况，取S₁＝0.3，局部极大值对应的区域如图13中的黑色区域所示。进行网格细化，对应于前述的步骤4，细化系数G₁取值为6。对局部极大值区域的细网格进行后向投影成像，对其他区域的网格进行成像结果的填充，得到细网格的成像结果如图14所示。细网格的成像结果满足退出条件，多尺度后向投影成像到此结束。在与实例1相同的计算机配置条件下，该成像算法所需的时间为20.9s。

运用传统的后向投影成像算法对相同的记录剖面进行处理，直接对成像区域进行精细网格划分，即将成像区域划分为(15×6)×(15×6)个网格，成像结果如图15所示。在相同的计算机配置情况下，所需的时间为248.5s。通过对比图14和图15可见，多尺度后向投影成像方法在目标所在位置周围进行了精细成像，与传统后向投影成像方法的成像结果相同。而在非目标区域，本专利所提的成像方法运用大尺度网格代替了小尺度网格，有效降低了运算量。成像结果保证了目标区域的分辨力，耗时降低为原成像算法的8.41％。

Claims (5)

1.一种探地雷达多尺度后向投影成像方法，其特征在于，包括以下步骤：

步骤1：设置成像区域，设定第一尺度的成像网格，运用后向投影成像方法进行成像处理，获得第一尺度的二维成像结果；尺度号i＝1；

步骤2：对第i尺度下的二维成像结果进行局部极大值的提取，确定处于局部极大值的成像区域；

步骤3：保持整个成像区域不变，设定第i+1尺度的成像网格，第i+1尺度的成像网格的尺寸小于第i尺度的成像网格，记第i尺度的成像网格和第i+1尺度的成像网格分别为粗成像网格和细成像网格；对处于局部极大值区域的细成像网格，再次运用后向投影成像方法进行成像处理；

步骤4：对非局部极大值区域的细成像网格进行成像结果的填充：这些网格处的成像值直接由与之对应的粗成像网格处的成像结果给出；从而获得第i+1尺度的二维成像结果；

步骤5：判断细成像网格是否满足设置的条件：

如不满足，则使i递增1，即i＝i+1，并返回步骤2；

如满足，则成像结束。

2.根据权利要求1所述的探地雷达多尺度后向投影成像方法，其特征在于，探地雷达进行一维沿线扫描，所获得的雷达记录剖面记为E₀(x，t)；E₀表示散射场的电场强度值，x是检测点的横向位置，t表示时间；该雷达记录剖面为一个二维矩阵，其尺寸为K×L，K表示时间维的采样点数，L表示空间维的横向采样点数；成像的目的便是通过探地雷达记录剖面数据E₀(x，t)获得整个探测区域的散射强度值；根据探地雷达空间扫描区域，确定横向维成像范围为[h_a，h_b]，纵向维成像范围为[g_a，g_b]；

第一尺度的成像网格的尺寸按照以下方法确定：

将该探地雷达空间扫描区域分别沿纵向深度维和横向维等间隔地划分为

和个网格；

和

按照如下方式加以确定： $\left\{\begin{array}{c}{L}_{z_1}=\mathrm{round}\left(\frac{1}{M}K\right)\\ L_{x_1}=\mathrm{round}\left(\frac{1}{N}L\right)\end{array}\right.;$ 其中round(·)函数表示取整运算，M和N取大于2的整数，同时保证 $\left\{\begin{array}{c}{L}_{z_1}\&GreaterEqual;3\\ L_{x_1}\&GreaterEqual;3\end{array}\right..$

3.根据权利要求1所述的探地雷达多尺度后向投影成像方法，其特征在于，步骤2的具体过程为：

提取第i尺度下的网格成像结果O₁(z_m，x_n)的最大值和最小值，分别记为

和

分以下三种情况提取局部极大值对应的成像区域：

(1)如果 $O_{1_{\max }}>0,O_{1_{\min }}<0,$

则设置阈值系数S₁和S₂，分别提取满足式

和式

的所有成像区域为局部极大值区域；

(2)如果 $O_{1_{\max }}>0,O_{1_{\min }}\&GreaterEqual;0,$

则设置阈值系数S₁，提取满足式

的所有成像区域为局部极大值区域；

(3)如果 $O_{1_{\max }}\&le;0,O_{1_{\min }}<0,$

则设置阈值系数S₂，提取满足式

的所有成像区域为局部极大值区域；

阈值系数S₁和S₂应满足 $\left\{\begin{array}{c}{S}_1\&Element;\left(\mathrm{0.3,0.9}\right)\\ S_2\&Element;\left(\mathrm{0.3,0.9}\right)\end{array}\right.;$ 局部极大值区域占据了部分成像区域，记为{(m，n)|z_m∈P；x_n∈Q}；其中，P和Q分别为纵向维成像范围[g_a，g_b]和横向维成像范围[h_a，h_b]的子集，分别表示了第i尺度下的成像结果局部极大值所对应的纵向和横向范围。

4.根据权利要求2所述的探地雷达多尺度后向投影成像方法，其特征在于，步骤5中所述的设置的条件为L_xf≥L，其中L_xf表示当前小尺度细网格成像结果的横向维网格划分数目。

5.根据权利要求1-4任一项所述的探地雷达多尺度后向投影成像方法，其特征在于，步骤3中，第i+1尺度的网格数量为第i尺度的网格数量的G₁×G₁倍，G₁为细化系数，G₁的取值为G₁∈[2，7]；在细化的过程中，整个成像区域的空间范围不变。

Images