The Multi-domain Union Clutter Suppression Algorithm Based on Robust Principal Component Analysis
-
摘要: 奇异值分解等传统算法在处理穿墙成像中的杂波抑制问题时，杂波消除不够彻底，目标成像质量不高，严重影响后续的目标检测与识别。为解决这一问题，该文基于鲁棒主成分分析理论，在回波域和图像域分别建立联合低秩稀疏模型，以光滑化快速交替线性化(SFAL)方法来求解模型，并对目标图像进行指数加权联乘多域图像融合处理，从而得到最终成像结果。仿真结果表明，该算法速度快、精度高，可有效改善目标成像质量，并能较好地满足穿墙成像的实时性和准确性要求。Abstract: In through-the-wall imaging, the clutter can not be eliminated completely through traditional algorithms, and affects seriously the subsequent target detection and recognition. To solve the problem, based on robust principal component analysis theory, a joint low-rank and sparse model is established in echo and image domain respectively. The models are solved by Smoothing Fast Alternating Linearization (SFAL) method. Then, the target images are dealt with exponentially weighted multiply multi-domain image fusion to obtain the final image. The simulation results indicate that the algorithm has great speed and accuracy with effective improvement on imaging quality of targets.
-
1. 引言
区别于自由空间中的其他雷达，穿墙成像雷达 (Through-the Wall Imaging Radar, TWIR)需要对墙后目标进行探测成像[1,2]。在这一过程中，由墙体反射造成的杂波信号会以较大的幅值来“遮蔽”目标信号，或者在时域上与目标信号交织重叠。因此，有效抑制墙体杂波[3]，是TWIR对墙后目标准确成像的重要前提。
诸如奇异值分解[4] (Singular Value Decomposition, SVD)一类的经典杂波抑制算法，仅能实现对主墙体杂波的消除，而残余杂波会显著降低目标成像质量。随着对TWIR的广泛应用和深入研究，目标检测、识别等工作对前期成像的要求越发严格，现有的杂波抑制算法难以达到穿墙成像的实时性和准确性条件。
近年来，机器学习理论被逐渐引入TWIR领域，如压缩感知[5,6](Compressive Sensing, CS)、矩阵补全[7](Matrix Completion, MC)以及鲁棒主成分分析[8](Robust Principal Component Analysis, RPCA)。其中，作为高光谱图像去噪及视频监控中的前景提取等研究方向的常用工具，RPCA[9]可以将数据矩阵唯一分解为低秩杂波矩阵与稀疏目标矩阵，即同时实现对杂波与目标的准确分离。对应于数据矩阵为穿墙回波矩阵或穿墙图像矩阵，RPCA可以分离得到其中的目标回波信号或目标像素分量，进而实现杂波抑制。所以，在本文中，基于RPCA理论，提出一种新的杂波抑制算法。
为满足穿墙成像的实时性和准确性要求，分别就算法的速度和精度进行设计改进：一是提出光滑化快速交替线性化方法来缩短RPCA问题的求解时间，从而加快算法的速度；二是多域联合处理，同时在回波域和图像域应用RPCA理论，通过对其得到的目标图像进行指数加权联乘融合处理来提高算法精度。
综上所述，本文的章节安排如下：第2节介绍RPCA原理；对回波域和图像域建模的工作在第3节完成；第4节对所提杂波抑制算法进行理论分析；在第5节设计穿墙场景仿真，通过与多种算法的性能对比，验证所提算法的速度和精度；最后在第6节总结全文。
2. RPCA原理
鲁棒主成分分析，是为增强主成分分析[10] (Principal Component Analysis, PCA)鲁棒性而产生的多项式时间算法。相较于经典的PCA，RPCA将低秩矩阵受到的轻微扰动扩展为任意大小的稀疏误差矩阵，突破了加性1维高斯噪声的限制。简单来讲，RPCA将数据矩阵
D∈Rl×s 表示为低秩矩阵A∈Rl×s 和稀疏矩阵E∈Rl×s 之和，并对其进行分解。作为一个优化问题，其数学表述为min (1) 其中，
{\text{rank(}}{\boldsymbol{A}}{\text{)}} 表示矩阵{\boldsymbol{A}} 的秩，{\left\| {\boldsymbol{E}} \right\|_0} 表示矩阵{\boldsymbol{E}} 的{l_0} 范数，\gamma > 0 为平衡低秩项和稀疏项的正则化参数。同时，为了保证矩阵的唯一分解，需满足如式(2)的条件[9]{\text{rank(}}{\boldsymbol{A}}{\text{)}} \le {{\rho }}\frac{s}{{{{(\lg l)}^2}}},{\text{ }}{\left\| {\boldsymbol{E}} \right\|_0} \le 0.1ls (2) 其中，
{{\rho }} > 0 为常数系数。然而，式(1)同时包含MC和
{l_0} 范数最小化这两个非确定性多项式(Non-deterministic Polynomial, NP)-hard的子问题，所以这一不连续的非凸优化问题并不存在理论上的有效解。不过，结合相关研究成果[11, 12]，以{l_1} 范数代替{l_0} 范数，以核范数代替矩阵的秩，便可以得到凸松弛后的优化问题为\mathop {\min }\limits_{{\boldsymbol{A}},{\boldsymbol{E}}} {\text{ }}{\left\| {\boldsymbol{A}} \right\|_ * } + \gamma {\left\| {\boldsymbol{E}} \right\|_1},{\text{ s.t. }}{\boldsymbol{D}} = {\boldsymbol{A}} + {\boldsymbol{E}} (3) 其中，
{\left\| {\boldsymbol{A}} \right\|_ * } 表示矩阵{\boldsymbol{A}} 的核范数，{\left\| {\boldsymbol{E}} \right\|_1} 表示矩阵{\boldsymbol{E}} 的{l_1} 范数。在这里，矩阵{\boldsymbol{A}} 和{\boldsymbol{E}} 要达到秩稀疏不相干的条件，以确保凸优化问题的正确求解。显然，在相当广泛的情况下，该条件都是成立的，即矩阵{\boldsymbol{A}} 和{\boldsymbol{E}} 可以被高概率地精确恢复[13]。对于式(3)中RPCA问题的求解，有以下几种常用的1阶方法：一是加速近端梯度[14](Accelerated Proximal Gradient, APG)方法，以线性化函数来实现对目标函数的局部逼近，使得单次迭代的计算量较小，但收敛速度较慢；二是精确增广拉格朗日乘子[15](Exact Augmented Lagrange Multipliers, EALM)方法，通过交替迭代增广拉格朗日函数来得到矩阵
{\boldsymbol{A}} 和{\boldsymbol{E}} ，可实现高阶线性收敛，不过每次迭代的SVD次数较多，大大影响了计算速度；三是非精确增广拉格朗日乘子[15](Inexact Augmented Lagrange Multipliers, IALM)方法，将EALM中的多轮交替最小化削减为1轮，可在保持原有收敛速度的基础上有效减少SVD次数。3. 联合低秩稀疏模型
3.1 回波域建模
首先，设定成像场景如下：TWIR以收发同置的均匀天线阵列对墙后目标进行探测成像，阵元数量为
N ，且天线阵列与匀质墙体表面平行，则第n 个天线阵元所接收的回波信号为u(n,t) = {u_{\text{w}}}(n,t) + {u_{{\text{tg}}}}(n,t) + {u_{{\text{no}}}}(n,t) + {u_{{\text{ant}}}}(n,t) (4) 其中，
{u_{\text{w}}}(n,t) 为墙体杂波信号，{u_{{\text{tg}}}}(n,t) 为目标回波信号，{u_{{\text{no}}}}(n,t) 为噪声信号，{u_{{\text{ant}}}}(n,t) 为天线耦合信号。在本文中，认为天线耦合波经预处理后被消除；且为便于建模分析，假设噪声信号对目标成像的影响相对有限，则回波信号可进一步简化为u(n,t) = {u_{\text{w}}}(n,t) + {u_{{\text{tg}}}}(n,t) (5) 随后，将式(5)的结果扩展至整个天线阵列，有2维回波矩阵
{\boldsymbol{U}} = {{\boldsymbol{U}}_{\text{w}}} + {{\boldsymbol{U}}_{{\text{tg}}}},{\text{ }}{\boldsymbol{U}} \in {{\boldsymbol{R}}^{T \times N}} (6) 其中，
T 为各天线通道信号的采样点数，{{\boldsymbol{U}}_{\text{w}}} 为杂波信号矩阵，{{\boldsymbol{U}}_{{\text{tg}}}} 为目标信号矩阵。在上述成像场景下，单就墙体杂波信号来讲，各天线阵元对应的墙体反射系数及信号传输路径是相同的，则其各自接收的墙体反射回波也是相同的[16-18]，即
{{\boldsymbol{U}}_{\text{w}}} 是低秩矩阵。从空间上看，墙后目标相对于整个探测区域是稀疏的；类似地，在回波域中，目标采样点相对于所有采样点也是稀疏的，即{{\boldsymbol{U}}_{{\text{tg}}}} 为稀疏矩阵。那么，根据式(3)，在回波域中建立联合低秩稀疏[19](Joint Low-Rank and Sparse, JLRS)模型，有\mathop {\min }\limits_{{{\boldsymbol{U}}_{\text{w}}},{{\boldsymbol{U}}_{{\text{tg}}}}} {\text{ }}{\left\| {{{\boldsymbol{U}}_{\text{w}}}} \right\|_ * } + \gamma {\left\| {{{\boldsymbol{U}}_{{\text{tg}}}}} \right\|_1}{\text{ subj }}{\boldsymbol{U}} = {{\boldsymbol{U}}_{\text{w}}} + {{\boldsymbol{U}}_{{\text{tg}}}} (7) 3.2 图像域建模
相较于回波域，图像域是对探测区域的直观表现，其应用RPCA理论的效果理应更好[20]。对于式(6)中的2维回波矩阵，以后向投影 (Back Projection, BP)算法进行成像处理，则图像中任一像素的幅度为
i(p,q) = {i_{\text{w}}}(p,q) + {i_{{\text{tg}}}}(p,q) (8) 其中，
{i_{\text{w}}}(p,q) 为位置(p,q) 处的像素杂波分量，{i_{{\text{tg}}}}(p,q) 为位置(p,q) 处的像素目标分量。将这一表示方法扩展至整个图像，有2维图像矩阵{\boldsymbol{I}} = {{\boldsymbol{I}}_{\text{w}}} + {{\boldsymbol{I}}_{{\text{tg}}}},{\text{ }}{\boldsymbol{I}} \in {{\boldsymbol{R}}^{P \times Q}} (9) 其中，
P ,Q 分别为图像中纵向与横向划分的像素点数量，{{\boldsymbol{I}}_{\text{w}}} 为杂波分量矩阵，{{\boldsymbol{I}}_{{\text{tg}}}} 为目标分量矩阵。实际上，BP算法是对经过时间补偿的各通道信号进行相干叠加。这样一来，各通道信号均以最大幅度在目标位置处叠加成像，目标的稀疏性大大增强，显然
{{\boldsymbol{I}}_{{\text{tg}}}} 是稀疏矩阵；而幅度不同的信号在墙体位置处叠加抵消，杂波的低秩性有所减弱，但仍可认为{{\boldsymbol{I}}_{\text{w}}} 是低秩矩阵。同样地，在图像域建立JLRS模型，有\mathop {\min }\limits_{{{\boldsymbol{I}}_{\text{w}}},{{\boldsymbol{I}}_{{\text{tg}}}}} {\text{ }}{\left\| {{{\boldsymbol{I}}_{\text{w}}}} \right\|_ * } + \gamma {\left\| {{{\boldsymbol{I}}_{{\text{tg}}}}} \right\|_1},{\text{ s.t. }}{\boldsymbol{I}} = {{\boldsymbol{I}}_{\text{w}}} + {{\boldsymbol{I}}_{{\text{tg}}}} (10) 4. 理论分析
4.1 光滑化快速交替线性化方法
根据上文对RPCA问题求解方法的分析，可以得到：线性化近似会减小单次迭代的计算量，而交替迭代会加速收敛，降低迭代复杂度。由此出发，为了加快RPCA问题求解速度，在这里提出一种光滑化快速交替线性化(Smoothing Fast Alternating Linearization, SFAL)方法。
以图像域JLRS模型为例，取定
{\left\| {{{\boldsymbol{I}}_{\text{w}}}} \right\|_ * } = f({\boldsymbol{x}}) ,\gamma {\left\| {{\boldsymbol{I}} - {{\boldsymbol{I}}_{\text{w}}}} \right\|_1} = g({\boldsymbol{x}}) ，则可将式(10)中的RPCA问题转化为凸函数之和的最小化问题[21]，即\mathop {\min }\limits_{{\boldsymbol{x}} \in {{\boldsymbol{R}}^{l \times s}}} F({\boldsymbol{x}}) = f({\boldsymbol{x}}) + g({\boldsymbol{x}}) (11) 然而，考虑到凸函数
f({\boldsymbol{x}}) ,g({\boldsymbol{x}}) 都是非光滑的，为了方便进行最小化求解，对其依次进行光滑化处理。区别于经典的光滑处理技术[22]，为了增强光滑化处理的可拓展性及减少光滑参数调整，优先选择光滑的迫近函数
c({\boldsymbol{v}}) 。c({\boldsymbol{v}}) 是{{\boldsymbol{R}}^{l \times s}} 上凸性参数{\kappa _b} = 1 的强凸函数，其梯度\nabla c({\boldsymbol{v}}) 满足Lipschitz连续，且Lipschitz指数{L_c} \ge {\kappa _c} 。令c({\boldsymbol{v}}) = \dfrac{1}{2}\left\| {\boldsymbol{v}} \right\|_2^2 ，并给定光滑参数\alpha > 0 ，则f({\boldsymbol{x}}) 的光滑近似函数[23]为\begin{split} {f_\alpha }({\boldsymbol{x}}) =& \mathop {\max }\limits_{{\boldsymbol{v}} \in {\text{dom}}({f^ * })} \left\{ {\left. { < {\boldsymbol{x}},{\boldsymbol{v}} > - {f^ * }({\boldsymbol{v}}) - \alpha c({\boldsymbol{v}})} \right\}} \right. \\ =& \frac{{\left\| {\boldsymbol{x}} \right\|_2^2}}{{2\alpha }} - {}^{{\alpha ^{ - 1}}}{f^ * }({\alpha ^{ - 1}}{\boldsymbol{x}}) \end{split} (12) 其中，
{f^ * } 为f 的Fenchel共轭函数，{\text{dom}}({f^ * }) 为函数{f^ * } 的定义域，{}^{{\alpha ^{ - 1}}}{f^ * }({\alpha ^{ - 1}}{\boldsymbol{x}}) 表示参数{\alpha ^{{{ - }}1}} 下凸函数{f^ * } 的Moreau包络，{\left\| \cdot \right\|_2} 表示{l_2} 范数。由于迫近函数的强凸性，式(12)中的优化问题有唯一的最优解[24]，即{J_\alpha }({\boldsymbol{x}}) = {\alpha ^{ - 1}}({\boldsymbol{x}} - {\text{pro}}{{\text{x}}_{\alpha f}}({\boldsymbol{x}})) (13) 其中，迫近算子
{\text{pro}}{{\text{x}}_{\alpha f}}({\boldsymbol{x}}) = \arg \mathop {\min }\limits_{\boldsymbol{r}} \Bigr\{ \left. f({\boldsymbol{r}}) + \dfrac{\alpha }{2}\left\| {{\boldsymbol{r}} - {\boldsymbol{x}}} \right\|_2^2:{\boldsymbol{r}} \in {\text{dom}}(f) \Bigr\} \right. 。对于光滑函数
{f_\alpha }({\boldsymbol{x}}) 来说，其梯度\nabla {f_\alpha }({\boldsymbol{x}}) = {J_\alpha }({\boldsymbol{x}}) 满足Lipschitz连续，且Lipschitz参数{L_{{f_\alpha }}} = {\alpha ^{ - 1}} 。同样地，对于函数g({\boldsymbol{x}}) 有光滑近似函数为{g_\beta }({\boldsymbol{x}}) = \frac{{\left\| {\boldsymbol{x}} \right\|_2^2}}{{2\beta }} - {}^{{\beta ^{ - 1}}}{g^ * }({\beta ^{ - 1}}{\boldsymbol{x}}) (14) 其中，
\beta > 0 为光滑参数，{g^ * } 为g 的Fenchel共轭函数。函数{g_\beta }({\boldsymbol{x}}) 的梯度\nabla {g_\beta }({\boldsymbol{x}}) = {\beta ^{ - 1}}({\boldsymbol{x}} - {\text{pro}}{{\text{x}}_{\beta g}}({\boldsymbol{x}})) 满足Lipschitz连续，且Lipschitz参数{L_{{g_\beta }}} = {\beta ^{ - 1}} 。此时，式(11)可转化为如下的光滑凸优化问题\mathop {\min }\limits_{{\boldsymbol{x}} \in {{\boldsymbol{R}}^{l \times s}}} F({\boldsymbol{x}}) = {f_\alpha }({\boldsymbol{x}}) + {g_\beta }({\boldsymbol{x}}) (15) 考虑到
{f_\alpha }({\boldsymbol{x}}) 和{g_\beta }({\boldsymbol{x}}) 是线性化可微函数，分别以其正则化形式构成目标函数F({\boldsymbol{x}}) 的2次近似，有\begin{split} G_{{\mu _g}}^{{g_\beta }}({\boldsymbol{x}},{{\boldsymbol{z}}^k}) =& {f_\alpha }({\boldsymbol{x}}) + {g_\beta }({{\boldsymbol{z}}^k}) + < \nabla {g_\beta }({{\boldsymbol{z}}^k}),{\boldsymbol{x}} - {{\boldsymbol{z}}^k} > \\ & + \frac{1}{{2{\mu _g}}}\left\| {{\boldsymbol{x}} - {{\boldsymbol{z}}^k}} \right\|_2^2\\[-18pt] \end{split} (16) \begin{split} G_{{\mu _f}}^{{f_\alpha }}({\boldsymbol{y}},{{\boldsymbol{x}}^k}) =& {g_\beta }({\boldsymbol{y}}) + {f_\alpha }({{\boldsymbol{x}}^k}) + < \nabla {f_\alpha }({{\boldsymbol{x}}^k}),{\boldsymbol{y}} - {{\boldsymbol{x}}^k} > \\ & + \frac{1}{{2{\mu _f}}}\left\| {{\boldsymbol{y}} - {{\boldsymbol{x}}^k}} \right\|_2^2 \\[-18pt] \end{split} (17) 其中，
{\mu _f} ,{\mu _g} 是惩罚项参数，{\boldsymbol{y}} ,{\boldsymbol{z}} 为迭代算子。忽略常数项，则式(16)、式(17)的极小值分别为\begin{split} & {{\boldsymbol{x}}^k} =\\ & \arg \mathop {\min }\limits_{\boldsymbol{x}} \left\{ {\left. {{f_\alpha }({\boldsymbol{x}}) + \frac{1}{{2{\mu _g}}}\left\| {{\boldsymbol{x}} - ({{\boldsymbol{z}}^k} - {\mu _g}\nabla {g_\beta }({{\boldsymbol{z}}^k}))} \right\|_2^2} \right\}} \right. \end{split} (18) \begin{split} & {{\boldsymbol{y}}^k} = \\ & \arg \mathop {\min }\limits_{\boldsymbol{y}} \left\{ {\left. {{g_\beta }({\boldsymbol{y}}) + \frac{1}{{2{\mu _f}}}\left\| {{\boldsymbol{y}} - ({{\boldsymbol{x}}^k} - {\mu _f}\nabla {f_\alpha }({{\boldsymbol{x}}^k}))} \right\|_2^2} \right\}} \right. \end{split} (19) 由于
{f_\alpha }({\boldsymbol{x}}) 和{g_\beta }({\boldsymbol{x}}) 为可分离函数，极小值的计算可简化为求解1维最小化问题，计算量较小。同时，为加速收敛，引入\eta 来更新算子{\boldsymbol{z}} ，有{\eta _{k + 1}} = \left(1 + \sqrt {1 + 4\eta _k^2} \right)/2 (20) {{\boldsymbol{z}}^{k{\text{ + }}1}} = {{\boldsymbol{y}}^k} + \frac{{{\eta _k} - 1}}{{{\eta _{k + 1}}}}({{\boldsymbol{y}}^k} - {{\boldsymbol{y}}^{k - 1}}) (21) 不难证明，上述方法的迭代复杂度为
O(\sqrt {L/\varepsilon } ) ，其中，L = {L_{{f_\alpha }}}{L_{{g_\beta }}}/({L_{{f_\alpha }}} + {L_{{g_\beta }}}) ，\varepsilon 表示最优解的阶数，且该迭代复杂度达到了1阶方法的理论极限[25]。同时，就单次迭代的计算复杂度而言，该算法同常见算法相近。最后，总结SFAL方法的流程如表1所示。表 1 SFAL方法输入：2维图像矩阵 {\mathbf{I}} \in {{\mathbf{R}}^{P \times Q}} ，凸函数f({\mathbf{x}}) = {\left\| {\mathbf{x}} \right\|_ * }，凸函数 g({\mathbf{x}}) = \gamma {\left\| {{\mathbf{I}} - {\mathbf{x}}} \right\|_1} ，正则化参数 \gamma = 1/\sqrt {\max (P,Q)} ； 输出：杂波分量矩阵{ {\mathbf{I} }_{\text{w} } } = { {\mathbf{x} }^{k{{ - } }1} }，目标分量矩阵{ {\mathbf{I} }_{ {\text{tg} } } } = {\mathbf{I} } - { {\mathbf{x} }^{k{{ - } }1} }。 (1) 初始化参数：\alpha = \beta = {10^{ - 6}}, {{\mathbf{x}}^0} = {{\mathbf{y}}^0} = {{\mathbf{z}}^1} = 0 , {\mu _f} = {\mu _g} = 1 , {\eta _1} = 1；k = 1。 (2) 根据式(12)和式(14)进行光滑化处理。 (3) 迭代解未收敛时执行步骤(4)到(7) (4) 根据式(18)和式(19)进行交替迭代； (5) 根据式(20)更新\eta ； (6) 根据式(21)更新{\mathbf{z}}； (7) k \leftarrow k + 1 (8) 结束循环 在以SFAL方法求解回波域中的JLRS模型时，仅需更改输入为2维回波矩阵
{\boldsymbol{U}} \in {{\boldsymbol{R}}^{T \times N}} ,\gamma = 1/\sqrt {\max (T,N)} ，其余参数无需调整。4.2 指数加权联乘多域图像融合处理
根据之前所建立的JLRS模型，不难发现，杂波和目标在各域中并非总是保持严格低秩或严格稀疏的。在回波域中，杂波的低秩性较好，目标的稀疏性较弱，在图像域中则恰恰相反。二者在性质上呈现出一种“互补关系”，而这种关系在各域得到的目标图像中体现得尤为明显。为了提高算法精度、改善目标成像质量，对目标图像进行指数加权联乘多域图像融合处理。
首先，对回波域和图像域中的目标图像作如下分析：在回波域中，JLRS模型中的目标信号矩阵为实矩阵(对原始回波数据取模求得)，隐藏了目标的相位信息，由此得到的目标图像相当于非相干BP成像结果，再加上残余杂波的影响，其方位向分辨率和聚焦效果都有所降低，但墙体主杂波抑制相对彻底；在图像域中，对原始回波成像直接进行分解，此时的目标图像为相干BP成像结果，其目标聚焦效果相对较好，但杂波分离不够彻底，且对目标形成一定的“遮蔽”效应。
经上述分析可得，相较于多通道或多角度的子图像融合，回波域和图像域中的目标图像存在一定的差异性，基本的联乘融合[26]对目标成像质量的改善很是有限。由此提出如下的指数加权联乘融合思路，即
{i_{{\text{mf}}}}(p,q) = i_1^a(p,q)*i_2^b(p,q),{\text{ }}a,b \ge 0 (22) 其中，
{i_{{\text{mf}}}}(p,q) 为融合图像中位置(p,q) 处的像素点幅度，{i_1}(p,q) ,{i_2}(p,q) 分别为回波域和图像域中目标图像在位置(p,q) 处的像素点幅度，a ,b 为对应子图像的加权指数。接下来，进一步讨论加权指数的取值。考虑到各域中目标图像的差异性，过大的加权指数会造成目标丢失，可初步限定
a,b \in [0,1] 。在此基础上，为了便于快速确定最佳加权指数取值，不妨以目标聚焦效果较好的图像域目标图像为融合主图像，通过调整回波域目标图像在联乘中的贡献来最大限度地抑制杂波和聚焦目标，即在b = 1 的情况下寻找合适的a \in [0,1] 。在这一过程中，综合考虑算法速度要求和联乘取值经验，a 的取值可精确到0.05 。最后，确定评价指标以表示不同加权指数取值对融合图像中目标成像质量的改善效果。一般来讲，目标杂波比[27] (Target to Clutter Ratio, TCR)是最常用的成像质量评价指标，其定义为
{R_{{\text{TC}}}} = 10\lg \left(\frac{{\dfrac{1}{{{n_{{\text{tg}}}}}}\displaystyle\sum\limits_{(p,q) \in {A_{{\text{tg}}}}} {{{\left| {{i_{{\text{mf}}}}(p,q)} \right|}^2}} }}{{\dfrac{1}{{{n_{{\text{cl}}}}}}\displaystyle\sum\limits_{(p,q) \in {A_{{\text{cl}}}}} {{{\left| {{i_{{\text{mf}}}}(p,q)} \right|}^2}} }}\right) (23) 其中，
{n_{{\text{tg}}}} ,{n_{{\text{cl}}}} 分别为图像中目标区域{A_{{\text{tg}}}} 和杂波区域{A_{{\text{cl}}}} 中的像素点数量。但是，该指标必须要获取目标位置等先验信息，这无疑让问题变得更加复杂。既然是对图像特性进行度量，不妨考虑像素均值这一概念，即
M = \frac{1}{{PQ}}\sum\limits_{p = 1}^P {\sum\limits_{q = 1}^Q {\left| {{i_{{\text{mf}}}}(p,q)} \right|} } (24) 随着加权指数
a 的增大，回波域目标图像在融合图像中的贡献相应增加，杂波被逐渐消除；在杂波基本消失后，继续增大加权指数，便会造成目标像素的损失，因此找到杂波消失的节点显得尤为重要。然而，在这一过程中，像素均值持续下降，且下降速度并没有出现明显的变化点。不过，受此启发，可考虑超均值像素数这一指标。相较于像素均值，超均值像素数能更加全面地描述图像中像素变化的整体情况和局部细节，尤其是杂波像素和目标像素的幅度变化过程。设图像中所有像素点的集合为
{\boldsymbol{S}} ，则该指标的定义为C = {\text{crad}}({\boldsymbol{B}}),{\text{ }}\forall (p,q) \in {\boldsymbol{B}},\left| {{i_{{\text{mf}}}}(p,q)} \right| > M (25) 其中，幅度超过像素均值的像素点的集合
{\boldsymbol{B}} \subseteq {\boldsymbol{S}} ,{\text{crad}}({\boldsymbol{B}}) 表示集合{\boldsymbol{B}} 的元素个数。同样是在加权指数增大的过程中，超均值像素数呈现出有规律的变化：从一开始，随着幅度较大的杂波像素被消除，超均值像素数以较快的速度持续减少；在杂波基本消失后，目标像素成为图像中待消除像素的“主体”，指标的减少速度会明显放缓。所以，以超均值像素数作为评价指标，并以其减速放缓的第一个点作为加权指数
a 。4.3 算法思路整合
结合上文内容，对所提杂波抑制算法的思路作进一步整合，其整体流程如图1所示。
同时，对算法的主要步骤可归纳如下：
第1阶段：利用SFAL方法得到各域中的目标图像，包括步骤1和步骤2。
步骤1 根据穿墙雷达回波信号，直接建立回波域JLRS模型；同时，以BP算法处理回波信号，得到原始图像，并以此建立图像域JLRS模型；
步骤2 利用SFAL方法依次求解各域中的JLRS模型，分别得到各域中的目标图像；
第2阶段：对各域中的目标图像进行指数加权联乘融合处理，包括步骤3和步骤4。
步骤3 给定
b = 1 ，并在a \in [0,1] 的范围内以0.05 的步长绘制超均值像素数的变化曲线，以曲线减速明显放缓的第一个点作为加权指数a ；步骤4 根据选定的加权指数
a 和b 进行指数加权联乘多域图像融合处理，得到最终的融合图像。5. 仿真验证
为了验证本文所提杂波抑制算法的速度与精度，利用MATLAB构建如下的穿墙场景：选择匀质混凝土墙体，其厚度为
0.2{\text{ m}} ，相对介电常数为5.0 ，磁导率为{{4\pi }} \times {\text{1}}{{\text{0}}^{ - 7}}{\text{ H/m}} ，电导率为5.0 \times {10^{ - 2}}{\text{ S/m}} 。天线阵列为均匀线阵，阵列孔径为2{\text{ m}} ，阵元间距为0.1{\text{ m}} ，阵元数量N = 21 ，其距离墙体的垂直距离为1.0{\text{ m}} 。天线工作在收发同置模式，各阵元依次发射信号并由自身接收回波，且各天线通道信号的采样点数T = 1024 。发射信号为单位幅度的高斯脉冲2阶导信号，脉冲形成因子为0.5{\text{ ns}} ，脉冲持续时间约为1.1{\text{ ns}} ，对应的–3 dB频谱范围为0.8～2.6 GHz。探测区域为墙后2×2.8 m的范围，共放置4个半径为3.5{\text{ cm}} 的金属小球，各金属球完全相同，在此场景的分辨率下可视作点目标处理。场景如图2所示。5.1 算法速度
在本文中，对算法速度的改进主要体现在对RPCA问题求解方法的优化。那么，不妨以图像域JLRS模型为例，验证SFAL方法的性能。
首先，在上述穿墙场景下进行目标探测，以BP算法对回波信号处理成像，成像网格设置为
256 \times 256 ，并建立JLRS模型。随后，依次以SFAL, APG, EALM, IALM方法求解该模型，得到目标图像。在参数设置方面，各方法采用默认参数，其迭代终止阈值统一设定为\zeta = {10^{ - 7}} 。在运行环境方面，所用联想台式机配置Intel Core i7-9700 3.00 GHz的8核中央处理器，其安装内存为8 GB，仿真软件为Windows 7系统下的MATLAB R2018b软件。最后，以迭代时间和目标杂波比等指标对各方法的速度及准确度进行对比，具体情况如表2所示。为了不失一般性，表中各数据均为多次独立重复试验结果的平均值。表 2 各方法性能对比算法类型 SFAL APG EALM IALM 目标杂波比(dB) 12.15 11.22 11.95 11.27 迭代次数(次) 8 163 13 39 迭代时间(s) 0.1327 1.4603 8.3925 0.5512 在准确度方面，各方法所得目标图像的目标杂波比基本持平，说明其目标成像质量相近；在速度方面，SFAL方法的迭代次数最少，迭代时间最短，相较于迭代次数较少的EALM方法及迭代时间较短的IALM方法仍有较大优势。总的来看，在准确度相当的前提下，SFAL方法以其速度优势展现出了良好的性能，而这也进一步加快了杂波抑制算法的整体速度。
5.2 算法精度
为了验证算法的精度，同样以图2的穿墙场景进行探测成像，其结果如图3所示。可以看到，其中只有墙体位置附近成像清晰，隐约可见墙后目标阴影。以本文所提算法进行杂波抑制，建立求解回波域和图像域JLRS模型，得到各自的目标图像，其结果分别如图4、图5所示。前者的墙体主杂波抑制较为彻底，但目标散焦严重，难以准确定位，且部分目标丢失；后者的目标聚焦准确，但不够清晰，且墙体位置附近有残余杂波。显然，单一域中的目标图像均未达到对目标成像质量的要求。
在进行多域联合处理之前，对加权指数
a 绘制超均值像素数的变化曲线。由图6可知，整体上该曲线呈下降趋势，初期的上升峰值可解释为杂波像素的快速消除和像素均值的急速下降所导致的超均值像素数的短暂增加，另有回波域加权指数a = 0.65 。同时，给定图像域加权指数b = 1 。根据上面得到的加权指数，进行指数加权联乘多域图像融合处理，最终成像结果如图7所示。相较于图8的背景对消成像和图9的SVD算法成像，多域联合成像中对墙体杂波尤其是残余杂波的抑制更为彻底，目标成像较为清晰，聚焦效果更好。
此外，以目标杂波比来定量分析算法的精度，相应数据如表3所示。可以看到，与SVD算法和背景对消相比，本文所提算法对原始成像中目标杂波比的提升异常明显。总的来看，本文所提算法精度较高，能够较为彻底地抑制杂波，有效提升目标杂波比。
表 3 各情况下的目标杂波比(dB)原始回波
成像多域联合
成像背景对消
成像SVD算法
成像目标杂波比 3.89 25.21 16.27 13.17 较原始成像改善 0 21.32 12.38 9.28 6. 结束语
本文针对穿墙雷达杂波抑制问题提出一种基于RPCA理论的多域联合抑制算法，该算法通过多域图像融合使目标成像准确聚焦，并利用SFAL方法来增强这一过程的实时性。经仿真证明，该算法具有良好的速度和精度，可实现对杂波的充分抑制，能够快速地为目标检测、识别等后续处理提供准确的目标信息。
实际上，本文算法仍限于BP成像和收发同置的工作模式。为了更广泛地适用于穿墙雷达领域，可考虑将所提算法推广至更一般的穿墙场景，如工作在收发分置模式下的天线阵列。而这意味着要解决一个关键问题，即杂波低秩性和目标稀疏性的增强问题。JLRS模型越标准，杂波和目标的分离就越彻底，目标成像质量就越高。同时，从降低单次迭代复杂度的角度入手来优化RPCA求解方法，可使算法的运算速度迈向新的台阶。
-
表 1 SFAL方法
输入：2维图像矩阵 {\mathbf{I}} \in {{\mathbf{R}}^{P \times Q}} ，凸函数f({\mathbf{x}}) = {\left\| {\mathbf{x}} \right\|_ * }，凸函数 g({\mathbf{x}}) = \gamma {\left\| {{\mathbf{I}} - {\mathbf{x}}} \right\|_1} ，正则化参数 \gamma = 1/\sqrt {\max (P,Q)} ； 输出：杂波分量矩阵{ {\mathbf{I} }_{\text{w} } } = { {\mathbf{x} }^{k{{ - } }1} }，目标分量矩阵{ {\mathbf{I} }_{ {\text{tg} } } } = {\mathbf{I} } - { {\mathbf{x} }^{k{{ - } }1} }。 (1) 初始化参数：\alpha = \beta = {10^{ - 6}}, {{\mathbf{x}}^0} = {{\mathbf{y}}^0} = {{\mathbf{z}}^1} = 0 , {\mu _f} = {\mu _g} = 1 , {\eta _1} = 1；k = 1。 (2) 根据式(12)和式(14)进行光滑化处理。 (3) 迭代解未收敛时执行步骤(4)到(7) (4) 根据式(18)和式(19)进行交替迭代； (5) 根据式(20)更新\eta ； (6) 根据式(21)更新{\mathbf{z}}； (7) k \leftarrow k + 1 (8) 结束循环 
下载: 导出CSV
表 2 各方法性能对比
算法类型 SFAL APG EALM IALM 目标杂波比(dB) 12.15 11.22 11.95 11.27 迭代次数(次) 8 163 13 39 迭代时间(s) 0.1327 1.4603 8.3925 0.5512
下载: 导出CSV
表 3 各情况下的目标杂波比(dB)
原始回波
成像多域联合
成像背景对消
成像SVD算法
成像目标杂波比 3.89 25.21 16.27 13.17 较原始成像改善 0 21.32 12.38 9.28
下载: 导出CSV
-
[1] 刘新, 阎焜, 杨光耀, 等. UWB-MIMO穿墙雷达三维成像与运动补偿算法研究[J]. 电子与信息学报, 2020, 42(9): 2253–2260. doi: 10.11999/JEIT190356XIN Liu, YAN Kun, YANG Guangyao, et al. Study on 3D imaging and motion compensation algorithm for UWB-MIMO through-wall radar[J]. Journal of Electronics &Information Technology, 2020, 42(9): 2253–2260. doi: 10.11999/JEIT190356 [2] SHAO Wenyi and MCCOLLOUGH T. Advances in microwave near-field imaging: Prototypes, systems, and applications[J]. IEEE Microwave Magazine, 2020, 21(5): 94–119. doi: 10.1109/MMM.2020.2971375 [3] ZHOU Yi, HUANG Chen, LIU Hongqing, et al. Front-wall clutter removal in through-the-wall radar based on weighted nuclear norm minimization[J]. IEEE Geoscience and Remote Sensing Letters, To be published. doi: 10.1109/lgrs.2020.3034568. [4] DOĞU S, AKINCI M N, ÇAYÖREN M, et al. Truncated singular value decomposition for through-the-wall microwave imaging application[J]. IET Microwaves, Antennas & Propagation, 2020, 14(4): 260–267. doi: 10.1049/iet-map.2019.0677 [5] YE Guodong, PAN Chen, DONG Youxia, et al. Image encryption and hiding algorithm based on compressive sensing and random numbers insertion[J]. Signal Processing, 2020, 172: 107563. doi: 10.1016/j.sigpro.2020.107563 [6] LEIGSNERING M, DEBES C, and ZOUBIR A M. Compressive sensing in through-the-wall radar imaging[C]. Proceedings of 2011 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Prague, Czech Republic, 2011: 4008–4011. doi: 10.1109/ICASSP.2011.5947231. [7] VAN HA T, BOUTERDOUM A, and PHUNG S L. A matrix completion approach for wall-clutter mitigation in compressive radar imaging of indoor targets[C]. 2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Calgary, Canada, 2018: 1608–1612. doi: 10.1109/ICASSP.2018.8462000. [8] TANG V H, BOUZERDOUM A, and PHUNG S L. Compressive radar imaging of stationary indoor targets with low-rank plus jointly sparse and total variation regularizations[J]. IEEE Transactions on Image Processing, 2020, 29: 4598–4613. doi: 10.1109/tip.2020.2973819 [9] CANDÈS E J, LI Xiaodong, MA Yi, et al. Robust principal component analysis?[J]. Journal of the ACM, 2011, 58(3): 1–37. doi: 10.1145/1970392.1970395 [10] TIVIVE F H C and BOUZERDOUM A. An improved SVD-based wall clutter mitigation method for through-the-wall radar imaging[C]. The 14th Workshop on Signal Processing Advances in Wireless Communications (SPAWC), Darmstadt, Germany, 2013: 430–434. doi: 10.1109/spawc.2013.6612086. [11] CANDES E J and TAO T. Near-optimal signal recovery from random projections: Universal encoding strategies?[J]. IEEE Transactions on Information Theory, 2006, 52(12): 5406–5425. doi: 10.1109/tit.2006.885507 [12] WEN Zaiwen, YIN Wotao, and ZHANG Yin. Solving a low-rank factorization model for matrix completion by a nonlinear successive over-relaxation algorithm[J]. Mathematical Programming Computation, 2012, 4(4): 333–361. doi: 10.1007/s12532-012-0044-1 [13] CHANDRASEKARAN V, SANGHAVI S, PARRILO P A, et al. Rank-sparsity incoherence for matrix decomposition[J]. SIAM Journal on Optimization, 2011, 21(2): 572–596. doi: 10.1137/090761793 [14] LIN Zhouchen, GANESH A, WRIGHT J, et al. Fast convex optimization algorithms for exact recovery of a corrupted low-rank matrix[R]. Coordinated Science Laboratory Report no. UILU-ENG-09-2214, DC-246, 2009: 1–18. [15] LIN Zhouchen, CHEN Minming, and MA Yi. The augmented Lagrange multiplier method for exact recovery of corrupted low-rank matrices[J]. arXiv preprint arXiv: 1009.5055, 2013. [16] 孙鑫. 超宽带穿墙雷达成像方法与技术研究[D]. [博士论文]. 国防科学技术大学, 2015.SUN Xin. Research on method and technique of ultra-wideband through-the-wall radar imaging[D]. [Ph. D. dissertation], National University of Defense Technology, 2015. [17] TANG V H, BOUZERDOUM A, and PHUNG S L. Multipolarization through-wall radar imaging using low-rank and jointly-sparse representations[J]. IEEE Transactions on Image Processing, 2018, 27(4): 1763–1776. doi: 10.1109/tip.2017.2786462 [18] TIVIVE F H C and BOUZERDOUM A. Joint low-rank and sparse based image reconstruction for through-the-wall radar imaging[C]. The 7th International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP), Curacao, Holland, 2017: 1–5. doi: 10.1109/camsap.2017.8313110. [19] TANG V H, BOUZERDOUM A, PHUNG S L, et al. Radar imaging of stationary indoor targets using joint low-rank and sparsity constraints[C]. 2016 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Shanghai, China, 2016: 1412–1416. doi: 10.1109/icassp.2016.7471909. [20] 韩银萍, 刘丽, 王冰洁, 等. 基于鲁棒主成分分析的混沌穿墙成像雷达杂波抑制方法[J]. 电子器件, 2020, 43(1): 142–146. doi: 10.3969/j.issn.1005-9490.2020.01.029HAN Yinping, LIU Li, WANG Bingjie, et al. Clutter removal using robust principal component analysis for chaos through-wall imaging radar[J]. Chinese Journal of Electron Devices, 2020, 43(1): 142–146. doi: 10.3969/j.issn.1005-9490.2020.01.029 [21] GOLDFARB D, MA Shiqian, and SCHEINBERG K. Fast alternating linearization methods for minimizing the sum of two convex functions[J]. Mathematical Programming, 2013, 141(1/2): 349–382. doi: 10.1007/s10107-012-0530-2 [22] NESTEROV Y. Smooth minimization of non-smooth functions[J]. Mathematical Programming, 2005, 103(1): 127–152. doi: 10.1007/s10107-004-0552-5 [23] TRAN-DINH Q. Adaptive smoothing algorithms for nonsmooth composite convex minimization[J]. Computational Optimization and Applications, 2017, 66(3): 425–451. doi: 10.1007/s10589-016-9873-6 [24] TRAN-DINH Q and CEVHER V. A primal-dual algorithmic framework for constrained convex minimization[J]. arXiv preprint arXiv: 1406.5403, 2015. [25] NESTEROV Y. Introductory Lectures on Convex Optimization: A Basic Course[M]. New York: Kluwer Academic, 2003: 45–125. [26] JIA Yong, CUI Guolong, KONG Lingjiang, et al. Multichannel and multiview imaging approach to building layout determination of through-wall radar[J]. IEEE Geoscience and Remote Sensing Letters, 2014, 11(5): 970–974. doi: 10.1109/lgrs.2013.2283778 [27] MCINTOSH B, VENKATARAMANAN S, and MAHALANOBIS A. Infrared target detection in cluttered environments by maximization of a target to clutter ratio (TCR) metric using a convolutional neural network[J]. IEEE Transactions on Aerospace and Electronic Systems, 2020, 57(1): 485–496. doi: 10.1109/TAES.2020.3024391 期刊类型引用(5)
1. 于浩，贾玮，昝继业，卞宇翔，刘金锁. 基于诱骗态的BB84协议量子秘密共享方案. 量子电子学报. 2019(03): 348-353 . 
百度学术2. CAO Dong，SONG Yaoliang，ZHU Cheng. A Novel Least-Entanglement-Assisted Asymmetric Quantum Codes Based on Sliding Grill. Chinese Journal of Electronics. 2014(03): 569-573 . 必应学术3. 王乐，邹丽，赵生妹. 一种含有安全可信任中心的量子秘密共享方案. 量子电子学报. 2014(05): 591-598 . 百度学术4. 曹东，宋耀良. 采用纠缠私钥实现多方量子隐蔽通信. 应用科学学报. 2012(01): 52-58 . 百度学术5. 袁建国，栗婵媛，黄胜，王永. 光通信中基于BIBD与循环矩阵分解的QC-LDPC码新颖构造方法. 光电子.激光. 2013(09): 1698-1701 . 百度学术其他类型引用(8)
-
下载:
百度学术
计量
- 文章访问数: 852
- HTML全文浏览量: 517
- PDF下载量: 122
- 被引次数: 13
下载:
