基于墙体补偿与多视角融合的穿墙雷达建筑物布局成像算法

尹子翔 杨小鹏 兰 天，2

（1.北京理工大学信息与电子学院雷达技术研究所，北京 100081；2.北京理工大学重庆创新中心，重庆 401120）

1 引言

近些年来，城市成为了灾害救援、反恐行动的主战场。这些作战往往需要进入密闭建筑物内部，如果能够在进入建筑物之前就能够获取到建筑物内部的结构信息，可以减少己方损失、提高作战行动成功率。因此，研究建筑布局成像，具有重要的现实意义和研究价值［1-2］。

对建筑布局成像，国内外许多科研机构都开展了相关研究。文献［3-4］利用边缘检测和Hough 变换在穿墙成像结果中突出显示建筑物墙壁，但是仅能显示墙体大致位置，没有对成像结果进行进一步处理。文献［5］提出了一种多视角建筑布局成像算法，使用M-N-K检测器融合成像，能够获取建筑物内部结构信息。文献［6］提出了一种基于多方位多尺度的建筑布局成像方法，使用非下采样Contourlet变换，将变换后的高频分量和低频分量采用相应的融合规则进行融合，可以增强建筑物墙体的轮廓和细节信息。以上两种算法完全基于图像处理的手段，并不考虑墙体参数，导致最终成像结果中墙体位置偏移、厚度展宽。文献［7］提出了一种适用于穿墙雷达的建筑布局成像算法，该算法通过Radon 变换提取墙体位置来进行墙体补偿，但是增强了多径鬼影等杂波的幅度，最终得到的建筑布局相对粗糙。

针对以上问题，本文提出了一种基于墙体补偿与多视角融合的建筑布局成像方法。在图像融合前进行墙体位置矫正，解决图像融合中出现的墙体展宽、位置偏移问题；在图像融合时采用最小值融合，能够在保留墙体的同时，抑制多径鬼影等杂波。首先，采用后向投影方法对多个视角的观测回波进行成像，对成像结果应用Radon变换，通过恒虚警率（CFAR）检测获取Radon平面峰值，提取墙体的位置信息；然后，以墙体的前后表面中心为分界线，对图像进行分割；进一步，对包含偏移墙体的子图像进行补偿，矫正墙体位置；最后，对多个视角墙体补偿后的图像进行最小值融合，抑制多径鬼影等杂波，得到正确清晰的建筑布局全景图像。

2 回波模型

对于一个穿墙雷达建筑布局成像场景，第k个合成孔径单元接收到的回波Yk(t)(k=1，2，…，L)由第i(i=1，2，…，N)堵墙的前表面反射波Yki，1(t)和后表面反射波Yki，2(t)组成。图1为三墙场景下电磁波传播示意图。在本文的成像场景中，成像目标为墙体而不是点目标，在成像中起主要作用的是垂直入射波的回波。因此可以忽略电磁波的折射，仅关注电磁波的垂直透射与反射。假设发射信号为S(t)，那么Yk(t)可以写为

其中，Ri为第i堵墙和第i-1堵墙之间的距离，R1为第一堵墙和收发天线之间的距离，Di为第i堵墙的厚度，c为光速，为电磁波在第i堵墙中的传播速度，其值为

其中，εi是第i堵墙的相对介电常数。

3 成像算法

后向投影法（Back Projection，BP）基本思想是首先对成像区域进行网格化，计算区域中像素点到收发天线的距离，进而得到传输时延，根据该时延寻址各个收发天线对应的回波信号并进行叠加处理。完成整个成像区域的遍历运算，得到最终的成像结果［8］。

以成像平面I(X，Y) 中任意一像素点Pq=(xq，yq)为例，该像素点对应的成像延时Tqk为

式中PT，k和PR，k(k=1，2，…，L)分别表示发射天线和接收天线的位置。将发射天线坐标PT，k=(xtk，ytk)，接收天线坐标PR，k=(xrk，yrk)代入公式（3），像素点Pq的成像聚焦时延可以写为

将该像素点Pq对应的所有回波信号叠加，得到像素点Pq的成像结果

最后遍历整个成像区域，得到最终的成像结果I(x，y)。

4 墙体补偿

4.1 墙体位置检测

实际应用中，墙体大多为平直的，利用这一特性采用直线检测方法检测原始图像I(X，Y)中的直线，获取墙体的位置信息，对偏移墙体进行矫正。本文采用Radon变换检测直线来获取墙体位置。

Radon变换定义如下

其中，I为成像平面，I(x，y)为BP 成像结果的像素值，ρ为直线到坐标原点的距离，θ为点(x，y)与x轴的夹角，δ为Dirac delta函数，其值为

这个函数使I(x，y)沿着直线

进行积分，如图2所示。

Radon 变换将图像投影到(ρ，θ)平面，(ρ，θ)平面中每一点对应(x，y)平面中一条直线的积分。因此，(x，y)中的墙体所在直线会在(ρ，θ)平面形成峰值，这就将原直线的检测转化为对峰值的检测。使用CFAR 检测得到峰值位置，提取峰值对应的(ρ，θ)值，得到(x，y)平面中该直线的斜率和与图像中心的斜距，从而确定该直线即墙体位置。

4.2 墙体位置矫正

由公式（4）计算得到的聚焦时延未考虑电磁波在墙体和在空气中传播的速度的差异，因此需在不同墙体及所在区域上矫正上述差异带来的位置偏移。

设有图1所示场景，墙体位置矫正的步骤如下：

1）对回波数据进行BP 成像得到原始图像I1(x，y)。

2）对I1(x，y)应用Radon 变换得到R1(ρ，θ)，通过CFAR 检测，得到的第一峰值和第二峰值分别对应I1(x，y)中第一堵墙的前表面和后表面。二者位置相减得到第一堵墙的成像厚度d1，其与真实厚度D1的关系为

3）以第一堵墙中心为分界线将图像I1(x，y)分割为两个子图像I11(x，y)和I12(x，y)。子图像I11(x，y)包含第一堵墙前表面及之前的部分，子图像I12(x，y)包含I1(x，y)中剩下的部分。其中，I11(x，y)中墙体位置正确，I12(x，y)中墙体位置偏移。

4）对第一堵墙造成的时延进行补偿，这个补偿时延为

补偿后I12(x，y)实际的成像结果应为

式（11）中对子图像每个像素点补偿一个常数时延Δτ1，等效于将该子图像向天线方向平移Δτ1·c=d1-D1的距离。因此无需对该子图像重新BP 成像，直接将该子图像向天线方向平移d1-D1，得到补偿后的图像I2(x，y)，即能达到与重新BP 成像相同的效果。对于公式（10）中d1可由步骤2）得到。D1为第一堵墙的真实厚度，可在真实场景中由相邻墙体厚度估计得到。通过这种方式，我们无需提前得知墙体的相对介电常数ε1便可对该墙体进行补偿。

5）对I2(x，y)应用Radon变换，得到第二堵墙的成像厚度d2，以第二堵墙中心为分界线将I2(x，y)分割为两个子图像I21(x，y) 和I22(x，y)。子图像I21(x，y)包含第一堵墙的后表面、第二堵墙的前表面以及它们之前的部分，子图像I22(x，y) 包含I2(x，y)剩下的部分。其中，I21(x，y)中墙体位置正确，I22(x，y)中墙体位置偏移。

6）对第二堵墙造成的时延进行补偿，这个补偿时延为

将I22(x，y)向天线方向平移(d1-D1)+(d2-D2)，得到补偿后的图像I3(x，y)。

7）对I3(x，y)应用Radon变换，以第三堵墙中心为分界线分割I3(x，y)，并对靠后的子图像进行补偿，得到I31(x，y)和补偿后的图像I4(x，y)。

8）将所有包含正确墙体位置的子图像I11(x，y)，I21(x，y)，I31(x，y)，I4(x，y)分别进行归一化并拼接，得到完整的图像I'(x，y)。

5 多视角融合

穿墙雷达探测的场景大多是封闭的室内环境，电磁波除了产生墙体前后表面回波外，还会产生墙体间多径回波、墙体内多径回波等。多径回波会导致成像中出现鬼影，严重影响对建筑布局的判断。

Q.Tan等研究发现多径鬼影对雷达探测位置具有方向依赖性（Aspect Dependence），即多径鬼影位置会随着雷达观测位置的变化而改变［9］。利用这一特性，从多个视角对建筑进行观测，通过对多个观测结果取最小值来抑制多径鬼影。

基于多视角融合的多径鬼影抑制算法原理如下：

（1）多视角观测下墙体成像位置不变

如图3 所示，从视角1 和视角2 对该墙体的左墙面进行成像。由于视角1 发射的电磁波仅在空气中传播，经过一次反射，接收到信号，因此墙体左表面成像位置就是真实位置。从视角2 观测，由于电磁波要穿透墙体，电磁波的传播速度发生改变，成像位置会变为图中左侧虚线位置。但通过第4 节墙体补偿后，墙体位置被矫正到图中右侧虚线位置，即真实位置。因此，对同一个墙面，从视角1 与视角2 进行观测得到的墙体成像结果一致且均为真实位置。此为多视角观测下目标墙体成像位置不变。

（2）多视角观测下墙体内多径鬼影成像位置不同

从视角1 与视角2 对墙体进行观测得到的墙体内多径鬼影如图4 所示。视角1 发射的电磁波在墙体内进行了三次反射，同时电磁波在墙体内的传播速度也会发生改变，因此该路径下多径回波产生的多径鬼影位置如图中右侧虚线所示。同理，由视角2 观测得到的墙体内多径鬼影如图中左侧虚线所示。综上可以发现，对同一个墙体由视角1 与视角2观测得到的墙体内多径鬼影位于墙体两侧。此为多视角观测下墙体内多径鬼影成像位置不同。

（3）多视角观测下墙体间多径鬼影成像位置不同

从视角1 与视角2 对墙体进行观测得到的墙体间多径鬼影如图5 所示。视角1 发射的电磁波在墙体间进行了三次反射，该路径下多径回波产生的多径鬼影位置如图中右侧虚线所示。同理，由视角2观测得到的墙体内多径鬼影如图中左侧虚线所示。因此，对同一组墙体，由视角1 与视角2 观测得到的墙体间多径鬼影位于该组墙体的两侧。此为多视角观测下墙体间多径鬼影成像位置不同。

（4）多径鬼影的方位依赖性

由前三点可以得出结论，在应用第4 节墙体补偿的前提下，多视角观测下墙体的成像位置不变，墙体内多径鬼影、墙体间多径鬼影成像位置不同，即多径鬼影的方位依赖性。基于这一特性，可以从多个视角对建筑进行观测，对多个观测结果取最小值来抑制多径鬼影，同时墙体信息得到保留。多视角观测示意图如图6所示。

假设有四个观测视角，每个观测视角得到的墙体补偿完毕后的图像为(xk，yk)(k=1，2，3，4)。由于每一个视角的成像结果都有自己对应的坐标系，因此在图像融合时需要统一各个视角的坐标系。假设全局坐标系为图6 所示的XOY坐标系，利用下式将局部成像结果转换为全局成像结果。

其中，Ik(x，y)为转换后的全局坐标系成像结果，rot为转换操作符，ψk=[90° 0° 270° 180°]为转换角度。

将多个视角的观测信息进行融合，得到

对每个像素点取所有图像对应位置的最小值，融合后成像结果仅在墙壁处有较高幅值。因为同一个多径回波不是所有观测视角都存在，所以取最小值后多径回波得到抑制，而且该方法不需要场景的先验信息即可有效的抑制多径鬼影。最终的I(x，y)为矫正墙体位置并抑制多径鬼影的建筑布局图像。

以三层墙体为例，墙体补偿与多视角融合方法流程图如图7所示。

6 仿真验证

下面通过仿真验证本文所提墙体补偿与多视角融合方法成像效果。使用gprMax［10-11］对一个房间模型进行仿真，房间结构布局如图8所示。建筑的墙体为0.24 m厚的砖砌成，相对介电常数为3.8。此外，模型中具有5 cm 厚的混凝土地板与天花板，相对介电常数为6.8。整个房间的大小为3 m×3 m×1 m。

发射信号为gaussiandot波形，中心频率1.0 GHz，通过偶极子天线发射。依次从四个视角对建筑物进行观测，每次观测时天线距离墙壁0.4 m，天线步进0.1 m，步进次数27 次，发射天线与接收天线距离0.2 m。

墙体补偿的仿真结果如图9 所示。图9 中，图9（a）为沿观测视角1 去掉直达波后的BP 成像。在此图中，第一堵墙前表面位置正确，其他墙体位置偏移。图9（b）为对图9（a）旋转90°后进行Radon变换后的结果，运用CFAR 检测得到第一堵墙的位置与厚度信息。图9（b）中，θ=90°列中的峰值对应图9（a）中的墙面。图9（c）和图9（d）分别为补偿第一堵墙与第二堵墙后的结果。在图9（c）中，第一堵墙后表面与第二堵墙前表面位置正确；在图9（d）中，第二堵墙后表面与第三堵墙前表面位置正确，其余墙体位置偏移。图9（e）为补偿第三堵墙后的图像，在此图中，第三堵墙后表面得到了矫正。图9（f）为所有墙体位置正确的子图像归一化后拼接而成的完整图像，在此图中，所有墙体的位置均正确。

多视角融合的仿真结果如图10 所示。图10中，图10（a）为对观测视角1 进行墙体补偿后的图像。图10（b）为对观测视角3 进行墙体补偿后的图像。图10（c）为图10（b）与图10（a）取最小值后的结果，即视角1 与视角3 融合后的结果。图10（d）为对图10（c）进行灰度线性变换增加对比度后的图像，在此图中，能够明显辨认三堵墙的位置，并且多径鬼影相比图10（a）、（b）显著减少。对观测视角2、4 进行相同的操作，得到图11 所示结果。

图12 给出了不同建筑布局成像方法的比较，图12（a）为对四个视角观测回波BP 成像后直接相加的结果，图12（b）为文献［5］所提方法得到的结果，图12（c）为文献［7］所提方法得到的结果，图12（d）为本文算法得到的结果。可以看出，与算数融合方法和文献［5］所提的M-N-K 方法相比，本文方法可以还原墙体的真实厚度与位置。文献［7］的方法虽然能还原墙体位置与厚度，但是增强了多径鬼影等杂波的强度，本文算法与之相比抑制了多径鬼影，改善了成像效果。

7 结论

本文针对穿墙雷达建筑布局成像中的墙体展宽、位置偏移和多径鬼影问题，提出了一种基于墙体补偿与多视角融合的穿墙雷达建筑布局成像方法。首先，建立了建筑布局成像的回波模型，推导了BP 成像算法；其次，针对墙体位置偏移问题，采用墙体位置检测并补偿聚焦时延的方法矫正墙体位置；之后，针对多径鬼影问题，采用多视角观测融合的方法抑制多径鬼影。仿真结果表明，与传统的建筑布局成像方法相比，本文所提方法能解决传统方法中的墙体展宽和位置偏移问题，并且能在图像融合的同时抑制多径鬼影，得到质量更高的建筑布局全景图。