Building Layout Imaging Algorithm of Wide-Beam Microwave Computerized Tomography

在微波层析透视探测成像领域，为了降低发射与接收天线的同步移动要求，提升系统探测效率，本文提出了宽波束微波CT探测模式，如图1所示。图中仅给出一组采样路径示例，可根据需求设计多组采样路径或不同形式的采样路径。其中，系统的发射天线数目为1，接收端采用包含若干天线的直线阵列或平面阵列接收信号。发射天线在有限个位置处发射扇形或锥形宽波束电磁波对未知场景进行探测，接收阵列沿若干规划路径移动，发射天线与接收阵列的中心天线间电磁传播直达路径与天线的移动路径相互垂直。在每组采样路径上均匀地采样信号功率信息，通过反演成像方法成像出未知场景内的建筑布局结构。

考虑二维场景的微波层析成像问题，通常情况下电磁波传播至位置${\boldsymbol{r}}$处的电场^[18]可表示为：

上式满足方程：

其近似解可表示为：

式中，${E_{{\rm{inc}}}}\left( {\boldsymbol{r}} \right)$表示位置${\boldsymbol{r}}$处的入射电场；相位$\phi \left( {\boldsymbol{r}} \right)$可表示为：

式中，$g\left( {{\boldsymbol{r}},{\boldsymbol{r}}'} \right)$表示自由空间中的格林函数；$I\left( {{\boldsymbol{r}}'} \right)$为位置${\boldsymbol{r}}'$处介质的电磁衰减率。式(3)包含了信号的幅度与相位信息，为获取信号的功率信息，对该式两边分别取共轭可得：

综合式(4)和式(5)，可得：

式中，${\text{Imag}}\left( \cdot \right)$表示对结果取虚部。对式(6)取对数得：

式中，$ {P_{\rm{RX}}}\left( {\boldsymbol{r}} \right)\;\left( {\rm{dBm}} \right) = 10{\rm{lg}}\left( {{{{{\left| {E\left( {\boldsymbol{r}} \right)} \right|}^2}} \mathord{\left/ {\vphantom {{{{\left| {E\left( {\boldsymbol{r}} \right)} \right|}^2}} {\left( {120{\text{π}} \times {{10}^{ - 3}}} \right)}}} \right. } {\left( {120{\text{π}} \times {{10}^{ - 3}}} \right)}}} \right) $为位置${\boldsymbol{r}}$处以dBm为单位的接收功率^[19]；$ {P_{{\rm{inc}}}}\left( {\boldsymbol{r}} \right) $为位置${\boldsymbol{r}}$处无物体时的入射功率。

为求解任意位置处的衰减特性参数$ I\left( {{\boldsymbol{r}}'} \right) $，将探测场景离散成$N$个等体积的单元格，用单元格中心位置向量${{\boldsymbol{r}}_n} (n \in \{ 1,2,\cdots,N\}) $表示，则式(4)可表示为：

在式(8)中，

式中，$\Delta V$表示每个单元格的体积，二维场景$D$中的高度维为1。令${{\boldsymbol{p}}_i}$、${{\boldsymbol{q}}_i}$分别表示第$i$次采样时，发射天线与接收天线的二维位置向量。

假设共采样$M$组数据，令：

第$i$个元素为：

式中，${E_{{\rm{inc}},{{\boldsymbol{p}}_i}}}\left( {{{\boldsymbol{r}}_n}} \right)$表示发射天线在$ {{\boldsymbol{p}}_i} $时，单元格${{\boldsymbol{r}}_n}$处的入射电场。综合式(10)和式(11)，令：

上式中${\boldsymbol{\varPhi }}$的任一元素可表示为：

式中，$g\left( {{{\boldsymbol{q}}_i},{{\boldsymbol{r}}_n}} \right)$表示自由空间的格林函数。位置$ {{\boldsymbol{p}}_i} $处的天线在${{\boldsymbol{r}}_n}$处产生的入射电场${E_{{\rm{inc}},{{\boldsymbol{p}}_i}}}\left( {{{\boldsymbol{r}}_n}} \right)$取决于天线的电磁波束覆盖范围，基于电磁传播直达路径的角度${\theta ^ * }$和天线的波束宽度${\theta _a}$，可得电磁波束的覆盖角度范围为$\left[ {{\theta ^ * } - {{{\theta _a}} \mathord{\left/ {\vphantom {{{\theta _a}} 2}} \right. } 2},{\theta ^ * } + {{{\theta _a}} \mathord{\left/ {\vphantom {{{\theta _a}} 2}} \right. } 2}} \right]$。由于倾斜角为90°的直线不存在斜率，因此根据电磁传播直达路径的角度${\theta ^ * }$是否为90°，用不同的限定范围来表示电磁波束的覆盖范围。

首先，当发射天线与接收阵列移动路径的角度$\theta \ne 0^\circ $，即发射天线与接收阵列的中心天线间电磁传播直达路径角度${\theta ^ * } \ne 90^\circ $时，电磁波的覆盖范围可直接由天线波束边界所在直线的斜率确定。在电磁波束方向图约束下，若${{\boldsymbol{r}}_n}$同时满足：

式中，${\theta _a}$表示发射与接收天线的波束宽度；$ K\left( {{{\boldsymbol{p}}_i},{{\boldsymbol{q}}_i};{{\boldsymbol{r}}_n}} \right) $表示发射天线${{\boldsymbol{p}}_i}$与单元格${{\boldsymbol{r}}_n}$连线所在直线的斜率，以及接收天线${{\boldsymbol{q}}_i}$与单元格${{\boldsymbol{r}}_n}$连线所在直线的斜率，则${E_{{\rm{inc}},{{\boldsymbol{p}}_i}}}\left( {{{\boldsymbol{r}}_n}} \right)$可表示为：

否则，${E_{{\rm{inc}},{{\boldsymbol{p}}_i}}}\left( {{{\boldsymbol{r}}_n}} \right) = 0$。

其次，当发射天线与接收阵列移动路径的角度$\theta = 0^\circ $，即电磁传播直达路径角度${\theta ^ * } = 90^\circ $时，由于倾斜角为$90^\circ $的直线不存在斜率，故将电磁波的覆盖范围分为3个区域表示，如式(16)所示。若${{\boldsymbol{r}}_n}$满足下式任一条件：

式中，$ {{\boldsymbol{r}}_n}\left( x \right) $、$ {{\boldsymbol{p}}_i}\left( x \right) $分别表示单元格${{\boldsymbol{r}}_n}$和发射天线${{\boldsymbol{p}}_i}$的$ x $坐标，则${E_{{\rm{inc}},{{\boldsymbol{p}}_i}}}\left( {{{\boldsymbol{r}}_n}} \right)$可表示为式(15)的形式，否则，${E_{{\rm{inc}},{{\boldsymbol{p}}_i}}}\left( {{{\boldsymbol{r}}_n}} \right) = 0$。

根据式(7)，将每一次的采样值代入式中，可得：

式中，向量${{\boldsymbol{P}}_{{\rm{RX}}}}\left( {{\rm{dBm}}} \right) = \left[ {P_{{\rm{RX}},{{\boldsymbol{p}}_1}}}\left( {{{\boldsymbol{q}}_1}} \right),{P_{{\rm{RX}},{{\boldsymbol{p}}_2}}}\left( {{{\boldsymbol{q}}_2}} \right),\cdots, {P_{{\rm{RX}},{{\boldsymbol{p}}_M}}}\left( {{{\boldsymbol{q}}_M}} \right) \right]^{\text{T}}$与向量${{\boldsymbol{P}}_{{\rm{inc}}}}\left( {{\rm{dBm}}} \right) = \left[ {P_{{\rm{inc}},{{\boldsymbol{p}}_1}}}\left( {{{\boldsymbol{q}}_1}} \right) , {P_{{\rm{inc}},{{\boldsymbol{p}}_2}}}\left( {{{\boldsymbol{q}}_2}} \right) , \cdots, {P_{{\rm{inc}},{{\boldsymbol{p}}_M}}}\left( {{{\boldsymbol{q}}_M}} \right) \right]^{\text{T}}$分别表示接收天线${{\boldsymbol{q}}_i}$处的接收信号功率与空场景时的接收信号功率。

结合式(12)与式(17)，可得：

式中，$ {{\boldsymbol{\varPhi }}_{\rm{R}}} $、$ {{\boldsymbol{\varPhi }}_{\rm{I}}} $、$ {{\boldsymbol{I}}_{\rm{R}}} $、$ {{\boldsymbol{I}}_{\rm{I}}} $分别为$ {\boldsymbol{\varPhi }} $与$ {\boldsymbol{I}} $的实部和虚部。由${{\boldsymbol{\varPhi }}_{\rm{R}}}{{\boldsymbol{I}}_{\rm{R}}} \geqslant {{\boldsymbol{\varPhi }}_{\rm{I}}}{{\boldsymbol{I}}_{\rm{I}}}$，则最终模型简化为：

式中，$ \Delta {\boldsymbol{P}}{\text{ = }}{\left[ {\Delta {{\boldsymbol{P}}_1},\Delta {{\boldsymbol{P}}_2},\cdots,\Delta {{\boldsymbol{P}}_M}} \right]^{\text{T}}} $为物体遮挡导致的衰减功率向量，$ M $为总采样点数；$ {{\boldsymbol{\varPhi }}_{\rm{R}}} $为$ M \times N $维系统矩阵，由式(13)求解可得；$ {{\boldsymbol{I}}_{\rm{R}}}{\text{ = }}{\left[ {{{\boldsymbol{I}}_1},{{\boldsymbol{I}}_2},\cdots,{{\boldsymbol{I}}_N}} \right]^{\text{T}}} $为建筑场景中各离散单元格的衰减特性参数。为简化后续符号体系，仍采用$\Delta {\boldsymbol{P}} = {\boldsymbol{\varPhi I}}$表示，其中系统矩阵$ {\boldsymbol{\varPhi }} $与成像矩阵$ {\boldsymbol{I}} $中的元素均为取实部后的数值。

针对宽波束微波CT模式的建筑布局层析成像问题，基于电磁波束方向图约束的Rytov信号模型，提出了一种改进型统计成像算法。在宽波束微波CT模式下，基于Rytov信号模型的建筑布局极大似然成像方法，即求解成像矩阵$ {\boldsymbol{I}} = {\left[ {{{\boldsymbol{I}}_1},{{\boldsymbol{I}}_2}, \cdots ,{{\boldsymbol{I}}_N}} \right]^{\text{T}}} $，使其对数似然函数达到最大值^[20-21]。将统计成像的迭代求解过程分为两部分。

1)在给定第$k$次迭代后成像结果的暂时解$ {{\boldsymbol{I}}^{\left( k \right)}} $及采样数据值$ \Delta {\boldsymbol{P}} $时，完备数据集$ {\boldsymbol{Z}} $的对数似然函数的条件期望可表示为：

2)对式(20)中条件期望最大化，从而得到成像矩阵$ {\boldsymbol{I}} $的更新解：

为求解成像矩阵的极大似然函数，设$ {h_i} > 0 $，$ {g_i}\left( {\boldsymbol{{\rm I}}} \right) = {\left( {{\boldsymbol{\varPhi I}}} \right)_i} = \displaystyle\sum\limits_{n = 1}^N {{{\boldsymbol{\varPhi }}_{i,n}}{{\boldsymbol{I}}_n}} $，则：

对于上式，存在：

忽略式(23)的常数项，替代函数$ {\Psi _{{\rm{EM}}}}\left( {{\boldsymbol{I}},{{\boldsymbol{I}}^{\left( k \right)}}} \right) $可表示为：

根据式(21)对式(24)求最大化，令$ {{\boldsymbol{I}}_n} $的偏导数为0，可得：

即得统计成像方程：

式中，$ {{\boldsymbol{I}}_n}^{\left( {k + 1} \right)} , {n = 1,2,\cdots,N} $为更新后成像矩阵中的第$n$个元素；$ {{\boldsymbol{\varPhi }}_{i,n}} $表示第$i$次采样时电磁波对第$n$个像素单元的影响；$M$和$N$分别表示总采样点数与建筑场景的离散单元格数；$ \Delta {{\boldsymbol{P}}_i} $为第$i$次采样时的衰减功率值。

在宽波束微波CT模式下，首先将所有采样值之和除以采样路径数与离散单元格数的积；然后，依采样角度数目$S$将采样数据值分为不同的子集合$\left( {{{{T}}_1},{{{T}}_2}, \cdots ,{{{T}}_S}} \right)$，分别利用每个子集合的数据更新成像矩阵，并将结果作为下次子迭代的初值，完成所有子集合对成像矩阵的校正，即为一轮统计成像过程；最后，将结果作为下一轮迭代的初值，重复上述操作，直至满足迭代终止条件，得到最终的建筑布局层析成像结果。具体步骤如下。

1) 将所采样建筑场景的平均值设定为初始成像矩阵的元素：

式中，${{\boldsymbol{I}}_n}, {n = 1,2, \cdots ,N} $表示建筑场景成像矩阵中任一元素；$ {N_q} $、$ {N_g} $分别表示不同角度采样路径的数目及单组路径的采样点数；$ \Delta {{\boldsymbol{P}}_{q,g}} $表示第$q$组采样路径中第$g$次采样时的衰减功率值；$N$为场景离散化后的单元格数目。

2) 对于第$s$个子集合，计算所有采样数据对应的估计值：

式中，${M_s}$为第$s$个子集合中采样数据的数目，$i = 1, 2, \cdots ,{M_s} $。对成像矩阵中第$n$个元素${{\boldsymbol{I}}_n} $进行校正，其中，${n = 1,2, \cdots ,N} $，其计算公式为：

3) 将下一子集合的数据重复式(28)和式(29)的更新过程，直至完成$S$个子集合的采样数据对成像矩阵的校正，即实现了一轮统计成像处理。

4) 将成像结果作为下一轮更新的初值，重复式(28)和式(29)的更新流程，直至满足终止条件得到建筑布局层析成像结果。所提方法的统计成像过程可归纳为：

在宽波束微波CT模式的建筑布局层析成像中，受限于探测场地及采集时间的要求，通常采用稀疏采样路径，因此将同一采样路径下的数据值作为一个子集合。

在宽波束微波CT探测模式下，为验证所提信号模型与成像方法的有效性，利用Matlab数值仿真软件开展了仿真实验，具体仿真参数如表1所示。

对于信号频率与带宽的选择，不同于雷达成像体制，微波层析成像体制中成像结果的分辨率与频率和带宽无关^[13]。为后续采用体积小、成本低的商用WiFi设备验证相关成像理论与方法，设定发射信号采用2.4 GHz正弦波。在建筑结构未知时，采样路径角度在0°～180°内均匀分布^[14]，成像效果较好，综合考虑采样复杂度、时间和成像效果，仿真与实测结果中均采用4组采样路径。为比较不同形式宽波束成像，仿真中增加了12组采样路径。微波层析体制的成像分辨率与采样间隔有关，视成像分辨率需求设为2 cm。采样间隔与单路径采样点数的乘积为单次采样路径的长度，需至少大于探测场景的边长，采样模式示意图如图2所示。

对于宽波束微波CT模式，系统利用多天线阵列接收信号，因此阵列天线数目对成像质量有较大影响。按照发射天线是否同步移动，图2中将宽波束微波CT模式分为：1)发射天线与接收阵列同步移动；2)发射天线与接收阵列均固定。仅给出$\theta = 0^\circ $和$\theta = 90^\circ $两组采样路径作为示意，实线为发射天线移动路径，虚线为接收阵列移动路径。

仿真场景如图3所示，场景尺寸为6.4 m×6.4 m，该场景由3个房间组成，其中墙体厚度均为0.16 m，场景离散单元格尺寸为0.02 m。通过数值仿真分析接收天线数分别为1、5和200时的结果，其中接收天线数为1时，即窄波束微波CT模式时，仿真结果如图4所示。接收天线数为5和200时，分别对应图2中两种宽波束微波CT模式，仿真结果如图5～图6所示。

为衡量成像结果的质量，引入图像处理领域常用的评价指标−结构相似性参数(structural similarity, SSIM)^[22]，利用成像结果与原始图像的像素间相关性表征其成像复原度。SSIM范围为−1～1，两幅图像完全一致时，其值为1，反之则为−1。

由图4～图6的仿真结果及表2不同接收天线数目时的成像结果参数对比可知，宽波束微波CT模式下接收天线数目为5和200时，均可实现建筑场景的层析成像，且稀疏采样条件下采样路径数较多时成像质量较好。接收天线数目为5时，即发射天线与接收阵列同步移动形式，4组采样路径即可重建出场景主要结构，得到较好的成像效果，SSIM为0.782；接收天线数目为200时，即发射天线与接收阵列均为固定形式，4组采样路径时仅可成像出部分结构，外部轮廓扩展严重，SSIM为0.241，成像质量较差。

在采样路径数目为12组，接收天线数目为1时，成像结果的SSIM为0.937，发射天线移动次数为2400次；接收天线数目为5与200时，成像结果的SSIM分别为0.915和0.816，发射天线移动次数分别为480次和12次。相较于窄波束微波CT模式，宽波束微波CT模式在牺牲部分成像复原度的情况下，减少了发射天线的移动次数，可以有效地节省对建筑场景进行数据采样的时间。

综上所述，宽波束微波CT模式可有效节省采样时间，在采样路径较少时，发射天线与接收阵列同步移动的宽波束微波CT模式具有明显优势。

为进一步验证所提信号模型与方法的有效性，基于电磁波束方向图约束的Rytov信号模型在发射天线与接收阵列同步移动的宽波束微波CT模式下开展了实验验证。根据3.1节中的仿真研究，兼顾系统复杂度与采样数据耗时，采用4组采样路径探测L形结构墙体，对应宽波束微波CT模式中接收天线数目为5，采样示意图如图7所示。

采用标准红砖搭建了L形结构墙体实验场景，其三维尺寸为2 m×2 m×1 m，墙体的厚度为0.23 m，具体实验参数如表3所示。4组采样路径的角度分别为0°、45°、90°和135°，发射天线与接收天线的移动间隔分别为10 cm和2 cm，发射天线移动1个位置时，接收天线需移动5个位置。

在成像复原度参数SSIM的基础上，均方根误差(root mean squared error, RMSE)被用于评估成像结果与原始场景的误差情况，其值越小，表示成像结果与原始场景的误差越小。

由图8宽波束微波CT模式下L形结构墙体的实测成像结果可知，不同算法均可成像出L形结构墙体的整体轮廓，验证了宽波束微波CT模式的可行性。在基于电磁波束方向图约束的Rytov信号模型下，相比于ART算法，本文方法可以减少成像伪影问题，抑制建筑结构边缘的成像扩展现象。

由表4不同算法的成像结果参数对比可知，相比于ART算法结果的SSIM为0.713，RMSE为0.314，本文方法成像结果的SSIM为0.808，RMSE为0.229，即所提方法的成像结果参数较优，成像质量更高，验证了宽波束微波CT模式下本文方法的有效性。

在建筑布局结构透视微波层析成像领域，针对窄波束微波CT模式收发天线同步移动要求高、采样耗时久等问题，本文设计了宽波束微波CT模式，建立了基于电磁波束方向图约束的Rytov信号模型，提出了一种基于宽波束信号的微波层析统计成像方法。通过数值仿真实验分析了不同形式宽波束微波CT模式的成像效果，在采样路径较少时，发射天线与接收阵列同步移动的的宽波束微波CT模式具有明显优势，最后在该探测模式下，通过实测数据验证了所提信号模型与方法的有效性。