0 引言

现代科技的迅速发展对光电经纬仪提出了更多要求，比如在保证观测精度的同时还要具有良好的机动性.在光电经纬仪中，反射镜的面形对成像质量有着很大的影响.高机动性需要减轻反射镜的支撑结构质量，而观测精度则需要保证反射镜镜面面形.一般来说，支撑结构质量越大其刚度越大，反射镜的面形越好.因此，这两种需求是相互矛盾的，可以用结构优化的方法，在减轻支撑结构质量的情况下寻找合理的结构来保证较好的反射镜面形.结构优化大致可分为参数优化、形状优化和拓扑优化.其中，拓扑优化因其对初始的拓扑没有限定、可能得到新颖而性能优良的优化结构而具独特的优势.常用的拓扑优化方法有带惩罚的各向同性材料(Simplified Isotropic Material with Penalization, SIMP)方法^[1]、水平集方法^[2]，以及双向渐进结构优化(Bi-Evolutionary Structural Optimization, BESO)方法^[3]，本文选择的是SIMP方法.针对反射镜的支撑结构的拓扑优化，LIU Shu-tian^[4]、HU Rui^[5]等以镜面的均方根差(Root Mean Square，RMS)为目标；WEI Lei^[6]等以支撑结构刚度和点位移为目标；HUANG Hai^[7]等以结构惯性矩为目标；PARK K S^[8]等寻找系统的Strehl比和镜面的RMS之间的联系并且以RMS的函数近似Strehl比作为目标.这些工作大都以反射镜的RMS、支撑结构刚度等为目标，优化这些目标可以对成像质量起到一定的作用.然而，要获取更好的镜面面形，对以上的目标进行优化是不够的.原因在于，RMS反映的是反射镜的光学表面的整体性质，没有局部的面形信息，不能确定具体的位移模式；而支撑结构刚度和成像质量之间没有明确的函数关系，刚度极大的结构不一定是成像质量最佳时对应的支撑结构.所以，要获取更好的面形，必须要对反射镜的面形进行更细致、更精准的描述.

本文选择Zernike多项式^[9]来描述变形的反射镜镜面.在光学设计中，Zernike多项式广泛应用于拟合曲面^[10]或波阵面^[11]来表达面形误差或者像差.相较于RMS，Zernike多项式能够描述镜面的局部形状，更加的精细、准确.以Zernike多项式的系数为目标的拓扑优化可以针对反射镜面形中特定的像差成分进行优化，能得到凭借现有的设计流程、工程经验所不能得到的新颖的、更优的结构，从而减小最终成像的像差.目前，可查找到的多局限于结构参数优化，其设计变量的个数相较于拓扑优化(一般在几千个设计变量以上)较少，直接以Zernike系数为优化目标或者约束的拓扑优化工作尚未见报导.在前述文献中，Zernike系数也仅仅被用来验证优化结构的有效性，而并没有作为目标或者约束参与到优化过程中.事实上，直接应用Zernike系数为目标或者约束进行优化有着很大的困难，依本文分析主要有两点：一方面，Zernike系数的数值计算较繁杂，通常需要在一个数据传递环境(如Sigfit)中集成光学计算软件(如CODEV、Zemax等)或其他第三方数据处理软件(如MATLAB)和有限元软件(如ANSYS)，先在有限元软件中计算需要拟合的实际面形的数据，再将面形数据传递到光学计算软件或第三方数据处理软件中计算Zernike系数，例如SAHU R^[12]等使用提取的有限元计算结果，借助MATLAB计算Zernike多项式系数.这种分步计算的方式流程复杂，而且在光学软件中需要对面形数据进行密集插值以保证精度，计算量很大.另一方面，将Zernike系数加入到拓扑优化的目标或者约束中，均需要求解Zernike系数对设计变量的敏度，而这是拓扑优化中较为关键的一个步骤.在上述的计算流程下，Zernike系数的敏度只能用差分的方式计算.但是，拓扑优化的设计变量和节点数相当，动辄数十万；如果采用差分的方式计算，需要以多软件协同的形式循环进行Zernike系数的求解，其计算量、计算时间显然是不可接受的.另外，商用有限元程序发展迅速，对于一些目标(如应变能、指定点位移)的敏度计算可以方便的在程序内获得.由于这些原因，在拓扑优化问题的建模过程中，含Zernike系数的目标或者约束常常会被一些易于在商用软件内求解敏度的物理量所等效替代^[12].然而，部分较低阶的Zernike多项式表征的位移模式相对简单，等效替代可以简单、有效的进行；其余较高阶的Zernike多项式表征的位移模式太过复杂，等效替代很难进行.为了对任意拟合的Zernike系数直接进行优化或者约束，寻找其对设计变量的敏度的计算方法是必要的.针对Zernike系数的计算，本文结合Zernike最小二乘拟合的算法以及有限元离散方法，实现了在有限元框架下快速求解Zernike系数，克服了数据传递以及密集再次插值的问题，化简计算流程的同时还能保证计算精度.更为重要的是，这样计算可以建立起Zernike系数和物理场的解之间明确的函数关系，这是以界面接口耦合形式的多软件协同计算的方式很难做到的，而这些函数关系又是求解Zernike系数敏度的重要基础.针对Zernike系数敏度的计算，本文选择较为成熟的SIMP方法，以已经建立的Zernike系数和物理场解的函数关系为基础，直接以Zernike系数为目标、约束，写出协调的拓扑优化数学模型，使用离散伴随的方法推导目标函数对设计变量的敏度.这种方式可以一次求解所有设计变量对应的敏度，尤其适合设计变量非常多的拓扑优化.本文实现了直接对Zernike系数进行优化或者约束.理论上，本文的算法可以对任意多项式拟合的Zernike系数的解析函数进行优化或者将其作为约束，具有一定的泛用性.此外，为了证明本文算法能够有效地优化Zernike系数，使用COMSOL二次开发编写程序，对水平放置和45°斜放双工况的反射镜镜面支撑进行了拓扑优化.

1 结构有限元分析与光学表面Zernike系数计算

1.1 反射镜模型及有限元计算

本文研究的模型及尺寸如图 1所示.简便起见，本文仅以平面镜为例进行优化理论方法的研究，但整个方法可以推广至任意形状的反射镜的支撑结构设计.整个模型分为反射镜以及其支撑结构两个区域，两者的半径均为R，反射镜区域厚度为t₁，支撑结构区域厚度为t₂.模型的载荷以及边界条件如图 2所示.为了便于描述，图 2(a)将模型沿着直径剖开，图 2(b)仅展示模型的支撑结构区域.反射镜区域记为Ω₁，受重力G作用，是拓扑优化中的非设计区域；支撑结构区域记为Ω₂，是设计区域，和底部三个固定点相连，并且固定点成120°圆周对称排列.在反射镜和支撑结构的交界面Γ上，本文假设位移场连续即u_Ω1=u_Ω2.这样，有限元计算仅考虑材质分界面带来的影响.

图 1. 研究模型及尺寸

Fig. 1. Study model and size

下载图片 查看所有图片

图 2. 载荷、边界条件与设计、非设计区域

Fig. 2. Load, boundary condition, design andnondsign domain

下载图片 查看所有图片

从图 1可以清晰地看到，整体结构呈偏平化，其横向尺度相较于厚度方向的尺度大了将近一个量级.特别是反射镜部分，相较于支撑结构更薄，结构的扁平化更严重，整体结构的行为可能会表现出部分厚板的特性.这种几何特征给高精度的精确有限元计算带来了困难.如果使用的网格较粗，往往不能准确地计算反射镜和支撑结构的变形，进而会影响镜面Zernike系数的计算.对此，本文选择使用六面体映射网格和三棱柱网格相结合，使用扫掠网格划分反射镜和支撑结构，通过增加映射网格、增加扫掠层数加密单元，三分之一的网格如图 3所示.这种网格划分可以使用更密的实体单元来比较精确地描述镜体和支撑结构的行为.此外，整体结构和支撑呈现出明显的120°圆周对称，因此网格的划分也呈现出同样的对称形式以提高计算精度.

图 3. 三分之一模型的网格

Fig. 3. Mesh of one third model

下载图片 查看所有图片

1.2 Zernike多项式及系数的计算

Zernike多项式是定义在单位圆上的一系列正交的多项式，利用这些多项式可以将镜面的位移表述成对应多项式所描述的位移模型的加和形式.同时，Zernike多项式和光学像差有着紧密的联系，特别是赛德尔初级五像差，可以直接用Zernike多项式对应的系数来描述.

根据资料^[1]，给定的表面ΔZ(r, θ)可以表述为

式中，θ为极角，r为极径，A_nm、B_nm为Zernike系数，R_n^m为Zernike多项式的生成表达式，其定义为

式(1)、(2)中，n和m分别表示径向和周向的波数，n-m必须是偶数且n≥m.以上的Zernike多项式定义在极坐标下，在本文中，θ、r的定义^[13]如图 4所示.

图 4. 角度θ的定义

Fig. 4. Definition of angle θ

下载图片 查看所有图片

Zernike系数的计算可以使用最小二乘拟合来实现.定义E为Zernike拟合后的面形与实际面形之间的拟合误差，有

式中，p为节点数，δ_i为实际面形在节点i处的位移，Z_i为Zernike拟合后的面形在节点i处的位移，W_i为节点i的权重.Z_i可以用Zernike多项式和系数的线性组合表述为

式中，q为Zernike拟合使用的项数，ϕ_ji为节点i处Zernike第j项的值，c_j为需要求解的Zernike系数.式(4)带入式(3)有

为了计算Zernike系数c_j，需使误差函数E取极小.在这里，可令E对Zernike系数的偏导数为零求解获得，即

这是一个线性代数方程组，可以表述为

其中，

1.3 有限元框架下的Zernike系数求解

由1.2节可以看到，Zernike多项式拟合是离散的，是基于镜面网格节点位移数据的，拟合的关键是计算H_jk和P_j，且有

式(10)、(11)中，W_i是权重，其意义是节点i在曲面上所占的面积.在实际工程应用中，镜面的变形相较于结构的尺寸是很小的.如果将变形后的曲面位移函数δ看作坐标的函数，δ_i为其在各个节点上的值，那么这个函数在曲面上是有限的.所以，可以认为W_i是曲面的一个分割，H_jk、P_j实际上是ϕ_jϕ_k、δϕ_j在曲面上的、以W_i为分割的黎曼积分^[14].分割W_i越精细，黎曼积分就越准确逼近真实的函数积分值；在较粗的网格中，单元内的δ、ϕ_j变化较大，用常数去近似会产生较大的误差，这也是需要密集插值计算Zernike系数的原因.如果节点个数趋向于∞，亦即W_i→0时，可以得到

式(12)、(13)中，两式的积分在曲面上进行，其中，积分号下的ϕ_i为Zernike多项式第i项，是解析而又连续的，给定曲面上的任意一点即可求出在该点的函数值.δ是曲面上节点的位移函数，应当是连续的，但δ是隐式的；有限元计算求解的只是δ在各个节点上的值δ_i，要想计算在任意一点的函数值需要进行单元内插值.

为了在有限元框架内求解Zernike系数，可以在曲面上进行ϕ_jϕ_k、δϕ_j的数值积分来计算∫ϕ_jϕ_kdΓ、∫δϕ_jdΓ，这种积分是使用高斯求积法^[15]来进行的.由于可以得到ϕ_i在各个节点上的准确数值，进而使用单元内插值的方式得到高斯积分点上的函数值；δ_i即是节点位移，也可以用形函数插值得到高斯积分点上的位移数值.这样，有

式(14)、(15)中，W_i^′、ϕ_ji^′、δ_i^′分别表示在高斯积分点i上的高斯积分权重、Zernike多项式第j项在高斯积分点i上的数值以及高斯积分点i上的位移.具有M个高斯积分点的高斯求解算法可以精确积分的多项式最高阶次为2M-1，这样可以用较少的积分点来得到较高的精度.

由于ϕ_i、δ在曲面上都是有限、连续的，二者在曲面上的积分也是有限、连续的，所以对于C₀连续函数，在同样比较逼近真实积分结果的情况下黎曼求积和高斯求积的数值积分结果是等价的.由此，得到以下定义

式(16)、(17)都可以在有限元计算中直接实现.这样，通过式(16)、(17)将Zernike系数的计算和有限元计算联系了起来，建立起Zernike系数和有限元计算结果之间明确的离散联系；此外，Zernike系数的计算可以直接在有限元数值计算框架中实现；而且，积分的计算在高斯积分点上进行，不再需要密集的插值，Zernike系数的计算网格和有限元网格重合，计算精确，计算效率也比较高.

考虑到高斯积分点上的函数值可以由单元节点插值表示，即

式(18)、(19)中，n_e为单元节点数目，N_i为单元节点i的形函数，δ_ei、ϕ_eji分别为δ、ϕ_j在网格节点上的值.由形函数的性质

有

式(21)、(22)中，ϕ_eji为Zernike多项式第j项在网格节点i上的值；δ_es为网格节点i上的位移；W_ij^′为包含网格节点i的单元中第j个高斯积分点的高斯积分权重；N_j为该高斯积分点的插值函数；m为网格节点数目；o为包含网格节点i的单元中的高斯积分点的数目.令

表示节点i的积分权重，那么有

这样，在高斯积分点上进行的计算转化为在网格节点上的计算，这和有限元计算在节点上进行是相吻合的，便于计算和推导.

2 基于Zernike系数的拓扑优化

2.1 拓扑优化模型

本文研究的优化模型如图 2所示.反射镜(Ω₁)为非设计域，支撑结构(Ω₂)为设计域.本文采用SIMP法，在设计域中基于网格节点进行设计变量的惩罚插值.在整个模型中，整体结构仅受一个弹性应力场作用，其载荷以及边界如1.1节中所述.

优化的数学模型为

上述数学模型中，设计变量为ρ，目标是Zernike系数C的函数；式(26)如1.2节中所述，其作用为计算Zernike系数；式(27)为弹性应力场离散的有限元方程；式(28)为各项Zernike系数的约束；式(29)为体积约束，V^*为体积分数上限，通常由实际需求确定；式(30)为设计变量取值约束.

2.2 敏度求解

为了求解2.1节中的优化问题，必须求得目标对设计变量的敏度.本文采用的是伴随敏度求解方法，其推导如下.

记目标为

在求解优化敏度问题前，式(26)、(27)均已经求解完毕，那么有

式(31)中，λ₁、λ₂分别为Zernike求解方程、弹性应力场方程的伴随变量，其维度分别和C，U对应.于是，目标对设计变量ρ的敏度为

为了消去场变量对设计变量的敏度，令含$\frac{\mathrm{d} \boldsymbol{C}}{\mathrm{d} \rho}$、$\frac{\mathrm{d} \boldsymbol{U}}{\mathrm{d} \rho}$的项为零得到伴随方程

求解式(33)、(34)得到伴随变量λ₁、λ₂，目标函数的敏度表达式化简为

2.3 数值实现

在2.1节中所述的优化问题是非线性的，求解该问题需要使用迭代的方式进行，在每个迭代步中求解近似的线性化子问题.具有单个约束(体积约束)的优化问题可以用很多算法求解，比如优化准则法(Optimality Criteria，OC)、进化算法等.然而，本文的约束更多，使用上述算法可能得不到一个满意的解.目标复杂、具有多个约束的优化问题可以由线性序列规划方法(Sequential Linear Programming, SLP)、CONLIN以及移动渐近线法(Method of Moving Asymptotes，MMA)等算法求解.本文使用MMA^[16]算法，这种算法将非凸的目标和约束作了凸近似，在每一个迭代步内求解一个线性化的子问题，用子问题的解近似原问题的解.

在实际的数值实现时，首先指定设计域的材料、边界条件、载荷等进行有限元计算初始化.每一个优化迭代步由以下三个部分组成：1)求解弹性应力场方程、Zernike系数方程；2)敏度分析，求解伴随变量、敏度；3)用MMA更新设计变量.对于有多个工况的优化问题，步骤1)、2)需要多次重复.当满足预设的收敛条件或达到最大迭代步时，迭代停止，获得优化的结构.

本文算法的实现使用了有限元软件COMSOL的二次开发功能.

3 算例及结果

3.1 模型设置

算例使用的反射镜、支撑几何尺寸如图 1所示，R=0.5 m，t₁=0.1 m，t₂=0.2 m.反射镜的材料为微晶玻璃，其杨氏模量为90.3 GPa，密度为2 530 kg/m³，泊松比为0.243；支撑结构的材料为超殷钢4J32，其杨氏模量为141 GPa，密度为8 100 kg/m³，泊松比为0.25.本文只考虑反射镜受重力G作用，忽略支撑结构的重力，固定边界如前所述.Zernike拟合项取Fringe列表前11项.同时，选择水平放置和倾斜放置两种情况下的算例来说明算法的有效性.

3.2 水平放置

反射镜水平放置，模型的载荷及边界条件如图 2所示.水平放置下，反射镜的面形Zernike系数中，轴向位移、场曲、初级球差和初级三叶像差占主要的部分，都需要优化，所以这是一个多目标的优化问题.在这里选择轴向位移为目标，其余的像差作为约束，得到如下的拓扑优化数学模型

式中，相关方程的含义同节2.1；c₄^*、c₉^*和c₁₀₁₁^*为约束的上限，可以根据实际的工程需求选取一个相对合理的数值，本文取10^-16；体积分数上限V^*取0.2；初始结构为设计变量0.3的均匀分布结构.

支撑结构的拓扑优化结果如图 5所示, 这里没有绘出支撑结构上面的反射镜，为完整的优化结构提供了两种视角.由于优化模型的固定边界条件(如图 5(a)中棕色圆锥形所示)和边界载荷具有三分之一圆周对称性，所以其拓扑优化结构也具有相应的对称性.从图 5(a)中看到，结构中存在较为细小的结构，如灰色圆圈、红色圆圈和紫色圆圈部分.这些细小结构的上表面和比较粗壮的六个支撑结构的上表面在一个平面上，和反射镜相连，也就是说细小结构也是实际支撑反射镜的.由于本模型是多目标优化，优化的约束中存在场曲、初级球差这类中心为位移模式曲面驻点的像散项，所以存在灰色圆圈、紫色圆圈内的结构来减小中心部分位移；圆周散布的六个较为粗壮的支撑结构结合红色圆圈内的结构可以减小三叶像散的烈度.图 5(b)、(c)分别展示了拓扑优化结构的俯视图和仰视图.

图 5. 支撑结构的完整优化结果

Fig. 5. Full optimal of supporting structure

下载图片 查看所有图片

图 6为初始结构与拓扑优化结构对应的合并后的Zernike系数，各项定义为

图 6. 各项系数优化前后对比

Fig. 6. Comparison of coefficients between before and after optimization

下载图片 查看所有图片

coef.1表征镜面的轴向位移位移模式；

coef.2表征镜面的倾斜位移模式；

coef.3表征镜面的初级慧差位移模式；

coef.4表征镜面的像散位移模式；

coef.5表征镜面的离焦位移模式；

coef.6表征镜面的初级球差位移模式；

coef.7表征镜面的三叶像散位移模式.

为了更加直观，图 6中纵坐标的函数值设为lgcoef.i+12，i=1, 2, …, 7.在图 6中，约束值对应的函数值为4，如红色虚线所示.图 6中带有黑色框的柱状图表示优化的目标与约束对应的数值.可以看到，轴向位移得到了较大的优化，达到了10^-6量级；同时，场曲、初级球差以及三叶像散均得到了很好的优化，基本都满足了约束.目标对迭代步的变化如图 7所示.优化结构在优化循环150步时得到，在150步时基本已经收敛.

图 7. c₁²随循环步的变化

Fig. 7. Changes of c₁² along with iteration steps

下载图片 查看所有图片

3.3 45度斜放

反射镜45°斜放，模型的载荷如图 8所示，固定边界条件仍然如图 2(b)所示.考虑到实际工程应用和模型的对称性，这里包含两个工况，且两个工况中的载荷是对称的.简洁起见，这里固定坐标轴以及系统结构，用旋转重力的方向来表示不同角度下的受力状态.如图 8所示，将模型沿着垂直于对称平面的直径剖开，图 8(a)和(b)分别展示了两种工况中的重力(G₁和G₂)方向.

图 8. 45°斜放情况下的载荷

Fig. 8. Loads in tilts of 45°

下载图片 查看所有图片

这种放置方式下，反射镜的面形变化中，轴向位移、倾斜、场曲、像散、初级慧差、初级球差以及初级三叶像散都较大，均需要优化，这也是一个多目标优化问题.这里选择轴向位移和倾斜项作为目标，其余项作为约束小于给定值的形式进行优化，得到的拓扑优化数学模型为

式中，c₄^*、c₅₆^*、c₇₈^*、c₉^*和c₁₀₁₁^*为约束的上限，可以根据实际的工程需求确定.在这里，可以假定一个相对合理的数值，比如都取10^-16.体积分数上限V^*取0.2，初始结构为设计变量0.2的均匀分布结构.

支撑结构的拓扑优化结果如图 9所示，除了图 9(b)、(c)以外，其余子图均未画出反射镜.图 9(a)为拓扑优化结构的正等轴测图，图中的蓝色圆锥表示固定点.可以明显地看到，拓扑优化结构有着明显的对称性，对称面如图 9(a)所示.这种结构的对称性是由两种工况中载荷的对称性造成的.图 9(b)、(c)为优化结构的左视图，同时还画出了反射镜，二者的视角如图 9(a)中粉色箭头所示，两种工况下的重力载荷如图中所示意并且和图 8相对应可以很明显的看到载荷的对称性.图 9(d)、(e)的3D视图视角分别对应图 9(a)中紫色和蓝色箭头.图 9(e)中蓝色的圆锥同样表示固定点.图 9(f)、(g)为优化结果的俯视图和仰视图.

图 9. 支撑结构的优化结果

Fig. 9. Optimized solution of supporting structure

下载图片 查看所有图片

图 10列出了初始结构与拓扑优化结构对应的合并后的Zernike系数，各项定义为

图 10. 各项系数优化前后对比

Fig. 10. Comparison of coefficients between before and after optimization

下载图片 查看所有图片

coef.1表征镜面的轴向位移和倾斜的“刚体位移”位移模式；

coef.2表征镜面的离焦位移模式；

coef.3表征镜面的像散位移模式；

coef.4表征镜面的初级慧差位移模式；

coef.5表征镜面的初级球差位移模式；

coef.6表征镜面的三叶像散位移模式.

为了更加直观，图 10中纵坐标的函数值设为lgcoef.i+12，i=1, 2, …, 6.在图 10中，约束值对应的函数值为4，如红色虚线所示.图 10中除了coef.1以外，其余所有的系数都是优化对象，可以明显的看到除了coef.1，其余所有的系数都减小到了10^-8左右的量级，基本满足了约束；coef.1本身也达到了10^-6量级.据此，本文的优化算法在这种工况下也适用.目标对迭代步的变化如图 11所示.优化结构在优化循环200步时得到，在200步时基本已经收敛.

图 11. c₁²+c₂²+c₃²随循环步的变化

Fig. 11. Changes of c₁²+c₂²+c₃² along with iteration steps

下载图片 查看所有图片

4 结论

本文将Zernike系数的计算结合到有限元框架中，建立了Zernike系数和有限元位移场的解之间的函数关系，实现了直接以Zernike系数为目标、约束的多目标反射镜支撑结构的拓扑优化.两个示意算例结果说明了本文算法的有效性.在本文理论的基础上，后续将通过增材制造等形式将优化结构制造出来，与理论结果相对照；同时可以考虑支撑结构的自重进行优化，并开发优化波像差的算法，即以整个系统在像空间内波面的Zernike系数为目标，对反射镜的支撑结构进行优化，从而实现系统级的拓扑优化.