阵列信号处理广泛地应用于雷达、声纳、无线通信、地震学等领域，波达方向估计(空间谱估计)作为阵列信号处理中的一个重要研究方向，近几十年来吸引了大批学者进行研究，以最大似然算法(ML)^[1]为代表的子空间拟合类算法估计精度最高，但该算法需要求解一个多维非线性优化问题，运算量相当大，不利于工程应用。

子空间类算法中MUSIC算法^[2]以分辨力强、稳定性好以及在一定条件下能够逼近克拉美-罗界^[7]而得到了很大的重视，但是需要谱峰搜索。而作为MUSIC算法的孪生算法，ESPRIT算法^[3]应用信号子空间的旋转不变特性得到波达方向估计的闭式解，避免了计算量大的问题。但是ESPRIT算法仅利用到信号子空间旋转不变关系矩阵的广义特征值，而由于噪声的影响导致自相关矩阵估计存在误差，这对特征值的估计造成了波动^[4]。

针对传统ESPRIT算法估计方差大的问题，本文利用波达方向与广义特征值向量的关系进行求解，通过关系矩阵与阵列流形的关系，进行信号波达方向估计。该方法计算量低，可以得到多个估计值，将其平均以达到对波动抑制的效果，最终达到与MUSIC算法逼近的估计性能。

1 信号模型

假设M个天线构成的均匀直线阵列，K个统计独立的远场窄带信号分别从方向${\theta _k}, k = 1, 2, \cdots, K$到达阵列，阵列在采样时刻t的接收向量为：

$ {\boldsymbol{x}}(t) = \sum\limits_{k = 1}^K {{{\boldsymbol{a}}_k}{s_k}(t)} +{\boldsymbol{v}}(t) $

(1)

式中，${{\boldsymbol{a}}_k} = {\left[{1, \exp ({\rm{j}}{\phi _k}), \cdot \cdot \cdot, \exp ({\rm{j}}(M-1){\phi _k})} \right]^{\rm{T}}}$是维度M×1的驱动向量，其中${\phi _k}{\rm{ = - 2 \mathsf{ π} }}d\sin \left( {{\theta _k}} \right)/\lambda $，d为阵元间距，λ为信号波长；${s_k}(t)$为对应${\theta _k}$方向的信号；${\boldsymbol{v}}(t) \in {C^{M \times 1}}$假设为高斯白噪声。假设快拍数为L，可以将式(1)等效为矩阵表达形式：

$\mathit{\boldsymbol{X}} = \mathit{\boldsymbol{AS}} + \mathit{\boldsymbol{V}}$

(2)

式中，$\mathit{\boldsymbol{X}} \in {C^{M \times L}}$；阵列流形$\mathit{\boldsymbol{A}} = [{\mathit{\boldsymbol{a}}_1},{\mathit{\boldsymbol{a}}_2}, \cdots ,{\mathit{\boldsymbol{a}}_K}] $；信号$\mathit{\boldsymbol{S}} \in {C^{K \times L}}$；噪声$\mathit{\boldsymbol{V}} \in {C^{M \times L}}$。假设信号与噪声独立，即$E(\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{V}}^{\rm{H}}}) = 0$，其中0是$K \times M$维的零矩阵。得到自相关矩阵：

$\mathit{\boldsymbol{R}} = E(\mathit{\boldsymbol{X}}{\mathit{\boldsymbol{X}}^{\rm{H}}}) = \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{R}}_S}{\mathit{\boldsymbol{A}}^{\rm{H}}} + \sigma _n^2{\mathit{\boldsymbol{I}}_M} $

(3)

式中，${\mathit{\boldsymbol{R}} _S}$为信号的自相关矩阵；$\sigma _n^2$为噪声的功率；${\mathit{\boldsymbol{I}} _M}$为$M \times M$的单位矩阵。但在实际中只能通过样本对自相关矩阵进行估计$\hat {\mathit{\boldsymbol{R}}} $：

$ \mathit{\boldsymbol{R}} \approx \mathit{\boldsymbol{\hat R}} = \frac{1}{{2L}}\sum\limits_{t = 1}^L {(\mathit{\boldsymbol{x}}(t){\mathit{\boldsymbol{x}}^{\rm{H}}}(t) + \mathit{\boldsymbol{y}}(t){\mathit{\boldsymbol{y}}^{\rm{H}}}(t))} $

(4)

式中，$\mathit{\boldsymbol{y}} (t)$是$\mathit{\boldsymbol{x}} (t)$的倒序共轭向量。

2 旋转不变特性

ESPRIT算法最基本的假设是存在两个完全相同的两个子阵，针对本文信号模型可选取阵列的前M-1个阵元为子阵1，后M-1个阵元为子阵2，两个子阵接收的数据分别为X₁和X₂：

$\left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{X}}_1} = {\mathit{\boldsymbol{A}}_1}\mathit{\boldsymbol{S}} + {\mathit{\boldsymbol{N}}_1}}\\ {{\mathit{\boldsymbol{X}}_2} = {\mathit{\boldsymbol{A}}_2}\mathit{\boldsymbol{S}} + {\mathit{\boldsymbol{N}}_2}} \end{array}} \right.$

(5)

由阵列流形的范德蒙特性易得到阵列流形旋转不变性关系：

$ {\mathit{\boldsymbol{A}}_2} = {\mathit{\boldsymbol{A}}_1}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} $

(6)

式中，$\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} = {\rm{diag}}(\exp ({\rm{j}}{\phi _1}), \exp ({\rm{j}}{\phi _2}), \cdots, \exp ({\rm{j}}{\phi _K}))$，是阵列流形旋转不变关系矩阵。${\rm{diag(}} \cdot {\rm{)}}$代表$\mathit{\boldsymbol{\varPhi}}$是一个对角矩阵，对角线上的元素为$\exp ({\rm{j}}{\phi _k})$，$k = 1, 2, \cdots, K$。

在信号不相关的假设下可知，阵列流形与自相关矩阵的信号子空间具有相同的值域。因此存在一个唯一的$K \times K$非奇异矩阵T使得下式成立：

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{U}}_{s1}}\mathit{\boldsymbol{T}} = {\mathit{\boldsymbol{A}}_1}}\\ {{\mathit{\boldsymbol{U}}_{s2}}\mathit{\boldsymbol{T}} = {\mathit{\boldsymbol{A}}_2} = {\mathit{\boldsymbol{A}}_1}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}} \end{array}} \right. $

(7)

式中，${\mathit{\boldsymbol{U}} _{s1}}$、${\mathit{\boldsymbol{U}} _{s2}}$分别对应子阵1和子阵2的信号子空间(也即对应式(8)中${\mathit{\boldsymbol{U}} _{s}}$前M-1行和后M-1行)：

$ \mathit{\boldsymbol{R}} = {\mathit{\boldsymbol{U}}_s}{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_s}\mathit{\boldsymbol{U}}_s^{\rm{H}} + {\mathit{\boldsymbol{U}}_n}{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_n}\mathit{\boldsymbol{U}}_n^{\rm{H}} $

(8)

式中，${\mathit{\boldsymbol{\varSigma}} _s} \in {C^{K \times K}}$和${\mathit{\boldsymbol{\varSigma}} _n} \in {C^{\left({M-K} \right) \times \left({M-K} \right)}}$均为对角矩阵，对角线元素为自相关矩阵的特征值，${\mathit{\boldsymbol{\varSigma}} _s}$中对应最大的K个值。求解式(7)有：

${\mathit{\boldsymbol{U}}_{s2}}\mathit{\boldsymbol{T}} = {\mathit{\boldsymbol{U}}_{s1}}\mathit{\boldsymbol{T \boldsymbol{\varPhi} }} $

(9)

由于T可逆，从而得到信号子空间旋转不变关系^[11]：

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{U}}_{s2}} = {\mathit{\boldsymbol{U}}_{s1}}\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}}\\ {\mathit{\boldsymbol{ \boldsymbol{\varPsi} }} = \mathit{\boldsymbol{T \boldsymbol{\varPhi} }}{\mathit{\boldsymbol{T}}^{ - 1}}} \end{array}} \right. $

(10)

3 基于特征向量的ESPRIT算法 3.1 基于广义特征向量的ESPRIT算法

可以将求解变换矩阵T和阵列流形旋转关系矩阵$ \boldsymbol{\varPhi} $的问题转换为求解$(\mathit{\boldsymbol{U}}_{s1}^{\rm{H}}{\mathit{\boldsymbol{U}}_{s2}}, \mathit{\boldsymbol{U}}_{s1}^{\rm{H}}{\mathit{\boldsymbol{U}}_{s1}})$的广义特征向量和广义特征值求解的问题：

$ \mathit{\boldsymbol{U}}_{s1}^{\rm{H}}{\mathit{\boldsymbol{U}}_{s2}}\mathit{\boldsymbol{T}} = \mathit{\boldsymbol{U}}_{s1}^{\rm{H}}{\mathit{\boldsymbol{U}}_{s1}}\mathit{\boldsymbol{T \boldsymbol{\varPhi} }} $

(11)

记由广义特征值组成的对角矩阵为$\hat {\boldsymbol{\varPhi} } $，传统ESPRIT算法基于广义特征值得到的波达方向估计为：

$ {\hat \theta _k} = {\rm{asin}}\left( { - \frac{{{\rm{angle}}({{\mathit{\boldsymbol{ \boldsymbol{\hat \varPhi} }}}_{kk}})\lambda }}{{2{\rm{ \mathsf{ π} }}d}}} \right)\;\;\;k = 1,2, \cdots ,K$

(12)

式中，angle(·)为相角运算符；${{{\mathit{\boldsymbol{ \boldsymbol{\hat \varPhi} }}}_{kk}}}$是第k个特征值。

但是，由式(7)可见，如果能够得到广义特征向量矩阵T，便可对子阵列流形进行估计：

$ {\mathit{\boldsymbol{A}}_1} \approx {\mathit{\boldsymbol{B}}_1} = {\mathit{\boldsymbol{U}}_{s1}}\mathit{\boldsymbol{TD}} $

(13)

式中，D是一个对角矩阵，反映了广义特征向量的参考相位是不确定的，但是该不确定性对波达方向的估计没有影响。为了减小噪声对阵列流形估计的影响，再将B₁中的所有元素幅度置1得到子阵列流形的最终估计${\hat {\mathit{\boldsymbol{A}} _1}}$：

$ \begin{array}{*{20}{l}} {\;\;\;\;\;{{\mathit{\boldsymbol{\hat A}}}_{1(i,j)}} = \frac{{{\mathit{\boldsymbol{B}}_{1(i,j)}}}}{{{\rm{abs}}({\mathit{\boldsymbol{B}}_{1(i,j)}})}}}\\ {i \in [1,2, \cdots ,M - 1],\;\;j \in [1,2, \cdots ,K]} \end{array} $

(14)

然后，由阵列流形的范德蒙特性有：

$ \begin{array}{*{20}{l}} {\;\;\;\;\frac{{{\mathit{\boldsymbol{A}}_{1(i,j)}}}}{{{\mathit{\boldsymbol{A}}_{1(i - 1,j)}}}} = \exp ({\rm{j}}{\phi _j})}\\ {i \in [2,3, \cdot \cdot \cdot ,M - 1],j \in [1,2, \cdot \cdot \cdot ,K]} \end{array} $

(15)

因此可以将阵列流形的估计带入式(15)，得到M-2个相位估计，将其平均，最终得到广义特征向量波达方向估计：

$ \left\{ {\begin{array}{*{20}{l}} {{{\hat \phi }_k} = \frac{{\mathit{\boldsymbol{\hat A}}_1^{\rm{H}}[1:M - 2,k]{{\mathit{\boldsymbol{\hat A}}}_1}[2:M - 1,k]}}{{M - 2}}}\\ {{{\hat \theta }_k} = {\rm{asin}}\left( { - \frac{{{\rm{angle}}({{\hat \phi }_k})\lambda }}{{2{\rm{ \mathsf{ π} }}d}}} \right)\;\;\;k = 1,2, \cdot \cdot \cdot ,K} \end{array}} \right. $

(16)

3.2 广义特征向量矩阵的求解

广义特征向量的求解可以采用QZ分解^[8]。在式(10)中给出了信号子空间旋转不变关系矩阵$ \boldsymbol{\varPsi }$同阵列流形旋转不变关系矩阵$\boldsymbol{\varPhi } $的关系：$\mathit{\boldsymbol{ \boldsymbol{\varPsi} }} = \mathit{\boldsymbol{T \boldsymbol{\varPhi} }}{\mathit{\boldsymbol{T}}^{ - 1}}$，从而可以将求解广义特征向量和广义特征值的问题转化为求解$ \boldsymbol{\varPsi }$，并通过$ \boldsymbol{\varPsi }$的特征分解得到特征向量和特征值^[12]。

对$ \boldsymbol{\varPsi }$的估计现在存在直接求解和迭代求解两重求解方式，其中直接求解方法包括最小二乘(LS)^[13]、总体最小二乘方法(TLS)^[9]。最小二乘算法的基本思想如下：

$ \min \left\| {\Delta {\mathit{\boldsymbol{U}}_{s2}}} \right\|_{\rm{F}}^2\;\;\;{\rm{s}}.{\rm{t}}.{\mathit{\boldsymbol{U}}_{s1}}\mathit{\boldsymbol{ \boldsymbol{\varPsi} }} = {\mathit{\boldsymbol{U}}_{s2}} + \Delta {\mathit{\boldsymbol{U}}_{s2}} $

(17)

在${\mathit{\boldsymbol{U}}_{s2}}$中引入范数平方最小扰动$\Delta {\mathit{\boldsymbol{U}}_{s2}}$的目的在于校正${\mathit{\boldsymbol{U}}_{s2}}$中存在的噪声。将目标函数对位置参量$ \boldsymbol{\varPsi }$求偏导很容易得到最小二乘估计式：

$ {\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}_{LS}} = {(\mathit{\boldsymbol{U}}_{s1}^{\rm{H}}{\mathit{\boldsymbol{U}}_{s1}})^{ - 1}}\mathit{\boldsymbol{U}}_{s1}^{\rm{H}}{\mathit{\boldsymbol{U}}_{s2}} $

(18)

于是对$(\mathit{\boldsymbol{U}}_{s1}^{\rm{H}}{\mathit{\boldsymbol{U}}_{s2}}, \mathit{\boldsymbol{U}}_{s1}^{\rm{H}}{\mathit{\boldsymbol{U}}_{s1}})$求广义特征向量的过程可以转化为对${(\mathit{\boldsymbol{U}}_{s1}^{\rm{H}}{\mathit{\boldsymbol{U}}_{s1}})^{ - 1}}\mathit{\boldsymbol{U}}_{s1}^{\rm{H}}{\mathit{\boldsymbol{U}}_{s2}}$求特征向量的过程。对比式(18)可以得到求广义特征值与最小二乘求解等效。

总体最小二乘算法针对本文内容指的是在${\mathit{\boldsymbol{U}}_{s1}}$和${\mathit{\boldsymbol{U}}_{s2}}$上分别增加扰动$\Delta {\mathit{\boldsymbol{U}}_{s1}}$、$\Delta {\mathit{\boldsymbol{U}}_{s2}}$对${\mathit{\boldsymbol{U}}_{s1}}$和${\mathit{\boldsymbol{U}}_{s2}}$中的误差进行同时校正。该算法需要先构造矩阵${\mathit{\boldsymbol{U}}_{s12}} = [{\mathit{\boldsymbol{U}}_{s1}}{\kern 1pt} {\kern 1pt} {\mathit{\boldsymbol{U}}_{s2}}] \in {\square}{^{(M -1) \times 2K}}$，再将$\mathit{\boldsymbol{U}}_{s12}^{\rm{H}}{\mathit{\boldsymbol{U}}_{s12}}$进行特征值分解得到$\mathit{\boldsymbol{U}}_{s12}^{\rm{H}}{\mathit{\boldsymbol{U}}_{s12}} = \mathit{\boldsymbol{W \boldsymbol{\varLambda} }}{\mathit{\boldsymbol{W}}^{\rm{H}}}$，并将W按式(18)分成4个$K \times K$的块矩阵，最终得到估计：

$ \mathit{\boldsymbol{W}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{W}}_{11}}}&{{\mathit{\boldsymbol{W}}_{12}}}\\ {{\mathit{\boldsymbol{W}}_{21}}}&{{\mathit{\boldsymbol{W}}_{22}}} \end{array}} \right] $

(19)

$ {{\mathit{\boldsymbol{ \varPsi}}\text{ }}_{\text{TLS}}}=-{{\mathit{\boldsymbol{W}}}_{12}}\mathit{\boldsymbol{W}}_{22}^{-1}\in {{\square }^{K\times K}} $

(20)

迭代求解算法包括STLS算法^[10]、SLS算法^[5]和CTLS算法^[6]等。但这类算法借助高斯迭代或逆迭代求解大大增加了运算量，因此对这类算法不做考虑。

3.3 算法流程及分析

根据上述分析可以得到基于广义特征向量的波达方向估计算法流程如下：

1) 将采样数据代入式(4)估计自相关矩阵R；

2) 对自相关矩阵进行奇异值分解，通过提取大特征值对应的特征向量划分信号子空间${U_s}$，如式(8)所示；

3) 将信号子空间按照子阵划分形式，分块为${\mathit{\boldsymbol{U}}_{s1}}$、${\mathit{\boldsymbol{U}}_{s2}}$；

4) 对式(18)或式(20)进行特征值分解或者利用QZ分解求解广义特征值分解，得到广义特征向量矩阵估计T；

5) 利用式(13)及式(14)，估计子阵1的阵列流形${\hat {\mathit{\boldsymbol{A}}}_1}$；

6) 利用式(16)对波达方向进行估计。

基于广义特征向量的ESPRIT算法具有传统ESPRIT算法计算量小的优点，增加的计算量集中体现在式(14)及式(16)的${\hat \phi _k}$上，仅仅增加了$2(M - 1)K$次乘法运算及$(M - 3)K$次加法运算。不考虑ESPRIT算法与MUSIC算法都存在的计算过程，MUSIC算法谱峰搜索需要计算量为$O[{N_\theta }(M(M-K) + M)]$，${N_\theta }$为网格划分数。而以QZ算法为例，需要$O(66{K^3})$^[8]，加上增加的计算量仍然远远小于MUSIC算法(一般有${N_\theta } {\square} M > K$)。

4 实验仿真

为了讨论基于广义特征向量的ESPRIT算法的性能，在同一环境下分别采用基于谱峰搜索的MUSIC算法、基于广义特征值的ESPRIT算法以及基于广义特征向量的ESPRIT算法进行波达方向估计的3个实验。

仿真条件为8个阵元组成的均匀线阵，阵元间隔为1/2个波长，噪声均为高斯白噪声。为了讨论在不同条件下的性能，每个试验都重复进行了500次独立实验。

4.1 实验设置

实验一：信噪比的影响

实验设置3个非相关的远场窄带高斯信号源分别从-20.1°、1°及29.2°入射到阵列。每次实验快拍数均取20。为了讨论信噪比对算法估计的影响，信噪比取值范围为-5 dB到20 dB，每隔3 dB取一个值，针对每个信噪比设置进行实验。结果如图 1所示。

实验二：快拍数的影响

信号源个数以及来波方向设置如实验一。快拍数取值范围为10~100，每隔5取一个值。针对每个设置的快拍数分别进行实验。信噪比均设置为10 dB。实验结果如图 2所示。

实验三：多目标估计效果

设置信噪比固定为10 dB，快拍数为20。为了考虑估计目标个数对算法性能的影响，分别进行了7次实验。第k($k = 1, 2, \cdots , 7$)次实验分别取[-50.35°, -35.75°, -19.6°, -4.2°, 10.41°, 24.3°, 43.8°]中的前k个值为入射信号的角度，入射信号均为非相关远场信号，实验结果如图 3所示。

4.2 实验结果分析

从3个实验看出，基于广义特征值分解的ESPRIT算法不管在任何情况下估计性能均优于基于特征值的ESPRIT算法，特别在信噪比小或快拍数小的情况下，基于广义特征向量的ESPRIT算法与MUSIC算法的均方根误差曲线接近重合。

为了度量基于广义特征向量的ESPRIT算法对基于广义特征值的ESPRIT算法的改进程度，引入改进率概念：

$ {\rm{imp}} = \frac{1}{{{\rm{Num}}}}\sum\limits_i {\frac{{({\rm{ev}}{{\rm{a}}_i} - {\rm{ev}}{{\rm{e}}_i})}}{{{\rm{ev}}{{\rm{a}}_i}}}} $

(21)

式中，${\rm{ev}}{{\rm{a}}_i}$为基于广义特征值的ESPRIT算法第i次试验下的误差；${\rm{ev}}{{\rm{a}}_i}$为基于广义特征向量的ESPRIT算法第i次试验下的误差；Num为设置的个数。

综合3次试验的结果，在不同信噪比、快拍数下，基于广义特征向量的ESPRIT的平均改进率达到了18%。在不同信源个数下的改进率分别为[15% 12% 10% 8% 5% 2% 0]。

5 结束语

旋转不变矩阵特征值含有波达方向信息，但是利用该信息估计波达方向误差较大。本文在ESPRIT算法框架下提出利用广义特征向量对波达方向进行估计。首先建立驱动矢量与广义特征向量的关系；然后利用阵列位置及信号波长的先验信息反演出波达方向，提高了ESPRIT算法的估计精度。