CN103970964B - 一种挠性卫星模态参数在轨辨识方法 - Google Patents

Abstract

一种挠性卫星模态参数在轨辨识方法，本发明涉及卫星模态参数在轨辨识领域，本发明要解决无法建立精确的动力学模型，频域法时效性低，难以辨识密频模态参数及获取阻尼信息使得估算、综合的参数存在误差，辨识精度降低以及时域法无法准确对模型定阶的问题而提出的一种挠性卫星模态参数在轨辨识方法。该方法是通过1、采集力矩和角速度信息；2、确定离散状态空间系统矩阵A、B、C和D；3、确定辨识矩阵A、B、C和D；4、建立动力学和帆板振荡方程；5、获得模态参数与力矩到角速度的传递函数；6、将矩阵A、B、C和D转换为传递函数等步骤实现的。本发明应用于卫星模态参数在轨辨识领域。

Description

一种挠性卫星模态参数在轨辨识方法

技术领域

本发明涉及一种挠性卫星模态参数在轨辨识方法。

背景技术

随着航天事业的发展，当代卫星结构越来越复杂，除了卫星主体部分外，增加了太阳能帆板以及各种用途的挠性附件。人们希望新一代的卫星能够拥有更高精度的指向能力，以完成复杂的空间任务，这就相应的增加了对控制系统的要求。

同时，随着挠性附件的复杂化以及对控制系统要求的提高，挠性卫星准确建模就变的十分重要。当前各类航天器需要承担的任务越来越多样化，构型也越来越复杂，比如机身、机翼以及太阳能帆板等都有成千上万个零件构成。同时出于很多因素的考虑，希望航天器的质量尽可能的轻，因此通常都选择低密度材料来制造太阳能帆板等挠性附件，这使得航天器的力学特性更加的复杂，同时对航天器的控制技术也提出了更高的要求。有限元分析法可以将一个大型复杂的结构分解为规模若干个小规模的结构，可以适应各种复杂的形状。因此该方法成为最初工程上广泛使用的方法。

当前阶段，我们有限元分析商用软件ANSYS,NASTRAN等可以方便的进行分析，利用有限元分析结合地面试验时，需测量结构参数，对结构运动学方程进行简化，在使用有限元分析法并结合地面试验的过程中无法对结构参数进行准确的测量，导致无法建立精确的动力学模型；同时还需要设法模拟太空的真实环境，然而在地面会由于地心引力、空气阻力等很多地面因素的影响，使得地面试验无法真实的模拟卫星在轨状态的太空环境，同时这其中还包括估算、综合的各种参数存在着误差的问题，因此识别精度受到了严重影响。比如JPL发射的探索木星的伽利略卫星在地面实验室中的频率误差为20％，振型误差为50％，

目前在轨辨识模态参数方法主要有频域法和时域法。频域法目前在研究工作状态下的辨识模态参数的频域方法较少，如峰值提取法、频域分解法以及频域多参考点模态识别法等方法。然而，峰值提取法无法辨识密集模态以及系统阻尼。频域分解法仅当输入激励为白噪声时，可以辨识成正交特性的模态参数，实际应用受到了极大的限制。频域多参考点模态识别法求解较为简单，计算量相对较小。然而却不可避免的仍存在着频域法辨识的共性问题：对噪声敏感，涉及复杂的复数运算，辨识精度以及时效性都较低，实际应用受到了严重影响。

当前阶段使用频域法实现在轨辨识模态参数时需要采样平均，耗时较长，使得该方法的时效性降低；并且将采集到的数据由时域变换到频域，意味着辨识精度相应降低。此外，频域法采用的峰值提取法难以辨识密频模态参数以及获取阻尼信息。

现有时域辨识方法在实际应用中相对比较普遍，主要包括Ibrahim时域法(Ibrahim Time Domain,ITD)，多个参考点复指数法(PRCE)以及特征系统实现算法(ERA)等方法。然而，Ibrahim时域法虽然有较为精确的辨识结果，却对噪声十分敏感，严重影响了对特征值以及特征向量的辨识精度；多个参考点复指数法(PRCE)虽然辨识精度较Ibrahim时域法有了较大的改善，却是基于提高模态辨识需要数据量的规模来实现的，时效性显著降低。特征系统实现算法(ERA)相比上述提到的两种时域辨识方法，具有辨识精度高、抗扰能力强的特点。然而在应用这种方法时，如果传感器数量配置过多，则会耗费较多的时间处理数据，时效性相应的降低；或者是由于干扰或传感器自身精度问题造成的测量值不准确时，辨识精度将大打折扣。

利用特征系统实现算法(ERA)辨识模态参数方法与频域法相比较，时域法在处理手段和估计密集模态、获取阻尼方面均具有一定的优势。同时，时域法不必将测试信号在时域与频域之间转换，在数据精度方面比频域法有优势，相应的辨识精度有所提高；此外，时域法可以不必需要激励数据，可以从响应数据中直接辨识模态参数。然而，时域法在获得一定成果的同时也面临着一定的重要问题：时域法没有自己独立的辨识方法，并且对噪声十分敏感，无法准确对模型定阶。因此频域法和时域法在实际应用中都受到很大的限制。

发明内容

本发明的目的是为了解决目前在工程中常使用有限元分析并结合地面试验的过程中无法建立精确的动力学模型，综合的各种参数存在着误差，影响识别精度、频域法实现在轨辨识模态参数时需要采样平均，耗时较长，时效性低；采用的峰值提取法难以辨识密频模态参数以及获取阻尼信息，辨识精度相应降低，以及时域法没有独立的辨识方法，对噪声十分敏感，无法准确对模型定阶等问题而提出的一种挠性卫星模态参数在轨辨识方法。

上述的发明目的是通过以下技术方案实现的：

步骤一、在外界诸多干扰不影响挠性卫星姿态稳定的前提下，于每个采样周期T平均准确的采集到执行机构施加到挠性卫星体上的力矩和挠性卫星体相对惯性坐标系的角速度信息；

步骤二、将采集到的力矩和角速度信息代入子空间辨识算法，确定离散状态空间系统矩阵A和C并且计算出离散状态空间系统矩阵B和D；

步骤三、假设当前所需辨识系统的系统状态方程，确定辨识离散状态空间系统矩阵A、B、C和D；

步骤四、基于挠性卫星姿态角的小角度假设，建立挠性卫星动力学方程和帆板振荡方程；

步骤五、将建立的挠性卫星动力学方程和帆板振荡方程解耦，通过拉式变换，获得模态参数与力矩到角速度的传递函数；

步骤六、将步骤三得到的辨识离散状态空间系统矩阵A、B、C和D转换为步骤五获得的模态参数与力矩到角速度的传递函数；即完成了一种挠性卫星模态参数在轨辨识方法。

发明效果

本发明提出的是一种基于子空间辨识算法的挠性卫星在轨辨识模态参数的方法，与参考文献挠性卫星姿态快速机动算法研究中提出的辨识系统矩阵的方法的不同之处在于，本发明创新性的利用输入输出矩阵获得相关扩展可控矩阵Γ_i和扩展可观测矩阵Γ_i行空间的正交补矩阵，获得辨识系统矩阵：

假设：

系统矩阵B,D则可以由式(123)计算得到：

然后将系统矩阵转化为传递函数，并获取模态参数。从而证明了《挠性卫星姿态快速机动算法研究》获得辨识系统矩阵的方法不是唯一的，可以利用输入输出矩阵获得相关扩展可控矩阵Γ_i和扩展可观测矩阵Γ_i行空间的正交补矩阵的方法，来获得辨识系统矩阵，为获得辨识系统矩阵提供了另一种方法。

本发明在卫星在轨状态下获得在轨期间有效信息进行模态参数辨识，从而建立精确的动力学模型。相比于在地面利用有限元分析试验时，无法模拟太空的真实环境，受地心引力、空气阻力等很多原因导致的辨识精度受到影响的情况，而本发明辨识模态参数的误差相对较小，辨识精度较高从而能够近似模拟太空的真实环境。

本发明不苛求小的采样时间，同时对输入信号的要求，仅仅是随机产生即可，无需苛刻的白噪声、各态历经或脉冲激励等苛刻条件，利用卫星自身可以测量得到的力矩和角速度信息，基于子空间辨识算法抗干扰能力强，准确对模型定阶(模态参数的阶数)相应的使得辨识精度也较高。

本发明核心步骤基于数学理论以及缜密的数学推导，编程求解过程较为简单，计算量也较小，使得辨识的时效性大幅度提升，可以为卫星在轨期间的短时内的姿态控制律重构、振荡隔离等时效性要求较强的控制需求提供便利。

本发明基于合理的传感器设置，基于现有的星上设备在轨获得角速度和力矩等相关信息，即可准确辨识挠性卫星模态振荡比较关键的前三阶模态的振荡频率、阻尼系数等模态信息，工程上实现容易，拥有很大的实用价值。

本发明分别给出了基于惯性空间坐标系以及轨道坐标系的仿真验证，每组仿真中均利用两组不同的模态参数进行了验证。前三阶模态振动频率的辨识相对误差平均为：5.42％,2.73％,0.92％；前三阶模态阻尼的辨识相对误差平均为：12.04％,1.11％,1,61％。而JPL发射的探索木星的伽利略卫星在轨辨识时，频率误差为5％，振型误差为20％。通过上述直观的对比，进一步证明了本方案具有较高的辨识精度，以及广阔的实用前景。

附图说明

图1是具体实施方式一提出的一种挠性卫星模态参数在轨辨识方法流程图；

图2是具体实施方式二提出的j维空间正交投影示意图；

图3是具体实施方式二提出的j维空间斜投影示意图；

图4是具体实施方式三提出的线性非时变确定性系统示意图，其中，Δ表示延时单元；

图5是具体实施方式二提出的矩阵输入输出方程几何示意图。

具体实施方式

具体实施方式一：本实施方式的一种挠性卫星模态参数在轨辨识方法，按照以下步骤制备的：

步骤一、在外界诸多干扰不影响挠性卫星姿态稳定的前提下，于每个采样周期T平均准确的采集到执行机构施加到挠性卫星体上的力矩和挠性卫星体相对惯性坐标系的角速度信息；

步骤二、将采集到的力矩和角速度信息代入子空间辨识算法，确定离散状态空间系统矩阵A和C并且计算出离散状态空间系统矩阵B和D；

步骤三、假设当前所需辨识系统的系统状态方程，确定辨识离散状态空间系统矩阵A、B、C和D；

步骤四、由于挠性卫星在相对轨道系三轴姿态较为稳定的情况下，挠性卫星保持相对稳定的姿态角，基于挠性卫星姿态角的小角度假设，建立挠性卫星动力学方程和帆板振荡方程；

步骤五、将建立的挠性卫星动力学方程和帆板振荡方程解耦，通过拉式变换，获得模态参数与力矩到角速度的传递函数；

步骤六、将步骤三得到的辨识离散状态空间系统矩阵A、B、C和D转换为步骤五获得的模态参数与力矩到角速度的传递函数如图1；即完成了一种挠性卫星模态参数在轨辨识方法。

本实施方式效果：

本实施方式在卫星在轨状态下获得在轨期间有效信息进行模态参数辨识，从而建立精确的动力学模型。相比于在地面利用有限元分析试验时，无法模拟太空的真实环境，受地心引力、空气阻力等很多原因导致的辨识精度受到影响的情况，而本实施方式辨识模态参数的误差相对较小，辨识精度较高。

本实施方式不苛求小的采样时间，同时对输入信号的要求，仅仅是随机产生即可，无需苛刻的白噪声、各态历经或脉冲激励等苛刻条件，利用卫星自身可以测量得到的力矩和角速度信息，基于子空间辨识算法抗干扰能力强，准确对模型定阶(模态参数的阶数)相应的使得辨识精度也较高。

本实施方式核心步骤基于数学理论以及缜密的数学推导，编程求解过程较为简单，计算量也较小，使得辨识的时效性大幅度提升，可以为卫星在轨期间的短时内的姿态控制律重构、振荡隔离等时效性要求较强的控制需求提供便利。

本实施方式基于合理的传感器设置，基于现有的星上设备在轨获得角速度和力矩等相关信息，即可准确辨识挠性卫星模态振荡比较关键的前三阶模态的振荡频率、阻尼系数等模态信息，工程上实现容易，拥有很大的实用价值。

本实施方式分别给出了基于惯性空间坐标系以及轨道坐标系的仿真验证，每组仿真中均利用两组不同的模态参数进行了验证。前三阶模态振动频率的辨识相对误差平均为：5.42％,2.73％,0.92％；前三阶模态阻尼的辨识相对误差平均为：12.04％,1.11％,1,61％。而JPL发射的探索木星的伽利略卫星在轨辨识时，频率误差为5％，振型误差为20％。通过上述直观的对比，进一步证明了本方案具有较高的辨识精度，以及广阔的实用前景。

具体实施方式二：本实施方式与具体实施方式一不同的是：步骤二中将采集到的力矩和角速度信息代入子空间辨识算法，确定离散状态空间系统矩阵A和C并且计算出离散状态空间系统矩阵B和D的实现步骤为：

利用“子空间辨识”算法由可以获得的力矩和角速度信息作为输入输出数据构成数据矩阵；该矩阵中列阵包含挠性卫星模态参数的相关信息；通过对该矩阵的行、列阵进行几何变换，辨识出相关可观(可控)矩阵和相邻采样周期的两个状态向量；在此基础上根据状态方程和奇异值分解等方法求出系统矩阵，进而完成辨识模态参数包括两个步骤：首先，利用扩展可观测矩阵Γi计算系统矩阵A和C；然后再计算系统矩阵B和D；

(一)通过输入输出数据构造Hankel块矩阵；其中输入块Hankel矩阵可以由输入数据构造，具体构成形式如下所示：

其中，

(1)i表示行块数，要求i至少应当比需要辨识的系统的阶数的平方要大；每一个行块含有m行(m为输入数据的个数)，U_0|2i-1有2mi行；j为d-2i+1(d表示输入输出数据的总数)，因为所有得到的数据都要求使用；

(2)输入块Hankel矩阵U_0|2i-1的脚标由该输入块Hankel的第一列中第一个行块u₀和最后一个行块u_2i-1确定，输入块Hankel矩阵U_0|i-1,U_i|2i-1的脚标的定义方式与之相同；输入块可以式中“p”表示“过去”，“f”表示“将来”；那么矩阵U_p表示由过去的输入数据构成的过去输入块Hankel矩阵，矩阵U_f表示由将来的输入数据构成的将来输入块Hankel矩阵；将矩阵U_0|2i-1分为两个矩阵块U_0|i-1,U_i|2i-1，并将U_0|i-1定义为矩阵U_p，将U_i|2i-1定义为矩阵U_f；

输出块Hankel矩阵Y_0|2i-1，Y_p，Y_f定义的方法类似，具体定义形式如下：

其中，矩阵Y_p表示由过去的输出数据构成的输出块Hankel矩阵，矩阵Y_f表示由将来的输出数据构成的输出块Hankel矩阵，输出块Hankel矩阵Y_0|2i-1的脚标由该输出块Hankel的第一列中第一个行块y₀和最后一个行块y_2i-1确定；

进一步，包含有过去输入U_p和过去输出Y_p的过去输入输出Hankel块矩阵W_p定义：

状态序列在子空间辨识算法中起着重要作用；状态序列X_i定义为：

x_i表示状态序列的第一个矩阵块,x_i+j-1表示状态序列的最后一个矩阵块，状态序列为n行j列，由脚标直观的确定矩阵维数，而表示状态序列X_i表示为实数集矩阵，维数为n行j列；

其中下标i表示状态序列的第一个元素的下标为每个矩阵行块数；同时有：

X_p＝X₀,X_f＝X_i (4)

其中，X_p表示过去的实数集构成的矩阵，X₀表示过去的实数集构成的矩阵第一个元素构成的矩阵，X_f表示将来的实数集构成的矩阵，X_i表示将来的实数集构成的矩阵第一个元素构成的矩阵，X_i表示实数集矩阵；

计算将来输出块Hankel矩阵Y_f的行空间沿着将来输入块Hankel矩阵U_f在过去输入输出的Hankel矩阵行空间上的斜投影Ο_i；

公式中：

U_f：定义为将来输入块Hankel矩阵；

Y_f：定义为将来输出块Hankel矩阵，其定义方法与U_f类似；

W_p：包含有过去输入输出的Hankel矩阵；

其中，

正交投影为：∏_B表示将一个矩阵的行空间投影到矩阵的行空间的几何操作如图2：

其中：

表示矩阵的伪逆，矩阵表示矩阵BB^T的伪逆矩阵；

表示将一个矩阵的行空间投影到矩阵的行空间的正交空间的几何操作，并且∏_B和存在式(7)的关系：

其中：I_j：为j阶单位方阵

A/B表示将矩阵的行空间投影到B上的几何操作：

那么

矩阵的行空间沿着的行空间在的行空间上的斜投影定

义为如图3：

其中，first r columns：表示矩阵的前r列；

此外斜投影同样可以按如下方法理解：将A的行空间向B和C的联合行空间正交投影，即：

(二)计算加权斜投影SVD：

对式(5)进行奇异值分解，可分解为：

W₁∈R^li×li：自定义权值矩阵，为满秩矩阵；W₂∈R^j×j：自定义权值矩阵，S为奇异值分解过程中的矩阵，U,V均为酉矩阵，酉矩阵U的列向量被称为Ο_i的左奇异向量，将U从前r列处分块为U＝(U₁ U₂)，酉矩阵V的列向量称为Ο_i的右奇异向量，V的前r列是Ο_i ^HΟ_i的r个非零特征值所对应的特征向量，将他们取为矩阵V₁，则V＝(V₁,V₂)，且满足rank(W_p)＝rank(W_PW₂)；

则

矩阵Ο_i满足：

Ο_i＝Γ_iX_f (13)

其中：Γ_i：扩展可观测矩阵，满秩；可表示为A,C：均为系统矩阵，而且系统矩阵A,C均可观；i：表示每个矩阵行块数；n为系统阶数；l表示输出量的个数；X_f为状态序列；需辨识的线性系统阶数n等于式(12)中非零奇异值的个数，式(12)中奇异值分解过程中的矩阵S的奇异值确定系统阶数n,并获得扩展可观矩阵为：

状态序列X_f中由W₂的列空间决定的部分可以表示为：

其中：A为相似矩阵

状态序列X_f等于：

(三)由矩阵U,S,W确定扩展可观测矩阵Γ_i和扩展可观测矩阵Γ_i行空间的正交补矩阵

公式中：

Γ_i：为扩展观测矩阵；

扩展可观测矩阵Γ_i行空间的正交补矩阵；

(四)由利用输入输出矩阵获得相关扩展可观测矩阵Γ_i和扩展可观测矩阵Γ_i行空间的正交补矩阵确定离散系统矩阵A,C；利用输入输出矩阵获得相关扩展可控矩阵Γ_i和扩展可观测矩阵Γ_i行空间的正交补矩阵获得辨识系统矩阵；

系统矩阵A由扩展可观测矩阵Γ_i的移位后矩阵计算得到：

那么可以得到：

公式中：

表示扩展可观测矩阵Γ_i去掉前l(l表示输出量的个数)行的矩阵；

表示扩展可观测矩阵Γ_i去掉后l(l表示输出量的个数)行的矩阵；

表示矩阵的伪逆，即表示矩阵的伪逆；

Γ_i的前l(l表示输出量的个数)行；

(五)辨识系统矩阵

矩阵输入输出方程定义为如图5：

Y_p＝Γ_iX_p+H_iU_p (22)

Y_f＝Γ_iX_f+H_iU_f (23)

X_f＝AⁱX_p+△_iU_p (24)

其中：

Γ_i：表示系统扩展可观测矩阵，定义为：

H_i：表示三角形式的Toeplitz矩阵，定义为：

△_i：表示系统扩展控矩阵，定义为：

其中，n表示确定系统的阶数；m表示输入的维数；

在式(23)左侧乘以可得：

为将系统矩阵B和D从中分离出来，在式(25)的右侧乘以可得

为将符号简化，将矩阵M代替式(26)中等号左侧部分即：

用矩阵L代替式(26)右侧的扩展可观测矩阵即

其中为将来输入块Hankel矩阵U_f的伪逆；

那么式(26)可表示为：

也可以写成如下形式：

系统矩阵B和D由式(29)(30)计算得到。其它步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：步骤三中假设当前所需辨识系统的系统状态方程，辨识离散状态空间系统矩阵A、B、C和D的过程为：

假设当前所需辨识系统的系统状态方程为：

x_k+1=Ax_k+Bu_k (31)

y_k-Cs_k+Du_k

式(31)线性非时变确定性系统包含输入u_k，输出y_k以及状态序列X_k，由系统状态矩阵A、B、C和D来描述；其中输入u_k可根据测量实际系统的输入得到，或者根据工程需求自行设定输入值，同样可以根据测量实际系统的输出得到y_k，由子空间辨识算法辨识系统矩阵A、B、C和D如图4；

根据可测的输入和输出解决如下问题：

(1)确定系统的阶数n；输入的维数为m行；

(2)确定系统矩阵系统的状态矩阵通过对相关矩阵的行向量空间进行几何操作得到；给定矩阵的行向量作为j维空间向量的坐标，A，B，C各自的行向量定义了一个j维线性向量子空间；j维线性向量子空间的几何操作通过RQ分解实现，通过RQ分解实现的结果，RQ分解就是奇异值分解；进而确定系统矩阵其它步骤及参数与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：步骤四中由于挠性卫星在相对轨道系三轴姿态较为稳定的情况下，挠性卫星保持相对稳定的姿态角，基于挠性卫星姿态角的小角度假设，建立挠性卫星动力学方程和帆板振荡方程的过程为：

选取卫星相对轨道系三轴姿态稳定的情况为辨识环境，基于挠性卫星姿态角的小角度假设建立挠性卫星动力学方程和帆板振荡方程：

公式中：

ω_s：挠性卫星转动角速度；

I_s：整星系统惯性并矢量；

为反对称矩阵，可表示为：

F_iS：帆板振动与星体转动耦合系数；

η_i：模态坐标；

T_e：挠性卫星在轨执行机动任务所受和力矩；

ξ_i：挠性附件的结构阻尼；

Ω_i:挠性附件振动频率。其它步骤及参数与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：步骤五中将建立的挠性卫星动力学方程和帆板振荡方程解耦，通过拉式变换，获得模态参数与力矩到角速度的传递函数的过程为：

基于小角度假设，可将式(32)解耦为三个独立的通道，那么式(32)简化为如下形式：

对式(33)进行拉式变换并进行化简，得到模态参数与力矩到角速度的传递函数的关系；

G(s)ω_s＝T_e (34)

其中：

公式中：

ε_i：挠性附件的结构阻尼；I_s：整星系统惯性并矢量；F_iS：帆板振动与星体转动耦合系数；T_e：挠性卫星在轨执行机动任务所受和力矩；Ω_i:挠性附件振动频率；ε_i为获得模态信息、s表示拉普拉斯变换算子。其它步骤及参数与具体实施方式一至四之一相同。

具体实施方式六：本实施方式与具体实施方式一至五之一不同的是：步骤六中将步骤三得到的辨识离散状态空间系统矩阵A、B、C和D转换为步骤五获得的模态参数与力矩到角速度的传递函数的过程为：

子空间辨识算法主要用于辨识离散状态空间系统矩阵，因此需将离散状态空间矩阵A、B、C和D转换为连续空间的系统矩阵A_c、B_c、C_c和D_c：

将为连续空间的系统矩阵A_c、B_c、C_c和D_c利用式(38)转换为从力矩到角速度的系统传递函数：

G_inv(s)＝C_c(sI-A_c)^-1B_c+D_c (38)

其中，I表示单位矩阵

公式中：

由G(s)函数分母中二次项中即可获得模态信息ε_i和Ω_i；I_sΓ示整星系统惯性并矢量。其它步骤及参数与具体实施方式一至五之一相同。

实施例一：

本实施例一种挠性卫星模态参数在轨辨识方法，具体是按照以下步骤制备的：

下面给出本发明基于子空间辨识算法的挠性卫星在轨模态参数辨识方法分别在惯性空间坐标系和轨道坐标系中的仿真验证与分析；

一、仿真说明；

卫星动力学以及运动学的相关参数选取如下：

帆板展开时整星系统相对Ob坐标系的转动惯量为：

左右帆板振动和卫星本体的耦合系数为：

系统用飞轮做力矩执行机构及力矩激励机构，三轴飞轮配置相同，设力矩幅值为：

T_max＝0.3Nm (43)

采用三轴稳定控制前提下，单轴激励的方法进行在轨辨识实验，其中激励通道为卫星偏航通道；偏航通道的激励力矩采用的是PD控制加随机信号和的形式，并且由simulink中random number模块产生，均值为0，方差为0.1；这样设计的目的是为了在做激励实验的同时，不至于使卫星姿态发生失稳，甚至翻转的恶性现象；PD控制以三轴的姿态欧拉角和欧拉角速度作为输入量，来计算PD控制力矩；随机信号的产生是为了能起到对偏航通道激励充分的目的；随机力矩信号形式可以任取，储存在计算机中，在轨辨识试验期间按时间顺序与PD控制力矩相加即可得出在每个采样周期期间加在偏航通道的力矩；三轴选用的PD控制器参数为：

姿态角和姿态角速度的初值为：

由于子空间辨识算法不依赖于系统初值，因此本实验在进行的初始时刻的星体角速度、姿态角、姿态角速度等均是存在一定的初值；ω_s为挠性卫星转动角速度即姿态角速度

设置仿真时间为400秒到500秒，采样周期为T＝0.2s；

利用采集到的数据，通过本发明算法，实现了离散系统状态矩阵计算、等价连续系统状态矩阵转换、传递函数建立等一系列过程；通过化简后的传递函数分母部分，可以较为准确的辨识出比较关键的前三阶模态信息；

二、相对惯性系姿态控制情况下的在轨辨识仿真分析；

第一组仿真实验为挠性卫星相对惯性系姿态稳定控制；这是因为相对轨道系的姿态控制，当三轴姿态达到稳定时，卫星本体的角速度将和轨道系角速度相同，因此卫星角速度将始终等于一个常值，而不是趋近于0，不能理想达到本发明算法小角速度假设的前提条件，出于保证辨识精度的原因考虑，采用相对惯性系稳定的姿态控制；

在本次仿真实验中，挠性卫星模态信息选取为：

Ω＝diag(0.421071.55711.92252.8435.8136)*2πrad/s (46)

ε＝diag(0.02620.02670.03970.02590.0178) (47)

400秒姿态稳定控制仿真结束后，利用本发明的核心算法，并结合matlab的d2c，ss2tf，tf2zp以及zpk等程序文件，计算得到了从T1到ωs1的传递函数相关的式(48)G(s)，化简后的零极点形式为：

分母部分的几个二次多项式的系数信息包含了模态信息2ε_iΩ_i及经过简单处理，辨识出的前三阶模态信息的结果如下所示：

误差分析：

(1)绝对误差分析：

一阶模态绝对误差：

△Ω₁＝|0.42107-0.4442|×2π＝0.02313×2πrad/s (51)

△ε₁＝|0.0262-0.0307|＝0.0045 (52)

二阶模态绝对误差：

△Ω₂＝|1.5571-1.5987|×2π＝0.0416×2πrad/s (53)

△ε₂＝|0.0267-0.0270|＝0.0003 (54)

三阶模态绝对误差：

△Ω₃＝|1.9225-1.9408|×2π＝0.0183×2πrad/s (55)

△ε₃＝|0.0397-0.0405|＝0.0008 (56)

(2)相对误差：

一阶模态相对误差：

二阶模态相对误差：

三阶模态相对误差：

由上述误差分析，可以看出本发明可以达到较高的辨识精度；

为了进一步证明辨识算法的有效性，将模型动力学部分挠性振荡的阻尼矩阵和频率矩阵修改如下：

Ω＝diag(0.62107 1.3571 1.8225 2.843 5.8136)*2πrad/s (63)

ε＝diag(0.0362 0.0367 0.0497 0.0259 0.0178) (64)

然后再一次用上述实验的条件和流程进行在轨辨识实验；结果为，从T₁到ω_s1的传递函数相关的G(s)化简后的零极点形式：

辨识出的前三阶模态信息的结果如下所示：

误差分析：

(1)绝对误差分析：

一阶模态绝对误差：

△Ω₁＝|0.62107-0.6543|×2π＝0.03323×2πrad/s (68)

△ε₁＝|0.0362-0.0387|＝0.0025 (69)

二阶模态绝对误差：

△Ω₂＝|1.3571-1.3949|×2π＝0.0378×2πrad/s (70)

△ε₂＝|0.0367-0.0371|＝0.0004 (71)

三阶模态绝对误差：

△Ω₃＝|1.8225-1.8389|×2π＝0.0164×2πrad/s (72)

△ε₃＝|0.0497-0.0503|＝0.0006 (73)

(2)相对误差：

一阶模态相对误差：

二阶模态相对误差：

三阶模态相对误差：

由以上误差分析，可以看出，在新的模态参数下，依然可以得到较高的辨识精度；因此证明了所设计的在轨辨识算法的有效性；

三、相对轨道系姿态控制情况下的在轨辨识仿真分析

考虑到目前一部分卫星，如气象卫星或观测卫星等，主要执行的是相对轨道系的姿态控制任务，因此第二组仿真实验给出卫星相对轨道系姿态稳定控制的情况；相对第一组仿真实验的主要修改是加入了轨道角速度，并在此基础上分析对模态参数在轨辨识的精度影响；

在本次仿真实验中，我们选取轨道参数为：

挠性卫星轨道半长轴为：

R_s＝6378120+1204480＝7582600m (80)

地心引力常数为：

M_u＝3.986015×10¹⁴ (81)

轨道周期为：

轨道角速度为：

在本次仿真分析中，卫星模态信息选取为：

Ω＝diag(0.42107 1.5571 1.9225 2.843 5.8136)*2πrad/s (168)

ε＝diag(0.0262 0.0267 0.0397 0.0259 0.0178) (169)

400秒姿态稳定控制仿真结束后利用本发明核心算法，并结合matlab的d2c，ss2tf，tf2zp以及zpk等程序计算得到了式G(s)，化简后的零极点形式：

分母部分的几个二次多项式的系数信息包含了模态信息2ε_iΩ_i及经过简单处理，辨识出的前三阶模态信息的结果如下所示：

误差分析：

(1)绝对误差分析：

一阶模态绝对误差：

△Ω₁＝|0.42107-0.4442|×2π＝0.02313×2πrad/s (173)

△ε₁＝|0.0262-0.0307|＝0.0045 (174)

二阶模态绝对误差：

△Ω₂＝|1.5571-1.5987|×2π＝0.0416×2πrad/s (175)

△ε₂＝|0.0267-0.0270|＝0.0003 (176)

三阶模态绝对误差：

△Ω₃＝|1.9225-1.9408|×2π＝0.0183×2πrad/s (177)

△ε₃＝|0.0397-0.0405|＝0.0008 (178)

(2)相对误差：

一阶模态相对误差：

二阶模态相对误差：

三阶模态相对误差：

可以看出即使引入了轨道角速度，即相对于轨道系姿态稳定下，同样可以达到较高的辨识精度；

为了进一步证明辨识算法的有效性，将模型动力学部分挠性振荡的阻尼矩阵和频率矩阵修改如下：

Ω＝diag(0.62107 1.3571 1.8225 2.843 5.8136)*2πrad/s (185)

ε＝diag(0.0362 0.0367 0.0497 0.0259 0.0178) (186)

然后再一次用上述实验的条件和流程进行在轨辨识实验，G(s)化简后的零极点形式：

辨识出的前三阶模态信息的结果如下所示：

误差分析：

(2)绝对误差分析：

一阶模态绝对误差：

△Ω₁＝|0.62107-0.6543|×2π＝0.03323×2πrad/s (190)

△ε₁＝|0.0362-0.0387|＝0.0025 (191)

二阶模态绝对误差：

△Ω₂＝|1.3571-1.3949|×2π＝0.0378×2πrad/s (192)

△ε₂＝|0.0367-0.0371|＝0.0004 (193)

三阶模态绝对误差：

△Ω₃＝|1.8225-1.8389|×2π＝0.0164×2πrad/s (194)

△ε₃＝|0.0497-0.0503|＝0.0006 (195)

(3)相对误差：

一阶模态相对误差：

二阶模态相对误差：

三阶模态相对误差：

可以看出，在新的模态参数下，依然可以得到较高的辨识精度；因此证明了本发明所设计的在轨辨识算法在考虑轨道角速度下的同样有效性；在卫星姿态控制方面本实施方式在进行在轨辨识时主要在姿态稳定控制基础上施加激励力矩进行，因此有必要介绍卫星姿态控制；获得并保持卫星在空间定向的技术叫做卫星的姿态控制，这种指向一般是指相对于某参考系的姿态，一般参考系主要是惯性系和轨道系；对于在轨的卫星，要求其姿态以给定的要求或规律变化；对于卫星的姿态控制可以分成两类，即被动和主动控制；主动控制则是指利用卫星自身的姿态确定环节、姿态控制器环节、执行机构等环节的联合作用，形成闭环反馈的控制方式；按稳定方式主要分为自旋稳定和三轴稳定；本实施方式主要在三轴稳定的基础上进行在轨辨识研究；在飞轮执行机构卫星的执行机构主要有：喷气、飞轮及磁力矩器；所谓飞轮执行机构就是应用角动量交换将卫星动量偏差转化为飞轮动量控制；国内外大多数卫星主要采用飞轮作为主要执行机构；在精度方面，飞轮为连续旋转工作类型，其控制精度比较高；在使用寿命方面，飞轮由于其能源由太阳帆板提供，理论上可以工作10年；本实施方式中采用飞轮作为姿态控制执行机构以及激励力矩产生机构，控制及激励数据都是基于考虑星载计算机采样周期进行设计的，更符合实际情况；在PID控制方面作为经典的姿态控制方法，PID控制仍不失为一种精确和具有先进性的控制规律；基于自动控制原理进行分析与设计，PID控制器设计时可以考虑系统动态特性和带宽等因素，并通过幅值裕度和相位裕度反映系统的鲁棒性能，因此至今仍为绝大多数的三轴稳定航天器所采用；在PID控制中，比例信号可增加系统通频带，但会使系统稳定性降低；微分信号给系统提供阻尼从而改善稳定性，但也使系统对噪声及干扰较为敏感；积分信号提高系统稳态精度；各控制参数物理意义明确，简单可靠，适当选择后可保证航天器具有较高的控制精度和良好的动态性能；本算法在轨辨识期间的卫星姿态稳定由PID控制保证。

Claims (5)

1.一种挠性卫星模态参数在轨辨识方法，其特征在于：一种挠性卫星模态参数在轨辨识方法具体是按照以下步骤进行的：

步骤一、在外界诸多干扰不影响挠性卫星姿态稳定的前提下，于每个采样周期T平均准确的采集到执行机构施加到挠性卫星体上的力矩和挠性卫星体相对惯性坐标系的角速度信息；

步骤二、将采集到的力矩和角速度信息代入子空间辨识算法，确定离散状态空间系统矩阵A和C并且计算出离散状态空间系统矩阵B和D；

步骤三、假设当前所需辨识系统的系统状态方程，确定离散状态空间系统矩阵A、B、C和D；

步骤四、基于挠性卫星姿态角的小角度假设，建立挠性卫星动力学方程和帆板振荡方程；

步骤五、将建立的挠性卫星动力学方程和帆板振荡方程解耦，通过拉式变换，获得模态参数与力矩到角速度的传递函数；

步骤六、将步骤三得到的离散状态空间系统矩阵A、B、C和D转换为步骤五获得的模态参数与力矩到角速度的传递函数；即完成了一种挠性卫星模态参数在轨辨识方法；

步骤二中将采集到的力矩和角速度信息代入子空间辨识算法，确定离散状态空间系统矩阵A和C并且计算出离散状态空间系统矩阵B和D的实现步骤为：

(一)通过输入输出数据构造Hankel块矩阵；其中输入块Hankel矩阵由输入数据构造，具体构成形式：

其中，

(1)i表示行块数，每一个行块含有m行，U_0|2i-1有2mi行；j为d-2i+1；

(2)输入块Hankel矩阵U_0|2i-1的脚标由该输入块Hankel的第一列中第一个行块u₀和最后一个行块u_2i-1确定，输入块式中“p”表示“过去”，“f”表示“将来”；那么矩阵U_p表示由过去的输入数据构成的过去输入块Hankel矩阵，矩阵U_f表示由将来的输入数据构成的将来输入块Hankel矩阵；将矩阵U_0|2i-1分为两个矩阵块U_0|i-1,U_i|2i-1，并将U_0|i-1定义为矩阵U_p，将U_i|2i-1定义为矩阵U_f；

输出块Hankel矩阵Y_0|2i-1，Y_p，Y_f具体定义形式如下：

其中，矩阵Y_p表示由过去的输出数据构成的输出块Hankel矩阵，矩阵Y_f表示由将来的输出数据构成的输出块Hankel矩阵，输出块Hankel矩阵Y_0|2i-1由该输出块Hankel的第一列中第一个行块y₀和最后一个行块y_2i-1确定；

包含过去输入U_p和过去输出Y_p的过去输入输出Hankel块矩阵W_p定义：

$W_p=\left(\begin{array}{c}{U}_P\\ Y_P\end{array}\right)=\left(\begin{array}{c}{U}_{0|i-1}\\ Y_{0|i-1}\end{array}\right)---\left(2\right)$

状态序列X_i定义为：

其中，x_i表示状态序列的第一个矩阵块,x_i+j-1表示状态序列的最后一个矩阵块，状态序列为n行j列，表示状态序列，状态序列X_i表示维数为n行j列实数集矩阵；

其中i表示状态序列的第一个元素的下标为每个矩阵行块数；同时有：

X_p＝X₀,X_f＝X_i (4)

其中，X_p表示过去的实数集构成的矩阵，X₀表示过去的实数集构成的矩阵第一个元素构成的矩阵，X_f表示将来的实数集构成的矩阵；

计算将来输出块Hankel矩阵Y_f的行空间沿着将来输入块Hankel矩阵U_f包含过去输入输出的Hankel矩阵行空间上的斜投影O_i；

$O_i=Y_f{/}_{U_f}{W}_p---\left(5\right)$

公式中：

U_f：定义为将来输入块Hankel矩阵；

Y_f：定义为将来输出块Hankel矩阵；

W_p：包含有过去输入输出的Hankel矩阵；

其中，

正交投影为：∏_B表示将一个矩阵的行空间投影到离散状态空间系统矩阵的行空间的几何操作：

其中：

表示矩阵的伪逆，矩阵表示矩阵BB^T的伪逆矩阵；

表示将一个矩阵的行空间投影到矩阵的行空间的正交空间的几何操作，并且∏_B和存在式(7)的关系：

${\mathrm{\Π}}_{B^{\⊥}}=I_j-{\mathrm{\Π}}_B---\left(7\right)$

其中：I_j：为j阶单位方阵

A/B表示将离散状态空间系统矩阵的行空间投影到B上的几何操作：

那么

$A/B^{\⊥}={\mathrm{A\Π}}_{B^{\⊥}}---\left(9\right)$

矩阵的行空间沿着的行空间在的行空间上的斜投影定义为：

$A{/}_{B}C=A\left(\begin{array}{cc}{C}^T& B^T\end{array}\right){\left(\begin{array}{cc}{\mathrm{CC}}^T& {\mathrm{CB}}^T\\ {\mathrm{BC}}^T& {\mathrm{BB}}^T\end{array}\right)}_{firstrcolumns}C---\left(10\right)$

其中，first r columns：表示矩阵的前r列；

将A的行空间向B和C的联合行空间正交投影，即：

$A{/}_{B}C=A/\left(\begin{array}{c}C\\ B\end{array}\right)---\left(11\right)$

(二)计算加权斜投影SVD：

对式(5)进行奇异值分解，分解为：

$\begin{array}{c}{W}_1{O}_i{W}_2={\mathrm{USV}}^T=\left(\begin{array}{cc}{U}_1& U_2\end{array}\right)\left(\begin{array}{cc}{S}_1& 0\\ 0& 0\end{array}\right)\left(\begin{array}{c}{V}_1^T\\ V_2^T\end{array}\right)\\ =U_1{S}_1{V}_1^T\end{array}---\left(12\right)$

W₁∈R^li×li：自定义权值矩阵，为满秩矩阵；W₂∈R^j×j：自定义权值矩阵，S为奇异值分解过程中的矩阵，U,V均为酉矩阵，酉矩阵U的列向量被称为O_i的左奇异向量，将U从前r列处分块为U＝(U₁ U₂)，酉矩阵V的列向量称为O_i的右奇异向量，V的前r列是O_i ^HO_i的r个非零特征值所对应的特征向量，则V＝(V₁,V₂)，且满足：

rank(W_p)＝rank(W_P W₂)；

则矩阵O_i满足：

O_i＝Γ_iX_f (13)

其中，Γ_i：扩展可观测矩阵，满秩；表示为A,C：均为离散状态空间系统矩阵，而且离散状态空间系统矩阵A,C均可观，i：表示每个矩阵行块数，n为系统阶数，l表示输出量的个数，X_f为状态序列，辨识的线性系统阶数n等于式(12)中非零奇异值的个数；式(12)中奇异值分解过程中的矩阵S的奇异值确定系统阶数n，并获得扩展可观矩阵为：

状态序列X_f中由W₂的列空间决定的部分表示为：

$X_f{W}_2=A^{-1}{S}_1^{1/2}{V}_1^T---\left(15\right)$

其中，A为离散状态空间系统矩阵；

状态序列X_f等于：

(三)由矩阵U,S,W确定扩展可观测矩阵Γ_i和扩展可观测矩阵Γ_i行空间的正交补矩阵

${\mathrm{\Γ}}_i=W_1^{-1}{U}_1{S}_1^{1/2}---\left(17\right)$

${\mathrm{\Γ}}_i^{\⊥}=U_2^T{W}_1---\left(18\right)$

其中，

Γ_i：为扩展观测矩阵；

扩展可观测矩阵Γ_i行空间的正交补矩阵；

(四)由利用输入输出矩阵获得相关扩展可观测矩阵Γ_i和扩展可观测矩阵Γ_i行空间的正交补矩阵确定离散状态空间系统矩阵A，C；

离散状态空间系统矩阵A由扩展可观测矩阵Γ_i的移位后矩阵计算得到：

$\underset{\‾}{{\mathrm{\Γ}}_i}A=\stackrel{\‾}{{\mathrm{\Γ}}_i}---\left(19\right)$

那么得到：

$C={\mathrm{\Γ}}_{i_{firstlrow}}---\left(21\right)$

其中，

表示扩展可观测矩阵Γ_i去掉前l行的矩阵；

Γ_i ：表示扩展可观测矩阵Γ_i去掉后l行的矩阵；

表示矩阵的伪逆，即表示矩阵Γ_i 的伪逆；

Γ_i的前l行；

(五)辨识系统矩阵

矩阵输入输出方程定义为：

Y_p＝Γ_iX_p+H_iU_p (22)

Y_f＝Γ_iX_f+H_iU_f (23)

X_f＝AⁱX_p+Δ_iU_p (24)

其中，

Γ_i：表示系统扩展可观测矩阵，定义为：

H_i：表示三角形式的Toeplitz矩阵，定义为：

Δ_i：表示系统扩展控矩阵，定义为：

其中，n表示确定系统的阶数；m表示输入的维数；

公式(23)左侧乘以得到公式(25)：

${\mathrm{\Γ}}_i^{\⊥}{Y}_f={\mathrm{\Γ}}_i^{\⊥}{H}_i.U_f---\left(25\right)$

在式(25)的右侧乘以得到公式(26)：

将矩阵M代替式(26)中即：

用矩阵L代替式(26)右侧的扩展可观测矩阵即

$L=\left(\begin{array}{cccc}{L}_1& L_2& ...& L_i\end{array}\right)={\mathrm{\Γ}}_i^{\⊥}---\left(28\right)$

其中，为将来输入块Hankel矩阵U_f的伪逆；

那么式(26)表示为：

$\begin{array}{c}\left(\begin{array}{cccc}{M}_1& M_2& ...& M_i\end{array}\right)=\left(\begin{array}{cccc}{L}_1& L_2& ...& L_i\end{array}\right)\\ \×\left(\begin{array}{ccccc}D& 0& 0& ...& 0\\ CB& D& 0& ...& 0\\ CAB& CB& D& ...& 0\\ ...& ...& ...& ...& ...\\ {\mathrm{CA}}^{i-2}B& {\mathrm{CA}}^{i-3}B& {\mathrm{CA}}^{i-4}B& ...& D\end{array}\right)\end{array}---\left(29\right)$

写成如下形式：

离散状态空间系统矩阵B和D由式(29)(30)计算得到。

2.根据权利要求1所述一种挠性卫星模态参数在轨辨识方法，其特征在于：步骤三中假设当前所需辨识系统的系统状态方程，辨识离散状态空间系统矩阵A、B、C和D的过程为：

假设当前所需辨识系统的系统状态方程为：

$\begin{array}{c}{x}_{k+1}={\mathrm{Ax}}_k+{\mathrm{Bu}}_k\\ y_k={\mathrm{Cx}}_k+{\mathrm{Du}}_k\end{array}---\left(31\right)$

式(31)线性非时变确定性系统包含输入u_k，输出y_k以及状态序列X_k，其中输入u_k根据测量实际系统的输入得到，根据测量实际系统的输出得到y_k，由子空间辨识算法辨识离散状态空间系统矩阵A、B、C和D。

3.根据权利要求1所述一种挠性卫星模态参数在轨辨识方法，其特征在于：步骤四中基于挠性卫星姿态角的小角度假设，建立挠性卫星动力学方程和帆板振荡方程为：

选取卫星相对轨道系三轴姿态稳定的情况为辨识环境，基于挠性卫星姿态角的小角度假设建立挠性卫星动力学方程和帆板振荡方程：

$\left\{\begin{array}{c}{I}_s{\stackrel{\·}{\mathrm{\ω}}}_s+{{\mathrm{\ω}}_s}^{\×}{I}_s{\mathrm{\ω}}_s+\underset{i}{\mathrm{\Σ}}\[F_{iS}{\stackrel{\·\·}{\mathrm{\η}}}_i+{{\mathrm{\ω}}_s}^{\×}{F}_{iS}{\stackrel{\·}{\mathrm{\η}}}_i\]=T_e\\ {\stackrel{\·\·}{\mathrm{\η}}}_i+2{\mathrm{\ξ}}_i{\mathrm{\Ω}}_i{\stackrel{\·}{\mathrm{\η}}}_i+{\mathrm{\Ω}}_i^2{\mathrm{\η}}_i+F_{iS}^T{\stackrel{\·}{\mathrm{\ω}}}_s=0\end{array}\right.---\left(32\right)$

公式中：

ω_s：挠性卫星转动角速度；

I_s：整星系统惯性并矢量；

ω_s ^×：为反对称矩阵，表示为：

F_iS：帆板振动与星体转动耦合系数；

η_i：模态坐标；

T_e：挠性卫星在轨执行机动任务所受和力矩；

ξ_i：挠性附件的结构阻尼；

Ω_i:挠性附件振动频率。

4.根据权利要求3所述一种挠性卫星模态参数在轨辨识方法，其特征在于：步骤五中将建立的挠性卫星动力学方程和帆板振荡方程解耦，通过拉式变换，获得模态参数与力矩到角速度的传递函数的过程为：

基于小角度假设，将式(32)解耦为三个独立的通道，那么式(32)简化为如下形式：

$\left\{\begin{array}{c}{I}_s{\stackrel{\·}{\mathrm{\ω}}}_s+\underset{i}{\mathrm{\Σ}}{F}_{iS}{\stackrel{\·\·}{\mathrm{\η}}}_i=T_e\\ {\stackrel{\·\·}{\mathrm{\η}}}_i+2{\mathrm{\ξ}}_i{\mathrm{\Ω}}_i{\stackrel{\·}{\mathrm{\η}}}_i+{\mathrm{\Ω}}_i^2{\mathrm{\η}}_i+F_{iS}^T{\stackrel{\·}{\mathrm{\ω}}}_s=0\end{array}\right.---\left(33\right)$

对式(33)进行拉式变换并进行化简，得到模态参数与力矩到角速度的传递函数的关系；

G(s)ω_s＝T_e (34)

其中：

$G\left(s\right)=I_{s}s+\underset{i}{\Σ}{F}_{is}{s}^2\frac{-F_{is}^{T}s}{s^2+2{\mathrm{\ϵ}}_i{\mathrm{\Ω}}_{i}s+{\mathrm{\Ω}}_i^2}---\left(35\right)$

公式中：

ξ_i：挠性附件的结构阻尼；I_s：整星系统惯性并矢量；F_iS：帆板振动与星体转动耦合系数；T_e：挠性卫星在轨执行机动任务所受和力矩；Ω_i:挠性附件振动频率；ε_i为获得模态信息、s表示拉普拉斯变换算子。

5.根据权利要求4所述一种挠性卫星模态参数在轨辨识方法，其特征在于：步骤六中将步骤三得到的离散状态空间系统矩阵A、B、C和D转换为步骤五获得的模态参数与力矩到角速度的传递函数的过程为：

将离散状态空间系统矩阵A、B、C和D转换为连续空间的系统矩阵A_c、B_c、C_c和D_c：

$A=e^{A_{c}T}---\left(36\right)$

$B={\∫}_0^T{e}^{A_{c}t}{\mathrm{dtB}}_c---\left(37\right)$

将为连续空间的系统矩阵A_c、B_c、C_c和D_c利用式(38)转换为从力矩到角速度的系统传递函数：

G_inv(s)＝C_c(sI-A_c)^-1B_c+D_c (38)

其中，I表示单位矩阵

公式中：

$G\left(s\right)=\frac{1}{G_{inv}\left(s\right)}=I_{s}s+\underset{i}{\Σ}{F}_{is}{s}^2\frac{-F_{is}^{T}s}{s^2+2{\mathrm{\ϵ}}_i{\mathrm{\Ω}}_{i}s+{\mathrm{\Ω}}_i^2}---\left(39\right)$

由G(s)函数分母中二次项中获得模态信息ε_i和Ω_i；I_s表示整星系统惯性并矢量。