一种针对机动目标的无人机鲁棒化跟踪方法与流程

本发明涉及一种鲁棒化的机动目标跟踪方法，属于联合递推贝叶斯估计理论和多源信息融合理论以增强无人机进行机动目标跟踪鲁棒性的技术领域。

背景技术：

无人机(unmannedaerialvehicle,uav)是一种在无载人驾驶情况下由自身动力系统驱动、以自主或受遥控方式完成特定动作与任务的飞行器。uav不仅在地形侦察、情报搜集、目标打击等军用领域大显身手，而且在地理测绘、防火监视、搜索救援等民用领域也广泛存在。随着材料科技、制造技术和计算机技术、通信技术等的迅速发展，uav本身也逐渐趋于微型化、轻量化、低成本和低噪声。总体来讲，uav的应用类型大致可划分为监视型、打击型和通信型，在大多数这三类属性的应用场合中，uav的任务都包含了目标跟踪这一关键环节，这一环节执行的效果好坏对后续的计算机系统的决策和飞行器任务的执行都有着较大的影响。

从本质上看，uav进行目标跟踪可以看作是一种移动传感器对目标的跟踪行为。一方面，区别于静态传感器进行的目标跟踪(被动型)，uav目标跟踪属于一种主动型目标跟踪，在跟踪方式的强度上要大大提升；另一方面，大多数静态传感器的目标跟踪单独存在，而uav目标跟踪则通常是uav整体任务的一个环节，它对跟踪就要求有更好的性能以维持整个系统的良好状态。

由于uav出现的初衷是为了避免人力浪费和保障人生安全，因此大部分uav进行目标跟踪常在复杂环境和约束条件下进行，这就要求其计算中心对实时信息的处理应具有更好的鲁棒性，以应对系统本身或观测过程中的突发情况和各类扰动。首先对uav配备的传感器装置而言，无论是具有逃逸策略的目标的电子对抗，还是恶劣天气抑或是复杂地形对无线电传播的干扰，都会使uav的传感器所探测到的观测信息的可靠性大大降低。其次，对uav本身系统而言，由于uav的飞行器属性，它自身便具有爬升、旋转、俯冲等强机动行为，这些潜在机动行为不仅会使配备的陀螺仪的输出误差随角速率的巨大变化而增大，进而给整个惯性系统的定位估计带来较大误差，同时会对机载传感器的量测过程造成较大干扰。再者，对目标而言，当被跟踪目标同样是飞行器，或者当地面灵活目标觉察到被侦察时选择迅速逃逸，目标本身也会具有uav的上述机动特性，这种机动性不仅会使计算中心对目标状态变迁服从高斯过程的假设失效，也会使计算中心内置默认的跟踪方法和目标的真实动力学特性出现模型失配现象，从而出现估计精度大幅度下滑，甚至发散以至于uav完全不能跟踪目标。

技术实现要素：

本发明针对uav在复杂环境下进行单目标跟踪时所遇到的实际问题，公开了一种鲁棒化的机动目标跟踪方法。该发明建立了跳变马尔可夫模型以描述uav进行机动目标跟踪时可能出现的模型失配问题，模型具有重尾非高斯属性的学生t分布替代传统高斯分布，以描述带有状态突变的目标系统过程和异常测量的uav传感器的观测过程，在交互式多模型框架下，按照递推贝叶斯估计方法设计了基于学生t分布概率密度函数的子模型广义非线性滤波方法，最后将各个子模型的鲁棒化估计结果进行信息融合以得到机动目标的跟踪结果。

该方法包括以下几个核心步骤：

步骤一：建立uav进行单机动目标跟踪时的学生t分布跳变马尔可夫模型；

步骤二：计算各个目标动力学子模型间的混合概率并基于学生t分布进行交互；

步骤三：对交互后的服从学生t分布的系统状态进行预测；

步骤四：对服从学生t分布的系统状态进行量测更新；

步骤五：基于学生t分布对子模型的模型概率进行更新；

步骤六：融合各个子模型的鲁棒估计结果并进行输出；

步骤七：再回到步骤三重新依次执行，直至目标跟踪结束。

进一步地，步骤一所述的学生t分布跳变马尔可夫模型包括系统状态转移方程xk＝f(xk-1,rk)+wk-1(rk)和传感器观测方程zk＝h(xk)+vk，其中下标k代表数据来源第k个离散采样时刻，zk为观测量，xk为目标状态变量，f(·)和h(·)分别为目标的系统函数和uav的传感器观测函数，rk为模型变量；

用p(·)代表概率密度函数，std(·)为学生t分布函数符号，过程噪声wk-1(rk)和观测噪声vk分别服从学生t分布，即p(wk-1(rk))＝std(wk|0,qk-1(rk),υw)和p(vk)＝std(vk|0,rk,υv)，其中qk-1(rk)和rk分别为过程噪声和观测噪声的尺度化矩阵，υw和υv分别为过程噪声和观测噪声的自由度参数；

此外，认为任意第r个子模型的状态的k时刻模型先验(即k-1时刻的模型后验)概率密度函数和k时刻模型后验概率密度函数p(xk|z1:k)＝std(xk|mk,pk,υ)均为学生t分布函数，其中和mk、和pk、υ^r和υ分别为各自学生t分布的均值、尺度化矩阵和自由度参数。

进一步地，步骤二包含以下子步骤：

2.1计算混合概率它代表k-1时刻子模型为模型r而k时刻为模型s的概率，可按下式计算：

其中πrs表征在k-1时刻到k时刻由模型r切换为模型s的概率，为第r个子模型在k时刻的模型先验概率(即k-1时刻的模型后验概率)，为一个归一化常量；

2.2计算混合的学生t分布的均值尺度化矩阵和自由度参数υ^0s分别如下：

υ^0s＝υ^r(4)

进一步地，步骤三需要计算学生t分布中均值的预测值尺度化矩阵的预测值和自由度参数的预测值υ^s，具体如下：

υ^s＝υ^0s(7)

进一步地，步骤四包含如下子步骤：

4.1计算第s个子模型的边缘化观测似然函数的均值自相关尺度化矩阵和互相关尺度化矩阵具体为：

4.2计算滤波增益和自由度参数即

4.3计算第s个子模型所对应的模型后验概率密度函数p(xk|rk＝s,z1:k)的均值和尺度化矩阵分别如下：

4.4用矩信息匹配法优化均值和尺度化矩阵得到新的均值和尺度化矩阵分别如下：

进一步地，步骤五更新子模型的模型后验概率

其中，似然函数的求取可通过计算第s个子模型的边缘化观测似然函数p^s(zk|z1:k-1)得到，即只需计算这一标准学生t分布概率密度函数。

进一步地，步骤六中，利用更新的模型后验概率与步骤四中的均值和尺度化矩阵可以得到k时刻融合的均值mk和尺度化矩阵pk如下：

本发明的有益效果是：

(1)解决uav跟踪机动目标时存在的模型失配问题，使估计器在目标发生机动行为后依然能快速而准确地跟踪上目标；

(2)较好应对因目标状态突变、uav观测异常等扰动引起的目标系统过程和传感器观测过程呈现的重尾非高斯特性问题，从而全面增强目标跟踪的鲁棒性；

(3)具有较好的一般性和可扩展性，可根据实际工作环境和适用情形进行适当修正，以得到特定条件下最优的跟踪方法，比如，目标有无机动行为可能、系统或观测模型非线性度的强弱、系统过程或观测过程发生突变或异常的可能性大小、多无人机或多目标等情形。

附图说明

图1是本发明方法的流程图。

图2是本发明中uav进行目标跟踪的场景示意图。

图3是投影在x-y平面内的uav目标跟踪效果图。

具体实施方式

下面将就具体实施对本发明作进一步的详细说明，方法整体流程可参考图1。

步骤一：建立uav进行单机动目标跟踪时的学生t分布跳变马尔可夫模型；

图2描述了本发明的一个实施例的工作场景设置，具体参数的设定和场景的假设并不限于该图描述的设置，该图展示了uav和运动目标在相邻两个采样时刻(k-1和k时刻)的空间位置信息情况。如图所示，一个uav在建有x-y-z坐标系的立体空间中飞行并对运动目标进行跟踪，被跟踪目标为一个地面目标，可认为其纵向高度不发生变化且高度为0，即目标仅在x-y平面上进行运动，此时uav也可保持在某一固定高度为h的平面上进行巡航和监视，无需发生纵向的位移。uav机身配有图像传感器和激光测距仪分别进行目标的测向和测距，gps模块和陀螺仪组合进行uav自身的定位、测速和测姿，默认在进行目标跟踪前uav精确已知自身在空间坐标系中的位置下标k代表数据来源于第k个离散采样时刻，分别代表uav的x方向位置、y方向位置、z方向位置，其中z方向位置即uav的高度可认为恒定不变。定义任意k时刻目标的状态变量为列向量目标状态维数n＝6，其中分别代表该时刻目标的x方向位置、y方向位置、x方向速度、y方向速度、x方向加速度、y方向加速度，此处默认目标z方向位置z方向速度z方向加速度进而建立uav跟踪机动目标的跳变马尔可夫模型如下：

xk＝f(xk-1,rk)+wk-1(rk)(1)

(1)为目标状态转移方程，它引入了模型变量rk，rk为目标动力学模型编号，认为目标状态变化规律可在匀速模型、匀加速模型、协调转弯模型、当前统计模型四个模型之间自由切换，即rk在有限集合δ＝{1，2，…,m}m＝4中取任意值，其切换规律可由一个维数为m×m的转移概率矩阵π＝[πrs]m×m(r,s∈δ)进行描述，矩阵中元素πrs＝pr(rk＝s|rk-1＝r)表征在k-1时刻到k时刻由模型r切换为模型s的概率，符号pr表示事件发生的概率。本发明中rk未出现在传感器观测方程(2)中，原因是uav的测向测距组合观测并不随目标运动模式的变化而发生改变，故可以将rk省去。f(·)为目标的系统函数，可以为任意运动目标的动力学函数，如前所述，本例限制其在匀速模型、匀加速模型、协调转弯模型、当前统计模型四个模型(模型数m＝4)之中切换，zk为uav获得的观测量，观测维数d＝2，如(2)中所示，h(·)在本例的实施中为一个测距测向组合的观测函数，wk-1(rk)和vk分别是具有重尾非高斯属性的过程噪声和观测噪声，用学生t分布对它们进行建模：

p(wk-1(rk))＝std(wk|0,qk-1(rk),υw)(3)

p(vk)＝std(vk|0,rk,υv)(4)

式中，p(·)为对应变量的概率密度函数，std(·)为学生t分布的符号，其中qk-1(rk)和rk分别为过程噪声和观测噪声的尺度化矩阵，υw和υv分别为过程噪声和观测噪声的自由度参数。自由度参数是学生t分布区别于高斯分布的重要特征之一，它控制了概率密度函数的峰态，因此可以描述因状态突变或观测异常等干扰导致的相关变量的概率密度函数呈现的重尾特性，这是高斯分布无法做到的。事实上，高斯分布是学生t分布在自由度参数趋于正无穷大时的一个特例。此外，为使后续表述方便，定义

进一步，由全概率公式可得k时刻目标状态的后验概率密度函数p(xk|z1:k)：

其中z1:k表示从第1个采样时刻到第k时刻的所有观测量。

根据(3)和(4)的设定，不难看出每个子模型的后验概率密度函数p(xk|rk＝s,z1:k)是一个学生t分布，因此p(xk|z1:k)可以看作是m个学生t分布的加权求和的结果，本发明用一个学生t分布std(xk|mk,pk,υ)来近似这一求和结果，即认为p(xk|z1:k)＝std(xk|mk,pk,υ)。以下步骤将阐述如何求解目标状态后验概率密度函数p(xk|z1:k)。

步骤二：计算各个目标动力学子模型间的混合概率并基于学生t分布进行交互；

由(5)可以看出，随着目标跟踪的进行，模型个数m的存在会使p(xk|z1:k)的计算量呈指数增长，于是考虑在交互式多模型框架下解决此问题。若k-1时刻所选的模型rk-1＝r(r∈δ)，则k时刻状态的模型先验概率密度函数p(xk-1|rk-1＝r,z1:k-1)表示为：

将各个子模型的模型概率密度函数进行交互处理，得到模型rk＝s(s∈δ)的模型概率密度函数p(xk-1|rk＝s,z1:k-1)：

其中定义为混合概率，它表征k-1时刻子模型为模型r而k时刻为模型s的概率，它可由贝叶斯公式计算得到：

上式中的分母为一个归一化常量。

由(7)可以看出，p(xk-1|rk＝s,z1:k-1)实际上相当于是m个学生t分布的加权求和，为了方便近似计算，此处仅用一个学生t分布来代替求和结果。为了尽可能保留原有m个学生t分布的信息，现用矩信息匹配法求取p(xk-1|rk＝s,z1:k-1)。令则混合后学生t分布的均值和尺度矩阵可计算如下：

混合后的学生t分布自由度参数υ^0s仍然沿用混合前的各个子模型的自由度参数υ^r，即：

υ^0s＝υ^r(11)

步骤三：对交互后的服从学生t分布的系统状态进行预测；

利用(9)-(11)对混合后的系统状态进行预测，预测的概率密度函数p(xk|rk＝s,z1:k-1)为：

其中，学生t分布中均值的预测值和尺度化矩阵的预测值可分别计算如下：

上式(14)中自由度参数υ^0s的预测值υ^s保持不变，即：

υ^s＝υ^0s(15)

步骤四：对服从学生t分布的系统状态进行量测更新；

根据递推贝叶斯估计公式，给出在k时刻获取观测量zk后模型s的模型后验概率密度函数p(xk|rk＝s,z1:k)：

为使后续表述简洁，给出同样服从学生t分布的模型s的边缘化观测似然函数该学生t分布的均值为：

自相关尺度化矩阵为：

互相关尺度化矩阵为：

经整理计算可得p(xk|rk＝s,z1:k)的均值

p(xk|rk＝s,z1:k)的尺度化矩阵

其中(21)中的滤波增益

自由度参数

其中，d为步骤一中给出的观测维数。注意(23)的更新法则会使自由度参数υ^s在每一个更新时刻递增，这样随着时间的推移和跟踪过程的进行，学生t分布就将失去它特有的重尾属性，从而失去对目标状态突变值或观测异常值等干扰的概率描述，进而影响跟踪方法的鲁棒性。为了保持跟踪方法较好的鲁棒性，自由度参数需处于一个相对较小的值。于是重新令先用自由度参数υ^s代替再利用矩信息匹配法优化均值和尺度化矩阵可得新的均值和尺度化矩阵

步骤五：基于学生t分布对子模型的模型概率进行更新；

再次利用贝叶斯公式，得到更新的模型s的模型后验概率

其中，该式中的似然函数就是模型s的边缘化观测似然函数p^s(zk|z1:k-1)，所以有：

要计算(27)，只需将υ^s代入学生t分布的概率密度函数表达式，具体如下：

其中γ(·)为伽玛函数(gammafunction)，π为圆周率，|·|为行列式运算符。

步骤六：融合各个子模型的鲁棒估计结果并进行输出；

为得到(5)中p(xk|z1:k)＝std(xk|mk,pk,υ)的最终融合估计结果，利用(29)和(30)分别融合在步骤四中得到的各个子模型估计的状态均值和尺度化矩阵融合后的状态均值mk和尺度化矩阵pk分别为：

为保持较好的鲁棒性，融合后的自由度参数可取为各个子模型的自由度参数的最小值：

υ＝min{υ¹,υ²,…,υ^s,…,υ^m}(31)

其中min{·}为最小值符号函数。

步骤七：再回到步骤三重新依次执行，直至目标跟踪结束。

若还需继续跟踪目标，则令k＝k+1进入下一采样时刻，回到步骤三并再依次顺序执行，反复以上过程直至目标跟踪结束。

需要指出的是，本发明具有较好的一般性和可扩展性。具体可解释如下：

1.可通过调节子模型个数m以获得更好的跟踪效果。若目标运动机动性和不可预测性较强，可增加m值以包含更多运动模型；若目标运动仅在固定运动模式切换，只需引入几个典型子模型即可。更特殊地，若目标确认无机动性，令m＝1时，此时仅利用步骤三和步骤四便可以完成目标跟踪。

2.当目标系统函数f(·)和观测函数h(·)非线性程度较弱时，涉及非线性函数积分运算的(13)、(14)、(17)、(18)和(19)可用泰勒级数展开进行近似求取；当非线性程度较为强烈，可考虑采用ut变换、球面-相径容积规则等数值方法进行积分运算。

3.若单个uav的观测方式发生变化，比如某一个时刻因故障或跟踪策略调整，只有测距信息而无测角信息，或者当有多个uav协同工作进行某一机动目标的跟踪，但并非所有uav都配有相同传感器而采用同一种观测方式，这些情形均可认为目标跟踪过程中还存在观测模型的切换问题，只需在(2)中对h(·)或vk添加对应的模型变量，联合目标运动的多模型一起构造马尔可夫链即可用本发明公开的方法进行后续跟踪处理。更特殊地，当目标为非机动目标时，只需考虑观测模型的切换而无需考虑目标运动的多模型。

本发明一个实施例的仿真效果见图3，该图为uav跟踪机动目标在x-y平面上投影的轨迹图。uav投影在x-y平面初始位置为(2000米，600米)，机动目标始于x-y平面的(2000米，1000米)处。机动目标出发后不久先以10度/秒的角速度大幅度顺时针方向转动，后立即以10度/秒的角速度持续作逆时针方向转动。设定整个目标状态变化过程和uav观测过程受到一定概率异常值的干扰。由图可见，面向非机动目标的鲁棒跟踪方法的轨迹曲线在目标发生第一次大幅度机动行为后便不能再跟踪上真实目标，轨迹随后处于发散状态，其跟踪性能也不再可靠，这个结果是由于模型失配问题导致的。面向机动目标的非鲁棒化跟踪方法相比面向非机动目标的鲁棒跟踪方法具有更高的可信度，它基本能跟踪上真实目标的运动轨迹，但它在目标发生机动行为时段的误差要明显大于本发明公开的目标跟踪方法，原因在于该方法设计并未考虑目标系统过程和观测过程中异常值对估计性能的影响，仅仅使用高斯分布进行设计会使跟踪方法鲁棒性变差。对照真实目标轨迹和本发明对应的轨迹，本发明公开的方法在全跟踪时段(包括目标发生机动行为的时间段)都能以相对较小的误差跟踪上真实目标，既能改善机动目标本身带来的模型失配问题，还具有较好的鲁棒性以应对系统过程和观测过程的异常发生，因此具有较为稳健和良好的跟踪效果。