互质面阵下的二维DOA跟踪方法与流程

本发明涉及阵列信号处理技术领域，尤其是一种互质面阵下的二维doa跟踪方法。

背景技术：

相对于传统面阵，互质面阵的主要思想是通过联合两个满足互质关系的子面阵进行空间谱估计，不仅能有效解决测向模糊问题，还可以获得较大的阵列孔径，显著提高波束形成和doa估计中的自由度，因此得到广泛应用。传统的经典超分辨率doa估计算法需要对接收信号的协方差矩阵进行特征值分解或者奇异值分解，计算量太大，难以满足doa实时跟踪的需求，因此研究doa跟踪算法变得很有实际意义。针对这个问题，国内外学者对其中的难点提出了诸多方法。其中比较有代表性的算法是子空间跟踪算法，包括past、pastd、双迭代最小二乘(bi-iterativeleastsquares，bi-ils)算法、双边迭代奇异值分解(bi-iterationsingularvaluedecomposition，bi-svd)算法等。本发明将互质面阵与pastd算法相结合，提出互质面阵下的doa实时跟踪算法，实现doa角度的自动配对与跟踪，由于无需进行协方差矩阵的构造和特征值分解，算法复杂度低。

技术实现要素：

本发明所要解决的技术问题在于，提供一种互质面阵下的二维doa跟踪方法，能够实现互质面阵下的任意动态目标的低复杂度的doa实时跟踪。

为解决上述技术问题，本发明提供一种互质面阵下的二维doa跟踪方法，包括如下步骤：

(1)通过联合互质面阵中的两个子面阵构建阵列接收信号的数据模型；

(2)利用pastd方法实时更新接收数据的信号子空间；

(3)基于阵列接收数据的信号子空间，利用互质面阵下的二维esprit算法得到doa估计，并通过联合阵列子阵间的角度关系实现解模糊。

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

阵元总数为m²+n²-1的互质面阵，由阵元数分别为n×n和m×m的两个均匀子面阵构成，且两个子阵的阵元间距分别为mλ/2和nλ/2；对于子面阵i，其在x轴和y轴上阵元的方向矩阵可以分别表示为其中，方向向量di为阵元间距；因此，子面阵i的接收信号xi(t)可以表示为

其中，s(t)为信源矢量，ni(t)为子面阵i的噪声，和о分别表示kronecker积和khatri-rao积。

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

(21)选择适当的初始值λk(0)，w(0)；

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

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

(24)当步骤(23)中的k＝k后，使得t＝t+1，再次从步骤(22)开始计算。

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

在互质面阵中，对于子面阵i，按下式构造ai1，ai2

其中，mi＝m或mi＝n，ai2＝ai1φiy，旋转因子

将信号子空间eis分解成eix＝es(1:mi(mi-1),:)，其中eix为eis的1到mi(mi-1)行，eiy为eis的mi+1到行，eix、eiy也可表达为eix＝ai1ti，eiy＝ai1φiyti，矩阵满秩，进而有

其中，则φiy的对角线的元素就是ψi的特征值，由最小二乘法则可以得到ψi的估计

将特征值分解可以获得然后根据的特征矢量，就能得到t的估计结果在不考虑噪声的条件下有

πi表示一个置换矩阵；因为和的特征值是一致的，所以将特征值分解就可以得到的估计为

其中为的第k个特征值；

同样对矩阵eis重构，可得类似地可以将e'is分解成矩阵e'ix＝e'is(1:mi(mi-1),:)和矩阵构造矩阵那么，在不考虑噪声的条件下有其中πi表示一个置换矩阵，类似于uik，可以得到的估计值

其中表示的是矩阵(e'ix)⁺e'iy第k个特征值；

通过联合两个子阵间的角度关系，可以消除互质面阵中角度估计的模糊值，即两个子阵相同的估计值即为真实角度，其余值为模糊值。

根据uik和vik的表达式可知，uik和vik的估计值有相同的列模糊，因此估计角度可以自动配对，角度配对完成后，根据以下公式即可得到信源仰角和方位角的估计值

本发明的有益效果为：本发明无需分解特征值和构造接收数据的协方差矩阵就能实现2d-doa信息的实时跟踪，复杂度低；对比互质阵下基于past进行doa跟踪的方法，本发明具有更好的跟踪性能。

附图说明

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

图2为本发明的pastd算法流程示意图。

图3为本发明实现的信源方位角的60s跟踪图。

图4为本发明实现的信源俯仰角的60s跟踪图。

图5为本发明实现的信源doa的60s跟踪图。

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

图7为本发明与past算法在不同信噪比下的跟踪性能对比示意图。

具体实施方式

一、数据模型

互质面阵的阵列结构如图1所示，与传统均匀面阵相比，互质面阵可以分成两个子阵，n表示第一个子阵在x轴和y轴方向的阵元数，m表示第二个子阵在x轴和y轴方向的阵元数。子阵一是由阵元数为n×n的均匀面阵构成，子阵二由阵元数为m×m的均匀面阵构成。子阵阵元在原点位置重合，故互质面阵的阵元总数为m²+n²-1。其中子阵一的阵元间距为d1＝mλ/2，子阵二的阵元间距为d2＝nλ/2，各阵元的位置可表示为如下集合

ds＝{(md1,nd1)|0≤m,n≤n-1}∪{(pd2,qd2)|0≤p,q≤m-1}

假设空间有k个非相干信源入射到上述互质面阵，θk,分别代表入射信源的俯仰角和方位角。考虑面阵中的子面阵i，其在x轴和y轴上阵元的方向矩阵可以分别表示为其中，方向向量因此，子面阵i的接收信号xi(t)可以表示为

其中，s(t)为信源矢量，ni(t)为子面阵i的噪声。和分别表示kronecker积和khatri-rao积。

二、利用pastd方法实时更新信号子空间

定义一个无约束代价函数

j(w)＝e{||x-ww^hx||²}

＝e{(x-ww^hx)^h(x-ww^hx)}

＝e{x^hx}-2e{x^hww^hx}+e{x^hww^hww^hx}

可以发现

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)

由binyang提出的定理可以知道当目标函数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)可以用自相关矩阵cyy(t)和互相关矩阵cxy(t)来表示，minj(w(t))最优解是wiener滤波器，即

其中，cyy(t)和cxy(t)更新公式为：

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

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

1)选择适当的初始值λk(0)，w(0)；

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

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

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

pastd的算法流程图如图2所示，算法最后一步通过xk(t)减去c(t)的第k个特征向量wk(t)达到算法紧缩。

三、利用2d-esprit得到doa信息

在互质面阵中，对于子面阵i，按下式构造ai1，ai2

其中，mi＝m或mi＝n，ai2＝ai1φiy，旋转因子

将由pastd得到的信号子空间eis分解成eix＝es(1:mi(mi-1),:)，其中eix为eis的1到mi(mi-1)行，eiy为eis的mi+1到行，eix、eiy也可表达为eix＝ai1ti，eiy＝ai1φiyti。矩阵满秩。进而有

其中，则φiy的对角线的元素就是ψi的特征值。由最小二乘法则可以得到ψi的估计

将特征值分解可以获得然后根据的特征矢量，就能得到t的估计结果在不考虑噪声的条件下有

πi表示一个置换矩阵。因为和的特征值是一致的，所以将特征值分解就可以得到的估计为

其中为的第k个特征值。

同样对矩阵eis重构，可得类似地可以将e'is分解成矩阵e'ix＝e'is(1:mi(mi-1),:)和矩阵构造矩阵那么，在不考虑噪声的条件下有其中πi表示一个置换矩阵。类似于uik，可以得到的估计值

其中表示的是矩阵(e'ix)⁺e'iy第k个特征值。

由于互质面阵中阵元间距大于半波长，因此会在非满秩面阵中会产生角度模糊。假设有空间角度为的非相干信源入射到互质面阵，考虑面阵中的子面阵i，子面阵i在x轴和y轴方向上相邻阵元的相位差与阵元间距之间的关系可分别表示为

其中kx和ky都为整数，取值范围分别为同时kx和ky还需满足当和确定时，存在一个或多个空间角度满足上式，当d＝mλ/2(m为大于1的整数)时，kx和ky的取值分别为m和m/2个，即使考虑到角度配对问题，仍有多个角度估计值满足上式，这些角度中只有一个为真实角度，其他均为模糊值。大量文献已证明，互质面阵的角度模糊消除可以通过联合两个子阵间的角度关系解决，两个子阵相同的估计值即为真实角度，其余值为模糊值。

根据uik和vik的表达式可知，uik和vik的估计值有相同的列模糊，因此估计角度可以自动配对。角度配对完成后，根据以下公式即可得到信源仰角和方位角的估计值

下面利用matlab仿真对本发明的算法性能进行分析。其中，采用求根均方误差(rootmeansquareerror，rmse)来评估算法doa估计性能，rmse定义如下

其中,k表示空间中的信源个数，l表示montecarlo试验次数。和分别表示第q次montecarlo试验时第k个信源仰角θk和方位角的估计值，θk,q和分别表示其精确值。

n和m分别表示互质面阵中的两个子阵在x轴和y轴方向的阵元数，阵列阵元总数为m²+n²-1。假设信源是以直线或曲线运动，方位角和仰角在时间t内线性变化，且整个跟踪过程中信源数保持不变。跟踪时间t设为60s。每隔1s跟踪一次，1s内快拍数j设为200。遗忘因子β设为0.97。pastd算法的初始参数设置如下：λk＝1(k＝1,…,k)，w(0)＝[ik,0]^t，ik为k×k的单位阵。

图3-5是snr为20db时基于pastd算法的doa跟踪结果，阵列参数m＝4，n＝5。从图3-5可以看出该算法能够有效的对信源的方位角和仰角进行跟踪。图5显示了两个信源的方位角和仰角在一个图中的变化轨迹，可看出本发明能够较精确的进行doa跟踪。

图6是本发明在不同阵元数下的doa跟踪性能比较图。在快拍数j为300的情况下，固定参数m为5，n＝[3,4,6]，可以看出随着阵元数的增加，本发明的跟踪性能变得越来越好。

图7显示了本发明与past算法在不同信噪比下的跟踪性能对比图，past算法同样基于图1的互质面阵模型。pastd算法随着信噪比提高性能逐步改善，而在信噪比升高后性能增长减缓，同时可以看出本发明性能优于past算法。本发明在高信噪比时仍可以较精确的实现信源二维doa跟踪。

本发明的复杂度分析：在图1所示的阵列模型下，pastd算法每次更新的复杂度为o(4(m²+n²)k+2k)，在互质面阵中得到信号子空间的情况下利用2d-esprit算法进行二维doa估计的复杂度为o(4k²m(m-1)+4k²n(n-1)+16k³)。由此可见，本发明的复杂度较低，利于实时跟踪doa信息。

符号说明：小写(大写)粗体字来表示向量(矩阵)。(·)^t、(·)^h分别表示矩阵或向量的转置、共轭转置。e(·)是统计期望。和分别表示kronecker积和khatri-rao积。diag(·)代表使用向量的元素作为对角元素的对角矩阵。angle(·)表示取相角。triu(·)表示取矩阵上三角元素。表示对值x的估计。