CN104777465A - 基于b样条函数任意扩展目标形状及状态估计方法 - Google Patents

Abstract

本发明公开了一种基B样条函数的扩展目标跟踪与形状估计方法，它属于制导和智能信息处理技术领域，主要解决对任意形状扩展目标跟踪与形状估计问题。该方法通过引入多时刻联合统计的方法，构建了扩展目标的伪量测集，并根据伪量测集更新目标形状信息，采用B样条函数估计目标形状，从而实现了对任意几何形状扩展目标的跟踪与形状估计，且具有良好的形状估计精确度和鲁棒性，可以满足实际工程***的设计需求，具有良好的工程应用价值。

Description

基于B样条函数任意扩展目标形状及状态估计方法

技术领域

本发明属于智能信息处理技术领域，涉及扩展目标的形状和状态的估计方法。具体地说是一种基于B样条函数任意扩展目标形状及状态估计方法，可用于各种交通管制、机器人导航和精确制导等***中的目标跟踪与形状估计。

背景技术

传统的低分辨雷达等探测***中，目标被当作单个点来处理，因为其相对于传感器分辨单元来说太小，仅占据一个分辨单元。但随着现代雷达等探测设备分辨率的不断提高，目标的回波信号可能分布在不同的距离分辨单元中，其探测场不再等效为一个点，即单个目标可能产生多个量测，本发明中称这样的目标为扩展目标。针对扩展目标跟踪，单个点状态已经难以充分描述扩展目标，需要综合考虑目标形状等信息进行检测与跟踪分析。

目前，扩展目标形状估计方法主要包括：Random Matrices(RM)方法和RandomHypersurface Models(RHMs)方法。其中，RM方法使用扩展目标量测集的离散方差矩阵作为目标的形状参数，在贝叶斯框架下，将Wishart分布产生随机的形状参数矩阵作为先验形状，并根据下一时刻的量测信息更新后验形状。该方法的优点是对单帧数据处理比较简单，适用于集群目标的形状估计。但是，RM方法只能获得目标的近似椭圆形状，如果目标不是椭圆型（如星型或其他不规则形状），则该方法获得的形状信息将不准确，直接影响对目标状态的分析。RHMs方法是随机集模型的一种，通过设定的形状方程，随机产生方程中的参量作为先验参数，并通过量测集和***噪声构造约束条件来筛选先验参数，从而拟合出目标的形状，可以看出，该方法过分依赖于目标的先验形状，并根据先验形状参数确定目标的形状方程，如果参数不合理，将直接影响形状估计的精度。

发明内容

针对上述问题，本发明提出一种基于B样条函数的扩展目标形状及状态估计方法，采用多帧统计技术，在极坐标系下选取样本控制点，并用B样条函数拟合出目标的形状函数，实现高噪声、低量测数、低传感器精度情况下对任意扩展目标的形状近似估计。

实现本发明的关键技术是：在极坐标系下引入B样条技术，采用贝叶斯滤波框架进行多时刻联合估计，统计不同角度上目标的轮廓长度，运用B样条函数拟合目标形状，实现对任意扩展目标的形状估计。

为实现上述目标，具体实现步骤如下：

（1）初始化参数：目标状态ξ₀＝{x₀,X₀,P₀,Δ₀}，其中，x₀表示目标位置信息，X₀为形状信息，P₀为运动噪声协方差，Δ₀为形状噪声；并假定Q和R分别表示状态噪声协方差和量测噪声协方差；设定参量d、和m，其中，d为角度划分宽度，为固定不变的角度集合，m为伪量测集的最大元素数。

（2）当k≥1，根据量测Y_k和状态x_k进行卡尔曼滤波。

（3）以滤波后目标的位置信息为原点建立极坐标系，记录量测相对于原点的坐标，添加到伪量测集形成新的伪量测集Z_k+1。

（4）根据伪量测集Z_k+1更新目标形状

(3a)在0到2π之间均匀产生n个角度，构成固定不变的角度集合并统计过原点且指向θ_i角度方向上射线附近的点；先构造以指向θ_i角度方向射线为对称轴，原点为底边中点的矩形，底边长度设为2d，将矩形之内的点划分到集合D_k+1,i中，并统计D_k+1,i中所有点到矩形底边距离的期望，则可近似地获得伪量测集轮廓的长度，其中，d为先验常数；

(3b)统计D_k+1,i获得k+1时刻目标的轮廓集合，根据轮廓集合进行卡尔曼滤波，更新目标的形状信息，更新后的形状信息为X_k+1。

（5）以X_k+1作为控制顶点，生成三次B样条函数，估计目标的形状。

（6）若下一时刻观测信息到达，转到步骤（2）进行迭代；否则，目标跟踪过程结束。

本发明具有以下优点：

（1）本发明采用了多帧统计的思想建立了伪量测集，不需要对量测率的分布做任何假设，从而可以在低量测率、高噪声的情况下准确估计目标的形状。

（2）本发明在高噪声、低量测率的极端情况下，引入B样条函数，可以对任意形状的扩展目标进行形状估计，为后续的目标的身份识别、航迹关联等提供了可靠的信息特征。

附图说明

图1是本发明的整体流程图；

图2是本发明方法对伪量测集为十字型目标的形状估计效果图；

图3是本发明方法估计的形状与目标真实形状的对比图；

图4是本发明方法对伪量测集为“Y”型目标的形状估计效果图；

图5是本发明方法估计的形状与目标真实形状的对比图；

图6是采用本发明方法与RM方法对十字型扩展目标前20帧的形状估计图；

图7是前20帧的量测图；

图8是采用本发明方法与RM方法的平均形状估计效果图；

图9是采用本发明方法与RM方法对“Y”型扩展目标的前20帧的形状估计图；

图10是前20帧的量测图；

图11是本发明方法与RM方法的平均形状估计效果图；

具体实施方式

一、基础理论介绍

1.卡尔曼滤波技术

假设单个目标的状态方程和量测方程分别表示为：

x_k+1＝Fx_k+Gw_k    (1)

y_k＝h(x_k)+v_k    (2)其中，x_k表示目标在k时刻的状态，F为一步转移矩阵，h(·)表示观测模型，w_k和v_k分别表示状态噪声和量测噪声，对应的协方差分别表示为Q_k和R_k。

假定已知k时刻目标的状态x_k和协方差P_k，则卡尔曼滤波步骤如下：

（1）预测下一时刻的目标状态

x_k+1|k＝Fx_k+Gw_k

（2）预测下一时刻的协方差矩阵

P_k+1|k＝FP_kF^T+GQG^T

（3）计算增益

K_k+1＝P_k+1|k[S_k+1|k]^-1

S_k+1|k＝P_k+1|k+R_k+1

（4）根据最新量测进行状态更新

x_k+1＝x_k+1|k+K_k+1(y_k+1-x_k+1|k)

（5）更新协方差矩阵

P_k+1＝[I-K_k+1]P_k+1|k

2.基于B样条函数形状估计方法

B样条函数方法：通过有限个控制点，形成一条光滑的曲线函数，并采用其拟合控制点的轮廓形状。若给定控制顶点集合U_k＝[μ₁,μ₂,…,μ_n]^T，将控制顶点代入B样条函数，则可获得形参为u的B样条函数

$S\left(u\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\μ}}_i\·N_{i,l}\left(u\right)$

其中，N_i,l(u)表示B样条曲线函数，u表示形式变量，l表示B样条函数的次数，当l＝3时，则为三次B样条函数，且可表示为：

二、本发明基于B样条函数的任意扩展目标形状及状态估计方法

参照图1，本发明的具体实施步骤如下：

步骤1.令初始时刻k＝0，初始化参数x₀、P₀、X₀、Δ₀、d、、R、Q和m。

步骤2.当k≥1，对目标运动状态进行卡尔曼滤波

（2.1）预测下一时刻目标运动状态与协方差矩阵：

x_k+1|k＝Fx_k

P_k+1|k＝FP_kF^T+GQG^T

（2.2）根据新获得的量测集合Y_k+1更新目标状态和协方差矩阵：

$x_{k+1}=x_{k+1|k}+K_{k+1}(\frac{1}{\left|Y_{k+1}\right|}\underset{i=1}{\overset{\left|Y_{k+1}\right|}{\mathrm{\Σ}}}{y}_{k+1,i}-{\tilde{x}}_{k+1|k})$

P_k+1＝[I-K_k+1H]P_k+1|k

其中，|·|表示集合中元素的个数，H为量测矩阵；

${\tilde{x}}_{k+1|k}={\mathrm{Hx}}_{k+1|k}$

K_k+1＝P_k+1|kH[S_k+1|k]^-1

S_k+1|k＝HP_k+1|kH^T+R

步骤3.量测集处理，生成新的伪量测集Z_k+1

(3.1)先将k+1时刻的量测集Y_k+1中每个元素减去更新后的质心坐标然后并入Z_k构成新的伪量测集

${\hat{Z}}_{k+1}=Z_k\∪{\{y_{k+1,l}-{\tilde{x}}_{k+1}|y_{k+1,l}\∈Y_{k+1}\}}_{l=1}^{\left|Y_{k+1}\right|};$

(3.2)假定伪量测集元素最大个数为m，判断是否大于m，若则令否则删除中时刻靠前的个元素，再设定 $Z_{k+1}={\hat{Z}}_{k+1};$

步骤4.根据伪量测集Z_k+1更新目标形状

(4.1)在0到2π均匀产生n个角度，构成固定不变的角度集合并将伪量测集按角度进行子集划分

$D_{k+1,i}=\left\{z_{k+1,l}\right|f_{L_i}\left(z_{k+1,l}\right)<\mathrm{dandC}\left(z_{k+1,l}\right)=\mathrm{ture}\}$

$f_{L_i}\left(z_{k+1,l}\right)=\frac{\left|\right|B_1{z}_{k+1,l}^{\left(1\right)}+B_2{z}_{k+1,l}^{\left(2\right)}\left|\right|}{\sqrt{B_1^2+B_2^2}}$

其中，d为划分宽度，和表示k+1时刻伪量测集中第l个元素，在笛卡尔坐标系x方向和y方向上的坐标；为点z_k+1,l到过原点沿θ_i角度方向直线L_i的距离，其中，B₁和B₂为直线L_i上满足的参数；C(·)表示z_k+1，l在角度正方向范围内的约束条件；表示过原点且垂直于L_i的直线，其中，A₁和A₂为直线上满足 $\frac{A_2}{A_1}=\tan \left({\mathrm{\&theta;}}_i\right)$ 的参数。

(4.2)计算集合D_k+1,i中元素到直线距离的期望，则目标的外轮廓可表示为：

${\mathrm{\&rho;}}_i=2\frac{1}{\left|D_{k+1,i}\right|}\underset{l=1}{\overset{\left|D_{k+1,i}\right|}{\mathrm{\&Sigma;}}}\frac{\left|\right|A_1{z}_{k+1,l}^{\left(1\right)}+A_2{z}_{k+1,l}^{\left(2\right)}\left|\right|}{\sqrt{A_1^2+A_2^2}}$

z_k+1,l∈D_k+1,i

$\frac{A_2}{A_1}=\tan \left({\mathrm{\&theta;}}_i\right);$

(4.3)根据k+1时刻得到的外轮廓集合，对目标的形状进行卡尔曼滤波。

预测：

X_k+1|k＝X_k

${\mathrm{\&Delta;}}_{k+1|k}={\mathrm{\&Delta;}}_k+\stackrel{\&OverBar;}{H}\left|\right|R\left|\right|$

更新：

$X_{k+1}=X_{k+1|k}+(X_{k|k}-{\hat{K}}_{k+1}\&CenterDot;2E\left(\right|\left|D_{k+1,i}\right|\left|\right))$

${\mathrm{\&Delta;}}_{k+1}=(I-{\hat{K}}_{k+1}){\mathrm{\&Delta;}}_{k+1|k}$

${\hat{K}}_{k+1}=\mathrm{diag}\left({\left\{\frac{{\mathrm{\&delta;}}_{k+1|k,i}}{{\mathrm{\&delta;}}_{k+1|k,i}+w_i^{\left(k\right)}}\right\}}_{i=1}^n\right)$

其中， ${\mathrm{\&Delta;}}_{k+1|k}={\left[{\mathrm{\&delta;}}_{k+1|k,i}\right]}_{n\&times;1}^T.$

步骤5.以X_k+1为控制顶点，生成三次B样条函数，估计目标形状：

(5.1)将极坐标下的点集映射到平面坐标系

(5.2)令为U_k+1集合末尾添加集合首前三个元素的扩展控制顶点集，可生成闭合的B样条函数：

$s\left(u\right)=\underset{i=1}{\overset{n+3}{\mathrm{\&Sigma;}}}{\underset{\&OverBar;}{\mathrm{\&mu;}}}_{k+1}\&CenterDot;N_{i,3}\left(u\right)$

其中，N_i,3(u)为三次B样条函数

步骤6.重复步骤2，继续对扩展目标进行形状和状态估计。

本发明可通过以下实验仿真进一步说明：

1.仿真条件及参数

假设多个目标在x-y平面上作匀速运动，目标运动状态表示为x＝[x,v_x,y,v_y]^T，其中，x和y分别为单个目标在笛卡尔坐标系中x方向和y方向上的位置，v_x和v_y分别为每个目标在x方向和y方向上的速度。目标的状态方程如式（1）所示，其中， $F=\left[\begin{array}{cccc}1& T& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& T\\ 0& 0& 0& 1\end{array}\right],G=\left[\begin{array}{cc}{T}^2/2& 0\\ T& 0\\ 0& T^2/2\\ 0& T\end{array}\right],$ T表示采样时间间隔。

量测方程为y_k＝Hx_k+v_k，其中， $H=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right].$ 仿真场景中过程噪声协方差为其中σ_w1＝σ_w2＝1，测噪声协方差为其中σ_v1＝σ_v2＝1，目标运动状态的初始协方差P₀=diag[5,1,5,1]，目标形状信息的初始方差向量角度集合为均匀分布在[0，2π]上的n个角度其中n＝20，每个角度的采样间隔d＝1.5。初始伪量测集Z₀为空集，伪量测集最大元素数量m＝45。本发明与RM方法做对比分析，RM方法中初始形状状态X₀＝diag(10,10)，其他的状态参数均与以上参数相同。

2.仿真内容及结果分析

仿真实验，将本发明方法与RM方法进行对比实验分析，主要从以下三个方面开展实验：

实验1：不同形状目标产生的伪量测集的形状估计

图2是本发明方法对伪量测集为十字型目标的形状估计效果图。可以看出，本发明方法可以较准确地从伪量测集中估计出目标的形状。

图3是本发明方法估计的形状与目标真实形状的对比图。

图4是本发明方法对伪量测集为“Y”型目标的形状估计效果图。

图5是本发明方法估计的形状与目标真实形状的对比图。

可以看出，本发明方法可以对不同形状的扩展目标进行形状估计。

实验2：十字型形状扩展目标形状的连续估计

图6是采用本发明方法与RM方法对十字型扩展目标前20帧跟踪的形状估计结果。可以看出，本发明方法估计的形状比RM方法的椭圆估计结果更加精确。

图7是前20帧的量测图。对比图6可以看出，在量测集噪声较大的情况下，本发明方法仍可以估计出正确的形状，具有较强的抗干扰性。

图8是本发明方法与RM方法的平均形状估计结果。可以看出，本发明方法可以较好地获得目标的形状特征。

实验3：“Y”型扩展目标形状连续估计

图9是采用本发明方法与RM方法对十字型扩展目标前20帧跟踪的形状估计结果。

图10是前20帧的量测图。对比图9可以看出，受噪声影响量测形状的改变，不影响本发明方法对扩展目标形状的正确估计，具有较强的抗干扰性能。

图11是本发明方法与RM方法的平均形状估计结果。可以看出，本发明方法能够更加准确地估计出目标的形状特征。

Claims (4)

1.基B样条函数的扩展目标跟踪与形状估计方法，包括： 

（1）初始化参数：目标状态ξ₀＝{x₀,X₀,P₀,Δ₀}，其中，x₀表示目标位置信息，X₀为形状信息，P₀为运动噪声协方差，Δ₀为形状噪声；并假定Q和R分别表示状态噪声协方差和量测噪声协方差；设定参量d、和m，其中，d为角度划分宽度，为固定不变的角度集合，m为伪量测集的最大元素数； 

（2）当k≥1，根据量测集Y_k和状态x_k进行卡尔曼滤波； 

（3）以滤波后目标的位置信息为原点建立极坐标系，记录量测相对于原点的坐标，添加到伪量测集形成新的伪量测集Z_k+1； 

（4）根据伪量测集Z_k+1更新目标形状； 

（5）以X_k+1作为控制顶点，生成三次B样条函数，估计目标的形状； 

（6）若下一时刻观测信息到达，转到步骤（2）进行迭代；否则，目标跟踪过程结束。 

2.根据权利要求书1所述的扩展目标跟踪与形状估计方法，其中，步骤(3)所述的量测集处理，生成新的伪量测集Z_k+1，按下述方法计算得到： 

(2.1)先将k+1时刻的量测集Y_k+1中每个元素减去更新后的质心坐标 

然后并入Z_k构成新的伪量测集，即 

(2.2)根据给定的伪量测集元素最大个数m，判断是否大于m，若 则令否则删除中时刻靠前的个元素，再设定

3.根据权利要求书1所述的扩展目标跟踪与形状估计方法，其中，步骤(4)所述根据伪量测集Z_k+1更新目标形状，按以下步骤实现： 

(3.1)在0到2π均匀产生n个角度，构成固定不变的角度集合并将伪量测集按角度进行子集划分， 

其中，d为划分宽度，和表示k+1时刻伪量测集中第l个元素，在笛卡尔坐标系x方向和y方向上的坐标；为点z_k+1,l到过原点沿θ_i角度方向直线L_i的距离，其中，B₁和B₂为直线L_i上满足的参数；C(·)表示z_k+1,l在角度正方向范围内的约束条件；表示过原点且垂直于L_i的直线，其中，A₁和A₂为直线上满足的参数； 

(3.2)计算集合D_k+1,i中元素到直线距离的期望，则目标的外轮廓可表示为： 

z_k+1,l∈D_k+1,i

(3.3)根据k+1时刻得到的外轮廓集合，对目标的形状进行卡尔曼滤波， 

预测： 

X_k+1|k＝X_k

更新： 

其中，

4.根据权利要求书1所述的扩展目标跟踪与形状估计方法，其中，步骤(5)所述的以X_k+1为控制顶点，生成三次B样条函数，估计目标形状，按下述方法计算得到： 

(4.1)将极坐标下的点集映射到平面坐标系 

(4.2)令为U_k+1集合末尾添加集合首前三个元素的扩展控制顶点集，可生成闭合的B样条函数： 

其中，N_i,3(u)为三次B样条函数 

。