ISSN 
-
在城市建筑区域中，采用低频电磁波进行透墙探测，包括对典型建筑场景内目标的探测定位[1-2]、跟踪[3-5]以及场景的布局结构成像[6-7]等，被广泛应用于地震救援、城市巷战等领域。由于城市中建筑构型的复杂性，感兴趣的目标通常会被建筑结构或室内家具等遮蔽，导致电磁波无法有效地探测到目标，影响目标定位和识别效果等。若能够获取建筑场景的结构作为已知先验信息，可进一步精确地得到探测目标的多维度信息，更有利于后续作战与救援行动的开展。
针对建筑布局结构成像问题，传统的成像思路是采用合成孔径雷达(synthetic aperture radar, SAR)从多个探测视角扫描场景，利用与发射天线同侧的接收天线获取反射和散射信号等，以实现建筑布局结构透视成像及建筑物内隐蔽目标的探测与成像[8-12]。2010年，美国陆军研究实验室搭建了具备2个发射天线和16个接收天线的车载雷达成像系统[8-9]，并从建筑物两个互相垂直的探测视角开展合成孔径扫描，采用后向投影(back projection, BP)方法将各视角下多个通道的数据进行相干成像，最后将各视角的成像结果非相干融合得到建筑结构成像结果。2013年，文献[10]基于车载SAR成像系统接收几种不同类型墙体的两个视角下信号回波，利用BP方法获取各视角的图像，并融合得到建筑成像结果。2018年，文献[11]将稀疏重构技术应用于穿墙建筑布局成像领域，提出了融合加权相干因子(coherence factor, CF)的稀疏成像算法，降低了旁瓣对物体成像的影响，增大了探测区域中弱散射物体的成像概率，通过实测结果对所提方法的有效性进行了分析与验证。以上研究采用超宽带(ultra wideband, UWB)雷达以合成孔径的形式探测并成像建筑布局[12]，为得到较高的准确度和分辨率，需要在系统体积、重量和功耗等方面付出代价，在系统复杂度和研发成本等方面存在一定的限制。
在建筑布局成像领域，文献[13-15]提出了一种新的解决思路，采用一对发射与接收天线分置的雷达系统在建筑场景外部以类CT形式探测场景，通过接收信号的强度信息，并利用场景的像素单元差异性对未知建筑布局进行成像。2018年，文献[16]基于微波CT采样模式，提出了改进型障碍物映射成像方法，采用逆面积椭圆传播模型表征信号传播路径穿过障碍物时的能量衰减，将方向信息引入空间相关矩阵，利用非负约束的改进型正则化算法提高了成像精度。相比于UWB雷达，这种新体制雷达具有一些独特优势：透射信号的单程接收形式，避免了信号的多程传播损耗；仅利用窄带信号的强度信息进行成像[17]，降低了系统复杂度等。然而，现阶段采用的微波CT模式仍存在不足之处，其采样模式借鉴于医学CT扫描，因此对发射与接收天线的同步移动要求较高，且扫描建筑场景的耗时较久。
在微波层析透视探测成像中，针对收发天线同步移动要求高、采样耗时久等问题，本文设计了宽波束微波CT模式，提出了一种基于宽波束信号的微波层析统计成像方法。通过求解建筑布局成像矩阵的最大似然函数，导出了统计成像求解方程，然后利用角度依赖的不同子集合数据更新成像矩阵，得到最终建筑布局图像，并通过实测结果验证了所提方法的有效性。
-
在微波层析透视探测成像领域，为了降低发射与接收天线的同步移动要求，提升系统探测效率，本文提出了宽波束微波CT探测模式，如图1所示。图中仅给出一组采样路径示例，可根据需求设计多组采样路径或不同形式的采样路径。其中，系统的发射天线数目为1，接收端采用包含若干天线的直线阵列或平面阵列接收信号。发射天线在有限个位置处发射扇形或锥形宽波束电磁波对未知场景进行探测，接收阵列沿若干规划路径移动，发射天线与接收阵列的中心天线间电磁传播直达路径与天线的移动路径相互垂直。在每组采样路径上均匀地采样信号功率信息，通过反演成像方法成像出未知场景内的建筑布局结构。
图 1 宽波束微波CT探测模式
考虑二维场景的微波层析成像问题，通常情况下电磁波传播至位置
${\boldsymbol{r}}$ 处的电场[18]可表示为：$$ E\left( {\boldsymbol{r}} \right) = {{\text{e}}^{{\text{j}}\psi \left( {\boldsymbol{r}} \right)}} $$ (1) 上式满足方程：
$$ \left[ {{\nabla ^2} + {k^2}\left( {\boldsymbol{r}} \right)} \right]E\left( {\boldsymbol{r}} \right) = 0 $$ (2) 其近似解可表示为：
$$ E\left( {\boldsymbol{r}} \right) = {E_{{\rm{inc}}}}\left( {\boldsymbol{r}} \right){{\text{e}}^{{\text{j}}\phi \left( {\boldsymbol{r}} \right)}} $$ (3) 式中，
${E_{{\rm{inc}}}}\left( {\boldsymbol{r}} \right)$ 表示位置${\boldsymbol{r}}$ 处的入射电场；相位$\phi \left( {\boldsymbol{r}} \right)$ 可表示为：$$ \phi \left( {\boldsymbol{r}} \right) = \frac{{ - {\text{j}}}}{{{E_{{\rm{inc}}}}\left( {\boldsymbol{r}} \right)}}\iiint\limits_{{R^3}} {g\left( {{\boldsymbol{r}},{\boldsymbol{r}}'} \right)I\left( {{\boldsymbol{r}}'} \right){E_{{\rm{inc}}}}\left( {{\boldsymbol{r}}'} \right){\text{d}}{\boldsymbol{r}}'} $$ (4) 式中，
$g\left( {{\boldsymbol{r}},{\boldsymbol{r}}'} \right)$ 表示自由空间中的格林函数；$I\left( {{\boldsymbol{r}}'} \right)$ 为位置${\boldsymbol{r}}'$ 处介质的电磁衰减率。式(3)包含了信号的幅度与相位信息，为获取信号的功率信息，对该式两边分别取共轭可得：$$ {E^*}\left( {\boldsymbol{r}} \right) = E_{{\rm{inc}}}^*\left( {\boldsymbol{r}} \right){{\text{e}}^{{\text{j}}{\phi ^*}\left( {\boldsymbol{r}} \right)}} $$ (5) 综合式(4)和式(5)，可得：
$$ {\left| {E\left( {\boldsymbol{r}} \right)} \right|^2} = {\left| {{E_{{\rm{inc}}}}\left( {\boldsymbol{r}} \right)} \right|^2}{{\text{e}}^{ - 2{\text{Imag}}\left( {\phi \left( {\boldsymbol{r}} \right)} \right)}} $$ (6) 式中，
${\text{Imag}}\left( \cdot \right)$ 表示对结果取虚部。对式(6)取对数得：$$ {P_{{\rm{RX}}}}\left( {\boldsymbol{r}} \right)\;\left( {{\rm{dBm}}} \right) = {P_{{\rm{inc}}}}\left( {\boldsymbol{r}} \right) + 10\lg \left( {{{\text{e}}^{ - 2}}} \right){\rm{Imag}}\left( {\phi \left( {\boldsymbol{r}} \right)} \right) $$ (7) 式中，
$ {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)可表示为：$$ \begin{split} & \phi \left( {\boldsymbol{r}} \right) = \frac{{ - {\text{j}}}}{{{E_{{\rm{inc}}}}\left( {\boldsymbol{r}} \right)}}\sum\limits_{n = 1}^N {I\left( {{{\boldsymbol{r}}_n}} \right){E_{{\rm{inc}}}}\left( {{{\boldsymbol{r}}_n}} \right)} \iiint\limits_{{R^3}} {g\left( {{\boldsymbol{r}},{\boldsymbol{r}}'} \right){\text{d}}{\boldsymbol{r}}'} \approx \\ &\qquad \frac{{ - {\text{j}}}}{{{E_{{\rm{inc}}}}\left( {\boldsymbol{r}} \right)}}\sum\limits_{n = 1}^N {g\left( {{\boldsymbol{r}},{{\boldsymbol{r}}_n}} \right)I\left( {{{\boldsymbol{r}}_n}} \right){E_{{\rm{inc}}}}\left( {{{\boldsymbol{r}}_n}} \right)} \Delta V \end{split} $$ (8) 在式(8)中，
$$ \iiint\limits_{{R^3}} {g\left( {{\boldsymbol{r}},{\boldsymbol{r}}'} \right)}{\text{d}}{\boldsymbol{r}}' \cong g\left( {{\boldsymbol{r}},{{\boldsymbol{r}}_n}} \right)\Delta V $$ (9) 式中，
$\Delta V$ 表示每个单元格的体积，二维场景$D$ 中的高度维为1。令${{\boldsymbol{p}}_i}$ 、${{\boldsymbol{q}}_i}$ 分别表示第$i$ 次采样时，发射天线与接收天线的二维位置向量。假设共采样
$M$ 组数据，令：$$ {\boldsymbol{\varGamma }} = {\left[ {{\phi _{{p_1}}}\left( {{{\boldsymbol{q}}_1}} \right),{\phi _{{p_2}}}\left( {{{\boldsymbol{q}}_2}} \right),\cdots,{\phi _{{p_M}}}\left( {{{\boldsymbol{q}}_M}} \right)} \right]^{\text{T}}} $$ (10) 第
$i$ 个元素为：$$ {\phi _{{p_i}}}\left( {{{\boldsymbol{q}}_i}} \right) = \frac{{ - {\text{j}}}}{{{E_{{\rm{inc}},{{\boldsymbol{p}}_i}}}\left( {{{\boldsymbol{q}}_i}} \right)}}\sum\limits_{n = 1}^N {g\left( {{{\boldsymbol{q}}_i},{{\boldsymbol{r}}_n}} \right)I\left( {{{\boldsymbol{r}}_n}} \right){E_{{\rm{inc}},{{\boldsymbol{p}}_i}}}\left( {{{\boldsymbol{r}}_n}} \right)\Delta V} $$ (11) 式中，
${E_{{\rm{inc}},{{\boldsymbol{p}}_i}}}\left( {{{\boldsymbol{r}}_n}} \right)$ 表示发射天线在$ {{\boldsymbol{p}}_i} $ 时，单元格${{\boldsymbol{r}}_n}$ 处的入射电场。综合式(10)和式(11)，令：$$ {\boldsymbol{\varGamma }} = - {\text{j}}{\boldsymbol{\varPhi I}} $$ (12) 上式中
${\boldsymbol{\varPhi }}$ 的任一元素可表示为：$$ {{\boldsymbol{\varPhi }}_{i,n}} = \frac{{g\left( {{{\boldsymbol{q}}_i},{{\boldsymbol{r}}_n}} \right){E_{{\rm{inc}},{{\boldsymbol{p}}_i}}}\left( {{{\boldsymbol{r}}_n}} \right)\Delta V}}{{{E_{{\rm{inc}},{{\boldsymbol{p}}_i}}}\left( {{{\boldsymbol{q}}_i}} \right)}} $$ (13) 式中，
$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}$ 同时满足：$$ \left\{ \begin{gathered} {\rm{tan}}\left[ {{{\left( {{\theta ^ * } - \frac{{{\theta _a}}}{2}} \right){\text{π}}} \mathord{\left/ {\vphantom {{\left( {{\theta ^ * } - \frac{{{\theta _a}}}{2}} \right){\text{π}}} {180^\circ }}} \right. } {180^\circ }}} \right] \leqslant K\left( {{{\boldsymbol{p}}_i},{{\boldsymbol{q}}_i};{{\boldsymbol{r}}_n}} \right) \\ K\left( {{{\boldsymbol{p}}_i},{{\boldsymbol{q}}_i};{{\boldsymbol{r}}_n}} \right) \leqslant {\rm{tan}}\left[ {{{\left( {{\theta ^ * } + \frac{{{\theta _a}}}{2}} \right){\text{π}}} \mathord{\left/ {\vphantom {{\left( {{\theta ^ * } + \frac{{{\theta _a}}}{2}} \right){\text{π}}} {180^\circ }}} \right. } {180^\circ }}} \right] \\ \end{gathered} \right. $$ (14) 式中，
${\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) = g\left( {{{\boldsymbol{r}}_n},{{\boldsymbol{p}}_i}} \right) = \frac{{{{\rm{e}}^{{\text{j}}{k_0}\left| {{{\boldsymbol{r}}_n} - {{\boldsymbol{p}}_i}} \right|}}}}{{4{\text{π}}\left| {{{\boldsymbol{r}}_n} - {{\boldsymbol{p}}_i}} \right|}} $$ (15) 否则，
${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}$ 满足下式任一条件：$$ \left\{ \begin{gathered} {\rm{tan}}\left[ {{{\left( {{\theta ^ * } - \frac{{{\theta _a}}}{2}} \right){\text{π}}} \mathord{\left/ {\vphantom {{\left( {{\theta ^ * } - \frac{{{\theta _a}}}{2}} \right){\text{π}}} {180^\circ }}} \right. } {180^\circ }}} \right] \leqslant K\left( {{{\boldsymbol{p}}_i},{{\boldsymbol{q}}_i};{{\boldsymbol{r}}_n}} \right) \\ K\left( {{{\boldsymbol{p}}_i},{{\boldsymbol{q}}_i};{{\boldsymbol{r}}_n}} \right) \leqslant {\rm{tan}}\left[ {{{\left( {{\theta ^ * } + \frac{{{\theta _a}}}{2}} \right){\text{π}}} \mathord{\left/ {\vphantom {{\left( {{\theta ^ * } + \frac{{{\theta _a}}}{2}} \right){\text{π}}} {180^\circ }}} \right. } {180^\circ }}} \right] \\ {{\boldsymbol{r}}_n}\left( x \right) = {{\boldsymbol{p}}_i}\left( x \right) \\ \end{gathered} \right. $$ (16) 式中，
$ {{\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)，将每一次的采样值代入式中，可得：
$$ \Delta {\boldsymbol{P}} = {\text{Imag}}\left( {\boldsymbol{\varGamma }} \right) = \frac{{{{\boldsymbol{P}}_{{\rm{RX}}}}\left( {{\rm{dBm}}} \right) - {{\boldsymbol{P}}_{{\rm{inc}}}}\left( {{\rm{dBm}}} \right)}}{{10\lg \left( {{{\text{e}}^{ - 2}}} \right)}} $$ (17) 式中，向量
${{\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)，可得：
$$ \Delta {\boldsymbol{P}} = {\text{Real}}({\boldsymbol{\varPhi I}}) = {{\boldsymbol{\varPhi }}_{\rm{R}}}{{\boldsymbol{I}}_{\rm{R}}} + {{\boldsymbol{\varPhi }}_{\rm{I}}}{{\boldsymbol{I}}_{\rm{I}}} $$ (18) 式中，
$ {{\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}} \approx {{\boldsymbol{\varPhi }}_{\rm{R}}}{{\boldsymbol{I}}_{\rm{R}}} $$ (19) 式中，
$ \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}} $ 的对数似然函数的条件期望可表示为：$$ {\Psi _{{\rm{EM}}}}\left( {{\boldsymbol{I}},{{\boldsymbol{I}}^{\left( k \right)}}} \right) = E\left[ {{\rm{log}}\left( {{\boldsymbol{Z}}\left| {\boldsymbol{I}} \right.} \right)\left| {\Delta {\boldsymbol{P}},{{\boldsymbol{I}}^{\left( k \right)}}} \right.} \right] $$ (20) 2)对式(20)中条件期望最大化，从而得到成像矩阵
$ {\boldsymbol{I}} $ 的更新解：$$ {{\boldsymbol{I}}^{\left( {k + 1} \right)}}{\text{ = }}\arg \max {\Psi _{\begin{subarray}{l} {\rm{EM}} \\ I \geqslant 0 \end{subarray} }}\left( {{\boldsymbol{I}},{{\boldsymbol{I}}^{\left( k \right)}}} \right) $$ (21) 为求解成像矩阵的极大似然函数，设
$ {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}} $ ，则：$$ L\left( {\boldsymbol{I}} \right) = \sum\limits_{i = 1}^M {\Delta {{\boldsymbol{P}}_i}\left[ {{\rm{log}}\left( {\sum\limits_{n = 1}^N {{{\boldsymbol{\varPhi }}_{i,n}}{{\boldsymbol{I}}_n} + {h_i}} } \right) - \left( {\sum\limits_{n = 1}^N {{\boldsymbol{\varPhi }}{}_{i,n}{{\boldsymbol{I}}_n} + {h_i}} } \right)} \right]} $$ (22) 对于上式，存在：
$$ \begin{split} & L\left( {\boldsymbol{I}} \right) \geqslant \sum\limits_{i = 1}^M {\left[ {\Delta {{\boldsymbol{P}}_i}\left\{ {\sum\limits_{n = 1}^N {\frac{{{{\boldsymbol{\varPhi }}_{i,n}}{{\boldsymbol{I}}_n}^{\left( k \right)}}}{{{g_i}\left( {{{\boldsymbol{I}}^{\left( k \right)}}} \right)}} \log \left( {\frac{{{{\boldsymbol{I}}_n}}}{{{{\boldsymbol{I}}_n}^{\left( k \right)}}} {g_i}\left( {{{\boldsymbol{I}}^{\left( k \right)}}} \right)} \right)} } \right.} \right.} + \\ &\qquad \left. {\left. {{\text{ }} \frac{{{h_i}}}{{{g_i}\left( {{{\boldsymbol{I}}^{\left( k \right)}}} \right)}}\log {g_i}\left( {{{\boldsymbol{I}}^{\left( k \right)}}} \right)} \right\} - \left( {\sum\limits_{n = 1}^N {{{\boldsymbol{\varPhi }}_{i,n}}{{\boldsymbol{I}}_n} + {h_i}} } \right)} \right] \end{split} $$ (23) 忽略式(23)的常数项，替代函数
$ {\Psi _{{\rm{EM}}}}\left( {{\boldsymbol{I}},{{\boldsymbol{I}}^{\left( k \right)}}} \right) $ 可表示为：$$ \begin{split} & {\Psi _{{\rm{EM}}}}\left( {{\boldsymbol{I}},{{\boldsymbol{I}}^{\left( k \right)}}} \right){\text{ = }}\sum\limits_{n = 1}^N {\left\{ {\left( {\sum\limits_{i = 1}^M {\frac{{{{\boldsymbol{\varPhi }}_{i,n}}{{\boldsymbol{I}}_n}^{\left( k \right)}}}{{{g_i}\left( {{{\boldsymbol{I}}^{\left( k \right)}}} \right)}}} \Delta {{\boldsymbol{P}}_i}} \right)\log {{\boldsymbol{I}}_n}} \right.} - \\ & \qquad\qquad\qquad \left. { \sum\limits_{i = 1}^M {{{\boldsymbol{\varPhi }}_{i,n}}} {{\boldsymbol{I}}_n}} \right\} \end{split} $$ (24) 根据式(21)对式(24)求最大化，令
$ {{\boldsymbol{I}}_n} $ 的偏导数为0，可得：$$ - \sum\limits_{i = 1}^M {{{\boldsymbol{\varPhi }}_{i,n}} + \frac{{{{\boldsymbol{I}}_n}^{\left( k \right)}}}{{{{\boldsymbol{I}}_n}}}} \sum\limits_{i = 1}^M {\frac{{\Delta {{\boldsymbol{P}}_i}{{\boldsymbol{\varPhi }}_{i,n}}}}{{\displaystyle\sum\limits_{n = 1}^N {{{\boldsymbol{\varPhi }}_{i,n}}{{\boldsymbol{I}}_n}^{\left( k \right)}} }}} = 0 $$ (25) 即得统计成像方程：
$$ {{\boldsymbol{I}}_n}^{\left( {k + 1} \right)} = \frac{{{{\boldsymbol{I}}_n}^{\left( k \right)}}}{{\displaystyle\sum\limits_{i = 1}^M {{{\boldsymbol{\varPhi }}_{i,n}}} }}\sum\limits_{i = 1}^M {\frac{{{{\boldsymbol{\varPhi }}_{i,n}}\Delta {{\boldsymbol{P}}_i}}}{{\displaystyle\sum\limits_{n = 1}^N {{{\boldsymbol{\varPhi }}_{i,n}}{{\boldsymbol{I}}_n}^{\left( k \right)}} }}} $$ (26) 式中，
$ {{\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} = \frac{{\displaystyle\sum\limits_{q = 1}^{{N_q}} {\displaystyle\sum\limits_{g = 1}^{{N_g}} {\Delta {{\boldsymbol{P}}_{q,g}}} } }}{{{N_q}N}} $$ (27) 式中，
${{\boldsymbol{I}}_n}, {n = 1,2, \cdots ,N} $ 表示建筑场景成像矩阵中任一元素；$ {N_q} $ 、$ {N_g} $ 分别表示不同角度采样路径的数目及单组路径的采样点数；$ \Delta {{\boldsymbol{P}}_{q,g}} $ 表示第$q$ 组采样路径中第$g$ 次采样时的衰减功率值；$N$ 为场景离散化后的单元格数目。2) 对于第
$s$ 个子集合，计算所有采样数据对应的估计值：$$ \Delta {{\boldsymbol{P}}_i} = \sum\limits_{n = 1}^N {{{\boldsymbol{\varPhi }}_{i,n}}{{\boldsymbol{I}}_n}} $$ (28) 式中，
${M_s}$ 为第$s$ 个子集合中采样数据的数目，$i = 1, 2, \cdots ,{M_s} $ 。对成像矩阵中第$n$ 个元素${{\boldsymbol{I}}_n} $ 进行校正，其中，${n = 1,2, \cdots ,N} $ ，其计算公式为：$$ {{\boldsymbol{I}}_n} = {{\boldsymbol{I}}_n}\frac{1}{{\displaystyle\sum\limits_{i = 1}^M {{{\boldsymbol{\varPhi }}_{i,n}}} }}\sum\limits_{i = 1}^M {\frac{{\Delta {{\boldsymbol{P}}_i}}}{{\displaystyle\sum\limits_{n = 1}^N {{{\boldsymbol{\varPhi }}_{i,n}}{{\boldsymbol{I}}_n}} }}{{\boldsymbol{\varPhi }}_{i,n}}} $$ (29) 3) 将下一子集合的数据重复式(28)和式(29)的更新过程，直至完成
$S$ 个子集合的采样数据对成像矩阵的校正，即实现了一轮统计成像处理。4) 将成像结果作为下一轮更新的初值，重复式(28)和式(29)的更新流程，直至满足终止条件得到建筑布局层析成像结果。所提方法的统计成像过程可归纳为：
$$ {{\boldsymbol{I}}_n}^{\left( {k + 1,s + 1} \right)} = \frac{{{{\boldsymbol{I}}_n}^{\left( {s,k} \right)}}}{{\displaystyle\sum\limits_{i = 1}^{{M_s}} {{{\boldsymbol{\varPhi }}_{i,n}}} }}\sum\limits_{i = 1}^{{M_s}} {\frac{{\Delta {{\boldsymbol{P}}_i}{{\boldsymbol{\varPhi }}_{i,n}}}}{{\displaystyle\sum\limits_{n = 1}^N {{{\boldsymbol{\varPhi }}_{i,n}}{{\boldsymbol{I}}_n}^{\left( {s,k} \right)}} }}} $$ (30) 在宽波束微波CT模式的建筑布局层析成像中，受限于探测场地及采集时间的要求，通常采用稀疏采样路径，因此将同一采样路径下的数据值作为一个子集合。
-
在宽波束微波CT探测模式下，为验证所提信号模型与成像方法的有效性，利用Matlab数值仿真软件开展了仿真实验，具体仿真参数如表1所示。
表 1 仿真参数表
仿真参数 数值 发射信号/GHz 2.4 发射天线数目/个 1 接收天线数目/个 1，5，200 单组路径采样点数/点 200 采样路径数目/组 4，12 采样间隔/m 0.02 对于信号频率与带宽的选择，不同于雷达成像体制，微波层析成像体制中成像结果的分辨率与频率和带宽无关[13]。为后续采用体积小、成本低的商用WiFi设备验证相关成像理论与方法，设定发射信号采用2.4 GHz正弦波。在建筑结构未知时，采样路径角度在0°～180°内均匀分布[14]，成像效果较好，综合考虑采样复杂度、时间和成像效果，仿真与实测结果中均采用4组采样路径。为比较不同形式宽波束成像，仿真中增加了12组采样路径。微波层析体制的成像分辨率与采样间隔有关，视成像分辨率需求设为2 cm。采样间隔与单路径采样点数的乘积为单次采样路径的长度，需至少大于探测场景的边长，采样模式示意图如图2所示。
图 2 不同形式的宽波束微波CT模式示意图
对于宽波束微波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所示。
图 3 仿真场景
图 4 接收天线数目为1时的仿真结果
图 5 接收天线数目为5时的仿真结果
图 6 接收天线数目为200时的仿真结果
为衡量成像结果的质量，引入图像处理领域常用的评价指标−结构相似性参数(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，成像质量较差。
表 2 不同接收天线数目时结果参数对比
接收天线数目/个 采样路径数目/组 SSIM 发射天线移动次数/次 1 4 0.804 800 12 0.937 2400 5 4 0.782 160 12 0.915 480 200 4 0.241 4 12 0.816 12 在采样路径数目为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所示。
图 7 宽波束微波CT模式采样示意图
采用标准红砖搭建了L形结构墙体实验场景，其三维尺寸为2 m×2 m×1 m，墙体的厚度为0.23 m，具体实验参数如表3所示。4组采样路径的角度分别为0°、45°、90°和135°，发射天线与接收天线的移动间隔分别为10 cm和2 cm，发射天线移动1个位置时，接收天线需移动5个位置。
表 3 实验参数表
实验参数 数值 发射信号/GHz 2.4 发射天线数目/个 1 等效接收天线数目/个 5 单组路径采样点数/点 121 采样路径数目/组 4 采样间隔/m 0.02 在成像复原度参数SSIM的基础上，均方根误差(root mean squared error, RMSE)被用于评估成像结果与原始场景的误差情况，其值越小，表示成像结果与原始场景的误差越小。
由图8宽波束微波CT模式下L形结构墙体的实测成像结果可知，不同算法均可成像出L形结构墙体的整体轮廓，验证了宽波束微波CT模式的可行性。在基于电磁波束方向图约束的Rytov信号模型下，相比于ART算法，本文方法可以减少成像伪影问题，抑制建筑结构边缘的成像扩展现象。
图 8 宽波束微波CT模式下L形结构墙体实测成像结果
由表4不同算法的成像结果参数对比可知，相比于ART算法结果的SSIM为0.713，RMSE为0.314，本文方法成像结果的SSIM为0.808，RMSE为0.229，即所提方法的成像结果参数较优，成像质量更高，验证了宽波束微波CT模式下本文方法的有效性。
表 4 宽波束微波CT模式下不同算法的成像结果对比
成像算法 SSIM RMSE ART算法 0.713 0.314 本文方法 0.808 0.229 -
在建筑布局结构透视微波层析成像领域，针对窄波束微波CT模式收发天线同步移动要求高、采样耗时久等问题，本文设计了宽波束微波CT模式，建立了基于电磁波束方向图约束的Rytov信号模型，提出了一种基于宽波束信号的微波层析统计成像方法。通过数值仿真实验分析了不同形式宽波束微波CT模式的成像效果，在采样路径较少时，发射天线与接收阵列同步移动的的宽波束微波CT模式具有明显优势，最后在该探测模式下，通过实测数据验证了所提信号模型与方法的有效性。
Building Layout Imaging Algorithm of Wide-Beam Microwave Computerized Tomography
-
摘要: 针对宽波束微波层析建筑布局成像问题，考虑实际应用中天线发射的电磁波具有一定波束宽度，建立了基于电磁波束方向图约束的Rytov信号模型，并提出了基于宽波束信号的微波层析统计成像方法。首先，通过构造建筑布局成像矩阵的最大似然函数，导出了统计成像求解方程；然后，通过稀疏角度依赖的各子集合数据更新成像矩阵，实现了宽波束微波层析建筑布局成像。该方法适用于稀疏采样角度下宽波束微波层析成像，且对噪声有较强抑制作用。通过仿真与实验数据验证了该信号模型与方法的有效性。Abstract: Aiming at the building layout imaging problem of wide-beam microwave computerized tomography (CT), considering that the electromagnetic waves emitted by antennas have a certain beamwidth in practical applications, the Rytov signal model is established based on the constraint of electromagnetic beam direction graph, and a statistical imaging method of microwave tomography based on the wide-beam signal is proposed. First, by constructing the maximum likelihood function of the imaging matrix of the building layout, the statistical imaging equation is derived. Then, the imaging matrix is updated by the sparse angle-dependent subset data, which realizes the microwave tomography of the building layout under wide-beam signal. The proposed method is suitable for wide-beam microwave CT imaging with sparse sampling angles, and has a strong suppression effect on noise. The validity of the proposed signal model and method are verified by simulation and experimental data.
-
表 2 不同接收天线数目时结果参数对比
接收天线数目/个 采样路径数目/组 SSIM 发射天线移动次数/次 1 4 0.804 800 12 0.937 2400 5 4 0.782 160 12 0.915 480 200 4 0.241 4 12 0.816 12 
下载: 导出CSV
-
[1] SU N Y, CHEN X L, GUAN J, et al. Maritime target detection based on radar graph data and graph convolutional network[J]. IEEE Geoscience and Remote Sensing Letters, 2022, 19: 1-5. [2] 崔国龙, 孔令讲, 杨建宇. 步进变频穿墙成像雷达中反投影算法研究[J]. 电子科技大学学报, 2008, 37(6): 864-867. doi: 10.3969/j.issn.1001-0548.2008.06.016 CUI G L, KONG L J, YANG J Y. Back-Projection algorithm to stepped-frequency through-the-wall radar imaging[J]. Journal of University of Electronic Science and Technology of China, 2008, 37(6): 864-867. doi: 10.3969/j.issn.1001-0548.2008.06.016 [3] KOCUR D, SVECOVA M, ROVNAKOVA J. Through-The-Wall localization of a moving target by two independent ultra wideband radar systems[J]. Sensors, 2013, 13(9): 11969-11997. doi: 10.3390/s130911969 [4] LIU Z Q, DENG Y J, MA F, et al. Target detection and tracking algorithm based on improved Mask RCNN and LMB[C]//International Conference on Control, Automation and Information Sciences. Xi'an: [s.n.], 2021: 1037-1041. [5] 孔令讲, 陈国浩, 崔国龙, 等. 基于均值漂移的穿墙雷达多目标跟踪[J]. 电子科技大学学报, 2019, 48(3): 321-325. doi: 10.3969/j.issn.1001-0548.2019.03.001 KONG L J, CHEN G H, CUI G L, et al. Mean-Shift based multiple targets tracking for TWR[J]. Journal of University of Electronic Science and Technology of China, 2019, 48(3): 321-325. doi: 10.3969/j.issn.1001-0548.2019.03.001 [6] CHEN J, ZHANG Y, LI H, et al. Ultrawideband tomographic imaging in multipath-rich environment[EB/OL]. (2021-07-26). https://ieeexplore.ieee.org/document/9494710. [7] 郭睿, 臧博, 景国彬. 基于机载InSAR干涉相位的高层建筑重构[J]. 电子科技大学学报, 2019, 48(5): 717-721. doi: 10.3969/j.issn.1001-0548.2019.05.009 GUO R, ZANG B, JING G B. High-Rise building reconstruction using airborne InSAR interferometric phase[J]. Journal of University of Electronic Science and Technology of China, 2019, 48(5): 717-721. doi: 10.3969/j.issn.1001-0548.2019.05.009 [8] LE C, DOGARU T, NGUYEN L, et al. Ultrawideband (UWB) radar imaging of building interior: Measurements and predictions[J]. IEEE Transactions on Geoscience and Remote Sensing, 2009, 47(5): 1409-1420. doi: 10.1109/TGRS.2009.2016653 [9] LE C, NGUYEN L, DOGARU T. Radar imaging of a large building based on near-field Xpatch model[C]//IEEE International Symposium on Antennas and Propagation and CNC-USNC/URSI Radio Science Meeting. Toronto: IEEE, 2010: 1-4. [10] SEVIGNY P, DIFILIPPO D J. A multi-look fusion approach to through-wall radar imaging[C]//2013 IEEE Radar Conference. Ottawa: IEEE, 2013: 1-6. [11] 金添, 宋勇平. 超宽带雷达建筑物结构稀疏成像[J]. 雷达学报, 2018, 7(3): 275-284. JIN T, SONG Y P. Sparse imaging of building structures by UWB radar[J]. Journal of Radars, 2018, 7(3): 275-284. [12] 金添, 宋勇平, 崔国龙, 等. 低频电磁波建筑物内部结构透视技术研究进展[J]. 雷达学报, 2021, 10(3): 342-359. doi: 10.12000/JR20119 JIN T, SONG Y P, CUI G L, et al. Research progress of low frequency electromagnetic wave perspective technology for internal structure of building[J]. Journal of Radars, 2021, 10(3): 342-359. doi: 10.12000/JR20119 [13] MOSTOFI Y. Compressive cooperative sensing and mapping in mobile networks[J]. IEEE Transactions on Mobile Computing, 2011, 10(12): 1769-1784. doi: 10.1109/TMC.2011.31 [14] GONZALEZ-RUIZ A, MOSTOFI Y. Cooperative robotic structure mapping using wireless measurements-a comparison of random and coordinated sampling patterns[J]. IEEE Sensors Journal, 2013, 13(7): 2571-2580. doi: 10.1109/JSEN.2013.2256122 [15] DEPATLA S, KARANAM C R, MOSTOFI Y. Robotic through-wall imaging: Radio-Frequency imaging possibilities with unmanned vehicles[J]. IEEE Antennas and Propagation Magazine, 2017, 59(5): 47-60. doi: 10.1109/MAP.2017.2731302 [16] XU S, LIU H, GAO F, et al. Experimental verification: Enabling obstacle mapping based on radio tomographic imaging[C]//IEEE 88th Vehicular Technology Conference. Chicago: IEEE, 2018: 1-5. [17] ZHANG Y, CHEN J, GUO S, et al. Building layout tomographic reconstruction via commercial WiFi signals[J]. IEEE Internet of Things Journal, 2021, 8(20): 15500-15511. doi: 10.1109/JIOT.2021.3073912 [18] CHEW W. Waves and fields in inhomogeneous media[M]. New York: IEEE Press, 1995: 447-488. [19] DEPATLA S, BUCKLAND L, MOSTOFI Y. X-Ray vision with only WiFi power measurements using Rytov wave models[J]. IEEE Transactions on Vehicular Technology, 2015, 64(4): 1376-1387. doi: 10.1109/TVT.2015.2397446 [20] ZHENG X, QU G, ZHOU J. Accelerated strategy for the MLEM algorithm[J]. Journal of X-Ray Science and Technology, 2021, 29(1): 135-149. doi: 10.3233/XST-200749 [21] MU X F, SHI R, LUO G, et al. High spatial resolution tomographic gamma scanning reconstruction with improved MLEM iterative algorithm based on split Bregman total variation regularization[J]. IEEE Transactions on Nuclear Science, 2021, 68(12): 2762-2770. doi: 10.1109/TNS.2021.3125001 [22] 佟雨兵, 张其善, 祁云平. 基于PSNR与SSIM联合的图像质量评价模型[J]. 中国图象图形学报, 2006(12): 1758-1763. doi: 10.11834/jig.2006012307 TONG Y B, ZHANG Q S, QI Y P. Image quality evaluation model based on the combination of PSNR and SSIM[J]. Journal of Image and Graphics, 2006(12): 1758-1763. doi: 10.11834/jig.2006012307 -
点击查看大图
图(8) / 表(4)
计量
- 文章访问数: 3748
- HTML全文浏览量: 1123
- PDF下载量: 80
- 被引次数: 0
