利用ADS-B数据的GM-EPHD滤波雷达系统误差联合估计方法

摘要

一种利用ADS-B数据的GM-EPHD滤波雷达系统误差联合估计方法。其首先将目标ADS-B观测值转换到以雷达站位置为中心的直角坐标系下，建立目标的ADS-B及雷达的观测方程，在线性高斯条件下建立目标及雷达系统误差的状态转移方程，将雷达系统误差扩展至目标状态中组成扩展状态，先使用坐标转换后的ADS-B观测数据对目标状态进行高斯混合-概率假设密度滤波，以获得较为准确的目标状态估计，再使用雷达观测数据对目标状态更新后的扩展状态进行高斯混合-扩展概率假设密度滤波，进而采用两步卡尔曼滤波器联合估计出目标状态和雷达系统误差，最后使用加权平均算法得到雷达系统误差的融合估计结果。本方法估计精度高、估计性能好。

说明书

技术领域

本发明属于传感器误差配准技术领域，特别是涉及一种利用ADS-B数据的 GM-EPHD滤波雷达系统误差联合估计方法。

背景技术

雷达误差一般分为两类：随机误差和系统误差。在多雷达融合跟踪系统中，雷达系 统误差估计已成为多雷达融合处理的先决条件，将会直接影响整个系统的工作性能。因 此需要对雷达系统误差进行估计，以此对雷达测量进行相应的补偿，这个过程也称为误 差配准。雷达系统误差包括以下几种：雷达站点坐标标定误差、空间坐标转换误差、雷 达测量误差。常用的雷达系统误差估计方法一般只考虑因测距、测角引起的测量误差。 现有的雷达系统误差估计方法可以归纳为离线估计方法和在线估计方法。离线估计方法 是通过对一段时间的雷达观测数据进行数据拟合，从而估计出雷达的系统误差，如最小 二乘法及最大似然法等。在线估计方法主要是利用滤波方法实现雷达系统误差的递推估 计，在线估计方法和离线估计方法相比具有实时误差估计的优点，因此得到国内外学者 的更多关注。2006年Bar-Shalom提出利用Kalman滤波进行雷达系统误差估计的方法。 2007年Herrero提出一种利用Kalman滤波实现目标状态与系统误差联合估计方法。

无论是基于数据拟合的离线估计方法还是基于滤波的在线估计方法，上述方法都需 要满足这样一个假设：目标状态和观测之间的关联关系预先已知。现有的方法通过最近 邻方法(NN)、联合概率数据关联方法(JPDA)等获得目标状态和观测之间的关联关 系。而对于多目标或密集杂波场景，想要获得准确的关联关系是十分困难的，而且错误 的目标状态与观测之间的关联关系将会严重影响雷达系统误差估计结果。2003年Mahler 在随机有限集理论框架下提出传递目标状态集合的后验概率密度的一阶统计矩的概率假 设密度(Probability Hypothesis Density，PHD)滤波理论。PHD滤波器将复杂的多目标状态 空间的运算转换为单目标状态空间内的运算，可有效避免多目标状态估计中复杂的数据 关联过程。2006年Vo提出了高斯混合模型PHD(Gaussian Mixture PHD,GM-PHD)滤波 器，给出了线性高斯条件的PHD滤波器的封闭解形式。2012年Wenling Li等将GM-PHD 滤波器用于多传感器多目标跟踪中，实现了目标状态和传感器系统误差的联合估计。

随着全球导航卫星系统(GNSS)和空空、空地数据链通信技术的发展，一种新型的 航空器运行监视技术—广播式自动相关监视技术(Automatic Dependent  Surveillance-Broadcast，ADS-B)正在航空器监视中广泛应用。航空器上的ADS-B机载收 发信机将本机GPS(Global Positioning System)导航设备获得的本机经度、纬度、速度、时 间、高度等数据通过数据链对外广播，ADS-B地面站通过接收有效空域内航空器所发送 的广播数据实现对航空器的监视，其定位精度即为机载GPS导航设备的定位精度，远远 优于雷达定位精度。因此利用ADS-B进行雷达系统误差估计成为该领域新的研究热点。 2009年Besada提出了利用ADS-B监视数据对空中交通管理系统航管雷达进行误差配准 的方法。2013年He You提出利用Kalman滤波及ADS-B监视数据进行雷达系统误差估 计的方法。但到目前为止尚未发现利用ADS-B的概率假设密度滤波雷达系统误差联合估 计方面的报道。

发明内容

为了解决上述问题，本发明的目的在于提供一种利用ADS-B数据的GM-EPHD滤波 雷达系统误差联合估计方法。

为了达到上述目的，本发明提供的利用ADS-B数据的GM-EPHD滤波雷达系统误差 联合估计方法包括按顺序进行的下列步骤：

1)建立目标的ADS-B及雷达的观测方程的S1阶段；

2)建立目标及雷达系统误差的状态转移方程的S2阶段；

3)利用GM-EPHD滤波器对雷达系统误差进行融合估计的S3阶段。

在步骤1)中，所述的建立目标的ADS-B及雷达的观测方程的方法是首先利用坐标 投影技术将目标的ADS-B观测值转换到以雷达站位置为中心的直角坐标系下，由于 ADS-B监视数据的定位精度远远优于雷达的定位精度，在不考虑ADS-B定位误差的基础 上，建立坐标转换后的目标的ADS-B观测方程，然后考虑雷达系统误差，建立目标的雷 达观测方程，为利用GM-EPHD滤波器对雷达系统误差进行融合估计做准备。

在步骤2)中，所述的建立目标及雷达系统误差的状态转移方程的方法是假设目标状 态及雷达系统误差均满足线性高斯特性，进而建立目标及雷达系统误差的状态转移方程， 为利用GM-EPHD滤波器对雷达系统误差进行融合估计做准备。

在步骤3)中，所述的利用GM-EPHD滤波器对雷达系统误差进行融合估计的方法 是首先将雷达系统误差扩展至目标状态中组成扩展状态，根据S2阶段中建立的目标及雷 达系统误差的状态转移方程分别对目标状态及雷达系统误差进行预测，然后使用S1阶段 中坐标转换后的ADS-B观测数据对目标状态进行高斯混合-概率假设密度滤波，以获取 较为准确的目标状态估计，接着使用S1阶段中的雷达观测数据对扩展状态进行高斯混合 -扩展概率假设密度滤波，进而采用两步卡尔曼滤波器联合估计出目标状态和雷达系统误 差，最后使用加权平均算法对雷达系统误差进行融合估计。

本发明提供的利用ADS-B数据的GM-EPHD滤波雷达系统误差联合估计方法首先将 目标的ADS-B观测值转换到以雷达站位置为中心的直角坐标系下，由于ADS-B监视数 据的定位精度远远优于雷达的定位精度，在不考虑ADS-B定位误差的基础上，建立目标 的ADS-B及雷达的观测方程，然后在线性高斯条件下建立目标及雷达系统误差的状态转 移方程，接着将雷达系统误差扩展至目标状态中组成扩展状态，先使用坐标转换后的 ADS-B观测数据对目标状态进行高斯混合-概率假设密度滤波，以获得较为准确的目标状 态估计，再使用雷达观测数据对目标状态更新后的扩展状态进行高斯混合-扩展概率假设 密度滤波，进而采用两步卡尔曼滤波器联合估计出目标状态和雷达系统误差，最后使用 加权平均算法得到雷达系统误差的融合估计结果。本发明方法具有估计精度高、估计性 能好等优点。

附图说明

图1为本发明提供的利用ADS-B数据的GM-EPHD滤波雷达系统误差联合估计方法 流程图。

图2为目标真实运动轨迹及雷达和ADS-B观测数据。

图3为目标真实运动轨迹及状态估计结果。

图4为目标状态估计的OSPA距离。

图5为目标数估计结果。

图6、图7分别为采用本发明方法和文献方法获得的雷达距离、方位误差估计比较图。

图8、图9分别为采用本发明方法和文献方法获得的雷达距离、方位误差估计的均方 根误差比较图。

图10为采用本发明方法对雷达系统误差校正后的雷达观测图。

具体实施方式

下面结合附图和具体实施例对本发明提供的利用ADS-B数据的GM-EPHD滤波雷达 系统误差联合估计方法进行详细说明。

图1为本发明提供的利用ADS-B数据的GM-EPHD滤波雷达系统误差联合估计方法 流程图。其中的全部操作都是在计算机系统中完成的，操作的主体均为计算机系统。

如图1所示，本发明提供的利用ADS-B数据的GM-EPHD滤波雷达系统误差联合估 计方法包括按顺序进行的下列步骤：

1)建立目标的ADS-B及雷达的观测方程的S1阶段：

本阶段是首先利用坐标投影技术将目标的ADS-B观测值转换到以雷达站位置为中心 的直角坐标系下，由于ADS-B监视数据的定位精度远远优于雷达的定位精度，在不考虑 ADS-B定位误差的基础上，建立坐标转换后的目标的ADS-B观测方程，接着考虑雷达系 统误差，建立目标的雷达观测方程，然后进入下一步S2阶段。

在此阶段中，所述的建立目标的ADS-B及雷达的观测方程的具体方法如下：由于 ADS-B使用WGS-84大地坐标系获得目标的经度、纬度及高度信息，而雷达使用以雷达 位置为中心的极坐标系获得目标的距离及方位信息。首先将目标的ADS-B观测值转换到 以雷达站位置为中心的直角坐标系下，利用ADS-B定位精度即GPS导航设备定位精度 远远优于雷达定位精度的优点，在不考虑ADS-B定位误差的基础上，建立坐标转换后的 目标的ADS-B观测方程，然后考虑雷达系统误差，建立目标的雷达观测方程。

a、建立坐标转换后的目标的ADS-B观测方程

首先将目标的ADS-B观测值转换到地心地固(ECEF)直角坐标系下，然后再转换成以 雷达站位置为中心的直角坐标系下。设ADS-B系统获得的k时刻第i个目标的观测值为 分别表示目标的经度、纬度及距离 海平面的高度，转换后的ECEF坐标系观测值为则：

$\left(\begin{array}{c}{\mathrm{Ex}}_{k,\mathrm{ADS}-B}^{\left(i\right)}=\left({C+H}_{k,\mathrm{ADS}-B}^{\left(i\right)}\right)\cos L_{k,\mathrm{ADS}-B}^{\left(i\right)}\cos {\lambda }_{k,\mathrm{ADS}-B}^{\left(i\right)}\\ {\mathrm{Ey}}_{k,\mathrm{ADS}-B}^{\left(i\right)}=(C+H_{k,\mathrm{ADS}-B}^{\left(i\right)})\cos L_{k,\mathrm{ADS}-B}^{\left(i\right)}\sin {\lambda }_{k,\mathrm{ADS}-B}^{\left(i\right)}\\ {\mathrm{Ez}}_{k,\mathrm{ADS}-B}^{\left(i\right)}=[C(1-e^2)+H_{k,\mathrm{ADS}-B}^{\left(i\right)}]\sin L_{k,\mathrm{ADS}-B}^{\left(i\right)}\end{array}\right)---\left(1\right)$

其中，e为地球的偏心率，E_q为赤道半径。

若雷达站所在位置在WGS-84坐标系下可表示为(L_RS,λ_RS,H_RS)，则其ECEF坐标系 下的位置(Ex_RS,Ey_RS,Ez_RS)同样可以由式(1)获得。

将ECEF坐标系下的目标ADS-B观测值通过式(2)转换 为以雷达站位置为中心的直角坐标系下的观测值

$\left(\begin{array}{c}{x}_{k,\mathrm{ADS}-B}^{\left(i\right)}\\ y_{k,\mathrm{ADS}-B}^{\left(i\right)}\\ z_{k,\mathrm{ADS}-B}^{\left(i\right)}\end{array}\right)=\left(\begin{array}{c}{\mathrm{Ex}}_{k,\mathrm{ADS}-B}^{\left(i\right)}-{\mathrm{Ex}}_{\mathrm{RS}}\\ {\mathrm{Ey}}_{k,\mathrm{ADS}-B}^{\left(i\right)}-{\mathrm{Ey}}_{\mathrm{RS}}\\ {\mathrm{Ez}}_{k,\mathrm{ADS}-B}^{\left(i\right)}-{\mathrm{Ez}}_{\mathrm{RS}}\end{array}\right)+T_{\mathrm{RS}}\left(\begin{array}{c}{\mathrm{Ex}}_{k,\mathrm{ADS}-B}^{\left(i\right)}\\ {\mathrm{Ey}}_{k,\mathrm{ADS}-B}^{\left(i\right)}\\ {\mathrm{Ez}}_{k,\mathrm{ADS}-B}^{\left(i\right)}\end{array}\right)---\left(2\right)$

其中，T_RS为旋转矩阵。

$T_{\mathrm{RS}}=\left(\begin{array}{ccc}-\sin {\lambda }_{\mathrm{RS}}\cos L_{\mathrm{RS}}& -\sin {\lambda }_{\mathrm{RS}}\sin L_{\mathrm{RS}}& \cos {\lambda }_{\mathrm{RS}}\\ \sin L_{\mathrm{RS}}& \cos L_{\mathrm{RS}}& 0\\ \cos {\lambda }_{\mathrm{RS}}\cos L_{\mathrm{RS}}& \cos {\lambda }_{\mathrm{RS}}\sin L_{\mathrm{RS}}& \sin {\lambda }_{\mathrm{RS}}\end{array}\right)$

由于ADS-B的定位精度远远优于雷达的定位精度，假设不考虑ADS-B定位误差， 并且忽略坐标转换误差，则由式(2)表示的以雷达位置为中心的直角坐标系下的ADS-B二 维空间测量值可以写成：

$\left(\begin{array}{c}{x}_{k,\mathrm{ADS}-B}^{\left(i\right)}\\ y_{k,\mathrm{ADS}-B}^{\left(i\right)}\end{array}\right)=\left(\begin{array}{c}{x}_k^i\\ y_k^i\end{array}\right)+\left(\begin{array}{c}{n}_{x,k}\\ n_{y,k}\end{array}\right)---\left(3\right)$

将式(3)用矩阵形式表示为：

$\left(\begin{array}{c}{z}_{k,\mathrm{ADS}-B}^i=h_{\mathrm{ADS}-B}\left(x_k^i\right)+n_k^{\mathrm{ADS}-B}\\ =H_k^{\mathrm{ADS}-B}+n_k^{\mathrm{ADS}-B}\end{array}\right)---\left(4\right)$

其中，表示以雷达为中心的直角坐标系下由目标的位置、速度 组成的状态矢量，$n_k^{\mathrm{ADS}-B}={[n_{x,k},n_{y,k}]}^T$表示随机噪声，且

$h_{\mathrm{ADS}-B}\left(x_k^i\right)=\left(\begin{array}{c}{x}_k^i\\ y_k^i\end{array}\right)---\left(5\right)$

由于监视区域内的所有目标均配备ADS-B机载设备，则k时刻可以获得N_k个目标的 ADS-B观测值，即$Z_k^{\mathrm{ADS}-B}=\{z_{k,\mathrm{ADS}-B,}^1{z}_{k,\mathrm{ADS}-B,}^2,...,z_{k,\mathrm{ADS}-B,}^{N_k}\}.$

b、建立目标的雷达观测方程

雷达使用以雷达位置为中心的极坐标系获得目标的距离和方位信息，雷达的观测模 型可以用式(6)表示：

$\left(\begin{array}{c}\widetilde{r_k^{\left(i\right)}}=r_k^{\left(i\right)}+\mathrm{\Delta r}+{\mathrm{nr}}_k^{\left(i\right)}\\ \widetilde{{\theta }_k^{\left(i\right)}}={\theta }_k^{\left(i\right)}+\mathrm{\Delta \theta }+n{\theta }_k^{\left(i\right)}\end{array}\right)---\left(6\right)$

其中，分别表示k时刻第i个目标的雷达距离、方位角测量值，分 别表示距离、方位对应的真实值，△r、△θ分别表示雷达距离、方位的系统误差，分别表示雷达测量距离、方位的观测噪声。

将式(6)用矩阵形式表示为：

$\left(\begin{array}{c}{z}_{k,\mathrm{radar}}^i=h_{\mathrm{radar}}\left(x_k^i\right)+b_k+n_k^{\mathrm{radar}}\\ =H_k^{\mathrm{radar}}+b_k+n_k^{\mathrm{radar}}\end{array}\right)---\left(7\right)$

其中，表示以雷达为中心的直角坐标系下由目标的位置、速度 组成的状态矢量，b_k＝[△r,△θ]^T表示由雷达距离、方位组成的系统误差， $n_k^{\mathrm{radar}}={[{\mathrm{nr}}_k^{\left(i\right)},{\mathrm{n\theta }}_k^{\left(i\right)}]}^T$表示雷达距离、方位的观测噪声，且

$h_{\mathrm{radar}}\left(x_k^i\right)=\left(\begin{array}{c}\sqrt{{\left(x_k^i\right)}^2+{\left(y_k^i\right)}^2}\\ \arctan \frac{y_k^i}{x_k^i}\end{array}\right)---\left(8\right)$

假设监视区域N_k个目标中有M_k个目标被雷达观测到(由于雷达有一定的检测概率， 一般情况下N_k≥M_k)，则k时刻可以获得M_k个目标的雷达观测值，即 $Z_k^{\mathrm{radar}}=\{z_{k,\mathrm{radar}}^1,z_{k,\mathrm{radar}}^2,...,z_{k,\mathrm{radar}}^{M_k}\}.$

2)建立目标及雷达系统误差的状态转移方程的S2阶段：

本阶段假设目标状态及雷达系统误差模型均满足线性高斯特性，进而建立目标及雷 达系统误差的状态转移方程，然后进入下一步S3阶段。

在此阶段中，所述的建立目标及雷达系统误差的状态转移方程的具体方法如下：假 设目标状态及雷达系统误差模型均满足线性高斯特性，考虑目标在二维平面内做线性运 动，则目标状态及雷达系统误差转移方程可以分别表示为

$x_k^i=F_{k-1}{x}_{k-1}^i+G{w}_{k-1}---\left(9\right)$

$b_k^i=b_{k-1}^i+v_{k-1}---\left(10\right)$

其中，F_k-1表示目标的状态转移矩阵，G为系统噪声系数矩阵，w_k-1为状态过程噪声， 且w_k-1～N(·,0,Q_k)，v_k-1为雷达系统误差过程噪声，且v_k-1～N(·,0,Q_b,k)。

$F_{k-1}=\left(\begin{array}{cccc}1& 0& \mathrm{\Delta t}& 0\\ 0& 1& 0& \mathrm{\Delta t}\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right)$

$G=\left(\begin{array}{cc}\frac{\Delta t^2}{2}& 0\\ 0& \frac{\Delta t^2}{2}\\ \mathrm{\Delta t}& 0\\ 0& \mathrm{\Delta t}\end{array}\right)$

$Q_k=\left(\begin{array}{cc}{\sigma }_x^2& 0\\ 0& {\sigma }_y^2\end{array}\right)$

$Q_{b,k}=\left(\begin{array}{cc}{\sigma }_r^2& 0\\ 0& {\sigma }_{\theta }^2\end{array}\right)$

3)利用GM-EPHD滤波器对雷达系统误差进行融合估计的S3阶段：

本阶段假设目标的状态模型及观测模型均满足线性高斯特性，目标的存活概率与目 标检测概率相互独立，并且假设雷达系统误差也满足线性高斯特性。为了实现雷达系统 误差和目标状态的联合估计，将雷达系统误差扩展至目标状态中组成扩展状态，使用高 斯混合-概率假设密度滤波器递归跟踪扩展状态。在此过程中，由于ADS-B的定位精度 远远优于雷达的定位精度，先使用坐标转换后的ADS-B观测数据对目标状态进行高斯混 合-概率假设密度滤波，以获得较为准确的目标状态估计，然后使用雷达观测数据对目标 状态更新后的扩展状态进行高斯混合-扩展概率假设密度滤波。由于雷达系统误差和目标 状态在雷达测量更新步骤通过似然函数耦合在一起，此时采用两步卡尔曼滤波器可同时 估计出雷达系统误差和目标状态。

在此阶段中，所述的利用GM-EPHD滤波器对雷达系统误差进行融合估计的具体方 法如下：

a、假设扩展状态新生、衍生过程的强度函数可以写成高斯混合形式

${\gamma }_k(x_k,b_k)=\underset{j=1}{\overset{J_{\gamma ,k}}{\Sigma }}{w}_{\gamma ,k}^{j}N\left(\right[x_k;b_k];[m_{\gamma ,k}^j;b_{\gamma ,k}^j],[P_{\gamma ,k}^j;B_{\gamma ,k}^j\left]\right)---\left(11\right)$

${\beta }_{k|k-1}(x_k,b_k|x_{k-1},b_{k-1})=\underset{i=1}{\overset{J_{\beta ,k}}{\Sigma }}{w}_{\beta ,k}^{i}N\left(\right[x_k;b_k];[x_{k-1};b_{k-1}],[Q_{x,k}^i;Q_{b,k}^i\left]\right)---\left(12\right)$

其中，J_γ,k表示新生目标的个数，分别表示第j个高斯混合形 式新生目标强度函数的均值、雷达系统误差均值、系统噪声协方差、误差协方差。J_β,k表 示衍生目标的个数，x_k-1、b_k-1、分别表示第i个高斯混合形式衍生目标强度 函数的均值、雷达系统误差均值、系统噪声协方差、误差协方差。

b、预测步：

假设扩展状态后验强度函数可以写成高斯混合形式

${\nu }_{k-1|k-1}=\underset{j=1}{\overset{J_{k-1}}{\Sigma }}{\omega }_{k-1}^{j}N\left(\right[x_{k-1};b_{k-1}];[m_{k-1|k-1}^j;b_{k-1|k-1}^j],[P_{k-1|k-1}^j;B_{k-1|k-1}^j\left]\right)---\left(13\right)$

则预测的强度函数也可以写成高斯混合形式

υ_k|k-1＝υ_S,k|k-1+υ_β,k|k-1+γ_k(x_k,b_k)           (14)

其中，υ_S,k|k-1表示扩展状态存活过程的强度函数，υ_β,k|k-1表示扩展状态衍生过程的强度函 数，γ_k(x_k,b_k)表示扩展状态新生过程的强度函数，且由式(14)给出。

${\nu }_{S,k|k-1}=p_S\underset{j=1}{\overset{J_{k-1}}{\Sigma }}{\omega }_{k-1}^{j}N\left(\right[x_k;b_k];[m_{S,k|k-1}^j;b_{S,k|k-1}^j],[P_{S,k|k-1}^j;B_{S,k|k-1}^j\left]\right)$

${\nu }_{\beta ,k|k-1}=\underset{j=1}{\overset{J_{k-1}}{\Sigma }}\underset{i=1}{\overset{J_{\beta ,k}}{\Sigma }}{\omega }_{k-1}^j{\omega }_{\beta ,k}^{i}N(x_k;m_{\beta ,k|k-1}^{j,i},P_{\beta ,k|k-1}^{j,i})N(b_k;b_{\beta ,k|k-1}^{j,i},B_{\beta ,k|k-1}^{j,i})$

$m_{S,k|k-1}^j=F_{k-1}{m}_{k-1|k-1}^j$

$P_{S,k|k-1}^j=F_{k-1}{P}_{k|k-1}^j{F}_{k-1}^T+Q_{x,k-1}$

$m_{\beta ,k|k-1}^{j,i}=m_{k-1|k-1}^j$

$m_{\beta ,k|k-1}^{j,i}=P_{k|k-1}^j+Q_{x,k}^i$

$b_{S,k|k-1}^j=b_{k-1|k-1}^j$

$B_{S,k|k-1}^j=B_{k-1|k-1}^j+Q_{b,k-1}$

$b_{\beta ,k|k-1}^{j,i}=b_{k-1|k-1}^j$

$B_{\beta ,k|k-1}^{j,i}=B_{k-1|k-1}^j+Q_{b,k}^i$

由于雷达系统误差和目标状态相互独立，在预测步中，采用卡尔曼滤波器分别传递 雷达系统误差和目标状态的高斯混合形式强度函数的权值、均值及协方差。

c、更新步：

假设预测强度函数也可以表示成如式(15)所示的高斯混合形式

${\upsilon }_{k|k-1}=\underset{j=1}{\overset{J_{k|k-1}}{\Sigma }}{\omega }_{k|k-1}^{j}N\left(\right[x_k;b_k];[m_{k|k-1}^j;b_{k|k-1}^j],[P_{k|k-1}^j;B_{k|k-1}^j\left]\right)---\left(15\right)$

则使用坐标转换后的ADS-B测量数据对预测强度函数进行更新，得到的后验强度函 数为

${\upsilon }_{k|k}^{\mathrm{ADS}-B}=\underset{z_k\in Z_k^{\mathrm{ADS}-B}}{\Sigma }{\upsilon }_{d,k}^{\mathrm{ADS}-B}(x_k;z_k)---\left(16\right)$

其中

${\upsilon }_{d,k}^{\mathrm{ADS}-B}(x_k;z_k)=\underset{j=1}{\overset{J_{k|k-1}}{\Sigma }}{\omega }_k^{j,\mathrm{ADS}-B}\left(z_k\right)N(x_k;m_{k|k}^{j,\mathrm{ADS}-B}\left(z_k\right),P_{k|k}^{j,\mathrm{ADS}-B})\times N(b_k;b_{k|k-1}^j,B_{k|k-1}^j)$

ADS-B测量数据更新后的高斯混合形式目标强度函数的权值、均值及协方差分别为

${\omega }_k^{j,\mathrm{ADS}-B}\left(z_k\right)=\frac{{\omega }_{k|k-1}^j{q}_k^{j,\mathrm{ADS}-B}\left(z_k\right)}{{\kappa }_k^{\mathrm{ADS}-B}\left(z_k\right)+{\Sigma }_{t=1}^{J_{k|k-1}}{\omega }_{k|k-1}^t{q}_k^{t,\mathrm{ADS}-B}\left(z_k\right)}$

$m_{k|k}^{j,\mathrm{ADS}-B}\left(z_k\right)=m_{k|k-1}^j+K_{k,x}^{j,\mathrm{ADS}-B}(z_k-{\hat{z}}_{k|k-1}^{j,\mathrm{ADS}-B})$

$P_{k|k}^{j,\mathrm{ADS}-B}=(I-K_{k,x}^{j,\mathrm{ADS}-B}{H}_k^{\mathrm{ADS}-B})P_{k|k-1}^j$

其中,

${\hat{z}}_{k|k-1}^{j,\mathrm{ADS}-B}=H_k^{\mathrm{ADS}-B}{m}_{k|k-1}^j$

$q_k^{j,\mathrm{ADS}-B}\left(z_k\right)=N(z_k;{\hat{z}}_{k|k-1}^{j,\mathrm{ADS}-B},H_k^{\mathrm{ADS}-B}{P}_{k|k-1}^j{\left(H_k^{\mathrm{ADS}-B}\right)}^T+R_k^{\mathrm{ADS}-B})$

$K_{k,x}^{j,\mathrm{ADS}-B}=P_{k|k-1}^j{\left(H_k^{\mathrm{ADS}-B}\right)}^T{[H_k^{\mathrm{ADS}-B}{P}_{k|k-1}^j{\left(H_k^{\mathrm{ADS}-B}\right)}^T+R_k^{\mathrm{ADS}-B}]}^{-1}$

由于ADS-B测量中不含杂波，故

当使用ADS-B测量数据对预测强度函数更新后，目标状态高斯项的权值、均值及协 方差通过卡尔曼滤波器得到传递。然后使用雷达测量数据再次更新强度函数，由于雷达 系统误差和目标状态通过似然函数耦合在一起，采用两步卡尔曼滤波器可以依次更新目 标状态和雷达系统误差高斯项的权值、均值及协方差。

雷达测量数据更新后的后验强度函数为

${\upsilon }_{k|k}^{\mathrm{radar}}=(1-p_D^{\mathrm{radar}}){\upsilon }_{k|k}^{\mathrm{ADS}-B}+\underset{z_k\in Z_k^{\mathrm{radar}}}{\Sigma }{\upsilon }_{d,k}^{\mathrm{radar}}(x_k;z_k)---\left(17\right)$

其中

${\upsilon }_{d,k}^{\mathrm{radar}}(x_k;z_k)=\underset{j=1}{\overset{J_{k,\mathrm{ADS}-B}}{\Sigma }}{\omega }_k^{j,\mathrm{radar}}\left(z_k\right)N(x_k;m_{k|k}^{j,\mathrm{radar}}\left(z_k\right),P_{k|k}^{j,\mathrm{radar}})N(b_k;b_{k|k}^j,B_{k|k}^j)$

雷达测量数据更新后的目标状态高斯项的权值、均值及协方差分别为

${\omega }_k^{j,\mathrm{radar}}\left(z_k\right)=\frac{{p_D^{\mathrm{radar}}\omega }_k^{j,\mathrm{ADS}-B}{q}_k^{j,\mathrm{radar}}\left(z_k\right)}{{\kappa }_k^{\mathrm{radar}}\left(z_k\right)+p_D^{\mathrm{radar}}{\Sigma }_{t=1}^{J_{k,\mathrm{ADS}-B}}{\omega }_k^{j,\mathrm{ADS}-B}{q}_k^{j,\mathrm{radar}}\left(z_k\right)}$

$m_{k|k}^{j,\mathrm{radar}}\left(z_k\right)=m_{k|k}^{j,\mathrm{ADS}-B}+K_{k,x}^{j,\mathrm{radar}}(z_k-{\hat{z}}_{k|k-1}^{j,\mathrm{radar}})$

$P_{k|k}^{j,\mathrm{radar}}=(I-K_{k,x}^{j,\mathrm{radar}}{H}_k^{\mathrm{radar}})P_{k|k}^{j,\mathrm{ADS}-B}$

雷达系统误差高斯项的均值、协方差分别为

$b_{k|k}^{j,\mathrm{radar}}\left(z_k\right)=b_{k|k-1}^j+K_{k,b}^{j,\mathrm{radar}}(z_k-{\hat{z}}_{k|k-1}^{j,\mathrm{radar}})$

$B_{k|k}^j=(I-K_{k,b}^{j,\mathrm{radar}})B_{k|k-1}^j$

其中，

${\kappa }_k^{\mathrm{radar}}\left(z_k\right)=\mathrm{\lambda c}\left(z_k\right)$

${\hat{z}}_{k|k-1}^{j,\mathrm{radar}}=H_k^{\mathrm{radar}}{m}_{k|k}^{j,\mathrm{ADS}-B}+b_{k|k-1}^j$

$q_k^{j,\mathrm{radar}}\left(z_k\right)=N(z_k;{\hat{z}}_{k|k-1}^{j,\mathrm{radar}},H_k^{\mathrm{radar}}{P}_{k|k}^{j,\mathrm{ADS}-B}{\left(H_k^{\mathrm{radar}}\right)}^T+R_k^{\mathrm{radar}}+B_{k|k-1}^j)$

$K_{k,x}^{j,\mathrm{radar}}=P_{k|k}^{j,\mathrm{ADS}-B}{\left(H_k^{\mathrm{radar}}\right)}^T{[H_k^{\mathrm{radar}}{P}_{k|k}^{j,\mathrm{ADS}-B}{\left(H_k^{\mathrm{radar}}\right)}^T+R_k^{\mathrm{radar}}+B_{k|k-1}^j]}^{-1}$

$K_{k,b}^{j,\mathrm{radar}}=B_{k|k-1}^j{[H_k^{\mathrm{radar}}{P}_{k|k}^{j,\mathrm{ADS}-B}{\left(H_k^{\mathrm{radar}}\right)}^T+R_k^{\mathrm{radar}}+B_{k|k-1}^j]}^{-1}$

λ为杂波数，c(z_k)为杂波分布函数，这里假设服从观测空间内的均匀分布。

d、对更新获得的扩展状态的强度函数的高斯项进行裁剪合并，获得目标状态和 雷达系统误差估计，方法同GM-PHD算法相同。

e、雷达系统误差融合估计

雷达系统误差融合估计的方法是使用加权平均算法对获得的多个雷达系统误差估计 值进行加权融合。

$\hat{b}=\frac{\underset{i=1}{\overset{J_k}{\Sigma }}{\omega }_k^{i,\mathrm{radar}}\left(z_k\right)b_{k|k}^{i,\mathrm{radar}}\left(z_k\right)}{\underset{i=1}{\overset{J_k}{\Sigma }}{\omega }_k^{i,\mathrm{radar}}\left(z_k\right)}---\left(18\right)$

本发明提供的利用ADS-B数据的GM-EPHD滤波雷达系统误差联合估计方法的效果 可以通过以下仿真结果进一步说明。

仿真数据描述：实验设置在一个[-1000,1000]×[-1000,1000]的监视区域，3个目标运 动，一部二坐标雷达处于监视区域的中心，观测获得目标距离、方位信息。3个目标均配 备机载ADS-B收发信机，地面站接收广播信息获得目标经度、纬度及高度信息，ADS-B 观测值转换到以雷达为中心的直角坐标系下。仿真雷达及ADS-B同步采样，周期均为 T＝1s，仿真100步，雷达检测概率为P_D,k＝P_D＝0.99。

雷达系统误差真值为(50m,0.01rad)。

雷达随机测量噪声协方差矩阵$R_k^{\mathrm{radar}}=\mathrm{diag}\{{\left(20m\right)}^2,{\left(0.005\mathrm{rad}\right)}^2\}.$

雷达观测杂波数为λ＝10，杂波分布函数c(z_k)服从[-1000,1000]×[-1000,1000]上的 均匀分布。

系统加速度扰动标准差为σ_v＝0.1m/s²，系统误差噪声协方差矩阵 B_k＝diag{(5m)²,(0.001rad)²}。

高斯项修剪门限T＝10^-5，合并门限U＝4，高斯项最大个数J_max＝100。

实验运行环境为Intel Core2 Quad CPU 2.66GHz，2GB内存，仿真软件为Matlab  R2010a。

为了比较本发明方法的估计结果，选取文献(Wenling Li,Jia Y,Du J,Yu F.Gaussian  mixture PHD filter for multi-sensor multi-target tracking with registration errors[J]. Signal Processing,2012)中的GM-EPHD联合估计方法，利用相同数据进行100次蒙特 卡罗比较试验，文献方法中的“sensor1”和“sensor2”都存在系统误差，利用“sensor1”和 “sensor2”的融合观测对“sensor1”进行系统误差估计，而本发明中使用ADS-B数据代替文 献方法中的“sensor2”，使用ADS-B观测和“sensor1”观测联合估计目标状态和“sensor1”的 系统误差。

图2为实验所设置的目标真实运动轨迹及雷达和ADS-B观测数据，‘—*—’为雷 达观测值，‘—o—’为ADS-B观测值，‘——’为目标位置的真实值。其中目标1从k＝1s 时刻出发，k＝100s结束，其初始位置为(250m,250m)；目标2从k＝1s时刻出发，k＝100s 结束，其初始位置为(-250m,-250m)；目标3在k＝30s时从目标1中分开，直到k＝100s 结束。

图3为目标真实运动轨迹及状态估计结果，‘—o—’为目标状态估计值，‘——’ 为目标位置的真实值；图4为目标状态估计的OSPA距离，实验中设置p＝2,c＝200；图 5为目标数估计结果，‘—△—’为目标数估计值，‘—o—’为目标数真实值。从图3、图 4可以看出，使用本发明方法进行目标状态估计的平均误差近似为5m，估计出的目标状 态和真实状态接近。从图5可以看出，采用本发明方法对目标数目估计结果和实际情况 相吻合。

图6、图7分别为采用两种方法获得的雷达距离、方位误差估计比较图，‘—*—’为 采用本发明方法估计结果，‘—△—’为采用文献方法的估计结果，‘——’为雷达系统误 差真值，对比估计结果可以看出，采用本发明方法估计稳定后雷达系统误差估计结果更 接近系统误差真值。图8、图9分别为采用两种方法获得的雷达距离、方位误差估计的均 方根误差比较图，‘—*—’为采用本发明方法估计的均方根误差，‘—△—’为采用文献方 法的均方根误差，对比结果可以看出，本发明方法获得的雷达系统误差估计的均方根误 差要远远小于文献方法。

图10为采用本发明方法对雷达系统误差校正后的雷达观测，‘—*—’为雷达系统误 差校正前的观测值，‘—o—’为采用本发明方法对雷达系统误差校正后的观测值，‘——’ 为目标状态真实值，对比校正前后的雷达观测值可以看出，使用本发明方法对雷达系统 误差进行校正后，雷达观测与真实运动轨迹已经非常接近了，证明本发明方法对雷达系 统误差估计是准确有效的。