CN104021293A - 一种基于结构最小二乘法的联合到达角-频率估计方法 - Google Patents

Abstract

本发明提供了一种基于结构最小二乘法的联合DOA-频率估计方法，解决了传统ESPRIT算法联合估计DOA和频率时，对相干信号失效的问题。不同于传统ESPRIT算法中使用的抽样协方差矩阵，本发明采用前后向平均协方差矩阵将其替换，规避了ESPRIT算法在处理相干信号存有的风险，同时也克服了最优时域因子难以获取的问题。考虑信号子空间的误差，再运用SLS求解旋转不变方程，获取精度更高的信号子空间，完成对信号DOA和频率的鲁棒估计。

Description

一种基于结构最小二乘法的联合到达角-频率估计方法

技术领域

本发明涉及阵列信号处理技术领域，尤其涉及一种联合到达角-频率估计方法。

背景技术

阵列信号处理是信号处理领域中的一个重要分支，经过几十年的发展已日趋成熟并且在雷达、生物医疗、勘探及天文等多个军事和国民经济领域都有着广泛的应用。其工作原理是将多个传感器组成传感器阵列，并利用这一阵列对空间信号进行接收和处理，目的是抑制干扰和噪声，提取信号的有用信息。与一般的信号处理方式不同，阵列信号处理是通过布置在空间的传感器组接收信号，并且利用信号的空域特性来滤波及提取信息。因此，阵列信号处理也常被称为空域信号处理。此外，阵列信号处理有着灵活的波束控制、很强的抗干扰能力与极高的空间超分辨能力等优点，因而受到了众多学者的关注，其应用范围也不断地增大。

近年来，波达方向、频率联合估计在雷达、移动通信等领域有着重要的应用背景并引起了广泛的讨论。在过去几十年间，波达方向与频率参数联合估计的研究发展迅速，两者的精确估计能够在提升链接性能的情况下，更好的保证信道信息。

自频率—波数的最小方差(MV)方法问世以来，方位—频率联合估计的线性预测方法、多维MUSIC方法、最大似然方法、ESPRIT方法等算法相继问世。这些方法中，ML(最大似然估计)作为理论最优，与最小二乘法估计在高斯白噪声的情况下是等价的。虽然最大似然估计有更出色的统计性能，但其需要计算庞大的多维优化。类ESPRIT的算法虽然折中了估计精度和计算复杂度之间存在的矛盾，但它们在处理相干信号时，性能却衰减严重。为克服前述的缺点，对抽样数据做时—空预处理，进而采用ESPRIT算法对DOA和频率作联合参数估计，然而实际应用中最优时域因子m₀难以获取，因为m₀是关于抽样点数N的线性函数，即m₀＝(3N+2)＝5，随着抽样点数的增大，其复杂度将增至O(M³)。

在信号频率接近或相干的情况下，如何对DOA和频率进行联合估计时，最大限度的提高估计精度，正是本发明需要阐述和解决的。

发明内容

本发明的目的是解决传统ESPRIT算法联合估计DOA和频率时，对相干信号失效的问题。不同于传统ESPRIT算法中使用的抽样协方差矩阵，本发明采用前后向平均协方差矩阵将其替换，规避了ESPRIT算法在处理相干信号存有的风险，同时也克服了最优时域因子难以获取的问题。考虑信号子空间的误差，再运用SLS求解旋转不变方程，获取精度更高的信号子空间，完成对信号DOA和频率的鲁棒估计。

本发明通过如下技术方案实现：

一种基于结构最小二乘法的联合到达角-频率估计方法，包括以下步骤：

1)获取M全向阵元的均匀线阵的P个窄带信号经下变频为基带信号后的M×1个观测向量x(t)的抽样数据矩阵

其中，P<M，F为抽样率，m为抽样个抽样子序列个数，N为每个序列的抽样数；

2)计算抽样协方差矩阵 表示X_m的抽样协方差矩阵，然后计算前后向平均矩阵R：

$R=\frac{1}{2}(\hat{R}+\mathrm{\&Pi;}\hat{R}*\mathrm{\&Pi;}),$

其中，Π表示交换矩阵，其反对角线上元素均为1，其余为0；

3)计算信号子空间U_s，U_s由的对应于P个最大的特征值的P个特征向量构成；

4)使用最小二乘LS算法求解下列方程，得到Φ_θ，Φ_f的初始估计值 该初始估计值和U_s都用于步骤6)迭代的初始化：

$\begin{array}{cc}{J}_{\mathrm{\&theta;}}^{\&UpArrow;}{U}_s{\mathrm{\&Psi;}}_{\mathrm{\&theta;}}=J_{\mathrm{\&theta;}}^{\&DownArrow;}{U}_s& {\mathrm{\&Psi;}}_{\mathrm{\&theta;}}=\mathrm{T\&Theta;}{T}^{-1}\end{array}$

$J_f^{\&UpArrow;}{U}_s{\mathrm{\&Psi;}}_f=J_f^{\&DownArrow;}{U}_s,{\mathrm{\&Psi;}}_f=\mathrm{T\&Theta;}{T}^{-1},$

其中， $\mathrm{\&Theta;}=\mathrm{diag}\{e^{j2\mathrm{\&pi;}\sin {\mathrm{\&theta;}}_{1}d/\mathrm{\&lambda;}},...,e^{j2\mathrm{\&pi;}\sin {\mathrm{\&theta;}}_{P}d/\mathrm{\&lambda;}}\},$ $\mathrm{\&Phi;}=\mathrm{diag}\{e^{j2\mathrm{\&pi;}{f}_1/F},...,e^{j2\mathrm{\&pi;}{f}_P/F}\},$ λ是信号波长，阵元间距 $d=\mathrm{\&lambda;}/2,J_{\mathrm{\&theta;}}^{\&UpArrow;}=I_m\&CircleTimes;\left[\begin{array}{cc}{I}_{M-1}& 0_1\end{array}\right],J_{\mathrm{\&theta;}}^{\&DownArrow;}=I_m\&CircleTimes;\left[\begin{array}{cc}{0}_1& I_{M-1}\end{array}\right],$ $J_f^{\&UpArrow;}=\left[\begin{array}{cc}{I}_{m-1}& 0_1\end{array}\right]\&CircleTimes;I_m,$ $J_f^{\&UpArrow;}=\left[\begin{array}{cc}{0}_1& I_{m-1}\end{array}\right]\&CircleTimes;I_M,$ 表示Kronecker积，0₁是一个(M-1)×1的零向量， I_M是M×M的单位矩阵，T是一个P×P的非奇异矩阵；

5)定义如下两个误差矩阵：

$E_0=J_0^{\&UpArrow;}{\stackrel{\&OverBar;}{U}}_s{\stackrel{\&OverBar;}{\mathrm{\&Psi;}}}_{\mathrm{\&theta;}}-J_{\mathrm{\&theta;}}^{\&DownArrow;}{\stackrel{\&OverBar;}{U}}_s$

$E_f=J_f^{\&UpArrow;}{\stackrel{\&OverBar;}{U}}_s{\stackrel{\&OverBar;}{\mathrm{\&Psi;}}}_f-J_f^{\&DownArrow;}{\stackrel{\&OverBar;}{U}}_s;$

6)对于第K次迭代，K>1，用高斯牛顿法求解下式，得到ΔΨ_θ，k，ΔΨ_f，k和ΔΨ_s，k：

$\underset{{\mathrm{\&Delta;\&Psi;}}_{\mathrm{\&theta;},k},{\mathrm{\&Delta;\&Psi;}}_{f,k}}{\min }\left|\right|{H_k\&CenterDot;\left[\begin{array}{c}\mathrm{vec}\left\{{\mathrm{\&Delta;\&Psi;}}_{\mathrm{\&theta;},k}\right\}\\ \mathrm{vec}\left\{{\mathrm{\&Delta;\&Psi;}}_{f,k}\right\}\\ \mathrm{vec}\left\{{\mathrm{\&Delta;U}}_{s,k}\right\}\end{array}\right]+\left[\begin{array}{c}\mathrm{vec}\left\{E_{\mathrm{\&theta;},k}\right\}\\ \mathrm{vec}\left\{E_{f,k}\right\}\\ \mathrm{\&kappa;}\&CenterDot;\mathrm{vec}\left\{{\mathrm{\&Delta;U}}_k\right\}\end{array}\right]\left|\right|}^2,$

其中，vec{·}表示矢量运算,和的改良矩阵分别是 ${\stackrel{\&OverBar;}{U}}_{s,k}=U_{s,k-1}+{\mathrm{\&Delta;U}}_{s,k-1},{\stackrel{\&OverBar;}{\mathrm{\&Psi;}}}_{\mathrm{\&theta;},k}={\mathrm{\&Psi;}}_{\mathrm{\&theta;},k-1}+{\mathrm{\&Delta;\&Psi;}}_{\mathrm{\&theta;},k-1}$ 和 ${\stackrel{\&OverBar;}{\mathrm{\&Psi;}}}_{f,k}={\mathrm{\&Psi;}}_{f,k-1}+{\mathrm{\&Delta;\&Psi;}}_{f,k-1},$ E_θ，k和E_f，k表示第k次迭代的误差矩阵，为第k次迭代的信号子空间的估计误差矩阵，

$H_k=\left[\begin{array}{ccc}{I}_P\&CircleTimes;\left(J_{\mathrm{\&theta;}}^{\&UpArrow;}{U}_{s,k}\right)& 0& {\mathrm{\&Psi;}}_{\mathrm{\&theta;},k}^T\&CircleTimes;J_{\mathrm{\&theta;}}^{\&UpArrow;}-I_P\&CircleTimes;J_{\mathrm{\&theta;}}^{\&DownArrow;}\\ 0& I_P\&CircleTimes;\left(J_f^{\&UpArrow;}{U}_{s,k}\right)& {\mathrm{\&Psi;}}_{f,k}^T\&CircleTimes;J_f^{\&UpArrow;}-I_P\&CircleTimes;J_f^{\&DownArrow;}\\ 0& 0& \mathrm{\&kappa;}{I}_{\mathrm{mMP}}\end{array}\right],$

$E_{\mathrm{\&theta;},k+1}\&ap;E_{\mathrm{\&theta;},k}+J_{\mathrm{\&theta;}}^{\&UpArrow;}{U}_{s,k}{\mathrm{\&Delta;\&Psi;}}_{\mathrm{\&theta;},k}+J_{\mathrm{\&theta;}}^{\&UpArrow;}{\mathrm{\&Delta;U}}_{s,k}{\mathrm{\&Psi;}}_{\mathrm{\&theta;},k}-J_{\mathrm{\&theta;}}^{\&DownArrow;}{\mathrm{\&Delta;U}}_{s,k},$

$E_{f,k+1}\&ap;E_{f,k}+J_f^{\&UpArrow;}{U}_{s,k}{\mathrm{\&Delta;\&Psi;}}_{f,k}+J_f^{\&UpArrow;}{\mathrm{\&Delta;U}}_{s,k}{\mathrm{\&Psi;}}_{f,k}-J_f^{\&DownArrow;}{\mathrm{\&Delta;U}}_{s,k};$

7)判断是否满足：如果满足则迭代终止，得到Ψ_θ、Ψ_f的终值；如果不满足，K＝K+1,然后回到步骤5；

8)根据Φ_θ，Φ_f的终值，得到角度和频率的估计，具体为：由Ψ_θ、Ψ_f的终值求得Θ和Φ；对Θ和Φ作特征分解，得到如下DOA和频率的估计:

${\widehat{\mathrm{\&theta;}}}_i={\sin }^{-1}\left(\frac{\mathrm{\&lambda;}\&CenterDot;\&angle;\left({\mathrm{\&alpha;}}_i\right)}{2\mathrm{\&pi;d}}\right)$

${\hat{f}}_i=\frac{F\&CenterDot;\&angle;\left({\mathrm{\&beta;}}_i\right)}{2\mathrm{\&pi;}},i=1,\&CenterDot;\&CenterDot;\&CenterDot;,P$

上式中，∠表示角运算，α_i和β_i分别表示Ψ_θ和Ψ_f的第i个特征值。

本发明的基于结构最小二乘法的联合DOA-频率估计方法，解决了传统ESPRIT算法联合估计DOA和频率时，对相干信号失效的问题。不同于传统ESPRIT算法中使用的抽样协方差矩阵，本发明采用前后向平均协方差矩阵将其替换，规避了ESPRIT算法在处理相干信号存有的风险，同时也克服了最优时域因子难以获取的问题。考虑信号子空间的误差，再运用SLS求解旋转不变方程，获取精度更高的信号子空间，完成对信号DOA和频率的鲁棒估计。

附图说明

图1是本发明的基于结构最小二乘法的联合DOA-频率估计方法的流程图；

图2是本发明的联合DOA-频率估计方法与传统算法的RMSE和SNR的关系曲线图；

图3是本发明的联合DOA-频率估计方法与传统算法的RMSE和角度间隔的关系曲线图；

图4是本发明的联合DOA-频率估计方法与传统算法的RMSE和频率间隔的关系曲线图。

具体实施方式

下面结合附图说明及具体实施方式对本发明进一步说明。

考虑一M全向阵元的均匀线阵，感兴趣信号带宽的中心频率f_c。假设有P(P<M)个窄带信号(d_p(t)},中心频率为f_c+f_p，p＝1，…，P，方向为{θ₁，…，θ_P}的远场信号，作用在均匀线阵上。下变频为基带信号后，M×1个观测向量为

$\begin{array}{c}x\left(t\right)=\underset{i=1}{\overset{P}{\mathrm{\&Sigma;}}}a\left({\mathrm{\&theta;}}_i\right)d_i\left(t\right)e^{j2\mathrm{\&pi;}{f}_{i}t/F}+n\left(t\right)\\ =A{\mathrm{\&Phi;}}^{t}d\left(t\right)+n\left(t\right)\end{array}---\left(1\right)$

公式(1)中，F为抽样率，是源信号向量，n(t)是均值为零、方差为的加性高斯白噪声过程，I_M是M×M的单位矩阵，阵列流型A＝[a(θ₁)，…，a(θ_P)]中第p个导向矢量为

$a\left({\mathrm{\&theta;}}_p\right)={[1,e^{-j2\mathrm{\&pi;}\sin {\mathrm{\&theta;}}_{p}d/\mathrm{\&lambda;}},\&CenterDot;\&CenterDot;\&CenterDot;,e^{-j2\mathrm{\&pi;}(M-1)\sin {\mathrm{\&theta;}}_{p}d/\mathrm{\&lambda;}}]}^T---\left(2\right)$

以及

$\mathrm{\&Phi;}=\mathrm{diag}\{e^{j2\mathrm{\&pi;}{f}_1/F},...,e^{j2\mathrm{\&pi;}{f}_P/F}\}.---\left(3\right)$

其中，λ是信号波长，阵元间距d＝λ/2。

假设上述窄带信号的信道是块衰落的，在短暂的抽样间隔内认为{d_p(k)}保持不变，即

$d\left(t\right)\&ap;d(t+\frac{1}{F})\&ap;...\&ap;d(t+\frac{m-1}{F}).---\left(4\right)$

这意味着，对于同样的抽样周期，开头的m(m＜＜F)个抽样是近似相同的。收集m个抽样子序列，其中每个序列包含N个抽样。最终，数据矩阵结构如下所示

将(1)和(4)代入公式(5)中，得

X≈A_mD_m+N_m              (6)

其中

$\begin{array}{c}{A}_m={\left[\begin{array}{cccc}{\left(A\right)}^T& {\left(\mathrm{A\&Phi;}\right)}^T& \&CenterDot;\&CenterDot;\&CenterDot;& {\left({\mathrm{A\&Phi;}}^{m-1}\right)}^T\end{array}\right]}^T\\ D_m=\left[\begin{array}{cccc}d\left(0\right)& \mathrm{\&Phi;d}\left(\frac{1}{F}\right)& \&CenterDot;\&CenterDot;\&CenterDot;& {\mathrm{\&Phi;}}^{N-1}d\left(\frac{N-1}{F}\right)\end{array}\right]\\ N_m=\left[\begin{array}{cccc}n\left(0\right)& n\left(\frac{1}{F}\right)& \&CenterDot;\&CenterDot;\&CenterDot;& n\left(\frac{N-1}{F}\right)\end{array}\right]\end{array}---\left(7\right)$

使

这样，公式(6)可写为

公式(9)中，表示Khatri-Rao积。

令表示X_m的抽样协方差矩阵。当源信号是相关时，将很难估计。因此，本发明采用前后向平均矩阵替换即

$R=\frac{1}{2}(\hat{R}+\mathrm{\&Pi;}\hat{R}*\mathrm{\&Pi;})---\left(10\right)$

式中，Π表示交换矩阵，其反对角线上元素均为1，其余为0。

注意到对应于P个最大的特征值的P个特征向量形成了信号子空间U_s，也就是，span{U_s}＝span{A_m}。

定义如下两个选择矩阵：

$J_{\mathrm{\&theta;}}^{\&UpArrow;}=I_m\&CircleTimes;\left[\begin{array}{cc}{I}_{M-1}& 0_1\end{array}\right]---\left(11a\right)$

$J_{\mathrm{\&theta;}}^{\&DownArrow;}=I_m\&CircleTimes;\left[\begin{array}{cc}{0}_1& I_{M-1}\end{array}\right]---\left(11b\right)$

以及

$J_f^{\&UpArrow;}=\left[\begin{array}{cc}{I}_{m-1}& 0_1\end{array}\right]\&CircleTimes;I_M---\left(12a\right)$

$J_f^{\&DownArrow;}=\left[\begin{array}{cc}{0}_1& I_{m-1}\end{array}\right]\&CircleTimes;I_M---\left(12b\right)$

式中表示Kronecker积，0₁是一个(M-1)×1的零向量。然后，对于DOA和频率的旋转不变方程可表示为

$J_{\mathrm{\&theta;}}^{\&UpArrow;}{U}_s{\mathrm{\&Psi;}}_{\mathrm{\&theta;}}=J_{\mathrm{\&theta;}}^{\&DownArrow;}{U}_s$

$J_f^{\&UpArrow;}{U}_s{\mathrm{\&Psi;}}_f=J_f^{\&DownArrow;}{U}_s---\left(13\right)$

其中

Ψ_θ=TΘT^-1

Ψ_f＝TΦT^-1                      (14)

T是一个P×P的非奇异矩阵，同时

$\mathrm{\&Theta;}=\mathrm{diag}\{e^{j2\mathrm{\&pi;}\sin {\mathrm{\&theta;}}_{1}d/\mathrm{\&lambda;}},\&CenterDot;\&CenterDot;\&CenterDot;,e^{j2\mathrm{\&pi;}\sin {\mathrm{\&theta;}}_{P}d/\mathrm{\&lambda;}}\}.---\left(15\right)$

通过求解公式(14)以及对Θ和Φ作特征分解，得到如下DOA和频率的估计

${\widehat{\mathrm{\&theta;}}}_i={\sin }^{-1}\left(\frac{\mathrm{\&lambda;}\&CenterDot;\&angle;\left({\mathrm{\&alpha;}}_i\right)}{2\mathrm{\&pi;d}}\right)---\left(16\right)$

${\hat{f}}_i=\frac{F\&CenterDot;\&angle;\left({\mathrm{\&beta;}}_i\right)}{2\mathrm{\&pi;}},i=1,\&CenterDot;\&CenterDot;\&CenterDot;,P---\left(17\right)$

式中，∠表示角运算，α_i和β_i分别表示Ψ_θ和Ψ_f的第i个特征值。

由于公式(13)是一个高度结构化和超定的方程，随着其中重叠元素的增加，LS(Least square，最小二乘)解的性能将会下降，这是由于假设在和中是不存在误差的，只需要最小化和中的误差即可。然而，实际情况是在公式(13)中每一项都存在误差，即U_s，Ψ_θ和Ψ_f中均存在误差。它们的精确估计可以表示为 和为了最小化ΔU_s，ΔΨ_θ和ΔΨ_f，需要用到一个迭代最小化过程。

定义如下两个误差矩阵：

$E_{\mathrm{\&theta;}}=J_{\mathrm{\&theta;}}^{\&UpArrow;}{\stackrel{\&OverBar;}{U}}_s{\stackrel{\&OverBar;}{\mathrm{\&Psi;}}}_{\mathrm{\&theta;}}-J_{\mathrm{\&theta;}}^{\&DownArrow;}{\stackrel{\&OverBar;}{U}}_s---\left(18a\right)$

$E_f=J_f^{\&UpArrow;}{\stackrel{\&OverBar;}{U}}_s{\stackrel{\&OverBar;}{\mathrm{\&Psi;}}}_f-J_f^{\&DownArrow;}{\stackrel{\&OverBar;}{U}}_s---\left(18b\right)$

在第k次迭代，让 ${\stackrel{\&OverBar;}{U}}_{s,k}=U_{s,k-1}+\mathrm{\&Delta;}{U}_{s,k-1},$ ${\stackrel{\&OverBar;}{\mathrm{\&Psi;}}}_{\mathrm{\&theta;},k}={\mathrm{\&Psi;}}_{\mathrm{\&theta;},k-1}+\mathrm{\&Delta;}{\mathrm{\&Psi;}}_{\mathrm{\&theta;},k-1}$ 和分别表示和的改良矩阵。让E_θ，k和E_f，k表示第k次迭代的误差矩阵，因此在第(k+1)次迭代，忽略二次项ΔU_s，kΔΨ_θ，k和ΔU_s，kΔΨ_f，k，可得到

$E_{\mathrm{\&theta;},k+1}\&ap;E_{\mathrm{\&theta;},k}+J_{\mathrm{\&theta;}}^{\&UpArrow;}{U}_{s,k}\mathrm{\&Delta;}{\mathrm{\&Psi;}}_{\mathrm{\&theta;},k}+J_{\mathrm{\&theta;}}^{\&UpArrow;}\mathrm{\&Delta;}{U}_{s,k}{\mathrm{\&Psi;}}_{\mathrm{\&theta;},k}-J_{\mathrm{\&theta;}}^{\&DownArrow;}\mathrm{\&Delta;}{U}_{s,k}---\left(19a\right)$

$E_{f,k+1}\&ap;E_{f,k}+J_f^{\&UpArrow;}{U}_{s,k}\mathrm{\&Delta;}{\mathrm{\&Psi;}}_{f,k}+J_f^{\&UpArrow;}\mathrm{\&Delta;}{U}_{s,k}{\mathrm{\&Psi;}}_{f,k}-J_f^{\&DownArrow;}\mathrm{\&Delta;}{U}_{s,k}---\left(19b\right)$

从上式可得

$\begin{array}{c}\mathrm{vec}\left\{E_{\mathrm{\&theta;},k+1}\right\}\&ap;\mathrm{vec}\left\{E_{\mathrm{\&theta;},k}\right\}+[I_P\&CircleTimes;\left(J_{\mathrm{\&theta;}}^{\&UpArrow;}{U}_{s,k}\right)]\&times;\mathrm{vec}\left\{\mathrm{\&Delta;}{\mathrm{\&Psi;}}_{\mathrm{\&theta;},k}\right\}\\ +[{\mathrm{\&Psi;}}_{\mathrm{\&theta;},k}^T\&CircleTimes;J_{\mathrm{\&theta;}}^{\&UpArrow;}-I_P\&CircleTimes;J_{\mathrm{\&theta;}}^{\&DownArrow;}]\&times;\mathrm{vec}\left\{{\mathrm{\&Delta;\&Psi;}}_{s,k}\right\}\\ \mathrm{vec}\left\{E_{f,k+1}\right\}\&ap;\mathrm{vec}\left\{E_{f,k}\right\}+[I_P\&CircleTimes;\left(J_f^{\&UpArrow;}{U}_{s,k}\right)]\&times;\mathrm{vec}\left\{\mathrm{\&Delta;}{\mathrm{\&Psi;}}_{f,k}\right\}\\ +[{\mathrm{\&Psi;}}_{f,k}^T\&CircleTimes;J_f^{\&UpArrow;}-I_P\&CircleTimes;J_f^{\&DownArrow;}]\&times;\mathrm{vec}\left\{\mathrm{\&Delta;}{U}_{s,k}\right\}\end{array}---\left(20\right)$

式中，vec{·}表示矢量运算。同时，定义为第k次迭代的信号子空间的估计误差矩阵。将(20)式整理为矩阵形式得到如下SLS问题

$\underset{\mathrm{\&Delta;}{\mathrm{\&Psi;}}_{\mathrm{\&theta;},k}\mathrm{\&Delta;}{\mathrm{\&Psi;}}_{f,k}}{\min }\left|\right|{H_k\&CenterDot;\left[\begin{array}{c}\mathrm{vec}\left\{\mathrm{\&Delta;}{\mathrm{\&Psi;}}_{\mathrm{\&theta;},k}\right\}\\ \mathrm{vec}\left\{{\mathrm{\&Delta;\&Psi;}}_{f,k}\right\}\\ \mathrm{vec}\left\{\mathrm{\&Delta;}{U}_{s,k}\right\}\end{array}\right]+\left[\begin{array}{c}\mathrm{vec}\left\{E_{\mathrm{\&theta;},k}\right\}\\ \mathrm{vec}\left\{E_{f,k}\right\}\\ \mathrm{\&kappa;}\&CenterDot;\mathrm{vec}\left\{\mathrm{\&Delta;}{U}_k\right\}\end{array}\right]\left|\right|}^2---\left(21\right)$

其中k>1是一个用户自定义的参数，其作用是保持ΔU_s的元素比E_θ和E_f中的大。

以及

$H_k=\left[\begin{array}{ccc}{I}_p\&CircleTimes;\left(J_{\mathrm{\&theta;}}^{\&UpArrow;}{U}_{s,k}\right)& 0& {\mathrm{\&Psi;}}_{\mathrm{\&theta;},k}^T\&CircleTimes;J_{\mathrm{\&theta;}}^{\&UpArrow;}-I_P\&CircleTimes;J_{\mathrm{\&theta;}}^{\&DownArrow;}\\ 0& I_P\&CircleTimes;\left(J_f^{\&UpArrow;}{U}_{s,k}\right)& {\mathrm{\&Psi;}}_{f,k}^T\&CircleTimes;J_f^{\&UpArrow;}-I_P\&CircleTimes;J_f^{\&DownArrow;}\\ 0& 0& \mathrm{\&kappa;}{I}_{\mathrm{mMP}}\end{array}\right]---\left(22\right)$

使用最小二乘解和U_s分别作为Ψ_θ ，Ψ_f和U_s的初始估计值，当K＝1时，ΔΨ_θ，k，ΔΨ_f，k和ΔU_s，k可以是任意的一个初始值。迭代终止的条件为

$\min \{{\left|\right|{\mathrm{\&Delta;\&Psi;}}_{\mathrm{\&theta;},k}\left|\right|}_F^2,{\left|\right|{\mathrm{\&Delta;\&Psi;}}_{f,k}\left|\right|}_F^2,{\left|\right|{\mathrm{\&Delta;u}}_{s,k}\left|\right|}_F^2\}\&le;\mathrm{\&epsiv;}---\left(23\right)$

式中，ε>0是一个预定义的小常量。当算法在第k次迭代收敛，那么和分别表示Ψ_θ和Ψ_f的最终估计值。DOA和频率的估计值最终可由公式(16)和(17)计算得出。

附图1是本发明提供的方法的算法流程图。基于结构最小二乘法的联合DOA-频率估计方法具体为：

下面，通过实验数据验证本发明的方法。考虑一阵元数M＝7的均匀线阵，阵列间距为d＝λ/2。噪声是均值为零、方差为1的白高斯噪声，且与信号互不相关。实验在RMSE(Root Mean Square Error，均方根误差)意义下，比较ESPIRT算法和JAFE算法。在JAFE算法中，L≈0.4M。在后续的实验中，设定L＝3。此外，抽样数m＝3，快拍数N＝64，收敛值ε＝10^-7，同时假设信号的个数是已知的，所有的仿真结果均由1000次蒙特卡洛实验获得。

实验1  RMSE与信噪比的关系

在这个仿真中，展示了RMSE的性能随信噪比所呈现的变化。固定噪声功率的情况下，从-6dB到30dB来变化信号功率。三个功率相同的信号作用到阵列上，入射角为θ₁＝10°，θ₂＝19°和θ₃＝30°。信号的中心频率分别为f₁＝2MHz，f₂＝2.06MHz和f₃＝2.2MHz。对于JAFE算法，时域和空域平滑因子都接近于3。从附图2中可以观察到，在图示的SNR变化范围内，本发明提供的方法(附图中的proposed)的性能，无论是DOA估计还是频率估计都优于ESPRIT和JAFE算法。

实验2  RMSE与角度间隔的关系

在该实验中，设法比较了DOA和频率的估计误差随角度间隔所呈现的变化趋势。此例中，只考虑两个信号，信噪比固定为SNR＝0dB，两个信号的中心频率分别为f₁＝2MHz和f₂＝2.15MHz，第一个信号的DOA设定为θ₁＝0°，第二个信号的DOA为θ₂＝0°+Δθ，其中Δθ∈[0°，16°]。从附图3可以看出，当角度间隔小于6°时，本发明提供的方法在DOA和频率估计这两方面的性能均优于其他两种算法。

实验3  RMSE与频率间隔的关系

在此仿真中呈现的是RMSE随频率间隔的变化。设定两个信号的DOA为θ₁＝0°和θ₂＝6°。第一个信号的频率为f₁＝2MHz，另外一个信号的频率为f₂＝(2+Δf)，Δf∈[0，200]kHz。在频率间隔很小时，如Δf<100kHz，两个信号是时域相关的。由于前后向平均以及SLS算法，很容易从附图4中看出本发明提供的方法的性能优于其他两种方法。当Δf足够大，JAFE稍逊ESPRIT方法，而本发明提供的方法和ESPRIT贴合到一起。

本发明提供的方法是一种基于SLS-ESPRIT的联合DOA-频率估计算法，与传统的ESPRIT算法不能很好处理相干信号不同，提及的算法相干信号的处理也有着相当的精度，在估计性能上有明显提高。

1.  对相干信号的处理。具体实施步骤如下：

首先，在抽样频率为F时，抽取个数为m(m＜＜F)的信号，此时认为抽样间隔内的信号是不变的，在快拍为N时，得到的阵列输出数据矩阵为其次，设定为X_m的抽样协方差矩阵，由于信号是相关的，用前后向平均矩阵替换估计性能不佳的抽样协方差矩阵降低算法对相干信号的敏感度，这一步所需的计算量为O(m²M²N)，加上后续对R的EVD运算所产生的计算量O(m³M³)，这一过程的计算量为O(m²M²N+m³M³)。

2.考虑了子空间里存在的误差，用搜索迭代求解信号子空间。具体实施步骤如下：

首先，在使用ESPRIT算法联合估计DOA和频率时，也考虑了子空间里存在的误差ΔU_s，构建如公式(21)所示的SLS问题。然后，利用SLS对不变方程进行搜索迭代，每一步迭代的计算量为O(10m³M³P³)，在 和的最小值均满足小于预定义值ε后，迭代将会终止。最终，信号子空间的估计值为DOA和频率的估计值也分别由 ${\mathrm{\&Psi;}}_{\mathrm{\&theta;}}^{\mathrm{SLS}}={\mathrm{\&Psi;}}_{\mathrm{\&theta;},k}+\mathrm{\&Delta;}{\mathrm{\&Psi;}}_{\mathrm{\&theta;},k}$ 和 ${\mathrm{\&Psi;}}_f^{\mathrm{SLS}}={\mathrm{\&Psi;}}_{f,k}+\mathrm{\&Delta;}{\mathrm{\&Psi;}}_{f,k}$ 给出，这个过程产生的计算量为O(10Km³M³P³)。

通过仿真分析，本发明提供的方法估计性能高于传统的ESPRIT算法和JAFE算法，不但完成了对DOA和频率的联合估计，同时也得到较为精确的信号子空间。

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。

Claims (2)

1.一种基于结构最小二乘法的联合到达角-频率估计方法，其特征在于：

所述方法包括以下步骤：

1)获取M全向阵元的均匀线阵的P个窄带信号经下变频为基带信号后的M×1个观测向量x(t)的抽样数据矩阵

其中，P<M，F为抽样率，m为抽样个抽样子序列个数，N为每个序列的抽样数；

2)计算抽样协方差矩阵 表示X_m的抽样协方差矩阵，然后计算前后向平均矩阵R：

$R=\frac{1}{2}(\hat{R}+\mathrm{\&Pi;}\widehat{R*}\mathrm{\&Pi;}),$

其中，Π表示交换矩阵，其反对角线上元素均为1，其余为0；

3)计算信号子空间U_s，U_s由的对应于P个最大的特征值的P个特征向量构成；

4)使用最小二乘LS算法求解下列方程，得到Φ_θ，Φ_f的初始估计值 该初始估计值和U_s都用于步骤6)迭代的初始化：

$J_{\mathrm{\&theta;}}^{\&UpArrow;}{U}_s{\mathrm{\&Psi;}}_{\mathrm{\&theta;}}=J_{\mathrm{\&theta;}}^{\&DownArrow;}{U}_s,{\mathrm{\&Psi;}}_{\mathrm{\&theta;}}=\mathrm{T\&Theta;}{T}^{-1}$

$J_f^{\&UpArrow;}{U}_s{\mathrm{\&Psi;}}_f=J_f^{\&DownArrow;}{U}_s,{\mathrm{\&Psi;}}_f=\mathrm{T\&Theta;}{T}^{-1},$

其中， $\mathrm{\&Theta;}=\mathrm{diag}\{e^{j2\mathrm{\&pi;}\sin {\mathrm{\&theta;}}_{1}d/\mathrm{\&lambda;}},\&CenterDot;\&CenterDot;\&CenterDot;,e^{j2\mathrm{\&pi;}\sin {\mathrm{\&theta;}}_{P}d/\mathrm{\&lambda;}}\},$ $\mathrm{\&Phi;}=\mathrm{diag}\{e^{j2\mathrm{\&pi;}{f}_1/F},\&CenterDot;\&CenterDot;\&CenterDot;,e^{j2\mathrm{\&pi;}{f}_P/F}\},$ λ是信号波长，阵元间距 $d=\mathrm{\&lambda;}/2,J_{\mathrm{\&theta;}}^{\&UpArrow;}=I_m\&CircleTimes;\left[\begin{array}{cc}{I}_{M-1}& 0_1\end{array}\right],J_{\mathrm{\&theta;}}^{\&DownArrow;}=I_m\&CircleTimes;\left[\begin{array}{cc}{0}_1& I_{M-1}\end{array}\right],$ $J_f^{\&UpArrow;}=\left[\begin{array}{cc}{I}_{m-1}& 0_1\end{array}\right]\&CircleTimes;I_M,J_f^{\&DownArrow;}=\left[\begin{array}{cc}{0}_1& I_{m-1}\end{array}\right]\&CircleTimes;I_M,$ 表示Kronecker积，0₁是一个(M-1)×1的零向量，I_M是M×M的单位矩阵，T是一个P×P的非奇异矩阵；

5)定义如下两个误差矩阵：

$E_{\mathrm{\&theta;}}=J_{\mathrm{\&theta;}}^{\&UpArrow;}{\stackrel{\&OverBar;}{U}}_s{\stackrel{\&OverBar;}{\mathrm{\&Psi;}}}_{\mathrm{\&theta;}}-J_{\mathrm{\&theta;}}^{\&DownArrow;}{\stackrel{\&OverBar;}{U}}_s$

$E_f=J_f^{\&UpArrow;}{\stackrel{\&OverBar;}{U}}_s{\stackrel{\&OverBar;}{\mathrm{\&Psi;}}}_f-J_f^{\&DownArrow;}{\stackrel{\&OverBar;}{U}}_s;$

6)对于第K次迭代，K>1，用高斯牛顿法求解下式，得到ΔΨ_θ，k，ΔΨ_f，k

和ΔU_s，k：

$\underset{\mathrm{\&Delta;}{\mathrm{\&Psi;}}_{\mathrm{\&theta;},k},\mathrm{\&Delta;}{\mathrm{\&Psi;}}_{f,k}}{\min }{\left|\right|H_k\&CenterDot;\left[\begin{array}{c}\mathrm{vec}\left\{\mathrm{\&Delta;}{\mathrm{\&psi;}}_{\mathrm{\&theta;},k}\right\}\\ \mathrm{vec}\left\{\mathrm{\&Delta;}{\mathrm{\&Psi;}}_{f,k}\right\}\\ \mathrm{vec}\left\{\mathrm{\&Delta;}{U}_{s,k}\right\}\end{array}\right]+\left[\begin{array}{c}\mathrm{vec}\left\{E_{\mathrm{\&theta;},k}\right\}\\ \mathrm{vec}\left\{E_{f,k}\right\}\\ k\&CenterDot;\mathrm{vec}\left\{{\mathrm{\&Delta;U}}_k\right\}\end{array}\right]\left|\right|}^2,$

其中，vec{·}表示矢量运算,和的改良矩阵分别是 ${\stackrel{\&OverBar;}{U}}_{s,k}=U_{s,k-1}+{\mathrm{\&Delta;U}}_{s,k-1},{\stackrel{\&OverBar;}{\mathrm{\&Psi;}}}_{\mathrm{\&theta;},k}={\mathrm{\&Psi;}}_{\mathrm{\&theta;},k-1}+{\mathrm{\&Delta;\&Psi;}}_{\mathrm{\&theta;},k-1}$ 和 ${\stackrel{\&OverBar;}{\mathrm{\&Psi;}}}_{f,k}={\mathrm{\&Psi;}}_{f,k-1}+{\mathrm{\&Delta;\&Psi;}}_{f,k-1},E_{\mathrm{\&theta;},k}$ 和E_f，k表示第k次迭代的误差矩阵，为第k次迭代的信号子空间的估计误差矩阵，

$H_k=\left[\begin{array}{ccc}{I}_P\&CircleTimes;\left(J_{\mathrm{\&theta;}}^{\&UpArrow;}{U}_{s,k}\right)& 0& {\mathrm{\&Psi;}}_{\mathrm{\&theta;},k}^T\&CircleTimes;J_{\mathrm{\&theta;}}^{\&UpArrow;}-I_P\&CircleTimes;J_{\mathrm{\&theta;}}^{\&DownArrow;}\\ 0& I_P\&CircleTimes;\left(J_f^{\&UpArrow;}{U}_{s,k}\right)& {\mathrm{\&Psi;}}_{f,k}^T\&CircleTimes;J_f^{\&UpArrow;}-I_P\&CircleTimes;J_f^{\&DownArrow;}\\ 0& 0& {\mathrm{kI}}_{\mathrm{mMP}}\end{array}\right]$

$E_{\mathrm{\&theta;},k}\&ap;E_{\mathrm{\&theta;},k-1}+J_{\mathrm{\&theta;}}^{\&UpArrow;}{U}_{s,k-1}{\mathrm{\&Delta;\&Psi;}}_{\mathrm{\&theta;},k-1}+J_{\mathrm{\&theta;}}^{\&UpArrow;}{\mathrm{\&Delta;U}}_{s,k-1}{\mathrm{\&Psi;}}_{\mathrm{\&theta;},k-1}-J_{\mathrm{\&theta;}}^{\&DownArrow;}{\mathrm{\&Delta;U}}_{s,k-1}$

$E_{f,k}\&ap;E_{f,k-1}+J_f^{\&UpArrow;}{U}_{s,k-1}{\mathrm{\&Delta;\&Psi;}}_{f,k-1}+J_f^{\&UpArrow;}{\mathrm{\&Delta;U}}_{s,k-1}{\mathrm{\&Psi;}}_{\mathrm{f\&theta;},k-1}-J_f^{\&DownArrow;}{\mathrm{\&Delta;U}}_{s,k-1};$

7)判断是否满足： $\min \{{\left|\right|{\mathrm{\&Delta;\&Psi;}}_{\mathrm{\&theta;},k}\left|\right|}_F^2,{\left|\right|{\mathrm{\&Delta;\&Psi;}}_{f,k}\left|\right|}_F^2,{\left|\right|{\mathrm{\&Delta;U}}_{s,k}\left|\right|}_F^2\}\&le;\&Element;$ 如果满足则迭代终止，得到Ψ_θ、Ψ_f的终值；如果不满足，K＝K+1,然后回到步骤5；

8)根据Φ_θ，Φ_f的终值，得到角度和频率的估计，具体为：由Ψ_θ、Ψ_f的终值求得Θ和Φ；对Θ和Φ作特征分解，得到如下DOA和频率的估计:

${\widehat{\mathrm{\&theta;}}}_i={\sin }^{-1}\left(\frac{\mathrm{\&lambda;}\&CenterDot;\&angle;\left({\mathrm{\&alpha;}}_i\right)}{2\mathrm{\&pi;d}}\right)$

${\hat{f}}_i=\frac{F\&CenterDot;\&angle;\left({\mathrm{\&beta;}}_i\right)}{2\mathrm{\&pi;}},i=1,...,P$

上式中，∠表示角运算，α_i和β_i分别表示Ψ_θ和Ψ_f的第i个特征值。

2.根据权利要求1所述的估计方法，其特征在于：所述ΔΨ_θ，k，ΔΨ_f，k和ΔU_s，k的迭代初始值是任给的一个值。