CN102694588B - 一种基于共轭扩展的波达方向估计方法 - Google Patents

Abstract

本发明公开了一种基于共轭扩展的波达方向估计方法，该方法先计算接收信号之间的二阶统计量，利用延迟自相关函数的性质得到导向矩阵的共扼矩阵，然后利用扩展矢量的延迟采样构造出一个新的伪快拍矩阵，最后对这个矩阵用四阶累积量的方法来估计信号的波达方向。通过阵列扩展的信号波达方向估计方法，本发明能够大大增加估计的信号波达方向数目，并且能够提高信号波达方向的估计精度。

Description

一种基于共轭扩展的波达方向估计方法

技术领域

本发明涉及一种信号波达方向(DOA)估计技术，利用智能天线对无线接收信号的波达方向进行估计，从而可以为网络定位、节点管理、商业服务与位置相关的信息服务提供技术支持。

背景技术

阵列信号处理中的波达方向(DOA)估计在雷达、无线通信等领域有着广泛的应用。其中，基于子空间的MUSIC算法对阵列信号的波达方向进行估计，开辟了阵列信号处理的新时代。该MUSIC算法是一种基于二阶统计量的阵列信号处理方法，二阶统计量的方法通常假定待估计的信号是高斯分布的，但是在实际应用中，信源产生的信号通常是非高斯分布的，使得估计的波达方向具有较大的误差。其中，高阶统计量是指高于二阶的统计量，它具有对任何形式的高斯过程的不可以有效地提取出信号的非高斯成份，最大程度地抑制高斯噪声。因此，高阶统计量常常被用来处理非高斯信号DOA估计问题。由于对称分布的随机信号的三阶累积量为零，因而阵列测向中常用四阶累积量。

目前对于波达方向具有多种估计方法，中国专利CN1523372A公开了一种基于TD-SCDMA标准实现圆-角定位的来波方向估计方法，提出了一种虚拟线阵四阶累量MUSIC算法以及虚拟线阵四阶累量修正MUSIC算法，其计算步骤如下：

1）将基于TD-SCDMA无线定位阵元空间均匀圆阵转换为模式空间虚拟均匀线阵；

2）构成虚拟线阵采样协方差矩阵；

3）对上述协方差矩阵进行特征值分解，求得协方差矩阵的特征值，并计算与各特征值对应的特征向量，从而构成信号子空间和噪声子空间；

3）采样MDL(Minimum Description Length Criterion,MDL)方法估计信号源个数；

4）计算噪声子空间的投影矩阵，并估计MUSIC的功率谱；

5）通过搜索MUSIC的功率谱的谱峰值，估计来波方向。

该方法不仅能够解决TD-SCDMA的圆阵智能天线的DOA估计问题，而且也能够解决TD-SCDMA对相关信号源的DOA估计问题。但是，该方法难以估计出更多的信号到达方向，而且其对信号到达方向的估计精度和分辨率也较低。

为了提高DOA的估计精度和分辨率，MUSIC-like算法利用阵列的输出来构造四阶累积量矩阵，并对该矩阵采用子空间方法进行方向估计。四阶累积量在阵列信号处理中的应用，明确提出了四阶累积量具有扩展阵列有效口径的能力，解释了阵列扩展的物理意义。由于该方法的阵元或虚拟阵元有限，因此其对信号波达方向的估计能力有限；而且，该方法的阵列有效口径不足，导致该方法对信号方向的估计精度也较低。

发明内容

针对现有技术的不足，本发明的目的是要提供一种基于共轭扩展的波达方向估计方法，该方法能够在有限的无线传感器数目下，估计出更多的信号到达方向，提高信号到达方向精度和提高分辨力。

为实现上述目的，本发明采用的技术方案如下：

一种基于共轭扩展的波达方向估计方法，本方法首先对接收的阵列信号进行二阶统计量处理，利用延迟自相关函数的性质得到导向矩阵的共扼矩阵，然后利用扩展矢量的延迟采样构造出一个新的伪快拍矩阵，最后对这个矩阵用四阶累积量的方法来估计信号的波达方向；

该方法具体步骤为：

1）根据第m个天线阵列的输出信号X_m(t)，利用延迟自相关的性质，按照公式：

$R_{x_m{x}_1}\left(\mathrm{\τ}\right)=\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{R}_{s_d{s}_d}\left(\mathrm{\τ}\right)e^{-j2\mathrm{\π}((a_m-a_1)\cos \left({\mathrm{\θ}}_d\right)+(b_m-b_1)\sin \left({\mathrm{\θ}}_d\right))/\mathrm{\λ}}+R_{n_m{n}_1}$

构造矩阵R_x(τ)，R′_x(τ)；

2）根据 $R\left(\mathrm{\τ}\right)=\left[\begin{array}{c}{R}_x\left(\mathrm{\τ}\right)\\ {\left(R_x^{\′}\left(\mathrm{\τ}\right)\right)}^{*}\end{array}\right],$ 对R(τ)在伪快拍数N_p个延迟τ(τ=T_s,2T_s,...N_pT_s)上采样，得到一个伪快拍矩阵

$\tilde{R}=[R\left(T_S\right),R\left(2T_S\right),...R\left(N_p{T}_S\right)];$

3）对所述伪快拍矩阵做四阶累积量计算得到一个四阶累积量矩阵C，其中：

$C=E\left\{(\tilde{R}\⊗{\tilde{R}}^{*}){(\tilde{R}\⊗{\tilde{R}}^{*})}^H\right\}-E\{\tilde{R}\⊗{\tilde{R}}^{*}\}\·E\left\{{(\tilde{R}\⊗{\tilde{R}}^{*})}^H\right\}-E\left\{\tilde{R}{\tilde{R}}^H\right\}\⊗E\left\{{\left(\tilde{R}{\tilde{R}}^H\right)}^{*}\right\};$

4）将所述四阶累积量矩阵C演化成以下形式：

$C=\stackrel{\‾}{B}{\mathrm{\Σ}}_s{\stackrel{\‾}{B}}^H;$

其中，

∑_s为对角线上的元素为R_s(τ)的峰值的对角矩阵；

得到导向矩阵

$\stackrel{\‾}{B}=[\stackrel{\‾}{b}\left({\mathrm{\θ}}_1\right),\·\·\·\stackrel{\‾}{b}\left({\mathrm{\θ}}_D\right)],$ i=1,2,3...；

5）对所述四阶累积量矩阵C进行特征值分解，得到特征值

$\tilde{C}=U_S{\mathrm{\Σ}}_S{U}_S^H+U_N{\mathrm{\Σ}}_N{U}_N^H;$

其中，

U_s是由大特征值对应的特征矢量张成的子空间，即信号子空间；

U_N是由小特征值对应的特征矢量张成的子空间，即噪声子空间；

6）对谱估计式P_MUSIC(θ)进行谱峰搜索，其中，所述谱估计式P_MUSIC(θ)为：

$P_{\mathrm{MUSIC}}\left(\mathrm{\θ}\right)=\frac{1}{b^H\left(\mathrm{\θ}\right){\hat{U}}_N{\hat{U}}_N^{H}b\left(\mathrm{\θ}\right)};$

该谱估计式P_MUSIC(θ)的谱峰对应的横坐标值即为信号的波达方向。

在上述技术方案中，步骤2）中所述导向矩阵还可以通过以下方式得到：

设该方法的初始条件为：信号数D，信号入射方向为θ_i(i＝1,2,…D)，阵元数为M，其中一阵元的坐标为(a_m,b_m)，得到：

$a\left({\mathrm{\θ}}_i\right)=\left[\begin{array}{c}1\\ \exp (-j\frac{2\mathrm{\π}}{\mathrm{\λ}}(a_2\cos \left({\mathrm{\θ}}_i\right)+b_2\sin \left({\mathrm{\θ}}_i\right)))\\ .\\ .\\ .\\ \exp (-j\frac{2\mathrm{\π}}{\mathrm{\λ}}(a_M\cos \left({\mathrm{\θ}}_i\right)+b_M\sin \left({\mathrm{\θ}}_i\right)))\end{array}\right]$

以及，

$a^{\′}\left({\mathrm{\θ}}_i\right)=\left[\begin{array}{c}\exp (-j\frac{2\mathrm{\π}}{\mathrm{\λ}}(a_2\cos \left({\mathrm{\θ}}_i\right)+b_2\sin \left({\mathrm{\θ}}_i\right)))\\ \exp (-j\frac{2\mathrm{\π}}{\mathrm{\λ}}(a_3\cos \left({\mathrm{\θ}}_i\right)+b_3\sin \left({\mathrm{\θ}}_i\right)))\\ .\\ .\\ .\\ \exp (-j\frac{2\mathrm{\π}}{\mathrm{\λ}}(a_M\cos \left({\mathrm{\θ}}_i\right)+b_M\sin \left({\mathrm{\θ}}_i\right)))\end{array}\right]$

根据公式：

$\stackrel{\‾}{b}\left({\mathrm{\θ}}_i\right)=\left[\begin{array}{c}a\left({\mathrm{\θ}}_i\right)\\ {\left(a^{\′}\left({\mathrm{\θ}}_i\right)\right)}^{*}\end{array}\right]\⊗{\left[\begin{array}{c}a\left({\mathrm{\θ}}_i\right)\\ {\left(a^{\′}\left({\mathrm{\θ}}_i\right)\right)}^{*}\end{array}\right]}^{*},$ (i=1,2,…D)

得到导向矩阵

$\stackrel{\‾}{B}=[\stackrel{\‾}{b}\left({\mathrm{\θ}}_1\right),\·\·\·\stackrel{\‾}{b}\left({\mathrm{\θ}}_D\right)].$

上述优选技术方案中，所述确定要估计的参数数目D，以及所述需要选择的参数伪快拍数N_p，均根据系统实际情况确定。

本发明的有益效果如下：

1、由于估计信号方向数与阵列形状密切相关，本发明所述的估计方法对于均匀线阵，M个实阵元可以扩展成4M-3个虚拟阵元。对于最佳阵列(如均匀圆阵)，M个实阵元可以扩展成M²-2M+1个虚拟阵元，因此本发明所述的估计方法具有更多的虚拟阵元，能够同时估计的信号方向数大大提高，极大地降低了对硬件的要求，有效地提高了对波达方向估计的准确性；

2、由于本发明所述的估计方法扩展出了更多的虚拟阵元，孔径增大，使得阵元间的互耦效应减弱，在信噪比比较低的情况下，本方明所述的估计方法依然可以得到很好的效果；

3、由于扩展出更多的虚拟阵元，阵列的孔径增大，孔径越大，角度估计的分辨率也就越高，所以估计的信号波达方向精度也更高；

4、通过实验分析，本发明所述的估计方法可以有效减少所需的快拍数。

为了充分地了解本发明的目的、特征和效果，以下将结合附图对本发明的构思、具体结构及产生的技术效果作进一步说明。

附图说明

图1是线性阵列结构示意图；

图2是本发明所述的估计方法的流程图；

图3是本发明所述估计方法的计算结果示意图；

具体实施方式

如图1和图2所示，本发明公开了一种基于共轭扩展的波达方向估计方法，首先对接收的阵列信号进行二阶统计量处理，利用延迟自相关函数的性质得到导向矩阵的共扼矩阵，然后利用扩展矢量的延迟采样构造出一个新的伪快拍矩阵，最后对这个矩阵用四阶累积量的方法来估计信号的波达方向。

假设天线阵列为均匀线阵，每个阵元都是一致全向，阵元数为M，阵元间距为d。假如空间有D个远场窄带非相干信号分别从θ_i(i＝1,2,…D)入射到天线阵列中。令(a_m,a_m)表示各阵元在坐标轴中的位置，其中(a₁,a₁)=(0,0)为参考阵元。其中，确定要估计的参数数目D，根据系统实际情况确定。由于空间信号相互独立信号与噪声统计独立，因此第m个阵列输出X_m(t)为：

$X_m\left(t\right)=\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{S}_d\left(t\right)\exp \{-j\frac{2\mathrm{\π}}{\mathrm{\λ}}(a_m\cos \left({\mathrm{\θ}}_d\right)+b_m\sin \left({\mathrm{\θ}}_d\right))\}+n_m\left(t\right);$ 其中，

X_m(t)为第m个天线的输出信号；S_d(t)为第d个信号；

n_m(t)为高斯白噪声；λ表示波长；

θ为第d个信号的波达方向

将上式改写成矢量形式为：

X(t)=A(θ)S(t)+N(t)；

其中A为导向矩阵。

本发明所述的估计方法具体步骤如下：

（1）根据第m个天线阵列的输出信号X_m(t)，构造矩阵R_x(τ)，R′_x(τ)：

如果第m个天线阵列的输出信号X_m(t)，则第m个阵元输出X_m(t)与第一个阵元输出X₁(t)间的延迟互相关函数为：

$R_{x_m{x}_1}\left(\mathrm{\τ}\right)=E\left\{x_m(t+\mathrm{\τ}/2)x_1^{*}(t-\mathrm{\τ}/2)\right\}=\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{R}_{s_d{s}_d}\left(\mathrm{\τ}\right)e^{-j2\mathrm{\π}((a_m-a_1)\cos \left({\mathrm{\θ}}_d\right)+(b_m-b_1)\sin \left({\mathrm{\θ}}_d\right))/\mathrm{\λ}}+R_{n_m{n}_1};$

其中，τ是延迟时间并且τ≠0；

如果用来表示s_d(t)的延迟自相关函数，来表示噪声的延迟互相关函数，根据上式有：

$R_{s_d{s}_d}\left(\mathrm{\τ}\right)=E\left\{s_d(t+\mathrm{\τ}/2)s_d^{*}(t-\mathrm{\τ}/2)\right\};$

$R_{n_m{n}_1}\left(\mathrm{\τ}\right)=E\left\{n_m(t+\mathrm{\τ}/2)n_1^{*}(t-\mathrm{\τ}/2)\right\}={\mathrm{\σ}}^2\mathrm{\δ}(m-1)$

其中，σ²是高斯噪声的方差，δ(·)是Dirac函数。

当τ≠0时，因此可以得到矢量形式：

R_x(τ)=AR_s(τ)，R′_x(τ)=A′R_s(τ)；

其中：

$R_x\left(\mathrm{\τ}\right)={[R_{s_1{s}_1}\left(\mathrm{\τ}\right),R_{s_2{s}_1}\left(\mathrm{\τ}\right),\·\·\·R_{s_M{s}_1}\left(\mathrm{\τ}\right)]}^T;$

$R_x^{\′}\left(\mathrm{\τ}\right)={[R_{x_2{x}_1}\left(\mathrm{\τ}\right),R_{x_3{x}_1}\left(\mathrm{\τ}\right),\·\·\·R_{s_M{s}_1}\left(\mathrm{\τ}\right)]}^T;$

A'=[a'(θ₁)，...a′(θ_D)]；

$a^{\′}\left({\mathrm{\θ}}_i\right)=\left[\begin{array}{c}\exp (-j\frac{2\mathrm{\π}}{\mathrm{\λ}}(a_2\cos \left({\mathrm{\θ}}_i\right)+b_2\sin \left({\mathrm{\θ}}_i\right)))\\ \exp (-j\frac{2\mathrm{\π}}{\mathrm{\λ}}(a_3\cos \left({\mathrm{\θ}}_i\right)+b_3\sin \left({\mathrm{\θ}}_i\right)))\\ .\\ .\\ .\\ \exp (-j\frac{2\mathrm{\π}}{\mathrm{\λ}}(a_M\cos \left({\mathrm{\θ}}_i\right)+b_M\sin \left({\mathrm{\θ}}_i\right)))\end{array}\right];$

（2）根据 $R\left(\mathrm{\τ}\right)=\left[\begin{array}{c}{R}_x\left(\mathrm{\τ}\right)\\ {\left(R_x^{\′}\left(\mathrm{\τ}\right)\right)}^{*}\end{array}\right],$ 对R(τ)在伪快拍数N_p个延迟τ(τ=T_s,2T_s，...N_pT_s)上采样，得到一个伪快拍矩阵

$\tilde{R}=[R\left(T_S\right),R\left(2T_S\right),...R\left(N_p{T}_S\right)].$

根据自相关函数的性质：

${\left(R_x^{\′}(-\mathrm{\τ})\right)}^{*}={\left(A^{\′}\right)}^{*}{R}_s^{*}(-\mathrm{\τ})={\left(A^{\′}\right)}^{*}{R}_s\left(\mathrm{\τ}\right);$

可以构造出一个新的矢量矩阵R(τ)：

$R\left(\mathrm{\τ}\right)=\left[\begin{array}{c}{R}_x\left(\mathrm{\τ}\right)\\ {\left(R_x^{\′}(-\mathrm{\τ})\right)}^{*}\end{array}\right]=\left[\begin{array}{c}A\\ {\left(A^{\′}\right)}^{*}\end{array}\right]R_s\left(\mathrm{\τ}\right)$

与对阵元输出X(t)在各个时刻t上做快拍类似，此处对R(τ)在N_p个延迟τ(τ=T_s,2T_s,…N_pT_s)上采样，得到一个伪快拍矩阵

$\tilde{R}=[R\left(T_s\right),R\left(2T_s\right),\·\·\·R\left(N_p{T}_s\right)].$

（3）对做四阶累积量计算得到一个四阶累积量矩阵C：根据上述伪快拍矩阵可以得到一个四阶累积量矩阵C，

$C=E\left\{(\tilde{R}\⊗{\tilde{R}}^{*}){(\tilde{R}\⊗{\tilde{R}}^{*})}^H\right\}-E\{\tilde{R}\⊗{\tilde{R}}^{*}\}\·E\left\{{(\tilde{R}\⊗{\tilde{R}}^{*})}^H\right\}-E\left\{\tilde{R}{\tilde{R}}^H\right\}\⊗E\left\{{\left(\tilde{R}{\tilde{R}}^H\right)}^{*}\right\};$

（4）将所述四阶累积量矩阵C演化成以下形式：

$C=\stackrel{\‾}{B}{\mathrm{\Σ}}_s{\stackrel{\‾}{B}}^H;$

其中，

∑_s为对角矩阵，其对角线上的元素为R_s(τ)的峰值；

因此可以得到导向矩阵为：

$\stackrel{\‾}{B}=[\stackrel{\‾}{b}\left({\mathrm{\θ}}_1\right),\·\·\·\stackrel{\‾}{b}\left({\mathrm{\θ}}_D\right)];$

住的注意的是，该导向矩阵也可以采用以下方式得到：

设该方法的初始条件为：信号数D，信号入射方向为θ_i(i＝1,2,…D)，阵元数为M，其中一阵元的阵元坐标为(a_m,b_m)，得到：

$a\left({\mathrm{\θ}}_i\right)=\left[\begin{array}{c}1\\ \exp (-j\frac{2\mathrm{\π}}{\mathrm{\λ}}(a_2\cos \left({\mathrm{\θ}}_i\right)+b_2\sin \left({\mathrm{\θ}}_i\right)))\\ .\\ .\\ .\\ \exp (-j\frac{2\mathrm{\π}}{\mathrm{\λ}}(a_M\cos \left({\mathrm{\θ}}_i\right)+b_M\sin \left({\mathrm{\θ}}_i\right)))\end{array}\right]$

以及，

$a^{\′}\left({\mathrm{\θ}}_i\right)=\left[\begin{array}{c}\exp (-j\frac{2\mathrm{\π}}{\mathrm{\λ}}(a_2\cos \left({\mathrm{\θ}}_i\right)+b_2\sin \left({\mathrm{\θ}}_i\right)))\\ \exp (-j\frac{2\mathrm{\π}}{\mathrm{\λ}}(a_3\cos \left({\mathrm{\θ}}_i\right)+b_3\sin \left({\mathrm{\θ}}_i\right)))\\ .\\ .\\ .\\ \exp (-j\frac{2\mathrm{\π}}{\mathrm{\λ}}(a_M\cos \left({\mathrm{\θ}}_i\right)+b_M\sin \left({\mathrm{\θ}}_i\right)))\end{array}\right]$

根据公式：

$\stackrel{\‾}{b}\left({\mathrm{\θ}}_i\right)=\left[\begin{array}{c}a\left({\mathrm{\θ}}_i\right)\\ {\left(a^{\′}\left({\mathrm{\θ}}_i\right)\right)}^{*}\end{array}\right]\⊗{\left[\begin{array}{c}a\left({\mathrm{\θ}}_i\right)\\ {\left(a^{\′}\left({\mathrm{\θ}}_i\right)\right)}^{*}\end{array}\right]}^{*},$ (i=1,2,…D)

得到导向矩阵

（5）对所述四阶累积量C进行特征值分解，得到特征值

$\tilde{C}=U_S{\mathrm{\Σ}}_S{U}_S^H+U_N{\mathrm{\Σ}}_N{U}_N^H;$

其中，

U_s是由大特征值对应的特征矢量张成的子空间，即信号子空间；

U_N是由小特征值对应的特征矢量张成的子空间，即噪声子空间；

从而得到噪声子空间U_N。

（6）根据信号子空间U_s与噪声子空间U_N是互相正交的，即信号子空间中的导向矢量也与噪声子空间正交，因此可以得到谱估计式P_MUSIC(θ)：

$P_{\mathrm{MUSIC}}=\frac{1}{{\stackrel{\‾}{b}}^H\left(\mathrm{\θ}\right){\hat{U}}_N{\hat{U}}_N^H\stackrel{\‾}{b}\left(\mathrm{\θ}\right)};$

对该谱估计式P_MUSIC(θ)进行谱峰搜索就可以计算出信号的波达角，该谱估计式P_MUSIC(θ)的谱峰对应的横坐标值即为信号的波达方向。

基于共轭矩阵的阵列扩展技术能够得到更多的虚拟矩阵，也就是说相同阵元数，利用本发明方法能够估计出更多的信号波达方向，而且精度也会更高。关于到底能够估计出多少信号波达方向，有以下的定理。对于M个实阵元，利用四阶累积量方法来进行阵列扩展，本发明所用的方法其扩展阵列阵元数最多为2M²-2M+1元、均匀线阵的扩展阵列元数为4M-3。

以下通过具体数据的计算过程，对本发明所述的估计方法进行进一步的阐述和验证。

在天线阵列为均匀线阵的情况下，假设初始条件为:阵元数M=3，阵元间距为半波长，信噪比为10dB，快拍数N_s=1000，伪快拍数N_p=200。

当有6个远场窄带信号（远场窄带信号和噪声都是随机产生）从40°、57°、85°、103°、123°和146°方向入射到天线阵列上，通过本发明所述的估计方法对信号的波达方向进行估计。由于本发明所述的估计方法可以将3个均匀线阵的实阵元扩展成4*3-3=9个阵元，所以最多能估计出8个信号波达方向。

1）根据延迟自相关函数，阵列输出可以构造出矩阵

$R_x\left(\mathrm{\τ}\right)=1.0e+003*{\left[\begin{array}{ccc}0& \·\·\·& 3.6152-0.0008i\\ 0& \·\·\·& -0.1074-0.1243i\\ 0& \·\·\·& -0.1965-0.3422i\end{array}\right]}_{3\×200};$

$R_x^{\′}\left(\mathrm{\τ}\right)=1.0e+003*{\left[\begin{array}{ccc}0& \·\·\·& -0.0991-0.1266i\\ 0& \·\·\·& -0.1936-0.3454i\end{array}\right]}_{2\×200};$

2）根据 $R\left(\mathrm{\τ}\right)=\left[\begin{array}{c}{R}_x\left(\mathrm{\τ}\right)\\ {\left(R_x^{\′}\left(\mathrm{\τ}\right)\right)}^{*}\end{array}\right],$ 对R(τ)在伪快拍数N_p个延迟τ(τ=T_s,2T_s，...N_pT_s)上采样，得到一个伪快拍矩阵：

$\tilde{R}=1.0e+003*{\left[\begin{array}{ccc}0& \·\·\·& 3.6152-0.0008i\\ 0& \·\·\·& -0.1074-0.1243i\\ 0& \·\·\·& -0.1965-0.3422i\\ 0& \·\·\·& -0.0991+0.1266i\\ 0& \·\·\·& -0.1936+0.3454i\end{array}\right]}_{5\×200}$

3）对所述伪快拍矩阵做四阶累积量计算得到一个四阶累积量矩阵C，其中：

$C=E\left\{(\tilde{R}\⊗{\tilde{R}}^{*}){(\tilde{R}\⊗{\tilde{R}}^{*})}^H\right\}-E\{\tilde{R}\⊗{\tilde{R}}^{*}\}\·E\left\{{(\tilde{R}\⊗{\tilde{R}}^{*})}^H\right\}-E\left\{\tilde{R}{\tilde{R}}^H\right\}\⊗E\left\{{\left(\tilde{R}{\tilde{R}}^H\right)}^{*}\right\};$

计算得到：

$C=1.0e+012*{\left[\begin{array}{ccc}-0.9747& \·\·\·& -0.1025+0.1302i\\ \·& \·\·\·& \·\\ -0.1025-0.1302i& \·\·\·& -1.0761\end{array}\right]}_{10\×10}$

4）计算导向矩阵根据 $\stackrel{\‾}{b}\left({\mathrm{\θ}}_i\right)=\left[\begin{array}{c}a\left({\mathrm{\θ}}_i\right)\\ {\left(a^{\′}\left({\mathrm{\θ}}_i\right)\right)}^{*}\end{array}\right]\⊗{\left[\begin{array}{c}a\left({\mathrm{\θ}}_i\right)\\ {\left(a^{\′}\left({\mathrm{\θ}}_i\right)\right)}^{*}\end{array}\right]}^{*},$ (i=1,2,…D)得到：

$\stackrel{\‾}{B}={\left[\begin{array}{ccc}0.1076-0.9942i& \·\·\·& 0.4710+0.8821i\\ \·& \·\·\·& \·\\ 1.0000& \·\·\·& 1.0000\end{array}\right]}_{10\×1800}$

5）对步骤3）中得到的四阶累积量C进行特征分解，得到对应于噪声子空间的特征向量，

$U_N={\left[\begin{array}{ccc}0.4382+0.0000i& \·\·\·& 0.1241-0.0000i\\ \·& \·\·\·& \·\\ 0.1975-0.2270i& \·\·\·& -0.1281-0.0615i\end{array}\right]}_{10\×4}$

6）构成谱估计式进行谱峰搜索得到波达方向，得到：

P_MUSIC(θ)=[-9.1245…-9.1245]_1×1800

请参阅图3，根据谱估计式的空间谱功率和波达方向作图，可以看出，在谱峰处对应的横坐标即为信号的波达方向。本实施例中，该谱估计式具有六个谱峰，其坐标值分别为(39.68.497)、(57.5 6.33)、(85.1 8.851)、(102.9 10.87)、(123.2 12.58)、(146.317.73)，其对应的横坐标分别为39.6°、57.5°、85.1°、102.9°、123.2°、146.3°。与理论值40°、57°、85°、103°、123°和146°相比较，可以看出该估计方法对此六个波达方向进行了很好的估计，而且准确度很高。

相对于现有技术，本发明所述的估计方法可以将实阵元可以扩展成多个虚拟阵元，对于最佳阵列(如均匀圆阵)，M个实阵元可以扩展成M²-2M+1个虚拟阵元，因此本发明所述的估计方法具有更多的虚拟阵元，极大降低了对硬件的要求，有效提高了对波达方向估计的准确性。

其次，提高了噪声免疫力，由于扩展出更多的虚拟阵元，孔径增大，阵元间的互耦效应减弱，在信噪比比较低的情况下，本方明所述的估计方法依然可以得到很好的效果。

再次，提高了角度分辨率，由于扩展出更多的虚拟阵元，阵列的孔径增大，孔径越大，角度估计的分辨率也就越高，所以估计的信号波达方向精度也更高；最后，所需快拍数，通过实验分析，本发明所述的估计方法可以有效减少所需的快拍数。

以上详细描述了本发明的较佳具体实施例，应当理解，本领域的普通技术无需创造性劳动就可以根据本发明的构思做出诸多修改和变化。因此，凡本技术领域中技术人员依本发明构思在现有技术基础上通过逻辑分析、推理或者根据有限的实验可以得到的技术方案，均应该在由本权利要求书所确定的保护范围之中。

Claims (3)

1.一种基于共轭扩展的波达方向估计方法，其特征在于，该方法首先对接收的阵列信号进行二阶统计量处理，进行延迟采样构造出一伪快拍矩阵，最后对该伪快拍矩阵用四阶累积量的方法来估计信号的波达方向； 

该方法具体步骤为： 

1)假设天线阵列为均匀线阵，每个阵元都是一致全向，阵元数为M，阵元间距为d，假如D个远场窄带非相干信号分别以θ_i，i＝1,2,...D投射到M个天线阵列上，根据第m个天线阵列的输出信号X_m(t)，其中 

其中X_m(t)为第m个天线的输出信号；S_d(t)为第d个信号；n_m(t)为高斯白噪声；λ表示波长；θ为第d个信号的波达方向； 

利用延迟自相关的性质，按照公式 

构造矩阵R_x(τ)，R'_x(τ)；其中表示第m个阵元输出X_m(t)与第一个阵元输出X₁(t)间的延迟互相关函数，表示s_d(t)的延迟自相关函数，表示噪声的延迟互相关函数； 

2)根据对R(τ)在伪快拍数N_p个延迟τ＝T_s,2T_s,...N_pT_s上采样，得到一个伪快拍矩阵

3)对所述伪快拍矩阵做四阶累积量计算得到一个四阶累积量矩阵C，其中： 

4)将所述四阶累积量矩阵C演化成以下形式： 

其中， 

∑_s为对角线上的元素为R_s(τ)的峰值的对角矩阵； 

得到导向矩阵

其中所述步骤4)中的导向矩阵通过以下方式得到： 

设该方法的初始条件为：信号数D，信号入射方向为θ_i，阵元数为M，其中一阵元的坐标为(a_m,b_m)，得到： 

以及， 

根据公式： 

得到导向矩阵

5)对所述四阶累积量矩阵C进行特征值分解，得到特征值

其中，U_S是由大特征值对应的特征矢量张成的子空间，即信号子空间； 

U_N是由小特征值对应的特征矢量张成的子空间，即噪声子空间； 

6)对谱估计式P_MUSIC(θ)进行谱峰搜索，其中，所述谱估计式P_MUSIC(θ)为： 

该谱估计式的P_MUSIC(θ)谱峰对应的横坐标值即为信号的波达方向。 

2.如权利要求1所述的基于共轭扩展的波达方向估计方法，其特征在于：所述信号数D，根据系统实际情况确定。 

3.如权利要求1所述的基于共轭扩展的波达方向估计方法，其特征在于：所述伪快拍数N_p，根据系统实际情况确定。