Thermal-structural-optical integrated analysis method based on the complete equations of rigid body motion

光机系统受到重力、热、结构力等内外载荷的作用发生结构变形，同时光学镜面面形也发生变化，导致其系统光学性能下降^[1]。光机热集成分析是评估这些载荷对系统光学性能影响的重要手段，已广泛应用于空间相机^[2-4]、航空相机^[5-6]、大型陆基望远镜^[7]、工业相机^[8]等光机系统设计阶段，特别是工作于空间恶劣环境的高精度光学相机，光机热集成分析结果已成为相机机械结构、光学系统、热系统等子系统设计甚至光校方案的重要评价依据。

目前，光机热集成分析的方法大致可分为两类：直接拟合法和分离刚体法，前者采用拟合方法直接对光学镜面变形误差进行曲面拟合，然后输入到光学系统以分析光学性能的变化^[9-10]。参考文献[11]采用非均匀有理B样条曲线并结合泽尼克（Zernike）多项式，对圆形光学元件面形误差曲面进行拟合和边缘延展。参考文献[12]应用泽尼克多项式拟合算法，构建光学面变形与光学分析数据的转换接口，获取了热不敏红外系统在高低温工况下的光学系统传递函数。对于结构变形相对很小的光机系统，直接拟合法可以得到较好的曲面表征效果，而在实际光机工程中，光机系统的结构变形量通常大于光学镜面的变形量，即光学镜片的刚体位移量大于其自身镜面的变形量，将光学镜面离焦、偏心、倾斜等刚体位移拟合在光学镜面弹性变形中，会使得光学镜面拟合结果不符合实际情况，从而导致较大的计算误差。此外，光学镜面的刚体位移量是光学系统光校的重要参考信息，因此，分离光学镜面变形误差中的刚体位移很有必要。

分离刚体法由美国学者Doyle K B等^[13]提出，该方法首先分离出由整体结构变形导致的光学镜面刚体位移，然后采用多项式拟合方法对光学镜面弹性变形进行曲面拟合，最后以刚体运动量作为光学镜面离焦、偏心、倾斜的输入，以拟合多项式系数作为光学镜面弹性变形的输入，将光学面的变形量转换为光学系统的参数以评价光学性能的变化。参考文献[14]在低热变形水冷镜的优化设计中，应用分离刚体法分析水冷镜在不同温度场分布情况下的面变形误差。参考文献[15]结合拓扑优化方法和分离刚体法，建立了多部件光机系统的拓扑优化方法。基于分离刚体法，参考文献[16]研究了微振动对航空相机光学性能的影响。自分离刚体法提出后，大多数研究集中在该方法的应用上，对其适用度和局限性的研究很少，而对于结构复杂、板件尺寸大、光学镜片多、环境载荷恶劣的光机系统，其光学镜面的变形模式复杂、刚体变形较大，该分离刚体法的计算结果与实际情况有较大误差。此外，光机热集成分析方法集成了结构力学、光学、热力学等学科的技术，涉及有限元求解、刚体运动计算、弹性变形表征、光学系统评价等内容，其中的有限元单元划分、光学镜面修正、刚体运动方程以及光学镜面面形等势必会影响最终的光学性能分析结果，但目前这些因素对于分析结果的影响程度尚不清楚。

文中提出基于刚体运动完备方程的光机热集成分析方法，旨在提高分离刚体法的计算精度。创新点及贡献主要体现在以下几个方面：（1）采用双立方样条插值法对光学镜面进行修正，避免了对光学面解析方程的依赖，弱化了光学镜面弹性变形的假设条件；（2）提出基于刚体运动完备方程的优化求解方法，提高了刚体运动量的计算精度；（3）明确光学镜面修正方法、刚体运动方程、点阵分布以及镜面面形对求解结果的影响；（4）展示所提方法在工程案例中的应用过程。

光学镜面的空间变形主要由光学面的刚体运动和弹性变形两部分组成，刚体运动主要包括旋转运动和平移运动，由光机系统的结构受载变形而带动镜室整体偏移导致，镜面的弹性变形是镜室变形挤压和镜片自身受温度载荷变形导致的。

记光学系统的全局坐标系为O，光学面顶点坐标系为S，且光学面顶点坐标系原点在全局坐标系中的坐标为(X_S, Y_S, Z_S)。为了方便描述光学面变形方程并消除光学面变形量的累计误差，选取全局坐标系作为所有光学面顶点坐标系的基准。将光学面顶点坐标系S原点移动至全局坐标系原点，并按照X、Y和Z轴的顺序围绕原点旋转到与全局坐标系重合，该过程的变换矩阵记为M_S→O。光学面点阵的顶点坐标系经过变换矩阵M_S→O操作后，最终与全局坐标系完全一致，光学面点阵的全局坐标变为顶点坐标系中的坐标，如光学面上第m个点在光学面顶点坐标系中的坐标可表示为：

式中：下标“Pm,O”代表第m个点在全局坐标系下的齐次坐标；下标“Pm,S”代表在顶点坐标系下的齐次坐标。

变换矩阵M_S→O可表示为：

式中：T_XYZ为平动变换矩阵；R_X、R_Y和R_Z分别为围绕X、Y和Z轴旋转的变换矩阵：

式中：θ_X为X轴旋转角；θ_Y为Y轴旋转角；θ_Z为Z轴旋转角。对于光学面顶点坐标系，可知$ {T_X}{\text{ = }} - {X_S} $，$ {T_Y}{\text{ = }} - {Y_S} $，$ {T_Z}{\text{ = }} - {Z_S} $。

系统结构发生变形后，光学面发生刚体运动和弹性变形，其中刚体运动可以用镜面顶点坐标系的运动表示，也即顶点坐标系的原点由(X_S, Y_S, Z_S)运动至(X_S’, Y_S’, Z_S’)，且围绕顶点坐标系X、Y和Z轴分别θ_X、θ_Y和θ_Z的顺序旋转，光学面变形前、后的顶点坐标系示意图如图1所示。为了保持与变形前的顶点坐标系有相同的运动基准，对变形后的光学面点阵进行M_S→O运动变换，变形后的光学面点阵在变形前的顶点坐标系中可表示为：

Figure 1.  Schematic diagram of optical surface vertex coordinate system

光学面上第m个点在顶点坐标系中的空间运动向量可以表示为：

光学面上第m个点的空间运动由镜面刚体运动和弹性变形导致，其刚体运动矩阵用M_R表示，弹性变形用d_m表示，则有以下关系式：

式中：${{\boldsymbol{d}}_m}{\text{ = }}{\left[ {d{X_{Pm,S}}\;d{Y_{Pm,S}}\;d{Z_{Pm,S}}\;{\text{0}}} \right]^{\rm{T}}}$。

对于有m个离散点的光学面，则可以通过矩阵表示：

式中：${{\boldsymbol{A}}_S} {\text{ = }} \left[ {\begin{array}{*{20}{c}} {{X_{P{\text{1}},S}}} & \cdots &{{X_{Pm,S}}} \\ {{Y_{P{\text{1}},S}}}& \cdots &{{Y_{Pm,S}}} \\ {{Z_{P{\text{1}},S}}}& \cdots &{{Z_{Pm,S}}} \\ 1& \cdots &1 \end{array}} \right]$和${{\boldsymbol{A'}}_S} {\text{ = }} \left[ {\begin{array}{*{20}{c}} {{{X'}_{P{\text{1}},S}}}& \cdots &{{{X'}_{Pm,S}}} \\ {{{Y'}_{P{\text{1}},S}}}& \cdots &{{{Y'}_{Pm,S}}} \\ {{{Z'}_{P{\text{1}},S}}}& \cdots &{{{Z'}_{Pm,S}}} \\ 1& \cdots &1 \end{array}} \right]$分别为光学面变形前、后的点阵在光学面顶点坐标系中的坐标；${{\boldsymbol{D}}_S}{\text{ = }}\left[ {{{\boldsymbol{d}}_1}\;{{\boldsymbol{d}}_2}\; \cdots \;{{\boldsymbol{d}}_m}} \right]$。

则公式（8）可以扩展为：

式中：${{\boldsymbol{P}}_S}{\text{ = }}\left[ {{{\boldsymbol{p}}_{1,S}}\;{{\boldsymbol{p}}_{2,S}}\; \cdots \;{{\boldsymbol{p}}_{m,S}}} \right]$。

将公式（11）代入公式（10）可得：

光学面点阵的弹性变形矩阵可表示为：

式中：${{\boldsymbol{\tilde P}}_S}{\text{ = }}{{\boldsymbol{M}}_{\text{R}}}{{\boldsymbol{A}}_S} - {{\boldsymbol{A}}_S}$为光学面点阵的刚体运动位移矩阵，此处的P_S表示光学面光线入射接触点的位移向量，而不是有限元法求解的光学面点阵的位移。如图2所示，光线与镜面接触区域由ab变为a'b'，矢向光线在光学面变形前的入射点为$ ({X_{Pm,S}},{Y_{Pm,S}},{Z_{Pm,S}}) $，该点在光学面变形后偏移至$ ({\bar X_{Pm,S}},{\bar Y_{Pm,S}},{\bar Z_{Pm,S}}) $，但光学面变形后的入射点其实是$ ({X'_{Pm,S}},{Y'_{Pm,S}},{Z'_{Pm,S}}) $，因而需要将光学面ab修正为光学面a'b'。

Figure 2.  Schematic diagram of light incident area before and after optical surface deformation

有限元法求解的节点位移为点$({\bar X_{Pm,S}},{\bar Y_{Pm,S}}, $$ {\bar Z_{Pm,S}}) $和点$ ({X_{Pm,S}},{Y_{Pm,S}},{Z_{Pm,S}}) $的差值（坐标转换后的结果）：

式中：${{{\bar {\boldsymbol{A}}}}_S}{\text{ = }}\left[ {\begin{array}{*{20}{c}} {{{\bar X}_{P{\text{1}},S}}}& \cdots &{{{\bar X}_{Pm,S}}} \\ {{{\bar Y}_{P{\text{1}},S}}}& \cdots &{{{\bar Y}_{Pm,S}}} \\ {{{\bar Z}_{P{\text{1}},S}}}& \cdots &{{{\bar Z}_{Pm,S}}} \\ 1& \cdots &1 \end{array}} \right]$。

光学面的刚体运动量可直接用于光机系统的光校调试，由于光学面的弹性变形无法用刚体运动表征，去除刚体运动后的镜面弹性变形量是最小值，即${{\boldsymbol{D}}_S}$的F范数${\left\| {{{\boldsymbol{D}}_S}} \right\|_{\text{F}}}$ 最小，可通过优化过程对刚体运动量进行求解。刚体运动量共有六个设计变量，其中三个平动变量为$ \Delta {X_S} $、$ \Delta {Y_S} $和$ \Delta {Z_S} $，三个转动变量为$ \Delta {\theta _X} $、$ \Delta {\theta _Y} $和$ \Delta {\theta _Z} $。为了方便描述刚体运动量的求解过程，采用${{\boldsymbol{D}}_S}$的F范数的平方作为目标最小值：

通过数值寻优法对刚体运动量进行求解，优化方程为：

Minimize: ${\left( {{{\left\| {{{\boldsymbol{D}}_S}} \right\|}_{\text{F}}}} \right)^{\text{2}}}$

Subject to:

在公式（16）中，首先需要确定光学面光线入射接触点的位移向量P_S，由于光机系统结构变形较小，位移向量P_S与有限元法求解的位移${{\boldsymbol{\bar P}}_S}$较为接近，可对${{\boldsymbol{\bar P}}_S}$进行双立方样条插值求解P_S。双立方样条插值方法用Green函数构建响应面与精确数据点之间的关系，响应面则由Green函数的线性组合构成，双立方样条插值法详细介绍见参考文献[17]。

在光学分析中，可通过光学面顶点坐标系或者辅助光学面顶点坐标系的离轴、厚度、偏转角等参数项定义刚体运动的顺序，而弹性变形是光学面上的微小波动，基于光学干涉检验的概念，可视之为加工或检测时的镜面畸变，需通过光学面顶点坐标系对其定义。目前，有多种光学面弹性变形的表征函数，其中应用最为广泛的是Zernike函数。以Fringe Zernike多项式为例，其表达式为：

式中：$ {A_i} $为第i项系数；$\; \rho $为归一化径向坐标；$ \varphi $为角度坐标。Fringe Zernike多项式的详细内容见参考文献[13]。

对于具有m个节点的光学面，其拟合方程为：

式中：${\boldsymbol{K}} {\text{ = }} \left[ {\begin{array}{*{20}{c}} {{F_1}\left( {{\rho _1},{\varphi _1}} \right)}&{{F_2}\left( {{\rho _1},{\varphi _1}} \right)}& \cdots &{{F_n}\left( {{\rho _1},{\varphi _1}} \right)} \\ {{F_1}\left( {{\rho _2},{\varphi _2}} \right)}&{{F_2}\left( {{\rho _2},{\varphi _2}} \right)}& \cdots &{{F_n}\left( {{\rho _2},{\varphi _2}} \right)} \\ \vdots & \vdots & \ddots & \vdots \\ {{F_1}\left( {{\rho _m},{\varphi _m}} \right)}&{{F_2}\left( {{\rho _m},{\varphi _m}} \right)}& \cdots &{{F_n}\left( {{\rho _m},{\varphi _m}} \right)} \end{array}} \right] $；${\boldsymbol{C}} = {\left[ {{A_1}\;{A_2}\; \cdots \;{A_n}} \right]^{\rm{T}}}$；${{\boldsymbol{{Z}}}_{\text{d}}}$为光学面顶点坐标系下的Z向弹性变形量，通过${\boldsymbol{M}}_{\text{R}}^{ - 1}{{\boldsymbol{D}}_S}$计算得到。

Fringe Zernike多项式的系数由公式（19）求解：

光学面误差分析结果可为光学系统在温度、重力、振动激励等载荷作用下的光学性能变化分析提供输入，若仅评价光学面变形误差对光学性能的影响，则不需要计算光学面刚体运动量，只需对变形后光学面进行表征。然而，刚体运动量可为光学面变形误差分析提供更深入的视角，区分相机结构变形和光学面自身变形，对光学系统的结构改进、材料替换等具有指导作用，并且刚体运动量可直接应用于光学系统的光校。

刚体运动求解是分离刚体法的关键步骤，在传统的分离刚体法中，对E求偏导，以最优解处为0构建方程组求解刚体运动。传统方法假设刚体运动很小，在方程组中，用转角代替其正弦值，用1代替转角的余弦值，并且只保留刚体运动变量的一次项，省略变量的二次及其以上高阶次项。如，对X方向平动变量$ \Delta {X_S} $求偏导，得到如下公式：

依次对$ \Delta {Y_S} $、$ \Delta {Z_S} $、$ \Delta {\theta _X} $、$ \Delta {\theta _Y} $和$ \Delta {\theta _Z} $求偏导，可以得到其他五组方程，联合六个方程组，可求解得到该六个刚体运动量。

在求解刚体运动之前，传统方法采用切线法或法线法对光学面进行修正。如图3所示，以矢向光线为例，需将$ c $点变形后的位置由$ \bar c $点修正为$ c' $点，在此过程中，计算$ c $点处的切线或法线需要光学面的解析方程作为输入。该切线或法线方法是基于以下假设：温度升高会引起光学面曲率半径的增大，图3中的光学面点阵将朝着−Z方向变动，但在复杂的光学系统中，光机结构及光学镜面受载变形并不完全满足该假设，此时用$ c $点处的切线或法线预测$ c' $点将会导致一定的误差。此外，传统方法用简化后方程组（如公式（20））求解刚体运动，如图3所示，$ c $点在刚体运动后的位置为$ {c^ * } $点，$ {c^ * } $点对应于修正后的$ d' $点，而传统方法依然以$ c' $点作为目标点，其求解的刚体运动量不是最优解。

Figure 3.  Schematic diagram of light incident area before and after optical surface deformation

与传统方法相比，文中所提方法采用双立方样条插值法对光学面变形进行修正，基于刚体运动完备方程，通过数值优化过程计算得到刚体运动的精确值，其求解过程只需要光学面点阵的坐标及其位移，而不需要光学面的解析方程。所提方法共有六个实施步骤：（a）在有限元软件中，建立对应于光学系统的相机结构有限元模型，施加载荷及边界约束条件，计算光学面上的节点位移；（b）按照有限元模型及输出结果文件格式，提取光学面点阵的坐标及其位移向量；（c）以光学面点阵的坐标及其位移向量作为输入，用双立方样条插值法对光学面变形误差进行修正；（d）利用优化方程（16）求解光学镜面的刚体运动量；（e）建立光学面弹性变形的Zernike多项式方程，利用公式（19）计算Zernike系数；（f）将光学面刚体运动量和弹性变形的Zernike系数，以光学分析软件的格式构建光学面变形误差模型，分析其对光学系统的性能影响。所提光学面变形误差计算方法的实施流程如图4所示。

Figure 4.  Flow chart of the proposed optical surface deformation error calculation method

对单个光学面施加理想刚体运动和弹性变形的形式定义光学面变形误差，用所提光学面变形误差计算方法和传统方法求解刚体运动量和弹性变形量，对比两种方法的计算精度，并分析刚体运动方程、变形修正方法等因素对计算结果的影响。

该光学面是标准球面，其矢高表达式为：

式中：c为曲率；R为径向坐标；k为圆锥曲线常数。在该案例中，不失一般地，设c为0.2，k为0.5，光学面径向最大值为1（即R_max为1），如图5所示。

Figure 5.  Optical surface - Standard sphere (c=0.2, k=0.5)

为了避免节点数量过少而影响光学面变形误差的计算精度，在X−Y平面径向和周向分别均匀采样101个节点，共有10001个节点。设置五种情形：刚体运动I（大位移）、刚体运动II（小位移）、刚体运动III（小位移和转角）、弹性变形、刚体运动+弹性变形，刚体运动量测试值如表1所示。由表1可知，所提方法准确地计算出刚体运动量，平动、转动分量与其测试值几乎完全相等。传统方法的计算精度相对偏低，在三个情形中，T_z的求解值与测试值有较大的偏差，在第三个情形中，$ {\theta _x} $、$ {\theta _y} $与测试值都有一定的偏差。

该案例中的光学面口径为1，为了分析光学面变形误差计算精度与光学面口径的相关性，分别将光学面口径扩大为10和100，相应地，c分别改为0.02和0.002，k仍然为0.5，计算结果如表2和表3所示。在三个情形中，所提方法都准确地计算出刚体运动量，而传统方法求得的T_z、$ {\theta _y} $与其测试值的偏差都相对较大，尤其在第三个情形中，T_z的求解值与其测试值的偏差很大，$ {\theta _x} $、$ {\theta _y} $也有较大的偏差。传统方法利用光学面的径向切线计算光学面的修正面，当存在刚体转动且光学面口径较大时，其光学面的修正面会产生较大的偏差，导致刚体运动量的计算误差偏大，如表1~3所示的第三情形的计算结果。此外，传统方法采用刚体运动的近似方程求解刚体运动量，当刚体转动量较大时，也会产生额外的误差。所提方法利用双立方样条插值法修正光学面，不需要计算光学面的径向或法向切线，并通过刚体运动完备方程求解刚体运动量，计算精度高且与光学面口径大小几乎没有关联。需要强调的是，所提算法和传统方法的算法在商业软件Matlab中的执行时间几乎相等。

在第四个情形中，光学面没有刚体运动，只有弹性变形。光学面弹性变形测试函数为Fringe Zernike多项式中的第四项（如图6所示），其表达式为：

Figure 6.  The fourth term of Fringe Zernike polynomial (Taking k₄=1 as an example)

式中：k₄为ZK₄项的系数；R为径向坐标；$ \phi $为角坐标。

由于光学面变形通常很小，该案例中的k₄设为0.001，也即光学面局部变形的最大值为0.001，计算结果如表1所示，相较于光学面的弹性变形量，两种方法求解的刚体运动量都很小且几乎相同，可以认为光学面无刚体运动，所提方法和传统方法都准确地识别出无刚体运动。

在第五个情形中，光学面同时存在刚体运动和弹性变形，刚体运动和弹性变形测试值都设为0.001，如表1所示，所提方法比传统方法更为精确地计算出刚体运动量，传统方法求解的T_z与测试值有较大的偏差，$ {\theta _x} $、$ {\theta _y} $都有一定的偏差。

为了进一步验证所提方法的有效性，随机生成150个测试样本，在每个测试样本中，刚体运动的六个分量在0~0.01之间随机生成，用统计误差分析指标评估光学面刚体运动量的计算精度。常用的误差分析指标有确定性系数R²、均方根误差（RMSE）、相对均方根误差（RRMSE）、最大绝对误差（MAE）和相对最大绝对误差（RMAE）^[18]，其数学表达式如下：

式中：$ {l_i} $为测试值；$ {\hat l_i} $为求解值，$ {\bar l_i} $为测试值的均值；std为测试值的标准差；m为测试样本的数量。

确定性系数和相对均方根误差反映了光学面变形误差模型的整体精度，相对最大绝对误差则表征光学面变形误差模型的局部精度。统计结果如表4所示，所提方法的统计误差指标值反映出其计算精度较高，六个刚体运动分量的确定性系数都等于1，相对均方根误差都很小，且相对最大绝对误差的最大值为7.7869 E-6。传统方法的统计误差指标值则反映出其计算精度相对较低，T_z、$ {\theta _x} $、$ {\theta _y} $的确定性系数相对较小，同时它们的相对均方根误差和相对最大绝对误差都较大，远大于所提方法的统计误差指标值。

刚体运动方程和光学面修正方法是影响光学面变形误差计算精度的关键因素，在所提方法中，刚体运动的优化求解过程是“动态的”，未变形前的光学面点阵按照刚体平动、X轴、Y轴、Z轴旋转的顺序运动，假设其运动向量为${\boldsymbol{L}}{\text{ = }}{\left[ {\Delta {X_S}\;\Delta {Y_S}\;\Delta {Z_S}\;\Delta {\theta _X}\;\Delta {\theta _Y}\;\Delta {\theta _Z}} \right]^{\rm{T}}}$，所提方法的刚体运动矩阵为：

传统方法中的刚体运动求解过程是“静态的”，为了分析传统方法中的刚体运动方程求解方法对刚体运动计算误差的影响，用该方程替换所提方法中的完备方程，保留所提方法中的光学面修正方法，对比两种方法的计算精度。

采用统计误差分析指标评估刚体运动量的计算误差，重新随机生成150个测试样本，每个测试样本中刚体运动的六个分量在0~0.01之间随机生成。统计结果如表5所示。所提方法和传统方程求解方法的整体计算精度都很高，确定性系统几乎都等于1，但传统方程求解方法的相对均方根误差值较大，且远远大于完备方程求解方法的相对均方根误差。此外，在表征局部计算精度的相对最大绝对误差方面，所提方法远小于传统方程求解方法，反映出所提方法在计算精度、结果稳定性等方面都优于传统方法中的运动方程求解方法。

传统方法的计算结果如表4所示，将传统方法中的面形修正方法改为双立方样条插值法，其计算结果如表5所示。表4对应的150个测试样本与表5的150个测试样本不同，但所提方法的误差指标值几乎相等，这反映了所提方法的计算结果具有较好的稳定性，同时也表明表4与表5中的统计结果可用于对比评价面形修正方法。由表4和表5可知，在整体计算精度方面，表5中的确定性系数几乎都等于1，而表4中T_z、$ {\theta _x} $、$ {\theta _y} $的确定性系数相对较小，表5中的相对均方根误差也都小于表4中的值；同样，在局部精度方面，表5中的相对最大绝对误差都小于表4中对应的指标值，双立方样条插值法的计算精度明显优于切线法。

光学面点阵由镜面有限元单元的划分方式及尺寸确定，目前光学镜面有限元模型的构建主要以经验性知识为主。以球形镜面为例，考虑到圆的对称性及六面体网格的划分，传统方法常用辐射状的网格划分方式，其节点分布如图7所示。为了降低矩阵的奇异性，节点的数量通常是待求参数项的2倍及以上，以Fringe Zernike多项式为例，则至少需要74个节点。

Figure 7.  Radial distribution of nodes on the optical surface. (a) 97 nodes；(b) 261 nodes；(c) 521 nodes；(d) 2191 nodes

为了分析节点数目对计算结果的影响，节点数分别设为97、261、521、1093及2191，如图7所示。刚体运动六个分量的测试值分别设为0.001、0.002、0.003、0.001、0.002、0.003，弹性变形测试函数选为Fringe Zernike多项式中的高阶次函数：第37项函数，其系数k₃₇设为0.001，表达式为：

用拟合残差RMS与弹性变形RMS的比值表征拟合效果，其表达式为：

式中：${{\boldsymbol{C}}^ * }$为光学面弹性变形的多项式拟合系数最优解；${\left\| {{\boldsymbol{Z}_{\text{d}}}} \right\|_{\text{2}}}$为${\boldsymbol{Z}_{\text{d}}}$的2范数。

计算结果如表6所示，Z向平动量的计算误差相对较大，其他刚体运动分量和弹性变形系数的误差都很小，面误差比也都很小，但随着节点数的增大，Z向平动量的计算值与其测试值越来越接近，但仍然有较大的误差。节点分布均匀性决定了节点的相对权重，聚簇处的节点权重大，稀疏处的节点权重小，而辐射状节点分布具有“外疏中间密”的特点，体现了中间部分权重大而外围节点的权重小，表6中的Z向平动量的误差是由高阶次变形和节点不均匀分布导致的。从方程组求解的角度看，节点数量多且分布均匀有利于提高光学面变形误差的计算精度；从有限元求解的角度，有限元单元选择四边形或六面体，形状越规整，网格尺寸越小，越有利于提高仿真精度。

利用商业有限元软件Altair Hypermesh建立近似均匀分布、四边形单元的点阵，如图8所示，节点数设为接近表6中的数目，分别为97、258、520、1092及2193，刚体运动和弹性变形的测试值保持不变。计算结果如表7所示，与辐射状节点分布情况的分析结论相似，但整体上的计算误差，特别是Z向平动量的计算误差比辐射状节点分布情况小，可知节点的分布均匀性降低了计算结果的误差。

Figure 8.  Approximately uniform distribution of nodes on the optical surface. (a) 97 nodes；(b) 258 nodes；(c) 520 nodes；(d) 2193 nodes

Fringe Zernike多项式中的前三项（离焦、X轴和Y轴偏转量）与刚体运动量的分量相同，为了表示更为复杂的光学面弹性变形，将第4项至第37项选为测试函数，系数k₄~k₃₇都设为0.001，刚体运动的测试值保持不变。在节点近似均匀分布情况下，刚体运动分量中的最大误差和弹性变形系数中的最大误差与节点数的关系如图9所示。由图9（a）可知，由于弹性变形模式极为复杂，刚体运动量的计算误差都很大，当节点数为97时，最大误差值为153.55%，但其随着节点数目的增大而不断减小；由图9（b）可知，弹性变形多项式系数的计算误差都很小，且随着节点数目的增大而不断减小，最大误差值为1.4%。

Figure 9.  Relationship between estimation errors and number of nodes. (a) Rigid body motion；(b) Elastic deformation coefficient

在实际光机工程中，主要关注Fringe Zernike多项式中的前15项（4阶次函数），其具有与光学相对应的物理含义。将测试函数设为Fringe Zernike多项式中的第4项至第15项，系数k₄~k₁₅都设为0.001，刚体运动的测试值保持不变，结果如图9所示。当节点数为97时，刚体运动量的计算误差为50.73%，随着节点数目的增大，计算误差先减小后增大，当节点数为1466（对应的单元平均尺寸为0.06）时，计算误差的最小值为4.64%；弹性变形的多项式系数的计算误差都很小，随着节点数目的增大，计算误差先减小后增大，当节点数为1789时，计算误差的最小值为0.05%。在节点近似均匀划分方法中，当节点数较小（如97）时，节点分布均匀性较好，但随着节点数目的增大，节点分布均匀性会出现恶化，如图8所示。这说明在该面形复杂度下，建立节点数约为1466的近似均匀网格可以得到较优的计算结果。

此外，对比表1~表7及图9可知，随着光学面变形复杂性的增加，如测试函数分别选为第4项、第37项、第4~15项和第4~37项，当节点数较小时，刚体运动量和弹性变形系数的计算误差不断增大，当节点数很大时，刚体运动量和弹性变形系数的计算误差呈先增大后减小的趋势。

光学面可能是标准球面的一部分，如X范围为[−a, a] (0<a<1)，此时Fringe Zernike多项式彼此之间是非正交的。为了评估所提方法对部分球形光学面的适用性，依次对a取值0.1、0.2、···、0.9，节点数取1466×a的近似值，节点分布如图10所示。测试函数分别选为第4项、第37项和第4~15项，并且系数都设为0.001，刚体运动的测试值保持不变。

Figure 10.  Approximately uniform distribution of nodes on the optical surface with X: (a) [−0.1, 0.1]；(b) [−0.3, 0.3]；(c) [−0.5, 0.5]；(d) [−0.7, 0.7]

刚体运动分量中的最大误差和弹性变形系数中的最大误差与面形范围系数a的关系如图11所示，由图11(a)可知，当测试函数为单个高阶函数时（第37项），刚体运动的计算误差整体上最小且基本不随a值变化，单个低阶函数其次，第4~15项函数对应的计算误差最大，后两者随着a值的增大而减小。由图11(b)可知，当测试函数为单个函数时，即简单变形模式，弹性变形系数计算误差值都很小，而当测试函数为第4~15项函数时，计算误差值整体上较大且随着a值的增大而快速减小，当a大于0.4时，计算误差都很小。由此可知，当面形范围系数a值较小时，在光学变形复杂度较高的情况下，所提方法的计算精度较低，这是由于光学面仅是标准球形面的较小部分，复杂面变形中的一部分被算法识别为该光学面的刚体运动，此时，光学面变形的多项式系数的计算误差也较大，无法反映出真实的物理含义。需要说明的是，尽管刚体运动量和光学面变形多项式系数的计算误差较大，最终的光学面拟合精度很高，它们的面误差比值都很小。在实际工程中，光学镜面在XY面上的投影可能是非标准圆，但可以用近似于图10中的光学面对其计算结果进行评价。

Figure 11.  Relationship between estimation errors and the surface range coefficient. (a) Rigid body motion；(b) Elastic deformation

该光学系统是基于史瓦西尔德构型的偏视场离轴两反光学系统，由主、次两块反射镜镜构成，其光路结构如图12(a)所示。平行光入射至主镜，经凸面主镜发散后至凹面次镜，由次镜聚焦成像到像面处，两块反射镜均为离轴二次曲线非球面。系统全视场角8.2°，线视场使用，系统焦距103.50 mm，入瞳直径42.29 mm，F/#为2.4。

Figure 12.  An optical system with two mirrors. (a) Optical design; (b) Mechanical structure; (c) Finite element model

基于工程实践经验设计简化的镜头光机结构，如图12(b)所示。光机结构件外形为呈字母“U”形的板件结构，两头设计安装镜框，中间段连接固定两端镜框，并设计镜头安装接口。反射材料为熔石英，结构件材料为铝合金6061-T6，反射镜使用硅橡胶粘结安装在镜框内。光机三维尺寸为248 mm×133 mm×160 mm，质量为2.19 kg。

光学系统的有限元模型如图12(c)所示，主镜和次镜的单元平均尺寸分别为1.8 mm和3.6 mm，节点数分别为1457和1432，整个有限元模型共有152148个单元、34807个节点。系统工况载荷为：重力载荷和−30°温度载荷，边界条件为：完全约束支撑板四个螺栓孔的6个自由度。利用商业有限元软件MSC NASTRAN求解该有限元模型，并提取两个光学面的节点坐标和位移信息，如图13所示。

Figure 13.  Finite element and displacement of optical surfaces. (a) Finite element-primary mirror; (b) Displacement-primary mirror；(c) Finite element-secondary mirror; (d) Displacement-secondary mirror

计算结果如表8所示，主镜和次镜的面误差比分别为0.44%和1.82%，弹性变形的波峰波谷值(peak-valley, PV)分别为3.4480 E-6 mm和1.9449 E-5 mm，远小于刚体位移值，光学面弹性变形和拟合残差如图14所示，拟合残差远小于弹性变形量，说明拟合效果较好。

Figure 14.  Elastic deformation and fitting residual. (a) Primary mirror deformation; (b) Deformation fitting residual of primary mirror；(c) Secondary mirror deformation; (d) Deformation fitting residual of secondary mirror

以表8中的变形误差参数和Fringe Zernike多项式系数作为输入，利用商业光学软件ZEMAX对该光学系统的光学性能进行评价。像面离焦量为2.1658 E-2 mm，光学面变形前、后的点阵图和调制传递函数（modulation transfer function, MTF）如图15所示，可知光学面变形对光斑和MTF的影响都较小，说明光学系统在该工况下的光学性能变化较小。

Figure 15.  Spot diagram and MTF curve. (a) Original spot diagram; (b) Original MTF；(c) Spot diagram of the deformed system; (d) MTF of the deformed system

恶劣环境下的复杂光机系统的光学性能集成分析一直是光学仪器领域的难题。文中提出基于刚体运动完备方程的光机热集成分析方法，利用双立方样条插值法对光学镜面进行修正，然后通过刚体运动完备方程的数值动态优化过程分离出刚体运动，最后采用最小二乘法求解出弹性变形多项式系数。该方法的求解过程不依赖光学面的解析方程，仅需要光学面点阵的坐标和位移。数值案例和工程案例的分析结果表明，所提方法在计算精度、结果稳定性等方面都优于传统的刚体分离方法，特别是在刚体位移相对较大、变形模式复杂的情况下，所提方法仍具有很高的计算精度，从而具有更大的适用范围；通过定性和定量分析刚体运动方程、光学面修正方法、光学面节点分布以及面形等因素对光机热分析结果的影响，提供了最佳的光机热集成分析条件。