-
无源定位方法在航空和电子领域广泛应用，与主动雷达相比具有体积小、隐蔽性高及低空探测能力等优点[1]，近年来成为各领域的研究热点。可用于定位的外辐射源包括调频广播FM[2] 信号、数字地面广播(digital video broadcasting-terrestrial，DVB-T)[3]信号、手机信号基站WiFi[4]等常见的民用信号。与多站系统相比，单站系统机动性强，不存在时间和数据同步的问题。因此研究高精度的单站外辐射源定位方法具有重要意义。
目标跟踪是典型的估计问题，关键是通过信号处理手段获取的观测信息估计得到目标的状态信息。根据利用的观测信息的不同，算法包括仅测方向(direction of arrival，DOA)[5]、仅测时差(time difference of arrival，TDOA)[6]、联合角度和时差[7]、联合角度和频差 (frequency difference of arrival，FDOA)[8]以及联合角度时差频差的定位跟踪算法[9]等。不同的观测信息使系统的复杂度不同，信息量越少越单一，系统越简单；信息量越多，系统越复杂。本文研究单站接收多外辐射源信号获取角度和时差信息构建定位跟踪模型，然后利用跟踪算法实现目标的高精度跟踪。
求解这类问题的常用算法包括扩展卡尔曼滤波(extended Kalman filter，EKF)算法、无迹卡尔曼滤波(unscented Kalman filter，UKF)算法和粒子滤波(particle filter，PF)算法[9]。EKF算法及其改进的算法因为其原理简单、计算量低而被长期研究和应用[10]。EKF算法将非线性估计问题通过泰勒展开线性化之后用传统的卡尔曼滤波算法实现[9]。但这类算法需要计算Hessian矩阵，对于比较复杂的非线性模型无疑增加了系统的复杂度。通过构建伪线性模型，利用卡尔曼滤波算法直接估计目标状态，减少了算法的复杂度，这种算法称为伪线性卡尔曼滤波PKF[4]。该算法最初被应用于仅测角的无源定位中，但需角度具有很高的精度。Koteswara将所有的观测量构造成一个向量，进而采用一种新的伪线性卡尔曼滤波(pseudo-linear Kalman filter，PKF)算法[11, 12]实现目标跟踪。但这种算法只能应用于仅有一个观测量的仅侧角模型中，且目标运动形式简单，算法使用限制高。
本文研究联合DOA和TDOA的单站外辐射源跟踪问题。首先根据几何关系构建观测方程，进一步构建基于角度和时差的定位跟踪模型，变形为线性形式。然后利用传统的卡尔曼滤波算法实现目标跟踪，滤波的初值由最小二乘(least squares，LS)算法获取，从而降低了整个系统的复杂度。进一步分析伪线性跟踪模型，通过迭代PKF(iterative PKF，IPKF)算法实现目标跟踪。
-
假设观测站位于坐标原点XR=[0,0,0]T，目标位置X=[x,y,z]T，辐射源位置XTrk=[xk,yk,zk]T，辐射源到观测站的距离${{d}^{k}}=\sqrt{{{\left( {{x}^{k}} \right)}^{2}}+{{\left( {{y}^{k}} \right)}^{2}}+{{\left( {{z}^{k}} \right)}^{2}}}$，目标到观测站的距离为$R=\sqrt{{{x}^{2}}+{{y}^{2}}+{{z}^{2}}}$，辐射源到目标的距离${{r}^{k}}=\sqrt{\left( {{\left( x-{{x}^{k}} \right)}^{2}}+{{\left( y-{{y}^{k}} \right)}^{2}}+{{\left( z-{{z}^{k}} \right)}^{2}} \right)}$，信号直达波与经目标反射的回波时间差为τk，其中k=1,…,K，K为辐射源个数。该系统如图 1所示。定义目标的平面距离${{R}_{2}}=\sqrt{{{x}^{2}}+{{y}^{2}}}$。
根据目标的角度，建立方位角的观测方程：
$$\theta =arctan\left( y/x \right)\Leftrightarrow xsin\left( \theta \right)=ycos\left( \theta \right)$$ (1) 由于实际测量方位角存在测量误差，所以有：
$${{\theta }_{m}}=\theta +{{v}_{\theta }}$$ (2) 将式(2)带入xsin(θm)－ycos(θm)，并将其在θ处进行泰勒展开，并省略误差高次项，可得：
$$\begin{align} & xsin\left( {{\theta }_{m}} \right)-ycos\left( {{\theta }_{m}} \right) \\ & =xsin\left( \theta \right)+xcos\left( \theta \right){{v}_{\theta }}-ycos\left( \theta \right)+ysin\left( \theta \right){{v}_{\theta }} \\ & =xsin\left( \theta \right)-ycos\left( \theta \right)+xcos\left( \theta \right)+ysin\left( \theta \right){{v}_{\theta }} \\ \end{align}$$ (3) 根据图 1的几何关系可知：
$${{R}_{2}}=\sqrt{{{x}^{2}}+{{y}^{2}}}=xcos\theta +ysin\theta $$ (4) 
图 1 单站PCL定位跟踪系统
Figure 1. Single-Observer PCL Positioning and Tracking System
将式(2)和式(3)带入式(4)并将之化简可得：
$$xsin({{\theta }_{m}})-ycos({{\theta }_{m}})={{R}_{2}}{{v}_{\theta }}$$ (5) 构建关于俯仰角的观测方程：
$$\varphi =arctan\left( z/\sqrt{{{x}^{2}}+{{y}^{2}}} \right)~$$ (6) 将式(4)代入式(6)得到R2sinφ=zcosφ。利用图 1的几何关系可知R=R2cosφ+zsinφ。对xcosθm+ysinθm在θ处泰勒展开可得：
$$\begin{align} & xcos\left( {{\theta }_{m}} \right)+ysin\left( {{\theta }_{m}} \right) \\ & =xcos\left( \theta \right)-xsin\left( \theta \right){{v}_{\theta }}+ysin\left( \theta \right)+ycos\left( \theta \right){{v}_{\theta }} \\ & ={{R}_{2}}+\left[ ycos\left( \theta \right)-xsin\left( \theta \right) \right]{{v}_{\theta }}={{R}_{2}}~ \\ \end{align}$$ (7) 式(6)可写为：
$$\begin{align} & xcos{{\theta }_{m}}sin{{\varphi }_{m}}+ysin{{\theta }_{m}}sin{{\varphi }_{m}}-zcos{{\varphi }_{m}} \\ & ={{R}_{2}}sin{{\varphi }_{m}}-zcos{{\varphi }_{m}} \\ \end{align}$$ (8) 对式(8)泰勒展开，可得：
$$\begin{align} & {{R}_{2}}sin({{\varphi }_{m}})-zcos({{\varphi }_{m}}) \\ & ={{R}_{2}}sin\left( \varphi \right)+{{R}_{2}}cos\left( \varphi \right){{v}_{\varphi }}-zcos\left( \varphi \right)+zsin\left( \varphi \right){{v}_{\varphi }} \\ & ={{R}_{2}}sin\left( \varphi \right)-zcos\left( \varphi \right)+\left[ {{R}_{2}}cos\left( \varphi \right)+zsin\left( \varphi \right) \right]{{v}_{\varphi }} \\ & =R{{v}_{\varphi }} \\ \end{align}$$ (9) 所以俯仰角的观测方程为：
$$R{{v}_{\varphi }}=xcos({{\theta }_{m}})sin({{\varphi }_{m}})+ysin({{\theta }_{m}})sin({{\varphi }_{m}})-zcos({{\varphi }_{m}})$$ (10) 根据时差的关系，构建方程：
$${{\tau }^{k}}=(R+{{r}^{k}}-{{d}^{k}})/c$$ (11) 根据R和rk的定义式，计算：
$$\begin{align} & {{R}^{2}}-{{({{r}^{k}})}^{2}} \\ & =(R+{{r}^{k}})(R-{{r}^{k}}) \\ & =2{{x}^{k}}x+2{{y}^{k}}y+2{{z}^{k}}z-{{({{x}^{k}})}^{2}}-{{({{y}^{k}})}^{2}}-{{({{z}^{k}})}^{2}} \\ \end{align}$$ (12) 将式(11)带入式(12)，并加上式(11)，可得：
$$R=\frac{{{x}^{k}}x+{{y}^{k}}y+{{z}^{k}}z}{c{{\tau }^{k}}+{{d}^{k}}}-\frac{{{({{d}^{k}})}^{2}}}{2(c{{\tau }^{k}}+{{d}^{k}})}+\frac{1}{2}(c{{\tau }^{k}}+{{d}^{k}})$$ (13) 由于cτk+dk>0，两边同时乘cτk+dk，得：
$$R(c{{\tau }^{k}}+{{d}^{k}})={{x}^{k}}x+{{y}^{k}}y+{{z}^{k}}z-12{{({{d}^{k}})}^{2}}+\frac{1}{2}{{(c{{\tau }^{k}}+{{d}^{k}})}^{2}}$$ (14) 根据时差观测值和真实值的关系，建立方程：
$${{\tau }^{k}}_{m}={{\tau }^{k}}+{{v}_{\tau }}$$ (15) 由于需对式(14)进行线性化建模，需将式(14)中的R用x、y、z和θ、φ来表示，然而实际观测方程中，只能得到θm、φm。
下面证明在省略误差高次项的条件下：
$$R=xcos\left( {{\theta }_{m}} \right)cos\left( {{\varphi }_{m}} \right)+ysin\left( {{\theta }_{m}} \right)cos\left( {{\varphi }_{m}} \right)+zsin({{\varphi }_{m}})$$ (16) 式(7)已证明在省略误差高次项条件下xcos(θm)+ysin(θm)=R2，所以式(16)的右边等于R2cos(φm)+zsin(φm)，即需证：
$$R={{R}_{2}}cos({{\varphi }_{m}})+zsin({{\varphi }_{m}})$$ (17) 式中，右边在φ处泰勒展开，省略高次误差项得：
$$\begin{align} & {{R}_{2}}cos\left( \varphi \right)-{{R}_{2}}sin\left( \varphi \right){{v}_{\varphi }}+zsin\left( \varphi \right)+zcos\left( \varphi \right){{v}_{\varphi }}= \\ & {{R}_{2}}cos\left( \varphi \right)+zsin\left( \varphi \right)+\left[ zcos\left( \varphi \right)-{{R}_{2}}sin\left( \varphi \right) \right]{{v}_{\varphi }} \\ \end{align}$$ (18) 根据已知几何关系R=R2cosφ+zsinφ;zcos(φ)=R2sin(φ)，所以式(18)化为：
$${{R}_{2}}cos({{\varphi }_{m}})+zsin({{\varphi }_{m}})=R$$ (19) 命题得证。
将式(15)带入式(14)，并将其在τk处泰勒展开，得到：
$$\begin{align} & {{x}^{k}}x+{{y}^{k}}y+{{z}^{k}}z-R\left( c{{\tau }^{k}}_{m}+{{d}^{k}} \right)-\frac{1}{2}{{\left( {{d}^{k}} \right)}^{2}}+ \\ & \frac{1}{2}{{({{c}^{k}}_{m}+{{d}^{k}})}^{2}}=({{c}^{2}}{{\tau }^{k}}_{m}+c{{d}^{k}}-cR){{v}^{k}}_{\tau }~ \\ \end{align}$$ (20) 令ρk=cτk+dk，则ρmk=cτmk+dk；将式(16)带入式(20)，构建关于时差的观测方程：
$$\begin{align} & [{{z}^{k}}-{{\rho }^{k}}_{m}sin{{\varphi }_{m}}]z-\frac{1}{2}{{({{d}^{k}})}^{2}}+\frac{1}{2}{{({{\rho }^{k}}_{m})}^{2}} \\ & =(c{{\rho }^{k}}_{m}-c{{r}_{1}}){{v}_{\tau }}~ \\ \end{align}$$ (21) 得到观测方程组为：
$$\left\{ \begin{array}{*{35}{l}} ~0=xsin({{\theta }_{m}})-ycos({{\theta }_{m}})-{{r}_{2}}{{v}_{\theta }} \\ 0=xcos{{\theta }_{m}}sin{{\varphi }_{m}}+ysin{{\theta }_{m}}sin{{\varphi }_{m}}-zcos{{\varphi }_{m}}-{{r}_{1}}{{v}_{\varphi }} \\ {{({{d}^{k}})}^{2}}-{{(\rho _{k}^{k})}^{2}}=2[{{x}^{k}}-{{\rho }^{k}}_{m}cos{{\theta }_{m}}cos{{\varphi }_{m}}]x+ \\ 2[{{y}^{k}}-\rho _{m}^{k}sin{{\theta }_{m}}cos{{\varphi }_{m}}]y+ \\ 2[{{z}^{k}}-\rho _{m}^{k}sin{{\varphi }_{m}}]z-2(c\rho _{m}^{k}-c{{r}_{1}}){{v}_{\tau }}, \\ k=1,2,\ldots ,K \\ \end{array} \right.$$ (22) 将式(22)构建成矩阵形式，得到n时刻的观测方程：
$${{Z}^{o}}_{m}\left( n \right)=H_{m}^{o}\left( n \right){{X}^{o}}\left( n \right)+v_{m}^{o}\left( n \right)$$ (23) 式中的各矩阵变量为：
$$\begin{align} & Z_{m}^{o}\left( n \right)= \\ & {{[0,0,{{({{d}^{1}})}^{2}}-{{(\rho _{m}^{1}\left( n \right))}^{2}},\ldots ,{{({{d}^{K}})}^{2}}-{{({{\rho }^{K}}_{m}\left( n \right))}^{2}}]}^{T}} \\ & {{v}^{o}}_{m}\left( n \right)= \\ & [-{{R}_{2}}\left( n \right){{v}_{\theta }},-R\left( n \right){{v}_{\varphi }},-2(c{{\rho }^{1}}_{m}\left( n \right)-cR\left( n \right))\centerdot \\ & v_{\tau }^{1},\ldots ,-2({{c}^{2}}\rho _{m}^{k}\left( n \right)-cR\left( n \right))v_{\tau }^{k}{{]}^{T}} \\ \end{align}$$ 那么，观测误差的协方差矩阵为：
$$\begin{align} & {{Q}^{o}}_{m}\left( n \right)=diag[{{R}^{2}}_{2}\left( n \right){{\sigma }^{2}}_{\theta },{{R}^{2}}\left( n \right){{\sigma }^{2}}_{\varphi },4(c{{\rho }^{1}}_{m}\left( n \right)- \\ & cR{{)}^{2}}{{\sigma }^{2}}_{{{\tau }^{1}}},\ldots ,4{{({{c}^{2}}{{\rho }^{k}}_{m}\left( n \right)-cR\left( n \right))}^{2}}{{\sigma }^{2}}_{{{\tau }^{K}}}{{]}^{T}} \\ \end{align}$$ 观测矩阵为：
$${{H}^{o}}_{m}\left( n \right)=\left[ \begin{array}{*{35}{l}} sin({{\theta }_{m}}) & -cos({{\theta }_{m}}) & 0 \\ cos{{\theta }_{m}}sin{{\varphi }_{m}} & sin{{\theta }_{m}}sin{{\varphi }_{m}} & -cos{{\varphi }_{m}} \\ 2({{x}^{1}}-{{\rho }^{1}}_{m}\left( n \right)cos{{\theta }_{m}}cos{{\varphi }_{m}}) & 2({{y}^{1}}-{{\rho }^{1}}_{m}\left( n \right)sin{{\theta }_{m}}cos{{\varphi }_{m}}) & 2({{z}^{1}}-{{\rho }^{1}}_{m}\left( n \right)sin{{\varphi }_{m}}) \\ \vdots & \vdots & \vdots \\ ~2({{x}^{K}}-{{\rho }^{K}}_{m}\left( n \right)cos{{\theta }_{m}}cos{{\varphi }_{m}}) & 2({{y}^{K}}-{{\rho }^{K}}_{m}\left( n \right)sin{{\theta }_{m}}cos{{\varphi }_{m}}) & 2({{z}^{K}}-{{\rho }^{K}}_{m}\left( n \right)sin{{\varphi }_{m}}) \\ \end{array} \right]$$ -
Kalman滤波算法是基于最小均方误差算法导出的，其思想是构建未知参数的线性矢量模型，通过递推处理寻求状态矢量的最小误差意义下的解。该算法要求相邻时刻的状态之间和状态与观测量之间满足线性关系。本文所研究的观测模型中，观测向量和观测矩阵均与含有误差的观测量有关，所以选择基于PKF展开。
根据目标运动形式不同来构建状态转移方程，主要包括匀速运动、匀加速运动和机动形式。本文主要研究滤波算法，并假设目标为匀速运动，且其状态矢量为X(n)=[x(n),y(n),z(n),(n),(n),(n)]T，目标运动状态转移方程为：
$$X\left( n+1 \right)=\Phi \left( n+1,n \right)X\left( n \right)+w\left( n \right)$$ (24) 其中，$\Phi \left( n+1,n \right)=\left[ \begin{array}{*{35}{l}} {{I}_{3}} & T{{I}_{3}} \\ 0 & {{I}_{3}} \\ \end{array} \right]$为状态转移矩阵；I3为3×3的单位矩阵；T为观测时间间隔。假设w(n)均值为零，方差为Q(n)过程噪声。
结合Kalman滤波的基本原理，本文PKF目标跟踪算法的过程如下。
1) 滤波初始化。根据伪线性模型，直接采用最小二乘算法获取目标的初始位置：
$$\begin{align} & {{{\hat{X}}}^{o}}\left( 1 \right)={{\left[ H_{m}^{o}\left( 1 \right){{\left[ H_{m}^{o}\left( 1 \right) \right]}^{T}} \right]}^{-1}}\centerdot \\ & {{\left[ H_{m}^{o}\left( 1 \right) \right]}^{T}}Z_{m}^{o}\left( 1 \right) \\ \end{align}$$ (25) 同理得下一时刻的目标位置估计${{{\hat{X}}}^{o}}$(2)，进而得到目标的初始速度为：
$${{{\hat{V}}}^{o}}\left( 1 \right)={{{\hat{X}}}^{o}}\left( 2 \right)-{{{\hat{X}}}^{o}}\left( 1 \right)$$ (26) 所以滤波的初值等于：
$$\hat{X}\left( 1 \right)={{[{{({{{\hat{X}}}^{o}}\left( 1 \right))}^{T}},{{({{{\hat{V}}}^{o}}\left( 1 \right))}^{T}}]}^{T}}$$ (27) 目标的状态误差和观测误差的协方差矩阵Qmo(1)、P(1)均为已知。预测误差相关矩阵$K\left( 1,0 \right)=E\left\{ \left[ \hat{X}\left( 1 \right)-\bar{X}\left( 1 \right) \right]{{\left[ \hat{X}\left( 1 \right)-\bar{X}\left( 1 \right) \right]}^{\text{T}}} \right\}$。将式(1)、式(6)和式(11)联合构造非线性观测方程：
$$Z=f({{X}^{o}})~$$ (28) 式中,Z=[θ,φ,τ1,…,τK]T，这里针对初始时刻分析，省略Z(1)中的(1)，对式(28)两边微分：
$$dZ={{\frac{\partial f({{X}^{o}})}{\partial {{X}^{o}}}}^{o}}d{{X}^{o}}~$$ (29) 令A${{A}^{o}}=\partial f\left( {{X}^{o}} \right)/\partial {{X}^{o}}$，那么$\text{d}{{X}^{o}}\text{=}{{\left( {{A}^{o}}{{\left( {{A}^{o}} \right)}^{\text{T}}} \right)}^{-1}}{{\left( {{A}^{o}} \right)}^{\text{T}}}\text{d}Z$，所以有：
$$\begin{align} & E\left[ d{{X}^{o}}{{\left( d{{X}^{o}} \right)}^{T}} \right]={{\left[ {{A}^{o}}{{\left( {{A}^{o}} \right)}^{T}} \right]}^{-1}}{{\left( {{A}^{o}} \right)}^{T}}{{Q}^{o}}{{A}^{o}} \\ & {{[{{A}^{o}}{{({{A}^{o}})}^{T}}]}^{-1}} \\ \end{align}$$ (30) 其中，${{Q}^{o}}=E\left[ \text{dZ}\cdot \text{d}{{\text{Z}}^{\text{T}}} \right]\text{=}\left[ \sigma _{\theta }^{2},\sigma _{\varphi }^{2},\sigma _{\tau }^{2},\cdots ,\sigma _{\tau }^{2} \right]$。将式(25)估计得到的${{{\hat{X}}}^{o}}$(1)计算Ao(1)，将Ao(1)带入式(30)，得到：
$$K\left( 1,0 \right)=\left[ \begin{array}{*{35}{l}} Ed{{X}^{o}}{{(d{{X}^{o}})}^{T}}] & 0 \\ 0 & U \\ \end{array} \right]$$ (31) 其中,U是初始时刻的速度和加速度误差的协方差矩阵，因难以得到理论值，给定一个初值，U=I6。
2) 状态预测。滤波算法开始n=1,2,3,…。
$$\hat{X}\left( n+1,n \right)=\Phi \hat{X}\left( n \right)$$ (32) 由于目标的状态不仅包括其位置，还包括速度，因此目标的观测矩阵应该为：
$$H\left( n \right)=\left[ H_{m}^{o}\left( n \right),{{0}^{\left( K+2 \right)\times 3}} \right]$$ (33) 3) 预测的协方差矩阵：
$$K\left( n,n-1 \right)=\Phi K\left( n-1,n-1 \right){{\Phi }^{H}}~$$ (34) 4) 计算增益矩阵：
$${{\left[ H\left( n \right)K\left( n,n-1 \right){{H}^{T}}\left( n \right)+Q_{m}^{o}\left( n \right) \right]}^{-1}}$$ (35) 5) 状态更新：
$$\begin{align} & \hat{X}~\left( n+1 \right)=\hat{X}\left( n+1,n \right)+ \\ & G\left( n \right)\left[ Z_{m}^{o}\left( n \right)-H\left( n \right)\hat{X}(n) \right] \\ \end{align}$$ (36) 6) 估计误差协方差矩阵更新：
$$\begin{align} & K\left( n+1,n \right)=\left[ \Phi -G\left( n \right)H\left( n \right) \right]K\left( n,n-1 \right) \\ & {{\left[ \Phi -G\left( n \right)H\left( n \right) \right]}^{T}}+P\left( n \right)+G\left( n \right)Q\left( n \right){{G}^{T}}\left( n \right)~ \\ \end{align}$$ (37) 7) 返回步骤2)。
-
本文§2推导了PKF算法的观测方法和观测误差的协方差矩阵，但是观测误差的协方差矩阵与目标到观测站的距离有关：
$$\begin{align} & {{Q}^{o}}_{m}\left( n \right)=diag[{{R}^{2}}_{2}\left( n \right){{\sigma }^{2}}_{\theta },{{R}^{2}}\left( n \right){{\sigma }^{2}}_{\varphi },4(c{{\rho }^{1}}_{m}\left( n \right)- \\ & cR{{)}^{2}}{{\sigma }^{2}}_{{{\tau }^{1}}},\ldots ,4{{({{c}^{2}}{{\rho }^{k}}_{m}\left( n \right)-cR\left( n \right))}^{2}}{{\sigma }^{2}}_{{{\tau }^{K}}}{{]}^{T}} \\ \end{align}$$ 由于目标当前时刻位置是未知的，实际处理过程中利用前一时刻的位置代替来计算误差的协方差矩阵，使滤波算法得以完成，但替代使协方差矩阵不够准确，进而使滤波的结果不够准确。针对此问题，本文提出通过迭代方法提高PKF算法中的观测误差协方差矩阵的准确性，进而提高跟踪误差，即IPKF目标跟踪算法。算法步骤如下。
1) 初始化；
2) 状态预测，计算预测的协方差矩阵，与PKF算法相同；
3) 迭代滤波，1≤i≤ L，L为迭代次数；
(1) 计算增益：
$$\begin{align} & {{G}^{i}}\left( n \right)=\Phi K\left( n,n-1 \right)\centerdot \\ & \Phi {{\left[ H\left( n \right)K\left( n,n-1 \right){{H}^{T}}\left( n \right)+{{[Q_{m}^{o}\left( n \right)]}^{i}} \right]}^{-1}}~ \\ \end{align}$$ (38) (2) 状态滤波：
$$\begin{align} & {{{\hat{X}}}^{i}}\left( n \right)=\hat{X}\left( n,n-1 \right)+ \\ & {{G}^{i}}\left( n \right)\left[ Z_{m}^{o}\left( n \right)-H\left( n \right)\hat{X}\left( n,n-1 \right) \right] \\ \end{align}$$ (39) (3) 更新协方差矩阵：
$$\begin{align} & {{K}^{i}}\left( n,n \right)=\left[ \Phi -{{G}^{i}}\left( n \right)H_{m}^{o}\left( n \right) \right]K\left( n,n-1 \right)\left[ \Phi - \right. \\ & G\left( n \right)H_{m}^{o}{{\left. \left( n \right) \right]}^{T}}+P\left( n \right)+G\left( n \right){{\left[ \left( n \right) \right]}^{i}}{{G}^{T}}\left( n \right) \\ \end{align}$$ (40) (4) i←i+1，返回第(1)步；
4) n←n+1，返回步骤(2)。
-
对于实时定位的模型，CRLB是目标定位误差的理想下界。CRLB由Fisher信息矩阵确定，假设矩阵为J，则矩阵的元素为：
$$\begin{align} & {{\left[ J\left( X \right) \right]}_{kl}}=E\left[ \frac{\partial lnp\left( Z,X \right)}{\partial k}\frac{\partial lnp\left( Z,X \right)}{\partial l} \right], \\ & \left( k,l=x,y,z \right) \\ \end{align}$$ (41) 式中,k、l代表 3×3矩阵J的下标；x、y、z为J的下标取值；该表达式的含义为，矩阵J对向量X=[z,y,z]T求导，即对每个元素求偏导，得到J的各个元素，即$J=\left[ \begin{array}{*{35}{l}} {{J}_{xx}} & {{J}_{xy}} & {{J}_{xz}} \\ {{J}_{yx}} & {{J}_{yy}} & {{J}_{yz}} \\ {{J}_{zx}} & {{J}_{zy}} & {{J}_{zz}} \\ \end{array} \right]$。
但是对于目标跟踪系统，定位的结果不仅与当前时刻的观测量有关，还与前一时刻的观测量有关，所以式(41)不能作为跟踪模型的CRLB。文献[16]推导了非线性滤波算法的CRLB。假设当前时刻的Fisher信息矩阵为Jn：
$${{J}_{n+1}}=D_{n}^{22}-D_{n}^{12}({{J}_{n}}+D_{n}^{11})D_{n}^{12}$$ (42) 式中，
$$D_{n}^{11}=E\left\{ -\Delta _{{{x}_{n}}}^{{{x}_{n}}}ln[p({{x}_{n+1}}|{{x}_{n}})] \right\}~$$ (43) $$D_{n}^{12}=E\left\{ -\Delta _{{{x}_{n}}}^{{{x}_{n+1}}}ln\left[ p({{x}_{n+1}}|{{x}_{n}}) \right] \right\}={{\left[ D_{n}^{21} \right]}^{T}}$$ (44) $$\begin{align} & D_{n}^{22}=E\left\{ -\Delta _{{{x}_{n+1}}}^{{{x}_{n+1}}}ln[p({{x}_{n+1}}|{{x}_{n}})] \right\}+ \\ & E\left\{ -\Delta _{{{x}_{n+1}}}^{{{x}_{n+1}}}ln[p({{z}_{n+1}}|{{x}_{n+1}})] \right\} \\ \end{align}$$ (45) 式中,Δab{·}为·对a和b求偏导。
根据式(24)可得：
$$\begin{align} & p({{X}_{n+1}}|{{X}_{n}})=\frac{1}{2\pi Q{{\left( n \right)}^{1/2}}}\cdot \text{ } \\ & ~exp\left\{ {{[{{X}_{n+1}}-\Phi {{X}_{n}}]}^{T}}{{Q}^{-1}}[{{X}_{n+1}}-\Phi {{X}_{n}}]~ \right\} \\ \end{align}$$ (46) 根据式(28)，可得：
$$\begin{align} & p\left( {{Z}_{n+1}}|{{X}_{n+1}} \right)=\frac{1}{{{\left| 2\pi R\left( n+1 \right) \right|}^{{1}/{2}\;}}}\cdot \text{ }\cdot exp \\ & \left\{ {{[{{Z}_{n+1}}-f({{X}_{n+1}})]}^{T}}{{R}^{-1}}\left( n+1 \right)[{{Z}_{n+1}}-f({{X}_{n+1}})] \right\} \\ \end{align}$$ (47) 将式(46)和式(47)带入式(43)~(45)计算得到：
$$\begin{align} & \left\{ {{[{{Z}_{n+1}}-f({{X}_{n+1}})]}^{T}}{{R}^{-1}}\left( n+1 \right)[{{Z}_{n+1}}-f({{X}_{n+1}})] \right\} \\ & D_{n}^{11}=-{{\Phi }^{T}}{{Q}^{-1}}\left( n \right)\Phi \\ \end{align}$$ (48) $$D_{n}^{12}=D_{n}^{21}={{\Phi }^{T}}{{Q}^{11}}\left( n \right)$$ (49) $$\begin{align} & D_{n}^{22}={{Q}^{-1}}\left( n+1 \right)+{{M}^{T}}\left( n+1 \right){{R}^{-1}}\left( n+1 \right) \\ & M\left( n+1 \right) \\ \end{align}$$ (50) 其中，
$$M\left( n+1 \right)=\frac{\partial f({{X}_{n+1}})}{\partial {{X}_{n+1}}}$$ (51) 将式(48)~(51)带入式(42)计算得到各时刻的Fisher信息矩阵，进而得到CRLB。那么n时刻目标跟踪误差满足：
$$E\left[ {{\left| \hat{X}\left( n \right)-X\left( n \right) \right|}_{i}} \right]\ge {{\left[ J_{n}^{-1}\left( i,i \right) \right]}^{1/2}}$$ (52) 其中，Jn－1(i,i)为n时刻的Fisher矩阵的逆矩阵的第i个对角元素。
-
为了证明算法的有效性，本文通过仿真实验对运动目标进行跟踪。由于文献[12, 13]中的PKF算法仅适用于只有一个观测量的系统，且目标运动形式简单，所以将本文提出的PKF、IPKF算法与传统的EKF算法、实时定位的CRLB和滤波的CRLB进行对比。
实验采用的外辐射源为5个，其位置坐标如表 1。为了反映跟踪误差的统计特性，跟踪位置偏差均为500次蒙特卡洛仿真的均值结果。参照文献[14]，实际的测向精度在0.5°~10°；根据文献[15]，选择时差估计误差为50~250 ns。
表 1 外辐射源的位置/km
Table 1. Location of the Passive Radars/km
XTr1 XTr2 XTr3 XTr4 XTr5 x 100 -100 0 0 100 y 0 0 100 -100 100 z 0.5 0.8 0.5 0.8 0.8 假设各个时刻的DOA和TDOA的测量误差相等，且误差分别为1°和100 ns。目标的初始位置在(150,150,5)km，以(200,100,50)m/s的速度作匀速直线运动。得到利用EKF、PKF和IPKF算法实现目标跟踪的偏差如图 2，并将跟踪的误差偏差与实时定位CRLB和跟踪CRLB对比。由图 2(a)可知，滤波跟踪算法的CRLB低于实时定位的CRLB，PKF跟踪算法的跟踪精度、收敛速度均高于EKF算法。经过迭代的PKF算法的跟踪精度和收敛速度均高于PKF算法，随着迭代次数的增加，跟踪的精度和稳定性得到提高，但是随着迭代次数的增加，计算复杂度也随之提高，且跟踪的偏差小于实时定位的CRLB。图 2(a)中所示的结果即为迭代3次之后的跟踪精度，已基本稳定，所以无需再增加迭代次数。EKF的跟踪误差比较大是由其跟踪算法不稳定造成的，对其稳定性的分析将在下文给出。
图 2 匀速运动目标跟踪偏差
Figure 2. Tracking Deviation of Uniform Motion Target
假设各时刻DOA和TDOA的测量误差分别为5°和500 ns，且各时刻测量误差相等。匀速运动目标的初始条件与图 2(a)的实验条件相同。利用EKF、IPKF算法实现目标跟踪的偏差和定位跟踪算法的CRLB如图 2(b)。由图 2(b)可知，IPKF算法跟踪的精度和收敛速度均高于EKF算法，随着迭代次数的增加，算法跟踪的性能得到提高。且迭代次数达到3次之后，算法性能基本稳定。对比图 2(a)和图 2(b)的跟踪结果可知，随着观测量误差的增大，跟踪的偏差也随之增大；且跟踪的误差偏离CRLB越远。这是因为PKF算法的观测方程是通过省略测量误差高次项的条件下构建的，当测量误差比较大时，观测方程就越不准确，进而导致算法跟踪的精度降低。
在相同的实验条件下，研究不同跟踪算法的稳定性。当跟踪最后时刻的位置误差大于6 km时算法发散。假设各个时刻的DOA和TDOA误差分别为3°和300 ns，得到不同算法1 000次仿真中发散的次数统计如表 2。由图 2和表 2所示的结果可以看出，IPKF算法的稳定性高于EKF算法，且迭代次数越高，IPKF算法的稳定性越高。
表 2 (σθ=σφ=3°,στ=300 ns)
Table 2. Statistics of the Diffuse Times of Different Tracking Algorithms (σθ=σφ=3°,στ=300 ns)
EKF PKF 2 IPKF 3 IPKF 4 IPKF 5 IPKF 发散次数 257 54 15 11 4 3 仿真实验证明了IPKF算法的目标跟踪精度、收敛速度和稳定性均高于EKF算法。虽然两种算法均存在泰勒展开线性化，但是IPKF算法是对测量误差省略的线性化，而EKF算法是对目标位置误差省略的线性化。那么当估计值存在较大误差时，EKF算法的线性化方法得到的方程将不准确，进而导致跟踪的精度得不到保证，而IPKF算法可以在较差初值时得到比较稳定的跟踪结果。
-
本文研究了基于角度和时差构建无源定位跟踪模型。通过变形，构建伪线性卡尔曼滤波模型。利用传统的卡尔曼滤波算法实现运动目标跟踪，滤波初值由最小二乘算法获取，进一步分析算法模型，提出IPKF目标跟踪算法。仿真实验说明，本文提出的IPKF算法的跟踪精度、收敛速度和算法问的稳定性均优于传统的EKF算法。从而证明了PKF在联合角度和时差的单站外辐射源定位跟踪应用领域的价值。
Single-Observer Passive Tracking Based on DOA and TDOA Using Iterative Pseudo-Liner Kalman Filter
-
摘要: 随着电子战、信息战在现代军事领域的地位日趋重要，基于外辐射源的定位跟踪方法成为现代雷达领域的研究热点。针对通过单站接收多外辐射源信号获取角度（direction of arrival，DOA）和时差（time difference of arrival，TDOA）信息对运动目标跟踪的问题，首先推导角度和时差的伪线性观测方程，在通过最小二乘（least squares，LS）算法获取初值的条件下，利用传统的卡尔曼滤波算法实现目标的跟踪，该方法称为伪线性卡尔曼滤波（pseudo-liner Kalman filter，PKF）算法。进一步分析观测方程，提出了利用迭代的IPKF（iterative PKF）目标跟踪算法，并推导其克拉美罗下界（Cramer-Rao lower bound，CRLB）。仿真实验分析说明，该IPKF算法的跟踪精度、收敛速度和稳定性均优于传统的扩展卡尔曼滤波（extended Kalman filter，EKF）算法，且迭代次数越多，性能越好，观测误差越小，跟踪误差越接近CRLB。Abstract: Given the growing importance of electronic and information warfare in modern military field, localization and tracking methods based on passive coherent radar have become a research hotspot. This paper studies the problems in passive tracking based on DOA and TDOA, obtained by single-observer receiving signals from multiple passive transmitters. First of all, we deduce the pseudo-liner observation equations for DOA and TDOA. A traditional Kalman filter is applied to track target, on the condition that the least squares (LS) algorithm is used to get the initial value, a pseudo-liner Kalman filter (PKF). Furthermore, after analyzing the observation equations, an Iterative PKF (IPKF) algorithm is used to track a target. We also deduce the CRLB of the proposed algorithm. Simulations show that the tracking precision, the rate of convergence, and the tracking stability of proposed algorithm are higher than extended Kalman filter (EKF). The higher the number of iterations, the better the performance and the tracking errors will be closer to CRLB with smaller observation error.
-
Key words:
- IPKF /
- single-observer target tracking /
- passive radar /
- joint DOA and TDOA /
- CRLB
-
表 1 外辐射源的位置/km
Table 1. Location of the Passive Radars/km
XTr1 XTr2 XTr3 XTr4 XTr5 x 100 -100 0 0 100 y 0 0 100 -100 100 z 0.5 0.8 0.5 0.8 0.8 
下载: 导出CSV
表 2 (σθ=σφ=3°,στ=300 ns)
Table 2. Statistics of the Diffuse Times of Different Tracking Algorithms (σθ=σφ=3°,στ=300 ns)
EKF PKF 2 IPKF 3 IPKF 4 IPKF 5 IPKF 发散次数 257 54 15 11 4 3
下载: 导出CSV
-
[1] Farina Alfonso, Kuschel Heiner. Guest Editorial Special Issue on Passive Radar (Part I)[J]. IEEE Aerospace and Electronic Systems Magazine, 2012, 27(10):5 doi: 10.1109/MAES.2012.6373906 [2] Belfiori F, Monni S, Van Rossum W, et al. Antenna Array Characterization and Signal Processing for an FM Radio-based Passive Coherent Location Radar System[J]. IET Radar Sonar & Navig., 2012, 6(8):687-696 http://ieeexplore.ieee.org/document/6323281/ [3] Poullin D, Flecheux M. Passive 3D Tracking of Low Altitude Targets Using DVB (SFN Broadcasters)[J].IEEE Aero-Space and Electronic Systems Magazine, 2012, 27(11):36-41 doi: 10.1109/MAES.2012.6380824 [4] Falcone P, Colone F, Lombardo P. Potentialities and Challenges of WiFi-based Passive Radar[J]. IEEE Aerospace and Electronic Systems Magazine, 2012, 27(11):15-26 doi: 10.1109/MAES.2012.6380822 [5] Lingren A G, Gong K F. Position and Velocity Estimation via Bearing Observations[J].IEEE Transactions on Aero-Space and Electronic Systems, 1978, 14(4):564-577 http://cn.bing.com/academic/profile?id=430ba5588610673d870ca5ff15bbcae0&encoded=0&v=paper_preview&mkt=zh-cn [6] Shen Junyang, Molisch A F, Salmi J M. Accurate Passive Location Estimation Using TOA Measurements[J]. IEEE Transactions on Wireless Communications, 2012, 11(6):182-192 http://cn.bing.com/academic/profile?id=eec0f34d5e2561d7b8cf56f5ef93c0e2&encoded=0&v=paper_preview&mkt=zh-cn [7] Li H W, Wang J. Particle Filter for Manoeuvring Target Tracking via Passive Radar Measurements with Glint Noise[J].IET Radar Sonar & Navigation, 2012; 6(3):180-189 http://cn.bing.com/academic/profile?id=f80a2005317a50111c25dcbd7a4930a1&encoded=0&v=paper_preview&mkt=zh-cn [8] Norouzi Y, Derakhshani M. Joint Time Difference of Arrival/Angle of Arrival Position Finding in Passive Radar[J]. IET Radar Sonar & Navigation, 2009, 3(2):167-176 http://cn.bing.com/academic/profile?id=bedb908b521e2351af2850537d81b136&encoded=0&v=paper_preview&mkt=zh-cn [9] 李书进,虞晖.基于扩展卡尔曼滤波的非线性带滑移滞变系统的实时估计[J].武汉大学学报·信息科学版,2004,29(1):89-92 http://ch.whu.edu.cn/CN/abstract/abstract4698.shtml Li Shujin, Yu Hui. Identification of Nonlinear Hysteretic Systems with Slip by Extended Kalman Filter[J].Geomatics and Information Science of Wuhan University, 2004, 29(1):89-92 http://ch.whu.edu.cn/CN/abstract/abstract4698.shtml [10] 王坚,刘超,高井祥,等.基于抗差EKF的GNSS/INS紧组合算法研究[J].武汉大学学报·信息科学版,2011,36(5):596-600 http://ch.whu.edu.cn/CN/abstract/abstract555.shtml Wang Jian, Liu Chao, Gao Jingxiang, et al. GNSS/INS Tightly Coupled Navigation Model Based on Robust EKF[J].Geomatics and Information Science of Wuhan University, 2011, 36(5):596-600 http://ch.whu.edu.cn/CN/abstract/abstract555.shtml [11] Koteswara R S. Pseudo-Linear Estimator for Bearings-only Passive Target Tracking[J]. IEE Radar Sonar and Navigation, 2001, 148(1):16-22 doi: 10.1049/ip-rsn:20010144 [12] Koteswara R S. Maneuvering Target Tracking Using Pseudo Linear Estimator with Active Sonar Measurements[C]. IEEE Conference on Signal Processing, Communications and Networking, India, 2007 [13] 辛云宏,杨万海.基于伪线性卡尔曼滤波的两站红外无源定位及跟踪技术[J].西安电子科技大学学报·自然科学版,2004,31(4):505-508 http://www.cnki.com.cn/Article/CJFDTOTAL-XDKD200404004.htm Xin Yunhong, Yang Wanhai. Pseudo-Linear Kalman Filter based Passive Location and Tracking Techniques by Two Infrared Stations[J]. Journal of Xidian University, 2004, 31(4):505-508 http://www.cnki.com.cn/Article/CJFDTOTAL-XDKD200404004.htm [14] Streetly M. Jean's Radar and Electronic Warfare Systems[R]. 14th Edn. Jean's Information Group Limited, UK, 2002 [15] Petropulu A, Nikias C L, Proakis J G. Cumulant Cepstrum of FM Signals and High-Resolution Time Delay Estimation[C]. IEEE Conference on Acoustics, Speech, and Signal Processing, New York, USA, 1988 [16] Tichavsk P, Muravchik C H, Arye Nehorai. Posterior Cramer-Rao Bounds for Discrete-Time Nonlinear Filtering[J]. IEEE Transactions on Signal Processing, 1998, 46(5):1386-1396 doi: 10.1109/78.668800 -
点击查看大图
图(2) / 表(2)
计量
- 文章访问数: 2573
- HTML全文浏览量: 127
- PDF下载量: 344
- 被引次数: 0
