DOA估计算法

阵列信号处理中的DOA估计算法 摘要：本文简要介绍了阵列信号处理的基本知识和其数学模型，并且对阵列信号处理中很重要的来波方向（DOA）估计方法进行了比较，主要包括古典谱估计方法、Capon最小方差法、多重信号分类（MUSIC）算法以及旋转不变因子空间（ESPRIT）算法。通过这些算法的介绍和比较，我们可以很方便地在不同的情况下选择不同的算法去对信号的来波方向进行估计。 关键词：阵列信号处理；来波方向（DOA）；MUSIC；自相关矩阵；特征分解；ESPRIT DOA Estimation Algorithms in Array Signal Processing Abstract：In this paper, we have introduced the basic knowledge and data model of array signal processing and have compared many DOA estimation methods in array signal processing，which included classical spectrum estimation method、Capon minimum variance method、MUSIC method and ESPRIT method。Through the introduction and comparison of these algorithms，we can choose different algorithm to estimate the DOA of signal in different situation，conveniently。 Key words：array signal processing；DOA；MUSIC；self-correction matrix；eigendecomposition；ESPRIT 1.引言 近几十年来，阵列信号处理作为信号处理的一个重要分支，在声纳、雷达、通信以及医学诊断等领域得到了相当广泛的应用和发展。阵列信号处理是指在一定大小空间的不同位置去设置传感器，组成传感器阵列，利用传感器阵列去接收空间中的信号并且通过一定的方法对接收的信号进行处理。阵列信号处理的目的是为了增强有用的信号，抑制无用的干扰和噪声，并且从接收的信号中提取出有用信号的特征以及信号所包含的信息。与传统的单个定向传感相比，传感器阵列具有比较高的信号增益、灵活的波束控制、很高的空间分辨率以及极强的干扰抑制能力。阵列信号处理研究的主要问题包括[5]：空间谱估计——对空间信号波达方向进行超分辨估计；零点形成技术——使天线的零点对准干扰方向；波束形成技术——使阵列方向图的主瓣指向所需的方向。其研究的三个主要方向分别在不同的时期进行了不同的主要研究，这三个阶段分别是：d 1、20世纪60年代主要集中在波束形成技术方面[1]，如自适应相控天线、自适应波束操控天线和自适应聚束天线等，主要目的是使阵列方向图的主瓣指向所需要的方向。 2、20世纪70年代主要集中在零点形成技术方面[2]，如自适应置零技术、自适应调零技术、自适应杂波抑制和自适应旁瓣相消等，可以提高信号输出的信噪比（SNR）。 3、20世纪80年代主要集中在空间谱估计方面[3]，如最大似然谱估计、最大熵谱估计、子空间谱估计等，它是现代谱估计理论与自适应阵列技术结合的产物，主要是研究在阵列处理带宽内空间信号的波达方向的估计问题，这标志着阵列信号处理研究的重大变化。 信号的波达方向（DOA）估计是阵列信号处理领域的一个非常重要的研究内容。信号的DOA估计算法大多是一种极值搜索法，即首先形成一个包含待估计参数的函数（一般是一个伪谱函数），然后通过对该函数进行峰值搜索，得到的极值就是信号的波达方向。这些算法主要包括：1965年Bartlett基于波束形成的思想提出的DOA估计算法，但是该算法不能分辨出两个空间距离小于波束宽度的信号源。1968年Schweppe首先研究了虽大似然估计算法（ML），但是比较重要的还是后来Capon提出的高进度的ML，该算法对于服从高斯分布的信源估计可以达到克劳—拉美界，但是需要对接收阵列数据的自相关矩阵进行求了逆运算，运算量相当大。1979年Schmidt提出了多重信号分类法[4]（Multiple Signal Classification，MUSIC）以及各种改进的MUSIC算法等，它们都需要进行特征值分解运算，可以得到比较高精度的参数估计，但是计算量太大。1985年Roy和Kailath提出了一种借助旋转不变技术的参数估计算法[6]（Estimating Signal Via Rotational Invariance Techniques，ESPRIT），它是利用阵列流行的某些特性形成一个可以直接求解的函数，能够比较方便的得到所需要的估计参数。在此之后，人们以MUSIC和ESPRIT为基础，提出了各种各样的算法，例如最小范数法[7]、ROOT-MUSIC[8]、TLS-ESPRIT[9]等。这些不同的算法是基于不同的理论提出的，并且建立在不同的约束条件之下，所以其特性和适用对象也会不同。 2.数据模型 2.1 平面波与阵列 在无线通信中我们通过天线对电磁波进行发射和接收。为了增加电磁波的利用率和电磁波的波束形状可控，一般采用阵列天线。在一般情况下，将一组传感器按一定的方式设置在空间不同的位置上组成传感器阵列，此传感器阵列能够接收空间的传播信号，然后对所接收到的信号经过适当的处理并提取所需的信号源和信号属性等信息，包括信号辐射源辐射信号的数目、方向、幅度等。一般来说，构成阵列的阵元可以按照任意的方式进行排列，但是通常是按照直线等距、圆周等距或平面等距排列的，并且取向相同。为了简化天线阵列的分析，通常作如下假设[10]： 1. 窄带假设：这样可以保证所有阵元几乎同时接收到该信号，即阵元接收之间的信号包络没有变化； 2. 信号的统计特性：假设入射到阵列的信号为平稳且各态历经，这样可以用时间平均来代替统计平均。噪声为互不相关的白噪声，方差为 。 3. 忽略阵元之间的互耦； 4. 信号的数目要小于阵元的数目，并且阵列接收到得所有信号的波达方向互不相同，信号之间互不相关； 5. 平面波假设：假设信源到阵列的距离远大于阵列的口径，从而所有入射到阵列的信号波前金额以近似为平面波。 假设在天线阵的原唱存在 个信号源，则所有到达阵列的波前可近似为平面波。若天线阵由 个全向天线组成，将第一个阵元设为参考阵元，则到达参考阵元的第 个信号为： （1） 式中， 为第 个信号的复包络，包含信号信息。 为空间信号的载波。由于信号满足窄带假设条件，则 ,那么经过传播延迟 后的信号可以示为： （2） 则理想情况下第 个阵元接收到的信号可以表示为: （3） 式中， 为第 个阵元到达第 个阵元时相对于参考阵元的时延， 为第 阵元上的加性噪声。根据式（2）和（3）可得，整个天线阵接收到得信号为： （4） 式中， 为信号 的方向向量， 为阵列流形， 为信号矩阵， 为加性噪声矩阵， 表示矩阵转置。 2.2 均匀线阵与均匀圆阵 在实际中一般使用均匀线阵和均匀圆阵等阵列结构。 （1）均匀线阵 均匀线阵（ULA：Uniform Linear Array）是一最简单常用的阵列形式，如图1所示，将 个阵元等距离排列成一直线，阵元间距为 。假定一信源位于远场，即其信号到达各阵元的波前为平面波，其波达方向（DOA）定义为与阵列法线的夹角 。 图1 ULA示意图 以第一个阵元为参考阵元，则各阵元相对参考阵元的时延为： （5） 由此可得等距线阵的方向向量为： （6） 当波长和阵列的几何结构确定时，该方向向量只与空间角 有关，因此等距线阵的方向向量记为 ，它与基准点的位置无关。若有 个信号源，其波达方向分别为 ， ，则阵列流形矩阵为： （7） 以上给出了等距线阵的方向向量的表示形式。实际使用的阵列结构方向向量 与空间角 一一对应，不能出现模糊现象。这里需要说明的是：阵元间距 是不能任意选定的，甚至有时需要非常精确的校准。假设 很大，相邻阵元的相位延迟就会超过 ，此时，阵列方向向量无法在数值上分辨出具体的相位延迟，就会出现相位模糊。可见，对于等距线阵来说，为了避免方向向量和空间叫之间的模糊，其阵元间距不能大于半波长 ，以保证阵列流形矩阵的各个列向量线性独立。天线阵列的输出为： （8） 其向量形式为： （9） 式中， 为权重向量。 （2）均匀圆阵 均匀圆周阵列简称均匀圆阵（UCA: Uniform Circular Array），是平面阵列，它的有效估计是二维的，能够同时确定信号的方位角和仰角。均匀圆阵由 个相同各向同性阵元均匀分布在 平面一个半径为 的圆周上，如图2所示。采用球面坐标系表示入射平面波的波达方向，坐标系的原点 位于阵列的中心，即圆心。信源俯角 是原点到信源的连线与 轴的夹角，方向角 则是原点到信源的连线在 平面上的投影与 轴之间的夹角。 阵列的第 个阵元与 轴之间的夹角为 ，则该处的位置向量为[5]： （10） 在某个时刻，原点和第 个阵元接收到得信号的复包络间的相位差为： （11） 式中， ， 。 均匀圆阵相对于波达方向为 的信号的方向向量为： （12） 图2 UCA示意图 3.DOA估计算法的特性比较 3.1 古典谱估计法 古典谱估计法是通过计算空间谱求取其局部最大值，从而估计出信号的波达方向。Bartlett波束形成方法是经典傅里叶分析对传感器阵列数据的一种自然推广。Bartlett方法使波束形成器的输出功率相对于某个输入信号最大。设希望来自 方向的输出功率为最大，则代价函数为： （13） 在白噪声方差 一定的情况下，权重向量的范数 不影响输出信噪比，故取权重向量的范数为1，用拉格朗日因子的方法求得上述最大优化问题的解为： （14） 从式（14）可以看出，阵列权重向量是使信号在各阵元上产生的延迟均衡，以便使它们各自的贡献最大限度地综合在一起。空间谱是以空间角为自变量分析到达波的空间分布，其定义为： （15） 将所有方向向量的集合 成为阵列流形。在实际应用中，阵列流形可以在阵列校准是确定或者利用接收的采样值计算得到。 从式（15）可知，利用空间谱的峰值就可以估计出信号的波达方向。当有 个信号存在时，对于不同的 ，利用式（15）计算得到不同的输出功率。最大输出功率对应的空间谱的峰值也就最大，而最大空间谱峰值对应的到的DOA值即为信号波达方向的估计值。古典谱估计方法将阵列所有可利用的自由度都用于在所需观测方向上形成一个波束。当只有一个信号时，这个方法是可行的。但是当存在来自多个方向的信号时，阵列的输出将包括期望信号和干扰信号，估计性能会急剧下降。而且该方法要受到波束宽度和旁瓣高度的限制，这是由于大角度范围的信号会影响观测方向的平均功率，因此，这种方法的空间分辨率比较低。我们可以通过增加天线阵列的阵元来提高分辨率，但是这样会增加系统的复杂度和算法对于空间的存储要求。 3.2 Capon最小方差法 为了解决Bartlett方法的一些局限性，Capon提出了最小方差法。该方法使部分(不是全部)自由度在期望观测方向形成一个波束，同时利用剩余的自由度在干扰信号方向形成零陷，可以使得输出功率最小，达到使非期望干扰的贡献最小的目的，同时增益在观测方向保持为常数，通常为1，如式 所示： （16） 其中 是接收信号 的协方差矩阵。 求解式（16）得到的权向量通常称为最小方差无畸变响应(MVDR，Mhmnum Variance Distortionless Response)波束形成器权值，因为对于某个观测方向，它使输出信号的方差(平均功率)最小，又能使来自观测方向的信号无畸变地通过(增益为1，相移为0)。这是个约束优化问题，可以利用拉格朗日乘子法求解。 令 ， 分别对 和 求偏导数可得： （17） 式（17）两端分别左乘 得： 上式两端分别右乘 得： 因此， （18） 对（18）式两端分别右乘 有： 所以， （19） 将式（19）带入（18）中，并对两边取共轭对称，最终得到： （20） 利用Capon波束形成法得到的空间功率谱公式如下： （21） 计算Capon谱并在全部 范围上搜索其峰值，就可估计出DOA。 虽然与古典谱估计法相比，Capon法能提供更佳的分辨率，但Capon法也有很多缺点。如果存在与感兴趣信号相关的其他信号，Capon法就不能再起作用，因为它在减小处理器输出功率时无意中利用了这种相关性，而没有为其形成零陷。换句话说，在使输出功率达到最小的过程中，相关分量可能会恶性合并。另外，Capon法需要对矩阵求逆运算，这样会使得计算量非常大。 3.3 MUSIC算法 Music算法是由R．O．Schmidt于1979年提出来，1986年重新发表的。它是最早的也是最经典的超分辨DOA估计方法，它利用了信号子空间和噪声子空间的正交性，构造空间谱函数，通过谱峰搜索，检测信号的DOA。它是建立在以下假设基础上的： (1)阵列形式为线性均匀阵，阵元间距不大于处理最高频率信号波长的二分之一； (2)处理器的噪声为加性高斯分布，不同阵元间距噪声均为平稳随机过程，且相互独立，空间平稳(各阵元噪声方差相等)； (3)空间信号为零均值平稳随机过程，它与阵元噪声相互独立，且信号间相互独立； (4)信号源少于阵元数，信号取样数大于阵元数。 在此假设基础上，Music算法对DOA的估计从理论上可以有任意高的分辨率。 Music算法原理如下： 由式（4）可得接收信号的协方差矩阵为： （22） 由于假设信号与噪声是不相关的，且噪声为平稳的加性高斯白噪声，因此式（22）中的二，三项为零，且有 。则式（22）简化为式（23）： （23） 式（23）中的 是有用信号的协方差矩阵。 由于假设信号源之间互不相关，因此 为满秩矩阵，其秩为 。而 为 维的矩阵，其秩也是 ，并且 是Hermite半正定矩阵，其秩也是 。因此，令 的特征值为 ，那么 的 个特征值为： 它们对应的特征向量分别为 ，其中前 个对应大特征值，后 个对应小特征值。 由此可以看出，协方差矩阵 经过特征值分解后可以产生 个较大的特征值和 个较小的特征值，并且这 个小特征值非常接近。所以当这些小特征值的重数 确定了，那么信号的个数就可以由式（24）估计出来： （24） 对于与 个最小特征值对应的特征向量，有： ， 即： ， 因为 满秩， 非奇异，因此： 或 这表明与 个最小特征值对应的特征向量，和 个信号特征值对应的方向向量正交，即信号子空间和噪声子空间正交。因此，我们构造 维的噪声子空间： 并定义Music空间谱为： （25） 或 （26） 由于信号子空间和噪声子空间正交，所以当 等于信号的入射角时，Music空间谱将产生极大值。因此当对Music空间谱搜索时，其 个峰值将对应 个信号的入射方向，这就是Music算法。 现将Music算法的步骤归纳如下： （1）收集信号样本 ， ，其中 为采样点数，估计协方差函数： （2）对 进行特征值分解： 式中 为特征值对角阵，且从大到小顺序排列 是对应的特征向量。 （3）利用最小特征值的重数 ，按照式（24）估计信号数 ，并构造噪声子空间 。 （4）按照式（25）搜索Music空间谱，找出 个峰值，得到DOA估计值。 尽管从理论上讲，Music算法可以达到任意精度分辨，但是也有其局限性。它在低信噪比的情况下不能分辨出较近的DOA，另外，当阵列流行存在误差时，对Music算法也有较大的影响。 3.4 Music算法的改进 人们对于Music算法提出了各种改进，以提高分辨率，减小计算的复杂度。其中一种改进方法是Barabell提出的求根-MUSIC算法[8]。这种方法根据多项式求根，可以提供更高的分辨能力，但是只适用于均匀线阵。Barabell提出的另一种是利用信号空间特征向量（主特征向量）的性质，定义了具有更加分辨率的有理谱函数。Schell在1989年也提出了利用信号的谱相干特性改善常规MUSIC算法性能的循环MUSIC算法[11]。下面分别介绍。 1. 求根-MUSIC算法 对于阵元间距为 的均匀线阵，方向向量 的第 个元素可以表示为： （27） 式（26）给出的MUSIC谱是一个全极点函数，即： （28） 式中， 。利用式（27），式（28）的分母可以写作： （29） 式中， 是 中第 行、第 列上的元素。将两个累加合并在一起，式（29）可以简化为： （30） 式中， 是 中第 条对角线上的元素之和。 定义如下多项式，即： （31） 则评价MUSIC谱 等价于评价单位圆上的多项式 ，因为 的根靠近单位圆，则MUSIC谱的峰值存在。在没有噪声的理想情况下，极点恰好落在单位圆上，位置由波达方向确定。换句话说， 在 处的一个极点，即在MUSIC谱产生峰值的位置上，故有： （32） 求根-MUSIC算法比MUSIC谱形式的算法具有更好的分辨率，而且在低信噪比SNR的情况下也能够很好的工作。该算法避免了传统MUSIC算法的谱峰搜索过程，大大节省了计算量，但是该算法只适用于等距线阵。 2. 循环MUSIC算法 循环MUSIC算法是一种利用接收信号的谱相干性和空间相干性的信号选择性的定向算法。将谱相关性和MUSIC结合起来，在相距很近的信号中只有一个感兴趣信号（Signal of Interest，SOI）且信号间隔小于阵列阙值时，能够分辨出该期望信号。循环MUSIC算法还不受入射到阵列上的信号数（包括SOI和干扰）必须小于阵元数这一要求的约束。 考察一个 元阵列，接收的 个信号在 频率上具有谱相关性，而干扰信号（数目任意）在该频率上不具有谱相关性。一个例子是，在严重的同信道干扰环境中检测一个具有特定谱相关和多径分量数的期望信号。令 ， 为期望信号， 为入射到阵列上的噪声和干扰分量。于是接收信号向量 可以写作： （33） 因为只有期望信号在 具有谱相关性，接收信号 的循环自相关矩阵 定义为： （34） 可以表示为： （35） 式中， 是期望信号的循环自相关矩阵，定义为： （36） 式中 （37） 显然，矩阵 的秩为 。对于 ， 的零空间由对应于零特征值的特征向量 张成为： （38） 如果信号不完全相关，则 满秩，等于 。因为 也满秩，由式（35）和式（38）知， 的零空间正交于期望信号的方向向量，即： ， （39） 将式（39）作为正交性的量度，式（25）的循环MUSIC谱可以定义如下： （40） 对于所有的 ，搜索出 的 个谱峰，就可以算出期望信号的波达方向。 3.5 ESPRIT算法 由于MUSIC算法需要进行谱峰搜索，计算量很大，因此在实际的应用中对于系统的计算速度要求较高。在MUSIC算法以后，人们开始研究各种不需要进行谱峰搜索的快速DOA算法。有Roy等人提出的旋转不变子空间（ESPRIT）算法是空间谱估计中的另一种经典算法[6]。ESPRIT算法的基本思想是利用旋转不变因子技术来估计信号参数，它把传感器阵列分解成两个完全相同的子阵列，在两个子阵列中每两个相对应的阵元有着相同的位移，即阵列具有平移不变性，每两个位移相同的传感器配对。在实际情况下，比如等间距的直线阵列或双直线阵列都可以满足ESPRIT算法对于阵列天线的要求。它同MUSIC算法一样，也需要对阵列接收数据自相关矩阵进行特征值分解，但是两者存在明显的不同，MUSIC算法利用了自相关矩阵信号子空间的正交性，而ESPRIT算法利用了自相关矩阵信号子空间的旋转不变特性。ESPRIT算法不需要知道阵列的几何结构，因此对于阵列的校准要求比较低，现在ESPRIT算法已经成为主要的DOA估计算法之一。 设由 个对偶极子组成的阵元数为 的天线阵列，两个子阵列对应元素具有相等的敏感度模式和相同的位移偏移量 ， 个独立的中心频率为 的窄带信号源入射到该阵列，两个子阵列第 组对应阵元的接收信号可以表示为： （41） （42） 式中， 表示第 个信号源的入射方向，将每个子阵列的接收信号表示成向量形式有： (43) （44） 式中， 是带噪声的数据向量， 表示两个阵列之间的相位延迟，也称为旋转不变因子， 和 为加性噪声向量。 定义整个阵列的接收向量为 ，用子阵列接收向量来表示： （45） 式中， （46） 天线阵列接收向量 的自相关矩阵为： （47） 设 ，则 的 个最小的广义特征值等于 ，而与 个最大广义特征值相对应的特征向量 满足： （48） 式中， 表示由矩阵中的向量张成的空间。 则存在唯一的非奇异矩阵 满足： （49） 利用阵列的旋转不变结构特性， 可以分解成为 和 。 （50） 由于 和 共享一个列空间， 的秩为 ，则： （51） 这表明存在一个唯一的秩为 的矩阵 可满足： （52） 定义： （53） 把式（53）带入式（52）可得： （54） 如果信号的入射方向不同，则阵列流行 是满秩的，则可以得到： （55） 显然， 的特征值必然等于对角矩阵 的对角元素，而 的列向量为 的特征向量。 ESPRIT算法避免了大多数DOA估计方法所固有的搜索过程，大大减小了计算量，并降低了对于硬件的存储要求。和MUSIC算法不同的是，ESPRIT算法不需要精确知道阵列的流行向量，因此，对阵列校正的要求不是很严格。 4.结论 本文对阵列信号的基本知识进行了简单的阐述和介绍。并且对阵列信号处理中很重要的DOA估计方法进行了介绍和探讨，比较了不同DOA估计方法的优缺点。特别是对于MUSIC算法，给出了两种改进的MUSIC算法，可以增加我们对于MISIC算法的理解。文中还指出了在不同的条件下这些方法的不同特性，为以后学习和研究DOA估计方法提供了比较有价值的参考。 参考文献 [1] Special issue on adaptive arrays, IEEE Trans. Antennas Propagat., 1964, vol. 12. [2] Special issue on adaptive arrays, IEEE Trans. Antennas Propagat., 1976, vol. 24. [3] Special issue on adaptive arrays, IEEE Trans. Antennas Propagat., 1986, vol. 34. [4] Schmidt R O. Multiple emitter location and signal parameter estimation. IEEE Trans. Antennas Propagat., 1986, 34（3）:276-280. [5] 张贤达，保铮. 通信信号处理. 北京：国防工业出版社，2000. [6] Roy R, Kailath T. Estimating of signal parameters via rotational invariance techniques. IEEE Trans. on Acoustics Speech, Signal Processing, 1989, 37（7）：984-995. [7] Kumaresan R, Tufts D W. Estimating the angle-of-arrival of multiple plane waves. IEEE Trans. Aerospace Electronic System, 1983, 19（1）：135-139. [8] Barael A J. Improving the resolution of eigen-structure based direction finding algorithms. ICASSP 1983, Boston, 1983, 336-339. [9] Zhang X D, Liang Y C. Prefiltering-based ESPRIT for estimating parameters of sinusoids in non-Gaussian ARMA noise. IEEE Trans. Signal Processing, 1995, 43（1）：349-353. [10] 王永良，陈辉，彭应宁，万群. 空间谱估计理论与算法. 北京：清华大学出版社, 2004. [11] Schell S V, Calabretta, Gardner W A, etal. Cyclic MUSCI Algorithms for Signal Selective DOA Estimation, Proc. of the IEEE Int’l. Conf. on Acoustics, Speech and Signal Processing, 1989:2278-2281.