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

本发明涉及阵列信号处理
技术领域：
，尤其是一种共形圆阵下的低复杂度二维doa估计方法。
背景技术：
：随着人们对目标定位精度要求的日益提高，传统的利用波束机械扫描的测向方法，在速度、精度和分辨率上都已无法满足实际应用的要求，而共形阵列的广泛使用不仅可以节省空间，最大限度地减少天线对飞行器空气动力学性能的影响，还可以扩展天线波束扫描范围，有效提高电磁兼容性能。与经典线阵、面阵不同，由于受到共形载体曲率的影响，共形阵列天线的单元方向图指向不一致，呈现出多极化特性，使得共形阵列天线的快拍数据建模时引入阵列入射信号的极化参数。因此，信源方位与极化状态的耦合是共形阵列下高分辨doa估计的难点。目前，对共形阵列高分辨doa估计方法的已有研究主要集中在导向矢量建模和简化模型条件下的经典算法移植。本发明针对共形圆阵这类特殊共形阵列结构，提出了该阵列下的低复杂度盲极化doa估计算法，实现了精确的doa估计。技术实现要素：本发明所要解决的技术问题在于，提供一种共形圆阵下的低复杂度二维doa估计方法，能够实现共形圆阵下的低复杂度doa精确估计。为解决上述技术问题，本发明提供一种共形圆阵下的低复杂度二维doa估计方法，包括如下步骤：(1)利用欧拉旋转实现阵列阵元极化方向图的全局旋转变换，解决共形阵列天线的多极化问题，从而构建共形圆阵下的阵列数据模型；(2)基于pastd方法得到阵列接收信号的信号子空间；(3)根据谱函数的结构特点，将四维music谱峰搜索降至二维music谱峰搜索，实现共形圆阵下的低复杂度二维doa估计。优选的，步骤(1)具体为：共形天线的导向矢量由入射信号参数(θ,φ,γ,η)共同决定，其中，γ,η为入射信号的极化参数，θ,φ为入射信号的俯仰角和方位角，对于一个由m个相同的全向阵列构成的均匀圆阵，其方向矩阵为：由于一般阵元方向图的定义和设计都是以本地局部坐标系作为参考，因此需要利用欧拉旋转实现阵元极化方向图的全局旋转变换，即将局部阵元方向图转换为全局阵元方向图gm(θ,φ)，每个阵元的变换步骤如下：(11)将全局坐标系内(θ,φ)处的单位矢量进行直角坐标表示：x＝sinθcosφ,y＝sinθsinφ,z＝cosθ(12)利用欧拉旋转变换将全局直角坐标变换至阵元的局部直角坐标并得到局部极坐标中对应的方位与阵列结构对应的旋转变换欧拉角及欧拉旋转变换矩阵分别定义为：dm＝2(m-1)π/m,em＝0,fm＝0(13)由阵元的局部极坐标响应得到其在局部直角坐标系下的表示：式中，为第m个阵元在局部坐标系下方向图的极化表示，并存在以下关系：(14)由阵元方向图的局部直角坐标系表示以及欧拉旋转逆变换得到阵元方向图的全局直角坐标表示：(15)最后将全局直角坐标下的阵元方向图转换成全局极坐标表示，得到gmθ,gmφ:gmθ(θm,φm)＝-gmz/sinθgmφ(θm,φm)＝-gmxsinφ+gmycosφ因此，阵列的接收信号模型为：式中，为接收数据矢量，为入射信号矢量；为阵列噪声矢量。优选的，步骤(2)具体为：(21)选择适当的初始值λn(0)，w(0)；(22)对每一个t＝1,2,…,j(j为快拍数)，使得x1(t)＝x(t)；(23)对每一个n＝1,2,…,n(n为信源数)，分别更新以下变量：阵列接收数特征值特征向量以及阵列接收数据xn+1(t)＝xn(t)-wn(t)yn(t)；(24)当步骤(23)中的n＝n后，使得t＝t+1，再次从步骤(22)开始计算；pastd算法最后一步通过xn(t)减去c(t)的第n个特征向量wn(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轴的夹角为：其位置矢量为rm＝(rcosδm,rsinδm,0)。本文采用球面坐标系表示入射平面波的波达方向，坐标系的原点o在阵列的中心即圆心。信源的仰俯角θ∈[0,π/2]是z轴与信源入射方向的夹角，方位角φ∈[0,2π]则是从x轴沿逆时针方向到信源入射方向在阵列入射平面上投影的夹角。设有n个窄带信号源从远场辐射到天线阵列，第n个信源到达角度为(θn,φn)。由于共形天线阵列中阵元坐标系的旋转关系和极化分量的旋转关系，使得共形天线的导向矢量由入射信号参数(θ,φ,γ,η)共同决定。其中，γ,η为入射信号的极化参数，具体定义为:对于传播方向上的任意一点，在该点的传播横截面上电场矢量的端点为随时间变化的一个极化椭圆，椭圆的形状、倾角和旋向取决于两个方向电磁场幅度比以及相位差。定义，tanγ＝ax/ay表示y方向的电场幅度与x方向的电场幅度比，η＝φy-φx表示y方向电场和x方向电场的相位差，取值范围为γ∈[0,π/2]，η∈[0,2π)。因此对于图1所示的共形圆阵，其方向矩阵为：其中极化分量p(γ,η)＝[p(γ1,η1),p(γ2,η2),…,p(γn,ηn)]，方向向量若用pmi(m＝1,…,m)表示第i个入射信号的极化矢量在第m个阵元的极化方向图上的投影，如图2所示，则p(γi,ηi)＝[p1i,p2i,…,pmi]t式中，pmi＝ui·gm，gm为第m个阵元方向图在全局坐标系中的正交极化分量表示。eθ和eφ为全局坐标系中正交极化基矢量；和为第i个入射信号极化在全局坐标系中的正交极化分解表示，一般阵元方向图的定义和设计都是以本地局部坐标系作为参考，即表示第m个阵元在局部坐标系内的方向图，因此需要利用欧拉旋转实现阵元极化方向图的全局旋转变换，即将局部阵元方向图转换为全局阵元方向图gm(θ,φ)。因此gm表示为：gm＝gmθeθ+gmφeφ式中，gmθ,gmφ为利用欧拉旋转实现的阵元极化方向图的全局旋转变换，表示第m个阵元方向图在全局坐标系中的正交极化分量。每个阵元的变换步骤如下：步骤1将全局坐标系内(θ,φ)处的单位矢量进行直角坐标表示：x＝sinθcosφ,y＝sinθsinφ,z＝cosθ步骤2利用欧拉旋转变换将全局直角坐标变换至阵元的局部直角坐标并得到局部极坐标中对应的方位与图1阵列结构对应的旋转变换欧拉角及欧拉旋转变换矩阵分别定义为：dm＝2(m-1)π/m,em＝0,fm＝0步骤3由阵元的局部极坐标响应得到其在局部直角坐标系下的表示：式中，为第m个阵元在局部坐标系下方向图的极化表示，并存在以下关系：步骤4由阵元方向图的局部直角坐标系表示以及欧拉旋转逆变换得到阵元方向图的全局直角坐标表示：步骤5最后将全局直角坐标下的阵元方向图转换成全局极坐标表示，得到gmθ,gmφ:gmθ(θm,φm)＝-gmz/sinθgmφ(θm,φm)＝-gmxsinφ+gmycosφ因此，阵列的接收信号模型为：式中，为接收数据矢量，为入射信号矢量；为阵列噪声矢量。二、利用pastd方法得到噪声子空间定义一个无约束代价函数可以发现e{xhwwhx}＝tr(e{whxxhw})＝tr(whcw)e{xhwwhwwhx}＝tr(e{whwxxhwhw})＝tr(whcwwhw)其中c表示接收数据x的自相关矩阵。假设w秩为n，则j(w)可以表示为：j(w)＝tr(c)-2tr(whcw)+tr(whcwwhw)binyang提出下列两个重要的定理：定理1：w是j(w)的一个平衡点，当且仅当w＝urq,ur是m×n的矩阵，q是n×n的任意酉矩阵。特征向量不在矩阵ur的特征值的和即为j(w)的值。定理2：只有ur是矩阵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)xh(t)，y(t)＝wh(t-1)x(t)代入上式，得到w(t)＝w(t-1)-μ[2x(t)yh(t)-x(t)yh(t)×wh(t-1)w(t-1)-w(t-1)wh(t-1)y(t)yh(t)]由于j(w)取得极小值时，wh(t-1)w(t-1)＝i，可得w(t)＝w(t-1)+μ[x(t)-w(t-1)y(t)]yh(t)由于w(t)跟踪时变子空间的能力较差，算法收敛慢。为了解决这个问题可以定义一个新的指数加权函数其中0<β≤1,表示遗忘因子，主要是保证在非稳定环境下将过去的数据降低权重，以保证跟踪的稳定性。当β＝1时对应常见的滑动窗口。进一步，可以认为y(i)＝wh(i-1)x(i)≈wh(t)x(i)因此得到已经修正的目标函数当目标函数全局最小时，w(t)可以用自相关矩阵cyy(t)和互相关矩阵cxy(t)来表示，minj(w(t))最优解是wiener滤波器，即其中，cyy(t)和cxy(t)更新公式为：pastd方法主要利用紧缩技术通过一定顺序来估计主分量，首先在目标数n为1时通过past方法估计出当前最主要的特征向量wn(t)，再把当前t时刻的数据xn(t)减去它在特征向量wn(t)上的投影得到xn+1(t)，然后重复上述步骤计算出wn+1(t),wn+2(t)…。利用pastd方法求解信源信号子空间的具体步骤如下：1)选择适当的初始值λn(0)，w(0)；2)对每一个t＝1,2,…,j(j为快拍数)，使得x1(t)＝x(t)；3)对每一个n＝1,2,…,n(n为信源数)，分别更新以下变量：阵列接收数yn(t)＝特征值特征向量以及阵列接收数据xn+1(t)＝xn(t)-wn(t)yn(t)；4)当步骤(3)中的n＝n后，使得t＝t+1，再次从步骤(2)开始计算；pastd的算法流程图如图3所示，算法最后一步通过xn(t)减去c(t)的第n个特征向量wn(t)达到算法紧缩。三、谱峰搜索得到doa信息由pastd算法得到噪声子空间un，由阵列信号处理的基础知识可知，阵列方向矢量张成的子空间与噪声子空间正交。因此定义阵列空间谱函数为将上式中的改写成令则谱函数可写成令显然有q(θi,φi)＝qh(θi,φi)，由于p(γi,ηi)为非零向量，p(γi,ηi)≠0，当参数(θ,φ,γ,η)＝(θi,φi,γi,ηi),i＝1,2,…,n时，有此时有ph(γi,ηi)q(θi,φi)p(γi,ηi)＝0；当参数(θ,φ,γ,η)≠(θi,φi,γi,ηi),i＝1,2,…,n时，ph(γ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仿真对本发明的算法性能进行分析。其中，采用求根均方误差(rootmeansquareerror,rmse)来评估算法doa估计性能，rmse定义如下：其中,n表示空间中的信源个数，l表示montecarlo试验次数。和分别表示第q次montecarlo试验时第n个信源仰角θn和方位角φn的估计值，和分别表示其精确值。表1是算法复杂度对比。其中，uv表示谱峰搜索次数。从表1可以看出本发明的复杂度低于基于降维music进行doa估计的复杂度。表1算法复杂度对比基于降维music的共形圆阵估计算法o{m3+m2j+(m3+m+1)uv}本发明提出的算法o{(4mn+n)j+(m2+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性能随着阵元数量的增加而提高。当前第1页12