CN111693962B - 一种基于交叉检验的目标运动模型估计方法 - Google Patents

Abstract

本发明涉及一种基于交叉检验的目标运动模型估计方法，属于传感器目标跟踪与数据融合领域。本发明为了克服拟合多项式描述目标运动模型时存在的过拟合现象，以及多维度下模型估计数量的指数级增长的问题，采用基于交叉检验的方法，依据s‑1个测量点数据计算多维度下目标可能存在的所有运动模型，将第s个测量点与旧有测量点特征值的交叉检验结果，作为基于s个测量点数据确定目标运动模型的依据，从而筛选出较少且有效目标运动模型。采用本发明确定的多维度下的目标运动模型估计数量大大减少，有效避免了过拟合现象，为在本系统中开展目标跟踪与数据融合工作提供了一种简单、有效且易于工程实现的目标运动模型估计方法。

Description

一种基于交叉检验的目标运动模型估计方法

技术领域

本发明属于传感器目标跟踪与数据融合领域，具体涉及一种基于交叉检验的目标运动模型估计方法。

背景技术

基于地面雷达测量数据估计飞机等空中目标的运动模型，是各类传感器数据融合及滤波算法提高目标状态估计准确性、降低计算复杂度的关键环节，也有助于预测目标趋势以及判别目标企图。一维空间中的目标运动模型可以描述为测量值相对于时间的函数，可以是一次函数，表示时间与测量值之间是直线关系，也可以是二次及二次以上曲线函数，通常使用拟合多项式表示。拟合多项式次数越高，对已有测量数据的拟合误差就越小，但会存在过拟合现象，即随着拟合次数的增高，新的测量点与拟合点之间误差可能会增大。当空中目标运动模型扩展到多维度时，各个维度上存在的过拟合现象会严重影响目标运动模型的估计与选择。比如：如果目标在x轴方向上的运动模型估计有a种，在y轴方向上的运动模型估计有b种，那么目标在平面坐标系中的运动模型估计就有a×b种。采用拟合多项式描述目标运动模型时存在的过拟合现象和多维度下空中目标运动模型估计数量的指数级增长，会严重影响传感器数据融合精度与效率。

发明内容

(一)要解决的技术问题

本发明要解决的技术问题是如何提供一种基于交叉检验的目标运动模型估计方法，以克服采用拟合多项式描述目标运动模型时存在的过拟合现象和多维度下空中目标运动模型估计数量的指数级增长的问题。

(二)技术方案

为了解决上述技术问题，本发明提出一种基于交叉检验的目标运动模型估计方法，该方法包括如下步骤：

步骤1：将某雷达针对某个空中目标的s个连续测量点数据(ρ_i,θ_i,h_i,t_i)转换为统一直角坐标{(X_i,Y_i,t_i)|i＝1,2,...,s}，其中，(ρ_i,θ_i,h_i,t_i)表示t_i时刻该雷达测得的目标距离ρ_i、方位θ_i和高度h_i，i＝1,2,…s；

步骤2：对目标X轴运动分量{(t_i,X_i)|i＝1,2,...,s-1}，使用基于最小二乘的多项式拟合方法，分别得到最高次数为1,2,...,m时的X轴拟合多项式系数

步骤3：对目标Y轴运动分量{(t_i,Y_i)|i＝1,2,...,s-1}，使用基于最小二乘的多项式拟合方法，分别得到最高次数为1,2,...,m时的Y轴拟合多项式系数

步骤4：依据s-1个测量点数据得到的目标在平面直角坐标系下可能存在的所有运动模型集合M′_m×m，用多项式系数表示为：

其中：

等价为：

j∈[1,m],i′∈[1,m]，表示目标运动模型在X轴方向上按照j次多项式拟合，在Y轴方向上按照i′次多项式拟合；

步骤5：计算第s个测量点(t_s,X_s,Y_s)按照运动模型集合M′_m×m中的每一个元素分别进行多项式拟合后的拟合误差矩阵D′_m×m为

步骤6：求前s-1个测量点中，按照运动模型集合M′_m×m中的每一个元素分别进行多项式拟合后拟合误差的最大值矩阵D″_m×m为

步骤7：通过D′_m×m和D″_m×m的均值构造模型筛选矩阵D_m×m为：

步骤8：查找所述模型筛选矩阵D_m×m中值最小的元素d_u,v(u∈[1,m],v∈[1,m])，裁剪并变换得到邻域矩阵B_3×3，根据邻域矩阵中所有值为0的元素对应的运动模型参数生成有效运动模型数组M″_L×2(L∈[1,9])；

步骤9：按照有效运动模型数组M″_L×2中的参数，依据s个测量点数据，在X、Y两个维度上分别采用基于最小二乘的多项式拟合方法求拟合多项式系数，最终得到目标运动模型估计集合M_L×2。

(三)有益效果

本发明针对拟合多项式描述目标运动模型时存在的过拟合现象，以及多维度下模型估计数量的指数级增长问题，采用基于交叉检验的方法，依据s-1个测量点数据计算多维度下目标可能存在的所有运动模型，将第s个测量点与旧有测量点特征值的交叉检验结果，作为基于s个测量点数据确定目标运动模型的依据，从而达到筛选出较少且有效目标运动模型的目的。通过对后续2个新测量点的实际检验用例表明，采用本发明确定的多维度下的目标运动模型估计数量大大减少，有效避免了过拟合现象，为在本系统中开展目标跟踪与数据融合工作提供了一种简单、有效且易于工程实现的目标运动模型估计方法。本发明模型估计方法科学、方案实施步骤合理，模型估计结果理想，对于提高后续目标跟踪与数据融合的准确性、降低计算复杂度具有重要意义，同时也有助于目标特征分析、目标趋势预测、目标企图判别等工作的开展。本发明所提供的方法的时间复杂度和空间复杂度都很低，可操作性和实用性很强。

附图说明

图1为本发明目标运动模型估计方法的主体实现步骤；

图2为本发明实施例中测量点在统一直角坐标系中的显示图；

图3为本发明实施例中模型筛选与新点验证误差数据相关关系图。

具体实施方式

为使本发明的目的、内容和优点更加清楚，下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。

针对拟合多项式描述目标运动模型时存在的过拟合现象，以及多维度下模型估计数量的指数级增长问题，提出了一种基于交叉检验的目标运动模型估计方法，采用基于交叉检验的方法，依据s-1个测量点数据计算多维度下目标可能存在的所有运动模型，将第s个测量点与旧有测量点特征值的交叉检验结果，作为基于s个测量点数据确定目标运动模型的依据，从而达到筛选出较少且有效目标运动模型的目的，尽可能地减小多维度下目标运动模型的估计数量，筛选出尽量少而且属于有效拟合的目标运动模型。通过对后续2个新测量点的检验表明，采用本发明确定的多维度下的目标运动模型估计数量大大减少，有效避免了过拟合现象，为在本系统中开展目标跟踪与数据融合工作提供了一种简单、有效且易于工程实现的目标运动模型估计方法。

为实现上述目的，本发明提供了一种基于交叉检验目标运动模型估计方法，所述估计方法应用于传感器目标跟踪与数据融合系统中。

所述估计方法包括如下步骤：

步骤1：将某雷达针对某个空中目标的s个连续测量点数据(ρ_i,θ_i,h_i,t_i)(表示t_i时刻该雷达测得的目标距离ρ_i、方位θ_i和高度h_i，i＝1,2,…s)转换为统一直角坐标{(X_i,Y_i,t_i)|i＝1,2,...,s}。

步骤2：对目标X轴运动分量{(t_i,X_i)|i＝1,2,...,s-1},使用基于最小二乘的多项式拟合方法，分别得到最高次数为1,2,...,m时的X轴拟合多项式系数

步骤3：对目标Y轴运动分量{(t_i,Y_i)|i＝1,2,...,s-1},使用基于最小二乘的多项式拟合方法，分别得到最高次数为1,2,...,m时的Y轴拟合多项式系数

步骤4：依据s-1个测量点数据得到的目标在平面直角坐标系下可能存在的所有运动模型集合M′_m×m，用多项式系数表示为：

其中：

等价为：

j∈[1,m],i′∈[1,m]，表示目标运动模型在X轴方向上按照j次多项式拟合，在Y轴方向上按照i′次多项式拟合。

步骤5：计算第s个测量点(t_s,X_s,Y_s)按照运动模型集合M′_m×m中的每一个元素分别进行多项式拟合后的拟合误差矩阵D′_m×m。

步骤6：求前s-1个测量点中，按照运动模型集合M′_m×m中的每一个元素分别进行多项式拟合后拟合误差的最大值矩阵D″_m×m。

步骤7：通过D′_m×m和D″_m×m的均值构造模型筛选矩阵D_m×m为：

步骤8：查找模型筛选矩阵D_m×m中值最小的元素d_u,v(u∈[1,m],v∈[1,m])，裁剪并变换得到邻域矩阵B_3×3。根据邻域矩阵中所有值为0的元素对应的运动模型参数生成有效运动模型数组M″_L×2(L∈[1,9])。

步骤9：按照有效运动模型数组M″_L×2中的参数，依据s个测量点数据，在X、Y两个维度上分别采用基于最小二乘的多项式拟合方法求拟合多项式系数，最终得到目标运动模型估计集合M_L×2。

所述步骤1包括如下步骤：

步骤1.1：获取雷达观测数据(ρ_i,θ_i,h_i,t_i)，i＝1,2,…s，为保证目标运动模型估计的即时性，一般选取不超过距当前时刻90秒的雷达测量点数据。比如雷达测量周期为10秒时，s取值7或8为宜。如果是两坐标雷达观测数据，则认为h_i＝0。θ_i取值范围[0,360)，单位为度，其中雷达正北方向为0度，正东方向为90度，计算过程中需要转换为弧度制。

步骤1.2：将雷达观测数据(ρ_i,θ_i,h_i,t_i)，i＝1,2,…s转换为以本站为中心的二维直角坐标(X_zi,Y_zi)：

步骤1.3：将(X_zi,Y_zi)，i＝1,2,…s转换为中心统一直角坐标(X_i,Y_i)：

X_i＝X_zicosδ_xz-Y_zisinδ_xz+X_zx

Y_i＝X_zisinδ_xz+Y_zicosδ_xz+Y_zx

其中：(X_zx,Y_zx)为雷达站址在中心统一直角坐标系中的坐标；δ_xz为雷达站址与直角坐标系中心点的经度差，单位为弧度。

所述步骤2包括如下步骤：

步骤2.1：将X轴运动分量{(t_i,X_i)|i＝1,2,...,s-1}简记为：{(x_i,y_i)||i＝1,2,...,n}，然后确定m值。m表示在X轴方向上目标运动模型以拟合多项式表示时的最高次数，代表了目标在该维度上依赖于时间的运动模型复杂度，m值越大，表示目标在给定时间内可能具有的机动性越高，m≥1。对于常规飞机目标，m取值为4或5即可满足要求。

步骤2.2：针对n个测量点数据{(x_i,y_i)||i＝1,2,...,n}，令j依次取值1,2,...,m，采用基于最小二乘的j次多项式拟合方法，求得X轴方向上m组拟合多项式系数，记为：

其中，基于最小二乘的j次多项式拟合方法实现原理如下：

以j次多项式形式描述的目标运动模型为：

y＝f_j(x)＝k₀+k₁x+k₂x²+...+k_px^j

其中：{k_p|p＝0,1,...,j}是待定的拟合多项式f_j(x)的系数。基于n个测量点数据{(x_i,y_i)||i＝1,2,...,n}，采用最小二乘法求{k_p|p＝0,1,...,j}值。即令：

当F(k_p|p＝0,1,...,j)取得最小值时，对于待定参数{k_p|p＝0,1,...,j}，有：

令J＝j+1，对于上式的求解，可以等价为：求J阶线性代数方程组E_J×JK_J＝C_J的解，其中：

r＝0,1,2,...,j×2，

p＝0,1,2,...,j。

n表示测量点个数，(x_i,y_i)表示第i个测量点的坐标，

表示x_i的p次方。我们采用全选主元高斯消去法可以成功求得J阶线性代数方程组E_J×JK_J＝C_J的解，记为：

所述步骤3包括如下步骤：

步骤3.1：将Y轴运动分量{(t_i,Y_i)|i＝1,2,...,s-1}简记为：{(x_i,y_i)||i＝1,2,...,n}，然后确定m值。m值意义如步骤2.1所述。通常认为目标在Y轴方向上的机动性能与X轴上相同。

步骤3.2：针对n个测量点数据{(x_i,y_i)||i＝1,2,...,n}，令j依次取值1,2,...,m，采用基于最小二乘的j次多项式拟合方法，求得Y轴方向上m组拟合多项式系数，记为：

求解原理同步骤2.2。

所述步骤5包括如下步骤：

步骤5.1：令j＝1，i′＝1，其中“＝”表示赋值为。

步骤5.2：取运动模型集合M′_m×m的元素m′_ji′，得到X轴上的拟合多项式系数

求第s个测量点(t_s,X_s,Y_s)在X轴向的拟合误差

步骤5.3：取运动模型集合M′_m×m的元素m′_ji′，得到Y轴上的拟合多项式系数

求第s个测量点(t_s,X_s,Y_s)在Y轴向的拟合误差

步骤5.4：则第s个测量点按照运动模型m′_ji′，在平面直角坐标系中的拟合误差为

其中j∈[1,m],i′∈[1,m]。

步骤5.5：i′值加1,如果i′≤m则转步骤5.2，否则转步骤5.6。

步骤5.6：i′＝1，j值加1,如果j≤m则转步骤5.2，否则转步骤5.7。

步骤5.7：则第s个测量点按照运动模型集合M′_m×m进行多项式拟合后的拟合误差矩阵为：

所述步骤6包括如下步骤：

步骤6.1：令i＝1，j＝1，i′＝1，w″＝0，其中“＝”表示赋值为。

步骤6.2：取运动模型集合M′_m×m的元素m′_ji′，得到X轴上的拟合多项式系数

求第i个测量点(X_i,Y_i,t_i)在X轴向的拟合误差

步骤6.3：取运动模型集合M′_m×m的元素m′_ji′，得到Y轴上的拟合多项式系数

求第i个测量点(X_i,Y_i,t_i)在Y轴向的拟合误差

步骤6.4：则第i个测量点按照运动模型m′_ji′，在平面直角坐标系中的拟合误差为

其中i∈[1,s-1],j∈[1,m],i′∈[1,m]。

步骤6.5：记录当前拟合误差最大值：如果

则

步骤6.6：i值加1,如果i≤s-1则转步骤6.2，否则转步骤6.7。

步骤6.7：d″_ji′＝w″，i＝1，w″＝0，i′值加1,如果i′≤m则转步骤6.2，否则转步骤6.8。

步骤6.8：i′＝1，j值加1,如果j≤m则转步骤6.2，否则转步骤6.9。

步骤6.9：前s-1个测量点中，按照运动模型集合M′_m×m中的每一个元素分别进行多项式拟合后拟合误差的最大值矩阵D″_m×m为：

所述步骤8包括如下步骤：

步骤8.1：找到模型筛选矩阵D_m×m中值最小的元素d_u,v(u∈[1,m],v∈[1,m])。

步骤8.2：裁剪D_m×m得到邻域矩阵B_3×3，并进行规则化变换：

其中：b_g,h(g∈[1,3],h∈[1,3])，取值变换规则如下：

b₂₂＝0，

其中α为邻域矩阵筛选系数，α≥1，一般根据经验给出，表示大于矩阵D_m×m中最小值元素α倍的元素所在位置对应的组合拟合多项式为非有效运动模型。

步骤8.3：根据邻域矩阵生成有效运动模型数组M″_L×2，L表示数组长度，L∈[1,9]。邻域矩阵B_3×3中值为0的元素对应的拟合多项式均为有效的目标运动模型，因此数组下标L取值最小为1，最大为9。二维数组M″_L×2的第一列表示目标在X轴方向上的有效拟合多项式次数，第二列表示目标在Y轴方向上的有效拟合多项式次数，均是大于等于1，小于等于m的自然数。例如：若b₃₂＝0，则对应的目标在X轴方向上作u+1次、在Y轴方向上作v次多项式拟合而构成的运动模型是目标可能存在的运动方式，此时数组元素M″_l1＝u+1，M″_l2＝v，l∈[1,L]。

所述步骤9包括如下步骤：

步骤9.1：令l＝1。

步骤9.2：令j＝M″_l1。

步骤9.3：将s个测量点的X轴运动分量{(t_i,X_i)|i＝1,2,...,s}简记为：{(x_i,y_i)||i＝1,2,...,n}，采用基于最小二乘的j次多项式拟合方法，求拟合多项式系数。包括如下步骤：

步骤9.3.1：求

r＝0,1,2,...,j×2，

p＝0,1,2,...,j。其中：n表示测量点个数，(x_i,y_i)表示第i个测量点的坐标，

表示x_i的p次方。

步骤9.3.2：令J＝j+1，构造矩阵E_J×J，C_J和K_J，其中：

步骤9.3.3：采用全选主元高斯消去法求J阶线性代数方程组E_J×JK_J＝C_J的解，记为：

步骤9.4：令j＝M″_l2。

步骤9.5：将s个测量点的Y轴运动分量{(t_i,Y_i)|i＝1,2,...,s}简记为：{(x_i,y_i)||i＝1,2,...,n}，采用基于最小二乘的j次多项式拟合方法，求拟合多项式系数。包括如下步骤：

步骤9.5.1：同步骤9.3.1。

步骤9.5.2：同步骤9.3.2。

步骤9.5.3：采用全选主元高斯消去法求J阶线性代数方程组E_J×JK_J＝C_J的解，记为：

步骤9.6：l值加1,如果l≤L则转步骤9.2，否则结束程序。

最终我们得到目标运动模型估计集合M_L×2表示为：

其中：L∈[1,9]，ω_l1是目标第l∈[1,L]种运动模型估计中在X轴方向上的拟合多项式系数，ω_l2是对应的其在Y轴方向上的拟合多项式系数，通常情况下ω_l1≠ω_l2，即第l种模型估计中目标在X轴和Y轴方向上分别按照不同的多项式次数和系数值进行拟合，不存在过拟合现象。

具体实施例1

本实施例具体描述本发明所提出的基于交叉检验的目标运动模型估计方法，所述方法应用于传感器目标跟踪与数据融合系统中。

所述估计方法包括如下步骤：

步骤1：将某雷达针对某个空中目标的s个连续测量点数据(ρ_i,θ_i,h_i,t_i)(表示t_i时刻该雷达测得的目标距离ρ_i、方位θ_i和高度h_i，i＝1,2,…s)转换为统一直角坐标{(X_i,Y_i,t_i)|i＝1,2,...,s}。

步骤1.1：为保证目标运动模型估计的即时性，根据雷达测量周期，本例中s取值为7。因为是两坐标雷达测量数据，h_i均取值为0。θ_i取值范围[0,360)，单位为度，其中雷达正北方向为0度，正东方向为90度，计算过程需要转换为弧度制。测量过程中t_i、ρ_i和θ_i均存在系统误差和随机误差。雷达站址统一直角坐标为(1485，2490)千米，定位系统误差为(0.2，-0.25)千米，雷达站址与直角坐标系中心点的经度差为π/180弧度。为验证实施效果，在当前7个测量点基础上，增加2个新点。雷达测量数据见表1.1。

表1.1：雷达测量数据

序号	时间t(秒)	距离ρ(km)	方位θ(度)	备注
					1	16	47.04	217.87	交叉验证点
2	26	47.20	221.43	交叉验证点
					3	36	48.61	223.83	交叉验证点
4	46	50.67	226.13	交叉验证点
					5	57	53.19	225.80	交叉验证点
6	67	55.55	226.65	交叉验证点
					7	77	58.21	225.08	当前时刻点
8	87	60.58	225.02	效果检验点
					9	97	63.10	224.01	效果检验点

步骤1.2：将雷达观测数据(ρ_i,θ_i,h_i,t_i)，i＝1,2,…9转换为以本站为中心的二维直角坐标(X_zi,Y_zi)：

计算结果见表1.2。

表1.2：以雷达为中心的二维直角坐标

步骤1.3：将(X_zi,Y_zi)，i＝1,2,…9转换为中心统一直角坐标(X_i,Y_i)：

X_i＝X_zicosδ_xz-Y_zisinδ_xz+X_zx

Y_i＝X_zisinδ_xz+Y_zicosδ_xz+Y_zx

其中：(X_zx,Y_zx)为雷达站址在中心统一直角坐标系中的修正坐标。δ_xz为雷达站址与直角坐标系中心点的经度差(单位为弧度)。计算结果见表1.3。所有测量点在统一直角坐标系中的显示如图2所示。

表1.3：测量点的中心统一直角坐标

步骤2：对目标X轴运动分量{(t_i,X_i)|i＝1,2,...,s-1},使用基于最小二乘的多项式拟合方法，分别得到最高次数为1,2,...,m时的X轴拟合多项式系数

包括如下步骤：

步骤2.1：将X轴运动分量{(t_i,X_i)|i＝1,2,...,7-1}简记为：{(x_i,y_i)||i＝1,2,...,6}，见表1.4所示。对于常规飞机目标，在X轴方向上的运动模型以拟合多项式表示时的最高次数m取值为5即可满足要求。

表1.4：X轴运动分量简记为(x_i,y_i)

步骤2.2：针对以上6个测量点数据{(x_i,y_i)||i＝1,2,...,6}，令j依次取值1,2,...,m，采用基于最小二乘的j次多项式拟合方法，求得X轴方向上m＝5组拟合多项式系数，记为：

我们以j＝4为例说明基于最小二乘的j次多项式拟合方法求解步骤如下。

计算：

令：

采用全选主元高斯消去法求得5阶线性代数方程组E_5×5K₅＝C₅的解为：

记为：

步骤3：对目标Y轴运动分量{(t_i,Y_i)|i＝1,2,...,6},使用基于最小二乘的多项式拟合方法，分别得到最高次数为1,2,...,5时的Y轴拟合多项式系数

包括如下步骤：

步骤3.1：将Y轴运动分量{(t_i,Y_i)|i＝1,2,...,6}简记为：{(x_i,y_i)||i＝1,2,...,6}，见表1.5所示。对于常规飞机目标，在Y轴方向上的运动模型以拟合多项式表示时的最高次数m取值为5即可满足要求。通常认为目标在Y轴方向上的机动性能与X轴上相同。

表1.5：Y轴运动分量简记为(x_i,y_i)

序号	x	y
			1	16	2452.12
2	26	2453.82
			3	36	2454.10
4	46	2454.01
			5	57	2452.01
6	67	2450.92

步骤3.2：针对以上6个测量点数据{(x_i,y_i)||i＝1,2,...,6}，令j依次取值1,2,...,5，采用基于最小二乘的j次多项式拟合方法，求得Y轴方向上5组拟合多项式系数，记为：

求解方法步骤2.2。

步骤4：依据6个测量点数据{(t_i,X_i,Y_i)|i＝1,2,...,6}(见表1.3)得到的目标在平面直角坐标系下可能存在的所有运动模型集合M′_5×5，用多项式系数表示为：

其中：m′_ji′，j∈[1,5],i′∈[1,5]，表示目标运动模型在X轴方向上按照j次多项式拟合，在Y轴方向上按照i′次多项式拟合，计算结果如下：

步骤5：计算第7个测量点(t₇,X₇,Y₇)＝(77,1444.71,2447.93)按照运动模型集合M′_5×5中的每一个元素分别进行多项式拟合后的拟合误差矩阵D′_5×5。

步骤5.1：以j＝3，i′＝4为例，计算过程如下。

步骤5.2：取运动模型集合M′_5×5的元素m′₃₄，得到X轴上的拟合多项式系数：

求第7个测量点(t₇,X₇,Y₇)在X轴向的拟合误差

步骤5.3：取运动模型集合M′_5×5的元素m′₃₄，得到Y轴上的拟合多项式系数：

求第7个测量点(t₇,X₇,Y₇)在Y轴向的拟合误差

步骤5.4：则第7个测量点(t₇,X₇,Y₇)按照运动模型m′₃₄，在平面直角坐标系中的拟合误差d′₃₄为：

最终计算得到第7个测量点(t₇,X₇,Y₇)按照运动模型集合M′_5×5进行多项式拟合后的拟合误差矩阵D′_5×5为：

步骤6：求前6个测量点中，按照运动模型集合M′_5×5中的每一个元素分别进行多项式拟合后拟合误差的最大值矩阵D″_5×5。

步骤6.1：令i＝1，j＝1，i′＝1，w″＝0。

以i＝2，j＝3，i′＝2为例，计算过程如下。

步骤6.2：取运动模型集合M′_5×5的元素m′₃₂，得到X轴上的拟合多项式系数：

求第i＝2个测量点(t₂,X₂,Y₂)＝(26,1454.59,2453.82)在X轴向的拟合误差

为：

步骤6.3：取运动模型集合M′_5×5的元素m′₃₂，得到Y轴上的拟合多项式系数：

求第i＝2个测量点(t₂,X₂,Y₂)＝(26,1454.59,2453.82)在Y轴向的拟合误差

为：

步骤6.4：则第i＝2个测量点(t₂,X₂,Y₂)＝(26,1454.59,2453.82)按照运动模型m′₃₂，在平面直角坐标系中的拟合误差

为：

步骤6.5：记录当前拟合误差最大值：如果

则

步骤6.6：i值加1,如果i≤6则转步骤6.2，否则转步骤6.7。

步骤6.7：d″₃₂＝w″，i＝1，w″＝0，i′值加1,如果i′≤5则转步骤6.2，否则转步骤6.8。

步骤6.8：i′＝1，j值加1,如果j≤5则转步骤6.2，否则转步骤6.9。

步骤6.9：前6个测量点中，按照运动模型集合M′_5×5中的每一个元素分别进行多项式拟合后拟合误差的最大值矩阵D″_5×5为：

步骤7：通过D′_5×5和D″_5×5的均值构造模型筛选矩阵D_5×5为：

步骤8：查找模型筛选矩阵D_5×5中值最小的元素d_u,v(u∈[1,m],v∈[1,m])，裁剪并变换得到邻域矩阵B_3×3。根据邻域矩阵中所有值为0的元素对应的运动模型参数生成有效运动模型数组M″_L×2(L∈[1,9])。

步骤8.1：找到模型筛选矩阵D_5×5中值最小的元素d_u,v＝d₃₂＝0.6065(u＝3,v＝2)。

步骤8.2：裁剪D_5×5得到邻域矩阵B_3×3，并进行规则化变换：

取邻域矩阵筛选系数α＝3，表示矩阵D_5×5中大于d₃₂×3＝0.6065×3＝1.8195的元素所在位置对应的组合拟合多项式为非有效运动模型，取值为1。按照取值变换规则，计算结果如下：

步骤8.3：根据邻域矩阵生成有效运动模型数组M″_L×2，L表示数组长度，L∈[1,9]。邻域矩阵B_3×3中值为0的元素(共4个)对应的拟合多项式均为有效的目标运动模型，此处数组下标L＝4。二维数组M″_4×2的第一列表示目标在X轴方向上的有效拟合多项式次数，第二列表示目标在Y轴方向上的有效拟合多项式次数。例如：若b₂₃＝0，则对应的目标在X轴方向上作u＝3次、在Y轴方向上作v+1＝3次多项式拟合而构成的运动模型是目标可能存在的运动方式，此时数组元素M″_l1＝u＝3，M″_l2＝v+1＝3，l∈[1,4]。最终生成的有效运动模型数组M″_4×2为：

步骤9：按照有效运动模型数组M″_4×2中的参数，依据s＝7个测量点数据(见表1.6)，在X、Y两个维度上分别采用基于最小二乘的多项式拟合方法求拟合多项式系数，最终得到目标运动模型估计集合M_4×2。

表1.6：当前测量点数据(s＝7)

序号	时间t(秒)	X(km)	Y(km)
				1	16	1456.98	2452.12
2	26	1454.59	2453.82
				3	36	1452.16	2454.10
4	46	1449.29	2454.01
				5	57	1447.72	2452.01
6	67	1445.48	2450.92
				7	77	1444.71	2447.93

步骤9.1：令l＝1。

步骤9.2：令j＝M″_l1。

步骤9.3：将s个测量点的X轴运动分量{(t_i,X_i)|i＝1,2,...,7}简记为：{(x_i,y_i)||i＝1,2,...,7}，采用基于最小二乘的j次多项式拟合方法，求拟合多项式系数。包括如下步骤：

步骤9.3.1：求

r＝0,1,2,...,j×2，

p＝0,1,2,...,j。其中：测量点个数n＝7，(x_i,y_i)表示第i个测量点的坐标，

表示x_i的p次方。

步骤9.3.2：令J＝j+1，构造矩阵E_J×J，C_J和K_J，其中：

以j＝M″₃₁＝3为例，E_J×J，C_J和K_J分别为：

步骤9.3.3：采用全选主元高斯消去法求J阶线性代数方程组E_J×JK_J＝C_J的解。以j＝M″₃₁＝3为例，解得：

记为：

步骤9.4：令j＝M″_l2。

步骤9.5：将7个测量点的Y轴运动分量{(t_i,Y_i)|i＝1,2,...,7}简记为：{(x_i,y_i)||i＝1,2,...,7}，采用基于最小二乘的j次多项式拟合方法，求拟合多项式系数。包括如下步骤：

步骤9.5.1：同步骤9.3.1。

步骤9.5.2：同步骤9.3.2。

步骤9.5.3：采用全选主元高斯消去法求J阶线性代数方程组E_J×JK_J＝C_J的解，记为：

步骤9.6：l值加1,如果l≤4则转步骤9.2，否则结束程序。

最终我们得到目标运动模型估计集合M_4×2表示为：

其中：ω_l1是目标第l∈[1,4]种运动模型估计中在X轴方向上的拟合多项式系数，ω_l2是对应的其在Y轴方向上的拟合多项式系数，通常情况下ω_l1≠ω_l2，即第l种模型估计中目标在X轴和Y轴方向上分别按照不同的多项式次数和系数值进行拟合。经过计算，M_4×2中各元素分别表示为：

步骤10：使用新的测量点(第8、9点，见表1.7)对本发明的实施效果进行验证。

表1.7：检验验证数据点(第8,9点)

序号	时间t(秒)	X(km)	Y(km)
				8	87	1443.10	2446.18
9	97	1442.16	2443.61

为验证目标运动模型估计集合M_4×2的有效性，我们首先根据已有的7个测量点(见表1.6)，计算所有可能的目标运动模型(同样取m＝5，计算过程与步骤2至步骤4相同)，得到的目标在平面直角坐标系下可能存在的所有运动模型集合

然后计算第8,9测量点(见表1.7)按照

进行拟合后的剩余误差矩阵

通过与步骤8得到的模型筛选矩阵D_5×5的比较，可以确定通过本发明方法得到的目标运动模型估计数量大大减少(本例中由25种减小到4种，对应虚线框部分)的同时，目标运动模型估计集合M_4×2中每一种模型对新测量点的拟合误差都在合理范围内，有效避免了过拟合现象的产生。

第8,9测量点按照

进行拟合后的剩余误差矩阵

的计算结果为：

其中：

的下标值对应本例计算结果目标运动模型估计集合M_4×2中的X轴和Y轴方向上的拟合多项式次数，而且

是

中的最小值，表示通过新点验证，目标在X轴上进行2次、在Y轴上进行3次多项式拟合次数是最佳运动模型估计结果。

将D_5×5和

中的数据转换为一维数组后，显示如图3所示。表明本例中使用当前7点进行模型筛选的误差数据与使用新点进行验证的剩余误差数据保持了非常一致的趋势，二者根据积差法定义的相关系数R为0.9659，属于高度相关。据此也可以说明本发明方法中的模型筛选方法有效。

本发明针对拟合多项式描述目标运动模型时存在的过拟合现象，以及多维度下模型估计数量的指数级增长问题，采用基于交叉检验的方法，依据s-1个测量点数据计算多维度下目标可能存在的所有运动模型，将第s个测量点与旧有测量点特征值的交叉检验结果，作为基于s个测量点数据确定目标运动模型的依据，从而达到筛选出较少且有效目标运动模型的目的。通过对后续2个新测量点的实际检验用例表明，采用本发明确定的多维度下的目标运动模型估计数量大大减少，有效避免了过拟合现象，为在本系统中开展目标跟踪与数据融合工作提供了一种简单、有效且易于工程实现的目标运动模型估计方法。本发明模型估计方法科学、方案实施步骤合理，模型估计结果理想，对于提高后续目标跟踪与数据融合的准确性、降低计算复杂度具有重要意义，同时也有助于目标特征分析、目标趋势预测、目标企图判别等工作的开展。本发明所提供的方法的时间复杂度和空间复杂度都很低，可操作性和实用性很强。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。例如但不限于以下几点：

(1)关于s的取值问题。

在本发明步骤1中，为保证目标运动模型估计的即时性(关注当前)，一般s取值7或8为宜。但在目标特征识别、目标运动趋势分析(关注历史)等工作中，采用本发明方法并增大s值，既能满足特殊功能的需要，同时也有助于提高目标运动模型估计的准确性。

(2)关于m的取值问题。

本发明步骤2中的m表征了目标在某个维度上依赖于时间的运动模型复杂度，对于不同维度，m值是不同的。本例中将X轴方向上和Y轴方向上的模型复杂度设置为相同，既符合飞机等空中目标的运动规律特征，也是为了方案叙述上的方便，实际使用中应根据不同维度的物理意义灵活掌握。

(3)拟合多项式系数求解方法的选择。

本发明采用了最小二乘法求拟合多项式系数，对于契比雪夫曲线拟合方法、最佳一致逼近的里米兹方法等其他拟合多项式系数求解方法同样适用，但要保证其在全过程中的一致性。

(4)合理考虑其他维度的运动模型。

本发明只给出了目标在X轴方向上和Y轴方向上的运动模型估计方法，在实际应用中，还可以考虑目标高度、速度、加速度等多个维度依赖于时间的运动模型估计，只要对相关步骤进行扩充即可。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。

Claims (10)

1.一种基于交叉检验的目标运动模型估计方法，其特征在于，该方法包括如下步骤：

步骤1：将某雷达针对某个空中目标的s个连续测量点数据(ρ_i,θ_i,h_i,t_i)转换为统一直角坐标{(X_i,Y_i,t_i)|i＝1,2,...,s}，其中，(ρ_i,θ_i,h_i,t_i)表示t_i时刻该雷达测得的目标距离ρ_i、方位θ_i和高度h_i，i＝1,2,…s；

步骤2：对目标X轴运动分量{(t_i,X_i)|i＝1,2,...,s-1}，使用基于最小二乘的多项式拟合方法，分别得到最高次数为1,2,...,m时的X轴拟合多项式系数

步骤3：对目标Y轴运动分量{(t_i,Y_i)|i＝1,2,...,s-1}，使用基于最小二乘的多项式拟合方法，分别得到最高次数为1,2,...,m时的Y轴拟合多项式系数

步骤4：依据s-1个测量点数据得到的目标在平面直角坐标系下可能存在的所有运动模型集合M′_m×m，用多项式系数表示为：

其中：

等价为：

表示目标运动模型在X轴方向上按照j次多项式拟合，在Y轴方向上按照i′次多项式拟合；

步骤5：计算第s个测量点(t_s,X_s,Y_s)按照运动模型集合M′_m×m中的每一个元素分别进行多项式拟合后的拟合误差矩阵D′_m×m为

步骤6：求前s-1个测量点中，按照运动模型集合M′_m×m中的每一个元素分别进行多项式拟合后拟合误差的最大值矩阵D″_m×m为

步骤7：通过D′_m×m和D″_m×m的均值构造模型筛选矩阵D_m×m为：

步骤8：查找所述模型筛选矩阵D_m×m中值最小的元素d_u,v(u∈[1,m],v∈[1,m])，裁剪并变换得到邻域矩阵B_3×3，根据邻域矩阵中所有值为0的元素对应的运动模型参数生成有效运动模型数组M″_L×2(L∈[1,9])；

步骤9：按照有效运动模型数组M″_L×2中的参数，依据s个测量点数据，在X、Y两个维度上分别采用基于最小二乘的多项式拟合方法求拟合多项式系数，最终得到目标运动模型估计集合M_L×2。

2.如权利要求1所述的基于交叉检验的目标运动模型估计方法，其特征在于，所述步骤1具体包括如下步骤：

步骤1.1：获取雷达观测数据(ρ_i,θ_i,h_i,t_i)，i＝1,2,…s；

步骤1.2：将雷达观测数据(ρ_i,θ_i,h_i,t_i)，i＝1,2,…s转换为以本雷达站为中心的二维直角坐标(X_zi,Y_zi)：

步骤1.3：将(X_zi,Y_zi)，i＝1,2,…s转换为中心统一直角坐标(X_i,Y_i)：

X_i＝X_zi cosδ_xz-Y_zi sinδ_xz+X_zx

Y_i＝X_zi sinδ_xz+Y_zi cosδ_xz+Y_zx

其中：(X_zx,Y_zx)为雷达站址在中心统一直角坐标系中的坐标；δ_xz为雷达站址与直角坐标系中心点的经度差，单位为弧度。

3.如权利要求1或2所述的基于交叉检验的目标运动模型估计方法，其特征在于，所述步骤2具体包括如下步骤：

步骤2.1：将X轴运动分量{(t_i,X_i)|i＝1,2,...,s-1}简记为：{(x_i,y_i)||i＝1,2,...,n}，然后确定m值，m表示在X轴方向上目标运动模型以拟合多项式表示时的最高次数；

步骤2.2：针对n个测量点数据{(x_i,y_i)||i＝1,2,...,n}，令j依次取值1,2,...,m，采用基于最小二乘的j次多项式拟合方法，求得X轴方向上m组拟合多项式系数，记为：

4.如权利要求1或2所述的基于交叉检验的目标运动模型估计方法，其特征在于，所述步骤3具体包括如下步骤：

步骤3.1：将Y轴运动分量{(t_i,Y_i)|i＝1,2,...,s-1}简记为：{(x_i,y_i)||i＝1,2,...,n}，然后确定m值,m表示在Y轴方向上目标运动模型以拟合多项式表示时的最高次数；

步骤3.2：针对n个测量点数据{(x_i,y_i)||i＝1,2,...,n}，令j依次取值1,2,...,m，采用基于最小二乘的j次多项式拟合方法，求得Y轴方向上m组拟合多项式系数，记为：

5.如权利要求1或2所述的基于交叉检验的目标运动模型估计方法，其特征在于，所述步骤5具体包括如下步骤：

步骤5.1：令j＝1，i′＝1；

步骤5.2：取运动模型集合M′_m×m的元素m′_ji′，得到X轴上的拟合多项式系数

求第s个测量点(t_s,X_s,Y_s)在X轴向的拟合误差

步骤5.3：取运动模型集合M′_m×m的元素m′_ji′，得到Y轴上的拟合多项式系数

求第s个测量点(t_s,X_s,Y_s)在Y轴向的拟合误差

步骤5.4：则第s个测量点按照运动模型m′_ji′，在平面直角坐标系中的拟合误差为

其中j∈[1,m],i′∈[1,m]；

步骤5.5：i′值加1，如果i′≤m则转步骤5.2，否则转步骤5.6；

步骤5.6：i′＝1，j值加1,如果j≤m则转步骤5.2，否则转步骤5.7；

步骤5.7：则第s个测量点按照运动模型集合M′_m×m进行多项式拟合后的拟合误差矩阵为：

6.如权利要求1或2所述的基于交叉检验的目标运动模型估计方法，其特征在于，所述步骤6具体包括如下步骤：

步骤6.1：令i＝1，j＝1，i′＝1，w″＝0；

步骤6.2：取运动模型集合M′_m×m的元素m′_ji′，得到X轴上的拟合多项式系数

求第i个测量点(X_i,Y_i,t_i)在X轴向的拟合误差

步骤6.3：取运动模型集合M′_m×m的元素m′_ji′，得到Y轴上的拟合多项式系数

求第i个测量点(X_i,Y_i,t_i)在Y轴向的拟合误差

步骤6.4：则第i个测量点按照运动模型m′_ji′，在平面直角坐标系中的拟合误差为

其中i∈[1,s-1],j∈[1,m],i′∈[1,m]；

步骤6.5：记录当前拟合误差最大值：如果

则

步骤6.6：i值加1，如果i≤s-1则转步骤6.2，否则转步骤6.7；

步骤6.7：d″_ji′＝w″，i＝1，w″＝0，i′值加1,如果i′≤m则转步骤6.2，否则转步骤6.8；

步骤6.8：i′＝1，j值加1,如果j≤m则转步骤6.2，否则转步骤6.9；

步骤6.9：前s-1个测量点中，按照运动模型集合M′_m×m中的每一个元素分别进行多项式拟合后拟合误差的最大值矩阵D″_m×m为：

7.如权利要求1或2所述的基于交叉检验的目标运动模型估计方法，其特征在于，所述步骤8具体包括如下步骤：

步骤8.1：找到模型筛选矩阵D_m×m中值最小的元素d_u,v(u∈[1,m],v∈[1,m])；

步骤8.2：裁剪D_m×m得到邻域矩阵B_3×3，并进行规则化变换：

其中：b_g,h(g∈[1,3],h∈[1,3])，取值变换规则如下：

其中α为邻域矩阵筛选系数，α≥1，表示大于矩阵D_m×m中最小值元素α倍的元素所在位置对应的组合拟合多项式为非有效运动模型；

步骤8.3：根据邻域矩阵生成有效运动模型数组M″_L×2，L表示数组长度，L∈[1,9]；二维数组M″_L×2的第一列表示目标在X轴方向上的有效拟合多项式次数，第二列表示目标在Y轴方向上的有效拟合多项式次数，均是大于等于1，小于等于m的自然数。

8.如权利要求1或2所述的基于交叉检验的目标运动模型估计方法，其特征在于，所述步骤9具体包括如下步骤：

步骤9.1：令l＝1；

步骤9.2：令j＝M_l″₁；

步骤9.3：将s个测量点的X轴运动分量{(t_i,X_i)|i＝1,2,...,s}简记为：{(x_i,y_i)||i＝1,2,...,n}，采用基于最小二乘的j次多项式拟合方法，求拟合多项式系数；

步骤9.4：令j＝M″_l2；

步骤9.5：将s个测量点的Y轴运动分量{(t_i,Y_i)|i＝1,2,...,s}简记为：{(x_i,y_i)||i＝1,2,...,n}，采用基于最小二乘的j次多项式拟合方法，求拟合多项式系数；

步骤9.6：l值加1,如果l≤L则转步骤9.2，否则结束程序；

最终我们得到目标运动模型估计集合M_L×2表示为：

其中：L∈[1,9]，ω_l1是目标第l∈[1,L]种运动模型估计中在X轴方向上的拟合多项式系数，ω_l2是对应的其在Y轴方向上的拟合多项式系数，即第l种模型估计中目标在X轴和Y轴方向上分别按照不同的多项式次数和系数值进行拟合，不存在过拟合现象。

9.如权利要求8所述的基于交叉检验的目标运动模型估计方法，其特征在于，所述步骤9.3具体包括如下步骤：

步骤9.3.1：求

其中：n表示测量点个数，(x_i,y_i)表示第i个测量点的坐标，

表示x_i的p次方；

步骤9.3.2：令J＝j+1，构造矩阵E_J×J，C_J和K_J，其中：

步骤9.3.3：采用全选主元高斯消去法求J阶线性代数方程组E_J×JK_J＝C_J的解，记为：

10.如权利要求8所述的基于交叉检验的目标运动模型估计方法，其特征在于，所述步骤9.5具体包括如下步骤：

步骤9.5.1：求

其中：n表示测量点个数，(x_i,y_i)表示第i个测量点的坐标，

表示x_i的p次方；

步骤9.5.2：令J＝j+1，构造矩阵E_J×J，C_J和K_J，其中：

步骤9.5.3：采用全选主元高斯消去法求J阶线性代数方程组E_J×JK_J＝C_J的解，记为：

Images