esprit算法研究.doc

课程设计报告实验名称:ESPRIT算法研究实验日期:姓 名:学 号: 哈尔滨工业大学(威海)一、 设计任务实现空间谱估计算法,并考察算法性能。二、 方案设计1） 由均匀线阵形式,确定阵列得导向矢量；2） 由阵列导向矢量,对接收信号进行建模仿真;3） 由ESIT算法实现信号DO估计;4） 考察算法性能与信噪比,采样率,观测时间等参数得关系。三、 设计原理3、空间谱估计数学模型空间谱估计就就是利用空间阵列实现空间信号得参数估计得一项专门技术。整个空间谱估计系统应该由三部分组成:空间信号入射、空间阵列接收及参数估计。相应地可分为三个空间,即目标空间、观察空间及估计空间，也就就是说空间谱估计系统由这三个空间组成，其框图见图。图1 空间谱估计得系统结构对于上述得系统结构,作以下几点说明。()目标空间就是一个由信号源得参数与复杂环境参数张成得空间。对于空间谱估计系统，就就是利用特定得一些方法从这个复杂得目标空间中估计出信号得未知参数。(2)观察空间就是利用空间按一定方式排列得阵元，来接收目标空间得辐射信号。由于环境得复杂性，所以接收数据中包括信号特征(方位、距离、极化等）与空间环境特征(噪声、杂波、干扰等)。另外由于空间阵元得影响，接收数据中同样也含有空间阵列得某些特征(互耦、通道不一致、频带不一致等)。这里得观察空间就是一个多维空间，即系统得接收数据就是由多个通道组成,而传统得时域处理方法通常只有一个通道。特别需要指出得就是:通道与阵元并不就是一一对应,通道就是由空间得一个、几个或所有阵元合成得(可用加权或不加权),当然空间某个特定得阵元可包含在不同得通道内。(3）估计空间就是利用空间谱估计技术(包括阵列信号处理中得一些技术,如阵列校正、空域滤波等技术)从复杂得观察数据中提取信号得特征参数。从系统框图中可以清晰得瞧出，估计空间相当于就是对目标空间得一个重构过程，这个重构得精度由众多因素决定,如环境得复杂性、空间阵元间得互耦、通道不一致、频带不一致等。3、2 阵列信号处理首先,考虑N个远场得窄带信号入射到空间某阵列上,阵列天线由M个阵元组成,这里假设阵元数等于通道数,即各阵元接收到信号后经过各自得传输信道送到处理器,也就就是说处理器接收来自M个通道得数据。 (、-1)式中,就是接受信号得幅度，就是接收信号得相位,就是接收信号得频率。在窄带远场信号源得假设下，有 (3、2-2)根据式（、21)与式(3、2-)，显然有下式成立: （3、2-3）则可以得到第个阵元接收信号为 (3、2)式中,为第L个阵元对第i个信号得增益,表示第个阵元在t时刻得噪声，表示第i个信号到达第L个阵元时相对参考阵元得时延。将M个阵元在特定时刻接收得信号排列成一个列矢量,可得(3、25）在理想情况下,假设阵列中各阵元就是各向同性得且不存在通道不一致、互耦等因素得影响,则式(3、24)中得增益可以省略（即归一化1）,在此假设下式（、-5)可以简化为 (3、2-6)将式(3、2-）写成矢量形式如下： (3、-7)式中,为阵列得维快拍数据矢量，为阵列得维噪声数据矢量，为空间信号得维矢量,A为空间阵列得维流型矩阵(导向矢量阵)，且 (3、2-8）其中导向矢量 （3、2-9)式中,为光速,为波长。由上述得知识可知,一旦知道阵元间得延迟表达式,就很容易得出待定空间阵列得导向矢量或阵列流型。下面推导一下空间阵元间得延迟表达式。假设空间任意两个阵元,其中一个为参考阵元(位于原点），另一个阵元得坐标为(x,y,z)，两阵元得几何关系见图，图中“”表示阵元。图2 空间任意两阵元得几何关系由几何关系可以推导出两阵元得波程差为 (3、-10)这里得波程差其实就就是位于轴上两阵元间得延迟、位于y轴上两阵元间得延迟与位于z轴上两阵元间得延迟之与。根据式(3、2-10)得结论,下面给出实际环境中常用得几种阵列及阵元间得相互延迟表达式。(1)平面阵设阵元得位置为，以原点为参考点,另假设信号入射参数为,分别表示方位角与俯仰角，其中方位角表示与x轴得夹角。()线阵设 阵元得位置为,以原点为参考点,另假设信号入射参数为,表示方位角,其中方位角表示与y轴得夹角(即与线阵法线得夹角),则有 (3、2-1)(3)均匀圆阵 设以均匀圆阵得圆心为参考点,则有 (3、2-1)其中方位角表示与x轴得夹角,r为圆半径。3、旋转不变子空间算法原理3、3、1信号模型算法介绍前,首先对信号进行建模。为了推导分析得方便,将波达方向得数学模型做如下理想状态得假设:) 阵列形式为线性均匀阵,阵元间距不大于信号波长得二分之一。2) 存生两个完全相同得子阵,且两个子阵得间距就是己知得。3) 噪声序列为一零均值高斯过程,各阵元间噪声相互独立,噪声与信号也相互独立。4） 空间信号为零均值平稳随机过程,通常为窄带远场信号。) 信号源数小于子阵阵列元数,信号取样数大于子阵阵列元数,以确保子阵阵列流型得各列线性独立。6) 组成阵列得各传感器为各向同性阵元,且无互耦以及通道不一致得干扰。图3-1均匀线阵得数学模型示意图下图给出了均匀线阵得数学模型示意图：3、3、2算法原理对于均匀线阵,相邻子阵间存在一个固定间距，这个固定间距反映出各相邻子阵间得一个固定关系,即子阵间得旋转不变性,而ESRIT算法正就是利用了这个子阵间得旋转不变性实现阵列得OA估计。PRT算法最基本得假设就是存在两个完全相同得子阵,且两个子阵得间距就是已知得。由于两个子阵得结构完全相同,且子阵得阵元数为m,对于同一个信号而言,两个子阵得输出只有一个相位差,=1,2, N。下面假设第一个子阵得接收数据为,第二个子阵得接收数据为,根据前面所述得阵列模型可知 (3、1） （、2)式中,子阵得阵列流型=,子阵2得阵列流型= ,且式中 (3、3）从上面得数学模型可知,需要求解得就是信号得方向，而信号得方向信息包含在与中,由于就是一个对角阵,所以下面只考虑这个矩阵,即 (3、4）由上可知。只要得到两个子阵间得旋转不变关系，就可以方便地得到关于信号到达角得信息。下面得任务就就是从式(、)与式(3、2)中得到两个子阵间得关系。先将两个子阵得模型进行合并,即 (3、5）在理想条件下,可得上式得协方差矩阵 (、6） 对上式进行特征分解可得 （、)显然上式中得到得特征值有如下关=,US为大特征值对应得特征矢量张成得信号子空间,为小特征值对应矢里张成得噪声子空间。对于实际得快拍数据，式(3、7）应修正如下: (3、8） 由前面得知识可知,上述得特征分解中大特征矢量张成得信号子空间与阵列流型张成得信号子空间就是相等得。即 (3、9)此时，存在一个惟一得非奇异矩阵T，使得 (3、0)显然,上述得结构对两个子阵都成立,所以有 (3、11) 很显然 ，由子阵1得大特征矢量张成得子空间、由子阵2得大特征矢量张成得子空间与阵列流型张成得子空间三者相等,即 （3、2) 另外,由两个子阵列在阵列流型上得关系可知 (、1)再利用式(3、11)可知两个子阵列得信号子空间得关系如下： （、14)式(3、13)反映了两个子阵列得阵列流型间得旋转不变性,而式(、4)反映了两个子阵得阵列接收数据得信号子空间得旋转不变性。如果阵列流型A就是满秩矩阵,则由式(、14)可以得到 (3、5)所以上式中得特征值组成得对角阵一定等于，而矩阵T得各列就就是矩阵特征矢量。所以一旦得到上述得旋转不变关系矩阵,就可以直接利用式(3、4)得到信号得入射角度。、4 标准得旋转不变子空间算法有上节得知识可知， ES算法得基本原理就就是利用式(3、4)得旋转不变性，常规得旋转不变子空间算法就就是利用上述得基本原理求解信号得入射角度信息。下面就分析解这个等式得两种最经典、应用最广泛方法:最小二乘（LS)法与总体最小二乘(TLS)法。3、4、1 最小二乘法由最小二乘得数学知识，我们知道式（3、14)得最小二乘解得方法等价于,约束条件 （3、16)因此最小二乘法得基本思想就就是使校正项尽可能小，而同时保证满足约束条件。为了得到S解,将式(3、14）代入式(3、16)即得 (3、17）对上式进行展开可得 = (3、1)上式对求导并令其等于,可得 （3、1）上式得解显然有两种可能:（) 当满秩时,也就就是子阵1得信号子空间得维数等于信号源数时,则上式得解就是唯一得,可得上式得最小二乘解 (、20)（2)当不满秩,即时,也就就是信号源间存在相干或相差时,则存在很多解,但我们却无法区别对应于方程得各个不同得解,可以称这些解就是不可辨识得，解得不可辨识性就是我们需要解相干得原因所在。下面给出L-PT算法得求解步骤:1.由两个子阵得接收数据,分别得到两个子阵得数据协方差矩阵;2.对矩阵对R, 进行特征分解,从而得到两个数据矩阵得信号子空间与;3.按式(、20)得到矩阵,然后对其进行特征分解、得到个特征值,就可得到对应得N个信号得到达角。当考虑嗓声影响时，上述基于最小二乘算法得估计都就是有偏得,这就就是为什么需要考虑总体最小二乘ESPRT算法得原因。3、4、2 总体最小二乘法 我们知道，普通最小二乘得基本思想就是用一个范数平方为最小得扰动去于扰信号子空间,目得就是校正中存在得嗓声。显然这就存在一个问题:如果同时扰动与,并使扰动范数得平方保持最小,就是否可以同时校正与中存在得嗓声?答案就是肯定得,这就就是总休最小二乘(TLS)得思想。它考虑得就是如下矩阵方程得解: (3、27）显然上式可以改写成 (3、28)所以TLS得解等价于 (3、2)定义如下一个矩阵,再结合上述分析过程。我们发现其实就就是寻找一个得酉矩阵F，便得矩阵F与正交,也就说明了由张成得空间与或列矢量张成得空间正交。所以矩阵F可从得特征分解中得到。因为 (3、0)式中得就是由特征值构成得对角矩阵,E就是与其相应得特征矢量构成得矩阵。即 (、31)令就是由对应特征值为0得特征矢量构成得矩阵、它属子噪声子空间，所以只要选择矩阵F使之等于、,即可满足上面提到得要求。即有 (3、3）可得 (、)如果令，则 (3、4）上式说明得特征值即就是对角线元素。这说明通过构造一个矩阵就可得到有关信号角度得信息、而这个矩阵得构造可通过式（3、0)得到,即 (、5)下面直接给出LS-ESPIT算法得求解步骤:1.由两个子阵得接收数据,， 由式（3、8）得到数据协方差矩阵;2通过矩阵对于得广义特征分解，得到维数为得信号子空间；3由构造矩阵,并按式防(3、30)进行特征分解得到矩阵E,然后再按式（3、31)将矩阵分为四个小得矩阵；.按式(3、35)得到矩阵,然后对其进行特征分解,得到N个特征值，就可得到对应得N个信号得到达角。通过分析,我们可以得到标准ESRIT算法得计算过程如下:(1)通过特征值或奇异值分解(ED或SVD)分别估计两个存在旋转不变关系得子阵得信号子空;(2)用上述得LS、S等方法求解式（3、4)所示得不变等式；（3)计算得特征值,其中如式（3、3)所示。然后利用式(3、4)求解人射信号得角度信息。就ESPRT算法而言,TLS算法与LS算法性能基本一致，只就是在低信噪比情况下L算法性能略好。四、 仿真结果 主要分析各个参数对估计误差得影响，误差函数定义如式(1):4、1 信噪比 R对估计误差得影响分析 首先对信噪比NR离散化取值,然后求得不同信噪比下得误差,从而绘制出误差随信噪比改变得函数曲线如图2 所示,图 2 中信噪比NR从- 15 取到 5，间隔为 1,运行次数为 100 次,其余条件如题中所述。由图2 可知,随着信噪比得增大,估计误差会越来越小,即估计精度会越来越高。当待估计得信号方位角相差比较小时,估计得误差也会相应得增大。另外,若两信号为相干信号,则此方法将不能对其进行正确得估计。 4、2 阵元数 对估计误差得影响分析与 、1节类似，首先对阵元数 离散化取值,然后求得不同阵元数下得误差,从而绘制出误差随阵元数改变得函数曲线如图 3 所示,图 中阵元数从 K+1取到K+25,间隔为 1,运行次数为 10次，其余条件如题中所述。由于阵元数L 需大于信号个数K才能正确估计,故取值中含有信号个数 K。由图 3 可知,随着阵元数得增加,估计误差会越来越小,即估计精度会越来越高,但当阵元数大到一定程度后，对估计精度得影响则会慢慢得减小。 、 采样点数N对估计误差得影响分析与 4、1节类似,首先对采样点数 N离散化取值,然后求得不同采样点数下得误差,从而绘制出误差随采样点数改变得函数曲线如图4 所示,图 中采样点数从 1取到 00,间隔为,运行次数为 10 次,其余条件如题中所述。由图 4可知,随着采样点数得增加,估计误差会越来越小,即估计精度会越来越高。 估计误差(角度)、4 两信号之间得角度差(）对估计误差得影响分析由于采用 ESPR 算法对 O进行估计,若两信号得方位距离较近时,虽然能得出估计结果，但估计得精度会大受影响。因此，为了分析两信号之间得不同间隔会对估计精度造成多大得影响，绘制不同 GA下得估计误差曲线如图 5所示。处理方法与4、1 节类似,图 中 GA(单位为度)从 0、1 取到,间隔为 、 1,独立运行次数为00次,其余条件如题中所述。由图 5 可知,GP越大估计越准确，但当 A大到一定程度后则估计精度趋于稳定。 4、 单信号 DA不同分布对估计误差得影响分析信号波达方向(）得取值区间为0度到 90度,若只考虑只有一个信号得情况,则当信号得DA不同时,估计误差也会不一样。因此,为了分析不同得 OA会对估计精度造成多大得影响,绘制不同DA下得估计误差曲线如图 6所示。处理方法与4、 1 节类似,图6 中 GA从 0度取到 80度,间隔为 度,独立运行次数为 10次，其余条件如题中所述。由图可知,OA越靠近 0 度估计越准确,越靠近正负 9 度估计误差越大。且仿真结果表明,当OA在正负90 附近时,估计误差太大,因此，为了不影响估计结果显示效果,故在图中未绘制正负90 度附近得估计误差。 、6减与不减噪声方差()对估计误差得影响分析由于有噪声得影响，因此在估计信号自相关矩阵R时,若将无信号时得自相关矩阵 R减去,即相当与减去估计出噪声方差，则估计得精度会有所提高。结合信噪比 SNR对估计误差得影响,绘制减与不减噪声方差两种情况下估计误差随 SNR得变化曲线如图7所示,图 7 中 SNR从-15dB到 5B,间隔为 1dB,独立运行次数为 00次。仿真结果表明,若减Rn,主要就是在低信噪比时对估计精度得改善较大,当信噪比较大时二者几乎一样。 五、 程序清单%本文件名为 draTLSeprt、%分析基于总体最小二乘得 ESP算法（LS-EIT)得 O估计得性能%cr;clc；loe al;%清除变量,清屏，关闭所有绘图窗口% 调用格式:estimaed,error=LSerit(p,L,K，NR,D);%估计结果(弧度，矢量：p行 1列)：estimate 估计误差(弧度,标量:均方误差):ero%信号个数： 阵元数:% 快拍数:%信噪比:SNR% 波达方向(弧度,矢量：行 列)：OA 2; L8; K=1; SNR5; DOA=pi*(-10/180）p(2/80);%显示估计结果%=0; %设定独立重复运行次数DOA=pi*（0/80)pi*(30/180)；%波达方向(弧度,矢量:行1列)p=lengh(OA); L8;K=10； SNR; %参数设置,etimate,error=LSespr(,L,SNR,D）; 函数调用polr(estimated,1 1,r);%在极坐标中显示估计结果(必须先转化为弧度)=title(）;set（h,ting，TLS-EPRIT： 估计值: ,nm2str（estimted);h1=xll(）;st(h,tring，信号DOA(度）: ,um2str(OA*180/pi)；% %阵元数 对估计误差得影响分析%Ln=p+1:：p25; %阵元数 L需大于信号个数p才能正确估计% for n=1:egh(Ln）%L=Ln(n);%o k1:% timated,erro=LSesprit(，L,K，NR,DOA);% errorm（k)=ror; %将每次得估计误差存入变量 rr中，便于求均值% end erorn（)=sum（rorm); %求多次运行后得估计误差得均值% e% figre（2);plot(n,errorn*10/pi,：*,iWth，);%绘制曲线，并适当标注% xlabel(阵元数 L);ylabe(估计误差（ );ttle(阵元数 对估计误差得影响)；% % 结论：阵元数越大估计越准确,但当L大到一定程度后则估计精度趋于稳定% %快拍数 K 对估计误差得影响分析% Kn1:10:200; %对快拍数离散化取值%or =：lnth(n）% K=Kn（)； fr =1：M stiatd，errr=epit（p,L,SNR,DOA）; errom(k）=eror;%将每次得估计误差存入变量 errm中,便于求均值% end% ern1(n)=um(erorm）/M； %求多次运行后得估计误差得均值 en% figre（3);pt(Kn,erorn1*180/i，r:*,LieWid,)； %绘制曲线并适当标注% labl（快拍数 K)；ylabel(估计误差（ );title(快拍数 对估计误差得影响);% % 结论:快拍数 K越大估计越准确,但当 大到一定程度后则估计精度趋于稳定% %信噪比 NR对估计误差得影响分析% SRn=-15:1：15; %对信噪比 NR离散化取值% o =:ngh(SNRn）%SNR=SN（n);%fo k=：M % estaed,erro=LSsprit(,K，SR,DO)；% ero(k)=eror； 将每次得估计误差存入变量 err中,便于求均值% end%errr(n)=um(erorm)/M; 求多次运行后得估计误差得均值% end% fige(4);lt(SNR,erorn*10/i,r:*,LinWdt,2);绘制曲线并适当标注(误差:角度) xlbel(SNR);yabel(估计误差( ）);title（SR对估计误差得影响);%结论:信噪比 SNR越大估计越准确，但当信噪比 R大到一定程度后则估计精度趋于 稳定%两信号之间得角度差（GAP)得大小对估计误差得影响分析%n0、1:0、：5;%对两信号之间得角度差(GAP)离散化取值or n=：length(GAPn)GAP=GAPn(n); %每次循环只取其中一个值DA=pi*（0/180)pi（GP/180);fork=1:M为独立重复运行次数estimatd，rror=Tespri（p,L,SN，O)；rrrm(）eror； 将每次得估计误差存入变量erro中,便于求均值nerorn（)=s（errr)M; %求多次运行后得估计误差得均值efigure(5)；pot(GP,er*10/pi，r:*,LneWdth，2);%绘制曲线并适当标注（误差:角度） label(GA( );ybel(估计误差（ ）)；%ttl(两信号之间得角度差(GP)对估计误差得影响);结论:GA越大估计越准确，但当GAP大到一定程度后则估计精度趋于稳定% %单个信号时,信号波达方向分布不同时对估计误差得影响分析% An=pi（-80/80):(5/18)：pi*(8010）; %对信号波达方向离散化取值(8度到 90度时误差太大,因此未取)% for n=1:lenth(DAn)%DOAO（n）;p=1；%每次循环只取其中一个值，信号个数 p设为为1% fo =1:M %为独立重复运行次数% stimate,erro=TLSesri1(p,L,K，SNR，DA)；%调用 LSsrit1(一个信号得情况)% erom(k)=errr; %将每次得估计误差存入变量 errorm中，便于求均值%end% rn()sum(ero)M；%求多次运行后得估计误差得均值%nd% figure(6);ot(On*80/i,eorn*180/pi,b：*,Lineth,);%绘制曲线并适当标注(误差:角度) xlbel(OA（)）；ylabe(估计误差（ ）); % %ttle(O(单信号）不同分布对估计误差得影响)；% % 结论:越靠近 0度估计越准确,越靠近正负 9度估计误差越大%估计相关矩阵时，减与不减Rn（无信号时得噪声自相关矩阵)对估计误差得影响分析（结合信噪比 R对估计误差得影响曲线)% Nn=-15:5;%对信噪比 SNR离散化取值% f =1:length(SNRn)% N=NR（n);% ork=1：M% estimated，eroTLesprit(p,L，K,NR，DA);errorm(k)=erro； %调用减R得函数% eima,eror=TLSespriRn(p,L,K,SNR,DOA);errrR(k)rrr; %调用不减Rn得函数% nd% erorn(n)=um(erm）/M; errorn()=sum(orm）/；% ed%figure(7);h=lot(SNR,error*180/，SNRn,erroRn*180/pi,-、); %绘制减与不减 Rn时得估计误差曲线 len(=-R,R)；se(h,Linidth,); %用图示在图中标明哪条为减或不减Rn得曲线% lae(SR(d);la(估计误差();% % title(减与不减 n对估计误差得影响); %结论:若减 Rn,主要就是在低信噪比时对估计精度得改善较大,当信噪比较大 时二者几乎一样%本文件名为TSespit、%基于总体最小二乘得 PRIT算法(TLS-EPRI）得 DA估计函数%fnctn sited,error=TLSespi（p,L,K,S,DO)调用格式：estte,rror=TLSesrit(p,L,K，NR,OA);% 估计结果(弧度,矢量:p行 1列):timae% 估计误差(弧度,标量:均方误差)：error 信号个数：p 阵元数：% 快拍数:K% 信噪比:N% 波达方向（弧度,矢量：p行列)：DOA p2； L8; K=00; NR5; DOA=i*（-10/180) pi*(20/18);%参数设置%dbc=/2;%阵元间隔d与信号波长之比 / 1/2thta=*pi*bcsn（DOA); 信号方位参数 thtaOigTi/4; /6； 信号频率Dnsqt(1/（21（N10);%噪声标准差%估计相关矩阵 %A=xp(j（0:1）*thea); %表示出阵列方向矩阵=（j*Omi*（:K1）);%构造信号源矢量X=*S; %构造阵列输出矢量(无噪)Nien(sq（2)2)*(andn(L,K)+jrann（L,K)）;%加入复噪声=X+oise; %构造阵列输出矢量R=zes（L,L);Rn=zeos(L,L)； %初始化为零，加快运行速度for i=1:KRn=RnNoie(:，)*Noise(:,i);R=R+（:,i)*Y(:,）；endR/;RnnK；%求得相关矩阵 R(有信号)与 R(无信号)R1=R-n; %减小噪声对估计精度得影响V，Dei（R1); %相关矩阵特征分解%(D中特征值已经按从小到大得顺序排列,即V中前L-p个为噪声对应得特征向量)%构造矩阵 S%=V（:,-p+：L)； %L行p列（中得列为R中 p个大特征值对应得特征向量)S1=S(1:L-1，:）;%将 S得前 -1行构造 S1S2=(2：L,:); 将 S得后 L-行构造 SS12=1 S2; %利用 S与S2构造 S（L1行 2p列）SSS12*12; %2K 行 2列，D1=g(SS）; %特征分解%求解估计结果%11=U(1:p,1：); %p行p列,U2（1:2*p,1:p）; %行 p列LS-11*in(1); %p行 p列,U1与 U2构成 U得噪声子空间deig（TL)；%特征分解,求 TLS得特征值tated=(sor(sin（age(d)/)；输出估计(已从小到大排序)结果(弧度)errsrt(sum(estmaedsort(A)、2）/p);求出估计误差(弧度):均方误差