CN113419209B - 一种锥面共形阵列盲极化波达方向估计方法 - Google Patents

Abstract

本发明提供一种锥面共形阵列盲极化波达方向估计方法，包括下列步骤：步骤一：对锥面共形阵列天线快拍进行数据建模；步骤二：实值盲极化波达方向估计方法。新算法将子空间估计过程实值化处理，通过子阵分割设计解决了协方差矩阵分裂引起的“镜像模糊”问题，并推导了参数估计的克拉美—罗界，所提算法在保证DOA估计精度的前提下降低了算法复杂度。

Description

一种锥面共形阵列盲极化波达方向估计方法

技术领域

本发明涉及信号与信息处理技术，具体涉及一种锥面共形阵列盲极化波达方向估计方法。

背景技术

在现代无线通信系统中，共形阵列天线由于能够与飞机、导弹以及卫星等高速运行的载体平台表面相共形，且并不破坏载体的外形结构及空气动力学等特性，成为天线领域的一个研究热点。其中，来波信号极化特征与方位信息的耦合是共形阵列天线估计的难点，共形阵列波达方向(direction of arrival,DOA)估计大都采用旋转不变子空间算法(Estimation of Signal Parameters via Rotational Invariance Techniques,ESPRIT)，将来波信号的二维角参数与极化参数解耦合，以达到不考虑极化参数，实现高分辨方位估计的效果，但此类算法的子空间估计涉及复值矩阵特的征值分解，计算量庞大的问题较为突出。高效超分辨算法致力于减少系统的复杂度和算法的运算量，以满足应用系统对多信源实时测向的需求，实值/半实值类算法就是其中的典型代表，实现超分辨算法实值化最具代表性的方法是基于阵列接收数据协方差阵的酉变换技术，文献充分利用中心对称阵列输出数据协方差矩阵为艾米特中心对称矩阵的特性，通过数据变换，将复数运算转换为实值运算，可将运算量缩减了75％。另外，协方差矩阵分裂估计对实值运算和阵列结构的任意性进行折衷，进而发展了两种全新的基于半实值运算的超分辨算法，摆脱了算法对阵列结构的依赖性，但还存在镜像方位的模糊问题，需要在算法应用过程中进一步处理。

发明内容

本发明的目的在于克服上述现有技术中的不足之处，提出一种锥面共形阵列盲极化波达方向估计方法。

本发明是通过如下方式实现的：

一种锥面共形阵列盲极化波达方向估计方法，包括下列步骤：

步骤一：对锥面共形阵列天线快拍进行数据建模；

步骤二：实值盲极化波达方向估计方法。

进一步，所述的步骤一具体包括下述步骤：

Step 1)利用锥面共形阵列的对称性，将锥面共形阵列划分为多个子阵，每一子阵的波达方向估计算法均相同；在锥面共形阵列的圆锥面的母线上划分l₁，l₂，l₃三个子阵，每个子阵上有m个阵元均匀分布，并且处于同一条母线上的相邻阵元间距为λ/2，其中，λ为入射信号波长；子阵l₁，l₂，l₃在全局坐标系XOY面上进行投影，三条母线与X轴顺时针的旋转方向的夹角分别为α₁，α₂，α₃，且α₁＝π-α₃，α₂∈(α₁，α₃)，以上子阵设置可实现方位角

的估计，将整个圆锥上的多个子阵综合起来，实现/>

的方位角估计；

Step 2)由共形阵列天线快拍数据建模方法，得锥面共形阵列的导向矢量为：

式(1)中，r_i为第i个阵元在全局坐标系中对单位强度来波信号的响应，G_i为第i个阵元在全局坐标系中的位置矢量，u为来波信号的方向矢量，θ和φ分别为来波信号相对于全局坐标系的俯仰角和方位角，j为虚数；

当有n个来波信号时，阵列的流形矩阵为：

A＝[a(θ₁，φ₁)，a(θ₂，φ₂)，…，a(θ_n，φ_n)] (2)

其中，此锥面共形阵列的快拍数据模型为：

X＝AS+N＝(A_θK_θ+A_φK_φ)S+N (3)

K_θ＝diag(k_1θ，k_2θ，…，k_nθ) (4)

S＝[s₁，s₂，…，s_n]^T (6)

N＝[n₁，n₂，…，n_n]^T (7)

式(3)-(7)中，K为对角阵，S为信号矢量，A_θ、A_φ分别为来波信号在极化正交基矢量上的流形矩阵，K_θ、K_φ分别为极化正交基矢量上投影分量构成的对角阵，K_nθ、K_nφ分别为第n个极化正交量上投影分量构成的对角阵，s_n为第n个辐射源，n_n为第n个加性噪声，N为噪声矢量，θ和

分别为入射信号在全局坐标系中的俯仰角和方位角。

进一步，所述的步骤二具体包括下述步骤：

Step 1)在已求出快拍数据矩阵X的基础上，按照式(8)得到此阵列的数据协方差矩阵R；

可将R-AOCM和I-AOCM表示为：

式(8)-(10)中，N为噪声矢量，X^H为快拍数据矩阵的共轭转置，Re(X)为快拍数据矩阵的实部，Im(X)为快拍数据矩阵的虚部，Im^T(X)为快拍数据矩阵的虚部转置，Re^T(X)为快拍数据矩阵的实部的转置；

Step 2)采用实值化算法，在R-ACOM上进行特征分解，得到实值的信号子空间，从而分别得到母线l₁、l₂、l₃对应的Φ₁、Φ₂、Φ₃，其中，Φ₁、Φ₂、Φ₃分别为三条母线上阵元构成子阵的相位差；

设原始的信号子空间span(S)及其共轭信号子空间span(S^*)的交集，假设此交集为span(G)，表示为：

式(11)中，G为信号子空间与共轭信号子空间交集的基矢量，S为信号子空间基矢量，S^*为信号子空间基矢量的共轭；

在对信号协方差矩阵R的实部或者虚部直接进行特征分解的处理后，提取信号子空间时，需要提取的特征向量应该翻倍，每个子阵对均由同一条母线上分布的m个阵元构造，即1～m-1阵元，构成第一个阵列，2～m阵元，构成另一个阵列，两个阵列之间的距离为λ/2.以此类推，可分别得到母线l₁、l₂、l₃上距离矢量各不相同的三个子阵对；对母线l₁上分布的阵元进行分割，子阵l₁₁的接收数据为：

子阵l₁₂的接收数据为：

其中，Φ₁为相位差矩阵：

Φ₁＝diag[exp(-jw₁₁)，…，exp(-jw_1n)] (14)

w_1i＝(2π/λ)dv₁·u_i＝πv₁·u_i (15)

展开即得：

X₁₁和X₁₂的最后一行组成了母线l₁的快拍数据X₁：

X₁＝A₁S+N₁＝[X₁₁；X₁₂(m-1，：)] (17)

式(12)-(17)中，N₁₁表示子阵l₁₁接收的噪声矢量，A₁为母线1上阵元构成阵列的流形矩阵，A₁₁，A₁₂分别为母线1上阵元构成的子阵1和子阵2的流形矩阵，S表示信号矢量，A_11θ、

分别为母线1上阵元构成的子阵1在极化正交基矢量上的流形矩阵，K_θ、/>

分别为极化正交基矢量上投影分量构成的对角阵，Φ₁为子阵相位差构成的对角阵，j表示虚部，w_1i为母线1上子阵对对第i个信号的响应相差，u_i为第i个信号方向矢量，θ_v1、/>

为母线1子阵对对应阵元距离矢量的表征参数，θ_i、/>

分别表示第i个来波的俯仰角与方位角，N₁为母线1上阵列的加性噪声矩阵，

N₁₂表示子阵l₁₂接收的噪声矢量，v₁为子阵l₁₁与子阵l₁₂距离的方向矢量，同理，将l₂、l₃划分为l₂₁、l₂₂和l₃₁、l₃₂，即可获得其快拍数据X₂和X₃以及Φ₂和Φ₃；

Step 3)对Φ₁、Φ₂、Φ₃进行特征分解，得到t_1i、t_2i、t_3i；

通过Φ₁，Φ₂，Φ₃这几个包含旋转不变关系的矩阵，结合式(15)求解信号的方向参数；

阵列的数据协方差矩阵为：

利用式(9)或(10)替换协方差矩阵，完成矩阵实值化处理，并进行特征分解：

U_S为阵列信号子空间的特征矢量矩阵，根据式(17)，U_S的第1～m-1行为子阵l₁₁对应的信号子空间，将其记为U_S11；U_S的第2～m行为子阵l₁₂对应的信号子空间，将其记为U_S1，同理得到U_S21，U_S22，U_S31，U_S32；

式(18)-(19)中，U_N为噪声子空间，

为信号子空间的共轭转置；/>

为噪声子空间的共轭转置；

对应ESPRIT里的最小二乘法：

式(20)中，

为母线1上子阵1的信号子空间的共轭转置，U_S1为母线1上子阵1的信号子空间；

且式(14)中的exp(-jw_1i)(i＝1，2，…，n)即为Φ₁的特征值t_1i，同理可得Φ₁的特征值t_1i，Φ₂的特征值t_2i，Φ₃的特征值t_3i；

Step 4)根据式(21)-式(23)，得到对应的N个接收信号的方向信息；

结合α₁＝π-α₃，α₂∈(α₁，α₃)分别求解t_1i，t_2i，t_3i，就可求出入射信号的方位角

入射信号的俯仰角θ_i：

z_ji＝angle(t_ji)/(2πd)，j＝1，2，3 (23)

式(21)-(23)中，z_ji为母线j上阵列对第i个信号的接收相位；

Step 5)通过限定子阵的空域覆盖范围，对估计出的信号俯仰角和方位角进行判定，排除虚假信号，从而得到最终准确的信号方位。

本发明的有益效果在于：提供了一种高效的锥面共形阵列天线盲极化DOA估计方法，新算法将子空间估计过程实值化处理，通过子阵分割设计解决了协方差矩阵分裂引起的“镜像模糊”问题，并推导了参数估计的克拉美—罗界，所提算法在保证DOA估计精度的前提下降低了算法复杂度。

附图说明

图1为锥面共形阵列示意图；

图2为含阵元锥面共形阵列示意图；

图3为锥面共形阵列天线俯视图；

图4为对母线上阵元的分割示意图；

图5为估计参数误差与信噪比分析图，极化幅度比为5，相差为30°，实验执行次数5000次；信噪比-5dB到15dB，步进1dB；阵元数30个，快拍次数100次；波达方向：波达1，俯仰角：65度，方位角：80度；波达2，俯仰角：80度，方位角：95度，图5a为本专利与文献[3]对于波达俯仰角的估计偏差，图5b为本专利与文献[3]对于波达方位角的估计偏差，图5c为本专利与文献[3]对于波达俯仰角的估计方差，图5d为本专利与文献[3]对于波达方位角的估计方差，图5e为本专利与文献[3]对于波达俯仰角的识别概率，图5f为本专利与文献[3]对于波达方位角的识别概率。

图6为估计参数误差与快拍次数分析图；极化幅度比为5，相差为30°，实验执行次数5000次；信噪比5dB；阵元数30个，快拍次数1000次起，步长100次，结束于5000次，波达方向：波达1，俯仰角：65度，方位角：80度；波达2，俯仰角：75度，方位角：90度；其中，图6a为本专利与文献[3]对于波达俯仰角的估计偏差，图6b为本专利与文献[3]对于波达方位角的估计偏差，图6c为本专利与文献[3]对于波达俯仰角的估计方差，图6d为本专利与文献[3]对于波达方位角的估计方差，图6e为本专利与文献[3]对于波达俯仰角的识别概率，图6f为本专利与文献[3]对于波达方位角的识别概率；

图7为估计参数误差与阵元数分析图，极化幅度比为5，相差为30°，实验执行次数5000次；信噪比0dB；快拍次数200次；阵元数150至600个，步进30个.波达方向：波达1，俯仰角：65度，方位角：80度；波达2，俯仰角：75度，方位角：90度；其中，图7a为本专利与文献[3]对于波达俯仰角的估计偏差，图7b为本专利与文献[3]对于波达方位角的估计偏差，图7c为本专利与文献[3]对于波达俯仰角的估计方差，图7d为本专利与文献[3]对于波达方位角的估计方差，图7e为本专利与文献[3]对于波达俯仰角的识别概率，图7f为本专利与文献[3]对于波达方位角的识别概率；

图8为估计参数的理论界限对比示意图。

图中文献[3]为齐子森，王布宏，郭英等，发表的《锥面共形阵列天线盲极化DOA估计算法》，发表在《电子学报》，2009，37(9)。

具体实施方式

下面结合附图和本发明的实例，对本发明作进一步的描述。

如图1-4所示，一种锥面共形阵列盲极化波达方向估计方法，包括下列步骤：步骤一：对锥面共形阵列天线快拍进行数据建模；

Step 1)利用锥面共形阵列的对称性，将锥面共形阵列划分为多个子阵，每一子阵的波达方向估计算法均相同；在锥面共形阵列的圆锥面的母线上划分l₁，l₂，l₃三个子阵，每个子阵上有m个阵元均匀分布，并且处于同一条母线上的相邻阵元间距为λ/2，其中，λ为入射信号波长；子阵l₁，l₂，l₃在全局坐标系XOY面上进行投影，三条母线与X轴顺时针的旋转方向的夹角分别为α₁，α₂，α₃，且α₁＝π-α₃，α₂∈(α₁，α₃)，以上子阵设置可实现方位角

的估计，将整个圆锥上的多个子阵综合起来，实现/>

的方位角估计；

Step 2)由共形阵列天线快拍数据建模方法，得锥面共形阵列的导向矢量为：

式(1)中，r_i为第i个阵元在全局坐标系中对单位强度来波信号的响应，G_i为第i个阵元在全局坐标系中的位置矢量，u为来波信号的方向矢量，θ和φ分别为来波信号相对于全局坐标系的俯仰角和方位角，j为虚数；

当有n个来波信号时，阵列的流形矩阵为：

A＝[a(θ₁，φ₁)，a(θ₂，φ₂)，…，a(θ_n，φ_n)] (2)

其中，此锥面共形阵列的快拍数据模型为：

X＝AS+N＝(A_θK_θ+A_φK_φ)S+N (3)

K_θ＝diag(k_1θ，k_2θ，…，k_nθ) (4)

S＝[s₁，s₂，…，s_n]^T (6)

N＝[n₁，n₂，…，n_n]^T (7)

式(3)-(7)中，K为对角阵，S为信号矢量，A_θ、A_φ分别为来波信号在极化正交基矢量上的流形矩阵，K_θ、K_φ分别为极化正交基矢量上投影分量构成的对角阵，K_nθ、K_nφ分别为第n个极化正交量上投影分量构成的对角阵，sn为第n个辐射源，n_n为第n个加性噪声，N为噪声矢量，θ和

分别为入射信号在全局坐标系中的俯仰角和方位角。

步骤二：实值盲极化波达方向估计方法，所述的步骤二具体包括下述步骤：

Step 1)在已求出快拍数据矩阵X的基础上，按照式(8)得到此阵列的数据协方差矩阵R；

可将R-AOCM和I-AOCM表示为：

式(8)-(10)中，N为噪声矢量，X^H为快拍数据矩阵的共轭转置，Re(X)为快拍数据矩阵的实部，Im(X)为快拍数据矩阵的虚部，Im^T(X)为快拍数据矩阵的虚部转置，Re^T(X)为快拍数据矩阵的实部的转置；

Step 2)采用实值化算法，在R-ACOM上进行特征分解，得到实值的信号子空间，从而分别得到母线l₁、l₂、l₃对应的Φ₁、Φ₂、Φ₃，其中，Φ₁、Φ₂、Φ₃分别为三条母线上阵元构成子阵的相位差；

设原始的信号子空间span(S)及其共轭信号子空间span(S^*)的交集，假设此交集为span(G)，表示为：

式(11)中，G为信号子空间与共轭信号子空间交集的基矢量，S为信号子空间基矢量，S^*为信号子空间基矢量的共轭；

在对信号协方差矩阵R的实部或者虚部直接进行特征分解的处理后，提取信号子空间时，需要提取的特征向量应该翻倍，每个子阵对均由同一条母线上分布的m个阵元构造，即1～m-1阵元，构成第一个阵列，2～m阵元，构成另一个阵列，两个阵列之间的距离为λ/2.以此类推，可分别得到母线l₁、l₂、l₃上距离矢量各不相同的三个子阵对；对母线l₁上分布的阵元进行分割，子阵l₁₁的接收数据为：

子阵l₁₂的接收数据为：

其中，Φ₁为相位差矩阵：

Φ₁＝diag[exp(-jw₁₁)，…，exp(-jw_1n)] (14)

w_1i＝(2π/λ)dv₁·u_iπv₁·u_i (15)

展开即得：

X₁₁和X₁₂的最后一行组成了母线l₁的快拍数据X₁：

X₁＝A₁S+N₁＝[X₁₁；X₁₂(m-1，：)] (17)

式(12)-(17)中，N₁₁表示子阵l₁₁接收的噪声矢量，A₁为母线1上阵元构成阵列的流形矩阵，A₁₁，A₁₂分别为母线1上阵元构成的子阵1和子阵2的流形矩阵，S表示信号矢量，A_11θ、A₁₁分别为母线1上阵元构成的子阵1在极化正交基矢量上的流形矩阵，K_θ、

分别为极化正交基矢量上投影分量构成的对角阵，Φ1为子阵相位差构成的对角阵，j表示虚部，w_1i为母线1上子阵对对第i个信号的响应相差，u_i为第i个信号方向矢量，θ_v1、/>

为母线1子阵对对应阵元距离矢量的表征参数，θ_i、/>

分别表示第i个来波的俯仰角与方位角，N1为母线1上阵列的加性噪声矩阵，N₁₂表示子阵l₁₂接收的噪声矢量，v₁为子阵l₁₁与子阵l₁₂距离的方向矢量，同理，将l₂、l₃划分为l₂₁、l₂₂和l₃₁、l₃₂，即可获得其快拍数据X₂和X₃以及Φ₂和Φ₃；

Step 3)对Φ₁、Φ₂、Φ₃进行特征分解，得到t_1i、t_2i、t_3i；

通过Φ₁，Φ₂，Φ₃这几个包含旋转不变关系的矩阵，结合式(15)求解信号的方向参数；

阵列的数据协方差矩阵为：

利用式(9)或(10)替换协方差矩阵，完成矩阵实值化处理，并进行特征分解：

U_S为阵列信号子空间的特征矢量矩阵，根据式(17)，U_S的第1～m-1行为子阵l₁₁对应的信号子空间，将其记为U_S11；U_S的第2～m行为子阵l₁₂对应的信号子空间，将其记为U_S12，同理得到U_S21，U_S22，U_S31，U_S32；

式(18)-(19)中，U_N为噪声子空间，

为信号子空间的共轭转置；/>

为噪声子空间的共轭转置；

对应ESPRIT里的最小二乘法：

式(20)中，

为母线1上子阵1的信号子空间的共轭转置，U_S1为母线1上子阵1的信号子空间；

且式(14)中的exp(-jw_1i)(i＝1，2，…，n)即为Φ₁的特征值t_1i，同理可得Φ₁的特征值t_1i，Φ₂的特征值t_2i，Φ₃的特征值t_3i；

Step 4)根据式(21)-式(23)，得到对应的N个接收信号的方向信息；

结合α₁＝π-α₃，α₂∈(α₁，α₃)分别求解t_1i，t_2i，t_3i，就可求出入射信号的方位角

入射信号的俯仰角θ_i：

z_ji＝angle(t_ji)/(2πd)，j＝1,2,3 (23)式(21)-(23)中，z_ji为母线j上阵列对第i个信号的接收相位；

Step 5)通过限定子阵的空域覆盖范围，对估计出的信号俯仰角和方位角进行判定，排除虚假信号，从而得到最终准确的信号方位。

实例：人体目标步态精细识别仿真实验

实验1：估计参数误差与信噪比分析

实验条件：阵列结构如图5a-图5f所示，极化幅度比为5，相差为30°，实验执行次数5000次；信噪比-5dB到15dB，步进1dB；阵元数30个，快拍次数100次；波达方向：波达1，俯仰角：65度，方位角：80度；波达2，俯仰角：80度，方位角：95度。

由仿真实验可知，所提算法可以解决多辐射源方位估计问题，随着信噪比的增加，本专利实值高效算法的DOA估计成功概率、估计方差以及估计偏差都逐步趋近文献[3]算法，当信噪比大于5dB时，对方位角的估计效果基本与文献[3]相同，对俯仰角的估计效果与文献[3]比存在明显差距，当信噪比大于5dB时，本专利算法的俯仰角估计偏差也在0.5°以内，满足空间目标的高分辨区分要求.

实验2：估计参数误差与快拍次数

实验条件：阵列结构如图6a-图6f所示，极化幅度比为5，相差为30°，实验执行次数5000次；信噪比5dB；阵元数30个，快拍次数1000次起，步长100次，结束于5000次，波达方向：波达1，俯仰角：65度，方位角：80度；波达2，俯仰角：75度，方位角：90度。由仿真实验可知，快拍次数的提升对本专利算法以及文献[3]算法的DOA估计效果影响并非十分明显，当快拍数在3000以上时，估计效果随着快拍的增加变化趋于平缓，整体估计精度的变化并不如信噪比对估计结果的影响.本专利算法对俯仰角与方位角的估计效果均未随快拍增加与文献[3]算法趋于一致，存在一定差距.这也说明协方差矩阵分裂造成了数据信息的丢失，影响了所采样本的信息含量，带来的差距不能单纯靠数据累计来补偿.

实验3：估计参数误差与阵元数分析

实验条件：阵列结构如图7a-图7f所示，极化幅度比为5，相差为30°，实验执行次数5000次；信噪比0dB；快拍次数200次；阵元数150至600个，步进30个.波达方向：波达1，俯仰角：65度，方位角：80度；波达2，俯仰角：75度，方位角：90度。由仿真实验可知，阵元个数的增加对参数估计效果的提升十分明显，成功概率、估计偏差以及估计方差均随着阵元个数增加快速与文献[4]算法的估计效果趋于一致，这也说明参数估计效果对天线口径的变化极其敏感.

实验4：参数估计理论界仿真

实验条件：阵列结构如图8所示，极化幅度比为5，相差为30°，实验执行次数5000次；信噪比-5dB到15dB，步进1dB；阵元数30个，快拍次数100次；波达方向：波达1，俯仰角：65度，方位角：80度；波达2，俯仰角：75度，方位角：90度.文献[3]算法的参数估计CRB推导详见该文献.由仿真结果可知：对协方差矩阵进行实值化处理后，参数的估计性能有损失，估计的理论线均大于实值化处理之前.但随着信噪比的增加，逐步趋于一致.仿真还呈现出俯仰角的估计效果好于方位角，主要原因是图1给出的阵列结构中母线上阵元数多于锥体横截面阵元数，所以对俯仰角分辨能力强，方位角差一些.

就运算复杂度而言，ESPRIT算法的运算复杂度主要集中在对数据协方差矩阵R的处理上，所以本专利仅对协方差矩阵R特征分解部分的计算量进行讨论，表1给出了锥面共形阵列天线盲极化DOA估计算法和锥面共形阵列天线实值盲极化DOA估计算法的运算复杂度对比情况其中，算法的运算复杂度均以各算法中包含的实数乘法运算次数来表示，单元4×O(M²L)给出了计算复数矩阵R的逆矩阵或EVD所包含的实数乘法次数，其中，M为阵元个数，L为仿真实验的快拍次数。

表1运算复杂度(实数乘法次数)对比

V算法	运算复杂度
		文献[3]算法	4×O(M2L)
本专利算法	O(M2L)

由表1可以看出，锥面共形阵列天线实值盲极化DOA估计算法在运算复杂度上相较于锥面共形阵列天线盲极化DOA估计算法，降低了约75％的运算量，按照此对比，节省时间大约在75％左右，但是在实际仿真中，由于本专利所采用的实值化算法是直接对数据协方差矩阵的实部进行特征分解，此操作导致我们为了准确提取到真实信号方位的信息，必须将真实信号和虚假信号的特征信息同时提取出来，进行一系列处理直至求出方位信息后，再加以甄别，最终得出准确的信号方位。这其中甄别真假信号的过程是原算法所没有的，所以，最终节省时间的百分比会小于75％，具体数值则与这部分操作所占时间的百分比有关。

实验5：算法仿真时间分析

实验条件：阵列结构如图1所示，极化幅度比为5，相差为30°，信噪比10dB；快拍次数100次；运行次数100次；波达方向：波达1，俯仰角：65度，方位角：80度；波达2，俯仰角：80度，方位角：95度；阵元数，150至600，步进150。

表2算法仿真时间对比

阵元数	150	300	450	600
					运算时间(文献[3])	0.2614	0.7047	1.7096	3.2995
运算时间(本专利)	0.1795	0.4366	1.0198	1.8644
					运算提高比	31.3％	38.0％	40.3％	43.5％

由计算仿真时间表2对比可以看出，算法经过实值化处理后可显著降低算法复杂度，总体运行时间是闫锋刚，沈毅，刘帅等的《高效超分辨波达方向估计算法综述》，发表在《系统工程与电子技术》，2015第3期上的算法的60％-70％之间，而且天线规模越大，阵元个数越多，复杂度下降比例越大，充分验证了算法的高效性。

本发明的有益效果在于：提供了一种高效的锥面共形阵列天线盲极化DOA估计方法，新算法将子空间估计过程实值化处理，通过子阵分割设计解决了协方差矩阵分裂引起的“镜像模糊”问题，并推导了参数估计的克拉美—罗界，所提算法在保证DOA估计精度的前提下降低了算法复杂度。

Claims (2)

1.一种锥面共形阵列盲极化波达方向估计方法，包括下列步骤：

步骤一：对锥面共形阵列天线快拍进行数据建模；

当有n个来波信号时，阵列的流形矩阵为：

A＝[a(θ₁,φ₁),a(θ₂,φ₂),…,a(θ_n,φ_n)] (2)

其中，此锥面共形阵列的快拍数据模型为：

X＝AS+N＝(A_θK_θ+A_φK_φ)S+N (3)

K_θ＝diag(k_1θ,k_2θ,…,k_nθ) (4)

K_φ＝diag(k_1φ,k_2φ,…,k_nφ) (5)

S＝[s₁,s₂,…,s_n]^T (6)

N＝[n₁,n₂,…,n_n]^T (7)

式(3)-(7)中，S为信号矢量，A_θ、A_φ

分别为来波信号在极化正交基矢量上的流形矩阵，k_θ、k_φ

分别为极化正交基矢量上投影分量构成的对角阵，k_nθ、k_nφ

分别为第n个极化正交量上投影分量构成的对角阵，s_n为第n个辐射源，n_n为第n个加性噪声，N为噪声矢量，θ和φ分别为来波信号在全局坐标系中的俯仰角和方位角；

步骤二：实值盲极化波达方向估计方法；

Step 1)在已求出快拍数据矩阵X的基础上，按照式(8)得到此阵列的数据协方差矩阵R；

可将R-AOCM和I-AOCM表示为：

式(8)-(10)中，N为噪声矢量，X^H为快拍数据矩阵的共轭转置，Re(X)为快拍数据矩阵的实部，Im(X)为快拍数据矩阵的虚部，Im^T(X)为快拍数据矩阵的虚部转置，Re^T(X)为快拍数据矩阵的实部的转置；

Step 2)采用实值化算法，在R-ACOM上进行特征分解，得到实值的信号子空间，从而分别得到母线l₁、l₂、l₃对应的Φ₁、Φ₂、Φ₃，其中，Φ₁、Φ₂、Φ₃分别为三条母线上阵元构成子阵的相位差；

设原始的信号子空间span(S)及其共轭信号子空间span(S^*)的交集，假设此交集为span(G)，表示为：

式(11)中，G为信号子空间与共轭信号子空间交集的基矢量，S为信号子空间基矢量，S^*为信号子空间基矢量的共轭；

在对信号协方差矩阵R的实部或者虚部直接进行特征分解的处理后，提取信号子空间时，需要提取的特征向量翻倍，每个子阵对均由同一条母线上分布的m个阵元构造，即1～m-1阵元，构成第一个阵列，2～m阵元，构成另一个阵列，两个阵列之间的距离为λ/2，以此类推，可分别得到母线l₁、l₂、l₃上距离矢量各不相同的三个子阵对；对母线l₁上分布的阵元进行分割，子阵l₁₁的接收数据为：

X₁₁＝A₁₁S+N₁₁＝(A_11θK_θ+A_11φK_φ)S+N₁₁ (12)

子阵l₁₂的接收数据为：

X₁₂＝A₁₂S+N₁₂＝(A_11θK_θ+A_11φK_φ)Φ₁S+N₁₂ (13)

其中，Φ₁为相位差矩阵：

Φ₁＝diag[exp(-jw₁₁),…,exp(-jw_1n)] (14)

w_1i＝(2π/λ)dv₁·u_i＝πv₁·u_i (15)

展开即得：

w_1i＝π[sin(θ_v1)cos(φ_v1)sin(θ_i)sin(φ_i)+sin(θ_v1)sin(φ_v1)sin(θ_i)sin(φ_i)+cos(θ_v1)cos(θ_i)] (16)

X₁₁和X₁₂的最后一行组成了母线l₁的快拍数据X₁：

X₁＝A₁S+N₁＝[X₁₁；X₁₂(m-1,:)] (17)

式(12)-(17)中，N₁₁表示子阵l₁₁接收的噪声矢量，A₁为母线1上阵元构成阵列的流形矩阵，A₁₁，A₁₂分别为母线1上阵元构成的子阵1和子阵2的流形矩阵，S表示信号矢量，A_11θ、A_11φ分别为母线1上阵元构成的子阵1在极化正交基矢量上的流形矩阵，K_θ、K_φ分别为极化正交基矢量上投影分量构成的对角阵，Φ₁为子阵相位差构成的对角阵，j表示虚部，w_1i为母线1上子阵对对第i个信号的响应相差，u_i为第i个信号方向矢量，θ_v1、φ_v1为母线1子阵对对应阵元距离矢量的表征参数，θ_i、φ_i分别表示第i个来波的俯仰角与方位角，N₁为母线1上阵列的加性噪声矩阵，N₁₂表示子阵l₁₂接收的噪声矢量，v₁为子阵l₁₁与子阵l₁₂距离的方向矢量，同理，将l₂、l₃划分为l₂₁、l₂₂和l₃₁、l₃₂，即可获得其快拍数据X₂和X₃以及Φ₂和Φ₃；

Step 3)对Φ₁、Φ₂、Φ₃进行特征分解，得到t_1i、t_2i、t_3i；

通过Φ₁，Φ₂，Φ₃这几个包含旋转不变关系的矩阵，结合式(15)求解信号的方向参数；

阵列的数据协方差矩阵为：

利用式(9)或(10)替换协方差矩阵，完成矩阵实值化处理，并进行特征分解：

U_S为阵列信号子空间的特征矢量矩阵，根据式(17)，U_S的第1～m-1行为子阵l₁₁对应的信号子空间，将其记为U_S11；U_S的第2～m行为子阵l₁₂对应的信号子空间，将其记为U_S12，同理得到U_S21，U_S22，U_S31，U_S32；

式(18)-(19)中，U_N为噪声子空间，

为信号子空间的共轭转置；/>

为噪声子空间的共轭转置；

对应ESPRIT里的最小二乘法：

式(20)中，

为母线1上子阵1的信号子空间的共轭转置，U_S11为母线1上子阵1的信号子空间；

且式(14)中的exp(-jw_1i),i＝1,2,…,n，即为Φ₁的特征值t_1i，同理可得Φ₁的特征值t_1i，Φ₂的特征值t_2i，Φ₃的特征值t_3i；

Step 4)根据式(21)-式(23)，得到对应的N个接收信号的方向信息；

结合α₁＝π-α₃，α₂∈(α₁，α₃)分别求解t_1i，t_2i，t_3i，就可求出来波信号的方位角φ_i：

来波信号的俯仰角θ_i：

z_ji＝angle(t_ji)/(2πd)，j＝1,2,3 (23)

式(21)-(23)中，z_ji为母线j上阵列对第i个信号的接收相位；

Step 5)通过限定子阵的空域覆盖范围，对估计出的信号俯仰角和方位角进行判定，排除虚假信号，从而得到最终准确的信号方位。

2.如权利要求1所述的一种锥面共形阵列盲极化波达方向估计方法，其特征在于：所述步骤一的具体步骤为：

步骤一：对锥面共形阵列天线快拍进行数据建模；

Step 1)利用锥面共形阵列的对称性，将锥面共形阵列划分为多个子阵，每一子阵的波达方向估计算法均相同；在锥面共形阵列的圆锥面的母线上划分l₁，l₂，l₃三个子阵，每个子阵上有m个阵元均匀分布，并且处于同一条母线上的相邻阵元间距为λ/2，其中，λ为来波信号波长；子阵l₁，l₂，l₃在全局坐标系XOY面上进行投影，三条母线与X轴顺时针的旋转方向的夹角分别为α₁,α₂,α₃,且α₁＝π-α₃，α₂∈(α₁，α₃)，以上子阵设置可实现方位角φ∈(α₁，α₃)的估计，将整个圆锥上的多个子阵综合起来，实现φ∈(0，2π)的方位角估计；

Step 2)由共形阵列天线快拍数据建模方法，得锥面共形阵列的导向矢量为：

式(1)中，r_i为第i个阵元在全局坐标系中对单位强度来波信号的响应，G_i为第i个阵元在全局坐标系中的位置矢量，u为来波信号的方向矢量，θ和φ分别为来波信号相对于全局坐标系的俯仰角和方位角，j为虚数；

得到锥面共形阵列的导向矢量。

Images