电磁信号的二维波达方向角估计方法

摘要

本发明公开了一种电磁信号的二维波达方向角估计方法，主要解决现有技术中二维方向角估计运算量大，造成目标侦察和无源定位反应速度慢和估计误差大的问题。其实现步骤为：采用天线接收机形成均匀平面阵；计算所有天线接收机回波信号；计算天线接收机回波信号的噪声子空间；令x轴夹角为固定值代入角度求解函数获得y轴夹角的方向向量，多次运算得到多组值；利用多组值生成幅度谱图；通过寻找幅度谱图中的峰值点获得x轴夹角；通过最小二乘法获得y轴夹角。本发明将二维波达方向角估计化简成一维波达方向角估计过程，大幅度降低了运算量，提高了目标侦察和无源定位反应速度，避免了因信息滞后引起的参数估计误差，可用于快速目标侦察和无源定位。

权利要求

1.一种电磁信号的二维波达方向角估计方法，包括以下步骤：

1)采用天线接收机形成均匀平面阵，其中x轴方向有N个天线接收机，y轴方向有M个天线接收机，且天线接收机间距均为d，N≥2，M≥2，0

2)将Q个不相关的窄带信号sq(t)以x轴夹角αq和y轴夹角βq入射到天线均匀平面阵上，其中，q＝1,2,...,Q，1≤Q

3)计算天线均匀平面阵上各天线接收机的回波信号Yn,m(t)；

其中， 表示矩阵叉乘运算，N(t)为噪声信号向量，an(αq)和am(βq)分别表示第q个入射的窄带信号的x轴夹角αq和y轴夹角βq相应的方向向量，展开形式表示为：T

an(αq)＝(1,exp(j2πd·cosαq/λ),...,exp[j2πd·(N-1)·cosαq/λ])Tam(βq)＝(1,exp(j2πd·cosβq/λ),...,exp[j2πd·(M-1)·cosβq/λ])T其中，j表示虚数单位，(·)表示矩阵转置运算；

4)根据天线接收机回波信号Yn,m(t)，计算噪声子空间Un；

5)构造角度求解函数为：

其中， 为中间变量，min(·)表示求最小值运H

算，(·)表示矩阵共轭转置运算，IM表示M×M维的单位矩阵；

6)设x轴夹角αq为(0°,180°]范围内的一个固定值，将x轴夹角αq代入中间变量函数G(αq)得到中间变量矩阵Gq，然后利用下式计算方向向量am(βq)：am(βq)＝Pmin(Gq)，

其中，Pmin(·)表示求解矩阵的最小特征值所对应的特征向量；

7)取x轴夹角αq为(0°,180°]范围内的其它固定值，重复执行步骤6)，获得相应的方向向量am(βq)；

8)将所取的多个x轴夹角αq和得到的多个方向向量am(βq)代入函数：中，获得相应的多个函数值Zq，再在平面中以(0,0)为原点坐标，以多组(αq,Zq)为x，y坐标绘制幅度值点，将各点相连得到幅度谱图；

9)从幅度谱图中按照从高到低的顺序寻找幅值较大的前Q个谱峰，将这些谱峰的峰值点对应的x轴坐标值作为求得的x轴夹角αq值；

10)根据每个求得的x轴夹角αq值，寻找对应的方向向量am(βq)值，再利用最小二乘法求得相应的y轴夹角βq值，完成电磁信号的二维波达方向角估计。

2.根据权利要求1所述的电磁信号的二维波达方向角估计方法，其中，步骤4)所述的计算天线接收机回波信号Yn,m(t)的噪声子空间Un，按如下步骤进行：H

4a)计算天线接收机信号的协方差矩阵R＝E[Yn,m(t)·Yn,m(t)]，其中E[·]表示求数学期望；

4b)对协方差矩阵R进行如下特征分解：

H

R＝U·Λ·U，

其中，Λ为协方差矩阵R的特征值矩阵，U为特征值所对应的特征向量矩阵；

4c)将特征值矩阵Λ中的特征值按从大到小排序，取后(NM-Q)个较小特征值对应的特征向量矩阵作为噪声子空间Un。

3.根据权利要求1所述的电磁信号的二维波达方向角估计方法，其中所述步骤6)中计算中间变量矩阵Gq的最小特征值所对应的特征向量，求得方向向量am(βq)，按如下步骤进行：

6a)对中间变量矩阵Gq进行特征分解：

其中，Λq为中间变量矩阵Gq的特征值矩阵，Uq为中间变量矩阵Gq的特征值所对应的特征向量矩阵；

6b)寻找中间变量矩阵Gq中最小的特征值Λqmin对应的特征向量Uqmin，将该特征向量Uqmin作为求得的方向向量am(βq)。

4.根据权利要求1所述的电磁信号的二维波达方向角估计方法，其中步骤10)所述的根据方向向量am(βq)，利用最小二乘法求得y轴夹角βq值，按如下步骤进行：

10a)根据方向向量am(βq)，利用下式计算角度向量g：g＝-angle(am(βq))

其中，angle(·)表示对向量中每个元素分别取相角；

10b)令矩阵 计算下式：

T T -1 T

[c0,c1]＝(P P) Pg，

-1 T

其中，c0和c1分别表示所求向量的两个元素，(·) 表示求逆矩阵运算，(·) 表示矩阵转置运算；

10c)求得y轴夹角βq值：

说明书

电磁信号的二维波达方向角估计方法

技术领域

[0001] 本发明属于电磁信号处理技术领域，特别涉及一种二维波达方向角估计方法，可用于目标侦察与无源定位。

背景技术

[0002] 波达方向角DOA估计是利用处于空间不同位置的天线信号阵列接收多个不同方向的信号源发出的电磁信号，运用现代信号处理方法快速准确的估计出信号源的方向，在雷达、声纳、无线通信等领域具有重要应用价值。二维波达方向角估计一般采用面阵或矢量传感器实现二维参数的估计，多数有效的二维波达方向角估计算法属于一维波达方向角估计算法的直接扩展，并没有充分利用阵列接收信号中携带的多维信息，往往存在运算量过大或存在角度配对等缺点。

[0003] 传统的二维波达方向角估计方法包括两种：一种是利用二维角度搜索同时获得二维角度；另一种是分别求取其中一维角度值，再进行角度配对。第一种方法需要使用二维搜索，其运算量是一维搜索的平方倍，因此运算量巨大；第二种方法在分别完成一维角度估计后还需要使用角度配对算法，该过程仍需要花费一定运算量，因此这两种方法均具有较大运算量。实际应用中，目标侦察与无源定位均需要在角度估计的基础上进行，如果角度估计运算量大将造成目标侦察和无源定位反应速度慢，甚至造成由于信息滞后引起的参数估计误差较大。

发明内容

[0004] 本发明的目的在于针对上述已有技术的不足，提出一种电磁信号的二维波达方向角估计方法，以大幅度减小估计的运算量，提高目标侦察和无源定位反应速度，避免因信息滞后引起的参数估计误差。

[0005] 为实现上述目的，本发明的实现步骤包括如下：

[0006] 1）采用天线接收机形成均匀平面阵，其中x轴方向有N个天线接收机，y轴方向有M个天线接收机，且天线接收机间距均为d，N≥2，M≥2，0

[0007] 2）将Q个不相关的窄带信号sq(t)以x轴夹角αq和y轴夹角βq入射到天线均匀平面阵上，其中，q=1,2,...,Q，1≤Q

[0008] 3）计算天线均匀平面阵上各天线接收机的回波信号Yn,m(t)；

[0009]

[0010] 其中， 表示矩阵叉乘运算，N(t)为噪声信号向量，an(αq)和am(βq)分别表示第q个入射的窄带信号的x轴夹角αq和y轴夹角βq相应的方向向量，展开形式表示为：

[0011] an(αq)=(1,exp(j2πd·cosαq/λ),...,exp[j2πd·(N-1)·cosαq/λ])T[0012] am(βq)=(1,exp(j2πd·cosβq/λ),...,exp[j2πd·(M-1)·cosβq/λ])T[0013] 其中，j表示虚数单位，(·)T表示矩阵转置运算；

[0014] 4）根据天线接收机回波信号Yn,m(t)，计算噪声子空间Un；

[0015] 5）构造角度求解函数为：

[0016] 其中， 为中间变量，min(·)表示求最H
小值运算，(·)表示矩阵共轭转置运算，IM表示M×M维的单位矩阵；

[0017] 6）设x轴夹角αq为(0°,180°]范围内的一个固定值，将x轴夹角αq代入中间变量函数G(αq)得到中间变量矩阵Gq，然后利用下式计算方向向量am(βq)：

[0018] am(βq)=Pmin(Gq)，

[0019] 其中，Pmin(·)表示求解矩阵的最小特征值所对应的特征向量；

[0020] 7）取x轴夹角αq为(0°,180°]范围内的其它固定值，重复执行步骤6），获得相应的方向向量am(βq)；

[0021] 8）将所取的多个x轴夹角αq和得到的多个方向向量am(βq)代入函数：

[0022] 中，获得相应的多个函数值Zq，再在平面中以(0,0)为原点坐标，以多组(αq,Zq)为x，y坐标绘制幅度值点，将各点相连得到幅度谱图；

[0023] 9）从幅度谱图中按照从高到低的顺序寻找幅值较大的前Q个谱峰，将这些谱峰的峰值点对应的x轴坐标值作为求得的x轴夹角αq值；

[0024] 10）根据每个求得的x轴夹角αq值，寻找对应的方向向量am(βq)值，再利用最小二乘法求得相应的y轴夹角βq值，完成电磁信号的二维波达方向角估计。

[0025] 本发明由于将二维波达方向角度估计过程化简成一维波达方向角度估计过程，只需要一维角度搜索即可同时获得二维波达方向角度，避免了二维角度搜索或角度配对过程，大大降低了运算复杂度。

[0026] 实验结果表明，本发明计算二维波达方向角度角度的运算量仅为[0027] Ο{LM2N2+M3N3+n[(M2N+M2)(MN-Q)+M3+M2]}，而传统2D-MUSIC方法计算二维波达方向角度角度的运算量为

[0028] Ο{LM2N2+M3N3+n2[MN(MN-Q)]}，其中求解αq与βq时的搜索点数均为n，n≥1，由于x轴方向接收机数量N、y轴方向接收机数量M、窄带信号数量Q和快拍数L均比n要小的多，因此本发明运算量的量级为n，而2D-MUSIC方法运算量的量级为n2，可见本发明方法的运算量远低于传统2D-MUSIC方法，能够提高目标侦察与无源定位的反应速度和参数估计准确性，避免了因信息滞后引起的参数估计误差。

附图说明

[0029] 图1是本发明的实现流程图；

[0030] 图2是本发明与现有2D-MUSIC算法的运算量对比图；

[0031] 图3是本发明与现有2D-MUSIC算法的x轴夹角估计均方根误差对比图；

[0032] 图4是本发明与现有2D-MUSIC算法的y轴夹角估计均方根误差对比图。

具体实施方式

[0033] 参照图1，本发明的电磁信号的二维波达方向角度估计，其实现步骤如下：

[0034] 步骤1：利用天线接收机形成均匀平面阵。

[0035] 在x轴方向每隔间距d放置1个天线接收机，共放置N个，将其称为1行天线接收机；在y轴方向每隔间距d放置1行天线接收机，共放置M行，形成天线接收机总个数为NM的天线均匀平面阵，其中N≥2，M≥2，0

[0036] 步骤2：将窄带信号入射到天线均匀平面阵上。

[0037] 设入射到天线均匀平面阵上的窄带信号为Q个，且信号之间不相关，信号形式用sq(t)表示；

[0038] 所述的窄带信号满足如下条件：

[0039] B·ΔTmax<<1，

[0040] 其中B表示信号带宽，ΔTmax表示信号到达任意两个天线接收机的时延差的最大值；

[0041] 所述窄带信号从空间任意方向射向天线均匀平面阵，该入射方向可被看作一条射线，该射线与天线均匀平面阵所在平面形成的角度称为二维波达方向角，用x轴夹角αq和y轴夹角βq表示，其中，q=1,2,...,Q，1≤Q

[0042] 本实例中要求窄带信号个数Q是已知的值，αq和βq是待估计的值。

[0043] 步骤3：计算天线均匀平面阵上各天线接收机的回波信号Yn,m(t)：

[0044]

[0045] 其中， 表示矩阵叉乘运算，N(t)为噪声信号向量，该向量是一个随机向量，在各天线接收机的回波信号Yn,m(t)中各不相同，信噪比越小，对波达方向角度的估计误差越大；

[0046] an(αq)表示第q个入射的窄带信号的x轴夹角αq的方向向量，am(βq)表示第q个入射的窄带信号的y轴夹角βq的方向向量，展开形式分别表示为：T

[0047] an(αq)=(1,exp(j2πd·cosαq/λ),...,exp[j2πd·(N-1)·cosαq/λ])T

[0048] am(βq)=(1,exp(j2πd·cosβq/λ),...,exp[j2πd·(M-1)·cosβq/λ])T

[0049] 其中，j表示虚数单位，(·)表示矩阵转置运算。

[0050] 步骤4：计算天线接收机回波信号Yn,m(t)的噪声子空间Un。

[0051] 现有技术中，获得天线接收机回波信号噪声子空间Un的方法有很多，例如特征分解法，奇异值分解法，多级维纳滤波法，双迭代法，幂迭代法等，本实例中使用了应用较为普遍的特征分解法，具体方法如下：H

[0052] 4a)计算天线接收机信号的协方差矩阵：R=E[Yn,m(t)·Yn,m(t)]，[0053] 其中E[·]表示求数学期望，(·)H表示矩阵共轭转置运算；

[0054] 4b)对协方差矩阵R进行特征分解，即：

[0055] R=U·Λ·UH，

[0056] 其中，Λ为协方差矩阵R的特征值矩阵，该矩阵的对角元素为特征值，且除对角元素以外的元素全为0，U为特征值所对应的特征向量矩阵；

[0057] 4c)将特征值矩阵Λ中的特征值按从大到小排序，取后(NM-Q)个较小特征值对应的特征向量矩阵作为噪声子空间Un。

[0058] 步骤5：构造角度求解函数

[0059] 5a)根据多信号分类特征分解MUSIC算法，构造出角度求解函数为：其中，函数V(αq,βq)为：

[0060]

[0061] 其中，min(·)表示求最小值运算；

[0062] 5b)将上述步骤5a)中的方程式变形为：

[0063]

[0064] 其中，IM表示M×M维的单位矩阵；

[0065] 5c)令中间变量 将上述步骤5b)

[0066] 中的方程式简化为：

[0067]

[0068] 5d）根据步骤5c)中的方程式V(αq,βq)，最终构造角度求解函数为：

[0069]

[0070] 步骤6：设x轴夹角αq为(0°,180°]范围内的一个固定值，计算方向向量am(βq)。

[0071] 6a）将x轴夹角αq代入式： 中，得到中间变量矩阵Gq；

[0072] 6b）对中间变量矩阵Gq进行特征分解：

[0073]

[0074] 其中，Λq为中间变量矩阵Gq的特征值矩阵，Uq为中间变量矩阵Gq的特征值所对应的特征向量矩阵；

[0075] 6c）寻找中间变量矩阵Gq中最小的特征值Λqmin对应的特征向量Uqmin，将该特征向量Uqmin作为求得的方向向量am(βq)。

[0076] 步骤7：取x轴夹角αq为(0°,180°]范围内的其它固定值，重复执行步骤6，获得相应的am(βq)值；αq在取值范围内按照从小到大的顺序，以固定角度间隔进行取值；取值间隔根据期望达到的角度估计精度进行设定，取值间隔越小，角度估计精度越高。

[0077] 步骤8：绘制幅度谱图。

[0078] 上述所取的多个x轴夹角αq和得到的多个方向向量am(βq)中，只有满足角度求解函数 的Q组αq和βq角度值才是所求的二维角度值，具体实现方法是将所取的多个x轴夹角αq和得到的多个方向向量am(βq)分别代入函数中，寻找极小值点对应的αq和βq角度值。

[0079] 人们一般习惯在幅度谱图中寻找极大值点，而不是极小值点，为了满足人们的观测习惯，将所取的多个x轴夹角αq和得到的多个方向向量am(βq)分别代入函数中，用寻找极大值点对应的αq和βq角度值的方法代替上述方法。

[0080] 具体步骤为：

[0081] 8a）将所取的多个角度αq和得到的多个方向向量am(βq)代入函数：中，获得相应的多个函数值Zq；

[0082] 8b）在平面中以(0,0)为原点坐标，以多组(αq,Zq)为x，y坐标，绘制幅度值点；

[0083] 8c）将各个幅度值点相连，得到幅度谱图。

[0084] 步骤9：完成x轴夹角αq的估计。

[0085] 从幅度谱图中按照从高到低的顺序寻找幅值较大的前Q个谱峰，本实例取Q=3；将这些谱峰的峰值点对应的x坐标值作为求得的x轴夹角αq值。

[0086] 步骤10：根据每个求得的x轴夹角αq值，寻找对应的方向向量am(βq)值，再利用最小二乘法求得相应的y轴夹角βq值，利用上述方法求解y轴夹角βq的次数为Q次，本实例取Q=3。

[0087] 求解y轴夹角βq的步骤为：

[0088] 10a)根据方向向量am(βq)利用下式求角度向量g：

[0089] g=-angle(am(βq))，

[0090] 其中，angle(·)表示对向量中每个元素取相角；

[0091] 10b)令矩阵 计算下式：

[0092] [c0,c1]=(PTP)-1PTg，

[0093] 其中，c0和c1表示所求向量的两个元素，(·)-1表示求逆矩阵运算；

[0094] 10c)求得y轴夹角βq值：

[0095]

[0096] 本发明的效果可通过以下仿真说明：

[0097] 1.仿真条件与方法：

[0098] 采用天线接收机形成均匀平面阵，其中x轴方向有8个天线接收机，y轴方向有8个天线接收机，天线接收机间距d等于λ/2。有三个窄带信号入射到均匀平面阵，其中窄带信号入射方向与x轴的夹角αq和窄带信号入射方向与y轴的夹角βq的取值分别为(60°,50°)、(90°,100°)和(145°,75°)。信噪比SNR为0dB，快拍数L等于64，角度搜索取值间隔为0.01度。

[0099] 为了进一步评价本发明的性能，对多次独立的实验结果进行了平均，并采用二维角度的均方根误差作为评价指标：

[0100]

[0101] 其中RMSEα和RMSEβ分别表示αq和βq的估计均方根误差，Num为实验次数，αq为第q个角度的真实值， 表示第q个角度的估计值。本仿真分别采用本发明和现有的2D-MUSIC算法进行。

[0102] 2.仿真内容与结果

[0103] 仿真1，利用本发明和现有的2D-MUSIC算法分别进行二维波达方向角度估计，并统计运算量，其结果如图2所示。图2中横坐标表示角度搜索次数，纵坐标表示运算量。

[0104] 从图2可以看出，本发明与现有2D-MUSIC算法相比，大幅度降低了角度估计的运算量，且随着角度搜索次数的增加，运算量降低幅度越大。

[0105] 仿真2，利用本发明和现有的2D-MUSIC算法分别进行100次独立的二维波达方向角度估计实验，分别计算x轴夹角的均方根误差和y轴夹角的均方根误差，其结果如图3和图4所示。其中图3为计算x轴夹角的均方根误差，图4为计算y轴夹角的均方根误差；图3和图4中横坐标表示信噪比，纵坐标表示均方根误差，信噪比变化范围为-4dB～8dB。

[0106] 从图3和图4可以看出，本发明与现有2D-MUSIC算法角度估计精度基本相当。

[0107] 综上，本发明在保证角度估计精度的同时，大幅度减小了运算量，保证了目标侦察和无源定位的快速反应，避免了因信息滞后引起的参数估计误差。