开放科学(资源服务)标志码(OSID):

0 概述

基于传感器阵列的波达方向(Direction of Arrival，DOA)估计被广泛应用于生产生活中^[1]，例如雷达、声呐、地震勘探、导航、声源跟踪等^[2-4]。现有DOA估计算法大多假设被测信号是窄带且噪声服从高斯分布^[5]，并利用多重信号分类^[6-7](Multiple Signal Classification，MUSIC)算法进行估计，因此，在处理宽带信号源时的性能受自然界声波、地震波等因素^[8-10]限制。

宽带信号DOA估计算法主要分为非相干信号子空间算法^[11-12]和相干信号子空间算法^[13-14]。非相干信号子空间算法仅分辨非相干信号，并且在低信噪比条件下估计性能不佳；相干信号子空间算法可用于相干信号，在低信噪比条件下估计准确度高，但其聚焦过程运算量较大且性能受波达角度预估计偏差的影响，波达角度预估计结果较差时会导致性能严重下降^[15]。文献[16]提出一种投影子空间正交性测试(TOPS)算法，利用多个频点子空间的正交性实现宽带信号DOA估计，但其估计精度不高且易出现伪峰。文献[17]提出修正的TOPS算法，利用信号子空间投影有效剔除伪峰，但其在低信噪比条件下估计性能不佳。文献[18]根据宽带阵列导向矢量在法线方向上的频率一致性，提出基于频域时延补偿的DOA估计算法，该算法性能优，但计算量较大。文献[19]采用接收到的数据构造聚焦矩阵以避免DOA预估计，但估计性能依赖聚焦频率的选取。文献[20]利用阵型中理想的低秩Toeplitz结构，在快拍数不足的情况下实现DOA估计，但该方法需要求解半定规划问题且计算复杂度较高。文献[21]提出信号子空间聚焦(Focusing Signal Subspace，FSS)算法，将参考频率的信号子空间特征向量和其他频率的信号子空间特征向量结合到Frobenius范数的约束中，实现聚焦矩阵的构造，该算法分辨率高且均方根误差低，无需进行初步的DOA估计，但在短快拍条件下性能不佳。

本文提出一种改进的信号子空间聚焦算法MFSS。根据子频带波长间隔与半波长的匹配度选取最佳参考频率及子频带，减少聚焦过程的运算量，同时将协方差矩阵平均处理为Hankel矩阵并进行奇异值分解重构，降低噪声及短快拍对协方差矩阵的影响。在此基础上，利用信号子空间聚焦构造最终的聚焦协方差矩阵，并通过Root-正交传播算子实现DOA估计。

1 宽带信号模型

均匀线阵模型如图 1所示。假设有$ P $个相互独立的远场宽带相干信号以角度$ {\theta }_{1}, {\theta }_{2}, \cdots , {\theta }_{P} $入射到$ M $个阵元的均匀线阵上，阵元间隔$ d $为信号最高频率对应的半波长。

第$ m $个阵元接收的信号如式(1)所示：

$ {x_m}\left( t \right) = \sum\limits_{p = 1}^P {{s_p}} (t - {\tau _{mp}}) + {n_m}\left( t \right),m = {\rm{1}},{\rm{2}}, \cdots ,M$

(1)

其中：$ {s}_{p}\left(t\right) $为阵元接收的第$ p $个信号；$ {n}_{m}\left(t\right) $为第$ m $个阵元接收的噪声；$ {\tau }_{mp} $为第$ p $个信号到达第$ m $个阵元较到达第1个阵元的延迟。$ {\tau }_{mp} $如式(2)所示：

$ {\tau }_{mp}=\frac{1}{c}(m-1)d\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{p} $

(2)

其中：$ c $为信号传播速度。

均匀线阵模型对接收数据$ {x}_{m}\left(t\right) $做离散傅里叶变换，将其划分为$ J $个子带，快拍数为$ K $，如式(3)所示：

$ {X}_{m}\left({f}_{j}\right)=\sum\limits _{p=1}^{P}{S}_{p}\left({f}_{j}\right){e}^{-j2\pi {\tau }_{mp}}+{N}_{m}\left({f}_{j}\right) $

(3)

矩阵形式的频域阵列信号接收模型如式(4)所示：

$ X\left( {{f_j}} \right) = \mathit{\boldsymbol{A}}({f_j},\theta )S\left( {{f_j}} \right) + N\left( {{f_j}} \right) $

(4)

其中：$ \mathit{\boldsymbol{A}}({f}_{j}, \theta ) $为阵列流型矩阵。列向量$ \mathit{\boldsymbol{\alpha }}({f}_{j}, {\theta }_{p}) $如式(5)所示：

$ \mathit{\boldsymbol{\alpha }}({f}_{j}, {\theta }_{p})=[1\mathrm{ }, {\mathrm{e}}^{-\mathrm{j}2\mathrm{\pi }{f}_{j}{\tau }_{2p}}, \mathrm{ }\cdots \mathrm{ }, {\mathrm{e}}^{-\mathrm{j}2\mathrm{\pi }{f}_{j}{\tau }_{Mp}}] $

(5)

$ \mathit{\boldsymbol{X}}\left({f}_{j}\right) $的协方差矩阵如式(6)所示：

$ \begin{aligned}\mathit{\boldsymbol{R}}\left({f}_{j}\right)= & E\left[\mathit{\boldsymbol{X}}\right({f}_{j}\left)\mathit{\boldsymbol{X}}^{\mathrm{H}}\right({f}_{j}\left)\right]=\\ & \mathit{\boldsymbol{A}}({f}_{j}, \theta )\mathit{\boldsymbol{R}}_{s}\left({f}_{j}\right)\mathit{\boldsymbol{A}}^{\mathrm{H}}({f}_{j}, \theta )+\mathit{\boldsymbol{R}}_{n}\left({f}_{j}\right)\end{aligned} $

(6)

其中：$ \mathit{\boldsymbol{R}}_{s}\left({f}_{j}\right) $、$ \mathit{\boldsymbol{R}}_{n}\left({f}_{j}\right) $分别为频率上信号与噪声的协方差矩阵。因此，$ \mathit{\boldsymbol{R}}_{s}\left({f}_{j}\right) $如式(7)所示：

$ \mathit{\boldsymbol{R}}_{s}\left({f}_{j}\right)=E\left[\mathit{\boldsymbol{S}}\right({f}_{j}\left)\mathit{\boldsymbol{S}}^{\mathrm{H}}\right({f}_{j}\left)\right] $

(7)

2 基于MFSS的宽带信号DOA估计 2.1 FSS算法

FSS算法的原理是结合参考频率和其他频率信号子空间的特征向量，在Frobenius范数约束下构造聚焦协方差矩阵，采用MUSIC算法实现DOA估计。

FSS算法对$ \mathit{\boldsymbol{R}}\left({f}_{j}\right) $进行特征分解，如式(8)所示：

$ \mathit{\boldsymbol{R}}\left({f}_{j}\right)=\mathit{\boldsymbol{U}}\left({f}_{j}\right)\mathrm{\Sigma }\left({f}_{j}\right)\mathit{\boldsymbol{U}}^{\mathrm{H}}\left({f}_{j}\right) $

(8)

其中：$ \mathit{\boldsymbol{U}}\left({f}_{j}\right) $为$ \mathit{\boldsymbol{R}}\left({f}_{j}\right) $的$ M\times M $维的特征向量矩阵；$ \mathit{\boldsymbol{U}}\left({f}_{j}\right)=\left[{e}_{1}\right({f}_{j}), {e}_{2}({f}_{j}), \cdots , {e}_{M}({f}_{j}\left)\right] $。由于信号与噪声相互独立，则信号子空间$ \mathit{\boldsymbol{U}}_{\mathrm{s}}\left({f}_{j}\right) $、噪声子空间$ \mathit{\boldsymbol{U}}_{n}\left({f}_{j}\right) $分别如式(9)、式(10)所示：

$ \mathit{\boldsymbol{U}}_{\mathrm{s}}\left({f}_{j}\right)=\left[{e}_{1}\right({f}_{j}), {e}_{2}({f}_{j}), \cdots , {e}_{P}({f}_{j}\left)\right] $

(9)

$ \mathit{\boldsymbol{U}}_{n}\left({f}_{j}\right)=\left[{e}_{P+1}\right({f}_{j}), {e}_{P+2}({f}_{j}), \cdots , {e}_{M}({f}_{j}\left)\right] $

(10)

FSS算法构造非奇异矩阵$ \mathit{\boldsymbol{T}}\left({f}_{j}\right) $，将其作为聚焦矩阵，$ \mathit{\boldsymbol{U}} $_s(f_j)如式(11)所示：

$ \mathit{\boldsymbol{U}}_{\mathrm{s}}\left({f}_{j}\right)=\mathit{\boldsymbol{T}}\left({f}_{j}\right)\mathit{\boldsymbol{U}}_{\mathrm{s}}\left({f}_{j}\right) $

(11)

为获得最小的聚焦误差，FSS算法使用Frobenius范数约束信号子空间，如式(12)所示：

$ \underset{\mathit{T}\left({f}_{j}\right)}{\mathrm{m}\mathrm{i}\mathrm{n}}\left|\right|\mathit{\boldsymbol{U}}_{\mathrm{s}}\left({f}_{0}\right)-\mathit{\boldsymbol{T}}\left({f}_{j}\right)\mathit{\boldsymbol{U}}_{\mathrm{s}}\left({f}_{j}\right)|{|}_{\mathrm{F}}^{2} $

(12)

其中：$ \mathit{\boldsymbol{T}}\left({f}_{j}\right) $为Hermitian矩阵。因此，聚焦矩阵$ \mathit{\boldsymbol{T}}\left({f}_{j}\right) $满足式(13)：

$ \mathit{\boldsymbol{T}}^{\mathrm{H}}\left({f}_{j}\right)\mathit{\boldsymbol{T}}\left({f}_{j}\right)=\mathit{\boldsymbol{I}} $

(13)

将式(12)重写为：

$ \begin{array}{l}\underset{\mathit{T}\left({f}_{\xi }\right)}{\mathrm{m}\mathrm{i}\mathrm{n}}\left|\right|\mathit{\boldsymbol{U}}_{\mathrm{s}}\left({f}_{0}\right)-\mathit{\boldsymbol{T}}\left({f}_{j}\right)\mathit{\boldsymbol{U}}_{\mathrm{s}}\left({f}_{j}\right)|{|}_{\mathrm{F}}^{2}=\\ \mathrm{t}\mathrm{r}\left(\mathit{\boldsymbol{U}}_{\mathrm{s}}\right({f}_{0}\left)\mathit{\boldsymbol{U}}_{\mathrm{s}}^{\mathrm{H}}\right({f}_{0}\left)\right)+\mathrm{t}\mathrm{r}\left(\mathit{\boldsymbol{U}}_{\mathrm{s}}\right({f}_{j}\left)\mathit{\boldsymbol{U}}_{\mathrm{s}}^{\mathrm{H}}\right({f}_{j}\left)\right)-\\ 2\mathrm{R}\mathrm{e}\left[\mathrm{t}\mathrm{r}\right(\mathit{\boldsymbol{U}}_{\mathrm{s}}\left({f}_{0}\right)\mathit{\boldsymbol{U}}_{\mathrm{s}}^{\mathrm{H}}\left({f}_{j}\right)\mathit{\boldsymbol{T}}^{\mathrm{H}}\left({f}_{j}\right)\left)\right]\end{array} $

(14)

式(14)前两项是固定值，若要获得最小的聚焦误差，应使$ \mathrm{R}\mathrm{e}\left[\mathrm{t}\mathrm{r}\right(\mathit{\boldsymbol{U}}_{\mathrm{s}}\left({f}_{0}\right)\mathit{\boldsymbol{U}}_{\mathrm{s}}^{\mathrm{H}}\left({f}_{j}\right)\mathit{\boldsymbol{T}}^{\mathrm{H}}\left({f}_{j}\right)\left)\right] $取最大值。

矩阵$ \mathit{\boldsymbol{C}} $如式(15)所示：

$ \mathit{\boldsymbol{C}} = {\mathit{\boldsymbol{U}}_{\rm{s}}}\left( {{f_0}} \right)\mathit{\boldsymbol{U}}_{\rm{s}}^{\rm{H}}\left( {{f_j}} \right) = \mathit{\boldsymbol{\tilde U}}\sum {{\mathit{\boldsymbol{\tilde V}}^{\rm{H}}}} $

(15)

其中：$ \mathit{\boldsymbol{\tilde U}} $和$ \mathit{\boldsymbol{\tilde V}} $分别为矩阵$ \mathit{\boldsymbol{C}} $奇异值分解后的左右奇异分量。则$ \mathrm{R}\mathrm{e}\left[\mathrm{t}\mathrm{r}\right(\mathit{\boldsymbol{U}}_{\mathrm{s}}\left({f}_{0}\right)\mathit{\boldsymbol{U}}_{\mathrm{s}}^{\mathrm{H}}\left({f}_{j}\right)\mathit{\boldsymbol{T}}^{\mathrm{H}}\left({f}_{j}\right)\left)\right] $满足：

$ \begin{array}{l} {\rm{Re}}[{\rm{tr}}({\mathit{\boldsymbol{U}}_{\rm{s}}}\left( {{f_0}} \right)\mathit{\boldsymbol{U}}_{\rm{s}}^{\rm{H}}\left( {{f_j}} \right){\mathit{\boldsymbol{T}}^{\rm{H}}}\left( {{f_j}} \right))] \le \\ |{\rm{tr}}({\mathit{\boldsymbol{U}}_{\rm{s}}}\left( {{f_0}} \right)\mathit{\boldsymbol{U}}_{\rm{s}}^{\rm{H}}\left( {{f_j}} \right){\mathit{\boldsymbol{T}}^{\rm{H}}}\left( {{f_j}} \right))| = \\ |{\rm{tr}}(\Sigma {{\mathit{\boldsymbol{\tilde V}}}^{\rm{H}}}{\mathit{\boldsymbol{T}}^{\rm{H}}}\left( {{f_j}} \right)\mathit{\boldsymbol{\tilde U}})| \end{array}$

(16)

令$ \mathit{\boldsymbol{Z}} = {\mathit{\boldsymbol{\tilde V}}^{\rm{H}}}{\mathit{\boldsymbol{T}}^{\rm{H}}}\left( {{f_j}} \right)\mathit{\boldsymbol{\tilde U}} $，$ \mathit{\boldsymbol{Z}} $是一个酉矩阵，$ \left|{Z}_{ii}\right|\le 1, i=\mathrm{1, 2}, \cdots , M $，带入式(16)中得：

$ \begin{array}{*{20}{l}} {{\rm{Re}}[{\rm{tr}}({\mathit{\boldsymbol{U}}_{\rm{s}}}\left( {{f_0}} \right)\mathit{\boldsymbol{U}}_{\rm{s}}^{\rm{H}}\left( {{f_j}} \right){\mathit{\boldsymbol{T}}^{\rm{H}}}\left( {{f_j}} \right))] \le }\\ {|{\rm{tr}}(\sum \mathit{\boldsymbol{Z}} )| \le \sum\limits_{i = 1}^M {\left| {{Z_{ii}}} \right|} {\sigma _{ii}} \le \sum\limits_{i = 1}^M {{\sigma _{ii}}} } \end{array} $

(17)

当$ \left|{Z}_{ii}\right|=1, i=\mathrm{1, 2}, \cdots , M $时，式(17)可取最大值，此时：

$ \begin{array}{*{20}{l}} {\mathit{\boldsymbol{Z}} = {{\mathit{\boldsymbol{\tilde V}}}^{\rm{H}}}{\mathit{\boldsymbol{T}}^{\rm{H}}}\left( {{f_j}} \right)\mathit{\boldsymbol{\tilde U}} = {{\mathit{\boldsymbol{\tilde V}}}^{\rm{H}}}\left( {\mathit{\boldsymbol{\tilde V}}{{\mathit{\boldsymbol{\tilde U}}}^{\rm{H}}}} \right)\mathit{\boldsymbol{\tilde U}} = I}\\ {{\rm{objectto}}\left| {{Z_{ii}}} \right| = 1,i = {\rm{1}},{\rm{2}}, \cdots ,M} \end{array} $

(18)

$ \mathit{\boldsymbol{T}}\left({f}_{j}\right) $如式(19)所示：

$ \mathit{\boldsymbol{T}}\left( {{f_j}} \right) = \mathit{\boldsymbol{\tilde U}}{\mathit{\boldsymbol{\tilde V}}^{\rm{H}}} $

(19)

则频率$ {f}_{j} $处的聚焦协方差矩阵如式(20)所示：

$ {\mathit{\boldsymbol{R}}_j}\left( {{f_0}} \right) = \mathit{\boldsymbol{T}}\left( {{f_j}} \right)\mathit{\boldsymbol{R}}\left( {{f_j}} \right){\mathit{\boldsymbol{T}}^{\rm{H}}}\left( {{f_j}} \right) $

(20)

最终的协方差矩阵$ {\mathit{\boldsymbol{R}}_{\rm{F}}} $是子频带的聚焦协方差矩阵的均值，如式(21)所示：

$ {\mathit{\boldsymbol{R}}_{\rm{F}}} = \frac{1}{J}\left( {\sum\limits_{j = 1}^J {{R_j}} \left( {{f_0}} \right)} \right)$

(21)

最后，通过MUSIC算法实现DOA估计。

2.2 MFSS算法

MFSS算法是将宽带信号分为J段并进行傅里叶变换，通过计算子频带波长间隔与半波长的插值选定参考频率，并筛选出3个子频带，将子频带的协方差矩阵处理为Hankel矩阵，采用奇异值分解去噪并重构协方差矩阵，利用信号子空间聚焦法构造聚焦协方差矩阵并通过Root-正交传播算子得到DOA估计值。MFSS算法流程如图 2所示。

2.2.1 参考频率及子频带选取

对于$ J $个子频带，MFSS算法计算频率$ {f}_{c}=\frac{c}{\lambda }=\frac{c}{2d} $，其中$ c $为速度，阵元间隔$ d=\frac{\lambda }{2} $。在子带$ \mathrm{1, 2}, \cdots , J $中，MFSS算法搜索波长间隔最接近正整数$ {j}_{0}=\frac{J}{{f}_{s}}\mathrm{m}\mathrm{o}\mathrm{d}\left(\frac{{f}_{c}}{f}\right) $的子带，其中$ \mathrm{m}\mathrm{o}\mathrm{d} $表示余数，$ {f}_{s} $为信号的采样频率。由于每个子频带信号的波长是不同的，因此MFSS算法需选取一个子频带与该波长间隔最为匹配。

MFSS算法选取该子频带的中心频率$ {f}_{0} $作为参考频率，使用该子频带及临近两条子频带进行后续的计算估计，临近的两条子频带频率为$ {f}_{\xi } $，$ \xi =\mathrm{1, 2} $。

2.2.2 协方差矩阵重构

在实际工程中，协方差矩阵由有限数量采样条件下获得接收数据的平均值构成，且噪声形式复杂。这些因素都会导致特征分解时信号子空间与噪声子空间划分模糊。针对该问题，本文采用Hankel矩阵奇异值分解法对协方差矩阵进行重构。

本文对选取子频带的数据协方差矩阵$ \mathit{\boldsymbol{R}}\left({f}_{\xi }\right) $次对角线及平行于次对角线直线上的元素进行平均处理，如式(22)所示：

$ \left\{\begin{array}{l}{\tilde{r}}_{mn}=\tilde{r}(m+n+1)=\frac{1}{M-k}\sum\limits _{i=1}^{M-k}{r}_{i, M-k+1-i}, 0\le k < M\\ {\tilde{r}}_{mn}=\tilde{r}(M-m-n+1)=\frac{1}{M+k}\sum\limits _{i=1-k}^{M+k}{r}_{i, M+k+1-i}, k < 0\end{array}\right. $

(22)

其中：$ k $为对角线；$ k=0 $为次对角线；$ 0\le k < M $为次对角线上方的对角线；-$ k $为次对角线下方的对角线；$ {r}_{mn} $和$ {\tilde{r}}_{mn} $分别为$ \mathit{\boldsymbol{R}}\left( {{f_\xi }} \right) $和平均处理后矩阵中第$ m $行、第$ n $列的元素。Hankel矩阵$ \mathit{\boldsymbol{\tilde R}}\left( {{f_\xi }} \right) $如式(23)所示：

$ \mathit{\boldsymbol{\tilde R}}\left( {{f_\xi }} \right) = \left[ {\begin{array}{*{20}{c}} {{r_{M - 1}}}&{{r_{M - 2}}}& \cdots &{{r_0}}\\ {{r_{M - 2}}}&{{r_{M - 3}}}& \cdots &{{r_{ - 1}}}\\ ?&?&{}&?\\ {{r_0}}&{{r_{ - 1}}}& \cdots &{{r_{ - M + 1}}} \end{array}} \right] $

(23)

本文对Hankel矩阵$ \mathit{\boldsymbol{\tilde R}}\left( {{f_\xi }} \right) $进行奇异值分解，如式(24)所示：

$ \mathit{\boldsymbol{\tilde R}}\left( {{f_\xi }} \right) = \mathit{\boldsymbol{UD}}{\mathit{\boldsymbol{V}}^{\rm{T}}}$

(24)

奇异值矩阵$ \mathit{\boldsymbol{D}} $中的大奇异值对应信号分量，而小奇异值对应噪声分量，因此，保留奇异值矩阵$ \mathit{\boldsymbol{D}}$的前$ P $个奇异值，后$ M-P $个奇异值为0，以实现去噪的目的，处理后的矩阵为$ \mathit{\boldsymbol{D}}_{P} $。奇异值矩阵利用$ \mathit{\boldsymbol{U}} $、$ {\mathit{\boldsymbol{V}}}^{\mathrm{T}} $和$ {\mathit{\boldsymbol{D}}}_{P} $构建新的矩阵$ {\mathit{\boldsymbol{R}}}_{P}\left({f}_{\xi }\right) $，如式(25)所示：

$ {\mathit{\boldsymbol{R}}}_{P}\left({f}_{\xi }\right)=\mathit{\boldsymbol{U}}{\mathit{\boldsymbol{D}}}_{P}{\mathit{\boldsymbol{V}}}^{\mathrm{T}} $

(25)

2.2.3 信号子空间聚焦

本文对$ {\mathit{\boldsymbol{R}}}_{P}\left({f}_{\xi }\right) $及$ \mathit{\boldsymbol{R}}\left({f}_{0}\right) $特征进行分解，求出$ {f}_{\xi } $及$ {f}_{0} $处的信号子空间$ {\mathit{\boldsymbol{U}}}_{\mathrm{s}}\left({f}_{\xi }\right) $及$ {\mathit{\boldsymbol{U}}}_{\mathrm{s}}\left({f}_{0}\right) $，对$ {\mathit{\boldsymbol{U}}}_{\mathrm{s}}\left({f}_{0}\right){\mathit{\boldsymbol{U}}}_{\mathrm{s}}^{\mathrm{H}}\left({f}_{\xi }\right) $的乘积进行奇异值分解，如式(26)所示：

$ {\mathit{\boldsymbol{U}}}_{\mathrm{s}}\left({f}_{0}\right){\mathit{\boldsymbol{U}}}_{\mathrm{s}}^{\mathrm{H}}\left({f}_{\xi }\right)={\tilde{\mathit{\boldsymbol{U}}}}_{\mathrm{F}}\mathrm{\Sigma }{\tilde{\mathit{\boldsymbol{V}}}}_{\bf{F}}^{\bf{H}} $

(26)

聚焦矩阵构造为$ \mathit{\boldsymbol{T}}\left({f}_{\xi }\right)={\tilde{\mathit{\boldsymbol{U}}}}_{\mathrm{F}}{\tilde{\mathit{\boldsymbol{V}}}}_{\bf{F}}^{\bf{H}} $，根据式(20)求得频率$ {f}_{\xi } $处的聚焦协方差矩阵$ {\mathit{\boldsymbol{R}}}_{\xi }\left({f}_{0}\right) $。最终的协方差矩阵$ {\mathit{\boldsymbol{R}}}_{\mathrm{M}\mathrm{F}} $如式(27)所示：

$ {\mathit{\boldsymbol{R}}}_{\mathrm{M}\mathrm{F}}=\frac{1}{3}\left(\mathit{\boldsymbol{R}}\left({f}_{0}\right)+\sum\limits _{\xi =1}^{2}{\mathit{\boldsymbol{R}}}_{\xi }\left({f}_{0}\right)\right) $

(27)

2.2.4 Root-正交传播算子

本文对协方差矩阵$ {\mathit{\boldsymbol{R}}}_{\mathrm{M}\mathrm{F}} $进行分块，如式(28)所示：

$ {\mathit{\boldsymbol{R}}}_{\mathrm{M}\mathrm{F}}=\left[\begin{array}{l}{R}_{1}\\ {R}_{2}\end{array}\right] $

(28)

其中：$ {R}_{1} $和$ {R}_{2} $分别为矩阵$ {\mathit{\boldsymbol{R}}}_{\mathrm{M}\mathrm{F}} $的前$ P $行及后$ M-P $行，假设$ {R}_{1} $满秩，则$ {R}_{2} $为$ {R}_{1} $的线性组合，故存在$ (M-P)\times P $维变换矩阵$ \hat{\mathit{\boldsymbol{P}}} $，如式(29)所示：

$ {R}_{2}={\hat{\mathit{\boldsymbol{P}}}}^{\mathrm{H}}{R}_{1} $

(29)

在有限快拍的情况下，式(29)并不成立，变换矩阵$ \hat{\mathit{\boldsymbol{P}}} $可以通过最小化代价函数求出，如式(30)、式(31)所示：

$ J\left(\hat{\mathit{\boldsymbol{P}}}\right)=\left|\right|{\hat{\mathit{\boldsymbol{P}}}}^{\mathrm{H}}{R}_{1}-{R}_{2}|{|}^{2} $

(30)

$ \hat{\mathit{\boldsymbol{P}}}=({\mathit{\boldsymbol{R}}}_{1}^{\mathrm{H}}{R}_{1}{)}^{-1}{\mathit{\boldsymbol{R}}}_{1}^{\mathrm{H}}{R}_{2} $

(31)

变换矩阵$ \hat{\mathit{\boldsymbol{P}}} $即为传播算子，根据传播算子获得噪声子空间估计$ \hat{\mathit{\boldsymbol{Q}}} $，如式(32)所示：

$ \hat{\mathit{\boldsymbol{Q}}}=[\hat{\mathit{\boldsymbol{P}}}, -{\mathit{\boldsymbol{I}}}_{M-P}] $

(32)

其中：$ {\mathit{\boldsymbol{I}}}_{M-P} $为$ M-P $维单位阵。Root-正交传播算子对$ \hat{\mathit{\boldsymbol{Q}}} $进行正交化，如式(33)所示：

$ {\hat{\mathit{\boldsymbol{Q}}}}_{0}=(\hat{\mathit{\boldsymbol{Q}}}{\hat{\mathit{\boldsymbol{Q}}}}^{\mathrm{H}}{)}^{-\frac{1}{2}}\hat{\mathit{\boldsymbol{Q}}} $

(33)

定义多项式为：

$ f\left(z\right)={z}^{M-1}\mathit{\boldsymbol{p}}\left(z\right){\hat{\mathit{\boldsymbol{Q}}}}_{0}{\hat{\mathit{\boldsymbol{Q}}}}_{0}^{\mathrm{H}}p\left(z\right) $

(34)

其中：$ z={\mathrm{e}}^{jw} $；$ \mathit{\boldsymbol{p}}\left(z\right)=[1, z, \cdots , {z}^{M-1}{]}^{\mathrm{T}} $。根据式(34)计算得到$ P $个接近于单位元上的根，对于均匀线阵，单位元上的根可通过式(35)求解：

$ {\theta }_{i}=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{i}\mathrm{n}\left(\frac{\lambda }{2\mathrm{\pi }d}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}\left\{{\hat{z}}_{i}\right\}\right), i=\mathrm{1, 2}, \cdots , P $

(35)

MFSS算法的步骤主要分为：1)对阵列接收到的宽带信号数据分段，并进行离散傅里叶变换；2)选取波长间隔最接近$ {j}_{0}=\frac{J}{{f}_{s}}\mathrm{m}\mathrm{o}\mathrm{d}\left(\frac{{f}_{c}}{f}\right) $子频带的中心频率作为参考频率，保留该子频带及其临近两条子频带；3)求得各频点处的协方差矩阵$ \mathit{\boldsymbol{R}}\left({f}_{\xi }\right) $，根据式(22)将其处理为Hankel矩阵并利用奇异值分解重构；4)对$ {\mathit{\boldsymbol{R}}}_{P}\left({f}_{\xi }\right) $及$ \mathit{\boldsymbol{R}}\left({f}_{0}\right) $进行特征分解，根据各频点处的信号子空间，构造聚焦矩阵$ \mathit{\boldsymbol{T}}\left({f}_{\xi }\right)={\tilde{\mathit{\boldsymbol{U}}}}_{\mathrm{F}}{{\tilde{\mathit{\boldsymbol{V}}}}_{\mathrm{F}}}^{\mathrm{H}} $，计算聚焦协方差矩阵$ {\mathit{\boldsymbol{R}}}_{\xi }\left({f}_{0}\right) $，与$ \mathit{\boldsymbol{R}}\left({f}_{0}\right) $进行平均操作得到最终的协方差矩阵$ {\mathit{\boldsymbol{R}}}_{\mathrm{M}\mathrm{F}} $；5)对协方差矩阵$ {\mathit{\boldsymbol{R}}}_{\mathrm{M}\mathrm{F}} $分块，根据式(32)求得噪声子空间估计$ \hat{\mathit{\boldsymbol{Q}}} $，计算式(35)的$ P $个接近于单位元上的根，得到信号的DOA估计。

3 算法复杂度分析

为评估MFSS算法的实用价值，本文分析MFSS、FSS、MTOPS^[17]及LR-MUSIC算法^[20]的复杂度。MFSS算法的时间复杂度主要由以下5项构成：1)选取参考频率及子频带时间复杂度$ O\left(J\right) $；2)构造协方差矩阵时间复杂度$ O\left(K{M}^{2}\right) $；3)重构协方差矩阵$ {\mathit{\boldsymbol{R}}}_{P} $时间复杂度$ O({M}^{2}+{M}^{3}) $；4)聚焦获得最终协方差矩阵$ {\mathit{\boldsymbol{R}}}_{MF} $时间复杂度$ O\left({M}^{3}\right) $；5)构造传播算子多项式求根时间复杂度$ O\left(P{M}^{2}\right) $。因此，MFSS算法总体时间复杂度约为$ O(J+P{M}^{2}+{M}^{3}) $。FSS算法需要$ J $个子带以实现聚焦协方差矩阵的构造，并采用MUSIC算法构造空间谱，时间复杂度约为$ O(JK{M}^{2}+J{M}^{3}) $。MTOPS算法构造协方差矩阵并特征分解的时间复杂度约为$ O(JK{M}^{2}+J{M}^{3}) $，计算判决矩阵的时间复杂度为$ O\left(2PM\right(M-P)+{P}^{2}(M-P\left)\right(J-1\left)\right) $，求解矩阵迹时间复杂度为$ O\left(2{P}^{3}\right(J-1\left)\right) $，MTOPS算法总体时间复杂度较高。LR-MUSIC算法构造聚焦矩阵时间复杂度为$ O\left(J{M}^{3}\right) $，计算聚焦协方差矩阵时间复杂度为$ O\left(JK{M}^{2}\right) $，此时的时间复杂度已接近FSS算法，在后续计算中最优Toeplitz矩阵通过半正定规划求出，计算量大，且总体复杂度远大于MFSS算法。与其他算法相比，MFSS算法在时间复杂度上具有较大优势，更利于工程应用。

4 仿真实验与分析

本文对MFSS、FSS、MTOPS、LR-MUSIC算法进行仿真对比，以验证MFSS算法在短快拍情况下的有效性。假设信号源数目已知，阵元数$ M=8 $，信号为频率$ 200~400 $ Hz的远场相干宽带信号，阵元间距为信号中心频率所对应的半波长，噪声为相互独立且与信号无关的高斯白噪声。DFT点数为128，子频带数目$ J=20 $。LR-MUSIC算法是通过CVX工具包解决SDP问题。

4.1 实验1

假定2个远场相干宽带信号入射角分别为$ 60° $和$ 80° $，快拍数为50，本文进行500次蒙特卡洛实验。DOA估计的均方根误差(R_RMSE)如式(36)所示：

$ {R}_{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}=\frac{1}{P}\sum\limits _{p=1}^{P}\sqrt{\frac{1}{C}\sum\limits _{c=1}^{C}({\tilde{\theta }}_{p, c}-{\theta }_{p}{)}^{2}} $

(36)

其中：$ C $为蒙特卡洛实验的数量。

在不同信噪比时4种算法的均方根误差对比如图 3所示，其信噪比范围从-20 dB以间隔2 dB升至5 dB。从图 3可以看出，在短快拍条件下，4种算法的均方根误差均随信噪比的增加而逐渐减小。在整个信噪比范围内，MFSS算法的均方根误差始终低于其他3种算法。因此，在短快拍低信噪比条件下，MFSS算法的估计误差最小。

估计成功率为$ {D}_{P}=\frac{{N}_{i}}{N} $，$ {N}_{i} $为进行$ C $次蒙特卡洛实验中满足$ |{\tilde{\theta }}_{i}-{\theta }_{i}| < 1° $的次数。在不同信噪比时4种算法的估计成功率对比如图 4所示。从图 4可以看出，MFSS、MTOPS、LR-MUSIC算法的估计成功率随信噪比的增大趋近1，MFSS算法将成功估计的信噪比门限降低至-4 dB。当信噪比低于-4 dB时，MFSS算法的估计成功率明显高于MTOPS、LR-MUSIC算法。当信噪比为0 dB时，这3种算法能完全成功估计两个信号的波达方向，在短快拍情况下，FSS算法始终无法完全成功估计。因此，MFSS算法在低信噪比条件下具有更高的估计成功率。

4.2 实验2

本文考虑2个入射角为$ {\theta }_{1}=60° $、$ {\theta }_{2}=60°+\mathrm{\Delta }\theta $的宽带相干信号，信噪比为5 dB，快拍数为50。若$ |{\tilde{\theta }}_{i}-{\theta }_{i}| < \frac{\mathrm{\Delta }\theta }{2} $，$ i=\mathrm{1, 2} $，则算法可以成功分辨两目标。不同角度间隔下4种算法的估计成功率如图 5所示。

从图 5可以看出，随着宽带信号源角度间隔增大，4种算法的分辨性能明显提高，当两信号源角度间隔为$ 5° $时，MFSS和MTOPS算法的估计成功率达到$ 100\mathrm{\%} $，而FSS和LR-MUSIC算法仍无法分辨两信号源。在角度间隔小于$ 5° $时，MFSS算法的估计成功率高于MTOPS算法，MFSS算法的成功估计角度间隔门限较FSS算法降低了$ 3° $。因此，MFSS算法在短快拍条件下能够分辨角度间隔更小的信号源。

4.3 实验3

本文考虑2个独立宽带相干信号源的入射角为$ {\theta }_{1}=60° $、$ {\theta }_{2}=80° $，信噪比为5 dB。在不同快拍数下4种算法的均方根误差对比如图 6所示。从图 6可以看出，快拍数从5上升至50时，MFSS算法的均方根误差最小。

在不同快拍数下4种算法的估计成功率对比如图 7所示。从图 7可以看出，在快拍数大于40时，MFSS、TOPS和LR-MUSIC算法的估计成功率较接近，逐渐趋于1。在短快拍条件下，MFSS算法的估计成功率始终高于其他3种算法。因此，在同条件下，MFSS算法具有更高的估计精度。

4.4 实验4

本文分别考虑2个独立宽带相干信号(入射角$ {\theta }_{1}=60° $、$ {\theta }_{2}=80° $)和3个独立宽带相干信号(入射角$ {\theta }_{1}=40° $、$ {\theta }_{2}=60° $、$ {\theta }_{3}=80° $)的情况，进行500次蒙特卡洛实验，得出算法的平均运算时间。

不同算法的运算时间对比如表 1所示。从表 1可以看出，MFSS算法的运算时间最短，且远小于MTOPS和LR-MUSIC算法的运算时间。相比FSS算法，在信源数为3时，MFSS算法的平均运算时间降低了21.14%，更具实用性。

5 结束语

本文提出一种无需角度预估计的信号子空间聚焦算法MFSS，利用奇异值的分布规律减少快拍数及噪声对估计性能的影响，通过波长间隔与阵元间距的匹配度选取最佳参考频点及子频带，降低运算量。仿真结果表明，MFSS算法在短快拍条件下能够有效提高估计精度。后续将提高算法在复杂噪声情况下的估计性能，使其适用于实际无线传感网络定位环境。