CN103941233A - 一种多平台主被动传感器协同跟踪的雷达间歇交替辐射控制方法 - Google Patents

Abstract

本发明公开了一种多平台主被动传感器协同跟踪的雷达辐射控制方法，对多个搭载主被动传感器的平台进行雷达辐射协同管控，在满足跟踪任务需求的条件下，控制雷达在时间上间歇开机、在空间上交替开机，并利用基于坐标交替转换的多平台主被动序贯滤波实现多平台辐射管控下的融合跟踪问题，该跟踪方法不需要多平台被动传感器量测时间同步，也不需要多平台被动传感器的交叉定位，能够适应被动量测不连续和无规律出现的情况，具有良好的适应性，适合工程应用，解决雷达长时间开机辐射易被侦察的问题，进一步提高雷达的隐蔽能力，进而提高多平台系统整体生存能力。

Description

一种多平台主被动传感器协同跟踪的雷达间歇交替辐射控制方法

一、技术领域

本发明隶属于传感器管理领域，适用于解决多平台主被动传感器协同跟踪的雷达辐射控制问题。 

二、背景技术

现代复杂电磁环境下，各种侦察和干扰设备日趋完善，为了尽可能提高作战平台的隐蔽性，防止已方雷达被对方侦察到进而受到干扰和攻击，就必须在满足跟踪精度的前提下尽可能减少雷达的电磁辐射，提高作战平台的隐蔽性和反侦察能力。 

现有的雷达辐射控制以及辐射控制下的主被动传感器融合跟踪技术主要集中在单平台，如单机载平台雷达、红外、电子支援措施的协同跟踪与雷达开关机管理等，这些方法只对单平台上的雷达在时间上进行间歇开机管理，而多平台环境下，不同平台被动传感器只能提供目标方位角、俯仰角，不能提供目标距离，无法将被动量测点坐标转换到同一平台坐标系下进行滤波跟踪，而双平台交叉定位方法虽然在一定程度上解决目标距离未知的问题，但一般需要被动量测时间同步，难以解决被动量测不同步或间歇出现的情况。可见，现有雷达辐射控制方法不能适用于多平台辐射管理下多平台主被动协同跟踪。 

三、发明内容

本发明所要解决的问题就是，针对现有方法不适用多平台协同雷达辐射管控的问题，提供一种多平台主被动传感器协同跟踪的雷达辐射控制方法，对多个搭载主被动传感器的平台进行雷达辐射协同管控，在满足跟踪任务需求的条件下，控制雷达时间上间歇开机，空间上交替开机，并实现了多平台雷达间歇交替辐射时的主被动融合跟踪，提高多平台探测系统的隐蔽性。 

本发明解决所述技术问题，采用技术方案步骤如下： 

(1)被动传感器引导雷达开机完成目标航迹起始； 

(2)利用协方差控制方法在时间上确定下一次需要雷达辐射的时刻； 

(3a)当预定开机时间到，没有得到新的被动量测，则跳到第(4a)步； 

(3b)当预定开机时间未到，有被动量测时，首先用所得到的被动量测点迹与融合航迹进行关联，然后将关联上的被动量测利用基于坐标交替转换的多平台主被动序贯滤波进行跟踪，回到第(2)步； 

(4a)在空间上选择不同平台交替开机辐射； 

(4b)雷达开机得到新的雷达量测后，先利用所得雷达量测点迹与融合航迹进行关联，然后将关联上的雷达量测利用基于坐标交替转换的多平台主被动序贯滤波算法进行跟踪，跳回到第(2)步。 

具体的，所述步骤(2)具体为： 

a1)假设雷达最小和最大辐射间隔分别用T_min、T_max表示，并令T_left=T_min，T_right=T_max，进入下一步； 

a2)如果T_right-T_left＜ε(ε为一小的常数)，则令下次雷达开机辐射间隔的时间T₁=T_left或T₁=T_right，跳到第a5)步，否则令T_test＝T_left+[(T_right-T_left)/2]； 

a3)令t₁=T_test，计算间隔时间为t₁时扩展卡尔曼滤波的预测协方差P(t₁|0) 

P(t₁|0)＝F(t₁)·P(0|0)·F′(t₁)+Γ(t₁)·σ_v ²·Γ′(t₁) 

其中，P(0|0)表示最近一次滤波协方差估计，F(t₁)表示状态转移矩阵，Γ(t₁)表示过程噪声分布矩阵，σ_v ²为零均值高斯白噪声的协方差； 

a4)判断预测协方差是否超过门限。假设跟踪精度门限为R_th，如果  $\sqrt{P_{\left(t_1\right|0)}\left(\mathrm{1,1}\right)+P_{\left(t_1\right|0)}\left(\mathrm{3,3}\right)+P_{\left(t_1\right|0)}\left(\mathrm{5,5}\right)}<R_{\mathrm{th}}$ (其中分别是预测协方差矩阵P(t₁|0)的对角线上X、Y、Z位置误差方差项)，则T_left=T_test，否则T_right=T_lest，跳回到第a2)步； 

a5)得到下次预定雷达开机辐射间隔时间为T₁。 

具体的，所述步骤(3b)具体为： 

b1)得到与融合航迹关联的被动量测后，首先将该被动量测到来前最近的融合航迹点 以及对应的协方差P(k|k)转换到获得该量测的传感器所在平台的载体坐标系下，得到对应的航迹点P₁ ^NED(k|k)； 

b2)利用P₁ ^NED(k|k)进行外推，得到状态预测 

${\hat{X}}_1^{\mathrm{NED}}(k+1|k)=F\&CenterDot;{\hat{X}}_1^{\mathrm{NED}}\left(k\right|k)$

预测协方差 

P₁ ^NED(k+1|k)＝F·P₁ ^NED(k|k)·F^T

b3)然后将载体坐标系下的预测航迹和对应的协方差P₁ ^NED(k+1|k)从X-Y-Z直角坐标下转换到ρ-θ-ε(距离-方位-俯仰)的极坐标系 

${X_i}^{\mathrm{SP}}(k+1|k)={U_k}^T{\hat{X}}_1^{\mathrm{NED}}{(k+1|k)}_{|(x,y,z)}{U}_k$

P_i ^SP(k+1|k)＝U_k ^TP₁ ^NED(k+1|k)_|(x，y，z)U_k

式中，U_k为直角坐标到球极坐标的转换矩阵。 

构造关联统计量 

g₁(k+1)＝d₁(k)^T[P₁+R₁]^-1d₁(k) 

其中，P₁为P_i ^SP(k|k+1)在方位和俯仰上的误差协方差矩阵，R₁为被动量测在方位俯仰上的量测误差矩阵，d₁(k)为在方位和俯仰上的新息。 

利用如下公式进行量测与航迹的关联判决 

g₁(k+1)＜th    关联 

g₁(k+1)＞th    不关联 

其中th是二自由度χ²分布的统计量g₁(k)的假设检验门限。 

b4)得到与融合航迹关联的被动量测Z(k+1)=[θ_PAS(k+1)，ε_PAS(k+1)]′后，先利用 计算被动量测方位和俯仰的预测值 

$h(k+1|k)={[{\mathrm{\&theta;}}_{\mathrm{PAS}}(k+1|k),{\mathrm{\&epsiv;}}_{\mathrm{PAS}}(k+1|k)]}^{\&prime;}={\left[\begin{array}{c}\arctan 2[-z(k+1|k),x(k+1|k)]\\ \arctan 2[y\left(k\right),\sqrt{x^2(k+1|k)+z^2(k+1|k)}]\end{array}\right]}^{\&prime;}$

其中坐标系以载体为坐标原点的北-天-东为X-Y-Z方向(以下均以此为坐标系)，arctan2表示matlab软件中对应的反正切函数； 

利用计算被动量测时扩展卡尔曼滤波的量测矩阵 

$H_{\mathrm{PAS}}(k+1)=\left[\begin{array}{ccccccccc}-\frac{\hat{z}(k+1|k)}{{\hat{r}}^2}& 0& 0& 0& -\frac{\hat{x}(k+1|k)}{{\hat{r}}^2}& 0& 0& 0& 0\\ -\frac{\hat{x}(k+1|k)\hat{z}(k+1|k)}{\hat{r}{\hat{R}}^2}& 0& \frac{\hat{r}}{{\hat{R}}^2}& 0& -\frac{\hat{y}(k+1|k)\hat{z}(k+1|k)}{\hat{r}{\hat{R}}^2}& 0& 0& 0& 0\end{array}\right]$

其中，H_PAS(k+1)对应的状态向量表示形式为 $\left[\begin{array}{ccccccccc}x& \stackrel{.}{x}& y& \stackrel{.}{y}& z& \stackrel{.}{z}& \stackrel{..}{x}& \stackrel{..}{y}& \stackrel{..}{z}\end{array}\right],$ $\hat{r}=\sqrt{{\left[\hat{x}\left(k\right|k-1)\right]}^2+{\left[\hat{z}\left(k\right|k-1)\right]}^2},\hat{R}=\sqrt{{\hat{r}}^2+{\left[\hat{y}\left(k\right|k-1)\right]}^2}.$

利用雷达测量精度性能指标，得到量测误差矩阵 

$R(k+1)=\left[\begin{array}{cc}{\widehat{\mathrm{\&sigma;}}}_{{\mathrm{\&theta;}}_{\mathrm{PAS}}}^2& 0\\ 0& {\widehat{\mathrm{\&sigma;}}}_{{\mathrm{\&epsiv;}}_{\mathrm{PAS}}}^2\end{array}\right]$

其中，分别表示被动量测的方位和俯仰的误差方差。 

将被动量测Z(k+1)=[θ_PAS(k+1)，ε_PAS(k+1)]′、被动量测时扩展卡尔曼滤波的量测矩阵H_PAS、量测误差矩阵R(k+1)、被动量测方位和俯仰的预测值h(k+1|k)分别带入扩展卡尔曼滤波公式： 

$K(k+1|k)={P_1}^{\mathrm{NED}}(k+1)\&CenterDot;H_{\mathrm{PAS}}^T(k+1)\&CenterDot;{[H_{\mathrm{PAS}}(k+1)\&CenterDot;{P_1}^{\mathrm{NED}}(k+1|k\left)\right)\&CenterDot;H_{\mathrm{PAS}}^T(k+1)+R(k+1)]}^T$

P₁ ^NED(k+1|k+1)＝[I-K(k+1)H_PAS(k+1)]P₁ ^NED(k+1|k) 

${\hat{X}}_1^{\mathrm{NED}}(k+1|k+1)={\hat{X}}_1^{\mathrm{NED}}(k+1|k)+K(k+1)\&CenterDot;[Z(k+1)-h(k+1|k)]$

其中，I表示与状态向量维数相同的单位矩阵。 

b5)然后将和P₁ ^NED(k+1|k+1)分别转换回原融合航迹的坐标系下，成为 原航迹的下一个滤波状态估计点及对应的协方差。 

本发明的有益效果是，可以解决现代战场复杂电场环境下雷达长时间开机辐射易被侦察的问题，控制雷达在时间上间歇开机，在空间上交替开机，提高雷达的隐蔽能力，进一步提高多平台系统整体生存能力。本发明利用基于坐标交替转换的多平台主被动序贯滤波实现多平台辐射管控下的融合跟踪问题，该方法不需要多平台被动传感器量测时间同步，也不需要多平台被动传感器的交叉定位，能够适应被动量测不连续和无规律出现的情况，适应性强，适合工程应用。 

四、附图说明

附图1本发明的方法步骤流程图； 

附图2是用本发明进行仿真实验时己方飞机和目标飞机的飞行轨迹图 

附图3是用本发明进行仿真实验时被动传感器量测数据率为1s一个时2架己方飞机雷达辐射情况 

附图4是用本发明进行仿真实验时被动传感器量测数据率为1s一个时目标跟踪误差 

附图5是用本发明进行仿真实验时被动传感器量测数据率为10s一个时雷达辐射情况 

附图6是用本发明进行仿真实验时被动传感器量测数据率为10s一个时跟踪误差 

附图7是用现有单平台方法进行仿真实验时被动传感器量测数据率为1s一个时己方飞机1雷达辐射情况 

附图8是用现有单平台方法进行仿真实验时被动传感器量测数据率为1s一个时目标跟踪误差 

附图9是用现有单平台方法进行仿真实验时被动传感器量测数据率为10s一个时己方飞机1雷达辐射情况 

附图10是用现有单平台方法进行仿真实验时被动传感器量测数据率为10s一个时目标跟踪误差 

五、具体实施方法 

下面参照附图对本发明创造做进一步详细说明 

参照图1，本发明的具体步骤如下 

步骤一：首先使多个平台的被动传感器一直处于探测状态，并利用被动传感器引导两个平台中的某一雷达进行开机探测，得到雷达测量后，将雷达量测利用逻辑法进行航迹起始，利用三点法完成滤波初始化，得到初始化的状态估计和协方差P(0|0)，并通过坐标换转，将状态估计和协方差转换到融合中心坐标系下，得到融合中心的初始状态和初 始协方差P₁ ^NED(0|0)； 

步骤二：根据目标跟踪精度的任务需求，预测下一次雷达开机的时间，具体步骤为： 

a1)假设雷达最小和最大辐射间隔分别用T_min、T_max表示，并令T_left=T_min，T_right=T_max，进入下一步； 

a2)如果T_right-T_left＜ε(ε为一小的常数)，则令下次雷达开机辐射间隔的时间T₁=T_left或T₁=T_right，跳到第a5)步，否则令T_test=T_left+[(T_right-T_left)/2]； 

a3)令t₁=T_test，计算间隔时间为t₁时扩展卡尔曼滤波的预测协方差P(t₁|0) 

P(t₁|0)＝F(t₁)·P(0|0)·F′(t₁)+Γ(t₁)·σ_v ²·Γ′(t₁) 

其中，P(0|0)表示最近一次滤波协方差估计，F(t₁)表示状态转移矩阵，Γ(t₁)表示过程噪声分布矩阵，σ_v ²为零均值高斯白噪声的协方差； 

a4)判断预测协方差是否超过门限。假设跟踪精度门限为R_th，如果  $\sqrt{P_{\left(t_1\right|0)}\left(\mathrm{1,1}\right)+P_{\left(t_1\right|0)}\left(\mathrm{3,3}\right)+P_{\left(t_1\right|0)}\left(\mathrm{5,5}\right)}<R_{\mathrm{th}}$ (其中分别是预测协方差矩阵P(t₁|0)的对角线上X、Y、Z位置误差方差项)，则T_left=T_test，否则T_right=T_test，跳回到第a2)步； 

a5)得到下次预定雷达开机辐射间隔时间为T₁。 

步骤三：1)当雷达预定开机时间到，没有得到新的被动量测，则跳到步骤四；否则，当预定开机时间未到，有被动量测时，先将被动量测点迹与融合航迹进行关联判决，关联判决的具体步骤为： 

b1)如果得到新的被动量测，首先将融合航迹和对应的协方差转换到该被动量测所在平台的载体坐标系下； 

b2)然后将载体坐标系下的航迹和对应的协方差从以该载体为原点的X-Y-Z直角坐标下转换到对应的ρ-θ-ε(距离-方位-俯仰)的极坐标系 

X_i ^SP(k+1|k)＝U_k ^TX_i(k+1|k)_|(x，y，z)U_k

P_i ^SP(k+1|k)＝U_k ^TP_i(k+1|k)_|(x，y，z)U_k

式中，U_k为直角坐标到极坐标的转换矩阵。 

b3)构造关联统计量 

g₁(k+1)=d₁(k)^T[P₁+R₁]^-1d₁(k) 

其中，P₁为P_i ^SP(k|k+1)在方位和俯仰上的误差协方差矩阵，R₁为被动量测在方位俯仰上的量测误差矩阵，d₁(k)为被动量测在方位和俯仰上的新息。 

b4)利用如下公式进行量测与航迹的关联判决 

g₁(k+1)＜th₁    关联 

g₁(k+1)＞th₁    不关联 

其中th₁是二自由度χ²分布的统计量g₁(k)的假设检验门限。 

2)将关联上的被动量测利用基于坐标交替转换的多平台主被动序贯滤波进行跟踪，然后跳回到步骤二。其中基于坐标交替转换的多平台主被动序贯滤波进行跟踪的具体步骤为： 

c1)得到新的被动量测后，首先将该被动量测到来前最近的融合航迹点以及对应的协方差P(k|k)转换到获得该量测的传感器所在平台的载体坐标系下，得到对应的航迹点 P₁ ^NED(k|k)； 

c2)利用P₁ ^NED(k|k)进行外推，得到状态预测 

${\hat{X}}_1^{\mathrm{NED}}(k+1|k)=F\&CenterDot;{\hat{X}}_1^{\mathrm{NED}}\left(k\right|k)$

预测协方差 

P₁ ^NED(k+1|k)＝F·P₁ ^NED(k|k)·F^T

利用计算被动量测方位和俯仰的预测值 

$h(k+1|k)={[{\mathrm{\&theta;}}_{\mathrm{PAS}}(k+1|k),{\mathrm{\&epsiv;}}_{\mathrm{PAS}}(k+1|k)]}^{\&prime;}={\left[\begin{array}{c}\arctan 2[-z(k+1|k),x(k+1|k)]\\ \arctan 2[y\left(k\right),\sqrt{x^2(k+1|k)+z^2(k+1|k)}]\end{array}\right]}^{\&prime;}$

其中坐标系以载体为坐标原点的北-天-东为X-Y-Z方向(以下均以此为坐标系)，arctan2表示matlab软件中对应的反正切函数； 

利用计算被动量测时扩展卡尔曼滤波的量测矩阵 

$H_{\mathrm{PAS}}(k+1)=\left[\begin{array}{ccccccccc}-\frac{\hat{z}(k+1|k)}{{\hat{r}}^2}& 0& 0& 0& -\frac{\hat{x}(k+1|k)}{{\hat{r}}^2}& 0& 0& 0& 0\\ -\frac{\hat{x}(k+1|k)\hat{z}(k+1|k)}{\hat{r}{\hat{R}}^2}& 0& \frac{\hat{r}}{{\hat{R}}^2}& 0& -\frac{\hat{y}(k+1|k)\hat{z}(k+1|k)}{\hat{r}{\hat{R}}^2}& 0& 0& 0& 0\end{array}\right]$

其中，状态向量表示形式为 $\left[\begin{array}{ccccccccc}x& \stackrel{.}{x}& y& \stackrel{.}{y}& z& \stackrel{.}{z}& \stackrel{..}{x}& \stackrel{..}{y}& \stackrel{..}{z}\end{array}\right],\hat{r}=\sqrt{{\left[\hat{x}\left(k\right|k-1)\right]}^2+{\left[\hat{z}\left(k\right|k-1)\right]}^2},$ $\hat{R}=\sqrt{{\hat{r}}^2+{\left[\hat{y}\left(k\right|k-1)\right]}^2}.$

利用雷达测量精度性能指标，得到量测误差矩阵 

$R(k+1)=\left[\begin{array}{cc}{\widehat{\mathrm{\&sigma;}}}_{{\mathrm{\&theta;}}_{\mathrm{PAS}}}^2& 0\\ 0& {\widehat{\mathrm{\&sigma;}}}_{{\mathrm{\&epsiv;}}_{\mathrm{PAS}}}^2\end{array}\right]$

其中，分别表示被动量测的方位和俯仰的误差方差。 

将被动量测Z(k+1)=[θ_PAS(k+1)，ε_PAS(k+1)]′、H_PAS、R(k+1)、H_PAS、h(k+1|k)带入扩展卡尔曼滤波公式： 

$K(k+1|k)={P_1}^{\mathrm{NED}}(k+1)\&CenterDot;H_{\mathrm{PAS}}^T(k+1)\&CenterDot;{[H_{\mathrm{PAS}}(k+1)\&CenterDot;{P_1}^{\mathrm{NED}}(k+1|k\left)\right)\&CenterDot;H_{\mathrm{PAS}}^T(k+1)+R(k+1)]}^T$

P₁ ^NED(k+1|k+1)＝[I-K(k+1)H_PAS(k+1)]P₁ ^NED(k+1|k) 

${\hat{X}}_1^{\mathrm{NED}}(k+1|k+1)={\hat{X}}_1^{\mathrm{NED}}(k+1|k)+K(k+1)\&CenterDot;[Z(k+1)-h(k+1|k)]$

其中，I表示与状态向量维数相同的单位矩阵。 

c3)然后将和P₁ ^NED(k+1|k+1)分别转换回原融合航迹的坐标系下，成为原航迹的下一个滤波状态估计点及对应的协方差。 

步骤四：1)在空间上选择不同的平台上的雷达进行交替辐射，选择方法可采用顺序交替或按平台的优先等级交替； 

2)当雷达开机辐射，得到新的雷达量测，先将雷达量测点迹与融合航迹进行关联判决，具体步骤为： 

d1)如果得到新的雷达量测，首先将融合航迹和对应的协方差转换到该雷达量测所在平台的载体坐标系下； 

d2)然后将载体坐标系下的航迹和对应的协方差从以载体为原点的X-Y-Z直角坐标下转换到ρ-θ-ε(距离-方位-俯仰)的极坐标系 

X_i ^SP(k+1|k)＝U_k ^TX_i(k+1|k)_|(x，y，z)U_k

P_i ^SP(k+1|k)＝U_k ^TP_i(k+1|k)_|(x，y，z)U_k

式中，U_k为直角坐标到极坐标的转换矩阵。 

d3)构造关联统计量 

g₂(k+1)=d₂(k)^T[P₂+R₂]^-1d₂(k) 

其中，P₂为P_i ^SP(k|k+1)在距离、方位和俯仰上的误差协方差矩阵，R₂为雷达量测在距离、方位、俯仰上的量测误差矩阵，d₂(k)为雷达量测在距离、方位和俯仰上的新息。 

d4)利用如下公式进行量测与航迹的关联判决 

g₂(k+1)＜th₂    关联 

g₂(k+1)＞th₂    不关联 

其中th₂是三自由度的χ²分布的统计量g₂(k)的假设检验门限。 

3)将关联上的雷达量测利用基于坐标交替转换的多平台主被动序贯滤波算法进行跟踪，跳回到步骤二，其中雷达量测利用基于坐标交替转换的多平台主被动序贯滤波具体步骤如下： 

e1)得到新的雷达量测后，首先将该被动量测到来前最近的融合航迹点以及对应的滤波协方差P(k|k)转换到获得该量测的传感器所在平台的载体坐标系下，得到对应的航迹点P₁ ^NED(k|k)； 

e2)利用P₁ ^NED(k|k)进行外推，得到一步预测 

${\hat{X}}_1^{\mathrm{NED}}(k+1|k)=F\&CenterDot;{\hat{X}}_1^{\mathrm{NED}}\left(k\right|k)$

预测协方差 

P₁ ^NED(k+1|k)＝F·P₁ ^NED(k|k)·F^T

利用计算雷达量测距离、方位和俯仰的预测值 

$h(k+1|k)=\left[\begin{array}{c}{\mathrm{\&rho;}}_{\mathrm{radar}}(k+1|k)\\ {\mathrm{\&theta;}}_{\mathrm{radar}}(k+1|k)\\ {\mathrm{\&epsiv;}}_{\mathrm{radar}}(k+1|k)\end{array}\right]=\left[\begin{array}{c}\sqrt{x^2+y^2+z^2}\\ \arctan 2[-z\left(k\right|k-1),x\left(k\right|k-1)]\\ \arctan 2[y\left(k\right),\sqrt{x^2\left(k\right|k-1)+z^2\left(k\right|k-1)]}\end{array}\right]$

其中坐标系以载体为坐标原点的北-天-东为X-Y-Z方向(以下均以此为坐标系)，arctan2表示matlab软件中对应的反正切函数； 

利用计算雷达量测时扩展卡尔曼滤波的量测矩阵 

$H_{\mathrm{radar}}=\begin{array}{}\end{array}\left[\begin{array}{ccccccccc}\frac{\hat{x}\left(k\right|k-1)}{\hat{r}}& 0& \frac{\hat{y}\left(k\right|k-1)}{\hat{r}}& 0& \frac{\hat{z}\left(k\right|k-1)}{\hat{r}}& 0& 0& 0& 0\\ -\frac{\hat{z}\left(k\right|k-1)}{{\hat{r}}^2}& 0& 0& 0& -\frac{\hat{x}\left(k\right|k-1)}{{\hat{r}}^2}& 0& 0& 0& 0\\ -\frac{\hat{x}\left(k\right|k-1)\hat{z}\left(k\right|k-1)}{\hat{r}{\hat{R}}^2}& 0& \frac{\hat{r}}{{\hat{R}}^2}& 0& -\frac{\hat{y}\left(k\right|k-1)\hat{z}\left(k\right|k-1)}{\hat{r}{\hat{R}}^2}& 0& 0& 0& 0\end{array}\right]$

其中，状态向量表示形式为 $\left[\begin{array}{ccccccccc}x& \stackrel{.}{x}& y& \stackrel{.}{y}& z& \stackrel{.}{z}& \stackrel{..}{x}& \stackrel{..}{y}& \stackrel{..}{z}\end{array}\right],\hat{r}=\sqrt{{\left[\hat{x}\left(k\right|k-1)\right]}^2+{\left[\hat{z}\left(k\right|k-1)\right]}^2},$ $\hat{R}=\sqrt{{\hat{r}}^2+{\left[\hat{y}\left(k\right|k-1)\right]}^2}.$

利用雷达测量精度性能指标，得到雷达量测误差矩阵 

$R(k+1)\left[\begin{array}{ccc}{\widehat{\mathrm{\&sigma;}}}_{{\mathrm{\&rho;}}_{\mathrm{radar}}}& 0& 0\\ 0& {\widehat{\mathrm{\&sigma;}}}_{{\mathrm{\&theta;}}_{\mathrm{radar}}}& 0\\ 0& 0& {\widehat{\mathrm{\&sigma;}}}_{{\mathrm{\&epsiv;}}_{\mathrm{radar}}}\end{array}\right]$

其中，分别表示雷达量测的距离、方位和俯仰的误差方差。 

将雷达量测Z(k+1)=[ρ_radar(k+1)，θ_radar(k+1)，ε_radar(k+1)]′、H_radar、R(k+1)、H_PAS、h(k+1|k)带入扩展卡尔曼滤波公式： 

$K(k+1|k)={P_1}^{\mathrm{NED}}(k+1)\&CenterDot;H_{\mathrm{PAS}}^T(k+1)\&CenterDot;{[H_{\mathrm{PAS}}(k+1)\&CenterDot;{P_1}^{\mathrm{NED}}(k+1|k\left)\right)\&CenterDot;H_{\mathrm{PAS}}^T(k+1)+R(k+1)]}^T$

P₁ ^NED(k+1|k+1)＝[I-K(k+1)H_PAS(k+1)]P₁ ^NED(k+1|k) 

${\hat{X}}_1^{\mathrm{NED}}(k+1|k+1)={\hat{X}}_1^{\mathrm{NED}}(k+1|k)+K(k+1)\&CenterDot;[Z(k+1)-h(k+1|k)]$

其中，I表示与状态向量维数相同的单位矩阵。 

e3)然后将和P₁ ^NED(k+1|k+1)分别坐标转换回到原融合航迹的坐标系下，成为原航迹的下一个滤波状态估计点及对应的协方差。 

本发明的效果可以通过以下matlab仿真结果进一步说明： 

仿真实验条件： 

设有2架己方作战飞机，初始地理坐标纬度、经度、高度分别为(36.14°，120°，10000m)、(36.23°，120°，10000m)，其纬度、经度、高度方向上初始速度分别为(0，-150m/s，0)、(0，-150m/s，0)，且均作匀速直线运动，仿真总时间110s，飞行的轨迹图如图2所示。每架飞机上分别搭载一个雷达、一个红外传感器和一个电子支援措施，雷达能测量目标的距离、方位、俯仰，红外传感器和电子支援措施能测量目标的方位、俯仰，目标飞机1架，其地理坐标纬度、经度、高度为(36.18°，119.5°，10000m)，其纬度、经度、高度方向上初始速度为(0，200m/s，0)；机载雷达最小辐射间隔为1秒，被动量测的数据率分两种情况：1)机载被动传感器量测每1秒有一个量测，2)机载被动传感器量测每10秒有一个量测；机载雷达测距误差为100m，方位角均方根误差为0.2°，俯仰角均方根误差为0.2°；机载红外传感器方位角均方根误差为0.1°，俯仰角均方根误差为0.1°；机载电子资源措施的方位角均方根误差为0.3°，俯仰角均方根误差为0.3°。跟踪任务需求精度为100米以内。 

实验一： 

用本发明方法对2架己方飞机探测目标过程中的雷达辐射进行管控，在被动传感器量测数据率为每1秒一个时，2架己方飞机雷达辐射情况如图3所示，融合跟踪误差如图4所示。在被动传感器量测数据率为每10秒一个时，2架己方飞机雷达辐射情况如图5所示，融合跟踪误差如图6所示。 

从图3可以看出，在需求精度限定条件下，己方飞机1和飞机2的雷达在时间上间歇开机，随着跟踪精度的提高开机时间间隔越来越大，并在空间上交替开机。同时，由于2架飞机的红外和电子支援措施等被动量测的辅助跟踪，雷达在一段时间后完全不需要开机就能满足跟踪需求。 

比较图3和图5可以看出，雷达辐射次数除了与跟踪精度需求有关，还与被动量测的数据率有关，即被动量测数据高时雷达辐射次数少，被动量测数据率低时雷达辐射次数多，雷达能自适应根据被动量测的数据率来调整辐射的时间。 

实验二： 

在同一仿真场景下，用现有的单平台方法来让己方飞机1完成对目标的跟踪和雷达辐射管控，对应的在被动传感器量测数据率为每1秒一个时，己方飞机1雷达辐射情况如图7所示，融合跟踪误差如图8所示。在被动传感器量测数据率为每10秒一个时，1架飞机搭载的 雷达辐射情况如图9所示，融合跟踪误差如图10所示。本发明的方法与传统的方法在相同的仿真场景下对应的雷达辐射次数比较如表1所示 

表1 现有方法与本发明方法对应的雷达辐射次数比较 

从表1可以看出现有单平台的雷达辐射方法中在被动传感器量测间隔为1s，本发明方法比现有方法飞机1雷达辐射次数减少约70％，在被动传感器量测间隔为10s时，飞机1雷达辐射次数减少约53％，极大的减少了飞机的雷达辐射时间，提高了平台的隐蔽性。 

Claims (3)

1.一种多平台主被动传感器协同跟踪的雷达间歇交替辐射控制方法，其特征在于包括以下步骤：

(1)被动传感器引导雷达开机完成目标航迹起始；

(2)利用协方差控制方法在时间上确定下一次需要雷达辐射的时刻；

(3a)当预定开机时间到，没有得到新的被动量测，则跳到第(4a)步；

(3b)当预定开机时间未到，有被动量测时，首先用所得到的被动量测点迹与融合航迹进行关联，然后将关联上的被动量测利用基于坐标交替转换的多平台主被动序贯滤波进行跟踪，回到第(2)步：

(4a)在空间上选择不同平台交替开机辐射；

(4b)雷达开机得到新的雷达量测后，先利用所得雷达量测点迹与融合航迹进行关联，然后将关联上的雷达量测利用基于坐标交替转换的多平台主被动序贯滤波算法进行跟踪，跳回到第(2)步。

2.根据权利要求1所述的一种多平台主被动传感器协同跟踪的雷达间歇交替辐射控制方法，其特征在于，所述步骤(2)具体为：

a1)假设雷达最小和最大辐射间隔分别用T_min、T_max表示，并令T_left=T_min，T_right=T_max，进入下一步；

a2)如果T_right-T_left＜ε(ε为一小的常数)，则令下次雷达开机辐射间隔的时间T₁＝T_left或T₁=T_right，跳到第a5)步，否则令T_test=T_lest+[(T_right-T_left)/2]；

a3)令t₁=T_lest，计算间隔时间为t₁时扩展卡尔曼滤波的预测协方差P(t₁|0)

P(t₁|0)＝F(t₁)·P(0|0)·F′(t₁)+Γ(t₁)·σ_v ²·Γ′(t₁)

其中，P(0|0)表示最近一次滤波协方差估计，F(t₁)表示状态转移矩阵，Γ(t₁)表示过程噪声分布矩阵，σ_v ²为零均值高斯白噪声的协方差；

a4)判断预测协方差是否超过门限。假设跟踪精度门限为R_th，如果 $\sqrt{P_{\left(t_1\right|0)}\left(\mathrm{1,1}\right)+P_{\left(t_1\right|0)}\left(\mathrm{3,3}\right)+P_{\left(t_1\right|0)}\left(\mathrm{5,5}\right)}<R_{\mathrm{th}}$ (其中分别是预测协方差矩阵P(t₁|0)的对角线上X、Y、Z位置误差方差项)，则T_left=T_test，否则T_right=T_lest，跳回到第a2)步；

a5)得到下次预定雷达开机辐射间隔时间为T₁。

3.根据权利要求1所述的一种多平台主被动传感器协同跟踪的雷达间歇交替辐射控制方法，其特征在于，所述步骤(3b)具体为：

b1)得到与融合航迹关联的被动量测后，首先将该被动量测到来前最近的融合航迹点以及对应的协方差P(k|k)转换到获得该量测的传感器所在平台的载体坐标系下，得到对应的航迹点P₁ ^NED(k|k)；

b2)利用P₁ ^NED(k|k)进行外推，得到状态预测

${\hat{X}}_1^{\mathrm{NED}}(k+1|k)=F\&CenterDot;{\hat{X}}_1^{\mathrm{NED}}\left(k\right|k)$

预测协方差

P₁ ^NED(k+1|k)＝F·P₁ ^NED(k|k)·F^T

b3)然后将载体坐标系下的预测航迹和对应的协方差P₁ ^NED(k+1|k)从X-Y-Z直角坐标下转换到ρ-θ-ε(距离-方位-俯仰)的极坐标系

${X_i}^{\mathrm{SP}}(k+1|k)={U_k}^T{\hat{X}}_1^{\mathrm{NED}}{(k+1|k)}_{|(x,y,z)}{U}_k$

P_i ^SP(k+1|k)＝U_k ^TP₁ ^NED(k+1|k)_|(x，y，z)U_k

式中，U_k为直角坐标到球极坐标的转换矩阵。

构造关联统计量

g₁(k+1)=d₁(k)^T[P₁+R₁]^-1d₁(k)

其中，P₁为P_i ^SP(k|k+1)在方位和俯仰上的误差协方差矩阵，R₁为被动量测在方位俯仰上的量测误差矩阵，d₁(k)为在方位和俯仰上的新息。

利用如下公式进行量测与航迹的关联判决

g₁(k+1)＜th    关联

g₁(k+1)＞th    不关联

其中th是二自由度χ²分布的统计量g₁(k)的假设检验门限。

b4)得到与融合航迹关联的被动量测Z(k+1)=[θ_PAS(k+1)，ε_PAS(k+1)]′后，先利用计算被动量测方位和俯仰的预测值

$h(k+1|k)={[{\mathrm{\&theta;}}_{\mathrm{PAS}}(k+1|k),{\mathrm{\&epsiv;}}_{\mathrm{PAS}}(k+1|k)]}^{\&prime;}={\left[\begin{array}{c}\arctan 2[-z(k+1|k),x(k+1|k)]\\ \arctan 2[y\left(k\right),\sqrt{x^2(k+1|k)+z^2(k+1|k)}]\end{array}\right]}^{\&prime;}$

其中坐标系以载体为坐标原点的北-天-东为X-Y-Z方向(以下均以此为坐标系)，arctan2表示matlab软件中对应的反正切函数。

利用计算被动量测时扩展卡尔曼滤波的量测矩阵

$H_{\mathrm{PAS}}(k+1)=\left[\begin{array}{ccccccccc}-\frac{\hat{z}(k+1|k)}{{\hat{r}}^2}& 0& 0& 0& -\frac{\hat{x}(k+1|k)}{{\hat{r}}^2}& 0& 0& 0& 0\\ -\frac{\hat{x}(k+1|k)\hat{z}(k+1|k)}{\hat{r}{\hat{R}}^2}& 0& \frac{\hat{r}}{{\hat{R}}^2}& 0& -\frac{\hat{y}(k+1|k)\hat{z}(k+1|k)}{\hat{r}{\hat{R}}^2}& 0& 0& 0& 0\end{array}\right]$

其中，H_PAS(k+1)对应的状态向量表示形式为 $\left[\begin{array}{ccccccccc}x& \stackrel{.}{x}& y& \stackrel{.}{y}& z& \stackrel{.}{z}& \stackrel{..}{x}& \stackrel{..}{y}& \stackrel{..}{z}\end{array}\right],$ $\hat{r}=\sqrt{{\left[\hat{x}\left(k\right|k-1)\right]}^2+{\left[\hat{z}\left(k\right|k-1)\right]}^2},\hat{R}=\sqrt{{\hat{r}}^2+{\left[\hat{y}\left(k\right|k-1)\right]}^2}.$

利用雷达测量精度性能指标，得到量测误差矩阵

$R(k+1)=\left[\begin{array}{cc}{\widehat{\mathrm{\&sigma;}}}_{{\mathrm{\&theta;}}_{\mathrm{PAS}}}^2& 0\\ 0& {\widehat{\mathrm{\&sigma;}}}_{{\mathrm{\&epsiv;}}_{\mathrm{PAS}}}^2\end{array}\right]$

其中，分别表示被动量测的方位和俯仰的误差方差。

将被动量测Z(k+1)=[θ_PAS(k+1)，ε_PAS(k+1)]′、被动量测时扩展卡尔曼滤波的量测矩阵H_PAS、量测误差矩阵R(k+1)、被动量测方位和俯仰的预测值h(k+1|k)分别带入扩展卡尔曼滤波公式：

$K(k+1|k)={P_1}^{\mathrm{NED}}(k+1)\&CenterDot;H_{\mathrm{PAS}}^T(k+1)\&CenterDot;{[H_{\mathrm{PAS}}(k+1)\&CenterDot;{P_1}^{\mathrm{NED}}(k+1|k\left)\right)\&CenterDot;H_{\mathrm{PAS}}^T(k+1)+R(k+1)]}^T$

P₁ ^NED(k+1|k+1)＝[I-K(k+1)H_PAS(k+1)]P₁ ^NED(k+1|k)

${\hat{X}}_1^{\mathrm{NED}}(k+1|k+1)={\hat{X}}_1^{\mathrm{NED}}(k+1|k)+K(k+1)\&CenterDot;[Z(k+1)-h(k+1|k)]$

其中，I表示与状态向量维数相同的单位矩阵。

b5)然后将和P₁ ^NED(k+1|k+1)分别转换回原融合航迹的坐标系下，成为原航迹的下一个滤波状态估计点及对应的协方差。