本发明涉及卫星实验技术领域，具体地，涉及平面变形绝对-相对欧拉角计算方法及系统。具体地，涉及结构变形引起的指向角度变化的计算方法，具体涉及一种利用平面节点或测点的坐标、位移、欧拉角转序计算平面坐标系指向变化欧拉角的方法，具体涉及一种平面变形绝对-相对欧拉角计算方法。
背景技术：
卫星上高精度设备对安装面的变形有严格的指标要求，卫星结构在设计阶段需要对各种工况下的安装面的指向变形进行量化分析，通常需要精确计算安装面在指定坐标系和转序下的变形欧拉角。
卫星在地面试验前后、经历运输和发射段力学环境前后，设备安装面可能发生变形，安装面指向可能发生变化。卫星在轨运行期间周期性受到太阳辐射等热载荷作用，设备安装面可能发生热变形，安装面指向变化可能影响设备的指向精度，进而影响重要设备甚至整星的正常工作。
申请号为cn201110380055.6的中国发明专利公开了一种基于角速度的欧拉角任意步长正交级数指数型近似输出方法，用于解决现有的飞行器机动飞行时欧拉角输出精度差的技术问题。技术方案是通过引入多个参数并将滚转、俯仰、偏航角速度采用改进的类似切比雪夫正交多项式的递推形式展开逼近，按照依次求解俯仰角、滚转角、偏航角，直接对欧拉角的表达式进行高阶逼近积分，使得欧拉角的求解按照超线性逼近，保证了确定欧拉角的时间更新迭代计算精度，从而提高了惯性设备输出飞行姿态的准确性。
申请号为cn201711194465.5的中国发明专利公开了一种曲面单层网壳梁单元欧拉角的计算方法，包括如下步骤：将曲面单层网壳进行三角化；建立节点与坐标文件node.txt及element.txt；找出与某一节点相连的节点集合，集合中的每一个节点与该节点组成了梁单元，对任意梁单元，求定向第三点坐标矢量；假定全局坐标系中z轴对应的矢量为z＝[0,0,1]，求局部坐标系各轴矢量；求得空间梁的欧拉角beta为局部坐标轴z与主平面为铅直面的空间梁之局部坐标轴z0之间的角度。本发明，可适用计算机编程运算，计算精度和计算效率均很高，可满足大规模曲面单层网壳设计的需求，适用于多种复杂曲面的处理，可处理大规模曲面单层网壳梁单元欧拉角的计算问题。
为了精确计算卫星上重要设备安装面在各种工况下的变形和指向变化，需要对整星进行力学、热变形分析，通过整星有限元模型和分析结果可以观察和评估安装面的位移、变形量等。对于在指定坐标系和转序下的安装面坐标系指向变化的欧拉角，则需要通过对安装面变形前后的节点坐标、位移数据进行处理后计算得到。
本发明的创新点为：利用平面变形前后节点坐标、位移计算参考向量，根据参考向量与平面坐标系之间的关系计算平面变形前后坐标系，避免直接利用节点数据计算坐标系造成坐标系不能反映平面整体变形的问题；将绝度欧拉角、相对欧拉角计算过程分开，将绝度欧拉角计算结果作为相对欧拉角计算输入，降低单独计算相对欧拉角的算法的复杂程度；计算相对欧拉角时引入参考坐标系，使算法简洁且具有模块化特征。
技术实现要素：
针对现有技术中的缺陷，本发明的目的是提供一种平面变形绝对及相对欧拉角计算方法及系统。
根据本发明提供的一种平面变形绝对及相对欧拉角计算方法，包括以下步骤：
步骤s1：根据平面节点的编号、坐标，寻找距离最大的两个节点x1、x2，寻找节点y2，满足向量
步骤s2：以y2到x1、x2两点的距离|x1y2|和|x2y2|中较小的值的三分之一为半径r，分别以节点x1、x2、y2为圆心，在平面内寻找半径r范围内的所有节点，分别组成3组节点，并根据3组节点的坐标和位移，计算平面变形前3个中心点
步骤s3：计算平面变形前2个参考向量
步骤s4：根据预先定义的平面初始坐标系参考点o0、x0、y0，计算平面变形前后坐标系的x轴向量；
步骤s5：使用平面所有节点变形前后的坐标数据和最小二乘法计算平面变形前后的法向量，即平面变形前后的z轴向量；
步骤s6：根据平面变形前后的x轴向量和法向量，计算平面变形前后的y轴向量；
步骤s7：计算平面变形前后的x、y、z轴向量的单位向量，并构建平面变形前后坐标系；
步骤s8：根据平面变形前后坐标系和指定转序计算相应的绝对欧拉角；
步骤s9：根据两个平面，第一平面和第二平面，同一工况下变形前后的坐标系的关系，计算第二平面的参考坐标系；
步骤s10：按照指定转序，计算第二平面变形后坐标系相对于第二平面参考坐标系的欧拉角，即得第二平面相对于第一平面变形的相对欧拉角。
优选地，所述y1是指一个点，满足x1、x2、y1、y2之间的向量关系，即四个点构成一个平行四边形，x1与x2、y1与y2分别为平行四边形两组对角点，x1、x2、y2为平面内搜寻找到的节点，而y1点为平面内满足向量关系的点；
所述预设角度范围为：30°～150°之间。
优选地,所述步骤s8：按照一定转序计算由一个坐标系到另一个坐标系的欧拉角，转序为绕坐标系3个不同坐标轴的转动顺序，按照转序定义转轴1、2、3与坐标系各轴的对应关系。
优选地,所述步骤s8：平面变形绝对欧拉角计算时，先计算平面变形后坐标系第3转轴向量在变形前坐标系第2、3转轴平面内的投影向量，投影向量与变形前坐标系第3转轴的夹角为绕第1转轴转动的欧拉角α。
优选地,所述步骤s8：平面变形绝对欧拉角计算时，平面变形前坐标系绕第1转轴转动α后得到中间过程坐标系，中间过程坐标系第3转轴与变形后坐标系第3转轴的夹角为绕第2转轴转动的欧拉角β，中间过程坐标系第2转轴与平面变形后坐标系第2转轴的夹角为绕第3转轴的欧拉角γ。
优选地，欧拉角的正负由相关向量计算结果的正负决定。
根据本发明提供的一种平面变形绝对及相对欧拉角计算系统，包括以下模块：
模块s1：根据平面节点的编号、坐标，寻找距离最大的两个节点x1、x2，寻找节点y2，满足向量
模块s2：以y2到x1、x2两点的距离|x1y2|和|x2y2|中较小的值的三分之一为半径r，分别以节点x1、x2、y2为圆心，在平面内寻找半径r范围内的所有节点，分别组成3组节点，并根据3组节点的坐标和位移，计算平面变形前3个中心点
模块s3：计算平面变形前2个参考向量
模块s4：根据预先定义的平面初始坐标系参考点o0、x0、y0，计算平面变形前后坐标系的x轴向量；
模块s5：使用平面所有节点变形前后的坐标数据和最小二乘法计算平面变形前后的法向量，即平面变形前后的z轴向量；
模块s6：根据平面变形前后的x轴向量和法向量，计算平面变形前后的y轴向量；
模块s7：计算平面变形前后的x、y、z轴向量的单位向量，并构建平面变形前后坐标系；
模块s8：根据平面变形前后坐标系和指定转序计算相应的绝对欧拉角；
模块s9：根据两个平面，第一平面和第二平面，同一工况下变形前后的坐标系的关系，计算第二平面的参考坐标系；
模块s10：按照指定转序，计算第二平面变形后坐标系相对于第二平面参考坐标系的欧拉角，即得第二平面相对于第一平面变形的相对欧拉角。
优选地，所述y1是指一个点，满足x1、x2、y1、y2之间的向量关系，即四个点构成一个平行四边形，x1与x2、y1与y2分别为平行四边形两组对角点，x1、x2、y2为平面内搜寻找到的节点，而y1点为平面内满足向量关系的点；
所述预设角度范围为：30°～150°之间。
优选地,所述模块s8：按照一定转序计算由一个坐标系到另一个坐标系的欧拉角，转序为绕坐标系3个不同坐标轴的转动顺序，按照转序定义转轴1、2、3与坐标系各轴的对应关系；
所述模块s8：平面变形绝对欧拉角计算时，先计算平面变形后坐标系第3转轴向量在变形前坐标系第2、3转轴平面内的投影向量，投影向量与变形前坐标系第3转轴的夹角为绕第1转轴转动的欧拉角α；
所述模块s8：平面变形绝对欧拉角计算时，平面变形前坐标系绕第1转轴转动α后得到中间过程坐标系，中间过程坐标系第3转轴与变形后坐标系第3转轴的夹角为绕第2转轴转动的欧拉角β，中间过程坐标系第2转轴与平面变形后坐标系第2转轴的夹角为绕第3转轴的欧拉角γ。
优选地，欧拉角的正负由相关向量计算结果的正负决定。
与现有技术相比，本发明具有如下的有益效果：
本发明通过多点平均法构造参考向量，根据平面变形前坐标系与参考向量之间的关系修正变形后平面坐标系，根据坐标系转序构建参考坐标系，针对不同转序下一个向量绕另一向量转动一定角度、坐标系绕某一向量转动一定角度编写统一算法，提高了平面变形绝对-相对欧拉角的计算精度和效率。
附图说明
通过阅读参照以下附图对非限制性实施例所作的详细描述，本发明的其它特征、目的和优点将会变得更明显：
图1是本发明单个平面变形前后绝对欧拉角计算流程示意图。
图2是本发明两个平面变形前后相对欧拉角计算流程示意图。
具体实施方式
下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变化和改进。这些都属于本发明的保护范围。
根据本发明提供的一种平面变形绝对及相对欧拉角计算方法，包括以下步骤：
步骤s1：根据平面节点的编号、坐标，寻找距离最大的两个节点x1、x2，寻找节点y2，满足向量
步骤s2：以y2到x1、x2两点的距离|x1y2|和|x2y2|中较小的值的三分之一为半径r，分别以节点x1、x2、y2为圆心，在平面内寻找半径r范围内的所有节点，分别组成3组节点，并根据3组节点的坐标和位移，计算平面变形前3个中心点
步骤s3：计算平面变形前2个参考向量
步骤s4：根据预先定义的平面初始坐标系参考点o0、x0、y0，计算平面变形前后坐标系的x轴向量；
步骤s5：使用平面所有节点变形前后的坐标数据和最小二乘法计算平面变形前后的法向量，即平面变形前后的z轴向量；
步骤s6：根据平面变形前后的x轴向量和法向量，计算平面变形前后的y轴向量；
步骤s7：计算平面变形前后的x、y、z轴向量的单位向量，并构建平面变形前后坐标系；
步骤s8：根据平面变形前后坐标系和指定转序计算相应的绝对欧拉角；
步骤s9：根据两个平面，平面1和平面2，同一工况下变形前后的坐标系的关系，计算平面2的参考坐标系；
步骤s10：按照指定转序，计算平面2变形后坐标系相对于平面2参考坐标系的欧拉角，即得平面2相对于平面1变形的相对欧拉角。
具体地，所述y1是指一个点，满足x1、x2、y1、y2之间的向量关系，即四个点构成一个平行四边形，x1与x2、y1与y2分别为平行四边形两组对角点，x1、x2、y2为平面内搜寻找到的节点，而y1点为平面内满足向量关系的点；
所述预设角度范围为：30°～150°之间。
具体地,所述步骤s8：按照一定转序计算由一个坐标系到另一个坐标系的欧拉角，转序为绕坐标系3个不同坐标轴的转动顺序，按照转序定义转轴1、2、3与坐标系各轴的对应关系。
具体地,所述步骤s8：平面变形绝对欧拉角计算时，先计算平面变形后坐标系第3转轴向量在变形前坐标系第2、3转轴平面内的投影向量，投影向量与变形前坐标系第3转轴的夹角为绕第1转轴转动的欧拉角α。
具体地,所述步骤s8：平面变形绝对欧拉角计算时，平面变形前坐标系绕第1转轴转动α后得到中间过程坐标系，中间过程坐标系第3转轴与变形后坐标系第3转轴的夹角为绕第2转轴转动的欧拉角β，中间过程坐标系第2转轴与平面变形后坐标系第2转轴的夹角为绕第3转轴的欧拉角γ。
具体地，欧拉角的正负由相关向量计算结果的正负决定。
根据本发明提供的一种平面变形绝对及相对欧拉角计算系统，包括以下模块：
模块s1：根据平面节点的编号、坐标，寻找距离最大的两个节点x1、x2，寻找节点y2，满足向量
模块s2：以y2到x1、x2两点的距离|x1y2|和|x2y2|中较小的值的三分之一为半径r，分别以节点x1、x2、y2为圆心，在平面内寻找半径r范围内的所有节点，分别组成3组节点，并根据3组节点的坐标和位移，计算平面变形前3个中心点
模块s3：计算平面变形前2个参考向量
模块s4：根据预先定义的平面初始坐标系参考点o0、x0、y0，计算平面变形前后坐标系的x轴向量；
模块s5：使用平面所有节点变形前后的坐标数据和最小二乘法计算平面变形前后的法向量，即平面变形前后的z轴向量；
模块s6：根据平面变形前后的x轴向量和法向量，计算平面变形前后的y轴向量；
模块s7：计算平面变形前后的x、y、z轴向量的单位向量，并构建平面变形前后坐标系；
模块s8：根据平面变形前后坐标系和指定转序计算相应的绝对欧拉角；
模块s9：根据两个平面，平面1和平面2，同一工况下变形前后的坐标系的关系，计算平面2的参考坐标系；
模块s10：按照指定转序，计算平面2变形后坐标系相对于平面2参考坐标系的欧拉角，即得平面2相对于平面1变形的相对欧拉角。
具体地，所述y1是指一个点，满足x1、x2、y1、y2之间的向量关系，即四个点构成一个平行四边形，x1与x2、y1与y2分别为平行四边形两组对角点，x1、x2、y2为平面内搜寻找到的节点，而y1点为平面内满足向量关系的点；
所述预设角度范围为：30°～150°之间。
具体地,所述模块s8：按照一定转序计算由一个坐标系到另一个坐标系的欧拉角，转序为绕坐标系3个不同坐标轴的转动顺序，按照转序定义转轴1、2、3与坐标系各轴的对应关系；
所述模块s8：平面变形绝对欧拉角计算时，先计算平面变形后坐标系第3转轴向量在变形前坐标系第2、3转轴平面内的投影向量，投影向量与变形前坐标系第3转轴的夹角为绕第1转轴转动的欧拉角α；
所述模块s8：平面变形绝对欧拉角计算时，平面变形前坐标系绕第1转轴转动α后得到中间过程坐标系，中间过程坐标系第3转轴与变形后坐标系第3转轴的夹角为绕第2转轴转动的欧拉角β，中间过程坐标系第2转轴与平面变形后坐标系第2转轴的夹角为绕第3转轴的欧拉角γ。
具体地，欧拉角的正负由相关向量计算结果的正负决定。
下面通过优选例，对本发明进行更为具体地说明。
优选例1：
本发明针对使用结构有限元法分析平面变形后无法直接得到平面变形指向变化欧拉角的缺陷，提出了一种根据结构有限元分析结果数据计算平面变形欧拉角的方法。本发明利用结构有限元分析结果中平面上节点的编号、初始坐标、位移数据，计算平面的绝对欧拉角；根据各个平面绝对欧拉角计算过程中产生的变形前后坐标系，计算多个平面相对变形的相对欧拉角。
根据本发明提供的一种平面变形绝对-相对欧拉角计算方法，包括以下步骤：
步骤一，根据平面节点的编号、坐标，寻找距离最大的两个节点x1、x2，寻找节点y2，满足向量
步骤二，分别以节点x1、x2、y2为圆心，以y2到x1、x2两点的距离|x1y2|和|x2y2|中较小的值的三分之一为半径，在平面内寻找3组节点(解释：1组节点举例：以x1为圆心，r为半径，所围成的圆内部的所有节点为1组节点。另两组节点则分别以x2、y2为圆心，r为半径，所围成的圆内的节点)，并根据3组节点的坐标和位移，计算平面变形前3个中心点
步骤三，计算平面变形前2个参考向量
步骤四，根据预先定义的平面初始坐标系参考点o0、x0、y0，计算平面变形前后坐标系的x轴向量。
步骤五，使用平面所有节点变形前后的坐标数据和最小二乘法计算平面变形前后的法向量，即平面变形前后的z轴向量。
步骤六，根据平面变形前后的x轴向量和法向量，计算平面变形前后的y轴向量。
步骤七，计算平面变形前后的x、y、z轴向量的单位向量，并构建平面变形前后坐标系。
步骤八，根据平面变形前后坐标系和指定转序计算相应的绝对欧拉角。
步骤九，根据两个平面(平面1和平面2)同一工况下变形前后的坐标系的关系，计算平面2的参考坐标系。
步骤十，按照指定转序，计算平面2变形后坐标系相对于平面2参考坐标系的欧拉角，即得平面2相对于平面1变形的相对欧拉角。通过步骤十中的两个坐标系计算两个平面变形相对欧拉角的计算流程与步骤八中通过两个坐标系计算绝对欧拉角的流程相似。
本发明对多个平面相对欧拉角的计算按照模块化设计思路拆分为单个平面的绝对欧拉角算法和多个平面相对欧拉角的快速算法。在计算单个平面的欧拉角时，已得到相对欧拉角计算所需的数据，避免计算相对欧拉角时再次处理平面节点编号、坐标、位移数据，提高计算效率。
优选例2：
下面对本发明的实施例作详细说明，本实施例以本发明技术方案为前提进行实施，给出了详细的实施方式和具体计算过程，但本发明的保护范围不限于下述的实施例。
本实施例的计算过程包括如图1所示的单个平面绝对欧拉角计算过程和图2所示的2个平面相对欧拉角计算过程。多个平面的相对欧拉角计算过程与图2相同，需要逐次计算两个平面相对欧拉角。
步骤一，根据平面节点的坐标，寻找3个节点x1、x2、y2，3个节点的坐标满足以下关系：
|x1x2|≥|xixj|,(xi,xj为平面内任意两个节点)
步骤二，计算平面变形前3个中心点
其中，
n1表示第1组节点的数量，平面内以x1为圆心，r为半径的圆内的节点数量。
n2表示第2组节点的数量，平面内以x2为圆心，r为半径的圆内的节点数量。
n3表示第3组节点的数量，平面内以y2为圆心，r为半径的圆内的节点数量。
式(2)中x、y、p表示节点，包含变形前坐标数据；u表示节点位移；上标0表示变形前，上标1表示变形后，(变形前后，节点编号相同，原始坐标相同，只是变形前无位移，变形后有位移)下同。
步骤三，计算平面变形前2个参考向量和变形后2个参考向量，变形前后参考向量满足以下关系：
步骤四，根据平面初始坐标系参考点o0(坐标系原点)、x0(x轴正向上一点)、y0(平面内一点，且∠x0o0y0∈(30°,150°))，计算平面变形前后坐标系的x轴向量x0和x1，计算过程如下：
步骤五，使用平面所有节点变形前后的坐标数据和最小二乘法拟合变形前后平面方程，并计算平面变形前后的法向量n0和n1，即平面变形前后的z轴向量z0和z1(最小二乘法是通用算法，可以计算拟合平面方程，由平面方程求平面法向量是常规算法)。
步骤六，根据平面变形前后的x轴向量x0和x1和z轴向量z0和z1，计算平面变形前后的y轴向量y0和y1，如下所示。
y0＝z0×x0
y1＝z1×x1(5)
步骤七，计算平面变形前后的x、y、z轴向量的单位向量，并构建平面变形前后坐标系s0和s1，计算过程如下：
其中，
vi表示变形前后与x、y、z轴平行的向量，下标0、1表示变形前、变形后。vi表示x0、x1、y0、y1、z0、z1(在式(5)及前文已求出)中任意一项。
步骤八，根据平面变形前后坐标系s0和s1和指定转序(例如先后绕z、x、y轴转动)，计算相应转序下的绝对欧拉角α、β、γ，计算过程如下：
1)调整平面变形坐标系的坐标轴向量排序。
其中，
2)计算平面变形后坐标系第3转轴在平面变形前坐标系第2、3转轴平面内的投影向量及其长度。
a，b，c分别表示方程的系数，解线性方程组可得。
l表示变形后轴3在变形前2、3轴平面内的投影向量的模或长度。
3)计算欧拉角α、β、γ。
如果l＝0，则按照下式计算：
β＝90°
γ＝0°(9)
如果l≠0，则按照下式计算：
将坐标系s0绕
欧拉角γ按照下式计算：
4)欧拉角α、β、γ单位转换。
式(9～12)所得欧拉角单位为度(°)，可根据角度单位关系转换为弧度、微弧度、角秒等。
步骤九，根据两个平面(平面1和平面2)同一工况下变形前后的坐标系的关系，计算平面2的参考坐标系，计算过程如下：
平面2和平面2变形前后的坐标系计算方法同式(1～6)。
步骤十，按照指定转序，计算平面2变形后坐标系
在本申请的描述中，需要理解的是，术语“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本申请和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本申请的限制。
本领域技术人员知道，除了以纯计算机可读程序代码方式实现本发明提供的系统、装置及其各个模块以外，完全可以通过将方法步骤进行逻辑编程来使得本发明提供的系统、装置及其各个模块以逻辑门、开关、专用集成电路、可编程逻辑控制器以及嵌入式微控制器等的形式来实现相同程序。所以，本发明提供的系统、装置及其各个模块可以被认为是一种硬件部件，而对其内包括的用于实现各种程序的模块也可以视为硬件部件内的结构；也可以将用于实现各种功能的模块视为既可以是实现方法的软件程序又可以是硬件部件内的结构。
以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变化或修改，这并不影响本发明的实质内容。在不冲突的情况下，本申请的实施例和实施例中的特征可以任意相互组合。