对于射电望远镜主动反射面调节问题的一点思路

题目

中国天眼——500米口径球面射电望远镜（Five-hundred-meter Aperture Spherical radio Telescope，简称FAST），是我国具有自主知识产权的目前世界上单口径最大、灵敏度最高的射电望远镜。它的落成启用，对我国在科学前沿实现重大原创突破、加快创新驱动发展具有重要意义。
FAST由主动反射面、信号接收系统（馈源舱）以及相关的控制、测量和支承系统组成（如图1所示），其中主动反射面系统是由主索网、反射面板、下拉索、促动器及支承结构等主要部件构成的一个可调节球面。主索网由柔性主索按照短程线三角网格方式构成，用于支承反射面板（含背架结构），每个三角网格上安装一块反射面板，整个索网固定在周边支承结构上。每个主索节点连接一根下拉索，下拉索下端与固定在地表的促动器连接，实现对主索网的形态控制。反射面板间有一定缝隙，能够确保反射面板在变位时不会被挤压、拉扯而变形。索网整体结构、反射面板及其连接示意图见图2和图3。

图1 FAST三维示意图

图2 整体索网结构

左边为(a) 右边为(b)
图3 反射面板、主索网结构及其连接示意图
主动反射面可分为两个状态：基准态和工作态。基准态时反射面为半径约 300米、口径为500米的球面（基准球面）；工作态时反射面的形状被调节为一个300米口径的近似旋转抛物面（工作抛物面）。图4是FAST在观测时的剖面示意图，C点是基准球面的球心，馈源舱接收平面的中心只能在与基准球面同心的一个球面（焦面）上移动，两同心球面的半径差为F=0.466R（其中R为基准球面半径，称F/R为焦径比）。馈源舱接收信号的有效区域为直径1米的中心圆盘。当FAST观测某个方向的天体目标S时，馈源舱接收平面的中心被移动到直线SC与焦面的交点P处，调节基准球面上的部分反射面板形成以直线SC为对称轴、以P为焦点的近似旋转抛物面，从而将来自目标天体的平行电磁波反射汇聚到馈源舱的有效区域。

图4 FAST剖面示意图
将反射面调节为工作抛物面是主动反射面技术的关键，该过程通过下拉索与促动器配合来完成。下拉索长度固定。促动器沿基准球面径向安装，其底端固定在地面，顶端可沿基准球面径向伸缩来完成下拉索的调节，从而调节反射面板的位置，最终形成工作抛物面。
本赛题要解决的问题是：在反射面板调节约束下，确定一个理想抛物面，然后通过调节促动器的径向伸缩量，将反射面调节为工作抛物面，使得该工作抛物面尽量贴近理想抛物面，以获得天体电磁波经反射面反射后的最佳接收效果。
请你们团队根据附录中的要求及相关参数建立模型解决以下问题：
1、当待观测天体𝑆位于基准球面正上方，即𝛼 = 0°, 𝛽 = 90°时，结合考虑反射面板调节因素，确定理想抛物面。
2、当待观测天体𝑆位于𝛼 = 36.795°, 𝛽 = 78.169°时，确定理想抛物面。建立反射面板调节模型，调节相关促动器的伸缩量，使反射面尽量贴近该理想抛物面。将理想抛物面的顶点坐标，以及调节后反射面300米口径内的主索节点编号、位置坐标、各促动器的伸缩量等结果按照规定的格式（见附件4）保存在“result.xlsx”文件中。
3、基于第2问的反射面调节方案，计算调节后馈源舱的接收比，即馈源舱有效区域接收到的反射信号与300米口径内反射面的反射信号之比，并与基准反射球面的接收比作比较。

摘要

本文针对“FAST”主动反射面形状调节的应用问题，通过旋转矩阵法，最小二乘法，建立理想抛物面模型，目标优化模型，旨在得出不同情况下调节主动反射面的最佳情况。
针对问题一，本文建立单目标优化模型，以接收区域面积s为目标函数，以促动器变化范围和原抛物线方程为限制条件，得出抛物面方程为：

在分析过程中，采用降维的方法将空间问题转化为平面问题，通过将立体图形投影至xoz平面，选取外围点作为拟合数据。采用最小二乘法拟合外围数据，得出主索节点所在平面方程，拟合优度 。以反射信号刚好能达到接收圆盘为极限情况，得出最大入射角应小于。依据最大入射角确立有效区间，对该区间内所有点的偏移不超过0.6作为限制条件，得出理想抛物面方程。
针对问题二，本文采用旋转矩阵法，以均方差误差为目标函数，以主索节点浮动范围为限制条件，建立目标优化函数模型为：

根据观测天体的方位角和仰角，通过旋转矩阵法，重新建立空间直角坐标系，并求出对应主节锁点坐标，在新空间直角坐标系中对有效接收范围曲面进行优化，得到理想抛物面方程。通过计算得出顶点坐标为（-49.34，-36.8645，-230.002）。具体对应编号的促动器伸缩量见附件四。
针对问题三，本文建立反射模型。通过对一个反射面板进行分析，建立以入射点为原点的空间直角坐标系，对入射信号和反射信号求解法向量，确立反射信号在馈源舱接收平面的点坐标。通过确立平面内三个点得出平面方程，并依据中心圆盘直径确立范围。

根据范围内的反射点坐标，反推工作抛物面内的有效主索节点，得出馈源舱的接收比为0.79，基准球面的接收比为0.025，优化结果较好。

问题分析

对于问题一的分析

问题一要求结合反射面调节等因素，确立理想抛物面。本文对反射面板进行分析，不难得出所有反射面板均处于主索节点构成的抛物面上。将主索节点坐标可视化，构建三维立体图，确立基准态时的近似抛物面方程。结合入射波和反射波，计算馈源舱可以接收到电磁波的有效接受区域，确立反射范围。通过调节促动器的顶端径向伸缩量，在伸缩范围内尽可能使得更多区域处于有效接受范围内，对抛物面采用降维分析抛物线最适开口大小，从而得到理想抛物面。

对于问题二的分析

问题二要求依据天体S，确立观测天体的理想抛物面。本文依据口径300米的抛物面为前提，以到截面距离为约束条件，确立需要调整的主索节点。对于偏转后的抛物面，本文采用旋转矩阵，以天体和馈源舱所在直线为z轴，重新建立空间直角坐标系。并依据旋转矩阵，得出对应主索节点的新坐标。以理想抛物面各点的均方差误差最小为目标函数，确立每个点的偏移距离，并确立顶点坐标。最后依据旋转矩阵的逆变换得出原坐标系中的顶点坐标。

对于问题三的分析.

问题二要求依据天体S，确立观测天体的理想抛物面。本文依据口径300米的抛物面为前提，以到截面距离为约束条件，确立需要调整的主索节点。对于偏转后的抛物面，本文采用旋转矩阵，以天体和馈源舱所在直线为z轴，重新建立空间直角坐标系。并依据旋转矩阵，得出对应主索节点的新坐标。以理想抛物面各点的均方差误差最小为目标函数，确立每个点的偏移距离，并确立顶点坐标。最后依据旋转矩阵的逆变换得出原坐标系中的顶点坐标。

模型假设

结合本题的实际，为确保模型求解的准确性和合理性，本文排除一些因素的干扰，提出以下几点假设：
1．假设所有促动器没有损坏；
2．假设入射波均平行射入主动反射面；
3．假设经过馈源舱的入射信号不被阻挡。

模型的建立与求解

经过以上的分析和准备，本文将逐步建立以下数学模型，进一步阐述模型的实际建立过程。

问题一的建立与求解

要确立理想抛物面，即可认为确立一抛物面，使得馈源舱接受的信号达到最大值。由题可得，馈源舱接受信号的有效区域为直径1米的中心圆盘，故首先需确立主动反射面的有效接受波范围。依据分析过程，确立以下过程：主索节点可视化——确立抛物面方程——计算有效接收区域——调节促动器伸缩量——确定理想抛物面。

确立基准态时的理想抛物面

主动反射面是由反射面板和主索网构成，而反射面板是铺设在主索网构成的三角网格上，故可以通过主索节点的坐标确立主动反射面板的大致曲面。通过Matlab对附件中的坐标进行可视化，得到主动反射面的三维立体图如下所示：

通过图像分析，得出主索节点由内而外呈五边形向外围螺旋扩散，且主索节点在空间分布上按照编号体现区域化，整体上可以大致视为抛物面。而抛物面可看作一条抛物线绕某一条轴旋转得到，故本文采用降维的方法，通过选取不同截面求出所对应的抛物线方程，进而得出对应抛物面表达式。
设抛物线的公式为，其中 为待定参数。要确立一条过坐标原点的抛物线，仅需选取两个图形上的点即可确立方程表达式。将抛物面视为由抛物线绕z轴旋转得到，用 代替x即可得到抛物面方程。
对于计算多组抛物面方程，本文采用最小二乘法进行函数拟合。根据上述拟合模型表达式，计算偏差平方和，对其进行求导，化简得到范德蒙德矩阵：

要得到截面数据，以投影在XOZ平面为例，以同一z值下的最大值和最小值作为拟合数据进行拟合，得到结果如图所示：

拟合方程为：

对于x的一次项，其量级远小于其他的参数量级，故本文将x的一次项式舍去，调整结果为：

为评价拟合结果，本文采用拟合优度 来评价拟合的结果。当
越接近1时，即说明误差平方和与0越接近，表明拟合的效果越好，计算公式如下：

通过Matlab对上式进行运算，得出 ，接近1，故认为拟合的效果很好，可以作为该抛物面的函数表达式。

有效信号接收面积

对于射电望远镜，要接收到入射波的信号，必须要确保入射波经反射后，可以打到馈源舱直径1米的中心圆盘上，这样便存在一个理论有效接受区域，如图所示。

以基准球面为坐标原点O，建立平面直角坐标系，假设入射波与法线的夹角为 ，入射波与主动反射面的交点为A，构建入射波反射图：

由焦径比为0.466得出p点坐标为（0，-160.2）。由题目可得，基准态时的半径为300m，且由基准态变为工作态时不对半径造成影响，可得OA=300m。根据口径为300m，得出AB=150m，OAB构成∠BAO=60°的直角三角形，进而得出B点的坐标为（0，-277.13），BP=116.93m。在计算入射角时，考虑到馈源舱的接收区域为直径1m的中心圆盘，在计算时还需减去圆盘宽度。设入射角为 ，根据正切定理可建立方程：

通过Matlab进行求解，可得 =25.6147°，所以入射角应小于25.6147°。

计算理想抛物面

在保证入射角不变且处于最大值的情况下，要使得接收信号波最多，意味着需要有更多的主动反射面处于可接收范围内。对于馈源舱来讲，接收区域始终为以P点上方一焦点，向两边各以38.7706°为夹角构成的扇形接收区。通过构建相似扇形，即延长OA距离，可以使得接收区域变多。
改变OA长度与促动器有关，要使得OA改变，即可看做促动器变化 的长度。以接收区域面积s为目标函数，以促动器变化范围和原抛物线方程为限制条件，建立目标规划方程：

由Matlab计算可得当时，s有最大值，此时的抛物线方程为：

即可确立理想抛物面的方程为：

问题二的建立与求解

确立需改变的主索节点

由题目可得共有2226个主索节点，而控制反射面是由主索节点对下拉索推进或者拉回实现的，现只需筛选出需要变化的主索节点，即可通过局部变化来得到最佳曲面，从而减少计算量。

假设改变后垂直于中心线SCP且截面圆直径为300m的平面为，设的表达式为：

其中， 。由题目可得基准球半径R=300，以基准球心作为坐标原点，根据勾股定理求出基准球心与平面 距离D的表达式为：

经计算得出，D约为1267.2131。对于每一个处于照明区域的主索节点 均满足以下目标规划模型：

通过Matlab计算，筛选出699个主索节点可以作为调整点

上图为选择作为的调整主索节点的对比图，其中红色点为筛选数据，蓝色点为原始数据，由图可知，筛选后的主索节点都位于西南方。

旋转前后点坐标关系

当天体S相对于球心的方位角为时，可以将所有节点变化到以z轴为顶点的抛物面上，即对空间直角坐标系进行旋转，以旋转后的SCP直线作为z轴，以C点作为原点建立空间直角坐标系。（S为天体，C为球心，P为馈源舱）
经分析，此时中心线的一个法向量为，建立SCP直线方程：

其中，d为常数项，和分别为方位角和仰角。结合问题一中计算所得馈源舱最大信号接收区域。假设原抛物线上存在一点（a，b，c），将该点绕y轴旋转 角度，用旋转矩阵表示该点可表示为：

将旋转矩阵与点的列向量相乘，得到绕y轴旋转角度后的坐标为：

再将该点绕原点O旋转 角度，用旋转矩阵表示为：

将旋转矩阵与点的列向量相乘，得到旋转前后的坐标关系可表示为：

其中W代表变化前后的参数，其计算公式为:

主索节点坐标优化

将空间直角坐标系旋转后，根据公式可得需要调整的主索节点相对于新空间直角坐标系的坐标。

由上图可以得出，对于任意主索节点A， 指向球心，假设调整后的偏差量为d，按照的方向径向进行，取沿球心方向为正。将基准态时的顶点坐标A记为（x，y，z），调整后的坐标N记为，本文建立如下方程组：

由于工作态反射面为旋转抛物面，故调整后的A点进过调整后得到的点N位于工作抛物面上，于是N满足理想抛物面方程，进而得到偏差量d的表达式为：

其中h为工作抛物面顶点与基准球面的相对位置。为使得工作状态的主索网能够最接近理想抛物面，设置优化的目标函数为主索点表面与理想抛物面的均方根误差最小[5]，计算均方根误差为：

其中n代表处于照明区域内的总节点个数，为照明区域内节点的偏差分布向量，满足：

对此，建立目标函数：

对于均方差超出最大上限时，取最大上限为0.6。当h为0时（对于这一点的h值，一应该是存在误差的，还有待优化），均方差误差最小，此时的各点位于最佳拟合抛物面上，顶点坐标为（-217.0118，195.7897，-78.6778）。接着根据旋转矩阵的逆变换，得出原空间直角坐标系的顶点坐标为（-49.34，-36.8645，-230.002）。

问题三的建立与求解

接收面积

要计算调节后馈源舱的接收比，通过计算馈源舱有效区域接收到的反射信号与300米口径内反射面的反射信号之比。下图为入射信号经反射板反射后，打到馈源舱接收信号的有效区域的示意图。

对于该问，本文选择达到圆盘上的信号占整个圆盘的面积的比值作为比较值。考虑到反射面是由三个主索节点构成的三角形，视信号沿直线传播，沿垂直方向射向反射面，在此条件下以反射点为坐标原点，构建空间直角坐标系。
对任一由构成的平面 ，分别计算两点构成的向量表达式为：

对这两个向量的叉乘积进行标准化，得到该平面的单位法向量。对于平面中的P1点，入射信号的单位向量为（0，0，-1），根据入射信号与反射信号关于法线对称，可以算出射出信号的单位向量。根据单位向量与实际坐标的关系，得出：

其中，根据题目可得，。根据上述方程组，得到，
的坐标分别为：

由此得出平面一个端点射到馈源舱接收平面的入射点的x坐标和y坐标。通过选取另外两组数据，再次计算反射信号在馈源舱接收平面的入射点坐标，从而确立接收平面 的方程。

以上为计算流程，即通过不断选取反射平面，计算对于该平面的反射信号入射向量，对叉乘积进行标准化得到该平面的单位法向量，进而得出射出法向量，计算平面内的点坐标。

接收信号比

依据上文计算所得的接收平面方程，以馈源舱接收信号的有效区域为直径1米的中心圆盘为限制条件，确立中心圆盘区域方程为:

结合上述约束条件，反推至口径300米的工作抛物面对应主索节点编号，计算有效主索节点个数。因此，可以将有效区域接收到的反射信号与300米口径内反射面的反射信号之比转化为筛选主索节点与有效主索节点之比。本文通过Matlab绘出挑选节点与有效节点的对比图。

经过计算，调节后的馈源舱接收比为0.79。选取改变前的主索节点作为样本数据，同样按照上述过程计算，以2226个主索节点作为样本，得出基准球面馈源舱的接收比为0.025。（关于上图的图例和显示不符合的原因，是视角问题造成的，取的步长小了）

总结

本篇文章都是自己的一些见解，在一些地方还存在瑕疵，还请大家海涵。
对于问题二的思路，我参考了一下该篇文章：https://zhuanlan.zhihu.com/p/433577623。对于文中matlab算法的实现，如有疑问可以私我。