[发明专利]基于张量积多胞鲁棒H2滤波的无陀螺卫星姿态确定方法

权利要求书

1.基于张量积多胞鲁棒H2滤波的无陀螺卫星姿态确定方法，其特征在于：包括以下步骤：

步骤1，建立卫星姿态的状态模型和星敏感器测量模型；

选取惯性系作为参考坐标系，确定卫星本体系相对于惯性系的方位；

建立卫星姿态动力学方程：

Jω·+ω×Jω=T]]>

其中，ω为在卫星本体系中相对于惯性系的角速度；J为卫星惯性阵，T为施加到卫星上的转动力矩，包含控制力矩和作用在卫星上的外部各种干扰力矩；

采用修正罗德里格参数σ作为姿态描述参数,卫星姿态运动学方程为：

σ·=M(σ)ω---(2)]]>式中，M(σ)=14[(1-σTσ)I3+2[σ×]+2σσT];]]>

采用星敏感器测量得到的测量模型为：

S_b=R(σ)S_i+ΔS

其中，S_i为恒星在惯性系中的单位方向矢量，ΔS为星敏感器测量噪声，R(σ)为惯性系到卫星本体系的姿态转换矩阵；

步骤2，采用雅可比线性化将步骤1建立的非线性系统变换为LPV系统，得到卫星姿态滤波状态误差模型和测量误差模型；

设系统的状态估计x^=σ^Tω^TT,]]>其中和满足：

ω^·=-J-1ω^×Jω^+J-1T^---(3)]]>

σ^·=M(σ^)ω^---(4)]]>

结合雅可比线性化，得到：

ω·-ω^·=Δω·=FωΔω+nω]]>

σ·-σ^·=Δσ·=-[ω^×]Δσ+14Δω]]>

式中：Fω=J-1[(Jω^)×]-J-1[ω^×]J,]]>n_ω=J^-1ΔT，ΔT为力矩误差；

选取误差状态变量Δx=Δσ^TΔω^TT,]]>则系统的滤波状态误差模型为

Δx·=FΔx+Gw---(5)]]>

式中：F=-[ω^×]14I3×303×3Fω,]]>G=03×3I3×3,]]>w＝n_ω

星敏感器测量残差为

ΔSb=Sb-S^b]]>

=R(σ)Si-R(σ^)Si+ΔS]]>

=[I-R(Δσ)]R(σ^)Si+ΔS]]>

Δσ为小量时，忽略二阶量，R(Δσ)≈I-4[Δσ×]，则星敏感器的测量残差方程为：

ΔSb=[4(R(σ^)Si)×]Δσ+ΔS]]>

采用m个星敏感器对卫星姿态进行测量，得到系统测量模型为：

z＝[S_b1，S_b2，...,S_bm]^T，

系统的测量误差模型为

Δz=H(σ^)Δx+v---(6)]]>

式中，H(σ^)=4(R(σ^)Si1)×03×34(R(σ^)Si2)×03×3······4(R(σ^)Sim)×03×3]]>v=ΔS1ΔS2···ΔSm,]]>S_bm为第m个星敏感器的测量模型，S_im为第m个星敏感器在惯性系中的单位方向矢量，ΔS_m为第m个星敏感器的测量噪声；

步骤3，根据步骤1中的卫星姿态动力学、运动学方程，得到变参数的上下界，代入仿射参数依赖的矩阵F，得到其多胞顶点F₁,F₂；

对于非仿射参数依赖的矩阵通过张量积模型转换方法获得LPV系统的测量误差模型的多胞描述；

系统滤波状态和测量误差方程的多胞描述形式为：

Δx·=FΔx+Gw]]>(8)

Δz=HΔx+v

其中，(F,H)={Σi=12λ1iFi,Σi=1sλ2iHi|Σi=12λ1i=1,Σi=1sλ2i=1},]]>(F_i,H_i)为多胞形系统的顶点；s为多胞顶点数目，为1模式、2模式、3模式下分别保留的奇异值个数相乘；

步骤4，结合鲁棒H2滤波获取姿态确定系统的状态估计校正量；

将多胞描述形式离散化，得到系统滤波状态方程和测量误差方程的离散多胞描述模型：

Δx_k+1=AΔx_k+Bw_k

Δz_k=CΔx_k+Dw_k

其中，A＝I+FT_s，B、C、D为多胞描述形式中相应参数的离散值，且(A,C)={Σi=12λ1iAi,Σi=1sλ2iCi|Σi=12λ1i=1,Σi=1sλ2i=1},]]>T_s为采样周期；Δx_k表示k时刻的状态误差，w_k为离散后的系统噪声，Δz_k表示k时刻的测量误差；

根据姿态确定滤波误差系统的离散多胞描述模型，并利用基于LMI技术的鲁棒H2滤波原理，得到姿态估计校正量的计算公式为：

x^mk+1=AFx^mk+BF(zk-z^k|k-1)]]>

Δx^k=CFx^mk+DF(zk-z^k|k-1)]]>

其中，是中间变量，是k-1时刻测量模型到k时刻的一步测量预测，是k时刻的状态估计值；滤波系数A_F、B_F、C_F、D_F计算公式为：AF=G2-1SA,BF=G2-1SB,]]>C_F=S_C、D_F=S_D

其中，S_A、S_B、S_C、S_D、G₂通过求解下述优化问题得到：

min trace(Z)

s.t.G11+G11T-P11jG2+G21T-P12jψ1j*G2+G21T-P22jψ2j**ψ3j******]]>

SA-F21TG11Bj+SBDjSA-λ2G2TG21Bj+SBDjψ4j-F11Bj-λ1SBDjP22j-λ2SA-λ2SAT-F21Bj-λ2SBDj*I>0,]]>

ZI-SDCj-SC-SDDj*P11jP12j0**P22j0***I>0,j=1,2,···,2s]]>

其中，λ₁,λ₂，G₁₁,G₂₁,G₂,F₁₁,F₂₁,S_A,S_B,S_C,S_D,P_11j,P_12j,P_22j为变量；

ψ1j=G11Aj+SBCj-F11T,]]>ψ_2j=G₂₁A_j+S_BC_j-λ₁G₂^T

ψ3j=P11j-F11Aj-λ1SBCj-AjTF11T-λ1CjTSBT]]>

ψ4j=P12j-λ1SA-AjTF21T-λ2CjTSBT]]>

步骤5，对EKF算法获取的系统姿态进行一步预测，得到：

σ^k|k-1=σ^k-1+M(σ^k-1)ω^k-1Ts]]>

ω^k|k-1=ω^k-1+(-J-1([ω^k-1×]Jω^k-1)+J-1T^)Ts]]>

式中，为k-1时刻角速度估计值和MRP估计值；分别为k时刻角速度预测值和MRP预测值；

步骤6，利用步骤4得到的估计校正量对步骤5得到的预测值进行校正，获得k时刻角速度和MRP状态估计为：

σ^k=Δσ^k⊗σ^k|k-1]]>

ω^k=ω^k|k-1+Δω^k]]>

步骤7，令k=k+1，将步骤6获得的k时刻状态估计值代入步骤5，得到k+1时刻一步预测；将该预测值代入步骤4，得到k+1时刻的姿态估计校正量；再进行步骤6，得到k+1时刻的角速度和MRP状态估计；

在无陀螺卫星运行过程中，按照上述方法，实现卫星姿态的实时获取。

2.根据权利要求1所述的基于张量积多胞鲁棒H2滤波的无陀螺卫星姿态确定方法，其特征在于：所述张量积模型转换方法的具体步骤为：

步骤3.1，提取张量

1）确定LPV系统的变参数空间P：由卫星姿态动力学和运动学方程得到变参数的变参数空间P；

2）对变参数空间P进行网格划分；

3）依据划分的网格分别求取矩阵在各个网格点的值，并保存在张量H′中；

步骤3.2，利用高阶奇异值分解简化张量模型，获得多胞顶点；

步骤3.2.1，对张量H′的1模式矩阵H₍₁₎进行奇异值分解：

H(1)=U1U1dD1D1dV1V1dT]]>

其中，对角矩阵D₁包含被保留的1模式非零奇异值中元素是被舍弃的零奇异值；U₁、V₁和分别为保留的奇异值对应的矩阵和舍弃的奇异值对应的矩阵；

令S(1)=D1V1T]]>则H(1)=U1D1V1T=U1S(1);]]>

对1模式矩阵H₍₁₎进行权系数标准化，具体过程为：

步骤（1），如果Σ((U1d)T)=0,]]>则令U‾1=U1φ1,]]>∑(X)=∑((U₁)^T)，φ₁=X₁，其中X₁为任意矩阵；如果Σ((U1d)T≠0,]]>则令U‾1=U1U1dΣ((U1d)T)φ1,]]>其中φ1=X1001;]]>得到的每行的和值为1；

步骤（2），取所包含元素的最小值α_min，令

ϵ=1,αmin≥-11αmin,else]]>

Z=1ncol+1(ϵI+1)]]>

式中n_col为的列数，I为单位矩阵；得到的中任意元素大于0小于1；

步骤(3)，结合步骤（1）和（2）得到：

U1·S(1)=U1φ1ZZ-1φ1-1S(1)=U‾1·S‾(1)]]>

将矩阵存储为张量则张量H′表示为张量积形式为：

H′=S‾1×1U‾1]]>

步骤3.2.2，采用步骤3.2.1的方法，对张量的2模式矩阵(S₁)₍₂₎进行奇异值分解，得到

(S1)(2)=U2D2V2T]]>

D₂中包含被保留的2模式非零奇异值；U₂、是相应的奇异值矩阵；将保存成张量S₂，并在此基础上进行权系数标准化得到

H′=S‾2×1U‾1×2U‾2]]>

对张量的3模式矩阵(S₂)₍₃₎进行奇异值分解，得到

(S2)(3)=U3D3V3T]]>

D₃中包含被保留的3模式非零奇异值；U₃、是相应的奇异值矩阵；将保存成张量S₃，并在此基础上进行权系数标准化得到

H′=S‾3×1U‾1×2U‾2×3U‾3]]>

最终得到高阶奇异值分解结果：

H′=S‾3⊗n=13U‾n,]]>且满足

||H(σ^)-S‾3⊗n=13Un||2≤Σn=13(Σin=In+1Rn(σin(n))2)---(7)]]>

其中为n模式奇异值分解下的第i_n个奇异值，R_n为n模式下奇异值的总个数,I_n为n模式下保留的奇异值个数；当奇异值之间差值比较大时，舍弃部分非零奇异值。