波达方向估计(direction of arrival, DOA)是阵列信号处理研究的热点问题, 在雷达、通信、电子对抗等领域都有着广泛的应用.传统的高分辨率测向方法如多重分类方法(multiple signal classification, MUSIC)、最小方差无畸变(minimum variance distortionless response, MVDR)方法等都是在入射信号频率相同且已知的情况下进行估计的.而随着现代电子战的信号环境日趋复杂, 测向面临着一个很宽的频段, 测向接收机的瞬时带宽朝着宽带发展, 在一个瞬时带宽内必然包含有很多不同频率的窄带信号且信号数可能会超过阵元数.因此, 如何实现对这些窄带信号的DOA进行准确估计是一个巨大的挑战.若直接利用传统的DOA估计方法会产生严重的估计误差.

近年来, 针对宽频段的DOA估计问题得到了广大学者的重视^[1-4].但现有成果都是建立在均匀线阵基础上进行研究的, 可估计信号数要少于阵元数, 严重损失了阵列的自由度.李鹏飞等^[5]提出基于空频域稀疏表示的宽频段DOA估计方法, 虽然可以用于欠定情况, 但不足在于通过对感兴趣的角度域进行均匀量化来构造完备字典, 未考虑实际情况下DOA信息与完备字典的匹配关系, 导致模型失配时估计性能会急剧下降.虽然更密集的网格有助于减小重构误差提高估计精度, 但过度密集的离散字典网格会使得完备字典的列相关性显著增加, 不能保证完备字典满足约束等容条件(restricted isometry property, RIP)^[6].

最近, 二级嵌套阵列^[7]作为一种新型的稀疏阵列结构被提出, 与相同物理阵元数(设为N)的均匀线阵相比, 二级嵌套阵列可通过虚拟变换使自由度扩展至O(N²).同时具有易于构造、物理阵元和虚拟阵元具有解析表达式等优点受到广泛关注, 为解决欠定DOA估计问题提供了一个新的思路, 一些有价值的研究成果相继出现^[8-14].而现有的研究成果并不能直接应用到宽频段欠定DOA估计问题中.

针对上述问题, 本文研究了一种基于二级嵌套阵列的宽频段欠定DOA估计方法.首先通过引入空间频率的概念, 对二级嵌套阵列接收数据进行处理, 将频率和DOA二维信息转化为一维空间频率信息；然后利用空间频率的空域稀疏性构造连续稀疏模型, 利用原始对偶方法以及多项式求根得到空间频率的高分辨估计；最后建立频域协方差矩阵大特征矢量之和稀疏模型对信号频率和空间频率进行配对, 最终实现信号频率与DOA的联合估计.该方法不需要对角度域进行离散化处理来构造完备字典, 避免了模型失配问题对估计性能的影响.因此, 具有更高的估计精度和分辨性能.

1 信号模型

二级嵌套阵列的物理结构如图 1所示, 二级嵌套阵列由2个不同阵元间距的均匀线阵组成, 其中第1级包含有N₁个阵元, 阵元间距为d₁, 第2级包含有N₂个阵元, 阵元间距为d₂=(N₁+1)d₁, 所有阵元位置集合为S={n₁d₁, n₁=0, 1, …, N₁-1}∪{n₂d₂-d₁, n₂=1, …, N₂}, 以第1个阵元为参考阵元.经过Khatri-Rao积, 二级嵌套阵列可以产生2N₂(N₁+1)-1个阵元的虚拟均匀线阵, 阵元间距为d₁, 各阵元位置为{nd₁, n=-N_i, …, N_i, N_i=N₂(N₁+1)-1}, 如图 2所示.

假设K个远场窄带平面波窄带信号分别从方向[θ₁, θ₂, …, θ_K]入射到二级嵌套阵列上, 每个信号的中心频率分别为[f₁, f₂, …, f_K].则第i个信号的导向矢量a(θ_i, f_i)可表示为

$ \begin{array}{l} \mathit{\boldsymbol{a}}({\theta _i}, {f_i}) = [1, {\rm{exp}}(-{\rm{j2 \mathsf{ π} }}{d_1}{f_i}{\rm{sin}}{\theta _i}/c), \ldots, \\ \;\;\;{\rm{exp}}(-{\rm{j2 \mathsf{ π} }}({N_1}-1){d_1}{f_i}{\rm{sin}}{\theta _i}/c), \ldots, \\ \;\;\;{\rm{exp}}( - {\rm{j2 \mathsf{ π} }}({N_2}{d_2} - {d_1}){f_i}{\rm{sin}}{\theta _i}/c){]^{\rm{T}}}. \end{array} $

(1)

c=3×10⁸ m/s.t时刻阵列接收数据可表示为

$ \mathit{\boldsymbol{X}}(t) = \sum\limits_{i = 1}^K {\mathit{\boldsymbol{a}}({\theta _i},{f_i}){s_i}(t) + \mathit{\boldsymbol{n}}(t)} = \mathit{\boldsymbol{As}}(t) + \mathit{\boldsymbol{n}}(t). $

(2)

式中:c为电磁传播速度, X(t)=[x₁(t), x₂(t), …, x_N1+N2(t)]^T表示所有阵元在采样时刻t时的输出信号；A=[a(θ₁, f₁), a(θ₂, f₂), …, a(θ_K, f_K)]表示所有入射信号DOA对应的(N₁+N₂)×K维阵列流型矩阵；s(t)=[s₁(t), s₂(t), …, s_K(t)]^T表示采样时刻t时的所有入射信号构成的K×1维信号向量；n(t)=[n₁(t), n₂(t), …, n_N1+N2(t)]^T表示阵元接收到的(N₁+N₂)×1维高斯白噪声向量, 其均值为0, 方差为σ².

假设K个入射信号之间互不相关, 则接收数据X(t)的协方差矩阵可表示为

$ \begin{array}{*{20}{l}} {{\mathit{\boldsymbol{R}}_X} = E[\mathit{\boldsymbol{X}}(t){\mathit{\boldsymbol{X}}^{\rm{H}}}(t)] = } \\ {\;\;\;\;\;\;\;\;\;\mathit{\boldsymbol{A}}{\rm{diag}}({P_1},{P_2}, \ldots ,{P_K}){\mathit{\boldsymbol{A}}^{\rm{H}}} + {\sigma ^2}\mathit{\boldsymbol{I}} \cdot } \end{array} $

(3)

式中：E[·]表示期望运算, (·)^H表示共轭转置, diag(·)表示求对角矩阵, P_i表示第i个入射信号的功率, I表示单位矩阵.

2 算法原理 2.1 数据预处理

由式(1) 可知, 宽频段DOA估计是一个二维参数估计问题, 直接求解计算量大.为了降低计算复杂度, 引入空间频率概念^[5], 空间频率定义为

$ {f_s} = {d_1}f{\rm{sin}}\theta /c. $

(4)

式中，f为入射信号的频率，θ为入射信号的方位角.f_sⁱ=d₁f_isinθ_i/c表示第i个入射信号的空间频率, 为了保证无模糊测向, 这里取d₁f_max/c≤0.5, 因此空间频率范围为-0.5≤f_s≤0.5.则式(1) 可表示为

$ \begin{array}{l} \mathit{\boldsymbol{a}}(f_s^i) = [1, {\rm{exp}}(-{\rm{j2 \mathsf{ π} }}f_s^i), \ldots, {\rm{exp}}(-{\rm{j2 \mathsf{ π} }}{N_1}f_s^i), \ldots, \\ \;\;\;\;\;\;\;\;\;\;\;{\rm{exp}}(-{\rm{j2 \mathsf{ π} }}({N_2}{N_1} + {N_2} - 1)f_s^i){]^{\rm{T}}}. \end{array} $

(5)

对协方差矩阵R_X进行向量化操作, 即将各列依次堆叠为一列可得(N₁+N₂)²×1维向量Z：

$ \mathit{\boldsymbol{Z}} = {\rm{vec}}\left( {{\mathit{\boldsymbol{R}}_\mathit{\boldsymbol{X}}}} \right) = \mathit{\Phi }\left( {f_s^1,f_s^2, \ldots ,f_s^K} \right)\mathit{\boldsymbol{P}} + {\sigma ^2}{{\bf{1}}_n}. $

(6)

$ \begin{array}{l} \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}(f_s^1, \ldots, f_s^K) = {\mathit{\boldsymbol{A}}^*} \odot \mathit{\boldsymbol{A}} = \left[{{\mathit{\boldsymbol{a}}^*}(f_s^1) \otimes \mathit{\boldsymbol{a}}(f_s^1), \ldots, } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {{\mathit{\boldsymbol{a}}^*}(f_s^K) \otimes \mathit{\boldsymbol{a}}(f_s^K)} \right]. \end{array} $

(7)

式中：vec (·)表示矩阵向量化操作, (·)^*表示共轭运算, ☉表示求Khatri Rao积操作, $ \otimes $表示求Kronecker积操作, 1_n=[e₁^T, e₂^T, …, e^T_N1+N2]^T, e_i表示除了第i个元素为1, 其余元素为0的(N₁+N₂)×1维列向量.对比式(2) 与(6) 可知, Φ的作用相当于式(2) 中的阵列流型矩阵A, 可看作一个扩展的虚拟阵列流型矩阵, Φ中共有2N₂(N₁+1)-1个不重复的行, 从而引起向量Z中某些元素不具有唯一性.因此, 剔除向量Z中重复的元素组成一个新向量$ \mathit{\boldsymbol{\tilde Z}} $, 则$ \mathit{\boldsymbol{\tilde Z}} $可表示为

$ \mathit{\boldsymbol{\tilde Z}} ＝ {\rm{ }}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_1}s + {\sigma ^2}{\mathit{\boldsymbol{w}}_1}. $

(8)

式中：$ \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{{\rm{ }}_1} = \left[ {\mathit{\boldsymbol{\tilde a}}\left( {f_s^1} \right),\mathit{\boldsymbol{\tilde a}}\left( {f_s^2} \right), \cdot \cdot \cdot ,\tilde a\left( {f_s^K} \right)} \right] $可看作为一个阵元数为2N₂(N₁+1)-1的虚拟线阵所对应的阵列流型矩阵, 其对应的导向矢量为

$ \begin{array}{l} \mathit{\boldsymbol{\tilde a}}({\eta _k}) = [{\rm{exp}}(-{\rm{j2 \mathsf{ π} }}({N_2}({N_1} + 1)-1)f_s^k), \ldots, \\ \;\;\;\;\;\;\;\;1, \ldots, {\rm{exp}}({\rm{j2 \mathsf{ π} }}({N_2}({N_1} + 1)-1)f_s^k){]^{\rm{T}}}. \end{array} $

(9)

阵元位置集合为{nd₁, n=-(N₂(N₁+1)-1), …, N₂(N₁+1)-1}. w₁表示除第N₂(N₁+1)-1元素为1外, 其余元素为0的(2N₂(N₁+1)-1)×1维列向量.通过比较式(9) 与(4) 可知, 经过向量化处理阵列孔径得到了扩展, 提高了阵列的自由度, 可实现对多于阵元数的信号进行DOA估计.

2.2 空间频率估计

由式(5) 可知, 通过引入空间频率, 宽频段DOA估计问题由一个频率和DOA二维信息估计问题降维成一维空间频率估计问题, 要想实现对DOA的高精度估计, 必须对空间频率具有较高的估计精度和分辨能力.本文引入超分辨理论^[15-16]用于空间频率的估计, 能够有效避免传统稀疏重构方法存在的模型失配问题对估计性能的影响, 提高估计精度和分辨性能.

将式两边同乘以exp (-jπn)并进行变换, 可得

$ \begin{array}{*{20}{l}} {y\left( n \right) = {\rm{exp}}\left( { - {\rm{j\pi }}n} \right)\left( {\tilde Z\left( n \right) - {\sigma ^2}{w_1}\left( n \right)} \right) = }\\ {{\rm{exp}}\left( { - {\rm{j\pi }}n} \right)\sum\limits_{k = 1}^K {{P_k}{\rm{exp}}} \left( {{\rm{j2\pi }}nf_s^k} \right) = }\\ {\sum\limits_{k = 1}^K {{P_k}{\rm{exp}}\left( { - {\rm{j2\pi }}n{\tau _k}} \right)} = \int_0^1 {P\left( \tau \right){\rm{exp}}\left( { - {\rm{j2\pi }}n{\tau _k}} \right){\rm{d}}\tau } \cdot } \end{array} $

(10)

式中:τ_k=(1-2f_s^k)/2, 0≤τ_k≤1, n=-(N₂(N₁+1)-1), …, N₂(N₁+1)-1, τ={τ_k, 1≤k≤K}.进一步式(10) 可写成矩阵形式

$ \mathit{\boldsymbol{y}} ＝ \mathit{\boldsymbol{FP}}. $

(11)

其中：F_{m, n}=exp(-j2πmτ_n), y=[y(-N²/2-N+2), …, y(N²/2+N-2)]^T, P(τ)在连续域上可表示为

$ P(\tau ) = \sum\limits_{k = 1}^K {{P_k}{\delta _{{\tau _k}}} \cdot } $

(12)

式中:δ_τk表示在τ_k的Dirac测量, P_k为第k个信号的功率.

为实现对τ_k的高精度估计, 式(12) 所表示的连续稀疏问题可以转化为相应的优化问题进行求解

$ \mathop {{\rm{min}}}\limits_P {\mkern 1mu} {\left\| \mathit{\boldsymbol{P}} \right\|_{{\rm{TV}}}}{\rm{s}}.{\rm{t}}.\mathit{\boldsymbol{y}} = \mathit{\boldsymbol{FP}} \cdot $

(13)

式中:$ {{\left\| \ \ \right\|}_{\text{TV}}} $表示全变差范数, 用于表示连续信号的稀疏性, $ {\left\| \mathit{\boldsymbol{P}} \right\|_{{\rm{TV}}}} $可表示为

$ {\left\| \mathit{\boldsymbol{P}} \right\|_{{\rm{TV}}}} = {\rm{sup}}\sum\limits_{j = 1}^\infty {\left\| {\mathit{\boldsymbol{P}}({\mathit{\boldsymbol{B}}_j})} \right\|} . $

(14)

式中:sup表示取上确界, P(B_j)表示集合[0, 1]的不重合子集, 当P是稀疏的, 此时$ {\left\| \mathit{\boldsymbol{P}} \right\|_{{\rm{TV}}}} = \sum\limits_{i = 1}^K {{P_i}} $等效于离散l₁范数.

在实际DOA估计中, 噪声功率σ²通常未知, 式(13) 可进一步等效为如下的优化问题：

$ \mathop {{\rm{min}}}\limits_{P, {\sigma ^2} \ge 0} {\left\| \mathit{\boldsymbol{P}} \right\|_{{\rm{TV}}}}{\rm{s}}{\rm{.t}}{\rm{.}}\left\| {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{FP}} - {\sigma ^2}{\mathit{\boldsymbol{w}}_2}} \right\|{{\rm{ }}_2} \le \mu \cdot $

(15)

式中：w₂(n)=w₁(n)e^-jπn, μ为正则化参数, 式(15) 对应的对偶问题可表示为

$ \begin{array}{l} \;\;\;\;\;\;\;\;\mathop {{\rm{max}}}\limits_u {\rm{Re}}[{\mathit{\boldsymbol{u}}^{\rm{H}}}y] - \mu {\left\| \mathit{\boldsymbol{u}} \right\|_2}\\ {\rm{s}}{\rm{.t}}{\rm{.}}{\left\| {{\mathit{\boldsymbol{F}}^{\rm{H}}}\mathit{\boldsymbol{u}}} \right\|_\infty } \le 1, Re[{\mathit{\boldsymbol{u}}^{\rm{H}}}{\mathit{\boldsymbol{w}}_2}] \le 0 \cdot \end{array} $

(16)

式中：Re()为求实部，u为对偶变量.式(16) 是一个无穷维的优化问题不能直接求解, 但第一个约束条件可等效为一个半正定矩阵约束^[15].因此, 该无穷维对偶问题可等效为相应的半正定规划问题

$ \begin{array}{l} \mathop {{\rm{max}}}\limits_{u, Q} {\rm{Re}}\left[{{\mathit{\boldsymbol{u}}^{\rm{H}}}y} \right] - \mu {\left\| \mathit{\boldsymbol{u}} \right\|_2}\\ {\rm{s}}.{\rm{t}}.\left[{\begin{array}{*{20}{l}} \mathit{\boldsymbol{Q}}\\ {{\mathit{\boldsymbol{u}}^{\rm{H}}}} \end{array}} \right.\left. {\begin{array}{*{20}{l}} \mathit{\boldsymbol{u}}\\ 1 \end{array}} \right] \ge 0, {\rm{Re}}\left[{{\mathit{\boldsymbol{u}}^{\rm{H}}}{\mathit{\boldsymbol{w}}_2}} \right] \le 0\\ \sum\limits_{i = 1}^{2{N_2}({N_1} + 1) - 1 - j} {{\mathit{\boldsymbol{Q}}_{i, i + j}}} = \left\{ \begin{array}{l} 1, j = 0;\\ 0, j = 1, 2, \ldots, {N_2}({N_1} + 1) - 1 \cdot \end{array} \right. \end{array} $

(17)

式中：Q为N₂(N₁+1)×N₂(N₁+1) 维埃尔米特矩阵, 对于式(17) 可利用CVX优化工具箱^[17]进行求解, 进一步对多项式$ 1 - {\left\| {{\mathit{\boldsymbol{F}}^{\rm{H}}}\mathit{\boldsymbol{u}}\left( \tau \right)} \right\|^{\rm{2}}} = 0 $求根得到τ的估计值τ_e, 从而得到空间频率f_s的估计值$ {\hat f_s} ＝ 1/2 - {\tau _e} $.

2.3 DOA与信号频率的匹配

由式(4) 空间频率的定义可知, 若知道所有入射信号频率并与相应的空间频率正确匹配, 就可得到入射信号的DOA信息.因此, DOA与信号频率的匹配等同于空间频率与信号频率的匹配.

取式(2) 的傅里叶变换, 可得阵列接收数据的频域模型：

$ \mathit{\boldsymbol{Y}}\left( {{f_k}} \right) = \mathit{\boldsymbol{A}}({f_k})\mathit{\boldsymbol{s}}({f_k}) + \mathit{\boldsymbol{N}}({f_k}), k = 1, 2, \ldots, K. $

(18)

式中：Y(f_k)=[y₁(f_k), y₂(f_k), …, y_N1+N2(f_k)]^T为阵列接收数据X(t)经过傅里叶变换在频率f_k处的阵列输出向量；s(f_k)=[s₁(f_k), s₂(f_k), …, s_K(f_k)]^T为K个信号在频率f_k处的输出向量；N(f_k)=[N₁(f_k), N₂(f_k), …, N_N1+N2(f_k)]^T为噪声在频率f_k处的输出向量.

空间频率与信号频率的配对过程进一步等效为寻找一个或多个空间频率与Y(f_k)相对应, 该问题可利用稀疏重构的思想进行求解.对于式(18), 如果直接构造频域稀疏模型容易受到噪声的影响, 导致配对的准确率下降.为了提高配对的准确性, 本文利用信号协方差矩阵大特征值对应的特征矢量为各信号方向矢量的一个线性组合的性质^[19]来建立配对函数.

假设信号与噪声相互独立, 中心频率为f_i的阵列协方差矩阵为

$ \begin{gathered} \mathit{\boldsymbol{R}}\left( {{f_i}} \right) = E\left[ {\mathit{\boldsymbol{Y}}\left( {{f_i}} \right){\mathit{\boldsymbol{Y}}^{\rm{H}}}\left( {{f_i}} \right)} \right] = \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\mathit{\boldsymbol{A}}\left( {{f_i}} \right){\mathit{\boldsymbol{R}}_s}\left( {{f_i}} \right){\mathit{\boldsymbol{A}}^{\rm{H}}}\left( {{f_i}} \right) + {\delta ^2}\mathit{\boldsymbol{I}} \cdot \hfill \\ \end{gathered} $

(19)

式中：R_s(f_i)为频率为f_i信号的协方差矩阵, 对R(f_i)进行特征分解可得

$ \mathit{\boldsymbol{R}}({f_i}) = \sum\limits_{j = 1}^{{N_1} + {N_2}} {{\lambda _j}{\mathit{\boldsymbol{e}}_j}\mathit{\boldsymbol{e}}_j^{\rm{H}}} . $

(20)

由于信号之间相互独立, 则信号协方差矩阵的秩为K, 噪声协方差矩阵R_N=δ²I, 则存在如下的线性关系：

$ {\mathit{\boldsymbol{R}}_N}{\mathit{\boldsymbol{e}}_k} = \sum\limits_{j = 1}^K {\mathit{\boldsymbol{a}}\left( {{f_i}, {\theta _j}} \right){h_k}\left( j \right).} $

(21)

式中：1≤k≤K, e_k为特征矢量, h_k(i)为线性组合因子.当噪声为理想高斯白噪声时, 式(21) 可进一步简化为

$ {\mathit{\boldsymbol{e}}_k} = \sum\limits_{j = 1}^K {\mathit{\boldsymbol{a}}({\theta _j}){h_k}(j).} $

(22)

式(22) 表明无论信号是否相干, 大特征值对应的特征矢量为各信号方向矢量的一个线性组合因子.令v为频域协方差矩阵的大特征值之和, 则

$ \mathit{\boldsymbol{v}} = \sum\limits_{i = 1}^K {\mathit{\boldsymbol{a}}\left( {{f_k}, {\eta _i}} \right)\mathit{\boldsymbol{h}}\left( i \right)} {\rm{ }} = \mathit{\boldsymbol{A}}\left( {{f_k}} \right)\mathit{\boldsymbol{h}} \cdot $

(23)

构造式(23) 的稀疏模型, 即

$ \mathit{\boldsymbol{v}} = \mathit{\boldsymbol{BH}} \cdot $

(24)

式中：完备字典B的形式为

$ \mathit{\boldsymbol{B}} = \left[{\mathit{\boldsymbol{a}}\left( {\hat f_s^1} \right), \mathit{\boldsymbol{a}}\left( {\hat f_s^2} \right), \cdot \cdot \cdot, \mathit{\boldsymbol{a}}\left( {\hat f_s^K} \right)} \right]. $

(25)

其原子是由2.2节估计得到空间频率构造而成, 则配对函数可定义为

$ {\rm{min}}{\left\| {v\left( {{f_i}} \right) - \mathit{\boldsymbol{BH}}} \right\|_{\rm{2}}}{\rm{ + }}\varepsilon {\left\| \mathit{\boldsymbol{H}} \right\|_{\rm{2}}}. $

(26)

式中：ε为常量.利用CVX优化工具箱进行求解, 进而根据H中大非零元素的位置得到与信号频率$ {\hat f_i} $相匹配的空间频率$ \hat f_s^i $.

当多个信号同时入射时需要进行多次配对, 为提高配对的效率采用可变完备字典形式, 即在第一次配对完成后剔除完备字典B中与信号频率$ {\hat f_1} $相对应的空间频率$ \hat f_s^1 $相对应的原子形成新的完备字典B₁, 重复上述配对过程将所有的信号频率与空间频率实现配对.

当信号频率$ \left\{ {{{\hat f}_j}} \right\}_{j = 1}^K $与空间频率$ \left\{ {\hat f_s^i} \right\}_{j = 1}^K $正确配对后, 入射信号DOA计算可得

$ {\hat \theta _i} = {\rm{arcsin}}(c\hat f_s^{ij}/{\hat f_j}). $

(27)

式中：$ \hat f_s^{ij} $表示与信号频率$ {\hat f_j} $相匹配的空间频率.

3 仿真实验及性能分析

为了验证本文方法的有效性, 本节中主要通过理论分析和仿真实验对本文方法的估计性能进行分析, 探讨相关因素对估计性能的影响, 并与文献[5]方法进行比较.仿真中本文方法选用阵元数为6的二阶嵌套阵列作为接收阵列, 其中：N₁=N₂=3, 阵元最小间距d₁=1 m.以均方根误差(root mean square errors, RMSE)作为衡量估计性能优劣的指标, 空间频率与角度的均方根误差分别定义为

$ {\rm{RMSE}}\left( {{f_s}} \right){\rm{ = }}\sqrt {\frac{1}{{KM}}\sum\limits_{m = 1}^M {\sum\limits_{j = 1}^K {{{\left( {\hat f_s^j - f_s^j} \right)}^2}} } } . $

(28)

$ {\rm{RMSE}}\left( \theta \right) = \sqrt {\frac{1}{{KM}}\sum\limits_{m = 1}^M {\sum\limits_{j = 1}^K {{{\left( {{{\hat \theta }_j} - {\theta _j}} \right)}^2}} } } . $

(29)

式中：M为蒙特卡罗实验次数；$ \hat f_s^j $为第j个空间频率f_s^j的估计值；$ {\hat \theta _i} $为第i个信号角度θ_i的估计值.

1) 有效性分析分析

为验证本文所提方法可以实现对多于阵元数目的信号进行DOA估计, 假设9个等功率远场窄带信号入射到二级嵌套阵列上.信号的DOA分别为[-60.12°、-51.20°、-40.14°、-30.03°、20.31°、70.09°、80.15°、10.02°、35.07°], 中心频率分别为[115、153]MHz, 快拍数为2 048, 信噪比SNR=10 dB, 如图 3所示为信号频率估计结果；如图 4所示为空间频率估计结果，A为幅值；如图 5所示为DOA与频率的联合估计结果；如表 1所示为实际值与估计值对比；由图 3~5及表 1的实验结果可知, 当信号数多于阵元数时, 本文方法能够对宽频段内中心频率不同的多个窄带信号进行频率和DOA的联合估计且具有较高的估计精度.这是因为本文方法利用二级嵌套阵列有效地扩展了阵列孔径增加了阵列自由度, 因此能够在欠定条件下进行DOA估计.

2) 频率误差对估计性能的影响

由式(27) 可知, 信号频率估计误差会对DOA估计的精度产生影响, 为了分析信号频率估计误差对DOA估计的影响, 对f求偏导可得

$ \frac{{{\rm{d}}\theta }}{{{\rm{d}}f}} = \frac{{{f_s}c}}{{{f^2}\sqrt {1 - f_s^2{c^2}/{f^2}} }}. $

(30)

由式(30) 可知, 当空间频率f_s一定时, 信号频率f越高, 频率估计误差对DOA的估计精度影响越小.

假设有2个等功率远场窄带信号入射到二级嵌套阵列上, 中心频率分别为[100、150] MHz, 信号的角度在0°~90°范围内随机产生, 分别分析绝对频率估计误差$ \Delta f = |f - \hat f| $以及相对频率误差rf(rf=Δf/f)对DOA估计精度的影响, 每个频率估计误差下进行500次蒙特卡罗实验, 具体实验结果如图 6所示；

由图 6的实验结果可知, 角度均方根误差随着信号频率估计误差的增加而增加.当信号频率越高时, 绝对频率估计误差对角度均方根误差影响越小, 与前面的理论分析相一致.而当采样点足够多时, 现有的基于FFT变换的测频方法在采样点足够多时, 测频误差可达到几Hz, 基本上可以忽略其对估计性能的影响.而在同等相对频率误差条件下, 角度均方根误差不随信号频率的大小而变化.

3) 空间频率估计精度比较

由空间频率的定义以及上述的分析可知, DOA估计精度主要取决于空间频率的估计精度.要想获得高精度的DOA估计就必须对空间频率有高的估计能力.这里将本文方法与传统稀疏重构方法Lasso方法^[18]、MUSIC方法^[7]以及W-SPSF方法^[10]作比较.

假设4个远场窄带信号入射到二级嵌套阵列上, 空间频率分别为[-0.095 1、0.250 7、0.101 5、0.381 1], 对于Lasso方法以及W-SPSF方法, 以采样间隔Δ=0.001对空间频率-0.5~0.5范围进行均匀采样来构造完备字典, 对于MUSIC方法假设预先知道入射信号数目.每个信噪比、快拍数下进行500次蒙特卡罗实验, 如图 7所示为快拍数T=512时, 信噪比SNR从-5~15 dB, 步长为5 dB, 4种方法的空间频率均方根误差比较；如图 8所示为信噪比SNR=5 dB时, 快拍数从500~900, 步长为100, 4种方法的空间频率均方根误差比较.

由图 7、8的实验结果可知, 在一定的信噪比、快拍数范围内, 3种方法的空间频率均方根误差随着信噪比、快拍数的增加而减小.本文方法在低信噪比、少快拍数情况下估计性能要优于Lasso方法、MUSIC方法以及W-SPSF方法.

4) 空间频率分辨性能比较

假设2个相近的等功率远场信号入射到二级嵌套阵列上, 空间频率分别为[0.342 8、0.390 8], 信号频率为[103、175] MHz, 信噪比SNR=5 dB, 快拍数T=512.将本文方法与MUSIC方法以及Lasso方法的分辨性能进行比较, 如图 9所示为空间频率分辨性能比较结果；

由图 9的实验结果可知, 当2个信号的空间频率相近时, MUSIC方法以及Lasso方法均无法进行分辨, 而本文方法仍然能够估计出2个信号的空间频率, 分辨性能更优, 进而可以获得更高的DOA估计精度.

5) DOA估计精度比较

假设2个等功率远场窄带信号入射到阵列上, 空间频率分别为[0.338 3、0.101 5], 信号频率为[103、175] MHz.将本文方法与文献[5]SRSF方法以及文献[2]的DW-MUSIC方法进行比较, 本文方法选用二级嵌套阵列, SRSF以及DW-MUSIC方法选用均匀线阵, Δ为完备字典的采样间隔, DW-MUSIC方法的搜索间隔为0.001, 每个信噪比、快拍数下进行500次蒙特卡罗实验, 分析信噪比、快拍数对估计性能的影响.如图 10所示为快拍数T=1 024, 信噪比SNR从-5~15 dB, 步长为5 dB时, 角度均方根误差随信噪比变化关系, 如图 11所示为信噪比SNR=5 dB, 快拍数从500~900, 步长为100时, 角度均方根误差随快拍数变化关系；

由图 10和11的实验结果可知, SRSF方法的估计精度与完备字典的间隔有着直接的关系, 当完备字典的采样间隔较大时, 空间频率与完备字典失配导致SRSF方法的估计性能明显下降, 而本文方法由于在整个空频域-进行重构, 无需构造字典从而避免了模型不匹配对估计精度的影响.另外, 由于本文方法采用二级嵌套阵列, 形成具有更多阵元数的虚拟均匀阵列扩展了阵列孔径, DOA估计精度得到了相应的提高.因此, 本文方法的估计精度要高于SRSF方法以及DW-MUSIC方法.

4 结语

在现代复杂电子战信号环境中, 对宽频段内的多个窄带信号进行频率和DOA联合估计具有重要的现实意义.本文针对宽频段欠定DOA估计问题, 研究了一种基于二级嵌套阵列的DOA估计方法, 该方法利用空间频率将频率与DOA的二维参数转化成一维参数信息, 结合二级嵌套阵列的自由度扩展特性和空间频率的空域稀疏性来构建连续稀疏模型；然后利用原始对偶方法以及多项式求根获得空间频率的高精度估计；最后, 利用频域协方差矩阵的大特征矢量为各信号的方向矢量的线性组合的性质来建立配对函数解决了空间频率与频率的配对问题.理论分析和仿真实验验证了本文方法对宽频段内的多个窄带信号具有较高的估计精度和分辨能力, 可估计的信号数要远大于实际阵元数.