1.   引言

天基探测雷达具有全天候、全天时的战略、战术预警能力，且具有不受地球曲率限制、不易受攻击等优势，在预警防御系统中具有巨大潜力。本文考虑面向大范围自主搜索发现任务的天基雷达。设计过程中，为了实现较好的最小可检测速度，通常方位向孔径尺寸较大；而综合考虑重量、收拢等因素，俯仰向口径尺寸较小，俯仰向测角能力较弱。因此，单个天基探测雷达对空中目标跟踪定位时，相比于目标径向距量测及方位角量测误差，目标俯仰角量测误差较大；仅利用单雷达径向距量测、方位角量测、俯仰角量测难以有效估计空中目标高度。空中目标高度信息的缺失严重放大了目标定位误差，同时不利于目标属性判别与威胁估计。天基探测雷达空中目标高度估计研究较为少见。Graves^[1]提出了一种基于多路径干扰模式的天基雷达空中目标高度估计方法，理论上难以适用于低空目标。而通过多个天基探测雷达组网，融合利用多个天基探测雷达径向距、方位角量测，可提高空中目标高度估计准确度，对改善目标定位精度及威胁估计能力等具有重要意义。

天基探测雷达目标定位涉及多次坐标转换，具有较强的非线性，需要采用非线性滤波算法估计目标运动状态。常见的非线性滤波方法有扩展卡尔曼滤波(Extended Kalman Filter, EKF)、无迹卡尔曼滤波(Unscented Kalman Filter, UKF)、转换量测的卡尔曼滤波(Converted Measurement Kalman Filter, CMKF)和粒子滤波(Particle Filter, PF)等。Kulikov等人^[2]基于全局误差控制常微分方程(Ordinary Differential Equations, ODE)求解器提出了一种精确的连续离散EKF方法，在高维雷达目标跟踪场景下对所提算法进行了验证，数值结果表明该方法可用于实际的目标跟踪，灵活度和鲁棒性更高。Kulikov等人^[3]和Zhang等人^[4]采用无迹变换(Unscented Transform, UT)来近似非线性系统的状态转移过程和量测方程，结果表明其滤波精度要优于EKF方法。然而，EKF和UKF滤波器忽略了状态协方差矩阵的欠估计问题，导致滤波精度下降。张连仲等人^[5]采用一种基于期望最大化的转换量测卡尔曼滤波方法，将球坐标系下的雷达量测值转换到直角坐标系下，利用期望最大化方法修正了量测噪声协方差，最后利用标准卡尔曼滤波器进行滤波，提高了野值干扰下的目标跟踪精度，但该方法需要利用真值计算量测转换误差矩阵，进而导致兼容性问题。Wu等人^[6]提出了一种重新加权的鲁棒粒子滤波方法来降低量测野值对目标跟踪精度的影响，结果表明该方法具有较高的跟踪精度，但是，由于使用蒙特卡罗采样，粒子滤波所需计算时间成本过高。Ait-El-Fquih等人^[7]提出了一种基于变分贝叶斯(Variational Bayes, VB)的粒子滤波方法，该方法将VB方法与随机抽样技术相结合以分割状态空间，从而推导出一种新的计算效率高的粒子滤波算法。闫文旭等人^[8]在VB框架下，将系统状态的后验分布近似于一个可解的变分分布，通过迭代优化证据下界，并考虑惩罚因子调整迭代的步长，在非线性基础上推导出闭合形式的迭代EKF。

近年来，概率图模型及运行其上的消息传递方法得到了极大发展，并被应用于目标跟踪定位问题中^[9-11]。概率图模型结合概率论和图论，利用图模型表示变量的联合概率分布，利用条件独立性将高维推断问题分解，从而构造高效的学习和推断算法。常见图模型包括贝叶斯网络、马尔可夫随机场以及因子图，且它们之间可以相互转化。同时，对于复杂推断问题，具体的消息传递方法近年来也备受关注^[12-15]，常见算法包括信度传播(Belief Propagation, BP)、环信度传播(Loopy BP, LBP)、变分消息传递(Variational Message Passing, VMP)、期望传播(Expectation Propagation, EP)、非参数信度传播(Nonparametric BP, NBP)、粒子信度传播(Particle BP, PBP)、广义信度传播(Generalized BP, GBP)、平均场近似与信度传播组合的消息传递等。

多天基雷达组网融合跟踪本质上是一个复杂的推断问题。考虑到天基平台计算资源的稀缺性等实际需求，需要准确、低计算复杂度的融合跟踪方法。鉴于因子图模型对复杂推断问题强大的表征能力，本文采用因子图对多部雷达目标跟踪问题进行建模，首先基于因子图描述传统的一致卡尔曼滤波，随后引入耦合参数，并考虑到目标跟踪模型的非线性特性，推导了分布式一致扩展卡尔曼滤波(Distributed Consistent Extended Kalman Filter, DCEKF)算法。随后，针对EKF线性化非线性模型可能导致误差变大的问题，在建立的因子图的基础上，利用分布式NBP (DCNBP)算法，有效解决了多雷达目标跟踪的非线性问题，相较于DCEKF算法跟踪精度明显提升。

2.   系统建模与问题描述

2.1   天基雷达目标跟踪坐标系转换

天基雷达提供的目标测量(径向距、方位角和俯仰角)位于雷达测量坐标系，而常见的空中目标经度、纬度、高度信息在大地坐标系中描述。为了量测信息的统一，将大地坐标系向雷达测量坐标系转换。首先是大地坐标系到地心(第4)赤道坐标系(简称ECEF坐标系)的转换，本文采用近似迭代方法^[16]，该方法具有更好的转换效率和精度，并且在地球半径50 km以内，不存在奇点和不收敛情况。下面介绍本文涉及的几种坐标系及坐标转换^[8]。

如图1所示，ECEF坐标系以地球球心 $O$ 为原点， ${{OX}}^{\rm{E}}$ 轴指向0经度方向， ${{OZ}}^{\rm{E}}$ 轴垂直赤道平面并指向地球自转角速度方向， ${{OX}}^{\rm{E}}$ , ${{OY}}^{\rm{E}}$ 与 ${{OZ}}^{\rm{E}}$ 轴构成右手法则坐标系。为了坐标转换需要，同样在图1中建立地心(第1)赤道坐标系(简称ECI，又称J2000.0坐标系)，其以地球球心 $O$ 为原点， ${{OX}}^{\rm{J}}$ 轴指向春分点， ${{OZ}}^{\rm{J}}$ 轴垂直赤道平面并指向地球自转角速度方向， ${{OX}}^{\rm{J}}$ , ${{OY}}^{\rm{J}}$ 与 ${{OZ}}^{\rm{J}}$ 轴构成右手法则坐标系。天基雷达通常安装在卫星上，故建立卫星轨道坐标系(简称LVLH坐标系)，其以卫星质心 $o$ 为原点， ${{ox}}^{\rm{L}}$ 轴为地心指向卫星质心的卫星矢径方向， ${{oz}}^{\rm{L}}$ 轴为卫星运动速度方向， ${{ox}}^{\rm{L}}$ , ${{oy}}^{\rm{L}}$ 和 ${{oz}}^{\rm{L}}$ 轴构成右手法则坐标系，下面介绍从ECEF坐标系到雷达量测坐标系的转换步骤，坐标转换的流程如图2所示。

图  1  坐标系示意图

Figure  1.  Coordinate systems

下载: 全尺寸图片 幻灯片

图  2  坐标系转换流程图

Figure  2.  Procedure of coordinate systems transformation

下载: 全尺寸图片 幻灯片

为了得到较为精确的坐标转换，考虑地球的极移、自转、章动、岁差的影响，从ECEF到ECI的坐标变换为

其中， $\left[{\boldsymbol{A}} \right]$ , $\left[ {\boldsymbol{B}} \right]$ , $\left[ {\boldsymbol{C}} \right]$ 和 $\left[ {\boldsymbol{D}} \right]$ 分别为 $t_0$ 时刻的极移、自转、章动和岁差矩阵，且

其中， $x_p$ 和 $y_p$ 为地极坐标， ${\theta _G}$ 为格林尼治恒星时角， ${\varepsilon _M}$ , $\Delta \varepsilon$ 和 $\Delta \phi$ 分别为平黄赤交角、交角章动和黄经章动， $Z_p$ , $\theta_p$ 和 $\xi_p$ 为岁差参数， $t_0$ 为历元起始时刻， ${{\boldsymbol{R}}_x}(\cdot)$ , ${{\boldsymbol{R}}_y}(\cdot)$ 和 ${{\boldsymbol{R}}_z}(\cdot)$ 分别为绕 $x$ 轴、 $y$ 轴和 $z$ 轴旋转的方向余弦矩阵。

从ECI坐标系到卫星LVLH坐标系的旋转矩阵为

其中，

其中， $\nu$ 为纬度幅角， $\varOmega$ 和 $i$ 为轨道倾角与升交点赤经。

卫星本体坐标系 ${{o}}{{{x}}_b}{{{y}}_b}{{{z}}_b}$ 以卫星质心为坐标原点 $o$ ，3轴为卫星的3个惯性主轴， ${{o}}{{{y}}_b}$ 轴垂直于星箭分离面， ${{o}}{{{z}}_b}$ 轴垂直指向星体对地面， ${{o}}{{{x}}_b}$ , ${{o}}{{{y}}_b}$ 与 ${{o}}{{{z}}_b}$ 构成右手直角坐标系。理想情况下，卫星入轨姿态角全为 $0{^\circ}$ ，卫星LVLH坐标系与卫星本体坐标系重合。卫星LVLH坐标系到卫星本体坐标系的转换为

其中， $\gamma_s$ , $\beta_s$ 和 $\alpha_s$ 分别为卫星在本体坐标系中姿态的滚动角、俯仰角和偏航角， ${\boldsymbol{M}}_1$ , ${\boldsymbol{M}}_2$ , ${\boldsymbol{M}}_3$ 分别为按照右手法则绕 $x$ , $y$ , $z$ 轴逆时针旋转一定角度的旋转矩阵， $x_c$ , $y_c$ 和 $z_c$ 为目标在卫星LVLH坐标系下的位置， ${x_b}$ , ${y_b}$ 和 ${z_b}$ 为目标在卫星本体坐标系中的位置。

如图3所示，天基雷达天线阵面坐标系 ${{o}}{{{x}}_a}{{{y}}_a}{{{z}}_a}$ 坐标原点 $o$ 在天线中心， ${{o}}{{{x}}_a}$ 轴指向天线横轴， ${{o}}{{{y}}_a}$ 轴指向天线纵轴， ${{o}}{{{z}}_a}$ 轴指向天线阵面法向， ${{o}}{{{x}}_a}$ , ${{o}}{{{y}}_a}$ 与 ${{o}}{{{z}}_a}$ 构成右手直角坐标系。一般天线安装在卫星的对地面(通常天线阵面坐标系 ${{o}}{{{x}}_a}$ 的轴与卫星本体坐标系的 ${{o}}{{{x}}_b}$ 轴重合)，主要用于对地任务，入轨后转动角度满足对地工作要求，当所有阵面姿态角均为 $0{^\circ}$ 时，天线阵面坐标系与卫星本体坐标系重合。卫星本体坐标系到天线阵面坐标系的转换为

图  3  天基雷达天线阵面坐标系

Figure  3.  Space-based radar antenna array coordinate system

下载: 全尺寸图片 幻灯片

其中， $\gamma_c$ , $\beta_c$ 和 $\alpha_c$ 分别为雷达天线阵面的滚动角、俯仰角和偏航角， $x_a$ , $y_a$ 和 $z_a$ 为目标在天线阵面坐标系下的位置， ${\boldsymbol{T}}_1$ , ${\boldsymbol{T}}_2$ , ${\boldsymbol{T}}_3$ 分别为按照右手法则绕 $x$ , $y$ 和 $z$ 轴逆时针旋转一定角度的旋转矩阵。

雷达阵面测量坐标系为以天线中心为原点的球坐标系，从天线阵面坐标系到雷达量测坐标系的转换为

其中， $r_z$ , $\theta_z$ 和 $\phi_z$ 分别为雷达阵面测量坐标系下的空中目标的径向距、方位角、俯仰角。

2.2   目标运动模型

不失一般性，考虑空中弱机动目标，此类目标的运动可以近似为匀速运动，其运动模型可以描述为

其中， ${{\boldsymbol{x}}}_k$ 为 $k$ 时刻目标在ECEF坐标系下的位置速度矢量， ${{\boldsymbol{x}}}_k = \left[ {{x_k},{y_k},{z_k},{{\dot x}_k},{{\dot y}_k},{{\dot z}_k}} \right]^{\mathrm{T}}$ ， ${\boldsymbol{F}}$ 为目标状态转移矩阵

其中， $T$ 为采样间隔， ${\boldsymbol{\omega}}_{k - 1}$ 为过程噪声，且假设 ${\boldsymbol{\omega}}_{k - 1} \sim {\cal{N}}\left({0,{\boldsymbol{Q}}} \right)$ 。

2.3   雷达量测模型

假设在 $k$ 时刻共有 $N$ 部天基雷达对目标进行测量，每部雷达得到的目标量测信息为径向距 $r_k^n$ 、方位角 $\theta _k^n$ 。为了方便描述多雷达组网分布式融合中每个雷达的局部目标跟踪，重写目标运动模型式(7)如下。在第 $n$ 个雷达节点，目标的运动模型可以近似为

其中， ${{\boldsymbol{x}}}_k^n$ 为第 $n$ 个雷达节点目标的局部状态向量， ${{\boldsymbol{x}}}_k^n = \left[ {x_k^n,y_k^n,z_k^n,\dot x_k^n,\dot y_k^n,\dot z_k^n} \right]^{\mathrm{T}}$ ， ${\boldsymbol{\omega}}_{k - 1}^n$ 为第 $n$ 个雷达节点的局部过程噪声，且假设 ${\boldsymbol{\omega}}_{k - 1}^n \sim {\cal{N}}\left( {0,{{\boldsymbol{Q}}^n}} \right)$ 。

由于俯仰角误差较大，仅利用雷达径向距和方位角量测估计目标高度及位置，第 $n$ 个雷达在 $k$ 时刻的量测方程表示为

其中， ${{{\boldsymbol{x}}}}_{a,k}^n = {\left[ {x_{a,k}^n,y_{a,k}^n,z_{a,k}^n} \right]^{\mathrm{T}}}$ 为目标在第 $n$ 个雷达量测坐标系下的位置矢量，且 ${{{\boldsymbol{x}}}}_{a,k}^n = g({{{\boldsymbol{x}}}}_k^n)$ ，其中， $g\left( \cdot \right)$ 为ECEF坐标系到雷达测量坐标系的转换函数， $\Delta r_k^n$ 和 $\Delta \theta_k^n$ 为第 $n$ 个雷达在 $k$ 时刻的测量误差，且假设 $\Delta r_k^n \sim {\cal{N}}\left( {0,\sigma _{r_k^n}^2} \right)$ , $\Delta \theta_k^n \sim {\cal{N}}\left( {0,\sigma _{\theta _k^n}^2} \right)$ , ${\boldsymbol{R}}_k^n = \text{diag}\{\sigma _{r_k^n}^2,\sigma _{\theta _k^n}^2\}$ 。定义 $k$ 时刻所有雷达的量测为 ${{\boldsymbol{y}}_k} = {\left[ {y_k^1, y_k^2, \cdots ,y_k^N} \right]^{\mathrm{T}}} = {\left[ {r_k^1,\theta_k^1, r_k^2,\theta_k^2, \cdots ,r_k^N,\theta_k^N} \right]^{\mathrm{T}}}$ ，则全局量测模型可以表示为

其中， ${{{\boldsymbol{v}}}_k} \sim {\cal{N}}\left( {0,{\boldsymbol{R}}_k} \right)$ , ${\boldsymbol{R}}_k = \text{diag}\{\sigma _{r_k^1}^2,\sigma _{\theta _k^1}^2, \sigma _{r_k^2}^2, \sigma _{\theta _k^2}^2, \cdots, \sigma _{r_k^N}^2,\sigma _{\theta _k^N}^2\}$ 。

2.4   基于因子图的天基雷达组网目标跟踪建模

因子图是概率图模型的一种。一般情况下，高维复杂随机系统的推断问题可以通过带有局部交互函数的概率图模型来建模，即利用因子分解方式和条件独立性假设紧凑表示多个变量的联合分布。因子图模型中包含变量节点和因子节点，其中每个变量节点表示一个随机变量，每个因子节点表示一个局部函数，连接变量节点和因子节点的边表示该因子是该变量的函数。因子图模型结合消息传递方法为解决复杂环境下的目标跟踪与多源信息融合问题提供了统一框架，将复杂高维推断问题转换为优化问题，在估计精度、计算复杂度和实现的灵活性方面具有显著优势，能够为目标跟踪问题提供效果好、效率高和可扩展的统一解决方案。为此，在上述空中目标与天基雷达量测模型的基础上，本文采用因子图方法对天基雷达组网目标跟踪进行建模。

记 ${\boldsymbol{X}} = \{{{\boldsymbol{x}}}_{0:K}^{n}\}_{n = 1}^N$ 表示从 $0$ 到 $K$ 时刻目标在 $N$ 部雷达下的状态矢量，记 ${\boldsymbol{Y}} = \{{\boldsymbol{y}}_{1:K}\}$ 表示从 $1$ 到 $K$ 时刻所有雷达的量测信息。则对于状态空间模型式(7)和量测模型式(11)， ${\boldsymbol{X}}$ 和 ${\boldsymbol{Y}}$ 的联合概率分布函数可分解为

其中， $p\left({{\boldsymbol{x}}}_k \mid {{\boldsymbol{x}}}_{k-1}\right)$ 为状态转移概率分布函数， $p\left({\boldsymbol{y}}_k \mid {{\boldsymbol{x}}}_k\right)$ 为似然函数。式(12)可以用因子图4表示，其中， $f_k$ 和 $h_k$ 分别表示 $k$ 时刻状态转移函数和似然函数对应的因子节点。

图  4  目标状态与量测联合概率分布函数的因子图分解

Figure  4.  Factor graph decomposition of the joint probability distribution function of target states and measurements

下载: 全尺寸图片 幻灯片

3.   基于因子图的分布式滤波算法

3.1   分布式一致扩展卡尔曼滤波

图4所示因子图模型中， $k$ 时刻，通过因子图上的递归计算得到状态 ${{{\boldsymbol{x}}}}_k$ 的后验边缘分布和状态 ${{{\boldsymbol{x}}}}_{k+1}$ 的预测分布，即

其中，定义 $p\left({{{\boldsymbol{x}}}}_k \mid {{\boldsymbol{y}}}_{1:k-1}\right) = N\left({ {\hat {{\boldsymbol{x}}}}}_k^{-},P_k^{-}\right)$ , $p\left({{{\boldsymbol{x}}}}_k \mid {{\boldsymbol{y}}}_{1:k}\right) = {\cal{N}}\left({ {\hat {{\boldsymbol{x}}}}}_k,{\boldsymbol{P}}_k\right)$ ，则 $k$ 时刻后验分布更新与传统卡尔曼滤波的状态估计相类似。

式(13)与式(14)为目标全局状态变量的更新与预测。根据2.3节的雷达观测模型，在局部状态变量中，引入一个耦合因子节点 $g_{ji}$ 来表示全局变量 ${{{\boldsymbol{x}}}}_k$ 的两个复制状态变量 ${{{\boldsymbol{x}}}}_k^i$ 和 ${{{\boldsymbol{x}}}}_k^j$ 在相邻节点 $i$ 和 $j$ 上的关系，其中， $i \in {D_j}$ ， ${D_j}$ 表示 $j$ 的邻居集合。将 ${g_{ji}}$ 定义为一个指数函数^[17]

其中， $\kappa > 0$ 为耦合参数。

图5为局部变量间的消息传递因子图，其中， $u,i \in {D_j}$ 。随着 $\kappa$ 的增加，相邻变量 ${{{\boldsymbol{x}}}}_k^i$ 和 ${{{\boldsymbol{x}}}}_k^j$ 的相关性越来越强，从而促进图模型中状态分布的一致性。当 $\kappa \to \infty$ 时，可以得到， ${g_{ji}}\left( {{{{\boldsymbol{x}}}}_k^j,{{{\boldsymbol{x}}}}_k^i} \right) \propto \delta \left( {{{{\boldsymbol{x}}}}_k^j - {{{\boldsymbol{x}}}}_k^i} \right)$ ，这时局部状态变量的分布在整个图模型中是相同的，即 $p\left( {{{{\boldsymbol{x}}}}_k^j|{{\boldsymbol{y}}}_{1:k}^1, {{\boldsymbol{y}}}_{1:k}^2, \cdots, {{\boldsymbol{y}}}_k^N} \right) = p\left( {{{{\boldsymbol{x}}}}_k^i|{{\boldsymbol{y}}}_{1:k}^1, {{\boldsymbol{y}}}_{1:k}^2, \cdots, {{\boldsymbol{y}}}_{1:k}^N} \right)$ ，且与全局变量分布 $p\left( {{{{{\boldsymbol{x}}}}_k}|{{\boldsymbol{y}}_{1:k}}} \right)$ 相同。

图  5  $k$ 时刻局部因子图模型

Figure  5.  A local factor graph model at time $k$

下载: 全尺寸图片 幻灯片

因此， $k$ 时刻局部变量 ${{{\boldsymbol{x}}}}_k^j$ 的后验边缘分布为

采用一阶线性化近似消息 ${m_{h_k^j \to {{{\boldsymbol{x}}}}_k^j}}$ ，由式(9)和式(10)可得

其中， ${\boldsymbol{H}}_k^j$ 为 $h_k^j\left( {{{{\boldsymbol{x}}}}_k^j} \right)$ 的雅可比矩阵，

变量节点 ${{{\boldsymbol{x}}}}_k^j$ 到其邻居因子节点 $g_{j i}$ 的消息为

因子节点 $g_{j i}$ 到变量节点 ${\boldsymbol{x}}_k^i$ 的消息为

将式(15)和式(21)代入到式(20)中，可以得到迭代方程。初始化 ${\boldsymbol{\mu}} _{ji}^0$ 和 ${\boldsymbol{\varLambda}} _{ji}^0$ 后，第 $l$ 次迭代时的方程为

经过 $L$ 次迭代之后，节点 ${{\boldsymbol{x}}}_k^j$ 接收邻居节点的 ${\boldsymbol{\mu}} _{ji}^L$ 和 ${\boldsymbol{\varLambda}} _{ji}^L$ 来计算自身的后验边缘分布的均值 $\hat{{{\boldsymbol{x}}}}_k^j$ 与方差 ${\boldsymbol{\varLambda}} _k^j$ ，即

至此，完成了DCEKF算法推导，具体流程总结如算法1所示。

1  天基雷达组网目标跟踪DCEKF算法

1.  DCEKF algorithm for target tracking in space-based radar networks

1: 初始化：令 $k = 0$, 对每个雷达节点 $j = 1,2,\cdots, N$，初始化 　 ${ {\hat {{\boldsymbol{x}}}}}_{0{{|}}0}^j$, ${\boldsymbol{P}}_{0|0}^j$；
2: for $k = 1:K$ do
3: 　 for $j = 1,2,\cdots, N$ do
4: 　　计算预测分布： ${ {\hat {{\boldsymbol{x}}}}}_{k + 1|k}^j = {{\boldsymbol{F}}_k}{ {\hat {{\boldsymbol{x}}}}}_k^j$, 　　　　 ${\boldsymbol{P}}_{k + 1|k}^j{{ = }}{{\boldsymbol{F}}_k}{\boldsymbol{P}}_k^j{\boldsymbol{F}}_k^{\mathrm{T}} + {{\boldsymbol{Q}}_k}$；
5: 　　根据式(18)、式(19)计算局部EKF更新 ${\boldsymbol{\mu}} _k^j$和 ${\boldsymbol{P}}_k^j$；
6: 　 end for
7: 　令 $l = 0$；
8: 　 for $j = 1,2,\cdots, N$ do
9: 　　初始化： ${\boldsymbol{\mu}} _{ji}^l,{\boldsymbol{\varLambda}} _{ji}^l,\forall i \in {D_j}$；
10: 　 end for
11: 　　 while $l \le L$ do
12: 　　 $l = l + 1$；
13: 　　 for $j = 1,2,\cdots, N$ do
14: 　　　对节点 $j$的每个邻居节点 $i \in {D_j}$，采用式(22)、式(23) 　　　　　 计算 ${\boldsymbol{\mu}} _{ji}^l$和 ${\boldsymbol{\varLambda}} _{ji}^l$；
15: 　　 end for
16: 　 end while
17: 　 for $j = 1,2,\cdots, N$ do
18: 　　采用式(24)、式(25)计算 ${\hat{\boldsymbol{x}}}_k^j$, ${\boldsymbol{\varLambda}} _k^j$, ${\boldsymbol{P}}_k^j = ({\boldsymbol{\varLambda}}_k^j)^{-1}$；
19: 　 end for
20: end for

下载: 导出CSV 

| 显示表格

3.2   分布式非参数信度传播算法

在DCEKF算法中，将量测方程式(10)中的非线性进行准线性化近似，以达到采用均值和方差表示消息的目的。但线性化近似可能产生较大损失，因此针对因子图模型中的非线性、非高斯概率分布，结合信度传播与粒子滤波，采用分布式非参数信度传播算法(DCNBP)计算因子图中的消息传递。在2.4节基础上加入耦合函数后， $k$ 时刻的因子图如图6所示。

图  6  $k$ 时刻因子图模型

Figure  6.  The factor graph at time $k$

下载: 全尺寸图片 幻灯片

从因子节点 $f_{k-1}$ 传递至变量节点 ${{{\boldsymbol{x}}}}_k^j$ 的消息 ${m_{{f_{k - 1}} \to {{{\boldsymbol{x}}}}_k^j}}$ 为

采用粒子化消息传递算法，根据重要性采样方法构建消息 $m_{{f_{k - 1}} \to {{{\boldsymbol{x}}}}_k^j}$ 的粒子化表示。不失一般性，选择提议分布 $q\left( {{{{\boldsymbol{x}}}}_k^j} \right) = p\left( {{{{\boldsymbol{x}}}}_k^j|{{{\boldsymbol{x}}}}_{k - 1}^j} \right)$ 。记从提议分布中采样的第 $s$ 个粒子状态为 ${{{\boldsymbol{x}}}}_{k,s}^j$ ，则消息 $m_{{f_{k - 1}} \to {{{\boldsymbol{x}}}}_k^j}$ 的粒子化表示为

其中，

将权重 $w_{{f_{k,s}}}^j$ 归一化，记为 $\bar{w}_{{f_{k,s}}}^j$ 。

从因子节点 ${g_{ij}}$ 传递至变量节点 ${{{\boldsymbol{x}}}}_k^j$ 的消息 ${m_{{g_{ij}} \to {{{\boldsymbol{x}}}}_k^j}}$ 可以写成

重要性采样的提议分布选择耦合函数 ${g_{ij}}\left( {{{{\boldsymbol{x}}}}_k^i,{{{\boldsymbol{x}}}}_k^j} \right)$ ，则消息 ${m_{{g_{ij}} \to {{{\boldsymbol{x}}}}_k^j}}$ 的粒子化表示为

其中，

从因子节点 $h_k^j$ 传递至变量节点 ${{{\boldsymbol{x}}}}_k^j$ 的消息为 ${m_{h_k^j \to {{{\boldsymbol{x}}}}_k^j}}$ 。重要性采样的提议分布可选择似然函数 $p({{\boldsymbol{y}}}_k^j | {{{\boldsymbol{x}}}}_k^j )$ ，但由于观测方程中 ${{\boldsymbol{y}}}_k^j$ 仅与目标的位置 ${{\boldsymbol{o}}}_k^j = [ {x_k^j,y_k^j,z_k^j} ]$ 有关，因此提议分布 $p({{\boldsymbol{y}}}_k^j | {{{\boldsymbol{x}}}}_k^j )$ 可以简化为 $p({{\boldsymbol{y}}}_k^j | {{\boldsymbol{o}}}_k^j )$ ，则目标位置 ${{\boldsymbol{o}}}_k^j$ 可以表示为

其中， ${g^{ - 1}}\left( \cdot \right)$ 表示从雷达阵面坐标系到ECEF坐标系的转换函数，且 $\phi_z \sim {\cal{U}}\left( {0,2\pi } \right)$ 。因此，从因子节点 $h_k^j$ 传递至变量节点 ${{{\boldsymbol{x}}}}_k^j$ 的消息为 $m_{h_k^j \to {{{\boldsymbol{x}}}}_k^j}$ ，其粒子化形式为

其中，

且

式(37)中， ${p_{\Delta r_k^n}}\left( {\Delta r_k^n} \right)$ 和 ${p_{\Delta \theta_k^n}}\left( {\Delta \theta_k^n} \right)$ 表示观测噪声的概率分布函数。由于 ${{\boldsymbol{o}}}_k^j$ 的粒子从式(37)中直接采样得到，因此， $S$ 个粒子的权重相等，均为 $1/S$ 。

从变量节点 ${{{\boldsymbol{x}}}}_k^j$ 传递给因子节点 $f_k$ 的消息可以写成

从变量节点 ${{{\boldsymbol{x}}}}_k^j$ 传递给因子节点 $g_{ji}$ 的消息可以写成

对于 ${m_{{{{\boldsymbol{x}}}}_k^j \to {g_{ji}}}}$ 和 ${m_{{{{\boldsymbol{x}}}}_k^j \to f_{k}}}$ 这种多个消息乘积的形式，采用核密度估计^[18]的方式进行粒子化采样。不失一般性，设 ${p\left({{\boldsymbol{x}}}\right)}$ 表示 ${D}$ 个输入消息的乘积

为了简化描述，式(40)去掉关于时间变量 $k$ 和节点 $j$ 的表示。第 $i$ 个输入的消息 ${m_i}\left( {{\boldsymbol{x}}} \right)$ 可表示为 ${m_i}\left( {{\boldsymbol{x}}} \right) = \left\{ {{{{\boldsymbol{x}}}_{i,s}},{w_{i,s}}} \right\}_{s = 1}^S$ ，因此，输入消息 ${m_i}\left( {{\boldsymbol{x}}} \right)$ 可通过核密度估计的方法非参数化表示为

其中， ${\cal{N}} \left( {{{\boldsymbol{x}}};{{{\boldsymbol{x}}}_{i,s}},{{\boldsymbol{P}}_i}} \right)$ 表示归一化的高斯密度函数，均值和协方差分别为 ${{\boldsymbol{x}}}_{i,s}$ 和 ${{\boldsymbol{P}}_i}$ 。式(41)中， $S$ 个高斯密度函数选取相同的协方差

其中，

因此， $D$ 个输入消息的乘积可以表示为

易知， $D$ 个高斯密度函数的乘积仍然是一个高斯密度函数

其中，

因此，式(45)所示的高斯密度函数 ${\cal{N}}\left( {{\boldsymbol{x}};{\bar{\boldsymbol{\mu}}} ,{\bar{\boldsymbol{P}}}} \right)$ 的权重 $\bar w$ 可以表示为

其中， $\left\{ {{w_i}} \right\}_{i = 1}^D$ 是 $D$ 个输入消息高斯混合式中高斯分量的权重。

由式(40)可知， $p\left( {{{\boldsymbol{x}}}} \right)$ 表示 $D$ 个高斯混合式的乘积，每个高斯混合又含有 $S$ 个高斯分量，则 $p\left( {{{\boldsymbol{x}}}} \right)$ 一共含有 $S^D$ 个高斯分量，直接采样的计算复杂度为 $O\left( {{S^D}} \right)$ 。因此，对 $D$ 个高斯混合式的乘积采用重要性采样^[18]，引入辅助变量用来表示抽取样本的高斯分量

且 ${\zeta _{1:D}} = \left\{ {{\zeta_1},{\zeta_2}, \cdots, {\zeta_D}} \right\}$ 。离散随机变量 ${\zeta_i} = {l_i}$ 表示从第 $i$ 个高斯混合式中的第 $l_i$ 个高斯分量中抽取一个样本，则式(40)可写成

令重要性采样的建议性分布为 $p\left( {{\zeta_{1:D}}} \right)$ ，且

从建议性分布 $p\left( {{\zeta_{1:D}}} \right)$ 中获得辅助变量 ${\zeta_{1:D}} \in {\left\{ {1,2, \cdots ,N} \right\}^D}$ 的一组样本 $l_{1:D}$ 。然后通过式(46)和式(47)计算出式(45)所示的 $D$ 个加权高斯密度函数乘积 $\displaystyle\prod\nolimits_{i = 1}^D {{w_{i,{l_i}}} {\cal{N}} \left( {{{\boldsymbol{x}}};{{{\boldsymbol{x}}}_{i,{l_i}}}, {{\boldsymbol{P}}_i}} \right)}$ 的均值 ${\bar{\boldsymbol{\mu}}}$ 、协方差 ${\bar{\boldsymbol{P}}}$ 和权重 $\bar w$ 。再从高斯密度函数 ${\cal{N}} \left( {{{\boldsymbol{x}}};{\bar{\boldsymbol{\mu}}} ,{\bar{\boldsymbol{P}}} } \right)$ 采样一个粒子 ${\hat{{\boldsymbol{x}}}} \sim {\cal{N}} \left( {{{\boldsymbol{x}}};{\bar{\boldsymbol{\mu}}} ,{\bar{\boldsymbol{P}}} } \right)$ ，计算权重 $\hat w = \bar w/\displaystyle\prod\nolimits _{i = 1}^D{w_{i,{l_i}}}$ ，对重要性权重 $\hat w$ 做归一化处理，使权重和为1。因此， $D$ 个输入消息乘积的粒子化形式为 $\left\{ {{{\hat{{\boldsymbol{x}}}}_s},{{\hat w}_s}} \right\}_{s = 1}^S$ ，该采样方法的计算复杂度仅为 $O\left( {DS} \right)$ 。基于重要性采样的消息乘积粒子化算法如算法2所示。

2  基于重要性采样的消息乘积粒子化算法

2.  Message product particleization algorithm based on importance sampling

Require: $D$个输入消息，且每个输入消息 ${m_i}\left( {{\boldsymbol{x}}} \right)$的粒子化形式 　为 $\left\{ {{{{\boldsymbol{x}}}_{i,s}},{w_{i,s}}} \right\}_{s = 1}^S$；
Ensure: $D$个输入消息乘积的粒子化形式 $\left\{ {{{{\hat{{\boldsymbol{x}}}}}_s},{{\hat w}_s}} \right\}_{s = 1}^S$；
1: for $i = 1,2, \cdots, D$ do
2: 　构造输入消息 ${m_i}\left( {{\boldsymbol{x}}} \right)$的高斯混合形式： 　　 ${m_i}\left( {{\boldsymbol{x}}} \right) = \displaystyle\sum\limits_{s = 1}^S {{w_{i,s}}{\cal{N}}\left( {{{\boldsymbol{x}}};{{{\boldsymbol{x}}}_{i,s}},{{\boldsymbol{P}}_i}} \right)}$；
3: 　 for $s = 1,2, \cdots, S$ do
4: 　　设置均值： ${{\boldsymbol{\mu}} _{i,s}} = {{\boldsymbol{x}}}_{i,s}$；
5: 　　设置协方差：
${{\boldsymbol{P}}_i} = {S^{\textstyle\frac{1}{6}}}\displaystyle\sum\limits_{s = 1}^S {\displaystyle\sum\limits_{l = 1}^S{{w_{i,s}}{w_{i,l}}}} \left({{{{\boldsymbol{x}}}_{i,s}} - {{{\bar{{\boldsymbol{x}}}}}_i}} \right) {\left( {{{\boldsymbol{x}}_{i,l}} - \bar{\boldsymbol{x}}_i} \right)^{\mathrm{T}}}$；
6: 　 end for
7: end for
8: for $s = 1,2, \cdots, S$ do
9: 　 for $i = 1,2, \cdots, D$ do
10: 　　从建议性分布为 $p\left( {{\zeta_{1:D}}} \right)$中获得辅助变量 　　　　 ${\zeta_{1:D}} \in {\left\{ {1,2, \cdots, S} \right\}^D}$的一组样本 $l_{1:D}$；
11: 　　根据式(46)和式(47)计算 $\prod\nolimits_{i = 1}^D {{w_{i,{l_i}}}{\cal{N}}\left( {{\boldsymbol{x}};{{\boldsymbol{\mu}}_{i,{l_i}}},{{\boldsymbol{P}}_i}} \right)}$ 的 　　　　 均值 ${\bar {\boldsymbol{\mu}}}$、协方差 ${\bar {\boldsymbol{P}} }$和权重 $\bar w$；
12: 　　采样一个粒子 ${\hat {\boldsymbol{x}}} \sim {\cal{N}}\left( {{\boldsymbol{x}};{\bar {\boldsymbol{\mu}}} ,{\bar {\boldsymbol{P}} }} \right)$；
13: 　　计算权重 $\hat w = \bar w/\prod {_{i = 1}^D{w_{i,{l_i}}}}$；
14: 　 end for
15: end for
16: 归一化权重得到 $D$个输入消息乘积的粒子化形式 　 $\left\{ {{{{\hat {\boldsymbol{x}}}}_s},{{\hat w}_s}} \right\}_{s = 1}^S$。

下载: 导出CSV 

| 显示表格

注意，粒子化消息 ${m_{{{\boldsymbol{x}}}_k^j \to {f_k}}}( {{{\boldsymbol{x}}}_k^j} )$ 和 ${m_{{{\boldsymbol{x}}}_k^j \to {g_{ji}}}} ( {{{\boldsymbol{x}}}_k^j} )$ 时，输入消息 ${m_{h_k^j \to {{\boldsymbol{x}}}_k^j}}$ 仅是关于 ${\boldsymbol{o}}_k^j$ 的函数，因此将变量节点 ${{\boldsymbol{x}}}_k^j$ 分为前3项 ${\boldsymbol{o}}_k^j = {[ {x_k^j,y_k^j,z_k^j} ]^{\mathrm{T}}}$ 和后3项 ${\boldsymbol{v}}_k^j = {[ {\dot x_k^j,\dot y_k^j,\dot z_k^j} ]^{\mathrm{T}}}$ 。然后，分别对消息 ${m_{{{\boldsymbol{x}}}_k^j \to {f_k}}} ( {{\boldsymbol{o}}_k^j} )$ , ${m_{{{\boldsymbol{x}}}_k^j \to {g_{ji}}}}( {{\boldsymbol{o}}_k^j} )$ 和 ${m_{{{\boldsymbol{x}}}_k^j \to {g_{ji}}}}( {{\boldsymbol{v}}_k^j} )$ 进行粒子化，最后将粒子 ${\boldsymbol{o}}_{k,s}^j$ 和粒子 ${\boldsymbol{v}}_{k,s}^j$ 进行合并，构成目标状态 ${\boldsymbol{x}}_k^j$ 的后验概率分布的粒子化形式 $p( {{{\boldsymbol{x}}}_k^j} ) = \{ {{\boldsymbol{x}}}_{k,s}^j, w_{k,s}^j \}_{s = 1}^S$ 。天基雷达组网目标定位NBP算法流程总结如算法3所示。

3  天基雷达组网目标跟踪DCNBP算法

3.  DCNBP algorithm for target tracking in space-based radar networks

1: 初始化：令 $k = 0$，对每个雷达节点 $j = 1,2,\cdots, N$，初始化 　信度 ${b_0}( {{{\boldsymbol{x}}}_0^j} ) = {p_0}( {{{\boldsymbol{x}}}_0^j} )$；
2: for $k = 1:K$ do
3: 　 for $j = 1,2,\cdots, N$ do
4: 　　根据式(27)—式(29)建立消息 ${m_{{f_{k - 1}} \to {{\boldsymbol{x}}}_k^j}}$的粒子化表示 　　　 $\{ {{{\boldsymbol{x}}}_{k,s}^j,\bar w_{{f_k},s}^j} \}_{s = 1}^S$；
5: 　　根据式(34)—式(37)建立消息 ${m_{h_k^j \to {{\boldsymbol{x}}}_k^j}}$的粒子化表示 　　　 $\{ {{\boldsymbol{x}}_{k,s}^j,w_{{h_k},s}^j} \}_{s = 1}^S$；
6: 　　根据算法2建立消息乘积 ${m_{{f_{k - 1}} \to {{\boldsymbol{x}}}_k^j}}{m_{h_k^j \to {{\boldsymbol{x}}}_k^j}}$的粒子化 　　　表示；
7: 　 end for
8: 　令 $l = 0$；
9: 　 for $j = 1,2,\cdots, N$ do
10: 　　粒子初始化 $m_{{g_{ji}} \to {{\boldsymbol{x}}}_k^j}^0 = \{ {{{\boldsymbol{x}}}_{ji,s}^0,w_{{g_{ji}},s}^0} \}_{s = 1}^S$；
11: 　 end for
12: 　 while $l \le L$ do
13: 　　 $l = l + 1$； 　14: 　　采用算法2计算 $m_{{{\boldsymbol{x}}}_k^j \to {g_{ji}}}^l = {m_{h_k^j \to {{\boldsymbol{x}}}_k^j}}\prod\limits_{u \in {D_j}\backslash i} {m_{{g_{uj}} \to {{\boldsymbol{x}}}_k^j}^{l - 1}}$， 　　　　粒子化形式为 $m_{{{\boldsymbol{x}}}_k^j \to {g_{ji}}}^l = \{ {\bar {\boldsymbol {x}}_{ji,s}^l,\bar w_{{g_{ji}},s}^l} \}_{s = 1}^S$；
15: 　　根据耦合关系更新 $m_{{g_{ji}} \to {{\boldsymbol{x}}}_k^j}^l = \{ {{{\boldsymbol{x}}}_{ji,s}^l,w_{{g_{ji}},s}^l} \}_{s = 1}^S$；
16: 　 end while 　17: 　根据算法2计算消息乘积 ${m_{{f_{k - 1}} \to {{\boldsymbol{x}}}_k^j}}{m_{h_k^j \to {{\boldsymbol{x}}}_k^j}}\prod\limits_{u \in {D_j}} {m_{{g_{uj}} \to {{\boldsymbol{x}}}_k^j}^L}$ 　　　 的粒子化形式以表示后验分布 $p( {{{\boldsymbol{x}}}_k^j} )$，根据后验分布得到 　　　 节点的估计值 $\hat {\boldsymbol {x}}_{k}^j$；
18: end for

下载: 导出CSV 

| 显示表格

综上，以线图拓扑下的天基雷达网络为例，基于DCNBP的天基雷达组网分布式目标融合跟踪结构图如图7所示。每个雷达节点包含利用本地量测的目标跟踪模块以及利用其他节点估计值的融合模块。初始化后，在 $k$ 时刻，每部雷达在目标跟踪模块基于本地量测进行本地跟踪。本地跟踪后，每部雷达再根据相邻雷达传递的本地跟踪值与本地量测进行消息融合更新，经过多次迭代更新，实现多部雷达目标跟踪的一致性融合。

图  7  DCNBP算法分布式实现结构图

Figure  7.  Block diagram of the distributed implementation of DCNBP

下载: 全尺寸图片 幻灯片

3.3   算法复杂度分析

本节通过计算执行的基本运算，比较分析DCEKF算法与DCNBP算法的计算复杂度。记目标的局部状态向量大小为 $s$ ，局部量测向量大小为 $M$ ，DCNBP算法中粒子数为 $S$ ，消息迭代次数为 $L$ 。

对于DCEKF算法，根据3.1节的算法1，利用式(18)、式(19)计算局部 ${\boldsymbol{\mu}} _k^j$ 和 ${\boldsymbol{P}}_k^j$ ，计算复杂度分别为 $O(M^3+sM^2+s^3)$ 和 $O(M^3+sM^2+s^2M+s^3)$ ；采用式(22)、式(23)通过 $L$ 次迭代计算节点 $j$ 的所有邻居节点的 ${\boldsymbol{\mu}} _{ji}^L$ 和 ${\boldsymbol{\varLambda}} _{ji}^L$ ，计算复杂度均为 $O(LM^3+ sLM^2+s^2LM+s^3L)$ ；通过式(24)、式(25)计算节点 $j$ 的 $\hat {\boldsymbol {x}}_k^j$ 和 ${\boldsymbol{\varLambda}} _k^j$ ，计算复杂度均为 $O(s^3)$ ；最后，计算下一时刻的预测分布 $\hat {\boldsymbol {x}}_{k + 1|k}^j$ 和 ${\boldsymbol{\varLambda}}_{k + 1|k}^j$ ，计算复杂度为 $O(s^2)$ 和 $O(s^3)$ 。综上，DCEKF算法的计算复杂度为 $O(LM^3+sLM^2+s^2LM+s^3L)$ 。

对于DCNBP算法，每个节点的消息均为单独计算，因此其状态向量和量测向量维度为 $s$ 和 $L$ 。每个时刻，首先采用式(28)和式(29)生成带权值粒子进行粒子初始化，生成粒子和权值的计算复杂度分别为 $s^2S$ 和 $s^2+sS+s$ ，从而该步每个节点的计算复杂度为 $O(s^2S)$ 。然后，采用算法2计算算法3中的所有消息乘积的粒子化形式。对于算法2，记节点的输入消息个数为 $D$ ，计算所有高斯混合消息的均值与协方差的复杂度为 $sS$ 和 ${s^2}DS^2$ 。根据设置辅助变量采样粒子，式(46)和式(47)的计算复杂度分别为 $s^3D+s^2D+s^2$ 和 $sD+s^2+2s$ 。因此， $D$ 个输入消息乘积的粒子化形式的计算复杂度为 $DS(s^3D+ s^2D+sD+2s^2+2s)$ 。从而，算法2的计算复杂度为 $O(s^2DS^2+s^3D^2S)$ 。因此，经过多次迭代计算每个节点消息的计算复杂度为 $O(s^2LDS^2+s^3LD^2S)$ 。最后通过算法2计算每个节点的后验分布。因此DCNBP算法的计算复杂度为 $O(s^2LDS^2+s^3LD^2S)$ 。

本文选择3颗天基雷达测量空中目标的径向距和方位角，单颗雷达的局部量测向量维度 $M$ 与消息输入个数大小 $D$ 均较小。因此，DCEKF和DCNBP算法的计算复杂度分别为 $O(s^3L)$ 和 $O(s^2S^2L+s^3LS)$ 。可以看出，DCNBP算法的粒子数越多，计算复杂度相较DCEKF算法越大。

3.4   滤波初始状态的选取与俯仰角先验信息

在实际跟踪中，滤波器初值的选取会对跟踪性能造成极大的影响，空中目标俯仰角误差过大是天基雷达目标定位的一个关键问题。

雷达得到的俯仰角误差越大，得到的滤波器初值与真实值偏差就越大。本文利用几何法^[19]先求解目标高度 $h$ ，再结合地球半径 ${r_e}$ 、径向距 ${r_z}$ 和方位角 ${\theta_z}$ 估计目标俯仰角 $\phi_z$ ，由雷达与目标空间位置关系可知

其中， $r_s$ 为雷达到地心的距离， $(x_c,y_c,z_c)$ 为目标在LVLH坐标系下的坐标。实际天基雷达天线较多安装在对地面，即阵面姿态角中仅滚动角变化较多，为简化推导，仅假设滚动角 ${{\gamma _c}}$ 非零，则

其中， $(x_a,y_a,z_a)$ 为目标在雷达阵面坐标系下的坐标，将式(52)代入到式(51)中，得

整理得

再将式(6)代入到式(54)，整理得

其中，

解方程，得

其中，

通过上述方法可估计出雷达量测的俯仰角信息，并估计出目标高度作为滤波算法初值。一方面提高了滤波算法初值精度，另一方面对于单部雷达将俯仰角估计值初值加入到量测方程中，可以部分解决单部雷达因为俯仰角缺失或误差太大无法跟踪目标的问题，但目标高度估计误差仍然较大。

4.   仿真与分析

本节通过仿真验证DCNBP算法的性能，与几何法、传统集中式EKF (CEKF)、DCEKF算法、集中式粒子滤波(CPF)算法进行对比，基于高度求解结果对俯仰角进行估计并将其加入到量测，仿真验证高度估计误差对目标定位精度的影响。

4.1   仿真场景设置

考虑一个空中匀速直线飞行的目标，高度为9 km，采用3颗天基探测雷达对此空中目标进行跟踪。3颗卫星轨道六根数，即轨道半长轴 $a$ ，偏心率 $e$ ，轨道倾角 $\eta$ ，近地点幅角 $\omega$ ，升交点赤经 $\varOmega$ , 平近点角 $M$ 参数设置如下。卫星1：6978.14 km, 1.369e–15, $20^{\circ}$ , $0^{\circ}$ , $0^{\circ}$ , $50^{\circ}$ ；卫星2：7078.14 km, 1.369e–15, $20^{\circ}$ , $0^{\circ}$ , $0^{\circ}$ , $60^{\circ}$ ；卫星3：7178.14 km, 1.369e–15, $20^{\circ}$ , $0^{\circ}$ , $0^{\circ}$ , $70^{\circ}$ ；3颗卫星姿态角均为 $0^\circ$ 。3部天基雷达阵面姿态角均为滚动角 $30^\circ$ 、俯仰角 $0^\circ$ 、偏航角 $0^\circ$ 。雷达量测噪声标准差相同，根据典型的雷达与目标参数(脉冲宽度、方位向波束宽度等)取值，设置为径向距0.09 km，方位角 $0.03^\circ$ 。

设置迭代次数 $L = 4$ ；系统噪声协方差均为 $Q = \text{diag}\{10^{-3},\;10^{-3},\; 10^{-3}, \;10^{-3}/T, \;10^{-3}/T, \;10^{-3}/T\}$ ， $T = 15\;{\mathrm{s}}$ ；仿真得到目标在ECEF坐标系下的位置，并将其转换到大地坐标系，得到目标高度估计。算法性能评估指标为空中目标高度的均方根误差，即

其中， $N$ 表示蒙特卡罗仿真次数，取 $N = 100$ ， $\tilde{h}_k$ 和 $h_k$ 分别表示 $k$ 时刻目标估计高度和真实高度。在CPF算法中，令粒子数 $S = 800$ 。

4.2   结果分析

首先令耦合参数 $\kappa = 500$ ，粒子数 $S = 800$ 。CEKF、DCEKF、DCNBP、几何法估计以及CPF的目标高度RMSE结果在图8中给出，从图中可以看出，几何法的目标高度估计误差最大，在1.5～3.0 km，这是因为仅采用了量测信息。集中式算法由于能够利用全局信息，一般情况下，相比于分布式估计方法，具有更高的估计精度。因此，可以看到，传统CEKF算法和CPF算法的精度最高，目标高度估计的RMSE为0.50 km和0.45 km。但实际工作时，集中式算法需要将所有量测信息与状态信息发送到一颗卫星上进行数据处理，对该单星的通信负载与计算速率要求高，并且如果该星发生故障后，可能会导致整个系统失效，从而，对于天基雷达应用，包括CEKF, CPF算法在内的传统集中式滤波器鲁棒性与灵活性较差。而DCEKF, DCNBP等分布式滤波算法能更好地克服此类问题。图8中，DCEKF算法、DCNBP算法高度估计RMSE分别为0.85 km和0.55 km。相比于DCEKF算法，DCNBP算法估计精度提升了35.3%，同时曲线波动更小，可见DCNBP算法不仅具有较高的估计精度，也具有较好的灵活性和鲁棒性。本文未实现分布式PF算法，理论上，与DCNBP算法相比，分布式PF算法具有更高的估计精度，但计算量也更大；天基平台计算资源较为稀缺，有必要在目标高度估计精度与计算量之间做出权衡。

图  8  空中目标高度估计RMSE

Figure  8.  RMSE of aerial target altitude estimation

下载: 全尺寸图片 幻灯片

DCNBP算法中，耦合参数 $\kappa$ 和采样粒子数 $S$ 为关键参数，能够影响算法估计性能。图9和图10分别为不同 $\kappa$ 和 $S$ 下DCNBP算法的目标高度估计结果，表1中列出了不同参数下算法估计的RMSE平均值。可以看出，随着 $\kappa$ 和 $S$ 增大，DCNBP算法高度估计的RMSE变小，精度也越高。这是由于随着 $\kappa$ 的增加，相邻变量的相关性越强，邻居节点量测信息利用越充分，估计精度越好。而根据3.1节，当 $\kappa \to \infty$ 时，局部变量与全局变量分布相同，高度估计精度并不会随着 $\kappa$ 的增加一直提高，其最终会趋于某一定值。同时当粒子数 $S$ 越多，对非线性分布的近似误差越小，最终的估计精度也更高，但 $S$ 的增加也会提高算法的计算复杂度与通信负载，在应用当中需要根据实际情况进行权衡。

图  9  空中目标高度估计RMSE(粒子数 $S = 400$ )

Figure  9.  RMSE of aerial target altitude estimation (particle number $S = 400$ )

下载: 全尺寸图片 幻灯片

图  10  空中目标高度估计RMSE(粒子数 $S = 800$ )

Figure  10.  RMSE of aerial target altitude estimation (particle number $S = 800$ )

下载: 全尺寸图片 幻灯片

表  1  不同参数下DCNBP算法估计的RMSE平均值(km)

Table  1.  Average RMSE obtained by DCNBP with different parameters (km)

耦合参数( $\kappa$)	$S = 400$	$S = 800$
$\kappa = 50$	0.8196	0.6596
$\kappa = 500$	0.7511	0.5506
$\kappa = 5000$	0.6389	0.4661

下载: 导出CSV 

| 显示表格

在天基雷达跟踪目标过程中，雷达测量误差与雷达到目标距离相关，因此分别取雷达对空中目标量测标准差为：径向距0.07 km，方位角 $0.02^\circ$ ；径向距0.09 km，方位角 $0.03^\circ$ ；径向距0.11 km，方位角 $0.04^\circ$ 。图11所示为不同雷达测量误差下DCNBP算法空中目标高度估计结果。由图可知，当目标与雷达距离越近时，雷达量测误差越小，目标最终的高度估计结果越精确。

图  11  不同测量误差空中目标高度估计RMSE

Figure  11.  RMSE of aerial target height estimation with different measurement errors

下载: 全尺寸图片 幻灯片

前文中推导了基于高度估计的雷达俯仰角计算方法，该计算结果可以被认为是一个雷达伪测量值，高度估计越准确，该俯仰角伪量测值误差就越小。为了进一步提高目标定位精度，将该伪测量值当作量测信息加入到滤波中，得到雷达径向距、方位角与俯仰角信息，对空中目标进行定位估计。图12为不同高度估计误差下空中目标定位结果的RMSE。可以看出，高度估计误差越大，俯仰角伪量测误差越大，目标定位结果的RMSE也就越大；当高度估计误差很大时，可等效为俯仰角信息缺失，此时，目标定位RMSE最大；随着高度估计误差的减小，俯仰角伪量测误差逐渐减小，目标定位RMSE也逐渐减小，可知对目标高度估计越准确，加入俯仰角估计后的目标位置估计精度越高。由此可得，优先对目标高度进行估计进而得到雷达俯仰角伪测量值能够弥补雷达在实际情况中俯仰角量测误差大的缺点，进一步提高目标定位精度。

图  12  不同高度误差空中目标定位估计RMSE

Figure  12.  RMSE of aerial target position estimation with different altitude errors

下载: 全尺寸图片 幻灯片

5.   结语

本文针对天基探测雷达空中目标定位时面临的俯仰角信息缺失、量测非线性等问题，提出了一种基于一致性的分布式天基雷达组网空中目标高度估计与定位方法。首先建立了空中目标运动模型与天基雷达量测模型，并基于因子图对传感器目标跟踪进行建模。考虑到雷达量测非线性，推导了基于因子图消息传递的DCEKF算法；针对DCEKF对高非线性问题近似误差大的问题，提出了基于采样的DCNBP算法。仿真验证表明，相较于DCEKFF、几何法，DCNBP算法对目标高度估计精度更高。同时，分析了DCNBP算法参数对估计性能的影响。最后，根据得到的目标高度值估计出目标的俯仰角，并将该估计值当作伪测量估计目标位置。仿真结果显示，目标高度估计越准确，加入俯仰角估计后的目标位置估计精度越高。