CN102980580A - 基于张量积多胞鲁棒h2滤波的无陀螺卫星姿态确定方法 - Google Patents

Abstract

本发明涉及一种基于张量积多胞鲁棒H2滤波的无陀螺卫星姿态确定方法，属于飞行器技术领域。本发明针对卫星姿态动力学方程和运动学方程的非线性特性，提出了基于张量积转换的多胞鲁棒H2滤波，将非线性滤波问题转化为线性滤波问题；首先建立姿态确定系统的状态方程和星敏感器测量方程，并利用雅可比线性化将非线性系统化为线性变参数误差系统，然后根据张量积模型转换建立LPV系统多胞模型描述，并结合鲁棒H2滤波获取姿态确定系统的状态估计校正量，最后利用估计校正量对EKF方法获取的姿态一步预测量进行校正，获取姿态估计值；避免了EKF方法中实时计算更新滤波增益，大大减小了滤波计算量。

Description

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

技术领域

本发明涉及一种基于张量积多胞鲁棒H2滤波的无陀螺卫星姿态确定方法，属于飞行器技术领域。

背景技术

星敏感器是一种高精度的测量仪器，在卫星的姿态测量系统中已得到了广泛的应用。以往受灵敏度和带宽限制，星敏感器不能作为卫星在轨运行的主要姿态敏感器，需要与陀螺组合在一起才能建立高精度的姿态参考，确保卫星姿态信息的可靠性。近年来，基于陀螺加姿态敏感器的姿态确定组合方案在工程中得到了广泛的应用。但是陀螺价格昂贵，长期运行后可能会降低性能或出现故障，而且一些卫星受到功率和重量的限制不能装载陀螺组件。同时星敏感器在技术上取得了很大进步，具有宽视场、高灵敏度、低噪声等效角等特性，其灵敏度和带宽都得到了显著改善，因此仅利用姿态角敏感器进行卫星姿态确定的技术受到了越来越多的关注。对于无陀螺的姿态确定，扩展卡尔曼滤波（Extended KalmanFiltering,EKF）是常用的方法。该方法由于简单和容易实现的优点得到了广泛的应用。但是EKF算法需要实时计算雅可比矩阵并更新滤波增益，计算量大。

发明内容

本发明的目的是为解决无陀螺卫星非线性姿态确定算法复杂、计算量大等问题，提出一种基于张量积多胞鲁棒H2滤波的无陀螺卫星姿态确定方法。

本发明方法针对卫星姿态动力学方程和运动学方程的非线性特性，提出了基于张量积转换的多胞鲁棒H2滤波，将非线性滤波问题转化为线性滤波问题。本发明首先建立姿态确定系统的状态方程和星敏感器测量方程，并利用雅可比线性化将非线性系统化为线性变参数（Linear Parameter Varing,LPV）误差系统，然后根据张量积(Tensor Product,TP)模型转换建立LPV系统多胞模型描述，并结合鲁棒H2滤波获取姿态确定系统的状态估计校正量，最后利用估计校正量对EKF方法获取的姿态一步预测量进行校正，获取姿态估计值。

本发明的技术解决方案如下：

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

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

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

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

$J\stackrel{\·}{\mathrm{\ω}}+\mathrm{\ω}\×\mathrm{J\ω}=T---\left(1\right)$

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

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

$\stackrel{\·}{\mathrm{\σ}}=M\left(\mathrm{\σ}\right)\mathrm{\ω}---\left(2\right)$

式中， $M\left(\mathrm{\σ}\right)=\frac{1}{4}[(1-{\mathrm{\σ}}^T\mathrm{\σ})I_3+2[\mathrm{\σ}\×]+{2\mathrm{\σ\σ}}^T].$

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

S_b=R(σ)S_i+ΔS

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

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

设系统的状态估计 $\hat{x}={\left[\begin{array}{cc}{\widehat{\mathrm{\σ}}}^T& {\widehat{\mathrm{\ω}}}^T\end{array}\right]}^T,$ 其中

和

满足：

$\stackrel{\·}{\widehat{\mathrm{\ω}}}={-J}^{-1}\widehat{\mathrm{\ω}}\×J\widehat{\mathrm{\ω}}+J^{-1}\hat{T}---\left(3\right)$

$\stackrel{\·}{\widehat{\mathrm{\σ}}}=M\left(\widehat{\mathrm{\σ}}\right)\widehat{\mathrm{\ω}}---\left(4\right)$

结合雅可比线性化，式(1)和(3)、式(2)和(4)分别相减，得到：

$\stackrel{\·}{\mathrm{\ω}}-\stackrel{\·}{\widehat{\mathrm{\ω}}}=\mathrm{\Δ}\stackrel{\·}{\mathrm{\ω}}=F_{\mathrm{\ω}}\mathrm{\Δ\ω}+n_{\mathrm{\ω}}$

$\stackrel{\·}{\mathrm{\σ}}-\stackrel{\·}{\widehat{\mathrm{\σ}}}=\mathrm{\Δ}\stackrel{\·}{\mathrm{\σ}}=-[\widehat{\mathrm{\ω}}\×]\mathrm{\Δ\σ}+\frac{1}{4}\mathrm{\Δ\ω}$

式中： $F_{\mathrm{\ω}}=J^{-1}[\left(J\widehat{\mathrm{\ω}}\right)\×]-J^{-1}[\widehat{\mathrm{\ω}}\×]J,$ n_ω=J^-1ΔT，ΔT为力矩误差。

选取误差状态变量 $\mathrm{\Δx}={\left[\begin{array}{cc}\mathrm{\Δ}{\widehat{\mathrm{\σ}}}^T& {\mathrm{\Δ}\widehat{\mathrm{\ω}}}^T\end{array}\right]}^T,$ 则系统的滤波状态误差模型为

$\mathrm{\Δ}\stackrel{\·}{x}=\mathrm{F\Δx}+\mathrm{Gw}---\left(5\right)$

式中： $F=\left[\begin{array}{cc}-[\widehat{\mathrm{\ω}}\×]& \frac{1}{4}{I}_{3\×3}\\ 0_{3\×3}& F_{\mathrm{\ω}}\end{array}\right],$ $G=\left[\begin{array}{c}{0}_{3\×3}\\ I_{3\×3}\end{array}\right],$ w＝n_ω

星敏感器测量残差为

${\mathrm{\ΔS}}_b=S_b-{\hat{S}}_b$

$=R\left(\mathrm{\σ}\right)S_i-R\left(\widehat{\mathrm{\σ}}\right)S_i+\mathrm{\ΔS}$

$=[I-R\left(\mathrm{\Δ\σ}\right)]R\left(\widehat{\mathrm{\σ}}\right)S_i+\mathrm{\ΔS}$

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

${\mathrm{\ΔS}}_b=[4\left(R\left(\widehat{\mathrm{\σ}}\right)S_i\right)\×]\mathrm{\Δ\σ}+\mathrm{\ΔS}$

采用m个星敏感器对卫星姿态进行测量，得到的系统测量模型z＝[S_b1,S_b2,..,S_bm]^T，系统的测量误差模型为

$\mathrm{\Δz}=H\left(\widehat{\mathrm{\σ}}\right)\mathrm{\Δx}+v---\left(6\right)$

式中， $H\left(\widehat{\mathrm{\σ}}\right)=\left[\begin{array}{cc}4\left(R\left(\widehat{\mathrm{\σ}}\right)S_{i1}\right)\×& 0_{3\×3}\\ 4\left(R\left(\widehat{\mathrm{\σ}}\right)S_{i2}\right)\×& 0_{3\×3}\\ \·\·\·& \·\·\·\\ 4\left(R\left(\widehat{\mathrm{\σ}}\right)S_{\mathrm{im}}\right)\×& 0_{3\×3}\end{array}\right]$ $v=\left[\begin{array}{c}{\mathrm{\ΔS}}_1\\ {\mathrm{\ΔS}}_2\\ \·\·\·\\ {\mathrm{\ΔS}}_m\end{array}\right],$ S_bm为第m个星敏感器的测量模型，S_im为第m个星敏感器在惯性系中的单位方向矢量，ΔS_m为第m个星敏感器的测量噪声。

步骤3，根据步骤1中的卫星姿态动力学、运动学方程，得到变参数

的上下界，代入仿射参数依赖的矩阵F，得到其多胞顶点F₁,F₂。

对于非仿射参数依赖的矩阵

通过张量积模型转换方法获得LPV系统的测量误差模型的多胞描述,具体包括如下步骤：

步骤3.1，提取张量

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

的变参数空间P。

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

3）依据划分的网格分别求取矩阵

在各个网格点的值，并保存在张量H′中。

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

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

$H_{\left(1\right)}=\left(\begin{array}{cc}{U}_1& U_1^d\end{array}\right)\left(\begin{array}{cc}{D}_1& \\ & D_1^d\end{array}\right){\left(\begin{array}{cc}{V}_1& V_1^d\end{array}\right)}^T$

其中，对角矩阵D₁包含被保留的1模式非零奇异值

中元素是被舍弃的零奇异值。U₁、V₁和

分别为保留的奇异值对应的矩阵和舍弃的奇异值对应的矩阵。

令 $S_{\left(1\right)}=D_1{V}_1^T$ 则 $H_{\left(1\right)}={U_{1}D}_1{V}_1^T=U_1{S}_{\left(1\right)};$

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

步骤（1），如果 $\mathrm{\Σ}\left({\left(U_1^d\right)}^T\right)=0,$ 则令 ${\stackrel{\‾}{U}}_1=U_1{\mathrm{\φ}}_1,$ ∑(X)=∑((U₁)^T)，φ₁=X₁，其中X₁为任意矩阵；如果 $\mathrm{\Σ}({\left(U_1^d\right)}^T\≠0,$ 则令 ${\stackrel{\‾}{U}}_1=\left[\begin{array}{cc}{U}_1& U_1^d\mathrm{\Σ}\end{array}\left({\left(U_1^d\right)}^T\right)\right]{\mathrm{\φ}}_1,$ 其中 ${\mathrm{\φ}}_1=\left(\begin{array}{cc}{X}_1& 0\\ 0& 1\end{array}\right).$ 从而

每行的和值为1。

步骤（2），取

所包含元素的最小值α_min，令

$\mathrm{\ϵ}=\left\{\begin{array}{cc}1& ,{\mathrm{\α}}_{\min }\≥-1\\ \frac{1}{{\mathrm{\α}}_{\min }}& ,\mathrm{else}\end{array}\right.$

$Z=\frac{1}{n_{\mathrm{col}+1}}(\mathrm{\ϵI}+1)$

式中n_col为

的列数，I为单位矩阵。从而

中的任意元素大于0小于1。

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

$U_1\·S_{\left(1\right)}=U_1{\mathrm{\φ}}_1{\mathrm{ZZ}}^{-1}{\mathrm{\φ}}_1^{-1}{S}_{\left(1\right)}={\stackrel{\‾}{U}}_1\·{\stackrel{\‾}{S}}_{\left(1\right)}$

将矩阵

存储为张量

则张量H′表示为张量积形式为：

$H^{\′}={\stackrel{\‾}{S}}_1{\×}_1{\stackrel{\‾}{U}}_1$

步骤3.2.2，采用步骤3.2.1的方法，对张量

的2模式矩阵(S₁)₍₂₎进行奇异值分解，得到

${\left(S_1\right)}_{\left(2\right)}=U_2{D}_2{V}_2^T$

D₂中包含被保留的2模式非零奇异值；U₂、是相应的奇异值矩阵。将

保存成张量S₂，并在此基础上进行权系数标准化得到

$H^{\′}={\stackrel{\‾}{S}}_2{\×}_1{\stackrel{\‾}{U}}_1{\×}_2{\stackrel{\‾}{U}}_2$

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

${\left(S_2\right)}_{\left(3\right)}=U_3{D}_3{V}_3^T$

D₃中包含被保留的3模式非零奇异值；U₃、

是相应的奇异值矩阵。将保存成张量S₃，并在此基础上进行权系数标准化得到

$H^{\′}={\stackrel{\‾}{S}}_3{\×}_1{\stackrel{\‾}{U}}_1{\×}_2{\stackrel{\‾}{U}}_2{\×}_3{\stackrel{\‾}{U}}_3$

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

$H^{\′}={\stackrel{\‾}{S}}_3\underset{n=1}{\overset{3}{\⊗}}{\stackrel{\‾}{U}}_n,$ 且满足

${\left|\right|H\left(\widehat{\mathrm{\σ}}\right)-{\stackrel{\‾}{S}}_3\underset{n=1}{\overset{3}{\⊗}}{U}_n\left|\right|}^2\≤\underset{n=1}{\overset{3}{\mathrm{\Σ}}}\left(\underset{i_n=I_n+1}{\overset{R_n}{\mathrm{\Σ}}}{\left({\mathrm{\σ}}_{i_n}^{\left(n\right)}\right)}^2\right)---\left(7\right)$

其中

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

的值在一定小的范围内。

根据张量定义可知张量

可以转换为用多胞顶点矩阵H₁,H₂...H_s（s为多胞顶点数目，为1模式、2模式、3模式下分别保留的奇异值个数相乘）表示，从而得到系统滤波状态和测量误差方程的多胞描述形式：

$\mathrm{\Δ}\stackrel{\·}{x}=\mathrm{F\Δx}+\mathrm{Gw}$ (8)

Δz=HΔx+v

其中， $(F,H)=\{\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{\mathrm{\λ}}_{1i}{F}_i,\underset{i=1}{\overset{s}{\mathrm{\Σ}}}{\mathrm{\λ}}_{2i}{H}_i|\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{\mathrm{\λ}}_{1i}=1,\underset{i=1}{\overset{s}{\mathrm{\Σ}}}{\mathrm{\λ}}_{2i}=1\},$ (F_i,H_i)为多胞形系统的顶点。

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

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

Δx_k+1=AΔx_k+Bw_k

Δz_k=CΔx_k+Dw_k

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

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

${\hat{x}}_{\mathrm{mk}+1}=A_F{\hat{x}}_{\mathrm{mk}}+B_F(z_k-{\hat{z}}_{k|k-1})$

${\mathrm{\Δ}\hat{x}}_k=C_F{\hat{x}}_{\mathrm{mk}}+D_F(z_k-{\hat{z}}_{k|k-1})$

其中，

是k时刻的中间变量，

是k-1时刻测量模型到k时刻的一步测量预测，

是k时刻的状态估计值。滤波系数A_F、B_F、C_F、D_F计算公式为：

$A_F=G_2^{-1}{S}_A,$ $B_F=G_2^{-1}{S}_B,$ C_F=S_C、D_F=S_D

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

min trace(Z)

$s.t.\left[\begin{array}{ccc}{G}_{11}+G_{11}^T-P_{11j}& G_2+G_{21}^T-P_{12j}& {\mathrm{\ψ}}_{1j}\\ *& G_2+G_{21}^T-P_{22j}& {\mathrm{\ψ}}_{2j}\\ *& *& {\mathrm{\ψ}}_{3j}\\ *& *& *\\ *& *& *\end{array}\right.$

$\left.\begin{array}{cc}{S}_A-F_{21}^T& G_{11}{B}_j+S_B{D}_j\\ S_A-{\mathrm{\λ}}_2{G_2}^T& G_{21}{B}_j+S_B{D}_j\\ {\mathrm{\ψ}}_{4j}& {-F}_{11}{B}_j-{\mathrm{\λ}}_1{S}_B{D}_j\\ P_{22j}-{\mathrm{\λ}}_2{S}_A-{\mathrm{\λ}}_2{S_A}^T& -F_{21}{B}_j-{\mathrm{\λ}}_2{S}_B{D}_j\\ *& I\end{array}\right]>0,$

$\left[\begin{array}{cccc}Z& I-S_D{C}_j& {-S}_C& -S_D{D}_j\\ *& P_{11j}& P_{12j}& 0\\ *& *& P_{22j}& 0\\ *& *& *& I\end{array}\right]>0,j=\mathrm{1,2},\·\·\·,2s$

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

${\mathrm{\ψ}}_{1j}=G_{11}{A}_j+S_B{C}_j-F_{11}^T,$ ψ_2j=G₂₁A_j+S_BC_j-λ₁G₂ ^T

${\mathrm{\ψ}}_{3j}=P_{11j}-F_{11}{A}_j-{\mathrm{\λ}}_1{S}_B{C}_j-A_j^T{F}_{11}^T-{\mathrm{\λ}}_1{C}_j^T{S}_B^T$

${\mathrm{\ψ}}_{4j}=P_{12j}-{\mathrm{\λ}}_1{S}_A-A_j^T{F}_{21}^T-{\mathrm{\λ}}_2{C}_j^T{S}_B^T$

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

${\widehat{\mathrm{\σ}}}_{k|k-1}={\widehat{\mathrm{\σ}}}_{k-1}+M\left({\widehat{\mathrm{\σ}}}_{k-1}\right){\widehat{\mathrm{\ω}}}_{k-1}{T}_s$

${\widehat{\mathrm{\ω}}}_{k|k-1}={\widehat{\mathrm{\ω}}}_{k-1}+({-J}^{-1}\left(\right[{\widehat{\mathrm{\ω}}}_{k-1}\×\left]J{\widehat{\mathrm{\ω}}}_{k-1}\right)+J^{-1}\hat{T})T_s$

式中， 为k-1时刻角速度估计值和MRP估计值。

分别为k时刻角速度预测值和MRP预测值。

步骤6，利用步骤4得到的估计校正量

对步骤5得到的预测值

进行校正，获得k时刻角速度和MRP状态估计为：

${\widehat{\mathrm{\σ}}}_k=\mathrm{\Δ}{\widehat{\mathrm{\σ}}}_k\⊗{\widehat{\mathrm{\σ}}}_{k|k-1}$

${\widehat{\mathrm{\ω}}}_k={\widehat{\mathrm{\ω}}}_{k|k-1}+\mathrm{\Δ}{\widehat{\mathrm{\ω}}}_k$

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

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

有益效果

本发明方法避免了EKF方法中忽略高阶项引入的模型误差；将线性变参数误差系统结合TP模型转换描述其多胞形式，代替了直接用边界值描述多胞系统，降低了多胞模型的保守性；同时结合鲁棒H2滤波思想，在EKF一步预测的基础上，进行预测校正，避免了EKF方法中实时计算更新滤波增益，大大减小了滤波计算量。

附图说明

图1为具体实施方式中采用本发明方法的MRP姿态误差；

图2为具体实施方式中采用本发明方法的姿态角速度误差；

图3为具体实施方式中EKF方法3-1-2转换欧拉角稳态误差；

图4为具体实施方式中采用本发明方法3-1-2转换欧拉角稳态误差；

图5为具体实施方式中EKF姿态速度稳态误差；

图6为具体实施方式中采用本发明方法的姿态角速度稳态误差；

图7为本发明基于张量积多胞鲁棒H2滤波的无陀螺卫星姿态确定方法流程图。

具体实施方式

为了更好的说明本发明的目的和优点，下面结合附图和实施例加以进一步说明。

（1）建立卫星姿态确定的状态模型和测量模型

卫星的姿态确定任务也即确定卫星本体系相对于参考坐标系的方位。选取惯性系作为参考坐标系，则卫星姿态动力学方程

$J\stackrel{\·}{\mathrm{\ω}}+\mathrm{\ω}\×\mathrm{J\ω}=T$

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

采用修正罗德里格参数(MRP)作为姿态描述参数，则卫星姿态运动学方程为：

$\stackrel{\·}{\mathrm{\σ}}=M\left(\mathrm{\σ}\right)\mathrm{\ω}$

式中， $M\left(\mathrm{\σ}\right)=\frac{1}{4}[(1-{\mathrm{\σ}}^T\mathrm{\σ})I_3+2[\mathrm{\σ}\×]+{2\mathrm{\σ\σ}}^T].$

星敏感器测量模型如下：

S_b=R(σ)S_i+ΔS

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

（2）利用雅可比线性化将非线性系统转化为LPV线性系统。

经雅可比线性化，可得：

$\mathrm{\Δ}\stackrel{\·}{\mathrm{\ω}}=F_{\mathrm{\ω}}\mathrm{\Δ\ω}+n_{\mathrm{\ω}}$

$\mathrm{\Δ}\stackrel{\·}{\mathrm{\σ}}=-[\widehat{\mathrm{\ω}}\×]\mathrm{\Δ\σ}+\frac{1}{4}\mathrm{\Δ\ω}$

式中： $F_{\mathrm{\ω}}=J^{-1}[\left(J\widehat{\mathrm{\ω}}\right)\×]-J^{-1}[\widehat{\mathrm{\ω}}\×]J,$ n_ω=J^-1ΔT

选取误差状态变量为 $\mathrm{\Δx}={\left[\begin{array}{cc}\mathrm{\Δ}{\widehat{\mathrm{\σ}}}^T& {\mathrm{\Δ}\widehat{\mathrm{\ω}}}^T\end{array}\right]}^T,$ 则系统的误差状态方程为

$\mathrm{\Δ}\stackrel{\·}{x}=\mathrm{F\Δx}+\mathrm{Gw}---\left(5\right)$

式中： $F=\left[\begin{array}{cc}-[\widehat{\mathrm{\ω}}\×]& \frac{1}{4}{I}_{3\×3}\\ 0_{3\×3}& F_{\mathrm{\ω}}\end{array}\right],$ $G=\left[\begin{array}{c}{0}_{3\×3}\\ I_{3\×3}\end{array}\right],$ w＝n_ω

星敏感器测量残差为

${\mathrm{\ΔS}}_b=S_b-{\hat{S}}_b$

$=R\left(\mathrm{\σ}\right)S_i-R\left(\widehat{\mathrm{\σ}}\right)S_i+\mathrm{\ΔS}$

$=[I-R\left(\mathrm{\Δ\σ}\right)]R\left(\widehat{\mathrm{\σ}}\right)S_i+\mathrm{\ΔS}$

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

${\mathrm{\ΔS}}_b=[4\left(R\left(\widehat{\mathrm{\σ}}\right)S_i\right)\×]\mathrm{\Δ\σ}+\mathrm{\ΔS}$

为得出三轴姿态，为提高观测精度，本实施例采用双星敏感器观测，即z=[S_b1,S_b2]^T，由此可得系统的测量残差方程为

$\mathrm{\Δz}=H\left(\widehat{\mathrm{\σ}}\right)\mathrm{\Δx}+v$

式中： $\mathrm{\Δz}=\left[\begin{array}{c}{\mathrm{\ΔS}}_{b1}\\ {\mathrm{\ΔS}}_{b2}\end{array}\right],$ $H\left(\widehat{\mathrm{\σ}}\right)=\left[\begin{array}{cc}4\left(R\left(\widehat{\mathrm{\σ}}\right)S_{i1}\right)\×& 0_{3\×3}\\ 4\left(R\left(\widehat{\mathrm{\σ}}\right)S_{i2}\right)\×& 0_{3\×3}\end{array}\right],$ $v=\left[\begin{array}{c}{\mathrm{\ΔS}}_1\\ {\mathrm{\ΔS}}_2\end{array}\right]$

（3）根据张量积理论用将LPV系统表示为多胞系统

矩阵F是仿射参数依赖的，可以通过变参数的上下界组合得到其多胞顶点F₁,F₂，而矩阵

不是仿射参数依赖的，需要通过张量积模型转换方法获得多胞模型。

不是仿射参数依赖的，其变参数为

由卫星姿态动力学和运动学方程可得到

的上下界，即变参数空间。利用平均网格划分的方法，在变参数空间中取10×10×10个网格点。

1)求取

在各个网格点的值并分别保存在张量H′中

2)对张量H′进行高阶奇异值分解，并在其1模式、2模式、3模式奇异值分解中分别保留1、1、2个奇异值。得到多胞形系统的2个顶点H₁,H₂。

$H_1=\left[\begin{array}{cccccc}0& -0.05596& 3.4486& 0& 0& 0\\ 0.0596& 0& -4.5081& 0& 0& 0\\ -3.4486& 4.5081& 0& 0& 0& 0\\ 0& -4.0566& 3.5199& 0& 0& 0\\ 4.0566& 0& -4.4006& 0& 0& 0\\ -3.5199& 4.4006& 0& 0& 0& 0\end{array}\right]$

$H_2=\left[\begin{array}{cccccc}0& -0.0363& 4.5595& 0& 0& 0\\ 0.0363& 0& -3.4909& 0& 0& 0\\ -4.5595& 3.4909& 0& 0& 0& 0\\ 0& -4.0034& 4.6161& 0& 0& 0\\ 4.0334& 0& -3.3744& 0& 0& 0\\ -4.6161& 3.3744& 0& 0& 0& 0\end{array}\right]$

将多胞描述系统离散化，得到离散的系统状态方程和测量方程：

$\mathrm{\Δ}{\hat{x}}_k=\mathrm{A\Δ}{\hat{x}}_k+\mathrm{B\Δ}{n}_k$

$\mathrm{\Δ}{\hat{y}}_k=\mathrm{C\Δ}{\hat{x}}_k+\mathrm{D\Δ}{n}_k$

其中， $\mathrm{\Ω}=\left\{(A,C)\right|(A,C)=\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{\mathrm{\α}}_i,(A^{\left(i\right)},C^{\left(i\right)}),$ ${\mathrm{\α}}_i\≥0,\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{\mathrm{\α}}_i=1\},$ A⁽ⁱ⁾,C⁽ⁱ⁾为多胞形系统顶点。

(4)获取姿态估计校正量

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

${\hat{x}}_{\mathrm{mk}+1}=A_F{\hat{x}}_{\mathrm{mk}}+B_F\mathrm{\Δ}{z}_k$

${\mathrm{\Δ}\hat{x}}_k=C_F{\hat{x}}_{\mathrm{mk}}+D_F\mathrm{\Δ}{z}_k$

其中，滤波系数A_F、B_F、C_F、D_F计算公式为：

$A_F=G_2^{-1}{S}_A,$ $B_F=G_2^{-1}{S}_B,$ C_F=S_C、D_F=S_D

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

min trace(Z)

$s.t.\left[\begin{array}{ccc}{G}_{11}+G_{11}^T-P_{11j}& G_2+G_{21}^T-P_{12j}& {\mathrm{\ψ}}_{1j}\\ *& G_2+G_{21}^T-P_{22j}& {\mathrm{\ψ}}_{2j}\\ *& *& {\mathrm{\ψ}}_{3j}\\ *& *& *\\ *& *& *\end{array}\right.$

$\left.\begin{array}{cc}{S}_A-F_{21}^T& G_{11}{B}_j+S_B{D}_j\\ S_A-{\mathrm{\λ}}_2{G_2}^T& G_{21}{B}_j+S_B{D}_j\\ {\mathrm{\ψ}}_{4j}& {-F}_{11}{B}_j-{\mathrm{\λ}}_1{S}_B{D}_j\\ P_{22j}-{\mathrm{\λ}}_2{S}_A-{\mathrm{\λ}}_2{S_A}^T& -F_{21}{B}_j-{\mathrm{\λ}}_2{S}_B{D}_j\\ *& I\end{array}\right]>0,$

$\left[\begin{array}{cccc}Z& I-S_D{C}_j& {-S}_C& -S_D{D}_j\\ *& P_{11j}& P_{12j}& 0\\ *& *& P_{22j}& 0\\ *& *& *& I\end{array}\right]>0,j=\mathrm{1,2,3,4}$

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

${\mathrm{\ψ}}_{1j}=G_{11}{A}_j+S_B{C}_j-F_{11}^T,$ ψ_2j=G₂₁A_j+S_BC_j-λ₁G₂ ^T

${\mathrm{\ψ}}_{3j}=P_{11j}-F_{11}{A}_j-{\mathrm{\λ}}_1{S}_B{C}_j-A_j^T{F}_{11}^T-{\mathrm{\λ}}_1{C}_j^T{S}_B^T$

${\mathrm{\ψ}}_{4j}=P_{12j}-{\mathrm{\λ}}_1{S}_A-A_j^T{F}_{21}^T-{\mathrm{\λ}}_2{C}_j^T{S}_B^T$

(5)利用估计校正量对EKF算法获取的姿态一步预测量进行校正，获取姿态估计值。

对系统状态进行一步预测，得到：

${\widehat{\mathrm{\σ}}}_{k|k-1}={\widehat{\mathrm{\σ}}}_{k-1}+M\left({\widehat{\mathrm{\σ}}}_{k-1}\right){\widehat{\mathrm{\ω}}}_{k-1}{T}_s$

${\widehat{\mathrm{\ω}}}_{k|k-1}={\widehat{\mathrm{\ω}}}_{k-1}+({-J}^{-1}\left(\right[{\widehat{\mathrm{\ω}}}_{k-1}\×\left]J{\widehat{\mathrm{\ω}}}_{k-1}\right)+J^{-1}\hat{T})T_s$

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

从而获得k时刻角速度和MRP状态估计为：

${\widehat{\mathrm{\σ}}}_k=\mathrm{\Δ}{\widehat{\mathrm{\σ}}}_k\⊗{\widehat{\mathrm{\σ}}}_{k|k-1}$

${\widehat{\mathrm{\ω}}}_k={\widehat{\mathrm{\ω}}}_{k|k-1}+\mathrm{\Δ}{\widehat{\mathrm{\ω}}}_k$

本实施例通过在matlab仿真平台上进行实验，以验证本发明提出的基于张量积的多胞H2滤波具有良好的性能。

在本实施例中，星敏感器测量噪声标准差为20″。姿态角初始值设定为[2°3°5°]，姿态角速度初始值为0，干扰力矩标准差为10^-3N·m，估计初始值与真实值初始值取值相等，估计误差协方差阵初始值为10^-8I，仿真总时长为200s，修正罗德里格参数(MRP)与姿态角之间的转换按照3-1-2的顺序进行折合。

下面与采用EKF无陀螺姿态确定的结果进行对比，图1和图2是本方法的MRP姿态误差和角速度姿态误差结果，最后误差在零值附近小范围波动，说明本方法能够准确的估计出MRP和姿态角速度。姿态角速度和3-1-2转换欧拉角跟踪误差如图3—图6所示，图中分别给出了采用EKF方法和基于张量积多胞鲁棒H2滤波方法的3-1-2转换欧拉角和角速度的稳态误差结果，。可以看出，两种方法都在30~40S内稳定，基于张量积多胞H2滤波姿态确定的精度和EKF姿态确定的精度相比稍差一些，但是稳态时误差变化比较稳定，更不容易发散，并且不需要实时计算雅可比矩阵并实时更新滤波增益计算量明显比EKF时减小。

Claims (2)

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

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

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

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

$J\stackrel{\·}{\mathrm{\ω}}+\mathrm{\ω}\×\mathrm{J\ω}=T$

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

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

$\stackrel{\·}{\mathrm{\σ}}=M\left(\mathrm{\σ}\right)\mathrm{\ω}---\left(2\right)$ 式中， $M\left(\mathrm{\σ}\right)=\frac{1}{4}[(1-{\mathrm{\σ}}^T\mathrm{\σ})I_3+2[\mathrm{\σ}\×]+{2\mathrm{\σ\σ}}^T];$

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

S_b=R(σ)S_i+ΔS

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

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

设系统的状态估计 $\hat{x}={\left[\begin{array}{cc}{\widehat{\mathrm{\σ}}}^T& {\widehat{\mathrm{\ω}}}^T\end{array}\right]}^T,$ 其中

和

满足：

$\stackrel{\·}{\widehat{\mathrm{\ω}}}={-J}^{-1}\widehat{\mathrm{\ω}}\×J\widehat{\mathrm{\ω}}+J^{-1}\hat{T}---\left(3\right)$

$\stackrel{\·}{\widehat{\mathrm{\σ}}}=M\left(\widehat{\mathrm{\σ}}\right)\widehat{\mathrm{\ω}}---\left(4\right)$

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

$\stackrel{\·}{\mathrm{\ω}}-\stackrel{\·}{\widehat{\mathrm{\ω}}}=\mathrm{\Δ}\stackrel{\·}{\mathrm{\ω}}=F_{\mathrm{\ω}}\mathrm{\Δ\ω}+n_{\mathrm{\ω}}$

$\stackrel{\·}{\mathrm{\σ}}-\stackrel{\·}{\widehat{\mathrm{\σ}}}=\mathrm{\Δ}\stackrel{\·}{\mathrm{\σ}}=-[\widehat{\mathrm{\ω}}\×]\mathrm{\Δ\σ}+\frac{1}{4}\mathrm{\Δ\ω}$

式中： $F_{\mathrm{\ω}}=J^{-1}[\left(J\widehat{\mathrm{\ω}}\right)\×]-J^{-1}[\widehat{\mathrm{\ω}}\×]J,$ n_ω=J^-1ΔT，ΔT为力矩误差；

选取误差状态变量 $\mathrm{\Δx}={\left[\begin{array}{cc}\mathrm{\Δ}{\widehat{\mathrm{\σ}}}^T& \mathrm{\Δ}{\widehat{\mathrm{\ω}}}^T\end{array}\right]}^T,$ 则系统的滤波状态误差模型为

$\mathrm{\Δ}\stackrel{\·}{x}=\mathrm{F\Δx}+\mathrm{Gw}---\left(5\right)$

式中： $F=\left[\begin{array}{cc}-[\widehat{\mathrm{\ω}}\×]& \frac{1}{4}{I}_{3\×3}\\ 0_{3\×3}& F_{\mathrm{\ω}}\end{array}\right],$ $G=\left[\begin{array}{c}{0}_{3\×3}\\ I_{3\×3}\end{array}\right],$ w＝n_ω

星敏感器测量残差为

${\mathrm{\ΔS}}_b=S_b-{\hat{S}}_b$

$=R\left(\mathrm{\σ}\right)S_i-R\left(\widehat{\mathrm{\σ}}\right)S_i+\mathrm{\ΔS}$

$=[I-R\left(\mathrm{\Δ\σ}\right)]R\left(\widehat{\mathrm{\σ}}\right)S_i+\mathrm{\ΔS}$

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

${\mathrm{\ΔS}}_b=[4\left(R\left(\widehat{\mathrm{\σ}}\right)S_i\right)\×]\mathrm{\Δ\σ}+\mathrm{\ΔS}$

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

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

系统的测量误差模型为

$\mathrm{\Δz}=H\left(\widehat{\mathrm{\σ}}\right)\mathrm{\Δx}+v---\left(6\right)$

式中， $H\left(\widehat{\mathrm{\σ}}\right)=\left[\begin{array}{cc}4\left(R\left(\widehat{\mathrm{\σ}}\right)S_{i1}\right)\×& 0_{3\×3}\\ 4\left(R\left(\widehat{\mathrm{\σ}}\right)S_{i2}\right)\×& 0_{3\×3}\\ \·\·\·& \·\·\·\\ 4\left(R\left(\widehat{\mathrm{\σ}}\right)S_{\mathrm{im}}\right)\×& 0_{3\×3}\end{array}\right]$ $v=\left[\begin{array}{c}{\mathrm{\ΔS}}_1\\ {\mathrm{\ΔS}}_2\\ \·\·\·\\ {\mathrm{\ΔS}}_m\end{array}\right],$ S_bm为第m个星敏感器的测量模型，S_im为第m个星敏感器在惯性系中的单位方向矢量，ΔS_m为第m个星敏感器的测量噪声；

步骤3，根据步骤1中的卫星姿态动力学、运动学方程，得到变参数

的上下界，代入仿射参数依赖的矩阵F，得到其多胞顶点F₁,F₂；

对于非仿射参数依赖的矩阵

通过张量积模型转换方法获得LPV系统的测量误差模型的多胞描述；

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

$\mathrm{\Δ}\stackrel{\·}{x}=\mathrm{F\Δx}+\mathrm{Gw}$ (8)

Δz=HΔx+v

其中， $(F,H)=\{\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{\mathrm{\λ}}_{1i}{F}_i,\underset{i=1}{\overset{s}{\mathrm{\Σ}}}{\mathrm{\λ}}_{2i}{H}_i|\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{\mathrm{\λ}}_{1i}=1,\underset{i=1}{\overset{s}{\mathrm{\Σ}}}{\mathrm{\λ}}_{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)=\{\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{\mathrm{\λ}}_{1i}{A}_i,\underset{i=1}{\overset{s}{\mathrm{\Σ}}}{\mathrm{\λ}}_{2i}{C}_i|\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{\mathrm{\λ}}_{1i}=1,\underset{i=1}{\overset{s}{\mathrm{\Σ}}}{\mathrm{\λ}}_{2i}=1\},$ T_s为采样周期；Δx_k表示k时刻的状态误差，w_k为离散后的系统噪声，Δz_k表示k时刻的测量误差；

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

${\hat{x}}_{\mathrm{mk}+1}=A_F{\hat{x}}_{\mathrm{mk}}+B_F(z_k-{\hat{z}}_{k|k-1})$

${\mathrm{\Δ}\hat{x}}_k=C_F{\hat{x}}_{\mathrm{mk}}+D_F(z_k-{\hat{z}}_{k|k-1})$

其中，是中间变量，

是k-1时刻测量模型到k时刻的一步测量预测，

是k时刻的状态估计值；滤波系数A_F、B_F、C_F、D_F计算公式为： $A_F=G_2^{-1}{S}_A,B_F=G_2^{-1}{S}_B,$ C_F=S_C、D_F=S_D

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

min trace(Z)

$s.t.\left[\begin{array}{ccc}{G}_{11}+G_{11}^T-P_{11j}& G_2+G_{21}^T-P_{12j}& {\mathrm{\ψ}}_{1j}\\ *& G_2+G_{21}^T-P_{22j}& {\mathrm{\ψ}}_{2j}\\ *& *& {\mathrm{\ψ}}_{3j}\\ *& *& *\\ *& *& *\end{array}\right.$

$\left.\begin{array}{cc}{S}_A-F_{21}^T& G_{11}{B}_j+S_B{D}_j\\ S_A-{\mathrm{\λ}}_2{G_2}^T& G_{21}{B}_j+S_B{D}_j\\ {\mathrm{\ψ}}_{4j}& {-F}_{11}{B}_j-{\mathrm{\λ}}_1{S}_B{D}_j\\ P_{22j}-{\mathrm{\λ}}_2{S}_A-{\mathrm{\λ}}_2{S_A}^T& -F_{21}{B}_j-{\mathrm{\λ}}_2{S}_B{D}_j\\ *& I\end{array}\right]>0,$

$\left[\begin{array}{cccc}Z& I-S_D{C}_j& {-S}_C& -S_D{D}_j\\ *& P_{11j}& P_{12j}& 0\\ *& *& P_{22j}& 0\\ *& *& *& I\end{array}\right]>0,j=\mathrm{1,2},\·\·\·,2s$

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

${\mathrm{\ψ}}_{1j}=G_{11}{A}_j+S_B{C}_j-F_{11}^T,$ ψ_2j=G₂₁A_j+S_BC_j-λ₁G₂ ^T

${\mathrm{\ψ}}_{3j}=P_{11j}-F_{11}{A}_j-{\mathrm{\λ}}_1{S}_B{C}_j-A_j^T{F}_{11}^T-{\mathrm{\λ}}_1{C}_j^T{S}_B^T$

${\mathrm{\ψ}}_{4j}=P_{12j}-{\mathrm{\λ}}_1{S}_A-A_j^T{F}_{21}^T-{\mathrm{\λ}}_2{C}_j^T{S}_B^T$

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

${\widehat{\mathrm{\σ}}}_{k|k-1}={\widehat{\mathrm{\σ}}}_{k-1}+M\left({\widehat{\mathrm{\σ}}}_{k-1}\right){\widehat{\mathrm{\ω}}}_{k-1}{T}_s$

${\widehat{\mathrm{\ω}}}_{k|k-1}={\widehat{\mathrm{\ω}}}_{k-1}+({-J}^{-1}\left(\right[{\widehat{\mathrm{\ω}}}_{k-1}\×\left]J{\widehat{\mathrm{\ω}}}_{k-1}\right)+J^{-1}\hat{T})T_s$

式中，

为k-1时刻角速度估计值和MRP估计值；

分别为k时刻角速度预测值和MRP预测值；

步骤6，利用步骤4得到的估计校正量

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

${\widehat{\mathrm{\σ}}}_k=\mathrm{\Δ}{\widehat{\mathrm{\σ}}}_k\⊗{\widehat{\mathrm{\σ}}}_{k|k-1}$

${\widehat{\mathrm{\ω}}}_k={\widehat{\mathrm{\ω}}}_{k|k-1}+\mathrm{\Δ}{\widehat{\mathrm{\ω}}}_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_{\left(1\right)}=\left(\begin{array}{cc}{U}_1& U_1^d\end{array}\right)\left(\begin{array}{cc}{D}_1& \\ & D_1^d\end{array}\right){\left(\begin{array}{cc}{V}_1& V_1^d\end{array}\right)}^T$

其中，对角矩阵D₁包含被保留的1模式非零奇异值

中元素是被舍弃的零奇异值；U₁、V₁和

分别为保留的奇异值对应的矩阵和舍弃的奇异值对应的矩阵；

令 $S_{\left(1\right)}=D_1{V}_1^T$ 则 $H_{\left(1\right)}={U_{1}D}_1{V}_1^T=U_1{S}_{\left(1\right)};$

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

步骤（1），如果 $\mathrm{\Σ}\left({\left(U_1^d\right)}^T\right)=0,$ 则令 ${\stackrel{\‾}{U}}_1=U_1{\mathrm{\φ}}_1,$ ∑(X)=∑((U₁)^T)，φ₁=X₁，其中X₁为任意矩阵；如果 $\mathrm{\Σ}({\left(U_1^d\right)}^T\≠0,$ 则令 ${\stackrel{\‾}{U}}_1=\left[\begin{array}{cc}{U}_1& U_1^d\mathrm{\Σ}\end{array}\left({\left(U_1^d\right)}^T\right)\right]{\mathrm{\φ}}_1,$ 其中 ${\mathrm{\φ}}_1=\left(\begin{array}{cc}{X}_1& 0\\ 0& 1\end{array}\right);$ 得到的

每行的和值为1；

步骤（2），取

所包含元素的最小值α_min，令

$\mathrm{\ϵ}=\left\{\begin{array}{cc}1& ,{\mathrm{\α}}_{\min }\≥-1\\ \frac{1}{{\mathrm{\α}}_{\min }}& ,\mathrm{else}\end{array}\right.$

$Z=\frac{1}{n_{\mathrm{col}+1}}(\mathrm{\ϵI}+1)$

式中n_col为

的列数，I为单位矩阵；得到的

中任意元素大于0小于1；

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

$U_1\·S_{\left(1\right)}=U_1{\mathrm{\φ}}_1{\mathrm{ZZ}}^{-1}{\mathrm{\φ}}_1^{-1}{S}_{\left(1\right)}={\stackrel{\‾}{U}}_1\·{\stackrel{\‾}{S}}_{\left(1\right)}$

将矩阵存储为张量

则张量H′表示为张量积形式为：

$H^{\′}={\stackrel{\‾}{S}}_1{\×}_1{\stackrel{\‾}{U}}_1$

步骤3.2.2，采用步骤3.2.1的方法，对张量

的2模式矩阵(S₁)₍₂₎进行奇异值分解，得到

${\left(S_1\right)}_{\left(2\right)}=U_2{D}_2{V}_2^T$

D₂中包含被保留的2模式非零奇异值；U₂、

是相应的奇异值矩阵；将

保存成张量S₂，并在此基础上进行权系数标准化得到

$H^{\′}={\stackrel{\‾}{S}}_2{\×}_1{\stackrel{\‾}{U}}_1{\×}_2{\stackrel{\‾}{U}}_2$

对张量

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

${\left(S_2\right)}_{\left(3\right)}=U_3{D}_3{V}_3^T$

D₃中包含被保留的3模式非零奇异值；U₃、

是相应的奇异值矩阵；将

保存成张量S₃，并在此基础上进行权系数标准化得到

$H^{\′}={\stackrel{\‾}{S}}_3{\×}_1{\stackrel{\‾}{U}}_1{\×}_2{\stackrel{\‾}{U}}_2{\×}_3{\stackrel{\‾}{U}}_3$

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

$H^{\′}={\stackrel{\‾}{S}}_3\underset{n=1}{\overset{3}{\⊗}}{\stackrel{\‾}{U}}_n,$ 且满足

${\left|\right|H\left(\widehat{\mathrm{\σ}}\right)-{\stackrel{\‾}{S}}_3\underset{n=1}{\overset{3}{\⊗}}{U}_n\left|\right|}^2\≤\underset{n=1}{\overset{3}{\mathrm{\Σ}}}\left(\underset{i_n=I_n+1}{\overset{R_n}{\mathrm{\Σ}}}{\left({\mathrm{\σ}}_{i_n}^{\left(n\right)}\right)}^2\right)---\left(7\right)$

其中

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

Images