CN104459751A - 基于gnss反射信号的双站雷达空间目标相对导航方法 - Google Patents

Abstract

基于GNSS反射信号的双站雷达空间目标相对导航方法，提出利用GNSS卫星作为照射源，探测在主动航天器的GNSS反射信号作空间目标的相对测量实现途径，利用运动方程和状态估计方法，确定航天器之间的相对位置和速度，得到目标航天器和主动航天器在不同距离时的运动方程模型。利用安装在主动航天器上的探测器接收直达和反射信号测量航天器之间的相对距离，建立基于GNSS反射信号双站雷达空间目标相对导航运动模型；采用卡尔曼滤波方法构造绝对状态估计器和相对状态估计器，从而获取航天器间的相对位置和相对速度信息。本发明方法不仅不需要单独的交会敏感器，并且主动航天器和目标航天器相对距离远近不同时，都能获得完备的观测信息。

Description

基于GNSS反射信号的双站雷达空间目标相对导航方法

技术领域

本发明属于通信领域，涉及一种空间目标的相对导航方法。

背景技术

空间目标相对测量，是指追踪飞行器基于其星载相对导航设备的测量，在一定条件下实时估计目标飞行器的相对状态，包括相对位置、相对速度或其它表征相对运动的信息。相对测量是飞行器轨道交会、伴飞与逼近过程制导和控制的输入，相对测量的性能将直接影响最终的制导和控制精度。目前空间目标常用的测量手段包括：微波交会对接雷达、激光雷达、卫星导航测量和CCD光学成像测量，这些测量手段在空间交会对接过程中得到了广泛应用，但上述传统的相对测量手段大多针对合作目标，即追踪飞行器与目标之间有信息交换。

随着人类探索、开发和利用外层空间的深入，对空间目标相对导航技术提出了更高的要求。诸如卫星由于故障、完全失效或任务结束而被放弃后，停留在空间将成为太空垃圾，不但占用了宝贵的轨道资源，还可能危及其它航天器的安全。为了尽可能挽回损失或净化轨道环境，各国正在研究以卫星维修、生命延长及太空垃圾清除为目的的在轨服务技术，即对故障卫星进行绕飞监测、跟踪接近、交会对接、在轨维修等操作。要实现这些目的，必须解决航天器间相对位姿的测量问题，与传统目标高精度相对导航问题不同，对于大多数已在轨服务的航天器和空间碎片等目标，其本身没有配备目标标识器和交会对接敏感器，并没有预先安装相对测量装置，空间目标之间没有直接的信息交流，这就使得空间目标相对导航问题能够得到的观测信息通常是不完备和不精确的，因而对相对导航信息的高精度确定造成困难。

针对空间目标的相对导航时空间交会中的难点问题，已经引起了相关研究部门的广泛关注，目前已经提出了一些相对导航位姿的测量方法，其中以基于立体视觉的航天器相对导航技术最为典型，其在应用时，需要在主动航天器上安装交会敏感器。安装在主动航天器上的交会敏感器在航天器相距较远时，很难获得完备的观测信息，导致自主导航算法不能满足相对导航需求，而只能依靠地面引导主动航天器进行交会对接，从而增加了地面测控的负担。

发明内容

本发明解决的技术问题是：克服现有技术的不足，提供了一种基于GNSS反射信号的双站雷达空间目标相对导航方法，该方法利用GNSS反射信号测量原理，建立了目标航天器和主动航天器的运动模型，在此基础上，提出了基于GNSS反射信号双站雷达空间目标相对导航模型，不仅不需要单独的交会敏感器，并且主动航天器和目标航天器相对距离远近不同时，都能获得完备的观测信息，自主导航算法能够满足相对导航需求，减轻了地面测控的负担，还可以实现全天候、全天时工作。

本发明的技术解决方案是：基于GNSS反射信号的双站雷达空间目标相对导航方法，包括如下步骤：

(1)在主动航天器上安装GNSS卫星导航接收机，利用卫星导航接收机同时接收GNSS卫星的直达信号，以及目标航天器反射的GNSS卫星的回波信号；

(2)主动航天器获取自身的位置和速度，作为惯性坐标系下主动航天器的绝对运动状态变量，同时估计目标航天器的位置和速度，作为惯性坐标系下目标航天器的绝对运动状态变量；将主动航天器的绝对运动状态变量和目标航天器的绝对运动状态变量进行差分，得到惯性坐标系下主动航天器和目标航天器之间的相对位置和相对速度，作为惯性坐标系下的相对运动状态变量；

(3)判断主动航天器和目标航天器之间的距离，当主动航天器和目标航天器之间距离超过d时，转步骤(4)；当主动航天器和目标航天器之间的距离不超过d时，转步骤(6)；所述d的大小取决于主动航天器的轨道高度s，满足d≤0.01s；

(4)将步骤(2)中获取的惯性坐标系下的全部绝对运动状态变量，以及主动航天器与GNSS卫星之间的钟差b作为绝对状态向量x，采用迭代的方式，根据k-1时刻绝对状态向量的估计值迭代更新k时刻绝对状态向量的估计值具体为：

(4-1)计算得到k-1时刻绝对状态向量中各绝对状态变量的导数值，其中：

$\left\{\begin{array}{c}{\stackrel{\·}{r}}_c=v_c\\ {\stackrel{\·}{v}}_c=-\frac{\mathrm{\μ}}{r_c^3}\·r_c+u_c\end{array}\right.$

$\left\{\begin{array}{c}{\stackrel{\·}{r}}_T=v_T\\ {\stackrel{\·}{v}}_T=-\frac{\mathrm{\μ}}{r_T^3}\·r_T+u_T\end{array}\right.$

$\stackrel{\·}{b}=0$

r_c和v_c分别为惯性坐标系下主动航天器的绝对位置和绝对速度，r_T和v_T分别为惯性坐标系下目标航天器的绝对位置和绝对速度，u_c和u_T分别为惯性坐标系下主动航天器和目标航天器所受外部合力，μ为地心引力常数；

(4-2)对步骤(4-1)计算得到的导数值，从k-1时刻到k时刻的时间区间进行积分，得到各绝对状态变量的变化量，在此基础上叠加k-1时刻相应绝对状态变量的估计值，作为k时刻绝对状态向量的预测值

(4-3)利用k时刻绝对状态向量的预测值预测得到k时刻所述直达信号和回波信号的预测值向量 $h_{\mathrm{abs}}\left({\hat{x}}_{k,k-1}\right)={[h_{\mathrm{\ρ}}\left({\hat{x}}_{k,k-1}\right),h_{\stackrel{\·}{\mathrm{\ρ}}}\left({\hat{x}}_{k,k-1}\right),h_{\mathrm{r\ρ}}\left({\hat{x}}_{k,k-1}\right),h_{r\stackrel{\·}{\mathrm{\ρ}}}\left({\hat{x}}_{k,k-1}\right)]}^T,$ 其中直达信号的预测值为：

$h_{\mathrm{\ρ}}\left({\hat{x}}_{k,k-1}\right)=|r_c-r_g|+b$

$h_{\stackrel{\·}{\mathrm{\ρ}}}\left({\hat{x}}_{k,k-1}\right)={(v_c-v_g)}^T\frac{r_c-r_g}{|r_c-r_g|}$

回波信号的预测值为：

$h_{\mathrm{r\ρ}}\left({\hat{x}}_{k,k-1}\right)=|r_c-r_T|+|r_T-r_g|-|r_g-r_c|$

$h_{r\stackrel{\·}{\mathrm{\ρ}}}\left({\hat{x}}_{k,k-1}\right)={(v_c-v_T)}^T\frac{r_c-r_T}{|r_c-r_T|}+{(v_T-v_g)}^T\frac{r_T-r_g}{|r_T-r_g|}$

r_g和v_g分别表示GNSS卫星在惯性坐标系下的绝对位置和绝对速度；

(4-4)将步骤(4-3)获取的直达信号和回波信号作为预测值，将步骤(1)中获取的直达信号和回波信号作为测量值，采用卡尔曼滤波的方法，计算得到绝对状态向量的修正值

(4-5)通过公式计算得到k时刻绝对状态向量的估计值

(5)从步骤(4)得到的绝对状态向量估计值中提取出对应的分量，分别得到惯性坐标系下主动航天器和目标航天器各自的绝对位置和绝对速度，再通过直接差分，得到惯性坐标系下主动航天器和目标航天器之间的相对位置和相对速度，作为相对导航的结果，并结束；

(6)将步骤(2)获取的惯性坐标系下的相对运动状态变量转换至轨道坐标系，将轨道坐标系下的全部相对运动状态变量作为相对状态向量X，采用迭代的方式，采用迭代的方式，根据k-1时刻相对状态向量的估计值迭代更新k时刻相对状态向量的估计值具体为：

(6-1)计算得到k-1时刻相对状态向量中各相对状态变量的导数值，其中：

$\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\Δ}}}_r={\mathrm{\Δ}}_v\\ {\stackrel{\·}{\mathrm{\Δ}}}_v=A\·{\mathrm{\Δ}}_r+B\·{\mathrm{\Δ}}_v+u_{\mathrm{\Δ}}\end{array}\right.$

Δ_r和Δ_v分别为轨道坐标系下主动航天器相对于目标航天器的相对位置和相对速度，u_Δ为轨道坐标系下主动航天器相对于目标航天器的相对外部合力，

$A=\left[\begin{array}{ccc}3{\mathrm{\ω}}^2& 0& 0\\ 0& 0& 0\\ 0& 0& -{\mathrm{\ω}}^2\end{array}\right]$

$B=\left[\begin{array}{ccc}0& 2\mathrm{\ω}& 0\\ -2\mathrm{\ω}& 0& 0\\ 0& 0& 0\end{array}\right]$

ω为主动航天器的平均轨道角速度；

(6-2)对步骤(6-1)计算得到的导数值，从k-1时刻到k时刻的时间区间进行积分，得到各相对状态变量的变化量，在此基础上叠加k-1时刻相应相对状态变量的估计值，作为k时刻相对状态向量的预测值

(6-3)利用k时刻相对状态向量的预测值预测得到k时刻所述直达信号和回波信号差值的预测值向量其中

$h_r\left({\hat{X}}_{k,k-1}\right)=\left|{\mathrm{\Δ}}_r\right|+|{\mathrm{\Δ}}_r-X_g|-\left|X_g\right|$

$h_v\left({\hat{X}}_{k,k-1}\right)={{\mathrm{\Δ}}_v}^T\frac{{\mathrm{\Δ}}_r}{\left|{\mathrm{\Δ}}_r\right|}+{({\mathrm{\Δ}}_v-V_g)}^T\frac{{\mathrm{\Δ}}_r-X_g}{|{\mathrm{\Δ}}_r-X_g|}$

X_g和V_g分别为GNSS卫星在轨道坐标系下的位置和速度；

(6-4)将步骤(6-3)获取的直达信号和回波信号的差值作为预测值，将步骤(1)中获取的直达信号和回波信号做差后的差值作为测量值，采用卡尔曼滤波的方法，计算得到相对状态向量的修正值

(6-5)通过公式计算得到k时刻相对状态向量的估计值

(7)从步骤(6)得到的相对状态向量估计值中提取出对应的分量，得到轨道坐标系下主动航天器和目标航天器的相对位置和相对速度，作为相对导航的结果并结束。

所述步骤(1)中利用卫星导航接收机的左旋天线接收GNSS卫星的直达信号，利用卫星导航接收机的右旋天线接收目标航天器反射的GNSS卫星的回波信号。所述的d的取值范围为1～100千米。

本发明与现有技术相比的优点在于：

(1)本发明方法兼顾了所有的GNSS反射信号的接收特性，以布鲁斯特角的大小划定采用左旋天线接收机接收GNSS反射信号，精确了反射信号的收集范围，为方法的正确实施奠定了基础；

(2)本发明方法中，在两航天器相距较远时，将两航天器相对运动采用绝对运动方程描述，在两航天器相距较近且运行在近圆轨道的条件下，将两航天器相对运动用Clohessy-Wiltshire方程进行近似描述，在同一方法中涵盖了相对导航距离的远近，方法具有普遍性和适用性；

(3)本发明方法中，基于GNSS反射信号双站雷达空间目标相对运动模型是相对导航方法的一个关键方程，针对绝对运动方程条件与C-W运动方程运动条件，分别在惯性参考系与轨道参考系中进行建立，弥补了空间目标相对导航模型在单一坐标系下参数难以转换的状况；

(4)全球卫星导航定位系统的反射信号相对导航技术作为一种新型的、低成本的、高机动性的导航技术，与其他相对导航优势互补，可以增加相对手段的多样性，弥补空间目标相对导航信息交流不足的状况，具有技术前沿性；并且本发明方法不需要单独的交会敏感器即能获得完备的观测信息，自主导航算法能够满足相对导航需求，减轻了地面测控的负担，还可以实现全天候、全天时工作。

附图说明

图1为本发明方法的流程框图；

图2为本发明GNSS反射信号探测器示意图；

图3为本发明GNSS信号传播时间示意图；

图4为本发明方法几何关系示意图。

具体实施方式

首先给出两个坐标系的定义。

(地心)惯性坐标系：坐标点位于地球质心，x轴指向春分点，z轴指向地球北极，y轴完成右手正交坐标系。其中，常用的地心惯性坐标系是J2000坐标系，它的坐标原点在地球质心，参考平面是J2000.0平赤道面，z轴向北指向平赤道面北极，x轴指向J2000.0平春分点，y轴与x轴和z轴组成直角右手系；

轨道坐标系：也即当地垂直水平坐标系，坐标原点固连于航天器所处轨道位置，x轴沿地心指向航天器质心的矢量方向，z轴指向航天器轨道角动量正方向，y轴完成右手正交坐标系。

下面将结合附图对本发明做进一步的详细说明。

本发明是一种基于GNSS反射信号双站雷达空间目标相对导航方法，其方法流程框图如图1所示，具体实施方式如下：

步骤1：主动航天器利用卫星导航接收机同时接收导航卫星的直达信号和目标航天器的回波反射信号，其中直达信号利用卫星导航接收机的左旋天线接收，目标航天器回波反射信号利用卫星导航接收机的右旋天线接收；

全球卫星导航系统(GNSS)信号是圆极化平面波，发射信号具有右旋圆极化(Right Hand Circular Polarization，RHCP)极性，其信号可以分解为水平和垂直的线性极化波。水平和垂直极化波的反射系数分别为：

${\mathrm{\Γ}}_H=\frac{\sin \mathrm{\θ}-\sqrt{\mathrm{\ϵ}-{\cos }^2\mathrm{\θ}}}{\sin \mathrm{\θ}+\sqrt{\mathrm{\ϵ}-{\cos }^2\mathrm{\θ}}}---\left(1\right)$

${\mathrm{\Γ}}_V=\frac{\mathrm{\ϵ}\sin \mathrm{\θ}-\sqrt{\mathrm{\ϵ}-{\cos }^2\mathrm{\θ}}}{\mathrm{\ϵ}\sin \mathrm{\θ}+\sqrt{\mathrm{\ϵ}-{\cos }^2\mathrm{\θ}}}---\left(2\right)$

式中，ε＝ε_r-j60λσ，ε_r为相对介电常数，σ为介质的导电率，λ为电磁波的波长，θ为入射角，即信号入射方向与入射点切面的夹角。

当GNSS信号经目标航天器反射后，其反射信号可以认为是由两部分组成的，一部分保持原来的极化方向(右旋)不变。一部分则与原来的极化方向相反，即左旋圆极化(Left Hand Circular Polarization，LHCP)，这两部分的反射系数可以分别表示为：

${\mathrm{\Γ}}_R=\frac{{\mathrm{\Γ}}_H+{\mathrm{\Γ}}_V}{2}---\left(3\right)$

${\mathrm{\Γ}}_R=\frac{{\mathrm{\Γ}}_H-{\mathrm{\Γ}}_V}{2}---\left(4\right)$

在GNSS信号在目标航天器发生反射时，当卫星仰角很小时，即入射角很大时，其反射波中右旋极化波占主要地位；当卫星仰角较大大时，即入射角较大时，其反射波中左旋极化波占主要地位；这种现象发生的临界角度为布鲁斯特角(brewster)。当入射角小于布鲁斯特角时，反射信号中RHCP波为主要成分；当入射角大于布鲁斯特角时，反射信号中LHCP波占主要成分。根据目标航天器的介质成份，一般可见，GNSS卫星高度角大于10°，因此采用LHCP天线可以接收到较强的目标航天器反射GNSS信号。

因此，安装在主动航天器上的GNSS探测器主要由一个指向目标航天器的右旋(RHCP)天线接收机和一个指向GNSS卫星的左旋(LHCP)天线接收机组成，如图2所示，分别接收GNSS卫星的直射信号和反射回波信号，主动航天器利用右旋天线接收GNSS卫星直射信号，用于导航定位解算，确定接收机的空间位置、计算反射信号的多普勒频移和估算时间延迟，反射信号由左旋天线接收，根据直射信号和反射信号的特性，输出相对于直射信号不同时间延迟时刻的反射信号的相关信息，解算相关信息。

步骤2：确定目标航天器和主动航天器在近距离和远距离时运动方程模型，利用运动方程和状态估计方法，确定航天器之间的相对位置和相对速度；

航天器的在轨运动规律可用下式近似描述

$\left\{\begin{array}{c}\stackrel{\·}{r}=v\\ \stackrel{\·}{v}=-\frac{{\mathrm{Gm}}_e}{r^3}\·r+u\end{array}\right.---\left(5\right)$

其中，r为航天器的位置矢量，v为航天器的速度矢量，G为万有引力常数，m_e为地球质量，u为航天器所受外部合力。记为地心引力常数，对位置矢量r、速度矢量v、外部合力u加以下标‘c’和‘T’，分别代表主动航天器和目标航天器的相关物理量，则两航天器各自的绝对运动方程为：

$\left\{\begin{array}{c}{\stackrel{\·}{r}}_c=v_c\\ {\stackrel{\·}{v}}_c=-\frac{\mathrm{\μ}}{r_c^3}\·r_c+u_c\end{array}\right.---\left(6\right)$

$\left\{\begin{array}{c}{\stackrel{\·}{r}}_T=v_T\\ {\stackrel{\·}{v}}_T=-\frac{\mathrm{\μ}}{r_T^3}\·r_T+u_T\end{array}\right.---\left(\text{7}\right)$

此外，主动航天器与GNSS卫星的钟差(以距离表示)记为b，因其变化缓慢，预测得到：

$\stackrel{\·}{b}=0---\left(8\right)$

构造绝对状态向量为x＝[r_c ^T,v_c ^T,r_T ^T,v_T ^T,b]^T，以函数f_abs等价表示式(6)～(8)，即有：

$\stackrel{\·}{x}=f_{\mathrm{abs}}(x,u)---\left(9\right)$

对上式离散化得到离散形式的运动方程为

x_k+1＝F_abs(x_k,u)                (10)

在两航天器相距较近，且运行在近圆轨道的条件下，两航天器相对运动可以用Clohessy-Wiltshire方程进行近似描述：

$\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\Δ}}}_r={\mathrm{\Δ}}_v\\ {\stackrel{\·}{\mathrm{\Δ}}}_v=A\·{\mathrm{\Δ}}_r+B\·{\mathrm{\Δ}}_v+u_{\mathrm{\Δ}}\end{array}\right.---\left(11\right)$

其中Δ_r、Δ_v、u_Δ是主动航天器相对于目标航天器的相对位置矢量、相对速度矢量、相对外部合力在目标航天器轨道坐标系下的表示，矩阵A和矩阵B由下式确定

$A=\left[\begin{array}{ccc}3{\mathrm{\ω}}^2& 0& 0\\ 0& 0& 0\\ 0& 0& -{\mathrm{\ω}}^2\end{array}\right]---\left(12\right)$

$B=\left[\begin{array}{ccc}0& 2\mathrm{\ω}& 0\\ -2\mathrm{\ω}& 0& 0\\ 0& 0& 0\end{array}\right]---\left(13\right)$

其中，ω为主动航天器圆轨道的平均轨道角速度。

构造相对状态向量为X＝[Δ_r ^T,Δ_v ^T]^T，以函数f_rel等价表示式(11)，即有

$\stackrel{\·}{X}=f_{\mathrm{rel}}(X,u_{\mathrm{\Δ}})---\left(14\right)$

对上式离散化得到离散形式的运动方程为

x_k+1＝F_abs(x_k,u_Δ)                      (15)

对式(10)和式(15)分别考虑噪声扰动向量w_abs和w_rel，将用于状态向量的传播预测。传播预测同时伴随有系统误差，其协方差矩阵分别表示为Q_abs和Q_rel，即

$Q_{\mathrm{abs}}\stackrel{\mathrm{\Δ}}{=}E(w_{\mathrm{abs}}\·{w_{\mathrm{abs}}}^T)---\left(16\right)$

$Q_{\mathrm{rel}}\stackrel{\mathrm{\Δ}}{=}E(w_{\mathrm{rel}}\·{w_{\mathrm{rel}}}^T)---\left(17\right)$

其中，E(·)表示数学期望，上标‘T’表示向量或矩阵的转置。当噪声向量的分量互相独立且服从高斯分布时，构造Q_abs、Q_rel的方法为

$Q_{\mathrm{abs}}=\left[\begin{array}{ccccc}{{\mathrm{\σ}}_{\mathrm{\ρc}}}^2{I}_{3\×3}& 0& 0& 0& 0\\ 0& {{\mathrm{\σ}}_{\stackrel{\·}{\mathrm{\ρ}}c}}^2{I}_{3\×3}& 0& 0& 0\\ 0& 0& {{\mathrm{\σ}}_{\mathrm{\ρT}}}^2{I}_{3\×3}& 0& 0\\ 0& 0& 0& {{\mathrm{\σ}}_{\stackrel{\·}{\mathrm{\ρ}}T}}^2{I}_{3\×3}& 0\\ 0& 0& 0& 0& {{\mathrm{\σ}}_b}^2\end{array}\right]---\left(18\right)$

$Q_{\mathrm{rel}}=\left[\begin{array}{cc}{{\mathrm{\σ}}_r}^2{I}_{3\×3}& 0\\ 0& {{\mathrm{\σ}}_v}^2{I}_{3\×3}\end{array}\right]---\left(19\right)$

其中，σ_ρc、分别为主动航天器绝对运动模型中位置、速度变化率的随机扰动的标准差；σ_ρT、分别为目标航天器绝对运动模型中位置、速度变化率的随机扰动的标准差；σ_b为主动航天器钟差变化率的标准差，σ_r、σ_v分别为相对运动模型中位置、速度变化率的随机扰动的标准差，I_3×3表示3阶单位矩阵。

步骤3：利用安装在目标航天器上的卫星导航接收机接收的直达信号和反射信号，测量目标航天器和主动航天器之间的相对距离，建立基于GNSS反射信号双站雷达空间目标相对导航运动模型。

GNSS系统中认为接收机与卫星间存在视距传输(Line of Sight，LoS)，接收机估计无线电波从发射到接收的传播延时，得到接收机到卫星的直线距离，从而得到测距的目的。

GNSS信号采用伪随机噪声(Pseudo Random Noise，PRN)序列作为扩频序列。这类序列具有类似随机二进制序列的频谱，但长度及内容固定。PRN序列的发送是周期性的，如GPS C/A码序列周期为1ms，P码周期为7天。假设接收机时钟理想同步，则卫星发射信号可通过测量当前接收信号与本地复现伪码间的时间差得到信号由卫星至接收机天线的传输时间。假设在t₂时刻卫星发射的信号在t₃时刻到达接收机，记光速为c，则卫星到接收机的距离可表示为r＝c(t₃-t₂)。受限于接收机时钟精度，接收机复现伪码与卫星实际发射信号并不同步。因此，在计算当前接收信号传输时间时，利用本地复现伪码表示的信号发射时间与实际信号发射时间存在误差。此时测量得到的距离ρ称作伪距，忽略电离层误差、GPS星历误差、接收机硬件误差以及多路径效应等其他误差，如图3所示，其定义简化如下：

ρ＝c(t₃-t₂)＝c(t₃-t₁)-c(t₂-t₁)＝|r_gc|+cδt

其中，δt表示GNSS系统时与接收机本地时间的偏差，根据步骤2的状态向量定义方法，采用函数的观点对上式进行等价描述，即当系统具有状态x时，主动航天器直接伪距由下式决定：

h_ρ(x)＝|r_c-r_g|+b

其中，r表示绝对位置矢量，下标‘c’和‘g’分别表示主动航天器和GNSS卫星，|·|表示向量的模长。其中GNSS卫星的绝对位置矢量通过卫星星历以及参考系转换的方式进行确定。

根据运动学关系，主动航天器的直接伪距率函数由下式确定

$h_{\stackrel{\·}{\mathrm{\ρ}}}\left(x\right)={(v_c-v_g)}^T\frac{r_c-r_g}{|r_c-r_g|}$

当目标航天器收到GNSS卫星信号后反射至主动航天器时，结合主动航天器的直接伪距，可以形成反射伪距差动测量。根据几何关系，如图4所示，当系统具有状态x时，反射伪距差动h_rρ(x)由下式决定：

h_rρ(x)＝|r_c-r_T|+|r_T-r_g|-|r_g-r_c|            (EQ-01)

其中下标‘T’表示目标航天器，其他符号定义如前。

根据运动学关系，反射伪距率差动函数由下式确定：

$h_{r\stackrel{\·}{\mathrm{\ρ}}}\left(x\right)={(v_c-v_T)}^T\frac{r_c-r_T}{|r_c-r_T|}+{(v_T-v_g)}^T\frac{r_T-r_g}{|r_T-r_g|}---(\mathrm{EQ}-02)$

根据上述函数，绝对状态观测的向量值函数h_abs(x)由下式进行构造

$h_{\mathrm{abs}}\left(x\right)={[h_{\mathrm{\ρ}}\left(x\right),h_{\stackrel{\·}{\mathrm{\ρ}}}\left(x\right),h_{\mathrm{r\ρ}}\left(x\right),h_{r\stackrel{\·}{\mathrm{\ρ}}}\left(x\right)]}^T---(\mathrm{EQ}-03)$

至此，在地心惯性系中考虑的绝对状态观测函数得到确定。

同理，需要建立在轨道参考系中考虑的相对状态观测函数。在轨道参考系中，式(EQ-01)可以表示为

h_r(X)＝|Δ_r|+|Δ_r-X_g|-|X_g|

其中h_r(X)是当系统相对状态为X时对应的反射伪距差动函数值，X_g为GNSS卫星在轨道参考系中的位置矢量，通过卫星星历和参考系转换的方式进行确定，其他符号定义如前。同样地，式(EQ-02)转换到轨道参考系中表示为

$h_v\left(X\right)={{\mathrm{\Δ}}_v}^T\frac{{\mathrm{\Δ}}_r}{\left|{\mathrm{\Δ}}_r\right|}+{({\mathrm{\Δ}}_v-V_g)}^T\frac{{\mathrm{\Δ}}_r-X_g}{|{\mathrm{\Δ}}_r-X_g|}$

其中，h_v(X)是当系统相对状态为X时对应的反射伪距率差动函数值，V_g为GNSS卫星在轨道参考系中的速度矢量，通过卫星星历和参考系转换的方式进行确定，其他符号定义如前。

根据上述函数，相对状态观测的向量值函数h_rel(X)由下式进行构造

h_rel(X)＝[h_r(X),h_v(X)]^T             (EQ-04)

式(EQ-03)和式(EQ-04)可以用于根据当前的状态向量预测测量结果。然而由于测量真值与测量函数之间存在的测量噪声η_abs和η_rel，对应于绝对状态观测关系和相对状态观测关系，这种预测将伴随有误差协方差矩阵R_abs和R_rel，即

$R_{\mathrm{abs}}\stackrel{\mathrm{\Δ}}{=}E({\mathrm{\η}}_{\mathrm{abs}}\·{{\mathrm{\η}}_{\mathrm{abs}}}^T)$

$R_{\mathrm{rel}}\stackrel{\mathrm{\Δ}}{=}E({\mathrm{\η}}_{\mathrm{rel}}\·{{\mathrm{\η}}_{\mathrm{rel}}}^T)$

当噪声向量的分量互相独立且服从高斯分布时，可以具体构造R_abs、R_rel为

$R_{\mathrm{abs}}=\left[\begin{array}{cccc}{{\mathrm{\σ}}_p}^2{I}_{m\×m}& 0& 0& 0\\ 0& {{\mathrm{\σ}}_{\stackrel{\·}{p}}}^2{I}_{m\×m}& 0& 0\\ 0& 0& {{\mathrm{\σ}}_{\mathrm{Dp}}}^2{I}_{m\×m}& 0\\ 0& 0& 0& {{\mathrm{\σ}}_{D\stackrel{\·}{p}}}^2{I}_{m\×m}\end{array}\right]$

$R_{\mathrm{rel}}=\left[\begin{array}{cc}{{\mathrm{\σ}}_{\mathrm{Dp}}}^2{I}_{m\×m}& 0\\ 0& {{\mathrm{\σ}}_{D\stackrel{\·}{p}}}^2{I}_{m\×m}\end{array}\right]$

式中，σ_Dp、分别为反射伪距差分测量信号误差的标准差；σ_p、分别为直接测量信号误差的标准差，I_m×m表示m阶单位矩阵，m是可用GNSS卫星的数量。

步骤4：采用卡尔曼滤波方法，对步骤1-3中所述的应用场合、系统组成和测量关系，构造绝对状态估计器和相对状态估计器，从而获取两航天器间的相对位置和相对速度信息。

步骤4-1

建立估计器的数据处理机制，并初始化估计器。具体流程为，将绝对状态向量和相对状态向量各自以时间标签进行标注，用下标符号‘k’表示第k个时间标签的数据，用下标符号‘k,k-1’表示第k-1与第k个时间标签之间的时间片段。初始时间标签为下标‘0’。

初始化绝对状态向量相对状态向量以及各自的误差协方差矩阵P₀(实际上分别为P_abs-0和P_rel-0)。其方法为：主动航天器采用传统的GNSS独立导航求解自身位置与速度，并根据轨道信息大致估计出目标航天器的位置与速度，将二者进行参考系转换表示到地心惯性参考系，结合钟差偏移初值确定同时将得到的两航天器位置和速度进行差分，进行参考系转换表示到目标航天器的轨道参考系，确定根据初值相对于真值的偏移程度，初始化P₀为较大的值。

步骤4-2

从时间标签0开始，递归计算估计状态，由以下①～⑤进行。

①状态预测。在时间标签为k-1的时刻，根据当前的利用式(10)和式(15)进行状态预测，得到

②测量预测。在时间标签为k-1的时刻，根据预测的状态利用式(EQ-03)和式(EQ-04)进一步预测测量值，得到

③准备更新。在时间标签为k-1的时刻，计算用于递归运算的状态转移矩阵Φ_k,k-1，观测变化率矩阵H_k，协方差预测矩阵P_k,k-1，增益矩阵K_k

${\mathrm{\Φ}}_{k,k-1}=e^{\mathrm{\Δt}\·\frac{\∂f}{\∂x}{|}_{{\hat{x}}_k}}$

$H_k=\frac{\∂h}{\∂x}{|}_{{\hat{x}}_k}$

P_k,k-1＝Φ_k,k-1P_k-1Φ_k,k-1 ^T+Q_k-1

K_k＝P_k,k-1H_k ^T(H_kP_k,k-1H_k ^T+R_k)^-1

其中，Δt为时间间隔，表示关于向量x的函数f在点处的梯度值。对于绝对状态估计器，其观测变化率矩阵为：

$H_k=\left[\begin{array}{ccccc}{\left(\frac{{\∂h}_{\mathrm{\ρ}}}{{\∂r}_c}\right)}^T& {\left(\frac{{\∂h}_{\mathrm{\ρ}}}{{\∂v}_c}\right)}^T& {\left(\frac{{\∂h}_{\mathrm{\ρ}}}{{\∂r}_T}\right)}^T& {\left(\frac{{\∂h}_{\mathrm{\ρ}}}{{\∂v}_T}\right)}^T& {\left(\frac{{\∂h}_{\mathrm{\ρ}}}{\∂b}\right)}^T\\ {\left(\frac{{\∂h}_{\stackrel{\·}{\mathrm{\ρ}}}}{{\∂r}_c}\right)}^T& {\left(\frac{{\∂h}_{\stackrel{\·}{\mathrm{\ρ}}}}{{\∂v}_c}\right)}^T& {\left(\frac{{\∂h}_{\stackrel{\·}{\mathrm{\ρ}}}}{{\∂r}_T}\right)}^T& {\left(\frac{{\∂h}_{\stackrel{\·}{\mathrm{\ρ}}}}{{\∂v}_T}\right)}^T& {\left(\frac{{\∂h}_{\stackrel{\·}{\mathrm{\ρ}}}}{\∂b}\right)}^T\\ {\left(\frac{{\∂h}_{\mathrm{r\ρ}}}{{\∂r}_c}\right)}^T& {\left(\frac{{\∂h}_{\mathrm{r\ρ}}}{{\∂v}_c}\right)}^T& {\left(\frac{{\∂h}_{\mathrm{r\ρ}}}{{\∂r}_T}\right)}^T& {\left(\frac{{\∂h}_{\mathrm{r\ρ}}}{{\∂v}_T}\right)}^T& {\left(\frac{{\∂h}_{\mathrm{r\ρ}}}{\∂b}\right)}^T\\ {\left(\frac{{\∂h}_{r\stackrel{\·}{\mathrm{\ρ}}}}{{\∂r}_c}\right)}^T& {\left(\frac{{\∂h}_{r\stackrel{\·}{\mathrm{\ρ}}}}{{\∂v}_c}\right)}^T& {\left(\frac{{\∂h}_{r\stackrel{\·}{\mathrm{\ρ}}}}{{\∂v}_T}\right)}^T& {\left(\frac{{\∂h}_{r\stackrel{\·}{\mathrm{\ρ}}}}{{\∂v}_T}\right)}^T& {\left(\frac{{\∂h}_{r\stackrel{\·}{\mathrm{\ρ}}}}{\∂b}\right)}^T\end{array}\right]{|}_{{\hat{x}}_k}$

其中：

$\frac{{\∂h}_{\mathrm{\ρ}}}{{\∂r}_c}=\frac{r_c-r_g}{|r_c-r_g|},$

$\frac{{\∂h}_{\mathrm{\ρ}}}{{\∂v}_c}=0_{3\×1},\frac{{\∂h}_{\mathrm{\ρ}}}{{\∂r}_T}=0_{3\×1},\frac{{\∂h}_{\mathrm{\ρ}}}{{\∂v}_T}=0_{3\×1},\frac{{\∂h}_{\mathrm{\ρ}}}{\∂b}=1$

$\frac{{\∂h}_{\stackrel{\·}{\mathrm{\ρ}}}}{{\∂r}_c}=(\frac{1}{|r_c-r_g|}-\frac{{(r_c-r_g)(r_c-r_g)}^T}{{|r_c-r_g|}^3})(v_c-v_g),\frac{{\∂h}_{\stackrel{\·}{\mathrm{\ρ}}}}{{\∂v}_c}=\frac{r_c-r_g}{|r_c-r_g|}$

$\frac{{\∂h}_{\stackrel{\·}{\mathrm{\ρ}}}}{{\∂r}_T}=0_{3\×1},\frac{{\∂h}_{\stackrel{\·}{\mathrm{\ρ}}}}{{\∂v}_T}=0_{3\×1},\frac{{\∂h}_{\stackrel{\·}{\mathrm{\ρ}}}}{\∂b}=0$

$\frac{{\∂h}_{\mathrm{r\ρ}}}{{\∂r}_c}=\frac{r_c-r_T}{|r_c-r_T|}-\frac{r_c-r_g}{|r_c-r_g|},\frac{{\∂h}_{\mathrm{r\ρ}}}{{\∂v}_c}=0_{3\×1}$

$\frac{{\∂h}_{\mathrm{r\ρ}}}{{\∂r}_T}=\frac{r_T-r_c}{|r_T-r_c|}+\frac{r_T-r_g}{|r_T-r_g|},\frac{{\∂h}_{\mathrm{r\ρ}}}{{\∂v}_T}=0_{3\×1},\frac{{\∂h}_{\mathrm{r\ρ}}}{\∂b}=0$

$\frac{{\∂h}_{\stackrel{\·}{\mathrm{\ρ}}}}{{\∂r}_c}=(\frac{1}{|r_c-r_g|}-\frac{(r_c-r_g){(r_c-r_g)}^T}{{|r_c-r_g|}^3})(v_c-v_g),\frac{{\∂h}_{\stackrel{\·}{\mathrm{\ρ}}}}{{\∂v}_c}=\frac{r_c-r_g}{|r_c-r_g|}$

$\frac{{\∂h}_{\stackrel{\·}{\mathrm{\ρ}}}}{{\∂r}_T}=0_{3\×1},\frac{{\∂h}_{\stackrel{\·}{\mathrm{\ρ}}}}{{\∂v}_T}=0_{3\×1},\frac{{\∂h}_{\stackrel{\·}{\mathrm{\ρ}}}}{\∂b}=0$

对于相对状态估计器，其观测变化率矩阵为：

$H_k=\left[\begin{array}{cc}{\left(\frac{{\∂h}_r}{{\∂\mathrm{\Δ}}_r}\right)}^T& {\left(\frac{{\∂h}_r}{{\∂\mathrm{\Δ}}_v}\right)}^T\\ {\left(\frac{{\∂h}_v}{{\∂\mathrm{\Δ}}_v}\right)}^T& {\left(\frac{{\∂h}_v}{{\∂\mathrm{\Δ}}_v}\right)}^T\end{array}\right]{|}_{{\hat{X}}_k}$

其中：

$\frac{{\∂h}_r}{{\∂\mathrm{\Δ}}_r}=\frac{{\mathrm{\Δ}}_r}{\left|{\mathrm{\Δ}}_r\right|}+\frac{{\mathrm{\Δ}}_r-X_g}{|{\mathrm{\Δ}}_r-X_g|},\frac{{\∂h}_r}{{\∂\mathrm{\Δ}}_v}=0_{3\×1}$

$\frac{{\∂h}_v}{{\∂\mathrm{\Δ}}_r}=(\frac{1}{\left|{\mathrm{\Δ}}_r\right|}-\frac{{\mathrm{\Δ}}_r{{\mathrm{\Δ}}_r}^T}{{\left|{\mathrm{\Δ}}_r\right|}^3}){\mathrm{\Δ}}_v+(\frac{1}{|{\mathrm{\Δ}}_r-X_g|}-\frac{({\mathrm{\Δ}}_r-X_g){({\mathrm{\Δ}}_r-X_g)}^T}{{|{\mathrm{\Δ}}_r-X_g|}^3})({\mathrm{\Δ}}_r-X_g)$

$\frac{{\∂h}_v}{{\∂\mathrm{\Δ}}_v}=\frac{{\mathrm{\Δ}}_r}{\left|{\mathrm{\Δ}}_r\right|}-\frac{{\mathrm{\Δ}}_r-X_g}{|{\mathrm{\Δ}}_r-X_g|}$

④状态更新。在时间标签为k的时刻，得到测量值z(z_abs和z_rel，分别对应于绝对状态和相对状态估计器)，利用测量值z_k，按下式进行状态更新

${\hat{x}}_k={\hat{x}}_{k,k-1}+K_k(z_k-h\left({\hat{x}}_{k,k-1}\right))$

P_k＝(I-K_kH_k)P_k,k-1(I-K_kH_k)^T+K_kR_kK_k ^T

⑤时间更新，递归计算。

将时间标签由k-1修改为k，重复步骤①～⑤

步骤4-3

由根据步骤4-2的结果，计算两个航天器间的相对位置和相对速度。具体方法为：在任意时间标签为k的时刻，从绝对状态向量中取出r_c,v_c,r_T,v_T，由相对状态向量中取出Δ_r,Δ_v。则当两航天器相距较远时，两航天器在地心惯性系中的相对位置确定为r_c-r_T，相对速度确定为v_c-v_T；当两航天器相距较近时，两航天器在目标航天器轨道参考系中的相对位置确定为Δ_r，相对速度确定为Δ_v。

这里的远近d的大小取决于主动航天器的轨道高度s，满足d≤0.01s，通常d的取值范围为1～100千米。

本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。

Claims (3)

1.基于GNSS反射信号的双站雷达空间目标相对导航方法，其特征在于包括如下步骤：

(1)在主动航天器上安装GNSS卫星导航接收机，利用卫星导航接收机同时接收GNSS卫星的直达信号，以及目标航天器反射的GNSS卫星的回波信号；

(2)主动航天器获取自身的位置和速度，作为惯性坐标系下主动航天器的绝对运动状态变量，同时估计目标航天器的位置和速度，作为惯性坐标系下目标航天器的绝对运动状态变量；将主动航天器的绝对运动状态变量和目标航天器的绝对运动状态变量进行差分，得到惯性坐标系下主动航天器和目标航天器之间的相对位置和相对速度，作为惯性坐标系下的相对运动状态变量；

(3)判断主动航天器和目标航天器之间的距离，当主动航天器和目标航天器之间距离超过d时，转步骤(4)；当主动航天器和目标航天器之间的距离不超过d时，转步骤(6)；所述d的大小取决于主动航天器的轨道高度s，满足d≤0.01s；

(4)将步骤(2)中获取的惯性坐标系下的全部绝对运动状态变量，以及主动航天器与GNSS卫星之间的钟差b作为绝对状态向量x，采用迭代的方式，根据k-1时刻绝对状态向量的估计值迭代更新k时刻绝对状态向量的估计值具体为：

(4-1)计算得到k-1时刻绝对状态向量中各绝对状态变量的导数值，其中：

$\left\{\begin{array}{c}{\stackrel{\·}{r}}_c=v_c\\ {\stackrel{\·}{v}}_c=-\frac{\mathrm{\μ}}{r_c^3}\·r_c+u_c\end{array}\right.$

$\left\{\begin{array}{c}{\stackrel{\·}{r}}_T=v_T\\ {\stackrel{\·}{v}}_T=-\frac{\mathrm{\μ}}{r_T^3}\·r_T+u_T\end{array}\right.$

$\stackrel{\·}{b}=0$

r_c和v_c分别为惯性坐标系下主动航天器的绝对位置和绝对速度，r_T和v_T分别为惯性坐标系下目标航天器的绝对位置和绝对速度，u_c和u_T分别为惯性坐标系下主动航天器和目标航天器所受外部合力，μ为地心引力常数；

(4-2)对步骤(4-1)计算得到的导数值，从k-1时刻到k时刻的时间区间进行积分，得到各绝对状态变量的变化量，在此基础上叠加k-1时刻相应绝对状态变量的估计值，作为k时刻绝对状态向量的预测值

(4-3)利用k时刻绝对状态向量的预测值预测得到k时刻所述直达信号和回波信号的预测值向量 $h_{\mathrm{abs}}\left({\hat{x}}_{k,k-1}\right)={[h_{\mathrm{\ρ}}\left({\hat{x}}_{k,k-1}\right),h_{\stackrel{\·}{\mathrm{\ρ}}}\left({\hat{x}}_{k,k-1}\right),h_{\mathrm{r\ρ}}\left({\hat{x}}_{k,k-1}\right),h_{r\stackrel{\·}{\mathrm{\ρ}}}\left({\hat{x}}_{k,k-1}\right)]}^T,$ 其中直达信号的预测值为：

$h_{\mathrm{\ρ}}\left({\hat{x}}_{k,k-1}\right)=|r_c-r_g|+b$

$h_{\stackrel{\·}{\mathrm{\ρ}}}\left({\hat{x}}_{k,k-1}\right)={(v_c-v_g)}^T\frac{r_c-r_g}{|r_c-r_g|}$

回波信号的预测值为：

$h_{\mathrm{r\ρ}}\left({\hat{x}}_{k,k-1}\right)=|r_c-r_T|+|r_T-r_g|-|r_g-r_c|$

$h_{r\stackrel{\·}{\mathrm{\ρ}}}\left({\hat{x}}_{k,k-1}\right)={(v_c-v_T)}^T\frac{r_c-r_T}{|r_c-r_T|}+{(v_T-v_g)}^T\frac{r_T-r_g}{|r_T-r_g|}$

r_g和v_g分别表示GNSS卫星在惯性坐标系下的绝对位置和绝对速度；

(4-4)将步骤(4-3)获取的直达信号和回波信号作为预测值，将步骤(1)中获取的直达信号和回波信号作为测量值，采用卡尔曼滤波的方法，计算得到绝对状态向量的修正值

(4-5)通过公式计算得到k时刻绝对状态向量的估计值 ${\hat{x}}_k;$

(5)从步骤(4)得到的绝对状态向量估计值中提取出对应的分量，分别得到惯性坐标系下主动航天器和目标航天器各自的绝对位置和绝对速度，再通过直接差分，得到惯性坐标系下主动航天器和目标航天器之间的相对位置和相对速度，作为相对导航的结果，并结束；

(6)将步骤(2)获取的惯性坐标系下的相对运动状态变量转换至轨道坐标系，将轨道坐标系下的全部相对运动状态变量作为相对状态向量X，采用迭代的方式，采用迭代的方式，根据k-1时刻相对状态向量的估计值迭代更新k时刻相对状态向量的估计值具体为：

(6-1)计算得到k-1时刻相对状态向量中各相对状态变量的导数值，其中：

$\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\Δ}}}_r={\mathrm{\Δ}}_v\\ {\stackrel{\·}{\mathrm{\Δ}}}_v=A\·{\mathrm{\Δ}}_r+B\·{\mathrm{\Δ}}_v+u_{\mathrm{\Δ}}\end{array}\right.$

Δ_r和Δ_v分别为轨道坐标系下主动航天器相对于目标航天器的相对位置和相对速度，u_Δ为轨道坐标系下主动航天器相对于目标航天器的相对外部合力，

$A=\left[\begin{array}{ccc}{3\mathrm{\ω}}^2& 0& 0\\ 0& 0& 0\\ 0& 0& -{\mathrm{\ω}}^2\end{array}\right]$

$B=\left[\begin{array}{ccc}0& 2\mathrm{\ω}& 0\\ -2\mathrm{\ω}& 0& 0\\ 0& 0& 0\end{array}\right]$

ω为主动航天器的平均轨道角速度；

(6-2)对步骤(6-1)计算得到的导数值，从k-1时刻到k时刻的时间区间进行积分，得到各相对状态变量的变化量，在此基础上叠加k-1时刻相应相对状态变量的估计值，作为k时刻相对状态向量的预测值

(6-3)利用k时刻相对状态向量的预测值预测得到k时刻所述直达信号和回波信号差值的预测值向量其中

$h_r\left({\hat{X}}_{k,k-1}\right)=\left|{\mathrm{\Δ}}_r\right|+|{\mathrm{\Δ}}_r-X_g|-\left|X_g\right|$

$h_v\left({\hat{X}}_{k,k-1}\right)={{\mathrm{\Δ}}_v}^T\frac{{\mathrm{\Δ}}_r}{\left|{\mathrm{\Δ}}_r\right|}+{({\mathrm{\Δ}}_v-V_g)}^T\frac{{\mathrm{\Δ}}_r-X_g}{|{\mathrm{\Δ}}_r-X_g|}$

X_g和V_g分别为GNSS卫星在轨道坐标系下的位置和速度；

(6-4)将步骤(6-3)获取的直达信号和回波信号的差值作为预测值，将步骤(1)中获取的直达信号和回波信号做差后的差值作为测量值，采用卡尔曼滤波的方法，计算得到相对状态向量的修正值

(6-5)通过公式计算得到k时刻相对状态向量的估计值 ${\hat{x}}_k.$

(7)从步骤(6)得到的相对状态向量估计值中提取出对应的分量，得到轨道坐标系下主动航天器和目标航天器的相对位置和相对速度，作为相对导航的结果并结束。

2.根据权利要求1所述的基于GNSS反射信号的双站雷达空间目标相对导航方法，其特征在于：所述步骤(1)中利用卫星导航接收机的左旋天线接收GNSS卫星的直达信号，利用卫星导航接收机的右旋天线接收目标航天器反射的GNSS卫星的回波信号。

3.根据权利要求1或2所述的基于GNSS反射信号的双站雷达空间目标相对导航方法，其特征在于：所述的d的取值范围为1～100千米。