CN102305630B - 基于扩展卡尔曼滤波的sar卫星自主定轨方法 - Google Patents

Abstract

基于扩展卡尔曼滤波的SAR卫星自主定轨方法，属于卫星自主定轨技术领域。现有卫星自主定轨技术主要采用GPS，由于GPS是一种半自主定轨方式，因此存在自主性和实时性差的问题。本发明中的定轨方法就可以解决这些问题。本发明的方法是：步骤一：建立基于轨道动力学的卫星运动方程，进而得到状态方程；步骤二：建立以SAR到地面标识点之间的距离和SAR与地面标识点之间的多普勒频移为观测量的观测方程；步骤三：建立扩展卡尔曼滤波的递推方程；步骤四：得到卫星的状态信息，即得到卫星位置和速度矢量

。本发明的方法提高了卫星的定轨精度，具有自主性和实时性的特点，能够实现低轨卫星的高精度实时自主定轨。

Description

基于扩展卡尔曼滤波的SAR卫星自主定轨方法

技术领域

本发明涉及一种SAR卫星自主定轨方法，属于卫星自主定轨技术领域。

背景技术

卫星的自主定轨是指不依赖地面系统的支持下，利用卫星自带的测量元件来确定卫星在惯性空间中的位置和速度，实现自主定轨。对于卫星来说，自主定轨有利于降低卫星对地面的依赖程度，尤其是在战争状态下，在地面系统遭受破坏时，仍能保持正常运行，这对军用卫星来讲，尤为重要；另外，自主定轨还可以减轻地面系统的负担，降低地面支持成本。

现有卫星自主定轨技术主要采用GPS，但是GPS严格来讲不是一种完全自主的定轨方式，它在地面有控制站点，工作人员通过这些控制站点对GPS星的轨道、时间等数据进行更新，它是一种半自主定轨方式，因此存在自主性和实时性差问题。

发明内容

本发明的目的是提供一种基于扩展卡尔曼滤波的SAR卫星自主定轨方法，现有卫星自主定轨技术主要采用GPS，由于GPS是一种半自主定轨方式，因此存在自主性和实时性差的问题。本发明中的定轨方法就可以解决这些问题。

本发明提出的基于扩展卡尔曼滤波的SAR卫星自主定轨方法，这种方法不需要地面测控站的支持，是一种完全自主的实时定轨方法。SAR为星载合成孔径雷达的缩写。

本发明解决上述问题采取的技术方案是：

基于扩展卡尔曼滤波的SAR  星自主定轨方法，所述方法包括以下步骤：

步骤一：建立基于轨道动力学的卫星运动方程，进而得到状态方程；

建立包含J₂项的卫星轨道动力学方程，在J2000惯性坐标系下，卫星的轨道动力学方程为：

$\stackrel{\·}{X}=F[X\left(t\right),t]---\left(1\right)$

写成状态方程，即：

$\left\{\begin{array}{c}\frac{\mathrm{dx}}{\mathrm{dt}}=v_x\\ \frac{\mathrm{dy}}{\mathrm{dt}}=v_y\\ \frac{\mathrm{dz}}{\mathrm{dt}}=v_z\\ \frac{{\mathrm{dv}}_x}{\mathrm{dt}}=-\mathrm{\μ}\frac{x}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-1.5)]+{\mathrm{\ΔF}}_x\\ \frac{{\mathrm{dv}}_y}{\mathrm{dt}}=-\mathrm{\μ}\frac{y}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-1.5)]+{\mathrm{\ΔF}}_y\\ \frac{{\mathrm{dv}}_z}{\mathrm{dt}}=-\mathrm{\μ}\frac{z}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-4.5)]+{\mathrm{\ΔF}}_z\end{array}\right.---\left(2\right)$

式中，

上述式(2)简写为

状态矢量X＝[x，y，z，v_x，v_y，v_z]^T，x，y，z，v_x，v_y，v_z分别为卫星在惯性坐标系的X，Y，Z方向上的位置和速度；μ为地心引力常数；r为卫星位置参数矢量；J₂为摄动系数；ΔF_x，ΔF_y，ΔF_z分别为X，Y，Z方向上的非球形摄动的高阶项及日、月引力摄动，太阳光压摄动以及大气阻力摄动未建模摄动的影响；上角标T为转置；R_e为地球半径；

在所述简化模型中，上述这些摄动的影响用系统噪声w(t)来表示，

上述式(2)为连续方程，将其离散化，得到状态方程，即：

X(k+1)＝Φ_k+1，kX(k)+W(k)    k＝1，2，3...(3)

式中，Φ_k+1，k为k至k+1时刻的一步状态转移矩阵；X(k)为第k时刻的状态；W(k)为系统噪声矩阵，对于白噪声有：

E[W(k)]＝0，E[W(k)W(j)]＝Q_kδ_kj    j＝1，2，3...(4)

式中，Q_k为系统噪声的方差矩阵；δ：表示δ函数；

步骤二：建立以SAR到地面标识点之间的距离和SAR与地面标识点之间的多普勒频移为观测量的观测方程；

当SAR的成像视场中出现地面标识点时，假设L_i为第i个地面标识点的经度，λ_i为第i个地面标识点的纬度；当i大于等于3时，就可以确定卫星的轨道；所述第i个地面标识点在地固坐标系下的坐标为：

$\left\{\begin{array}{c}{x}_{\mathrm{di}}=R_e\cos {\mathrm{\λ}}_i\cos L_i\\ y_{\mathrm{di}}=R_e\cos {\mathrm{\λ}}_i\sin L_i\\ z_{\mathrm{di}}=R_e\sin {\mathrm{\λ}}_i\end{array}\right.---\left(5\right)$

式中，R_e为地球半径；根据地固坐标系和惯性坐标系的转换关系，得到第i个地面标识点在惯性坐标系下的坐标，即：

$\left[\begin{array}{c}{x}_{\mathrm{mi}}\\ y_{\mathrm{mi}}\\ z_{\mathrm{mi}}\end{array}\right]=\left[\begin{array}{ccc}\cos \mathrm{\Ω}& -\sin \mathrm{\Ω}& 0\\ \sin \mathrm{\Ω}& \cos \mathrm{\Ω}& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_{\mathrm{di}}\\ y_{\mathrm{di}}\\ z_{\mathrm{di}}\end{array}\right]---\left(6\right)$

第i个地面标识点在惯性坐标系下的速度，即：

$\left[\begin{array}{c}{v}_{\mathrm{mxi}}\\ v_{\mathrm{myi}}\\ v_{\mathrm{mzi}}\end{array}\right]=\left[\begin{array}{ccc}0& -{\mathrm{\ω}}_e& 0\\ {\mathrm{\ω}}_e& 0& 0\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}{x}_{\mathrm{mi}}\\ y_{\mathrm{mi}}\\ z_{\mathrm{mi}}\end{array}\right]---\left(7\right)$

式中：ω_e为地球自转角速度；Ω＝Ω₀+ω_e(t-t₀)为t时刻的恒星时；Ω₀为t₀时刻的恒星时；

根据上述式(6)和(7)，得到卫星到第i个地面标识点的距离方程，即：

${\mathrm{\ρ}}_i=\sqrt{{(x-x_{\mathrm{mi}})}^2+{(y-y_{\mathrm{mi}})}^2+{(z-z_{\mathrm{mi}})}^2}+V_1,i=\mathrm{1,2,3}---\left(8\right)$

利用上述式(8)得到多普勒频移方程，即：

$f_i=-\frac{2}{{\mathrm{\λ\ρ}}_i}[(v_x-v_{\mathrm{mxi}})\left({x-x}_{\mathrm{mi}}\right)+(v_y-v_{\mathrm{myi}})(y-y_{\mathrm{mi}})+(v_z-v_{\mathrm{mzi}})(z-z_{\mathrm{mi}})]+V_2---\left(9\right)$

式中，i＝1，2，3

利用上述式(8)和(9)得到观测方程，即：

$G=\left[\begin{array}{c}\mathrm{\ρ}\\ f\end{array}\right]=h[X\left(t\right),t]+V\left(t\right)---\left(10\right)$

式(10)中的ρ＝[ρ₁，ρ₂，ρ₃]^T，f＝[f₁，f₂，f₃]^T，V＝[V₁，V₂]^T

式中：V为测量噪声；G为观测量；

步骤三：建立扩展卡尔曼滤波的递推方程；

设定k时刻的滤波值为将步骤一中的式(3)中的非线性向量函数F围绕滤波值

展开成泰勒级数，并略去二次以上项，得到线性化方程：

$X(k+1)=f[\hat{X}\left(k\right|k),k]+\frac{\∂F[X\left(k\right|k),k]}{\∂X\left(k\right)}{|}_{X\left(t\right)=\hat{X}\left(k\right|k)}[X\left(k\right)-\hat{X}\left(k\right|k)]+\mathrm{\Γ}(X\left(k\right),k)W\left(k\right)---\left(11\right)$

式中，Γ为系数矩阵；

将步骤二中的式(10)中的非线性向量函数h离散化后，并围绕滤波值

展开成泰勒级数，并略去二次以上项，得到离散化的观测方程，即：

$G(k+1)=h[\hat{X}(k+1|k),k+1]+\frac{\∂h[\hat{X}(k+1),k+1]}{\∂X(k+1)}{|}_{X\left(t\right)=\hat{X}(k+1|k)}[X(k+1)-\hat{X}(k+1|k)]+V(k+1)---\left(12\right)$

将上述式(11)和(12)代入标准的卡尔曼滤波方程中，得到扩展卡尔曼滤波的递推方程，即：

$\hat{X}(k+1|k+1)=F[\hat{X}\left(k\right|k),k]+K(k+1)\{G(k+1)-h[\hat{X}(k+1|k),k+1\left]\right\}---\left(13\right)$

其中：

K(k+1)＝P(k+1|k)H^T(k+1)[H(k+1)P(k+1|k)H^T(k+1)+R(k+1)]^-1    (14)

P(k+1|k)＝Φ(k+1，k))P(k|k)Φ^T(k+1，k)+Q(k)                (15)

P(k+1|k+1)＝[I-K(k+1)H(k+1)]P(k+1|k)                       (16)

滤波初始值为：

$\hat{X}\left(0\right|0)=E\left[X\left(0\right)\right]$

P(0|0)＝P_x(0)

上述式(15)中一步状态转移矩阵

$\mathrm{\Φ}(k+1,k)=\frac{\∂F[X\left(k\right),k}{\∂X\left(k\right)}{|}_{X\left(k\right)=\hat{X}\left(k\right|k)}$

上述式(16)中的观测矩阵

$H(k+1)=\frac{\∂h[X(k+1,k+1)}{\∂X(k+1)}{|}_{X\left(k\right)=\hat{X}(k+1|k)}---\left(17\right)$

式中，K为增益矩阵；P(k|k)为k时刻的协方差矩阵；P(k+1|k)为协方差矩阵的一步预测；Q为系统噪声的协方差矩阵；R为测量噪声的协方差矩阵；I为单位矩阵；

步骤四：得到卫星的状态信息；

利用步骤三建立的扩展卡尔曼滤波的递推方程得到卫星的状态信息，所述卫星的状态信息即为卫星的状态估计

所述卫星的状态估计

包括卫星位置和速度矢量

本发明的有益效果是：

(1)本发明的方法是一种不依赖地面设备的完全自主定轨方法，这种方法利用合成孔径雷达天线发送的脉冲信号测得卫星到地面标识点的距离，同时根据多普勒测速原理测得卫星相对地面标识点的相对速度，从而确定卫星的轨道。其基本原理是利用合成孔径雷达测量信息结合卫星轨道动力学方程，利用最优估计的方法估计卫星的位置和速度。本发明提出的方法只需人为建立地面标识点，不需要人为参与维护，是一种完全自主定轨方法。

(2)本发明中的地面标识点是事先设计好，对其在地固坐标系中的位置信息建库并存在星载计算机上，一旦SAR观测到地面标识点，即可进行识别，然后利用本发明提出的方法对卫星的位置和速度进行估计，所以本发明中提出的自主定轨方法具有实时定轨的能力。

(3)本发明提出的自主定轨方法还可以兼顾应用。对于有效载荷为合成孔径雷达的卫星，如对面进行侦查的军用卫星，除了对地面进行之外，不需要其他的辅助手段就可以完成卫星的自主定轨，达到兼顾应用的目的。

综上，本发明的方法适应近地卫星的自主实时定轨，能够实现低轨卫星的高精度实时自主定轨。

附图说明

图1是本发明的流程图；图2是本发明中的测量信息(SAR到地面标识点的距离和SAR与地面标识点之间的多普勒频移)示意图；图2中(x_m1，y_m1，z_m1)，(z_m2，y_m2，z_m2)，(x_m3，y_m3，z_m3)为地面标识点的坐标，ρ₁，ρ₂，ρ₃为SAR卫星到地面标识点的距离。

具体实施方式

具体实施方式一：结合图1说明；本实施方式的基于扩展卡尔曼滤波的SAR卫星自主定轨方法，所述方法包括以下步骤：

步骤一：建立基于轨道动力学的卫星运动方程，进而得到状态方程；

建立包含J₂项的卫星轨道动力学方程，在J2000惯性坐标系下，卫星的轨道动力学方程为：

$\stackrel{\·}{X}=F[X\left(t\right),t]---\left(1\right)$

写成状态方程，即：

$\left\{\begin{array}{c}\frac{\mathrm{dx}}{\mathrm{dt}}=v_x\\ \frac{\mathrm{dy}}{\mathrm{dt}}=v_y\\ \frac{\mathrm{dz}}{\mathrm{dt}}=v_z\\ \frac{{\mathrm{dv}}_x}{\mathrm{dt}}=-\mathrm{\μ}\frac{x}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-1.5)]+{\mathrm{\ΔF}}_x\\ \frac{{\mathrm{dv}}_y}{\mathrm{dt}}=-\mathrm{\μ}\frac{y}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-1.5)]+{\mathrm{\ΔF}}_y\\ \frac{{\mathrm{dv}}_z}{\mathrm{dt}}=-\mathrm{\μ}\frac{z}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-4.5)]+{\mathrm{\ΔF}}_z\end{array}\right.---\left(2\right)$

式中，

上述式(2)简写为

状态矢量X＝[x，y，z，v_x，v_y，v_z]^T，x，y，z，v_x，v_y，v_z分别为卫星在惯性坐标系的X，Y，Z方向上的位置和速度；μ为地心引力常数；r为卫星位置参数矢量；J₂为摄动系数；ΔF_x，ΔF_y，ΔF_z分别为X，Y，Z方向上的非球形摄动的高阶项及日、月引力摄动，太阳光压摄动以及大气阻力摄动未建模摄动的影响；上角标T为转置；R_e为地球半径；

在所述简化模型中，上述这些摄动的影响用系统噪声w(t)来表示，

上述式(2)为连续方程，将其离散化，得到状态方程，即：

X(k+1)＝Φ_k+1，kX(k)+W(k)    k＝1，2，3...(3)

式中，Φ_k+1，k为k至k+1时刻的一步状态转移矩阵；X(k)为第k时刻的状态；W(k)为系统噪声矩阵，对于白噪声有：

E[W(k)]＝0，E[W(k)W(j)]＝Q_kδ_kj    j＝1，2，3...(4)

式中，Q_k为系统噪声的方差矩阵；δ：表示δ函数；

步骤二：建立以SAR到地面标识点之间的距离和SAR与地面标识点之间的多普勒频移为观测量的观测方程(见图2)；

当SAR的成像视场中出现地面标识点时，假设L_i为第i个地面标识点的经度，λ_i为第i个地面标识点的纬度；当i大于等于3时，就可以确定卫星的轨道；所述第i个地面标识点在地固坐标系(所述地固坐标系是指：Z轴指向北极，X轴指向格林尼治子午面和赤道平面的交点，Y轴与X轴和Z轴满足右手定则)下的坐标为：

$\left\{\begin{array}{c}{x}_{\mathrm{di}}=R_e\cos {\mathrm{\λ}}_i\cos L_i\\ y_{\mathrm{di}}=R_e\cos {\mathrm{\λ}}_i\sin L_i\\ z_{\mathrm{di}}=R_e\sin {\mathrm{\λ}}_i\end{array}\right.---\left(5\right)$

式中，R_e为地球半径；根据地固坐标系和惯性坐标系的转换关系，得到第i个地面标识点在惯性坐标系下的坐标，即：

$\left[\begin{array}{c}{x}_{\mathrm{mi}}\\ y_{\mathrm{mi}}\\ z_{\mathrm{mi}}\end{array}\right]=\left[\begin{array}{ccc}\cos \mathrm{\Ω}& -\sin \mathrm{\Ω}& 0\\ \sin \mathrm{\Ω}& \cos \mathrm{\Ω}& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_{\mathrm{di}}\\ y_{\mathrm{di}}\\ z_{\mathrm{di}}\end{array}\right]---\left(6\right)$

第i个地面标识点在惯性坐标系下的速度，即：

$\left[\begin{array}{c}{v}_{\mathrm{mxi}}\\ v_{\mathrm{myi}}\\ v_{\mathrm{mzi}}\end{array}\right]=\left[\begin{array}{ccc}0& -{\mathrm{\ω}}_e& 0\\ {\mathrm{\ω}}_e& 0& 0\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}{x}_{\mathrm{mi}}\\ y_{\mathrm{mi}}\\ z_{\mathrm{mi}}\end{array}\right]---\left(7\right)$

式中：ω_e为地球自转角速度；Ω＝Ω₀+ω_e(t-t₀)为t时刻的恒星时；Ω₀为t₀时刻的恒星时；

根据上述式(6)和(7)，得到卫星到第i个地面标识点的距离方程，即：

${\mathrm{\ρ}}_i=\sqrt{{(x-x_{\mathrm{mi}})}^2+{(y-y_{\mathrm{mi}})}^2+{(z-z_{\mathrm{mi}})}^2}+V_1,i=\mathrm{1,2,3}---\left(8\right)$

利用上述式(8)得到多普勒频移方程，即：

$f_i=-\frac{2}{{\mathrm{\λ\ρ}}_i}[(v_x-v_{\mathrm{mxi}})\left({x-x}_{\mathrm{mi}}\right)+(v_y-v_{\mathrm{myi}})(y-y_{\mathrm{mi}})+(v_z-v_{\mathrm{mzi}})(z-z_{\mathrm{mi}})]+V_2---\left(9\right)$

式中，i＝1，2，3

利用上述式(8)和(9)得到观测方程，即：

$G=\left[\begin{array}{c}\mathrm{\ρ}\\ f\end{array}\right]=h[X\left(t\right),t]+V\left(t\right)---\left(10\right)$

式(10)中的ρ＝[ρ₁，ρ₂，ρ₃]^T，f＝[f₁，f₂，f₃]^T，V＝[V₁，V₂]^T

式中：V为测量噪声；G为观测量；

步骤三：建立扩展卡尔曼滤波的递推方程；

设定k时刻的滤波值为

将步骤一中的式(3)中的非线性向量函数F围绕滤波值

展开成泰勒级数，并略去二次以上项，得到线性化方程：

$X(k+1)=f[\hat{X}\left(k\right|k),k]+\frac{\∂F[X\left(k\right|k),k]}{\∂X\left(k\right)}{|}_{X\left(t\right)=\hat{X}\left(k\right|k)}[X\left(k\right)-\hat{X}\left(k\right|k)]+\mathrm{\Γ}(X\left(k\right),k)W\left(k\right)---\left(11\right)$

式中，Γ为系数矩阵；

将步骤二中的式(10)中的非线性向量函数h离散化后，并围绕滤波值

展开成泰勒级数，并略去二次以上项，得到离散化的观测方程，即：

$G(k+1)=h[\hat{X}(k+1|k),k+1]+\frac{\∂h[\hat{X}(k+1),k+1]}{\∂X(k+1)}{|}_{X\left(t\right)=\hat{X}(k+1|k)}[X(k+1)-\hat{X}(k+1|k)]+V(k+1)---\left(12\right)$

将上述式(11)和(12)代入标准的卡尔曼滤波方程中，得到扩展卡尔曼滤波的递推方程，即：

$\hat{X}(k+1|k+1)=F[\hat{X}\left(k\right|k),k]+K(k+1)\{G(k+1)-h[\hat{X}(k+1|k),k+1\left]\right\}---\left(13\right)$

其中：

K(k+1)＝P(k+1|k)H^T(k+1)[H(k+1)P(k+1|k)H^T(k+1)+R(k+1)]^-1    (14)

P(k+1|k)＝Φ(k+1，k)P(k|k)Φ^T(k+1，k)+Q(k)                 (15)

P(k+1|k+1)＝[I-K(k+1)H(k+1)]P(k+1|k)                       (16)

滤波初始值为：

$\hat{X}\left(0\right|0)=E\left[X\left(0\right)\right]$

P(0|0)＝P_x(0)

上述式(15)中一步状态转移矩阵

$\mathrm{\Φ}(k+1,k)=\frac{\∂F[X\left(k\right),k}{\∂X\left(k\right)}{|}_{X\left(k\right)=\hat{X}\left(k\right|k)}$

上述式(16)中的观测矩阵

$H(k+1)=\frac{\∂h[X(k+1,k+1)}{\∂X(k+1)}{|}_{X\left(k\right)=\hat{X}(k+1|k)}---\left(17\right)$

式中，K为增益矩阵；P(k|k)为k时刻的协方差矩阵；P(k+1|k)为协方差矩阵的一步预测；Q为系统噪声的协方差矩阵；R为测量噪声的协方差矩阵；I为单位矩阵；

步骤四：得到卫星的状态信息；

利用步骤三建立的扩展卡尔曼滤波的递推方程得到卫星的状态信息，所述卫星的状态信息即为卫星的状态估计

所述卫星的状态估计

包括卫星位置和速度矢量

本实施方式中，在步骤一完成后，步骤二运行之前，在卫星在轨运行时，需要利用合成孔径雷达对地面成像(属于现有技术)。

在步骤一中建立卫星的轨道动力学方程之前，首先要人为的建造多个位置已知的地面标识点，并将所述多个地面标识点的位置信息建成库，储存在星载计算机中。

本发明说明书中未作详细解释的内容属于本领域专业技术人员公知的现有技术。

Claims (1)

1.基于扩展卡尔曼滤波的SAR卫星自主定轨方法，其特征在于：所述方法包括以下步骤：

步骤一：建立基于轨道动力学的卫星运动方程，进而得到状态方程；

建立包含J₂项的卫星轨道动力学方程，在J2000惯性坐标系下，卫星的轨道动力学方程为：

$\stackrel{\·}{X}=F[X\left(t\right),t]---\left(1\right)$

写成状态方程，即：

$\left\{\begin{array}{c}\frac{\mathrm{dx}}{\mathrm{dt}}=v_x\\ \frac{\mathrm{dy}}{\mathrm{dt}}=v_y\\ \frac{\mathrm{dz}}{\mathrm{dt}}=v_z\\ \frac{{\mathrm{dv}}_x}{\mathrm{dt}}=-\mathrm{\μ}\frac{x}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-1.5)]+\mathrm{\Δ}{F}_x\\ \frac{{\mathrm{dv}}_y}{\mathrm{dt}}=-\mathrm{\μ}\frac{x}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-1.5)]+\mathrm{\Δ}{F}_y\\ \frac{{\mathrm{dv}}_z}{\mathrm{dt}}=-\mathrm{\μ}\frac{x}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-4.5)]+\mathrm{\Δ}{F}_z\end{array}\right.---\left(2\right)$

式中，

上述式(2)简写为状态矢量X＝[x，y，z，v_x，v_y，v_z]^T，x，y，z，v_x，v_y，v_z分别为卫星在惯性坐标系的X，Y，Z方向上的位置和速度；μ为地心引力常数；r为卫星位置参数矢量；J₂为摄动系数；ΔF_x，ΔF_y，ΔF_z分别为X，Y，Z方向上的非球形摄动的高阶项及日、月引力摄动，太阳光压摄动以及大气阻力摄动未建模摄动的影响；上角标T为转置；R_e为地球半径；

在所述简化模型中，上述这些摄动的影响用系统噪声w(t)来表示，

上述式(2)为连续方程，将其离散化，得到状态方程，即：

X(k+1)＝Φ_k+1，kX(k)+W(k)k＝1，2，3...(3)

式中，Φ_k+1，k为k至k+1时刻的一步状态转移矩阵；X(k)为第k时刻的状态；W(k)为系统噪声矩阵，对于白噪声有：

E[W(k)]＝0，E[W(k)W(j)]＝Q_kδ_kj j＝1，2，3...(4)

式中，Q_k为系统噪声的方差矩阵；δ：表示δ函数；

步骤二：建立以SAR到地面标识点之间的距离和SAR与地面标识点之间的多普勒频移为观测量的观测方程；

当SAR的成像视场中出现地面标识点时，假设L_i为第i个地面标识点的经度，λ_i为第i个地面标识点的纬度；当i大于等于3时，就可以确定卫星的轨道；所述第i个地面标识点在地固坐标系下的坐标为：

$\left\{\begin{array}{c}{x}_{\mathrm{di}}=R_e\cos {\mathrm{\λ}}_i\cos L_i\\ y_{\mathrm{di}}=R_e\cos {\mathrm{\λ}}_i\sin L_i\\ z_{\mathrm{di}}=R_e\sin {\mathrm{\λ}}_i\end{array}\right.---\left(5\right)$

式中，R_e为地球半径；根据地固坐标系和惯性坐标系的转换关系，得到第i个地面标识点在惯性坐标系下的坐标，即：

$\left[\begin{array}{c}{x}_{\mathrm{mi}}\\ y_{\mathrm{mi}}\\ z_{\mathrm{mi}}\end{array}\right]=\left[\begin{array}{ccc}\cos \mathrm{\Ω}& -\sin \mathrm{\Ω}& 0\\ \sin \mathrm{\Ω}& \cos \mathrm{\Ω}& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_{\mathrm{di}}\\ y_{\mathrm{di}}\\ z_{\mathrm{di}}\end{array}\right]---\left(6\right)$

第i个地面标识点在惯性坐标系下的速度，即：

$\left[\begin{array}{c}{v}_{\mathrm{mxi}}\\ v_{\mathrm{myi}}\\ v_{\mathrm{mzi}}\end{array}\right]=\left[\begin{array}{ccc}0& -{\mathrm{\ω}}_e& 0\\ {\mathrm{\ω}}_e& 0& 0\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}{x}_{\mathrm{mi}}\\ y_{\mathrm{mi}}\\ z_{\mathrm{mi}}\end{array}\right]---\left(7\right)$

式中：ω_e为地球自转角速度；Ω＝Ω₀+ω_e(t-t_o)为t时刻的恒星时；Ω_o为t₀时刻的恒星时；

根据上述式(6)和(7)，得到卫星到第i个地面标识点的距离方程，即：

${\mathrm{\ρ}}_i=\sqrt{{(x-x_{\mathrm{mi}})}^2+{(y-y_{\mathrm{mi}})}^2+{(z-z_{\mathrm{mi}})}^2}+V_1,i=\mathrm{1,2,3}---\left(8\right)$

利用上述式(8)得到多普勒频移方程，即：

$f_i=-\frac{2}{{\mathrm{\λ\ρ}}_i}[(v_x-v_{\mathrm{mxi}})(x-x_{\mathrm{mi}})+(v_y-v_{\mathrm{myi}})(y-y_{\mathrm{mi}})+(v_z-v_{\mathrm{mzi}})(z-z_{\mathrm{mi}})]+V_2---\left(9\right)$

式中，i＝1，2，3

利用上述式(8)和(9)得到观测方程，即：

$G=\left[\begin{array}{c}\mathrm{\ρ}\\ f\end{array}\right]=h[X\left(t\right),t]+V\left(t\right)---\left(10\right)$

式(10)中的ρ＝[ρ₁，ρ₂，ρ₃]^T，f＝[f₁，f₂，f₃]^T，V＝[V₁，V₂]^T

式中：V为测量噪声；G为观测量；

步骤三：建立扩展卡尔曼滤波的递推方程；

设定k时刻的滤波值为

将步骤一中的式(1)中的非线性向量函数F围绕滤波值

展开成泰勒级数，并略去二次以上项，得到线性化方程：

$X(k+1)=f[\hat{X}\left(k\right|k),k]+\frac{\∂F[X\left(k\right|k),k]}{\∂X\left(k\right)}{|}_{X\left(t\right)=\hat{X}\left(k\right|k)}[X\left(k\right)-\hat{X}\left(k\right|k)]+\mathrm{\Γ}(X\left(k\right),k)W\left(k\right)---\left(11\right)$

式中，Γ为系数矩阵；

将步骤二中的式(10)中的非线性向量函数h离散化后，并围绕滤波值

展开成泰勒级数，并略去二次以上项，得到离散化的观测方程，即：

$G(k+1)=h[\hat{X}(k+1|k),k+1]+\frac{\∂h[X(k+1),k+1]}{\∂X(k+1)}{|}_{X\left(t\right)=\hat{X}(k+1|k)}[X(k+1)-\hat{X}(k+1|k)]+V(k+1)---\left(12\right)$

将上述式(11)和(12)代入标准的卡尔曼滤波方程中，得到扩展卡尔曼滤波的递推方程，即：

$\hat{X}(k+1|k+1)=F[\hat{X}\left(k\right|k),k]+K(k+1)\{G(k+1)-h[\hat{X}(k+1|k),k+1\left]\right\}---\left(13\right)$

其中：

K(k+1)＝P(k+1|k)H^T(k+1)[H(k+1)P(k+1|k)H^T(k+1)+R(k+1)]^-1(14)

P(k+1|k)＝Φ(k+1，k)P(k|k)Φ^T(k+1，k)+Q(k)(15)

P(k+1|k+1)＝[I-K(k+1)H(k+1)]P(k+1|k)(16)

滤波初始值为：

$\hat{X}\left(0\right|0)=E\left[X\left(0\right)\right]$

P(0|0)＝P_x(0)

上述式(15)中一步状态转移矩阵

$\mathrm{\Φ}(k+1,k)=\frac{\∂F[X\left(k\right),k]}{\∂X\left(k\right)}{|}_{X\left(k\right)=\hat{X}\left(k\right|k)}$

上述式(16)中的观测矩阵

$H(k+1)=\frac{\∂h\left[X(k+1,k+1)\right]}{\∂X(k+1)}{|}_{X\left(k\right)=\hat{X}(k+1|k)}---\left(17\right)$

式中，K为增益矩阵；P(k|k)为k时刻的协方差矩阵；P(k+1|k)为协方差矩阵的一步预测；Q为系统噪声的协方差矩阵；R为测量噪声的协方差矩阵；I为单位矩阵；

步骤四：得到卫星的状态信息；

利用步骤三建立的扩展卡尔曼滤波的递推方程得到卫星的状态信息，所述卫星的状态信息即为卫星的状态估计

所述卫星的状态估计

包括卫星位置和速度矢量 $\hat{X}={[x,y,z,v_x{,v}_y,v_z]}^T.$

Images