基于动态规划和双变量容积卡尔曼滤波的空基外辐射源定位算法

卢雨，衣晓

(海军航空大学，山东 烟台 264001)

0 引言

空基外辐射源定位(APCL)是指利用空基外辐射源对目标进行定位跟踪的一种技术[1-6]。区别于地基外辐射源如数字广播电视[2-3]、移动手机信号[4]等，以及天基外辐射源如通信卫星[5]、全球定位导航信号[6]等，空基外辐射源信号则是指由空中移动平台所辐射的电磁信号。这种外辐射源信号不具备地基、天基两种信号所拥有的稳定性和可靠性，而且外辐射源信号是非合作的，空基外辐射源平台的运动状态未知且难测，这就极大地影响了系统对目标的定位跟踪精度。因此，如何在复杂的杂波环境下提取目标状态，并降低外辐射源不确定性对系统定位跟踪精度的影响是提升APCL系统实用效能的重要环节。

针对无源定位系统中存在传感器状态不确定性的影响，文献[7]提出一种迭代最大似然概率数据关联(IML-PDA)算法，该算法采用迭代预测-更新框架，有效提高了目标的定位精度。针对外辐射源定位系统中存在外辐射源状态不确定的影响，文献[8]给出了全面详细的应对方案，该方案采用联合状态估计的方法，提出一种递推的航迹起始/保持算法，将外辐射源状态和目标状态同时进行更新，当目标和外辐射状态相互独立时，可以取得较好的效果。实际应用背景下，过程和量测噪声并非服从高斯分布且统计特性不精确，针对这种情况，文献[9]提出一种基于最大熵方法的变分贝叶斯自适应卡尔曼滤波，进而提高了系统的滤波精度和稳定性。当环境中存在强杂波干扰时，则需要考虑目标是否存在，如何提取出目标状态，文献[10]提出了一种基于序贯蒙特卡洛与概率假设密度滤波相结合的目标跟踪方法，在主动分布式声纳系统的目标定位应用中取得了较好的效果，但是该方法在被动跟踪场景中的应用效果还有待研究。文献[11]基于极值理论确定航迹回溯的阈值，利用误差传播理论构建动态规划(DP)的状态转移范围，并通过幅值积累的方式提取出了目标的状态，但是采用的状态扩维(SA)方法计算量较大，且没有考虑目标和外辐射源状态的相互独立性对目标和外辐射源状态估计造成的影响。针对APCL系统在实际应用中可能出现的不可观测问题，文献[12]基于条件数定义的可观测度对其进行了详细的分析。

考虑到APCL系统的应用环境相对复杂，目标回波功率相对微弱，相比于外辐射源的直达波而言，容易出现被杂波淹没的问题。本文研究的主要问题就是如何提取杂波中的目标量测，并降低外辐射源状态不确定性的影响，实现对目标和外辐射源的状态估计。针对杂波环境下的空基外辐射源定位问题，本文提出一种基于DP和双变量容积卡尔曼滤波(DP-BVCKF)的定位跟踪算法。首先基于误差传播理论，确定出目标的状态转移范围；其次，选取回波信号幅度值作为 DP的值函数，在每个候选测量值的状态转移范围中寻找相邻时刻之间的最大幅度值进行递推积累；再次，基于杂波的统计特性和极值理论来解算DP的目标提取阈值，逆向提取来自目标的测量值序列[11]；然后利用双变量容积卡尔曼滤波(BVCKF)算法同时对目标和外辐射源进行状态估计，从原理上降低外辐射源状态不确定性对目标状态估计造成的影响。

1 问题描述

Xk+1=FkXk+Gkuk+Vk，

(1)

强杂波环境下的APCL系统测量值将会包含源于目标、外辐射源以及杂波的信息。假设k时刻观测站接收到Nk个量测数据，其数据集合为

{zk(i)=

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

幅度值Ak(i)的测量方程为

(10)

(11)

2 DP-BVCKF算法描述

当目标回波较弱时，难以区分杂波测量值和目标测量值，再加上外辐射源状态不确定性的影响，直接采用一般的滤波算法很容易形成虚假目标航迹。针对APCL系统的上述问题，本文提出DP-BVCKF算法，该算法主要分两个步骤：一是对一组已获取的含杂波量测值，进行DP[16]处理，逆向提取目标的量测序列；二是采取BVCKF算法，同时对目标和外辐射源进行状态估计，降低外辐射源状态不确定性对目标状态估计的影响。具体算法流程如图2所示。图2中k+1|k、Px,k+1|k分别为目标的状态向量和协方差，k+1|k、Py,k+1|k分别为外辐射源的状态向量和协方差，k+1|k为量测的一步预测值，h为量测转移矩阵，Kx、Ky分别为目标和外辐射源的卡尔曼滤波增益。

图2 DP-BVCKF算法流程图

2.1 构建DP的状态转移范围

在APCL系统的非线性测量空间中，当前时刻的目标在下一时刻可能出现的状态集合记为目标的状态转移范围记为Ω(zk+1)。考虑到外辐射源不确定性的影响，Ω(zk+1)的大小主要取决于θk(i)、φk(i)以及rk(i)，其范围参数可以由目标方位角和到达时间差的变化量来表示，即Ω(zk+1)=[Δθk,Δrk]。因此，根据APCL系统的非线性测量信息可以构建目标的状态转移范围[11]为

(Δθk)2=(Δxk)2(∂θk/∂xk)2|xk=α(zk(i))+

(Δyk)2(∂θk/∂yk)2|yk=β(zk(i))，

(12)

(Δrk)2=(Δxk)2(∂rk/∂xk)2|xk=α(zk(i))+

(Δyk)2(∂rk/∂yk)2|yk=β(zk(i))+

(Δxt,k)2(∂rk/∂xt,k)2|xt,k=t,k+

(Δyt,k)2(∂rk/∂yt,k)2|yt,k=t,k，

(13)

α(zk(i))=

(14)

β(zk(i))=

(15)

2.2 幅值积累和航迹逆向提取

从杂波中提取目标量测数据的核心在于杂波的幅值大小取决于其统计特性，而目标的幅值受多种因素影响，在利用DP方法选取幅值Ak(i)作为值函数[16]对相邻时刻中的量测数据进行幅值积累时，目标的状态序列所积累的值函数往往是最大的。根据极值理论，可以求取一定检测概率下的杂波幅度极值，并将其作为幅值积累的门限阈值，逆向提取目标的状态序列。记DP的值函数为I(zk(j))，j表示用于值函数累积的第j个测量值。目标状态序列逆向提取的具体步骤如下：

1)确定值函数积累的门限阈值。根据极值理论，杂波幅值的极大值服从Gumbel分布，其函数表示式为

(16)

式中：FA(·)为Gumbel分布函数；PFT为虚警率；a和b为FA(·)参数，可通过(17)式、(18)式[17]迭代估计得到：

(17)

(18)

2)幅值积累初始化。当k=1时，遍历当前时刻的所有量测数据，令值函数I(z1(i))=A1(i)，记η0(z1(i))=0为目标状态序列。

3)幅值积累，阈值判断，终止积累。当k>1时，遍历k时刻的所有量测数据，进行幅值积累，过程为

(19)

(20)

式中：Ωk-1(zk(i))为第k-1时刻量测数据可能转移到zk(i)的状态转移范围；ηk-1(zk(i))为转移到zk(i)的第k-1时刻测量值。再根据步骤1求得的门限阈值S，提取出满足判定条件的所有状态序列集合M(i)，判定条件为

{M(i)}={z1∶M(i)∶I(zM(i))>S},

(21)

式中：z1∶M(i)表示满足判定条件的M个状态序列。

2.3 基于BVCKF算法的状态估计

BVCKF算法的核心思想是将量测一步预测过程中对单变量的多维积分扩展为对双变量的二重多维积分，并在互协方差计算过程中对单变量的概率密度函数求解转换为对双变量联合概率密度函数的边缘概率密度函数的求解。由于非线性函数的多维积分不易求解[18]，该算法采用了球面-径向容积准则[19-20]来对非线性函数的多维积分进行近似计算，既降低了计算量，又保证了较高的近似精度，具体算法流程描述如下：

(22)

(23)

(24)

式中：n为状态向量的维数；ξi为初始容积点集；[1]i表示第i个初始容积点，初始容积点集[1]的具体形式为

[1]=

2)状态一步预测。首先计算经过系统状态方程f(·)传递的容积点：

(25)

(26)

k+1|k(k+1|k)T+Γk，

(27)

式中：ωi为权重，ωi=1/(2n)。同理可得Y在k+1时刻的状态一步预测k+1|k和一步预测协方差矩阵Py,k+1|k.

3)重新构造容积点。根据球面-径向容积准则，当得到状态一步预测k+1|k及其对应的协方差矩阵后，需要重新构造容积点，以便更好地近似变量的概率密度函数：

(28)

(29)

4)双变量联合的量测一步预测。将对单一变量的多重积分扩展为对双变量的二重积分，再对二重积分进行数值近似，得到双变量联合的量测一步预测及其协方差为

(30)

(31)

(32)

5)增益的计算。首先计算状态变量X的互协方差：

(33)

(34)

再计算状态变量的卡尔曼增益：

(35)

同理可得状态变量Y对应的卡尔曼增益Ky.

6)状态更新。利用计算得到的卡尔曼增益和最新时刻的量测值就可以对状态变量值和协方差矩阵进行更新：

k+1|k+1=k+1|k+Kx(Zk+1-k+1)，

(36)

(37)

式中：Zk+1对应3.2节中提取的M(i)中的目标量测数据。同理可得Y在k+1时刻的状态估计值k+1|k+1及其对应的协方差矩阵Py,k+1|k+1.

由上述推导可以看出，BVCKF算法考虑了滤波过程中出现两个随机变量的情况，为同时处理目标与外辐射源的状态估计提供了可靠的理论工具。同时采用球面-径向容积准则来处理非线性函数的多维向量积分。这种理论方法上的改进可以从原理上降低APCL系统中外辐射源状态不确定性对目标状态估计造成的影响。相比于以往对单变量的SA方法[9,11]和新息协方差修正算法[11]而言，在算法精度和运算速度上都有着显著的提升。

3 仿真分析

为了充分验证所提算法的有效性，本文选取目标状态估计的均方根误差(RMSE)和算法的单次运行耗时作为性能评价指标，对比分析了DP-SA算法[11]和本文所提的DP-BVCKF算法在定位跟踪性能上的差异。由于文献[11]已经对DP-SA算法中DP阶段提取目标量测的性能进行了较为详细的仿真分析，本文不再赘述。

3.1 仿真参数设置

目标、外辐射源以及观测站的初始状态如表1所示。

表1 仿真场景的初始状态

目标和外辐射源共同的过程噪声为

Γk=diag{(5 m)2,(5 m)2,(0.5 m/s)2,(0.5 m/s)2}。

量测噪声的协方差表示为

Rk=diag{(0.005 rad)2,(0.005 rad)2,(50 m)2,(2 m/s)2}。

目标等效椭球体的长短轴半径B=1 m、C=0.5 m，外辐射源辐射功率密度E=1 W/m3. 杂波的密度参数λ=0.01 m-1·rad-1，幅值参数ζ=0.5，样本数c=6 000，系统的检测概率Pd=0.95,PFT=0.05. 探测时长M=50 s，探测间隔ΔT为1 s和在21～25 s时，目标和外辐射源做机动转弯运动，转弯率ω设定为0.1 rad/s. 其余时间目标和外辐射源的转弯率为0 rad/s，即做匀速直线运动；观测站在整个探测周期内作ω=0.1 rad/s的匀速圆周运动。蒙特卡洛仿真次数设定为200次。基于DP-BVCKF算法的仿真场景如图3所示。

图3 基于DP-BVCKF算法的仿真场景示意图

3.2 算法跟踪性能对比分析

为了更好地验证分析DP-SA算法和本文所提的DP-BVCKF算法在性能上的差异，本文分别基于不同的初始协方差和不同的转弯率给出了对比仿真的结果，如图4和图5所示。

图4 不同初始协方差下算法性能对比图

从图4和图5中可以看出：在目标做匀速运动时，两种算法都能取得很好的收敛效果，且无论处于何种情形下，DP-BVCKF算法的收敛精度都更高。这就很好地验证了本文的理论分析，因为在DP阶段两种算法的性能效果一致;而在状态估计阶段，BVCKF在原理上具有根本性的优势，既能避免SA增加的额外计算量，又能实时对外辐射源状态进行同步的更新估计，极大地降低了外辐射源状态不确定性的影响。但是DP-SA算法和本文所提的DP-BVCKF算法都存在一个明显的缺陷，即对目标机动时的跟踪效果不够理想，容易出现滤波发散的现象。而且，由于DP-SA算法中以扩展卡尔曼滤波作为核心滤波工具，相比于DP-BVCKF算法而言，滤波发散现象更加严重。

对比图4(a)、图4(b)、图4(c)、图4(d)可以看出，当初始协方差增大后，两种算法对目标和外辐射源的位置估计精度受影响较小。但是DP-BVCKF算法的状态估计效果要明显好于DP-SA算法。这是因为DP-SA算法进行SA后，目标状态和外辐射源状态之间并不独立，对协方差进行滤波处理时，二者状态的不确定性会互相影响。

由图5可以看出，当目标和外辐射源的转弯率增大后，两种算法对目标和外辐射源的位置估计精度都明显下降。但是DP-BVCKF算法受转弯率的影响相对较小，这是因为两种算法选用的核心滤波算法在原理上存在差异。仿真结果表明在处理非线性滤波跟踪问题时，基于球面-容积准则的确定性采样近似方法性能要明显优于1阶泰勒展开式的线性化处理方法。

表2给出了两种算法在不同杂波密度下的运行耗时。由表2可以看出：杂波密度的增大使得目标状态提取过程的耗时增加；由于DP-SA算法在SA时增加了一定的计算量，导致其在状态估计环节产生的耗时要高于DP-BVCKF算法。

4 结论

本文针对杂波环境下空基外辐射源系统的定位问题，提出了一种基于DP-BVCKF的定位跟踪算法。该方法首先利用DP算法从含有杂波的量测数据中提取出源于目标的量测信息，然后提出BVCKF算法来降低外辐射源状态不确定性对目标状态估计的影响。通过与DP-SA算法的仿真对比分析表明，本文所提算法能有效地提取出目标的量测信息，并极大地提高目标状态估计的精度，在兼顾实时性的同时保持较高的准确性。但是通常情况下杂波的统计特性是复杂多变的，本文仅考虑了杂波统计特性已知的情况下对目标数据的提取问题，当杂波统计特性未知时，该方法需要进一步改进。而且在非线性非高斯环境下，目标状态估计面临的问题将会进一步复杂化，尤其是在应对机动目标跟踪，以及多目标跟踪时，需要对算法做出调整，这将是后续的研究重点。