esprit算法研究

精品文档课程设计报告实验名称：ESPRIT算法研究实验日期：姓 名：学 号： 哈尔滨工业大学（威海）一、 设计任务实现空间谱估计算法，并考察算法性能。二、 方案设计1） 由均匀线阵形式，确定阵列的导向矢量；2） 由阵列导向矢量，对接收信号进行建模仿真；3） 由ESPRIT算法实现信号DOA估计；4） 考察算法性能与信噪比，采样率，观测时间等参数的关系。三、 设计原理3.1空间谱估计数学模型空间谱估计就是利用空间阵列实现空间信号的参数估计的一项专门技术。整个空间谱估计系统应该由三部分组成：空间信号入射、空间阵列接收及参数估计。相应地可分为三个空间，即目标空间、观察空间及估计空间，也就是说空间谱估计系统由这三个空间组成，其框图见图1。图1 空间谱估计的系统结构对于上述的系统结构，作以下几点说明。(1)目标空间是一个由信号源的参数与复杂环境参数张成的空间。对于空间谱估计系统，就是利用特定的一些方法从这个复杂的目标空间中估计出信号的未知参数。(2)观察空间是利用空间按一定方式排列的阵元，来接收目标空间的辐射信号。由于环境的复杂性，所以接收数据中包括信号特征(方位、距离、极化等)和空间环境特征(噪声、杂波、干扰等)。另外由于空间阵元的影响，接收数据中同样也含有空间阵列的某些特征(互耦、通道不一致、频带不一致等)。这里的观察空间是一个多维空间，即系统的接收数据是由多个通道组成，而传统的时域处理方法通常只有一个通道。特别需要指出的是：通道与阵元并不是一一对应，通道是由空间的一个、几个或所有阵元合成的(可用加权或不加权)，当然空间某个特定的阵元可包含在不同的通道内。(3)估计空间是利用空间谱估计技术(包括阵列信号处理中的一些技术，如阵列校正、空域滤波等技术)从复杂的观察数据中提取信号的特征参数。从系统框图中可以清晰的看出，估计空间相当于是对目标空间的一个重构过程，这个重构的精度由众多因素决定，如环境的复杂性、空间阵元间的互耦、通道不一致、频带不一致等。3.2 阵列信号处理首先，考虑N个远场的窄带信号入射到空间某阵列上，阵列天线由M个阵元组成，这里假设阵元数等于通道数，即各阵元接收到信号后经过各自的传输信道送到处理器，也就是说处理器接收来自M个通道的数据。 (3.2-1)式中，是接受信号的幅度，是接收信号的相位，是接收信号的频率。在窄带远场信号源的假设下，有 (3.2-2)根据式(3.2-1)和式(3.2-2)，显然有下式成立： (3.2-3)则可以得到第L个阵元接收信号为 (3.2-4)式中，为第L个阵元对第i个信号的增益，表示第L个阵元在t时刻的噪声，表示第i个信号到达第L个阵元时相对参考阵元的时延。将M个阵元在特定时刻接收的信号排列成一个列矢量，可得(3.2-5)在理想情况下，假设阵列中各阵元是各向同性的且不存在通道不一致、互耦等因素的影响，则式(3.2-4)中的增益可以省略(即归一化1)，在此假设下式(3.2-5)可以简化为 (3.2-6)将式(3.2-6)写成矢量形式如下： (3.2-7)式中，为阵列的维快拍数据矢量，为阵列的维噪声数据矢量，为空间信号的维矢量，A为空间阵列的维流型矩阵(导向矢量阵)，且 (3.2-8)其中导向矢量 (3.2-9)式中，c为光速，为波长。由上述的知识可知，一旦知道阵元间的延迟表达式，就很容易得出待定空间阵列的导向矢量或阵列流型。下面推导一下空间阵元间的延迟表达式。假设空间任意两个阵元，其中一个为参考阵元(位于原点)，另一个阵元的坐标为(x ,y, z)，两阵元的几何关系见图，图中“”表示阵元。图2 空间任意两阵元的几何关系由几何关系可以推导出两阵元的波程差为 (3.2-10)这里的波程差其实就是位于x轴上两阵元间的延迟、位于y轴上两阵元间的延迟和位于z轴上两阵元间的延迟之和。根据式(3.2-10)的结论，下面给出实际环境中常用的几种阵列及阵元间的相互延迟表达式。(1)平面阵设阵元的位置为，以原点为参考点，另假设信号入射参数为，分别表示方位角与俯仰角，其中方位角表示与x轴的夹角。(2)线阵设 阵元的位置为，以原点为参考点，另假设信号入射参数为，表示方位角，其中方位角表示与y轴的夹角（即与线阵法线的夹角），则有 (3.2-11)(3)均匀圆阵 设以均匀圆阵的圆心为参考点，则有 (3.2-12)其中方位角表示与x轴的夹角，r为圆半径。3.3旋转不变子空间算法原理3.3.1信号模型 算法介绍前，首先对信号进行建模。为了推导分析的方便，将波达方向的数学模型做如下理想状态的假设：1) 阵列形式为线性均匀阵，阵元间距不大于信号波长的二分之一。2) 存生两个完全相同的子阵，且两个子阵的间距是己知的。 3) 噪声序列为一零均值高斯过程，各阵元间噪声相互独立，噪声与信号也相互独立。4) 空间信号为零均值平稳随机过程，通常为窄带远场信号。5) 信号源数小于子阵阵列元数，信号取样数大于子阵阵列元数，以确保子阵阵列流型的各列线性独立。6) 组成阵列的各传感器为各向同性阵元，且无互耦以及通道不一致的干扰。图3-1均匀线阵的数学模型示意图下图给出了均匀线阵的数学模型示意图:3.3.2 算法原理对于均匀线阵，相邻子阵间存在一个固定间距，这个固定间距反映出各相邻子阵间的一个固定关系，即子阵间的旋转不变性，而ESPRIT算法正是利用了这个子阵间的旋转不变性实现阵列的DOA估计。ESPRIT算法最基本的假设是存在两个完全相同的子阵，且两个子阵的间距是已知的。由于两个子阵的结构完全相同，且子阵的阵元数为m，对于同一个信号而言，两个子阵的输出只有一个相位差，=1，2， N。下面假设第一个子阵的接收数据为，第二个子阵的接收数据为，根据前面所述的阵列模型可知 (3.1) (3.2)式中，子阵1的阵列流型=A，子阵2的阵列流型= ，且式中 （3.3）从上面的数学模型可知，需要求解的是信号的方向，而信号的方向信息包含在A和中，由于是一个对角阵，所以下面只考虑这个矩阵，即 （3.4）由上可知。只要得到两个子阵间的旋转不变关系，就可以方便地得到关于信号到达角的信息。下面的任务就是从式（3.1）和式（3.2）中得到两个子阵间的关系。先将两个子阵的模型进行合并，即 （3.5）在理想条件下，可得上式的协方差矩阵 (3.6) 对上式进行特征分解可得 （3.7）显然上式中得到的特征值有如下关=，US为大特征值对应的特征矢量张成的信号子空间，为小特征值对应矢里张成的噪声子空间。对于实际的快拍数据，式(3.7)应修正如下： (3.8) 由前面的知识可知，上述的特征分解中大特征矢量张成的信号子空间与阵列流型张成的信号子空间是相等的。即 (3.9)此时，存在一个惟一的非奇异矩阵T，使得 （3.10）显然，上述的结构对两个子阵都成立，所以有 （3.11） 很显然 ，由子阵1的大特征矢量张成的子空间、由子阵2的大特征矢量张成的子空间与阵列流型A张成的子空间三者相等，即 （3.12） 另外，由两个子阵列在阵列流型上的关系可知 （3.13）再利用式（3.11）可知两个子阵列的信号子空间的关系如下: （3.14）式（3.13）反映了两个子阵列的阵列流型间的旋转不变性，而式（3.14）反映了两个子阵的阵列接收数据的信号子空间的旋转不变性。如果阵列流型A是满秩矩阵，则由式（3.14）可以得到 （3.15）所以上式中的特征值组成的对角阵一定等于，而矩阵T的各列就是矩阵特征矢量。所以一旦得到上述的旋转不变关系矩阵，就可以直接利用式（3.4）得到信号的入射角度。3.4 标准的旋转不变子空间算法有上节的知识可知， ESPRIT算法的基本原理就是利用式（3.14）的旋转不变性，常规的旋转不变子空间算法就是利用上述的基本原理求解信号的入射角度信息。下面就分析解这个等式的两种最经典、应用最广泛方法:最小二乘（LS）法和总体最小二乘（TLS）法。3.4.1 最小二乘法由最小二乘的数学知识，我们知道式（3.14）的最小二乘解的方法等价于，约束条件 （3.16）因此最小二乘法的基本思想就是使校正项尽可能小，而同时保证满足约束条件。为了得到LS解，将式（3.14）代入式（3.16）即得 （3.17）对上式进行展开可得 = （3.18）上式对求导并令其等于0，可得 （3.19）上式的解显然有两种可能:(1) 当满秩时，也就是子阵1的信号子空间的维数等于信号源数时，则上式的解是唯一的，可得上式的最小二乘解 （3.20） (2) 当不满秩，即时，也就是信号源间存在相干或相差时，则存在很多解，但我们却无法区别对应于方程的各个不同的解，可以称这些解是不可辨识的，解的不可辨识性是我们需要解相干的原因所在。下面给出LS-ESPRIT算法的求解步骤：1由两个子阵的接收数据,，分别得到两个子阵的数据协方差矩阵；2对矩阵对R, 进行特征分解，从而得到两个数据矩阵的信号子空间和；3按式（3.20）得到矩阵，然后对其进行特征分解.得到N个特征值，就可得到对应的N个信号的到达角。当考虑嗓声影响时，上述基于最小二乘算法的估计都是有偏的，这就是为什么需要考虑总体最小二乘ESPRIT算法的原因。3.4.2 总体最小二乘法 我们知道，普通最小二乘的基本思想是用一个范数平方为最小的扰动去于扰信号子空间，目的是校正中存在的嗓声。显然这就存在一个问题：如果同时扰动和，并使扰动范数的平方保持最小，是否可以同时校正和中存在的嗓声?答案是肯定的，这就是总休最小二乘（TLS）的思想。它考虑的是如下矩阵方程的解: (3.27)显然上式可以改写成 （3.28）所以TLS的解等价于 （3.29）定义如下一个矩阵，再结合上述分析过程。我们发现其实就是寻找一个的酉矩阵F，便得矩阵F与正交，也就说明了由F张成的空间与或列矢量张成的空间正交。所以矩阵F可从的特征分解中得到。因为 （3.30）式中的是由特征值构成的对角矩阵，E是与其相应的特征矢量构成的矩阵。即 （3.31）令是由对应特征值为0的特征矢量构成的矩阵.它属子噪声子空间，所以只要选择矩阵F使之等于、，即可满足上面提到的要求。即有 (3.32)可得 （3.33）如果令，则 （3.34）上式说明的特征值即是对角线元素。这说明通过构造一个矩阵就可得到有关信号角度的信息.而这个矩阵的构造可通过式（3.30）得到，即 （3.35）下面直接给出TLS-ESPRIT算法的求解步骤：1由两个子阵的接收数据，, 由式(3.8)得到数据协方差矩阵R；2通过矩阵对于的广义特征分解，得到维数为的信号子空间；3由构造矩阵，并按式防（3.30）进行特征分解得到矩阵E，然后再按式(3.31)将矩阵分为四个小的矩阵；4按式（3.35）得到矩阵，然后对其进行特征分解,得到N个特征值，就可得到对应的N个信号的到达角。通过分析，我们可以得到标准ESPRIT算法的计算过程如下：(1)通过特征值或奇异值分解（EVD或SVD）分别估计两个存在旋转不变关系的子阵的信号子空；(2)用上述的LS、TLS等方法求解式（3.14）所示的不变等式；(3)计算的特征值，其中如式（3.3）所示。然后利用式(3.4)求解人射信号的角度信息。就ESPRIT算法而言，TLS算法与LS算法性能基本一致，只是在低信噪比情况下TLS算法性能略好。四、 仿真结果 主要分析各个参数对估计误差的影响，误差函数定义如式(1)：4.1 信噪比 SNR对估计误差的影响分析 首先对信噪比 SNR离散化取值，然后求得不同信噪比下的误差，从而绘制出误差随信噪比改变的函数曲线如图 2 所示，图 2 中信噪比 SNR从- 15 取到 15，间隔为 1，运行次数为 100 次，其余条件如题中所述。由图 2 可知，随着信噪比的增大，估计误差会越来越小，即估计精度会越来越高。当待估计的信号方位角相差比较小时，估计的误差也会相应的增大。另外，若两信号为相干信号，则此方法将不能对其进行正确的估计。 4.2 阵元数 L对估计误差的影响分析与 4.1节类似，首先对阵元数 L离散化取值，然后求得不同阵元数下的误差， 从而绘制出误差随阵元数改变的函数曲线如图 3 所示，图 3 中阵元数从 K+1 取到K+25，间隔为 1，运行次数为 100 次，其余条件如题中所述。由于阵元数 L 需大于信号个数 K才能正确估计，故取值中含有信号个数 K。由图 3 可知，随着阵元数的增加，估计误差会越来越小，即估计精度会越来越高，但当阵元数大到一定程度后，对估计精度的影响则会慢慢的减小。 4.3 采样点数 N对估计误差的影响分析与 4.1节类似，首先对采样点数 N离散化取值，然后求得不同采样点数下的误差，从而绘制出误差随采样点数改变的函数曲线如图 4 所示，图 4 中采样点数从 10 取到 200，间隔为 5，运行次数为 100 次，其余条件如题中所述。由图 4可知，随着采样点数的增加，估计误差会越来越小，即估计精度会越来越高。 估计误差（角度）4.4 两信号之间的角度差（GAP）对估计误差的影响分析由于采用 ESPRI T算法对 DOA进行估计，若两信号的方位距离较近时，虽然能得出估计结果，但估计的精度会大受影响。因此，为了分析两信号之间的不同间隔会对估计精度造成多大的影响，绘制不同 GAP下的估计误差曲线如图 5所示。处理方法与 4. 1 节类似，图 5 中 GAP（单位为度）从 0. 1 取到 5，间隔为 0. 1，独立运行次数为 100 次，其余条件如题中所述。由图 5 可知，GAP越大估计越准确，但当 GAP大到一定程度后则估计精度趋于稳定。 4.5 单信号 DOA不同分布对估计误差的影响分析信号波达方向（DOA）的取值区间为-90度到 90度，若只考虑只有一个信号的情况，则当信号的 DOA不同时，估计误差也会不一样。因此，为了分析不同的 DOA会对估计精度造成多大的影响，绘制不同 DOA下的估计误差曲线如图 6所示。处理方法与 4. 1 节类似，图 6 中 GAP从- 80 度取到 80 度，间隔为 5 度，独立运行次数为 100 次，其余条件如题中所述。由图 6 可知，DOA越靠近 0 度估计越准确，越靠近正负 90 度估计误差越大。且仿真结果表明，当 DOA在正负 90 附近时，估计误差太大，因此，为了不影响估计结果显示效果，故在图中未绘制正负 90 度附近的估计误差。 2.6 减与不减噪声方差(Rn)对估计误差的影响分析由于有噪声的影响，因此在估计信号自相关矩阵 R时，若将无信号时的自相关矩阵 Rn减去，即相当与减去估计出噪声方差，则估计的精度会有所提高。结合信噪比 SNR对估计误差的影响，绘制减与不减噪声方差两种情况下估计误差随 SNR的变化曲线如图 7所示，图 7 中 SNR从-15dB到 5dB，间隔为 1dB，独立运行次数为 100次。仿真结果表明，若减 Rn，主要是在低信噪比时对估计精度的改善较大，当信噪比较大时二者几乎一样。 五、 程序清单%本文件名为 drawTLSesprit.m %分析基于总体最小二乘的 ESPRIT算法（TLS-ESPRIT）的 DOA估计的性能%clear;clc;close all; %清除变量，清屏，关闭所有绘图窗口% 调用格式:estimated,error=TLSesprit(p,L,K,SNR,DOA);% 估计结果(弧度，矢量：p行 1列)：estimated% 估计误差(弧度，标量：均方误差)：error% 信号个数：p% 阵元数：L% 快拍数：K% 信噪比：SNR% 波达方向(弧度，矢量：p行 1列)：DOA% p=2; L=8; K=100; SNR=5; DOA=pi*(-10/180) pi*(20/180);%显示估计结果%M=100; %设定独立重复运行次数DOA=pi*(0/180) pi*(30/180); %波达方向(弧度，矢量：p行 1列)p=length(DOA); L=8; K=100; SNR=5; %参数设置，estimated,error=TLSesprit(p,L,K,SNR,DOA); %函数调用polar(estimated,1 1,r*); %在极坐标中显示估计结果(必须先转化为弧度)h=title();set(h,string,TLS-ESPRIT: 估计值: ,num2str(estimated);h1=xlabel();set(h1,string,信号 DOA(度): ,num2str(DOA*180/pi);% %阵元数 L对估计误差的影响分析% Ln=p+1:1:p+25; %阵元数 L需大于信号个数 p才能正确估计% for n=1:length(Ln)% L=Ln(n);% for k=1:M% estimated,error=TLSesprit(p,L,K,SNR,DOA);% errorm(k)=error; %将每次的估计误差存入变量 errorm中，便于求均值% end% errorn(n)=sum(errorm)/M; %求多次运行后的估计误差的均值% end% figure(2);plot(Ln,errorn*180/pi,r:*,LineWidth,2); %绘制曲线，并适当标注% xlabel(阵元数 L);ylabel(估计误差( );title(阵元数 L对估计误差的影响);% % 结论：阵元数 L越大估计越准确，但当 L大到一定程度后则估计精度趋于稳定% %快拍数 K 对估计误差的影响分析% Kn=10:10:200; %对快拍数离散化取值% for n=1:length(Kn)% K=Kn(n);% for k=1:M% estimated,error=TLSesprit(p,L,K,SNR,DOA);% errorm(k)=error; %将每次的估计误差存入变量 errorm中，便于求均值% end% errorn1(n)=sum(errorm)/M; %求多次运行后的估计误差的均值% end% figure(3);plot(Kn,errorn1*180/pi,r:*,LineWidth,2); %绘制曲线并适当标注% xlabel(快拍数 K);ylabel(估计误差( );title(快拍数 K 对估计误差的影响);% % 结论：快拍数 K越大估计越准确，但当 K 大到一定程度后则估计精度趋于稳定% %信噪比 SNR对估计误差的影响分析% SNRn=-15:1:15; %对信噪比 SNR离散化取值% for n=1:length(SNRn)% SNR=SNRn(n);% for k=1:M % estimated,error=TLSesprit(p,L,K,SNR,DOA);% errorm(k)=error; %将每次的估计误差存入变量 errorm中，便于求均值% end% errorn(n)=sum(errorm)/M; %求多次运行后的估计误差的均值% end% figure(4);plot(SNRn,errorn*180/pi,r:*,LineWidth,2);%绘制曲线并适当标注(误差：角度)% xlabel(SNR);ylabel(估计误差( );title(SNR对估计误差的影响);% %结论：信噪比 SNR越大估计越准确，但当信噪比 SNR大到一定程度后则估计精度趋于 稳定%两信号之间的角度差（GAP）的大小对估计误差的影响分析%GAPn=0.1:0.1:5; %对两信号之间的角度差（GAP）离散化取值for n=1:length(GAPn)GAP=GAPn(n); %每次循环只取其中一个值DOA=pi*(0/180) pi*(GAP/180);for k=1:M %M为独立重复运行次数estimated,error=TLSesprit(p,L,K,SNR,DOA);errorm(k)=error; %将每次的估计误差存入变量 errorm中，便于求均值enderrorn(n)=sum(errorm)/M; %求多次运行后的估计误差的均值endfigure(5);plot(GAPn,errorn*180/pi,r:*,LineWidth,2);%绘制曲线并适当标注(误差：角度) xlabel(GAP( );ylabel(估计误差( );% title(两信号之间的角度差（GAP）对估计误差的影响);%结论：GAP越大估计越准确，但当 GAP大到一定程度后则估计精度趋于稳定% %单个信号时，信号波达方向分布不同时对估计误差的影响分析% DOAn=pi*(-80/180):(5/180):pi*(80/180); %对信号波达方向离散化取值（80度到 90度时误差太大，因此未取）% for n=1:length(DOAn)% DOA=DOAn(n);p=1; %每次循环只取其中一个值,信号个数 p设为为 1% for k=1:M %M为独立重复运行次数% estimated,error=TLSesprit1(p,L,K,SNR,DOA); %调用 TLSesprit1（一个信号的情况）% errorm(k)=error; %将每次的估计误差存入变量 errorm中，便于求均值% end% errorn(n)=sum(errorm)/M; %求多次运行后的估计误差的均值% end% figure(6);plot(DOAn*180/pi,errorn*180/pi,b:*,LineWidth,2);%绘制曲线并适当标注(误差：角度)% xlabel(DOA( );ylabel(估计误差( ); % % title(DOA（单信号）不同分布对估计误差的影响);% % %结论：越靠近 0度估计越准确，越靠近正负 90度估计误差越大% %估计相关矩阵 R时，减与不减 Rn（无信号时的噪声自相关矩阵）对估计误差的 影响分析(结合信噪比 SNR对估计误差的影响曲线)% SNRn=-15:1:5; %对信噪比 SNR离散化取值% for n=1:length(SNRn)% SNR=SNRn(n);% for k=1:M% estimated,error=TLSesprit(p,L,K,SNR,DOA);errorm(k)=error; %调用减Rn的函数% estimated,error=TLSespritRn(p,L,K,SNR,DOA);errormRn(k)=error; %调用不减 Rn的函数% end% errorn(n)=sum(errorm)/M; errornRn(n)=sum(errormRn)/M;% end% figure(7);h=plot(SNRn,errorn*180/pi,SNRn,errornRn*180/pi,r-.); %绘制减与不减 Rn时的估计误差曲线% legend(R=R-Rn,R); set(h,LineWidth,2); %用图示在图中标明哪条为减或不减 Rn的曲线% xlabel(SNR(dB);ylabel(估计误差( );% % title(减与不减 Rn对估计误差的影响);% % %结论：若减 Rn，主要是在低信噪比时对估计精度的改善较大，当信噪比较大 时二者几乎一样%本文件名为 TLSesprit.m %基于总体最小二乘的 ESPRIT算法（TLS-ESPRIT）的 DOA估计函数% function estimated,error=TLSesprit(p,L,K,SNR,DOA)% 调用格式:estimated,error=TLSesprit(p,L,K,SNR,DOA);% 估计结果(弧度，矢量：p行 1列)：estimated% 估计