基于径向最大面截距优化下的FAST反射面形状调节

黄建平，朱德凯，李星星，杨秀丽

（台州学院 电子与信息工程学院，浙江 临海 317000）

0 引言

500 m口径球面射电望远镜（Five-hundred-meter Aperture Spherical radio Telescope，FAST）的工作原理是：FAST在反射面板调节约束下，确定一个理想抛物面，将反射面通过促动器径向伸缩量的调节变为工作抛物面,获得天体电磁波经反射面反射后的最佳接收效果，从而将来自目标天体的平行电磁波反射汇聚到馈源舱的有效区域。为了获得尽可能理想的抛物面，提高馈源舱的信号接收效果，本文主要考虑解决以下3个问题：

（1）当待观测天体S位于基准球面正上方，即α=0°,β=90°时，结合考虑反射面板调节因素，确定理想抛物面。

（2）当待观测天体S位于α=36.795°,β=78.169°时，确定理想抛物面。建立反射面板调节模型，调节相关促动器的伸缩量，使反射面尽量贴近该理想抛物面。

（3）基于问题2的反射面调节方案，计算调节后馈源舱的接收比，即馈源舱有效区域接收到的反射信号与300 m口径内反射面的反射信号之比，并与基准反射球面的接收比作比较。

1 模型的建立和求解

1.1 建模前的准备

1.1.1 理想凹面镜的方程建立

从几何光学得到理想凹面镜的方程如下：

其定义域D为：x2+y2≤R2。进而确定该凹透镜的焦点坐标F为（0，0，p/2）。

1.1.2 相邻节点之间的变化

由于制程界限的限制，主索节点在径向上可以发生±0.6 m的变化，而这种变化可能会使相邻节点之间的距离发生一定的变化，从而导致最大形变量不会超过0.07%（这是一个硬性的技术指标）。经验证，±0.6 m的变化是安全的。具体方法是对一块反射面板所在的三角网格进行形变分析，其形变情况如图1所示。图1左边为三角网格的立体形变情况图［1］，右边为三角网格一条边的最大形变情况图。

图1 三角网格形变情况图

为了后续能够充分利用主索节点±0.6 m的径向变化，我们需要保证最大形变量能控制在0.07%内。最大形变量的计算公式为：

式（2）中η即为形变量，其中l、l′已在图1中给出几何意义，Δr为径变化量。对于l的计算方法，采用每两个主索节点之间距离的平均数，具体计算公式如下：

式（3）中的P代表节点坐标，N取2226。通过Matlab计算，可以得到最大形变量η=0.053%。

从结果可以发现，主索节点在径向上变化±0.6 m，其形变量并不会超过0.07%，因此后述中不需要再考虑形变量过大所导致的一系列问题。

1.1.3 空间坐标的变换矩阵

由于基准球面具有各向同性，并且为了在后文中方便计算，需要利用到空间坐标的变换矩阵，即

1.2 问题1的建模与求解

1.2.1 FAST的理想抛物面方程

通过下拉索和促动器的配合工作，能够保证主索节点充分地靠近理想抛物面，在此先不考虑径向伸缩范围的限制，重点研究FAST理想抛物面方程的构建。由于方位角为0°，仰角为90°，这说明星体是处于FAST的中心正上方，于是可以由式（1）得到此时FAST的理想抛物面方程应为：

其定义域D为：x2+y2≤R12。这里的R1取150 m。由于需要保证射向照明区域的平行光经过反射面板反射能够聚焦在焦面上，则必须要满足方程的焦点落在该焦面上，其中焦点坐标应为(0,0,-0.534R)，焦距为p/2，于是可以得到p和h之间的关系式为：

由式（6）可知p可以直接由h表示，所以只需要通过确定最优的h，便可确定最理想的抛物面方程。

1.2.2 径向最大面截距的定义

设散点集合为Q={Q1,Q2,…,Qn}，中心点为C，其中对∀i∈{1,2,…,n}，满足Qi∈ ℝ3，并且C∈ℝ3。现有一函数f也定义在ℝ3并且满足f:ℝ3↦ℝ，我们称散点集合Q到函数f的径向（中心点方向）最大面截距为Dm，其几何意义如图2所示。

图2 最大面截距几何示意图

首先，通过 ∀i∈{1,2,…,n}，可以得到Qi与C的距离，设Qi的坐标为(xi,yi,zi)，C的坐标为(xc,yc,zc)，于是可以得到直线方程：

将式（7）与f联立求解，并选取方向上的解作为节点，记为，最后便得到径向最大面截距的计算公式为：

1.2.3 问题1模型求解

由式（6）表明p受h影响，因此确定最优的h是问题的关键所在。为此，这里采用有限元分析的方法［2］，将以最小化径向最大面截距为目标函数进行求解，得

下面从两个方面确定最优解所在区间。

（1）对解域进行划分。

由现实情况可以粗略地确定-B≤h≤B,其中B取0.6 m，也就是说h的最优解域不可能超出这个范围，下面将解域进行离散化。

设解域为A，在空间维度上对其进行均匀划分，即将A划分为N份，每段的长度为Δh=2B/N，于是得到离散点集H={-B+kΔh|k=1,2,…,N}，与此同时，解域A被分割成子解域。

（2）确定最优解所在的子解域。

将所得到的Hi(i∈ 1,2,…,N)代入式（5）中，得到其定义域D为：x2+y2≤R21。根据所有主索节点的坐标，将这些散点构成的集合记为Q，于是利用式（8）便可以得到与Hi相对应的径向最大面截距Dm(Hi)，据此找到min{Dm(Hi)|i=1,2,…,N} 所对应的下标，将其记为b，最后确定最优解所在子解域为[Hb-1,Hb+1]。

通过上述（1）（2）两点的分析可知，在确定了最优解所在区间之后，只需将N变得足够大就可以确定其近似最优解，最终得到的近似最优解为：hbest=-300.2860m。经过计算发现，其所对应的径向最大面截距Dm(hbest)为0.3893 m，也就是说每个促动器的伸缩量均在径向伸缩范围内。最后，再将hbest代入f的方程，便可以得到理想抛物面方程为：

其定义域D为：。该理想抛物面图像如图 3所示。

图3 理想抛物面图像（深色部分）

1.3 问题2的建模与求解

1.3.1 反射面板和理想抛物面的贴合方式分析

由于本文默认主索节点均分布在旋转抛物面上，下面对三种贴合方式进行分析，具体情形如图4所示［3］。

图4 各种反射面板的贴合方式

根据图4得出板间离差公式为：

经过计算发现，板间离差的最大值为54.9 μm，说明图4中第3种采用的贴合方式是较为可靠的，可以在后文中使用。

1.3.2 平行光线的垂直入射板

平行光线从方位角α=36.795°、仰角β=78.169°的方向射向基准球面，由此可知经过圆心C的光线必定会和基准球面垂直，联立球面方程和光线方程求得目标解的坐标为G(-49.2544,-36.8403,-293.6270)。由于每个主索节点只能在径向上进行变化，故垂直入射点并不会因为反射面板的移动而发生改变。

根据垂直入射点不变的性质，再针对某一块板或某一个主索节点计算散点集Q中每个点到G的距离，至少筛选出其中最近的3个点，分别记下它们的坐标和编号，具体如表1所示。

表1 垂直入射点信息

通过相关的数据分析发现，所有主索节点到达G的距离均比R大了0.4 m，再结合表1可得出星光的垂直入射点就是编号为“D27”的主索节点，也就是理想抛物面的顶点，其坐标见表1第2列。

1.3.3 理想抛物面以及各工程参数的求法

由于平行光线是从方位角α=36.795°、仰角β=78.169°的方向射向基准球面的，这会给计算带来很大的不便。根据FAST光学面以及主索节点所在半球面所具有的球面对称性（各向同性），采用式（7）将原散点集Q变换成为新的散点集合，这样平行光线能够从正上方射入，变换矩阵为：

通过新的散点集即可求出最佳的hbest=-300.3360m，其所对应的Dm(hbest)为0.3359 m。再将其代入新坐标系下的抛物面方程，便可计算每个促动器的伸缩量大小，以及新坐标系下理想抛物面上各节点的坐标′及其编号，最后通过逆变换求出原坐标系下理想抛物面上各节点的坐标Q′以及原坐标系下抛物面的方程。变换坐标系下的抛物面方程为：

其定义域D为：。原坐标系下理想抛物面的方程为：

理想抛物面在原坐标系下的图像如图5所示。

将藜麦种子磨成粉末，过100目筛，取筛下物放入烘箱烘干备用。称取适量预处理后的藜麦粉并将其与0.20%的NaOH溶液按照1：5的比例混合，搅拌3 h后，静置18 h。将混合物于4000 r/min离心10 min，去上层黄色蛋白质沉淀和下层灰色沉淀，留中间白色淀粉反复离心、洗涤，调pH至中性，继续洗涤至上层无黄色沉淀为止。最后将产物置于 40 ℃烘箱中干燥48 h，粉碎，过100目筛后封袋保存。

图5 理想抛物面方程图像

经过分析发现，调节后反射面300 m口径内的主索节点一共有692个，部分主索节点编号及其位置坐标和促动器伸缩量如表2所示。

表2 300 m口径内的主索节点编号、坐标和促动器伸缩量

1.4 问题3的建模与求解

1.4.1 基准球面的接收比

对于基准球面，其光学面是一个球面，由于球面具有空间对称性和各向同性的特点，故只需考虑α=0°，β=90°的情况即可。光学面的光路反射图如图6所示。

图6 基准球面的光路反射图

对图6分析后可知，光线进入半球面时和地面的夹角为0°，经过反射后其夹角变为θout=2θ-π/2，再加上空间角的限制，可得出以下结果：

式（16）中的范围根据对称性只考虑了半边。对于出射光线能否射入馈源舱，我们只需要考虑出射光线lout和直线y=F-R的交点坐标x(θ)即可。如果它能够在馈源舱的有效半径内，就说明这部分光线经过反射能够被接收到。根据入射点Pin的坐标(-Rcosθ,-Rsinθ)，再结合出射角，便可得到反射光线所在的直线方程为：

将式（17）与y=F-R联立即可得到：

函数图像如图7所示。

图7 函数x的图像

对式（18）求驻点，发现不存在解析解，但是能够给出其数值解为1.3522，其对应的x=-3.0581，这说明反射光线一开始不能射入馈源舱，接着一个小角度内能射入馈源舱，然后又不能射入馈源舱，最后再次能够射入馈源舱，共4种状态。这4种状态的临界θ值如表3所示。

表3 4种状态的临界θ值

由于光线是垂直入射的，即在面板上单位面积接收到的能量大小是不一样的，所以这里需要进行相应的调整，平行光能流通量如图8所示。

图8 平行光能流通量图

因为光具有“沿直线传播”的特性，而它同时也是一种电磁波信号，所以可以将其看做“无旋场”；而对于来自遥远星体的平行光线，其能流密度Φ在空间上是均匀的，于是将图8的“半球壳+顶部”看做一个封闭的曲面，得到方程：

这就说明整个半球面的光能流通量等于半球壳顶部的光能流通量，根据式（19）可知：

通过前面所得到的θi便可以求出其所对应的Xi，其计算方法为：

将所得到的Xi记入，如表4所示。

表4 4种状态的临界X值

通过表4中X的临界值，可以确定能够打到馈源舱上的入射点X所在的区间范围为：X∈[1 02.4737,107.8680]∪[0 ,7.3782]，可求得基准反射球面的接收比为：

从式（22）中可以看出，接收比很小。由于未考虑板间缝隙以及馈源舱遮挡情况所带来的误差，实际上基准反射球面只能够将少部分的光汇聚到馈源舱内。

1.4.2 调节后馈源舱的接收比

由于问题2中说到光线的入射方向为α=36.795°，β=78.169°，说明光线并不是垂直顶面入射的。为了避免光线斜射带来的问题，首先可以采用问题2中的变换坐标系来解决；接着利用Monte Carlo算法，通过随机产生入射光线得到能够进入馈源仓光线的数量；最后计算接受比。但是因为所有的反射面板是离散分布的，所以确定入射面板是解决本问题的关键，从以下三方面展开讨论：

（1）确定入射面板。

由于结果是在变换后的坐标系中进行分析的，这里采用的也是变换坐标系下的坐标以及变换坐标系下的抛物面方程，将变换坐标系下的主索节点向z=0平面进行投影得到的点集记为。通过相关数据，可以确定每个反射面板的三个主索节点对应的编号，利用该编号便可以找到对应的主索节点坐标。将该坐标按照T矩阵进行坐标变换便可得到变换坐标系下的主索节点坐标，也就是将编号与进行了绑定。

采用Monte Carlo算法随机产生Xi，其坐标满足，并且其在z轴方向上的投影为0。计算Xi与点集中每个元素的距离，取3个距离最近的点，并找到每个点所对应的，再找到每个点对应的编号，由此可以确定反射面板。下面给出1个入射点与最近的3个点关系的具体例子，结果如表5所示。

表5 入射点与最近的3个点关系的一个具体实例

（2）反射光线的求解。

利用各节点的编号及其所绑定的，便可计算每个面板的法向量［4］，并选取在z轴方向投影大于0的法向量，得到

反射原理图如图9所示。

图9 反射原理图

结合施密特正交分解的方法容易得到反射光线的方程为：

（3）Monte Carlo算法的实现。

将之前得到的反射光线与z=F-R联立，得到交点坐标为(xhi,yhi)。通过判断该点是否满足，便可确定该反射光线能不能被馈源器接收。也就是说，只要给定入射光线在顶面生成的坐标，就可以得到该光线是否能通过馈源仓。此处会发现“枚举法”失效，原因在于顶面半径为150 m，如果想让误差小于0.1 m，则至少需要进行200万次循环迭代，从而使时间拉得很长。因此，我们考虑了Monte Carlo算法，并将初值N设为100万，利用MATLAB运行，得到馈源舱有效区域接收到的反射信号与300 m口径内反射面的反射信号之比为=29.50%。这就说明基准球面经过调整成为类抛物面的时候，来自星体的平行光线能够尽可能多地被馈源舱接收到。

2 结语

针对FAST工作原理，本文对FAST主动反射面板的形状调节进行了研究与探讨。为了提高其工作效率，必须将反射面尽可能贴近理想抛物面。因此文中提出一种基于径向最大面截距的优化模型，通过采用有限元分析、蒙特卡洛等算法模拟得到一系列工程参数在该模型下的标准取值；另外将该模型下建立的工作态反射面与基准态反射面对馈源舱信号的接收比进行了比较。结果表明，本文所提出的方法具有一定的优越性。