相干信号空间谱估计测向Matlab仿真研究

相干信号空间谱估计测向Matlab仿真研究1.引言由于多径传播、电磁干扰等因素的影响，相干信源存在的电磁环境是经常碰到的。当空间存在相干源时，经典的超分辨DOA估计方法：MUSIC算法和ESPRIT算法，已经失去了其高分辨性能优势，有时甚至不能正确地估计出信源的真实方位。新MUSIC算法在空间不存在相干源时，其估计性能基本上是和MUSIC算法是接近的，但若有相干源存在时，其估计性能也是大大降低。因此，若将其用于相干源，必须和经典的MUSIC算法一样，首先对阵列输出的协方差矩阵进行各种去相干处理，然后再采用新MUSIC算法实现对相干信源的DOA估计。基于加权空间平滑的MUSIC算法，该算法充分利用了子阵输出的自相关信息和互相关信息，将阵列协方差矩阵的所有子阵阵元数阶子矩阵进行加权平均，而权矩阵的选取以平滑后等价的信源协方差矩阵与对角阵的逼进为约束条件，以期对相干信源最大限度地去相干，改进常规空间平滑算法对相干源的分辨力。基于此本文提出的基于加权空间平滑的新MUSIC算法，以实现对相干源最大限度的去相干，实现相干源的高分辨DOA估计。2.窄带阵列相干源的数学模型和空间平滑算法2.1窄带阵列相干源的数学模型对于M元均匀线阵，阵源间距为d，且假设均为各向同性阵元。阵列远场中在以线阵轴线法线为参考的(1,2,...)kPkθ处有P个窄带点源以平面波入射，以阵列第一阵元为参考点，某一特定信号到达线阵时，各阵元接收信号间仅仅存在因波程差引起的相位差。阵列接收的快拍数据可由下式表示为:tttX()=A()S()+N()(2-1)(2-1)式中()tX为M1快拍数据矢量，是阵元输出信号，1()[(),......()]TMtxtxtX。()tN为M1阵列阵元噪声矢量，且是与信号源不相关的白噪声，均值为0，方差为2n，1()[(),......()]TMtntntN。()tS为输入信号矢量1()[(),......()]TMtststS。()A为阵列的流形矩阵,如式(2-2)所示。1()[(),......()]PAaa=11()(1)()1...jjMee22()(1)()1...jjMee.........()(1)()1...ppjjMee(2-2)（2()sin()iid,1,2,...iM）向量()ia称为第i个信号源的方向向量。矩阵A列向量即是阵列方向向量的集合，它表示所有信源的方向，称为阵列的方向矩阵。阵列的协方差矩阵R定义为2[()()]HHEttSRXXARAI,式中[()()]HEttSRSS为信源的协方差矩阵。当空间信源互不相干时，对协方差矩阵R进行特征分解，构造信号子空间和噪声子空间，利用其正交性直接采用MUSIC或ESPRIT算法进行DOA估计，令其特征值分解为HsRUU,式中，221(,......,)Mdiag，如果信号不相干，则A为列满秩，()rankPHSARA，于是221(,...,,0,...0)PdiagHHSUARAU，将R同时左乘HU和右乘U，可以得到：22221(,...,,0,...0)HHPdiagHHSURU=UARAUUU=I，所以自相关矩阵R的特征值为2222221,...,,...P，前P个特征值为信号特征值,后M-P个特征值为噪声特征值,由信号特征值对应的所有特征向量形成的矢量空间称为信号子空间,信号子空间同时可以用信号的方向矢量()ia表示。由噪声特征值对应的特征向量形成的矢量空间称为噪声子空间,二者组成信号自相关矩阵的信息空间,利用信号矢量与噪声空间的所有矢量都是正交的这一特性来估计来波信号的到达角（如式2-3所示）,即MUSIC法。1=()HHP（）（）aUUa（2-3）代码如下：clearall;clc;p=2;%入射信号数目M=4;%阵元个数fc=1e9;%入射信号中频为1GDOA=[-20,40]/180*pi;%信号入射DOAfs=3*fc;%采样频率N=512;%采样个数snr=10;%信噪比T=1/fs;%采样时间间隔c=3e8;%波速d=c/fc*0.5;%阵元间距t=0:T:(N-1)*T;%采样时间区间s1=sqrt(2)*cos(2*pi*fc*t);%信号数据s2=sqrt(2)*cos(2*pi*(fc+5e8)*t);ss=[s1;s2];s=ss(1:p,:);%%%%%%%%计算阵列流形矩阵AA=zeros(M,p);fork=1:pforkk=1:MA(kk,k)=exp(-j*2*pi*fc*(kk-1)*d*sin(DOA(k))/c);endend%%%%%%%%阵列接收数据y=A*s;y=awgn(y,snr);R=y*y'/N;%%%%%%%%计算噪声子空间[v,dd]=eig(R);if(dd(1,1)dd(2,2))Un=v(:,p+1:M);elseUn=v(:,1:(M-p));enddo=-90:90;pu=zeros(1,length(do));kg=1;fork=-90:90a=zeros(M,1);forkk=1:Ma(kk,1)=exp(-j*2*pi*fc*(kk-1)*d*sin(k/180*pi)/c);endpu(1,kg)=1/(a'*Un*Un'*a);kg=kg+1;endplot(do,10*log10(abs(pu)),'-r','linewidth',2);gridon;title('MUSIC测向');xlabel('波达方向');ylabel('MUSIC谱');disp('MUSIC测向结果:');fork=1:p[k1,k2]=max(pu);DOA_guji(k)=(k2-1)-90;pu(k2)=0;endDOA_guji仿真结果如下：图1无相干信号源时的MUSIC测向-100-80-60-40-20020406080100-1001020304050MUSIC测向波达方向MUSIC谱但若存在相干信源时，阵列输出信号协方差的秩()rankPR，对信号协方差矩阵进行特征值分解后，得到的较大的特征值个数小于P，而特征值为2的个数将大于M−P。与此相对应的信号子空间的向量也少于P，即特征向量展开的信号子空间的维数少于1()[(),......()]PAaa的列数。对某些相干源的方向矢量()ia，1,2,...iP将不正交于噪声子空间，不出现零点，所以，有些源在空间谱曲线中将不呈现峰值，造成谱估计的漏报。在这里我们不妨将上面程序的信号源是S1和S2设置成相干信号，S2=2S1；仿真结果如下：图2有相干信号源时的MUSIC测向因此，我们要对阵列输出的协方差矩阵首先进行预处理，使其阵列协方差矩阵的秩恢复为信号元数P，然后再采用MUSIC或ESPRIT算法。-100-80-60-40-20020406080100-6-5-4-3-2-1012MUSIC测向波达方向MUSIC谱2.2空间平滑算法（spatialsmoothing）2.2.1前向空间平滑算法将M个阵元的均匀线阵，分成相互交错的P个子阵，每个子阵包含的阵元数为m个，即满足M＝p+m-1。信号源数为N。图3前向空间算法原理图如图3所示，取第一个子阵(最左边的子阵)为参考子阵，那么各个子阵的输出矢量分别为:11222311[,,...,][,,...,]...[,,...,]fmfmfpppMxxxxxxxxxXXX（2-4）对于第k个子阵有：(1)11()[,,...,]()()()fkkkkkmmktxxxttXADsn（2-5）其中：12sin()0...0djeD22sin()0...0dje............2sin()00...Ndje（2-6）那么该子阵的数据协方差矩阵为：(1)(1)2()(())kkHkmsmRRADADI（2-7）其中，mA是一个m×p的参考子阵(通常取第一个子阵)的导向矢量矩阵，1()[(),......()]mmmNAaa，22sin()sin()()[1,,....,]kkddjjTmkeea，sR为信号的协方差矩阵，HsEssR。前向空间平滑技术是通过求各个子阵协方差矩阵的均值来实现的，即取前向平滑修正的协方差矩阵为:11pfkpkRR（2-8）可以证明，当满足mN,pN时，前向空间平滑数据协方差矩阵fR是满秩的。即可以通过特征分解求得相应的信号子空间和噪声子空间。2.2.2前后向空间平滑算法如果按照图4划分阵列，即称为后向平滑的方法划分子阵，那么各个子阵的输出矢量为:图4后向空间平滑算法原理图11121211[,,...,][,,...,]...[,,...,]bMMMmbMMMmbpmmxxxxxxxxxXXX(2-9)那么，第k个子阵的数据矢量为：*12()[,,...,]bkMkMkMmktxxxX（2-10）比较前向平滑和后向平滑的数据矢量，可以得到前向平滑中第k个子阵与后向平滑中第p-k+1个子阵之间存在如下关系：**(1)**1()(())()()bfkpkkmktJXtJstJntXAD（2-11）其中J为m的交换矩阵。00...1J0......0...1...010...0，所以后向平滑第p-k+1个子阵的数据协方差矩阵为:*(1)*(1)*21()()bkkHHpkmmsRJADRDAJI(2-12)那么后向空间平滑修正的数据矩阵为:111pbbpkkpRR(2-13)取前向平滑和后向平滑数据协方差矩阵的平均，即前后向空间平滑的数据矩阵，即2bffbRRR(2-14)同样可以证明，当满足mN,pN时，后向空间平滑数据协方差矩阵bR是满秩的。3.Matlab仿真程序代码clearall;clc;p=2;%入射信号数目M=8;%阵元个数L=5;%将阵列划分为相互重叠的L个子阵m=M-L+1;%每个子阵中的阵元个数fc=1e9;%入射信号中频为1GDOA=[-20,40]/180*pi;%信号入射DOAfs=3*fc;%采样频率N=100;%采样个数snr=10;%信噪比T=1/fs;%采样时间间隔c=3e8;%波速d=c/fc*0.5;%阵元间距t=0:T:(N-1)*T;%采样时间区间s1=sqrt(2)*cos(2*pi*fc*t);%信号数据s2=2*s1;s3=sqrt(2)*cos(2*pi*(fc+5e8)*t);ss=[s1;s2];s=ss(1:p,:);%%%%%%%%计算阵列流形矩阵AA=zeros(M,p);fork=1:pforkk=1:MA(kk,k)=exp(-j*2*pi*fc*(kk-1)*d*sin(DOA(k))/c);endendy=A*s;y=awgn(y,snr);%阵列接收数据R=y*y'/N;[v,dd]=eig(R);%对协方差矩阵进行特征值分解%%%%%计算噪声子空间if(dd(1,1)dd(2,2))Un=v(:,p+1:m);elseUn=v(:,1:(L-p));enddo=-90:90;pu1=zeros(1,length(do));kg=1;fork=-90:90a=zeros(M,1);forkk=1:Ma(kk,1)=exp(-j*2*pi*fc*(kk-1)*d*sin(k/180*pi)/c);endpu1(1,kg)=abs(1/(a'*Un*Un'*a));kg=kg+1;endholdon;rf=zeros(L,L);rb=zeros(L,L);Z=zeros(L,N);X=zeros(L,N);Rf=zeros(L,L);Rb=zeros(L,L);%%%%