作者简介:陈 昊(1992-),男(汉族),安徽省滁州市人,西北工业大学硕士研究生. E-mail: 710276588@mail.nwpu.edu.cn
中文责编:方 圆; 英文责编:雨 辰
西北工业大学自动化学院, 西安 710129
Chen Hao, Yang Feng, Wang Yongqi, and Pan Quan School of Automation, Northwestern Polytechnical University, Xi'an 710129, P.R.China
information fusion; sky-wave over-the-horizon radar; multipath; Gaussian mixture probability hypothesis density(GMPHD); unscented Kalman filter; complex environment
DOI: 10.3724/SP.J.1249.2014.02130
针对天波超视距雷达(over-the-horizon radar, OTHR)多目标跟踪所面临的多路径和低检测概率问题,结合混合高斯概率假设密度(Gaussian mixture probability hypothesis density, GMPHD)滤波器免数据关联以及计算耗费低的优点,提出多路径不敏卡尔曼混合高斯概率假设密度(multipath unscented Kalman-Gaussian mixture probability hypothesis density, MPUK-GMPHD)融合框架.该框架将多路径效应等效为多传感器,构建面向OTHR的多传感器意义下的概率假设密度(probability hypothesis density, PHD)融合算法. 通过多路径信息的有效融合避免多路径独立量测更新带来PHD目标数过估问题,并采用不敏卡尔曼滤波处理量测模型非线性问题.仿真结果表明,在OTHR多目标跟踪的复杂环境下,MPUK-GMPHD融合算法能够较准确地估计目标状态和目标数,缓解了直接利用GMPHD滤波器处理带来的目标数过估和较大计算量的问题.
Sky wave over-the-horizon radar(OTHR)usually suffers from multipath propagation effects and low detection probabilities when tracking multiple targets. In this paper we present a multipath unscented Kalman-Gaussian mixture probability hypothesis density(MPUK-GMPHD)fusion framework to overcome these defects. The GMPHD filter is capable of avoiding data associations while maintaining low computational complexity. A probability hypothesis density fusion algorithm is proposed for OTHR, wherein multipath propagation effects are analyzed using a multi-sensor model. Particularly, the target number over-estimation problem is effectively avoided by fusing the multipath information, and an unscented Kalman filter is utilized to solve the nonlinear problem in measurement model. Simulation results show that the proposed MPUK-GMPHD framework estimates target state and number more accurately than conventional methods in complex environments, e.g. OTHR multi-target tracking. It overcomes the target number over-estimation and high computational complexity problems in existing GMPHD filter based algorithms.
高频雷达是远程或超远程目标探测系统中一种重要的探测手段. 其中,天波超视距雷达(over-the-horizon radar, OTHR)利用电磁波在电离层与地/海面之间的反射作用传输高频能量[1],其作用距离不受地球曲率的限制,可实现对隐身战斗机、轰炸机、无人机、洲际导弹、巡航导弹和大中型舰船等高价值目标的远程预警,具有典型的四抗性能(侦察-反侦察、干扰-反干扰、摧毁-反摧毁和隐身-反隐身).但相比常规雷达,OTHR数据处理面临低数据率、低检测概率和低测量精度以及多路径传播效应的问题[2]. 在低检测概率及密集杂波环境下进行目标跟踪往往会导致跟踪精度低、失跟率高的现象,同时多路径效应会导致同一目标产生多个量测,从而形成多条航迹[3-4]. 例如,多路径Viterbi数据关联跟踪算法(multipath Viterbi data association,MVDA)[5]及修正概率数据关联(modified probabilistic data association, MPDA)等[6-10]传统处理算法都涉及复杂的数据关联过程,计算量巨大且难以处理目标数时变问题.Mahler基于有限集统计(finite set statistics,FISST)理论,提出概率假设密度(probability hypothesis density,PHD),即目标状态后验密度的一阶矩,实现对目标状态和目标数的估计.PHD滤波器免数据关联的特性使其成为学界研究热点[11-12].对其积分是监视区域内的目标期望数,峰值点则对应目标状态[13]. PHD有两种实现方法:序贯蒙特卡洛PHD(sequential Monte Carlo PHD,SMC-PHD)[14-15]和混合高斯PHD(Gaussian mixture PHD,GM-PHD)[16-17].由于PHD滤波器在递归过程中存在函数的积分运算,所以不存在一般意义下的解析解,仅在线性高斯情况下,可利用高斯分布的性质得其解析表达形式.而混合高斯模型可以无限逼近任一分布,这也使GM-PHD成为研究热点.
OTHR多路径效应使目标通过多个传播路径产生多个对应量测,使用标准PHD滤波算法进行处理,得到的目标个数无法反映真实目标个数,且易形成目标同源航迹.此外,OTHR雷达量测方程存在严重的非线性,直接采用GM-PHD滤波器,其后验强度不能表示为混合高斯形式.传统算法多采用拓展卡尔曼滤波(extended Kalman filter,EKF)算法[18].EKF算法在非线性方程的线性化过程中仅使用了一次项,忽略了二次以上的高阶项,由此引入的误差降低了目标估计精度.
不敏卡尔曼滤波(unscented Kalman filter,UKF)算法首先对状态变量进行不敏变换得到新的扩维变换向量,然后由变换向量实现滤波过程.理论证明,对于任意非线性问题均可达到三阶以上精度.本研究将UKF算法与多路径融合算法相结合,给出一种面向OTHR多目标跟踪的多路径不敏卡尔曼GMPHD融合框架. 该框架将多路径效应等效为多传感器,构建面向OTHR的多传感器意义下的PHD滤波器,将每条路径的PHD作为一个子单元进行独立的预测和更新,构建多路径高斯分量矩阵,利用多传感器高斯分量合并算法对多路径信息进行有效融合[19].采用UKF算法[20]解决滤波过程中量测模型存在的非线性问题.针对PHD滤波器只能输出独立状态信息,并未给出目标航迹信息,采用二维指派方法给出每个目标的独立航迹.仿真结果表明,该滤波框架能有效处理复杂环境下多目标跟踪过程中的多路径问题,形成稳定系统航迹.
图1为OTHR多路径示意图. 由发射机传送的电磁信号经电离层反射到达目标,从目标传回的信号再次通过电离层反射回到接收机.假设存在2层分离的电离层,电离层虚高分别为hE和hF, 则至多存在4条传播路径:① EE路径.传送信号E层反射,接受信号E层反射; ② EF路径.传送信号E层反射,接受信号F层反射; ③ FE路径.传送信号F层反射,接受信号E层反射; ④ FF路径.传送信号F层反射,接受信号F层反射. 假设目标的动态模型为
(k+1)=F(k)x(k)+ν(k)(1)
zm(k)=Hm(x(k))+wm(k),
m=1, 2, 3, 4(2)
其中,m为传播路径; Hm(·)为量测函数. ν(k)和wm(k)是相互独立的零均值高斯白噪声,其协方差分别是Qk和Rmk; k时刻的目标状态空间x(k)=[ρ(k),(·overρ)(k), b(k), b·(k)]T包含目标在参考坐标系下距离接收机的地面径向距离ρ、 径向距离变化率(·overρ)、 方位角b及方位角变化率b·; 第m条路径的雷达量测 zm(k)=[Rmg(k), Rmr(k), Amz(k)]T包含斜距Rmg、 多普勒Rmr以及视方位角Amz, 多普勒Rgr是Rmg的变化率. 假设检测概率Pmd是不随时间和状态变化的常数定值,具体地理坐标系与雷达坐标系间的变换关系详见文献[21].
设密度融合算法
虽然PHD滤波器在进行状态估计中可免数据关联,但由于方程中存在积分与递归运算,故并无一般意义下的解析解,SMC执行方式虽然能处理一定数量的目标,但在杂波密度较高等复杂场景下,计算量会急剧上升.而利用不敏卡尔曼滤波处理量测模型非线性问题,同时假设检测概率和存活概率是与状态和时间无关的常量,可得非线性高斯条件下PHD迭代更新解析形式.针对多路径效应,本研究给出多路径不敏卡尔曼混合高斯概率假设密度融合算法的具体处理流程.
假设Dk-1(x)=∑Jk-1i=1ω(i)k-1N(x; m(i)k-1, P(i)k-1)为k-1时刻多目标PHD,由于量测方程存在一定的非线性,本研究用不敏卡尔曼滤波进行处理.首先计算2nx+1个δ采样点 ξj及对应权值Wj
{ξ(i)j(k-1|k-1)=m(i)k-1 j=0
ξ(i)j(k-1|k-1)=m(i)k-1+
(((nx+κ)P(i)k-1)1/2 )j j=1,…, nx
ξ(i)j+nx(k-1|k-1)=m(i)k-1-
(((nx+κ)P(i)k-1)1/2 )j j=1,…, nx(3)
{W(i)j=κ/(nk+κ)j=0
W(i)j=1/[2(nk+κ)] j=1,…, nx
W(i)j+nx=1/[2(nk+κ)] j=1,…, nx(4)
其中,κ是一个尺度参数,其取值为满足(nk+κ)≠0的任意值;(((nx+κ)P(i)k-1)1/2 )j是(nx+κ)P(i)k-1均方根矩阵的第j行或第j列; nx为状态向量的维数.根据状态方程(1),可得δ点的1步预测为
ξ(i)j(k|k-1)=Fk·ξ(i)j(k-1|k-1)(5)
其中,Fk为状态转移矩阵. 利用1步预测δ点 ξ(i)j(k|k-1)以及权值Wi, 可得状态预测估计和状态预测协方差
m(i)k|k-1=∑2nxj=0W(i)jξ(i)j(k|k-1)(6)
P(i)k|k-1=∑2nxj=0W(i)jΔX(i)j(k|k-1)
ΔX(i)j'(k|k-1)+Q(k)(7)
其中,ΔX(i)j(k|k-1)=ξ(i)j(k|k-1)-m(i)k|k-1.
根据量测方程(2),可得预测量测点
ζ(i)j(k|k-1)=Hm(ξ(i)j(k|k-1))(8)
预测量测和相应的协方差为
η(i)k|k-1=∑2nxj=0W(i)jζ(i)j(k|k-1)(9)
P(i)zz=R(k)+∑2nxj=0W(i)jΔZ(i)j(k|k-1)
ΔZ(i)j'(k|k-1)(10)
其中,ΔZ(i)j(k|k-1)=ζ(i)j(k|k-1)-η(i)k|k-1.
同样可得测量和状态向量的交互协方差
P(i)xz=∑2nxj=0W(i)jΔX(i)j(k|k-1)
ΔZ(i)j'(k|k-1)(11)
预测PHD仍可写成高斯混合形式
Dk|k-1(x)=∑Jk|k-1j=0ω(i)k|k-1N(x; m(i)k|k-1, P(i)k|k-1)(12)
利用量测信息对PHD进行更新
Dk(x)=(1-pD, k)Dk|k-1(x)+∑z∈ZkDD, k(x; Z)
DD,k(x; Z)=∑Jk|k-1j=1ω(j)k(z)N(x; m(j)k|k(z), P(j)k|k)(13)
对应估计状态、协方差及高斯项权值计算如下:
ω(j)k(z)=(pD, kω(i)k|k-1q(j)k(z))/(κk(z)+pD, k∑Jk|k-1=1ω()k|k-1q()k(z))(14)
m(j)k|k(z)=m(j)k|k-1+K(j)k(z-η(i)k|k-1)(15)
P(j)k|k=P(j)k|k-1-K(j)kS(j)kK(j)k'(16)
K(j)k=P(i)zz(17)
K(j)k=P(i)xz·(P(i)zz)-1(18)
综上,可将更新PHD表示为混合高斯形式
Dk(x)=∑Jki=1ω(j)k(z)N(x; m(i)k, P(i)k)(19)
将每条路径的PHD作为一个子单元按照上述两步进行独立的预测和更新.构建如下的PHD矩阵表达式
[D(11)k,1 D(12)k,1 … D(1j)k,1 … D(1M1k)k,1
D(21)k,1 D(22)k,1 … D(2j)k,1 … D(2M1k)k,1
…
D(i1)k,1 D(i2)k,1 … D(ij)k,1 … D(iM1k)k,1
…
D(Jk1)k,1 D(Jk2)k,1 … D(Jk j)k,1 … D(JkM1k)k,1 …|D(11)k,L D(12)k,L … D(1j)k,L … D(1MLk)k,L
D(21)k,L D(22)k,L … D(2j)k,L … D(2MLk)k,L
…
D(i1)k,L D(i2)k,L … D(ij)k,L … D(iMLk)k,L
…
D(Jk1)k,L D(Jk 2)k,L … D(Jk j)k,L … D(JkMLk)k,L][]+|D(10)k
D(20)k
D(i0)k
D(Jk0)k](20)
其中, “[]”表示在这个矩阵表达式中的高斯分量经过列行合并后所有的高斯分量相加; “”表示不考虑交集部分的集合并运算. 共有L条路径,其中,
D(ij)k, τ=ω(ij)k, τN(x; m(ij)k, τ, P(i)k, τ),
τ=1,2,…, L(21)
D(i0)k=ω(i0)kN(x; m(i)k|k-1, P(i)k|k-1)(22)
ω(i0)k=(1-pD, k)ω(i)k|k-1(23)
多路径PHD融合分为3个步骤: ① 根据预先设置的舍弃门限去除权值小于舍弃门限的PHD; ② 在PHD矩阵中对每一列进行PHD合并,即假设在同一条路径的量测集合中源自同一目标的真实量测最多只有一个,因此在进行列合并时,只有属于同一列的PHD才能合并; ③ 在PHD矩阵中进行PHD的行合并,得到最终的融合PHD,然后进行提取得到目标数和状态估计. 详细过程如下.
假设PHD的舍弃门限为Tth, 可得新的PHD权值矩阵W^-k, τ=[ω^-(ij)k, τ]Jk×Mτ,(τ=1,2,…,L). 其中,
ω^-(ij)k, τ={ω(ij)k, τ, ω(ij)k, τ>Tth
0, 其他(24)
这样每个PHD就变成D(ij)k, τ=ω^-(ij)k, τN(x; m(ij)k, τ, P(i)k, τ).
假设PHD列合并的合并门限为Uc, {D(ij)k, τ}Jki=1为PHD矩阵Dk=[D(ij)k, τ]Jk×Mτ的第j列. 对于每个τ=1,2,…,L, j=1,2,…,Mτ,
πj=argmax1≤i≤Jkω^-(ij)k, τ(25)
Ωj, τ={1≤i≤Jk|(m(πj j)k, τ-m(ij)k, τ)T
(P(ij)k, τ)-1(m(πj j)k, τ-m(ij)k, τ)≤Uc}(26)
ω~(πj j)k, τ=∑i∈Ωj, τω^-(ij)k, τ(27)
m~(πj j)k, τ=1/(ω~(πj j)k, τ)∑i∈Ωj, τω^-(ij)k, τm(ij)k, τ(28)
P~(πj j)k, τ=∑i∈Ωj, τ(ω^-(ij)k, τ)/(ω~(πj j)k, τ)[P(ij)k, τ+
(m(πj j)k, τ-m(ij)k, τ)(m(πj j)k, τ-m(ij)k, τ)T](29)
于是,进一步得到新的PHD权值矩阵
W^k, τ=[Aω^G1(ij)k, τ]Jk×Mτ, τ=1,2,…,L
其中,Aω^G1(ij)k, τ={ω~(πj j)k, τ, i=πj
0, i∈Ωj, τ, j≠0
ω(ij)k, τ, 其他(30)
此时,每个PHD为D(ij)k, τ=Aω^G1(ij)k, τN(x; Am^G1(ij)k, τ, P^(ij)k, τ). 其中,
Am^G1(ij)k, τ={m~(πj j)k, τ, i=πj
m(ij)k, τ, 其他
P^(ij)k, τ={P~(πj j)k, τ, i=πj
P(ij)k, τ, 其他(31)
记Aω^G1(j)k, τ=∑Jki=1Aω^G1(ij)k, τ,(τ=1, 2, …, L), 有0<Aω^G1(j)k, τ≤1.
假设PHD行合并的合并门限为Ur. 记Σi={D^)(l)k, i}Jk, il=1=Lτ=1{D(ij)k, i}Mτj=1, i=1,2,…,Jk,且D^)(l)k, i=w^)(l)k, iN(x; m^)(l)k, i, P^)(l)k, i), 因此有Jk, i=∑Lτ=1Mτ
Ii={l=1,…, Jk, i|ω^)(l)k, i>0}(32)
η=0
循环
η=η+1
ξi=argmaxl∈Ii ω^)(l)k, i(33)
Ξk, i={l∈Ii|(m^)(l)k, i-m^)(ξi)k, i)T(P^)(ξi)k, i+
P^)(l)k, i)-1(m^)(l)k, i-m^)(ξi)k, i)≤Ur}(34)
ω~i=∑l∈Ξk, iω^)(l)k, i(35)
ω~(η)k, i=1-∏l∈Ξk, i(1-ω^)(l)k, i)(36)
m~(η)k, i=1/(ω^-i)∑l∈Ξk, iω^)(l)k, im(l)k, i(37)
P~(η)k, i=∑l∈Ξk, i(ω^)(l)k, i)/(ω^-i)[P^)(l)k, i+
(m^)(l)k, i-m^)(ξi)k, i)(m^)(l)k, i-m^)(ξi)k, i)T](38)
Ii=Ii/Ξk, i(39)
直到Ii=.于是得Ψk, i={w^)(l)k, iN(x; m~(l)k, i, P~(l)k, i)},(i=1,2,…,Jk), 输出Ψk=∪Jki=1Ψk, i.
取权值大于等于0.5的高斯分量的均值作为目标状态的估计.
考虑电离层为E(高度hE=100 km)和F(高度hF=220 km)两层,存在EE、EF、FE和FF四种传播模式,雷达监视区域为1 000~1 400 km,径向距速率为0.013 889~0.222 2 km/s,方向角为0.409 81~0.074 53 rad,量测噪声标准差:径向距σRg=5 km,多普勒σRr=0.001 km/s,方向角σAz=0.003 rad; 扫描区域内杂波数服从泊松分布,杂波在整个监视区域内服从均匀分布,平均杂波个数为50个.采样时间Ts=10 s,存活概率ek|k-1=0.99.具体参数如下:
Fk=[1 Ts 0 0
0 1 0 0
0 0 1 Ts
0 0 0 1],
R(m)k=[25 0 0
0 1×10-6 0
0 0 9×10-6],
Qk=[25 0 0 0
0 1×10-6 0 0
0 0 9×10-6 0
0 0 0 1.8×10-11].
Tth=1×10-4, Uc=Ur=4.
为验证算法有效性,采用最优次模式分配(optimal subpattern assignment,OSPA)[23]指标作为性能评判准则.
对于任意k∈N={1,2,…},Ψk为{1,2,…, k}的全排列组合, p∈[1,∞],c>0, 给定2个多目标状态集合X={x1, x2, …, xm}和Y={y1, y2, …, yn},m, n∈N0={0,1,2,…}, 如果m≤n, 两集合间的OSPA指标为
dcp(X,Y)=(1/nminΨ∈Ψn∑mi=1d(c)(xi, yΨ(i))p+
cp(n-m))1/p(40)
否则dcp(X,Y):=dcp(Y, X)(41)
其中, d(c)(x, y):=min(c, dp(x, y))为x, y∈Rn之间的距离; dp为p范数; c为截断距离. 仿真过程取p=2, c=10 km.
考虑单目标在监测区域内做匀速直线运动,将本研究提出的MPUK-GMPHD融合算法与OTHR跟踪处理的经典算法MPDA算法[24]进行对比分析,目标的初始状态为x=(1 100 km, 0.15 km/s, 0.157 2 rad,0.000 087 2 rad/s)
考虑多目标情况,验证所提出MPUK-GMPHD算法的有效性. 3个目标在监视区域内做匀速直线运动,x1=(1 100 km, 0.15 km/s, 0.465 28 rad, 0.000 087 2 rad/s),x2=(1 150 km, 0.15 km/s, 1.043 43 rad, -0.000 087 2 rad/s)以及x3=(1 200 km, 0.21 km/s, 0.8 rad, -0.000 067 2 rad/s)分别为3个目标的初始状态,忽略衍生目标出现情况,整个仿真时间持续600 s.
图2为单目标情况下MPUK-GMPHD算法和MPDA算法的估计航迹图. 可见,两种算法都能对目标进行有效跟踪.图3为50次仿真下两种算法的平均OSPA距离图,
由仿真结果可知,本研究提出的MPUK-GMPHD在精度上优于MPDA算法. MPDA需要考虑数据关联,这存在极大的风险性,而MPUK-GMPHD具有免数据关联的优点.MPDA在进行状态估计时,先将多路径量测信息进行合并,这导致部分信息的缺失,而MPUK-GMPHD在进行更新时,使用了所有的量测信息,然后进行多路径信息的融合,对雷达传感器采集数据的使用率更高,跟踪效果更好.图4为目标的多路径量测图. 可见,由于OTHR多路径效应,一个目标产生了多个同源量测,且不同路径的量测之间可能存在交叉和重叠,增加了环境的复杂程度,用传统方法进行处理易造成目标过估,形成多条同源航迹.图5为地理坐标系中3目标航迹图,仿真结果显示,本研究提出的MPUK-GMPHD融合算法可以实现对多目标的有效跟踪,通过融合算法将多路径信息融合之后形成稳定航迹.图6为MPUK-GMPDH多目标跟踪的OSPA距离图. 可见,MPUK-GMPHD融合算法具有较高的估计精度.图7为目标数估计图,结果显示MPUK-GMPHD融合算法有效避免了因多路径效应产生的目标数过估问题.
本研究提出一种面向天波超视距雷达多目标跟踪的MPUK-GMPHD融合算法,能有效避免复杂的数据关联,从而降低了计算量,增加了算法的执行效率.将多路径问题等效为多传感器,独立预测更新每条路径的PHD值并借助多传感器高斯分量合并算法对多路径信息进行有效融合,避免了目标数过估等问题.仿真结果表明,在复杂环境下, MPUK-GMPHD算法性能优异,仿真精度高,能够形成稳定航航迹.
深圳大学学报理工版
JOURNAL OF SHENZHEN UNIVERSITY SCIENCE AND ENGINEERING
(1984年创刊 双月刊)
主 管 深圳大学
主 办 深圳大学
编辑出版 深圳大学学报理工版编辑部
主 编 阮双琛
国内发行 深圳市邮电局
国外发行 中国国际图书贸易集团有限公司(北京399信箱)
地 址 北京东黄城根北街16号
邮 编 100717
电 话 0755-26732266
0755-26538306
Email journal@szu.edu.cn
标准刊号 ISSN 1000-2618
CN 44-1401/N