2021年高教社杯全国大学生数学建模A题——基于几何模型的“FAST”主动反射面的形状调节

题目简述

2021年高教社杯全国大学生数学建模竞赛A题是关于“FAST”（500米口径球面射电望远镜）主动反射面形状调节的问题。

赛题要解决的问题是：在反射面板调节约束下，确定一个理想抛物面，然后通过调节促动器的径向伸缩量，将反射面调节为工作抛物面，使得该工作抛物面尽量贴近理想抛物面，以获得天体电磁波经反射面反射后的最佳接收效果。

具体问题

1. 确定理想抛物面：当待观测天体位于基准球面正上方（$\alpha =0°$ , $\beta =90°$）时，结合反射面板调节因素，确定理想抛物面的形状。

2. 调节反射面：当待观测天体位于特定位置（$\alpha =36.795°$ , $\beta =78.169°$）时，确定理想抛物面，并建立反射面板调节模型，通过调节促动器的伸缩量，使反射面尽量贴近理想抛物面。将理想抛物面的顶点坐标，以及调节后反射面 300 米口径内的主索节点编号、位置坐标、各促动器的伸缩量等结果按照规定的格式（见附件 4）保存在“result.xlsx”文件中。

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

题目附录

1. 主动反射面共有主索节点 2226 个，节点间连接主索 6525 根，不考虑周边支承结构连接的部分反射面板，共有反射面板 4300 块。基准球面的球心在坐标原点，附件 1 给出了所有主索节点的坐标和编号，附件 2 给出了促动器下端点（地锚点）坐标、基准态时上端点（顶端）的坐标，以及促动器对应的主索节点编号，附件 3 给出了 4300 块反射面板对应的主索节点编号。
2. 基准态下，所有主索节点均位于基准球面上。
3. 每一块反射面板均为基准球面的一部分。反射面板上开有许多直径小于 5 毫米的小圆孔，用于透漏雨水。由于小孔的直径小于所观察的天体电磁波的波长，不影响对天体电磁波的反射，所以可以认为面板是无孔的。
4. 电磁波信号及反射信号均视为直线传播。
5. 主索节点调节后，相邻节点之间的距离可能会发生微小变化，变化幅度不超过 0.07%。
6. 将主索节点坐标作为对应的反射面板顶点坐标。
7. 通过促动器顶端的伸缩，可控制主索节点的移动变位，但连接主索节点与促动器顶端的下拉索的长度保持不变。促动器伸缩沿基准球面径向趋向球心方向为正向。假设基准状态下，促动器顶端径向伸缩量为 0，其径向伸缩范围为-0.6～+0.6 米。
8. 天体 S 的方位可用方位角𝛼和仰角𝛽来表示（见图 5）。
   

问题分析

理想抛物面的确定

理想抛物面的方程由顶点坐标和焦点坐标决定。观测天体位于基准球面正上方（$\alpha =0°$ , $\beta =90°$），其距离FAST很远，故可以认为天体的电磁波均平行主轴射入天眼，可将馈源舱设定为焦点。当待观测天体位于基准球面正上方时，抛物面顶点位于Z轴上，通过主索节点坐标可确定顶点坐标，进而得到抛物线方程，再经旋转得到理想抛物面方程。

反射面板调节模型的建立

问题 2 分为两部分，第一部分要求在天体 S 位于位于$\alpha =36.795°$ , $\beta =78.169°$ 时，建立理想抛物面关系式，并求出顶点坐标；第二部分要求在第一部分的基础上，计算调节后反射面 300 米口径内的主索节点被调节后的位置坐标、各促动器的伸缩量等数据

第一部分

目标天体的位置发生变化，馈源仓的中心在一定范围内随之移动，始终保持“天体-球心-馈源舱”三点一线。理想抛物面的焦点坐标和顶点坐标发生变化，但是焦距始终保持不变，故抛物面整体形状保持不变｡因此，对于不同观测角度的理想抛物面，可利用坐标系变换和布尔莎模型求解，类同问题一中的理想抛物面，即可得出所求理想抛物面关系式｡

第二部分

通过促动器顶端沿基准球面径向的伸缩，可控制主索节点的移动变位，且促动器顶端径向伸缩范围为-0.6～+0.6 米。在确定理想抛物面后，为调节反射面尽量贴近理想抛物面，则需要计算基准态下主索节点距理想抛物面的径向距离，来确定促动器的伸缩量和主索节点的调节坐标。

需要注意的是，假如存在节点距离理想抛物面的径向距离超过0.6米，则此时促动器的伸缩量不足以满足要求，则需要调整理想抛物面。具体的调整方式为——不断调整顶点的位置，重新计算理想抛物面和径向距离，直至满足要求；或者，采用最优化的方法，遍历所有的顶点位置（顶点所对应促动器伸缩-0.6～+0.6 米），寻找所有促动器移动量最少的情况。

馈源舱接收比计算

通过微元对每个三角反射板进行调整，得到调整后的三角反射面板主索节点坐标，从而确定调整后三角板的位置。接着通过几何模型求取每个三角发射板的法向量，利用光的反射定理和几何公式计算待测天体电磁波经过三角面板反射之后的信号，通过投影统计反射信号经过馈源舱有效区域的量，进而计算调节后馈源舱的接收比，并与基准面的接收比进行比较｡

根据反射面调节后主索节点的坐标，确定三角反射面板的位置和偏转角度，利用光的反射定律计算反射信号方向，通过平行投影法和微元法计算馈源舱接收比。

模型假设

1. 计算过程中所有数据保留五位小数，假设导致的误差在合理范围内，对结果影响可忽略不计。
2. 收集的数据真实有效。

符号说明

模型建立与求解

（一）基于光学性质的理想抛物面确定模型

1. 抛物线方程的确定
   • 顶点坐标：当待观测天体S位于基准球面正上方时，抛物面顶点位于 $Z$ 轴上，设为$(0,0,z)$。通过Excel计算附件1中节点 $Z$ 坐标的平均值，得到顶点A的坐标为$(0,0,-300.4)$。
   • 焦点坐标：根据旋转抛物面的光学性质，将馈源舱P设定为焦点，半焦距$$\frac{1}{2}p=0.466R\text{即}p=0.932R=279.6$$
     则焦点的$Z$ 坐标为$-300.4+0.466R=-160.6$，焦点坐标为P $(0,0,-160.6)$。
   • 抛物线方程：将顶点与焦点坐标代入抛物线方程$x^2=2py$，可得理想抛物线方程为$$y=0.0017857x^2-300.4$$
2. 抛物面方程确定：理想抛物面由抛物线绕Z轴旋转180°得到，令$y=z$，$x^2=x^2+y^2$，得到理想抛物面方程$$z=0.0017857(x^2+y^2)-300.4$$
   
3. 理想抛物面的检验：考虑反射面板调节约束，促动器顶端径向伸缩范围为-0.6 ~ +0.6米。在基准球面上随机取点 $B$ 与基准面圆心 $O$ 进行连线，设定与 $Z$ 轴的夹角为$\theta$，工作态时反射面口径与基准球面半径构成等边三角形，$\theta$范围为$150^{\circ }\sim 180^{\circ }$。
   改变 $\theta$ 的范围，计算基准面和理想抛物面之间的偏移量，发现当$\theta$在$155^{\circ }\sim 165^{\circ }$时，偏移量超过0.6米，不满足条件。
   
4. 理想抛物面的调整：将理想抛物面顶点定义在$[-299.8,301.0]$之间移动，通过MATLAB计算不同顶点位置的最大、最小偏移量，得到顶点的z轴坐标位于$[-300.9,-300.5]$之间时满足条件。
   
   选取平均中心点$(0,0,-300.7)$作为顶点，得到优化后的抛物面方程$$z=0.00178207(x^2+y^2)-300.7$$
   

（二）基于布尔莎模型的反射面板调节

在天体 S 位于位于$\alpha =36.795°$ , $\beta =78.169°$ 时，建立理想抛物面关系式，并求出顶点坐标；并计算调节后反射面 300 米口径内的主索节点被调节后的位置坐标、各促动器的伸缩量等数据

布尔莎七参数

布尔莎七参数（Bursa-Wolf 参数）是一种用于坐标转换的数学模型，主要用于将一个坐标系统中的点转换到另一个坐标系统中。

布尔莎七参数的定义

布尔莎七参数包括三个平移参数、三个旋转参数和一个缩放参数，具体如下：

1. 三个平移参数（Translation Parameters）
   用于调整两个坐标系的原点位置差异

   • $T_x$：沿X轴的平移量
   • $T_y$：沿Y轴的平移量
   • $T_z$：沿Z轴的平移量

2. 三个旋转参数（Rotation Parameters）
   用于调整两个坐标系的轴向差异。这些旋转角度通常以弧度表示，且在实际应用中通常非常小

   • ${\epsilon }_x$：绕X轴的旋转角度
   • ${\epsilon }_y$：绕Y轴的旋转角度
   • ${\epsilon }_z$：绕Z轴的旋转角度

3. 一个缩放参数（Scale Factor）
   用于调整两个坐标系的尺度差异。缩放因子通常是一个接近于1的值（如1 + s），其中$s$ 表示尺度变化量

   • $s$：缩放因子

坐标转换公式

假设点 $P$ 在源坐标系中的坐标为$(X,Y,Z)$，在目标坐标系中的坐标为 $(X^{\prime },Y^{\prime },Z^{\prime })$，则布尔莎七参数的转换公式为：

$$\left[\begin{array}{c}{X}^{\prime }\\ Y^{\prime }\\ Z^{\prime }\end{array}\right]=R\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]+\left[\begin{array}{c}{T}_x\\ T_y\\ T_z\end{array}\right]+\left[\begin{array}{c}sX\\ sY\\ sZ\end{array}\right]$$

其中：

• 第一个矩阵是旋转矩阵，用于调整旋转差异。
• 第二个向量是平移向量，用于调整原点位置差异。
• 第三个向量是缩放向量，用于调整尺度差异。

$R$是欧拉角旋转矩阵，具体可参考https://blog.csdn.net/m0_49842530/article/details/144094059
假设旋转顺序为 绕Z轴旋转 ${\epsilon }_z$ → 绕Y轴旋转 ${\epsilon }_y$ → 绕X轴旋转 ${\epsilon }_x$，则旋转矩阵可以表示为：
$$R=R_x({\epsilon }_x)\cdot R_y({\epsilon }_y)\cdot R_z({\epsilon }_z)$$
其中：
$$R_x({\epsilon }_x)=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos ({\epsilon }_x)& -\sin ({\epsilon }_x)\\ 0& \sin ({\epsilon }_x)& \cos ({\epsilon }_x)\end{array}\right]$$
$$R_y({\epsilon }_y)=\left[\begin{array}{ccc}\cos ({\epsilon }_y)& 0& \sin ({\epsilon }_y)\\ 0& 1& 0\\ -\sin ({\epsilon }_y)& 0& \cos ({\epsilon }_y)\end{array}\right]$$
$$R_z({\epsilon }_z)=\left[\begin{array}{ccc}\cos ({\epsilon }_z)& -\sin ({\epsilon }_z)& 0\\ \sin ({\epsilon }_z)& \cos ({\epsilon }_z)& 0\\ 0& 0& 1\end{array}\right]$$

参数的求解

布尔莎七参数通常通过已知的控制点（即在两个坐标系中都有坐标的点）来求解。具体步骤如下：

1. 选择多个已知的控制点（至少需要三个不共线的点）。
2. 建立坐标转换方程组。
3. 使用最小二乘法或其他优化方法求解七个参数。

参数的求解过程详细解释

1. 选择多个已知的控制点
   控制点是那些在两个坐标系中都有已知坐标的点。这些点作为转换的参考，用于建立坐标转换方程组。通常，选择的控制点数量越多，求解的精度越高，但至少需要三个不共线的点来求解布尔莎七参数。

2. 建立坐标转换方程组
   对于每个控制点 $P_i$，假设其在源坐标系中的坐标为 $(X_i,Y_i,Z_i)$，在目标坐标系中的坐标为 $(X_i^{\prime },Y_i^{\prime },Z_i^{\prime })$。根据布尔莎七参数的坐标转换公式，可以写出以下方程组：
   $$\left[\begin{array}{c}{X}_i^{\prime }\\ Y_i^{\prime }\\ Z_i^{\prime }\end{array}\right]=R\left[\begin{array}{c}{X}_i\\ Y_i\\ Z_i\end{array}\right]+\left[\begin{array}{c}{T}_x\\ T_y\\ T_z\end{array}\right]+\left[\begin{array}{c}s{X}_i\\ s{Y}_i\\ s{Z}_i\end{array}\right]$$
   对于 $n$ 个控制点，可以得到 $3n$ 个方程。这些方程可以写成矩阵形式：
   $$\mathbf{Ax}=\mathbf{b}$$
   其中

   • $\mathbf{A}$ 是系数矩阵
   • $\mathbf{x}$ 是参数向量（包含 $T_x,T_y,T_z,{\epsilon }_x,{\epsilon }_y,{\epsilon }_z,s$）
   • $\mathbf{b}$是观测向量（包含所有控制点的目标坐标与源坐标的差值）

3. 使用最小二乘法或其他优化方法求解七个参数
   由于方程组通常是超定的（即方程数量多于未知数数量），直接求解可能没有精确解。因此，通常使用最小二乘法来求解参数向量 $\mathbf{x}$，使得残差平方和最小。最小二乘解可以通过以下公式求得：
   $$\mathbf{x}=({\mathbf{A}}^T\mathbf{A}{)}^{-1}{\mathbf{A}}^T\mathbf{b}$$
   其中，${\mathbf{A}}^T$ 是 $\mathbf{A}$ 的转置矩阵，$({\mathbf{A}}^T\mathbf{A}{)}^{-1}$ 是 ${\mathbf{A}}^T\mathbf{A}$ 的逆矩阵。

4. 注意事项

   • 旋转矩阵的非线性
     • 旋转矩阵 $R$ 是非线性的，因为它依赖于旋转角度 ${\epsilon }_x,{\epsilon }_y,{\epsilon }_z$的三角函数。因此，直接使用最小二乘法可能需要迭代求解。
     • 通常需要将非线性问题转化为线性问题，或者使用非线性最小二乘法（如 Levenberg-Marquardt 算法）来求解。
   • 初始值的重要性
     对于大角度旋转，初始值的选择对求解结果有很大影响。如果初始值不合适，可能导致求解失败或收敛到错误的解。

理想抛物面求解模型的建立

1. 理想抛物面求解模型的分析：由于目标天体平行电磁波入射时焦距不变，理想抛物面形状不变。在坐标系合理变换下，利用问题一所得方程，通过布尔莎模型将新坐标系中的方程转换为原坐标系中的方程。

2. 基于坐标转换的新坐标系的建立：将圆心C设为原点，通过偏转原坐标系使抛物面对称轴落在(z)轴上，且保证抛物面顶点距圆心距离不变，这样问题一所得方程中的常数项无需改变。
   考虑旋转顺序为 绕Z轴旋转 ${\epsilon }_z$ → 绕Y轴旋转 ${\epsilon }_y$ → 绕X轴旋转 ${\epsilon }_x$，则相对应的布尔莎模型的参数为

   • $T_x=0$
   • $T_y=0$
   • $T_z=0$
   • ${\epsilon }_x=0$
   • ${\epsilon }_y=\beta -90^{\circ }$
   • ${\epsilon }_z=-\alpha$
   • $s=0$
     

3. 理想抛物面坐标的求解：以圆心C为原点建立坐标系，根据$\alpha$、$\beta$确定问题二中抛物面顶点坐标。通过MATLAB计算角度关系和转换矩阵，得到$\alpha =36.759^{\circ }$，$\beta =78.169^{\circ }$时的理想抛物面表达式，并代入问题一方程得到在原坐标系中的方程。
   假设问题一中的理想抛物线的坐标为$(X,Y,Z)$，问题二中的理想抛物线的坐标为$(X^{\prime },Y^{\prime },Z^{\prime })$，我们尝试将问题二中的理想抛物线通过旋转与问题一中理想抛物线重合，我们根据旋转，获取$(X,Y,Z)$与$(X^{\prime },Y^{\prime },Z^{\prime })$之间转换关系

$$\begin{gathered} \left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos ({\epsilon }_x)& -\sin ({\epsilon }_x)\\ 0& \sin ({\epsilon }_x)& \cos ({\epsilon }_x)\end{array}\right]\left[\begin{array}{ccc}\cos ({\epsilon }_y)& 0& \sin ({\epsilon }_y)\\ 0& 1& 0\\ -\sin ({\epsilon }_y)& 0& \cos ({\epsilon }_y)\end{array}\right]\left[\begin{array}{ccc}\cos ({\epsilon }_z)& -\sin ({\epsilon }_z)& 0\\ \sin ({\epsilon }_z)& \cos ({\epsilon }_z)& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{X}^{\prime }\\ Y^{\prime }\\ Z^{\prime }\end{array}\right]+\left[\begin{array}{c}{T}_x\\ T_y\\ T_z\end{array}\right]+\left[\begin{array}{c}s{X}^{\prime }\\ s{Y}^{\prime }\\ s{Z}^{\prime }\end{array}\right] \\ =\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}\cos (-11.831^{\circ })& 0& \sin (-11.831^{\circ })\\ 0& 1& 0\\ -\sin (-11.831^{\circ })& 0& \cos (-11.831^{\circ })\end{array}\right]\left[\begin{array}{ccc}\cos (-36.795^{\circ })& -\sin (-36.795^{\circ })& 0\\ \sin (-36.795^{\circ })& \cos (-36.795^{\circ })& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{X}^{\prime }\\ Y^{\prime }\\ Z^{\prime }\end{array}\right] \\ =\left[\begin{array}{ccc}0.784& 0.587& -0.198\\ -0.599& 0.800& 0\\ 0.158& 0.119& 0.980\end{array}\right]\left[\begin{array}{c}{X}^{\prime }\\ Y^{\prime }\\ Z^{\prime }\end{array}\right] \end{gathered}$$

将所得$(X,Y,Z)$ 带入问题一所得方程：$z=0.00178207(x^2+y^2)-300.7$，得到 $\alpha =36.759^{\circ }$，$\beta =78.169^{\circ }$ 时的理想抛物面表达式：
$$z=3.823x+2.859y-6674.7\times \sqrt{-1.916\times 10^{-7}{x}^2+5.122\times 10^{-7}xy-3.424\times 10^{-7}{y}^2+0.0008753y+1.048+6352.9}$$

4. 理想抛物面点坐标对应自变量的范围求解：通过对理想抛物面在XOY坐标系投影椭圆方程的求解，确定$x$、$y$的范围。
5. 检验所得关系式准确程度：将所得理想抛物面顶点坐标代入表达式，验证在不同坐标系中的坐标误差小于$0.0001\%$，证明关系式准确。

基于几何模型的相关调节数据的计算

• 利用几何模型计算主索节点调节后的坐标：通过联立直线方程与抛物面方程，计算主索节点调节后的坐标。
  
  求解公式，获取最终调节量
  

• 判断调节后主索节点是否在为得到理想抛物面需调节区域内：根据理想抛物面点坐标对应自变量的范围进行初步判断，再结合主索节点到圆心距离与调节后点到圆心距离的差值是否满足(\leq 0.6)米的条件，确定最终需调节的主索节点。

（三）基于光反射定律和几何模型的馈源舱接收比的计算

1. 主索节点的调节与三角反射面板转动的关系：反射面调节通过下拉索与促动器配合，根据主索节点调节后的坐标确定三角反射面板的准确位置和调节后的角度。
2. 三角反射板法向量计算：已知三角板上三个主索点坐标，利用向量叉乘计算平面法向量。
3. 馈源舱接收比的计算
   • 光的反射定律：反射光线与入射光线、法线满足“三线共面，两线分居，两角相等”，光具有可逆性。
   • 空间中关于直线的对称问题：通过建立方程组求解空间中直线关于平面法线对称的直线方向向量。
   • 利用微元法、平行投影法和几何模型计算馈源舱接收比：利用平行投影法将三维问题转化为二维问题，对馈源舱接收信号的有效区域投影，将反射信号投射到同一平面。采用微元法在特定区域撒点，判断点是否在三角形与圆的重合部分，计算重合面积，进而得到馈源舱接收比。最终得到调整后抛物面的接收比为(1.55805%)，基准球面的接收比为(0.00042594%)。

七、模型评价

1. 模型优点：综合考虑多种约束条件，寻得最优参数，保证了结果的正确性和合理性，且对模型结果进行检验，误差小于(0.0001%)。程序运行时间快，空间占用小。
2. 模型缺点：使用大量数学公式，计算难度较高。

八、参考文献

[1] 王振德.现代科技百科全书:[M].桂林:广西师范大学出版社
[2] 南仁东, 姜鹏. 500m口径球面射电望远镜(FAST)[J]. 机械工程学报, 2017(17):1 - 3.
[3] 张立云. 500m口径球面射电望远镜对贵州发展推动[J]. 中国新技术新产品, No.311(01)
[4] 韦钺,马文双,李明君,等.工程测量中平面坐标转换软件设计及应用[J].测绘工程,2013,22(4):76 - 79.
[5] 郭英起,唐彬,张秋江,等.基于空间直角坐标系的高精度坐标转换方法研究[J].大地测量与地球动力学,2012,32(3):125 - 128.