留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

以欠采样速率实现盲谱感知及二维DOA估计

张展 魏平 高林 张花国

张展, 魏平, 高林, 张花国. 以欠采样速率实现盲谱感知及二维DOA估计[J]. 电子科技大学学报, 2022, 51(3): 357-362. doi: 10.12178/1001-0548.2021386
引用本文: 张展, 魏平, 高林, 张花国. 以欠采样速率实现盲谱感知及二维DOA估计[J]. 电子科技大学学报, 2022, 51(3): 357-362. doi: 10.12178/1001-0548.2021386
ZHANG Zhan, WEI Ping, GAO Lin, ZHANG Huaguo. Achieving Blind Spectrum Sensing and Two-Dimensional DOA Estimation with Sub-Nyquist Sampling Rate[J]. Journal of University of Electronic Science and Technology of China, 2022, 51(3): 357-362. doi: 10.12178/1001-0548.2021386
Citation: ZHANG Zhan, WEI Ping, GAO Lin, ZHANG Huaguo. Achieving Blind Spectrum Sensing and Two-Dimensional DOA Estimation with Sub-Nyquist Sampling Rate[J]. Journal of University of Electronic Science and Technology of China, 2022, 51(3): 357-362. doi: 10.12178/1001-0548.2021386

以欠采样速率实现盲谱感知及二维DOA估计

doi: 10.12178/1001-0548.2021386
基金项目: 国家自然科学基金(61971103)
详细信息
    作者简介:

    张展(1994 − ),男,博士生,主要从事阵列信号处理技术、sub-Nyquist技术、谱感知技术等方面的研究

    通讯作者: 魏平,E-mail:pwei@uestc.edu.cn
  • 中图分类号: TN95

Achieving Blind Spectrum Sensing and Two-Dimensional DOA Estimation with Sub-Nyquist Sampling Rate

图(6)
计量
  • 文章访问数:  4789
  • HTML全文浏览量:  1692
  • PDF下载量:  39
  • 被引次数: 0
出版历程
  • 收稿日期:  2021-12-14
  • 修回日期:  2022-02-15
  • 刊出日期:  2022-05-25

以欠采样速率实现盲谱感知及二维DOA估计

doi: 10.12178/1001-0548.2021386
    基金项目:  国家自然科学基金(61971103)
    作者简介:

    张展(1994 − ),男,博士生,主要从事阵列信号处理技术、sub-Nyquist技术、谱感知技术等方面的研究

    通讯作者: 魏平,E-mail:pwei@uestc.edu.cn
  • 中图分类号: TN95

摘要: 提出了两种以欠采样率实现盲谱感知与二维波达方向(DOA)估计的算法。这两种算法利用了MWC的周期信号所对应的傅里叶级数系数矩阵来确定未知信号频谱所在的子频带位置,然后利用二维阵列的信号模型估计得到包含二维DOA的空间相位。通过空间相位信息与子频带位置信息能重构出子带谱信息,以此能得到频率的高精度估计,最后通过空间相位分别估计得到信号的方位角与俯仰角。仿真验证了所提欠采样接收机与相应算法的有效性,以及在低信噪比环境下依然有较好的鲁棒性。

English Abstract

张展, 魏平, 高林, 张花国. 以欠采样速率实现盲谱感知及二维DOA估计[J]. 电子科技大学学报, 2022, 51(3): 357-362. doi: 10.12178/1001-0548.2021386
引用本文: 张展, 魏平, 高林, 张花国. 以欠采样速率实现盲谱感知及二维DOA估计[J]. 电子科技大学学报, 2022, 51(3): 357-362. doi: 10.12178/1001-0548.2021386
ZHANG Zhan, WEI Ping, GAO Lin, ZHANG Huaguo. Achieving Blind Spectrum Sensing and Two-Dimensional DOA Estimation with Sub-Nyquist Sampling Rate[J]. Journal of University of Electronic Science and Technology of China, 2022, 51(3): 357-362. doi: 10.12178/1001-0548.2021386
Citation: ZHANG Zhan, WEI Ping, GAO Lin, ZHANG Huaguo. Achieving Blind Spectrum Sensing and Two-Dimensional DOA Estimation with Sub-Nyquist Sampling Rate[J]. Journal of University of Electronic Science and Technology of China, 2022, 51(3): 357-362. doi: 10.12178/1001-0548.2021386
  • 无线电频谱资源匮乏,引入角度信息可以增加频谱利用率[1],这使得频谱感知问题变为二维联合谱(空间谱与频谱)的感知问题。然而,以奈奎斯特率进行宽频谱感知会带来较大的数据处理压力。利用欠采样技术可以有效地降低采样率,实现以低速率进行盲谱感知与波达方向(direction of arrival, DOA)估计。

    基于欠采样技术的谱感知与DOA估计算法主要分为两大类:基于多陪集采样器的算法[2-6]和基于调制宽带转换器(modulated wideband coverter, MWC)的算法[7-9]。然而,上述算法解决的都是盲谱感知与一维DOA估计的问题。文献[10]在文献[8]所提的L形阵列的基础上,提出了一个三维阵列。该阵列能拆分为两个L形阵列,并根据文献[8]所提出的算法,提出了一种基于ESPRIT技术的算法,实现以欠采样速率进行频率与二维DOA的联合估计。随后,文献[11]在文献[4]所提算法的基础上,在混合了多陪集与MWC结构的L形阵列上,提出了一种能在欠采样速率下进行频率与二维DOA联合估计的算法。然而,目前这样的算法有着复杂的硬件结构,且在低信噪比下性能受限。

    本文提出了一种基于单通道MWC的任意二维阵列的欠采样接收机结构,基于该结构提出了两种新的算法。这两种算法能充分利用MWC采样器的特性[3],在欠采样条件下实现盲谱感知与二维DOA估计。提出的接收机结构与算法能大幅度降低硬件资源的消耗以及复杂程度,并在低信噪比环境下,仍能保持较好的感知性能。

    • 三维空间中,$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$个阵元引入的空间相位差可以表示为:

      $$ \varphi _k^m = - 2{\text π} {f_k} \dfrac{{{{\boldsymbol{p}}^m}{{\boldsymbol{r}}_k}}}{{{c}}} $$ (1)

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

      $$ {u_m}(t) = \sum\limits_{k = 1}^K {{s_k}(t){{\text{e}}^{{\rm{j}}\varphi _k^m}}} + {w_m}(t) $$ (2)

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

      $$ {U_m}(f) = \sum\limits_{k = 1}^K {{S_k}(f){{\text{e}}^{{\rm{j}}\varphi _k^m}}} + {W_m}(f) $$ (3)

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

      $$ {\boldsymbol{U}}(f) = {\boldsymbol{AS}}(f) + {\boldsymbol{W}}(f) $$ (4)

      式中,${\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}$为间隔的频谱搬移,因此相乘后的接收信号的频谱可以表示为:

      $$ U_m^p(f) = \sum\limits_{l = - \infty }^\infty {{b_l}{U_m}(f - l {f_p})} $$ (5)

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

      $$ {U_m}({{\text{e}}^{{\rm{j}}2{\text{π}}f{T_s}}}) = \sum\limits_{l = - N + 1}^{L - 1} {{b_{ - l}}{U_m}(f + l {f_p})} \;\;\;\;f \in {{\mathcal{F}}_s} $$ (6)

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

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

      $$ {\bar U_{m,l}}(f) \triangleq {U_m}(f + l {f_p})\;\;\;\;f \in {{\mathcal{F}}_p} $$ (7)

      式中,$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}}})$定义其子带:

      $$ {Z_{m,n}}({{\text{e}}^{{\rm{j}}2{\text{π }}f{T_s}}}) \triangleq {U_m}({{\text{e}}^{{\rm{j}}2{\text{π }}(f + n{f_p}){T_s}}})\;\;\;\;f \in {{\mathcal{F}}_p} $$ (8)

      式中,$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]$以及选择矩阵组:

      $$ {{\boldsymbol{\varXi}} _n} = \left( {\begin{array}{*{20}{c}} {{{\boldsymbol{0}}_{\left( {N - n - 1} \right) \times L}}} \\ {{{\boldsymbol{I}}_{L \times L}}} \\ {{{\mathbf{0}}_{n \times L}}} \end{array}} \right) $$ (9)

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

      $$ {Z_{m,n}}({{\text{e}}^{{\rm{j}}2{\text{π }}f{T_s}}}) = {\boldsymbol{b}}{{\boldsymbol{\varXi}} _n}{\bar {\boldsymbol{U}}_m}(f) $$ (10)

      进而可以得到:

      $$ {{\boldsymbol{Z}}_m}({{\text{e}}^{{\rm{j}}2{\text{π }}f{T_s}}}) = {\boldsymbol{V}}{\bar {\boldsymbol{U}}_m}(f) $$ (11)

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

      $$ {\boldsymbol{V}} = \left( {\begin{array}{*{20}{c}} {{\boldsymbol{b}}{{\boldsymbol{\varXi}} _0}} \\ {{\boldsymbol{b}}{{\boldsymbol{\varXi}} _1}} \\ \vdots \\ {{\boldsymbol{b}}{{\boldsymbol{\varXi}} _{N - 1}}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {{{\boldsymbol{v}}_0}}&{{{\boldsymbol{v}}_1}}& \cdots &{{{\boldsymbol{v}}_{L - 1}}} \end{array}} \right) $$ (12)
    • 根据式(7)和式(8)对子带的定义,同样可以定义信号谱${S_k}(f)$的子带谱:

      $$ {\bar S_{k,l}}(f) \triangleq {S_k}(f + l {f_p})\;\;\;\;f \in {{\mathcal{F}}_p} $$ (13)

      可以将$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 U_{m,l}}(f) = {{\boldsymbol{A}}^m}{\tilde {\boldsymbol{S}}_l}(f) + {\bar W_{m,l}}(f) $$ (14)

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

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

      $$ {\bar {\boldsymbol{U}}_m}(f) = \left( {{{\boldsymbol{I}}_L} \otimes {{\boldsymbol{A}}^m}} \right)\tilde {\boldsymbol{S}}(f) + {\bar {\boldsymbol{W}}_m}(f) $$ (15)

      式中,$ \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$个阵元的欠采样信号模型:

      $$ {{\boldsymbol{Z}}_m}({{\text{e}}^{{\rm{j}}2{\text{π }}f{T_s}}}) = \left( {{{\boldsymbol{A}}^m} \otimes {\boldsymbol{V}}} \right)\bar {\boldsymbol{S}}(f) + {\boldsymbol{V}}{{\bar {\boldsymbol{W}}}_m}(f) $$ (16)

      式中,$\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$个阵元通道的输出进行组合,得到:

      $$ \begin{split} {\boldsymbol{Z}}({{{\rm{e}}} ^{{\rm{j}}2 {\text π} f{T_s}}}) = & ({\boldsymbol{A}} \otimes {\boldsymbol{V}} ){\bar {\boldsymbol{S}}} (f) + ({{\boldsymbol{I}}_M} \otimes {\boldsymbol{V}} ){\bar {\boldsymbol{W}}} (f)\\ = & {\boldsymbol{G}}\bar {\boldsymbol{S}}(f) + {{\boldsymbol{I}}_{\boldsymbol{V}}}\bar {\boldsymbol{W}}(f) \end{split} $$ (17)

      式中,$\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)$是一个长度为KLK稀疏向量。因此,信号源所在子带位置索引集合可表示为$\varOmega = \left\{ {{l_k}:0 \leqslant {l_k} \leqslant L - 1,k = 1,2, \cdots ,K} \right\}$,其中${l_k}$代表了第$k$个信号源所在子带的位置,且它们互不相同。这个稀疏向量的支撑集$\varSigma$满足:

      $$ \varSigma = \left\{ {{L_k} \triangleq \left( {k - 1} \right)L + \left( {{l_k} + 1} \right):{l_k} \in \varOmega ,k = 1,2, \cdots ,K} \right\} $$ (18)

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

      $$ {\boldsymbol{Z}}({{\text{e}}^{{\rm{j}}2{\text{π}}f{T_s}}}) = {{\boldsymbol{G}}_\varSigma }{\bar {\boldsymbol{S}}^\varSigma }(f) + {{\boldsymbol{I}}_{\boldsymbol{V}}}\bar {\boldsymbol{W}}(f) $$ (19)

      式中,${{\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$个元素可以表示为:

      $$ {Z_{m,n}}({{\text{e}}^{{\rm{j}}2{\text{π }}f{T_s}}}) = \sum\limits_{k = 1}^K {{A_{mk}}{V_{n,{l_k}}}{{\bar S}_{k,{l_k}}}(f)} + {\boldsymbol{b}}{{\boldsymbol{\varXi}} _n}{\bar {\boldsymbol{W}}_m}(f) $$ (20)

      式中,${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}}={\displaystyle \sum _{k=1}^{K}{{\boldsymbol{a}}}_{k}\circ {{\boldsymbol{v}}}_{{l}_{k}}}\circ {\bar{S}}_{k,{l}_{k}}(f) $$ (21)

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

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

      $$ \mathop {\min }\limits_{{\boldsymbol{A}},{{\boldsymbol{V}}_\varOmega },{{\boldsymbol{S}}_\varSigma }} \left\| {{\mathcal{Z}} - \sum\limits_{k = 1}^K {{{\boldsymbol{a}}_k} \circ {{\boldsymbol{v}}_{{l_k}}} \circ {{\bar S}_{k,{l_k}}}(f)} } \right\|_{\rm{F}}^2 $$ (22)

      理论上,分解得到的因子向量分别可以组成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}$包含了信号的空间相位信息组,因此可以通过下式,利用极大似然估计的方式确定其空间相位组:

      $$ ( {\hat \phi _k^x,\hat \phi _k^y} ) = \arg \mathop {\max }\limits_{\phi _k^x,\phi _k^y \in \left[ { - {\text π} ,{\text π} } \right]} \left| {{{\boldsymbol{a}}^{\text{H}}}(\phi _k^x,\phi _k^y){{\hat {\boldsymbol{a}}}_k}} \right| $$ (23)

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

      $$ \upsilon _k^{(1)} = \dfrac{{{{\boldsymbol{a}}^{\text{H}}}(\hat \phi _k^x,\hat \phi _k^y){{\hat {\boldsymbol{a}}}_k}}}{{{{\boldsymbol{a}}^{\text{H}}}(\hat \phi _k^x,\hat \phi _k^y){\boldsymbol{a}}(\hat \phi _k^x,\hat \phi _k^y)}} $$ (24)

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

      $$ {\hat l_k} = \arg \mathop {\max }\limits_{j \in \left[ {0,L - 1} \right]} \left| {\dfrac{{{{\left( {{{\hat {\boldsymbol{v}}}_k}} \right)}^H}{{\boldsymbol{v}}_j}}}{{\left\| {{{\hat {\boldsymbol{v}}}_k}} \right\|\left\| {{{\boldsymbol{v}}_j}} \right\|}}} \right| $$ (25)

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

      $$ \upsilon _k^{(2)} = \dfrac{{{\boldsymbol{v}}_{{{\hat l}_k}}^{\text{H}}{{\hat {\boldsymbol{v}}}_k}}}{{{\boldsymbol{v}}_{{{\hat l}_k}}^{\text{H}}{{\boldsymbol{v}}_{{{\hat l}_k}}}}} $$ (26)

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

      $$ {\hat {\bar S}_{k,{l_k}}}(f) = \dfrac{{{{\hat {\bar S}}_k}(f)}}{{\upsilon _k^{(3)}}} $$ (27)

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

      $$ {\hat f_k} = {\hat l_k} {f_p} + {\hat {\bar f}_k} $$ (28)

      式中,折叠后的信号中心频率${\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))$,信号源对应的方位角、俯仰角、以及重构的信号频谱可通过下式分别估计得到:

      $$ {\hat \theta _k} = \arctan \left( {\frac{{\hat \phi _k^x}}{{\hat \phi _k^y}}} \right) $$ (29)
      $$ {\hat \alpha _k} = \arcsin \left( {\frac{{\hat \phi _k^x {\text{c}}}}{{2{\text{π}} {{\hat f}_k} \sin ({{\hat \theta }_k}) d}}} \right) $$ (30)
      $$ {\hat S_k}(f) = {\hat {\bar S}_{k,{l_k}}}(f - {\hat l_k} {f_p}) $$ (31)

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

    • 若忽略式(19)的噪声项,其协方差矩阵为:

      $$ {\boldsymbol{R}} = {\text{E}}\left[ {{\boldsymbol{Z}}({{\text{e}}^{{\rm{j}}2{\text{π }}f{T_s}}}){{\boldsymbol{Z}}^{\text{H}}}({{\text{e}}^{{\rm{j}}2{\text{π }}f{T_s}}})} \right] = {{\boldsymbol{G}}_\varSigma }{{\boldsymbol{R}}_{{\boldsymbol{\bar S}}}}{\boldsymbol{G}}_\varSigma ^{\text{H}} $$ (32)

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

      $$ p\left( {l,{\phi ^x},{\phi ^y}} \right) = \dfrac{1}{{{{\left\| {{\boldsymbol{g}}_l^{\text{H}}({\phi ^x},{\phi ^y}){{\boldsymbol{u}}_N}} \right\|}^2}}} $$ (33)

      式中,${{\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 {\bar {\boldsymbol{S}}}^\varSigma }(f) = \hat {\boldsymbol{G}}_\varSigma ^\dagger {\boldsymbol{Z}}({{\text{e}}^{{\rm{j}}2{\text{π }}f{T_s}}}) $$ (34)

      式中,$\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来描述,即:

      $$ {\text{RMSE}} \triangleq \sqrt {\frac{1}{{{N_m}K}}\sum\limits_{{n_m} = 1}^{{N_m}} {\sum\limits_{k = 1}^K {{{\left( {\mu _k^{{n_m}} - \hat \mu _k^{{n_m}}} \right)}^2}} } } $$ (35)

      式中,$\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

      图23中可以看出,当系统设置满足条件时,能够完成对信号源所在子带的正确估计,因而能够估计得到信号的子带频谱,并由此重构出完整的信号谱。信号子带频谱的正确估计以及所在子带位置的正确估计保证了载波频率的高精度估计,因而对信号源的方位角和俯仰角也能正确估计。此外,仿真环境对应的奈奎斯特采样率为$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。图46分别表示了估计的载波频率、俯仰角、方位角的RMSE曲线。

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

      图  3  信号三维参数估计图

      图  4  载波频率的估计性能

      图  5  俯仰角的估计性能

      图  6  方位角的估计性能

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

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

参考文献 (12)

目录

    /

    返回文章
    返回