CN104237862A - 基于ads-b的概率假设密度滤波雷达系统误差融合估计方法 - Google Patents

Abstract

一种基于ADS-B的概率假设密度滤波雷达系统误差融合估计方法。其包括建立基于ADS-B的雷达系统误差观测方程；建立雷达系统误差状态方程；构造雷达系统误差状态和观测有限集；概率假设密度滤波雷达系统误差融合估计等阶段。本发明提供的估计方法首先将目标的ADS-B和雷达的观测转换到以雷达为中心的直角坐标系下，由于ADS-B监视精度远远优于雷达定位精度，在不考虑ADS-B定位误差的基础上，建立基于ADS-B的雷达系统误差测量方程和状态方程，然后使用概率假设密度滤波器对坐标统一后的ADS-B与雷达的测量差值进行融合估计，进而得到雷达系统误差的估计结果。本发明方法具有估计精度高、估计性能好等优点。

Description

基于ADS-B的概率假设密度滤波雷达系统误差融合估计方法

技术领域

本发明属于传感器误差配准技术领域，特别是涉及一种基于ADS-B的概率假设密度滤波雷达系统误差融合估计方法。

背景技术

雷达观测误差主要分为随机误差和系统误差两类。在多雷达融合跟踪系统中，雷达系统误差估计已成为多雷达融合处理的先决条件，其会直接影响整个系统的工作性能。因此需要对雷达系统误差进行估计，以此对雷达测量进行相应的补偿，这个过程也称为误差配准。现有的雷达误差估计方法可以归纳为离线估计方法和在线估计方法两类。其中离线估计方法是通过对一段时间的雷达观测数据进行数据拟合，从而估计系统误差，如最小二乘方法及最大似然方法等。在线估计方法主要是利用滤波方法实现系统误差的递推估计。在线估计方法与离线方法相比具有实现实时误差校正的优点，因此得到的国内外学者的更多关注。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滤波器的封闭解形式。2009年Vo利用类似扩展Kalman滤波的方法给出了了弱非线性条件下的GM-PHD滤波。

随着全球导航卫星系统(GNSS)和空空、空地数据链通信技术的发展，一种新型的航空器运行监视技术—广播式自动相关监视技术(Automatic DependentSurveillance-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的概率假设密度滤波雷达系统误差融合估计方法。

为了达到上述目的，本发明提供的基于ADS-B的概率假设密度滤波雷达系统误差融合估计方法包括按顺序进行的下列步骤：

1)建立基于ADS-B的雷达系统误差观测方程的S1阶段；

2)建立雷达系统误差状态方程的S2阶段；

3)构造雷达系统误差状态和观测有限集的S3阶段；

4)概率假设密度滤波雷达系统误差融合估计的S4阶段。

在步骤1)中，所述的建立基于ADS-B的雷达系统误差观测方程的方法是利用坐标投影技术首先将目标ADS-B观测值转换到以雷达位置为中心的直角坐标系下，然后将目标雷达极坐标的观测值转换到以雷达站位置为中心的直角坐标系下，在不考虑ADS-B定位误差的基础上，建立基于ADS-B的雷达系统误差观测方程，为概率假设密度滤波雷达系统误差融合估计做准备。

在步骤2)中，所述的建立雷达系统误差状态方程的方法是将雷达的系统误差和目标的状态看作两个相互独立的过程，并且假设雷达的系统误差服从马尔可夫过程，进而建立雷达系统误差状态转移方程，为概率假设密度滤波雷达系统误差融合估计做准备。

在步骤3)中，所述的构造雷达系统误差状态和观测有限集的方法是认为每一个目标均可获取一组雷达系统误差估计值，然后将所有目标的雷达系统误差估计值组成雷达系统误差状态有限集，且将坐标统一后的任意的雷达观测值与ADS-B观测值的差值作为雷达系统误差观测值，以此组成雷达系统误差观测有限集。

在步骤4)中，所述的概率假设密度滤波雷达系统误差融合估计的方法是使用概率假设密度滤波器对雷达系统误差进行融合估计。

本发明提供的基于ADS-B的概率假设密度滤波雷达系统误差融合估计方法首先将目标的ADS-B和雷达的观测转换到以雷达为中心的直角坐标系下，由于ADS-B监视精度远远优于雷达定位精度，在不考虑ADS-B定位误差的基础上，建立基于ADS-B的雷达系统误差测量方程和状态方程，然后使用概率假设密度滤波器对坐标统一后的ADS-B与雷达的测量差值进行融合估计，进而得到雷达系统误差的估计结果。本发明方法具有估计精度高、估计性能好等优点。

附图说明

图1为本发明提供的基于ADS-B的概率假设密度滤波雷达系统误差融合估计方法流程图。

图2为目标运动轨迹及观测值。

图3为雷达斜距估计随时间变化曲线图。

图4为雷达方位角估计随时间变化曲线图。

图5为雷达俯仰角估计随时间变化曲线图。

图6为雷达斜距估计均方根误差随时间变化曲线图。

图7为雷达方位角估计均方根误差随时间变化曲线图。

图8为雷达俯仰角估计均方根误差随时间变化曲线图。

图9为雷达系统误差校正前后的雷达观测对比图。

具体实施方式

下面结合附图和具体实施例对本发明提供的基于ADS-B的概率假设密度滤波雷达系统误差融合估计方法进行详细说明。

图1为本发明提供的基于ADS-B的概率假设密度滤波雷达系统误差融合估计方法流程图。其中的全部操作都是在计算机系统中完成的，操作的主体均为计算机系统。

如图1所示，本发明提供的基于ADS-B的概率假设密度滤波雷达系统误差融合估计方法包括按顺序进行的下列步骤：

1)建立基于ADS-B的雷达系统误差观测方程的S1阶段：

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

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

a、将目标ADS-B观测值转换到以雷达位置为中心的直角坐标系下

首先将目标的ADS-B观测值转换到地心地固(ECEF)坐标系下，然后再转换为以雷达站位置为中心的坐标下的观测值。设ADS-B系统获得的k时刻第i个目标的观测为 $(L_{k,\mathrm{ADS}-B}^{\left(i\right)},{\mathrm{\λ}}_{k,\mathrm{ADS}-B}^{\left(i\right)},H_{k,\mathrm{ADS}-B}^{\left(i\right)}),$ $L_{k,\mathrm{ADS}-B}^{\left(i\right)},{\mathrm{\λ}}_{k,\mathrm{ADS}-B}^{\left(i\right)},H_{k,\mathrm{ADS}-B}^{\left(i\right)}$ 分别表示目标的经度、纬度及距离海平面的高度。则转换后的ECEF坐标系系观测值为 $({\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)}),$ 则：

$\left\{\begin{array}{c}{\mathrm{Ex}}_{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)}\cos {\mathrm{\λ}}_{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 {\mathrm{\λ}}_{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)$

其中， $C=\frac{E_q}{{(1-e^2{\sin }^2{L}_{k,\mathrm{ADS}-B}^{\left(i\right)})}^{\frac{1}{2}}},$ e为地球的偏心率，E_q为赤道半径。

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

将ECEF坐标系下的目标ADS-B观测值 $({\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)})$ 通过式(2)转换为以雷达站位置为中心的直角坐标下的观测值 $({\mathrm{Cx}}_{k,\mathrm{ADS}-B}^{\left(i\right)},{\mathrm{Cy}}_{k,\mathrm{ADS}-B}^{\left(i\right)},{\mathrm{Cz}}_{k,\mathrm{ADS}-B}^{\left(i\right)}):$

$\left[\begin{array}{c}{\mathrm{Cx}}_{k,\mathrm{ADS}-B}^{\left(i\right)}\\ {\mathrm{Cy}}_{k,\mathrm{ADS}-B}^{\left(i\right)}\\ {\mathrm{Cz}}_{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 {\mathrm{\λ}}_{\mathrm{RS}}\cos L_{\mathrm{RS}}& -\sin {\mathrm{\λ}}_{\mathrm{RS}}\sin L_{\mathrm{RS}}& \cos {\mathrm{\λ}}_{\mathrm{RS}}\\ \sin L_{\mathrm{RS}}& \cos L_{\mathrm{RS}}& 0\\ \cos {\mathrm{\λ}}_{\mathrm{RS}}\cos L_{\mathrm{RS}}& \cos {\mathrm{\λ}}_{\mathrm{RS}}\sin L_{\mathrm{RS}}& \sin {\mathrm{\λ}}_{\mathrm{RS}}\end{array}\right)$

b、将目标雷达极坐标的观测值转换到以雷达站位置为中心的直角坐标系下

雷达观测数据使用以雷达位置为中心的极坐标系获得目标的斜距、方位及俯仰信息。雷达的观测模型可以用式(3)表示：

其中，分别表示k时刻第i个目标的斜距、方位及俯仰的观测值， 分别表示斜距、方位及俯仰对应的真实值，分别表示斜距、方位及俯仰的系统误差，分别表示斜距、方位及俯仰的观测噪声。

设雷达k时刻获得的第i个目标的斜距、方位及俯仰分别为转换后以雷达站位置为中心的直角坐标系下的目标观测值为

以雷达斜距、方位及俯仰的系统误差作为状态向量，以雷达站位置为中心的直角坐标系下的雷达观测值与ADS-B观测值的差值 $z={[\widetilde{{\mathrm{Cx}}_{k,\mathrm{RADAR}}^{\left(i\right)}}-{\mathrm{Cx}}_{k,\mathrm{ADS}-B}^{\left(i\right)},\widetilde{{\mathrm{Cy}}_{k,\mathrm{RADAR}}^{\left(i\right)}}-{\mathrm{Cy}}_{k,\mathrm{ADS}-B}^{\left(i\right)},\widetilde{{\mathrm{Cz}}_{k,\mathrm{RADAR}}^{\left(i\right)}}-{\mathrm{Cz}}_{k,\mathrm{ADS}-B}^{\left(i\right)}]}^T$ 作为雷达系统误差的观测量，构建观测方程，将式(3)误差模型转换到以雷达站位置为中心的直角坐标系下，可以表示为：

由于ADS-B监视精度远远优于雷达定位精度，在不考虑ADS-B系统误差的情况下，假设以ADS-B观测值作为目标位置的真值，即则有：

其中，同理有：

从而得到雷达系统误差的观测方程如下：

其中，

表示斜距、方位及俯仰的观测噪声，且n_k～N(·,0,R_k)。

2)建立雷达系统误差状态方程的S2阶段：

本阶段是将雷达的系统误差独立于目标，并且服从马尔可夫过程，从而建立雷达系统误差状态转移方程，然后进入下一步S3阶段。

在此阶段中，所述的建立雷达系统误差状态方程的具体方法如下。假设雷达系统误差独立于目标，且服从马尔可夫过程，则雷达系统误差状态转移方程可以表示为：

其中， $F_{k-1}=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right],$ 表示系统噪声，且w_k-1～N(·,0,Q_k-1)。

3)构造雷达系统误差状态和观测有限集的S3阶段：

将雷达的系统误差估计和目标的状态估计看作两个相互独立的过程，并且假设雷达的系统误差估计服从马尔可夫过程，进而建立雷达系统误差状态转移方程，然后进入下一步S4阶段。

在此阶段中，所述的构造雷达系统误差状态和观测有限集的具体方法如下。假设监视区域内存在N_k个目标，均配备ADS-B机载设备。其中M_k个目标被雷达观测到(由于雷达有一定的检测概率，一般情况下N_k≥M_k)。

a、构造雷达系统误差状态有限集合

以雷达斜距、方位及俯仰的系统误差构造误差状态量。对于M_k个被雷达观测到的目标，每一个目标均可获得一组雷达系统误差估计值，因此可以构造一个含有M_k个元素的系统误差状态集合，即：

$X_k=\{x_k^{\left(1\right)},...,x_k^{\left(M_k\right)}\}\∈F\left(\mathrm{\χ}\right)---\left(9\right)$

其中，F(χ)表示系统误差状态的有限集合。

b、构造雷达系统误差观测有限集合

假设k时刻在监视区域内存在N_k个配备ADS-B机载设备的目标，则可以获得N_k个目标个ADS-B观测，即：

$Z_{k,\mathrm{ADS}-B}=\{z_{k,\mathrm{ADS}-B}^{\left(1\right)},z_{k,\mathrm{ADS}-B}^{\left(2\right)}...,z_{k,\mathrm{ADS}-B}^{\left(N_k\right)}\}---\left(10\right)$

其中， $z_{k,\mathrm{ADS}-B}^{\left(i\right)}={[{\mathrm{Cx}}_{k,\mathrm{ADS}-B}^{\left(i\right)},{\mathrm{Cy}}_{k,\mathrm{ADS}-B}^{\left(i\right)},{\mathrm{Cz}}_{k,\mathrm{ADS}-B}^{\left(i\right)}]}^T.$

同时利用雷达对目标进行观测，获得M_k个目标的雷达观测值，即：

$Z_{k,\mathrm{RADAR}}=\{z_{k,\mathrm{RADAR}}^{\left(1\right)},z_{k,\mathrm{RADAR}}^{\left(2\right)}...,z_{k,\mathrm{RADAR}}^{\left(M_k\right)}\}---\left(11\right)$

其中， $z_{k,\mathrm{RADAR}}^{\left(i\right)}={[\widetilde{{\mathrm{Cx}}_{k,\mathrm{RADAR}}^{\left(i\right)}},\widetilde{{\mathrm{Cy}}_{k,\mathrm{RADAR}}^{\left(i\right)}},\widetilde{{\mathrm{Cz}}_{k,\mathrm{RADAR}}^{\left(i\right)}}]}^T.$

将M_k个雷达观测值与N_k个ADS-B观测值任意两者之间的差值作为雷达系统误差的观测值，构成M_k×N_k个误差观测值集合Z_k：

$Z_k=\{z_k^{\left(1\right)},z_k^{\left(2\right)}...,z_k^{(N_k\×M_k)}\}\∈F\left(Z\right)---\left(12\right)$

其中，F(Ζ)表示系统误差观测的有限集合。

对于雷达系统误差来说，在M_k×N_k个误差观测值中，只有M_k个来自同一目标的雷达观测值与ADS-B观测值的差值才能作为正确的误差观测来更新雷达系统误差的估计，而其余的M_k×(N_k-1)则应将其视为杂波观测。这样避免了传统滤波方法需利用数据关联技术获得同一目标的雷达观测与ADS-B观测对的过程。

4)概率假设密度滤波雷达系统误差融合估计的S4阶段：

本阶段使用概率假设密度滤波器对雷达系统误差进行融合估计。

在此阶段中，所述的概率假设密度滤波雷达系统误差融合估计的具体方法如下。

a、以雷达斜距、方位及俯仰的系统误差构造误差状态向量，假设k-1时刻的误差强度函数为υ_k-1(x)：

${\mathrm{\υ}}_{k-1}\left(x\right)=\underset{i=1}{\overset{J_{k-1}}{\mathrm{\Σ}}}{\mathrm{\ω}}_{k-1}^{\left(i\right)}N(x;m_{k-1}^{\left(i\right)},P_{k-1}^{\left(i\right)})---\left(13\right)$

b、预测步：

预测目标强度函数得到υ_k|k-1(x)，由于假设雷达系统误差始终存在且不变，因此不存在新生状态，且存活概率p_S,k＝1，即：

${\mathrm{\υ}}_{k|k-1}\left(x\right)={\mathrm{\υ}}_{S,k|k-1}\left(x\right)=\underset{i=1}{\overset{J_{k-1}}{\mathrm{\Σ}}}{\mathrm{\ω}}_{k-1}^{\left(j\right)}N(x;m_{S,k|k-1}^{\left(i\right)},P_{S,k|k-1}^{\left(j\right)})---\left(14\right)$

其中， $\begin{array}{c}{m}_{S,k|k-1}^{\left(j\right)}=F_{k-1}{m}_{|k-1}^{\left(j\right)}\\ P_{S,k|k-1}^{\left(j\right)}=Q_{k-1}+F_{k-1}{P}_{k-1}^{\left(j\right)}{F}_{k-1}^T\end{array}$

c、更新步：

以雷达观测值与ADS-B观测值的差值构建观测有限集合进行雷达系统误差估计更新：

$v_{k|k}\left(x\right)=(1-P_{D,k})v_{k|k-1}\left(x\right)+\underset{z\∈Z_k}{\mathrm{\Σ}}{v}_{D,k}(x;z)---\left(15\right)$

其中，

$v_{D,k}(x;z)=\underset{j=1}{\overset{J_{k|k-1}}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^{\left(i\right)}\left(z\right)N(x;m_{k|k}^{\left(j\right)}\left(z\right),P_{k|k}^{\left(j\right)})$

${\mathrm{\ω}}_k^{\left(i\right)}\left(z\right)=\frac{p_{D,k}{\mathrm{\ω}}_{k|k-1}^{\left(j\right)}{q}_k^{\left(j\right)}\left(z\right)}{{\mathrm{\κ}}_k\left(z\right)+p_{D,k}\underset{l=1}{\overset{J_{k|k-1}}{\mathrm{\Σ}}}{\mathrm{\ω}}_{k|k-1}^{\left(l\right)}{q}_k^{\left(l\right)}\left(z\right)}$

$\begin{array}{c}{\mathrm{\κ}}_k\left(z\right)=\mathrm{\λc}\left(z\right)\\ q_k^{\left(j\right)}\left(z\right)=N(z;H_k{m}_{k|k-1}^{\left(j\right)},H_k{P}_{k|k-1}^{\left(j\right)}{H}_k^T+R_k)\\ m_{k|k}^{\left(j\right)}\left(z\right)=m_{k|k-1}^{\left(j\right)}+K_k^{\left(j\right)}(z-H_k{m}_{k|k-1}^{\left(j\right)})\\ P_{k|k}^{\left(j\right)}=[I-K_k^{\left(j\right)}{H}_k]P_{k|k-1}^{\left(j\right)}\\ K_k^{\left(j\right)}=P_{k|k-1}^{\left(j\right)}{H}_k^T{(H_k{P}_{k|k-1}^{\left(j\right)}{H}_k^T+R_k)}^{-1}\end{array}$

λ为杂波数，可由下式估计获得：

$\widehat{\mathrm{\λ}}=(N_k-1)\×M_k$

其中，N_k为k时刻获得的目标ADS-B观测数，M_k为k时刻获得的目标雷达观测数。c(z)为杂波分布函数，这里假设服从观测空间的均匀分布。

d、对更新获得的系统误差状态的强度函数υ_k|k(x)的高斯项进行裁剪合并获得目标状态估计，方法与GM-PHD算法相同。

e、对获得的多个雷达系统误差估计值进行加权融合：

$\hat{x}=\frac{\underset{i=1}{\overset{J_k}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^{\left(i\right)}\left(z\right)x_k^{\left(i\right)}}{\underset{i=1}{\overset{J_k}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^{\left(i\right)}\left(z\right)}---\left(16\right)$

本发明提供的基于ADS-B的概率假设密度滤波雷达系统误差融合估计方法的效果可以通过以下仿真结果进一步说明。

仿真数据描述：实验设置在一个[-10000,10000]×[-10000,10000]的监视区域，5个目标运动(如图2所示)，采样周期T＝1s，仿真100步，系统误差状态向量为其中分别表示斜距、方位、俯仰对应的误差值，雷达检测概率为P_D,k＝P_D＝0.95，误差存活概率为P_S,k＝P_S＝1。观测噪声n_k～N(.；0,R_k)，其中，1σ_n＝[50,0.005,0.005]^T。系统噪声w_k～N(.；0,Q_k)，其中，方差σ_w＝[20,0.002,0.002]^T。高斯项修剪阈值为T＝10^-5，合并门限U＝4，高斯项最大个数J_max＝100，实验运行环境为Intel Core2 Quad CPU 2.66GHz，2GB内存，仿真软件为MatlabR2010a。

图2为5个目标的运动轨迹及观测，其中，‘——’为目标真实运动轨迹，‘—*—’为目标ADS-B观测，‘—△—’为目标雷达观测。

图3为雷达斜距估计随时间变化曲线图，图4为雷达俯仰角估计随时间变化曲线图，图5为雷达方位角估计随时间变化曲线图。其中，‘——’为误差真值，‘—*—’为采用本发明方法的估计结果。

图6为雷达斜距估计均方根误差随时间变化曲线图，图7为雷达俯仰角估计均方根误差随时间变化曲线图，图8为雷达方位角估计均方根误差随时间变化曲线图。其中，‘—*—’为采用本发明方法估计的均方根误差。

图9为雷达系统误差校正前后的雷达观测对比图，从图中可以看出，误差校正后的雷达观测接近目标真实的运动轨迹。

Claims (5)

1.一种基于ADS-B的概率假设密度滤波雷达系统误差融合估计方法，其特征在于：所述的一种基于ADS-B的概率假设密度滤波雷达系统误差融合估计方法包括按顺序进行的下列步骤：

1)建立基于ADS-B的雷达系统误差观测方程的S1阶段；

2)建立雷达系统误差状态方程的S2阶段；

3)构造雷达系统误差状态和观测有限集的S3阶段；

4)概率假设密度滤波雷达系统误差融合估计的S4阶段。

2.根据权利要求1所述的基于ADS-B的概率假设密度滤波雷达系统误差融合估计方法，其特征在于：在步骤1)中，所述的建立基于ADS-B的雷达系统误差观测方程的方法是利用坐标投影技术首先将目标ADS-B观测值转换到以雷达位置为中心的直角坐标系下，然后将目标雷达极坐标的观测值转换到以雷达站位置为中心的直角坐标系下，在不考虑ADS-B定位误差的基础上，建立基于ADS-B的雷达系统误差观测方程，为概率假设密度滤波雷达系统误差融合估计做准备。

3.根据权利要求1所述的基于ADS-B的概率假设密度滤波雷达系统误差融合估计方法，其特征在于：在步骤2)中，所述的建立雷达系统误差状态方程的方法是将雷达的系统误差和目标的状态看作两个相互独立的过程，并且假设雷达的系统误差服从马尔可夫过程，进而建立雷达系统误差状态转移方程，为概率假设密度滤波雷达系统误差融合估计做准备。

4.根据权利要求1所述的基于ADS-B的概率假设密度滤波雷达系统误差融合估计方法，其特征在于：在步骤3)中，所述的构造雷达系统误差状态和观测有限集的方法是认为每一个目标均可获取一组雷达系统误差估计值，然后将所有目标的雷达系统误差估计值组成雷达系统误差状态有限集，且将坐标统一后的任意的雷达观测值与ADS-B观测值的差值作为雷达系统误差观测值，以此组成雷达系统误差观测有限集。

5.根据权利要求1所述的基于ADS-B的概率假设密度滤波雷达系统误差融合估计方法，其特征在于：在步骤4)中，所述的概率假设密度滤波雷达系统误差融合估计的方法是使用概率假设密度滤波器对雷达系统误差进行融合估计。