Building Layout Tomography Method Based on Joint Multidomain Direct Wave Estimation
-
摘要: 在进入陌生建筑物之前掌握其内部结构信息,能够为反恐作战、灾害救援、监视管控等多种应用提供支持,具有重要的现实意义和研究价值。为实现建筑布局结构信息获取,该文开展了基于多域联合直达波估计的建筑布局层析成像方法研究。首先,建立了线性近似模型,实现了直达波信号传播时延与未知建筑布局图像之间的映射关系;在此模型基础上,分析了在层析成像模式下直达波信号与多径信号在快时间域、慢时间域与多普勒域中的分布特性,提出了一种基于多域联合的直达波估计算法,实现了多径干扰抑制与直达波信号精确估计;此外,提出了一种总变分约束的投影矩阵自适应修正代数重建算法,提升了有限数据下的建筑布局反演质量;最后,电磁仿真与实测实验结果证明了所提出的建筑布局层析成像方法的有效性,其重建结果的结构相似性指标分别可达到91.2%和81.7%,显著优于现有建筑布局层析成像方法。Abstract: Obtaining internal layout information before entering unfamiliar buildings is crucial for various applications, such as counter-terrorism operations, disaster relief, and surveillance, highlighting its great practical significance and research value. To enable the acquisition of the building layout information, this paper presents a building layout tomography method based on joint multidomain direct wave estimation. First, a linear approximation model is established to map the relationship between the propagation delay of direct wave signals and the layout of the unknown building. Using this model, the distribution characteristics of direct wave and multipath signals in the fast-time, slow-time, and Doppler domains are analyzed in the tomographic imaging mode. A joint multidomain direct wave estimation algorithm is then proposed to achieve the suppression of multipath interference and precise estimation of direct wave signals. Additionally, a projection matrix adaptive correction algebraic reconstruction algorithm with total variation constraints is proposed, which enhances building layout inversion quality under limited data scenarios. Finally, electromagnetic simulation and experimental results demonstrate the effectiveness of the proposed building layout tomography method, with structural similarity indices of 91.2% and 81.7% for the reconstructed results, significantly outperforming existing building layout tomography methods.
-
1. 引言
城市环境下的反恐作战、灾害救援等行动,作战、救援人员往往需要深入陌生建筑物。缺乏建筑内部结构信息将严重影响人员安全和行动顺利。因此,利用微波信号透视墙体等障碍物,探测未知建筑盲区具有重要的社会意义和研究价值[1−8]。
针对建筑布局重建问题,传统的解决方案多采用超宽带(Ultra-Wideband, UWB)雷达从多个视角对建筑场景进行合成孔径成像,通过反射、散射信号反演建筑布局结构及建筑物内隐蔽目标。2008年,美国陆军研究实验室Nguyen等人[9]公布了具备两个发射天线和16个接收天线的车载穿墙雷达系统。次年,Le等人[10]利用上述穿墙雷达系统从建筑物两个互相垂直的探测视角开展合成孔径扫描,采用后向投影(Back Projection, BP)方法将各视角下多个通道的数据进行相干成像,最后将各视角的成像结果非相干融合得到建筑结构成像结果。为了提升建筑布局重建质量,金添等人[11]将稀疏重构技术应用于穿墙成像领域,提出了基于加权相干因子(Coherence Factor, CF)的稀疏成像算法,并建立多层墙体位置矫正模型,实现了墙体的高精度定位与成像,并通过实测数据验证了所提方法的有效性。由于电磁波在穿透墙体时伴随着速度变化、能量衰减以及多次反射,基于合成孔径成像模式的穿墙雷达成像方法存在内墙偏移、内墙被淹没和虚假墙体等现象[12],依赖墙体参数信息估计和墙体位置检测等复杂流程进行校正。
不同于合成孔径成像技术,美国加州大学Mostofi等人[13−16]受医学计算机断层成像技术(Computer Tomography, CT)启发,提出了一种新的穿墙成像技术,后文将该类技术称为射频层析成像技术(Radio Tomographic Imaging, RTI)。RTI采用分布式系统,将发射节点和接收节点放置在建筑物的对侧以采样穿透建筑物的透射信号。随后,利用不同测量位置处接收信号强度信息(Receive Signal Strength Information, RSSI)的差异并结合压缩感知技术实现了未知建筑结构的反演重构。相比于传统UWB合成孔径雷达成像方案,RTI具有一些独特优势:透射信号的单程接收形式,减少了信号传播损耗;仅利用窄带信号的强度信息进行成像,降低了系统要求等[17]。然而,RSSI是对接收信号功率的时间平均,难以避免多径效应带来的干扰,致使成像质量恶化。为解决这一问题,文献[18]提出了利用宽带信号分离多径信号与直达波信号的思路,并根据最大幅度假设(Maximum Amplitude Estimation, MAE)对直达波信号的传播时延进行估计。尽管如此,MAE在多径丰富的环境中并不可靠,Chen等人[19]在此基础上提出了一种基于多径甄别的直达波估计算法,考虑了多径干扰主要由侧墙一次反射引起,通过BP成像中侧墙的散射点位置甄别多径信号。虽然该方法能够有效抑制多径影响进而提升成像性能,然而存在处理流程复杂、依赖多径散射点集中等问题。
为了实现建筑布局的精确重构,围绕多径干扰、稀疏成像问题,本文提出了一种基于多域联合直达波估计的建筑布局层析成像方法。首先,通过线性近似模型映射直达波时延与建筑布局的关系;然后,提出了一种基于多域联合的直达波估计算法(Direct wave estimation algorithm based on multi domain joint, MD-DE),该方法能够有效提升多径环境下直达波信号的估计精度;此外,针对病态反演问题,提出了一种总变分约束的投影矩阵自适应修正代数重建算法(Projection Matrix Adaptive Modification Algebraic Reconstruction Technique algorithm with Total Variation constrained, PMAM-ART-TV),实现了稀疏数据下的建筑布局高质量重建。最后,通过仿真结果和实验结果验证了所提方法的有效性。
2. 基于线性近似的层析成像模型
以二维场景探测问题为例,层析成像工作模式二维示意图如图1所示。发射天线和接收天线分别被分置在场景的两侧,发射天线在一侧发射雷达信号,接收端在另一侧接收穿过场景的雷达信号。随后,收发天线沿着预设的路径同步移动,以获取不同位置处经未知场景调制的透射信号。完成多位置的扫描后,采集的透射数据将用于反演扫描区域内的介质参数分布情况,进而实现建筑布局成像。
在RTI中,每次采样得到的信号可被认为是直达波信号与多径信号的叠加,令发射信号为s(t),第m个采样位置的接收信号r\left( {{{{\boldsymbol{X}}}_m},t} \right)可表示为
r\left( {{{{\boldsymbol{X}}}_m},t} \right) = {\sigma _m}s\left( {t - {\tau _m}} \right) + \sum\limits_{l = 1}^{{L_{{m}}}} {\sigma _m^ls\left( {t - \tau _m^l} \right)} + n\left( t \right) (1) 其中,{{{\boldsymbol{X}}}_m}表示第m个采样位置处发射天线和接收天线的坐标,{\sigma _m}和{\tau _m}表示直达波的衰减系数和传播时延,\sigma _m^l和\tau _m^l表示多径信号的衰减系数和传播时延,{L_m}为第m次采样时多路径的数量,n\left( t \right)为噪声。
根据电磁传播理论[18],电磁波在不同介质中的传播速度不同,因此根据电磁波传播时延可反推传播路径上介质的物理参数信息。具体地,电磁波传播时延与介质的介电常数以及传播路径的长度有关,表示为
\tau = \int_{{{\mathbb{L}}_{\overrightarrow {{\text{TX}}} \to \overrightarrow {{\text{RX}}} }}} {\frac{{\sqrt {\varepsilon \left( {r} \right)} }}{{\mathrm{c}}}} {\text{ d}}{r} (2) 其中, {{\mathbb{L}}_{\overrightarrow {{\text{TX}}} \to \overrightarrow {{\text{RX}}} }} 表示以收发天线为两端点的直达波路径,如图1所示; \varepsilon \left( {r} \right) 表示直达波路径上位置r处介质的相对介电常数,c表示光速。根据式(2)可知,电磁波传播时延与传播路径上不同位置处的介电常数有关,通过计算反演介电常数空间分布,即可表征未知场景的建筑布局信息。为了便于描述场景中建筑布局的分布情况,将整个二维场景离散为N = P \times Q个尺寸为\Delta d \times \Delta d的网格。此时,令每个网格的坐标为{{{\boldsymbol{r}}}_n},n \in \left\{ {1,2,\cdots,N} \right\},假设各个网格内的物理参数相同,则二维场景可表示为
{O} = {\left[ {\Delta d\sqrt {\varepsilon \left( {{{{{r}}}_1}} \right)} /{\text{c }}\Delta d\sqrt {\varepsilon \left( {{{{{r}}}_2}} \right)} /{{\mathrm{c}}}{\text{ }} \cdots {\text{ }}\Delta d\sqrt {\varepsilon \left( {{{{{r}}}_N}} \right)} /{\text{c}}} \right]^{\mathrm{T}}} (3) 通过对场景离散化,式(2)可以被重写为
\tau = \int_{{{\mathbb{L}}_{\overrightarrow {{\text{TX}}} \to \overrightarrow {{\text{RX}}} }}} {\frac{{\sqrt {\varepsilon \left( {{\boldsymbol{r}}} \right)} }}{{\text{c}}}} {\text{ d}}{r} \approx \sum\limits_{n \in {{\mathbb{L}}_{\overrightarrow {{\text{TX}}} \to \overrightarrow {{\text{RX}}} }}} {\frac{{\sqrt {\varepsilon \left( {{{{\boldsymbol{r}}}_n}} \right)} }}{{\text{c}}}} \Delta d (4) 当收发天线沿着预设路径以固定步长获取多个采样位置的数据后,通过对各个采样位置处的直达波时延进行估计,可以得到一个列向量{{\boldsymbol{P}}} = \left[ {{{\bar \tau }_1},{{\bar \tau }_2},\cdots,{{\bar \tau }_M}} \right],其中,M为总采样点个数。将所有的观测值和离散网格式代入式(4),即可得到如下方程式:
\begin{split} \left[ {\begin{array}{*{20}{c}} {{{\bar \tau }_1}} \\ {{{\bar \tau }_2}} \\ \vdots \\ {{{\bar \tau }_M}} \end{array}} \right] = \;&\left[ {\begin{array}{*{20}{c}} {{A_{1,1}}}&{{A_{1,2}}}& \cdots &{{A_{1,N}}} \\ {{A_{2,1}}}&{{A_{2,2}}}& \cdots &{{A_{2,N}}} \\ \vdots & \vdots & \ddots & \vdots \\ {{A_{M,1}}}&{{A_{M,2}}}& \cdots &{{A_{M,N}}} \end{array}} \right]\\ \;&\times\left[ {\begin{array}{*{20}{c}} {\Delta d\sqrt {\varepsilon ({{{\boldsymbol{r}}}_1})} /{\text{c}}} \\ {\Delta d\sqrt {\varepsilon ({{{\boldsymbol{r}}}_2})} /{\text{c}}} \\ \vdots \\ {\Delta d\sqrt {\varepsilon ({{{\boldsymbol{r}}}_N})} /{\text{c}}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{\eta _1}} \\ {{\eta _2}} \\ \vdots \\ {{\eta _M}} \end{array}} \right] \end{split} (5) 使用向量和矩阵代替,式(5)可以表示为
{{\boldsymbol{P}}={\boldsymbol{AO}}}+{\boldsymbol{\eta}} (6) 其中, \eta 表示观测噪声,A为M \times N阶的线性层析投影矩阵,用于描述观测数据与离散化场景之间的映射关系,具体取值为
{A_{m,n}} = \left\{ {\begin{array}{*{20}{c}} {1,}&{{{{\boldsymbol{r}}}_n} \in {{{\boldsymbol{L}}}_m}} \\ {0,}&{{{{\boldsymbol{r}}}_n} \notin {{{\boldsymbol{L}}}_m}} \end{array}} \right. (7) 其中, {{{\boldsymbol{L}}}_m} 表示第m次测量时收发天线间的直达波路径,若第m次观测时,第n个网格在收发设备间的连线上,则{A_{m,n}} = 1,否则{A_{m,n}} = 0,如图2所示。
本节基于线性近似层析投影矩阵A,实现了高维空间矩阵与低维观测数据之间的映射。由于未知场景中多径信号的传播路径难以被估计与建模,线性近似层析模型仅考虑了收发天线之间的直接路径信息。因此,需要从原始接收信号中精确估计直达波时延,从而构建准确的模型方程。
3. 基于多域联合的直达波估计算法
电磁波在未知区域内传播过程中,遇到障碍物会发生散射现象,导致接收信号是直达波信号与多路径信号的叠加。本文根据层析成像扫描特点提出了一种基于多域联合的直达波估计算法,实现了多径信号与直达波信号的分离。
令发射信号为超宽带步进频信号,在0 \le t \le T的脉冲周期内,其形式为
\begin{split} S(t) =\;& \sum\limits_{k = 0}^{K - 1} {\text{exp}}\left( {{\mathrm{j}}2\pi \left( {{f_0} + k\Delta f} \right)t} \right) \\ & \times {\text{rect}}\left( {\frac{{t - k{T_0} - {T_0}/2}}{{{T_0}}}} \right) \end{split} (8) 其中,{f_0}为起始频率,K为频点数, \Delta f 为步进频率,每个频点的持续时间为{T_0}, \text{rect}(\cdot) 表示中心为0、宽度为1的矩形窗。
在第m次采样位置处的接收信号 {R_m}(t) 可以被认为是发射信号经过多条传播路径调制后不同幅度、不同时延的信号之和,即
\begin{split} R(m,t) =\;& \sum\limits_{i = 1}^I {\sum\limits_{k = 0}^{K - 1} {{\alpha _{m,i}}{\text{exp}}\left( {{\mathrm{j}}2\pi \left( {{f_0} + k\Delta f} \right)\left( {t - {\tau _{m,i}}} \right)} \right)} } \\ & \times {\text{rect}}\left( {\frac{{t - \tau - k{T_0} - {T_0}/2}}{{{T_0}}}} \right)\\[-1pt] \end{split} (9) 其中,I为总传播路径数,{\alpha _{m,i}}和{\tau _{m,i}}分别表示第m次采样时第i条传播路径的幅度衰减比和传播时延。对回波信号执行混频、采样、补零以及傅里叶变换操作,即可得到第m次采样信号的脉冲压缩结果[3]:
\begin{split} d\left( {m,t} \right) =\;& \sum\limits_{i = 1}^I {{\alpha _{m,i}}B{\text{sinc}}\left[ {B\left( {t - {\tau _{m,i}}} \right)} \right]} \\ & \times {\text{exp}}\left( {{\mathrm{j}}2\pi {f_{\mathrm{c}}}\left( {t - {\tau _{m,i}}} \right)} \right) \end{split} (10) 其中,B为发射信号等效脉冲的带宽,{f_{\mathrm{c}}}为等效脉冲的中心频率。取第i条路径信号所处的慢时间序列,可得到多普勒频率项:
{f_i}\left( m \right) = {\text{exp}}\left( {{\mathrm{j}}2\pi {f_{\mathrm{c}}}\Delta {\tau _i}\left( m \right)} \right) (11) 其中,\Delta {\tau _i}\left( m \right)表示第i条路径信号随着采样位置m变化的到达时延差。图3(b)给出了图3(a)仿真场景下的0°路径扫描信号的脉冲压缩结果,其中黑色虚线表示直达波信号,其余回波均为多径信号。从图3中可以看出,多径信号和直达波信号传播路径不同,导致其传播时延也不同。此外,由于建筑结构在空间中的分布通常存在一定的连续性,相邻采样点的直达波的到达时延差接近0;而扫描过程中障碍物与收发天线间存在相对运动,多径信号的到达时延差明显不为0。结合式(11),对于一小段区域内的连续采样信号,可以得到如下结论:
(1) 在墙体结构连续部分,直达波路径信号的多普勒频率接近0;
(2) 多径信号的多普勒频率不为0。
图4给出了第96个采样点的脉冲压缩结果,从图4中可以看出接收信号中存在明显的多径干扰,且其幅度相较于直达波更强,导致直接使用MAE[18]对直达波进行估计时无法得到正确的结果。
结合上述结论,对第91~101个采样点的脉冲压缩结果在采样点维度进行傅里叶变换,得到距离-多普勒(Range-Doppler, R-D)图像,如图5所示。从R-D谱可以看出,直达波在多普勒零频附近实现了积累,而多径信号积累在非零部分,和理论分析结论相符合。最后,取出零频处的一维距离像,如图6所示,可以看到多径信号得到有效抑制。此时,通过寻找最大峰值即可得到直达波时延。
综上,基于多域联合的直达波估计方法的主要处理步骤介绍如下:
(1) 窗函数截取
对某一采样点的接收信号进行直达波估计时,我们假设与其距离越近的采样点的直达波具有更高的相似性。为此,本文采用三角窗在采样点维度对该采样点周围的信号进行截取。通过应用三角窗,在保留局部信号直达波相似特征的同时,减轻墙体结构突变引起的不利影响。三角窗的表达式如下:
w\left( {m,{m_0},t} \right) = \left\{ \begin{aligned} & {1 - \frac{{\left| {m - {m_0}} \right|}}{{{w_d}}}{\text{ , }}\left| {m - {m_0}} \right| \le {w_d}} \\ & {0{\text{ , }}\left| {m - {m_0}} \right| > {w_d}} \end{aligned} \right. (12) 其中,{m_0}表示窗函数的中心位置, {w_d} 表示窗函数单侧宽度。
随后,利用三角窗对采样信号加窗处理,即可得到:
y\left( {m,t} \right) = d\left( {m,t} \right) \cdot w\left( {m,{m_0},t} \right),\;\; m{{ = 1,\,2,}}\, \cdots {\text{,}}\,M (13) (2) 采样点傅里叶变换
完成窗函数截取后,对加窗截取后的信号在采样点维度进行傅里叶变换得到R-D谱:
Y\left( {\omega ,t} \right) = {\text{FFT[}}y\left( {m,t} \right){\text{]}} (14) 其中, Y\left( {\omega ,t} \right) 表示R-D谱,\omega 表示多普勒频率。
(3) 直达波时延估计
取R-D谱的0频信号,并寻找最大峰值的索引即可得到直达波时延估计结果:
{\bar \tau _{{m_0}}} = \mathop {{\text{argmax}}}\limits_t \left[ {Y\left( {0,t} \right)} \right] (15) 至此,实现了采样点{m_0}处的直达波时延估计。以\Delta m的步长修改窗函数中心位置,重复上述操作即可完成对各个采样位置处的直达波时延估计。
4. 总变分约束的投影矩阵自适应修正代数重建算法
受限于实际环境与探测效率,建筑布局层析成像系统往往只能从有限的角度进行探测,获得的观测数据数量不可避免地远少于未知区域的网格数量,导致式(6)为欠定方程,直接对其求解较为困难。为此,本文提出了一种总变分约束的投影矩阵自适应修正代数重建算法。具体地,首先利用代数重建技术[20]使得重建结果的正投影值逼近观测数据;其次,施加总变分约束[21]以平滑图像;然后,根据图像的暂时解更新投影矩阵[22],用于下一次迭代成像;最终,通过反复迭代实现精确的反演成像。
4.1 代数重建技术
代数重建技术(Algebraic Reconstruction Technique, ART)的基本思想是将层析投影矩阵与当前成像结果相乘以获取投影值。然后,根据投影值与观测数据之间的差值计算修正值。最后,根据修正值更新成像网格的灰度值,以得到逼近观测数据的成像结果。因此,代数重建技术可被总结为如下步骤:
步骤1 计算当前图像的投影值
P_m^k = \sum\limits_{n = 1}^N {{A_{m,n}}O_n^k}, \quad m = 1,2, \cdots ,M (16) 其中, O_n^k 表示第k次迭代时第n个网格的灰度值, {A_{m,n}} 表示第m次测量时第n个网格的投影系数, P_m^k 表示第k次迭代时第m个投影矩阵对应的投影值。
步骤2 更新差值
\Delta C_{n}^k = \sum\limits_{m = 1}^M {\left( {{P_{om}} - P_m^k} \right)\frac{{{A_{m,n}}}}{{\displaystyle\sum\limits_{n = 1}^N {{A_{m,n}}} }}}, \quad n = 1,2, \cdots ,N (17) 其中, {P_{om}} 表示第m个采样位置处的观测数据。
步骤3 更新图像
O_n^{k + 1} = O_n^k + \Delta C_n^k, \quad n = 1,2, \cdots ,N (18) 经过上述3个步骤即完成了一次代数重建迭代。
4.2 基于总变分约束的建筑布局图像平滑
对于建筑布局重构任务,期望得到边缘完整、结构平滑的图像,以清晰地反映墙体信息。本文采用了分裂Bregman TV算法[23],通过引入分裂变量和Bregman迭代,有效施加总变分约束,实现了快速收敛、适应大规模数据的图像平滑与去噪。总变分模型表示如下:
\mathop {\min }\limits_u {\left\| {{\nabla _x}{{\boldsymbol{u}}}} \right\|_1} + {\left\| {{\nabla _y}{{\boldsymbol{u}}}} \right\|_1} + \frac{\mu }{2}\left\| {{{\boldsymbol{u}}} - {{\boldsymbol{O}}}} \right\|_2^2 (19) 其中, {\nabla _x} , {\nabla _y} 分别表示水平与垂直方向的差分算子,初次迭代中O为代数重建结果。为了便于对式(19)求解,引入松弛变量将式(19)转化为
\begin{split} & \mathop {\min }\limits_u {\left\| {{{{\boldsymbol{d}}}_x}} \right\|_1} + {\left\| {{{{\boldsymbol{d}}}_y}} \right\|_1} + \frac{\mu }{2}\left\| {{{\boldsymbol{u}}} - {{\boldsymbol{O}}}} \right\|_2^2 \\ & {\text{ s}}{\text{.t}}{\text{. }}{{{\boldsymbol{d}}}_x} = {\nabla _x}{{\boldsymbol{u}}},{{{\boldsymbol{d}}}_y} = {\nabla _y}{{\boldsymbol{u}}} \end{split} (20) 根据Bregman迭代,引入辅助向量{{\boldsymbol{b}}}_x^k, {{\boldsymbol{b}}}_y^k可以得到:
\begin{split} & \mathop {\min }\limits_{{{\boldsymbol{d}}_x},{{\boldsymbol{d}}_y},u} {\left\| {{{{\boldsymbol{d}}}_x}} \right\|_1} + {\left\| {{{{\boldsymbol{d}}}_y}} \right\|_1} + \frac{\mu }{2}\left\| {{{\boldsymbol{u}}} - {{\boldsymbol{O}}}} \right\|_2^2 \\ & +\frac{\lambda }{2}\left\| {{{{\boldsymbol{d}}}_x} - {\nabla _x}{{\boldsymbol{u}}} - {{\boldsymbol{b}}}_x^k} \right\| + \frac{\lambda }{2}\left\| {{{{\boldsymbol{d}}}_y} - {\nabla _y}{{\boldsymbol{u}}} - {{\boldsymbol{b}}}_y^k} \right\| \end{split} (21) 为了获取u的最优解,需求解以下子问题:
\begin{split} {u^{k + 1}} = \;&\mathop {\arg \min }\limits_{\boldsymbol{u}} \frac{\mu }{2}\left\| {{{\boldsymbol{u}}} - {{\boldsymbol{O}}}} \right\|_2^2 + \frac{\lambda }{2}\left\| {{{{\boldsymbol{d}}}_x} - {\nabla _x}{{\boldsymbol{u}}} - {{\boldsymbol{b}}}_x^k} \right\| \\ \;& + \frac{\lambda }{2}\left\| {{{{\boldsymbol{d}}}_y} - {\nabla _y}{{\boldsymbol{u}}} - {{\boldsymbol{b}}}_y^k} \right\| \\[-1pt] \end{split} (22) 此问题通过Gauss-Siedel方法求解[23],表达式如下:
\begin{split} u_{i,j}^{k + 1} = \;&\frac{\lambda }{{\mu + 4\lambda }}\bigr(u_{i + 1,j}^k + u_{i - 1,j}^k + u_{i,j + 1}^k + u_{i,j - 1}^k \\ \;& +d_{x,i - 1,j}^k - d_{x,i,j}^k + d_{y,i,j - 1}^k - d_{y,i,j}^k - b_{x,i - 1,j}^k \\ \;& + b_{x,i,j}^k - b_{y,i,j - 1}^k + b_{y,i,j}^k \bigr) + \frac{\lambda }{{\mu + 4\lambda }}{O_{i,j}} \\[-1pt] \end{split} (23) 其中,i,j分别表示当前计算的像素点在二维图像x,y方向上的索引值。
此外,式(22)中{{\boldsymbol{b}}}_x^k, {{\boldsymbol{b}}}_y^k根据式(28)进行迭代更新:
\left.\begin{aligned} & {{\boldsymbol{b}}}_x^{k + 1} = {{\boldsymbol{b}}}_x^k + {\nabla _x}{{{\boldsymbol{u}}}^{k + 1}} - {{\boldsymbol{d}}}_x^{k + 1} \\ & {{\boldsymbol{b}}}_y^{k + 1} = {{\boldsymbol{b}}}_y^k + {\nabla _y}{{{\boldsymbol{u}}}^{k + 1}} - {{\boldsymbol{d}}}_y^{k + 1} \end{aligned}\right\} (24) 进一步,利用广义收缩公式[23]迭代更新 {{{\boldsymbol{d}}}_x} 与 {{{\boldsymbol{d}}}_y} :
\left.\begin{aligned} & {{\boldsymbol{d}}}_x^{k + 1} = {\text{shrink}}\left( {{\nabla _x}{{{\boldsymbol{u}}}^{k + 1}} + {{\boldsymbol{b}}}_x^k,\frac{1}{\lambda }} \right) \\ & {{\boldsymbol{d}}}_y^{k + 1} = {\text{shrink}}\left( {{\nabla _y}{{{\boldsymbol{u}}}^{k + 1}} + {{\boldsymbol{b}}}_y^k,\frac{1}{\lambda }} \right) \end{aligned}\right\} (25) 通过反复迭代式(23)至式(25),即可实现基于总变分约束的图像平滑与去噪。
4.3 投影矩阵更新
投影矩阵的精确表征是高质量层析成像的关键。由于空气的相对介电常数约等于1,代入式(4)可以得到空场景下的直达波的传播时延为
{\tau _{{\text{empty}}}} = \int\nolimits_{{{\mathbb{L}}_{\overrightarrow {{\mathrm{TX}}} \to \overrightarrow {{\mathrm{RX}}} }}} {\frac{1}{{\mathrm{c}}}} {\text{ d}}{{\boldsymbol{r}}} (26) 将直达波时延估计结果减去空场景直达波时延后可以得到由障碍物引起的额外传播时延,表示如下:
{\tau _{{\text{extra}}}} = \bar \tau - {\tau _{{\text{empty}}}} = \int\nolimits_{{{\mathbb{L}}_{\overrightarrow {{\text{TX}}} \to \overrightarrow {{\text{RX}}} }}} {\frac{{\sqrt {\varepsilon \left( {{\boldsymbol{r}}} \right)} - 1}}{\mathrm{{c}}}} {\text{ d}}{{\boldsymbol{r}}} (27) 将观测向量中的直达波估计结果\bar \tau 替换为额外传播时延 {\tau _{{\text{extra}}}} 。此时,额外传播时延被认为是由传播路径中的障碍物所引起的,避免观测数据被投影至空气网格能够有效提升建筑布局图像重建的精度。具体地,通过遍历图像暂时解的各个像素点,找到灰度值接近0的网格,并根据网格位置修正投影矩阵系数,数学表达式如下:
{A_{mn}} = 0\quad {\text{s}}.{\text{t}}.{\text{ }}{O_n} < \eta ,\;\; {m = 1,2, \cdots ,M} (28) 其中, \eta 为判断网格属性为空气的阈值。该阈值采用大津算法(Otsu’s Method)[24],通过最大化类间方差的方式自动确定。随后,更新的投影矩阵将被用于下一次迭代成像。
综上,所提算法的处理流程如算法1所示。
1 PMAM-ART-TV算法流程1. PMAM-ART-TV algorithm flow输入: P, A,初始化 {\boldsymbol{O}}=0 ,外部停止标准{\varepsilon _o},内部停止标准
{\varepsilon _i},外部迭代次数t = 0,外部最大迭代次数 {T_o} ,内部最大迭代
次数{T_i}输出: \tilde{{\boldsymbol{O}}}={\boldsymbol{O}}^{t} repeat 1、代数重建迭代: k = 0 求解式(16)、式(17)更新 \Delta C_{}^{k + 1} 求解式(18)更新 {{O}}_n^{k + 1} k = k + 1 直到k = {T_i}或者 \Vert {{\boldsymbol{O}}}^{k+1}-{{\boldsymbol{O}}}^{k}\Vert \le {\varepsilon }_{i} ,输出 {{{\boldsymbol{O}}}_{{\mathrm{ART}}}} 2、总变分约束迭代: k = 0, {{\boldsymbol{O}}} = {{{\boldsymbol{O}}}_{{\mathrm{ART}}}} , {{{\boldsymbol{u}}}^k} = 0 求解式(23)更新 {{{\boldsymbol{u}}}^{k + 1}} 求解式(24)更新{{\boldsymbol{b}}}_x^{k + 1}, {{\boldsymbol{b}}}_y^{k + 1} 求解式(25)更新 {{\boldsymbol{d}}}_x^{k + 1} , {{\boldsymbol{d}}}_y^{k + 1} {{\boldsymbol{O}}} = {{{\boldsymbol{u}}}^{k + 1}} , k = k + 1 直到k = {T_i}或者 \Vert {{\boldsymbol{O}}}^{k+1}-{{\boldsymbol{O}}}^{k}\Vert \le {\varepsilon }_{i} ,输出O 3、更新投影矩阵 求解式(28)更新A, {{{\boldsymbol{O}}}^{t + 1}} , t = t + 1 until直到t = {T_o}或者 \Vert {{\boldsymbol{O}}}^{t+1}-{{\boldsymbol{O}}}^{t}\Vert \le {\varepsilon }_{o} 5. 仿真验证
为验证所提成像模型与方法的有效性,利用有限时域差分电磁仿真软件Gprmax[25]开展了层析成像仿真实验,仿真场景和仿真参数如图7和表1所示,其中收发天线沿着0°, 45°, 90°和135° 4组采样路径对场景进行扫描。
表 1 仿真参数Table 1. Simulation parameters类型 仿真参数 数值 信号参数 发射信号 步进频信号 中心频率 1.5 GHz 带宽 600 MHz 单频点持续时间 100 μs 扫描参数 采样路径长度 10 m 采样路径数目 4组 采样间隔 0.05 m 场景参数 场景尺寸 2 m×2 m 墙体厚度 0.20 m 相对介电常数 4 电导率 0.01 S/m 首先,仿真场景0°角下的观测数据估计结果如图8所示。为了体现所提直达波估计算法的优越性,对比了MAE算法的直达波估计结果与RSSI。上述估计结果均被归一化,以确保不同类型估计结果之间的可比性。从图8中可以看出,RSSI和MAE估计结果受到严重的多径干扰,导致其估计结果存在明显的抖动,与理想值相差较大。值得注意的是,所提算法充分考虑了多径信号在时域和多普勒域的分布特性,通过保留零频分量,实现了多径干扰抑制,使得其时延估计结果与理论值最为匹配。
为了进一步客观评价观测数据的估计精度,引入了均方根误差(Root Mean Squared Error, RMSE)指标用于评估观测数据与理想值的误差情况,其计算公式如(29)所示。
{\text{RMSE}} = \sqrt {\frac{1}{M}\sum\limits_{m = 1}^M {{{\left( {{{\tilde P}_m} - {\text{ }}{P_m}} \right)}^2}} } (29) 其中, {\tilde P_m} 表示归一化的观测数据估计结果, {P_m} 表示理想值。由式(29)可知,当RMSE越小,表示观测数据与理想值的误差越小。仿真所获取的观测数据的RMSE计算结果如表2所示,可以发现RSSI观测数据误差最大;而MAE算法在一定程度上分离了多径和直达波信号,其估计结果误差更小;所提直达波估计方法能够有效抑制多径从而得到误差最小的观测数据估计结果。
基于RSSI, MAE和MD-DE方法估计的观测数据,本文进一步开展了成像性能评估,并引入主流的ART算法、Tikhonov算法[26]以及TV算法[15]进行了比较,成像结果如图9所示。图9(a)给出了不同成像算法基于RSSI观测数据的成像结果,可以发现观测数据受到严重的多径干扰,所有算法的重建结果均存在墙体边缘模糊、伪影严重的问题。相比较而言,图9(b)中基于MAE观测数据的成像结果,其墙体的结构和边缘更为清晰,且成像结果中噪声较少;但是由于MAE估计结果依然无法完全去除多径的影响,部分直达波时延估计不正确,导致成像结果中出现了明显的墙体展宽和结构缺失问题。由此可知,获取精确的观测数据是实现高质量层析成像的必要条件。
图9(c)给出了基于所提出直达波估计算法得到的成像结果,由于观测数据远小于成像网格数量,ART算法和Tikhonov算法无法精确求解该问题,成像结果中墙体伪影依然严重,无法获取正确的墙体结构;而总变分约束与建筑结构分布的先验特性较为吻合,使得TV算法显著改善了成像质量,得到了墙体平滑、边缘清晰的建筑布局图像;所提出的成像算法在总变分约束的基础上,根据图像暂时解对投影矩阵进行了修正,使得观测数据被更精确地投影到成像结果中,进一步抑制了图像中的伪影,获得了更为清晰的建筑布局图像,证明了所提算法的优越性。
为了进一步客观评价各种算法的成像质量,本文引入了结构相似性参数(Structural Similarity, SSIM)[27], SSIM指标计算公式如下:
{\text{SSIM}} = \frac{{(2{\mu _1}{\mu _2} + {c_1})(2{\sigma _{1,2}} + {c_2})}}{{(\mu _1^2 + \mu _2^2 + {c_1})(\sigma _1^2 + \sigma _2^2 + {c_2})}} (30) 其中, {\mu _i} , {\sigma _i} , {\sigma _{12}} 和 {c_i}\left( {i = 1,2} \right) 分别为平均标准差、协方差和比例因子。该指标用于评估成像结果与原始图像间的相关性,若图像间相似度越高,则SSIM评价指标越接近1。不同成像结果的SSIM指标如表3所示,可以看出利用所提方法获得的成像结果与原图具有最高的相似度,表明了所提方法的优越性。
表 3 仿真成像结果SSIM指标对比(%)Table 3. Comparison of SSIM indicators in simulation imaging results (%)方法 ART Tikhonov TV PMAM-ART-TV RSSI 30.8 28.7 48.9 68.5 MAE 34.1 51.4 72.8 86.2 MD-DE 40.6 58.9 88.9 91.2 此外,为了比较层析成像模式和传统合成孔径成像模式的差异,给出了相同仿真场景和信号参数条件下的合成孔径成像结果。由于倾斜探测视角对合成孔径成像模式增益较为有限,改用了0°, 90°, 180°, 270° 4条平行路径对仿真场景进行探测扫描。图10(a)展示了经典后向投影算法(Back Projection, BP)得到的多视角融合成像结果,可以发现由于电磁波在穿透墙体时能量损失较大,导致成像结果中杂波显著、伪影严重。尽管基于低秩的TV正则化算法[28]能够有效地抑制杂波并恢复墙体结构,仍然存在墙体偏移和伪影问题,导致多视角融合结果中墙体边缘模糊、结构难以清晰识别。这些问题需要通过墙体参数估计和墙体补偿等技术手段来解决。相比之下,本文所提成像方法在不依赖墙体参数信息的情况下,能够获得较为准确的建筑结构信息,显示出明显的技术优势。
最后,为了探究噪声对所提建筑布局层析成像方法的影响,在接收信号中加入了以直达波信号功率为基准的高斯白噪声,并分别在–10~10 dB的信噪比下对不同的层析成像方法进行了对比分析,如图11(a)所示。由于RSSI信息取自信号的平均功率,且成像过程依赖于不同位置处RSSI信息的相对变化情况,故噪声对基于RSSI信息的成像结果的质量并未产生较大影响,其SSIM指标维持在一个较为平稳的水平。而对于MAE算法而言,当信噪比≤0 dB时,直达波信号被淹没在噪声当中,无法准确估计直达波信号,导致基于MAE算法的成像结果的SSIM指标在信噪比为0 dB左右时发生断崖式下跌,建筑布局重建质量较差。然而,所提出的MD-DE算法得益于对相邻多个采样点信号进行了相干处理,直达波信号的信噪比得到提升,使得直达波信号在–6 dB的信噪比情况下依然能被正确估计。此外,图11(b)给出了信噪比为–6 dB时,所提方法的建筑布局层析成像结果,可以看到即使在低信噪比的情况下,墙体的结构与布局依然能被清晰重构。
6. 实测验证
为验证所提建筑布局成像方法在真实场景的有效性,本文基于一发一收步进频雷达系统开展了相关实验。实验场景如图12所示,其中实验场景的尺寸为3.2{\text{ m}} \times {\text{3}}{\text{.2 m}},墙体厚度为0.12{\text{ m}},内部菱形障碍物为0.7{\text{ m}} \times 0.7{\text{ m}}。两个微带天线被置于两条移动滑轨上,分别发射步进频信号和接收透过实验场景的信号,具体雷达参数如表4所示。在实验过程中,收发天线将沿着0°, 45°, 90°和135°4条路径完成扫描,并基于透射信号数据实现探测场景重构。实测层析成像结果与SSIM指标计算结果分别如图13和表5所示。
表 4 雷达系统参数Table 4. Radar system parameters参数名称 数值 中心频率 1.9 GHz 带宽 600 MHz 频率步进 2 MHz 频点持续时间 100 μs 表 5 实测成像结果SSIM指标对比(%)Table 5. Comparison of SSIM indicators in actual imaging results (%)方法 ART Tikhonov TV PMAM-ART-TV RSSI 31.4 30.6 31.2 31.7 MAE 43.9 49.8 73.6 75.5 MD-DE 46.3 56.6 78.6 81.7 从成像结果可以发现,基于RSSI得到的成像结果质量均较差,无法分辨墙体和障碍物结构。直达波时延信息可以在一定程度上抑制多径的影响,并更精确地通过线性模型表征观测数据与场景中物理参数的关系,从而提升成像质量。然而,MAE无法准确估计直达波,导致其成像结果中存在墙体与障碍物成像结果不平滑、结构与伪影难以区分等问题,如图13(b)所示。本文所提方法在建筑布局重构性能上表现更佳,得到了伪影更少、建筑结构更清晰的层析成像结果。同时,所提算法重建结果的SSIM指标高于其他结果,也证明了所提算法在实测数据下的有效性与优越性。
7. 结语
本文提出了一种基于多域联合直达波估计的建筑布局层析成像方法。相较于以往的层析成像方法,该方法实现了观测数据的准确估计和稀疏数据下的高质量重构。首先,分析和建模了层析成像模式下采样信号中直达波信号与多径信号在时域和多普勒域的分布差异;然后,提出了基于多域联合的直达波估计算法,有效抑制了多径干扰,提升了直达波估计精度;此外,提出了总变分约束的投影矩阵自适应修正代数重建算法,能够在平滑图像同时提升代数迭代反演的精度。仿真结果和实验结果验证了所提方法的性能。
-
1 PMAM-ART-TV算法流程
1. PMAM-ART-TV algorithm flow
输入: P, A,初始化 {\boldsymbol{O}}=0 ,外部停止标准{\varepsilon _o},内部停止标准
{\varepsilon _i},外部迭代次数t = 0,外部最大迭代次数 {T_o} ,内部最大迭代
次数{T_i}输出: \tilde{{\boldsymbol{O}}}={\boldsymbol{O}}^{t} repeat 1、代数重建迭代: k = 0 求解式(16)、式(17)更新 \Delta C_{}^{k + 1} 求解式(18)更新 {{O}}_n^{k + 1} k = k + 1 直到k = {T_i}或者 \Vert {{\boldsymbol{O}}}^{k+1}-{{\boldsymbol{O}}}^{k}\Vert \le {\varepsilon }_{i} ,输出 {{{\boldsymbol{O}}}_{{\mathrm{ART}}}} 2、总变分约束迭代: k = 0, {{\boldsymbol{O}}} = {{{\boldsymbol{O}}}_{{\mathrm{ART}}}} , {{{\boldsymbol{u}}}^k} = 0 求解式(23)更新 {{{\boldsymbol{u}}}^{k + 1}} 求解式(24)更新{{\boldsymbol{b}}}_x^{k + 1}, {{\boldsymbol{b}}}_y^{k + 1} 求解式(25)更新 {{\boldsymbol{d}}}_x^{k + 1} , {{\boldsymbol{d}}}_y^{k + 1} {{\boldsymbol{O}}} = {{{\boldsymbol{u}}}^{k + 1}} , k = k + 1 直到k = {T_i}或者 \Vert {{\boldsymbol{O}}}^{k+1}-{{\boldsymbol{O}}}^{k}\Vert \le {\varepsilon }_{i} ,输出O 3、更新投影矩阵 求解式(28)更新A, {{{\boldsymbol{O}}}^{t + 1}} , t = t + 1 until直到t = {T_o}或者 \Vert {{\boldsymbol{O}}}^{t+1}-{{\boldsymbol{O}}}^{t}\Vert \le {\varepsilon }_{o} 表 1 仿真参数
Table 1. Simulation parameters
类型 仿真参数 数值 信号参数 发射信号 步进频信号 中心频率 1.5 GHz 带宽 600 MHz 单频点持续时间 100 μs 扫描参数 采样路径长度 10 m 采样路径数目 4组 采样间隔 0.05 m 场景参数 场景尺寸 2 m×2 m 墙体厚度 0.20 m 相对介电常数 4 电导率 0.01 S/m 表 2 观测数据误差对比
Table 2. Comparison of observation data errors
表 3 仿真成像结果SSIM指标对比(%)
Table 3. Comparison of SSIM indicators in simulation imaging results (%)
方法 ART Tikhonov TV PMAM-ART-TV RSSI 30.8 28.7 48.9 68.5 MAE 34.1 51.4 72.8 86.2 MD-DE 40.6 58.9 88.9 91.2 表 4 雷达系统参数
Table 4. Radar system parameters
参数名称 数值 中心频率 1.9 GHz 带宽 600 MHz 频率步进 2 MHz 频点持续时间 100 μs 表 5 实测成像结果SSIM指标对比(%)
Table 5. Comparison of SSIM indicators in actual imaging results (%)
方法 ART Tikhonov TV PMAM-ART-TV RSSI 31.4 30.6 31.2 31.7 MAE 43.9 49.8 73.6 75.5 MD-DE 46.3 56.6 78.6 81.7 -
[1] 金添, 宋勇平, 崔国龙, 等. 低频电磁波建筑物内部结构透视技术研究进展[J]. 雷达学报, 2021, 10(3): 342–359. doi: 10.12000/JR20119.JIN Tian, SONG Yongping, CUI Guolong, et al. Advances on penetrating imaging of building layout technique using low frequency radio waves[J]. Journal of Radars, 2021, 10(3): 342–359. doi: 10.12000/JR20119. [2] AHMAD F and AMIN M G. Noncoherent approach to through-the-wall radar localization[J]. IEEE Transactions on Aerospace and Electronic Systems, 2006, 42(4): 1405–1419. doi: 10.1109/TAES.2006.314581. [3] 崔国龙. 超宽带穿墙雷达合成孔径成像算法研究与实现[D]. [硕士论文], 电子科技大学, 2008: 6–16.CUI Guolong. Research and implementation of synthetic aperture imaging algorithm for ultra-wideband through-the-wall radar[D]. [Master dissertation], University of Electronic Science and Technology of China, 2008: 6–16. [4] NKWARI P K M, SINHA S, and FERREIRA H C. Through-the-wall radar imaging: A review[J]. IETE Technical Review, 2018, 35(6): 631–639. doi: 10.1080/02564602.2017.1364146. [5] GUO Shisheng, CUI Guolong, KONG Lingjiang, et al. An imaging dictionary based multipath suppression algorithm for through-wall radar imaging[J]. IEEE Transactions on Aerospace and Electronic Systems, 2018, 54(1): 269–283. doi: 10.1109/TAES.2017.2756298. [6] CHEN Jiahui, GUO Shisheng, LUO Haolan, et al. Non-line-of-sight multi-target localization algorithm for driver-assistance radar system[J]. IEEE Transactions on Vehicular Technology, 2023, 72(4): 5332–5337. doi: 10.1109/TVT.2022.3227971. [7] CHEN Jiahui, LI Nian, GUO Shisheng, et al. Enhanced 3-D building layout tomographic imaging via tensor approach[J]. IEEE Transactions on Geoscience and Remote Sensing, 2024, 62: 5105614. doi: 10.1109/TGRS.2024.3391282. [8] LI Nian, CHEN Jiahui, GUO Shisheng, et al. Building layout tomographic imaging based on MIMO-UWB radar[J]. IEEE Transactions on Geoscience and Remote Sensing, 2024, 62: 5108713. doi: 10.1109/TGRS.2024.3435831. [9] NGUYEN L, RESSLER M, and SICHINA J. Sensing through the wall imaging using the Army Research Lab ultra-wideband synchronous impulse reconstruction (UWB SIRE) radar[C]. SPIE 6947, Radar Sensor Technology XII, Orlando, USA, 2008: 69470B. doi: 10.1117/12.776869. [10] 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. [11] 金添, 宋勇平. 超宽带雷达建筑物结构稀疏成像[J]. 雷达学报, 2018, 7(3): 275–284. doi: 10.12000/JR18031.JIN Tian and SONG Yongping. Sparse imaging of building layouts in ultra-wideband radar[J]. Journal of Radars, 2018, 7(3): 275–284. doi: 10.12000/JR18031. [12] ZHANG Yang, CHEN Jiahui, GUO Shisheng, 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. [13] MOSTOFI Y. Cooperative wireless-based obstacle/object mapping and see-through capabilities in robotic networks[J]. IEEE Transactions on Mobile Computing, 2013, 12(5): 817–829. doi: 10.1109/TMC.2012.32. [14] DEPATLA S, BUCKLAND L, and 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. [15] DEPATLA S, KARANAM C R, and 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] KARANAM C R and MOSTOFI Y. 3D through-wall imaging with unmanned aerial vehicles using WiFi[C]. The 16th ACM/IEEE International Conference on Information Processing in Sensor Networks, Pittsburgh, USA, 2017: 131–142. doi: 10.1145/3055031.3055084. [17] 张扬. 建筑布局结构透视微波层析成像理论与方法[D]. [博士论文], 电子科技大学, 2022: 20–82.ZHANG Yang. Theory and method of perspective microwave tomography of building layout structure[D]. [Ph.D. dissertation], University of Electronic Science and Technology of China, 2022: 20–82. [18] GUO Qichang, LI Yanlei, LIANG Xingdong, et al. A novel CT-mode through-the-wall imaging method based on time delay estimation[J]. IEEE Geoscience and Remote Sensing Letters, 2021, 18(8): 1381–1385. doi: 10.1109/LGRS.2020.3000423. [19] CHEN Jiahui, ZHANG Yang, LI Huquan, et al. Ultrawideband tomographic imaging in multipath-rich environment[J]. IEEE Geoscience and Remote Sensing Letters, 2022, 19: 8016305. doi: 10.1109/LGRS.2021.30961830. [20] YASWANTH K, BHATTACHARYA S, and KHANKHOJE U K. Algebraic reconstruction techniques for inverse imaging[C]. 2016 International Conference on Electromagnetics in Advanced Applications, Cairns, Australia, 2016: 756–759. doi: 10.1109/ICEAA.2016.7731509. [21] ANDREU F, CASELLES V, DÍAZ J I, et al. Some qualitative properties for the total variation flow[J]. Journal of Functional Analysis, 2002, 188(2): 516–547. doi: 10.1006/jfan.2001.3829. [22] LI Nian, ZHAN Yang, CHEN Jiahui, et al. Building layout tomographic imaging method with sparse scan angles[C]. The 2022 IEEE 8th International Conference on Computer and Communications, Chengdu, China, 2022: 2128–2132. doi: 10.1109/ICCC56324.2022.10066042. [23] GOLDSTEIN T and OSHER S. The split bregman method for l1-regularized problems[J]. SIAM Journal on Imaging Sciences, 2009, 2(2): 323–343. doi: 10.1137/080725891. [24] OTSU N. A threshold selection method from gray-level histograms[J]. IEEE Transactions on Systems, Man, and Cybernetics, 1979, 9(1): 62–66. doi: 10.1109/TSMC.1979.4310076. [25] WARREN C, GIANNOPOULOS A, and GIANNAKIS I. GprMax: Open source software to simulate electromagnetic wave propagation for ground penetrating radar[J]. Computer Physics Communications, 2016, 209: 163–170. doi: 10.1016/j.cpc.2016.08.020. [26] HAMILTON B R, MA Xiaoli, BAXLEY R J, et al. Radio frequency tomography in mobile networks[C]. 2012 IEEE Statistical Signal Processing Workshop, Ann Arbor, USA, 2012: 508–511. doi: 10.1109/SSP.2012.6319745. [27] CHEN Jiahui, GUO Shisheng, CUI Guolong, et al. Building layout reconstruction with transmissive and reflective signals[J]. IEEE Transactions on Geoscience and Remote Sensing, 2022, 60: 5116612. doi: 10.1109/TGRS.2022.3199270. [28] ZHAO Jifang, JIN Liangnian, and LIU Qinghua. Through-the-wall radar sparse imaging for building walls[J]. The Journal of Engineering, 2019, 2019(21): 7403–7405. doi: 10.1049/joe.2019.0541. -