ISSN 
Fast Near-field Millimeter-Wave Imaging Algorithm via Generalized Approximate Message Passing
-
摘要: 传统近场毫米波均匀采样成像由于扫描时间长和计算代价大等问题无法实现实时成像。为此，该文构建了近场毫米波压缩采样成像模型及相应的观测矩阵，提出了一种基于广义近似消息传递的近场毫米波压缩采样成像快速算法。该算法将广义近似消息传递有效嵌入到期望最大化框架，加快了收敛速度；并利用快速傅里叶变换、小波滤波等方式构造了观测矩阵的快速算子，避免了大型观测矩阵的构造、存储与计算，进一步提升了算法运算速度。实验结果表明该算法可以快速、有效地从压缩采样数据中重建近场毫米波二维图像，并在重建效果与运算时间上都优于主流的快速迭代阈值收缩算法。Abstract: The conventional near-field millimeter-wave imaging from uniform sampling measurements does not work well in realizing real-time imaging due to high computational cost and long scan time. To this end, we present a near-field millimeter-wave imaging model from under-sampled measurements and the corresponding sensing matrix. Meanwhile, a fast millimeter-wave imaging algorithm based on generalized approximate message passing is proposed. The algorithm accelerates the convergence rate by embedding the generalized approximate message passing into the expectation maximization framework. In order to circumvent the construction and storage of large-scale observation matrix, we construct a fast implementation of the observation matrix using fast Fourier transform and wavelet filtering, which leads to an improved computing speed at the same time. Our experiments results confirm that the algorithm can reconstruct 2-D millimeter-wave images rapidly from under-sampled measurements, and it exhibits better performance in terms of both reconstruction quality and execution time compared with the fast iterative shrinkage-thresholding algorithm.
-
Keywords:
- AMP /
- compressed sensing /
- millimeter-wave imaging /
- nondestructive examination /
- SAR
-
近场毫米波合成孔径雷达(synthetic aperture radar, SAR)[1]成像系统可以有效检查各种复杂结构与复合材料，在无损检测与评价(nondestructive testing and evaluation, NDT&E)[2]中有着广泛应用。天线探头发射毫米波(1~10 mm)质询信号，穿透介质材料后与试验样品(SUT)的内部结构作用产生反向散射数据，成像系统利用重建算法将采集的反向散射数据重建为试验样品的高分辨率图像。
但是，传统近场毫米波SAR成像系统要获得被测物体高质量重建图像，其空间和频率采样都必须遵循奈奎斯特采样定理，要求成像系统进行均匀采样，由此带来扫描时间过长和计算量过大等问题，给实时成像带来极大困难。
压缩感知(compressed sensing, CS)理论[3]的提出给近场毫米波实时成像提供了有力工具。压缩感知将数据采样和压缩过程合二为一，对于稀疏信号或可压缩信号，只要观测矩阵满足有限等距特性(restricted isometry property, RIP)，就可以通过远低于奈奎斯特采样率的速率进行采样并精确重建出原始信号[4]。将压缩感知理论用于近场毫米波SAR成像可以有效解决扫描时间过长问题，为实时成像提供了有效解决方法。
基于压缩采样的近场毫米波成像问题成为近年来的研究热点[5-6]，许多压缩感知重建算法被应用到此问题，主要分为3类：贪婪类算法[7-8]、凸优化方法[9-10]以及稀疏贝叶斯学习方法[11-12]等。其中，稀疏贝叶斯学习算法有良好的性能，但它在每次迭代过程中都涉及到矩阵求逆运算，计算复杂度随信号维度指数量级增长。如此高的计算代价严重影响它在实际问题中的应用，尤其是高维大数据规模实际问题。
广义近似消息传递(generalized approximate message passing, GAMP) [13]是近年来新出现的消息传递框架下的贝叶斯迭代算法，由Rangan对近似消息传递算法(approximate message passing, AMP)[14]进行推广与泛化，使其适用于更一般的观测矩阵。
本文针对近场毫米波压缩采样成像问题提出了一种复数域下的基于GAMP的快速成像算法。该算法采用稀疏贝叶斯方法中的分层稀疏先验模型(hierarchical prior model)[17]，通过GAMP迭代得到待重建信号真实后验的近似，规避了稀疏贝叶斯学习算法中矩阵求逆的过程，并利用期望最大化(expectation-maximization, E-M)方法[15]更新超参数。此外，通过构造观测矩阵的快速实现方式，避免了传统贝叶斯方法中大规模观测矩阵的构造，有效缩减了内存代价并降低了算法运算时间。
1. 近场毫米波SAR成像模型
在笛卡尔坐标系内，2-D近场毫米波SAR成像系统通过探头对被测物进行扫描并采集反向散射数据，天线探头位于$(x', y', 0)$，被测物上任一点位于$(x, y, {z_0})$，定义反射场与入射场的比值为反射率函数，用$f(x, y, {z_0})$表示。天线探头采集的反向散射数据$s(x', y', \omega )$是被测物上所有点的反射叠加:
$$ s(x', y', \omega ) = \iint f(x, y, {z_0}){{\text{e}}^{ - {\text{j}}2k\sqrt {{{(x - x')}^2} + {{(y - y')}^2} + z_0^2} }}{\mkern 1mu} {\text{d}}x {\text{d}}y $$ (1) 式中，$k = \omega /c$为波数；$\omega $为工作角频率；c为光速。
指数函数表示从$(x', y', 0)$传出的球面波，可表示为一系列的平面波的叠加：
$$ {{\text{e}}^{ - {\text{j}}2k\sqrt {{{(x - x')}^2} + {{(y - y')}^2} + z_0^2} }} = \iint {{\text{e}}^{{\text{j}}{{k'}_x}(x' - x) + {\text{j}}{{k'}_y}(y' - y) + {\text{j}}{k_z}{z_0}}}{\kern 1pt} {\text{d}}{k'_x}{\text{d}}{k'_y} $$ 将它带入式(1)并去掉撇号，得到：
$$ s(x, y, \omega ) = F_{2{\text{D}}}^{ - 1}\{ {F_{2D}}\{ f(x, y, {z_0})\} {{\text{e}}^{{\text{j}}{k_z}{z_0}}}\} $$ (2) 式中，${F_{2{\text{D}}}}\{ f(x, y)\} = \iint f(x, y){{\text{e}}^{ - {\text{j}}({k_x}x + {k_y}y)}}{\text{d}}x{\text{d}}y$表示2-D傅里叶变换；$F_{2{\rm{D}}}^{ - 1}\left( \cdot \right)$表示2-D傅里叶反变换。
由电磁平面波的散射关系$k_x^2 + k_y^2 + k_z^2 = {(2k)^2}, $并结合式(2)，得到2-D SAR成像正向变换：
$$ f(x, y, {z_0}) = F_{2{\text{D}}}^{ - 1}\{ {F_{2{\text{D}}}}\{ s(x, y, \omega )\} {{\text{e}}^{ - {\text{j}}{z_0}\sqrt {4{k^2} - k_x^2 - k_y^2} }}\} $$ (3) 此正向变换将天线探头采集到的反向散射数据$s(x, y, \omega )$转换为2-D被测物体图像$f(x, y, {z_0})$。相反，2-D SAR成像反向变换将被测物体图像转换为散射数据。
为解决实时成像问题，本文主要研究近场毫米波压缩采样成像，即天线探头以欠采样方式随机采点，此过程完整观测模型可以表示为：
$$ y = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}f + n = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{\mathit{\Psi }^{\rm{T}}}x + n = \mathit{\boldsymbol{A}}x + n $$ (4) $y \in {\mathbb{C}^{\mathit{\boldsymbol{M}}}}$为压缩采样得到的随机采样数据；$f \in {\mathbb{C}^{\mathit{\boldsymbol{N}}}}$为被测物待恢复图像；$n \in {\mathbb{C}^{\mathit{\boldsymbol{N}}}}$为复高斯噪声；$\mathit{\Psi } \in {\mathbb{C}^{{\mathit{\boldsymbol{N}}} \times {\mathit{\boldsymbol{N}}}}}$为正交的小波变换矩阵；将f变换为小波域下的稀疏系数x，$\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} = {\mathit{\boldsymbol{U}}_n}F_{2{\rm{D}}}^{ - 1}\{ I[{F_{2{\rm{D}}}}( \cdot )]\} $为观测矩阵，表示随机采样过程，其中I是相位补偿项，Un是随机采样的二进制矩阵。
2. 基于GAMP的压缩成像快速算法
本压缩成像快速算法采用稀疏贝叶斯学习中的分层稀疏先验模型，通过该模型中超参数的引入，可以更好地表征信号的稀疏性，并利用E-M方法从数据中学习更新超参数，可以有效地提升图像的重建效果；另外本算法利用GAMP更新真实后验概率的近似解，规避了传统SBL方法中的矩阵求逆过程，将算法复杂度由$\mathcal{O}({N^3})$降低至$\mathcal{O}(MN)$，提升了算法运算速度；此外还通过2-D FFT、小波滤波等方式构造了式(4)中的观测矩阵A的快速实现算子，避免了大规模观测矩阵的构造、存储与运算，进一步提升了整个图像重建的运行速度(对于128×128的数据，此大矩阵维度为16 384×16 384)，从而实现了快速实时成像。
2.1 分层稀疏先验
GAMP算法要求信号先验信息已知，本算法采用一种分层稀疏先验[17]。
在第一层，假设xn服从零均值、方差为$\theta _n^{ - 1}$的高斯分布，且各个分量xn相互独立，因此有：
$$ p(x|\theta ) = \prod\limits_{n = 1}^N p ({x_n}|{\theta _n}) = \prod\limits_{n = 1}^N \mathcal{N} ({x_n}|0, \theta _n^{ - 1}) $$ (5) 式中，${\theta _n}$为控制小波系数xn稀疏性的超参数。
在第二层，给超参数${\theta _n}$赋Gamma分布(高斯分布的对偶分布)：
$$ p(\theta ) = \prod\limits_{n = 1}^N {\text{G}} {\text{amma}}({\theta _n}|a, b) = \prod\limits_{n = 1}^N {{\mathit{\Gamma} ^{ - 1}}} (a){b^a}\theta _n^{a - 1}{{\text{e}}^{ - b{\theta _n}}} $$ (6) 式中，Gamma函数为$\mathit{\Gamma }(a) = \int_{{\rm{ }}0}^{{\rm{ }}\infty } {{t^{a - 1}}} {{\rm{e}}^{ - t}}{\rm{d}}t$。
由式(5)和式(6)通过$p({x_n}) = \int p ({x_n}|{\theta _n})p({\theta _n}){\text{d}}{\theta _n}$得到x的边缘分布，文献[17]指出此积分结果为关于xn的Student-t分布，该分布以高概率在坐标原点附近取值，说明该分层先验能够很好地表征xn的稀疏性。
高斯噪声n均值为0、方差为${\sigma ^2}$，为方便估计，为噪声方差倒数$\beta = 1/{\sigma ^2}$赋予Gamma分布：$p(\beta ) = {\text{Gamma}}(\beta |c, d) = {\mathit{\Gamma} ^{ - 1}}(c){d^c}{\beta ^{c - 1}}{{\text{e}}^{ - d\beta }}$。
2.2 广义近似消息传递
令超参数集合为$\eta \buildrel \Delta \over = \{ \theta, \beta \} $，在已知先验信息和噪声前提下，假定超参数$\eta $已知，GAMP可以从观测值y中得到x后验分布$p(x|y, {\theta ^{(t)}}, {\beta ^{(t)}})$的近似估计。
根据文献[13]，GAMP主要存在两个重要的近似：
一是在输入端通过下式近似真实后验分布$p({x_n}|y, \eta )$：
$$ \hat p({x_n}|y, \mu _n^r, \phi _n^r, \eta ) = \frac{{p({x_n}|\eta )\mathcal{N}({x_n}|\mu _n^r, \phi _n^r)}}{{\int_x p ({x_n}|\eta )\mathcal{N}({x_n}|\mu _n^r, \phi _n^r)}} $$ (7) 式中，$\mu _n^r$和$\phi _n^r$是GAMP算法迭代中间变量，为了简洁表示，在推导过程中将迭代次数标号k略去。将式(5)带入式(7)，根据高斯分布乘积法则可得到$\hat p({x_n}|y, \mu _n^r, \phi _n^r$)服从高斯分布，其均值和方差分别为：$\mu _n^x = \mu _n^r/(1 + {\theta _n}\phi _n^r)$，$\phi _n^x = \phi _n^r/(1 + {\theta _n}\phi _n^r)$。
另一个近似是关于无噪声输出量${z_m} = a_m^{\text{T}}x$，其中$a_m^{\text{T}}$表示矩阵A的第m行。GAMP在输出端通过下式近似真实后验分布$p({z_m}|y, \eta )$，有：
$$ \hat p({z_m}|y, \mu _m^p, \phi _m^p, \eta ) = \frac{{p({y_m}|{z_m}, \eta )\mathcal{N}({z_m}|\mu _m^p, \phi _m^p)}}{{\int_{{\text{ }}z} p ({y_m}|{z_m}, \eta )\mathcal{N}({z_m}|\mu _m^p, \phi _m^p)}} $$ (8) 式中，$\mu _m^p$和$\phi _m^p$是GAMP算法迭代中间变量，同样在推导过程中将标号k略去。由高斯噪声假设有$p({y_m}|{z_m}, \eta ) = \mathcal{N}({y_m}|{z_m}, 1/\beta )$，再次利用高斯分布乘积法则可得到$\hat p({z_m}|y, \mu _m^p, \phi _m^p, \eta )$也服从高斯分布，其均值和方差分别为：$\mu _m^z = (\phi _m^p\beta {y_m} + \mu _m^p)/(1 + \beta \phi _m^p)$，$\phi _m^z = \phi _m^p/(1 + \beta \phi _m^p)$。
通过上述两个近似，可以推导出GAMP算法中最为关键的尺度估计函数${g_{{\rm{in}}}}( \cdot )$和${g_{{\text{out}}}}( \cdot )$。
在MMSE-GAMP形式下，根据尺度估计函数定义[13]，${g_{{\text{in}}}}( \cdot )$可以表示为：
$$ {g_{{\text{in}}}}(\mu _n^r, \phi _n^r, \eta ) = \mu _n^x = \mu _n^r/(1 + {\theta _n}\phi _n^r) $$ (9) 同时${g_{{\text{in}}}}( \cdot )$关于$\mu _n^r$的偏导数满足：
$$ \phi _n^r\frac{\partial }{{\partial \mu _n^r}}{g_{{\text{in}}}}(\mu _n^r, \phi _n^r, \eta ) = \phi _n^x = \phi _n^r/(1 + {\theta _n}\phi _n^r) $$ (10) ${g_{{\text{out}}}}( \cdot )$可以表示为：
$$ \begin{gathered} {g_{{\text{out}}}}(\mu _m^p, \phi _m^p, \eta ) = \frac{1}{{\phi _m^p}}(\mu _m^z - \mu _m^p) = \\ \frac{1}{{\phi _m^p}}\left( {\frac{{\phi _m^p\beta {y_m} + \mu _m^p}}{{1 + \beta \phi _m^p}} - \mu _m^p} \right) \\ \end{gathered} $$ (11) ${g_{{\text{out}}}}( \cdot )$关于$\mu _m^p$的偏导数满足：
$$ \phi _m^p\frac{\partial }{{\partial \mu _m^p}}{g_{{\text{out}}}}(\mu _m^p, \phi _m^p, \eta ) = \frac{{\phi _m^z - \phi _m^p}}{{\phi _m^p}} = \frac{{ - \beta \phi _m^p}}{{1 + \beta \phi _m^p}} $$ (12) 根据文献[13]，GAMP迭代完整步骤如下：
1) $\forall m \in \{ 1, 2, \cdots, M\} $:
$$ \mu _m^z(k) = \sum\limits_n {{a_{mn}}} \mu _n^x(k), \phi _m^p(k) = \sum\limits_n {a_{mn}^2} \phi _n^x(k) $$ $$ \mu _m^p(k) = \mu _m^z(k) - \phi _m^p(k)\mu _m^s(k - 1) $$ 2) $\forall m \in \{ 1, 2, \cdots, M\} $:
$$ \mu _m^s(k) = {g_{{\text{out}}}}(\mu _m^p(k), \phi _m^p(k), {\eta ^{(t)}}) $$ $$ \phi _m^s(k) = - \frac{\partial }{{\partial \mu _m^p}}{g_{{\text{out}}}}(\mu _m^p(k), \phi _m^p(k), {\eta ^{(t)}}) $$ 3) $\forall n \in \{ 1, 2, \cdots, N\} $:
$$ \phi _n^r(k) = {(\sum\limits_m {a_{mn}^2} \phi _m^s(k))^{ - 1}} $$ $$ \mu _n^r(k) = \mu _n^x(k) + \phi _n^r(k)\sum\limits_m {{a_{mn}}} \mu _m^s(k) $$ 4) $\forall n \in \{ 1, 2, \cdots, N\} $:
$$ \hat x(k + 1){\text{ = }}\mu _n^x(k + 1) = {g_{{\text{in}}}}(\mu _n^r(k), \phi _n^r(k), {\eta ^{(t)}}) $$ $$ \phi _n^x(k + 1) = \phi _n^r(k)\frac{\partial }{{\partial \mu _n^r}}{g_{{\text{in}}}}(\mu _n^r(k), \phi _n^r(k), {\eta ^{(t)}}) $$ (13) 2.3 超参数估计
期望最大化(E-M)算法可以实现在估计信号x的同时学习超参数$\eta $，主要分为E-step和M-step两个步骤[16]：
E-step：计算$Q(\eta |{\eta ^{(t)}}) = {{\text{E}}_{x|y, {\eta ^{(t)}}}}[{\text{log}}p(x, y|\eta )]$,
M-step：更新${\eta ^{(t + 1)}} = \mathop {\arg {\text{max}}}\limits_\eta {\mkern 1mu} Q(\eta |{\eta ^{(t)}})$。
利用E-M算法进行超参数估计，将x看作隐变量迭代求解Q函数，即：
$$ \begin{gathered} {\eta ^{(t + 1)}} = \mathop {\arg {\text{max}}}\limits_\eta Q(\eta |{\eta ^{(t)}}) = \\ \mathop {\arg {\text{max}}}\limits_\eta {{\text{E}}_{x|y, {\eta ^{(t)}}}}[{\text{log}}p(x, y|\eta )] \\ \end{gathered} $$ 首先对超参数$\theta $进行估计，Q函数关于${\theta _n}$求偏导数有：
$$ \frac{\partial }{{\partial {\theta _n}}}Q(\eta |{\eta ^{(t)}}) = \frac{\partial }{{\partial {\theta _n}}}{{\text{E}}_{x|y, {\eta ^{(t)}}}}[{\text{log}}p({x_n}|{\theta _n})p({\theta _n};a, b)] = \\\frac{1}{{2{\theta _n}}} - \frac{{\langle x_n^2\rangle }}{2} + \frac{{a - 1}}{{{\theta _n}}} - b $$ (14) 式中，$\langle \cdot \rangle $表示关于概率$\hat p({x_n}|y, \mu _n^r(k), \phi _n^r(k), {\eta ^{(t)}})$求期望，此概率为GAMP算法中关于$p(x|y, {\eta ^{(t)}})$的近似。于是有$\langle x_n^2\rangle = \frac{{{{(\mu _n^r(k))}^2}}}{{{{(1 + \theta _n^t\phi _n^r(k))}^2}}} + \frac{{\phi _n^r(k)}}{{1 + \theta _n^t\phi _n^r(k)}}$，令式(14)等于0，得到${\theta _n}$的估计表达式：
$$ \theta _n^{(t + 1)} = (2a - 1)/(2b + \langle x_n^2\rangle )\;\;\;\forall n \in 1, 2, \cdots , N $$ (15) 接着估计超参数$\beta $，Q函数关于$\beta $求偏导数有：
$$ \frac{\partial }{{\partial \beta }}Q(\eta |{\eta ^{(t)}}){\text{ = }}\frac{\partial }{{\partial \beta }}{{\text{E}}_{z|y, {\eta ^{(t)}}}}[{\text{log}}p(y|z, \beta )p(\beta ;c, d)] = \\\frac{M}{{2\beta }} - \frac{1}{2}\sum\limits_{m = 1}^M {\langle (} {y_m} - {z_m}{)^2}\rangle + \frac{{c - 1}}{\beta } - d $$ (16) 式中，$\langle \cdot \rangle $表示关于概率$\hat p({z_m}|y, \mu _m^p(k), \phi _m^p(k), {\eta ^{(t)}})$求期望，此概率为GAMP算法中关于$p({z_m}|y, {\eta ^{(t)}})$的近似。于是有$\langle {({y_m} - {z_m})^2}\rangle = {({y_m} - \mu _m^z)^2} +\phi _m^z $，令式(16)偏导数等于0，得到$\beta $的估计表达式：
$$ {\beta ^{(t + 1)}} = \frac{{M + 2c - 2}}{{2d + \sum\limits_m {\langle (} {y_m} - {z_m}{)^2}\rangle }} $$ (17) 2.4 算法总结
本文提出的近场毫米波压缩采样成像快速算法步骤总结如下：
1) 初始化：
k=0，初始化$\mu _m^s( - 1) = 0$，$\forall m \in \{ 1, 2, \cdots, M\} $；${\theta ^{(0)}}$和${\beta ^{(0)}}$；$\{ \mu _n^x(k)\} _{n = 1}^N$和$\{ \phi _n^x(k)\} _{n = 1}^N$；
2) GAMP迭代：
$t \geqslant 0$，给定$\theta _n^{(t)}$和$\beta _n^{(t)}$，进行式(13)中的GAMP迭代；
3) E-M超参数更新：
利用GAMP迭代得到的结果和中间量按照式(15)、式(17)中的规则更新超参数$\theta _n^{(t + 1)}$和$\beta _n^{(t + 1)}$；
重复步骤2)、3)直到${\sum\limits_n {\left| {\hat x\left( {k + 1} \right) - \hat x\left( k \right)} \right|} ^2} \le \epsilon, \epsilon $为预设定的误差容限，输出$\hat x({k_t})$(kt为最后一次迭代索引)，再通过小波逆变换${\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}^{\rm{T}}}$将重建小波系数$\hat x({k_t})$变换为重建图像$\hat f$。
3. 实验
3.1 实验方法
为验证本文所提出的基于广义近似消息传递的近场毫米波压缩成像快速算法(以下简称GAMP)的有效性，以金属物检测为实例，在不同采样率下进行二维图像重建实验，并将结果与原始的AMP[14]和主流的快速迭代阈值收缩算法[18](fast iterative shrinkage-threshtlding algorithm, FISTA)的重建结果进行对比。
检测实物如图 1所示，从右上至左下依次为小螺帽、垫圈和滚珠。
此近场毫米波成像系统工作在150 GHz。天线探头由开口的矩形波导构成，其位于被测物上方40 mm处，在X方向和Y方向以0.5 mm的步进对物体进行光栅扫描，扫描面积为64 mm×64 mm，通过采集反向散射数据得到大小为128×128的全采样数据，实验中通过对全采样数据随机采样得到不同采样率$\sigma $(10%、20%、30%、40%、50%)下的欠采样数据，之后利用重建算法将其重建为二维图像。
本实验通过欠采样数据重构图像与全采样数据重构图像两者之间的结构相似性(structural similarity, SSIM)来评价算法重建质量[19]，两图像f1、f2的SSIM值定义如下：
$$ {\text{SSIM}}({f_1}, {f_2}) = \frac{{(2{\mu _1}{\mu _2} + {c_1})(2{\sigma _{12}} + {c_2})}}{{(\mu _1^2{\text{ + }}\mu _2^2 + {c_1})(\sigma _1^2 + \sigma _2^2 + {c_2})}} $$ 式中，${\mu _1}$、${\mu _2}$分别为图像${f_1}$、${f_2}$的平均值；${\sigma _1}$、${\sigma _2}$分别为图像${f_1}$、${f_2}$的方差；${\sigma _{12}}$为图像${f_1}$、${f_2}$的协方差。由定义可知两相同图像的SSIM值为1，因此若SSIM值越大则算法重构效果越好。
实验全部算法均在MATLAB R2015b上运行实现，计算机的CPU为3.6 GHz的Intel(R) Core(TM) i7-4790，内存为8 GB。
3.2 实验结果及分析
在本文提出的基于GAMP的快速算法中参数设定为：$a = c = 1$，$b = d = {10^{ - 6}}$。在对比实验中，为保证对比有效，对于3种重构算法，采用相同的小波基(8阶消失矩的Symmlet)，附加相同的零均值复高斯噪声(输入信噪比为60 dB)，采用相同的break判别条件(误差容限$\epsilon$为10-4)。
第一个实验为利用本文提出的快速算法对不同采样率下的欠采样数据进行重构，不同采样率下的重建图像如图 2所示，其中图 2a为全采样数据重建图像，将其作为基准来评价欠采样数据重建图像；图 2b为20%采样率的重建图像，SSIM为0.602 0，可以大致分辨出3个物体，但图像较为模糊且存在很多阴影；图 2c为30%采样率的重建图像，SSIM为0.716 8，图像较清晰，阴影减少，可以较好地分辨出3个物体；图 2d为40%采样率的重建图像，SSIM为0.798 6，图像质量进一步提升，阴影进一步减少，可以很好地分辨出3个物体。由此验证了本算法可以有效地从欠采样数据中重建图像。
第二个实验为3种算法的重建结果(SSIM值、PSNR值与运行时间Time)对比，3种算法在不同采样率下分别进行100次重复试验，得到的平均结构相似性值SSIM、峰值信噪比PSNR和平均运行时间Time如表 1所示，由表中数据对比分析可知不同采样率下，本文提出的基于GAMP的快速成像算法在运算速度方面与原始的AMP算法基本一致，优于目前主流的FISTA；在成像质量上略优于AMP和FISTA，进一步验证了本算法的有效性。
表 1 3种重构算法不同采样率下的重建结果对比算法 采样率 10% 20% 30% 40% 50% SSIM PSNR Time/s SSIM PSNR Time/s SSIM PSNR Time/s SSIM PSNR Time/s SSIM PSNR Time/s GAMP 0.480 8 20.58 2.86 0.602 1 23.08 2.75 0.716 7 25.17 2.49 0.798 5 26.83 2.22 0.864 5 28.67 2.08 AMP 0.457 6 19.36 2.55 0.565 1 21.77 2.76 0.679 6 23.65 2.40 0.761 1 25.14 2.31 0.835 8 26.69 2.28 FISTA 0.450 6 19.28 6.22 0.561 2 21.67 5.56 0.683 9 23.56 4.11 0.769 1 24.95 3.63 0.845 0 26.58 3.12 4. 结束语
针对近场毫米波成像系统中的高计算量与实时成像问题，本文提出了一种基于GAMP的快速成像算法。该算法通过在期望最大化框架中引入GAMP，获得较低的算法复杂度，并通过构建观测矩阵A的快速实现算子进一步加快算法运行速度。实验结果表明该算法可以从欠采样数据中快速、精确重建近场毫米波二维图像，并且在重建质量与运算速度上都优于主流的FISTA。
-
表 1 3种重构算法不同采样率下的重建结果对比
算法 采样率 10% 20% 30% 40% 50% SSIM PSNR Time/s SSIM PSNR Time/s SSIM PSNR Time/s SSIM PSNR Time/s SSIM PSNR Time/s GAMP 0.480 8 20.58 2.86 0.602 1 23.08 2.75 0.716 7 25.17 2.49 0.798 5 26.83 2.22 0.864 5 28.67 2.08 AMP 0.457 6 19.36 2.55 0.565 1 21.77 2.76 0.679 6 23.65 2.40 0.761 1 25.14 2.31 0.835 8 26.69 2.28 FISTA 0.450 6 19.28 6.22 0.561 2 21.67 5.56 0.683 9 23.56 4.11 0.769 1 24.95 3.63 0.845 0 26.58 3.12 
下载: 导出CSV
-
[1] 皮亦鸣, 杨建宇, 付毓生, 等.合成孔径雷达成像原理[M].成都:电子科技大学出版社, 2007. PI Yi-ming, YANG Jian-yu, FU Yu-sheng, et al. Synthetic aperture radar imaging principle[M]. Chengdu:University of Electronic Science and Technology of China Press, 2007.
[2] KHARKOVSKY S, ZOUGHI R. Microwave and millimeter wave nondestructive testing and evaluation - overview and recent advances[J]. Instrumentation & Measurement Magazine IEEE, 2017, 10(2):26-38. http://cn.bing.com/academic/profile?id=5de2dc58b12addb71bae0791c9b1b0f0&encoded=0&v=paper_preview&mkt=zh-cn
[3] DONOHO D L. Compressed sensing[J]. IEEE Transactions on Information Theory, 2006, 52(4):1289-1306. DOI: 10.1109/TIT.2006.871582
[4] CANDES E J. The restricted isometry property and its implications for compressed sensing[J]. Comptes Rendus Mathematique, 2008, 346(9-10):589-592. DOI: 10.1016/j.crma.2008.03.014
[5] BI Dong-jie, XIE Yong-le, MA Lan, et al. Multifrequency compressed sensing for 2-D near-field synthetic aperture radar image reconstruction[J]. IEEE Transactions on Instrumentation & Measurement, 2017, 66(4):777-791. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=3a74264e67d261dcdd4d13dc4025bb58
[6] ZHANG X, HUANG P, LI X. 2D SAR imaging scheme based on compressive sensing[J]. Electronics Letters, 2014, 50(2):114-116. DOI: 10.1049/el.2013.3253
[7] SOUSSEN C, GRIBONVAL R, IDIER J, et al. Joint K-step analysis of orthogonal matching pursuit and orthogonal least squares[J]. IEEE Transactions on Information Theory, 2013, 59(5):3158-3174. DOI: 10.1109/TIT.2013.2238606
[8] NEEDELL D, TROPP J A. CoSaMP:Iterative signal recovery from incomplete and inaccurate samples[J]. Applied & Computational Harmonic Analysis, 2009, 26(3):301-321. http://d.old.wanfangdata.com.cn/Periodical/jsjgc201306005
[9] ZHANG Z, RAO B D. Iterative reweighted algorithms for sparse signal recovery with temporally correlated source vectors[C]//IEEE International Conference on Acoustics, Speech and Signal Processing.[S.l.]: IEEE, 2011: 3932-3935.
[10] SADEGHI M, BABAIE Z M. Iterative sparsification-projection:Fast and robust sparse signal approximation[J]. IEEE Transactions on Signal Processing, 2016, 64(21):5536-5548. DOI: 10.1109/TSP.2016.2585123
[11] JI S, XUE Y, CARIN L. Bayesian compressive sensing[J]. IEEE Transactions on Signal Processing, 2008, 56(6):2346-2356. DOI: 10.1109/TSP.2007.914345
[12] 叶娅兰, 何文文, 程云飞, 等.面向压缩感知的基于相关性字典学习算法[J].电子科技大学学报, 2017, 46(5):703-708. DOI: 10.3969/j.issn.1001-0548.2017.05.011 YE Ya-lan, HE Wen-wen, CHENG Yun-fei, et al. Correlation based dictionary learning algorithm for compressed sensing[J]. Journal of University of Electronic Science and Technology of China, 2017, 46(5):703-708. DOI: 10.3969/j.issn.1001-0548.2017.05.011
[13] RANGAN S. Generalized approximate message passing for estimation with random linear mixing[C]//IEEE International Symposium on Information Theory Proceedings.[S.l.]: IEEE, 2011: 2168-2172.
[14] DONOHO D L, MALEKI A, MONTANARI A. Messagepassing algorithms for compressed sensing[J]. Proceedings of the National Academy of Sciences of the United States of America, 2009, 106(45):18914. DOI: 10.1073/pnas.0909892106
[15] DEMPSTER A P, LAIRD N M, RUBIN D B. Maximum likelihood from incomplete data via the EM algorithm[J]. Journal of the Royal Statistical Society, 1977, 39(1):1-38. http://cn.bing.com/academic/profile?id=09c6ef77a32f13f8c2d65f1f902e1d29&encoded=0&v=paper_preview&mkt=zh-cn
[16] SUN Y, BABU P, PALOMAR D P. Majorization-minimization algorithms in signal processing, communications, and machine learning[J]. IEEE Transactions on Signal Processing, 2016, 65(3):794-816. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=7f9896f7e524ade132ab15e3addc018e
[17] TIPPING M E. Sparse bayesian learning and the relevance vector machine[J]. Journal of Machine Learning Research, 2001, 1(3):211-244. http://d.old.wanfangdata.com.cn/Periodical/dsjs201802008
[18] BECK A, TEBOULLE M. A fast iterative shrinkage-thresholding algorithm for linear inverse problems[J]. Siam Journal on Imaging Sciences, 2009, 2(1):183-202. DOI: 10.1137/080716542
[19] WANG Z, BOVIK A C, SHEIKH H R, et al. Image quality assessment:from error visibility to structural similarity[J]. IEEE Transactions on Image Processing, 2004, 13(4):600-612. DOI: 10.1089-fpd.2009.0394/
-
期刊类型引用(2)
1. 李海，谢瑞杰，谢伶莉，孟凡旺. 复杂地形环境下基于GAMP-STAP的低空风切变风速估计方法. 电子与信息学报. 2023(02): 576-584 . 
百度学术
2. 帅萍，高乐，吕珏，毕东杰，李西峰，谢永乐. 毫米波近场合成孔径二维成像的快速算法研究. 科技与创新. 2021(18): 33-36 . 百度学术
其他类型引用(2)
计量
- 文章访问数: 6571
- HTML全文浏览量: 1863
- PDF下载量: 130
- 被引次数: 4