天文光学望远镜主镜轴支撑的一种变网格复合优化法

摘要

本发明涉及天文光学望远镜主镜轴支撑优化法，尤其涉及天文光学望远镜主镜轴支撑的一种变网格复合优化法。本发明要解决的技术问题是现有的天文光学望远镜主镜轴支撑优化法速度过慢，无法满足高效快速计算的需要。为解决上述技术问题，本发明采用的技术方案是：天文光学望远镜主镜轴支撑的一种变网格复合优化法，本发明的优点是：其通过人为预先将镜面分为密网格区和疏网格区，再通过零阶优化法检验划分的正确与否，对正确划分的密网格区单独进行一阶优化法计算，快速得出高精度的支撑点位置，按这种方法，在中等性能PC(个人计算机)下，主镜支撑的每一种设计在一个小时内得到一个最佳方案。

说明书

技术领域

本发明涉及天文光学望远镜主镜轴支撑优化法，尤其涉及天文光学望远镜主镜轴支撑的一种变网格复合优化法。

背景技术

天文望远镜(AstronomicalTelescope)是观测天体的重要工具，可以毫不夸张地说，没有望远镜的诞生和发展，就没有现代天文学。随着望远镜在各方面性能的改进和提高，天文学也正经历着巨大的飞跃，迅速推进着人类对宇宙的认识。

由于望远镜的集光能力随着口径的增大而增强，望远镜的集光能力越强，就能够看到更暗更远的天体，因此，天体物理的发展需要更大口径的望远镜。

但是，随着望远镜口径的增大，一系列的技术问题接踵而来。一方面，望远镜的自重过大会使镜头变形相当明显，另一方面，镜体温度不均也令镜面产生畸变，进而影响成像质量。

如何设计合理的主镜轴支撑以保证在该支撑下，镜头的变形能够尽量小，现有的方法是先根据镜面的口径大小人工设定主镜轴支撑的目标函数，设计变量与状态变量；通过纯数学推导的方式计算出大概的支撑位置所在，然后，利用一阶优化法对整个镜面经行优化，得到具体的支撑位置所在，观测优化得到的各个支撑点和纯数学方式计算得出的各个支撑点之间的距离是否在合理的误差范围内，如果都在，优化结束。采用一阶优化法的优点是精度非常高，然而，其缺点也很明显，运算速度太慢，按现有的方法，在使用高性能计算机的条件下，主镜支撑的每一种设计在一个月内得到一个最佳方案都属于较快的，而主镜轴支撑的设计方案一般都不止一种。此外，其还有一个缺点就是得到的最佳值是局部最优值。

发明内容

本发明要解决的技术问题是现有的天文光学望远镜主镜轴支撑优化法速度过慢，无法满足高效快速计算的需要。

为解决上述技术问题，本发明采用的技术方案是：天文光学望远镜主镜轴支撑的一种变网格复合优化法，包括以下步骤：1)对主镜面建立目标函数，该目标函数对每种主镜轴支撑设计都返回一个目标值，根据目标值所在的位置将主镜面划分成疏网格区和密网格区，目标值所在位置为密网格区，其余位置为疏网格区；2)以步骤1中得到的疏网格步长为基准设置零阶优化法的容差，使得容差等于疏网格步长，然后，对主镜支撑进行零阶优化，得到一个低精度全局值；3)对步骤2中得到的优化结果中支撑点位置返回值进行判断，判断各个支撑点位置是否落在步骤1中划分的密网格区域内；如果是，确定该密网格划分正确，进行步骤4；如果否，则重复步骤1，对网格进行重新划分设置，直至步骤2中得到的优化支撑点位置均落在密网格区域内，确定该密网格划分正确；4)以步骤3得到的划分正确的密网格步长为基准设置一阶优化法的容差，使得容差等于密网格步长，然后，对密网格区进行局部的一阶优化，得到高精度支撑点值。

采用零阶优化法进行优化，其优点是速度非常快，然而，缺点更明显，其精度甚至低于单纯通过算法求得的目标函数的目标值，因此，在现有的优化法中不会采用零阶优化法进行优化工作，然而，在本发明方法中，先在电脑中通过算法预先设定目标值，将整个主镜面划分疏网格区和密网格区，再采用零阶优化法得到全局最优值，对密网格区的划分经行校验，是的密网格区的位置能够确定，再通过高精度的一阶优化法仅仅正对各个密网格进行局部优化，获得最优值，这样设置密网格区域得到精度相当整个镜子都是密网格情况，而速度却比整个镜子都是密网格快。因为网格越密，网格数量就越多，要解的方程就越多，速度就越慢。以上所有步骤的实现都是基于大型有限元软件ANSYS平台与APDL语言程序编程实现。

本发明的优点是：其通过人为预先将镜面分为密网格区和疏网格区，再通过零阶优化法检验划分的正确与否，对正确划分的密网格区单独进行一阶优化法计算，快速得出高精度的支撑点位置，按这种方法，在中等性能PC(个人计算机)下，主镜支撑的每一种设计在一个小时内得到一个最佳方案。

附图说明

图1是WFST主镜剖面示意图。

图2是WFST主镜支撑点分布示意图。

图3是WFST主镜参数化模型。

图4是WFST主镜1/72之一网格划分图。

图5是WFST整个主镜网格划分图。

图6是WFST主镜轴支撑载荷施加示意图。

图7是支撑点为N_P＝27时目标函数RMSe与设计变量R1,R2和R3关系四维图。

图8是支撑点为N_P＝27且目标函数RMSe最佳值28.78nm时主镜反射面相对于原坐标系形变图。

图9是支撑点为N_P＝27且目标函数RMSe最佳值28.78nm时最小半光程误差分布图。

图10是支撑点为N_P＝39时目标函数RMSe与设计变量R1,R2和R3关系四维图。

图11是支撑点为N_P＝39且目标函数RMSe最佳值9.32nm时主镜反射面相对于原坐标系形变图。

图12是支撑点为N_P＝39且目标函数RMSe最佳值9.32nm时最小半光程误差分布图。

图13是支撑点为N_P＝54时目标函数RMSe与设计变量R1,R2和R3关系四维图。

图14是支撑点为N_P＝54且目标函数RMSe最佳值5.29nm时主镜反射面相对于原坐标系形变图。

图15是支撑点为N_P＝54且目标函数RMSe最佳值5.29nm时最小半光程误差分布图。

具体实施方式

本发明包括以下步骤：

1)对主镜面建立目标函数，该目标函数Π对每种主镜轴支撑设计都返回一个目标值，根据目标值所在的位置将主镜面划分成疏网格区和密网格区，目标值所在位置为密网格区，其余位置为疏网格区；

2)以步骤1中得到的疏网格步长为基准设置零阶优化法的容差，使得容差等于疏网格步长，然后，对主镜支撑进行零阶优化，得到一个低精度全局值；

3)对步骤2中得到的优化结果中支撑点位置返回值进行判断，判断各个支撑点位置是否落在步骤1中划分的密网格区域内；如果是，确定该密网格划分正确，进行步骤4；如果否，则重复步骤1，对网格进行重新划分设置，直至步骤2中得到的优化支撑点位置均落在密网格区域内，确定该密网格划分正确；

4)以步骤3得到的划分正确的密网格步长为基准设置一阶优化法的容差，使得容差等于密网格步长，然后，对密网格区进行局部的一阶优化，得到高精度支撑点值。

步骤1的优化法数学原理如下

一般结构优化问题数学公式描述如下：

$\begin{array}{cc}\hat{\Pi }=\min \Pi \left(\gamma \right)& \gamma =[{\gamma }_1,{\gamma }_2,L,{\gamma }_m]\\ s.t.& {\gamma }_i\le {\gamma }_i\le {\gamma }_i,(i=1,2,3,L,m)\\ & g_j\le g_j\left(\gamma \right)\le g_j,(j=1,2,3,L,n_1)\\ & h_j\le h_j\left(\gamma \right)\le h_j,(j=1,2,3,L,n_2)\\ & w_j\le w_j\left(\gamma \right)\le w_j,(j=1,2,3,L,n_3)\end{array}---\left(1\right)$

目标函数(Π):对每种可能的设计，(Π)都返回一个目标值，一般通过各种优化法寻求(Π)最小值。设计变量(γ):一个函数或矢量数组，目标函数随着设计变量变化而变化。状态变量(g,h,w):对给定结构，每一组设计变量对应一组状态变量。

步骤2中涉及的零阶优化法数学原理如下：

利用最小二乘法拟合，可以把目标函数写成以下形式:

$\hat{\Pi }=a_0+\underset{i}{\overset{m}{\Sigma }}{a}_i{\gamma }_i+\underset{i}{\overset{m}{\Sigma }}\underset{j}{\overset{m}{\Sigma }}{b}_{ij}{\gamma }_i{\gamma }_j---\left(2\right)$

通过罚函数，把设计变量和状态变量约束问题转化成无约束问题，可以得到：

$F\left(\gamma ,p_k\right)=\hat{\Pi }+{\Pi }_0{p}_k[\underset{i=1}{\overset{m}{\Sigma }}\gamma \left({\gamma }_i\right)+\underset{i=1}{\overset{n_1}{\Sigma }}G\left(g_i\right)+\underset{i=1}{\overset{n_2}{\Sigma }}H\left(h_i\right)+\underset{i=1}{\overset{n_3}{\Sigma }}W\left(w_i\right)]---\left(3\right)$

这里γ_i是设计变量,g_i,h_i,w_i是状态变量,Υ,G,H,W分别是它们的罚函数。Π₀是目标参考值，p_k表示响应面参数。当设计变化或状态变量接近它们的约束值时候，惩罚函数值急剧增大。

步骤4中涉及的一阶优化法数学原理如下：

一阶优化法无约束问题方程如下：

$Q(\gamma ,q)=\Pi /{\Pi }_0+\underset{i=1}{\overset{m}{\Sigma }}{P}_{\gamma }\left({\gamma }_i\right)+q[\underset{i=1}{\overset{n_1}{\Sigma }}{P}_g\left(g_i\right)+\underset{i=1}{\overset{n_2}{\Sigma }}{P}_h\left(h_i\right)+\underset{i=1}{\overset{n_3}{\Sigma }}{P}_w\left(w_i\right)]---\left(4\right)$

这里Q(γ,q)是无量纲无约束目标函数；P_γ,P_g,P_h,P_w分别为设计变量和状态变量罚函数。Π₀是目标函数Π参考值.在整个设计空间内，对目标函数和状态变量罚函数进行微分。对每一个迭代(j)，引入优化搜索方向d⁽ⁱ⁾..则下一步(j+1)的设计变量如式(5)所示.在该方程中S_j为线搜索参数，对应于搜索方向d^(j)上最小的Q值。

γ^(j+1)＝γ^(j)+S_jd^(j)(5)

当第j步迭代目标函数和j+1步迭代和最佳值(b)满足一些步骤方程(6),则收敛。这里τ为目标函数容差。

|Π^(j)-Π^(j-1)|≤τand|Π^(j)-Π^(b)|≤τ(6)

主镜最小半光程误差数学描述

对于小变形弹性体，各应变分量与位移本构关系如下：

$\begin{array}{c}{ϵ}_{xx}=\frac{\partial u}{\partial x};{ϵ}_{xy}={ϵ}_{yx}=\frac{1}{2}\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\\ {ϵ}_{yy}=\frac{\partial v}{\partial y};{ϵ}_{yz}={ϵ}_{zy}=\frac{1}{2}\left(\frac{\partial w}{\partial y}+\frac{\partial v}{\partial z}\right)\\ {ϵ}_{zz}=\frac{\partial w}{\partial z};{ϵ}_{xz}={ϵ}_{zx}=\frac{1}{2}\left(\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}\right)\end{array}---\left(7\right)$

这里以反射面为抛物面主镜为例，假如主镜在没有重力时，反射面具有理想抛物面形状：

X²+Y²＝4f(Z+c).(8)

这里f为焦距，c为顶点。在重力的作用下主镜结构将发生变形，反射面表面节点将偏离原来抛物面表面的位置。假设这时存在一个反射面的最佳吻合抛物面，设这个最佳的抛物面在新的坐标系上表示为：

$X_1^2+Y_1^2=4f_1(Z_1+c).$

考虑抛物面面上的任一点i,则通过该点并垂直于抛物面表面的法线的方向余弦为：

$2{\mathrm{cos\alpha }}_1=-\frac{X_i}{\sqrt{f(f+Z_i+c)}}$

$2{\mathrm{cos\alpha }}_2=-\frac{Y_i}{\sqrt{f(f+Z_i+c)}}$

$2{\mathrm{cos\alpha }}_3=-\frac{2f}{\sqrt{f(f+Z_i+c)}}$

假设在该点在外力下的位移是(u_i,v_i,w_i)，变形后的该点与最佳抛物面的距离为Δ_i，则有

X-(X_i+u_i)＝±Δ_icosα₁

Y-(Y_i+v_i)＝±Δ_icosα₂

Z-(Z_i+w_i)＝±Δ_icosα₃

对于N个节点，变形表面对于理想抛物面的均方根距离偏差为

${\mathrm{RMS}}_g=\sqrt{\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}{\Delta }_i^2}$

上面所得均方根误差并不能代表表面有效偏差，真正表面均方根误差应该是通过优化最小半光程得到。均方根最小半光程误差如下：

${\mathrm{RMS}}_e=\sqrt{\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}{e}_i^2}$

e_i＝Δ_icosβ_i

这里β是表面法线和抛物面轴线的夹角。

实施例1

以2.5米大视场巡天望远镜(WFST)主镜轴支撑为例子

目标函数，设计变量与状态变量的设定

如图1所示，WFST观测的波长范围为320nm到1000nm,所以在重力载荷下轴支撑设计均方根最小半光程误差不能超过10nm。WFST主镜剖面是新月形，上表面和下表面都是焦距f＝2.13抛物面，主镜厚度h＝120.00mm,内径φ₁＝1000.00mm，外径φ₂＝2500.00mm。

如图2所示，WFST主镜一共由N_P力线性促动器支撑，这些支撑点分为3圈，由里到外，它们半径分别为R₁,R₂和R₃。这些支撑点位置分别在以R₁,R₂和R₃为半径外接正多边形的顶点上。

综上所述，WFST主镜轴支撑优化的目标函数，设计变量和状态变量，可设为：

均方根最小半光程误差为目标函数，它的设计变量为支撑点3个外接圆半径R₁,R₂,R₃；每圈起始支撑点与x轴角度θ₁,θ₂,θ₃；每个支撑点力F_J。为了简化支撑结构优化参数，所有支撑点的力和每圈起始角度设为定值：

这里N_P1、N_P2和N_P3分别第一圈，第二圈和第三圈支撑点数量。所以，目标函数只剩下R₁、R₂和R₃三个变量，则方程(14)可以简化成以下的形式

$\begin{array}{c}\hat{\Pi }={\mathrm{RMS}}_e=\sqrt{\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}{e}_i^2}\\ [{\gamma }_1,{\gamma }_2,L,{\gamma }_m]=[R_1,R_2,R_3]\\ g={\mathrm{RMS}}_g=\sqrt{\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}{\Delta }_i^2}\end{array}----\left(16\right)$

这里RMS_e和RMS_g都是通过APDL编程提取ANSYS计算结果并通过最佳抛物面拟合计算得到。未来有需要，本发明还可以对变量θ₁,θ₂,θ₃和F_J等进行优化。

有限元模型的建立

这里主镜材料采用的是微晶玻璃，如表1所示。

表1微晶玻璃材料性能

主镜参数化模型由APDL编程所建。

如图3-6所示，主镜参数化模型有两种扇形体组成：厚的扇形体和薄的扇形体。在网格生成的过程中，薄的扇形体区域生成密网格区域，厚的扇形体区域生成疏网格区域，内圈三个点分别约束住它们θ和z方向(柱坐标下)自由度，这样主镜是静定支撑，不会产生由约束引起计算误差；其他每一个点施加相同的力mg/N_P；整个主镜施加加速度g。

在优化的过程中，首先以较大的容差(该容差等于疏网格步长)采用零阶优化方法对主镜支撑进行优化，得到的是一个精度较低全局值。其次，对零阶优化结果中支撑点位置返回值进行判断，看这些支撑点位置是否落在密网格区域：如果是，则一阶优化法开始基于零阶优化法结果和之前网格设置在密网格区域以较小容差进行优化；如果不，则根据零阶优化结果对网格进行重新划分设置，保证之前零阶优化支撑点位置落在密网格区域，之后再启动一阶优化法基于零阶优化结果和重划网格设置在密网格区域以较小容差进行优化。这样，我们较快地得到的是全局的高精度最佳值。按这种方法，在中等性能PC(个人计算机)下，主镜支撑的每一种设计在一个小时内得到一个最佳方案。

有限元结果

WFST主镜支撑有以下三种Np方案。第一种是总支撑点Np＝27,内圈N_P1＝6点，第二圈N_P2＝9点，外圈N_P3＝12点；第二种是总支撑点Np＝39,内圈N_P1＝9点，第二圈N_P2＝12点，外圈N_P3＝18点；第三种是总支撑点Np＝54,内圈N_P1＝12点，第二圈N_P2＝18点，外圈N_P3＝24点。

图7、图10和图13分别表示N_P＝27,39和54时目标函数RMS_e随设计变量R₁,R₂和R₃变化而变化的四维图。

可见三种支撑点设计方案的目标函数相应最佳值分别为28.78nm,9.32nm和5.29nm。

目标函数RMS_e对设计变量R₁,R₂和R₃变化非常敏感，所以支撑点在密网格区域有助于计算收敛和提高计算精度。

图8、图11和图14分别表示三种支撑点设计方案最佳结果时主镜反射面相对原坐标变形云图，其形变云图是对称的。

图9、图12和图15分别表示这三种支撑点设计方案最佳结果时相对于新坐标系的最小半光程误差分布图。该分布图也是对称分布的，且每个支撑点所在的位置一般比别的区域变形要大，而且这些支撑点位置也是我们比较关注的区域和它们精度要求也比较高。所以，这也是我们要求支撑点所在区域网格要比其他区域要密的原因。在同样的计算机条件下，我们在支撑点所在位置对网格进行了细化，在别的区域对网格采取了相反做法，这样可以使我们得到较快计算速度和较高的计算精度。

如表2所示，为目标函数RMSe最佳值与此时相应设计变量R₁,R₂,R₃以及峰峰值PV。从表中可以看出，支撑点从27点增加到39点，RMS_e和PV有了急剧减少。从支撑点39增加到54点，RMSe和PV减小幅度很少。这说明，支撑点为54点时，RMS_e和PV接近了极限值。

依据WFST设计要求，在重力载荷下轴支撑设计均方根最小半光程误差一定要小于10nm。

从图7-15和表2可以看出,当N_P＝39和54它们最佳均方根最小半光程误差分别为9.32nm和5.29nm都满足WFST对主镜轴支撑设计要求。虽然N_P＝39时RMS_e＝0.029λ_min(λ_min＝320nm)大于N_P＝54时RMS_e＝0.017λ_min，但N_P＝39时支撑结构设计更为简单。所以在相同满足设计要求的前提下，N_P＝39时结构设计更接近于结构简化的极限值。

从以上的结果可以看出，变网格混合优化法是一种高效率伴随这高精度的优化方法。

表2最佳优化结果列表。