CN109856622A - 一种约束条件下的单雷达直线航迹线目标状态估计方法 - Google Patents

Abstract

本发明属于传感器目标跟踪与数据融合技术领域，具体涉及一种约束条件下的单雷达直线航迹线目标状态估计方法。该方法针对目标由匀速圆周运动转变为匀速直线运动后单雷达直线航迹线，将目标在匀速圆周运动模型中的先验信息作为当前匀速直线运动模型的约束条件，在统一直角坐标系中对当前有限个测量点进行垂直距离加权估计得到与时间无关的目标运动方程和航向估值，通过目标时间‑路程加权匀速估计模型得到线速度和分量速度估值，进而最终确定目标状态位置估计结果。对比试验表明，该方法不但易于工程实现，而且比传统卡尔曼滤波方法具有更优的目标状态估计效果和收敛速度。

Description

一种约束条件下的单雷达直线航迹线目标状态估计方法

技术领域

本发明属于传感器目标跟踪与数据融合技术领域，具体涉及在复杂数据环境中，且系统噪声和观测噪声均未知的情况下，针对目标由匀速圆周运动转变为匀速直线运动的一种约束条件下的单雷达直线航迹线目标状态估计方法。

背景技术

目标跟踪是指利用传感器所获得的对运动目标的时序离散观测数据，对目标的运动状态(位置、航向、速度等)进行估计和跟踪的方法。卡尔曼滤波方法是现阶段雷达中最常用的、占据主导地位的目标跟踪算法。它能够针对不同运动模型，通过累积测量对目标运动状态进行最小方差估计和跟踪，具有较好的跟踪精度，但通常收敛速度较慢，从开始滤波器建立到稳定跟踪一般所用时间为将近10个采样周期。而在通常情况下，当目标处于匀速直线飞行状态时，采用适合线性运动模型的标准卡尔曼滤波方法(KF)；当目标沿圆弧航迹线匀速飞行时，采用适合非线性运动模型的扩展卡尔曼滤波方法(EKF)。在每次针对不同运动模型的卡尔曼滤波器建立初期(一分钟以内)，都是目标运动状态的不稳定跟踪期，而对于飞机等高速、高机动目标，滤波收敛快慢直接影响到目标跟踪的稳定度和对目标的识别速度。

另外，在多雷达组网探测与目标跟踪系统(以下简称“系统”)中还存在实际系统数据环境复杂(测量数据采样间隔较大、测量误差、坐标转换误差以及数据标定与传输误差同时存在)与卡尔曼滤波方法应用条件比较苛刻的矛盾；卡尔曼滤波方法应用中过程噪声协方差Q和量测噪声协方差R在参数值选取上存在困难，不能兼顾稳态性和瞬时性的问题；系统测量中不可避免的野值出现，导致卡尔曼滤波发散或目标状态估计在短时间内失效，表现为估值航迹线出现较大锯齿、目标航向或速度估计出现突变等。

发明内容

(一)要解决的技术问题

本发明要解决的技术问题是：如何提供一种约束条件下的单雷达直线航迹线目标状态估计方法。

(二)技术方案

为解决上述技术问题，本发明提供一种约束条件下的单雷达直线航迹线目标状态估计方法，所述估计方法应用于单雷达数据跟踪与多雷达数据融合系统的前期数据预处理过程中；所述估计方法包括如下步骤：

步骤1：将某雷达针对某个由匀速圆周运动转变为匀速直线运动后直线段上目标的n个观测点数据(ρ_i,θ_i,h_i,t_i)转换为统一直角坐标(x_i,y_i,t_i)；其中，(ρ_i,θ_i,h_i,t_i)表示t_i时刻该雷达测得的目标距离ρ_i、方位θ_i和高度h_i，i＝1,2,…n；

步骤2：对n个观测点数据{(x_i,y_i),i＝1,2,…n}使用约束条件下的不加权直线参数估计模型粗略估计目标运动状态方程y-k₁x-d₁＝0；

步骤3：将方程y-k₁x-d₁＝0作为初始状态，对n个观测点数据{(x_i,y_i),i＝1,2,…n}使用约束条件下的加权直线参数估计模型进行迭代估计，得到最佳目标运动状态方程y-kx-d＝0，并使用取点定向法计算出t_n时刻的目标航向K_n；

步骤4：使用时间-路程加权匀速估计模型迭代估计目标最佳线速度方程S-Vt-S₀＝0，得到t_n时刻的目标速度V_n并计算出目标在X、Y方向上的分量速度v_xn和v_yn：

步骤5：计算n个观测点数据{(x_i,y_i,t_i),i＝1,2,…n}到直线y-kx-d＝0的垂足点的X轴坐标均值和时间均值并令则t_n时刻目标在统一直角坐标系中的位置估值(P_xn,P_yn)为：

由此得到目标由匀速圆周运动转变为匀速直线运动后，在t_n时刻的位置(P_xn,P_yn)、速度V_n和航向K_n。

(三)有益效果

与现有技术相比较，本发明提供一种适应多雷达组网探测与目标跟踪系统复杂数据环境，针对目标由匀速圆周运动转变为匀速直线运动后单雷达直线航迹线的目标状态估计方法。该方法将目标在匀速圆周运动模型中的先验信息作为当前匀速直线运动模型的约束条件，在统一直角坐标系中对当前有限个测量点进行垂直距离加权估计得到与时间无关的目标航向估值，通过目标时间-路程加权匀速估计得到线速度和分量速度估值，进而最终确定目标状态位置估计结果。对比试验表明，该方法不但易于工程实现，而且比传统卡尔曼滤波方法具有更优的目标状态估计效果和收敛速度。

附图说明

图1为本发明技术方案中目标状态估计方法的流程图。

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

图3为本发明实施例2中目标真实位置、卡尔曼估值位置和本发明估值位置三者比较显示图；

图4为本发明实施例2中目标测量位置、卡尔曼估值位置和本发明估值位置三者比较显示图；

图5为本发明实施例2中使用卡尔曼方法和使用本发明方法得到的所有时刻的时间配准点间距显示图；

图6为本发明实施例2中使用卡尔曼方法和使用本发明方法得到的所有时刻的平均时间配准点间距显示图；

图7为本发明实施例2中目标真实速度、卡尔曼估计速度和本发明估计速度三者比较显示图；

图8为本发明实施例2中使用卡尔曼方法和使用本发明方法得到的各点速度差显示图；

图9为本发明实施例2中使用卡尔曼方法和使用本发明方法得到的平均速度差随测量点数变化情况显示图；

图10为本发明实施例2中目标真实航向、卡尔曼估计航向和本发明估计航向三者比较显示图；

图11为本发明实施例2中使用卡尔曼方法和使用本发明方法得到的各点航向差显示图；

图12为本发明实施例2中使用卡尔曼方法和使用本发明方法得到的平均航向差随测量点数变化情况显示图；

图13为本发明实施例2中含有野值时目标真实位置与测量位置在统一直角坐标系中的显示图；

图14为本发明实施例2中含有野值时目标真实位置、卡尔曼估值位置和本发明估值位置三者比较显示图；

图15为本发明实施例2中含有野值时目标测量位置、卡尔曼估值位置和本发明估值位置三者比较显示图；

图16为本发明实施例2中含有野值时使用卡尔曼方法和使用本发明方法得到的所有时刻的时间配准点间距显示图；

图17为本发明实施例2中含有野值时使用卡尔曼方法和使用本发明方法得到的所有时刻的平均时间配准点间距显示图；

图18为本发明实施例2中含有野值时目标真实速度、卡尔曼估计速度和本发明估计速度三者比较显示图；

图19为本发明实施例2中含有野值时使用卡尔曼方法和使用本发明方法得到的各点速度差显示图；

图20为本发明实施例2中含有野值时使用卡尔曼方法和使用本发明方法得到的平均速度差随测量点数变化情况显示图；

图21为本发明实施例2中含有野值时目标真实航向、卡尔曼估计航向和本发明估计航向三者比较显示图；

图22为本发明实施例2中含有野值时使用卡尔曼方法和使用本发明方法得到的各点航向差显示图；

图23为本发明实施例2中含有野值时使用卡尔曼方法和使用本发明方法得到的平均航向差随测量点数变化情况显示图；

图24为本发明技术方案中用到的取点定向法实现流程图。

具体实施方式

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

针对以上现有技术的问题，本发明通过对目标由匀速圆周运动转变为匀速直线运动初期时运动规律的研究，将目标在匀速圆周运动模型中的先验信息作为当前匀速直线运动模型的约束条件，最终目的就是基于系统复杂数据环境，在无须获取系统噪声和观测噪声等先验知识的情况下，有效提高约束条件下的单雷达直线航迹线目标状态估计精度和收敛速度。通过与传统卡尔曼滤波算法在目标位置、航向和速度估值上进行对比实验表明，本发明取得了更好的估值效果和跟踪收敛速度，为在本系统中开展目标跟踪与数据融合工作提供了一种稳定、有效且易于工程实现的滤波方法。

具体而言，为解决现有技术的问题，本发明提供一种约束条件下的单雷达直线航迹线目标状态估计方法，如图1所示，所述估计方法应用于单雷达数据跟踪与多雷达数据融合系统的前期数据预处理过程中；所述估计方法包括如下步骤：

步骤1：将某雷达针对某个由匀速圆周运动转变为匀速直线运动后直线段上目标的n个观测点数据(ρ_i,θ_i,h_i,t_i)转换为统一直角坐标(x_i,y_i,t_i)；其中，(ρ_i,θ_i,h_i,t_i)表示t_i时刻该雷达测得的目标距离ρ_i、方位θ_i和高度h_i，i＝1,2,…n；此处要保证n个观测点数据均处于目标匀速直线运动期间；

步骤2：对n个观测点数据{(x_i,y_i),i＝1,2,…n}使用约束条件下的不加权直线参数估计模型粗略估计目标运动状态方程y-k₁x-d₁＝0；

步骤3：将方程y-k₁x-d₁＝0作为初始状态，对n个观测点数据{(x_i,y_i),i＝1,2,…n}使用约束条件下的加权直线参数估计模型进行迭代估计，得到最佳目标运动状态方程y-kx-d＝0，并使用取点定向法计算出t_n时刻的目标航向K_n；

步骤4：使用时间-路程加权匀速估计模型迭代估计目标最佳线速度方程S-Vt-S₀＝0，得到t_n时刻的目标速度V_n并计算出目标在X、Y方向上的分量速度v_xn和v_yn：

步骤5：计算n个观测点数据{(x_i,y_i,t_i),i＝1,2,…n}到直线y-kx-d＝0的垂足点的X轴坐标均值和时间均值并令则t_n时刻目标在统一直角坐标系中的位置估值(P_xn,P_yn)为：

由此得到目标由匀速圆周运动转变为匀速直线运动后，在t_n时刻的位置(P_xn,P_yn)、速度V_n和航向K_n。

其中，所述步骤1包括如下步骤：

步骤1.1：选择某雷达观测数据集{(ρ_i,θ_i,h_i,t_i),i＝1,2,…n}，是针对某个由匀速圆周运动转变为匀速直线运动后直线段上目标的时序观测数据；如果是三坐标雷达观测数据，一般认为h_i保持不变；如果是两坐标雷达观测数据，则认为h_i＝0。θ_i取值范围[0,360)，单位为度，其中雷达正北方向为0度，正东方向为90度。

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

步骤1.3：将(X_zi,Y_zi)，i＝1,2,…n转换为中心统一直角坐标(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.1中，t_i＜t_i+1，且n不大于N，N表示参与目标状态估计的最近时间内有限个观测数据点的个数，N取4～15。

其中，所述步骤2包括如下步骤：

步骤2.1：用所有观测点{(x_i,y_i),i＝1,2,…n}到待求直线的距离l_i的平方和最小作为估计模型，用待求直线航迹线必须与已知圆相切作为约束条件，计算在此条件下的这条直线的最佳参数(k,d)；其中，k为直线的斜率，d为该直线在x轴上的截距；

其中的已知圆就是目标在进入匀速直线运动前所处的圆周运动轨迹，为已知变量，其圆心为半径为r；

所述约束条件下的不加权直线参数估计模型如下：

由(2)式解得：

(3)式代入(1)式得：

对于(4)式，应有f(k)对于k的导数等于零，即：

记

则(5)式化简为：

为了消去“根号”，将(6)改写成：

(7)式两边平方后得：

(a₃k²+(a₂-b₂)k-a₃)²＝(a₁r+b₁rk)²(1+k²) (8)

(8)式两边展开并整理后得：

记：

则(9)式可化简为：

k⁴+bk³+ck²+fk+e＝0 (10)

在此采用求实系数代数方程全部根的牛顿-下山法(《C常用算法程序集》-徐士良编著，清华大学出版社，1996)对该式进行求解，得到k₁、k₂、k₃和k₄四个根；将k₁、k₂、k₃、k₄带入(3)式，对应得到d₁、d₂、d₃、d₄；

步骤2.2：按照距离最小原则确定模型的合理解；

(k₁,d₁)、(k₂,d₂)、(k₃,d₃)和(k₄,d₄)都是约束条件下的不加权直线参数估计模型的实根；按照点{(x_i,y_i)，i＝1,2,…n}到所求直线的距离的平方和最小原则，确定合理的直线参数值；该问题也可以简化为：计算点(x_n,y_n)分别到直线y＝k₁×x+d₁、y＝k₂×x+d₂、y＝k₃×x+d₃和y＝k₄×x+d₄的距离l₁,l₂,l₃,l₄：

取|l₁|、|l₂|、|l₃|和|l₄|中最小的|l_i|所对应的(k_i,d_i)，作为所求直线的合理参数，记为(k₁,d₁)。

其中，所述步骤3包括如下步骤：

步骤3.1：计算各点(x_i,y_i)到直线y-k₁x-d₁＝0的距离|l_i|之和；

n为数据点数；

步骤3.2：求各点(x_i,y_i)到直线y-k₁x-d₁＝0的距离l_i；

式中m表示迭代次数，n表示数据点数；m初始值为1，即：k⁽¹⁾＝k₁，d⁽¹⁾＝d₁；

步骤3.3：求|l_i|的倒数；

步骤3.4：求各点的权值w_i；

步骤3.5：用所有观测点{(x_i,y_i),i＝1,2,…n}到待求直线的加权距离w_il_i的平方和最小作为估计模型，用待求直线航迹线必须与已知圆相切作为约束条件，计算在此条件下的这条直线的最佳参数(k,d)；

其中的已知圆就是目标在进入匀速直线运动前所处的圆周运动轨迹，为已知变量，其圆心为半径为r；

所述约束条件下的加权直线参数估计模型如下：

由(12)式得：

(13)式代入(11)得：

对于(14)式，应有f(k)对于k的导数等于零，即：

记： 则(15)式化简为：

为了消去“根号”，将(16)改写成：

(17)式两边平方后得：

(a'₃k²+(a'₂-b'₂)k-a'₃)²＝(a'₁r+b'₁rk)²(1+k²) (18)

(18)式两边展开并整理后得：

记：

则(19)式可化简为：

k⁴+b'k³+c'k²+f'k+e'＝0 (20)

在此采用求实系数代数方程全部根的牛顿-下山法(《C常用算法程序集》-徐士良编著，清华大学出版社，1996)对该式进行求解，得到k'₁、k'₂、k'₃和k'₄四个根；将k'₁、k'₂、k'₃、k'₄带入(13)式，对应得到d'₁、d'₂、d'₃、d'₄；

步骤3.6：按照距离最小原则确定方程的合理解；

(k'₁,d'₁)、(k'₂,d'₂)、(k'₃,d'₃)和(k'₄,d'₄)都是模型的实根；我们按照点{(x_i,y_i)，i＝1,2,…n}到所求直线的距离的平方和最小原则，确定合理的直线参数值；该问题也可以简化为：计算点(x_n,y_n)分别到直线y＝k'₁×x+d'₁、y＝k'₂×x+d'₂、y＝k'₃×x+d'₃和y＝k'₄×x+d'₄的距离l'₁,l'₂,l'₃,l'₄：

m加1，并取|l'₁|、|l'₂|、|l'₃|和|l'₄|中最小的|l'_i|所对应的(k'_i,d'_i)，作为所求直线的合理参数，记为(k^(m),d^(m))，m表示迭代次数；

步骤3.7：计算所有数据点到新直线y-k^(m)x-d^(m)＝0的加权距离之和f^(m)(k^(m),d^(m))；

式中m表示迭代次数，n表示数据点数；

步骤3.8：判别是否为“最佳”解；

若f^(m)(k^(m),d^(m))≥f^(m-1)(k^(m-1),d^(m-1))或f^(m)(k^(m),d^(m))≤10^-6，则输出解(k^(m-1),d^(m-1))，并简记为(k,d)；否则重复步骤3.2-步骤3.8；式中m表示迭代次数；

步骤3.9：使用取点定向法计算出t_n时刻的目标航向K_n。

其中，所述步骤3.9中取点定向法的实现过程为：设(k,d)是估计得到的雷达观测航迹线参数，(x₁,y₁)和(x_n,y_n)是雷达对目标的首末测量点，经过滤波后，这两个点的坐标变为(x₁,y′₁)和(x_n,y′_n)，其中：y′₁＝k×x₁+d，y′_n＝k×x_n+d；令：Δx＝x_n-1-x₁，Δy＝y′_n-1-y′₁，π为圆周率，接下来依次进行如下判断和计算；

①如果Δy等于0，转②，否则转③；

②如果Δx大于0，航向K取值为0度，否则，航向K取值为180度，程序结束；

③如果Δx等于0，转④，否则转⑤；

④如果Δy大于等于0，航向K取值为90度，否则，航向K取值为270度，取点定向结束；

⑤如果Δy大于0，航向K取值为度，否则，航向K取值为度，取点定向结束。

取点定向法计算机实现流程如图24所示。航向K_n正北方向为0度，正东方向为90度。

其中，所述步骤4包括如下步骤：

步骤4.1：选取n个测量点{(x_i,y_i,T_i),i＝1,2,…n}作为时间-路程加权匀速估计模型的计算点；

(1)当直线航迹段的实际测量点数m′小于步骤1中确定的目标观测点数n时，需要使用目标处于圆弧运动期的n-m′个测量点，此时这n个点分别对应的时刻按照递增顺序表示为：t_m'-n+1,t_m'-n+2,…,t_-1,t₀,t₁,t₂,…,t_m'，其中：t_m'-n+1,…,t_-1,t₀为圆弧航迹段上的测量点时刻，t₁,t₂,…,t_m'为直线航迹段上的测量点时刻；

(2)当直线航迹段的实际测量点数m′不小于步骤1中确定的目标观测点数n时，这n个点均为直线航迹段上的测量点，分别对应的时刻按照递增顺序表示为t_m'-n+1,t_m'-n+2,…,t_m'；

令：

T₁＝t_m'-n+1，

T₂＝t_m'-n+2，

……

T_n＝t_m'；

并取T_i时刻对应的测量点坐标(x_i,y_i)，构成{(x_i,y_i,T_i),i＝1,2,…n}作为时间-路程加权匀速估计模型的计算点；

步骤4.2：计算T_i时刻目标在航迹线上的垂点(cx_i,cy_i),i＝1,2,…n；

(1)当(x_i,y_i,T_i)为直线航迹段(y-kx-d＝0)上的测量点时；

令：A＝k，B＝-1，C＝d，则：T_i时刻(x_i,y_i)在直线y-kx-d＝0上的垂点(cx_i,cy_i)为：

(2)当(x_i,y_i,T_i)为圆心为半径为r的圆弧航迹段上的测量点时；

圆弧航迹线标准方程为：

式中圆心和半径r已知；

过观测点(x_i,y_i)与圆心的直线方程为：

解(22)式得：

将(23)式代入(21)式得：

解(24)式得：

代入(23)式得：

取(x₁,y₁)和(x₂,y₂)中距离测量点(x_i,y_i)最近的点作为合理解，记为(cx_i,cy_i)；

步骤4.3：计算从T₁时刻起，到T_i时刻，目标经过的路程S_i,i＝1,2,…n；

步骤4.3.1：求T_i时刻与上一计算点T_i-1时刻间的路程ΔS_i，i＝1时，ΔS_i＝0；

(1)当(x_i,y_i,T_i)为直线航迹段上的测量点时；

(2)当(x_i,y_i,T_i)为圆心为半径为r的圆弧航迹段上的测量点时；

求两相邻计算点之间的弦长：

求圆心到弦的距离：

求圆心角：

求圆心角θ_i的对应的弧长即为ΔS_i：

ΔS_i＝rθ_i

步骤4.3.2：计算从T₁时刻起，到T_i时刻，目标经过的路程S_i；

步骤4.4：使用时间-路程加权匀速估计模型迭代估计目标最佳线速度方程S-Vt-S₀＝0；

该步骤4.4包括：

步骤4.4.1：对n个计算点的数据{(T_i,S_i),i＝1,2,…n}使用时间-路程不加权匀速估计模型粗略估计目标线速度方程S-V₀₁t-S₀₁＝0；

时间-路程不加权匀速估计模型为：

式中S_i是通过步骤4.3.2计算出来的对应T_i时刻目标经过的路程，S₀为待求的调节路程，v为待求的匀速度；

接下来对模型(25)进行求解。

f(S₀,v)分别对S₀和v求偏导数并令其等于零，以t_i表示T_i，则有：

(26)、(27)式化简得：

由(28)式解得：

(30)式代入(29)式得：

解(31)式得：

(32)式代入(30)式得：

由以上推导可知，不加权匀速度估计模型求解过程如下；

令：则有：

并将(v,S₀)记为(V₀₁,S₀₁)；

步骤4.4.2：计算各点(T_i,S_i)到S-V₀₁t-S₀₁＝0的距离|l_i|之和；

n为数据点数，V⁽¹⁾＝V₀₁，

步骤4.4.3：求各点(T_i,S_i)到S-V₀₁t-S₀₁＝0的距离l_i；

式中m表示迭代次数，n表示数据点数；m初始值为1，即：V⁽¹⁾＝V₀₁，

步骤4.4.4：求|l_i|的倒数；

步骤4.4.5：求各点的权值w_i；

步骤4.4.6：对n个计算点{(T_i,S_i),i＝1,2,…n}使用时间-路程加权匀速估计模型估计目标线速度方程S-Vt-S₀＝0；

时间-路程加权匀速估计模型为：

式中S_i是通过步骤4.3.2计算出来的对应T_i时刻目标经过的路程，w_i是通过步骤4.4.5计算出来的权值，S₀为待求的调节路程，v为待求的匀速度；

接下来对模型(34)进行求解。

f(S₀,v)分别对S₀和v求偏导数并令其等于零，则有：

(35)、(36)式化简得：

由(37)式解得：

(39)式代入(38)式得：

解(40)式得：

令：

则：

将(v,S₀)记为m表示迭代次数；

步骤4.4.7：m加1，计算所有数据点{(T_i,S_i),i＝1,2,…n}到新直线的加权距离之和

式中m表示迭代次数，n表示数据点数；

步骤4.4.8：判别是否为“最佳”解；

若或则输出解并简记为(V,S₀)；否则重复步骤4.4.3至步骤4.4.8；

步骤4.5：计算t_n时刻目标在X、Y方向上的分量速度v_xn和v_yn；

式中n表示数据点数，t_n表示末点时刻；V_n表示末点速度，即步骤4.4.8求得的V；k表示直线航迹线的斜率，在步骤3.8中已求出；

其中，所述步骤5包括如下步骤：

步骤5.1：求步骤4.2得到的垂点{(T_i,cx_i,cy_i),i＝1,2,…n}的X轴坐标均值和时间均值

步骤5.2：求垂点的Y轴均值

其中：(k,d)表示直线航迹线的斜率和截距，在步骤3.8中已求出；

步骤5.3：求t_n时刻目标在统一直角坐标系中的位置坐标(P_xn,P_yn)；

至此得到了目标由匀速圆周运动转变为匀速直线运动后，在t_n时刻的位置(P_xn,P_yn)(步骤5.3)、速度V_n(步骤4.4.8)和航向K_n(步骤3.9)估值。

实施例1

本实施例具体描述本发明所提出的一种约束条件下的单雷达直线航迹线目标状态估计方法，所述估计方法应用于单雷达数据跟踪与多雷达数据融合系统的前期数据预处理过程中。

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

步骤1：将某雷达针对某个由匀速圆周运动转变为匀速直线运动后直线段上目标的n个观测点数据(ρ_i,θ_i,h_i,t_i)(表示t_i时刻该雷达测得的目标距离ρ_i、方位θ_i和高度h_i，i＝1,2,…n。)转换为统一直角坐标(x_i,y_i,t_i)。此处要保证n个观测点数据均处于目标匀速直线运动期间。

步骤1.1：选择某雷达观测数据集{(ρ_i,θ_i,h_i,t_i),i＝1,2,…n}，是针对某个由匀速圆周运动转变为匀速直线运动后直线段上目标的时序观测数据，要求t_i＜t_i+1。本例中n取值为7，因为是两坐标雷达观测数据，认为h_i＝0。相关数据见表1.1～表1.2。

表1.1：基本参数

表1.2：雷达观测数据

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

计算结果见表1.3。

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

序号	时间t(秒)	X<sub>z</sub>(km)	Y<sub>z</sub>(km)
				1	849	-30.78	-356.85
2	860	-28.99	-355.85
				3	871	-23.94	-354.57
4	881	-19.82	-353.49
				5	890	-17.2	-352.16
6	899	-12.58	-351.21
				7	911	-10.22	-349.66

步骤1.3：将(X_zi,Y_zi)，i＝1,2,…7转换为中心统一直角坐标(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.4。

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

序号	时间t(秒)	x(km)	y(km)
				1	849	1469.22	2193.15
2	860	1471.01	2194.15
				3	871	1476.06	2195.43
4	881	1480.18	2196.51
				5	890	1482.80	2197.84
6	899	1487.42	2198.79
				7	911	1489.78	2200.34

步骤2：对n个观测点数据{(x_i,y_i),i＝1,2,…n}使用“约束条件下的不加权直线参数估计模型”粗略估计目标运动状态方程y-k₁x-d₁＝0。包括如下步骤：

步骤2.1：用所有观测点{(x_i,y_i),i＝1,2,…n}到待求直线的距离平方和最小作为估计模型，用待求直线航迹线必须与已知圆相切作为约束条件，计算在此条件下的这条直线的最佳参数(k,d)。其中，k为直线的斜率，d为该直线在x轴上的截距；其中已知圆圆心半径r如表1.1所示。对(1)、(2)式的求解步骤如下。

①计算a₁,a₂,a₃,b₁,b₂(n＝7)。

②计算b,c,f,e。

③采用求实系数代数方程全部根的牛顿-下山法解方程k⁴+bk³+ck²+fk+e＝0，得到所有解k₁、k₂、k₃和k₄。将k₁、k₂、k₃、k₄带入(3)式，对应得到d₁、d₂、d₃、d₄。计算结果见表1.5。

表1.5：方程的所有解

步骤2.2：按照距离最小原则确定模型的合理解。

计算测量点(x₇,y₇)＝(1489.78，2200.34)分别到直线y＝k₁×x+d₁、y＝k₂×x+d₂、y＝k₃×x+d₃和y＝k₄×x+d₄的距离l₁,l₂,l₃,l₄：

取|l₁|、|l₂|、|l₃|和|l₄|中最小的|l₄|所对应的(k₄,d₄)＝(0.345903，1684.723792)，作为所求直线的合理参数，记为(k₁,d₁)。

步骤3：将方程y-k₁x-d₁＝0作为初始状态，对7个观测点数据{(x_i,y_i),i＝1,2,…7}使用“约束条件下的加权直线参数估计模型”进行迭代估计，得到最佳目标运动状态方程y-kx-d＝0，并使用“取点定向法”计算出t₇时刻的目标航向K₇。包括如下步骤：

步骤3.1：计算各点(x_i,y_i)到直线y-k₁x-d₁＝0的距离|l_i|之和。

步骤3.2：求各点(x_i,y_i)到直线y-k₁x-d₁＝0的距离l_i。

式中m表示迭代次数，n表示数据点数。m初始值为1，即：k⁽¹⁾＝k₁，d⁽¹⁾＝d₁。计算结果见表1.6。

步骤3.3：求|l_i|的倒数。计算结果见表1.6。

步骤3.4：求各点的权值w_i。计算结果见表1.6。

表1.6:距离和权值计算结果

序号	l<sub>i</sub>	P<sub>i</sub>	w<sub>i</sub>
				1	0.207934	4.809229	0.158118
2	0.567144	1.763220	0.057971
				3	0.125351	7.977578	0.262288
4	-0.206875	4.833839	0.158928
				5	0.200250	4.993769	0.164186
6	-0.414357	2.413377	0.079347
				7	0.275911	3.624357	0.119162

步骤3.5：用所有观测点{(x_i,y_i),i＝1,2,…n}到待求直线的加权距离的平方和最小作为估计模型，用待求直线航迹线必须与已知圆相切作为约束条件，计算在此条件下的这条直线的最佳参数(k,d)。其中已知圆圆心半径r如表1.1所示。对(11)、(12)式的求解步骤如下。

①计算a'₁,a'₂,a'₃,b'₁,b'₂(n＝7)。

②计算b',c',f',e'。

③采用求实系数代数方程全部根的牛顿-下山法解方程k⁴+b'k³+c'k²+f'k+e'＝0，得到所有解k'₁、k'₂、k'₃和k'₄。将k'₁、k'₂、k'₃、k'₄带入(3)式，对应得到d'₁、d'₂、d'₃、d'₄。计算结果见表1.7。

表1.7：方程的所有解

序号	k'<sub>i</sub>	d'<sub>i</sub>
			1	0.61831927677849863	1282.5138982275448
2	0.61831927677849863	1282.5138982275448
			3	0.53180308969301637	1410.5992309474163
4	0.34911833477857185	1679.9972604501922

步骤3.6：按照距离最小原则确定模型的合理解。

计算测量点(x₇,y₇)＝(1489.78，2200.34)分别到直线y＝k'₁×x+d'₁、y＝k'₂×x+d'₂、y＝k'₃×x+d'₃和y＝k'₄×x+d'₄的距离l'₁,l'₂,l'₃,l'₄：

m值加1，并取|l'₁|、|l'₂|、|l'₃|和|l'₄|中最小的|l'₄|所对应的(k'₄,d'₄)＝(0.349118，1679.997260)，作为所求直线的合理参数，记为(k^(m),d^(m))。

步骤3.7：计算所有数据点到新直线y-k^(m)x-d^(m)＝0的加权距离之和f^(m)(k^(m),d^(m))，m表示迭代次数。

(当前m＝2)

步骤3.8：判别是否为“最佳”解。

因f⁽¹⁾(k⁽¹⁾,d⁽¹⁾)＝1.997822，f⁽²⁾(k⁽²⁾,d⁽²⁾)＝0.221203,f⁽²⁾(k⁽²⁾,d⁽²⁾)＜f⁽¹⁾(k⁽¹⁾,d⁽¹⁾)，则重复步骤3.2至步骤3.8。经计算，m＝27、执行步骤3.8时，f⁽²⁷⁾(k⁽²⁷⁾,d⁽²⁷⁾)＝0.000000617345<10^-6，因此输出解(k⁽²⁶⁾,d⁽²⁶⁾)＝(0.361058，1662.438159)，并简记为(k,d)。

步骤3.9：使用“取点定向法”计算出t₇时刻的目标航向K₇。

取首末点坐标(x₁,y₁)和(x₇,y₇)，并计算y′₁和y′₇得：

x₁＝1469.22，y₁＝2193.15，y′₁＝k×x₁+d＝2192.912703；

x₇＝1489.78，y₇＝2200.34，y′₇＝k×x₇+d＝2200.335525。

按照图24所示流程计算得K₇＝70.147454度。

步骤4：使用“时间-路程加权匀速估计模型”迭代估计目标最佳线速度方程S-Vt-S₀＝0，得到t_n时刻的目标速度V_n并计算出目标在X、Y方向上的分量速度v_xn和v_yn。包括如下步骤：

步骤4.1：选取7个测量点{(x_i,y_i,T_i),i＝1,2,…7}作为“时间-路程加权匀速估计模型”的计算点。

因当前直线航迹段的实际测量点数为7，分别对应的时刻按照递增顺序表示为t₁,t₂,…,t₇，令T_i＝t_i，i＝1,2,…7，并取T_i时刻对应的测量点坐标(x_i,y_i)，构成{(x_i,y_i,T_i),i＝1,2,…7}作为“时间-路程加权匀速估计模型”的计算点，见表1.8。

步骤4.2：计算T_i时刻目标在航迹线上的垂点(cx_i,cy_i),i＝1,2,…n。

因{(x_i,y_i,T_i),i＝1,2,…7}为直线航迹段y-kx-d＝0上的测量点,令：A＝k＝0.361058，B＝-1，C＝d＝1662.438159，则：T_i时刻(x_i,y_i)在直线y-kx-d＝0上的垂点(cx_i,cy_i)为：

计算结果见表1.8。

表1.8：速度估计模型的计算点、垂点和路程

步骤4.3：计算从T₁时刻起，到T_i时刻，目标经过的路程S_i,i＝1,2,…7。

步骤4.3.1：求T_i时刻与上一计算点T_i-1时刻间的路程ΔS_i(i＝1时，ΔS_i＝0)。

因{(x_i,y_i,T_i),i＝1,2,…7}为直线航迹段上的测量点，计算结果见表1.8。

步骤4.3.2：计算从T₁时刻起，到T_i时刻，目标经过的路程计算结果见表1.8。

步骤4.4：使用“时间-路程加权匀速估计模型”迭代估计目标最佳线速度方程S-Vt-S₀＝0。

步骤4.4.1：对7个计算点的数据{(T_i,S_i),i＝1,2,…7}使用“时间-路程不加权匀速估计模型”粗略估计目标线速度方程S-V₀₁t-S₀₁＝0。对模型(25)的求解步骤如下。

①计算

②计算v,S₀，并记为(V₀₁,S₀₁)。

(单位：km/秒)，

(单位：km)。

步骤4.4.2：计算各点(T_i,S_i)到S-V₀₁t-S₀₁＝0的距离|l_i|之和。

其中：V⁽¹⁾＝V₀₁＝0.375466。，

步骤4.4.3：求各点(T_i,S_i)到S-V₀₁t-S₀₁＝0的距离l_i。

式中m表示迭代次数，n表示数据点数。m初始值为1，即：V⁽¹⁾＝V₀₁，计算结果见表1.9。

步骤4.4.4：求|l_i|的倒数。计算结果见表1.9。

步骤4.4.5：求各点的权值w_i。计算结果见表1.9。

表1.9：权值计算过程表

序号	l<sub>i</sub>	P<sub>i</sub>	w<sub>i</sub>
				1	0.857912	1.165621	0.062245
2	-1.251272	0.799187	0.042677
				3	-0.198914	5.027293	0.268461
4	0.292688	3.416602	0.182449
				5	-0.170360	5.869917	0.313458
6	1.114817	0.897008	0.047901
				7	-0.644872	1.550696	0.082808

步骤4.4.6：对7个计算点{(T_i,S_i),i＝1,2,…7}使用“时间-路程加权匀速估计模型”估计目标线速度方程S-Vt-S₀＝0。对模型(34)的求解步骤如下。

①计算w₀，wSt，wS，wt₁，wt₂(n＝7)。

②计算v,S₀，并记为m表示迭代次数。

步骤4.4.7：m值加1，计算所有数据点{(T_i,S_i),i＝1,2,…7}到新直线的加权距离之和m表示迭代次数。

(当前m为2)

步骤4.4.8：判别是否为“最佳”解。

因 则重复步骤4.4.3至步骤4.4.8。经计算，m＝3时，则输出解并简记为(V,S₀)。

步骤4.5：计算t₇时刻目标在X、Y方向上的分量速度v_x7和v_y7。

v_x7＝V₇cosk＝0.370678×cos(0.361058)＝0.346778，

v_y7＝V₇sink＝＝0.370678×sin(0.361058)＝0.130947。

式中V₇表示末点速度，即步骤4.4.8求得的V；k表示直线航迹线的斜率，在步骤3.8中已求出。

步骤5：计算7个观测点数据{(x_i,y_i,t_i),i＝1,2,…7}到直线y-kx-d＝0的垂足点的X轴坐标均值和时间均值并令最终求出t₇时刻目标在统一直角坐标系中的位置估值(P_x7,P_y7)。包括如下步骤：

步骤5.1：求表1.8中垂点{(T_i,cx_i,cy_i),i＝1,2,…7}的X轴坐标均值和时间均值

步骤5.2：求垂点的Y轴均值

其中：(k,d)表示直线航迹线的斜率和截距，在步骤3.8中已求出。

步骤5.3：求t₇＝911秒时，目标在统一直角坐标系中的位置坐标(P_x7,P_y7)。

至此我们就得到了目标由匀速圆周运动转变为匀速直线运动后，在t₇＝911秒时：

位置(P_x7,P_y7)＝(1490.189003，2200.660633)km，

速度V₇＝0.370678公里/秒＝1334.439229公里/小时，

航向K₇＝70.147454度。(正北方向为0度，正东方向为90度)

实施例2

本实施例具体描述本发明所提出的一种约束条件下的单雷达直线航迹线目标状态估计方法与传统卡尔曼滤波方法的对比试验过程，以及二者在执行效果上的比较。

所述对比试验过程包括如下步骤：

步骤1：产生模拟数据的真实值和测量值。

按照表2.1中的模拟目标参数值产生匀速圆周运动后转入匀速直线飞行状态的目标的60个位置点(前30个处于匀速圆弧航迹段，后30个处于匀速直线航迹段)的真实值和测量值，结果见表2.2。

表2.1：基本参数

表2.2：雷达观测数据的真实值和测量值

步骤2：将目标点数据{(ρ_i,θ_i,h_i,t_i)，i＝1,2,…60。}转换为统一直角坐标{(X_i,Y_j,t_i),i＝1,2,…60。}。

X_i＝X_zicosδ_xz-Y_zisinδ_xz+X_zx

Y_i＝X_zisinδ_xz+Y_zicosδ_xz+Y_zx

其中：所有测量点的高度h_i均简化为0，(X_zx,Y_zx)为雷达站址在中心统一直角坐标系中的坐标。δ_xz为雷达站址与直角坐标系中心点的经度差(单位为弧度)。坐标计算结果见表2.3。目标点真实位置和测量位置在中心统一直角坐标系中的显示如图2所示。

表2.3：目标点的中心统一直角坐标

步骤3：对表2.3中的直线航迹段测量值{(t_i,X_i,Y_j),i＝31,32,…60}使用本发明所述估计方法逐点进行目标位置(wP_xi,wP_yi)、速度wV_i和航向wK_i估计。其中：参数n取值为7。表2.4是n＝7时每次参与估值计算的点序号表。表2.5是应用本发明的估计结果。

表2.4：n＝7时每次参与估值计算的点

表2.5：应用本发明(n＝7)逐点进行目标状态估计

步骤4：对表2.3中的直线航迹段测量值{(t_i,X_i,Y_j),i＝31,32,…60}使用传统卡尔曼滤波方法逐点进行目标位置(kP_xi,kP_yi)、速度kV_i和航向kK_i估计。其中：Q、R参数取值见表2.1。表2.6是应用卡尔曼滤波方法的逐点估计结果。

表2.6：应用卡尔曼方法逐点进行目标状态估计

目标点真实位置、测量点位置、本发明估计位置和卡尔曼滤波估计位置在中心统一直角坐标系中的显示如图3、图4所示。图上可见，两种滤波方法得到的航迹线几乎重合，均优于测量航迹线，接下来我们对时间配准点、速度和航向等指标进一步统计、比较。

步骤5：时间配准点比较

为了进一步比较位置估值效果，我们引入时间配准点概念，即：相同时刻目标真实位置点与估计位置点称为一对时间配准点。t_i时刻一对时间配准点之间的距离称为时间配准点间距d_i，d₁～d_m的平均值定义为t_m时刻的平均时间配准点间距vd_m。根据表2.3～2.6的计算结果，可以统计使用本发明方法和使用卡尔曼方法得到时间配准点间距与平均时间配准点间距值，见表2.7和表2.8，其中(rP_x，rP_y)表示目标真实位置，wvd₆表示使用本发明方法得到的滤波器开始工作后1分钟以内估值点的平均时间配准点间距，wvd₃₀表示使用本发明方法得到的所有估值点的平均时间配准点间距；kvd₆表示使用卡尔曼方法得到的滤波器开始工作后1分钟以内估值点的平均时间配准点间距，kvd₃₀表示使用卡尔曼方法得到的所有估值点的平均时间配准点间距。

表2.7：本发明时间配准点间距wd与平均间距wvd

表2.8：卡尔曼时间配准点间距kd与平均间距kvd

图5和图6是本发明方法和卡尔曼方法时间配准点间距与平均时间配准点间距的对比图。经过统计，全程有21个点的时间配准点间距wd_i≤kd_i，占全部点数的70.00％；航迹起始1分钟内，总共得到6个测量点，其中有5个点的时间配准点间距wd_i≤kd_i，占比83.33％；所有时刻的平均时间配准点间距均有wvd_i≤kvd_i；全程平均时间配准点间距wvd₃₀占kvd₃₀的68.32％；航迹起始1分钟内平均时间配准点间距wvd₆占kvd₆的24.46％。本发明方法在测量初期、点数较少时能获得更加准确的目标位置估值，优势非常明显。于是可得到以下结论：在本次试验中，以目标真实位置为基准，使用本发明方法得到的目标位置点优于卡尔曼方法位置点。

步骤6：速度估值比较

我们定义t_i时刻目标估计速度与真实速度差的绝对值为速度差dv_i，dv₁～dv_m的平均值定义为t_m时刻的平均速度差vdv_m。根据表2.3～2.6的计算结果，可以统计使用本发明方法和使用卡尔曼方法得到的所有时间点上的速度差和平均速度差值，见表2.9，其中rV表示目标真实速度；wV_i、wdv_i、wvdv_i表示使用本发明方法得到的t_i时刻目标速度估值、速度差和平均速度差；kV_i、kdv_i、kvdv_i表示使用卡尔曼方法得到的t_i时刻目标速度估值、速度差和平均速度差。

表2.9：速度、速度差和平均速度差

图7、图8和图9是本发明方法和卡尔曼方法分别在速度、速度差和平均速度差的对比图。经过统计，全程有27个估值点的速度差wdv_i≤kdv_i，占全部点数的90.00％；航迹起始1分钟内得到的总共6个估值点的速度差均有wdv_i≤kdv_i，占比100.00％；所有时刻的平均速度差均有wvdv_i≤kvdv_i；全程平均速度差wvdv₃₀只占kvdv₃₀的25.54％；航迹起始1分钟内平均速度差wvdv₆只占kvdv₆的6.36％。图上可以看出，本发明得到的速度估值更加准确、稳定，而且本发明方法在测量初期、点数较少时就能获得更加准确的速度估值，优势非常明显。

步骤7：航向估值比较

我们定义t_i时刻目标估计航向与真实航向差的绝对值为航向差dk_i，dk₁～dk_m的平均值定义为t_m时刻的平均航向差vdk_m。根据表2.3～2.6的计算结果，可以统计使用本发明方法和使用卡尔曼方法得到的所有时间点上的航向差和平均航向差值，见表2.10，其中rK表示目标真实航向；wK_i、wdk_i、wvdk_i表示使用本发明方法得到的t_i时刻目标航向估值、航向差和平均航向差；kK_i、kdk_i、kvdk_i表示使用卡尔曼方法得到的t_i时刻目标航向估值、航向差和平均航向差。

表2.10：航向、航向差和平均航向差

图10、图11和图12是本发明方法和卡尔曼方法分别在航向、航向差和平均航向差的对比图。经过统计，全程有30个点的航向差wdk_i≤kdk_i，占全部点数的100.00％；航迹起始1分钟内得到的总共6个估值点的航向差均有wdk_i≤kdk_i，占比100.00％；所有时刻的平均航向差均有wvdk_i≤kvdk_i；全程平均航向差wvdk₃₀只占kvdk₃₀的40.60％；航迹起始1分钟内平均航向差wvdk₆占kvdk₆的49.40％。图上可以看出，本发明得到的航向估值更加准确、稳定，而且本发明方法在测量初期、点数较少时就能获得更加准确的航向估值，优势非常明显。

步骤8：出现野值后的估值效果比较

为了模拟野值出现，我们在表2.3的第36个(t＝899秒)测量点上人为加入较大误差。图13～23分别显示了出现野值后，本发明方法和卡尔曼方法在航迹线、时间配准点、速度和航向等指标上的比较效果，表2.11是测量中包含野值与不包含野值时各项指标的统计结果。

表2.11：包含/不包含野值时各项指标统计

从以上图、表可以看出，卡尔曼滤波方法在野值出现后，位置、速度、和航向估值出现较大突变。相比之下，本发明方法对野值抑制效果较为明显，无论是在野值点处，还是其后几个点上，估值效果都优于卡尔曼方法。

本发明基于复杂数据环境，在无须获取系统噪声和观测噪声等先验知识的情况下，通过对目标由匀速圆周运动转变为匀速直线运动初期时运动规律的研究，将目标在匀速圆周运动模型中的先验信息作为当前匀速直线运动模型的约束条件，在统一直角坐标系中对当前有限个测量点进行垂直距离加权估计得到与时间无关的目标运动方程和航向估值，通过目标时间-路程加权匀速估计模型得到线速度和分量速度估值，进而最终确定目标状态位置估计结果。通过与传统卡尔曼滤波方法在时间配准点间距、速度差和航向差等指标的对比试验表明，本发明不但易于工程实现，而且比传统卡尔曼滤波方法具有更优的起始和全程目标状态估计效果，以及野值抑制能力。该方法有效解决了传统卡尔曼滤波方法应用条件比较苛刻、过程噪声协方差和量测噪声协方差难以确定和系统测量中出现野值，可能导致滤波发散或目标状态估计突变等问题，极大地改善了单雷达数据质量，为在本系统中开展目标跟踪与数据融合工作提供了一种稳定、有效且易于工程实现的滤波方法。本发明估计方法科学、方案实施步骤合理，目标估计效果理想，对于提高了多雷达目标状态估计的一致性和准确性具有重要意义。本发明所提供的方法的时间复杂度和空间复杂度都很低，可操作性和实用性很强。

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

(1)为了节省计算时间，在步骤1中，可以存储前n－1个测量点的统一直角坐标，每次只对第n个测量点进行极坐标到统一直角坐标转换。

(2)关于n的取值问题，我们认为一般取4～15为宜。但是，在相同数据环境中，在保证所取的n个测量点均处于目标匀速直线飞行状态的前提下，n越大，目标状态估值精度就越高。

(3)在步骤4.3.1中，求T_i时刻与上一计算点T_i-1时刻间的路程ΔS_i时，当i＝1时，如果ΔS_i设置为与上一估值点(如：圆弧航迹段末点)的距离，效果就更为理想。

(4)在步骤2.1和步骤3.5中，也可以采用费拉里法等其它方法对一元四次方程(10)、(20)进行求解。

Claims (8)

1.一种约束条件下的单雷达直线航迹线目标状态估计方法，其特征在于，所述估计方法应用于单雷达数据跟踪与多雷达数据融合系统的前期数据预处理过程中；所述估计方法包括如下步骤：

步骤1：将某雷达针对某个由匀速圆周运动转变为匀速直线运动后直线段上目标的n个观测点数据(ρ_i,θ_i,h_i,t_i)转换为统一直角坐标(x_i,y_i,t_i)；其中，(ρ_i,θ_i,h_i,t_i)表示t_i时刻该雷达测得的目标距离ρ_i、方位θ_i和高度h_i，i＝1,2,…n；

步骤2：对n个观测点数据{(x_i,y_i),i＝1,2,…n}使用约束条件下的不加权直线参数估计模型粗略估计目标运动状态方程y-k₁x-d₁＝0；

步骤3：将方程y-k₁x-d₁＝0作为初始状态，对n个观测点数据{(x_i,y_i),i＝1,2,…n}使用约束条件下的加权直线参数估计模型进行迭代估计，得到最佳目标运动状态方程y-kx-d＝0，并使用取点定向法计算出t_n时刻的目标航向K_n；

步骤4：使用时间-路程加权匀速估计模型迭代估计目标最佳线速度方程S-Vt-S₀＝0，得到t_n时刻的目标速度V_n并计算出目标在X、Y方向上的分量速度v_xn和v_yn：

步骤5：计算n个观测点数据{(x_i,y_i,t_i),i＝1,2,…n}到直线y-kx-d＝0的垂足点的X轴坐标均值和时间均值并令则t_n时刻目标在统一直角坐标系中的位置估值(P_xn,P_yn)为：

由此得到目标由匀速圆周运动转变为匀速直线运动后，在t_n时刻的位置(P_xn,P_yn)、速度V_n和航向K_n。

2.如权利要求1所述的约束条件下的单雷达直线航迹线目标状态估计方法，其特征在于，所述步骤1包括如下步骤：

步骤1.1：选择某雷达观测数据集{(ρ_i,θ_i,h_i,t_i),i＝1,2,…n}，是针对某个由匀速圆周运动转变为匀速直线运动后直线段上目标的时序观测数据；

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

步骤1.3：将(X_zi,Y_zi)，i＝1,2,…n转换为中心统一直角坐标(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为雷达站址与直角坐标系中心点的经度差，单位为弧度。

3.如权利要求2所述的约束条件下的单雷达直线航迹线目标状态估计方法，所述步骤1.1中，t_i＜t_i+1，且n不大于N，N表示参与目标状态估计的最近时间内有限个观测数据点的个数，N取4～15。

4.如权利要求2所述的约束条件下的单雷达直线航迹线目标状态估计方法，其特征在于，所述步骤2包括如下步骤：

步骤2.1：用所有观测点{(x_i,y_i),i＝1,2,…n}到待求直线的距离l_i的平方和最小作为估计模型，用待求直线航迹线必须与已知圆相切作为约束条件，计算在此条件下的这条直线的最佳参数(k,d)；其中，k为直线的斜率，d为该直线在x轴上的截距；

其中的已知圆就是目标在进入匀速直线运动前所处的圆周运动轨迹，为已知变量，其圆心为半径为r；

所述约束条件下的不加权直线参数估计模型如下：

由(2)式解得：

(3)式代入(1)式得：

对于(4)式，应有f(k)对于k的导数等于零，即：

记

则(5)式化简为：

为了消去“根号”，将(6)改写成：

(7)式两边平方后得：

(a₃k²+(a₂-b₂)k-a₃)²＝(a₁r+b₁rk)²(1+k²) (8)

(8)式两边展开并整理后得：

记：

则(9)式可化简为：

k⁴+bk³+ck²+fk+e＝0 (10)

在此采用求实系数代数方程全部根的牛顿-下山法对该式进行求解，得到k₁、k₂、k₃和k₄四个根；将k₁、k₂、k₃、k₄带入(3)式，对应得到d₁、d₂、d₃、d₄；

步骤2.2：按照距离最小原则确定模型的合理解；

(k₁,d₁)、(k₂,d₂)、(k₃,d₃)和(k₄,d₄)都是约束条件下的不加权直线参数估计模型的实根；按照点{(x_i,y_i)，i＝1,2,…n}到所求直线的距离的平方和最小原则，确定合理的直线参数值；该问题也可以简化为：计算点(x_n,y_n)分别到直线y＝k₁×x+d₁、y＝k₂×x+d₂、y＝k₃×x+d₃和y＝k₄×x+d₄的距离l₁,l₂,l₃,l₄：

取|l₁|、|l₂|、|l₃|和|l₄|中最小的|l_i|所对应的(k_i,d_i)，作为所求直线的合理参数，记为(k₁,d₁)。

5.如权利要求4所述的约束条件下的单雷达直线航迹线目标状态估计方法，其特征在于，所述步骤3包括如下步骤：

步骤3.1：计算各点(x_i,y_i)到直线y-k₁x-d₁＝0的距离|l_i|之和；

n为数据点数；

步骤3.2：求各点(x_i,y_i)到直线y-k₁x-d₁＝0的距离l_i；

式中m表示迭代次数，n表示数据点数；m初始值为1，即：

k⁽¹⁾＝k₁，d⁽¹⁾＝d₁；

步骤3.3：求|l_i|的倒数；

步骤3.4：求各点的权值w_i；

步骤3.5：用所有观测点{(x_i,y_i),i＝1,2,…n}到待求直线的加权距离w_il_i的平方和最小作为估计模型，用待求直线航迹线必须与已知圆相切作为约束条件，计算在此条件下的这条直线的最佳参数(k,d)；

其中的已知圆就是目标在进入匀速直线运动前所处的圆周运动轨迹，为已知变量，其圆心为半径为r；

所述约束条件下的加权直线参数估计模型如下：

由(12)式得：

(13)式代入(11)得：

对于(14)式，应有f(k)对于k的导数等于零，即：

记： 则(15)式化简为：

为了消去“根号”，将(16)改写成：

(17)式两边平方后得：

(a'₃k²+(a'₂-b'₂)k-a'₃)²＝(a'₁r+b'₁rk)²(1+k²) (18)

(18)式两边展开并整理后得：

记：

则(19)式可化简为：

k⁴+b'k³+c'k²+f'k+e'＝0 (20)

在此采用求实系数代数方程全部根的牛顿-下山法对该式进行求解，得到k'₁、k'₂、k'₃和k'₄四个根；将k'₁、k'₂、k'₃、k'₄带入(13)式，对应得到d'₁、d'₂、d'₃、d'₄；

步骤3.6：按照距离最小原则确定方程的合理解；

(k'₁,d'₁)、(k'₂,d'₂)、(k'₃,d'₃)和(k'₄,d'₄)都是模型的实根；我们按照点{(x_i,y_i)，i＝1,2,…n}到所求直线的距离的平方和最小原则，确定合理的直线参数值；该问题也可以简化为：计算点(x_n,y_n)分别到直线y＝k'₁×x+d'₁、y＝k'₂×x+d'₂、y＝k'₃×x+d'₃和y＝k'₄×x+d'₄的距离l'₁,l'₂,l'₃,l'₄：

m加1，并取|l'₁|、|l'₂|、|l'₃|和|l'₄|中最小的|l'_i|所对应的(k'_i,d'_i)，作为所求直线的合理参数，记为(k^(m),d^(m))，m表示迭代次数；

步骤3.7：计算所有数据点到新直线y-k^(m)x-d^(m)＝0的加权距离之和f^(m)(k^(m),d^(m))；

式中m表示迭代次数，n表示数据点数；

步骤3.8：判别是否为“最佳”解；

若f^(m)(k^(m),d^(m))≥f^(m-1)(k^(m-1),d^(m-1))或f^(m)(k^(m),d^(m))≤10^-6，则输出解(k^(m-1),d^(m-1))，并简记为(k,d)；否则重复步骤3.2至步骤3.8。式中m表示迭代次数；

步骤3.9：使用取点定向法计算出t_n时刻的目标航向K_n。

6.如权利要求5所述的约束条件下的单雷达直线航迹线目标状态估计方法，所述步骤3.9中取点定向法的实现过程为：设(k,d)是估计得到的雷达观测航迹线参数，(x₁,y₁)和(x_n,y_n)是雷达对目标的首末测量点，经过滤波后，这两个点的坐标变为(x₁,y₁′)和(x_n,y_n′)，其中：y₁′＝k×x₁+d，y_n′＝k×x_n+d；令：Δx＝x_n-1-x₁，Δy＝y′_n-1-y₁′，π为圆周率，接下来依次进行如下判断和计算；

①如果Δy等于0，转②，否则转③；

②如果Δx大于0，航向K取值为0度，否则，航向K取值为180度，程序结束；

③如果Δx等于0，转④，否则转⑤；

④如果Δy大于等于0，航向K取值为90度，否则，航向K取值为270度，取点定向结束；

⑤如果Δy大于0，航向K取值为度，否则，航向K取值为度，取点定向结束。

7.如权利要求5所述的约束条件下的单雷达直线航迹线目标状态估计方法，其特征在于，所述步骤4包括如下步骤：

步骤4.1：选取n个测量点{(x_i,y_i,T_i),i＝1,2,…n}作为时间-路程加权匀速估计模型的计算点；

(1)当直线航迹段的实际测量点数m′小于步骤1中确定的目标观测点数n时，需要使用目标处于圆弧运动期的n-m′个测量点，此时这n个点分别对应的时刻按照递增顺序表示为：t_m'-n+1,t_m'-n+2,…,t_-1,t₀,t₁,t₂,…,t_m'，其中：t_m'-n+1,…,t_-1,t₀为圆弧航迹段上的测量点时刻，t₁,t₂,…,t_m'为直线航迹段上的测量点时刻。

(2)当直线航迹段的实际测量点数m′不小于步骤1中确定的目标观测点数n时，这n个点均为直线航迹段上的测量点，分别对应的时刻按照递增顺序表示为t_m'-n+1,t_m'-n+2,…,t_m'。

令：

T₁＝t_m'-n+1，

T₂＝t_m'-n+2，

……

T_n＝t_m'；

并取T_i时刻对应的测量点坐标(x_i,y_i)，构成{(x_i,y_i,T_i),i＝1,2,…n}作为时间-路程加权匀速估计模型的计算点；

步骤4.2：计算T_i时刻目标在航迹线上的垂点(cx_i,cy_i),i＝1,2,…n；

(1)当(x_i,y_i,T_i)为直线航迹段(y-kx-d＝0)上的测量点时；

令：A＝k，B＝-1，C＝d，则：T_i时刻(x_i,y_i)在直线y-kx-d＝0上的垂点(cx_i,cy_i)为：

(2)当(x_i,y_i,T_i)为圆心为半径为r的圆弧航迹段上的测量点时；

圆弧航迹线标准方程为：

式中圆心和半径r已知；

过观测点(x_i,y_i)与圆心的直线方程为：

解(22)式得：

将(23)式代入(21)式得：

解(24)式得：

代入(23)式得：

取(x₁,y₁)和(x₂,y₂)中距离测量点(x_i,y_i)最近的点作为合理解，记为(cx_i,cy_i)；

步骤4.3：计算从T₁时刻起，到T_i时刻，目标经过的路程S_i,i＝1,2,…n；

步骤4.3.1：求T_i时刻与上一计算点T_i-1时刻间的路程ΔS_i，i＝1时，ΔS_i＝0；

(1)当(x_i,y_i,T_i)为直线航迹段上的测量点时；

(2)当(x_i,y_i,T_i)为圆心为半径为r的圆弧航迹段上的测量点时；

求两相邻计算点之间的弦长：

求圆心到弦的距离：

求圆心角：

求圆心角θ_i的对应的弧长即为ΔS_i：

ΔS_i＝rθ_i

步骤4.3.2：计算从T₁时刻起，到T_i时刻，目标经过的路程S_i；

步骤4.4：使用时间-路程加权匀速估计模型迭代估计目标最佳线速度方程S-Vt-S₀＝0；

该步骤4.4包括如下步骤：

步骤4.4.1：对n个计算点的数据{(T_i,S_i),i＝1,2,…n}使用时间-路程不加权匀速估计模型粗略估计目标线速度方程S-V₀₁t-S₀₁＝0；

时间-路程不加权匀速估计模型为：

式中S_i是通过步骤4.3.2计算出来的对应T_i时刻目标经过的路程，S₀为待求的调节路程，v为待求的匀速度；

接下来对模型(25)进行求解；

f(S₀,v)分别对S₀和v求偏导数并令其等于零，以t_i表示T_i，则有：

(26)、(27)式化简得：

由(28)式解得：

(30)式代入(29)式得：

解(31)式得：

(32)式代入(30)式得：

由以上推导可知，不加权匀速度估计模型求解过程如下：

令：则有：

并将(v,S₀)记为(V₀₁,S₀₁)。

步骤4.4.2：计算各点(T_i,S_i)到S-V₀₁t-S₀₁＝0的距离|l_i|之和；

其中，n为数据点数，V⁽¹⁾＝V₀₁，

步骤4.4.3：求各点(T_i,S_i)到S-V₀₁t-S₀₁＝0的距离l_i；

式中m表示迭代次数，n表示数据点数；m初始值为1，即：V⁽¹⁾＝V₀₁，

步骤4.4.4：求|l_i|的倒数；

步骤4.4.5：求各点的权值w_i；

步骤4.4.6：对n个计算点{(T_i,S_i),i＝1,2,…n}使用时间-路程加权匀速估计模型估计目标线速度方程S-Vt-S₀＝0；

时间-路程加权匀速估计模型为：

式中S_i是通过步骤4.3.2计算出来的对应T_i时刻目标经过的路程，w_i是通过步骤4.4.5计算出来的权值，S₀为待求的调节路程，v为待求的匀速度。接下来对模型(34)进行求解。

f(S₀,v)分别对S₀和v求偏导数并令其等于零，则有：

(35)、(36)式化简得：

由(37)式解得：

(39)式代入(38)式得：

解(40)式得：

令：

则：

将(v,S₀)记为m表示迭代次数。

步骤4.4.7：m加1，计算所有数据点{(T_i,S_i),i＝1,2,…n}到新直线的加权距离之和

式中m表示迭代次数，n表示数据点数。

步骤4.4.8：判别是否为“最佳”解；

若或则输出解并简记为(V,S₀)；否则重复步骤4.4.3至步骤4.4.8；

步骤4.5：计算t_n时刻目标在X、Y方向上的分量速度v_xn和v_yn；

式中n表示数据点数，t_n表示末点时刻；V_n表示末点速度，即步骤4.4.8求得的V；k表示直线航迹线的斜率，在步骤3.8中已求出。

8.如权利要求7所述的约束条件下的单雷达直线航迹线目标状态估计方法，其特征在于，所述步骤5包括如下步骤：

步骤5.1：求步骤4.2得到的垂点{(T_i,cx_i,cy_i),i＝1,2,…n}的X轴坐标均值和时间均值

步骤5.2：求垂点的Y轴均值

其中：(k,d)表示直线航迹线的斜率和截距，在步骤3.8中已求出；

步骤5.3：求t_n时刻目标在统一直角坐标系中的位置坐标(P_xn,P_yn)；

至此得到了目标由匀速圆周运动转变为匀速直线运动后，在t_n时刻的位置(P_xn,P_yn)、速度V_n和航向K_n估值。