2. 武汉理工大学 智能交通系统研究中心, 湖北 武汉 430063;
3. 水路公路交通安全控制与装备教育部工程研究中心, 湖北 武汉 430063
2. Intelligent Transportation System Research Center, Wuhan University of Technology, Wuhan 430063, China;
3. Engineering Research Center for Transportation Safety, Ministry of Education, Wuhan 430063, China
随着内河水运信息化的飞速发展,为了保证航运业安全有效地运行,雷达被广泛运用于追踪船舶运行状态方面,从而避免船舶发生碰撞造成不必要的经济损失。特别是长江干线,现在已经成为了世界上运输量最大并且最繁忙的通航河流。因此为了防止航行事故的发生,对船舶进行交通管制显得尤为重要。普通船用雷达一般通过人工标绘来获取船舶的航速、航向等航行数据。随着技术的不断发展,现在的船用雷达普遍具有自动标绘功能(automatic radar plotting aid,ARPA),可以自动获取上述数据,对运动的目标进行提取与跟踪。然而,近年来,随着跨临河建筑物的快速增加、船舶交通流量的不断增加,雷达所采集到的ARPA目标数据含有大量噪声,给目标跟踪带来了很大干扰。如何去除这些干扰,以达到更好的跟踪效果,成为了一个亟需解决的问题。
目标跟踪的含义是利用远程测量数据,通过使用数据关联算法和Kalman滤波器来估计运动目标的状态[1]。根据跟踪算法所采用的结构,可以采用不同类型关联程序[2-5],一旦数据关联完成,线性的或扩展的卡尔曼滤波器便对目标的位置和速度给出进一步的估计值。雷达数据处理过程中,数据关联是目标跟踪的核心部分。文献[2]通过采用Kalman滤波器与最近邻算法相结合实现了对车辆轨迹的追踪。文献[6]通过将聚类分析方法运用在概率数据关联算法中,以降低目标在复杂环境下的关联错误率,提高了目标量测的稳定性。最近邻算法[7]原理简单,计算量小;概率关联算法[8]与标准卡尔曼滤波在最大存储量上几乎相等,甚至基本不变,因此其算法容易实现。但在目标密度较大时,这两种方法都容易出现丢失目标或误跟目标的问题。
本文根据获取的APRP目标,提出一种基于目标的速度和航向对目标数据进行关联,去除真实目标中的噪声干扰,并将速度和航向作为BRB(Belief rule-based)系统的前置属性,建立置信规则库,通过ER算法来演绎推理,以此验证该方法的可靠性,最后采用Kalman滤波实现对目标的跟踪算法。
1 基于速度和航向的目标航迹检测本文所采集的目标数据均为ARPA雷达返回的雷达数据。虽然ARPA设备具备一定的自动识别目标能力,但是在内河水域中它会将礁石、桥梁等目标都识别为船只。为此,在滤波之前,必须将这些错误的目标点剔除。本文提出将选取的目标点先进行速度和航向的检测,如果满足条件,则将此目标点作为录取点来进行Kalman滤波,从而形成一个航迹上的跟踪点。
1) 速度检测。常见的船舶速度值为0~20m/s,其转弯速度或加速度为2(°)/s。因此,在选取目标点时,首先对其速度进行检测,如果目标点的速度太小或超过20m/s,则认为此目标点为伪目标点,并将其舍去,只有速度在最大速度和最小速度之间的目标点,才给予录取。但是仅仅依靠速度上的检测是不够的。如图 1所示,满足速度检测条件的点迹有A、B。
根据最近邻域准则A点离预测目标点最近,因此A点应该是所要关联的点迹。但因为船舶不可能在很短的时间内做出接近90°的转弯到达A点,因而A点不是真实目标点迹,而相对较远的B点才是所要选取的目标点。因此还需要对目标点迹做出符合航向上实际转向的关联,也就是航向检测。
2) 航向检测。当观测点符合速度检测后,再对其进行航向检测,如果观测点也符合航向检测,最后将符合两者条件的点作为最终可以相关联的目标点,进行目标点迹关联。航向检测图如图 2所示。
图 2中,A、B、C为真实目标点迹,Z1、Z2为目标观测点迹;CF为航迹BC的延长线,以CF为中线,两边分别有角度为θ的航向限制角,即CD、CE与CF的夹角。设目标观测点与CF的夹角为α,当α在θ范围内时,则满足航向检测,否则不满足。图 2中可以看到,观测点Z1与CF的夹角在θ范围内,故满足航向检测条件,可以进行目标点关联;而观测点Z2与CF的夹角大于θ,故不满足航向检测条件,不予关联,将其舍弃。
对于目标航迹关联的整体流程图如图 3所示。
2 BRB推理验证航迹检测法本文采用BRB(belief rule-based)推理方法来验证前文所提航迹检测方法的可靠性。它是传统IF-THEN专家知识库的扩充,使用证据推理(ER)算法演绎推理过程,可以被用来处理多种形式的定性与定量不确定信息[9]。数学上,一个置信规则可以被定义为[10]:IFA1∧A2k∧…∧ATkTHEN {(C1,βk1),(C2,βk2),…,(CN,βkN)} with a rule weight(规则权重)θk and attribute weights(属性权重)δ1,δ2,…,δT,
上面提到的置信规则可以这样来理解,IF(A1isA1k,A2isA2k,…,andATisATk)THEN the consequence is {(C1,βk1),(C2,βk2),…,(CN,βkN)},其中A1,A2,…,AT是置信规则的前提属性,C1,C2,…,CN是用于评价结果的评价等级,Ajk∈{Cj1,Cj2,…,CjTj}是相对于前提属性Aj,j∈{1,2,…,N}的一个评价等级,βkl是第k条规则下结果在评价等级Cl,l∈{1,2,…,N}中的置信度,Tl是第k条规则下前提属性的总数量。如果
样本 | 速度/(m·s-1) | 航向/(°) |
1 | 2.09 | 3.89 |
2 | 4.02 | 5.36 |
3 | 4.32 | 5.46 |
4 | 3.04 | 2.26 |
5 | 7.08 | 3.08 |
6 | 7.45 | 4.11 |
7 | 2.68 | 2.14 |
8 | 3.25 | 5.32 |
9 | 1.72 | 2.48 |
10 | 2.40 | 5.04 |
在BRB系统中,其推理过程包括:多种形式信息的输入转换、规则激活权重的计算以及使用ER聚集BRB规则得到最终置信度。
2.1 BRB输入转换一个前提属性值的输入转换指的是将此值转换成不同的置信度,并将这些置信度分配给此前提属性的不同参考值。这就相当于将一个输入值转换成一个对应于前提属性参考值的置信度分布。特别地,一个前提属性Pi的输入值(它的置信度假设为xi)可以通过以下置信分布来转换[10]:
$S\left( {{P}_{i}},{{x}_{i}} \right)=\left\{ ({{h}_{in}},{{x}_{in}}),n=1,2,\cdots ,{{n}_{i}} \right\},i=1,2,\cdots ,{{T}_{l}}$ | (1) |
式中:S表示估计的置信度分配到前提属性的输入值的分布情况,hin(第i个值)是输入前提属性Pi的第n个参考值,αin是对应于参考值hin的置信度(αin≥0)。其中,
本文中,选取目标点的航向和速度作为BRB框架中的前提属性,它们的输入值可以通过历史的APRA雷达数据来获取,并将输入值转换为相应的前提属性参考值的隶属度[16]。对应于速度和航向的参考值,有:速度——Fast(速度快)、Low(速度慢)、Normal(正常速度);航向——Large(角度大)、Small(角度小)、Normal(正常角度)。此时,可以通过专家评价分配置信度αin给这些语言值(评价等级),该分配置信度随后被分布在前提属性不同参考值hin[High(H),Medium(M),Low(L))]的置信度αin的项中。上述输入转换的过程可以通过下面的方程式来阐述:
$\begin{align} & H\ge {{x}_{i}}\ge M,Medium=\frac{H-{{x}_{i}}}{H-M}, \\ & High=\left( 1-Medium \right),Low=0.00 \\ \end{align}$ | (2) |
$\begin{align} & M>{{x}_{i}}\ge L,Low=\frac{M-{{x}_{i}}}{H-L}, \\ & Medium=(1-Low),High=0.00 \\ \end{align}$ | (3) |
式中:H代表最大值,xi代表输入值,M代表中间值,L代表最小值。
在进行输入转换之前,需要设定前提属性的参考值,假设两个前提属性的参考值分布如下:速度={(Fast,10.4),(Normal,5.3),(Low,0.4)};航向={(Large,10.8),(Normal,4.7),(Small,1.2)}。然后结合式(8)、(9)完成对置信度αin的计算。比如,选择一个样本,速度=6m/s,航向=4°,则速度隶属于{(Fast,10.4),(Normal,5.3),(Low,0.4)}的程度为(0.1473,0.8527,0),即Medium=(10.4-6)/(10.4-5.3)=0.862,High=1-0.8627,Low=0;同理,航向隶属于(Large,10.8),(Normal,4.7),(Small,1.2)的程度为(0,0.800,0.200)。
2.2 计算激活权重一般地,在置信规则中,连接符号“∧”通常被用来代表前置属性的逻辑关系。这就意味着只有当规则中所有前因被激活时,得到的结果才具有可信度。根据上述置信分布的计算,第k条规则下激活权重ωk可以由式(4)来计算[10]。
${{\omega }_{k}}=\frac{{{\theta }_{k}}\prod\limits_{i-1}^{{{T}_{k}}}{{{\left( {{\alpha }_{ik}} \right)}^{{{{\bar{\delta }}}_{i}}}}}}{\sum\nolimits_{l=1}^{L}{[{{\theta }_{l}}\prod\limits_{i-1}^{{{T}_{i}}}{{{\left( {{\alpha }_{il}} \right)}^{{{{\bar{\delta }}}_{i}}}}]}}},{{\bar{\delta }}_{i}}=\frac{{{\delta }_{i}}}{{{\max }_{i=1,2,\cdots ,{{T}_{k}}}}\left\{ {{\delta }_{i}} \right\}}$ | (4) |
式中:αik(i=1,2,…,Tk)作为第i个前置属性的输入置信度,代表个体匹配度,被用来评定第k条规则下的参考值Aik。
ER方法是用来集合所有L规则下的前提属性数据分组,从而通过给定的前提属性Pi输入值获得结果属性中每一个参考值的信任程度。本文采用分析的ER算法[11],其产生的输出结果O(Y)由结果属性的参考值组成,如式(5)所示。
$O\left( Y \right)=S\left( {{P}_{i}} \right)=\left\{ \left( {{C}_{j}},{{\beta }_{j}} \right),j=1,2,\cdots N \right\}$ | (5) |
式中:βj代表的是对应于一个结果参考值Cj的置信度。βj的计算由ER算法的分析方程式(6)来求得。
最终的合并结果或通过ER推理所得到的输出表达式为:{(C1,β1),(C2,β2),…,(CN,βN)},其中βj是隶属于结果属性中第j个参考值Cj的最终置信度。本文中选取的输出结果分布分别是真实目标、虚假目标,根据表 1的样本数据,结合式(1)~(6),其训练后最终的置信规则库如表 2所示。
规则数 | 规则权值 | 速度 | 航向 | 结果 |
1 | 0.76 | Fast | Large | [0.6046,0.3954] |
2 | 0.90 | Fast | Normal | [0.8125,0.1875] |
3 | 0.90 | Fast | Small | [0.8170,0.1830] |
4 | 0.80 | Normal | Large | [0.7805,0.2195] |
5 | 1.00 | Normal | Normal | [0.9349,0.0651] |
6 | 1.00 | Normal | Small | [0.9844,0.0156] |
7 | 0.90 | Low | Large | [0.7603,0.2397] |
8 | 0.90 | Low | Normal | [0.8751,0.1249] |
9 | 0.59 | Low | Small | [0.5191,0.4809] |
${{\beta }_{j}}=\frac{\prod\limits_{k=1}^{L}{({{w}_{k}}{{\beta }_{jk}}+1-{{\omega }_{k}}\sum\limits_{j=1}^{N}{{{\beta }_{jk}})-\prod\limits_{k=1}^{L}{(1-{{\omega }_{k}}\sum\limits_{{}}^{N}{j=1{{\beta }_{jk}})}}}}}{\sum\limits_{j=1}^{N}{\prod\limits_{k=1}^{L}{({{\omega }_{k}}{{\beta }_{jk}}+1-{{\omega }_{k}}\sum\limits_{j=1}^{N}{{{\beta }_{jk}})-(N-1)\prod\limits_{k=1}^{L}{(1-{{w}_{k}}\sum\limits_{j=1}^{N}{{{\beta }_{jk}})-\prod\limits_{k=1}^{L}{(1-{{\omega }_{k}})}}}}}}}$ | (6) |
表 2中的置信度规则反映了速度、航向与真实目标之间的因果关系,以规则5为例,其含义是:在当前时刻,速度的取值在正常范围内,而航向也在其规定的正常范围内时,其确定的真实目标的置信度为93.49%,虚假目标的置信度为6.51%,且该规则成立的置信度为1,也就是说完全相信此规则下所得到的结果是成立的。其他的规则可以用同样的方法来解释。因此,文中所提出的以速度和航向来检测目标航迹的方法具有较高可靠性。
3 目标航迹Kalman滤波Kalman滤波器[14]输入的是笛卡尔坐标系下的目标位置,在目标跟踪过程中,通过时间tn-1下的系统状态,结合时间tn下的测量参数,利用Kalman滤波进行当前时刻tn下系统状态估计。该滤波器以递归的方式来计算它所得到的解[5]。
当前时间标记下的系统状态包含的是笛卡尔坐标系下目标点的位置和速度[1]:
${{\beta }_{j}}=\frac{\prod\limits_{k=1}^{L}{({{w}_{k}}{{\beta }_{jk}}+1-{{\omega }_{k}}\sum\limits_{j=1}^{N}{{{\beta }_{jk}})-\prod\limits_{k=1}^{L}{(1-{{\omega }_{k}}\sum\limits_{{}}^{N}{j=1{{\beta }_{jk}})}}}}}{\sum\limits_{j=1}^{N}{\prod\limits_{k=1}^{L}{({{\omega }_{k}}{{\beta }_{jk}}+1-{{\omega }_{k}}\sum\limits_{j=1}^{N}{{{\beta }_{jk}})-(N-1)\prod\limits_{k=1}^{L}{(1-{{w}_{k}}\sum\limits_{j=1}^{N}{{{\beta }_{jk}})-\prod\limits_{k=1}^{L}{(1-{{\omega }_{k}})}}}}}}}$ | (7) |
假设目标是匀速运动,则基于先前状态下的时刻t的状态预测方程可以由(8)式来表述:
${{X}_{t}}=A\cdot {{X}_{{{t}_{n-1}}}}+{{W}_{{{t}_{n-1}}}}$ | (8) |
式(8)等价于:
${{\left[ \begin{matrix} x \\ y \\ {{v}_{x}} \\ {{v}_{y}} \\ \end{matrix} \right]}_{{{t}_{n}}}}=\left[ \begin{matrix} 1 & 0 & T & 0 \\ 0 & 1 & 0 & T \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{matrix} \right]\cdot {{\left[ \begin{matrix} x \\ y \\ {{v}_{x}} \\ {{v}_{y}} \\ \end{matrix} \right]}_{{{t}_{n-1}}}}+{{\left[ \begin{matrix} 0 \\ 0 \\ {{w}_{x}} \\ {{w}_{y}} \\ \end{matrix} \right]}_{{{t}_{n-1}}}}$ | (9) |
式中:Wtn-1是期望值为零和协方差矩阵为Ex的高斯噪声状态向量,A是状态转移矩阵,T代表的是测量周期。
本文中,对实验的结果只考虑目标位置的估计,因此测量矩阵只包含了目标位置的信息:
$Z=\left[ \begin{matrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ \end{matrix} \right]$ | (10) |
因此观测方程可以表示为
${{Y}_{{{t}_{n}}}}=Z\cdot {{X}_{{{t}_{n}}}}+{{n}_{{{t}_{n}}}}$ | (11) |
式中n为高斯零均值和协方差矩阵Ey的测量噪声。
Kalman方程描述当前的状态和状态估计更新的预测为
$\left\{ \begin{align} & {{X}_{{{t}_{n}}}}=A\cdot {{X}_{{{t}_{n}}}} \\ & \bar{P}=A\cdot {{P}_{{{t}_{n-1}}}}\cdot {{A}^{T}}+{{E}_{x}} \\ & {{G}_{{{t}_{n}}}}={{{\bar{P}}}_{{{t}_{n}}}}\cdot {{Z}^{T}}\cdot (Z\cdot {{{\bar{P}}}_{{{t}_{n}}}}\cdot {{Z}^{T}}+{{E}_{y}})-1 \\ & {{X}_{{{t}_{n}}}}={{{\bar{X}}}_{{{t}_{n}}}}+{{G}_{{{t}_{n}}}}\cdot ({{Y}_{{{t}_{n}}}}-Z\cdot {{{\bar{X}}}_{{{t}_{n}}}}) \\ & {{P}_{{{t}_{n}}}}=(I-{{G}_{{{t}_{n}}}}\cdot Z)\cdot {{{\bar{P}}}_{{{t}_{n}}}} \\ \end{align} \right.$ | (12) |
式中:G代表的是Kalman增益,P是预测的误差协方差矩阵,P是状态估计误差的协方差矩阵。
协方差为测量噪声Ey和状态的误差估计Ex具有以下含义:当测量值受噪声的影响较大时,可以给测量噪声Ey赋予较大的权值。这样处理即默认更多地依赖于基于模型的滤波器的预测值,而不是噪声Ey的预测值,同时也假设所依靠的模型有较高的精确度,因此协方差矩阵具有更小的值。
4 仿真结果及分析使用真实的测量数据在MATLAB环境下进行目标点航迹关联以及Kalman滤波跟踪的仿真。航迹起始位置为(150m,100m),初始速度为(40m/s,30m/s),采样间隔T=2s,Monte Carlo仿真次数为150次。仿真流程图如图 4所示。
图 4中p代表观测点的位置,vp表示径向速度,θ表示航向。在每一次时间戳t,对上述这些数据进行检测关联,速度大于20m/s或小于0的点舍去;不在航向θ以内的点舍去,在这里选取的航向范围为:0°<θ<60°;最后结合最小邻域法选取最合适的目标点,并将其送入Kalman滤波器中进行目标跟踪,最终显示目标的航迹。通过获取的ARPA雷达数据得到船舶的位置信息,其位置绘制如图 5所示。
图 5中有900个ARPA目标点,可以明显看到目标的轨迹位于笛卡尔坐标的平面内,但是在这些目标点中,并不全是真实的目标点,还有些点属于杂波干扰所导致的虚假目标点,这些目标点需要被剔除。接下来的步骤是数据关联,在一定的范围内,不符合速度和航向判断条件的目标被去除掉,只有有效的目标点被显示出来。其数据关联后的目标点迹如图 6所示。
图 6中可以看到在进行数据的关联后,之前的干扰点以及不符合判断条件的目标点都被很好地剔除掉,为后续的目标跟踪及航迹显示提供了更加真实可靠的目标点。在进行Kalman滤波跟踪后所得到的目标轨迹平滑度取决于测量以及状态噪声协方差矩阵的选择。通过实验可知,要想得到平滑度较好的轨迹,需要对测量噪声矩阵分配较大的权值。Kalman滤波后的目标轨迹如图 7所示。
从以上实验结果看出,利用一个简单的线性Kalman滤波器,结合基于目标点径向速度、航向的检测(结合最近邻关联算法)的算法简单可行,且能得到令人满意的跟踪效果。图 8给出三种算法对目标跟踪的RMSE均方根误差,其中(a)为概率数据互联算法(PDAF),它主要用于解决杂波环境下的目标跟踪问题,并且误跟和丢失目标的概率较小,计算量相对较小。可以看到本文所提出的算法具有更高的跟踪精度。此外,表 3给出了150次Monte Carlo仿真实验中,三种算法所用的耗时,进一步论证相对于其他两种算法,本文算法的计算复杂度更小,跟踪实时性高。
1) 建立了基于BRB的雷达目标识别系统,该系统与传统的IF-THEN规则库相比较,其在知识的描述上更加丰富,应用上更具有现实意义。结果表明通过速度和航向来检测目标航迹具有较高的置信度。
2) 为了实现真实目标的轨迹跟踪,构建了两个重要的信号处理结构单元:数据关联单元和Kalman滤波跟踪单元。数据关联阶段主要完成更新被测目标现有的轨迹,通过对目标点径向速度和航向的判断,来实现虚假目标点的删除,其目的是尽可能地减少错误目标轨迹的生成。Kalman滤波器在本文中被用来消除测量噪声,使得真实目标轨迹更加平滑。
3) 将本文算法与其他传统算法进行了对比分析,实验结果表明,所提算法不仅在跟踪精度和鲁棒性上更好,而且用时也相对更少。 总之,通过仿真实验分析,本文所提算法无论是在雷达目标跟踪精度上,还是在可操作性上,均符合实际工程的要求,可以为雷达目标的跟踪与实现提供一种新的途径。
[1] | BAR-SHALOM Y, LI X R, KIRUBARAJAN T. Estimation with applications to tracking and navigation:theory, algorithms and software[M]. New York: John Wiley & Sons,2001 . |
[2] | MACAVEIU A, CÂMPEANU A. Automotive radar target tracking by Kalman filtering[C]//Proceedings of the 201311th International Conference on Telecommunication in Modern Satellite, Cable and Broadcasting Services. Nis, Serbia:IEEE, 2013, 2:553-556. |
[3] |
陈志敏, 薄煜明, 吴盘龙, 等. 拟蒙特卡罗粒子滤波改进算法及其在雷达目标跟踪中的应用[J].
应用科学学报, 2012, 30(6): 607–612.
CHEN Zhimin, BO Yuming, WU Panlong, et al. Improved quasi-monte-carlo particle filtering and its application to radar target tracking[J]. Journal of applied sciences, 2012, 30(6): 607–612. |
[4] | KUHN H W. The Hungarian method for the assignment problem[J]. Naval research logistics, 2005, 52(1): 7–21. |
[5] | SIMON D. Kalman filtering with state constraints:a survey of linear and nonlinear algorithms[J]. IET control theory and application, 2010, 4(8): 1303–1318. |
[6] |
李秀良. 基于聚类分析的雷达航迹跟踪数据关联算法[J].
指挥信息系统与技术, 2010, 1(3): 62–65.
LI Xiuliang. Data association algorithm in target tracking process of radar based on clustering analysis[J]. Command information system & technology, 2010, 1(3): 62–65. |
[7] | ABBAS M A, Shoukry A A. CMUNE:a clustering using mutual nearest neighbors algorithm[C]//Proceedings of the 201211th international conference on information science, signal processing and their applications (ISSPA). Montreal, QC:IEEE, 2012:1192-1197. |
[8] |
叶晓雪, 王新民, 李俨, 等. 修正概率关联算法在航迹关联中的应用[J].
计算机仿真, 2011, 28(3): 19–21.
YE Xiaoxue, WANG Xinmin, LI Yan, et al. Application of modified PDA algorithm in dense clutter[J]. Computer simulation, 2011, 28(3): 19–21. |
[9] |
李彬, 王红卫, 杨剑波, 等. 基于置信规则推理的库存控制方法[J].
华中科技大学学报:自然科学版, 2011(7): 76–79.
LI Bin, WANG Hongwei, YANG Jianbo, et al. Belief rule-based inference method for inventory control[J]. Journal of Huazhong University of Science and Technology:natural science edition, 2011(7): 76–79. |
[10] | YANG Jianbo, LIU Jun, WANG Jin, et al. Belief rule-base inference methodology using the evidential reasoning approach-RIMER[J]. IEEE transactions on systems, man, and cybernetics-part a:systems and humans, 2006, 36(2): 266–285. |
[11] | YANG Jianbo, XU Dongling. Evidential reasoning rule for evidence combination[J]. Artificial intelligence, 2013, 205: 1–29. |