一种高分辨率小卫星颤振响应的确定方法.pdf

一种高分辨率小卫星颤振响应的确定方法，本发明采用单框控制力矩陀螺转子动不平衡产生的干扰力矩的确定方法，解决了此类高分辨率小卫星的颤振干扰确定问题；本发明通过适当的状态变量设置和替换，实现卫星颤振动力学方程自由度间的解耦和降阶，将其转化为一阶微分方程，利用数值解法，确定卫星在颤振干扰源作用下产生的姿态响应，方法简明规范；本发明对于卫星在颤振干扰作用下引起的姿态响应，结合卫星姿态控制规律消除颤振姿态响应中的可控部分，最终确定卫星在颤振干扰源作用下的颤振响应；本发明的确定方法还可以通过调节陀螺的安装方式和陀螺的加工精度来调整控制力矩陀螺的干扰力矩，最终达到提高图像质量的目的。

1.  一种高分辨率小卫星颤振响应的确定方法，其特征在于通过以下步骤实现：
第一步，根据公式(1)建立高分辨率小卫星的颤振动力学方程，
Ibω·b+Σk=1KFskτ··k=TdFskTω·b+τ··k+2μΩkτ·k+Ωk2τk=0(k=1,2,···,K)---(1)]]>
其中，I_b为星体转动惯量；ω_b为星体转动角速度；F_sk为第k个太阳翼附件的柔性对卫星转动的耦合系数，τ_k＝[τ_k1 τ_k2…τ_k6]^T，为第k个太阳翼附件前6阶模态坐标列阵，μ为太阳翼附件的阻尼系数，Ω_k＝diag(ω_k1ω_k2…ω_k6)，为第k个太阳翼附件前6阶模态频率列阵，ω_kj(j＝1，2，…，6)为第k个太阳翼第j阶模态频率，K为太阳翼附件的总数，T_d为卫星控制系统执行机构单框架控制力矩陀螺群动不平衡对卫星本体产生的干扰力矩；
第二步，确定卫星单框架控制力矩陀螺群动不平衡对卫星本体产生的干扰力矩T_d；
第三步，根据第二步确定的干扰力矩T_d，利用第一步中建立的颤振动力学方程，确定高分辨率小卫星在颤振干扰力矩作用下的姿态响应，
(3.1)将公式(1)的颤振动力学方程改写为公式(9)的形式，
Mx··+Cx·+Kx=Q---(9)]]>
其中，

广义坐标变量其中θ、ψ为卫星中心体的滚动角、俯仰角和偏航角；
(3.2)将公式(9)的二阶常微分方程转化为公式(10)的一阶状态方程形式：
y·=Ay+P---(10)]]>
其中，y=xx·,]]>A=0n×nIn×n-M-1K-M-1C,]]>P=0n×1M-1Q;]]>n是卫星系统的颤振动力学方程组自由度数，I_n×n是n×n的单位阵；
(3.3)利用一阶常微分方程初值问题的求解算法解算公式(10)得到卫星在颤振干扰力矩作用下的滚动、俯仰和偏航姿态响应，即滚动角、俯仰角和偏航角；
第四步，利用卫星姿态控制规律确定第一步建立的卫星柔性动力学方程中的可控姿态响应；
第五步，将第三步中确定的卫星在颤振干扰力矩作用下的姿态响应，减去第四步确定的可控姿态响应，获得卫星在颤振干扰力矩作用下引起的颤振响应。

2.  根据权利要求1所述的一种高分辨率小卫星颤振响应的确定方法，其特征在于：所述第二步中干扰力矩T_d根据以下步骤确定，
(2.1)利用转子质量特性测试方法，获得度量陀螺转子的主惯性轴相对旋转轴线倾斜程度的量J_xz、J_yz；
(2.2)根据公式组(2)得到陀螺转子的动不平衡度值和初始相位角，
I0=Jxz2+Jyz2cosφ=-Jyz/I0sinφ=Jxz/I0---(2)]]>
其中，I₀为陀螺转子的动不平衡度值，φ为陀螺转子转不平衡度的初始相位角；
(2.3)卫星控制系统执行机构测量部分测得每个陀螺转子相对自身框架的转角即转子转角γ_i和框架相对陀螺基座的转角即框架转角ζ_i，以及转子绕旋转轴旋转的角速度i＝1，2…N，N是陀螺的总数；
(2.4)利用步骤(2.2)得到I₀和φ和步骤(2.3)得到的角速度根据公式(3)得到陀螺转子在转子坐标系ox_my_mz_m中干扰力矩，

其中T_dm为陀螺转子在转子坐标系ox_my_mz_m下产生的干扰力矩，为转子坐标系ox_my_mz_mx_m向、y_m向的单位矢量；
(2.5)利用步骤(2.3)得到的转子转角γ_i和框架转角ζ_i，根据公式(4)、(5)得到框架坐标系ox_ry_rz_r相对陀螺基座坐标系ox_sy_sz_s的转换矩阵A_rsi和转子坐标系ox_my_mz_m相对框架坐标系ox_ry_rz_r的转换矩阵A_mri，
Arsi=1000cosζisinζi0-sinζicosζi---(4)]]>
Amri=cosγisinγi0-sinγicosγi0001---(5)]]>
A_rsi是第i个陀螺的框架坐标系相对陀螺基座坐标系的转换矩阵，A_mri是第i个陀螺的转子坐标系相对框架坐标系的转换矩阵；
(2.6)通过公式(6)得到由N个陀螺组成的正锥形陀螺群的陀螺基座坐标系ox_sy_sz_s到星体坐标系o_bx_by_bz_b的转换矩阵，陀螺群中的每一个陀螺的框架轴与卫星星体坐标系o_bx_by_bz_b的z_b轴共面，并与z_b轴的夹角为β，框架轴在卫星星体坐标平面x_by_b的投影与x_b轴夹角α_i，
Msbi=sinβcosαisinβsinαicosβ-cosβcosαi-cosβsinαisinβsinαi-cosαi0---(6)]]>
其中，M_sbi是第i个陀螺的陀螺基座坐标系到星体坐标系的转换矩阵；
(2.7)利用步骤(2.5)中得到的转换矩阵A_rsi、A_mri和步骤(2.6)得到的转换矩阵M_sbi，将公式(3)经过三次坐标转换，从转子坐标系转化到卫星星体坐标系，得到每i个陀螺转子在卫星星体坐标系内的干扰力矩T_di，
Tdi=MsbiT·ArsiT·AmriT·Tdm---(7);]]>
(2.8)根据公式(8)得到陀螺群在卫星星体坐标系内的干扰力矩T_d，Td=Σi=1NTdi]]>

一种高分辨率小卫星颤振响应的确定方法
技术领域
本发明涉及一种高分辨率小卫星颤振响应的确定方法，特别是采用控制力矩陀螺执行机构的高分辨率小卫星的颤振响应确定方法，属于卫星动力学技术领域，可用于卫星姿态干扰、颤振分析。
背景技术
卫星成像期间，活动部件运动、卫星结构的柔性等会使卫星主体以及CCD相机产生扰动，其中幅值较小、频率较高、控制系统无法测控的部分为颤振。高分辨率小卫星与大卫星相比，惯性质量轻，更易受到干扰影响。
2003年于登云、王跃宇等在“遥感卫星颤振响应分析技术研究”和“遥感卫星颤振响应分析的实用方法”文中，就颤振动力学分析技术进行了探讨，给出了基于拉格朗日法的遥感卫星动力学模型及其扰动源扰动力矩模型，给出了考虑太阳阵和相机柔性时在太阳阵步进、磁带机和动量轮运动综合干扰下卫星星体颤振响应的数值分析结果。
2008年NASA戈达德飞行中心Kuo-Chia Liu和Peiman Maghami等人在“Reaction Wheel Disturbance Modeling，Jitter Analysis，and Validation Tests forSolar Dynamics Observatory”一文中介绍了太阳动力观测计划SDO针对在轨主要干扰源反作用轮所作的干扰建模和颤振分析，详细介绍了干扰模型的验证修正和反作用颤振分析的地面验证试验以及为了满足观测仪器的颤振要求对反作用轮转速限速的措施等。以上研究都主要针对反作用轮干扰等进行卫星颤振分析，不能满足高分辨率小卫星颤振确定的需要。
本发明针对目前的技术空白，首次提出了一种利用单框控制力矩陀螺动不平衡产生的干扰力矩来确定卫星的颤振响应，利用本发明确定的干扰力矩，可以对敏捷小卫星进行姿态扰动分析、颤振分析，从而评判卫星成像质量，测试卫星的分辨率，还可以通过调整控制力矩陀螺的干扰力矩最终使卫星的分辨率达到设计要求。
高分辨率小卫星的敏捷机动性日益得到重视和发展，此类高分辨率小卫星主要采用控制力矩陀螺作为执行机构，目前尚未见到对于此类高分辨率小卫星的主要干扰源控制力矩陀螺的干扰建模和颤振分析。
发明内容
本发明的技术解决问题是：克服现有技术的不足，提供一种高分辨率小卫星颤振响应的确定方法，填补现有技术空白，本发明方法为卫星成像质量的影响分析及最终改进卫星成像质量提供依据。
本发明的技术解决方案是：一种高分辨率小卫星颤振响应的确定方法，通过以下步骤实现：
第一步，根据公式(1)建立高分辨率小卫星的颤振动力学方程，

其中，I_b为星体转动惯量；ω_b为星体转动角速度；F_sk为第k个太阳翼附件的柔性对卫星转动的耦合系数，τ_k＝[τ_k1 τ_k2…τ_k6]^T，为第k个太阳翼附件前6阶模态坐标列阵，μ为太阳翼附件的阻尼系数，Ω_k＝diag(ω_k1 ω_k2…ω_k6)，为第k个太阳翼附件前6阶模态频率列阵，ω_kj(j＝1，2，…，6)为第k个太阳翼第j阶模态频率，K为太阳翼附件的总数，T_d为卫星控制系统执行机构单框架控制力矩陀螺群动不平衡对卫星本体产生的干扰力矩；
第二步，确定卫星单框架控制力矩陀螺群动不平衡对卫星本体产生的干扰力矩T_d；
第三步，根据第二步确定的干扰力矩T_d，利用第一步中建立的颤振动力学方程，确定高分辨率小卫星在颤振干扰力矩作用下的姿态响应，
(3.1)将公式(1)的颤振动力学方程改写为公式(9)的形式，
Mx··+Cx·+Kx=Q---(9)]]>
其中，

广义坐标变量，其中θ、ψ为卫星中心体的滚动角、俯仰角和偏航角；
(3.2)将公式(9)的二阶常微分方程转化为公式(10)的一阶状态方程形式：

其中，n是卫星系统的颤振动力学方程组自由度数，I_n×n是n×n的单位阵；
(3.3)利用一阶常微分方程初值问题的求解算法解算公式(10)得到卫星在颤振干扰力矩作用下的滚动、俯仰和偏航姿态响应，即滚动角、俯仰角和偏航角；
第四步，利用卫星姿态控制规律确定第一步建立的卫星柔性动力学方程中的可控姿态响应；
第五步，将第三步中确定的卫星在颤振干扰力矩作用下的姿态响应，减去第四步确定的可控姿态响应，获得卫星在颤振干扰力矩作用下引起的颤振响应。
所述第二步中干扰力矩T_d根据以下步骤确定，
(2.1)利用转子质量特性测试方法，获得度量陀螺转子的主惯性轴相对旋转轴线倾斜程度的量J_xz、J_yz；
(2.2)根据公式组(2)得到陀螺转子的动不平衡度值和初始相位角，

其中，I₀为陀螺转子的动不平衡度值，φ为陀螺转子转不平衡度的初始相位角；
(2.3)卫星控制系统执行机构测量部分测得每个陀螺转子相对自身框架的转角即转子转角γ_i和框架相对陀螺基座的转角即框架转角ζ_i，以及转子绕旋转轴旋转的角速度i＝1，2…N，N是陀螺的总数；
(2.4)利用步骤(2.2)得到I₀和φ和步骤(2.3)得到的角速度，根据公式(3)得到陀螺转子在转子坐标系ox_my_mz_m中干扰力矩，

其中T_dm为陀螺转子在转子坐标系ox_my_mz_m下产生的干扰力矩，为转子坐标系ox_my_mz_mx_m向、y_m向的单位矢量；
(2.5)利用步骤(2.3)得到的转子转角γ_i和框架转角ζ_i，根据公式(4)、(5)得到框架坐标系ox_ry_rz_r相对陀螺基座坐标系ox_sy_sz_s的转换矩阵A_rsi和转子坐标系ox_my_mz_m相对框架坐标系ox_ry_rz_r的转换矩阵A_mri，


A_rsi是第i个陀螺的框架坐标系相对陀螺基座坐标系的转换矩阵，A_mri是第i个陀螺的转子坐标系相对框架坐标系的转换矩阵；
(2.6)通过公式(6)得到由N个陀螺组成的正锥形陀螺群的陀螺基座坐标系ox_sy_sz_s到星体坐标系o_bx_by_bz_b的转换矩阵，陀螺群中的每一个陀螺的框架轴与卫星星体坐标系o_bx_by_bz_b的z_b轴共面，并与z_b轴的夹角为β，框架轴在卫星星体坐标平面x_by_b的投影与x_b轴夹角α_i，

其中，M_sbi是第i个陀螺的陀螺基座坐标系到星体坐标系的转换矩阵；
(2.7)利用步骤(2.5)中得到的转换矩阵A_rsi、A_mri和步骤(2.6)得到的转换矩阵M_sbi，将公式(3)经过三次坐标转换，从转子坐标系转化到卫星星体坐标系，得到每i个陀螺转子在卫星星体坐标系内的干扰力矩T_di，
Tdi=MsbiT·ArsiT·AmriT·Tdm---(7);]]>
(2.8)根据公式(8)得到陀螺群在卫星星体坐标系内的干扰力矩T_d，
Td=Σi=1NTdi]]>

本发明设计原理：
对于高分辨率小卫星，考虑太阳翼的柔性，采用柔性航天器动力学建模方法，建立卫星颤振动力学方程，然后建立主要干扰源单框架控制力矩陀螺动不平衡产生的干扰力矩模型，进行卫星姿态响应分析，最后结合姿态控制结果确定卫星的颤振响应。包括下列步骤：
1、建立卫星颤振动力学方程。根据柔性航天器动力学建模方法，结合小卫星部件的刚度特点，主要考虑太阳翼的柔性，建立小卫星柔性航天器颤振动力学方程；
2、建立主要干扰源单框架控制力矩陀螺产生的干扰力矩数学模型。由定轴转动部件的动不平衡对其转轴产生的干扰力矩出发，结合单框架控制力矩陀螺转子的两自由度转动，即转子相对于框架的转动和框架相对于陀螺基座的转动，通过连续的坐标变换获得单框架控制力矩陀螺转子动不平衡对卫星产生的干扰力矩数学模型。
(1)建立分析中用到的坐标系
定义转子坐标系ox_my_mz_m，为转子本体坐标系，与转子固连；
定义框架坐标系ox_ry_rz_r，坐标轴分别沿框架轴、力矩轴和转子角动量轴，与框架固连，初始时刻ox_ry_rz_r和ox_my_mz_m重合；
定义陀螺基座坐标系ox_sy_sz_s，当框架角为零时，它与框架坐标系ox_ry_rz_r重合。根据陀螺群的构形和每个陀螺框架轴的方向布置，可确定每个陀螺基座的安装矩阵M。各坐标系如图5所示。
(2)各坐标系之间的转换
控制力矩陀螺转子相对框架的转角记为γ，框架相对陀螺基座的转角记为ζ。陀螺基座坐标系到转子坐标系的变换可由两次旋转得到，即先绕x_r轴转过ζ，再绕z_m轴转过γ。
框架坐标系ox_ry_rz_r相对陀螺基座坐标系ox_sy_sz_s的转换矩阵为

转子坐标系ox_my_mz_m相对框架坐标系ox_ry_rz_r的转换矩阵为

对于锥形构形的单框控制力矩陀螺群，每只陀螺的框架轴与星体某一体轴(例如z轴)共面，与该轴夹角β，框架轴在星体坐标平面x_by_b的投影与x_b轴夹角α，框架角零位的角动量平行于x_by_b平面，可确定每个陀螺基座坐标系到星体坐标系得转换矩阵，即陀螺基座的安装矩阵为

(3)确定转子动不平衡在转子坐标系下产生的干扰力矩根据达朗贝尔原理，定轴转动部件旋转时对转动轴产生的惯性力矩为：

特别是当等速旋转时

Ud=k→×(Jxzi→+Jyzj→)]]>
其中为转动部件绕定轴旋转的角速度，U_d为动不平衡，由转子主惯性轴和旋转轴线不重合引起，J_xz，J_yz是度量转子的主惯性轴相对旋转轴线倾斜程度的量，由转子质量特性测试获得。
对于单框架控制力矩陀螺，其转子绕角动量轴作匀速转动，因此转子动不平衡产生的干扰力矩T_dm在转子坐标系o_mx_my_mz_m中为

令I₀为动不平衡大小度量，则上式可以表示为

(4)确定转子动不平衡在星体坐标系下产生的干扰力矩
根据转子动不平衡在转子坐标系下产生的干扰力矩表达式，通过连续的坐标变换将其转换到星体坐标系下，最终确定转子动不平衡在星体坐标系下产生的干扰力矩。
在星体o_bx_by_bz_b坐标系中，T_d为
Td=MsbT·ArsT·AmrT·Tdm]]>
将上文中的各转换矩阵代入，得转子动不平衡对卫星本体产生的干扰力矩为

3、确定柔性航天器在颤振干扰力矩作用下的姿态响应。颤振动力学方程为二阶常微分方程组，方程自由度数高，各自由度间耦合严重。为便于分析，首先定义系统广义坐标x＝[x₁ x₂…x_n]^T(n为系统广义自由度数)，将颤振动力学方程改写成
Mx··+Cx·+Kx=Q]]>
其中M为航天器系统的广义质量矩阵，C为系统的广义阻尼矩阵，K为系统的广义刚度矩阵，Q为系统所受的广义外力。
设状态变量y=xx·]]>，通过状态变量替换，将系统的颤振动力学方程转化为一阶常微分方程
y·=Ay+P]]>
其中

从而可利用一阶常微分方程初值问题的求解算法获得卫星在颤振干扰力矩作用下的姿态响应；
4、确定航天器的颤振响应。对于卫星在颤振干扰力矩作用下的姿态响应，其中一部分为可控姿态，按照1建立的小卫星柔性动力学方程，根据卫星姿态控制规律确定可以控制的姿态。将3确定的卫星在颤振干扰力矩作用下的姿态响应，减去卫星姿态控制系统可控的姿态部分，最终获得卫星在颤振干扰力矩作用下引起的颤振响应。卫星的颤振响应可用于成像质量分析评估，也为卫星成像质量改进提供了技术基础。
本发明与现有技术相比有益效果为：
(1)本发明采用单框控制力矩陀螺转子动不平衡产生的干扰力矩的确定方法，解决了此类高分辨率小卫星的颤振干扰确定问题；
(2)本发明通过适当的状态变量设置和替换，实现卫星颤振动力学方程自由度间的解耦和降阶，将其转化为一阶微分方程，利用数值解法，确定卫星在颤振干扰源作用下产生的姿态响应，方法简明规范；
(3)本发明对于卫星在颤振干扰作用下引起的姿态响应，结合卫星姿态控制规律消除颤振姿态响应中的可控部分，最终确定卫星在颤振干扰源作用下的颤振响应；
(4)本发明的确定方法还可以通过调节陀螺的安装方式和陀螺的加工精度来调整控制力矩陀螺的干扰力矩，最终达到提高图像质量的目的。
附图说明
图1为本发明单框架控制力矩陀螺坐标系示意图；
图2为本发明实施方式中采用的正五面锥形控制力矩陀螺构形示意图；
图3为利用本发明后某敏捷高分辨率小卫星本体滚动颤振角；
图4为利用本发明后某敏捷高分辨率小卫星本体俯仰颤振角；
图5为利用本发明后某敏捷高分辨率小卫星本体偏航颤振角。
具体实施方式
对于高分辨率小卫星，主要考虑太阳翼的柔性，采用柔性航天器动力学建模方法，建立卫星颤振动力学方程，一方面通过建立主要颤振干扰源单框架控制力矩陀螺产生的干扰力矩模型，进行卫星姿态响应分析，另一方面根据卫星姿态控制系统控制规律获取卫星的可控姿态响应，最后从卫星在颤振干扰源作用下的姿态响应中消除可以控制的姿态部分确定卫星的颤振响应，为卫星成像质量评估提供依据。
下面根据具体实例来说明本发明的确定方法。
以一个固定安装3个太阳翼，采用单框架控制力矩陀螺为控制执行机构的小卫星为例，其颤振响应确定的具体步骤为：
一、建立卫星颤振动力学方程
根据柔性航天器动力学建模方法，将卫星中心体作为刚体，太阳翼作为附件，考虑太阳翼的柔性，取前6阶模态，建立卫星系统的颤振动力学方程：

式中，I_b为星体转动惯量；ω_b为星体转动角速度；F_sk为附件k柔性对卫星转动的耦合系数；τ_k＝[τ_k1 τ_k2…τ_k6]^T，为附件k前6阶模态坐标列阵；μ为太阳翼附件的阻尼系数；Ω_k＝diag(ω_k1 ω_k2…ω_k6)，为附件k前6阶模态频率列阵。
二、建立主要颤振干扰源单框架控制力矩陀螺动不平衡产生的干扰力矩模型
由定轴转动部件的动不平衡对其转轴产生的干扰力矩出发，结合单框架控制力矩陀螺转子的两自由度转动，即转子转轴相对于框架转动和框架相对于陀螺基座的转动，通过连续的坐标变换获得单框架控制力矩陀螺转子动不平衡对卫星产生的干扰力矩数学模型。
1、陀螺群中的每一个陀螺的框架轴与卫星星体坐标系o_bx_by_bz_b的z_b轴共面，并与z_b轴的夹角为β，框架轴在卫星星体坐标平面x_by_b的投影与x_b轴夹角α_i，i＝1，2…5，则每个陀螺基座坐标系ox_sy_sz_s到星体坐标系o_bx_by_bz_b的转换矩阵M_sbi为

如图2所示，五面锥形单框架陀螺群，陀螺群中的每一个陀螺的框架轴在卫星星体坐标平面x_by_b的投影与x_b轴夹角分别为α₁＝90°，α₂＝π-18°，α₃＝-(π-54°)，α₄＝-54°，α₅＝18°，代入公式(6)中可得五个陀螺基座的安装矩阵。
2、利用转子质量特性测试方法，获得度量每个陀螺转子的主惯性轴相对旋转轴线倾斜程度的量J_xz、J_yz。
3、根据公式组(2)得到陀螺转子的动不平衡度值和初始相位角，

其中，I₀为陀螺转子的动不平衡度值，φ为陀螺转子动不平衡度的初始相位角。
4、卫星控制系统执行机构测量部分测得每个陀螺转子相对自身框架的转角即转子转角γ_i和框架相对陀螺基座的转角即框架转角ζ_i，以及转子绕旋转轴旋转的角速度
5、根据公式(3)得到每个陀螺转子在转子坐标系ox_my_mz_m中干扰力矩，

其中T_dm为陀螺转子的在转子坐标系ox_my_mz_m中干扰力矩，为转子坐标系ox_my_mz_mx_m向、y_m向的单位矢量。
6、根据公式(4)、(5)得到框架坐标系ox_ry_rz_r相对陀螺基座坐标系ox_sy_sz_s的转换矩阵A_rsi和转子坐标系ox_my_mz_m相对框架坐标系ox_ry_rz_r的转换矩阵A_mri


A_rsi是第i个陀螺的框架坐标系相对陀螺基座坐标系的转换矩阵，A_mri是第i个陀螺的转子坐标系相对框架坐标系的转换矩阵；
7、利用公式(6)的转换矩阵M_sbi和公式(4)、(5)的转换矩阵A_rsi、A_mri，将公式(3)经过三次坐标转换，从转子坐标系转化到卫星星体坐标系，得到第i个陀螺转子在卫星星体坐标系内的干扰力矩T_di，
Tdi=MsbiT·ArsiT·AmriT·Tdm---(7).]]>
8、根据公式(8)得到陀螺群在卫星星体坐标系内的干扰力矩T_d，
Td=Σi=1NTdi]]>

对于五面锥形构形的单框力矩陀螺群，忽略转子动不平衡度的初始相位角φ_i，可得其转子动不平衡对星体产生的总干扰力矩为



式中，ζ_i(i＝1，2，…，5)分别为五个框架的转角；γ为转子的转角，五个控制力矩陀螺转子转角相同，由卫星控制系统执行机构测量部分给出。
三、确定柔性航天器在颤振干扰力矩作用下的姿态响应
上述卫星系统的颤振动力学方程组公式(1)自由度数为n＝21，设置广义坐标变量其中θ、ψ为中心体滚动角、俯仰角和偏航角，将系统动力学方程改写为
Mx··+Cx·+Kx=Q]]>
其中


设状态变量y=xx·]]>，将二阶常微分方程转化为状态方程(一阶方程)形式：
y·=Ay+P,]]>A=0n×nIn×n-M-1K-M-1C,]]>P=0n×1M-1Q]]>
利用一阶常微分方程初值问题的求解算法获得卫星在颤振干扰力矩作用下的姿态响应。
四、确定航天器的颤振姿态响应
对于步骤三中确定的姿态响应，其中一部分为可以控制的姿态，按照步骤一建立的小卫星柔性动力学方程，根据卫星姿态控制规律(此为本领域公知常识，具体分析参见《卫星姿态动力学与控制》，宇航出版社，主编屠善澄，1999年第1版，P239～P296)确定可控姿态响应。将步骤三确定的卫星在颤振干扰力矩作用下的姿态响应，减去卫星姿态控制系统可控的姿态响应，最终获得卫星在颤振干扰力矩作用下引起的颤振响应。
根据上述卫星颤振响应的确定方法，确定了该卫星本体在控制力矩陀螺动不平衡干扰下的颤振姿态角，图2-图4分别为该卫星本体的滚动、俯仰、偏航颤振角，由图可见卫星本体沿滚动、俯仰、偏航三个方向的最大颤振角位移分别为：0.0057角秒、0.0071角秒和0.0252角秒，从而为后续进行颤振对成像质量的影响评估提供了输入依据。
本发明未详细说明部分属本领域技术人员公知常识。