CN103344940B - 低复杂度的doa估计方法及系统 - Google Patents

Abstract

本发明提供了一种低复杂度的DOA估计方法及系统，在低复杂度的DOA估计方法中对样本数据做酉变换并将变换得到的数据分解为两部分，计算的自相关矩阵以及与的互相关矩阵，并分别提取出它们的实部。本发明的有益效果是只需要计算两个实数域的子样本协方差矩阵，通过Nystr？m方法构建出信号子空间，避免了构造整个样本协方差矩阵以及其特征值分解，从而进一步降低了复杂度。

Description

低复杂度的DOA估计方法及系统

技术领域

本发明涉及信号处理领域，一种低复杂度的DOA估计方法及系统。

背景技术

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

在阵列信号处理领域，最重要的两个研究方向是自适应滤波和空间谱估计，其中自适应滤波技术先于空间谱估计产生，而且其应用在工程系统中已十分广泛。然而，对于空间谱估计虽然在近30年中得到了快速的发展，相关研究内容十分广泛，但其工程应用系统却不多见。这里，空间谱是阵列信号处理领域中的一个重要概念，表示信号在空间各个方向上的能量分布。如果可以获取信号的空间谱，就能估计信号的波达方向（direction-of-arrival,DOA），因此，空间谱估计也常被称作“DOA估计”或者“超分辨谱估计”。所谓DOA是指空间信号到达天线阵列的方向，即信号的入射角。

在传统的超分辨信号波达方向估计算法中，以MUSIC算法和ESPRIT算法为代表的子空间方法受到最为广泛的研究。这类算法的特点是将阵列接收数据的协方差矩阵分解为相互正交的信号子空间和噪声子空间，并利用信号和噪声子空间直接的关系来估计信号的到达角。然而，此类算法均需对协方差矩阵做特征值分解，所需的运算量较大，约为其中M是阵元个数。随着技术的发展，天线尺寸也越来越庞大，子空间类算法的实际应用也受到了巨大的阻碍。

发明内容

为了解决现有技术中的问题，本发明提供了一种低复杂度的DOA估计方法。

本发明提供了一种低复杂度的DOA估计方法，包括如下步骤：

A.对样本数据做酉变换并将变换得到的数据分解为两部分：和其中，K为用户自定义参数且P≤K≤min{M,N}，P为信源个数；

B.计算的自相关矩阵以及与的互相关矩阵，并分别提取出它们的实部，分别用和表示；

C.定义一个新的矩阵 $Z=\left[\begin{array}{c}{\tilde{R}}_{11}\\ {\tilde{R}}_{21}\end{array}\right]{\tilde{R}}_{11}^{-1/2},$ 对Z^HZ做特征值分解同时构造矩阵Π=ZU_Z，信号子空间由Π的前P列组成，即U_S=Π(:,1:P)。

作为本发明的进一步改进，在所述步骤A中包括如下步骤：

A1.考虑一M阵元的均匀线阵，令X表示其接收的样本数据矩阵，令 $\tilde{X}=Q_M^{H}X;$

A2.分解数据样本为 $\tilde{X}=\left[\begin{array}{c}{\tilde{X}}_1\\ {\tilde{X}}_2\end{array}\right],$ 其中 这里N为样本数，K为用户定义参数满足P≤K≤min{M,N}，P为信源个数。

作为本发明的进一步改进，在所述步骤B中包括如下步骤：

B1.分别计算的自相关矩阵R₁₁与和的互相关矩阵R₂₁；

B2.令和分别表示R₁₁与R₂₁的实部。

作为本发明的进一步改进，在所述步骤C中包括如下步骤：

C1.利用公式 $Z=\left[\begin{array}{c}{\tilde{R}}_{11}\\ {\tilde{R}}_{21}\end{array}\right]{\tilde{R}}_{11}^{-1/2}$ 定义矩阵Z，令表示Z^HZ的特征值分解，其中Λ_Z=diag{λ₁,…,λ_K}表示由特征值组成的对角矩阵，并且所有特征值按降序排列，即λ₁≥…≥λ_K，U_Z=[u₁,…,u_K]表示相应的特征矢量矩阵，这里u_i表示第i个特征值λ_i所对应的特征矢量；

C2.通过公式Π=ZU_Z构造矩阵Π，则信号子空间由Π的前P列构成，即U_S=Π(:,1:P)。

作为本发明的进一步改进，该DOA估计方法还包括如下步骤：

D.按如下公式定义稀疏酉矩阵

$Q_M=\left\{\begin{array}{cc}\frac{1}{\sqrt{2}}\left[\begin{array}{cc}{I}_l& j{I}_l\\ J_l& -j{J}_l\end{array}\right]& \mathrm{for\; M}=2l\\ \frac{1}{\sqrt{2}}\left[\begin{array}{ccc}{I}_l& 0_{l\×1}& j{I}_l\\ 0_{l\×1}^T& \sqrt{2}& 0_{l\×1}^T\\ J_l& 0_{l\×1}& -j{J}_l\end{array}\right]& \mathrm{for\; M}=2l+1\end{array}\right.$

式中，J_l为l×l的交换矩阵，其反对角线上的元素为1其余均为0；

E.定义矩阵 $K_1=\mathrm{Re}\left\{Q_{M-1}^H{J}_{s2}{J}_M\right\},$ $K_2=\mathrm{Im}\left\{Q_{M-1}^H{J}_{s2}{J}_M\right\},$ 其中J_s1和J_s1分别定义为：J_s1=[I_(M-1)×M,0_(M-1)×1]，J_s2=[0_(M-1)×1,I_(M-1)×M]，表示选择矩阵；

F.利用步骤C2中得到的信号子空间结合步骤E中的K₁与K₂，令其中表示矩阵伪逆运算；

G.对Ψ做特征值分解，得其中ψ_i表示特征根，e_i表示相对于的特征矢量，则第i个信号的到达角表示为

${\widehat{\mathrm{\θ}}}_i={\sin }^{-1}\left(\frac{\mathrm{\λ}\·{\tan }^{-1}\left({\widehat{\mathrm{\ψ}}}_i\right)}{\mathrm{\πd}}\right),i=\mathrm{1,2},\·\·\·,P$

其中，和分别表示θ_i与ψ_i的估计值。

本发明还提供了一种低复杂度的DOA估计系统，包括：

复数据实数化模块，用于令X表示其接收的样本数据矩阵，令 $\tilde{X}=Q_M^{H}X;$

分解数据样本模块，用于分解数据样本为 $\tilde{X}=\left[\begin{array}{c}{\tilde{X}}_1\\ {\tilde{X}}_2\end{array}\right],$ 其中 这里N为样本数，K为用户定义参数满足P≤K≤min{M,N}，P为信源个数；

估计协方差矩阵模块，用于分别计算的自相关矩阵R₁₁与和的互相关矩阵R₂₁，令和分别表示R₁₁与R₂₁的实部。

作为本发明的进一步改进，该DOA估计系统还包括：

定义矩阵Z模块，用于定义一个新的矩阵 $Z=\left[\begin{array}{c}{\tilde{R}}_{11}\\ {\tilde{R}}_{21}\end{array}\right]{\tilde{R}}_{11}^{-1/2};$

对Z^HZ做特征值分解模块，用于令表示Z^HZ的特征值分解，

其中Λ_Z=diag{λ₁,…,λ_K}表示由特征值组成的对角矩阵，并且所有特征值按降序排列，即λ₁≥…≥λ_K，U_Z=[u₁,…,u_K]表示相应的特征矢量矩阵，这里u_i表示第i个特征值λ_i所对应的特征矢量；

构造矩阵模块，用于通过公式Π=ZU_Z构造矩阵Π；

得到信号子空间模块，用于信号子空间由Π的前P列构成，即U_S=Π(:,1:P)。

作为本发明的进一步改进，该DOA估计系统还包括：

定义矩阵模块，用于定义矩阵 $K_1=\mathrm{Re}\left\{Q_{M-1}^H{J}_{s2}{J}_M\right\},$ $K_2=\mathrm{Im}\left\{Q_{M-1}^H{J}_{s2}{J}_M\right\};$ 定义选择矩阵模块，用于其中J_s1和J_s1分别定义为：J_s1=[I_(M-1)×M,0_(M-1)×1]，J_s2=[0_(M-1)×1,I_(M-1)×M]，表示选择矩阵；

最小二乘解模块，用于利用得到的信号子空间结合定义矩阵模块中的K₁与K₂，令其中表示矩阵伪逆运算；

对Ψ做特征值分解及DOA估计值模块，用于对Ψ做特征值分解，得其中ψ_i表示特征根，e_i表示相对于的特征矢量，则第i个信号的到达角表示为

${\widehat{\mathrm{\θ}}}_i={\sin }^{-1}\left(\frac{\mathrm{\λ}\·{\tan }^{-1}\left({\widehat{\mathrm{\ψ}}}_i\right)}{\mathrm{\πd}}\right),i=\mathrm{1,2},\·\·\·,P$

其中，和分别表示θ_i与ψ_i的估计值。

本发明的有益效果是：本发明只需要计算两个实数域的子样本协方差矩阵和通过简单的矩阵变换就可以直接获得信号子空间，避免了构造整个样本协方差矩阵以及其特征值分解，从而进一步降低了复杂度。

附图说明

图1是DOA估计均方误差与信噪比的关系图（M=15）。

图2是DOA估计均方误差与信噪比的关系图（M=25）。

图3是DOA估计均方误差与快拍数的关系图（M=15）。

图4是DOA估计均方误差与快拍数的关系图（M=25）。

图5是算法时间复杂度随阵元数增加的关系图。

图6是本发明DOA估计方法流程图。

具体实施方式

本发明公开了一种低复杂度的DOA估计方法，包括步骤W1至W3：

在步骤W1中，对样本数据做酉变换并将变换得到的数据分解为两部分：和其中，K为用户自定义参数且P≤K≤min{M,N}，P为信源个数；

在步骤W2中，计算的自相关矩阵以及与的互相关矩阵，并分别提取出它们的实部，分别用和表示；

通过步骤W1和步骤W2，避免了构造一个完整的样本协方差矩阵。同时，在后续的信号子空间构造中，该算法只涉及对和的矩阵运算。这一步所需的复杂度为O(MNK)，而传统的方法则需要O(MN²)，通常在天线大阵列的环境下，K的取值将远小于min{M,N}，因此复杂度O(MNK)《O(MN²)。

在步骤W3中，定义一个新的矩阵 $Z=\left[\begin{array}{c}{\tilde{R}}_{11}\\ {\tilde{R}}_{21}\end{array}\right]{\tilde{R}}_{11}^{-1/2},$ 对Z^HZ做特征值分解同时构造矩阵Π=ZU_Z，信号子空间由Π的前P列组成，即U_S=Π(:,1:P)。

通过步骤W3，避免对整个样本协方差矩阵的特征值分解，直接获得信号子空间步骤W3所需的复杂度为O(MK²)，而传统的方法则需要O(M³)，当K的取值远小于min{M,N}时，可以大大降低算法的复杂度。

通过仿真分析，改算法的DOA估计性能高于传统的ESPRIT算法，并且与传统的酉ESPRIT算法相当，但是，其所需的复杂度远低于其他两种算法，在天线尺寸很大的情况下，优势尤为明显。

如图6所示，作为本发明DOA估计方法的一具体实施例，该方法是基于方法的，具体包括以下步骤：

1)考虑一M阵元的均匀线阵，令X表示其接收的样本数据矩阵，令同时将其分解为 $\tilde{X}=\left[\begin{array}{c}{\tilde{X}}_1\\ {\tilde{X}}_2\end{array}\right],$ 其中 这里N为样本数，K为用户定义参数满足P≤K≤min{M,N}，P为信源个数；

2)分别计算的自相关矩阵R₁₁与和的互相关矩阵R₂₁；

3)令和分别表示R₁₁与R₂₁的实部；

4)利用公式 $Z=\left[\begin{array}{c}{\tilde{R}}_{11}\\ {\tilde{R}}_{21}\end{array}\right]{\tilde{R}}_{11}^{-1/2}$ 定义矩阵Z，令表示Z^HZ的特征值分解，其中Λ_Z=diag{λ₁,…,λ_K}表示由特征值组成的对角矩阵，并且所有特征值按降序排列，即λ₁≥…≥λ_K，U_Z=[u₁,…,u_K]表示相应的特征矢量矩阵，这里u_i表示第i个特征值λ_i所对应的特征矢量；

5)通过公式Π=ZU_Z构造矩阵Π，则信号子空间可以由Π的前P列构成，即U_S=Π(:,1:P)；

6)按如下公式定义稀疏酉矩阵

$Q_M=\left\{\begin{array}{cc}\frac{1}{\sqrt{2}}\left[\begin{array}{cc}{I}_l& j{I}_l\\ J_l& -j{J}_l\end{array}\right]& \mathrm{for\; M}=2l\\ \frac{1}{\sqrt{2}}\left[\begin{array}{ccc}{I}_l& 0_{l\×1}& j{I}_l\\ 0_{l\×1}^T& \sqrt{2}& 0_{l\×1}^T\\ J_l& 0_{l\×1}& -j{J}_l\end{array}\right]& \mathrm{for\; M}=2l+1\end{array}\right.$

式中，J_l为l×l的交换矩阵，其反对角线上的元素为1其余均为0；

7)定义矩阵 $K_1=\mathrm{Re}\left\{Q_{M-1}^H{J}_{s2}{J}_M\right\},$ $K_2=\mathrm{Im}\left\{Q_{M-1}^H{J}_{s2}{J}_M\right\},$ 其中J_s1和J_s1分别定义为：J_s1=[I_(M-1)×M,0_(M-1)×1]，J_s2=[0_(M-1)×1,I_(M-1)×M]，表示选择矩阵；

8)利用步骤5得到的信号子空间结合步骤7中的K₁与K₂，令其中表示矩阵伪逆运算；

9)对Ψ做特征值分解，得其中ψ_i表示特征根，e_i表示相对于的特征矢量，则第i个信号的到达角可以表示为

${\widehat{\mathrm{\θ}}}_i={\sin }^{-1}\left(\frac{\mathrm{\λ}\·{\tan }^{-1}\left({\widehat{\mathrm{\ψ}}}_i\right)}{\mathrm{\πd}}\right),i=\mathrm{1,2},\·\·\·,P$

这里，和分别表示θ_i与ψ_i的估计值。

本方法通过酉变换Q_M将整个DOA估计所涉及的矩阵运算从复数域变换到实数域，通过方法构造Π，直接得到信号子空间U_S，避免了对整个样本协方差矩阵的计算以及对其做特征值分解来获得U_S，复杂度很低。

本发明还公开了一种低复杂度的DOA估计系统，包括：

复数据实数化模块，用于令X表示其接收的样本数据矩阵，令 $\tilde{X}=Q_M^{H}X;$

分解数据样本模块，用于分解数据样本为 $\tilde{X}=\left[\begin{array}{c}{\tilde{X}}_1\\ {\tilde{X}}_2\end{array}\right],$ 其中 这里N为样本数，K为用户定义参数满足P≤K≤min{M,N}，P为信源个数；

估计协方差矩阵模块，用于分别计算的自相关矩阵R₁₁与和的互相关矩阵R₂₁，令和分别表示R₁₁与R₂₁的实部。

该DOA估计系统还包括：

定义矩阵Z模块，用于定义一个新的矩阵 $Z=\left[\begin{array}{c}{\tilde{R}}_{11}\\ {\tilde{R}}_{21}\end{array}\right]{\tilde{R}}_{11}^{-1/2};$

对Z^HZ做特征值分解模块，用于令表示Z^HZ的特征值分解，

其中Λ_Z=diag{λ₁,…,λ_K}表示由特征值组成的对角矩阵，并且所有特征值按降序排列，即λ₁≥…≥λ_K，U_Z=[u₁,…,u_K]表示相应的特征矢量矩阵，这里u_i表示第i个特征值λ_i所对应的特征矢量；

构造矩阵模块，用于通过公式Π=ZU_Z构造矩阵Π；

得到信号子空间模块，用于信号子空间由Ξ的前P列构成，即U_S=Π(:,1:P)。

该DOA估计系统还包括：

定义矩阵模块，用于定义矩阵 $K_1=\mathrm{Re}\left\{Q_{M-1}^H{J}_{s2}{J}_M\right\},$ $K_2=\mathrm{Im}\left\{Q_{M-1}^H{J}_{s2}{J}_M\right\};$

定义选择矩阵模块，用于其中J_s1和J_s1分别定义为：J_s1=[I_(M-1)×M,0_(M-1)×1]，J_s2=[0_(M-1)×1,I_(M-1)×M]，表示选择矩阵；

最小二乘解模块，用于利用得到的信号子空间结合定义矩阵模块中的K₁与K₂，令其中表示矩阵违逆运算；

对Ψ做特征值分解及DOA估计值模块，用于对Ψ做特征值分解，得其中ψ_i表示特征根，e_i表示相对于的特征矢量，则第i个信号的到达角表示为

${\widehat{\mathrm{\θ}}}_i={\sin }^{-1}\left(\frac{\mathrm{\λ}\·{\tan }^{-1}\left({\widehat{\mathrm{\ψ}}}_i\right)}{\mathrm{\πd}}\right),i=\mathrm{1,2},\·\·\·,P$

其中，和分别表示θ_i与ψ_i的估计值。

本发明的目的是解决传统的旋转不变子空间类DOA估计算法高复杂度问题，提出了一种新的基于方法的快速酉ESPRIT算法，该方法只需要计算两个实数域的子样本协方差矩阵和通过简单的矩阵变换就可以直接获得信号子空间，避免了构造整个样本协方差矩阵以及其特征值分解，从而进一步降低了复杂度。

考虑一M阵元的均匀线阵，阵元间距d=λ/2，λ表示信号波长。假设P个独立的远场窄带信号入射到该阵列，入射角分别为θ₁,…,θ_P；噪声是空时均白的高斯随机过程，方差为σ²，且噪声与信号互不相关，则阵列接收数据可以表示为

X(t)=AS(t)+N(t)fort=1,…,N

式中，A=[a₁(θ),…,a_P(θ)]为阵列流型，S(t)=[s₁(t),…,s_P(t)]^T为信号矩阵，N(t)为噪声， $a_i\left({\mathrm{\θ}}_i\right)=[1,e^{-j2\mathrm{\π}d\sin \left({\mathrm{\θ}}_i\right)/\mathrm{\λ}},\·\·\·,e^{-j2\mathrm{\π}(M-1)d\sin \left({\mathrm{\θ}}_i\right)/\mathrm{\λ}}]$ 表示第i个信号的导向矢量，则阵列接收协方差矩阵可以表示为

R=E{X(t)X^H(t)}=AR_SA^H+σ²I

3.1信号子空间估计

将阵列接收数据X分解为

$X=\left[\begin{array}{c}{X}_1\\ X_2\end{array}\right]$

令 为了得到完整的信号子空间，K必须大于信号个数P并且小于min(M,N)，即{K|P≤K≤min(M,N)}。

信号子空间由如下定理给出：

定理1：令 $Z=\left[\begin{array}{c}{R}_{11}\\ R_{21}\end{array}\right]R_{11}^{-1/2},$ 且为Z^HZ的EVD分解，其中Λ_Z=diag[λ₁,…,λ_K]为特征值矩阵，U_z=[u₁,…,u_K]为相应的特征矢量矩阵，这里特征根以降序排列，即λ₁≥…≥λ_K。则信号子空间可以表示为

U_s=Π(:,1:P)

式中，Π=ZU_Z。

证明：证明见本章第3.3节。

利用阵列的旋转不变性，可知

J_s1AΦ=J_s2A(3-1)

我们的目的是利用酉变换Q_M将复阵列流行A转换为实值阵列流行即

$\tilde{A}=Q_M^{H}A---(3-2)$

其中

$Q_M=\left\{\begin{array}{cc}\frac{1}{\sqrt{2}}\left[\begin{array}{cc}{I}_l& j{I}_l\\ J_l& -j{J}_l\end{array}\right]& \mathrm{for\; M}=2l\\ \frac{1}{\sqrt{2}}\left[\begin{array}{ccc}{I}_l& 0_{l\×1}& j{I}_l\\ 0_{l\×1}^T& \sqrt{2}& 0_{l\×1}^T\\ J_l& 0_{l\×1}& -j{J}_l\end{array}\right]& \mathrm{for\; M}=2l+1\end{array}\right.---(3-3)$

是稀疏酉矩阵，J_l为交换矩阵，其反对角线上的元素为1其余均为0。

为了找到一个与实值阵列流型对应的实值信号子空间我们令同时将其分解为

$\tilde{X}=\left[\begin{array}{c}{\tilde{X}}_1\\ {\tilde{X}}_2\end{array}\right]---(3-4)$

令 ${\tilde{R}}_{11}=\mathrm{Re}\left\{E\right[{\tilde{X}}_1{\tilde{X}}_1^H\left]\right\},$ ${\tilde{R}}_{21}=\mathrm{Re}\left\{E\right[{\tilde{X}}_2{\tilde{X}}_1^H\left]\right\},$ 这里，Re{·}表示取实部。将与代入定理1中，就可以得到

3.2波达方向估计

将式（3-1）改写为

$J_{s1}{Q}_M\tilde{A}\mathrm{\Φ}=J_{s2}{Q}_M\tilde{A}---(3-5)$

由于将式（3-5）左乘经过简单的数学变换，我们可以得到

$K_1\tilde{A}\mathrm{\Ω}=K_2\tilde{A}---(3-6)$

其中，

$\mathrm{\Ω}=\mathrm{diag}\{\tan \left(\frac{{\mathrm{\φ}}_1}{2}\right),\·\·\·,\tan \left(\frac{{\mathrm{\φ}}_P}{2}\right)\}---(3-7)$

$K_1=\mathrm{Re}\left\{Q_{M-1}^H{J}_{s2}{J}_M\right\}---(3-8)$

$K_2=\mathrm{Im}\left\{Q_{M-1}^H{J}_{s2}{J}_M\right\}---(3-9)$

由定理1可知，和张成同一个空间，因此，必定存在一个非奇异矩阵T使得

$\tilde{A}={\tilde{U}}_S{T}^{-1}---(3-10)$

将式（3-10）代入式（3-6），可得

$K_1{\tilde{U}}_S\mathrm{\Ψ}=K_2{\tilde{U}}_S---(3-11)$

式中，Ψ=T^-1ΦT。这里需要强调Φ和Ψ是相似矩阵，因此，他们有相同的特征根。通过最小二乘，我们很容易得到

对Ψ做特征值分解，即

$\mathrm{\Ψ}=\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{\mathrm{\ψ}}_i{e}_i{e}_i^H---(3-13)$

因为Ψ和Φ拥有相同的特征根，则信号到达角可以表示为

${\widehat{\mathrm{\θ}}}_i={\sin }^{-1}\left(\frac{\mathrm{\λ}\·{\tan }^{-1}\left({\widehat{\mathrm{\ψ}}}_i\right)}{\mathrm{\πd}}\right),i=\mathrm{1,2},\·\·\·,P---(3-14)$

3.3证明定理1

假设U_S是理想的信号子空间，则必定存在一个非奇异矩阵T使得U_S=AT。事实上，Π可以改写为

$\mathrm{\Π}={\mathrm{ZU}}_Z=\left[\begin{array}{c}{R}_{11}\\ R_{21}\end{array}\right]Q=\left[\begin{array}{c}{R}_{11}\\ R_{21}\end{array}\right]R_{11}^{-1/2}{U}_Z---(3-15)$

式中，R₁₁与U_Z均为满秩矩阵，也为满秩矩阵。因此，

$\left[\begin{array}{c}{R}_{11}\\ R_{21}\end{array}\right]=\left[\begin{array}{c}{A}_1{R}_s{A}_1^H\\ A_2{R}_s{A}_1^H\end{array}\right]+\left[\begin{array}{c}{\mathrm{\σ}}_n^2{I}_K\\ 0_{(M-K)\×K}\end{array}\right]$

$=A(R_s{A}_1^H+{\left(A^{H}A\right)}^{-1}{A}^H\left[\begin{array}{c}{\mathrm{\σ}}_n^2{I}_K\\ 0_{(M-K)\×K}\end{array}\right])---(3-16)$

$=A(R_s+{\mathrm{\σ}}_n^2{\left(A^{H}A\right)}^{-1})A_1^H$

将式（3-16）代入（3-15），易得

$\mathrm{\Π}=A(R_s+{\mathrm{\σ}}_n^2{\left(A^{H}A\right)}^{-1})A_1^{H}Q---(3-17)$

$=\mathrm{AB}{A}_1^{H}Q$

式中，

又因为

$R_x=U_S{\mathrm{\Λ}}_S{U}_S^H+{\mathrm{\σ}}_n^2{U}_n{U}_n^H---(3-18)$

$={\mathrm{AR}}_S{A}^H+{\mathrm{\σ}}_n^2{I}_M$

上式右乘U_S得

${\mathrm{AR}}_S{A}^H{U}_S=U_S({\mathrm{\Λ}}_S-{\mathrm{\σ}}_n^2{I}_M)---(3-19)$

因为 $U_S{U}_S^H=A{\left(A^{H}A\right)}^{-1}{A}^H,$ 由式(3-19)可得

又因为 $B=R_S+{\mathrm{\σ}}_n^2{\left(A^{H}A\right)}^{-1},$ 得

将U_S=AT代入上式，即可得

B=TΛ_ST^H(3-22)

由于T和Λ_S都是满秩矩阵，且他们的秩均为P，所以B也是秩为P的满秩矩阵。

令则Π=AH。又因为具有范德蒙结构，矩阵H的前P列是线性独立的。因此，必定存在一个非奇异矩阵使得

$U_S=\mathrm{\Π}(:,1:P)=A\tilde{T}---(3-23)$

式中,定理1得证。

考虑一均匀线阵，阵列间距为d=λ/2。两个互不相干的窄带信号入射到该阵列，信号入射方向为θ₁=1°和θ₂=3°。此外，噪声是均值为零的白高斯噪声，且与信号互不相关。实验中，我们比较传统的ESPIRT算法、酉ESPRIT算法和本发明所提出的快速酉ESPRIT算法。所有的仿真结果均由500次蒙特卡洛实验获得。

实验1算法统计性能与信噪比的关系。

在这个仿真中，我们固定快拍数N=40，信噪比从-10dB增加到10dB。图1为阵元数M=15的情况下的DOA估计均方根误差图；图2为M=25时的DOA估计均方根误差图。对比图1和图2可以看出本发明所提的算法总是可以获得和酉ESPRIT算法相同的性能，无论K是大还是小。这也说明该算法对K的选择并不敏感，倘若K≥P。在低信噪比时，本发明的性能始终优于传统的ESPRIT算法。当信噪比大于5dB时，所有的三种算法都获得了相同的性能。所以，当M和N固定且信噪比较大时，我们可以选择一个相对小的K来进一步降低复杂度。从图中还可以发现，当K取值变大时，本发明的性能也逐渐提升。

实验2算法统计性能与快拍数的关系。

在此仿真中，我们固定信噪比为8dB。图3中，K=M=15时，本发明所提的算法退化为酉ESPRIT算法。从图4中可以看出，本发明所提出的算法性能始终优于ESPRIT算法。在低信噪比的情况下，所有的酉ESPRIT算法的性能都远好于ESPRIT算法。与酉ESPRIT算法相比，当K=12时，本发明所提的算法性能稍差；但是，随着K逐渐增大，例如K=15，该算法的性能与酉ESPRIT算法相当。

实验3时间复杂度

假设DOA个数为2，快拍数为100，仿真结果如图5。随着M增大，本发明所提的算法比其他两种算法所需的时间更少，尤其是当M很大时。

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

Claims (4)

1.一种低复杂度的DOA估计方法，其特征在于，包括如下步骤：

A.对样本数据做酉变换并将变换得到的数据分解为两部分：和其中，K为用户自定义参数且P≤K≤min{M,N}，P为信源个数；

B.计算的自相关矩阵以及与的互相关矩阵，并分别提取出它们的实部，分别用和表示；

C.定义一个新的矩阵 $Z=\left[\begin{array}{c}{\tilde{R}}_{11}\\ {\tilde{R}}_{21}\end{array}\right]{\tilde{R}}_{11}^{-1/2},$ 对Z^HZ做特征值分解同时构造矩阵Π＝ZU_Z，信号子空间由Π的前P列组成，即U_S＝Π(:,1:P)；

D.按如下公式定义稀疏酉矩阵

$Q_M=\left\{\begin{array}{cc}\frac{1}{\sqrt{2}}\left[\begin{array}{cc}{I}_l& {\mathrm{jI}}_l\\ J_l& -{\mathrm{jJ}}_l\end{array}\right]& forM=2l\\ \frac{1}{\sqrt{2}}\left[\begin{array}{ccc}{I}_l& 0_{l\×1}& {\mathrm{jI}}_l\\ 0_{l\×1}^T& \sqrt{2}& 0_{l\×1}^T\\ J_l& 0_{l\×1}& -{\mathrm{jJ}}_l\end{array}\right]& forM=2l+1\end{array}\right.$

式中，J_l为l×l的交换矩阵，其反对角线上的元素为1其余均为0；

E.定义矩阵 $K_1=\mathrm{Re}\left\{Q_{M-1}^H{J}_{s2}{J}_M\right\},K_2=\mathrm{Im}\left\{Q_{M-1}^H{J}_{s2}{J}_M\right\},$ 其中J_s1和J_s2分别定义为：

J_S1＝[I_(M-1)×(M-1),O_(M-1)×1]

J_S2＝[O_(M-1)×1,I_(M-1)×(M-1)],表示选择矩阵；

F.利用步骤C中得到的信号子空间结合步骤E中的K₁与K₂，令其中表示矩阵违逆运算；

G.对Ψ做特征值分解，得其中ψ_i表示特征根，e_i表示相对于的特征矢量，则第i个信号的到达角表示为

${\widehat{\mathrm{\θ}}}_i={\sin }^{-1}\left(\frac{\mathrm{\λ}\·{\tan }^{-1}\left({\widehat{\mathrm{\ψ}}}_i\right)}{\mathrm{\π}d}\right),i=1,2,...,P$

其中，和分别表示θ_i与ψ_i的估计值，C为复数域，N为样本数，M为阵元个数，U_Z为Nystrom分解得到的特征矢量矩阵，Λ_Z为与U_Z对应的特征值矩阵；

其中，l是为了说明Q_M如何产生而引入的一个临时变量；J_M是交换矩阵即反对角线元素为1，其余元素均为0；0_(M-1)×1阵为(M-1)×1的零向量,I_(M-1)×(M-1)是一个(M-1)×(M-1)的单位阵；U_s为信号子空间；d＝λ/2为阵列相邻两个阵元之间的间距,λ为发射信号波长。

2.根据权利要求1所述的DOA估计方法，其特征在于，在所述步骤A中包括如下步骤：

A1.考虑一M阵元的均匀线阵，令X表示其接收的样本数据矩阵，令 $\tilde{X}=Q_M^{H}X;$

A2.分解数据样本为 $\tilde{X}=\left[\begin{array}{c}{\tilde{X}}_1\\ {\tilde{X}}_2\end{array}\right],$ 其中这里N为样本数，K为用户定义参数满足P≤K≤min{M,N}，P为信源个数；

Q_M ^H定义为： $Q_M=\left\{\begin{array}{cc}\frac{1}{\sqrt{2}}\left[\begin{array}{cc}{I}_l& {\mathrm{jI}}_l\\ J_l& -{\mathrm{jJ}}_l\end{array}\right]& forM=2l\\ \frac{1}{\sqrt{2}}\left[\begin{array}{ccc}{I}_l& 0_{l\×1}& {\mathrm{jI}}_l\\ 0_{l\×1}^T& \sqrt{2}& 0_{l\×1}^T\\ J_l& 0_{l\×1}& -{\mathrm{jJ}}_l\end{array}\right]& forM=2l+1\end{array}\right.$ 其用途是获得实值样本协方差矩阵，即将复数域的样本协方差矩阵转换到实数域。

3.根据权利要求2所述的DOA估计方法，其特征在于，在所述步骤B中包括如下步骤：

B1.分别计算的自相关矩阵R₁₁与和的互相关矩阵R₂₁；

B2.令和分别表示R₁₁与R₂₁的实部。

4.根据权利要求3所述的DOA估计方法，其特征在于，在所述步骤C中包括如下步骤：

C1.利用公式 $Z=\left[\begin{array}{c}{\tilde{R}}_{11}\\ {\tilde{R}}_{21}\end{array}\right]{\tilde{R}}_{11}^{-1/2}$ 定义矩阵Z，令表示Z^HZ的特征值分解，其中Λ_Z＝diag{λ₁,…,λ_K}表示由特征值组成的对角矩阵，并且所有特征值按降序排列，即λ₁≥…≥λ_K，U_Z＝[u₁，…,u_K]表示相应的特征矢量矩阵，这里u_i表示第i个特征值λ_i所对应的特征矢量；

C2.通过公式Π＝ZU_Z构造矩阵Π，则信号子空间由Π的前P列构成，即U_S＝Π(:,1:P)。