一定约束下主动反射面的形状调节

Shape Adjustment of Active Reflector under Certain Constraints

摘要: 利用抛物面型射电望远镜观测天体时，往往需要调节主动反射面变为工作抛物面，并让该工作抛物面尽量贴近理想抛物面，以使馈源仓接收到天体反射电磁波的效果最佳。本文详细讨论了在一个特殊角度下，建立单目标最优化模型来确定理想抛物面，再通过空间坐标的旋转求出任意角度下的理想抛物面。分析了在一定约束下，通过调节促动器的伸缩量，使工作抛物面逼近理想抛物面。具体讲述了馈源仓接收反射信号的接收方式、反射信号及有效信号接收量的计算、定义了馈源仓的接收比的计算方法。

Abstract: When using a parabolic radio telescope to observe celestial bodies, it is often necessary to adjust the active reflector to become a working paraboloid, and make the working paraboloid as close as pos-sible to the ideal paraboloid, so that the feed bin can receive the best effect of the reflected electro-magnetic waves from the celestial body. This paper discusses in detail the establishment of a sin-gle-objective optimization model under a special angle to determine the ideal paraboloid, and then obtains the ideal paraboloid at any angle through the rotation of the space coordinates. It is ana-lyzed that the working paraboloid can approach the ideal paraboloid by adjusting the expansion and contraction of the actuator under certain constraints. Specifically, the receiving method of the reflected signal received by the feed bin, the calculation of the reflected signal and the effective sig-nal received, and the calculation method of the receiving ratio of the feed bin are defined.

1. 引言

建立在中国西南部贵州平塘县天然喀斯特洼地台址的FAST是目前国际上测量最精密、口径最大、观测精度、灵敏度最高的单口径射电望远镜，是我国自主研发且具有自主知识产权的500米口径宽球面射电望远镜，在推动我国原始科技创新领域的发展，天文学、基础科学的发展等方面具有重大意义 [1]。抛物面型射电望远镜FAST主要是由主动反射面与馈源舱两大部分组成，其中主动反射面系统是由主索网、反射面板、下拉索、促动器及支承结构等主要部件构成的一个可调节球面，并且有基准态和工作态两个状态 [2]。基准态时，FAST为反射面为半径约300米、口径为500米的基准球面；工作态时反射面的形状被调节为一个300米口径的近似旋转抛物面称为工作抛物面 [3]。由于受反射面板的调节约束，需要确定一个理想抛物面，通过调节促动器的径向伸缩量，将让该工作抛物面尽量贴近理想抛物面 [4]，以使馈源仓接收到的天体反射电磁波效果最佳，馈源舱接收平面的中心只能在与基准球面同心的球面即焦面上移动 [5]。对理想抛物面的确定，在已有的研究中，李明辉通过建立主索节点变位数值模型，分析了不同形变策略下组索节点径向、经向、维向的位移量 [6]；薛建兴对FAST瞬时抛物面拟合精度进行了研究，为反射单元初始面形确定及瞬时抛物面拟合精度调节提供了参考 [7]；杜敬利提出的主索网主动反射面的形状调节算法对抛物面顶点位置及非照明部分的节点位移调整进行了研究 [8] 等这些研究都为FAST的科学发展奠定了基础，但都未具体分析在任意角度下FAST瞬时抛物面变化及其表达式。

本文主要从理想抛物面的确定、工作抛物面的调节、馈源仓接收比的计算三个方面对FAST的工作原理进行讨论。首先，通过建立多目标最优化模型，确定了特定角度下的理想抛物面，再通过欧拉旋转定理求出了任意角度下的理想抛物面。然后在一定约束条件下，调节工作抛物面使其逼近理想抛物面，具体分析如何求出调节后各主索节点的坐标。最后，根据该坐标确定了所有参与变换的反射面板，并建立相关模型得出馈源仓的接收比。

2. 理想抛物面的确定

2.1. 模型的建立

假设观测天体T位于方位角和仰角处，为求出一般理想抛物面的方程，先确定，即观测天体位于FAST正上方时理想抛物面方程，再通过欧拉旋转定理得出一般情况下的抛物面方程。下面建立优化模型求出观测天体T位于正上方时的理想抛物面方程。

基准球面的球心C为原点坐标，馈源舱接收面的中心与直线TC的交点P (馈源仓平面中心)为主动反射面工作态时工作抛物面的焦点。基准球面与焦球面的半径差为 $F=0.466R$， $R=300$ 米，h为理想工作抛物面下顶点A到基准球面下顶点的距离。平面工作抛物线的示意图如下所示：

Figure 1. Parabola diagram

图1. 抛物线示意图

由图1可知，抛物线焦距 $\left|AP\right|=F+h$，顶点坐标 $A\left(0,-\left(R+h\right)\right)$，易求出抛物线方程为：

$f_1\left(x,z\right):x^2=4\left(F+h\right)z+4\left(F+h\right)\left(R+h\right)$ (1)

由上式可知，只需确定h的值即可求出抛物线方程。为使由基准球面到工作抛物面的变化较小，要考虑到理想抛物面应该尽量贴近基准球面，即使抛物线与圆弧围成的面积越小越好。

工作抛物线与基准圆弧在 $\left[-150,150\right]$ 范围内围成的面积为S，圆弧的方程为：

$f_2\left(x,z\right):x^2+z^2=R^2$

联立上式与式(1)可求出S的表达式为：

$S={\displaystyle {\int }_{-150}^{150}\left|f_1\left(x,z\right)-f_2\left(x,z\right)\right|\text{d}x}$ (2)

为保证基准球面内，边缘的抛物面成型的完整性，抛物面边缘与基准球面要满足平滑过渡 [6]，影响抛物面与基准球面相切的平滑程度的因素有边缘相切位置处主索节点的移动变化量与相切处的导数差值。

假设位于口径300米的基准球圆弧面上有一主索节点 $B_1\left(x_{B1},y_{B1}\right)$。抛物线上与之对应的一点 $B_2\left(x_{B2},y_{B2}\right)$， $x_{B2}=-150$，则基准圆弧在 $B_1$ 点的切线斜率 $k_1$，抛物线在 $B_2$ 点的切线斜率 $k_2$ 为：

$k_1{}_1=\frac{x_{B1}}{\sqrt{R^2+x_{B1}^2}},k_2=\frac{x_{B2}}{2\left(F+h\right)}$

两切线的斜率差为：

$\Delta k=\left|k_2-k_1\right|$ (3)

由于促动器是径向位移，所以需考虑边沿主索节点的径向位移量 $\Delta L$ 最小。直线 $C{B}_2$ 与基准圆弧的交点即为 $B_1$，容易求出 $B_1$ 的坐标为：

$x_{B1}=-\sqrt{\frac{R^2{x}_{B2}^2}{y_{{}_{B2}}^2+x_{B2}^2}},y_{B1}=-\frac{y_{B2}}{x_{B2}}\sqrt{\frac{R^2{x}_{B2}^2}{y_{{}_{B2}}^2+x_{B2}^2}}$

又因为主索节点的位移是通过促动器的伸缩来完成的，而促动器的伸缩量不能超过0.6米，由距离公式得到位移量 $\Delta L$ 的表达式为：

$\Delta L=\sqrt{{\left(x_{B1}-x_{B2}\right)}^2+{\left(y_{B1}-y_{B2}\right)}^2}\le 0.6$ (4)

由于出现多个目标函数，需要对各因素进行加权计算 [9]，通过分析得出影响抛物面的三个主要因素的权重比为： $S:\Delta k:\Delta L=0.9:0.02:0.08$，由此可得到关于h的单目标优化模型为：

$\min W=0.9S+0.02\Delta k+0.08\Delta L$

$\text{s}\text{.t}\text{.}\{\begin{array}{l}S={\displaystyle {\int }_{-150}^{150}\left|f_1\left(x,z\right)-f_2\left(x,z\right)\right|}\\ \Delta k=\left|k_2-k_1\right|\\ \Delta L=\sqrt{{\left(x_{B1}-x_{B2}\right)}^2+{\left(y_{B1}-y_{B2}\right)}^2}\le 0.6\end{array}$

2.2. 模型的求解

基于遍历求解法通，过MATLAB软件来进行求解并找出最优值。为确保结果的精度，以0.001为步长对h进行遍历 [10]，绘制出目标函数和h的关系图为图2：

Figure 2. The image of the objective function with respect to h

图2. 目标函数关于h的图像

求得 $h=0.368$ 米时为最优解，将h代入(1)式并绕Z轴旋转一周可得理想抛物面方程为：

$x^2+y^2=560.6720z+168407.9273$ (5)

式(5)即为 $\alpha =0˚$ 、 $\beta =90˚$ 时的理想抛物面方程，其顶点坐标为 $\left(0,0,-300.368\right)$，下面根据欧拉旋转定理求解一般情况下的抛物面方程。

定理1 [11] 欧拉旋转定理阐述为在三维空间中，任意两个共原点的三维空间坐标系，两坐标之间可以进行绕含原点的固定轴旋转，旋转后的结果为旋转前的坐标与相应旋转矩阵的乘积，绕各轴旋转的旋转矩阵如下：

$R_Y\left(\alpha \right)=\left[\begin{array}{ccc}\cos \alpha & 0& \sin \alpha \\ 0& 1& 0\\ -\sin \alpha & 0& \cos \alpha \end{array}\right]\begin{array}{cc}& R_Z\left(\alpha \right)=\left[\begin{array}{ccc}\cos \alpha & -\sin \alpha & 0\\ \sin \alpha & \cos \alpha & 0\\ 0& 0& 1\end{array}\right]\end{array}\begin{array}{cc}& R_X\left(\alpha \right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \alpha & -\sin \alpha \\ 0& \sin \alpha & \cos \alpha \end{array}\right]\end{array}$

现在天体T位于方位角 $\alpha$ 和仰角 $\beta$ 处，如图3所示。则此时的抛物面可以看成原抛物面( $\alpha =0˚$ 、 $\beta =90˚$ 时的理想抛物面)先绕Z轴逆时针旋转角，再绕Y轴逆时针旋转 $\omega =\pi /2-\beta$ 角得到。设原抛物面上一点 $P_0\left(x_0,y_0,z_0\right)$，旋转后抛物面对应点为 $P_0\left(x_0,y_0,z_0\right)$，则有：

$C{P}_0=R_Z\left(\alpha \right)R_Y\left(\omega \right)C{P}_1$

Figure 3. Azimuth and elevation diagram

图3. 方位角与仰角示意图

由上式可以求得旋转后抛物面的方程，化简后可得：

${\phi }_1^2\left(x,y,z\right)+{\phi }_2^2\left(x,y,z\right)-{\rho }_1={\rho }_2{\phi }_3\left(x,y,z\right)$ (6)

方程中 ${\rho }_1=161631.873$， ${\rho }_2=537.394$，其它 ${\phi }_1\left(x,y,z\right)$ 、 ${\phi }_2\left(x,y,z\right)$ 、 ${\phi }_3\left(x,y,z\right)$ 分别为：

$\{\begin{array}{l}{\phi }_1\left(x,y,z\right)=x\cos \alpha \cos \beta -z\sin \beta +y\cos \beta \sin \alpha \\ {\phi }_2\left(x,y,z\right)=x\sin \alpha -y\cos \alpha \\ {\phi }_3\left(x,y,z\right)=x\cos \alpha \sin \beta +z\cos \beta +y\sin \beta \sin \alpha \end{array}$

式(6)从形式上看，与式(1.5)比较接近，在一定条件下也可以看作是由抛物线旋转得到。此式即为一般情况下，观测天体T位于方位角 $\alpha$ 和仰角 $\beta$ 处时理想工作抛物面方程，该抛物面顶点坐标求出为：

$-300.368\left(\cos \alpha \cos \beta ,\cos \beta \sin \alpha ,\sin \beta \right)$

理想抛物面的确定为后面实际工作抛物面的调节奠定了基础。确定理想工作抛物面的方程后，需要通过调节促动器使实际工作抛物面尽量逼近该理想抛物面。

3. 工作抛物面的调节

3.1. 主索节点的确定

工作抛物面的调节是由下拉索和促动器配合完成的，下拉索长度固定，促动器沿基准球面径向安装，底端固定但顶端可伸缩，从而可调节反射面板，形成工作抛物面 [10]。基准态时的主索节点、促动器上下端点坐标已知，并且500米口径的基准球面上共有主索节点2226个，反射面板4300块，而工作抛物面的口径为300米，所以为了能够调节工作抛物面，我们需要确定所有位于工作抛物面上主索节点，再通过每个组索节点所对应的促动器伸缩，使抛物面得到调节。

在图3中，由于抛物面的顶点坐标在直线TC上，故直线TC方向向量容易求出为： $n=\left(\cos \alpha \cos \beta ,\cos \beta \sin \alpha ,\sin \beta \right)$，任取基准球上一点 $Q\left(x_1,y_1,z_1\right)$，易知直线QC的方向向量为： $n_1=\left(x_1,y_1,z_1\right)$，则有点Q到SC的距离为：

$d=\frac{n{n}_1}{\left|n\right|}$

将所有组索节点的坐标点代入上式，且由 $d\le 150$，即可确定工作抛物面上的主索节点。

3.2. 促动器伸缩量的计算

抛物面调节的本质是促动器的伸缩，因此需要计算出促动器的伸缩量。设促动器上端点为 $O_1\left(x_{O_1},y_{O_1},z_{O_1}\right)$，下端点为 $O_2\left(x_{O_2},y_{O_2},z_{O_2}\right)$，主索节为 $O_3\left(x_{O_3},y_{O_3},z_{O_3}\right)$，则可求出三点所在直线的方程为：

$\frac{x-x_{O_3}}{x_{O_1}-x_{O_2}}=\frac{y-y_{O_3}}{y_{O_1}-y_{O_2}}=\frac{z-z_{O_3}}{z_{O_1}-z_{O_2}}$

将上式与(1.6)式理想抛物面方程联立，即可求出该主索节点经促动器径向伸缩后落在理想抛物面上的坐标为 $O_4\left(x_{O_4},y_{O_4},z_{O_4}\right)$，伸缩量为：

$L=\sqrt{{\left(x_{O_4}-x_{O_3}\right)}^2+{\left(y_{O_4}-y_{O_3}\right)}^2+{\left(z_{O_4}-z_{O_3}\right)}^2}$

若求出 $\left|L\right|\le 0.6$，则满足要求；若 $\left|L\right|>0.6$，则不满足要求，此时需将伸缩量固定为0.6，再根据上式反求出变换后的主索节点坐标，该坐标即为工作抛物面上的该组索节点的坐标，与理想情况有一定差距。

4. 馈源仓接收比计算

4.1. 反射直线的确定

由于三维空间具有各向同性的特点，而上面求出的工作抛物面是以图3中TC为中心轴，抛物面上的主索节点坐标可认为是在OXYZ标系下的坐标。通过旋转变换将这些坐标转化为 $O{X}^{\prime }{Y}^{\prime }{Z}^{\prime }$ 坐标系下的坐标，即此时抛物面以 $Z^{\prime }$ 轴为中心轴。如此，可使信号馈源仓接收比的计算更加方便。在 $O{X}^{\prime }{Y}^{\prime }{Z}^{\prime }$ 坐标系下信号的入射方向为负 $Z^{\prime }$ 轴，设变换后任意一点坐标为 $H_1\left(x_{H_1},y_{H_1},z_{H_1}\right)$，即可求出过该点的入射直线方程为：

$\frac{x-x_{H_1}}{0}=\frac{y-y_{H_1}}{0}=\frac{z-z_{H_1}}{1}$

$H_1$ 为反射面板上的一点，设与 $H_1$ 同一块反射板上另外两点的坐标为 $H_2\left(x_{H_2},y_{H_2},z_{H_2}\right)$ 、 $H_3\left(x_{H_3},y_{H_3},z_{H_3}\right)$，由于反射面板相对整个FAST系统较小，其板面弯曲程度可忽略，因此可将反射面板可近似为平面板，从而可求出该平面板的法向量为：

$n_2=H_1{H}_2\cdot H_1{H}_3=\left|\begin{array}{ccc}{e}_x& e_y& e_z\\ x_{H_2}-x_{H_1}& y_{H_2}-y_{H_1}& z_{H_2}-z_{H_1}\\ x_{H_3}-x_{H_1}& y_{H_3}-y_{H_1}& z_{H_3}-z_{H_1}\end{array}\right|=a{e}_x-b{e}_y+c{e}_z=\left(a,b,c\right)$

将 $n_2=\left(a,b,c\right)$ 单位化后的到 $n_3=\left(a^{\prime },b^{\prime },c^{\prime }\right)$，其中 $a^{\prime }$ 、 $b^{\prime }$ 、 $c^{\prime }$ 分别为：

$a^{\prime }=\frac{a}{\sqrt{a^2+b^2+c^2}},b^{\prime }=\frac{b}{\sqrt{a^2+b^2+c^2}},c^{\prime }=\frac{c}{\sqrt{a^2+b^2+c^2}}$

反射面板上，入射直线与反射直线示意图如下所示：

Figure 4. Schematic diagram of the incident line and the reflection line

图4. 入射直线与反射直线示意图

图4中，若 $\left|MH\right|=1$，则 $MH=\left(0,0,1\right)$， $\cos \theta =l_1{n}_3/\left|l_1\right|\left|n_3\right|$， $HR=\left(a^{\prime }\cos \theta ,b^{\prime }\cos \theta ,c^{\prime }\cos \theta \right)$，由向量运算法则容易求出：

$HN=MH+2HR=\left(2a^{\prime }\cos \theta ,2b^{\prime }\cos \theta ,2c^{\prime }\cos \theta +1\right)$

式中， $HN$ 即为反射直线的方向向量，由此即可求出反射直线的方程为：

$\frac{x-x_{H_1}}{2a^{\prime }\cos \theta }=\frac{y-y_{H_1}}{2b^{\prime }\cos \theta }=\frac{z-z_{H_1}}{2c^{\prime }\cos \theta +1}$ (7)

4.2. 接收比计算

在 $O{X}^{\prime }{Y}^{\prime }{Z}^{\prime }$ 坐标系下，过馈源舱中心点与焦面所切的平面与 $O{X}^{\prime }{Y}^{\prime }$ 平面平行，故可得该切平面的方程为： $Z^{\prime }=R-F$，联立该方程与式(7)即可求出在焦平面上的反射点，反射面板三个顶点所对应的反射线与切平面相交于三点，这三个反射点可围成一个三角形，该三角形即为信号经过反射面板后反射到切平面方向上的投影。故该三角形的面积可称为反射信号的面积，而该三角形与馈源仓舱接收信号的有效区域相交的面积即为馈源舱收到的有效信号面积。

馈源仓舱接收信号的有效区域，是中心点坐标在 $P_0\left(0,0,R-F\right)$ 处，且半径r为0.5的圆盘。由于三角形三个点的坐标可以求出，所以三边直线方程也可求出，用点到直线的距离公式求出中心点到三边的距离分别为 $d_1$ 、 $d_2$ 、 $d_3$，比较其与r的大小关系，确定圆与三角形的位置。考虑到投影三角形与馈源舱的有效接收区域圆的交面积存在多种情况，且圆的半径较小而投影三角形的面积较大，故仅对以下几种情况进行讨论。

1) 接收圆与反射投影三角形无相交

若 $d_1>r$ 、 $d_2>r$ 、 $d_3>r$，说明接收圆包含在三角形内，相交面积为： $S_1=\pi r^2$ ；

若 $d_1\le r$ 、 $d_2\le r$ 、 $d_3\le r$ (至多一个取等号)，说明接收圆在三角形外，此时相交面积为： $S_1=0$ ；

若 $d_1$ 、 $d_2$ 、 $d_3$ 至少有一个大于该边所对应的三角形的高与圆半径r之和 $\Delta H$，则 $S_1=0$。

2) 接收圆与反射投影三角形相交一边

若 $d_1$ 、 $d_2$ 、 $d_3$ 中只有 $d_1\le r$ 时，且 $d_2<\Delta H$ 、 $d_3<\Delta H$，说明接收圆与三角形一边相交，如图5左图所示。联立直线与圆的方程即可求出弦长 $\Delta I$，再由余弦定理可求得：

$\cos \beta =\frac{2r^2-\Delta I^2}{2r^2},S_{扇}=\frac{\beta }{\text{2}}{r}^2$

此时，若圆心 $P_0$ 在三角形内部时，有效接收面积为： $S_1=\pi r^2-S_{扇}+d_1\Delta I/2$ ；若圆心在三角形外部，则有效接收面积为： $S_1=\pi r^2-S_{扇}+d_1\Delta I/2$。

Figure 5. Intersection diagram

图5. 相交示意图

3) 接收圆与反射投影三角形相交两边

若 $d_1$ 、 $d_2$ 、 $d_3$ 中有两个小于r，另一个小于 $\Delta H$，说明接收圆与三角形相交两条边。联立圆与两条直线的方程即可求出圆所截三角形两条边 $I_1$ 、 $I_2$ 的长度。相交部分可近似看成以三角形顶点为圆心的扇形，其半径近似为 $r^{\prime }=\left(I_1+I_2\right)/2$，而角 $\omega$ 可由余弦定理求出，故可得有效接收面积为： $S_1=\omega {\left(r^{\prime }\right)}^2/2$。

上面分析了馈源仓接收到的有效面积，而反射的总面积即为投影三角形的面积 $S_2$，由于三角形三点坐标已求出，每条边的长度也可求出，根据海伦公式容易求出 $S_2$。

第二节中已经确定了所有参与工作的组索节点，将这些点坐标与每块反射面板的三个顶点坐标进行比对，则可确定参与工作的反射面板，假设所有工作的反射面板共有N块，累加这N块反射面板反射三角形的面积即为反射信号总的信号量，累加这N块反射面板所反射到馈源仓有效区域的面积即为接收到的总的信号量，故接收比可定义为两者的比值，计算式为

$\eta =\frac{{\displaystyle \underset{i=1}{\overset{N}{\sum }}{S}_1^i}}{{\displaystyle \underset{i=1}{\overset{N}{\sum }}{S}_2^i}}$ (8)

5. 实例分析

当天体的方向为 $\alpha =36.795˚$ 、 $\beta =78.169˚$ 时，由式(6)可求出理想工作抛物面的方程为：

${\phi }_1^2\left(x,y,z\right)+{\phi }_2^2\left(x,y,z\right)-{\rho }_1={\rho }_2{\phi }_3\left(x,y,z\right)$

$\{\begin{array}{l}{\phi }_1\left(x,y,z\right)=0.1641x-0.9788z+0.1228y\\ {\phi }_2\left(x,y,z\right)=0.5990x-0.8008y\\ {\phi }_3\left(x,y,z\right)=0.7838x+0.2050z+0.5862y\\ {\rho }_1=161631.873,{\rho }_2=537.394\end{array}$

顶点坐标也可求得为： $\left(-49.2904,-36.8852,-294.0001\right)$，在各组索节点坐标已知的情况下，根据第二节的具体分析，确定了工作抛物面上的组索节点共有692个，通过MATLAB可绘制出分布示意图如图6：

Figure 6. Workgroup index node distribution diagram

图6. 工作组索节点分布示意图

通过促动器对这692个组索节点进行调节，使工作抛物面逼近理想抛物面，可得到调节后组索节点坐标及相应促动器伸缩量见表1 (有代表性的少量数据)：

Table 1. Original coordinates, adjusted coordinates, actuator telescopic amount

表1. 原坐标、调节后坐标、促动器伸缩量

表1中可看出，节点1的伸缩量为正值，节点2的伸缩量为负值，节点3的伸缩量刚好为−0.6，说明按理想抛物面调节其伸缩量可能超过−0.6，而将其固定在−0.6，这几组数据代表了可能出现的几种情况。

抛物面调节好后，通过第三节的分析可以计算出馈源仓的接收比为 ${\eta }_1=5.59\%$，同理可以算出未调节即基准状态下馈源仓的接收比为 ${\eta }_2=2.5\%$，提升了2.2326倍。

6. 结束语

本文对于FAST理想工作抛物面的确定，考虑了工作态和基准态的贴合程度以及边缘光滑程度两方面因素，使得求解出的抛物面方程更为准确。在理想抛物面的确定过程中，利用抛物面和球面各向同性的特点，将其由三维空间计算简化为二维平面计算，计算较为方便。同时，对于馈源仓接收比的计算通过坐标选择的方法，减轻了计算难度。在对于实际模型的变换策略中，可通过对于多个主索节点的变换形式进行分析，建立单个节点的变换对周围节点的影响模型，从而得出更为准确的主索节点位置坐标以及对应促动器的伸缩量。