基于高斯-牛顿迭代的三星时差定位融合算法

2020-12-09徐英杰郭福成

航天电子对抗 2020年5期

徐英杰，郭福成

（国防科技大学电子科学学院，湖南长沙410000）

0 引言

三星时差定位系统通过三个卫星构成编队接收辐射源信号到达各颗卫星的时间差，结合辐射源在地球表面这一先验知识即可求出其位置[1]。它在电磁目标监视、海上救援搜索等方面具有重要的应用价值。

对于地球上的一个静止辐射源目标，卫星编队可能会多次接收到其辐射的信号，在过顶期间对辐射源进行多次观测，有效利用多次观测的数据进行融合定位可以提高定位精度并节省增加卫星数量的成本。Zhang 等人[2]提出了一种仅基于TDOA 测量和不同阶段结果的数据融合模型，该模型由Chan 算法[3]计算得出，实验表明当辐射源估计被非视线误差和不同的测量误差干扰时，其定位精度优于Chan 算法。吕明等人[4]针对传统时差定位会出现定位模糊的问题，提出了应用独立测量数据融合的加权最小二乘算法，将多个观测副站测量的独立辐射源位置信息进行融合，充分利用了各个观测副站提供的冗余量测信息，实验证明其算法可提高整个定位系统的定位精度。许丞梁[5]将到达时间差和到达频率差（FDOA）、到达角（AOA）、高程等多个量测融合来进行辐射源的定位，实验表明这些冗余参数可以提高对辐射源的定位精度，扩大高精度定位的范围。

对于低轨卫星而言，卫星在不同时刻下的定位精度是不同的，在工程应用中经常将整轨数据的所有定位点进行合批处理来提高定位精度。目前常用到的位置合批处理方法有选取主星星下点离辐射源最近时的定位结果；所有单次定位结果的直接平均和所有单次定位结果的加权平均。然而这些方法从理论上来讲都不是最优的，选取主星星下点离辐射源最近时的定位结果没有用到整轨的星历点数据；直接平均的结果没有考虑到卫星在不同时刻下的定位精度差异；加权平均方法所用到的加权矩阵是单次定位得到的结果，理论上不是最优估计。为此本文提出一种将多次观测的时差数据与高程信息融合的三星时差定位算法，在使用相同数据的条件下，获得比现有位置合批处理方法更高的定位精度。

1 定位数学模型

使用三颗卫星对地球上的辐射源目标uo=[ x,y,z]T分别在N 个不同观测时刻下进行时差和高程融合定位。卫星位置会受到随机误差的影响，设n时刻时的三星位置为si,n=+Δsi,n(i=1,2,3;n=1,2,…,N)，其中= []T是卫星i 真实但未知的位置，Δsi,n是相应的随机误差。由于三星在N个不同时刻下对辐射源目标进行了观测，令s=[]T表示不同时刻下的卫星位置，卫星位置误差为Δs=s-so，其中so=[]T，定义Δs 为零均值，协方差矩阵为Qs的高斯随机向量。

以卫星1 作为参考卫星进行观测，则三星在n 时刻下观测的卫星i(i=2,3)与1 之间的TDOA 定义为di1,n，即：

式中，c 为已知的信号传播速度，ri1,n=cdi1,n是此时刻下两星间的到达距离差（range-difference-of-arrival，RDOA）为此时刻下真实的RDOA,ωi1,n/c 为此时刻的TDOA 测量噪声为此时刻辐射源目标uo与卫星1 之间的真实距离为此时刻辐射源与卫星i之间的真实距离，为：

为便于计算，将获得的TDOA 与c 相乘得到距离差测量，即n 时刻下观测的距离差量测为mt,n=[r21,n,r31,n]T，其测量误差为Δmt,n=mt,n-，其中= []T，定义Δmt,n为零均值，协方差矩阵为Qt,n的高斯随机向量。则N 次观测下的距离差量测表达为mt= []T=+Δmt，是距离差观测的真值，假设不同时刻之间的观测相互独立，则Δmt是均值为零，协方差为Qt=diag{Qt,1,…,Qt,n,…,Qt,N}的高斯随机向量。本文将辐射源高程信息也融合进观测方程，其表示为mh=+Δmh，其中mh= (R+Δh)2，R 为地球半径和辐射源真实高程之和，高程误差Δh 为零均值，方差为Qh的高斯随机变量。已知正球面模型下的地球方程为：

忽略二阶误差项，则Δmh=2RΔh 为零均值，方差为4R2Qh的高斯随机变量。因此其距离差和高程信息融合的量测表达为m= [,mh]T=mo+Δm，mo为观测的真值，假设时差观测和高程观测相互独立，则Δm 是均值为零，协方差为的高斯随机向量，本文假设卫星位置误差Δs 和观测误差Δm 相互独立。

2 位置合批的融合算法

2.1 CRLB

融合算法的CRLB 为：

具体推导详见文献[6]，这里不再赘述。

2.2 现有的位置合批处理

目前工程上广泛应用的合批处理方法主要有三种：一是取主星星下点距离辐射源位置最近时的定位结果。原理为取N 个不同时刻的单次定位结果，将n时刻的单次定位结果设为ûn(n=1,2,…,N )，将n 时刻下的主星星下点到辐射源的距离设为dn(n=1,2,…,N )，则主星星下点距离辐射源位置最近时的定位结果表达式为：

二是取N 个不同时刻的单次定位结果直接平均,其表达式为：

三是取N 个不同时刻的单次定位结果进行加权平均，其表达式为：

式中，Wn为n 时刻的定位误差协方差矩阵，n=1,2,…,N。由于单次定位时辐射源真实位置未知，无法计算每次的定位误差协方差矩阵，因此可以根据单次定位估计的结果,来求得CRLB 近似作为辐射源定位误差协方差矩阵Wn来代入计算，其CRLB 计算原理同本文2.1 部分，只不过用代替辐射源真实位置uo，用带误差的卫星位置s 代替真实卫星位置so。

2.3 时差定位融合算法

由于辐射源和卫星位置都是未知的，因此需要联合估计uo和so，为简化这一问题，文献[7]给出了一种仅需估计uo的加权最小二乘算法，其表达式为：

加权矩阵W =Qm+DoQsDoT，其中Do的形式与∂mo/∂so相同，只不过要用带误差的卫星位置s 代替真实卫星位置so。根据文献[8]可知，求解非线性最小二乘问题可采用Gauss-Newton 迭代，假设第k 次迭代辐射源的位置估计为，则第k+1 次迭代的结果为：

迭代初值可选取N 个观测时刻中的任一时刻得到的单次定位结果，其单次定位求解可参考文献[9]。

以上求解都是建立在正球面模型下，当辐射源位于高纬度时，采用正球面会有较大的误差，因此可以通过球面迭代得到WGS-84 椭球模型下的结果。计算步骤如下：

1) 根据正球面模型下的定位方法得到辐射源的位置估计，并转换为经纬高大地坐标，从而得到辐射源的纬度初始值；

2)利用公式(12)更新当地地球半径radius

式中,a=6 378 137 m 为地球半长轴，e 为第一偏心率，e2=0.006 694 379 990 13，B 为辐射源的纬度；

3) 重复前2 步，直到|radiusi-radiusi-1| ＜ε（ε 表示一个足够小的值）。

3 仿真分析

本节通过实验仿真验证上述理论。假设辐射源位置的经度和纬度为117.8°W 和19.6°N，辐射源位于海平面，高程为零。用于实验的卫星星历点数据基于STK 软件生成，三星的轨道高度设为1 100 km，轨道倾角为63.4°，星间距为50 ～110 km，轨道历元时刻从2020 年8 月29 日4 时开始，将采样间隔设置为1 s，整轨取960 个星历点数据作为样本。

令单次观测的时差误差协方差阵Qt,n为：

式中，δt表示TDOA 噪声的标准差，δh表示高程误差，c=3×108m/s，由此可得 Qm=diag{Qt,1,…,Qt,n,…,Qt,N,4R2Qh}为2N +1 维方阵。卫星位置误差协方差阵 Qs为 Qs=diag{Qs,1,…,Qs,n,…,Qs,N}，其中Qs,n为：

式中，δs为卫星位置噪声的标准差，I 为9×9 的单位矩阵。

令δt=0.2 μs，δs=50 m，δh=200 m，利用采样所得到的星历点数据进行定位。在计算本文所提算法时，选择采样数据中的第一个星历点数据进行单次定位来作为迭代初值，其余三种位置合批处理均使用采样所得星历点数据进行融合定位。为方便表示，称主星星下点距离辐射源位置最近时的定位结果为方法1；位置直接平均结果为方法2；位置加权平均结果为方法3。

接着通过计算定位误差的均方根误差（RMSE）来比较本文提出的算法和工程上常用的合批处理方法之间的性能。这里融合次数N 分别取5、30、60、120、240、320、480 和960 次，即对得到的960 个星历样本每间隔960/N s 进行一次采样，用采样所得星历点数据进行融合定位，比较不同融合次数下各种合批方法定位精度，实验结果如图1 所示。

通过实验结果可以得出，在相同条件下，本文算法的性能最好，其均方根误差达到了CRLB，方法3 的性能其次，随着融合次数的增加，方法2 的性能逐渐优于方法1，这与理论上的分析一致。本文算法是利用原始时差量测，通过非线性最小二乘估计得到辐射源位置的最优估计。综上所述，本文提出的算法定位精度要优于现有的各种位置合批处理方法。