一种基于贝叶斯方法的目标方位跟踪方法与流程

本发明属于阵列信号处理领域。

背景技术：

在雷达、通信、声呐、气象等领域，阵列信号处理有广泛而重要的应用。波达方向角(directionofarrival,doa)估计是阵列信号处理中的一种侧向技术，而实际情况中，目标多处于运动状态，doa跟踪是对运动目标doa进行实时估计，高精度的doa跟踪在个人定位业务、战场移动通信和语音信号处理等领域中都有广泛应用。然而，目标的运动通常导致空间谱扩展，而且运动目标的快拍信号难以长时间积累，这导致传统静止目标的doa估计方法不再适用。因此，高精度的doa跟踪成为了研究重点。

目前，最具有代表性的doa跟踪方法为：子空间类方法和贝叶斯类方法。

(一)子空间类方法。子空间类方法是分段估计出信号波达方向作为瞬时值，通过跟踪或者更新子空间实现doa跟踪。在yangb.anextensionofthepastdalgorithmtobothrankandsubspacetracking[j].signalprocessingletters,ieee,1995,2(9):179-182(一种扩展的实现秩和子空间跟踪的pastd算法).中提出的紧缩投影近似子空间跟踪(projectionapproximationsubspacetrackingdeflation,pastd)算法由于计算复杂度较低而被引入到doa跟踪问题当中。在sanchez-araujoj,marcoss.anefficientpastd-algorithmimplementationformultipledirectionofarrivaltracking[j].ieeetransactionsonsignalprocessing,1999,47(8):2321-2324(一种高效的用于多目标到达角度跟踪的pastd算法).中用pastd算法实时估计信号子空间和噪声子空间，并结合卡尔曼滤波实现doa跟踪。在低信噪比的情况下，这类方法的跟踪精度较差。

(二)贝叶斯类方法。基于贝叶斯跟踪方法是以贝叶斯原理为基础，利用未知参数的先验信息和样本信息得出后验信息，然后根据后验信息进行参数估计，可以在相干信源和低信噪比的情况下跟踪。典型的方法有序列蒙特卡洛方法，即粒子滤波(particlefilter,pf)方法。在ortonm,fitzgeraldw.abayesianapproachtotrackingmultipletargetsusingsensorarraysandparticlefilters[j].signalprocessing,ieeetransactionsonsignalprocessing,2002,50(2):216-223(使用传感阵列和粒子滤波进行多目标跟踪的贝叶斯方法).中，将粒子滤波算法用于一维阵列的doa跟踪。但是粒子滤波算法需要用到状态转移概率和观测噪声概率函数，在这些条件未知的情况下，可能出现较大误差。

技术实现要素：

本发明是为了解决低信噪比情况下，基于空间类doa跟踪方法精度差，以及粒子滤波doa跟踪方法需要的已知条件过多的问题，本发明提供了一种基于贝叶斯方法的目标方位跟踪方法。

一种基于贝叶斯方法的目标方位跟踪方法，它包括如下步骤：

步骤一：初始化当前t时刻的信号到达角θk(t)、信号功率pk(t)、噪声功率σ²(t)和角度变化量△θk(t)，使得θk(t)＝θk(t-1)、pk(t)＝pk(t-1)、σ²(t)＝σ²(t-1)、△θk(t)＝0；

其中，θk(t-1)、pk(t-1)和σ²(t-1)均已知，k表示第k个信号源，且k＝1,2,3,...,k，k为空间信号源个数，θk(t-1)为t-1时刻的信号到达角，pk(t-1)为t-1时刻的信号功率，σ²(t-1)为t-1时刻的噪声功率，角度变化量△θk(t)为当前时刻相对于上一时刻的角度变化量；

其中，t＝1,2,3,...,t′-1,t′，t′为总观测时间；

步骤二：在θk(t-1)处，进行一阶泰勒展开，从而获得阵列流形矩阵a；

步骤三：采用贝叶斯方法对步骤一中的信号功率pk(t)、噪声功率σ²(t)和阵列流形矩阵a进行处理，获得后验均值μ和后验方差∑；

步骤四：根据后验均值μ和后验方差∑更新信号功率pk(t)、噪声功率σ²(t)和角度变化量△θk(t)；

步骤五：采用em算法追踪m次更新后的信号功率pk(t)、噪声功率σ²(t)和角度变化量△θk(t)，判断m次更新后的信号功率pk(t)、噪声功率σ²(t)和角度变化量△θk(t)的数值是否收敛，结果为是，执行步骤六，结果为否，令m＝m+1，执行步骤一；

其中，m为重复执行步骤一至步骤四的次数；

步骤六：收敛状态下，第m次更新后的信号功率pk(t)为当前t时刻信号功率pk(t)的估计值，更新后的噪声功率σ²(t)为当前t时刻噪声功率σ²(t)的估计值，更新后的角度变化量△θk(t)为当前t时刻角度变化量θk(t-1)的估计值，

根据当前t时刻角度变化量△θk(t)的估计值，获得当前t时刻的信号到达角θk(t)的估计值，其中，θk(t)＝θk(t-1)+△θk(t)；

步骤七：判断t是否等于t′，结果为否，令t＝t+1，执行步骤一；结果为是，执行步骤八；

步骤八：完成了总观测时间t′内，各时刻的信号功率pk(t)、噪声功率σ²(t)和信号到达角θk(t)的估计，从而实现对目标方位的实时追踪。

所述的步骤二中，在θk(t-1)处，进行一阶泰勒展开，从而获得阵列流形矩阵a的具体过程为：

步骤二一：根据泰勒展开式，将t时刻的流型向量a(θk(t))，在θk(t-1)处进行一阶泰勒展开得到：

a(θk(t))＝a(θk(t-1))+a'(θk(t-1))△θk(t)(公式一)，

其中，

a(θk(t-1))为t-1时刻的流型向量，a′(θk(t-1))表示t-1时刻的流型向量的导数，e为自然常数，j为虚数单位，d为阵元间距，λ为信号波长，m为阵元个数，[]^t表示矩阵转置；

步骤二二：根据流型向量a(θk(t))，获得t时刻的阵列流形矩阵a；

a＝a1+bδ(公式二)，

其中，

a1＝[a(θ1(t-1)),a(θ2(t-1)),a(θ3(t-1)),...,a(θk-1(t-1)),a(θk(t-1))]，

b＝[a'(θ1(t-1)),a'(θ2(t-1)),a'(θ3(t-1)),...,a'(θk-1(t-1)),a'(θk(t-1))]，

δ＝diag([△θ1(t),△θ2(t),△θ3(t),...,△θk-1(t),△θk(t)])，

a1为t-1时刻的阵列流形矩阵，b为t-1时刻流形向量的导数构成的矩阵，δ为角度校正量，diag()表示对角矩阵函数。

所述的步骤三中，采用贝叶斯方法对步骤一中的信号功率pk(t)、噪声功率σ²(t)和阵列流形矩阵a进行处理，获得后验均值μ和后验方差∑的具体过程为：

μ＝pa^h(σ²(t)i+apa^h)^-1x(公式三)，

∑＝(p^-1+σ^-2(t)a^ha)^-1(公式四)，

其中，

p＝diag([p1(t),p2(t),p3(t),...,pk-1(t),pk(t)])，

p为信号的协方差矩阵，diag()表示对角矩阵函数，i为单位矩阵。

所述的步骤四中，根据后验均值μ和后验方差∑更新信号功率pk(t)、噪声功率σ²(t)和角度变化量△θk(t)的具体过程为：

[△θ1(t),△θ2(t),△θ3(t),...,△θk-1(t),△θk(t)]^t＝u^-1v(公式七)，

其中，

(∑)(n,n)表示∑的第n行第n列元素，

(μ)n表示μ的第n行，n为正整数；

x＝as+n，

其中，上标h表示共轭转置，表示取实部，⊙表示哈达玛乘积，l表示快拍数，diag()表示表示对角矩阵函数，v和u均表示中间变量，x∈c^m×l为阵列接收数据，s∈c^k×l为入射信号，n∈c^m×l为观测噪声，c为复数集合，tr()表示矩阵的迹。

本发明考虑的是角度慢变模式下的doa跟踪，即相邻时刻信号到达角度值及信号功率的变化很小。本发明在信号到达角和功率缓慢变化的情形下，根据前一时刻的信号到达角、信号功率和噪声功率的估计值，以及当前时刻的观测值，利用em算法，对前一时刻的信息进行校正，进而实现对信号到达角度的追踪。

本发明带来的有益效果是，

在本发明中，利用泰勒展开公式，将角度慢变模式下的doa跟踪建模为一个动态模型，并基于贝叶斯理论，将角度的跟踪转化为概率模型中的参数估计问题，采用em算法来求解。本发明提出的跟踪方法的优势在于能够同时实时跟踪信号的到达角度、信号功率以及噪声功率，且无需知道状态转移概率等先验信息，而且与稀疏贝叶斯学习方法相比，计算量小。

本发明所述一种基于贝叶斯方法的目标方位跟踪方法与子空间类doa跟踪方法相比，本发明跟踪方法精度提高30％以上，且在低信噪比的情况下尤为明显。

本发明所述一种基于贝叶斯方法的目标方位跟踪方法与粒子滤波doa跟踪方法相比，本发明所述方法无需知道状态转移概率和观测噪声概率函数。

本发明在进行doa跟踪的同时，本发明还可以实现信号功率和噪声功率的多参数联合跟踪。

附图说明

图1为本发明所述的一种基于贝叶斯方法的目标方位跟踪方法的流程图；

图2为有两个入射信号情形下的由本发明方法得到的信号到达角(doa)追踪图；

图3为利用本发明所述的一种基于贝叶斯方法的目标方位跟踪方法，得到的信号功率追踪图；

图4为利用本发明所述的一种基于贝叶斯方法的目标方位跟踪方法，得到的噪声功率追踪图；

图5为本发明方法与pastd-kalman(紧缩投影近似子空间跟踪-卡尔曼滤波)算法进行doa跟踪的均方根误差与信噪比关系曲线的对比图。

具体实施方式

具体实施方式一：参见图1说明本实施方式，本实施方式所述的一种基于贝叶斯方法的目标方位跟踪方法，它包括如下步骤：

步骤一：初始化当前t时刻的信号到达角θk(t)、信号功率pk(t)、噪声功率σ²(t)和角度变化量△θk(t)，使得θk(t)＝θk(t-1)、pk(t)＝pk(t-1)、σ²(t)＝σ²(t-1)、△θk(t)＝0；

其中，θk(t-1)、pk(t-1)和σ²(t-1)均已知，k表示第k个信号源，且k＝1,2,3,...,k，k为空间信号源个数，θk(t-1)为t-1时刻的信号到达角，pk(t-1)为t-1时刻的信号功率，σ²(t-1)为t-1时刻的噪声功率，角度变化量△θk(t)为当前时刻相对于上一时刻的角度变化量；

其中，t＝1,2,3,...,t′-1,t′，t′为总观测时间；

步骤二：在θk(t-1)处，进行一阶泰勒展开，从而获得阵列流形矩阵a；

步骤三：采用贝叶斯方法对步骤一中的信号功率pk(t)、噪声功率σ²(t)和阵列流形矩阵a进行处理，获得后验均值μ和后验方差∑；

步骤四：根据后验均值μ和后验方差∑更新信号功率pk(t)、噪声功率σ²(t)和角度变化量△θk(t)；

步骤五：采用em算法追踪m次更新后的信号功率pk(t)、噪声功率σ²(t)和角度变化量△θk(t)，判断m次更新后的信号功率pk(t)、噪声功率σ²(t)和角度变化量△θk(t)的数值是否收敛，结果为是，执行步骤六，结果为否，令m＝m+1，执行步骤一；

其中，m为重复执行步骤一至步骤四的次数；

步骤六：收敛状态下，第m次更新后的信号功率pk(t)为当前t时刻信号功率pk(t)的估计值，更新后的噪声功率σ²(t)为当前t时刻噪声功率σ²(t)的估计值，更新后的角度变化量△θk(t)为当前t时刻角度变化量θk(t-1)的估计值，

根据当前t时刻角度变化量△θk(t)的估计值，获得当前t时刻的信号到达角θk(t)的估计值，其中，θk(t)＝θk(t-1)+△θk(t)；

步骤七：判断t是否等于t′，结果为否，令t＝t+1，执行步骤一；结果为是，执行步骤八；

步骤八：完成了总观测时间t′内，各时刻的信号功率pk(t)、噪声功率σ²(t)和信号到达角θk(t)的估计，从而实现对目标方位的实时追踪。

本实施方式，在本发明中，利用泰勒展开公式，将角度慢变模式下的doa跟踪建模为一个动态模型，并基于贝叶斯理论，将角度的跟踪转化为概率模型中的参数估计问题，采用em算法来求解。本发明提出的跟踪方法的优势在于能够同时实时跟踪信号的到达角度、信号功率以及噪声功率，且无需知道状态转移概率等先验信息，而且与稀疏贝叶斯学习方法相比，计算量小。

em算法的英文全称为：expectationmaximization。

原理分析：首先根据角度慢变的假设，根据一阶泰勒展开式，得到与上一时刻角度和角度更新量有关的阵列流形矩阵，以便有效利用上一时刻的信息。然后进行em迭代，每次迭代中根据上一轮求得的信号功率、噪声功率及角度校正量的估计值计算信号的后验概率密度，根据该后验概率密度，对信号与噪声的联合概率密度求均值并使之最大化从而更新参数的估计值。获得角度校正量后，即可计算出实时的到达角度。

具体实施方式二：参见图1说明本实施方式，本实施方式与具体实施方式一所述的一种基于贝叶斯方法的目标方位跟踪方法的区别在于，所述的步骤二中，在θk(t-1)处，进行一阶泰勒展开，从而获得阵列流形矩阵a的具体过程为：

步骤二一：根据泰勒展开式，将t时刻的流型向量a(θk(t))，在θk(t-1)处进行一阶泰勒展开得到：

a(θk(t))＝a(θk(t-1))+a'(θk(t-1))△θk(t)(公式一)，

其中，

a(θk(t-1))为t-1时刻的流型向量，a′(θk(t-1))表示t-1时刻的流型向量的导数，e为自然常数，j为虚数单位，d为阵元间距，λ为信号波长，m为阵元个数，[]^t表示矩阵转置；

步骤二二：根据流型向量a(θk(t))，获得t时刻的阵列流形矩阵a；

a＝a1+bδ(公式二)，

其中，

a1＝[a(θ1(t-1)),a(θ2(t-1)),a(θ3(t-1)),...,a(θk-1(t-1)),a(θk(t-1))]，

b＝[a'(θ1(t-1)),a'(θ2(t-1)),a'(θ3(t-1)),...,a'(θk-1(t-1)),a'(θk(t-1))]，

δ＝diag([△θ1(t),△θ2(t),△θ3(t),...,△θk-1(t),△θk(t)])，

a1为t-1时刻的阵列流形矩阵，b为t-1时刻流形向量的导数构成的矩阵，δ为角度校正量，diag()表示对角矩阵函数。

具体实施方式三：参见图1说明本实施方式，本实施方式与具体实施方式一所述的一种基于贝叶斯方法的目标方位跟踪方法的区别在于，所述的步骤三中，采用贝叶斯方法对步骤一中的信号功率pk(t)、噪声功率σ²(t)和阵列流形矩阵a进行处理，获得后验均值μ和后验方差∑的具体过程为：

μ＝pa^h(σ²(t)i+apa^h)^-1x(公式三)，

∑＝(p^-1+σ^-2(t)a^ha)^-1(公式四)，

其中，

p＝diag([p1(t),p2(t),p3(t),...,pk-1(t),pk(t)])，

p为信号的协方差矩阵，diag()表示对角矩阵函数，i为单位矩阵。

具体实施方式四：参见图1说明本实施方式，本实施方式与具体实施方式一所述的一种基于贝叶斯方法的目标方位跟踪方法的区别在于，所述的步骤四中，根据后验均值μ和后验方差∑更新信号功率pk(t)、噪声功率σ²(t)和角度变化量△θk(t)的具体过程为：

[△θ1(t),△θ2(t),△θ3(t),...,△θk-1(t),△θk(t)]^t＝u^-1v(公式七)，

其中，

(∑)(n,n)表示∑的第n行第n列元素，

(μ)n表示μ的第n行，n为正整数；

x＝as+n，

其中，上标h表示共轭转置，表示取实部，⊙表示哈达玛乘积，l表示快拍数，diag()表示表示对角矩阵函数，v和u均表示中间变量，x∈c^m×l为阵列接收数据，s∈c^k×l为入射信号，n∈c^m×l为观测噪声，c为复数集合，tr()表示矩阵的迹。

具体实施方式五：参见图1说明本实施方式，本实施方式与具体实施方式二所述的一种基于贝叶斯方法的目标方位跟踪方法的区别在于，所述的步骤二二中，a＝[a(θ1(t)),a(θ2(t)),a(θ3(t)),...,a(θk-1(t)),a(θk(t))]。

仿真验证：

利用本发明所述的一种基于贝叶斯方法的目标方位跟踪方法进行仿真试验：

步骤一一，仿真参数设置如下：阵元个数m＝12，快拍数l＝20，阵元间距d＝λ/2，噪声功率σ²(t)＝1，观测时间t′＝100；

步骤一二：建立仿真信号运动状态模型：

xk(t)＝fxk(t-1)+gvk(t)

式中，vk(t)为零均值的高斯白噪声，且方差为0.001，其中表示θk(t)的速度；

步骤一三：设置入射信号个数k＝2，初始角度分别为40°和30°，信噪比的变化分别为8db至-3db和3db至8db，得到的doa、信号功率和噪声功率跟踪曲线分别如图2、图3和图4所示，其中，doa为信号到达角；

步骤一四：为了更好地分析本发明方法的doa跟踪性能，将其与pastd算法结合卡尔曼滤波的doa跟踪方法进行比较。进行100次蒙特卡洛试验，得到两种算法的均方根误差随信噪比的关系曲线如图5所示。

仿真结果证明：

由图2、图3和图4可以看出，本发明的方法可以很好的进行信号到达角度、信号功率和噪声功率的实时跟踪；由图5可以看出，本发明的方法与pastd算法结合卡尔曼滤波的doa跟踪方法相比具有更高的跟踪精度，在低信噪比的情况下优势更加明显。