一种针对单雷达直线航迹线的目标状态估计方法与流程

本发明属于传感器目标跟踪与数据融合
技术领域：
，具体涉及一种在复杂数据环境中，且系统噪声和观测噪声均未知的情况下，针对单雷达直线航迹线的目标状态估计方法。
背景技术：
：目标跟踪是指利用传感器(雷达、声纳、红外等)所获得的对运动目标(飞机、坦克、舰艇等)的时序离散观测数据，对目标的运动状态(位置、速度、航向、加速度等)进行估计和跟踪的方法。卡尔曼滤波是目标跟踪与数据融合应用中一种经典的最优线性无偏估计方法，当目标状态空间和测量空间均为线性高斯系统时，可以在干扰噪声背景中对目标运动状态进行最小方差估计和跟踪，而且它是一种递推估计方法，根据每一时刻的观测数据进行滤波更新，大大减少了计算量和存储量，容易满足实时计算要求，因而得到了广泛的实际应用并成为目标跟踪领域的理论基础。在多雷达组网探测与目标跟踪系统(以下简称“系统”)中，虽然传感器类型相对单一，但是雷达跨代使用、部署分散且地域跨度大，微波、短波、光缆、电缆等多种数据传输手段并存，探测环境易受电磁、气象、季节等因素影响，因此，传统卡尔曼滤波方法在本系统应用中存在以下难点：1、传统卡尔曼滤波应用条件比较苛刻，与实际系统中复杂的应用数据环境之间存在矛盾。实际系统中测量数据采样间隔较大、测量误差、坐标转换误差以及数据标定与传输误差同时存在。而卡尔曼滤波要求系统模型精确以及系统误差模型和观测误差模型已知，这在实际应用中是很难满足的，或者在系统工作过程中，任一误差分量模型发生变化，都将导致传统卡尔曼滤波方法发散或精度降低。2、过程噪声协方差q和量测噪声协方差r在参数值选取上存在困难，且不能兼顾稳态性和瞬时性。传统卡尔曼滤波算法中需要预先估计过程噪声协方差q和量测噪声协方差r，且在整个滤波过程中保持不变，所以当这两个估计噪声协方差矩阵参数不准确时，在滤波过程中将产生累计误差，导致滤波发散。同时，当q取值较小，r取值较大时，滤波后形成的目标状态估值稳态性较好，但瞬时性变差，表现在目标位置估计上就是时间配准点(同一时刻目标真实位置点与估计位置点)间距较大，有时可能落后一个测量周期。反之，当q取值较大，r取值较小时，滤波后形成的目标状态估值瞬时性较好，但稳态性较差，表现为：目标位置估计上更加靠近测量点位置，但是目标航向和速度估值变化较大，这样滤波也就失去了意义。3、系统测量中不可避免的野值出现，对传统卡尔曼滤波器的工作产生严重影响。由于雷达故障、人为因素或外部观测条件发生改变，就有可能使观测数据发生突变，即出现野值(也称为飞点)。野值可能导致滤波发散或目标状态估计在短时间内失效，表现为估值航迹线出现较大锯齿、目标航向和速度估计突变。技术实现要素：(一)要解决的技术问题本发明要解决的技术问题是：如何提供一种在复杂数据环境中，且系统噪声和观测噪声均未知的情况下，针对单雷达直线航迹线的目标状态估计方法。(二)技术方案为解决上述技术问题，本发明提供一种针对单雷达直线航迹线的目标状态估计方法，所述估计方法应用于单雷达数据跟踪与多雷达数据融合系统的前期数据预处理过程中；所述修正方法包括如下步骤：步骤1：将某雷达针对某个处于匀速直线飞行期的目标的n个观测点数据(ρi,θi,hi,ti)转换为统一直角坐标(xi,yj,ti)；其中，(ρi,θi,hi,ti)表示ti时刻该雷达测得的目标距离ρi、方位θi和高度hi，i＝1,2,…n；步骤2：对目标x轴运动分量{(ti,xi),i＝1,2,…n}使用加权直线参数估计模型进行迭代估计，得到tn时刻目标在x方向上的位置pxn、速度vxn和方向kxn；步骤3：对目标y轴运动分量{(ti,yj),i＝1,2,…n}使用加权直线参数估计模型进行迭代估计，得到tn时刻目标在y方向上的位置pyn、速度vyn和方向kyn；步骤4：则tn时刻目标在统一直角坐标系中的位置为(pxn,pyn)，此时，目标速度vn和航向kn分别为：(三)有益效果针对上述现有技术中的三个问题，本发明提供一种基于复杂数据环境，且无须获取系统噪声和观测噪声等先验知识的情况下，针对单雷达直线航迹线观测数据的目标状态估计方法，通过与传统卡尔曼滤波算法在目标位置、航向和速度估值上进行对比实验表明，本发明取得了更好的估值效果，为在本系统中开展目标跟踪与数据融合工作提供了一种稳定、有效且易于工程实现的滤波方法。与现有技术相比较，本发明提出的适应多雷达组网探测与目标跟踪系统复杂数据环境，针对单雷达直线航迹线的目标状态估计方法，充分利用目标处于匀速直线运动状态这一特点，分别对当前有限测量点的x、y分量进行相对于的测量时间的垂直距离加权估计，然后根据得到的分量直线航迹线模型参数确定最终的目标状态估计(位置、速度、航向)结果。对比试验表明，该方法不但易于工程实现，而且比传统卡尔曼滤波方法具有更优的目标状态估计效果。附图说明图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中含有野值时使用卡尔曼方法和使用本发明方法得到的平均航向差随测量点数变化情况显示图。具体实施方式为使本发明的目的、内容、和优点更加清楚，下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。为解决现有技术的问题，本发明提供一种针对单雷达直线航迹线的目标状态估计方法，如图1所示，所述估计方法应用于单雷达数据跟踪与多雷达数据融合系统的前期数据预处理过程中；所述修正方法包括如下步骤：步骤1：将某雷达针对某个处于匀速直线飞行期的目标的n个观测点数据(ρi,θi,hi,ti)转换为统一直角坐标(xi,yj,ti)；其中，(ρi,θi,hi,ti)表示ti时刻该雷达测得的目标距离ρi、方位θi和高度hi，i＝1,2,…n；步骤2：对目标x轴运动分量{(ti,xi),i＝1,2,…n}使用加权直线参数估计模型进行迭代估计，得到tn时刻目标在x方向上的位置pxn、速度vxn和方向kxn；步骤3：对目标y轴运动分量{(ti,yj),i＝1,2,…n}使用加权直线参数估计模型进行迭代估计，得到tn时刻目标在y方向上的位置pyn、速度vyn和方向kyn；步骤4：则tn时刻目标在统一直角坐标系中的位置为(pxn,pyn)，此时，目标速度vn和航向kn分别为：其中，所述步骤1包括如下步骤：步骤1.1：选择某雷达观测数据集{(ρi,θi,hi,ti),i＝1,2,…n}，是针对某个处于匀速直线飞行期的目标的时序观测数据；如果是三坐标雷达观测数据，一般认为hi保持不变；如果是两坐标雷达观测数据，则认为hi＝0。θi取值范围[0,360)，单位为度，其中雷达正北方向为0度，正东方向为90度。步骤1.2：将雷达观测数据(ρi,θi,hi,ti)，i＝1,2,…n转换为以本站为中心的二维直角坐标(xzi,yzi)：步骤1.3：将(xzi,yzi)，i＝1,2,…n转换为中心统一直角坐标(xi,yi)：xi＝xzicosδxz-yzisinδxz+xzxyi＝xzisinδxz+yzicosδxz+yzx其中：(xzx,yzx)为雷达站址在中心统一直角坐标系中的坐标；δxz为雷达站址与直角坐标系中心点的经度差，单位为弧度。其中，所述步骤1.1中，ti<ti+1，且n不大于n，n表示参与目标状态估计的最近时间内有限个观测数据点的个数，一般n取4～15为宜。其中，所述步骤2包括如下步骤：步骤2.1：对目标x轴运动分量{(ti,xi),i＝1,2,…n}采用不加权直线参数估计模型粗略估计在x方向上的目标运动状态方程x-k1t-d1＝0；不加权直线参数估计模型的运算方法包括如下步骤：步骤2.1.1：用目标x轴运动分量{(ti,xi),i＝1,2,…n}(简记为：{(xi,yi)，i＝1,2,…n})的所有点到某条直线的距离li的平方和最小作为条件构造直线，计算在此条件下的这条直线的最佳参数(k,d)；其中，k为该直线的斜率，d为该直线在x轴上的截距；对于(1)式，应有f(k,d)分别对k和d求偏导数，并等于零，即有下式成立：于是有：化简得：-kd2+(-a1k2+2b1k+a1)d+c0k2+(a2-b2)k-c0＝0(2)其中：化简得：b1-a1k-d＝0(3)其中：由(3)式解得：d＝b1-a1k(4)(4)式代入(2)式得：整理后得：记a＝c0-a1b1，c＝a1b1-c0，则(5)式化简为：ak2+bk+c＝0(6)解(6)式得：(7)，(8)分别代入(4)式得：d1＝b1-a1k1d2＝b1-a1k2步骤2.1.2：按照距离最小原则确定方程的合理解；(k1,d1)和(k2,d2)都是方程(2)的实根，且k1×k2＝-1，即解得的两条直线相互垂直；按照点{(xi,yi)，i＝1,2,…n}到所求直线的距离的平方和最小原则，确定合理的直线参数值；该问题也可以简化为：计算点(x1,y1)分别到直线y＝k1×x+d1和直线y＝k2×x+d2的距离l1,l2：若|l1|<|l2|，则取(k1,d1)，否则取(k2,d2)作为所求直线的合理参数，记为(k1,d1)；步骤2.2：采用加权直线参数估计模型迭代估计目标在x方向上的运动状态方程x-kxt-dx＝0；加权直线参数估计模型的运算方法包括如下步骤：步骤2.2.1：计算各点(xi,yi)到直线y-k1x-d1＝0的距离|li|之和；n为数据点数；步骤2.2.2：求各点(xi,yi)到直线y-k1x-d1＝0的距离li；式中m表示迭代次数，n表示数据点数；m初始值为1，即：k(1)＝k1，d(1)＝d1；步骤2.2.3：求|li|的倒数；步骤2.2.4：求各点的权值vi；步骤2.2.5：求解加权直线参数估计模型；用所有数据点{(xi,yi)，i＝1,2,…n}到某条直线的加权距离(vi×li)的平方和最小作为条件构造直线，计算在此条件下的这条直线的最佳参数(k,d)；对于(9)式，应有f(k,d)分别对k和d求偏导数，并等于零，即有下式成立：于是有：记：则(10)式化简为：c0'-b2'k+b1'kd+a2'k-c0'k2-a1'k2d+a1'd+b1'kd-kd2＝0(11)由(12)式解得：d＝b′1-a′1k(13)将(13)式代入(11)式得：(-c0'-a1'b1')k2+(b′12+a2'-a′12-b2')k+(c0'+a1'b1')＝0(14)记a'＝-c0'-a1'b1'，b'＝b′12+a′2-a′12-b'2，c'＝c0'+a1'b1'，则(14)式化简为：a'k+b'k+c'＝0(15)解(15)式得：将(16)，(17)分别代入(13)式得：d1＝b1'-a1'k1d2＝b1'-a1'k2步骤2.2.6：按照距离最小原则确定方程的合理解；按照点{(xi,yi)，i＝1,2,…n}到所求直线的距离的平方和最小原则，确定合理的直线参数值；该问题也可以简化为：计算点(x1,y1)分别到直线y＝k1×x+d1和直线y＝k2×x+d2的距离l1,l2：m加1，若|l1|+|l2|>lmin(lmin初值为106)，则输出(k(m-1),d(m-1))作为所求直线的合理参数，迭代过程退出；否则lmin＝|l1|+|l2|；若|l1|<|l2|，则取(k1,d1)，否则取(k2,d2)作为所求直线的合理参数，记为记为(k(m),d(m))，m表示迭代次数；步骤2.2.7：计算所有数据点到新直线y-k(m)x-d(m)＝0的加权距离之和f(m)(k(m),d(m))；式中m表示迭代次数，n表示数据点数；步骤2.2.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))，并简记为(kx,dx)；否则重复步骤2.2.2至步骤2.2.8；式中m表示迭代次数；步骤2.3：计算tn时刻目标在x方向上的位置pxn、速度vxn和方向kxn：pxn＝kxtn+dx，vxn＝kx，kxn＝kx。其中，所述步骤3包括如下步骤：步骤3.1：对目标y轴运动分量{(ti,yi),i＝1,2,…n}采用步骤2.1的不加权直线参数估计模型粗略估计在y方向上的目标运动状态方程y-k1t-d1＝0；步骤3.2：采用步骤2.2的加权直线参数估计模型迭代估计目标在y方向上的运动状态方程y-kyt-dy＝0；步骤3.3：计算tn时刻目标在y方向上的位置pyn、速度vyn和方向kyn：pyn＝kytn+dy，vyn＝ky，kyn＝ky。实施例1本实施例具体描述本发明所提出的一种针对单雷达直线航迹线的目标状态估计方法，所述估计方法应用于单雷达数据跟踪与多雷达数据融合系统的前期数据预处理过程中。所述估计方法包括如下步骤：步骤1：将某雷达针对某个处于匀速直线飞行期的目标的n个观测点数据(ρi,θi,hi,ti)(表示ti时刻该雷达测得的目标距离ρi、方位θi和高度hi，i＝1,2,…n。)转换为统一直角坐标(xi,yj,ti)。如果n＝1，则输出目标位置(x1,y1,t1)、目标速度和航向均为0，程序退出。步骤1.1：选择某雷达针对某个处于匀速直线飞行期的目标的n个时序观测数据{(ρi,θi,hi,ti),i＝1,2,…n}，其中：ti<ti+1，且n取值为7。因为是两坐标雷达观测数据，认为hi＝0。相关数据见表1.1～表1.2。表1.1：基本参数表1.2：雷达观测数据序号时间t(秒)距离ρ(km)方位θ(度)高度h(km)130213.67205.970240209.76205.970351205.81205.680461201.68205.10570199.27205.960681194.32205.560790191.82204.80步骤1.2：将雷达观测数据(ρi,θi,hi,ti)，i＝1,2,…29转换为以本站为中心的二维直角坐标(xzi,yzi)：计算结果见表1.3。表1.3：以雷达为中心的二维直角坐标序号时间t(秒)xz(km)yz(km)130-93.56-192.10240-91.85-188.59351-89.17-185.49461-85.54-182.64570-87.24-179.16681-83.84-175.30790-80.47-174.13步骤1.3：将(xzi,yzi)，i＝1,2,…29转换为中心统一直角坐标(xi,yi)：xi＝xzicosδxz-yzisinδxz+xzxyi＝xzisinδxz+yzicosδxz+yzx其中：(xzx,yzx)为雷达站址在中心统一直角坐标系中的坐标。δxz为雷达站址与直角坐标系中心点的经度差(单位为弧度)。计算结果见表1.4。表1.4：测量点的中心统一直角坐标序号时间t(秒)x(km)y(km)1301503.102334.162401504.882337.713511507.612340.854611511.292343.765701509.652347.226811513.122351.137901516.512352.37步骤2：对目标x轴运动分量{(ti,xi),i＝1,2,…7}使用“加权直线参数估计模型”进行迭代估计，得到tn时刻目标在x方向上的位置pxn、速度vxn和方向kxn。步骤2.1：对目标x轴运动分量{(ti,xi),i＝1,2,…7}采用“不加权直线参数估计模型”粗略估计在x方向上的目标运动状态方程x-k1t-d1＝0；所述不加权直线参数估计模型的运算过程包括如下步骤：步骤2.1.1：用所有点{(ti,xi),i＝1,2,…7}(简记为：{(xi,yi)，i＝1,2,…7})到某条直线的距离li的平方和最小作为条件构造直线，计算在此条件下的这条直线的最佳参数(k,d)。求得直线参数。对(1)式的求解步骤如下。表1.5：目标在x轴的运动分量序号t＝xx＝y1301503.102401504.883511507.614611511.295701509.656811513.127901516.51①计算a1,a2,b1,b2,c0(n＝7)。表1.6：参数计算结果1a160.428571428571431a24054.7142857142858b11509.4479365716902b22278451.8354599578c091298.2384983306②计算a,b,c。a＝c0-a1b1c＝a1b1-c0表1.7：中间参数计算结果2a84.4560454984603b384.33980139158666c-84.4560454984603③解方程，计算所有解。d1＝b1-a1k1d2＝b1-a1k2表1.8：x方向观测航迹线参数解粗略计算结果k10.21004805826508055k2-4.7608152546594873d11496.755032479386d21797.1372012461136步骤2.1.2：按照距离最小原则确定方程的合理解。计算测量点(x1,y1)分别到直线y＝k1×x+d1和直线y＝k2×x+d2的距离l1,l2：因|l1|<|l2|，则取(k1,d1)＝(0.210048，1496.755032)作为所求直线的合理参数，记为(k1,d1)。步骤2.2：采用“加权直线参数估计模型”迭代估计目标在x方向上的运动状态方程x-kxt-dx＝0；所述加权直线参数估计模型的运算过程包括如下步骤：步骤2.2.1：计算各点(xi,yi)到直线y-k1x-d1＝0的距离|li|之和(n＝7)。步骤2.2.2：求各点(xi,yi)到直线y-k1x-d1＝0的距离li。式中：m表示迭代次数，初始值为1，k(1)＝k1，d(1)＝d1。计算结果见表1.9。步骤2.2.3：求|li|的倒数。计算结果见表1.9。步骤2.2.4：求各点的权值vi。计算结果见表1.9表1.9:距离和权值计算结果序号lipivi10.04098624.398420.6193492-0.275443.6305620.09216130.1345857.4302210.18861541.6803210.5951240.0151075-1.773070.5639950.0143176-0.636981.5699050.03985270.8295951.2054080.030599步骤2.2.5：求解加权直线参数估计模型。用所有数据点{(xi,yi)，i＝1,2,…7}到某条直线的加权距离(vi×li)的平方和最小作为条件构造直线，计算在此条件下的这条直线的最佳参数(k,d)。对(9)式的求解步骤如下。①计算a0,a1,a2,b1,b2,c0(n＝7)。表1.10：中间参数计算结果3a0’0.43062036934622938a1’32.286170046601832a2’1094.2668062910109b1'1503.5792870236251b2'2260752.9983791234c0-48555.779283781565②计算a',b',c'。a'＝-c0'-a1'b1'，b'＝b′12+a'2-a′12-b'2，c'＝c0'+a1'b1'；表1.11：中间参数计算结果4a’10.962744388460123b'49.54401736240834c’-10.962744388460123③解方程，计算所有解。d1＝b1'-a1'k1d2＝b1'-a1'k2表1.12:x方向观测航迹线参数解第m次计算结果k10.21138550191927813k2-4.7306934057467718d11496.7544587632729d21656.3152587599034步骤2.2.6：按照距离最小原则确定方程的合理解。计算测量点(x1,y1)分别到直线y＝k1×x+d1和直线y＝k2×x+d2的距离l1,l2：m值加1，因|l1|+|l2|＝2.338488<lmin＝106，则lmin＝|l1|+|l2|＝2.338488。因|l1|<|l2|，则取(k1,d1)＝(0.211386，1496.754459)作为所求直线的合理参数，记为(k(m),d(m))。步骤2.2.7：计算所有数据点到新直线y-k(m)x-d(m)＝0的加权距离之和f(m)(k(m),d(m))，m表示迭代次数。步骤2.2.8：判别是否为“最佳”解。因f(1)(k(1),d(1))＝5.370975，f(2)(k(2),d(2))＝0.146692,f(2)(k(2),d(2))<f(1)(k(1),d(1))，则重复步骤2.2.2-步骤2.2.8。经计算，m＝3、执行步骤2.2.6时，|l1|+|l2|＝775.635331，大于上一次迭代时的lmin＝2.338488，因此输出解(k(2),d(2))，并简记为(kx,dx)＝(0.211386，1496.754459)。步骤2.3：计算tn时刻目标在x方向上的位置pxn、速度vxn和方向kxn(n＝7)：pxn＝kxtn+dx＝0.211386×90+1496.754459＝1515.779154，vxn＝kx＝0.211386，kxn＝kx＝0.211386。步骤3：对目标y轴运动分量{(ti,yi),i＝1,2,…n}使用“加权直线参数估计模型”进行迭代估计，得到tn时刻目标在y方向上的位置pyn、速度vyn和方向kyn(n＝7)。表1.13：目标在y轴的运动分量序号t＝xy＝y1302334.162402337.713512340.854612343.765702347.226812351.137902352.37步骤3.1：对目标y轴运动分量{(ti,yi),i＝1,2,…n}采用“不加权直线参数估计模型”粗略估计在y方向上的目标运动状态方程y-k1t-d1＝0，具体过程与步骤2.1类同。解得(k1,d1)＝(0.312850，2324.980758)作为所求直线的合理参数。步骤3.2：采用“加权直线参数估计模型”迭代估计目标在y方向上的运动状态方程y-kyt-dy＝0，具体过程与步骤2.2类同，解得(ky,dy)＝(0.314830，2324.824389)。步骤3.3：计算tn时刻目标在y方向上的位置pyn、速度vyn和方向kyn(n＝7)：pyn＝kytn+dy＝0.314830×90+2324.824389＝2353.159116，vyn＝ky＝0.314830，kyn＝ky＝0.314830。步骤4：则tn时刻目标在统一直角坐标系中的位置(pxn,pyn)、目标速度vn和航向kn分别为(n＝7)：(pxn,pyn)＝(1515.779154，2353.159116)公里实施例2本实施例具体描述本发明所提出的一种针对单雷达直线航迹线的目标状态估计方法与传统卡尔曼滤波方法的对比试验过程，以及二者在执行效果上的比较。所述对比试验过程包括如下步骤：步骤1：产生模拟数据的真实值和测量值。按照表2.1中的模拟目标参数值产生匀速直线飞行状态下目标的29个位置点的真实值和测量值，结果见表2.2。接下来将按照本发明所述方法对这29个测量点逐一进行位置、速度和航向估计。表2.1：基本参数表2.2：雷达观测数据的真实值和测量值步骤2：将目标点数据{(ρi,θi,hi,ti)，i＝1,2,…29。}转换为统一直角坐标{(xi,yi,ti),i＝1,2,…29。}。xi＝xzicosδxz-yzisinδxz+xzxyi＝xzisinδxz+yzicosδxz+yzx其中：(xzx,yzx)为雷达站址在中心统一直角坐标系中的坐标。δxz为雷达站址与直角坐标系中心点的经度差(单位为弧度)。计算结果见表2.3。目标点真实位置和测量位置在中心统一直角坐标系中的显示如图2所示。表2.3：目标点的中心统一直角坐标序号时刻t真实值x真实值y测量值x测量值y101496.56462326.20471497.01742326.1440281498.25052328.58321498.00452328.77403201500.62312331.93061500.09692332.30214311502.94572335.20761503.13972334.81185401504.98332338.08231504.98852337.97386501507.00222340.93071507.23812340.76067591508.76742343.42121509.84562343.16638711511.29772346.99111510.48782347.24609811513.35492349.89351513.13172350.412710911515.45612352.85801516.08662352.6654111001517.36632355.55311516.79062355.9697121081519.10122358.00071519.99872357.4387131191521.42872361.28451521.92052361.1152141281523.36142364.01131523.51742364.1355151391525.66742367.26481526.75572366.6431161521528.36232371.06691528.70042370.8347171591529.70632372.96311528.87592373.2831181721532.47372376.86751532.64522376.8696191791533.92152378.91011533.95682378.7162201891536.07172381.94381536.40342381.7430211991538.18762384.92911538.51912384.7892222091540.28492387.88811540.28002387.6327232191542.37932390.84301541.57052391.3966242291544.48642393.81581544.07292393.7185252421547.12352397.53641547.23352397.4669262501548.80152399.90391549.22112399.6019272611551.04092403.06331550.23792403.2995282701552.85412405.62151552.37592405.7809292811555.26002409.01591555.12912409.1223步骤3：对表2.3中的测量值{(ti,xi,yi),i＝1,2,…29}使用本发明所述估计方法逐点进行目标位置(wpxi,wpyi)、速度wvi和航向wki估计。其中：参数n取值为7。表2.4是n＝7时每次参与估值计算的点序号表。表2.5是应用本发明的估计结果。表2.4：n＝7时每次参与估值计算的点表2.5：应用本发明(n＝7)逐点进行目标状态估计步骤4：对表2.3中的测量值{(ti,xi,yi),i＝1,2,…29}使用传统卡尔曼滤波方法逐点进行目标位置(kpxi,kpyi)、速度kvi和航向kki估计。其中：q、r参数取值见表2.1。表2.6是应用卡尔曼滤波方法的逐点估计结果。表2.6：应用卡尔曼方法逐点进行目标状态估计目标点真实位置、测量点位置、本发明估计位置和卡尔曼滤波估计位置在中心统一直角坐标系中的显示如图3、图4所示。图上可见，两种滤波方法得到的航迹线几乎重合，均优于测量航迹线，接下来我们对时间配准点、速度和航向等指标进一步统计、比较。步骤5：时间配准点比较为了进一步比较位置估值效果，我们引入时间配准点概念，即：相同时刻目标真实位置点与估计位置点称为一对时间配准点。ti时刻一对时间配准点之间的距离称为时间配准点间距di，d1～dm的平均值定义为tm时刻的平均时间配准点间距vdm。根据表2.3～2.6的计算结果，可以统计使用本发明方法和使用卡尔曼方法得到时间配准点间距与平均时间配准点间距值，见表2.7和表2.8，其中(rpx，rpy)表示目标真实位置，wvd29表示使用本发明方法得到的所有估值点的平均时间配准点间距，kvd29表示使用卡尔曼方法得到的所有估值点的平均时间配准点间距。表2.7：本发明时间配准点间距wd与平均间距wvd表2.8：卡尔曼时间配准点间距kd与平均间距kvd图5和图6是本发明方法和卡尔曼方法时间配准点间距与平均时间配准点间距的对比图。经过统计，有23个点的时间配准点间距wdi≤kdi，占全部点数的79.31％；所有时刻的平均时间配准点间距均有wvdi≤kvdi；全程平均时间配准点间距wvd29只占kvd29的30.15％。本发明方法在测量初期、点数较少时就能获得更加准确的目标位置估值，优势非常明显。于是可得到以下结论：在本次试验中，以目标真实位置为基准，使用本发明方法得到的目标位置点优于卡尔曼方法位置点。步骤6：速度估值比较我们定义ti时刻目标估计速度与真实速度差的绝对值为速度差dvi，dv1～dvm的平均值定义为tm时刻的平均速度差vdvm。根据表2.3～2.6的计算结果，可以统计使用本发明方法和使用卡尔曼方法得到的所有时间点上的速度差和平均速度差值，见表2.9，其中rv表示目标真实速度；wvi、wdvi、wvdvi表示使用本发明方法得到的ti时刻目标速度估值、速度差和平均速度差；kvi、kdvi、kvdvi表示使用卡尔曼方法得到的ti时刻目标速度估值、速度差和平均速度差。表2.9：速度、速度差和平均速度差图7、图8和图9是本发明方法和卡尔曼方法分别在速度、速度差和平均速度差的对比图。经过统计，有28个点的速度差wdvi≤kdvi，占全部点数的96.55％；所有时刻的平均速度差均有wvdvi≤kvdvi；全程平均速度差wvdv29只占kvdv29的8.01％。图上可以看出，本发明得到的速度估值更加准确、稳定，而且本发明方法在测量初期、点数较少时就能获得更加准确的速度估值，优势非常明显。步骤7：航向估值比较我们定义ti时刻目标估计航向与真实航向差的绝对值为航向差dki，dk1～dkm的平均值定义为tm时刻的平均航向差vdkm。根据表2.3～2.6的计算结果，可以统计使用本发明方法和使用卡尔曼方法得到的所有时间点上的航向差和平均航向差值，见表2.10，其中rk表示目标真实航向；wki、wdki、wvdki表示使用本发明方法得到的ti时刻目标航向估值、航向差和平均航向差；kki、kdki、kvdki表示使用卡尔曼方法得到的ti时刻目标航向估值、航向差和平均航向差。表2.10：航向、航向差和平均航向差图10、图11和图12是本发明方法和卡尔曼方法分别在航向、航向差和平均航向差的对比图。经过统计，有19个点的航相差wdki≤kdki，占全部点数的65.52％；所有时刻的平均航向差均有wvdki≤kvdki；全程平均航向差wvdk29只占kvdk29的50.89％。图上可以看出，本发明方法在测量初期、点数较少时就能获得更加准确的航向估值，优势较为明显。步骤8：出现野值后的估值效果比较为了模拟野值出现，我们在表2.3的第22个测量点上人为加入较大误差。图13～23分别显示了出现野值后，本发明方法和卡尔曼方法在航迹线、时间配准点、速度和航向等指标上的比较效果，表2.11是测量中包含野值与不包含野值时各项指标的统计结果。表2.11：包含/不包含野值时各项指标统计从以上图表可以看出，加入野值后，卡尔曼滤波方法在野值点处，位置、速度和航向估值均出现较大突变；另外，野值还导致了卡尔曼方法滤波发散，影响了其后几个点的目标状态估值效果。相比之下，本发明方法对野值抑制效果较为明显，无论是在野值点处，还是其后几个点上，估值效果都优于卡尔曼方法。本发明基于复杂数据环境，在无须获取系统噪声和观测噪声等先验知识的情况下，充分利用目标处于匀速直线运动状态这一特点，分别对当前n个有限测量点的x、y分量进行相对于的测量时刻的垂直距离加权迭代估计，然后根据得到的分量直线航迹线模型参数确定最终的目标状态估计(位置、速度、航向)结果。通过与传统卡尔曼滤波方法在时间配准点间距、速度差和航相差等指标的对比试验表明，本发明不但易于工程实现，而且比传统卡尔曼滤波方法具有更优的目标状态估计效果和野值抑制能力。该方法有效解决了传统卡尔曼滤波方法应用条件比较苛刻、过程噪声协方差和量测噪声协方差难以确定和系统测量中出现野值，可能导致滤波发散或目标状态估计突变等问题，极大地改善了单雷达数据质量，为在本系统中开展目标跟踪与数据融合工作提供了一种稳定、有效且易于工程实现的滤波方法。本发明估计方法科学、方案实施步骤合理，目标估计效果理想，对于提高了多雷达目标状态估计的一致性和准确性具有重要意义。本发明所提供的方法的时间复杂度和空间复杂度都很低，可操作性和实用性很强。以上所述仅是本发明的优选实施方式，应当指出，对于本
技术领域：
的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。例如但不限于以下几点：(1)为了节省计算时间，在步骤1中，可以存储前n－1个测量点的统一直角坐标，每次只对第n个测量点进行极坐标到统一直角坐标转换。(2)关于n的取值问题，我们认为一般取4～15为宜。但是，在相同数据环境中，在保证目标处于匀速直线飞行状态的前提下，n越大，目标状态估值精度就越高。(3)本发明中涉及到的目标航向也称为格式化或标准化航向，是指从正北(y轴)算起，顺时针到目标行进方向的夹角，正北为0度，正东为90度；航向取值范围：大于等于0度，小于360度。通过步骤4的计算公式得到的航向单位为弧度，且需要考虑kyn→0时的情况。因此，在实际工程应用中，存在航向kn计算时取极值和归一化到[0,360)度的问题，在此给出步骤4中更加详细的kn(单位为弧度)计算公式：当前第1页12