采用MVFSA算法的射电望远镜曲面优化研究

沈世云，王扬帆，罗吉润，姜诲名

(1.重庆邮电大学 理学院，重庆 400065；2.重庆邮电大学 通信与信息工程学院，重庆 400065)

0 引言

中国天眼FAST射电望远镜是国家“十一·五”立项的国家大科学工程。作为目前世界上单口径最大、灵敏度最高的射电望远镜，其落成启用对中国加快建设科技强国、实现高水平科技自立自强有重大意义。中国天眼FAST射电望远镜的控制系统主要包括望远镜总控制、馈源支撑控制和反射面系统三大系统。主反射面设计的目的是调节反射面的形状和方向，使观测天体发出的平行方向电磁波可以通过反射面集中在馈送舱的有效区域。将主动反射面调整为工作抛物面是主反射面技术的重要过程，该过程通过将主索结点上的促动装置与下拉索结合来实现。

目前，已有众多学者对大口径射电望远镜的主动反射面调节优化或智能算法优化进行了研究。朱丽春[1]通过主反射面控制系统控制策略，建立自适应、自学习机制节点位移控制模型，采用开闭环结合控制策略实现反射面高精度实时成型；薛建兴等[2]基于反射面单元动态面形精度分析研究FAST瞬时抛物面的拟合精度；王从思[3]根据促动器支撑节点位置确定该点促动器所在点的法线方向及促动器所在直线方程，计算出面向抛物面的大型天线主动面促动器最佳调整量；陈聪聪等[4]研究了一种确定网状反射面天线原理误差计算的方法；张建春等[5]基于传统模拟退火算法提出了更高求解效率的非常快速模拟退火算法；王星童等[6]利用快速模拟退火算法加快了路径规划问题的迭代收敛速度；陈华根等[7]提出了一种改进的非常快速模拟退火算法，相比传统算法有了更高的稳健性和求解效率。

由于传统方法全优化的时间冗长，本文基于现有研究，依据FAST射电望远镜的主动反射面相关结构和其主动反射面工作原理，针对反射面对不同方向天体观察时工作面需要调整的问题，建立以工作面与理想抛物面拟合程度最高为目标的优化模型，提出了采用改进的非常快速模拟退火算法来优化射电望远镜工作面，并对比不同算法的求解效率。

1 天眼FAST的组成与反射面状态

1.1 FAST组成

FAST 主要是由反射面、望远镜总控、馈源支撑三大控制系统、信号接收系统以及支承系统组成。其中，主动反射面系统结构如图1所示，是由背架、反射面板、主索网、促动器及下拉索等主要部件构成的一个可调节球面。每个主索节点连接一根与固定在地表的促动器连接的下拉索，通过下拉索和促动器的调节实现对主索网的形态控制。

图1 主动反射面系统结构示意图

1.2 主动反射面的基准态与工作态

主动反射面包括基准态与工作态。如图2的FAST剖面所示，基准球面与焦面是以点C为圆点的部分同心球面，其中两面半径差为F=0.466R(R=139.8 m)，基准面可视作半径R=300 m、口径D=500 m的球面。当FAST进入观测状态并需要接收天体S发出的平行电磁波信号时，主动反射面处于工作态。为了反射面可以将来自目标观测天体的平行电磁波汇聚到馈源舱上直径为1 m的有效圆盘区域，馈源舱接受平面的中心将在焦面上移动，具体位置为直线lSC与焦面的交点P。同时，主动反射面上有1个口径为300 m的区域需要调整为以点P为焦点的反射抛物面。

图2 FAST剖面示意图

2 理想抛物面的确定

2.1 瞬时过程的假设

通过文献[8]得知，FAST望远镜从基准态转变为工作态的过程是不断反馈、优化和调整的过程，通过馈源舱对平行电磁波的接收结果，控制系统对没有将光束射入馈源舱的反射面单元进行调整和优化。整个系统完全优化完毕共耗时约11 000 s，调整优化的过程十分庞杂。故在求取理想抛物面时，选取FAST即将从基准态转变为工作态时的瞬间进行分析，尽可能忽略变化过程的影响。

2.2 抛物面的准面定义

如图3所示，L1、L2、L3绕对称轴旋转360°得到同样的抛物面，其上任一抛物线的准线共平面称这个平面为抛物面的准面。对任意一条抛物线而言，抛物线上的任意一点到焦点的距离等于该点到准线的距离。可以推得，抛物面上的任意一点到焦点的距离等于该点到准面的距离。

图3 抛物面与抛物线的关系

2.3 理想抛物面方程式

2.3.1直线lSC方程和馈源舱焦点的确定

以图2中点C作为坐标空间系的原点。如图4所示，设点M为直线lSC上的任意一点，直线方向α=36.795°，β=78.169°。设m为直线SC的方向向量，可以求得直线lSC的方向向量为：

图4 直线SC位置示意图

m=(-cosβ·cosα，-cosβ·sinα，-sinβ)

lSC直线参数方程为：

求得直线lSC方程的参数形式可以减少变量数量，在仅知道t的情况下就可以得出直线上点的坐标，简化了求取结点坐标和理想抛物面的过程。

结合上文可知，由于点P在直线lSC上，所以考虑使用直线参数方程的形式来求解点P坐标。

根据图2可得CP=(1-0.466)R=160.2 m，所以tsp=160.2。将其代入式(1)可得点P的坐标为(-26.3，-19.7，-156.8)。

2.3.2理想抛物面方程的确定

图5是基准面和理想抛物面的剖面图。从剖面看，其抛物面的任意一条母线上的点到焦点的距离等于其到准线的距离。由2.2节可知，抛物面上任意一点到焦点的距离等于其到准面的距离，所以可根据此关系建立方程。

图5 理想抛物面与基准面示意图(剖面)

设抛物面上任意一点B坐标为(x，y，z)，已知点P的坐标为(-26.3，-19.7，-156.8)。根据空间中两点距离公式求出点B到点P距离dBP。

求取lSC与准面的交点。设直线lSC与准面交于点T。由抛物面性质可知，直线lSC垂直于准面，直线lSC的方向向量即为准面的法向量。联立直线lSC的方向向量m与点T坐标，得到准面方程为：

x+0.75y+5968z+2682.19=0

点B到准面的距离k以及点B到焦点P的距离dBP表示为:

(1)

(2)

令式(1)与式(2)相等，得到一个抛物隐式方程。

理想抛物面顶点坐标如表1所示。

表1 理想抛物面顶点坐标

2.4 工作面区域的确定

确定一个理想抛物面后，可以通过调节促动器的径向伸缩量将反射面调节为工作抛物面，使得该工作抛物面尽量拟合理想抛物面，以此使得平行电磁波尽可能多地反射汇聚到馈源舱上的有效圆盘区域。由于工作区口径只有300 m，所以需要先确定索网中促动器需要伸缩改变的主索结点。

图6 工作态口径范围剖面示意图

计算得到各个夹角θ。若θ∈[0°，ε+30°]，则主索节点处在口径范围内；若θ在此范围内，则主索节点不属于口径范围。为了保证工作面完整性，引入1个任意小的数ε。待调节的主索结点如图7所示。通过计算得到，共有695个主索结点需要。

图7 待调整的主索结点示意图

3 优化模型的建立

3.1 衡量拟合效果的面形精度

如图8所示，由于促动器只能沿径向运动，所以主索结点的运动轨迹是沿直线lAC运动的。P1为主索结点的最理想位置，位于理想抛物面上；P2为通过促动器调节后的主索结点的新位置；点P3为基准球面上的主索结点，即主索结点的初始状态；点P4为促动器；Δx1为主索节点与促动器相连的距离；Δx2为初始位置与新位置的径向差值；Δx3为新位置与最理想位置的径向差值。

图8 主索结点与理想抛物面的位置关系示意图(剖面)

当Δx3越小时，证明经过促动器调节后的主索结点新位置更靠近理想抛物面上的位置，此时调节效果最好。为此，引入面形精度RMS来衡量反射面形精度。反射面面形精度是衡量射电望远镜综合性能的重要指标之一，直接影响望远镜的工作效率。本文中用反射面主索结点实际坐标与理想抛物面拟合误差来衡量反射面的RMS：

其中：Xi、Yi、Zi分别表示第i个节点调节成理想抛物面时的理想坐标，位于理想抛物面上；xi、yi、zi表示为第i个节点调整后的实际坐标。计算695个实际坐标与理想坐标的距离平方和的平均值，然后再开平方得到RMS。

3.2 优化模型的约束条件

3.2.1约束条件1：伸缩长度限制

各结点促动器伸缩调整长度为ti，其中i=1，2，3，…，695。根据实际意义，认为趋向球心时ti取正，远离球心时ti取负。长度调整需要在一定范围内，因此可以得到第1个关于促动器伸缩调整长度的约束条件：

-0.6≤ti≤0.6

3.2.2约束条件2：伸缩方向限制

促动器仅可以在径向方向上变化，所以采用直线方程的参数形式进行径向变化量表示：

3.3 优化模型的目标函数

设调整后节点之后，第i个结点原坐标调整为ui、vi、wi。在调整坐标时，通过调整各促动器的伸缩量来实现各结点向理想坐标位置逼近，从而实现工作面向理想抛物面逼近。综合考虑RMS精度值越小越好的目标，以及需要满足调整长度的范围和调整方向的约束条件，建立如下优化模型：

695个主索结点对应关系见表2。其中，单位方向向量为原坐标指向球心的向量单位化后所得向量。

表2 695个主索结点对应关系

通过求解上述优化模型可以得到最为拟合理想抛物面的工作面和各个主索结点的变化量ti。

4 优化模型的求解

4.1 模拟退火算法

1953年，Metropolis提出模拟退火算法(simulated annealing，SA)；Kirkpatrick在1983年成功使用SA解决了组合优化问题。SA作为一种有效的非线性组合寻优算法，在理论上已被证明是严格的，且有大量的实际应用证明该算法有效，可靠，可处理规模较大的背包问题，且能够鲁棒地得到满意解[9-10]。SA算法的关键思想是：算法会有一定概率接受一个劣解并跳出某个局部最优解空间并尝试寻找更优解。从全局搜索的角度来看，算法会尽力朝最优解进行搜索寻优。

SA算法原理和金属固体退火的原理相似，就是将热力学的理论运用到统计学中。从机理角度解释算法思想就是先将金属固体加温到一个足够高的温度，再逐渐降温冷却；在降温的过程中粒子会自发地进行热运动以使得固体状态朝着自由能减小，以此朝着能达到自由能最小的平衡态进行单点寻优。具体步骤如下：

步骤1给定一个较高且合适的初始温度T0，在符合模型给定变量区间内随机选择一个初始模型m0，并计算初始模型对应的目标函数值E(m0)；

步骤2采用合适的扰动方法，在符合模型给定变量变化区间内，扰动当前模型产生一个新模型m，并计算新模型的目标函数值E(m)，同时计算目标函数差值ΔE=E(m)-E(m0)；

步骤3在当前温度T下，若ΔE<0，则接受新模型m；否则，按照概率P接受当前较劣模型；当模型被接受时，用新模型m替换初始模型m0；

步骤4在温度T下，重复进行多次扰动和接受判别，即重复步骤2与步骤3共Lk次；

步骤5通过一定的降温退火计划下降温度T；

步骤6重复步骤2—5，直到温度下降到终止温度或者搜索达到收敛条件，认为本次已经搜索到一个较优模型解。

4.2 改进的非常快速模拟退火算法

Ingher在1989年提出了一种模拟退火算法的改进算法，即非常快速模拟退火算法(very fast simulated annealing，VFSA)。VFSA算法较SA算法具有更高的效率，使得模拟退火理论真正可实际应用并能解决更多的实际问题。VFSA算法的流程与SA算法流程一致，只在模型扰动、接受概率或者退火计划上作出了改进。相比传统算法，该算法在高温时有着更高的全局搜索效率。但由于在低温下仍然进行全局搜索，搜索的扰动幅度过大，难以继续搜寻更优模型，因此在低温状态下搜索方式可以作出改进以增加低温搜索效率。本文基于改进的非常快速模拟退火算法(modified very fast simulated annealing，MVFSA)[7]，利用Matlab对模型进行求解。

4.2.1模型扰动

模拟退火算法是在当前温度下扰动当前模型产生一个新模型，算法扰动由一个随机分布函数控制。VFSA算法使用依赖温度的似Cauchy分布扰动产生新模型。模型扰动规则如下：

ri=T·sgn(γ-0.5)·((1+T-1)|2γ-1|-1)

(3)

mi′=(0.5+ri)·(mimax-mimin)

(4)

其中：mi′表示扰动后的模型，mimax和mimin分别表示mi取值范围的上、下限；ri是扰动因子；γ表示(0，1)区间上的随机数；sgn(x)为符号函数。

4.2.2接受概率

模拟退火算法中新模型根据目标函数的差值ΔE，按照Metropolis准则来接受新模型，目的是跳出局部最优解。该准则是传统模拟退火的重要环节。

1)若ΔE<0(新模型较优)，则新模型无条件被接受；

2)若ΔE>0(新模型较劣)，计算接受概率P：

P=exp(-ΔE/T)

(5)

在ΔE>0的情况下，产生一个分布在[0，1]内的随机数ξ，若P>ξ则接受新模型，否则舍弃新模型。

4.2.3退火计划

诸多实践证明，采取温度T以指数变化方式降温的计划较符合固体退火的实质：

T=T0·exp(-ν(j-1)0.5)

(6)

其中：ν为模拟退火算法中的温度衰减系数；j为算法迭代次数。

4.2.4局部搜索阶段

在4.2.1与4.2.3节中的扰动与退火目的是进行全局搜索以找到最优区间解。当完成全局搜索算法达到收敛时，进行回火升温操作，此时从一较低温度开始退火，继续在已找到的最优解区间内进行局部搜索寻优。本阶段中改进的模型扰动与退火计划表示为：

mi′=mi+(ri-0.5)·(mimax-mimin)/L(j)

(7)

T=T0·exp(-ν(j-j1/ρ)0.5)

(8)

式(7)和(8)由式(5)和式(6)改进而来，其中L(j)为算法搜索空间限制因子，与算法迭代次数有关。当迭代次数j增大时，L(j)也增大，这意味着搜索范围将变小，通过搜索空间限制因子来达到算法在一个更小的搜索区间内寻找更优模型的目的。j1为上一阶段迭代次数；ρ为降温形状因子，与当前温度T成反比关系，意味着当T越小时ρ越大。算法的改进目的是使得搜索初期即在高温时有较大的搜索范围和劣解接受概率，在低温时可以在已找到的最优解区间继续以较小的搜索范围寻找更优解。

4.3 求解待调整主索结点位置和促动器伸缩量

用Matlab语言编写MVFSA算法对3.3节中的模型进行全局搜索和局部搜索。设置MVFSA算法的温度衰减系数ν=0.90，初始温度T0=200，终止温度为Tend=10-7。

通过随机函数的生成，选取1个在模型取值区间内的合法初始模型m0(t1，t2，…，t695)，并根据面形精度计算该初始模型的函数值RMS0。每次达到平衡后都会进行降温，即采用4.2.3的退火计划，直到T

在相同初始解下，进行50次的MVFSA算法搜索，取平均值得到RMS=0.167 533。695个主索结点对应关系见表3。

表3 695个主索结点对应关系

5 结果分析与算法比较

5.1 结果分析

如图9所示，在MVFSA算法下，算法迭代13次找到1个全局最优解区间，在温度很低时再经过4次迭代在该区间内找到更优解，提高了结果的精度和准确性。但是解的优化效率并没有提高很多，原因是模型变量多达695个，且每个变量取值区间相对搜索范围较狭窄，很难大幅度提高寻优结果。

图9 MVFSA算法各代最佳函数值曲线

5.2 算法比较

为了验证算法的高效性，将MVFSA算法、VFSA算法和传统的SA算法进行比较。结果如表4与图10所示。

表4 不同算法寻优结果

图10 不同算法各代最佳函数值曲线

比较在同一初解、相同初始温度T0=200、终止温度Tend=10-7和温度衰减系数ν=0.90时3种算法的搜索效率。在迭代开始时，VFSA算法和MVFSA算法迭代曲线斜率的绝对值远大于传统SA算法，且在迭代约8次时就快速找到了1个较优解区间。而在相同迭代次数下，传统SA算法仍在进行低效或者无效的搜索。

MVFSA算法在适当回火至一个较低温度时进行局部搜索，相比VFSA算法多进行了5次迭代搜索，收敛也十分迅速。二者总搜索时间几乎一致，但采用MVFSA算法得到的RMS值较原来下降了2.78%，意味着拟合理想抛物面程度更高。

上述结果表明：求解本文中建立的优化模型，采用MVFSA算法有更高的搜索效率和更优的搜索结果。

6 结论

针对FAST射电望远镜的主动反射面调节优化问题进行研究。通过机理建模的方式从下拉索的调节量出发，基于改进的非常快速模拟退火的智能算法对衡量拟合效果的面形精度优化模型进行求解。以695个待调整主索结点的下拉索促动器伸缩量为研究变量。该算法在非常快速模拟退火算法全局快速搜索的基础上，继续在局部进行快速搜索，在几乎相同且短暂的总搜索时间下得到更优的搜索结果，使得主动反射面的工作面更加拟合理想抛物面。

实际中，采用该算法或智能算法来解决射电望远镜的主动反射面调节案例较少。本文中提出的优化模型可以推广到其余口径天线或者射电望远镜的主动反射面调节优化中，提供了一种较好的优化思路，有利于馈源舱更好地接受天体发出的电磁波信号。后续研究将考虑更多影响因素，进一步提高算法的效率和准确性，提高主动反射面的工作效率。