林建军，王骧予涵，班晓军*，黄显林
（1.哈尔滨工业大学 航天学院，黑龙江 哈尔滨 150000；
2.哈尔滨工程大学 航天与建筑工程学院，黑龙江 哈尔滨 150001）
与有源定位相比，无源定位不需要向对方发射探测信号，具有隐蔽性良好、生存能力强的优点［1］。近年来，机载单站快速定位技术是发展最快的无源侦察定位技术之一，因其作用距离远、机动性强的特点，在军事上存在着巨大的应用价值。根据定位原理的不同，该技术可以分为到达角（Angle of Arrival，AOA），DOA，TDOA 等。随着阵列天线及其信号处理技术的发展，AOA测量的硬件平台成本降低，这使得基于AOA 的目标跟踪技术具有更广的适用范围［2］。
三维到达角（Three-Dimension Angle of Arrival，3D-AOA）目标跟踪是通过测量目标与运动站之间的方位角、俯仰角，进而对目标进行运动状态估计的技术，其测量方程具有非线性。常见的基于角度信息的目标运动分析方法可以分为伪线性滤波法与非线性滤波法［3］，其中，伪线性滤波法由于其计算量小、不易发散的优点而倍受关注。文献［4］在二维条件下对方位角进行伪线性处理，提出了适用于匀速运动目标的伪线性滤波器（Pseudo-Linear Estimator，PLE），相比于扩展卡尔曼滤波器（EKF）具有更高的稳定性，但该方法是一种有偏估计。为了解决PLE 的有偏性问题，文献［5］提出了一种辅助变量法（IV），但没有闭式解且对初值敏感；
文献［6］考虑了观测方程噪声特性提出了WIV，获得了接近最大似然估计（MLE）的估计性能。通过添加约束条件求解伪线性方程的极小值，文献［7］给出了渐近无偏的估计器，文献［8］将该方法应用到二维匀速目标跟踪上，文献［9-10］则进一步考虑观测平台的位置误差，对该方法进行了优化。在单站条件下，目标运动的可观测信息量少，易出现目标运动信息不完全可观的问题，若不能有效机动无法完成对目标的有效跟踪［11］，可以通过增加传感器数量［12］、引入速率先验条件［13］、固定角度探测范围［14］等方式，使传感器在静止条件下也能对目标进行有效跟踪。然而，基于伪线性模型的目标运动状态估计研究主要集中在小距离以及二维条件下，对远距离三维场景下的目标定位跟踪问题研究较少，且要求观测平台具有较为复杂的机动特性。在实际环境中，侦察机若处于安全状态，通常会以巡航状态工作。而借助于目标识别技术，能够对目标型号进行有效识别，并给出相应的先验速率参数。
本文在众多伪线性测量研究的基础上，针对巡航状态下的目标，提出了一种观测平台非机动条件下基于三维到达角的目标运动状态无偏估计法。该方法借助目标速率先验的条件，以相对运动模型对系统进行建模，将伪线性最小二乘方法应用到该模型中；
针对伪线性最小二乘存在有偏性的问题，通过添加约束条件对伪线性最小二乘法进行优化求解，达到消除模型有偏性的作用，从而得到一种具有渐进无偏特性的目标运动状态估计方法。为了保证算法的运行效率，算法的解为闭式解，而根据测量噪声特性对算法提出约束条件，使得问题转换为带约束的最优化问题，能够在牺牲伪线性最小二乘方法一定的运行效率后大大提高算法的目标跟踪精度。
2.1 运动模型
当目标以巡航状态工作时，运动状态近似呈匀速直线运动。即使目标出现小机动其状态方程成为时变系统，根据分段线性定常系统（Piece-Wise Const Systerm，PWCS）理论，仍可以认为目标在短时间内速度不变［15］。当观测平台以匀速运动形式巡航时，二者之间的相对运动仍保持匀速。在该场景下，定义观测平台的位置为po=[pox，poy，poz]T，速度为vo=[vox，voy，voz]T，目标的位置为pt=[ptx，pty，ptz]T，速度为vt=[vtx，vty，vtz]T。以目标与观测平台之间的相对运动作为研究对象，假设观测平台固定于点O，目标做匀速直线运动，其先验速率为vc。如图1 所示，以正东方向为X轴，正北方向为Y轴，天向为Z轴，建立右手系直角坐标系。目标与观测站连线与Y轴夹角为方位角θ，与该连线的X-Y平面投影夹角为俯仰角φ，目标的位置与速度分别用p=pt-po=[px，py，pz]T，v=vt-vo=[vx，vy，vz]T表示。
图1 东北天坐标系Fig.1 North-East-Up（ENU）coordinate system
tk表示k时刻（k=0，1，…），目标与观测站的距离为rk，那么初始时刻距离为r0。忽略过程噪声的影响，任意时刻tk的目标位置可表示为：
2.2 三维AOA 量测模型
三维空间中，k时刻方位角、俯仰角的测量方程为：
式中：假设nk，mk为服从高斯分布的测量噪声，均值为0，方差分别为，测量方位角与俯仰角之间相互独立，由此可得测量方程协方差矩阵
不失一般性，假设采样时间间隔为T，将式（1）代入式（3）中进行代数变换：
然而，实际应用过程中无法获得真值，因此用测量值代替，而会产生新的误差向量。为了形式简洁，假设初始时刻的角测量为真值，用·表示含有噪声影响的变量。
经过推导可以得到噪声向量εk为：
2.3 可观性分析
不考虑测量噪声，对方位角的测量方程进行伪线性处理。当k=1，2，…，k，k≥2 时，有：
定义伪观测矩阵：
式中N=[I2×2，02×1]T。
令M=HTθ Hθ，则有：
经推导得到M的行列式值为：
3.1 伪线性最小二乘法
利用观测信息重新构建新观测量的方法称为伪线性化，而基于伪线性化观测方程得到最小二乘解的方法称为伪线性最小二乘法（Pseudolinear Least Square，PLE）。在上一节中，忽略了方位角与俯仰角的测量噪声，证明了系统的可观测性。然而，伪线性化处理使得量测方程中的噪声特性发生变化，这会影响最小二乘法的估计精度。
将1 到k时刻的伪线性观测方程写成矩阵形式：
根据先验速率vc，进一步可以得到目标的初始距离、位置矢量、目标与观测平台距离：
式中e0=[cosφ0sinθ0，cosφ0cosθ0，sinφ0]T。
与ε均受角测量噪声的影响，从而具有相关性，这使得该估计方法并不是无偏估计，估计期望偏差为E{PLE-b}=-E{()-1}≠0。
3.2 约束总体最小二乘法
3.2.1 算法原理
基于伪线性观测方程的PLE 方法是一种有偏估计。当角测量噪声增大或者距离增大时算法性能会恶化，因此，本文提出了一种具有无偏性的约束总体最小二乘法（Constrain Total Least Square，CTLS），使得算法能够适用于远距离场景。PLE 本质上是一种以观测方程均方误差最小的方法，其损失函数表示为：
为了解决其有偏性问题，定义增广矩阵，
那么，可将损失函数J重新表示为：
由于观测噪声的存在，测量矩阵与观测量存在误差，可将其分解为真值与误差部分。于是，增广矩阵可以表示为：
更进一步，将损失函数J展开，并对等式两边取期望，得到：
从理论上看，若E［J］能够取最小值，即残差均方差E[εTε]取最小时，v为最优解。但是，由于测量误差的存在，v出现偏差。定义：
将噪声项约束为常数，即设vTΩv为任意常数c，此时将不会影响E［J］取到极值时的v，选择c=1。
于是，问题可重新被描述为：
拉格朗日极值法可处理带约束条件的极值问题，构造函数：
对方程求导并令其导数为0，得：
于是，v*，λ*即为矩阵（P，Ω）的最小广义特征向量与最小特征值，改进的渐进无偏估计b*为：
为了提高算法的运行效率，P，Ω可以由递推形式获得。
k时刻记矩阵P，Ω为P（k），Ω(k)，那么k+1时刻有：
同理可得：
算法步骤如表1 所示。
表1 约束总体最小二乘法算法步骤Tab.1 Steps of constrained total least squares algorithm
3.2.2 算法无偏性证明
将约束条件代入式（18）可得：
从理论上看，f(v，λ)的最小值也就是E[J]的最小值为1。对式（23）两边同时左乘vT，此时λ=1。将λ=1 代入式（23），并将Pv展开：
根据随机矩阵理论，随机矩阵的极限也是它的渐进期望，于是当测量数据量足够大时有：
此时向量v∈null（A），等价于：
定义b的真值为b0，显然其满足Hb0=-y。根据可观性理论可知，H为列满秩矩阵，因此b只有唯一解，即：b=b0。那么，当测量数据量足够大时，b依概率收敛于真值。
4.1 性能指标
克拉美罗下界（Cramer-Rao Lower Bound，CRLB）能够给出状态估计的理论最小方差。基于非机动条件下的观测平台对非机动目标进行运动状态估计时，角度测量误差与观测平台的位置误差是目标估计精度的重要影响因素，而本文不考虑观测平台位置误差的影响。定义状态观测量κ=[θ，φ]T，观测站与目标之间距离为r，待估计向量p=[px，py，pz]T。
克拉美罗下界等于费舍尔矩阵I的逆，即：
于是，状态估计的位置误差与距离误差的理论下界相同，可由式（35）得到：
其中：M表示蒙特卡洛仿真次数，，xk分别表示k时刻状态后验估计与真值。
时均误差可用于衡量算法在一段时间的估计性能，其定义如下：
式中：U=N-L+1，N为跟踪过程的总时刻数，L为开始记录误差的时刻数。
η用于衡量算法的提升精度，其方程如下：
式中：Index表示待评价指标参数，下标new，old分别表示初始方法与改进方法。
4.2 3D-AOA 的非机动单站运动仿真
4.2.1 不同角测量噪声下的仿真对比
仿真场景设计如下：采样时间为100 s，采样周期为0.2 s，观测站起始点坐标原点O，速度为[0，250，0]T，目标初始位置为[100，100，10]T，速度为[-340，80，-1.21]T，相对先验速率vc为380 m/s。蒙特卡洛仿真次数为200，忽略过程噪声。在光电跟踪系统的角测量误差水平能达到6 mrad［16］。因此，将观测角噪声分别假设为方差为0.1°～0.3°的高斯白噪声，且独立不相关。以相对距离误差（Relative Distance Error，RDE）与绝对位置误差（Absolue Position Error，APE）作为评价算法的定位精度的评价指标，分别为：
仿真场景下，观测站与目标运动的位置态势如图2 所示。仿真对比了PLE，IV，WIV，EKF 以及约束总体最小二乘法。其中，EKF 设置初值b0|-1=03×1，初始状态误差的协方差P0|-1=10I3×3。
图2 观测站-目标运动态势Fig.2 Observation station -target motion position
图3 比较了各算法在100 s 时，不同角度测量标准差条件下的RDE 与APE。表2 与表3 比较了不同算法相比PLE 的精度。从图3（a）可以看出，当测量噪声为0.1°时，各算法均能有效收敛到10%以内的相对距离误差。其中，约束总体最小二乘法与IV，WIV 方法误差接近。随着测量噪声逐渐增大到0.2°，各算法性能受到不同程度的影响。PLE 与CTLS 的估计误差均保持平稳上升，但PLE 的收敛精度明显低于CTLS；
EKF，IV，WIV 均出现较大的波动，性能不够稳定，EKF 的估计精度反而有所提升，这可能与初值P的选择是否恰当有关。当测量噪声增大至0.3°时，EKF 与PLE 方法在100 s 时的相对距离误差达到30% 以上，此时无法有效工作，而IV 与WIV 则出现发散倾向，RDE 与APE 均高于其他算法。其中，WIV 的绝对位置误差超过100 km，但是CTLS 在100 s 时的RDE 仍能达到10%以内。根据表2 与表3，在100 s 时CTLS 的提升效果最为稳定与明显，在角测量标准差为0.1°，0.15°，0.2°，0.25°，0.3°时，相对距离的估计精度相比于PLE 分别提高了70%，80.8%，82.9%，85.5%，85.9%；
绝对位置的估计精度相比于PLE 分别提高了 65.4%，75.5%，78.3%，81.3%，81.8%，精度的提升效果随着角测量标准差的增大逐步提升，算法的优势逐渐明显。
表2 不同角测量标准差下的RDE 及性能提升效果Tab.2 RDE and performance improvement effects under different angle measurement standard deviations（%）
表3 不同角测量标准差下的APE 及性能提升效果Tab.3 APE and performance improvement effects under different angle measurement standard deviations
图3 不同角度测量标准差下的算法性能比较Fig.3 Comparison of algorithm performance at different angles of measurement standard deviation
为了进一步精确评估算法的精度与收敛性能，计算50～100 s 的平均时间偏差，如表4 所示。可以看出，在0.1°，0.2°，0.3°时，CTLS 的相对距离误差与绝对距离误差均为最低，这说明CTLS 相比于其他算法具有较高的收敛速度与定位精度。
表4 不同角测量误差RDE 的平均时间偏差Tab.4 Time-ave deviation of RDE under different angle measurement standard deviations
4.2.2 不同初始距离条件下的仿真对比
在速度一定的条件下，目标起始点会影响角度观测角的变化，因此进一步对初始距离的影响进行探究。取单位起始点位置矢量为[100，100，5]T（km），设置不同起始点，分别为起始点位置矢量的0.5～2 倍，间隔为0.25 倍，其他条件相同，角测量标准差均设置为0.1°。图4 为算法在不同初始距离条件下的RDE 与APE。由图4 可知，当取0.5 倍单位初始距离时，除了EKF外，其他算法的相对距离误差与绝对距离误差均十分接近；
当初始距离逐渐提高至一倍初始距离时，各算法性能逐渐下降并有所区别，100 s 时按RDE 排序为PLE＞EKF＞IV＞WIV＞CTLS，APE 为EKF＞PLE＞IV＞WIV＞CTLS，各算法100 s 时均能达到10%以内的相对距离误差；
当初始距离增加到2 倍初始距离时，各种算法误差大大提高，其中IV，WIV 出现较大波动，当初始距离为1.75～2 倍单位初始距离之间时，二者的APE 均大于PLE 与EKF，已超过100 km，对目标的运动信息估计作用已失效。而PLE，EKF，CTLS 的误差均以较小幅度增加，CTLS 的误差曲线明显低于PLE，EKF。表5 与表6 比较了不同初始距离条件下不同算法相比PLE 的精度。CTLS 在初始距离分别为0.5 倍、1 倍、1.5 倍、2倍初始距离时，相对距离的估计精度相比PLE 分别提高33.3%，71.2%，80.9%，80.8%，绝对位置估计精度相比PLE 分别提高20%，65.4%，75.6%，75.8%，精度提升效果随着初始距离的增大逐步提升，且提升效果始终最佳，性能稳定。
表5 不同初始距离下的RDE 及性能提升效果Tab.5 RDE and performance improvement effects under different initial distances（%）
表6 不同初始距离下的APE 及性能提升效果Tab.6 APE and performance improvement effects under different initial distances
图4 不同初始距离的算法性能比较Fig.4 Comparison of algorithm performance for different initial distances
在50～100 s 内，不同初始距离条件下的时均RDE 与APE 如表7 所示。通过比较可以发现，当初始距离较近，取到0.5 倍单位初距时，各算法的性能接近，IV，WIV 与CTLS 能够收敛到时均RDE 为1%，APE 为1 km 以内。而初始距离增加到单位初始距离的1 倍与2 倍时，其他算法均在CTLS 性能参数的1.5 倍以上。其中，在初始距离为单位初始距离的2 倍时，其他算法的50～100 s 时均RDE 均在60%以上。
表7 不同初距RDE 平均时间偏差Tab.7 Time-ave deviation of RDE under different initial distances
最后，以EKF 的一次蒙特卡洛仿真时间作为单位时间，各算法的相对运行时间如表8 所示。可以看出，CTLS 的算法运行速度与IV 相当，而WIV 的运行效率远低于其他算法，且实验过程中其定位并不明显优于IV，与IV 一样受测量噪声与初始距离误差的影响较大，具有性能不稳定的缺点。
表8 算法相对运行时间对比Tab.8 Comparison of algorithm"s relative running time（s）
通过对比算法在不同测量误差与不同初始距离下的仿真性能，可知本文提出的CTLS 相比于PLE，EKF 具有更高的定位精度，估计效果均较为稳定，不易出现较大波动；
而相比于同样具有无偏特性的IV，WIV，当测量噪声较小或初始距离较近时，CTLS 均能达到十分接近克拉美罗下界的估计效果，但是当测量噪声较高或初始距离较远时，则具有更强的适应性。
本文研究了百公里级别的单站无源定位问题，在运动平台非机动条件下，利用先验速率的已知条件，提出了一种渐进无偏的基于3D-AOA的目标跟踪算法。以相对运动模型为研究对象，对系统的量测方程进行了伪线性处理，并对该系统的可观性进行了分析。然后，介绍了基于伪线性测量方程的PLE 方法，针对PLE 的有偏问题，提出了一种具有渐进无偏性能的CTLS，对该算法原理进行了推导，并证明了算法的无偏性。最后，通过仿真分析在0.1°，0.2°，0.3°测量标准差时发现，CTLS 在50～100 s 内的平均RDE 分别为6%，12%，21%，平均APE 为9，19，35 km。同时，还比较了在不同初始距离条件下的算法性能，在50 km 左右的仿真场景中各算法性能接近，时均RDE 均能达到5%以内；
而当距离逐渐提高到200 km 级别时，其他算法逐渐失效，而CTLS 在100 s 内已能明显看出收敛趋势，时均APE 为30%，明显优于其他算法，100 s 的APE已能达到10%以内；
比较各算法相对PLE 方法的精度提升效果发现，CTLS 的精度提升效果最好、最稳定。仿真实验结果表明，CTLS 与其他具有无偏性的算法相比，对角测量误差与初始距离的抗干扰能力最强，且定位精度高，性能稳定，不易发生较大波动，同时兼具收敛速度快的优点。而且算法的运行速度与EKF 以及PLE 同一个量级，运行效率较高，因此特别适合用于远距离场景、角测量误差较大时的单站定位跟踪问题。