二维DOA估计算法与对比实验

Projectreportabouttwo-dimensionalDOAestimation题目：考虑一个20阵元数的双线性均匀线阵，现有三个信源入射，它们的波达方向⑴OA）分别是（10o,10o）,（20o,20o）和（30o,30。），请用2D-MUSIC算法，2D-ESPRIT算法，2D-Capon算法,2D-PM算法以及DOA矩阵方法来估计这些信源的波达方向。1.信号接收模型如图1，考虑N个不同二维DOA（0〃,。〃）,n=1,2,・..,N的窄带远场信号s”（t）,在离散时间t入射有2M个传感器的双平行均匀线阵时。x轴和y轴上信源的方向矢量分别为：1Z其中X为波长A=[a（0「q）,11'asourceyM/'-其中X为波长A=[a（0「q）,11'asourceyM/'-X,•aSubarray2Subarray1ej2兀dsinO/os虬'入£j2兀(M-1)dsinO/os虬项T(1)（0@）=[1ej2兀dsin0^sin。「'Xd是阵元间距，x轴M个阵元对应方向矩阵为（0,。），…，ax（0n,中）]，具体表示为：(2)1J

ej2兀dsinOjcosQ'Xej2兀dsin0Ncos^N入eCmxn(3)£j2兀(M-1)dsin01cos虹入£j2兀(M-1)dsin02cos。2小...ej2兀(M-1)dsin0NcosON项y轴2个阵元对应方向矩阵为A=y体表示为:a（0,。），a（0,。），.•.，a（0,。）]，具y11y22yNNJeC2xN(4)TOC\o"1-5"\h\z1eC2xN(4)ej2兀dsinqsin虹人ej2兀dsin92sin虹'人ej2兀dsin©Nsin巾/人双平行线阵中子阵列1的接收信号为x(t)=As(t)+n(t)(5)子阵列2的接收信号为x2(t)=A?s(t)+n2(t)(6)其中n1(t)和n2(t)分别表示子阵列1和2的与信号不相干的加性高斯白噪声，鱼=diagej2兀dsin鱼=diagej2兀dsin01sin虹入，...，ej2兀dsin9Nsin虹入，s(t)=[s1(t),表示信源矢量。可得整个双平行线阵的接收信号为x(t)=:?)卜2s(t)+n，"2」s(t)]TeCnx1(7)式(7)也可表示为(t)(t)=[AOA]s(t)+n(t)xi*-yx(8)n(t)=[气(t)n2(t)]T；(9)yAoA表示A和A的Khatri-Rao积。由Khatri-Rao积的定义，接收信号可以表示为(t)=「a(9,。)®a(9,。n(t)=[气(t)n2(t)]T；(9)yAoA表示A和A的Khatri-Rao积。由Khatri-Rao积的定义，接收信号可以表示为(t)=「a(9,。)®a(9,。)，.•.，a(9,。)®a(9,。)]s(t)+n(t)=U£Uh+U£UhSSSNNN(10)2.二维DOA估计算法2D-MUSIC算法2D-MUSIC算法是通过二谱峰搜索来确定声源波达方向，假设噪声是独立同分布为零均值的高斯过程，且信号和噪声不相干。R为阵列协方差矩阵，可以R=E{(t)x(t)h}=(AOA)R(AOAA+E{(t)n(t)h}yxsyx式中旧「n(t)n(t)H序212MeC2Mx2M表示噪声协方差矩阵，E是期望操作，上标H表示共轭转置，信号协方差矩阵Rs=E〈(t)s(t)htEcnxn。Us表示(10)式特征分解后N个大特征值对应的特征向量张成的空间，称作信号子空间；UN表示

（2M-N）个小特征值对应的特征向量张成的空间，称作噪声子空间，且与信号子空间正交。在空间谱估计算法中，2D-MUSIC算法通过谱估计式的峰值二维搜索，来确定入射的方向信息，2D-MUSIC算法空间谱函数为P2D-MUSIC（11）'"）fa（0,©）®a（0,8）］hUUh「a（0,时®a（0,时一P2D-MUSIC（11）1—y式中U州是噪声子空间，经过式（10）特征分解可以得到。则DOA角度估计可以通过二维谱峰搜索使得式（11）输出值最大化的角度信息得到。现考虑20阵元数的双平行均匀线阵，有三个信源入射，它们的波达方向（DOA）分别是（10o,10o）,（20o,20o）和（30o,30o），快拍数L=200。在信噪比相对较高的噪声环境中，2D-MUSIC谱函数能够以图2（b）所示的等高线来表示声源DOA，可以看到波峰的大概位置。然而，当信噪比较低时，峰值可能被扭曲，估计的DOA可能偏离真实值，如图2（a）所示。4030201010204030201010203040elevation（degree）图2（a）信噪比SNR=10dB图2（b）信噪比SNR=20dB2D-Capon算法在空间谱估计算法中，2D-Capon算法也是通过谱估计式的峰值二维搜索，来确定入射的方向信息，2D-Capon算法空间谱函数为P2D-Capon（0’0）=7=―7-—X/a（0,。）®a（0,P2D-Capon（0’0）=7=―7-—X/a（0,。）®a（0,0）yxR-ia（0,。）®a（0,。）Lyx」（12）4030201010203040elevation(degree)图4030201010203040elevation(degree)图3(a)信噪比SNR=10dB2D-ESPRIT算法基于空间旋转不变性。ESPRIT算法作为一种子空间角估计方法，需要对接收数据的交叉谱矩阵进行特征值分解（EVD）或奇异值分解（SVD）。该算法要求接收信号具有一定的不变性结构，而不变性结构通常依赖于已知的两束相同子阵列之间的平移。该算法不需要对谱峰进行搜索，可以获得自动配对的参数估计。该算法为DOA估计提供了封闭解。我们分别构造矩阵A1和A2:（13）（14）A1=[AD1（Aj,A2=[AD（13）（14）A1和A2之间相差了一个旋转因子也，即A/A=，其中①=diagej2ndsin01sin8「A,•••,ej2冗dsin0^sin巾川通过R.特征分解得到信号子空间U,，由于阵列的平移不变性，U,可以分解为两部分，构造矩阵E=气（1：M,:），E=气（M+1:2M,:）。其中E表示取U,的1〜〃行；E表示取Us的必+1〜2〃行。E和E可以表示为TOC\o"1-5"\h\zE=A1T（15）E=A2T=A1QT（16）其中T为KxK维的满秩矩阵。由式（15）和（16）可得E=ET-10T=E呼（17）式中，呼=T-1皿T。由此可知，E和E张成相似的子空间，且矩阵鱼的对角元素为呼的特征值。根据最小二乘法，呼的估计值可由式（18）得出：人、呼=E+E（18）

对W进行特征值分解，得到鱼的估计值击，利用W的特征向量，我们得

yy一■一一-A到了矩阵T的估计值T。由于W和鱼的特征值相同，对W进行特征值分解得到ej2Wsinensi吐爪，yu=sinesin。的估计值可通过式（19）得到，即U=angle以永/2兀dn=1,2,...,N（19）式中，广是矩阵W的第n个特征值，angle（•）表示取复数的相位角。之后对信号n子空间US进行重构，U；=UST-1,即（20子空间US进行重构，U；=UST-1,（20）AD（A）yMx-由U构造矩阵E'=U（1:2（M-1），:），E=U'（3:2M,:）。利用与前面相似TOC\o"1-5"\h\zSxSyS方法，我们得到v=sinvcosv的估计值\o"CurrentDocument"v=angle（8）X/2兀dn=1,2,...，N（21）nn式中，七是矩阵（exAEy的第n个特征值。根据式（19）和式（21）可知，七和七的估计值由相同的列模糊，我们能够得到自动配对的方位角和俯仰角，（七,七）配对完成，得到二维DOA估计值:U+。2）e=arctan（U/v）现考虑20阵元数的双平行均匀线阵，有三个信源入射，它们的波达方向^=arcsin（n（22）（23）（DOA）分别是（10o,10o）,（20oU+。2）e=arctan（U/v）现考虑20阵元数的双平行均匀线阵，有三个信源入射，它们的波达方向^=arcsin（n（22）（23）2D-PM算法对矩阵A,分块可得(24)(24)A=a、1

XA1-x2」在假定阵列无空间模糊（即A^列满秩）的情形下，A广Cn、n为非奇异矩阵，Ax2CCW-N）阵，Ax2CCW-N）xN，可以由A线性变换得到。X1因此该阵列方向矩阵A可以写成AX1D1?ytAX2D2ytAX1D1(AytAx2D2(A「AX1PHAAX1（25）X1式中，P为传播算子，IeCnxn为单位矩阵。式（10）中阵列接收信号协方差矩阵R*可以写成两部分R=[GH]（26）k其中G其中G是由Rk的第1列到第N列组成,H是Rk剩下的（2M-N）列组成。定义传播矩阵P=（GHG）-1GHHeCNx（2M-N），定义矩阵（27）AD(A)AD(A（27）构造矩阵E=E（1:M,：），E广E（M+1:2M,:）。其中E表示取E的1〜"行;E^表示取E的M+1~2M行。用类似2D-ESPRIT的算法，求得u=sin0sin。的估计值：U-angleGL/2兀dn=1,2,...,N（28）TOC\o"1-5"\h\z然后对信号子空间E进行重构，有E'=ET-1，同理我们得到v=sinvcosu的估nnn计值\o"CurrentDocument"v=angle（£）X/2丸dn=1,2,...,N（29）nn式中，£n是矩阵（E：>E；的第n个特征值。根据式（28）和式（29）可知，七和七的估计值由相同的列模糊，我们能够得到自动配对的方位角和俯仰角，化,与）配对完成，得到二维DOA估计值：—）U2—）U2+v2/

kk0=arctan（U/V）j=arcsin(n（30）（31）现考虑20阵元数的双平行均匀线阵，有三个信源入射，它们的波达方向（DOA）分别是（10。,10。）,（20。,20。）和（30o,30o），快拍数L=200。2D-PM算法的估计效果如图5所示。3.算法性能比较实验3.算法性能比较实验在以下仿真实验中，估计性能的评价方式为均方根误差（RMSE），定义为(32)RMSE=1£_L"一"+«Q

n=1j=1

(32)其中0,4是第j次蒙特卡洛试验的方位角和俯仰角估计值，酊，『是角度k,jk,jk,jk,j真实值。实验一：信噪比SNR=5dB，快拍数从50到500变化，变化步长50。运行500次蒙特卡洛（MC）试验，上述几种算法的对比仿真效果如图6所示。实验二：快拍数L=200不变，信噪比SNR从5dB变化到20dB，变化步长为5dB，运行500次蒙特卡洛（MC）试验，上述几种算法的对比仿真效果如图7所示。图7快拍数L=2004.DOA矩阵算法4.1阵列结构和信号模型信号接收阵列由如图8所示的两个相互平行的线性子阵列组成，两个子阵列分别为子阵X。和子阵七。每个子阵列由M个传感器，两相邻传感器沿X方向的间距为刁，两个子阵列的间距也为刁。假设空间有N个非相关的远场窄带信号s(t)(1<n<N)入射到此阵列，其与X轴的夹角为0，与Y轴的夹角为甲。子阵列X和Y在t时刻对应输出信号分别为子阵,■图8线性平行阵列结构图TOC\o"1-5"\h\zx(t)=As(t)+n(t)(33)y(t)=A职(t)+n(t)(34)式中，％(t)，n,(t)为两个子阵列的加性高斯白噪声向量，与信号s(t)相互独立;x(t)=「X(t),...,X(t汁，y(t)=「,(t),...,,(t)V，s(t)=「s(t),...,s(t)V，1-1MA1-1MA1-1N」A=、,...,an],其中a=1,ej2Wcos(0n)/^，…，ej2兀(m-1)dcos(0n)口T,n=1,…，K(35)鱼=diagCj2兀d*&])/入，…，ej2膈cos(^N)/入｝(36)式中，人为波长。4.2DOA矩阵方法设接收到的数据x(t)的自相关矩阵为R了其表达式为TOC\o"1-5"\h\zR=E[x(t)xh(t牛APAh+b21m(37)式中，P=E「s(t)sh(t)]为信号源的协方差矩阵，I为MxM的单位矩阵，6为1-JM加性高斯白噪声的方差。考虑到噪声自身的独立性，且独立于信号，假设y(t)和X(t)的互相关矩阵为R弟，则R=E「y(t)xh(t)]=A鱼PAh(38)对R进行特征值分解(EVD)，令七…点n表示矩阵R的N个大特征值，在白噪声的假设下，可以由(M-K)个小特征值的平均得到噪声方差6的估计值。然后通过去除噪声的影响，可以得到C=APAh=R-a21m(39)根据DOA矩阵算法的思想，可以定义DOA矩阵为R=RC+(40)定理1如果A与P满秩，①无相同对角元素，则DOA矩阵的N个非零特征值等于①中的N个对角元素，而这些值对应的特征向量等于相应的信号方向矢量，即RA=A鱼(41)因此只要对DOA矩阵R进行EVD，就可以得到矩阵A和鱼，进而通过式(35)和式(36)得到DOA的估计，这种方法成为DOA矩阵方法。现考虑20阵元数的双平行均匀线阵，有三个信源入射，它们的DOA分别是(10o,15o),(25o,30o)和(30o,40o)，信噪比SNR=20dB，快拍数2=500。DOA矩阵算法的估计效果如图9所示。图9DOA矩阵算法实验一：信噪比SNR=20dB，快拍数从50到500变化，变化步长50。运行1000次蒙特卡洛（MC）试验，算法仿真