球面射电望远镜主动反射面的形状调节算法

Shape Regulation Algorithm for the Active Reflection Surface of the Spherical Radio Telescope

摘要: 为获得天体电磁波经反射面反射后的最佳接收效果，中国科学院研制了500米口径球面射电望远镜。该望远镜在反射面板调节约束条件下，通过调节促动器的径向伸缩量，将反射面调节为工作抛物面，使得该工作抛物面尽量贴近理想抛物面。本文以基准球面的球心为原点建立空间直角坐标系，只对被观测天体位于Z轴上方(简称为情形1)和方位角为36.7953˚、仰角为78.169˚ (简称为情形2)进行建模与分析。对于情形1，该文设计了一种算法。该算法基于馈源舱到理想抛物面与基准球面交点的长度选取照明区域的主索节点。继之以非线性最小二乘法得到被选取的主索节点与理想抛物面方程的拟合，其拟合优度R² = 0.9995，均方根误差RMSE = 0.2085，求得的拟合抛物面较为理想。对于情形2，因天体位置变化，得到的馈源舱到理想抛物面与基准球面交点的长度有误差，该文将选取照明区域的主索节点的判断条件，即馈源舱到基准球面的所有主索节点的距离是否小于等于PM (馈源舱到理想抛物面与基准球面交点的长度)进行放缩，得到与状况1相同个数的主索节点坐标。然后用非线性最小二乘法进行拟合，得到近似理想抛物面。其拟合优度R² = 0.9983，均方根误差RMSE = 0.6965，拟合效果较好。

Abstract: In order to obtain the best receiving effect of celestial electromagnetic waves after being reflected by the reflector, the Chinese Academy of Sciences has developed a 500-meter caliber sphere radio telescope. The telescope adjusts the reflection surface to a working parabolic surface by adjusting the radial elongation of the actuator, making the working parabolic surface as close to the ideal parabolic surface as possible. In this paper, the spatial right-angle coordinate system is established from the center of the benchmark sphere as the origin, which only models and analyzes the observed objects located above the Z axis (referred to as situation 1) and an azimuth of 36.7953˚ and an elevation of 78.169˚ (situation 2). For situation 1, the algorithm is designed in the paper. The algorithm selects the main cable node of the illumination area based on the length of the feed compartment to the intersection between the ideal parabolic plane and the reference sphere. Following the nonlinear least squares method to obtain the fit of the selected main cable node and the ideal parabolic equation, its goodness of fit R² = 0.9995, root mean square error RMSE = 0.2085. The obtained fitted parabola is relatively ideal. For situation 2, due to the change of the object position, the selection of the judgment condition of the main cable node of the lighting area is whether the distance of the source compartment from the main cable node to the reference space is less than or equal to the length of the PM (feed compartment to the intersection between the ideal parabolic and the reference sphere) to obtain the same number of main cable node coordinates as situation 1. It is then fitted by nonlinear least squares to obtain an approximate ideal parabola. Its goodness-of-fit R² = 0.9983, root mean square error RMSE = 0.6965, and fitting effect is good.

1. 引言

500米口径球面射电望远镜又被称为“中国天眼”，简称为“FAST”，是中国科学院国家天文台的一座射电望远镜。FAST主体工程不仅是目前世界上最大的填充口径射电望远镜 [1]；而且还是世界第二大单一口径射电望远镜，仅次于俄罗斯RATAN-600环状射电望远镜 [2]。

FAST由主动反射面、信号接收系统(馈源舱)以及相关的控制、测量和支承系统组成，其中主动反射面系统是由主索网、反射面板、下拉索、促动器及支承结构等主要部件组成的一个可调节球面。在不考虑周边支承结构连接的部分反射面板的情况下，主动反射面共有2226个主索节点，6525根节点间连接主索，4300块反射面板 [3]。

主动反射面系统有两种状态：基准态和工作态。主动反射面处于基准态时反射面为半径约300米、口径500米的球面；而工作态时通过促动器、下拉索和节点调节反射面的径向位置，形成一个口径为300米且不断移动的瞬时抛物面，此抛物面近似旋转抛物面 [4]。馈源舱接收平面的中心只能在与基准球面同心的一个球面(焦面)上移动，两同心球面的半径差为 $F=0.466R$。为使目标观测天体的平行电磁波反射汇聚到馈源舱的有效区域，FAST通过促动器的顶端沿基准球面径向伸缩调节下拉索改变反射面板的位置形成工作抛物面。且反射面调节为工作抛物面时是在反射面板调节约束下，通过调节促动器的径向伸缩量而成的 [5]。

由于小孔的直径小于所观察的天体电磁波的波长，不影响对天体电磁波的反射，因此本文可以认为面板是无孔的，同时电磁波信号及反射信号均视为直线传播。本文只对两种状况进行建模与分析：状况1：被观测天体位于Z轴上方；状况2：被观测天体位于方位角为36.7953˚，仰角为78.169˚，以基准球面的球心为原点建立空间直角坐标系。为确定一个理想的抛物面，使工作抛物面尽量贴近理想抛物面，给出了瞬时抛物面的拟合算法，以获得天体电磁波经反射面反射后的最佳接收效果，同时也给出了状况2反射面贴近理想抛物面时反射面板的调节模型以及调节相关促动器的伸缩量。

2. 理想抛物面的确定

2.1. 情形1被观测天体位于Z轴上方

2.1.1. 照明区域主索节点的提取

以基准球面的球心C为坐标原点，建立 $X,Y,Z$ 三维空间直角坐标系，球心C向被观测天体的方向为正方向。当待观测天体S位于基准球面正上方，即 $\alpha =0˚$， $\beta =90˚$，此时“FAST”天眼剖面图如图1所示：

Figure 1. Schematic diagram of section for FAST eye located at α = 0˚, β = 90˚

图1. FAST天眼位于α = 0˚，β = 90˚时剖面示意图

通过查找资料显示，中国天眼主动反射面馈源舱的照明角度为110˚~120˚。又因为反射面在工作状态时是一个300米口径的近似旋转抛物面，当 $\theta =120˚$ 时照明区域内的工作抛物面与基准球面刚好贴合，可以使天体目标的平行电磁波在反射面板上反射的有效区域面积最大，获得天体电磁波经反射面反射后的最佳接受效果。而当照明区域 ${110}^{\circ }\le \theta <{120}^{\circ }$ 时，工作抛物面在照明区域内的反射区域面积要远小于当 $\theta =120˚$ 时的反射区域面积。

因此，我们提取照明区域主索节点时取 $\theta =120˚$ 时的工作抛物面，此时它更贴近于理想抛物面，但仍存有一定的误差值。

$\begin{array}{l}\because R=300\text{m},F=0.466\text{m},MN=300\text{m}\\ \therefore OP垂直平分NM\\ \therefore PM=\frac{OP}{\cos \left(30˚\right)}=100\sqrt{3}\text{m}\\ \because CP=\left(1-0.466\right)R=160.2\text{m}\\ \therefore 馈源仓的坐标P为\left(0,0,160.2\right)\end{array}$ (1)

由式(1)便可得到提取照明区域主索节点约束条件为：

${\left(x_i-p_x\right)}^2+{\left(y_i-p_y\right)}^2+{\left(z_i-p_z\right)}^2\le P{M}^2$ (2)

其中， $x_i,y_i,z_i$ 为主索节点， $p_x,p_y,p_z$ 为馈源舱的坐标，当主索节点到馈源舱P点的距离小于等于PM的距离时，我们便认为该主索节点为照明区域的点。为提取该照明区域的主索节点，我们设计了一种算法，其流程如图2所示：

Figure 2. Algorithmic process for extracting the main cable node coordinates of the illumination area

图2. 提取照明区域的主索节点坐标算法流程

使用python实现提取照明区域的主索节点坐标算法，最终我们共提取586个主索节点坐标，照明区域如图3、图4和图5所示，图6、图7和图8为该照明区域针状图。

Figure 3. View of the illumination areas of the observed objects located at α = 0˚, β = 90˚

图3. 被观测天体位于α = 0˚，β = 90˚时的照明区域视图

Figure 4. X-Y graph of the illumination areas of the observed objects located at α = 0˚, β = 90˚

图4. 被观测天体位于α = 0˚，β = 90˚时的照明区域X-Y视图

Figure 5. X-Z graph of the illumination areas of the observed objects located at α = 0˚, β = 90˚

图5. 被观测天体位于α = 0˚，β = 90˚时的照明区域X-Z视图

Figure 6. Needle plots of the illumination areas of the observed objects located at α = 0˚, β = 90˚

图6. 被观测天体位于α = 0˚，β = 90˚时的照明区域针状图

Figure 7. Needle-like X-Z views of the illumination areas of the observed objects located at α = 0˚, β = 90˚

图7. 被观测天体S位于α = 0˚，β = 90˚时的照明区域针状X-Z视图

Figure 8. Needle-like Y-Z views of the illumination areas of the observed objects located at α = 0˚, β = 90˚

图8. 被观测天体S位于α = 0˚，β = 90˚时的照明区域针状图Y-Z视图

2.1.2. 抛物面的确定

为找到最贴近于理想抛物面的工作抛物面，我们通过抛物线方程： $y^2=2pz$ 绕Z轴旋转后形成抛物面： $x^2+y^2=2pz$，建立目标模型方程： $x^2+y^2=2pz+c$，其中c是工作抛物面贴近理想抛物面之间存在的误差 [6]。为使工作抛物面更贴近理想抛物面，我们将提取出来的照明区域主索节点与带有未知参数的抛物面进行拟合，经运算可得到参数c的值，此时所得到的方程认为最理想的工作抛物面 [7]。

建立抛物线方程

$y^2=2pz$ (3)

贴近理想抛物线方程

$y^2=2pz+c^{\prime }$ (4)

其中 $c^{\prime }$ 是贴近理想抛物面需填补的误差值。

式(4)绕Z轴旋转后得到旋转抛物面：

$x^2+y^2=2pz$ (5)

我们根据工作抛物面建立目标函数方程：

$x^2+y^2=2pz+c$ (6)

经转换得：

$z=\frac{x^2+y^2+c}{2p}$ (7)

由于 $2p=559.2$，所以：

$z=\frac{x^2+y^2+c}{559.2}$ (8)

注：工作状态时的主动反射面——300米口径的近似旋转抛物面(工作抛物面)是近似由抛物线旋转而来，所以与理想抛物面存在一定误差c。

使用MATLABcftool拟合箱中非线性最小二乘法将提取出来的主索节点与式(8)拟合得到 $c=-168300$。将c代入式(8)中可得抛物面的方程：

$z=\frac{x^2+y^2-168300}{559.2}$ (9)

2.2. 情形二被观测天体位于方位角为36.7953˚，仰角为78.169˚

2.2.1. 照明区域主索节点提取

根据对状况二的分析，我们根据题中给出的条件 $R=300\text{m}$， $F=0.466R$， $\alpha =36.7953˚$， $\beta =78.169˚$ 利用旋转坐标变换和几何对称性得到观测天体 $S\left(-x,-y,-z\right)$ 、抛物面顶点 $A\left(x,y,z\right)$ 和馈源舱 $P\left(p_x,p_y,p_z\right)$，且 $PM=100\sqrt{3}\text{m}$ (如图9所示)。

由已知条件我们便可得到照明区域主索节点的约束条件为：

$\sqrt{{\left(x_i+p_x\right)}^2+{\left(y_i+p_y\right)}^2+{\left(z_i+p_z\right)}^2}\le 100\sqrt{3}$ (10)

使用情形1提取照明区域主索节点的算法，共得到567个主索节点。由于天体位置变化后，照明区域的主索节点个数不变。且馈源舱的坐标是近似的，得到的馈源舱到理想抛物面与基准球面焦点的长度有误差，因此我们将判断条件，即馈源舱到基准球面的所有主索节点的距离是否小于等于PM进行放缩，增长PM的长度，其照明区域主索节点的约束条件为：

$\sqrt{{\left(x_i+p_x\right)}^2+{\left(y_i+p_y\right)}^2+{\left(z_i+p_z\right)}^2}\le 101\sqrt{3}$ (11)

经搜索判断我们最终得到586个主索节点，其照明区域如图10、图11和图12所示，图13、图14和图15为该照明区域针状图。

Figure 9. Schematic diagram of section for FAST eye located at α = 36.7953˚, β = 78.169˚

图9. 观测天体位于α = 36.7953˚，β = 78.169˚时FAST的剖面示意图

Figure 10. View of the illumination areas of the observed objects located at α = 36.7953˚, β = 78.169˚

图10. 被观测天体位于α = 36.7953˚，β = 78.169˚时的照明区域视图

Figure 11. X-Y graph of the illumination areas of the observed objects located at α = 36.7953˚, β = 78.169˚

图11. 被观测天体位于α = 36.7953˚，β = 78.169˚的照明区域X-Y视图

Figure 12. X-Z graph of the illumination areas of the observed objects located at α = 36.7953˚, β = 78.169˚

图12. 被观测天体位于α = 36.7953˚，β = 78.169˚的照明区域X-Z视图

Figure 13. Needle plots of the illumination areas of the observed objects located at α = 36.7953˚, β = 78.169˚

图13. 被观测天体位于α = 36.7953˚，β = 78.169˚的照明区域针状图

Figure 14. Needle-like X-Z views of the illumination areas of the observed objects located at α = 36.7953˚, β = 78.169˚

图14. 被观测天体位于α = 36.7953˚，β = 78.169˚的照明区域针状图X-Z视图

Figure 15. Needle-like Y-Z views of the illumination areas of the observed objects located at α = 36.7953˚, β = 78.169˚

图15. 被观测天体S位于α = 36.7953˚，β = 78.169˚的照明区域针状图Y-Z视图

2.2.2. 坐标系变换

为方便观察与计算，我们以C点为原点建立三维空间直角坐标系，如图16所示。通过坐标系的变换，我们可以得到观测天体 $S\left(x,y,z\right)$ 、抛物面顶点 $A\left(x,y,z\right)$ 和馈源舱 $P\left(p_x,p_y,p_z\right)$ 的坐标。

当 $H=160$ 米时，观测天体S与馈源舱P关于原点对称；当 $H\approx 300.97$ 米时，观测天体S与抛物面顶点A关于原点对称。通过式(12)坐标系变换得到抛物面顶点A和馈源舱P的坐标。

$\begin{array}{l}x=-300\cos \beta \cos \alpha \hfill \\ y=-300\sin \alpha \cos \beta \hfill \\ z=-300\sin \beta \hfill \\ P=R-0.466R=160.2\hfill \\ P_x=-P\cos \beta \cos \alpha \hfill \\ P_y=-P\sin \alpha \cos \beta \hfill \\ P_z=-P\sin \beta \hfill \end{array}$ (12)

将已知参数代入式(12)，我们得到被观测天体S的坐标为(26.3018, 19.6727, 156.7968)、抛物面顶点A的坐标为(−49.4129, −36.9588, −294.5721)和馈源舱P的坐标为(−26.3018, −19.6727, −156.7968)，保留4位小数。

Figure 16. The coordinate system of observed objects located at α = 36.7953˚, β = 78.169˚

图16. 被观测天体位于α = 36.7953˚，β = 78.169˚的坐标系

2.2.3. 抛物面的确定

对于情形二抛物面的求解，我们依然建立带有未知参数的方程(式(13))，将提取出来的586个主索节点与式(13)拟合，求解未知参数c。

$z=\frac{x^2+y^2+c}{559.2}$ (13)

使用MATLABcftool拟合箱中非线性最小二乘法将提取出来的主索节点与式(13)拟合得到c = −168,300。将c代入式(13)中可得抛物面的方程：

$z=\frac{x^2+y^2-168,100}{559.2}$ (14)

2.2.4. 带有约束条件的非线性规划方程的建立

为使反射面尽量贴近理想抛物面，我们建立了反射面板调节模型。通过调节相关促动器的伸缩量，可控制主索节点的移动变位，此时工作状态下的抛物面与理想抛物面之间的距离最短，我们认为该调节状态是最理想的 [8]。

由于连接主索节点与促动器顶端的下拉索的长度不能变，且促动器伸缩沿基准球面径向趋向球心方向为正向。在基准状态下，促动器顶端的径向伸缩量为0，其径向伸缩范围为−0.6~+0.6米。我们将此条件作为调节模型的约束条件，主索节点的移动变位作为目标方程，我们得到了带有约束条件的非线性规划方程 [9]，其方程如式(15)：

$\begin{array}{l}\min \Delta d_i=\sqrt{\Delta x_i^2+\Delta y_i^2+\Delta z_i^2}\\ \text{s}\text{.t}\text{.}\Delta d_i=\sqrt{\Delta x_i^2+\Delta y_i^2+\Delta z_i^2\le 0.6}\\ z_i+\Delta z_i=\frac{{\left(x_i+\Delta x_i\right)}^2+{\left(y_i+\Delta y_i\right)}^2-168,100}{559.2}\\ i=1,2,\cdots ,586\\ \Delta x_i,\Delta y_i,\Delta z_i\ge 0\end{array}$ (15)

$\Delta d_i$ 为第i个主索节点到理想抛物面的距离，即促动器的伸缩量； $\Delta x_i$ 、 $\Delta y_i$ 、 $\Delta z_i$ 分别为X轴、Y轴和Z轴上的位移。

我们利用python依次求得照明区域的主索节点到理想抛物面的距离，即伸缩量和位移 $\Delta x$ 、 $\Delta y$ 、 $\Delta z$。

3. 抛物面的拟合算法

非线性指的是 $f_i\left(x\right)$ 无法表示为x的线性关系，即 $f_i\left(x\right)$ 无法表示为线性关系，而是非线性关系。当 $f\left(x\right)$ 是x的非线性关系时，此时的最小二乘为非线性最小二乘。

由于非线性最小二乘问题无法构造出一个矩阵线性方程组，就没有办法求解析解，因此本问使用莱文贝格–马夸特方法(Levenberg-Marquardt algorithm)求解非线性最小二乘的解析解。

莱文贝格–马夸特方法能提供数非线性最小化(局部最小)的数值解。该算法不仅可以在执行时修改参数达到结合高斯–牛顿算法和梯度下降法的优点，而且还会对这两种方法的不足进行改善。

该方法是一个迭代的方法。我们根据泰勒展开式能把 $f\left(p+{\delta }_p\right)$ 写为下面的近似，如式(16)，这样写有两个优点：一是线性；二是只需要一阶微分。

$f\left(p+{\delta }_p\right)\approx f\left(p\right)+J{\delta }_p$ (16)

其中J是f的雅可比矩阵。对于每次的迭代我们假设这次迭代的点是 $p_k$，要找到一个 ${\delta }_{p,k}$ 让 $\left|x-f\left(p+{\delta }_{p,k}\right)\right|\approx \left|x-f\left(p\right)-J{\delta }_{p,k}\right|=\left|{ϵ}_k-{\delta }_{p,k}\right|$ 最小。根据投影公式我们知道当式(17)被满足的时候能有最小误差：

$\left(J^{\text{T}}J\right){\delta }_P=J^{\text{T}}{ϵ}_k$ (17)

将这个公式略加修改得到式(18)：

$\left[\mu I+\left(J^{\text{T}}J\right)\right]{\delta }_p=J^{\text{T}}{ϵ}_k$ (18)

当 $\mu$ 大的时候这种算法会接近最速下降法，小的时候会接近高斯–牛顿方法。为了确保每次 $ϵ$ 长度的减少，我们先采用一个小的 $\mu$，如果 $ϵ$ 长度变大就增加 $\mu$。

当这个算法达到以下条件时结束迭代：如果发现 $ϵ$ 长度变化小于特定的给定值就结束；发现 ${\delta }_p$ 变化小于特定的给定值就结束；到达了迭代的上限设定就结束 [10]。

4. 算法的精度评估

4.1. 情形1理想抛物面拟合精度评估

使用非线性最小二乘法将提取出的586个主索节点与所得到的理想抛物面做拟合，其误差平方和 $SSE=25.43$，拟合优度 $R^2=0.9995$，均方误差 $RMSE=0.2085$。从拟合结果来看，我们建立的模型取得了好的拟合效果，其拟合效果如图17、图18和图19所示。该模型有较好的稳定性与可行性，可获得天体电磁波经反射面反射后的最佳接收效果。

Figure 17. Fitted ideal parabolic view of the 586 main cable nodes extracted

图17. 提取出的586个主索节点拟合理想抛物面视图

Figure 18. Fitted ideal parabolic X-Y views of the 586 main cable nodes extracted

图18. 提取出的586个主索节点拟合理想抛物面X-Y视图

Figure 19. Fitted ideal parabolic X-Z views of the 586 main cable nodes extracted

图19. 提取出的586个主索节点拟合理想抛物面X-Z视图

Figure 20. Fitted ideal parabolic view of the 586 main cable nodes extracted

图20. 提取出的586个主索节点拟合理想抛物面视图

Figure 21. Fitted ideal parabolic X-Y views of the 586 main cable nodes extracted

图21. 提取出的586个主索节点拟合理想抛物面X-Y视图

Figure 22. Fitted ideal parabolic X-Z views of the 586 main cable nodes extracted

图22. 提取出的586个主索节点拟合理想抛物面X-Z视图

4.2. 情形二的理想抛物面拟合精度评估

使用与状况1相同的方法将提取出的586个主索节点与所得到的理想抛物面做拟合，其误差平方和 $SSE=283.8$，拟合优度 $R^2=0.9983$，均方误差 $RMSE=0.6965$。从拟合效果较好，其拟合效果如图20、图21和图22所示。我们所得到的理想抛物面可靠性较强，可获得天体电磁波经反射面反射后的最佳接收效果。

5. 结论

为获得天体电磁波经反射面反射后的最佳接收效果，本文根据FAST的结构和主动反射面板的调节机制，用提取照明区域主索节点算法将提取出的主索节点与带有未知参数的方程拟合，确定了情形1和情形2的理想抛物面。基于Python的FAST主动反射面的形状调节算法，建立了状况二反射面贴近理想抛物面时反射面板的调节模型，同时也计算出了反射面板调节后相关促动器的伸缩量。

在反射面板的调节模型中，本文只考虑了在基准状态下，促动器顶端径向伸缩量为0，其径向伸缩范围为−0.6米~+0.6米。但在实际的调节中，主索节点调节后，相邻节点之间的距离可能会发生微小变化。并且通过促动器顶端的伸缩，可控制主索节点的移动变位，连接主索节点与促动器顶端的下拉索的长度也保持不变。由于考虑的反射面板约束条件与实际调节有差距，会造成抛物面的拟合精度不够。在后续的工作中，会继续完善反射面板调节模型，提高抛物面的拟合精度，以获得更好的天体电磁波经反射面反射后的接受效果。

基金项目

本研究由湖南省衡阳市科技创新平台计划(重点实验室)项目“大数据交通应急管理实验室”资助，批准号202010041588。本研究由国家自然科学基金项目11271370资助。

NOTES

^*通讯作者。