一种利用双阵列天线确定微小卫星相对姿态的方法与流程

本发明属于微小卫星编队相对姿态确定技术领域，具体涉及一种利用双阵列天线确定微小卫星相对姿态的方法。

背景技术：

微小卫星具有重量轻、体积小、技术含量高和研制周期短等一系列优点，除此之外，还可以采用模块化设计技术及标准化星体，在流水线上批量生产并储存、便于机动发射。还可以利用分布式星座、编队或者引入人工智能等新技术成果、用智能星群完成更加复杂的应用任务，而这种任务利用传统的大卫星是不能完成的。所以微小卫星在军用领域或者民用领域，包括商业通信、空间科研、国防军事、行星探测等方面都有非常广泛的应用。微小卫星编队星间相对姿态确定是实现编队自主运行和各种空间任务的基础，如何设计适用于微小卫星的相对姿态确定方法变得尤为重要。

微小卫星体积小、重量轻，传统的姿态确定敏感器如星敏感器、不仅体积大，而且非常昂贵，增加了微小卫星设计成本，并不适合在微小卫星上使用。此外还有利用gps卫星进行相对姿态确定的方法，这种方法需要在微小卫星上安装gps接收机，违背了微小卫星精简设计的原则，另外gps信号容易受到外界干扰，并且当微小卫星在远地轨道运行的时候，也就是说当微小卫星高于gps卫星时，接收不到gps信号，这种方法就失效了。texasa&muniversity开发了一种基于视觉的导航传感器系统(nisnav)，这种系统在主从航天器上分别安装光学传感器和光源，利用光学敏感器感知光源，得到航天器到光源之间的视线向量。应为这种方法需要额外的为航天器安装敏感器及光源，所以不可避免的增加了航天器的体积，不适用与微小卫星，并且这种方法在主从航天器之间距离非常大的时候将不再适用。此外还有利用相机识别特征点进行姿态确定的方法，这种方法需要用到非常复杂和好事的图像处理方法，这种方法对光照要求也很高，在航天器位于阴影区的时候这种方法失效。

技术实现要素：

本发明所要解决的技术问题在于针对上述现有技术中的不足，提供一种利用双阵列天线确定微小卫星相对姿态的方法，运算复杂度低，实时性好，设计简单、应用范围广、利用微小卫星现有通信模块即可进行微小卫星编队相对姿态确定。

本发明采用以下技术方案：

一种利用双阵列天线确定微小卫星相对姿态的方法，使用相对姿态动力学方程作为系统方程通过探测信号组成探测帧，然后设计发射接收时序，主航天器c的两个阵列天线ant1和ant2接收来自从航天器阵列天线的发射信号，得到输出信号，根据输出信号与探测信号做滑动相关，得到阵列冲激响应矩阵；结合波达角空间谱公式，进行二维谱峰搜索得到两组波离角和波达角；利用两组波达角和波离角得到四个单位视线向量uk1，uk2，u'k1和u'k2；利用双向单位视线向量与相对姿态旋转矩阵构造测量方程z，求出测量噪声方差阵rz；利用triad算法得到初始时刻的姿态四元数qk0；然后根据系统方程和测量方程z，结合测得的初值xk0以及测量噪声方差阵进行扩展卡尔曼滤波得到逐渐收敛的相对姿态和相对旋转角速度。

具体的，确定系统方程具体为：主航天器上安装两个阵列天线作为接收端，从航天器上安装天线阵列作为发射端，建立主从航天器本体坐标系及测量坐标系，从航天器dk的本体坐标系相对于主航天器c的本体坐标系的相对姿态四元数为qk＝[q1k,q2k,q3k,q4k]^t，相对旋转角速度为[wxk,wyk,wzk]^t，根据姿态动力学方程得到从航天器dk与主航天器c之间的系统方程

进一步的，系统方程具体为：

其中，第k个从航天器相关的状态量xk为xk＝[q1k,q2k,q3k,q4k,wxk,wyk,wzk]^t，[ixk,iyk,izk]^t为第k个从航天器的转动惯量。

具体的，结合波达角空间谱公式，进行二维谱峰搜索得到波离角和波达角具体为：根据不同发射信号，利用码分多址技术对各个从航天器进行辨识，得到各个从航天器的信道冲激响应向量，根据冲激响应向量分别确定各个从航天器与主航天器的两个阵列天线之间通信链路直达径的波达角和波离角。

进一步的，设接收端收到包含了s条多径的信号，经过通道识别得到接收端的第n个振元接收到的来自发射天线第m个发射天线的输出信号i表示主航天器的第i个接收天线阵列；分别用与k个从航天器使用的基带探测信号a1(τ)，a2(τ)，a3(τ)…ak(τ)进行相关，对于每条多径信号得到相对应的冲激响应的幅值和相位，其中功率最大的与视线向量对应，对于m＝1,2...m，n＝1,2...n，得到两个与主航天器c上两个阵列天线ant1、ant2相关的m*n的冲激响应矩阵，利用m*n矩阵分别得到ant1与发射天线通信链路之间的波达角和波离角以及ant2与发射天线通信链路之间的波达角和波离角

具体的，利用双向单位视线向量与相对姿态旋转矩阵的关系构造测量方程z，求出测量噪声方差阵rz；

测量方程z如下：

其中，和分别为uk1和uk2的测量误差，和分别为u'k1和u'k2的测量误差；

测量噪声方差阵rz如下：

其中，为四个单位向量uk1，uk2，u'k1和u'k2的测量噪声方差矩阵。

进一步的，令uk1＝[αk1,βk1,λk1]^t为从主航天器第一个阵列天线到第k个从航天器的阵列天线之间的单位视线向量；uk2＝[αk2,βk2,λk2]^t为从主航天器第二个阵列天线到第k个从航天器的阵列天线之间的单位视线向量；u'k1＝[α'k1,β'k1,λ'k1]^t为从第k个从航天器的阵列天线到主航天器第一个阵列天线之间的单位视线向量；u'k2＝[α'k2,β'k2,λ'k2]^t为从第k个从航天器的阵列天线到主航天器第二个阵列天线之间的单位视线向量，四个单位视线向量计算方程具体为：

四个单位向量uk1，uk2，u'k1和u'k2的测量噪声方差矩阵表达式如下：

其中，和分别为uk1和uk2的测量误差，和分别为u'k1和u'k2的测量误差。

具体的，利用triad算法得到扩展卡尔曼滤波的初始状态具体为：令s＝[t1,t2,t3]，b＝[b1b2b3]，利用四个单位向量构造正交基底，计算得到利用得到初始时刻的相对四元数qk0，记为[q1k0,q2k0,q3k0,q4k0]^t。

进一步的，将初始时刻用阵列向量测得的四个单位向量记为和和构造一组正交基t1,t2,t3如下：

和构造一组正交基b1b2b3如下：

具体的，分别对各个从航天器的相对状态进行迭代，得到逐步收敛的各个从航天器的相对姿态，将系统方程测量方程z和滤波初值xk0以及测量噪声方差阵rz代入扩展卡尔曼滤波器中得到估计的第k个航天器的相对状态xk，令k＝1,2,lk，得到所有k个从航天器的相对姿态和相对旋转角速度。

与现有技术相比，本发明至少具有以下有益效果：

本方法一种利用双阵列天线确定微小卫星相对姿态的方法，在主航天器上安装两个阵列天线，构造两组旋转矩阵方程作为测量方程，增加系统可观性，通过加入los径的波达角和波离角信息，采用天线阵列估计出的信号方向解决航天器间相对姿态估计问题，利用伪随机序列的基带探测信号组成探测帧，并在其两端分别加上保护头和保护尾，使用bpsk对基带探测信号进行调制得到多径信号，方法简单，测量精度高，利用航天器自带的通信数传模块即可进行los向量的测量，不需额外姿态测量模块，特别能够减轻微小卫星负载，提高空间使用率。特别适用于微小卫星自主相对姿态的确定。

进一步的，根据姿态动力学方程得到从航天器dk与主航天器c之间的系统方程，作为滤波器的系统方程，用于预测下一时刻的姿态。

进一步的，利用码分多址技术能够尽可能地减少外界干扰，得到更加真实的信道冲激响应，减小视线向量的估计误差，提高姿态估计精度。

进一步的，利用冲激响应确定直达径的波达角和波离角可以减小多径对直达径带来的干扰，提高估计精度。

进一步的，利用双向视线向量与相对姿态旋转矩阵关系构造测量方程和测量方差作为卡尔曼滤波器的测量方程，对姿态估计值进行更新。

进一步的，利用triad算法得到滤波器的初始状态让滤波初值与真实状态更加接近，可以加快滤波器收敛的速度，提高估计精度。

进一步的，只需要利用阵列天线测得los向量即可进行相对姿态的估计，使滤波器设计简单，设计复杂度降低。

综上所述，本发明运算复杂度低，估计精度高，适用范围广，不过分依赖外界环境条件，利用cdma技术对多个从航天器进行辨识，实现同时对多个从航天器进行相对姿态确定的目的。

下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。

附图说明

图1为本发明的编队构形示意图；

图2为本发明的双天线阵列面板示意图

图3为本发明的其中一个从航天器与主航天器两个阵列天线之间波达角、波离角和单位视线向量示意图；

图4为本发明俯仰角和方位角示意图。

具体实施方式

本发明提供了一种利用双阵列天线确定微小卫星相对姿态的方法，先在主航天器上安装两个阵列天线，每个从航天器上安装接收阵列天线，建立坐标系，根据相对姿态动力学方程构造系统方程；设计各个从航天器的阵列天线发射信号及发射接收时序；利用码分多址技术对各个航天器进行辨识，据输出信号与对应探测信号做滑动相关，分别得到某时刻各个从航天器的los径冲激响应向量，分别确定各个从航天器到主航天器两个阵列天线的波达角和波离角；利用波达角和波离角得到相应的单位视线向量；利用单位视线向量与相对姿态旋转矩阵的关系构造测量方程，并求得测量噪声方差矩阵；利用triad算法得到初始时刻的相对姿态四元数。根据系统方程，测量方程，测量噪声方差阵以及初始时刻的姿态四元数初值进行扩展卡尔曼滤波得到逐渐收敛的相对姿态和相对旋转角速度。

利用航天器现有的通信模块进行相对姿态的估计，尽可能地减少微小卫星重量，提高微小卫星的空间使用率。利用阵列天线测得的双向单位视线向量与姿态旋转矩阵的关系构造测量方程，可以增加可观测性。利用冲激响应得到相对距离和波达角，运算复杂度低，估计精度高。利用码分多址技术对各个从航天器进行辨识，可以同时进行针对多个从航天器的相对姿态的确定。

本发明一种利用双阵列天线确定微小卫星相对姿态的方法，包括以下步骤：

s1、航天器测量系统设计，即在主航天器上安装两个阵列天线作为接收端，从航天器上安装天线阵列作为发射端，建立主从航天器本体坐标系及测量坐标系，不同的从航天器设计不同的发射信号；

请参阅图1所示，编队构形包括一个主航天器和k个从航天器。主航天器c上安装两个接收天线阵列和接收系统接收探测信号，k个从航天器上分别安装发射天线系统发射探测信号，天线阵列的测量坐标系与主航天器的本体坐标系的相对关系已知。

请参阅图2所示，主航天器c上安装两个阵列天线，记作ant1和ant2，这两个天线阵列在主航天器本体坐标系中的位置已知。从航天器dk到两个阵列天线分别有一个视线向量。

请参阅图2所示，以c和dk为例介绍发射系统和接收系统。主航天器c上安装两个接收天线阵列，包含n个振元，从航天器dk上安装发射天线，包含m个振元。其中，主航天器c的本体坐标系与其上安装的两个天线阵列的测量坐标系重合。

主航天器c的本体坐标系的原点在航天器几何中心，y轴(yb)的方向与天线阵列所在的平面垂直并指向该平面，z轴(zb)指向主航天器c的正上方，x轴(xb)由右手准则给出，主航天器所安装的两个天线阵列的坐标系原点在天线阵列的几何中心，其x轴(xs)，y轴(ys)和z轴(zs)分别与xb，yb和zb平行，指向也相同；

主航天器c的轨道坐标系原点在航天器几何中心，x轴(xo)的方向与地心到航天器的方向相同，y轴(yo)在轨道平面内，与速度指向相同，垂直于x轴，z轴(xb)垂直于轨道平面，由右手准则给出。

长度为x的伪随机序列作为从航天器dk的基带探测信号ak(τ)，其表达式为：

其中，τ表示时间，表示宽度为tb的矩形脉冲信号，x为序列长度，n为伪随机序列ak(τ)的长度序号。

对于不同的从航天器，选择不同的伪随机序列作为各个从航天器的基带探测信号，即a1(τ)，a2(τ)…ak(τ)，进行发射，各个探测信号互相正交。

对于从航天器dk，l个伪随机序列组成一个探测帧uk(τ)，其表达式为：

其中tp＝xtb，l为pn序列个数。

在探测帧u(τ)两端分别加上保护头和保护尾作为本方法中的基本探测信号，该探测信号经过bpsk调制后经m个发射端天线发射出去。在接收端，两个阵列天线的n个天线振元同时接收信号。

s2、根据步骤s1中主从航天器的本体坐标系构建主航天器与各个从航天器的相对姿态动力学方程作为扩展卡尔曼滤波器的系统方程；

请参阅图3所示，从航天器dk的本体坐标系相对于主航天器的本体坐标系的相对姿态四元数为qk＝[q1k,q2k,q3k,q4k]^t，为相对姿态四元数相对与时间的导数，第k个从航天器相对与主航天器的相对旋转角速度为[wxk,wyk,wzk]^t，[ixk,iyk,izk]^t为第k个从航天器的转动惯量，根据姿态动力学方程，可得从航天器dk与主航天器c之间的系统方程。

所述系统方程具体为：

其中，xk＝[q1k,q2k,q3k,q4k,wxk,wyk,wzk]^t，

s3、根据步骤s1中设计的不同发射信号，利用码分多址技术对各个从航天器进行辨识，得到各个从航天器的信道冲激响应向量，根据冲激响应向量分别确定各个从航天器与主航天器的两个阵列天线之间通信链路直达径的波达角和波离角；

分别用从航天器本地伪随机序列与主航天器接两个天线收到的信号进行滑动相关，对于主航天器上的每个天线阵列，可以得到k个从航天器的信道冲激响应向量，第k个从航天器相对应的冲激响应中功率最大的径即为los径。

假设接收端收到了包含了s条多径的信号，经过通道识别，可以得到接收端的第n个振元接收到的来自发射天线第m个发射天线的输出信号其中下标i表示主航天器的第i个接收天线阵列。

分别用与k个从航天器使用的基带探测信号，即a1(τ)，a2(τ)，a3(τ)…ak(τ)进行相关，由与伪随机序列很强的自相关特性，对于每条多径信号都可以得到相对应的冲激响应的幅值和相位，其中功率最大的与视线向量对应，对于m＝1,2...m，n＝1,2...n，可以得到两个与ant1和ant2有关的m*n的los径复数阵列冲激响应矩阵，通过提取复数阵列冲激响应阵列的某一行1*n构成的向量，利用空间谱进行二维谱峰搜索得到波达角和下标k表示第k个从航天器，1和2分别表示主航天器上的两个天线阵列。、

通过提取所述复数阵列冲激响应阵列的某一列m*1构成的向量，利用空间谱进行二维谱峰搜索得到波达角和下标中k表示第k个从航天器，1和2分别表示主航天器上的两个天线阵列。

利用这两个m*n矩阵可以分别得到ant1与发射天线通信链路之间的波达角和波离角，如图2所示的波达角和波离角以及ant2与发射天线通信链路之间的波达角和波离角，如图2所示的波达角和波离角

举例如下：主航天器c上安装两个接收阵列天线和接收机，从航天器d1，d2…dk上安装发射阵列天线和发射机，发射机发射信号，接收机接收信号。主航天器c本体坐标系与其上的两个阵列天线测量坐标系之间旋转矩阵为i3×3，主航天器的本体坐标系与其上的两个阵列天线的测量坐标系重合。本实例中航天器c上两个阵列天线的振元数量n＝4，从航天器上的阵列天线的振元数为m＝4，阵列天线的方向图在暗室中测得为4×180×180的矩阵q。

dk发射机使用长度为1023的m序列作为伪随机序列，基带探测信号ak(t)的码速率为62.5兆比特/秒，也即式(1)中的tb＝16ns，其中ns表示纳秒。一个探测帧u(t)由两个伪随机序列连接组成，也即式(2)中k＝2。探测帧通过bpsk调制，载波频率为2.6ghz。调制后的探测帧表示为u′(t)，航天器dk的天线发射u′k(t)。

其余从航天器d2，d3…dk使用的伪随机序列与d1使用的伪随机序列相互正交。本实例中令k＝3，即有三个从航天器。

s4、假设主航天器本体坐标系与其上安装的两个阵列天线的测量坐标系重合，利用步骤s3中确定的各个从航天器的波达角和波离角得到分别在各个阵列天线测量坐标系下的单位视线向量，利用各个视线向量与相对姿态旋转矩阵之间的关系构造扩展卡尔曼滤波器的测量方程，计算测量噪声方差阵；

请参阅图3所示，主航天器上阵列天线ant1到从航天器dk的单位视线向量在ant1测量坐标系中的表示为uk1，从航天器dk到主航天器上阵列天线ant1之间的单位视线向量在从航天器测量坐标系中的表示为u'k1，这两个向量之间的关系由两个坐标系之间的姿态旋转矩阵来进行转换，即

主航天器上阵列天线ant2到从航天器dk的单位视线向量在ant2测量坐标系中的表示为uk2，从航天器dk到主航天器上阵列天线ant2之间的单位视线向量在dk测量坐标系中的表示为u'k2，这两个向量之间的关系由两个坐标系之间的姿态旋转矩阵来进行转换，即

其中，为第k个从航天器测量坐标系到主航天器第一个阵列天线之间的相对姿态转换矩阵，为第k个从航天器测量坐标系到主航天器第二个阵列天线之间的相对姿态转换矩阵，因为主航天器的两个天线的测量坐标系与主航天器的本体坐标系重合，从航天器测量坐标系与其本体坐标系重合，所以第k个从航天器本体坐标系到主航天器本体坐标系的相对姿态转换矩阵。

这四个单位视线向量可以由s3中得到的波达角和波离角得到，请参阅图4所示，θ和分别为单位视线向量在坐标系中的俯仰角和方位角。

令uk1＝[αk1,βk1,λk1]^t为从主航天器第一个阵列天线到第k个从航天器的阵列天线之间的单位视线向量，该向量在主航天器第一个阵列天线测量坐标系中表示。

uk2＝[αk2,βk2,λk2]^t为从主航天器第二个阵列天线到第k个从航天器的阵列天线之间的单位视线向量，该向量在主航天器第二个阵列天线测量坐标系中表示。

u'k1＝[α'k1,β'k1,λ'k1]^t为从第k个从航天器的阵列天线到主航天器第一个阵列天线之间的单位视线向量，该向量在第k个从航天器的阵列天线测量坐标系中表示。

u'k2＝[α'k2,β'k2,λ'k2]^t为从第k个从航天器的阵列天线到主航天器第二个阵列天线之间的单位视线向量，该向量在第k个从航天器的阵列天线测量坐标系中表示。

所述单位视线向量计算方程具体为：

本发明中利用阵列天线得到的波达角和波离角的噪声为高斯白噪声，波达角和波离角的真实值与测量值之间存在测量误差，即：

其中，为角度的测量值，为角度的真实值，为测量误差，并且测量误差的协方差矩阵由下式得到：

其中δθ与为已知量，本实例中δθ与设为0.001弧度。

若不考虑噪声，单位视线向量之间的转换关系，即扩展卡尔曼滤波器的测量方程可以表示如下：

其中，表示姿态旋转矩阵由姿态四元数表示，具体为：

进一步的，如果考虑噪声，利用单位视线向量之间的转换关系，则扩展卡尔曼滤波器的测量方程可以表示如下：

其中和分别为uk1和uk2的测量误差，和分别为u'k1和u'k2的测量误差。

根据公式(6)(7)(8)(9)可以得到四个单位向量uk1，uk2，u'k1和u'k2的测量噪声方差矩阵，具体表达式如下：

测量方程(17)的测量噪声矩阵可以由下式表示：

s5、利用triad算法得到扩展卡尔曼滤波的初始状态；

对于初始时刻用阵列向量测得的四个单位向量，记为和

用和构造一组正交基t1,t2,t3,具体构造方式由下式得到：

用和构造一组正交基a1,a2,a3,具体构造方式由下式得到：

令s＝[t1t2t3]，b＝[b1b2b3]，根据公式和可以得到则根据可以反解得到初始姿态四元数qk0＝[q1k0,q2k0,q3k0,q4k0]^t。

s6、利用步骤s2、步骤s4和步骤s5中构造的扩展卡尔曼滤波器的系统方程和测量方程以及初始状态，分别对各个从航天器的相对状态进行迭代，得到逐步收敛的各个从航天器的相对姿态。

第k个从航天器相关的状态量具体为：

xk＝[q1k,q2k,q3k,q4k,wxk,wyk,wzk]^t(26)

系统方程具体为：

其中，

测量方程具体为：

系统敏感性矩阵具体为：

滤波初值为其中相对姿态四元数初值qk0＝[q1k0,q2k0,q3k0,q4k0]^t由步骤s5得到，wk0＝[wxk0,wyk0,wzk0]^t＝[0,0,0]^t。

给定初始状态和均方误差初值后，将上述系统方程、测量方程和滤波初值以及测量噪声方差阵rz代入扩展卡尔曼滤波器中得到估计的第k个航天器的相对状态xk，令k＝1,2,lk即可得到所有k个从航天器的相对姿态和相对旋转角速度。

本发明设计简单，运算复杂度较低，利用航天器现有模块进行相对姿态的确定，有效减小了微小卫星的重量，提高了微小卫星的空间使用率，实现了微小微型结构多功能。

以上内容仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明权利要求书的保护范围之内。