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

摘要

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

说明书

技术领域

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

背景技术

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

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

发明内容

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

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

1）采用天线接收机形成均匀平面阵，其中x轴方向有N个天线接收机，y 轴方向有M个天线接收机，且天线接收机间距均为d，N≥2，M≥2， 0<d≤λ/2，λ为入射窄带信号波长；

2）将Q个不相关的窄带信号s_q(t)以x轴夹角α_q和y轴夹角β_q入射到天线 均匀平面阵上，其中，q=1,2,...,Q，1≤Q<NM，α_q表示待求的第q个窄带信号 的入射方向与x轴的夹角，β_q表示待求的第q个窄带信号的入射方向与y轴的 夹角，且α_q∈(0°,180°]，β_q∈(0°,180°]；

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

$Y_{n,m}\left(t\right)=\underset{q=1}{\overset{Q}{\Sigma }}[a_n\left({\alpha }_q\right)\otimes a_m\left({\beta }_q\right)]·s_q\left(t\right)+N\left(t\right)$

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

a_n(α_q)=(1,exp(j2πd·cosα_q/λ),...,exp[j2πd·(N-1)·cosα_q/λ])^T

a_m(β_q)=(1,exp(j2πd·cosβ_q/λ),...,exp[j2πd·(M-1)·cosβ_q/λ])^T

其中，j表示虚数单位，(·)^T表示矩阵转置运算；

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

5）构造角度求解函数为：$\min \left(a_m^H\left({\beta }_q\right)G\left({\alpha }_q\right)a_m\left({\beta }_q\right)\right),$

其中，$G\left({\alpha }_q\right)={[a_n\left({\alpha }_q\right)\otimes I_M]}^H{U}_n{U}_n^H[a_n\left({\alpha }_q\right)\otimes I_M]$为中间变量，min(·)表示 求最小值运算，(·)^H表示矩阵共轭转置运算，I_M表示M×M维的单位矩阵；

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

a_m(β_q)=P_min(G_q)，

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

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

8）将所取的多个x轴夹角α_q和得到的多个方向向量a_m(β_q)代入函数：

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

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

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

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

实验结果表明，本发明计算二维波达方向角度角度的运算量仅为

Ο{LM²N²+M³N³+n[(M²N+M²)(MN-Q)+M³+M²]}，而传统2D-MUSIC 方法计算二维波达方向角度角度的运算量为

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

附图说明

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

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

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

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

具体实施方式

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

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

在x轴方向每隔间距d放置1个天线接收机，共放置N个，将其称为1行 天线接收机；在y轴方向每隔间距d放置1行天线接收机，共放置M行，形成 天线接收机总个数为NM的天线均匀平面阵，其中N≥2，M≥2，0<d≤λ/2， λ为入射窄带信号波长。

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

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

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

B·ΔT_max<<1，

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

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

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

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

$Y_{n,m}\left(t\right)=\underset{q=1}{\overset{Q}{\Sigma }}[a_n\left({\alpha }_q\right)\otimes a_m\left({\beta }_q\right)]·s_q\left(t\right)+N\left(t\right)$

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

a_n(α_q)表示第q个入射的窄带信号的x轴夹角α_q的方向向量，a_m(β_q)表示 第q个入射的窄带信号的y轴夹角β_q的方向向量，展开形式分别表示为：

a_n(α_q)=(1,exp(j2πd·cosα_q/λ),...,exp[j2πd·(N-1)·cosα_q/λ])^T

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

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

步骤4：计算天线接收机回波信号Y_n,m(t)的噪声子空间U_n。

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

4a)计算天线接收机信号的协方差矩阵：R=E[Y_n,m(t)·Y_n,m(t)^H]，

其中E[·]表示求数学期望，(·)^H表示矩阵共轭转置运算；

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

R=U·Λ·U^H，

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

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

步骤5：构造角度求解函数$\min \left(a_m^H\left({\beta }_q\right)G\left({\alpha }_q\right)a_m\left({\beta }_q\right)\right).$

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

$V({\alpha }_q,{\beta }_q)={[a_n\left({\alpha }_q\right)\otimes a_m\left({\beta }_q\right)]}^H{U}_n{U}_n^H[a_n\left({\alpha }_q\right)\otimes a_m\left({\beta }_q\right)],$

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

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

$V({\alpha }_q,{\beta }_q)=a_m^H\left({\beta }_q\right){[a_n\left({\alpha }_q\right)\otimes I_M]}^H{U}_n{U}_n^H[a_n\left({\alpha }_q\right)\otimes I_M]a_m\left({\beta }_q\right),$

其中，I_M表示M×M维的单位矩阵；

5c)令中间变量$G\left({\alpha }_q\right)={[a_n\left({\alpha }_q\right)\otimes I_M]}^H{U}_n{U}_n^H[a_n\left({\alpha }_q\right)\otimes I_M],$将上述步骤5b)

中的方程式简化为：

$V({\alpha }_q,{\beta }_q)=a_m^H\left({\beta }_q\right)G\left({\alpha }_q\right)a_m\left({\beta }_q\right);$

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

$\min \left(a_m^H\left({\beta }_q\right)G\left({\alpha }_q\right)a_m\left({\beta }_q\right)\right).$

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

6a）将x轴夹角α_q代入式：$G\left({\alpha }_q\right)={[a_n\left({\alpha }_q\right)\otimes I_M]}^H{U}_n{U}_n^H[a_n\left({\alpha }_q\right)\otimes I_M]$中， 得到中间变量矩阵G_q；

6b）对中间变量矩阵G_q进行特征分解：

$G_q=U_q·{\Lambda }_q·U_q^H$

其中，Λ_q为中间变量矩阵G_q的特征值矩阵，U_q为中间变量矩阵G_q的特征 值所对应的特征向量矩阵；

6c）寻找中间变量矩阵G_q中最小的特征值Λ_qmin对应的特征向量U_qmin， 将该特征向量U_qmin作为求得的方向向量a_m(β_q)。

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

步骤8：绘制幅度谱图。

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

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

具体步骤为：

8a）将所取的多个角度α_q和得到的多个方向向量a_m(β_q)代入函数： $Z_q({\alpha }_q,a_m\left({\beta }_q\right))=1/\left[a_m^H\left({\beta }_q\right)G\left({\alpha }_q\right)a_m\left({\beta }_q\right)\right]$中，获得相应的多个函数值Z_q；

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

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

步骤9：完成x轴夹角α_q的估计。

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

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

求解y轴夹角β_q的步骤为：

10a)根据方向向量a_m(β_q)利用下式求角度向量g：

g=-angle(a_m(β_q))，

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

10b)令矩阵$P=\left(\begin{array}{cc}1& 0\\ 1& 1\\ .& .\\ .& .\\ .& .\\ 1& M-1\end{array}\right),$计算下式：

[c₀,c₁]=(P^TP)^-1P^Tg，

其中，c₀和c₁表示所求向量的两个元素，(·)^-1表示求逆矩阵运算；

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

${\beta }_q=\arccos \left(\frac{c_1\lambda }{2\mathrm{\pi d}}\right).$

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

1.仿真条件与方法：

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

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

$\left(\begin{array}{c}{\mathrm{RMSE}}_{\alpha }=\frac{1}{Q}\underset{q=1}{\overset{Q}{\Sigma }}\sqrt{\frac{1}{\mathrm{Num}}\underset{\mathrm{num}=1}{\overset{\mathrm{Mum}}{\Sigma }}{({\hat{\alpha }}_q-{\alpha }_q)}^2}\\ {\mathrm{RMSE}}_{\beta }=\frac{1}{Q}\underset{q=1}{\overset{Q}{\Sigma }}\sqrt{\frac{1}{\mathrm{Mum}}\underset{\mathrm{num}=1}{\overset{\mathrm{Num}}{\Sigma }}{({\hat{\beta }}_q-{\beta }_q)}^2}\end{array}\right),$

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

2.仿真内容与结果

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

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

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

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

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