1.   引言

压缩感知^[1,2]是近些年来兴起的一种新的信号处理技术，可以应用在压缩成像和远场诊断等众多领域^[3,4]。该理论采用非自适应线性投影来保持原始信号的结构，然后以远低于奈奎斯特频率的采样频率对信号进行采样，最后通过数值优化问题准确重构出原始信号。测量矩阵设计是压缩感知理论的核心问题之一，测量矩阵分为随机测量矩阵和确定性测量矩阵。已有研究表明，随机测量矩阵是一类性能优良的测量矩阵，可以很大概率地恢复信号。然而随机测量矩阵在实际应用时所需存储量巨大，硬件难以实现，并且随机测量矩阵不能保证每次生成的矩阵都具备较好的性能。为解决这些实际问题，学者提出了确定性测量矩阵。目前国内外学者提出了许多构造确定性测量矩阵的方法，如文献[5]基于有限域得到一个基矩阵，通过嵌入矩阵得到一类稀疏测量矩阵；文献[6]基于Kasami码提出一种构造测量矩阵的方法；文献[7]利用有限域多项式的相关性构造确定性测量矩阵；文献[8]使用平衡不完全块设计测量矩阵。然而这些确定性测量矩阵维度较为固定，构造具有更加灵活维度的测量矩阵成为研究的难题。

卷积压缩感知是近年来提出的一类新型压缩感知技术。卷积压缩感知测量矩阵由于是从一个确定性测量矩阵中随机抽取若干行构成的，因此具有矩阵维度灵活的优点，这类矩阵也称为半确定性测量矩阵。目前卷积压缩感知测量矩阵的相关成果并不多，文献[9]通过抽取Legendre序列提出一种卷积压缩感知测量矩阵的构造方法；文献[10]提出一种基于抽样Sidelnikov序列的卷积压缩感知测量矩阵构造方法；Li等人^[11]通过扩展Frank-Zadoff-Chu(FZC)序列以及扩展Golay序列构造卷积压缩感知测量矩阵，但其相关性参数较大。为解决这一问题，进一步优化卷积压缩感知测量矩阵，本文将分圆类理论引入压缩感知测量矩阵设计中，构造的测量矩阵相关性比文献[11]低。仿真实验表明本文得到的矩阵不仅可以用于重构1维信号，而且可以用于重构2维图像以及信道估计。

2.   基本概念

2.1   卷积压缩感知

自然界存在的信号一般不是绝对稀疏的，信号$x$通过变换域成为具有稀疏性的信号，压缩感知中信号的压缩采样过程可以表述为

$x \in {R^N}$是长度为$N$的原始信号，${{\varPhi}}$是大小为$M \times N$的测量矩阵，${{\varPsi}}$是大小为$N \times N$的稀疏基，$y \in {R^M}$为压缩采样后得到的信号。

卷积压缩感知的核心思想是对信号进行卷积和滤波，然后对输出的信号进行2次采样，最终完成信息的获取。Tropp^[12]最早提出了卷积压缩感知的思想。卷积压缩感知理论经过众多研究者的研究与实验证明，提出了一种使用N-tap随机滤波的循环卷积理论^[13]，可以用式(2)表示部分循环矩阵

${R_\Omega }$是一个采样算子，$N \times N$的方阵${{A}}$通过${{a}} = [{a_0}\;{a_1}\; ···· \;{a_{N - 1}}]^{\rm{T}}$循环移位构成。矩阵${{A}}$可以利用快速傅里叶变换(Fast Fourier Transform, FFT)进行对角化，如式(3)所示

${{F}}$是一个$N \times N$的快速傅里叶变换矩阵，${{F}}$的第$i$行和第$j$列元素${F_{i,j}} = \omega _N^{ij}$, ${\omega _N} = \exp ({{{\rm{2}}\pi \sqrt {{\rm{ - 1}}} }/ N})$, $0 \le i \le N - 1$, ${\rm{0}} \le j \le N - 1$, “${( \cdot )^{\rm{T}}}$”表示矩阵的转置，${{V}} = {\rm{diag}}({{\sigma}} )$表示一个$N \times N$的对角矩阵，${{\sigma}} = {[{\sigma _0}\;{\sigma _1}\; ···· \;{\sigma _{N - 1}}]^{\rm{T}}}$。式(3)表明循环运算是可以快速计算的，${{a}}$可以通过取对角线序列${{\sigma}}$的快速傅里叶逆变换(Inverse Fast Fourier Transform, IFFT)得到，如式(4)所示

一般来说，对于傅里叶光学以及编码孔径成像等领域要求其测量矩阵是全实数矩阵^[11]。为了使矩阵${{A}}$是一个全实数矩阵，对角线序列${{\sigma}}$需要满足共轭对称性条件^[14]：${{{\sigma}} _k} = {{\sigma}} _{N - k}^{\rm{H}}$, $1 \le k \le N - 1$, “${( \cdot )^{\rm{H}}}$”表示序列的共轭对称性；同时为保证测量矩阵的良好性质，序列${{a}}$需要满足条件：$\mu ({{a}})\! =\!\! \mathop {\max}\limits_{0 \le t \le N - 1} \left| {{a_t}} \right| \le \delta$。

2.2   限制等距性

在压缩感知理论中，限制等距性(Restricted Isometry Property, RIP)^[15]是判断测量矩阵性能优劣的重要指标。为了重构稀疏信号，Candes等人^[15]证明了测量矩阵${{\varPhi}}$必须满足限制等距性条件。

定义1　对于给定的测量矩阵${{\varPhi}}$和所有$k$稀疏信号$x$，找到最小的常数${\delta _k}$使之满足

如果所有$k$阶稀疏信号$x$均满足式(5)的最小常数${\delta _k}$，那么称测量矩阵${{\varPhi}}$具备$k$阶RIP性质，称${\delta _k}$为$k$阶限制等距常量。

定义2^[11]　对于$M \times N$的矩阵${{A}}$, ${A_{i,j}}$($0 \le i \le M - 1$, $0 \le j \le N - 1$)代表矩阵${{A}}$的第$i$行和第$j$列元素。对于$N \times N$的矩阵${{A}}$, ${{A}}$的最大相关性参数$\mu ({{A}})$定义为

文献[16]推导出RIP性质与相关性之间的关系为：${\delta _k} \le (k - 1)\mu ({{A}})$，因此可知当矩阵相关性足够小时，矩阵满足限制等距常数的RIP性质。一般来说，测定一个测量矩阵的RIP特性是非常困难的，测量矩阵通常采用容易计算的相关性作为理论支撑。可以通过构造具有较低相关性的测量矩阵进而保证其良好的RIP，因此相关性成为确定性测量矩阵优劣的主要评判标准。

2.3   分圆类

定义3^[17]　设$p = ef + 1$, $p$为奇素数，$e$为正整数，${F_p}$为$p$阶有限域，$F_p^* = {F_p}\backslash \{ 0\}$。设$\alpha$为${F_p}$的本原元，则$F_p^* = \{ {\alpha ^0},{\alpha ^1}, ···· ,{\alpha ^{p - 2}}\}$。定义子集

称$D_\lambda ^{(e,p)}$为有限域${F_p}$的$e$阶分圆类，则有${F}_{p}^{*}= {D}_{0}^{(e,p)}\cup {D}_{1}^{(e,p)}\cup ···· \cup {D}_{e-1}^{(e,p)}$。

相应的$e$阶分圆数定义为$(\lambda ,j) \,=\, \left| ({D_\lambda } \,+\, 1) \cap {D_j} \right|$, ${\rm{0}} \le \lambda ,j \le e - 1$。

3.   卷积压缩感知测量矩阵的构造

本节利用分圆类构造卷积压缩感知测量矩阵，所得到的测量矩阵元素均为实数，适用于2维图像的重构。与文献[10]和文献[11]的方法相比，本节得到的矩阵具有更优的RIP特性。仿真实验表明，新矩阵可以应用于重构1维信号、2维图像以及信道估计。

3.1   构造方法

3.1.1   基于2阶分圆类的测量矩阵

设$p = 2f + 1$, $f$为偶数。基于2阶分圆类构造序列${{s}} = [{s_0}\;{s_1}\; ··· \;{s_{p - 1}}]$, ${s_t}(0 \le t \le p - 1)$为

定理1　构造的序列$s$满足共轭对称性和低相关性，即(1) ${s_t} = {s_{p - t}}$, $1 \le t \le p - 1$; (2) $\left| {{{\hat s}_t}} \right| \le 1 + {{\rm{1}} /{\sqrt p }}$, $\hat {{s}}$表示序列 ${{s}}$的$p$点傅里叶逆变换。

证明　证明满足共轭对称性：$p = 2f + 1$，因为$f$为偶数，设$f = 2r$, $\alpha$为${F_p}$的本原元，则${\alpha ^{\frac{{p - 1}}{2}}} = {\rm{ - }}1$即$- 1 = {\alpha ^{2r}} \in D_0^{(2,p)}$。假设$t \in D_\lambda ^{(2,p)}$，那么$t$表示为$t = {\alpha ^{2m + \lambda }}$，则${\rm{ - }}t = {\alpha ^{2r}}{\alpha ^{2m + \lambda }} = {\alpha ^{2(r + m) + \lambda }} \in D_\lambda ^{(2,p)}$, $t$和${\rm{ - }}t$都属于集合$D_\lambda ^{(2,p)}$，由于都是模$p$运算，所以$p - t \in D_\lambda ^{(2,p)}$，即${s_t} = {s_{p - t}}$, $1 \le t \le p - 1$。

证明相关性：$p = {\rm{2}}f{\rm{ + 1}} = {\rm{4}}r{\rm{ + 1}} \equiv {\rm{1mod4}}$, 2阶分圆数为

${{s}} = [{s_0}\;{s_1}\; ··· \;{s_{N - 1}}]$, $N = p = 2f + 1$, $\hat {{s}}$表示序列${{s}}$的$N$点IFFT，即

(1) 当$k = 0$时，有

(2) 当$k \!\ne\! {\rm{0}}$时，令${s_t} \!=\! {\rm{1}} \!-\! 2{b_t}$，即${s_t} \!=\! {{\rm{( - 1)}}^{^{{b_t}}}}$，则有

即$\left| {{{\hat s}_k}} \right| = \dfrac{{\rm{2}}}{{\sqrt N }}\left| {{{\hat b}_k}} \right|$, $k \ne 0$。

推导出${\hat b_k} = \dfrac{{{\rm{ - }}1 \pm \sqrt N }}{2}$，进而可得$\mathop {\max }\limits_{0 \le k \le N - 1} \left| {{{\hat s}_k}} \right| = \dfrac{2}{{\sqrt N }}\dfrac{{\sqrt N + 1}}{2} = 1 + \dfrac{1}{{\sqrt N }}$。证毕

根据定理1和定义2可以得到定理2。

定理2　将序列${{s}}$作为对角向量${{\sigma}}$，按照式(4)得到向量${{a}}$，对${{a}}$进行循环移位得到矩阵${{A}}$，则${{A}}$的相关性参数为$\mu ({{A}}) \le 1 + {1 / {\sqrt N }}$。

3.1.2   基于$e$$(e > 2)$阶分圆类的测量矩阵

设$p = ef + 1$, $e$和$f$均为偶数，$p$为奇素数。基于$e$阶分圆类构造序列${{s}} = [{s_0}\;{s_1}\; ··· \;{s_{p - 1}}]$, ${s_t}(0 \le t \le p - 1)$为

定理3　构造的序列$s$满足共轭对称性和低相关性，即：(1)${s_t} = {s_{p - t}}$, $1 \le t \le p - 1$; (2)$\left| {{{\hat s}_t}} \right| \le {\rm{2}}$, $\hat {{s}}$表示序列${{s}}$的$p$点IFFT。

证明　证明满足共轭对称性：$p = ef + 1$，因为$f$为偶数，设$f = 2r$, $\alpha$为${F_p}$的本原元，则${\alpha ^{\frac{{p - 1}}{{\rm{2}}}}} = {\rm{ - }}1$即$- 1 = {\alpha ^{er}} \in D_0^{(e,p)}$。假设$t \in D_\lambda ^{(e,p)}$，那么$t$表示为$t = {\alpha ^{em + \lambda }}$，则${\rm{ - }}t = {\alpha ^{er}}{\alpha ^{em + \lambda }} = {\alpha ^{e(r + m){\rm{ + }}\lambda }} \in D_\lambda ^{(e,p)}$, $t$和${\rm{ - }}t$都属于集合$D_\lambda ^{(e,p)}$，由于都是模$p$运算，所以$p - t \in D_\lambda ^{(e,p)}$，即${s_t} = {s_{p - t}}$, $1 \le t \le p - 1$。

当$e > 2$时，求解所得的卷积压缩感知测量矩阵的相关性较为困难，然而通过计算机分析发现，测量矩阵的相关性满足低相关性条件。图1给出基于$e$阶分圆类构造矩阵的相关性$\mu ({{A}})$的分布图。

图 1  相关性分布图

下载: 全尺寸图片 幻灯片

图1分别给出了基于2阶、4阶、6阶、8阶分圆类构造矩阵的相关性$\mu ({{A}})$的分布图，截取长度$N$为$1 \le N \le 1600$，从图中可以看出基于分圆类构造矩阵的相关性参数$\mu ({{A}}) \le 2$。证毕

根据定理3和定义2可以得到定理4。

定理4　将序列${{s}}$作为对角向量${{\sigma}}$，按照式(4)得到向量${{a}}$，对${{a}}$进行循环移位得到矩阵${{A}}$，则${{A}}$的相关性参数为$\mu ({{A}}) \le 2$。

表1列举了已有序列构造测量矩阵的不同方法，并对其进行了对比。本文基于$e$阶分圆类构造的测量矩阵比已有序列构造的矩阵相关性参数更加接近理论值，因此具有更好的RIP性质。

表 1  与已有序列构造的矩阵相关性比较

	对角向量$\sigma$	序列长度$N$	相关性参数$\mu ({{A}})$
文献[10]	抽样Sidelnikov序列	$N = \dfrac{ { {p^m} - 1} }{c}$, $c$为偶数	$\sqrt {c + {{\rm{1}} / N}} + {1 / {\sqrt N }}$, $N$为偶数
			$\sqrt {c + {1 / N}}$, $N$为奇数
文献[11]	Extended Frank-Zadoff-Chu(扩展FZC)序列	$N$为偶数	$4 + {4 / {\sqrt N }}$
		$N$为奇数	$2.69 + {{8.15} / {\sqrt N }}$
	Extended Golay(扩展Golay)序列	$N = {2^{{k_1}}}{10^{{k_2}}}{26^{{k_3}}}$, $N$为偶数，${k_1},{k_2},{k_3}$为整数	$2 + {2 / {\sqrt N }}$
		$N = {2^{{k_1}}}{10^{{k_2}}}{26^{{k_3}}} \pm 1$, $N$为奇数，${k_1},{k_2},{k_3}$为整数	$2 + {1 / {\sqrt N }}$
本文	由2阶分圆类得到的序列	$N = p$为奇素数，$p \equiv 1({\rm{mod}}4)$	$1 + {1 / {\sqrt N }}$
	由$e > 2$阶分圆类得到的序列	$N = p$为奇素数	$\mu ({{A} }) \le 2$

下载: 导出CSV 

| 显示表格

3.2   仿真实验分析

3.2.1   1维信号仿真

在模拟中使用的$k$稀疏信号的非零坐标位服从高斯分布，本文模拟在无噪声条件和有噪声条件下信号$x$的恢复性能，使用正交匹配跟踪(Orthogonal Matching Pursuit, OMP)算法^[18]作为恢复算法，若重构信号估计值$\hat x$的平方误差相当小，即满足${\Vert \widehat{x}-x\Vert }_{2}<{10}^{-6}$，则认为原始信号恢复成功。1维信号的仿真结果如图2、图3、图4、图5所示。

图 2  不同稀疏度下的重构百分比

下载: 全尺寸图片 幻灯片

图 3  不同稀疏度下的重构百分比

下载: 全尺寸图片 幻灯片

图 4  不同稀疏度下的输出信噪比

下载: 全尺寸图片 幻灯片

图 5  不同稀疏度下的输出信噪比

下载: 全尺寸图片 幻灯片

图2和图4中$e=2, f=128,\, {R}_{\Omega }=32$，生成大小为32×257的测量矩阵，图3和图5中$e{=4},\, f{=130}, {R}_{\Omega }{=64}$，生成大小为64×521的测量矩阵。图2和图3在无噪声条件下基于本文提出的测量矩阵恢复信号的性能优于随机高斯矩阵以及文献[11]构造的测量矩阵；图4和图5在被采样信号的信噪比(Signal-to-Noise, SNR)为30 ${\rm{dB}}$的条件下给出了不同稀疏度下的输出信噪比。仿真结果表明，信号在有噪声的情况下，本文提出的测量矩阵恢复信号的性能优于同等条件下随机高斯矩阵以及文献[11]构造的测量矩阵。

3.2.2   2维图像仿真

本文以大小为256×256的Lena图像作为原始的2维图像信号，以小波变换基对信号进行稀疏变换，使用Tvl1-l2^[19]算法作为恢复算法。恢复图像质量的标准采用峰值信噪比(Peak Signal to Noise Ratio, PSNR)来进行评价。2维图像的仿真结果如图6所示。

图 6  不同测量矩阵重构2维图像

下载: 全尺寸图片 幻灯片

图6中$e=2, \,f{=260},\, {R}_{\Omega }{=256}$，随机抽取256列生成大小为256×256的测量矩阵，图中显示不同测量矩阵对应的PSNR分别为：(a)原始图像；(b)随机高斯测量矩阵：13.0 dB；(c)扩展Golay序列构造的测量矩阵：13.0 dB；(d)基于分圆类构造的测量矩阵：13.1 dB。由实验结果可知，使用同一种重构算法与稀疏基的情况下，针对同一个2维图像信号，基于分圆类构造的测量矩阵恢复图像的效果优于文献[11]提出的扩展Golay序列构造的测量矩阵恢复图像的效果。

3.2.3   信道估计

考虑一般MIMO系统中的信道估计，其中导频序列从$t$个发射机传输并由$r$个接收机观察。由于每个接收机的信道估计方案都是相同的，因此可以假设系统中有一个接收机，而不会丧失通用性。假设每个发射机发射一个导频序列${{{\phi}} _i} = {[{\phi _i}(1)\;{\phi _i}(2)\; ··· \;{\phi _i}(M)]^{\rm{T}}}$其中${\rm{1}} \le i \le t$。从第$i$个发射机到接收机的信道脉冲响应(Channel Impulse Response, CIR)定义为${{{h}}_i} = {[{h_i}(1)\;{h_i}(2)\; ··· \;{h_i}(L)]^{\rm{T}}}$，其中$L$是所有信道脉冲响应上的最大信道长度。假设$M \ge L$，如果只激活第$i$个发射机，接收机观察导频序列和相应信道脉冲之间的线性卷积为$y_i^0 = {{{\phi}} _i} * {{{h}}_i} \triangleq {{\varPhi}} _i^0{{{h}}_i}$，当考虑所有$t$个发射机时，接收机的总体观测值可以表示为

$y \in {\mathbb{C}^M}$, $w$代表采样的AWGN噪声，${{\varPhi}} \in {\mathbb{C}^{M \times N}}$, ${{h}} \in {\mathbb{C}^N}$是组合信道脉冲响应。

模拟一个长度为$M=257, N=257$的多输入多输出系统，同时估计发射机$t = 2$和一个接收机之间的信道，并进行测量SNR从10 ${\rm{dB}}$到20 ${\rm{dB}}$。在每个SNR处，将组合多路径的数目从60改变到140，对应于每个信道的平均多路径数目从6增加到14。使用FPC-AS算法^[20]作为恢复算法，在每次模拟中进行500次实验，实验结果如图7、图8所示。

图 7  信道脉冲响应的实值及其估计

下载: 全尺寸图片 幻灯片

图 8  均方误差性能与组合多路径数

下载: 全尺寸图片 幻灯片

图7显示了在SNR为30 ${\rm{dB}}$的情况下，单个信道脉冲响应的实值及其估计，仿真结果表明，基于分圆类构造的测量矩阵应用在信道估计中估计值与实值完美重合。图8描述在不同SNR情况下均方误差(Mean Squared Error, MSE)性能与组合多路径数的关系，从图中可以看出，SNR将对均方误差性能产生影响，SNR越大，MSE越小，对信道估计的估计值更加精确；对应同一个SNR，基于分圆类构造的测量矩阵的MSE接近文献[11]中基于扩展FZC序列以及扩展Golay序列构造的测量矩阵的均方误差。

4.   结束语

本文基于分圆类构造了一类卷积压缩感知测量矩阵，序列的相关性从理论上确保了测量矩阵的恢复性能。实验仿真结果表明，本文构造的全实数系数矩阵不仅在1维信号的恢复中效果优于随机高斯矩阵以及已有序列构造的测量矩阵，还可以用于2维图像的重构以及信道估计。因此本文构造的测量矩阵在压缩感知以及无线通信方向具有广泛的应用前景。