随着小卫星功能模块化、一箭多星和应急发射等技术的发展，在轨卫星的数量急剧增加。而传统地面测控站受布站位置的约束，无法提供连续的测控能力。一旦小卫星的导航星座信号受到干扰，而自身又不具备自主导航功能，则需要利用天基测量信息为其提供一种应急导航能力。因此，有必要研究利用多星测距对空间合作目标提供协同导航的方法^[1-2]。

本文研究的协同导航方法在本质上属于多源信息融合滤波问题，主要有集中式和分布式2种。集中式滤波要求融合中心汇集所有节点的测量信息，进行集中式数据处理，一般会对融合中心造成较高的通信和计算压力，而且一旦该中心出现故障，则容易导致系统崩溃。为了分散通信和计算压力，降低节点故障对系统的影响，分布式滤波逐渐成为研究的热点问题^[3]。Olfati-Saber^[4]建立了一种分布式卡尔曼滤波(Distributed Kalman Filter，DKF)计算框架，但该方法在实际工程应用中容易受到系统非线性因素的限制。进而，Battistelli和Chisci^[5]提出了适用于非线性系统的分布式扩展卡尔曼滤波(Distributed Extended Kalman Filter，DEKF)。然而，DEKF对非线性系统的一阶线性化近似精度较低，并且需要计算系统的雅可比矩阵。为了提高非线性系统滤波精度，学者们相继提出了无迹卡尔曼滤波(Unscented Kalman Filter，UKF)^[6]和容积卡尔曼滤波(Cubature Kalman Filter，CKF)^[7]。其中，CKF将非线性高斯加权积分分解为球面积分和径向积分，并采用一组等权值的容积点计算非线性高斯加权积分，具有比EKF和UKF更高的滤波精度和数值计算稳定性^[8]，已在工程中取得广泛的应用^[9-10]。Wang等^[11]采用正则单形变换群计算球面积分，提出了球面单形-径向容积卡尔曼滤波(Spherical Simplex-Radial Cubature Kalman Filter，SSRCKF)，并且指出球面单形准则具有比CKF中的球面准则更高的精度。然而，CKF和SSRCKF 2种方法均采用矩匹配法计算径向积分，而该方法无法保证径向积分计算的最优性。为了提高径向积分精度，Shovan和Swati^[12]提出一种容积求积分卡尔曼滤波(Cubature Quadrature Kalman Filter，CQKF)，其采用高阶高斯-拉盖尔求积分准则计算径向积分，具有比矩匹配法更高的精度。然而，上述非线性滤波需要向分布式滤波中推广。

本文将球面单形准则与高斯-拉盖尔求积分准则相结合，提出了一种新的球面单形-径向容积求积分准则，将该准则嵌入DKF框架，得到了一种适用于非线性系统的分布式球面单形-径向容积求积分卡尔曼滤波(Distributed Spherical Simplex-Radial Cubature Quadrature Kalman Filter，DSSRCQKF)算法，并将该算法应用于解决多颗卫星对空间合作目标的协同导航问题，仿真实验结果验证了本文算法的有效性。

1 分布式协同导航数学模型

用导航星表示自身具备自主导航能力的卫星，用目标星表示空间合作目标，协同导航示意图如图 1所示。s个导航星处在一个分布式通信网络中，使用无向图G=(V, E)对该通信网络进行建模^[13]。其中，V={1, 2, …, s}为导航星集合，E={(d, q)|d, q∈V}为通信链路集合。当导航星d和q之间可以通信时，(d, q)∈E，此时称导航星d和q互为邻居星。导航星d的邻居星集合用N_d表示，同时记J_d=N_d∪(d)为邻居星集与自身的并集。假设每个导航星均具备对目标星的测距能力，通过与其邻居星的分布式数据融合实现对目标星的协同一致定轨。

1.1 协同导航状态方程

地球J₂项非球形摄动是卫星所受到的最主要的摄动力。在J2000地心惯性坐标系中，考虑J₂项摄动影响，卫星的轨道动力学模型为^[14]

式中：X=(x y z v_x v_y v_z)^T为卫星的轨道状态；J₂为带谐项系数；μ为地球引力常数；R_e为地球半径；(f_x f_y f_z)^T为地球高阶非球形摄动、三体引力摄动和太阳光压摄动等摄动之和在3个坐标轴上的分量，在研究中可以等效成零均值的高斯白噪声。

用四阶龙格库塔法可以将式(1)写成如下离散状态方程的形式：

式中：X_k=(x_k y_k z_k v_x,k v_y,k v_z,k)^T为k时刻的轨道状态；w_k为系统噪声。

1.2 协同导航量测方程

使用微波测距技术可以直接测量目标星与导航星d之间的径向距离R_d，设导航星d的轨道位置为(x_d y_d z_d)^T，则测量值与轨道状态X有如下几何关系：

将式(3)的几何关系写成如下离散量测方程的形式：

式中：z_d,k为k时刻的量测值；v_d,k为量测噪声。

2 分布式协同导航滤波算法 2.1 球面单形-径向容积求积分准则

考虑积分，令x=ry，y满足y^Ty=1为单位球体表面，r≥0为球体半径，则I(f)可以分解为如下球面积分S(r)和径向积分R^[7]：

式中：U_n为单位球体表面；σ(y)为面积分元素。

一般难以得到积分S(r)和R的解析解，因此考虑采用数值积分的方法对其进行近似。由文献[11, 15]可知，可以采用如下由2n+2个积分点构成的三阶球面单形准则对球面积分进行近似：

式中：为n维单位球面的表面积，为Gamma函数；a_l=[a_l,1 a_l,2 … a_l,n]^T，l=1, 2, …, n+1为n维单形的顶点，其元素定义为

对径向积分，令r²=t，得到，进而得，令，则，采用如下高斯-拉盖尔求积分准则对R中的积分项进行近似处理^[12]：

式中：t_j为求积分点；A_j为对应的权值。求积分点可以由如下p阶切比雪夫-拉盖尔多项式^[12]的解求得：

相应的权值可以由式(11)解得：

该准则的近似精度取决于求积分点的个数，当p=2时可以得到

由式(10)、式(11)解得t₁、A₁、t₂和A₂的值为

非线性卡尔曼滤波的核心是计算高斯加权积分，N(·)为高斯概率密度函数，且变量x的均值为x，协方差矩阵为P_x。该积分具有如下恒等形式^[7]：

结合式(12)和t₁、t₂、A₁、A₂，可以得到如下计算I_N的球面单形-径向容积求积分准则：

特别的，当式(9)中的p=1时，可以解得t₁=n/2，A₁=Γ(n/2)，进而可推导出文献[11, 15]中的球面单形-径向容积准则。由此可见，本文提出的球面单形-径向容积求积分准则的精度要高于球面单形-径向容积准则。

用a=[a₁ a₂ … a_n+1]表示由a_l组成的矩阵，用[·]_i表示矩阵的第i列，利用准则式(15)构造容积点和权值ω_i分别为

2.2 DSSRCQKF算法

由状态方程式(2)和量测方程式(4)构成如下滤波系统方程，为了使算法研究更具有一般性，考虑状态向量x_k∈Rⁿ，量测向量z_k∈R^c，同时为了使公式的表述更为简洁，在没有歧义的情况下省略下标d。

式中：w_k和v_k为不相关的零均值高斯白噪声，其协方差矩阵分别为Q_k和R_k。利用统计线性误差传播方法可以得到交叉协方差矩阵P_xz的等价表达式P_xz≈P_k^- H_k^T，从而定义伪观测矩阵，将式(16)和式(17)所示的容积点和权值嵌入DEKF^[5]的更新过程可以得到如下DSSRCQKF算法，具体计算步骤如下：

步骤1  滤波初始化。

循环k=1, 2, …，完成以下步骤。

步骤2  时间更新。分别利用和代替式(16)中的x和P_x计算容积点，并利用f(·)计算其非线性传递：

按照式(17)中ω_i的取值加权合并X_k⁽ⁱ⁾得到先验状态估计值：

计算先验状态协方差矩阵P_k^-：

步骤3  量测更新。分别利用和代替式(16)中的x和P_x计算容积点，并利用h(·)计算其非线性传递：

按照式(17)中ω_i的取值加权合并Z_k⁽ⁱ⁾得到量测预测值ẑ_k：

计算交叉协方差矩阵P_xz：

步骤4  与邻居星进行信息交互。向邻居星发送信息m_k=(i_k I_k )，并从邻居星接收信息，其中，I_k和i_k定义为

式中：ε_k=z_k-ẑ_k为滤波新息；R_k为量测误差协方差矩阵。

步骤5  对接收到的信息进行融合。

步骤6  状态更新。计算后验协方差矩阵P_k⁺：

计算后验状态估计：

式中：，ε为一个小量，为矩阵的Frobenius范数。

从式(18)可以看出，本文算法并没有对非线性函数的具体形式进行约束，因此DSSRCQKF算法不仅适用于本文所述的多星对合作目标的协同导航问题，同样可以应用于传感器网络目标跟踪、编队卫星协同导航和地基测控协同定轨等非线性系统协同滤波问题。

3 仿真实验

图 2为某型号卫星的地面实验模拟器，其上运行的高精度轨道预报算法经过实际飞行任务的检验，是可靠的仿真实验数据来源。仿真中考虑4颗导航星和1颗目标星，其轨道六根数见表 1。仿真时间为1 Jul 2016 12 : 00 : 00/UTC至1 Jul 2016 13 : 30 : 00/UTC，4颗导航星间的通信拓扑关系如图 3所示，假设星间测距精度为20 m。

滤波初值为

初始协方差矩阵为

对比本文提出的DSSRCQKF算法与DEKF算法，参数ε取为0.01，用均方根误差(RMSE)描述导航星对目标星的实时定轨精度，运行200次Monte Carlo仿真，仿真结果如图 4~图 7所示。统计平均定轨RMSE，并分别列于表 2和表 3。从仿真结果可以看出，本文DSSRCQKF算法实时定位精度约为19 m，定速精度约为1.71 m/s，而DEKF算法定位精度约为30 m，定速精度约为1.73 m/s，从而表明本文算法具有更高的导航精度。对于多星协同导航问题，2种算法通过导航星间的分布式通信与数据融合，均可实现对目标星轨道状态的一致估计。每颗导航星的估计结果基本相同，结果间细微的差别主要由系统的非线性引起，因为在非线性卡尔曼滤波中，假设后验概率密度服从高斯分布，本质上是一种次优滤波方法，无法像线性卡尔曼滤波一样得到理论上的最优估计。而且为了将非线性卡尔曼滤波嵌入分布式滤波中而引入的伪观测矩阵同样会带来一些误差，但每颗导航星估计值间的差别在应用中是可以接受的。从导航星间的通信拓扑结构可以看出，每颗导航星仅与其邻居星通信，滤波中间数据在整个无线网络中分布式流动，从而避免了将所有数据发送到融合中心的集中式处理，提高了系统的灵活性。

4 结论

1) 本文将球面单形准则与高斯-拉盖尔求积分准则相结合，提出一种新的球面单形-径向容积求积分准则，并将该准则嵌入DKF框架，得到一种适用于多星对空间目标协同导航的DSSRCQKF算法。该算法要求每颗导航星仅与其邻居星进行通信和数据融合，通过数据的分布式流动实现对目标星轨道状态的一致估计，从而避免了传统集中式处理中较高的通信和计算压力。

2) 仿真实验结果表明，本文DSSRCQKF算法中各导航星对目标星的实时定位精度约为19 m，定速精度约为1.71 m/s。DEKF算法中各导航星对目标星的定位精度约为30 m，定速精度约为1.73 m/s，相比之下本文算法将定位精度提高了11 m，定速精度提高了0.02 m/s，从而验证了本文算法的有效性。

3) 本文提出的DSSRCQKF算法还可以进一步应用于其他非线性系统协同滤波问题。