一种雷达组网目标状态与系统误差联合估计算法的制作方法

本发明涉及雷达组网信息融合目标跟踪领域，特别是涉及量测高精度、系统高维数、复杂强非线性条件下的机载雷达多平台协同跟踪领域，具体涉及一种雷达组网目标状态与系统误差联合估计算法。

背景技术：

雷达组网信息融合所能带来的巨大效益已经得到了世界各国的公认。尽管雷达组网信息融合技术依然在不断飞速发展，但是由于实际系统中各雷达探测目标误差的存在，系统实时融合效果的保障已成为雷达组网信息融合技术领域长期以来十分棘手的问题。实际应用表明，在多雷达组网跟踪系统中，雷达系统误差的存在会导致目标跟踪均方根误差比理论值要大。当雷达系统误差太大时，就会出现多部雷达融合跟踪甚至不如单部雷达目标跟踪效果的情况。最恶劣情况下，雷达系统误差会导致来自同一轨迹的多雷达量测互联失败，产生相对同一目标的多条航迹，这样本来是同一个目标的航迹，却由于相互偏差较大而可能被认为是不同的目标，从而给航迹关联及融合带来模糊和困难，使融合得到的系统航迹的性能下降。尤其是在目标密集、编队飞行等复杂的场景中更容易造成航迹关联混乱、融合精度降低，进而使整体系统融合失去意义，丧失了雷达组网系统本身应有的优点。因此，雷达组网系统误差配准技术是确保雷达网性能稳定所不可缺少的基础和关键技术，其目的是准确估计并消除组网雷达的系统误差。

从坐标转换的角度来看，常用的雷达组网系统误差配准算法主要可以分为两类：其中一类是基于球(极)投影的二维平面系统误差配准算法。该类算法是把各雷达的量测投影到与地球正切的局部载体坐标系上，然后再变换到二维平面来估计各雷达的系统误差。此类方法在工程上应用较多，具有算法简单、便于实现等优点，但是此类雷达系统误差配准技术也存在以下缺陷：

①球(极)投影法虽然利用高阶近似来提高精度，但由于地球是椭球而不是圆球，所以在投影时会给测量引入误差；

②球(极)投影法会使数据变形。球(极)保角投影，只保证方位角不变形，并不能保证斜距不变形，这样会导致系统误差不再是常数，而且与量测有关；

③在二维公共坐标系中只能估计方位角和径向距离的系统误差，无法估计俯仰角系统误差。

所以基于球(极)投影法的系统误差配准技术通常用于短距离雷达组网系统误差配准，而对于远距离雷达系统误差的配准，我们多采用基于地心地固(ecef)坐标系的雷达系统误差配准技术。

从数据处理的角度来看，雷达组网系统误差配准算法又可分为离线处理法和在线处理法。离线处理法主要包括最小二乘算法(ls)、广义最小二乘算法(gls)、最大似然算法(ml)和精确极大似然估计算法(eml)。该类算法需要对一段时间内的数据进行集中处理，因此计算量相对比较大。在线处理法主要是基于扩展卡尔曼滤波(ekf)和无迹卡尔曼滤波(ukf)的实时估计算法，计算量相对较小。但是ekf需要计算雅克比矩阵，在量测高精度、系统高维数、复杂强非线性的情况下往往不能很好的近似系统状态的后验概率密度函数，导致滤波精度大大降低，甚至造成滤波发散。虽然ukf避免了ekf的缺点，估计精度有一定的提高，但当ukf应用于高维非线性系统时，其确定性采样点的权值容易出现负值的情况，导致其矩积分中引入截断误差，无法保证数值稳定性和状态协方差矩阵的半正定性，造成滤波精度下降。

技术实现要素：

本发明的目的在于提供一种雷达组网目标状态与系统误差联合估计算法，克服或减轻现有技术的至少一个上述缺陷。

本发明的目的通过如下技术方案实现：一种雷达组网目标状态与系统误差联合估计算法，包括如下步骤，

步骤一：根据机载雷达测量原理，结合机载平台的地理位置及姿态角，构建出含有系统误差、测量噪声以及目标状态的每部机载雷达量测数学模型；

步骤二：根据目标运动状态，联合每部机载雷达的系统误差，构建出每一时刻含有过程噪声的扩维后目标运动状态数学模型；

步骤三：根据步骤一、步骤二构建的数学模型，设置测量噪声、过程噪声、扩维后目标运动状态初值及其估计误差协方差矩阵初值，应用ckf滤波方法，实现对目标运动状态及每部机载雷达系统误差的实时同步估计。

优选地是，所述机载雷达为两部。

优选地是，所述机载雷达测量原理通过如下公式实施：

x、y、z表示为目标相对于每部机载雷达的位置在ecef坐标系下每个坐标轴上的投影。

优选地是，所述机载平台的地理位置为机载平台所在大地坐标系下的经度、纬度及海拔高度；所述机载平台的姿态角为偏航角、纵摇角及横摇角。

优选地是，所述步骤一所构建的含有系统误差、测量噪声以及目标状态的每部机载雷达量测数学模型公式：

ri(k)为目标距离，θi(k)为目标方位角，ηi(k)为目标俯仰角；其中φi(k),αi(k)分别表示k时刻机载雷达i所在载体平台的偏航角、纵摇角及横摇角；a为由机载直角坐标系转换到end坐标系的坐标旋转矩阵表示为：

k时刻机载雷达i在大地坐标系的坐标为xisp(k)＝[li(k)λi(k)hi(k)]^t，li(k),λi(k),hi(k)分别表示k时刻机载雷达i所在的纬度、经度及海拔高度；t为由end坐标系转换到ecef坐标系的坐标旋转矩阵矩阵表示为：

xt(k)＝[xt(k)yt(k)zt(k)]^t表示为k时刻目标在ecef坐标系下的坐标；xis(k)表示为k时刻机载雷达i在ecef坐标系下的位置表示为：

其中，eq表示赤道半径，e表示地球偏心率；bi(k)为系统误差，ni(k)为测量噪声。

优选地是，所述步骤二所构建的每一时刻含有过程噪声的扩维后目标运动状态数学模型公式：

xa(k+1)＝fa(k)xa(k)+wa(k)

其中，fa(k)＝diag(f(k),i6×6)，f(k)＝diag(fx(k),fy(k),fz(k))，i＝x,y,z；wa(k)为零均值高斯白噪声；xa(k)为联合目标状态x(k)和两部雷达的系统误差b1(k),b2(k)扩维后的目标运动状态，表示为xa(k)＝[x(k)^tb1(k)^tb2(k)^t]^t。

优选地是，所述测量噪声表示为分别表示目标距离、方位角和俯仰角测量误差的协方差；过程噪声表示为qa＝diag(q,06×6)，q＝diag(qx,qy,qz)，i＝x,y,z，q表示噪声的功率谱密度。

优选地是，所述ckf滤波方法包括时间更新与量测更新两个步骤，

a、时间更新：

ⅰ、通过cholesky分解系统状态估计误差协方差矩阵p(k|k)；

p(k|k)＝s(k|k)s^t(k|k)

ⅱ、根据spherical-radialcubature准则计算cubature点；

ⅲ、通过状态方程传播cubature点；

ⅳ、计算系统状态的先验均值及先验协方差矩阵；

ⅴ、计算系统状态的先验协方差矩阵；

b、量测更新：

ⅰ、通过cholesky分解系统状态先验协方差矩阵p(k+1|k)；

p(k+1|k)＝s(k+1|k)s^t(k+1|k)

ⅱ、根据spherical-radialcubature准则计算cubature点；

ⅲ、通过量测方程传播cubature点；

zi(k+1|k)＝ha(xi(k+1|k))

ⅳ、计算量测的先验均值；

ⅴ、计算量测的先验协方差矩阵；

ⅵ、计算量测和状态向量的互相关协方差矩阵；

ⅶ、计算卡尔曼滤波增益；

ⅷ、计算系统状态的后验均值；

ⅸ、计算系统状态的后验协方差矩阵；

p(k+1|k+1)＝p(k+1|k)-w(k+1)pzz(k+1|k)w^t(k+1)。

本发明所提供的一种雷达组网目标状态与系统误差联合估计算法的有益效果在于，首先采用基于ecef坐标系的雷达系统误差配准技术，解决了传统的基于坐标投影的二维平面系统误差配准算法中存在的诸如数据变形、无法估计俯仰角系统误差及不适于远距离系统误差配准等固有问题。其次，该方法对系统状态进行扩维处理，把系统误差作为系统的未知、待估计状态，并引入容积卡尔曼滤波(ckf)算法，实现了对目标状态和系统误差的联合估计。本发明所提出的方法无需计算雅克比矩阵，而且在递推运算中具有很强的数值稳定性，解决了传统非线性高斯滤波器在量测高精度、系统高维数、复杂强非线性情况下对系统状态后验概率密度函数近似精度不高这一问题，提高了雷达组网条件下对系统误差的估计精度及目标的跟踪精度，有效地实现了雷达组网系统误差的实时配准，同时提升了目标跟踪的可靠性和稳定性。同时本发明具有良好的扩展性和适应性，可广泛应用于火控、监视、预警等战术功能，对稳定目标跟踪有较高要求的多平台有源/无源雷达信息融合跟踪系统，应用前景广阔，应用价值巨大。

附图说明

图1是本发明雷达组网目标状态与系统误差联合估计算法的流程图；

图2是根据本发明一个实施例中的仿真环境经纬图；

图3是根据本发明一个实施例中的第一雷达距离系统误差估计效果图；

图4是根据本发明一个实施例中的第二雷达距离系统误差估计效果图；

图5是根据本发明一个实施例中的第一雷达方位角系统误差估计效果图；

图6是根据本发明一个实施例中的第二雷达方位角系统误差估计效果图；

图7是根据本发明一个实施例中的第一雷达俯仰角系统误差估计效果图；

图8是根据本发明一个实施例中的第二雷达俯仰角系统误差估计效果图；

图9是根据本发明一个实施例中的第一雷达径向距离估计效果图；

图10是根据本发明一个实施例中的第二雷达径向距离估计效果图；

图11是根据本发明一个实施例中的第一雷达方位角估计效果图；

图12是根据本发明一个实施例中的第二雷达方位角估计效果图；

图13是根据本发明一个实施例中的第一雷达俯仰角估计效果图；

图14是根据本发明一个实施例中的第二雷达俯仰角估计效果图。

具体实施方式

为使本发明实施的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行更加详细的描述。在附图中，自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。所描述的实施例是本发明一部分实施例，而不是全部的实施例。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

下面结合附图对本发明的雷达组网目标状态与系统误差联合估计算法做进一步详细说明。

一种雷达组网目标状态与系统误差联合估计算法，其核心在于机载雷达系统的数学建模及非线性滤波算法的设计。首先完成对机载雷达系统数学模型的构建；然后采用基于ecef坐标系的系统误差配准技术，以解决传统基于坐标投影的二维平面系统误差配准算法中存在的诸如数据变形、无法估计俯仰角系统误差及不适于远距离系统误差配准等固有问题。最后对系统状态进行扩维处理，把系统误差作为系统未知、待估计的状态，并引入ckf滤波算法，实现对目标状态和系统误差的联合估计。见图1所示，具体包括以下几个步骤：

1)步骤一：根据机载雷达测量原理，结合机载平台的地理位置及姿态角，构建出含有系统误差、测量噪声以及目标状态的每部机载雷达量测数学模型。

本发明中的系统选择由两部三坐标机载雷达组成。机载雷达i(i＝1,2)在极坐标系下对同一目标进行同步测量，存在距离系统误差方位角系统误差俯仰角系统误差且假设它们均为常量加性系统误差，表示为

机载雷达i在k时刻测量得到含有系统误差bi(k)和测量噪声ni(k)的目标距离ri(k)、目标方位角θi(k)和目标俯仰角ηi(k)。其中ni(k)为零均值高斯白噪声，对应的噪声协方差矩阵表示为分别表示目标距离、目标方位角和目标俯仰角测量误差的协方差，且系统误差与测量噪声之间是相互独立的。

假设不含系统误差的真值量测表示为ri'(k),θi'(k),ηi'(k)，则k时刻机载雷达i在极坐标系下的测量方程为：

将极坐标系下的测量转化到局部直角载体坐标系下，可得：

式中h^-1(r,θ,η)＝[rsinθcosη,rcosθcosη,rsinη]^t。

假设k时刻机载雷达i所在的载体平台姿态角表示为其中φi(k),αi(k)分别表示k时刻机载雷达i所在载体平台的偏航角、纵摇角及横摇角。则根据载体坐标系和end坐标系间的转换关系，可把机载雷达i在载体坐标系的测量转换到end坐标系中：

式中，

假设k时刻机载雷达i在大地坐标系的坐标为xisp(k)＝[li(k)λi(k)hi(k)]^t，li(k),λi(k),hi(k)分别表示k时刻机载雷达i所在的纬度、经度及海拔高度。则根据大地坐标系和ecef坐标系间的转换关系，可把k时刻机载雷达i在大地坐标系的位置转换到ecef坐标系中：

式中，eq表示赤道半径，e表示地球偏心率。根据end坐标系和ecef坐标系间的转换关系，将k时刻机载雷达i在end坐标系的量测转换到ecef坐标系中：

式中

因此，k时刻机载雷达i在ecef坐标系下的量测可以表示为：

假设k时刻目标在ecef坐标系下的坐标为xt(k)＝[xt(k)yt(k)zt(k)]^t，如果机载雷达i不存在系统误差，且不含测量噪声，那么此时机载雷达i在ecef坐标系下的量测就应该等同于目标此时在ecef坐标系下的坐标，即：

xis(k)+t(xisp(k))a(vi(k))h^-1(ri'(k),θi'(k),ηi'(k))＝xt(k)

经过整理，可得机载雷达的量测方程为：

其中，x、y、z表示为目标相对于每部机载雷达的位置在ecef坐标系下每个坐标轴上的投影。

2)步骤二：根据目标运动状态，联合每部机载雷达的系统误差，构建出每一时刻含有过程噪声的扩维后目标运动状态数学模型。

首先将系统的状态进行扩维，把系统误差作为系统未知、待估计的状态。联合目标状态x(k)和两部雷达的系统误差b1(k),b2(k)，构建新的系统状态(即扩维后目标运动状态)为：

xa(k)＝[x(k)^tb1(k)^tb2(k)^t]^t

由于假设机载雷达的系统误差为不变常量，因此可得：

bi(k+1)＝i3×3bi(k)

扩维后，新的系统状态xa(k)的状态转移方程可表示为：

xa(k+1)＝fa(k)xa(k)+wa(k)

其中，fa(k)＝diag(f(k),i6×6)，f(k)＝diag(fx(k),fy(k),fz(k))，i＝x,y,z。wa(k)＝diag(w(k),06×6)，wa(k)为零均值高斯白噪声，其噪声协方差矩阵为qa＝diag(q,06×6)，q＝diag(qx,qy,qz)，i＝x,y,z，q表示噪声的功率谱密度。

3)根据步骤一、步骤二构建的数学模型，设置测量噪声、过程噪声、扩维后目标运动状态初值及其估计误差协方差矩阵初值，应用ckf滤波方法，实现对目标运动状态及每部机载雷达系统误差的实时同步估计。

首先将步骤二中构建的新的系统状态和机载雷达1、2的量测方程联合，可得：

za(k)＝ha(xa(k))+na(k)

式中，na(k)仍为零均值高斯噪声，对应的噪声协方差矩阵表示为ra(k)＝diag(r1(k),r2(k))。

接下来对扩维后的系统状态方程和量测方程应用ckf滤波方法，进而对系统状态及系统误差进行联合估计。其中，ckf滤波方法包括时间更新与量测更新两个步骤。

a、时间更新：

ⅰ、通过cholesky分解系统状态估计误差协方差矩阵p(k|k)；

p(k|k)＝s(k|k)s^t(k|k)

ⅱ、根据spherical-radialcubature准则计算cubature点；

ⅲ、通过状态方程传播cubature点；

ⅳ、计算系统状态的先验均值及先验协方差矩阵；

ⅴ、计算系统状态的先验协方差矩阵；

b、量测更新：

ⅰ、通过cholesky分解系统状态先验协方差矩阵p(k+1|k)；

p(k+1|k)＝s(k+1|k)s^t(k+1|k)

ⅱ、根据spherical-radialcubature准则计算cubature点；

ⅲ、通过量测方程传播cubature点；

zi(k+1|k)＝ha(xi(k+1|k))

ⅳ、计算量测的先验均值；

ⅴ、计算量测的先验协方差矩阵；

ⅵ、计算量测和状态向量的互相关协方差矩阵；

ⅶ、计算卡尔曼滤波增益；

ⅷ、计算系统状态的后验均值；

ⅸ、计算系统状态的后验协方差矩阵；

p(k+1|k+1)＝p(k+1|k)-w(k+1)pzz(k+1|k)w^t(k+1)。

需要说明的是，假设在加性噪声条件下，考虑如下离散形式的系统状态空间模型：

其中，表示k时刻系统状态向量(nx为状态维数)，表示k时刻外部量测向量(nz为量测维数)。fa(·)表示状态转移函数，ha(·)表示量测函数。与分别表示系统噪声和量测噪声，二者互不相关且均为零均值高斯白噪声，噪声协方差矩阵分别为qa(k)和ra(k)。

下面结合一个实施例通过其数值仿真详细说明该发明所提出的方法。

假设系统由两部机载雷达1、2和一个空中目标构成。机载雷达平台设置为直升机平台，飞行速度较慢，飞行海拔高度为1km，目标设置为战斗机目标，飞行速度较快，飞行海拔高度为2km。机载雷达扫描周期为1s，目标沿经线飞行，具体仿真场景如图2所示。

各机载雷达的量测噪声ni(k)均为零均值高斯白噪声，目标距离、目标方位角和目标俯仰角测量误差的协方差分别为雷达1机载平台姿态角变化规律可以描述为v1(k)＝[0.002k,0.01+0.002k,0.01+0.002k]^t，雷达2机载平台姿态角变化规律可以描述为v2(k)＝[0.002k,0.001k,0.001k]^t，各个雷达的系统误差都为bi＝[1000m,0.0087rad,0.0087rad]^t。下面对ecef-ckf-asr算法的估计性能进行仿真分析。对ecef-ckf-asr算法进行50次蒙特卡洛仿真结果如图3～图14所示。可以看出，随着ecef-ckf-asr算法估计的收敛，ecef-ckf-asr算法给出的目标状态估计基本消除了系统误差的影响。

表1ecef-ckf-asr算法系统误差估计精度

从表1可以看出，ecef-ckf-asr算法对各个系统误差的估计精度基本上都达到95％以上，可实现算法收敛。因此ecef-ckf-asr算法对各个雷达的系统误差具有良好的估计效果，很好地解决了机动雷达的误差配准问题。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。