三维空间中，$K$个分布在宽频谱范围$ \mathcal{F}=\left[0, f_{n}\right) $内且频谱信息未知(即载频未知且基带信息未知)的独立远场窄带信号${s_k}(t)$，分别以${\varTheta _k} = ({\theta _k},{\alpha _k})$方向平行入射到接收机中，其中${\theta _k}$和${\alpha _k}$分别为方位角和俯仰角，${f_n}$为频谱范围${\mathcal{F}}$对应的奈奎斯特率。本文的主要目的是，结合阵列模型和欠采样技术，以欠采样率实现盲谱感知与二维DOA估计。此外，本文假设所有信号的信息带宽${B_k}$≤$B$，即${S_k}(f - {f_k}) \in $$ [ - B/2,{\text{ }}B/2]$，其中${S_k}(f)$和${f_k}$分别表示第$k$个信号${s_k}(t)$的连续时间傅里叶变换和对应的载波频率。

单一线阵无法实现对信号俯仰角和方位角的估计，因此至少需要将阵列的维度扩展到两维，用${{\boldsymbol{p}}^m} = [p_x^m,p_y^m,p_z^m] d$表示第$m$个阵元的位置，其中$d$代表阵列中相邻阵元的最小间距，满足$d \leqslant {{c}}/(2{f_n})$，其中${{c}}$代表光速。因此，可以表示出第$k$个信号的方向向量${{\boldsymbol{r}}_k} = - {[\sin ({\alpha _k})\sin ({\theta _k}),\sin ({\alpha _k})\cos ({\theta _k}),{\rm{cos}}({\alpha _k})]^{\rm T}}$。所以，相对参考阵元，第$k$个信号向第$m$个阵元引入的空间相位差可以表示为：

该阵元接收到的信号模型可以表示为：

式中，$m = 1,2, \cdots ,M$；${w_m}(t)$表示第m个阵元上接收到的白噪声。将上式转化到频域进行分析，可得：

式中，${U_m}(f)$和${W_m}(f)$分别表示第m个阵元接收信号和白噪声的频域形式。因此，可以将式(3)表示为阵列形式：

式中，${\boldsymbol{U}}(f) = {\left[ {{U_1}(f),{U_2}(f), \cdots ,{U_M}(f)} \right]^{\text{T}}}$代表所有阵元接收到的信号；${\boldsymbol{S}}(f) = {\left[ {{S_1}(f),{S_2}(f), \cdots ,{S_K}(f)} \right]^{\text{T}}}$代表所有的信号源；${\boldsymbol{W}}(f) = {\left[ {{W_1}(f),{W_2}(f), \cdots ,{W_M}(f)} \right]^{\text{T}}}$代表所有阵元接收到的噪声；${\boldsymbol{A}}$代表阵列的方向矩阵；第$\left( {m,k} \right)$个元素表示为${A_{mk}} = {{\text{e}}^{{\rm{j}}\varphi _k^m}}$。

图1给出了单通道的MWC的组成结构，包含一个乘法器、一个周期信号$p(t)$、一个低通滤波器$h(t)$以及一个模数转换器(analog to digital converter, ADC)，其中周期信号$p(t)$的周期间隔为${T_p}$，ADC的采样间隔为${T_s}$，低通滤波器$h(t)$通带范围为${{\mathcal{F}}_s} = [0,1/{T_s})$。此外，${T_s}$与${T_p}$满足${T_p} = N {T_s}$，其中$N$代表大于等于1的整数。因此，基带频谱范围为${{\mathcal{F}}_p} = [0,1/{T_p})$。

图  1  单通道的调制宽带转换器

MWC结构的物理本质为：当接收信号${u_m}(t)$在时域与周期信号$p(t)$相乘时，接收信号的模拟频谱${U_m}(f)$会以${f_p} = 1/{T_p}$为间隔进行连续平移，且乘上对应的尺度因子。该尺度因子为周期信号的傅里叶级数系数。当通过低通滤波器后，会保留通带范围${{\mathcal{F}}_s}$内的频谱。该频谱为${\mathcal{F}}$范围内以${f_p}$为间隔的子频带谱与对应的尺度因子相乘后组合叠加的结果。

本文提出的欠采样接收机结构由任意二维阵列和图1所示的单通道MWC组成，其中每个阵元连接相同的单通道MWC。根据式(1)得到的空间相位差可知，二维阵列的方向矩阵${\boldsymbol{A}}$可以由一组空间相位表示为$(\phi _k^x,\phi _k^y)$，其中$\phi _k^x = 2{\text{π }} {f_k}d \sin ({\alpha _k}) \sin ({\theta _k})/ {{c}}$，$\phi _k^y = 2{\text{π}}{f_k}d\sin ({\alpha _k})\cos ({\theta _k})/{{c}}$。

接收信号与周期信号相乘实际上是对频谱进行了以$1/{T_p}$为间隔的频谱搬移，因此相乘后的接收信号的频谱可以表示为：

式中，${b_l}$为周期信号的傅里叶级数系数。当该信号通过低通滤波器$h(t)$后，其频谱范围被约束在${{\mathcal{F}}_s}$范围内。经过低速率ADC采样后，得到离散序列${u_m}[n]$，其离散时间傅里叶变换可以表示为：

式中，$L = \left\lceil {{f_n}/{f_p}} \right\rceil $表示了频谱${\mathcal{F}}$可以被划分出宽度为${f_p}$的子带个数。

整个观测频谱可以划分出$L$个子带，因此可以由${U_m}(f)$定义其宽度为${f_p}$的子带：

式中，$l = 0,1,2, \cdots ,L - 1$，这$L$个子带频谱能排成向量形式${\bar {\boldsymbol{U}}_m}(f) = [{\bar U_{m,0}}(f),{\bar U_{m,1}}(f),{\bar U_{m,2}}(f), \cdots , {\bar U_{m,L - 1}}(f)]^{\text{T}}$。同理，由于通过通带范围为${{\mathcal{F}}_s}$的低通滤波器后，其通带内含有$N$个宽度为${f_p}$的子带，因此可以由${U_m}({{\text{e}}^{{\rm{j}}2{\text{π }}f{T_s}}})$定义其子带：

式中，$n = 0,1,2, \cdots ,N - 1$，同样也能排成向量形式${{\boldsymbol{Z}}_m}({{\text{e}}^{{\rm{j}}2{\text{π }}f{T_s}}}) = {[{Z_{m,0}}({{\text{e}}^{{\rm{j}}2{\text{π }}f{T_s}}}),{Z_{m,1}}({{\text{e}}^{{\rm{j}}2{\text{π }}f{T_s}}}), \cdots ,{Z_{m,N - 1}}({{\text{e}}^{{\rm{j}}2{\text{π }}f{T_s}}})]^{\text{T}}}。$

为方便描述，定义${\boldsymbol{b}} = \left[ {{b_{N - 1}},{b_{N - 2}}, \cdots ,{b_{ - L + 1}}} \right]$以及选择矩阵组：

式中，$n = 0,1,2, \cdots ,N - 1$。因此，可以由式(6)得到：

进而可以得到：

式中，系数矩阵${\boldsymbol{V}}$可以表示为：

根据式(7)和式(8)对子带的定义，同样可以定义信号谱${S_k}(f)$的子带谱：

可以将$K$个信号源对应位置的子带谱排列为${\tilde {\boldsymbol{S}}_l}(f) = {[{\bar S_{1,l}}(f),{\bar S_{2,l}}(f), \cdots ,{\bar S_{K,l}}(f)]^{\rm T}}$。$L$个子带对应的向量可得到$\tilde {\boldsymbol{S}}(f) = {[\tilde {\boldsymbol{S}}_0^{\rm T}(f),\tilde {\boldsymbol{S}}_1^{\rm T}(f), \tilde {\boldsymbol{S}}_2^{\rm T}(f), \cdots ,\tilde {\boldsymbol{S}}_{L - 1}^{\rm T}(f)]^{\rm T}}$。

因此，由式(4)和子带谱定义式(7)可以得到：

式中，$ {\bar W_{m,l}}(f) \triangleq {W_m}(f + l {f_p}) $表示第$m$个天线接收噪声的第$l$个子带频谱。

因此，可以将第$m$个天线的接收信号的子频带谱按前文定义的方式进行重排，得到：

式中，$ \otimes $表示Kronecker积；${{\boldsymbol{A}}^m}$为方向矩阵$ {\boldsymbol{A}} $的第$m$行；${\bar {\boldsymbol{W}}_m}(f) = {[{\bar W_{m,0}}(f),{\bar W_{m,1}}(f),{\bar W_{m,2}}(f), \cdots ,{\bar W_{m,L - 1}}(f)]^{\rm T}}$。

将上式带入式(10)中，可得第$m$个阵元的欠采样信号模型：

式中，$\bar {\boldsymbol{S}}(f) = {[\bar {\boldsymbol{S}}_1^{\rm T}(f),\bar {\boldsymbol{S}}_2^{\rm T}(f), \cdots ,\bar {\boldsymbol{S}}_K^{\rm T}(f)]^{\rm T}}$是区别于$\tilde {\boldsymbol{S}}(f)$的另一种重排所有信号源的子带谱的方式，因此${\bar {\boldsymbol{S}}_k}(f) = {[{\bar S_{k,0}}(f),{\bar S_{k,1}}(f),{\bar S_{k,2}}(f), \cdots ,{\bar S_{k,L - 1}}(f)]^{\rm T}}$为第$k$个信号源的所有子带谱排成的向量形式。进而可以将$M$个阵元通道的输出进行组合，得到：

式中，$\bar {\boldsymbol{W}}(f) = {[\bar {\boldsymbol{W}}_1^{\rm T}(f),\bar {\boldsymbol{W}}_2^{\rm T}(f), \cdots ,\bar {\boldsymbol{W}}_M^{\rm T}(f)]^{\rm T}}$。

由于信号${s_k}(t)$被认为是窄带信号，若假定每一个信号仅存在于一个子带内，则意味着$\bar S(f)$是一个长度为KL的K稀疏向量。因此，信号源所在子带位置索引集合可表示为$\varOmega = \left\{ {{l_k}:0 \leqslant {l_k} \leqslant L - 1,k = 1,2, \cdots ,K} \right\}$，其中${l_k}$代表了第$k$个信号源所在子带的位置，且它们互不相同。这个稀疏向量的支撑集$\varSigma$满足：

欠采样接收模型式(17)可以进一步简化为：

式中，${{\boldsymbol{V}}_\varOmega } = [{{\boldsymbol{v}}_{{l_1}}},{{\boldsymbol{v}}_{{l_2}}} \cdots ,{{\boldsymbol{v}}_{{l_K}}}]$；${{\boldsymbol{G}}_\varSigma } = {\boldsymbol{A}} \odot {{\boldsymbol{V}}_\varOmega }$，$ \odot $代表Khatri-Rao积；${\bar {\boldsymbol{S}}^\varSigma }(f) = {[{\bar S_{1,{l_1}}}(f),{\bar S_{2,{l_2}}}(f), \cdots ,{\bar S_{K,{l_K}}}(f)]^{\rm T}}$。

根据式(19)可知，${\boldsymbol{Z}}({{\text{e}}^{{\rm{j}}2{\text{π}}f{T_s}}})$是一个$MN$个元素的向量，该向量的第$\left( {m-1} \right) \times N+n$个元素可以表示为：

式中，${V_{n,{l_k}}}$代表矩阵${\boldsymbol{V}}$的第$\left( {n,{l_k}} \right)$个元素。实际上，${Z_{m,n}}({{\text{e}}^{{\rm{j}}2{\text{π }}f{T_s}}})$是一个$T$点的DFT序列。因此，忽略上式噪声项后可以表示为一个秩为$K$的三阶张量：

式中，∘代表张量的外积；${\mathcal{Z}} \in {{\mathbb{C}}^{M \times N \times T}}$。

张量的秩一分解需要满足唯一分解条件，该条件的分析可以参考文献[9]。假设该条件成立，此时${\mathcal{Z}}$的唯一分解可以通过求解下式得到：

理论上，分解得到的因子向量分别可以组成3个因子矩阵，具体为：${\boldsymbol{A}} = \left[ {{{\boldsymbol{a}}_{\text{1}}},{{\boldsymbol{a}}_2}, \cdots ,{{\boldsymbol{a}}_K}} \right]$、${{\boldsymbol{V}}_\Omega } = [{{\boldsymbol{v}}_{{l_1}}},{{\boldsymbol{v}}_{{l_2}}}, \cdots , $$ {{\boldsymbol{v}}_{{l_K}}}]$、${{\boldsymbol{S}}_\varSigma } = [{\bar S_{1,{l_1}}}(f),{\bar S_{1,{l_2}}}(f), \cdots ,{\bar S_{K,{l_K}}}(f)]$。具体的解法可以参考文献[12]中的多种迭代最小二乘算法。

虽然${\mathcal{Z}}$唯一分解，却仍存在排列和尺度上的不确定性，但这并不影响配对分解得到的因子向量，即$({\hat {\boldsymbol{a}}_k},{\hat {\boldsymbol{v}}_k},{{\hat {\bar s}}_k}(f))$，其中分解得到的因子向量与理论上的因子向量有如下关系：${\hat {\boldsymbol{a}}_k} = \upsilon _k^{(1)}{{\boldsymbol{a}}_k}$；${\hat {\boldsymbol{v}}_k} = \upsilon _k^{(2)}{{\boldsymbol{v}}_{{l_k}}}$；${{\hat {\bar S}}_k}(f) = \upsilon _k^{(3)}{\bar S_{k,{l_k}}}(f)$；$\displaystyle\prod _{i = 1}^3\upsilon _k^{(i)} = 1$，其中$\upsilon _k^{(i)}$为因子向量对应的尺度因子。

由于分解得到的向量${\hat {\boldsymbol{a}}_k}$包含了信号的空间相位信息组，因此可以通过下式，利用极大似然估计的方式确定其空间相位组：

其对应的缩放因子$\upsilon _k^{(1)}$可以通过下式计算得到：

相似的是，分解得到的因子向量${\hat v_k}$包含了信号源所在频带的位置信息，因此其频带位置${l_k}$可以通过在傅里叶系数矩阵${\boldsymbol{V}}$中寻找与之最相似的列来得到，即通过向量投影的方式：

其对应的缩放因子$\upsilon _k^{(2)}$可以通过下式计算得到：

根据缩放因子间的关系，可以计算得到$\upsilon _k^{(3)} = 1/(\upsilon _k^{(1)}\upsilon _k^{(2)})$。因此，信号源所在子带的频谱可以通过下式得到：

根据折叠频率${\bar f_k}$的定义${\bar f_k} \triangleq \boldsymbolod \left( {{f_k},{f_p}} \right)$可知，如果能够得到折叠后的信号中心频率${\hat {\bar f}_k}$，即能够通过下式求得信号的载波频率：

式中，折叠后的信号中心频率${\hat {\bar f}_k}$可以通过多种方式从重构的子带谱$ {\hat {\bar S}_{k,{l_k}}}(f) $中得到，如通过滑动平均(类似CFAR检测过门限)的方式得到子带谱$ {\hat {\bar S}_{k,{l_k}}}(f) $的频率起始点和截止点，由此计算出折叠频率。

通过已配对估计的参数组$(\hat \phi _k^x,\hat \phi _k^y,{\hat f_k},{\hat {\bar S}_{k,{l_k}}}(f))$，信号源对应的方位角、俯仰角、以及重构的信号频谱可通过下式分别估计得到：

此时，已经实现了对信号源的盲谱感知及二维DOA估计，将该方法命名为JSS-2DOA-CP算法。

若忽略式(19)的噪声项，其协方差矩阵为：

式中，${{\boldsymbol{R}}_{{\boldsymbol{\bar S}}}}$是${\bar {\boldsymbol{S}}^\varSigma }(f)$的协方差矩阵。因此，可以通过对${\boldsymbol{R}}$使用子空间分解的方法，构造其伪谱：

式中，${{\boldsymbol{g}}_l}({\phi ^x},{\phi ^y}) = {\boldsymbol{a}}({\phi ^x},{\phi ^y}) \odot {{\boldsymbol{v}}_l}$是式(19)的导向矢量，其中包含了信号的空间相位信息和信号源所在子带的位置信息；${{\boldsymbol{u}}_N}$为${\boldsymbol{R}}$分解得到的噪声子空间。因此，伪谱中最大的$K$个峰值的索引信息对应了信号的空间相位信息和频带位置信息$ ({\hat l_k},\hat \phi _k^x,\hat \phi _k^y) $。

当得到信号源的相位信息和频带位置信息后，导向矢量${{\boldsymbol{g}}_{{{\hat l}_k}}}(\hat \phi _k^x,\hat \phi _k^y)$也就能估计得到，因此其阵列流行矩阵${\hat {\boldsymbol{G}}_\varSigma }$也能得到。根据最小二乘原理，信号所在频带的频谱能通过下式重构得到：

式中，$\hat {\boldsymbol{G}}_\varSigma ^\dagger = {\left( {\hat {\boldsymbol{G}}_\varSigma ^{\text{H}}{{\hat {\boldsymbol{G}}}_\varSigma }} \right)^{ - 1}}\hat {\boldsymbol{G}}_\varSigma ^{\text{H}}$。此时，信号源的相位信息、频带位置信息以及对应的子带频谱信息都已经得到且完成配对，接下来可以通过式(28)～(31)实现盲谱感知及二维DOA估计，最后将该方法命名为JSS-2DOA-MUSIC算法。

本节通过两个仿真实验分别验证了所提的接收机结构以及对应的盲感知算法的有效性，和它们在不同信噪比下的感知性能。这两个仿真实验的部分参数设置相同：奈奎斯特率为$10$ GHz；信号源的带宽为$50$ MHz；MWC周期信号的基频为$ f_{p}=100 $ MHz，因此共有$L = 100$个子频带；相邻阵元的最小间距为$d = 0.5 \times {{c}}/{f_n}$；ADC的采样率为${f_s} = 400$ MHz，因此$N = 4$。此外，信噪比为${\text{SNR}} \triangleq {\text{E[|}}s(t){{\text{|}}^2}{\text{]}}/{\sigma ^2}$，且参数估计的性能通过均方根误差RMSE来描述，即：

式中，$\mu $和$\hat \mu $分别表示估计量的真实值和估计值；下标$k$表示第$k$个信号源；上标${n_m}$表示第${n_m}$次蒙特卡洛随机试验。

第一个仿真实验验证所提采样接收机结构以及对应盲感知算法的有效性。该仿真设置了11个$50$ MHz带宽的复信号源。为了满足文献[9]中分析的系统参数选择条件，二维阵列设置为包含$M = 9$个天线的L形阵列，即每个子阵列分别包含5个阵元。信噪比SNR设置为$0$ dB。图2描述了原始信号频谱以及JSS-2DOA-CP算法和JSS-2DOA-MUSIC算法重构的信号频谱，并给出了其中一个信号源的重构频谱放大图。图3给出了信号源的三维参数的估计图，图中的3个坐标轴分别代表方位角$\theta $、俯仰角$\alpha $、载波频率f。

从图2～3中可以看出，当系统设置满足条件时，能够完成对信号源所在子带的正确估计，因而能够估计得到信号的子带频谱，并由此重构出完整的信号谱。信号子带频谱的正确估计以及所在子带位置的正确估计保证了载波频率的高精度估计，因而对信号源的方位角和俯仰角也能正确估计。此外，仿真环境对应的奈奎斯特采样率为$10$ GHz，而该系统用到了9个天线，每个RF电路的采样率为$400$ MHz，因此系统总采样率为$3.6$ GHz，远低于奈奎斯特率。该仿真也说明了所提出的采样结构以及对应的算法能够保证用以远低于奈奎斯特率的总采样率，在噪声环境下，实现盲谱感知以及二维DOA估计。

第二个实验是为了验证所提算法的性能，选择了文献[10]所提出的基于CaSCADE^[8]的衍生算法和文献[11]所提出的基于PABSS^[4]的衍生算法作为本仿真实验的对比算法。该仿真设置了两个带宽为50 MHz的复信号源。所提的接收机结构是包含7个天线的L形阵列，即$M = 7$。为了保证公平性，即系统的总采样率一致，设文献[10]中算法对应的三维阵列，每个子阵列包含10个天线，总计28个天线；设文献[11]所提的L阵列，每个子阵包含15个天线，总计29个天线。因此，本文所提的系统总采样率为$2.8$ GHz，文献[10]和文献[11]算法的总采样率分别为$2.8$ GHz和$2.9$ GHz。仿真的信噪比SNR变化范围以3 dB为间隔，设置为$ - 10$～$20$ dB。图4～6分别表示了估计的载波频率、俯仰角、方位角的RMSE曲线。

图  2  原始信号频谱及重构的信号频谱

图  3  信号三维参数估计图

图  4  载波频率的估计性能

图  5  俯仰角的估计性能

图  6  方位角的估计性能

从图中可以看出，所提算法的性能相比于另外两种算法，在谱重构以及频率估计的性能上有着显著优势，这是因为本文算法的载波频率估计依赖于对信号源所处子频带位置的估计，因此频率的估计精度会大幅度提高，进而带来谱重构的优势。此外，方位角和俯仰角的估计性能在低信噪比情况下也有显著优势。不仅如此，相比于对比算法的系统设计，本文所提出的系统结构有更少的硬件资源开销，因此本文所提的接收机结构以及对应的盲感知算法，能在降低资源开销的同时仍保证感知性能的鲁棒性。

本文提出的基于压缩采样的盲谱感知及二维DOA估计算法，充分利用了MWC中的周期信号产生对应的傅里叶级数系数矩阵，通过估计出信号源频谱所处的子带位置，以此提高载波频率的估计性能，进而使得在低信噪比环境下，信号源的方位角和俯仰角的估计精度大幅度提高。这种区别于传统的利用空间相位信息估计载频与二维度DOA的思路，大幅度提高了感知性能以及估计的鲁棒性。此外，设计的接收机系统以及对应的感知算法能有效降低系统硬件资源的开销，对认知无线电中多维联合谱感知问题具有重要意义。