随着内河水运信息化的飞速发展，为了保证航运业安全有效地运行，雷达被广泛运用于追踪船舶运行状态方面，从而避免船舶发生碰撞造成不必要的经济损失。特别是长江干线，现在已经成为了世界上运输量最大并且最繁忙的通航河流。因此为了防止航行事故的发生，对船舶进行交通管制显得尤为重要。普通船用雷达一般通过人工标绘来获取船舶的航速、航向等航行数据。随着技术的不断发展，现在的船用雷达普遍具有自动标绘功能(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为真实目标点迹，Z₁、Z₂为目标观测点迹；CF为航迹BC的延长线，以CF为中线，两边分别有角度为θ的航向限制角，即CD、CE与CF的夹角。设目标观测点与CF的夹角为α，当α在θ范围内时，则满足航向检测，否则不满足。图 2中可以看到，观测点Z₁与CF的夹角在θ范围内，故满足航向检测条件，可以进行目标点关联；而观测点Z₂与CF的夹角大于θ，故不满足航向检测条件，不予关联，将其舍弃。

对于目标航迹关联的整体流程图如图 3所示。

2 BRB推理验证航迹检测法

本文采用BRB(belief rule-based)推理方法来验证前文所提航迹检测方法的可靠性。它是传统IF-THEN专家知识库的扩充，使用证据推理(ER)算法演绎推理过程，可以被用来处理多种形式的定性与定量不确定信息^[9]。数学上，一个置信规则可以被定义为^[10]：IFA₁∧A₂^k∧…∧A_T^kTHEN {(C₁,β_k1),(C₂,β_k2),…,(C_N,β_kN)} with a rule weight(规则权重)θ_k and attribute weights(属性权重)δ₁,δ₂,…,δ_T，$\sum\limits_{l=1}^{N}{{{\beta }_{kl}}\le 1}$。

上面提到的置信规则可以这样来理解，IF(A₁isA₁^k,A₂isA₂^k,…,andA_TisA_T^k)THEN the consequence is {(C₁,β_k1),(C₂,β_k2),…,(C_N,β_kN)}，其中A₁,A₂,…,A_T是置信规则的前提属性，C₁,C₂,…,C_N是用于评价结果的评价等级，A_j^k∈{C_j1,C_j2,…,C_jTj}是相对于前提属性A_j，j∈{1,2,…,N}的一个评价等级，β_kl是第k条规则下结果在评价等级C_l,l∈{1,2,…,N}中的置信度，T_l是第k条规则下前提属性的总数量。如果$\sum\limits_{l=1}^{N}{{{\beta }_{kl}}\le 1}$，那么置信规则被认为具有完整性，否则就不具有完整性。据此，为了建立一个BRB系统，离散的评价等级，比如：Hight、Medium、Low，就需要被定义到每一个前提属性和结果上。BRB系统结构中的几种类型参数：置信度、规则权重和属性权重，可以通过对收集到的样本数据(航向、速度)进行训练来获取，样本数据如表 1所示。经过训练后的BRB系统此时可以用来预测目标的真实性。

在BRB系统中，其推理过程包括：多种形式信息的输入转换、规则激活权重的计算以及使用ER聚集BRB规则得到最终置信度。

2.1 BRB输入转换

一个前提属性值的输入转换指的是将此值转换成不同的置信度，并将这些置信度分配给此前提属性的不同参考值。这就相当于将一个输入值转换成一个对应于前提属性参考值的置信度分布。特别地，一个前提属性P_i的输入值(它的置信度假设为x_i)可以通过以下置信分布来转换^[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表示估计的置信度分配到前提属性的输入值的分布情况，h_in(第i个值)是输入前提属性P_i的第n个参考值，α_in是对应于参考值h_in的置信度(α_in≥0)。其中，$\sum\limits_{n=1}^{{{n}_{i}}}{{{\alpha }_{in}}\le 1\left( i=1,2,\cdots ,{{T}_{l}} \right)}$，n_i是参考值的数量。

本文中，选取目标点的航向和速度作为BRB框架中的前提属性，它们的输入值可以通过历史的APRA雷达数据来获取，并将输入值转换为相应的前提属性参考值的隶属度^[16]。对应于速度和航向的参考值，有：速度——Fast(速度快)、Low(速度慢)、Normal(正常速度)；航向——Large(角度大)、Small(角度小)、Normal(正常角度)。此时，可以通过专家评价分配置信度α_in给这些语言值(评价等级)，该分配置信度随后被分布在前提属性不同参考值h_in[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代表最大值，x_i代表输入值，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,…,T_k)作为第i个前置属性的输入置信度，代表个体匹配度，被用来评定第k条规则下的参考值A_i^k。${{\theta }_{k}}\prod\limits_{i=1}^{{{T}_{k}}}{{{\left( {{\alpha }_{ik}} \right)}^{{{{\bar{\delta }}}_{i}}}}}$被称作联合匹配度，它反映的是在第k条规则下，整个前提属性输入值和参考值之间的匹配度。δ_i(或δ_i)代表属性权重，需要注意的是，如果$^{{{{\bar{\delta }}}_{i}}}$，那么${{\left( \alpha _{ik}^{k} \right)}^{{{{\bar{\delta }}}_{i}}}}=1$，这就意味着一个零重要性的前提属性对激活权重不会有任何影响；如果${{\bar{\delta }}_{i}}=1$，那么${{(\alpha _{ik}^{k})}^{{{{\bar{\delta }}}_{i}}}}=\alpha _{ik}^{k}$，这就意味着大多数的前提属性对激活权重具有重大影响。随着激活权重的计算完成，接下来对于特定输入向量下的估计输出可以用ER理论推理得到。

2.3 ER推理输出

ER方法是用来集合所有L规则下的前提属性数据分组，从而通过给定的前提属性P_i输入值获得结果属性中每一个参考值的信任程度。本文采用分析的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代表的是对应于一个结果参考值C_j的置信度。β_j的计算由ER算法的分析方程式(6)来求得。

最终的合并结果或通过ER推理所得到的输出表达式为：{(C₁,β₁),(C₂,β₂),…,(C_N,β_N)}，其中β_j是隶属于结果属性中第j个参考值C_j的最终置信度。本文中选取的输出结果分布分别是真实目标、虚假目标，根据表 1的样本数据，结合式(1)～(6)，其训练后最终的置信规则库如表 2所示。

${{\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]输入的是笛卡尔坐标系下的目标位置，在目标跟踪过程中，通过时间t_n-1下的系统状态，结合时间t_n下的测量参数，利用Kalman滤波进行当前时刻t_n下系统状态估计。该滤波器以递归的方式来计算它所得到的解^[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)

式中：W_tn-1是期望值为零和协方差矩阵为E_x的高斯噪声状态向量，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为高斯零均值和协方差矩阵E_y的测量噪声。

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是状态估计误差的协方差矩阵。

协方差为测量噪声E_y和状态的误差估计E_x具有以下含义：当测量值受噪声的影响较大时，可以给测量噪声E_y赋予较大的权值。这样处理即默认更多地依赖于基于模型的滤波器的预测值，而不是噪声E_y的预测值，同时也假设所依靠的模型有较高的精确度，因此协方差矩阵具有更小的值。

4 仿真结果及分析

使用真实的测量数据在MATLAB环境下进行目标点航迹关联以及Kalman滤波跟踪的仿真。航迹起始位置为(150m,100m)，初始速度为(40m/s,30m/s)，采样间隔T=2s，Monte Carlo仿真次数为150次。仿真流程图如图 4所示。

图 4中p代表观测点的位置，v_p表示径向速度，θ表示航向。在每一次时间戳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仿真实验中，三种算法所用的耗时，进一步论证相对于其他两种算法，本文算法的计算复杂度更小，跟踪实时性高。

5 结论

1) 建立了基于BRB的雷达目标识别系统，该系统与传统的IF-THEN规则库相比较，其在知识的描述上更加丰富，应用上更具有现实意义。结果表明通过速度和航向来检测目标航迹具有较高的置信度。

2) 为了实现真实目标的轨迹跟踪，构建了两个重要的信号处理结构单元：数据关联单元和Kalman滤波跟踪单元。数据关联阶段主要完成更新被测目标现有的轨迹，通过对目标点径向速度和航向的判断，来实现虚假目标点的删除，其目的是尽可能地减少错误目标轨迹的生成。Kalman滤波器在本文中被用来消除测量噪声，使得真实目标轨迹更加平滑。

3) 将本文算法与其他传统算法进行了对比分析，实验结果表明，所提算法不仅在跟踪精度和鲁棒性上更好，而且用时也相对更少。 总之，通过仿真实验分析，本文所提算法无论是在雷达目标跟踪精度上，还是在可操作性上，均符合实际工程的要求，可以为雷达目标的跟踪与实现提供一种新的途径。