CN101915943A - 均匀背景介质的介电常数和隐蔽目标参数的联合反演方法 - Google Patents

Abstract

本发明提出了一种均匀背景介质的介电常数和隐蔽目标参数的联合反演方法，实现方式是首先设定一个介电常数的估计区间，然后对该区间中的每个估计值计算生成字典矩阵，再将这些字典矩阵按行排列生成联合字典矩阵。以此为基础运用基于压缩感知的成像算法进行优化求解，获得联合目标矩阵。再对联合目标矩阵进行分割提取，同时获得背景介质的介电常数值和探测区域的成像结果，从而较大程度提高了探地雷达的实时探测成像能力。

Description

均匀背景介质的介电常数和隐蔽目标参数的联合反演方法

技术领域

本发明属于探地雷达的高分辨成像技术领域，涉及一种均匀背景介质的介电常数和隐蔽目标参数的联合反演方法。

技术背景

探地雷达是一种有效的无损探测技术。它通过空域扫描向探测区域发射电磁波并接收散射回波，可实现对未知区域内部的成像处理，获得未知区域中的隐蔽目标参数，即目标分布信息和散射强度信息，有效应用于市政工程、考古、地雷探测、反恐等多种场合。探地雷达的成像方法有多种，基于压缩感知技术的成像算法基于回波信号的稀疏性、测量矩阵的随机性和优化算法完成对探测区域的高分辨成像处理，很大程度上放宽了探地雷达成像对空间维和时间维高密度采样的要求。以二维均匀无耗介质区域中的点散射型目标为例，对给定的成像区域x∈[x₁，x₂]；z∈[z₁，z₂]，设定成像分辨单元为δ_xz，其x方向和z方向的网格数目分别记为L_x和L_z，则可将目标函数表征为矢量形式

O＝[O_2D(1：L_z，1)；O_2D(1：L_z，2)；L；O_2D(1：L_z，L_x)]  (1)

该矢量的长度为L_o＝L_xL_z。其中O_2D为二维目标散射强度的矩阵描述形式，O_2D(m，n)＝σ(m+(n-1)×L_z)，m∈[1，L_z]，n∈[1，L_x]，σ(q)表示第q个分辨单元的散射强度值。若该处无目标，则σ(q)＝0。对该区域沿x方向进行空域扫描即可获得多个孔径点处的散射回波，进而可以进行成像处理。空域采样道数记为L_a，每道数据的时间采样点数记为N_t。当均匀背景介质的介电常数ε_d已知时，可以根据电磁散射机理建立从目标函数矢量到散射数据矢量之间的字典矩阵，记为Ψ，其尺寸为L_aN_t×L_o。该矩阵由L_a个子矩阵按列排列，可记为：

$\mathrm{\Ψ}=[{\mathrm{\Ψ}}_1;L;{\mathrm{\Ψ}}_s;L;{\mathrm{\Ψ}}_{L_a}]---\left(2\right)$

该矩阵中，第s个孔径位置处对应的子阵记为Ψ_s，其尺寸为N_t×L_o。根据电磁散射理论，该子阵的元素值可表示为

${\mathrm{\Ψ}}_s={\left[{\mathrm{\Ψ}}_s\right]}_{N_t,q}=\frac{s(t_n-\mathrm{\τ}(s,q))}{{\left|\right|s(t_n-\mathrm{\τ}(s,q))\left|\right|}^2}---\left(3\right)$

其中，

d(s，q)表示第s个孔径位置处到第q个分辨单元处的空间距离，

为真空中的介电常数值，c为光速。根据压缩感知理论，当选择的随机采样矩阵Φ满足ΦΨ具有有限等距性质这一条件时，就可以通过优化算法求解如下方程实现对目标矢量O的求解

$\hat{O}=\arg \min {\left|\right|O\left|\right|}_1$ s.t.x_d＝ΦΨO    (4)

其中x_d表示目标回波随机采样向量。

对探测区域的成像处理需要预先知道均匀无耗背景介质的介电常数ε_d。若介电常数未知，则无法进行区域的成像处理。现有的处理方法是设定一个介电常数估计区间，分别运用该区间中的每一个介电常数估计值进行成像处理，再根据成像结果的聚焦程度来确定介质的介电常数。这种处理方法是基于分步反演的，运算量大，估计精度不高，无法实现介质电磁参数和目标参数的联合反演，制约了探地雷达的应用。

发明内容

为了克服现有的介质介电常数和隐蔽目标参数分步反演的不足，提高探地雷达对未知介质中目标的成像反演能力，本发明提供了一种均匀背景介质的介电常数和隐蔽目标参数的联合反演方法。以探地雷达的压缩感知成像算法为基础，通过设计生成联合字典矩阵，建立了散射回波所要满足的约束方程。然后基于优化算法对该约束方程进行求解，即可同时获得介质的电磁参数和目标的参数信息，大大提高了探地雷达对电磁参数未知介质中目标的成像反演能力。

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

一种均匀背景介质的介电常数和隐蔽目标参数的联合反演方法，包括以下步骤：

步骤1：设定初始介电常数估计矢量：根据预估的介电常数区间设定一个初始介电常数估计矢量

初始介电常数估计矢量中的D个值均匀分布；

步骤2：生成联合字典矩阵：对初始介电常数估计矢量中的每个估计值ε_m，m∈[1，D]分别计算生成字典矩阵

m∈[1，D]；然后将该D个字典矩阵按行排列生成联合字典矩阵

联合字典矩阵Ψ_JOINT对应的目标函数矢量为O_JOINT，目标函数矢量的形式为：

步骤3：求解目标函数的估计值：根据

s.t.x＝ΦΨ_JOINTO_JOINT求取目标函数的估计值

x表示原始的探地雷达记录剖面经过随机采样后的数据；其中Φ为选择的随机采样矩阵，且满足ΦΨ_JOINT具有有限等距性质这一条件；的形式为：

步骤4：获得探测区域的成像结果和该区域的介电常数值：由于

的第s个矢量

为非零矢量，由此反推，初始介电常数估计矢量

中对应的第s个数值ε_s即为真实的背景介质的介电常数值；

中仅有的一个非零子矩阵

即为探测区域的成像结果。

D值为4～7之间的整数。

原始的探地雷达记录剖面经过随机采样后的数据如图2所示，横向维是空间维，即空间的各孔径采样点；纵向维是时间维，即回波数据的时间维变量。在运用压缩感知理论进行处理时，通过选定合适的随机采样矩阵Φ对原始的探地雷达剖面进行处理，处理后的剖面图如图9所示。图9所示的二维矩阵即为此处的x，它包含了目标函数矢量的解联合反演就是从图9所示的二维数据中按照本专利提出的方法对背景介质的介电常数和隐蔽目标参数(即散射点在二维空间的位置参数和各散射点的散射强度参数)进行联合反演。

即为约束方程s.t.x＝ΦΨ_JOINTO_JOINT的解，当然也就是我们要得到的目标函数的估计值。本专利中，先建立起约束方程，目的就是求得其中的

求得该矩阵后，也就实现了联合反演，也就获得了目标函数的估计值。

表示的是矢量，是为了便于构造联合字典矩阵。此处将其称为“非零子矩阵”，可以视其为是多行一列的矩阵。实际上表示的是一个二维空间的目标分布情况，即所需的成像结果。

是

的估计值，其实质也是成像结果的模拟。对一个二维成像结果而言，可表示为一个二维矩阵。横向和纵向分别按照分辨单元进行网格划分。当某个网格处没有目标存在时，对应的二维矩阵中该位置处就是零。有目标存在时，这个位置处就有值，值的大小即为该目标的散射强度值。本专利中，为了便于构造联合字典矩阵，将这样一个二维矩阵按列排列，即如式(1)的表述。采用本专利的方法实现了联合反演后，实际上获得的是一个按列排列的矢量。再将这个矢量转换为二维矩阵，就变成了二维成像结果了。

有益效果：

基于压缩感知的基本理论，设计生成目标函数矩阵和联合字典矩阵，通过对约束方程的优化求解，可以实现均匀背景介质的介电常数和探测区域中目标参数的联合反演，很大程度上降低了分步反演的运算量，提高了探地雷达对探测区域进行成像探测的实时处理能力。

探地雷达对某个探测区域进行扫描探测，可以获得一个雷达记录剖面。若要获得该探测区域的成像结果，就需要对记录剖面进行成像处理。而成像处理需要预知背景介质的介电常数值。当背景介质的介电常数值未知时，就需要对其进行估计。介电常数的估计精度会影响到成像结果的聚焦质量。传统的处理方法中，有基于原始记录剖面中的衍射曲线的形式进行估计的，也有基于成像运算进行估计的。前者的估计是实时估计，但精度较差。后者的估计精度较高，但运算量很大。基于成像运算进行估计时，首先是设定一个介电常数估计矢量，包括多个可能的介电常数值，然后运用每个介电常数估计值进行一次成像运算，再对多个成像结果分别计算其聚焦值，选择聚焦值最大(聚焦效果最好)的那个成像结果作为探测区域的成像结果，对应的介电常数值作为背景介质的介电常数估计(见参考文献2：修志杰，陈洁，方广有，李芳。基于F-K偏移及最小熵技术的探地雷达成像法，电子与信息学报，2007，29(4)：827-830)。这种方法的运算量很大，针对每个介电常数都要进行一次成像运算，然后计算各自成像结果的聚焦值。这种方法只能离线运行，通过事后处理获得成像结果和介电常数估计。

本发明基于压缩感知的基本理论，仅需若干个稀疏的时间采样和空间采样值，即可运用优化求解的方法实现背景介质的介电常数值和探测区域的成像结果的联合反演，即运用约束方程在设计生成的联合目标矩阵中对成像结果和背景介质的介电常数值进行联合寻优估计。该方法避开了原始处理方法中对每个介电常数估计值分别进行成像运算和成像结果的聚焦值运算的大运算量处理，极大程度降低了联合反演的运算量，可以在线实时运行，实时获得成像结果和介电常数估计。

附图说明

图1示出了二维成像区域中的三个目标；图中，横坐标表示横向维的空间采样矢量，纵坐标表示纵向维(也就是深度维)的空间采样矢量。右边的竖直坐标是颜色标尺，表示了图中所示的颜色对应的数值。此图中，从右边的标尺可见，黑色对应0值，浅灰色对应0.25，白色对应0.5。说明图中仅有三个分辨单元非零，其值分别为0.25，0.25和0.5。此处，该值实际上表示的是各散射点的散射强度值。

图2示出了图1所示的目标的空域散射回波；

图3示出了图1所示的成像区域对应的联合字典矩阵；

图4示出了图3所示的联合字典矩阵中某一介电常数估计值对应的子矩阵；

图5示出了图4所示的子矩阵中某一空间孔径点对应的子矩阵；

图6示出了图5所示的子矩阵中空间孔径点对应某一分辨单元的传递函数；

图7示出了所有空间孔径点对应的随机采样矩阵；

图8示出了图7所示的随机采样矩阵中某一空间孔径点对应的子矩阵；

图9示出了采用图7所示的随机采样矩阵对图2所示的散射回波进行随机采样后的散射数据；

图10示出了运用优化算法获得的联合目标矩阵；

图11示出了图10所示的联合目标矩阵进行分割后的真实成像结果。

具体实施方式

以下将结合图和具体实施过程对本发明做进一步详细说明。

首先根据预估的介电常数区间设定一个介电常数的估计矢量

此处的估值区间是根据原始雷达记录剖面中目标的衍射曲线形式和探测环境介电常数的经验值确定的。该矢量的数目，即D的值一般可取为4～7个。在该估计区间中，D个值是均匀分布的。数目越多，后续的运算量也越大，估计值也越精确，一般取为4到7个即可；然后对估计矢量中的每个估计值ε_m，m∈[1，D]按照文献[1]中的方法分别计算生成字典矩阵

(文献1：AliC.Gurbuz，James H.McClellan，Waymond R.Scott Jr..Compressive sensing tor subsurface imaging using ground penetrating radar.Signal Processing，2009(89)：1959-1972：)；【该文献中，给出了当背景介质的介电常数已知的情况下的成像反演方法。此处，是对介电常数估计矢量中的某个估计值运用文献中的方法进行成像反演，文献中的方法可以直接用在此处。中文标题：Ali C.Gurbuz，James H.McClellan，Waymond R.Scott Jr..基于压缩感知的探地雷达地下成像.Signal Processing，2009(89)：1959-1972。】然后将这D个字典矩阵按行排列生成联合字典矩阵

最后基于压缩感知理论进行目标函数矢量O_JOINT的求解。

联合字典矩阵Ψ_JOINT的尺寸为L_aN_t×DL_o。L_a为空域采样道数记，N_t为每道数据的时间采样点数记，L_o＝L_xL_z，L_x和L_z分别是成像分辨单元在x方向和z方向的网格数目。介质参数的不同，导致字典矩阵中各分辨单元到各孔径点传递函数的不同，也就是时延值τ(s，q)的不同。以二维成像为例，是沿着一维空间方向对地下探测区域进行扫描的。扫描方式有多种，如合成孔径扫描、实孔径扫描、阵列天线扫描等。孔径点即为该一维扫描方向上的空间采样点，一般为均匀分布形式。为获得精细的成像结果，空间采样要满足奈奎斯特采样定理，即两个孔径点的间距要小于某个数值。空间采样越密，成像结果越好，当然所采集的数据量也就越大，数据处理的运算量也越大。

对某个区域进行成像探测，一般是首先确定所要探测目标的尺寸，也就是确定成像的分辨单元；然后根据该分辨单元确定空间采样的密度，即孔径点间距；再进行数据采集和成像处理。

二维成像中，整个图像可以用一个二维矩阵来表示。图像中各像素所占据的空间单元即为分辨单元，各分辨单元处各自的数值是不同的，分别表示了二维空间中各位置处的散射强度值。分辨单元是图像中能分辨相邻两个目标的最小单元。在一个分辨单元内的所有空域区间中，认为其散射强度是相同的。Ψ_JOINT矩阵中某一子矩阵

的尺寸为L_aN_t×L_o。该子矩阵又包括L_a个孔径点处各自对应的子矩阵中，第s个孔径位置处的子阵记为

其尺寸为N_t×L_o，Ψ_JOINT矩阵包括D个子阵，其中第m个子阵记为

而

矩阵包括了所有孔径点遍历所有分辨单元的信息。而孔径点一共有L_a个，因此

矩阵就包括了L_a个子阵。在该L_a个子阵中，第s个子阵

表示第s个孔径点遍历所有分辨单元的矩阵，其尺寸为N_t×L_o。上文中的“第s个孔径位置处”即为“第s个孔径点处”。

根据电磁散射理论，该子阵的元素值可表示为

${\mathrm{\Ψ}}_{s,{\mathrm{\ϵ}}_m}={\left[{\mathrm{\Ψ}}_{s,{\mathrm{\ϵ}}_m}\right]}_{N_t,q}=\frac{s(t_n-\mathrm{\τ}(s,q))}{{\left|\right|s(t_n-\mathrm{\τ}(s,q))\left|\right|}^2}---\left(5\right)$

其中，

d(s，q)表示第s个孔径点处到第q个分辨单元处的空间距离，孔径点是空间采样的一个采样点，分辨单元是图像中能够分辨相邻两个目标的最小单元。此处所说的孔径点到分辨单元的距离是指孔径点到分辨单元这个区域的中心的距离。c为光速，ε_m即为介电常数的估计矢量

中对应的介电常数值。O_JOINT矢量为D个尺寸为L_o×1(与矢量O相同)的矢量按列排列而成的，即：

其中仅有介电常数为真实值对应的子矩阵

非零，该子矩阵

即为式(1)的形式，表征了成像区域中的目标分布和散射系数信息。从以上的分析可见：介电常数的估计是隐含在O_JOINT这一矢量中的，目标位置矢量和散射强度矢量是包含在O_JOINT的某一非零子矩阵中的，而O_JOINT是典型的稀疏矢量。

根据压缩感知理论，当选择的随机采样矩阵Φ(对具有类似尖峰、正弦、子波、Gabor函数等信号形式的有限支集的

矩阵而言，具有独立同分布的Bernoulli随机变量或Gaussian随机变量都可以用来构造随机采样矩阵Φ，如文献1中所示。即随机采样矩阵的每个数值都是一个随机变量，各随机变量满足独立同分布的性质。)满足ΦΨ_JOINT具有有限等距性质这一条件时，就可以通过优化算法求解如下方程实现对目标矢量O_JOINT的求解

${\hat{O}}_{\mathrm{JOINT}}=\arg {\min \left|\right|O_{\mathrm{JOINT}}\left|\right|}_1$ s.t.x＝ΦΨ_JOINTO_JOINT    (7)

上式的意义是：在等式x＝ΦΨ_JOINTO_JOINT的约束下，选择具有最小1-范数的矢量O_JOINT的估计值。其中，argmin表示使目标函数取最小值时的变量值，‖O_JOINT‖₁表示矢量O_JOINT的1-范数。1-范数的定义如下：令矢量x＝(x₁，x₂，...，x_n)，则该矢量的1-范数为：‖x‖₁＝|x₁|+|x₂|+...+|x_n|。图9所示的即为本例中的x。

根据以上的联合矩阵的构造形式，上式的反演结果有如下形式：

此处的意义是：在公式6确定的模型下，运用该专利所提的方法进行估计，获得的估计值和原始模型的初始值相同。该公式也可表示为

即：在D个长度为L_o的矢量构成的联合矢量

中，仅有第s个矢量为真实的目标矢量O，其他D-1个矢量均为零。而初始介电常数估计矢量

中对应的第s个数值ε_s即为真实的背景介质的介电常数值。获得了目标矢量的估计值

后，该估计值中仅有的一个非零子矩阵

即为探测区域的成像结果。

本发明的最终结果是可以同时获得探测区域的成像结果和该区域的介电常数值。运用探地雷达对探测区域进行扫描探测可以获得一个记录剖面。某些应用中，是以获得该探测区域的成像结果为目的的。某些应用中，是以获得该背景介质的介电常数为目的的。还可能有其他的应用目的。不同的应用需求下，对记录剖面的处理方法也不同。本专利是针对隐蔽目标参数和背景介质的介电常数联合反演(隐蔽目标参数的反演也就是对雷达记录剖面的成像处理)这个问题而提出的。)而ε_d即为背景介质的介电常数值，这样就实现了均匀背景的介电常数和隐蔽目标参数的联合反演。

实施例1：

本实例是针对二维探测成像的，但本技术并不局限于二维区域，对三维区域的探测成像也适用。

首先设定一成像区域，设该区域的介电常数为16ε₀，将该区域沿横向和纵向分割为400个成像单元，横向和纵向的单元数目都为20。图1给出了该区域中三个散射点的空域分布情况。各散射点都为全向散射且互不影响，即散射回波满足线性叠加关系。各散射点的灰度值表示其散射强度值，这三个散射点的散射强度分别为0.25，0.25和0.5。对该区域进行成像探测时，空域采样为随机采样的形式，沿x轴横向进行采样，采样点数为30。在每个采样点处，向下发射电磁波并接收散射回波，时域采样点数为256。通过射线寻迹仿真计算，整个孔径的散射回波如图2所示。

实际探测时，通过空间维x方向采样即可获得类似图2所示的记录剖面。雷达成像便是对该记录剖面进行处理获得类似图1所示的目标分布和各目标的散射强度。当背景介质的介电常数已知时，可采用基于压缩感知理论的成像算法。而当背景介质的介电常数未知时，就需要采样本专利中的联合字典矩阵进行介质的介电常数和目标散射函数的联合反演。

首先定义介电常数的估计矢量为

该估计矢量是根据原始记录剖面的曲线形式以及对探测区域的经验判断选定估计区间。在估计区间中，均匀分布了5个数值，然后基于射线寻迹生成联合字典矩阵Ψ_JOINT，该矩阵的二维显示结果如图3所示。该矩阵有五个子矩阵按行排列，每个子矩阵对应于介电常数估计矢量中的一个介电常数估计值，此处给出介电常数估计值为12ε₀的子矩阵

的二维显示结果，如图4所示。该子矩阵

包含了所有横向孔径点与成像区域每一分辨单元间的传递函数，各孔径点对应的传递函数矩阵按列排列。此处给出第一个孔径点对应的块矩阵的二维显示结果，如图5所示。该矩阵共有400列，每一列表示了当前扫描点到任一分辨单元的传递函数，该传递函数是一时域信号，时间采样点数为256点。此处给出第一个孔径点到第一个分辨单元的传递函数形式，如图6所示。

根据压缩感知理论，选择随机采样矩阵Φ为服从N(0，1)分布的随机矩阵，如图7所示。本例以服从N(0，1)分布的随机矩阵为例，但本技术并不局限于这种形式的随机矩阵，只要满足ΦΨ_JOINT具有有限等距性质的随机矩阵Φ都适用。该矩阵为一广义对角阵，即每一孔径点处对应的子随机矩阵Φ_s，s∈[1，L_a]占据全随机矩阵对角线上的某一块。此处给出第一个孔径点对应的子随机矩阵Φ₁，如图8所示。该矩阵的大小为10×256，将原始散射回波(256点的时域采样)运用该矩阵转换为10点随机采样的时域数据。对每一孔径的散射回波都运用随机矩阵进行转换，得到随机采样的回波数据，如图9所示。基于压缩感知的雷达成像处理便是从如图9所示的随机采样回波数据中反演目标参数信息。通过建立联合字典矩阵，成像处理可转化为式(7)的优化问题进行求解。运用优化算法对式(7)进行迭代处理，可获得联合目标矩阵的估计值该估计值为一列向量，对其进行二维重排，结果如图10所示。从图中可见，共获得了五个介电常数值对应的二维成像结果，只有第三个介电常数值对应的成像结果非零，其他四个成像结果均为零。则估计矢量

中的第三个估计值即为真实的介电常数值，而对应的成像结果即为目标参数的反演结果，如图11所示。成像结果与图1所示的真实目标分布相同。

Claims (2)

1.一种均匀背景介质的介电常数和隐蔽目标参数的联合反演方法，其特征在于，包括以下步骤：

步骤1：设定初始介电常数估计矢量：根据预估的介电常数区间设定一个初始介电常数估计矢量

初始介电常数估计矢量中的D个值均匀分布；

步骤2：生成联合字典矩阵：对初始介电常数估计矢量中的每个估计值ε_m，m∈[1，D]分别计算生成字典矩阵

然后将该D个字典矩阵按行排列生成联合字典矩阵联合字典矩阵Ψ_JOINT对应的目标函数矢量为O_JOINT，目标函数矢量的形式为：

步骤3：求解目标函数的估计值：根据

s.t.x＝ΦΨ_JOINTO_JOINT求取目标函数的估计值

x表示原始的探地雷达记录剖面经过随机采样后的数据；其中Φ为选择的随机采样矩阵，且满足ΦΨ_JOINT具有有限等距性质这一条件；

的形式为：

步骤4：获得探测区域的成像结果和该区域的介电常数值：由于

的第s个矢量

为非零矢量，由此反推，初始介电常数估计矢量

中对应的第s个数值ε_s即为真实的背景介质的介电常数值；

中仅有的一个非零子矩阵即为探测区域的成像结果。

2.根据权利要求1所述的均匀背景介质的介电常数和隐蔽目标参数的联合反演方法，其特征在于，D值为4～7之间的整数。

Images