CN105549005B - 一种基于网格划分的动态目标波达方向跟踪方法 - Google Patents

Abstract

本发明公开了一种基于网格划分的动态目标波达方向跟踪方法，首先在模拟的环境中布置传感器阵列；然后建立远场传感器阵列的窄带信号测量模型；构建空间含噪声的源信号模型；把平面环境平均划分网格，每个网格的信号入射角的范围大小相同；再构建目标转移模型，预测目标下一刻可能的状态；结合目标的状态转移模型和传感器测量模型生产卡尔曼滤波方程，求解出划分的每个网格中的信号状态值；最后根据网格中的状态值，利用合适的方法确定目标处于哪个网格，从而确定目标的位置。该方法利用网格划分构建目标转移和信号测量的线性方程，再结合卡尔曼滤波方法实现目标跟踪，克服了传统方法的计算复杂，计算量大，并且最后估计精度不足的缺点。

Description

一种基于网格划分的动态目标波达方向跟踪方法

技术领域

本发明涉及一种基于网格划分的动态目标波达方向跟踪方法，广泛应用于基于无线传感器的动态目标实时跟踪。

背景技术

动态目标跟踪是利用监测到移动物体发射或反射的信号通过一定算法对移动的物体运动轨迹进行评估。首先需要使用特定的接收信号的设备不断地接收信号，并对信号进行滤波等预处理，得到需要的信号信息；然后对得到的信号通过特定的一些算法进行滤波或信号变换等操作直接或间接求解出需要的值。

目前来看，目标追踪技术在民用和军事两个领域都被广泛的应用着，在有的研究当中分为单目标追踪和多目标追踪两大类。不管什么样的分类，目的都是一个，把获取到的目标信息传输给终端电脑，然后进行各种分析。单目标追踪主要有空间追踪和时间-空间追踪。目前提高追踪精度和对追踪的数据进行融合已成为目标追踪的主要研究方向。当前单目标的追踪过程主要还是进行着简单环境的研究，利用计算技术和通信技术更高效的信息融合获取有效信息，以及如何以较少的传感器节点能量消耗来提高测量精度和传感器网络寿命。

波达方向(Direction Of Arrive,DOA)估计作为阵列信号处理中一项重要的研究内容，近年来得到了广泛的关注，并取得了一系列的研究成果。在雷达、声纳等领域的工程实际应用中，由于目标信号源通常是移动的，因此需要对运动目标的DOA进行准确的跟踪估计。DOA代表性算法有MUSIC，ESPRIT等，然而这些算法往往只用来分析当目标处于静止状态的情况。如果目标信号方向在不断变化，那么这些算法需要反复的分解接收数据的协方差矩阵，无论是从实时性还是计算量上来说都是不容易实现的。基于DOA的跟踪方法主要有三类：

1)子空间跟踪算法

子空间跟踪算法是一种自适应估计阵列协方差矩阵的最小特征值所对应的特征向量的方法。这类算法通过各种途径更新子空间(信号子空间和噪声子空间)信息，避免数据协方差矩阵不断重复分解或奇异值分解，有效降低了计算量。这类算法具有计算量小、实时性好的优点，但其在低信噪比情况下估计性能欠佳，而且本身不具备解相干能力。

2)以分段估计代替实时估计

这类方法将传统的DOA估计方法直接推广到DOA跟踪过程，通过时间分段，利用每一段时间的DOA估计值代替DOA瞬时值。这类方法思路简单，在目标DOA变化较小或几乎 不变的条件下效果较好，其缺点是实时性和精度不高，不能处理DOA变化较快的场合。

3)滤波类DOA跟踪方法

这类的有效跟踪方法有很多。如经典的扩展卡尔曼滤波(Extend Kalman Filter,EKF)是使用较多的算法，但EKF仅仅利用了非线性函数的泰勒级数展开式的一阶项，而当系统高度非线性或非高斯时，EKF方法将导致滤波发散。还有无迹卡尔曼滤波(unscentedKalman Filter,UKF)、粒子滤波(Particle Filter,PF)等。

发明内容

本发明的目的是针对现有技术存在的不足，提出一种基于网格划分的动态目标波达方向跟踪方法，该方法克服了传统DOA跟踪技术的计算复杂度高，计算量大，计算时间长以至于难以满足实际应用的需求等问题。达到计算量小，计算复杂度低，可以很好的满足实际应用的实时需求。

为了达到上述目的，本发明的构思是：在平面环境布置传感器阵列，首先建立远场传感器阵列的窄带信号测量模型；构建空间含噪声的源信号模型；然后把平面环境平均划分网格，每个网格的信号入射角的范围大小相同，构建阵列信号的网格划分平面模型；再构建目标转移模型，预测目标下一刻可能的状态；最后结合目标的状态转移模型和传感器测量模型生产卡尔曼滤波方程，可以求解出目标最佳估计位置。

根据上述发明构思，本发明采用的技术方案是：

一种基于网格划分的动态目标波达方向跟踪方法，包括以下几个步骤：

(1)在监测区域建立直角坐标系，在坐标系的一个坐标轴上等间隔距离的布置传感器；

(2)构建远场传感器阵列的窄带信号测量模型及含噪声的源信号模型；

(3)把监测平面环境平均划分网格G等份，每个网格的信号入射角的范围大小相同，即，π/G；

(4)在平面环境中构建目标合理的转移模型，为了实现卡尔曼滤波跟踪，构建卡尔曼方程，所以这个转移模型应当是线性的；

(5)结合目标的转移模型和信号的测量模型，再根据卡尔曼滤波的原理，快速的低复杂度的估计出的动态目标当前位置。

所述步骤(2)中构建远场传感器阵列的窄带信号测量模型和含噪声的源信号模型，具体如下：

传感器阵列在直角坐标系的横轴上等间隔距离布置，一般间隔为信号的一半波长，即λ/2(λ为信号波长)，假设第一个传感器阵元的坐标为(0,0)，则其它阵元的坐标为(λ/2,0)、 (λ,0)、…、(N*λ/2,0)；传感器的测量模型表示为：

X(t)＝A*S(t)+R(t) (1)

其中，t为时刻，S(t)为t时刻在区域的M个目标所发射的信号，是M*1维向量；R(t)为t时刻N个传感器接收噪声，是N*1维向量，A为阵列的N*M维的导向矢量阵，X(t)为N个传感器接收的信号值，是N*1向量；

上述中的A＝[a(θ₁),a(θ₂),...,a(θ_M)]^T

其中表达式中的d为传感器位置间隔距离，λ为信号波长，θ_i表示第i个目标到达传感器的角度，i＝1,2,…,M，而α(θ_i)表示N个传感器对第i个目标的导向向量。

所述步骤(3)中把监测平面环境平均划分网格，使得测量方程为线性的，具体如下：

网格划分是以原点为基点，把坐标轴的上半区域平均划分G等份，即每份的范围大小为π/G，划分网格后，式(1)中的A就变成有固定的维数的N*G维的导向矩阵，不会根据在监测区域中的目标个数M变化而变化，因此A由如下等式表示：

A＝[a(θ₁),a(θ₂),...,a(θ_G)]^T

其中的θ_i(i＝1,2,…G)在划分网格时是已知的，例如θ₁＝π/G，θ₂＝2π/G，θ_G-1＝(G-1)π/G，因此也需要相应的改变式(1)中的S(t)，把S(t)变成G*1维的向量，但是S(t)是稀疏的，即在某个网格中有目标，在S(t)对应的位置上才有值，其他的位置上都为0。(1)式中的X(t)和R(t)都将变成G*1维向量。

所述步骤(4)中构建目标状态概率转移模型，具体如下：

假设目标以一种概率在三个相邻的网格间移动，即目标当前网格，与其相邻的下一个网格，与其相邻的上一个网格，

假设目标以1/3的概率停留在当前网格，或移动到相邻的两个网格，则目标的状态移动方程可以表示如下：

S(t)＝F*S(t-1)+Q(t) (2)

其中S(t-1)表示前一时刻的信号的值，F是G*G维的转移概率矩阵，Q(t)为当前时刻t的 转移噪声

其中f_ij是F的元素，i,j分别是F的行和列。

所述步骤(5)中利用卡尔曼滤波实现动态目标的位置实时估计，具体如下：

构建卡尔曼滤波方程

上式中的符号X(t)、S(t)、S(t-1)分别简化成X_t、S_t、S_t-1。根据卡尔曼滤波的计算步骤：

S_t|t-1＝F*S_t-1

P_t|t-1＝F*P_t-1F^T+Q

K＝P_t|t-1*A*(A*P_t|t-1*A^T+R)^-1

P_t＝(E_G*G-K*A)*P_t|t-1

S_t＝S_t|t-1+K*(X_t-A*S_t|t-1)

其中K为卡尔曼增益，F为转移概率矩阵，S_t为当前信号的值，S_t-1为前一时刻信号的值，S_t|t-1为当前时刻的信号预测值，P_t为当前时刻S_t的协方差矩阵，P_t-1为前一时刻S_t-1的协方差矩阵，P_t|t-1为当前时刻预测信号值的协方差矩阵，E_G*G为G*G维的单位矩阵；X_t为当前时刻的信号测量值。Q、R分别为高斯噪声Q(t)和R(t)的协方差。

根据上面的卡尔曼滤波计算步骤得到当前时刻的G个网格中的信号的值，需要根据这G个元素值计算出目标最后的位置，有两种可行的方法：

(1)把计算出的G个元素值中前M个最大信号值所处网格坐标作为M个目标的位置；

(2)根据当前时刻t的状态估计值S_t，根据聚类算法把S_t的G个元素分成M类，即假设监测区域有M个目标，然后求每类所有元素网格坐标的平均值，当作目标所处的位置。

与现有技术相比，本发明具有如下显而易见的优势：

本方法采用网格划分的策略，使得信号测量方程具有线性的特性，更加有利于计算；而动态目标的状态移动可以规定在网格间移动，这样目标的下一刻的状态就有更加明确的可预测性；由于卡尔曼滤波是最小均方误差意义上的最优估计，结合线性的测量方程和状态转移方程，可以达到最佳估计，从而提高精确度；同时，利用卡尔曼滤波可以大大的减少计算量和计算复杂度，减少计算时间，可以更好的满足实际需求。可广泛应用于基于DOA动态目标跟踪等领域。

本方法构建了线性的测量方程和状态转移方程，可以完美的结合卡尔曼滤波，提高目标定位跟踪精确，最主要的是可以大量的减少计算量和计算复杂的，提高计算速度，为把基于波达方向的目标跟踪技术应用于实际打下基础。

附图说明

图1为本方法的流程图。

图2为本发明坐标系统远场信号模型。

图3为平面划分网格及目标转移示意图。

具体实施方式

为了更好地理解本发明的技术方案，以下结合附图对本发明的具体实施例做进一步详细描述。

如图1所示，一种基于网格划分的动态目标波达方向跟踪方法，利用网格划分后的传感器测量特性，结合构建的线性状态转移模型，利用卡尔曼滤波算法对目标当前状态进行估计，具体实施如下：

(1)在监测区域建立直角坐标系，在坐标系的一个坐标轴上等间隔距离的的布置传感器。构建阵列接收信号。具体如下：

如图2所示，基于波达方向的信号处理主要依赖于传感器阵列，在需要监视的范围建立直角坐标系，在横轴上均匀等间隔的放置传感器节点，但是传感器节点的间隔距离是有要求的，一般为λ/2，其中λ为目标发射信号的波长。由于使用窄带远场信号，目标发射的信号在到达传感器时可以近似认为是平行的，目标信号到达不同传感器的值是不同的，根据图2，得到目标的角度位置。

τ*c＝(λ/2)sinθ θ＝arcsin(2τc/λ)

其中τ是两个传感器接收同一信号的时延，c为光速。θ为信号入射角度

因此得到目标与原点的连线的夹角θ，即目标的位置。

(2)构建远场传感器阵列的窄带信号测量模型及含噪声的源信号模型：

如图2所示，设有N个窄带远场信号到达M个阵元组成天线阵列上。假设阵元间不存在互藕以及相位不一致性，则接收信号表达式可以写为：

式中，u_i(t)是接收信号幅度，是接收信号的相位，τ是接收信号短暂延时，ω₀是信号角频率。远场窄带信号可以写成下面形式：

那么，由式(2)和(3)，可以得到下式:

则第k个阵元接收到的总N个远场信号的综合为：

式中，g_ki为第k个阵元对第i个信号的增益，n_k(t)为第k个阵元在t时刻的噪声，τ_ki表示第i个信号到达第k个阵元时相对于参考阵元的时延。

在理想情况下，式(5)中的增益g_ki可以省略(归一化为1)，因此可将M个阵元在同一时刻接收的信号排列成一个如下矩阵

写成矢量形式：

X(t)＝AS(t)+N_x(t) (7)

设采样数为L，则式中，X(t)为阵列的M*L维矢量，N_x(t)为阵列的M*L维噪声矢量，S(t)为空间信号的N*L维矢量，A为M*N维的空间阵列流型矩阵，且

A＝[a₁(ω₀) a₂(ω₀) ... a_N(ω₀)]

由此可知，只要能够求出阵元间的延迟的表达式，就能够写出空间阵列流形表达式矩阵。

对于线阵传感器阵列，设阵元的位置表示为X_k(k＝1,2,,,M)，参考点位于原点，入射信号的一维波达方向为θ_i(i＝1,2,…,N)，且该角是入射信号与Y轴的夹角，则延迟表达式变为：

因此，式(7)可表示如下矩阵形式

(3)把平面环境平均划分网格G等份，每个网格的信号入射角的范围大小相同。

在平面阵中，使用一种新颖的测量来使得式(7)能满足线性表示，由于式(7)中的A是随着检测区域发射的信号数量变化而变化的，这种特性不利于计算，因此，可以如图3所示，把检测区域划分固定的G个格子，每个格子的范围大小相同。这样式(10)中的θ_i(i＝1,2,…,G)具有固定的值，θ₁＝π/G，θ₂＝2π/G，θ_G-1＝(G-1)π/G。因此式(7)中的A就是固定的值，A＝[a(θ₁),a(θ₂),...,a(θ_G)]^T。

再来看S(t)，由于划分了G个网格，S(t)就是G*1维的矢量，有固定的维数，而在这些网格中，有些网格有目标，他们发射信号，而另外些网格没有目标信号，没有目标信号的网格在S(t)对应的位置的值为0。例如，只有在第1和2个网格有目标，则S(t)可表示为

S(t)＝[s₁(t) s₂(t) 0 ... 0]^T

因此，S(t)是根据网格中目标的个数和状态不同而变化的。X(t)和N_x(t)也是一样的。

(4)在平面环境中构建目标的转移模型。

如图3所示，在上述的步骤(3)中，构建平面网格模型，通过网格划分把测量过程线性化。目标在这些网格中转移，在一定条件下也可以把目标在网格的状态转移过程线性化。假设在检测平面划分网格足够多，使得一个网格在一个时刻最多只能容纳一个目标，也就是在同一时刻，两个目标物体不能在同一网格中。这样，在多目标跟踪时，不会因为两个目标在同一个网格而结果计算却只有一个目标，从而造成严重的误差。

本发明规定目标在下一个时刻最多只能移动一个网格的距离，而且只能移动到它当前所在网格的相邻的网格或停留在它当前所处的网格，即假设目标当前所处网格为第i个网格，那么它下一个时刻所处的位置有三个可能：第i-1个网格，第i个网格或第i+1个网格。而且规定目标移动三个网格的概率相等，即1/3。因此有：

S(t)＝F*S(t-1)+N_s(t)

(11)

式中，k表示时刻，S(t)表示发射信号，N_s(t)为转移噪声，F是G*G维的转移矩阵。

即

其中f_ij是F的元素，i,j分别是F的行和列

式(11)即为所需要的状态转移方程，它具有线性特性。

(5)结合目标的转移模型和信号的测量模型，再根据卡尔曼滤波的原理，可以快速的且低 复杂度的估计出动态目标当前位置。

卡尔曼滤波是一个最优化自回归数据处理算法，是最小均方误差为原则的最佳线性滤波器。根据前一个估计值和最近一个观察数据来估计信号的当前值，需要已知系统的状态方程和量测方程。

根据上面步骤(1)～(4)的分析，可以知道动态目标跟踪的状态方程和测量方程，即式(7)和式(11)，构建卡尔曼滤波方程：

简化上式符号，用X_t表示X(t)、S_t表示S(t)、S_t-1表示S(t-1)。卡尔曼滤波的计算过程可以如下所示：

S_t|t-1＝F*S_t-1

P_t|t-1＝F*P_t-1F^T+Q

K＝P_t|t-1*A*(A*P_t|t-1*A^T+R)^-1

P_t＝(E_G*G-K*A)*P_t|t-1

S_t＝S_t|t-1+K*(X_t-A*S_t|t-1)

其中K为卡尔曼增益，F为转移概率矩阵，S_t为当前信号的值，S_t-1为前一时刻信号的值，S_t|t-1为当前时刻的信号预测值，P_t为当前时刻S_t的协方差矩阵，P_t-1为前一时刻S_t-1的协方差矩阵，P_t|t-1为当前时刻预测信号值的协方差矩阵，E_G*G为G*G维的单位矩阵；X_t为当前时刻的信号测量值。Q、R分别为高斯噪声N_x(t)和N_s(t)的协方差。

根据状态转移方程和各个时刻的测量值，通过卡尔曼滤波，可以知道目标当前时刻的状态估计值。因此，需要根据状态估计值来确定目标的位置。可以有如下两种测量来估计目标的最终位置。

(1)因为当前时刻的信号估计值S_t为G*1维的向量，根据S_t的G个元素，取其中前M个最大值的元素，把这M个元素的位置，当中在监测区域内M个目标的位置。

(2)求平均法，根据S_t的G个元素，通过聚类算法分成M个类，然后把每个类中所有元素的网格坐标求平均，即可求出M个位置坐标Pⁱ。

式中，i是类的序号，k是类中的第k个元素，Gⁱ为第i类中所有的元素集合，g_k为第k个元素网格的坐标。

通过上面其中一种方法确定的M个目标的位置后，还需要通过关联算法把M个值关联到具体的目标上才会很好的显示目标的轨迹，简单常用关联算法有KNN关联算法。从而实现动态目标的轨迹跟踪。

Claims (4)

1.一种基于网格划分的动态目标波达方向跟踪方法，其特征在于，包括以下几个步骤：

(1)在监测区域建立直角坐标系，在坐标系的一个坐标轴上等间隔距离的布置传感器；

(2)构建远场传感器阵列的窄带信号测量模型及含噪声的源信号模型；

(3)把监测平面环境平均划分网格G等份，每个网格的信号入射角的范围大小相同，即π/G；

(4)在平面环境中构建运动目标合理的转移模型，为了实现卡尔曼滤波跟踪，构建卡尔曼方程，所以这个转移模型应当是线性的，具体步骤如下：

假设目标以一种概率在三个相邻的网格间移动，即目标当前网格，与其相邻的下一个网格，与其相邻的上一个网格，

假设目标以1/3的概率停留在当前网格，或移动到相邻的两个网格，则目标的状态移动方程表示如下：

S(t)＝F*S(t-1)+Q(t) (2)

其中S(t-1)表示前一时刻的信号的值，F是G*G维的转移概率矩阵，Q(t)为当前时刻转移噪声，

其中，f_ij是F的元素，i,j分别是F的行和列；

(5)结合目标的转移模型和信号的测量模型，再根据卡尔曼滤波的原理，估计出的动态目标当前位置。

2.根据权利要求1所述的基于网格划分的动态目标波达方向跟踪方法，其特征在于，所述步骤(2)中构建远场传感器阵列的窄带信号测量模型和含噪声的源信号模型，具体如下：

传感器阵列在直角坐标系的横轴上等间隔距离布置，一般间隔为信号的一半波长，即λ/2，其中λ为信号波长，假设第一个传感器阵元的坐标为(0,0)，则其它阵元的坐标依次为(λ/2,0)、(λ,0)、…、(N*λ/2,0)，传感器的测量模型可表示为：

X(t)＝A*S(t)+R(t) (1)

其中，t为时刻，S(t)为t时刻在区域中的M个目标所发射信号的值，是M*1维向量；R(t)为t时刻N个传感器接收的噪声，是N*1维向量；A为传感器阵列的N*M维的导向矩阵，X(t)为N个传感器接收的信号值，是N*1维向量；

上述中的A＝[a(θ₁),a(θ₂),...,a(θ_M)]^T

其中表达式中的d为传感器位置间隔距离，λ为信号波长，θ_i表示第i个目标到达传感器的角度，i＝1,2,…,M，而α(θ_i)表示N个传感器对第i个目标的导向向量。

3.根据权利要求1所述的基于网格划分的动态目标波达方向跟踪方法，其特征在于，所述步骤(3)中把监测平面环境平均划分网格，使得测量方程为线性的，具体如下：

网格划分是以原点为基点，把坐标轴的上半区域平均划分G等份，即每份的角度范围大小为π/G，划分网格后，式(1)中的A就变成有固定维数的N*G维的导向矩阵，不会根据在监测区域中的目标个数M变化而变化，因此A由如下等式表示：

A＝[a(θ₁),a(θ₂),...,a(θ_G)]^T

其中的θ_i，i＝1,2,…,G在划分网格时是已知的，θ₁＝π/G，θ₂＝2π/G，θ_G-1＝(G-1)π/G，因此也需要相应的改变式(1)中的S(t)，把S(t)变成G*1维的向量，但是S(t)是稀疏的，即在某个网格中有目标，在S(t)对应的位置上才有值，其他的位置上都为0，(1)式中的X(t)和R(t)都将变成G*1维向量。

4.根据权利要求1所述的基于网格划分的动态目标波达方向跟踪方法，其特征在于，所述步骤(5)中结合目标的转移模型和信号的测量模型，再根据卡尔曼滤波的原理，估计出的动态目标当前位置，具体如下：

构建卡尔曼滤波方程

<mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mi>X</mi> <mo>(</mo> <mi>t</mi> <mo>)</mo> <mo>=</mo> <mi>A</mi> <mo>*</mo> <mi>S</mi> <mo>(</mo> <mi>t</mi> <mo>)</mo> <mo>+</mo> <mi>R</mi> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mtd> </mtr> <mtr> <mtd> <mi>S</mi> <mo>(</mo> <mi>t</mi> <mo>)</mo> <mo>=</mo> <mi>F</mi> <mo>*</mo> <mi>S</mi> <mo>(</mo> <mi>t</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> <mo>+</mo> <mi>Q</mi> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mtd> </mtr> </mtable> </mfenced>

简化上式符号，用X_t表示X(t)、S_t表示S(t)、S_t-1表示S(t-1)，根据卡尔曼滤波的计算步骤：

S_t|t-1＝F*S_t-1

P_t|t-1＝F*P_t-1F^T+Q

K＝P_t|t-1*A*(A*P_t|t-1*A^T+R)^-1

P_t＝(E_G*G-K*A)*P_t|t-1

S_t＝S_t|t-1+K*(X_t-A*S_t|t-1)

其中K为卡尔曼增益，F为转移概率矩阵，S_t为当前信号的值，S_t-1为前一时刻信号的值，S_t|t-1为当前时刻的信号预测值，P_t为当前时刻S_t的协方差矩阵，P_t-1为前一时刻S_t-1的协方差矩阵，P_t|t-1为当前时刻预测信号值的协方差矩阵，E_G*G为G*G维的单位矩阵；X_t为当前时刻的信号测量值，Q、R分别为高斯噪声Q(t)和R(t)的协方差；

根据上面的卡尔曼滤波计算步骤得到当前时刻的G个网格中的信号的值，然后根据这G个元素值计算出目标最后的位置，有两种方法：

(1)把计算出的G个元素值中前M个最大信号值所处网格坐标作为M个目标的位置；

(2)根据当前时刻t的状态估计值S_t，根据聚类算法把S_t的G个元素分成M类，即假设监测区域有M个目标，然后求每类所有元素网格坐标的平均值，当作目标所处的位置。