CN109188342A - 共形圆阵下的低复杂度二维doa估计方法 - Google Patents

Abstract

本发明公开了一种共形圆阵下的低复杂度二维DOA估计方法，利用PASTd方法得到接收数据的信号子空间，然后基于降维MUSIC算法将共形圆阵的四维谱峰搜索降至二维谱峰搜索，实现共形圆阵下的低复杂度DOA估计。本发明的有益效果为：对比共形圆阵下基于降维MUSIC的DOA估计算法，本发明算法复杂度低，利于硬件实现；对比共形圆阵下基于PAST的DOA估计算法，本发明具有更好的DOA估计性能。

Description

共形圆阵下的低复杂度二维DOA估计方法

技术领域

本发明涉及阵列信号处理技术领域，尤其是一种共形圆阵下的低复杂度二维DOA估计方法。

背景技术

随着人们对目标定位精度要求的日益提高，传统的利用波束机械扫描的测向方法，在速度、精度和分辨率上都已无法满足实际应用的要求，而共形阵列的广泛使用不仅可以节省空间，最大限度地减少天线对飞行器空气动力学性能的影响，还可以扩展天线波束扫描范围，有效提高电磁兼容性能。与经典线阵、面阵不同，由于受到共形载体曲率的影响，共形阵列天线的单元方向图指向不一致，呈现出多极化特性，使得共形阵列天线的快拍数据建模时引入阵列入射信号的极化参数。因此，信源方位与极化状态的耦合是共形阵列下高分辨DOA估计的难点。目前，对共形阵列高分辨DOA估计方法的已有研究主要集中在导向矢量建模和简化模型条件下的经典算法移植。本发明针对共形圆阵这类特殊共形阵列结构，提出了该阵列下的低复杂度盲极化DOA估计算法，实现了精确的DOA估计。

发明内容

本发明所要解决的技术问题在于，提供一种共形圆阵下的低复杂度二维DOA估计方法，能够实现共形圆阵下的低复杂度DOA精确估计。

为解决上述技术问题，本发明提供一种共形圆阵下的低复杂度二维DOA估计方法，包括如下步骤：

(1)利用欧拉旋转实现阵列阵元极化方向图的全局旋转变换，解决共形阵列天线的多极化问题，从而构建共形圆阵下的阵列数据模型；

(2)基于PASTd方法得到阵列接收信号的信号子空间；

(3)根据谱函数的结构特点，将四维MUSIC谱峰搜索降至二维MUSIC谱峰搜索，实现共形圆阵下的低复杂度二维DOA估计。

优选的，步骤(1)具体为：

共形天线的导向矢量由入射信号参数(θ,φ,γ,η)共同决定，其中，γ,η为入射信号的极化参数，θ,φ为入射信号的俯仰角和方位角，对于一个由M个相同的全向阵列构成的均匀圆阵，其方向矩阵为：

由于一般阵元方向图的定义和设计都是以本地局部坐标系作为参考，因此需要利用欧拉旋转实现阵元极化方向图的全局旋转变换，即将局部阵元方向图转换为全局阵元方向图g_m(θ,φ)，每个阵元的变换步骤如下：

(11)将全局坐标系内(θ,φ)处的单位矢量进行直角坐标表示：

x＝sinθcosφ,y＝sinθsinφ,z＝cosθ

(12)利用欧拉旋转变换将全局直角坐标变换至阵元的局部直角坐标并得到局部极坐标中对应的方位与阵列结构对应的旋转变换欧拉角及欧拉旋转变换矩阵分别定义为：

D_m＝2(m-1)π/M,E_m＝0,F_m＝0

(13)由阵元的局部极坐标响应得到其在局部直角坐标系下的表示：

式中，为第m个阵元在局部坐标系下方向图的极化表示，并存在以下关系：

(14)由阵元方向图的局部直角坐标系表示以及欧拉旋转逆变换得到阵元方向图的全局直角坐标表示：

(15)最后将全局直角坐标下的阵元方向图转换成全局极坐标表示，得到g_mθ,g_mφ:

g_mθ(θ_m,φ_m)＝-g_mZ/sinθ

g_mφ(θ_m,φ_m)＝-g_mXsinφ+g_mYcosφ

因此，阵列的接收信号模型为：

式中，为接收数据矢量，为入射信号矢量；为阵列噪声矢量。

优选的，步骤(2)具体为：

(21)选择适当的初始值λ_n(0)，W(0)；

(22)对每一个t＝1,2,…,J(J为快拍数)，使得x₁(t)＝X(t)；

(23)对每一个n＝1,2,…,N(N为信源数)，分别更新以下变量：阵列接收数 特征值特征向量 以及阵列接收数据x_n+1(t)＝x_n(t)-W_n(t)y_n(t)；

(24)当步骤(23)中的n＝N后，使得t＝t+1，再次从步骤(22)开始计算；PASTd算法最后一步通过x_n(t)减去C(t)的第n个特征向量W_n(t)达到算法紧缩。

优选的，步骤(3)具体为：

根据阵列信号处理的基础知识，定义阵列空间谱函数为

将上式中的改写成

令则谱函数可写成

令易证Q(θ,φ)为半正定矩阵，当且仅当(θ,φ)＝(θ_i,φ_i)，i＝1,2,…,N时，Q(θ,φ)为奇异矩阵，即Q(θ,φ)在信源的真实入射方向处奇异，此时对矩阵Q(θ,φ)有det(Q(θ_i,φ_i))＝0,i＝1,...,N，由此可得共形圆阵的DOA估计谱函数

利用上式进行二维搜索即可得到入射信号的DOA信息。

本发明的有益效果为：对比共形圆阵下基于降维MUSIC的DOA估计算法，本发明算法复杂度低，利于硬件实现；对比共形圆阵下基于PAST的DOA估计算法，本发明具有更好的DOA估计性能。

附图说明

图1为本发明的阵列结构示意图。

图2为本发明的入射信号极化在阵元极化方向图的投影图。

图3为本发明的PASTd算法流程示意图。

图4为本发明在不同信噪比下的各个算法的RMSE性能对比示意图。

图5为本发明在不同快拍数下的各个算法的RMSE性能对比示意图。

图6为本发明的算法在不同阵元数下的RMSE性能对比示意图。

具体实施方式

符号：小写(大写)粗体字来表示向量(矩阵)。(·)^T、(·)^H分别表示矩阵或向量的转置、共轭转置。E(·)是统计期望。⊙表示Hadamard积运算。diag(·)代表使用向量的元素作为对角元素的对角矩阵。angle(·)表示取相角。triu(·)表示取矩阵上三角元素。表示对值x的估计。

本发明定义如下坐标系：[x,y,z]表示阵列全局直角坐标，表示阵元局部直角坐标，[r,θ,φ]表示阵元全局极坐标，表示阵元局部极坐标。

一、数据模型

共形圆阵的阵列结构如图1所示，M个相同的全向阵列均匀分布在平面X-Y上一个半径为R的圆周上。相邻两阵元间的夹角为β＝2π/M，阵元序号m按逆时针方向排序，第m个阵元同圆心之间的连线与X轴的夹角为：

其位置矢量为r_m＝(Rcosδ_m,Rsinδ_m,0)。

本文采用球面坐标系表示入射平面波的波达方向，坐标系的原点O在阵列的中心即圆心。信源的仰俯角θ∈[0,π/2]是Z轴与信源入射方向的夹角，方位角φ∈[0,2π]则是从X轴沿逆时针方向到信源入射方向在阵列入射平面上投影的夹角。设有N个窄带信号源从远场辐射到天线阵列，第n个信源到达角度为(θ_n,φ_n)。

由于共形天线阵列中阵元坐标系的旋转关系和极化分量的旋转关系，使得共形天线的导向矢量由入射信号参数(θ,φ,γ,η)共同决定。其中，γ,η为入射信号的极化参数，具体定义为:对于传播方向上的任意一点，在该点的传播横截面上电场矢量的端点为随时间变化的一个极化椭圆，椭圆的形状、倾角和旋向取决于两个方向电磁场幅度比以及相位差。定义，tanγ＝A_x/A_y表示Y方向的电场幅度与X方向的电场幅度比，η＝φ_y-φ_x表示Y方向电场和X方向电场的相位差，取值范围为γ∈[0,π/2]，η∈[0,2π)。因此对于图1所示的共形圆阵，其方向矩阵为：

其中极化分量P(γ,η)＝[p(γ₁,η₁),p(γ₂,η₂),…,p(γ_N,η_N)]，方向向量 若用p_mi(m＝1,…,M)表示第i个入射信号的极化矢量在第m个阵元的极化方向图上的投影，如图2所示，则

p(γ_i,η_i)＝[p_1i,p_2i,…,p_Mi]^T

式中，p_mi＝u_i·g_m，g_m为第m个阵元方向图在全局坐标系中的正交极化分量表示。E_θ和E_φ为全局坐标系中正交极化基矢量；和为第i个入射信号极化在全局坐标系中的正交极化分解表示，

一般阵元方向图的定义和设计都是以本地局部坐标系作为参考，即表示第m个阵元在局部坐标系内的方向图，因此需要利用欧拉旋转实现阵元极化方向图的全局旋转变换，即将局部阵元方向图转换为全局阵元方向图g_m(θ,φ)。因此g_m表示为：

g_m＝g_mθE_θ+g_mφE_φ

式中，g_mθ,g_mφ为利用欧拉旋转实现的阵元极化方向图的全局旋转变换，表示第m个阵元方向图在全局坐标系中的正交极化分量。每个阵元的变换步骤如下：

步骤1将全局坐标系内(θ,φ)处的单位矢量进行直角坐标表示：

x＝sinθcosφ,y＝sinθsinφ,z＝cosθ

步骤2利用欧拉旋转变换将全局直角坐标变换至阵元的局部直角坐标并得到局部极坐标中对应的方位与图1阵列结构对应的旋转变换欧拉角及欧拉旋转变换矩阵分别定义为：

D_m＝2(m-1)π/M,E_m＝0,F_m＝0

步骤3由阵元的局部极坐标响应得到其在局部直角坐标系下的表示：

式中，为第m个阵元在局部坐标系下方向图的极化表示，并存在以下关系：

步骤4由阵元方向图的局部直角坐标系表示以及欧拉旋转逆变换得到阵元方向图的全局直角坐标表示：

步骤5最后将全局直角坐标下的阵元方向图转换成全局极坐标表示，得到g_mθ,g_mφ:

g_mθ(θ_m,φ_m)＝-g_mZ/sinθ

g_mφ(θ_m,φ_m)＝-g_mXsinφ+g_mYcosφ

因此，阵列的接收信号模型为：

式中，为接收数据矢量，为入射信号矢量；为阵列噪声矢量。

二、利用PASTd方法得到噪声子空间

定义一个无约束代价函数

可以发现

E{x^HWW^Hx}＝tr(E{W^Hxx^HW})＝tr(W^HCW)

E{x^HWW^HWW^Hx}＝tr(E{W^HWxx^HW^HW})＝tr(W^HCWW^HW)

其中C表示接收数据x的自相关矩阵。假设W秩为N，则J(W)可以表示为：

J(W)＝tr(C)-2tr(W^HCW)+tr(W^HCWW^HW)

Bin Yang提出下列两个重要的定理：定理1：W是J(W)的一个平衡点，当且仅当W＝U_rQ,U_r是M×N的矩阵，Q是N×N的任意酉矩阵。特征向量不在矩阵U_r的特征值的和即为J(W)的值。定理2：只有U_r是矩阵C的N个主特征向量组成的情况下，目标函数J(W)取得极小值。在其他任意情况下，J(W)的平衡点都是鞍点。由以上定理可以知道当目标函数J(W)取极小值时，W的列空间等价于信号子空间。

在现实情况中，为了更新得到t时刻的子空间W(t)，需要利用t-1时刻的子空间W(t-1)以及t时刻的阵元数据x(t)。在此我们选择梯度下降法，选取下降梯度为

所以

其中，μ>0表示需要适当选择的步长值。将C(t)＝x(t)x^H(t)，y(t)＝W^H(t-1)x(t)代入上式，得到

W(t)＝W(t-1)-μ[2x(t)y^H(t)-x(t)y^H(t)

×W^H(t-1)W(t-1)-W(t-1)W^H(t-1)y(t)y^H(t)]

由于J(W)取得极小值时，W^H(t-1)W(t-1)＝I，可得

W(t)＝W(t-1)+μ[x(t)-W(t-1)y(t)]y^H(t)

由于W(t)跟踪时变子空间的能力较差，算法收敛慢。为了解决这个问题可以定义一个新的指数加权函数

其中0<β≤1,表示遗忘因子，主要是保证在非稳定环境下将过去的数据降低权重，以保证跟踪的稳定性。当β＝1时对应常见的滑动窗口。进一步，可以认为

y(i)＝W^H(i-1)x(i)≈W^H(t)x(i)

因此得到已经修正的目标函数

当目标函数全局最小时，W(t)可以用自相关矩阵C_yy(t)和互相关矩阵C_xy(t)来表示，minJ(W(t))最优解是Wiener滤波器，即

其中，C_yy(t)和C_xy(t)更新公式为：

PASTd方法主要利用紧缩技术通过一定顺序来估计主分量，首先在目标数N为1时通过PAST方法估计出当前最主要的特征向量W_n(t)，再把当前t时刻的数据x_n(t)减去它在特征向量W_n(t)上的投影得到x_n+1(t)，然后重复上述步骤计算出W_n+1(t),W_n+2(t)…。

利用PASTd方法求解信源信号子空间的具体步骤如下：

1)选择适当的初始值λ_n(0)，W(0)；

2)对每一个t＝1,2,…,J(J为快拍数)，使得x₁(t)＝X(t)；

3)对每一个n＝1,2,…,N(N为信源数)，分别更新以下变量：阵列接收数y_n(t)＝特征值特征向量 以及阵列接收数据x_n+1(t)＝x_n(t)-W_n(t)y_n(t)；

4)当步骤(3)中的n＝N后，使得t＝t+1，再次从步骤(2)开始计算；

PASTd的算法流程图如图3所示，算法最后一步通过x_n(t)减去C(t)的第n个特征向量W_n(t)达到算法紧缩。

三、谱峰搜索得到DOA信息

由PASTd算法得到噪声子空间U_N，由阵列信号处理的基础知识可知，阵列方向矢量张成的子空间与噪声子空间正交。因此定义阵列空间谱函数为

将上式中的改写成

令则谱函数可写成

令显然有Q(θ_i,φ_i)＝Q^H(θ_i,φ_i)，由于p(γ_i,η_i)为非零向量，p(γ_i,η_i)≠0，当参数(θ,φ,γ,η)＝(θ_i,φ_i,γ_i,η_i),i＝1,2,…,N时，有此时有p^H(γ_i,η_i)Q(θ_i,φ_i)p(γ_i,η_i)＝0；当参数(θ,φ,γ,η)≠(θ_i,φ_i,γ_i,η_i),i＝1,2,…,N时，p^H(γ_i,η_i)Q(θ_i,φ_i)p(γ_i,η_i)>0。因此Q(θ,φ)为半正定矩阵，当且仅当(θ,φ)＝(θ_i,φ_i)，i＝1,2,…,N时，Q(θ,φ)为奇异矩阵，即Q(θ,φ)在信源的真实入射方向处奇异，此时对矩阵Q(θ,φ)有det(Q(θ_i,φ_i))＝0,i＝1,...,N。由此可得共形圆阵的DOA估计谱函数

利用上式进行二维搜索即可得到入射信号的DOA信息。

下面利用MATLAB仿真对本发明的算法性能进行分析。其中，采用求根均方误差(Root Mean Square Error,RMSE)来评估算法DOA估计性能，RMSE定义如下：

其中,N表示空间中的信源个数，L表示Monte Carlo试验次数。和分别表示第q次Monte Carlo试验时第n个信源仰角θ_n和方位角φ_n的估计值，和分别表示其精确值。

表1是算法复杂度对比。其中，uv表示谱峰搜索次数。从表1可以看出本发明的复杂度低于基于降维MUSIC进行DOA估计的复杂度。

表1算法复杂度对比

基于降维MUSIC的共形圆阵估计算法	O{M<sup>3</sup>+M<sup>2</sup>J+(M<sup>3</sup>+M+1)uv}
		本发明提出的算法	O{(4MN+N)J+(M<sup>2</sup>+M+1)uv}

图4是不同信噪比下的各个算法的RMSE性能对比图。其中，阵元数M＝7，快拍数J＝300，信源数K＝2，入射信号的方向分别为(20°,5°)、(40°,25°)。从图4可以看出，本发明的RMSE性能随着信噪比的增加而提高，且优于基于PAST的共形圆阵DOA估计算法。

图5是不同快拍数下的各个算法的RMSE性能对比图。其中，阵元数M＝7，信噪比为20dB，信源数K＝2，入射信号的方向分别为(20°,5°)、(40°,25°)。从图5可以看出，本发明的RMSE性能随着快拍数的增加而提高，且优于基于PAST的共形圆阵DOA估计算法。

图6是本发明在不同阵元数下的性能对比图。其中，信噪比为20dB，信源数K＝2，快拍数J＝300，入射信号的方向分别为(20°,5°)、(40°,25°)。图6表明，本发明的RMSE性能随着阵元数量的增加而提高。

Claims (4)

1.共形圆阵下的低复杂度二维DOA估计方法，其特征在于，包括如下步骤：

(1)利用欧拉旋转实现阵列阵元极化方向图的全局旋转变换，解决共形阵列天线的多极化问题，从而构建共形圆阵下的阵列数据模型；

(2)基于PASTd方法得到阵列接收信号的信号子空间；

(3)根据谱函数的结构特点，将四维MUSIC谱峰搜索降至二维MUSIC谱峰搜索，实现共形圆阵下的低复杂度二维DOA估计。

2.如权利要求1所述的共形圆阵下的低复杂度二维DOA估计方法，其特征在于，步骤(1)具体为：

共形天线的导向矢量由入射信号参数(θ,φ,γ,η)共同决定，其中，γ,η为入射信号的极化参数，θ,φ为入射信号的俯仰角和方位角，对于一个由M个相同的全向阵列构成的均匀圆阵，其方向矩阵为：

由于一般阵元方向图的定义和设计都是以本地局部坐标系作为参考，因此需要利用欧拉旋转实现阵元极化方向图的全局旋转变换，即将局部阵元方向图转换为全局阵元方向图g_m(θ,φ)，每个阵元的变换步骤如下：

(11)将全局坐标系内(θ,φ)处的单位矢量进行直角坐标表示：

x＝sinθcosφ,y＝sinθsinφ,z＝cosθ

(12)利用欧拉旋转变换将全局直角坐标变换至阵元的局部直角坐标并得到局部极坐标中对应的方位与阵列结构对应的旋转变换欧拉角及欧拉旋转变换矩阵分别定义为：

D_m＝2(m-1)π/M,E_m＝0,F_m＝0

(13)由阵元的局部极坐标响应得到其在局部直角坐标系下的表示：

式中，为第m个阵元在局部坐标系下方向图的极化表示，并存在以下关系：

(14)由阵元方向图的局部直角坐标系表示以及欧拉旋转逆变换得到阵元方向图的全局直角坐标表示：

(15)最后将全局直角坐标下的阵元方向图转换成全局极坐标表示，得到g_mθ,g_mφ:

g_mθ(θ_m,φ_m)＝-g_mZ/sinθ

g_mφ(θ_m,φ_m)＝-g_mXsinφ+g_mYcosφ

因此，阵列的接收信号模型为：

式中，为接收数据矢量，为入射信号矢量；为阵列噪声矢量。

3.如权利要求1所述的共形圆阵下的低复杂度二维DOA估计方法，其特征在于，步骤(2)具体为：

(21)选择适当的初始值λ_n(0)，W(0)；

(22)对每一个t＝1,2,…,J(J为快拍数)，使得x₁(t)＝X(t)；

(23)对每一个n＝1,2,…,N(N为信源数)，分别更新以下变量：阵列接收数 特征值特征向量 以及阵列接收数据x_n+1(t)＝x_n(t)-W_n(t)y_n(t)；

(24)当步骤(23)中的n＝N后，使得t＝t+1，再次从步骤(22)开始计算；PASTd算法最后一步通过x_n(t)减去C(t)的第n个特征向量W_n(t)达到算法紧缩。

4.如权利要求1所述的共形圆阵下的低复杂度二维DOA估计方法，其特征在于，步骤(3)具体为：

根据阵列信号处理的基础知识，定义阵列空间谱函数为

将上式中的改写成

令p(γ_i,η_i)＝[p(γ₁,η₁) … p(γ_M,η_M)]^T，则谱函数可写成

令易证Q(θ,φ)为半正定矩阵，当且仅当(θ,φ)＝(θ_i,φ_i)，i＝1,2,…,N时，Q(θ,φ)为奇异矩阵，即Q(θ,φ)在信源的真实入射方向处奇异，此时对矩阵Q(θ,φ)有det(Q(θ_i,φ_i))＝0,i＝1,...,N，由此可得共形圆阵的DOA估计谱函数

利用上式进行二维搜索即可得到入射信号的DOA信息。