CN104298850A - 信源数未知的相干信号测向方法及系统 - Google Patents

Abstract

本发明提供了一种信源数未知的相干信号测向方法及系统，该相干信号测向方法包括初始化步骤、四阶累积量矩阵计算步骤、构造Toeplitz矩阵步骤、构造中间变量步骤、功率谱计算步骤、搜索步骤。本发明的有益效果是：本发明通过构造一组基于四阶累积量的Toeplitz矩阵对相干信号进行解相关，由于Toeplitz矩阵具有联合对角化结构，所以不需要任何有关信源数目的先验知识，本发明设计出了一种新的代价函数，能有效地估计出波达方向。

Description

信源数未知的相干信号测向方法及系统

技术领域

本发明涉及信号处理领域，尤其涉及信源数未知的相干信号测向方法及系统。 

背景技术

利用空间分布的传感器阵列测向是雷达、声呐、无线通信等应用中的重要任务。目前，提出了许多波达方向估计算法，其中，基于子空间的算法，例如ESPRIT和MUSIC算法在假设条件满足的情况下为估计波达方向提供了一个很好的解决方案。其假设条件为，信源数目已知，信源互不相关或部分相关且噪声是空间不相关的白噪声，即协方差矩阵与单位矩阵成正比。当任何一个假设条件不成立时，基于子空间算法的性能将严重下降。 

实际情况中，信源数目往往是未知的。最常用的信源数估计方法是赤池信息量准则(AIC)和最小描述准则(MDL)，然而，当采样数太小或信噪比太低时，这两种估计方法并不适用。尽管后来提出了各种改进的算法，但在一些极端的条件下，尤其是未知噪声性质时，其正确检测的概率仍然很低。 

实际传输过程中由于多路径传播的影响，接收的数据中包含了许多相干信号而导致信源协方差矩阵是秩亏的，这意味着基于子空间的方法中信源互不相关的假设条件无法成立。在这种情况下，空间平滑技术(SS)及其改进算法被提出来解决了信号相干的问题，这些算法对信源协方差矩阵做了预处理，将阵列分成多个子阵并求平均使得输出的协方差矩阵满秩，满足了子空间算法的假设条件。 

众所周知，大多数假设噪声为空间不相关的DOA估计算法对噪声模型很敏感。尽管如此，当已知协方差矩阵时，空间相关的噪声很容易通过预白化处理来满足不相关条件。但实际情况下，有限的观测值计算得到的噪声协方差矩阵并不精确。 

发明内容

为了解决现有技术中的问题，本发明提供了一种信源数未知的相干信号测向方法。 

本发明提供了一种无需信源数估计的相干信号波达方向估计方法，包括如下步骤： 

初始化步骤：考虑一具有(2M+1)个阵元的均匀线阵，假设有P(P≤M+1)个窄带远场信号从不同方向{θ₁,…,θ_P}入射到该阵列，前K个信号是互相关的，其余(P-K)个信号是互不相关的并且独立于前K个信号，x(t)表示接收样本数据； 

四阶累积量矩阵计算步骤：计算样本x(t)的(2M+1)(M+1)个四阶累积量矩阵，C_i,j,-M≤i≤j≤M； 

构造Toeplitz矩阵步骤：对每个累积量矩阵的前(M+1)行构造(M+1)个Toeplitz矩阵

构造中间变量步骤：计算矩阵和 

其中 

$a\left({\mathrm{\θ}}_p\right)={[e^{j2\mathrm{\π}M\sin {\text{\θ}}_p\mathrm{\Δ}/\mathrm{\λ}},...,1,...,e^{-j2\mathrm{\π}M\sin {\mathrm{\θ}}_p\mathrm{\Δ}/\mathrm{\λ}}]}^T$ 是导向矢量； 

功率谱计算步骤：通过公式 

计算功率谱； 

搜索步骤：搜索θ，找到P(θ)的峰值对应的角度即为DOA的估计值。 

作为本发明的进一步改进，在所述初始化步骤中：考虑一具有(2M+1)阵元的均匀线阵(ULA)，假设有P(P≤M+1)个远场窄带信号从不同方向{θ₁,…,θ_P}入射到该阵列，前K个信号是互相干的，其余(P-K)信号均是互不相关的且独立于前K个信号，令第一个信号d₁(t)作为参考信号，则第k个相干信号可以表示为： 

$s_k\left(t\right)={\mathrm{\β}}_k{e}^{\mathrm{j\δ}{\mathrm{\φ}}_k}{s}_1\left(t\right),k=2,...,K---\left(1\right)$

其中，β_k幅度衰减因子，δφ_k是相位变化。实际上，δφ_k不会影响信号之间的相干性，不失一般性，令δφ_k＝0,k＝2,…,K，第m个阵元所接收的信号可以表示为 

$\begin{array}{c}{x}_m\left(t\right)=\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{s}_i\left(t\right)e^{-j2\mathrm{\π}m\sin {\mathrm{\θ}}_i\mathrm{\Δ}/\mathrm{\λ}}+n_m\left(t\right)\\ =s_1\left(t\right)\underset{i=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\β}}_i{e}^{-j2\mathrm{\π}m\sin {\mathrm{\θ}}_i\mathrm{\Δ}/\mathrm{\λ}}\\ +\underset{i=K+1}{\overset{P}{\mathrm{\Σ}}}{s}_i\left(t\right)e^{-j2\mathrm{\π}m\sin {\mathrm{\θ}}_i\mathrm{\Δ}/\mathrm{\λ}}+n_m\left(t\right)\end{array}---\left(2\right)$

其中，s_i(t)是第i个信号的复包络，β₁＝1,λ是载波波长，Δ＝λ/2是阵元间间距，(2)式可改写为向量形式： 

x(t)＝[x_-M(t)，…，x₀(t)，…，x_M(t)]^T

                                   (3) 

＝Ad(t)+n(t) 

其中s(t)＝[s₁(t),…,s_P(t)]^T是源信号向量，A＝[a(θ₁),…,a(θ_P)]是阵列流行，其中，第p个导向矢量可以表示为： 

$a\left({\mathrm{\θ}}_p\right)={[e^{j2\mathrm{\π}M\sin {\text{\θ}}_p\mathrm{\Δ}/\mathrm{\λ}},...,1,...,e^{-j2\mathrm{\π}M\sin {\mathrm{\θ}}_p\mathrm{\Δ}/\mathrm{\λ}}]}^T.---\left(4\right)$

作为本发明的进一步改进，在所述四阶累积量矩阵计算步骤中,零均值平稳信号的四阶累积量定义如下： 

$\begin{array}{c}\mathrm{cum}(x_{k_1}\left(t\right),x_{k_2}^{*}\left(t\right),x_{l_1}\left(t\right),x_{l_2}^{*}\left(t\right))\\ =E\left\{x_{k_1}\left(t\right)x_{k_2}^{*}\left(t\right)x_{l_1}\left(t\right)x_{l_2}^{*}\left(t\right)\right\}\\ -E\left\{x_{k_1}\left(t\right)x_{k_2}^{*}\left(t\right)\right\}E\left\{x_{l_1}\left(t\right)x_{l_2}^{*}\left(t\right)\right\}\\ -E\left\{x_{k_1}\left(t\right)x_{l_1}\left(t\right)\right\}E\left\{x_{k_2}^{*}\left(t\right)x_{l_2}^{*}\left(t\right)\right\}\\ -E\left\{x_{k_1}\left(t\right)x_{l_2}^{*}\left(t\right)\right\}E\left\{x_{k_2}^{*}\left(t\right)x_{l_1}\left(t\right)\right\}\end{array}\text{---}\left(8\right)$

其中，-M≤k₁,k₂,l₁,l₂≤M； 

定义一个M×M的累积量矩阵它的第(l₁,l₂)个元素为 

$C_{k_1,k_2}(l_1,l_2)=\mathrm{cum}(x_{k_1}\left(t\right),x_{k_2}^{*}\left(t\right),x_{l_1}\left(t\right),x_{l_2}^{*}\left(t\right)).---\left(9\right)$

累积-M≤l₁,l₂≤M范围内的元素得到 

$\begin{array}{c}{C}_{k_1,k_2}=\mathrm{cum}(x_{k_1}\left(t\right),x_{k_2}^{*}\left(t\right),x\left(t\right),x^{*}\left(t\right))\\ =A\left(\mathrm{\θ}\right){\mathrm{\Γ}}_{k_1,k_2}{A}^H\left(\mathrm{\θ}\right)\end{array}\text{---}\left(10\right)$

其中 

是第p个信源的四阶累积量。在-M≤l₁,l₂≤M范 围内，总共可以得到(2M+1)²个累积量矩阵，然而，由是中心厄米特矩阵可知 

$C_{k_1,k_2}=C_{k_2,k_1}^H---\left(12\right)$

即和包含相同的数据信息，对应的特征向量和特征值也相同。 

作为本发明的进一步改进，在所述构造Toeplitz矩阵步骤中，考虑(2M+1)均匀线阵，利用四阶累积量技术构造(2M+1)(M+1)²个Toeplitz矩阵；在所述构造中间变量步骤中，利用(2M+1)(M+1)²个Toeplitz矩阵构造出中间变量F和G(θ)。 

本发明还提供了一种信源数未知的相干信号测向系统，包括： 

初始化单元：用于考虑一具有(2M+1)个阵元的均匀线阵，假设有P(P≤M+1)个窄带远场信号从不同方向{θ₁,…,θ_P}入射到该阵列，前K个信号是互相关的，其余(P-K)个信号是互不相关的并且独立于前K个信号，x(t)表示接收样本数据； 

四阶累积量矩阵计算单元：用于计算样本x(t)的(2M+1)(M+1)个四阶累积量矩阵，C_i,j,-M≤i≤j≤M； 

构造Toeplitz矩阵单元：用于对每个累积量矩阵的前(M+1)行构造(M+1)个Toeplitz矩阵

构造中间变量单元：用于计算矩阵和 

其中 

$a\left({\mathrm{\θ}}_p\right)={[e^{j2\mathrm{\π}M\sin {\text{\θ}}_p\mathrm{\Δ}/\mathrm{\λ}},...,1,...,e^{-j2\mathrm{\π}M\sin {\mathrm{\θ}}_p\mathrm{\Δ}/\mathrm{\λ}}]}^T$ 是导向矢量； 

功率谱计算单元：用于通过公式 

计算功率谱； 

搜索单元：用于搜索θ，找到P(θ)的峰值对应的角度即为DOA的估计值。 

作为本发明的进一步改进，在所述初始化单元中：考虑一具有(2M+1) 阵元的均匀线阵(ULA)，假设有P(P≤M+1)个远场窄带信号从不同方向{θ₁,…,θ_P}入射到该阵列，前K个信号是互相干的，其余(P-K)信号均是互不相关的且独立于前K个信号，令第一个信号d₁(t)作为参考信号，则第k个相干信号可以表示为： 

$s_k\left(t\right)={\mathrm{\β}}_k{e}^{\mathrm{j\δ}{\mathrm{\φ}}_k}{s}_1\left(t\right),k=2,...,K---\left(1\right)$

其中，β_k幅度衰减因子，δφ_k是相位变化。实际上，δφ_k不会影响信号之间的相干性，不失一般性，令δφ_k＝0,k＝2,…,K，第m个阵元所接收的信号可以表示为 

$\begin{array}{c}{x}_m\left(t\right)=\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{s}_i\left(t\right)e^{-j2\mathrm{\π}m\sin {\mathrm{\θ}}_i\mathrm{\Δ}/\mathrm{\λ}}+n_m\left(t\right)\\ =s_1\left(t\right)\underset{i=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\β}}_i{e}^{-j2\mathrm{\π}m\sin {\mathrm{\θ}}_i\mathrm{\Δ}/\mathrm{\λ}}\\ +\underset{i=K+1}{\overset{P}{\mathrm{\Σ}}}{s}_i\left(t\right)e^{-j2\mathrm{\π}m\sin {\mathrm{\θ}}_i\mathrm{\Δ}/\mathrm{\λ}}+n_m\left(t\right)\end{array}---\left(2\right)$

其中，s_i(t)是第i个信号的复包络，β₁＝1,λ是载波波长，Δ＝λ/2是阵元间间距，(2)式可改写为向量形式： 

x(t)＝[x_-M(t),…,x₀(t),…,x_M(t)]^T

                                  (3) 

＝Ad(t)+n(t) 

其中s(t)＝[s₁(t),…,s_P(t)]^T是源信号向量，A＝[a(θ₁),…,a(θ_P)]是阵列流行，其中，第p个导向矢量可以表示为： 

$a\left({\mathrm{\θ}}_p\right)={[e^{j2\mathrm{\π}M\sin {\text{\θ}}_p\mathrm{\Δ}/\mathrm{\λ}},...,1,...,e^{-j2\mathrm{\π}M\sin {\mathrm{\θ}}_p\mathrm{\Δ}/\mathrm{\λ}}]}^T.---\left(4\right)$

作为本发明的进一步改进，在所述四阶累积量矩阵计算步骤中,零均值平稳信号的四阶累积量定义如下： 

$\begin{array}{c}\mathrm{cum}(x_{k_1}\left(t\right),x_{k_2}^{*}\left(t\right),x_{l_1}\left(t\right),x_{l_2}^{*}\left(t\right))\\ =E\left\{x_{k_1}\left(t\right)x_{k_2}^{*}\left(t\right)x_{l_1}\left(t\right)x_{l_2}^{*}\left(t\right)\right\}\\ -E\left\{x_{k_1}\left(t\right)x_{k_2}^{*}\left(t\right)\right\}E\left\{x_{l_1}\left(t\right)x_{l_2}^{*}\left(t\right)\right\}\\ -E\left\{x_{k_1}\left(t\right)x_{l_1}\left(t\right)\right\}E\left\{x_{k_2}^{*}\left(t\right)x_{l_2}^{*}\left(t\right)\right\}\\ -E\left\{x_{k_1}\left(t\right)x_{l_2}^{*}\left(t\right)\right\}E\left\{x_{k_2}^{*}\left(t\right)x_{l_1}\left(t\right)\right\}\end{array}\text{---}\left(8\right)$

其中，-M≤k₁,k₂,l₁,l₂≤M； 

定义一个M×M的累积量矩阵它的第(l₁,l₂)个元素为 

$C_{k_1,k_2}(l_1,l_2)=\mathrm{cum}(x_{k_1}\left(t\right),x_{k_2}^{*}\left(t\right),x_{l_1}\left(t\right),x_{l_2}^{*}\left(t\right)).---\left(9\right)$

累积-M≤l₁,l₂≤M范围内的元素得到 

$\begin{array}{c}{C}_{k_1,k_2}=\mathrm{cum}(x_{k_1}\left(t\right),x_{k_2}^{*}\left(t\right),x\left(t\right),x^{*}\left(t\right))\\ =A\left(\mathrm{\θ}\right){\mathrm{\Γ}}_{k_1,k_2}{A}^H\left(\mathrm{\θ}\right)\end{array}\text{---}\left(10\right)$

其中 

是第p个信源的四阶累积量。在-M≤l₁,l₂≤M范围内，总共可以得到(2M+1)²个累积量矩阵，然而，由是中心厄米特矩阵可知 

$C_{k_1,k_2}=C_{k_2,k_1}^H---\left(12\right)$

即和包含相同的数据信息，对应的特征向量和特征值也相同。 

作为本发明的进一步改进，在所述构造Toeplitz矩阵单元中，考虑(2M+1)均匀线阵，利用四阶累积量技术构造(2M+1)(M+1)²个Toeplitz矩阵；在所述构造中间变量单元中，利用(2M+1)(M+1)²个Toeplitz矩阵构造出中间变量F和G(θ)。 

本发明的有益效果是：本发明通过构造一组基于四阶累积量的Toeplitz矩阵对相干信号进行解相关，由于Toeplitz矩阵具有联合对角化结构，所以不需要任何有关信源数目的先验知识，本发明设计出了一种新的代价函数，能有效地估计出波达方向。 

附图说明

图1是对称均匀线阵图。 

图2是空间谱比较图(垂直线代表真实DOA)。 

图3是空间谱比较图(垂直线代表真实DOA)。 

图4是DOA估计均方根误差与信噪比的关系图。 

图5是DOA分辨概率与信噪比的关系图。 

图6是DOA估计均方根误差与快拍数的关系图。 

图7是DOA分辨概率与快拍数的关系图。 

图8是DOA估计均方根误差与相关系数的关系图。 

图9是DOA分辨概率与相关系数的关系图。 

具体实施方式

本发明公开了一种信源数未知的相干信号测向方法，包括如下步骤： 

考虑一具有(2M+1)阵元的均匀线阵(ULA)。假设有P(P≤M+1)个远场 窄带信号从不同方向{θ₁,…,θ_P}入射到该阵列，前K个信号是互相干的，其余(P-K)信号均是互不相关的且独立于前K个信号。令第一个信号d₁(t)作为参考信号，则第k个相干信号可以表示为 

$s_k\left(t\right)={\mathrm{\β}}_k{e}^{\mathrm{j\δ}{\mathrm{\φ}}_k}{s}_1\left(t\right),k=2,...,K---\left(1\right)$

其中，β_k幅度衰减因子，δφ_k是相位变化。实际上，δφ_k不会影响信号之间的相干性，不失一般性，令δφ_k＝0,k＝2,…,K，第m个阵元所接收的信号可以表示为 

$\begin{array}{c}{x}_m\left(t\right)=\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{s}_i\left(t\right)e^{-j2\mathrm{\π}m\sin {\mathrm{\θ}}_i\mathrm{\Δ}/\mathrm{\λ}}+n_m\left(t\right)\\ =s_1\left(t\right)\underset{i=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\β}}_i{e}^{-j2\mathrm{\π}m\sin {\mathrm{\θ}}_i\mathrm{\Δ}/\mathrm{\λ}}\\ +\underset{i=K+1}{\overset{P}{\mathrm{\Σ}}}{s}_i\left(t\right)e^{-j2\mathrm{\π}m\sin {\mathrm{\θ}}_i\mathrm{\Δ}/\mathrm{\λ}}+n_m\left(t\right)\end{array}---\left(2\right)$

其中，s_i(t)是第i个信号的复包络，β₁＝1,λ是载波波长，Δ＝λ/2是阵元间间距，(2)式可改写为向量形式： 

x(t)＝[x_-M(t),…,x₀(t),…,x_M(t)]^T

                                 (3) 

＝Ad(t)+n(t) 

其中s(t)＝[s₁(t),…,s_P(t)]^T是源信号向量，A＝[a(θ₁),…,a(θ_P)]是阵列流行，其中，第p个导向矢量可以表示为： 

$a\left({\mathrm{\θ}}_p\right)={[e^{j2\mathrm{\π}M\sin {\text{\θ}}_p\mathrm{\Δ}/\mathrm{\λ}},...,1,...,e^{-j2\mathrm{\π}M\sin {\mathrm{\θ}}_p\mathrm{\Δ}/\mathrm{\λ}}]}^T.$

这里，假设n(t)是圆对称零均值高斯分布，其二阶距为 

并且 

其中R_n是N×N的未知正定厄尔米特矩阵，代表噪声的空间相关性，另外 

$\mathrm{\δ}(t_1-t_2)=\left\{\begin{array}{cc}0,& t_1\≠t_2\\ 1,& t_1=t_2\end{array}\right.---\left(7\right)$

四阶累积量： 

传统的阵列处理技术只利用了采样数据的二阶累积量，当信号服从高斯分布时，二阶统计量是充分的，然而，大多数通信系统中通常使用非高 斯信号，例如AQM和BPSK，这时，二阶统计量不能完全描述其统计性质。因此，高阶统计量能提供更多的数据信息。本发明中，用四阶累积矩阵去代替传统的采样协方差矩阵(SCM)，从而消除空间有色高斯噪声。零均值平稳信号的四阶累积量定义如下： 

$\begin{array}{c}\mathrm{cum}(x_{k_1}\left(t\right),x_{k_2}^{*}\left(t\right),x_{l_1}\left(t\right),x_{l_2}^{*}\left(t\right))\\ =E\left\{x_{k_1}\left(t\right)x_{k_2}^{*}\left(t\right)x_{l_1}\left(t\right)x_{l_2}^{*}\left(t\right)\right\}\\ -E\left\{x_{k_1}\left(t\right)x_{k_2}^{*}\left(t\right)\right\}E\left\{x_{l_1}\left(t\right)x_{l_2}^{*}\left(t\right)\right\}\\ -E\left\{x_{k_1}\left(t\right)x_{l_1}\left(t\right)\right\}E\left\{x_{k_2}^{*}\left(t\right)x_{l_2}^{*}\left(t\right)\right\}\\ -E\left\{x_{k_1}\left(t\right)x_{l_2}^{*}\left(t\right)\right\}E\left\{x_{k_2}^{*}\left(t\right)x_{l_1}\left(t\right)\right\}\end{array}\text{---}\left(8\right)$

其中，-M≤k₁,k₂,l₁,l₂≤M。 

定义一个M×M的累积量矩阵它的第(l₁,l₂)个元素为 

$C_{k_1,k_2}(l_1,l_2)=\mathrm{cum}(x_{k_1}\left(t\right),x_{k_2}^{*}\left(t\right),x_{l_1}\left(t\right),x_{l_2}^{*}\left(t\right)).---\left(9\right)$

累积-M≤l₁,l₂≤M范围内的元素得到 

$\begin{array}{c}{C}_{k_1,k_2}=\mathrm{cum}(x_{k_1}\left(t\right),x_{k_2}^{*}\left(t\right),x\left(t\right),x^{*}\left(t\right))\\ =A\left(\mathrm{\θ}\right){\mathrm{\Γ}}_{k_1,k_2}{A}^H\left(\mathrm{\θ}\right)\end{array}\text{---}\left(10\right)$

其中 

是第p个信源的四阶累积量。在-M≤l₁,l₂≤M范围内，总共可以得到(2M+1)²个累积量矩阵，然而，由是中心厄米特矩阵可知 

$C_{k_1,k_2}=C_{k_2,k_1}^H---\left(12\right)$

即和包含相同的数据信息，对应的特征向量和特征值也相同。因此，不需要计算所有(2M+1)²个累积量矩阵，(2M+1)(M+1)个不同的累积量矩阵包含了所有用于DOA估计的数据信息。 

去相干过程： 

当信号高度相关或相干时，秩亏的将导致传统子空间技术性能下降，为了解决这个问题，估计DOA之前进行去相干处理。 

论点1：给定一个四阶累积量其第(l₁,l₂)个元素可以表示为 

$C_{k_1,k_2}(l_1,l_2)=\underset{n=1}{\overset{P}{\mathrm{\Σ}}}{\mathrm{\φ}}_{l_1,n}\·e^{j2\mathrm{\π}{l}_2\sin {\mathrm{\θ}}_n/d}.---\left(13\right)$

其中 

${\mathrm{\φ}}_{l_1,n}=\left\{\begin{array}{c}{\mathrm{\β}}_n^{*}{\mathrm{\γ}}_{4,s_1}\underset{p=1}{\overset{K}{\mathrm{\Σ}}}\underset{q=1}{\overset{K}{\mathrm{\Σ}}}\underset{m=1}{\overset{K}{\mathrm{\Σ}}}{a}_{k_{1}p}{a}_{k_{2}q}^{*}{a}_{l_{1}m}{\mathrm{\β}}_p{\mathrm{\β}}_q^{*}{\mathrm{\β}}_m,n=1,...,K\\ {\mathrm{\γ}}_{4,s_n}\·{\left|{\mathrm{\β}}_n\right|}^4{a}_{k_{1}n}{a}_{k_{2}n}^{*}{a}_{l_{1}n},n=K+1,...,P\end{array}\right.$

证明：详细证明见附录。 

构造如下Toeplitz矩阵： 

其中 

$\tilde{A}=[\tilde{a}\left({\mathrm{\θ}}_1\right),...,\tilde{a}\left({\mathrm{\θ}}_P\right)]---\left(15\right)$

$\begin{array}{c}\tilde{a}\left({\mathrm{\θ}}_P\right)={[1,e^{-j2\mathrm{\π}\sin {\mathrm{\θ}}_p\mathrm{\Δ}/\mathrm{\λ}},...,e^{-j2\mathrm{\π}M\sin {\mathrm{\θ}}_p\mathrm{\Δ}/\mathrm{\λ}}]}^T\\ i=1,...,P.\end{array}---\left(17\right)$

观察(14)式可以发现，是范德蒙矩阵，其向量线性独立。因此，是秩为P的行满秩矩阵，又由论点1得，φ_m,i≠0，i＝1,…,P。所以，无论源信号是否为相干信号，Φ都是满秩对角矩阵。 

DOA估计: 

矩阵可以写成 

${\tilde{C}}_m=\tilde{A}{\widetilde{\mathrm{\Φ}}}_m{\tilde{A}}^H=\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{\mathrm{\φ}}_{m,i}\tilde{a}\left({\mathrm{\θ}}_i\right)\tilde{a}{\left({\mathrm{\θ}}_i\right)}^H---\left(18\right)$

明显地，(18)具有联合对角化结构，和张成相同的列空间 

$\mathrm{pan}\left\{{\tilde{C}}_m\right\}=\mathrm{span}\left\{\tilde{A}\right\}---\left(19\right)$

由C的第-m到m行是共轭对称的，可得其中J是反对角线全为1，其它元素全为0的交换矩阵。因为C_-m和C_m包含相同的数据信息，所以不需要对所有(2M+1)行做Toeplitz转换。这里，我们选择C的前(M+1)行对应的Toeplitz矩阵，即包含了所有的统计信息。最终得到的满秩，(M+1)个Toeplitz矩阵确定阵列流行的列空间并估计对应DOA参数。对于第p个信源，总是存在向量与其余(P-1)个导向矢量张成的列空间正交 

$b_p\⊥\mathrm{range}\{\tilde{a}\left({\mathrm{\θ}}_1\right),...,\tilde{a}\left({\mathrm{\θ}}_{p-1}\right),\tilde{a}\left({\mathrm{\θ}}_{p+1}\right),...,\tilde{a}\left({\mathrm{\θ}}_P\right)\}.---\left(20\right)$

同样可得 

${\tilde{a}}^H\left({\mathrm{\θ}}_i\right)b_p=\left\{\begin{array}{cc}{\tilde{a}}^H\left({\mathrm{\θ}}_i\right)b_p,& i=p\\ 0,& i\≠p\end{array}\right..---\left(21\right)$

将(21)带入(18) 

${\tilde{C}}_m{b}_p=\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{s}_{m,i}\tilde{a}\left({\mathrm{\θ}}_i\right){\tilde{a}}^H\left({\mathrm{\θ}}_i\right)b_p=g_m\tilde{a}\left({\mathrm{\θ}}_p\right).---\left(22\right)$

由(22)式可以确定，如果θ是真实的DOA之一，则总是存在一个标量g(i)使和平行，例如 

$\begin{array}{cc}{\tilde{C}}_{m}b=g_m\tilde{a}\left(\mathrm{\θ}\right),& -M\≤m\≤0.---\left(23\right)\end{array}$

由此引出如下优化问题 

$\begin{array}{cc}\underset{\mathrm{\θ}}{\min }& J(\mathrm{\θ},g,b)=\underset{m=-M}{\overset{0}{\mathrm{\Σ}}}{\left|\right|{\tilde{C}}_{m}b-g_m\widetilde{a\left(\mathrm{\θ}\right)}\left|\right|}^2\\ s.t.& \left|\right|g\left|\right|=1\end{array}---\left(24\right)$

其中a(θ)是参数θ对应的导向矢量，b是(M+1)×1向量，

由于b和g是未知参数,很难对(24)直接搜索DOA，下面简化问题(24)，使之不依赖于b和g，代价函数拓展为 

令 

又因为 $\underset{m=-M}{\overset{0}{\mathrm{\Σ}}}{g}_m={\left|\right|g\left|\right|}^2=1,{\tilde{a}}^H\left(\mathrm{\θ}\right)\tilde{a}\left(\mathrm{\θ}\right)=M+1,$ (25)改写为 

J(θ,g,b)＝b^HFb-b^HG(θ)g-g^HG^H(θ)b+M+1    (28) 

固定θ和g，对b一阶求导如下 

$\frac{\∂J(\mathrm{\θ},g,b)}{\∂b}=2(\mathrm{Fb}-G\left(\mathrm{\θ}\right)g)=0---\left(29\right)$

从而使 

把(30)带入(24)，优化问题简化为 

最小化(31)等效于最大化的特征值分解为 $\underset{i=1}{\overset{M+1}{\mathrm{\Σ}}}{\mathrm{\λ}}_i{u}_i{u}_i^H,$ 其中λ₁≥…≥λ_M+1代表特征值，对应特征向量为于是可得 

最后一个等式成立当且仅当g是最大特征值对应的特征向量，即g＝u₁，λ₁是对应的最大特征值。(31)可进一步改写为 

其中maxeig(…)代表最大特征值 

(2M+1)(M+1)个累积量矩阵包含所有不同的数据信息，每个累积量矩阵可以转换成(M+1)个Toeplitz矩阵。因为所有的Toeplitz矩阵有相同的对角化结构，所以可将代价函数写成 

相干信号的伪输出功率谱为 

本发明还公开了一种信源数未知的相干信号测向系统，包括： 

初始化单元：用于考虑一具有(2M+1)个阵元的均匀线阵，假设有P(P≤M+1)个窄带远场信号从不同方向{θ₁,…,θ_P}入射到该阵列，前K个信号是互相关的，其余(P-K)个信号是互不相关的并且独立于前K个信号， x(t)表示接收样本数据； 

四阶累积量矩阵计算单元：用于计算样本x(t)的(2M+1)(M+1)个四阶累积量矩阵，C_i,j,-M≤i≤j≤M； 

构造Toeplitz矩阵单元：用于对每个累积量矩阵的前(M+1)行构造(M+1)个Toeplitz矩阵

构造中间变量单元：用于计算矩阵和 其中  $a\left({\mathrm{\θ}}_p\right)={[e^{j2\mathrm{\π}M\sin {\text{\θ}}_p\mathrm{\Δ}/\mathrm{\λ}},...,1,...,e^{-j2\mathrm{\π}M\sin {\mathrm{\θ}}_p\mathrm{\Δ}/\mathrm{\λ}}]}^T$ 是导向矢量； 

功率谱计算单元：用于通过公式 计算功率谱； 

搜索单元：用于搜索θ，找到P(θ)的峰值对应的角度即为DOA的估计值。 

在所述初始化单元中：考虑一具有(2M+1)阵元的均匀线阵(ULA)，假设有P(P≤M+1)个远场窄带信号从不同方向{θ₁,…,θ_P}入射到该阵列，前K个信号是互相干的，其余(P-K)信号均是互不相关的且独立于前K个信号，令第一个信号d₁(t)作为参考信号，则第k个相干信号可以表示为： 

$s_k\left(t\right)={\mathrm{\β}}_k{e}^{\mathrm{j\δ}{\mathrm{\φ}}_k}{s}_1\left(t\right),k=2,...,K---\left(1\right)$

其中，β_k幅度衰减因子，δφ_k是相位变化。实际上，δφ_k不会影响信号之间的相干性，不失一般性，令δφ_k＝0,k＝2,…,K，第m个阵元所接收的信号可以表示为 

$\begin{array}{c}{x}_m\left(t\right)=\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{s}_i\left(t\right)e^{-j2\mathrm{\π}m\sin {\mathrm{\θ}}_i\mathrm{\Δ}/\mathrm{\λ}}+n_m\left(t\right)\\ =s_1\left(t\right)\underset{i=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\β}}_i{e}^{-j2\mathrm{\π}m\sin {\mathrm{\θ}}_i\mathrm{\Δ}/\mathrm{\λ}}\\ +\underset{i=K+1}{\overset{P}{\mathrm{\Σ}}}{s}_i\left(t\right)e^{-j2\mathrm{\π}m\sin {\mathrm{\θ}}_i\mathrm{\Δ}/\mathrm{\λ}}+n_m\left(t\right)\end{array}---\left(2\right)$

其中，s_i(t)是第i个信号的复包络，β₁＝1,λ是载波波长，Δ＝λ/2是阵元间间距，(2)式可改写为向量形式： 

x(t)＝[x_-M(t),…,x₀(t),…,x_M(t)]^T

                                 (3) 

＝Ad(t)+n(t) 

其中s(t)＝[s₁(t),…,s_P(t)]^T是源信号向量，A＝[a(θ₁),…,a(θ_P)]是阵列流行，其中，第p个导向矢量可以表示为： 

$a\left({\mathrm{\θ}}_p\right)={[e^{j2\mathrm{\π}M\sin {\text{\θ}}_p\mathrm{\Δ}/\mathrm{\λ}},...,1,...,e^{-j2\mathrm{\π}M\sin {\mathrm{\θ}}_p\mathrm{\Δ}/\mathrm{\λ}}]}^T.---\left(4\right)$

在所述四阶累积量矩阵计算步骤中,零均值平稳信号的四阶累积量定义如下： 

$\begin{array}{c}\mathrm{cum}(x_{k_1}\left(t\right),x_{k_2}^{*}\left(t\right),x_{l_1}\left(t\right),x_{l_2}^{*}\left(t\right))\\ =E\left\{x_{k_1}\left(t\right)x_{k_2}^{*}\left(t\right)x_{l_1}\left(t\right)x_{l_2}^{*}\left(t\right)\right\}\\ -E\left\{x_{k_1}\left(t\right)x_{k_2}^{*}\left(t\right)\right\}E\left\{x_{l_1}\left(t\right)x_{l_2}^{*}\left(t\right)\right\}\\ -E\left\{x_{k_1}\left(t\right)x_{l_1}\left(t\right)\right\}E\left\{x_{k_2}^{*}\left(t\right)x_{l_2}^{*}\left(t\right)\right\}\\ -E\left\{x_{k_1}\left(t\right)x_{l_2}^{*}\left(t\right)\right\}E\left\{x_{k_2}^{*}\left(t\right)x_{l_1}\left(t\right)\right\}\end{array}\text{---}\left(8\right)$

其中，-M≤k₁,k₂,l₁,l₂≤M； 

定义一个M×M的累积量矩阵它的第(l₁,l₂)个元素为 

$C_{k_1,k_2}(l_1,l_2)=\mathrm{cum}(x_{k_1}\left(t\right),x_{k_2}^{*}\left(t\right),x_{l_1}\left(t\right),x_{l_2}^{*}\left(t\right)).---\left(9\right)$

累积-M≤l₁,l₂≤M范围内的元素得到 

$\begin{array}{c}{C}_{k_1,k_2}=\mathrm{cum}(x_{k_1}\left(t\right),x_{k_2}^{*}\left(t\right),x\left(t\right),x^{*}\left(t\right))\\ =A\left(\mathrm{\θ}\right){\mathrm{\Γ}}_{k_1,k_2}{A}^H\left(\mathrm{\θ}\right)\end{array}\text{---}\left(10\right)$

其中 

是第p个信源的四阶累积量。在-M≤l₁,l₂≤M范围内，总共可以得到(2M+1)²个累积量矩阵，然而，由是中心厄米特矩 阵可知 

$C_{k_1,k_2}=C_{k_2,k_1}^H---\left(12\right)$

即和包含相同的数据信息，对应的特征向量和特征值也相同。 

在所述构造Toeplitz矩阵单元中，考虑(2M+1)均匀线阵，利用四阶累积量技术构造(2M+1)(M+1)²个Toeplitz矩阵；在所述构造中间变量单元中，利用(2M+1)(M+1)²个Toeplitz矩阵构造出中间变量F和G(θ)。 

在本发明中，考虑一均匀线阵如图1所示，有5个全向性传感器，阵列间距d＝λ/2，信噪比定义为所有源信号功率与每个传感器接收的加性噪声之比。在实验中，我们假设信源是4QAM调制信号，空间相关噪声服从高斯分布。噪声协方差矩阵的第(k,l)个元素为 

$R_n(k,l)={\mathrm{\σ}}_n^2{\mathrm{\γ}}^{|k-l|}{e}^{\mathrm{j\π}(k-l)/2}---\left(36\right)$

其中，功率电平用于调节信噪比的值，回归系数γ用于调节噪声之间的空间相关性。γ值越大，噪声相关性越大，γ＝0时噪声为高斯白噪声。 

实验1空间谱比较 

在这个实验中，三个功率相同的信号分别从21°,-30°,-10°入射，后两个信号相干。快拍数L＝500。信噪比设为8dB。分别讨论γ＝0.7和γ＝0两种情况，即空间相关高斯噪声和白噪声的情况。主要比较本专利所提出算法的性能和FBSS算法性能，同时，引入同样不需要已知信源数的四阶累积量算法做对比。在FBSS算法中，子阵大小为4。图2为归一化空间谱，通过减去最小的输出功率值再除以最大谱值来实现归一化。图中，黑色的虚线表示真实DOA，可以看出，在噪声空间相关的环境下，本专利所提出的算法正确找出了所有的信源，FOC方法只成功找出一个峰值，无法区分两个相干信号，FBSS算法也可得到三个可见峰值，但与实际DOA值偏差太大。图3为白噪声环境下算法性能的比较，这种情况下，FBSS算法的偏差减小并达到了最佳角度分辨率，本专利算法依然正确找到三个DOA，FOC算法仍然不能正确估计出两个相干信号。 

实验2RMSE与PR随信噪比的关系 

在这个实验中，比较了算法的RMSE随信噪比的关系，同时画出算法对应的分辨率概率(PR)曲线作为参考。对于FBSS算法，ESPRIT-like 算法和EVESPA算法，需要信源数作为先验知识。考虑有三个信源，一个入射角度为20°，另外两个相干信源分别从-42°和-15°入射，快拍数N＝800。设置γ＝0.9，使噪声环境更为严峻，令信噪比从-5dB到30dB增长，FBSS算法中子阵大小为4。通过1000次蒙特卡洛仿真计算出RMSE估计值，定义如下 

$\mathrm{RMSE}=\sqrt{\frac{1}{1000P}\underset{i=1}{\overset{P}{\mathrm{\Σ}}}\underset{j=1}{\overset{1000}{\mathrm{\Σ}}}{({\widehat{\mathrm{\θ}}}_{i,j}-{\mathrm{\θ}}_i)}^2}.---\left(37\right)$

观测图4可得，在低信噪比情况下，本专利算法性能最优，在空间相关噪声和信号相干的情况下，ESPRIT-like算法和FOC算法都不能正常工作，EVESPA算法性能优于FBSS算法却不能随SNR增长而提高。SNR<15dB时，本专利算法性能优于FBSS算法，反之，本专利算法性能低于FBSS算法，这是由于FBSS算法中天线孔径更大。同时，在高信噪比情况下，空间相关噪声对性能的扰动变小，孔径的影响更为重要。当SNR大于20dB，FBSS算法性能超过本专利算法并达到最优。图5是经验概率分辨率与SNR的关系，可以看出，本专利算法的概率分辨率最大，当SNR>15dB时，本专利算法与FBSS算法的概率分辨率都达到最大值。然而，EVESPA算法的概率分辨率较低，通过图3、4中可看出，这是由于其RMSE不随SNR的增长而增大。 

实验3RMSE与PR随快拍数的关系 

固定SNR为10dB，快拍数从100增加到1000，其他参数与实验2相同。图6为RMSE与样本数的关系，图7为相应的概率分辨率与样本数的关系。观察图6，本专利算法性能最优，其次是EVESPA算法，ESPRIT-like算法性能仍旧最差，因为该算法只利用了协方差矩阵中的一小部分信息，并且不能消除噪声的影响。FOC算法在空间相关噪声环境下具有鲁棒性，但由于不能处理相干信号，其性能依然很差。 

实验4RMSE与PR随相关系数的关系 

在这个实验中，考虑有三个信号分别从-10°，8°，35°入射，信噪比和快拍数固定为12dB和800dB。令第一个信号不相关，后两个信号相关。图8是估计DOA的RMSE随第二个信号与第三个信号间相关系数的变化，其中采样信号的相关系数由一阶自回归过程生成： 

$s_3\left(i\right)=\mathrm{\&rho;}{s}_2\left(i\right)+\sqrt{1-{\left|\mathrm{\&rho;}\right|}^2}\&CenterDot;e\left(i\right),i=1,...,N.---\left(38\right)$

本发明和FBSS方法的性能独立于两个信源的相关性，然而，FOSS随ρ的增长而恶化。对比其它四种估计算法，无论信号是否完全相干，本发明在估计精度上都有很大提升。观察图9中概率分辨率与相关系数的关系，当信号不相关或部分相关，本发明与FBSS方法都能成功得到所有的DOA。然而，当信号高度相关甚至相干时，FBSS估计失败的概率高达50％，本发明解决概率仍然可以达到100％。 

下面，我们介绍几个累积量的性质： 

CP1：如果是常数，是随机变量，则 

$\mathrm{cum}({\mathrm{\&alpha;}}_1{x}_1,...,{\mathrm{\&alpha;}}_n{x}_n)=\left(\underset{i=1}{\overset{n}{\mathrm{\&Pi;}}}{\mathrm{\&alpha;}}_1\right)\mathrm{cum}(x_1,...,x_n).---\left(39\right)$

CP2：累积量可加 

cum(x₁+y₁,x₂,…,x_n) 

＝cum(x₁,…,x_n)+cum(y₁,…,y_n).            (40) 

CP3：如果随机变量独立于随机变量

cum(x₁+y₁,…,x_n+y_n) 

＝cum(x₁,…,x_n)+cum(y₁,…,y_n).            (41) 

为了简化符号，令为第n个导向矢量的第m个元素。下面开始证明，第(k,l)个元素的(N×N)四阶累积量矩阵估计值为 

其中 ${\mathrm{\&gamma;}}_{4,s_p}=\mathrm{cum}(s_p\left(t\right),s_p^{*}\left(t\right),s_p\left(t\right),s_p^{*}\left(t\right)).$

定义如下中间变量 

${\mathrm{\&phi;}}_{l_1,n}=\left\{\begin{array}{c}{\mathrm{\&beta;}}_n^{*}{\mathrm{\&gamma;}}_{4,s_1}\underset{p=1}{\overset{K}{\mathrm{\&Sigma;}}}\underset{q=1}{\overset{K}{\mathrm{\&Sigma;}}}\underset{m=1}{\overset{K}{\mathrm{\&Sigma;}}}{a}_{k_{1}p}{a}_{k_{2}q}^{*}{a}_{l_{1}m}{\mathrm{\&beta;}}_p{\mathrm{\&beta;}}_q^{*}{\mathrm{\&beta;}}_m,n=1,...,K\\ {\mathrm{\&gamma;}}_{4,s_n}\&CenterDot;{\left|{\mathrm{\&beta;}}_n\right|}^4{a}_{k_{1}n}{a}_{k_{2}n}^{*}{a}_{l_{1}n},n=K+1,...,P\end{array}\right.$

对于-M≤l₁,l₂≤M，我们有 

论点1证明完毕。 

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

Claims (8)

1.一种信源数未知的相干信号测向方法，其特征在于，包括如下步骤：

初始化步骤：考虑一具有(2M+1)个阵元的均匀线阵，假设有P(P≤M+1)个窄带远场信号从不同方向{θ₁,…,θ_P}入射到该阵列，前K个信号是互相关的，其余(P-K)个信号是互不相关的并且独立于前K个信号，x(t)表示接收样本数据；

四阶累积量矩阵计算步骤：计算样本x(t)的(2M+1)(M+1)个四阶累积量矩阵，C_i,j,-M≤i≤j≤M；

构造Toeplitz矩阵步骤：对每个累积量矩阵的前(M+1)行构造(M+1)个Toeplitz矩阵

构造中间变量步骤：计算矩阵和

其中

$a\left({\mathrm{\&theta;}}_p\right)={[e^{j2\mathrm{\&pi;}M\sin {\text{\&theta;}}_p\mathrm{\&Delta;}/\mathrm{\&lambda;}},...,1,...,e^{-j2\mathrm{\&pi;}M\sin {\mathrm{\&theta;}}_p\mathrm{\&Delta;}/\mathrm{\&lambda;}}]}^T$ 是导向矢量；

功率谱计算步骤：通过公式

计算功率谱；

搜索步骤：搜索θ，找到P(θ)的峰值对应的角度即为DOA的估计值。

2.根据权利要求1所述的相干信号测向方法，其特征在于，在所述初始化步骤中：考虑一具有(2M+1)阵元的均匀线阵(ULA)，假设有P(P≤M+1)个远场窄带信号从不同方向{θ₁,…,θ_P}入射到该阵列，前K个信号是互相干的，其余(P-K)信号均是互不相关的且独立于前K个信号，令第一个信号d₁(t)作为参考信号，则第k个相干信号可以表示为：

$s_k\left(t\right)={\mathrm{\&beta;}}_k{e}^{\mathrm{j\&delta;}{\mathrm{\&phi;}}_k}{s}_1\left(t\right),k=2,...,K---\left(1\right)$

其中，β_k幅度衰减因子，δφ_k是相位变化。实际上，δφ_k不会影响信号之间的相干性，不失一般性，令δφ_k＝0,k＝2,…,K，第m个阵元所接收的信号可以表示为

$\begin{array}{c}{x}_m\left(t\right)=\underset{i=1}{\overset{P}{\mathrm{\&Sigma;}}}{s}_i\left(t\right)e^{-j2\mathrm{\&pi;}m\sin {\mathrm{\&theta;}}_i\mathrm{\&Delta;}/\mathrm{\&lambda;}}+n_m\left(t\right)\\ =s_1\left(t\right)\underset{i=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_i{e}^{-j2\mathrm{\&pi;}m\sin {\mathrm{\&theta;}}_i\mathrm{\&Delta;}/\mathrm{\&lambda;}}\\ +\underset{i=K+1}{\overset{P}{\mathrm{\&Sigma;}}}{s}_i\left(t\right)e^{-j2\mathrm{\&pi;}m\sin {\mathrm{\&theta;}}_i\mathrm{\&Delta;}/\mathrm{\&lambda;}}+n_m\left(t\right)\end{array}---\left(2\right)$

其中，s_i(t)是第i个信号的复包络，β₁＝1,λ是载波波长，Δ＝λ/2是阵元间间距，(2)式可改写为向量形式：

x(t)＝[x_-M(t),…,x₀(t),…,x_M(t)]^T

                                 (3)

＝Ad(t)+n(t)

其中s(t)＝[s₁(t),…,s_P(t)]^T是源信号向量，A＝[a(θ₁),…,a(θ_P)]是阵列流行，其中，第p个导向矢量可以表示为：

$a\left({\mathrm{\&theta;}}_p\right)={[e^{j2\mathrm{\&pi;}M\sin {\text{\&theta;}}_p\mathrm{\&Delta;}/\mathrm{\&lambda;}},...,1,...,e^{-j2\mathrm{\&pi;}M\sin {\mathrm{\&theta;}}_p\mathrm{\&Delta;}/\mathrm{\&lambda;}}]}^T.---\left(4\right)$

3.根据权利要求2所述的相干信号测向方法，其特征在于，在所述四阶累积量矩阵计算步骤中,零均值平稳信号的四阶累积量定义如下：

$\begin{array}{c}\mathrm{cum}(x_{k_1}\left(t\right),x_{k_2}^{*}\left(t\right),x_{l_1}\left(t\right),x_{l_2}^{*}\left(t\right))\\ =E\left\{x_{k_1}\left(t\right)x_{k_2}^{*}\left(t\right)x_{l_1}\left(t\right)x_{l_2}^{*}\left(t\right)\right\}\\ -E\left\{x_{k_1}\left(t\right)x_{k_2}^{*}\left(t\right)\right\}E\left\{x_{l_1}\left(t\right)x_{l_2}^{*}\left(t\right)\right\}\\ -E\left\{x_{k_1}\left(t\right)x_{l_1}\left(t\right)\right\}E\left\{x_{k_2}^{*}\left(t\right)x_{l_2}^{*}\left(t\right)\right\}\\ -E\left\{x_{k_1}\left(t\right)x_{l_2}^{*}\left(t\right)\right\}E\left\{x_{k_2}^{*}\left(t\right)x_{l_1}\left(t\right)\right\}\end{array}\text{---}\left(8\right)$

其中，-M≤k₁,k₂,l₁,l₂≤M；

定义一个M×M的累积量矩阵它的第(l₁,l₂)个元素为

$C_{k_1,k_2}(l_1,l_2)=\mathrm{cum}(x_{k_1}\left(t\right),x_{k_2}^{*}\left(t\right),x_{l_1}\left(t\right),x_{l_2}^{*}\left(t\right)).---\left(9\right)$

累积-M≤l₁,l₂≤M范围内的元素得到

$\begin{array}{c}{C}_{k_1,k_2}=\mathrm{cum}(x_{k_1}\left(t\right),x_{k_2}^{*}\left(t\right),x\left(t\right),x^{*}\left(t\right))\\ =A\left(\mathrm{\&theta;}\right){\mathrm{\&Gamma;}}_{k_1,k_2}{A}^H\left(\mathrm{\&theta;}\right)\end{array}\text{---}\left(10\right)$

其中

是第p个信源的四阶累积量。在-M≤l₁,l₂≤M范围内，总共可以得到(2M+1)²个累积量矩阵，然而，由是中心厄米特矩阵可知

$C_{k_1,k_2}=C_{k_2,k_1}^H---\left(12\right)$

即和包含相同的数据信息，对应的特征向量和特征值也相同。

4.根据权利要求1所述的相干信号测向方法，其特征在于，在所述构造Toeplitz矩阵步骤中，考虑(2M+1)均匀线阵，利用四阶累积量技术构造(2M+1)(M+1)²个Toeplitz矩阵；在所述构造中间变量步骤中，利用(2M+1)(M+1)²个Toeplitz矩阵构造出中间变量F和G(θ)。

5.一种信源数未知的相干信号测向系统，其特征在于，包括：

初始化单元：用于考虑一具有(2M+1)个阵元的均匀线阵，假设有P(P≤M+1)个窄带远场信号从不同方向{θ₁,…,θ_P}入射到该阵列，前K个信号是互相关的，其余(P-K)个信号是互不相关的并且独立于前K个信号，x(t)表示接收样本数据；

四阶累积量矩阵计算单元：用于计算样本x(t)的(2M+1)(M+1)个四阶累积量矩阵，C_i,j,-M≤i≤j≤M；

构造Toeplitz矩阵单元：用于对每个累积量矩阵的前(M+1)行构造(M+1)个Toeplitz矩阵

构造中间变量单元：用于计算矩阵和

其中

$a\left({\mathrm{\&theta;}}_p\right)={[e^{j2\mathrm{\&pi;}M\sin {\text{\&theta;}}_p\mathrm{\&Delta;}/\mathrm{\&lambda;}},...,1,...,e^{-j2\mathrm{\&pi;}M\sin {\mathrm{\&theta;}}_p\mathrm{\&Delta;}/\mathrm{\&lambda;}}]}^T$ 是导向矢量；

功率谱计算单元：用于通过公式

计算功率谱；

搜索单元：用于搜索θ，找到P(θ)的峰值对应的角度即为DOA的估计值。

6.根据权利要求5所述的相干信号测向系统，其特征在于，在所述初始化单元中：考虑一具有(2M+1)阵元的均匀线阵(ULA)，假设有P(P≤M+1)个远场窄带信号从不同方向{θ₁,…,θ_P}入射到该阵列，前K个信号是互相干的，其余(P-K)信号均是互不相关的且独立于前K个信号，令第一个信号d₁(t)作为参考信号，则第k个相干信号可以表示为：

$s_k\left(t\right)={\mathrm{\&beta;}}_k{e}^{\mathrm{j\&delta;}{\mathrm{\&phi;}}_k}{s}_1\left(t\right),k=2,...,K---\left(1\right)$

其中，β_k幅度衰减因子，δφ_k是相位变化。实际上，δφ_k不会影响信号之间的相干性，不失一般性，令δφ_k＝0,k＝2,…,K，第m个阵元所接收的信号可以表示为

$\begin{array}{c}{x}_m\left(t\right)=\underset{i=1}{\overset{P}{\mathrm{\&Sigma;}}}{s}_i\left(t\right)e^{-j2\mathrm{\&pi;}m\sin {\mathrm{\&theta;}}_i\mathrm{\&Delta;}/\mathrm{\&lambda;}}+n_m\left(t\right)\\ =s_1\left(t\right)\underset{i=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_i{e}^{-j2\mathrm{\&pi;}m\sin {\mathrm{\&theta;}}_i\mathrm{\&Delta;}/\mathrm{\&lambda;}}\\ +\underset{i=K+1}{\overset{P}{\mathrm{\&Sigma;}}}{s}_i\left(t\right)e^{-j2\mathrm{\&pi;}m\sin {\mathrm{\&theta;}}_i\mathrm{\&Delta;}/\mathrm{\&lambda;}}+n_m\left(t\right)\end{array}---\left(2\right)$

其中，s_i(t)是第i个信号的复包络，β₁＝1,λ是载波波长，Δ＝λ/2是阵元间间距，(2)式可改写为向量形式：

x(t)＝[x_-M(t),…,x₀(t),…,x_M(t)]^T

                                    (3)

＝Ad(t)+n(t)

其中s(t)＝[s₁(t),…,s_P(t)]^T是源信号向量，A＝[a(θ₁),…,a(θ_P)]是阵列流行，其中，第p个导向矢量可以表示为：

$a\left({\mathrm{\&theta;}}_p\right)={[e^{j2\mathrm{\&pi;}M\sin {\text{\&theta;}}_p\mathrm{\&Delta;}/\mathrm{\&lambda;}},...,1,...,e^{-j2\mathrm{\&pi;}M\sin {\mathrm{\&theta;}}_p\mathrm{\&Delta;}/\mathrm{\&lambda;}}]}^T.---\left(4\right).$

7.根据权利要求6所述的相干信号测向系统，其特征在于，在所述四阶累积量矩阵计算步骤中,零均值平稳信号的四阶累积量定义如下：

$\begin{array}{c}\mathrm{cum}(x_{k_1}\left(t\right),x_{k_2}^{*}\left(t\right),x_{l_1}\left(t\right),x_{l_2}^{*}\left(t\right))\\ =E\left\{x_{k_1}\left(t\right)x_{k_2}^{*}\left(t\right)x_{l_1}\left(t\right)x_{l_2}^{*}\left(t\right)\right\}\\ -E\left\{x_{k_1}\left(t\right)x_{k_2}^{*}\left(t\right)\right\}E\left\{x_{l_1}\left(t\right)x_{l_2}^{*}\left(t\right)\right\}\\ -E\left\{x_{k_1}\left(t\right)x_{l_1}\left(t\right)\right\}E\left\{x_{k_2}^{*}\left(t\right)x_{l_2}^{*}\left(t\right)\right\}\\ -E\left\{x_{k_1}\left(t\right)x_{l_2}^{*}\left(t\right)\right\}E\left\{x_{k_2}^{*}\left(t\right)x_{l_1}\left(t\right)\right\}\end{array}\text{---}\left(8\right)$

其中，-M≤k₁,k₂,l₁,l₂≤M；

定义一个M×M的累积量矩阵它的第(l₁,l₂)个元素为

$C_{k_1,k_2}(l_1,l_2)=\mathrm{cum}(x_{k_1}\left(t\right),x_{k_2}^{*}\left(t\right),x_{l_1}\left(t\right),x_{l_2}^{*}\left(t\right)).---\left(9\right)$

累积-M≤l₁,l₂≤M范围内的元素得到

$\begin{array}{c}{C}_{k_1,k_2}=\mathrm{cum}(x_{k_1}\left(t\right),x_{k_2}^{*}\left(t\right),x\left(t\right),x^{*}\left(t\right))\\ =A\left(\mathrm{\&theta;}\right){\mathrm{\&Gamma;}}_{k_1,k_2}{A}^H\left(\mathrm{\&theta;}\right)\end{array}\text{---}\left(10\right)$

其中

是第p个信源的四阶累积量。在-M≤l₁,l₂≤M范围内，总共可以得到(2M+1)²个累积量矩阵，然而，由是中心厄米特矩阵可知

$C_{k_1,k_2}=C_{k_2,k_1}^H---\left(12\right)$

即和包含相同的数据信息，对应的特征向量和特征值也相同。

8.根据权利要求5所述的相干信号测向系统，其特征在于，在所述构造Toeplitz矩阵单元中，考虑(2M+1)均匀线阵，利用四阶累积量技术构造(2M+1)(M+1)²个Toeplitz矩阵；在所述构造中间变量单元中，利用(2M+1)(M+1)²个Toeplitz矩阵构造出中间变量F和G(θ)。