2021年高教社杯全国大学生数学建模A题——基于几何模型的“FAST”主动反射面的形状调节
文章目录
- 题目简述
- 具体问题
- 题目附录
- 问题分析
- 理想抛物面的确定
- 反射面板调节模型的建立
- 第一部分
- 第二部分
- 馈源舱接收比计算
- 模型假设
- 符号说明
- 模型建立与求解
- (一)基于光学性质的理想抛物面确定模型
- (二)基于布尔莎模型的反射面板调节
- 布尔莎七参数
- 布尔莎七参数的定义
- 坐标转换公式
- 参数的求解
- 参数的求解过程详细解释
- 理想抛物面求解模型的建立
- 基于几何模型的相关调节数据的计算
- (三)基于光反射定律和几何模型的馈源舱接收比的计算
- 主索节点调节与三角反射面板转动的关系
- 三角反射板法向量计算
- 馈源舱接收比的计算
- 可以改进的部分
- 参考文献
题目简述
2021年高教社杯全国大学生数学建模竞赛A题是关于“FAST”(500米口径球面射电望远镜)主动反射面形状调节的问题。
赛题要解决的问题是:在反射面板调节约束下,确定一个理想抛物面,然后通过调节促动器的径向伸缩量,将反射面调节为工作抛物面,使得该工作抛物面尽量贴近理想抛物面,以获得天体电磁波经反射面反射后的最佳接收效果。
具体问题
-
确定理想抛物面:当待观测天体位于基准球面正上方( α = 0 ° \alpha = 0 \degree α=0° , β = 90 ° \beta= 90 \degree β=90°)时,结合反射面板调节因素,确定理想抛物面的形状。
-
调节反射面:当待观测天体位于特定位置( α = 36.795 ° \alpha = 36.795 \degree α=36.795° , β = 78.169 ° \beta= 78.169 \degree β=78.169°)时,确定理想抛物面,并建立反射面板调节模型,通过调节促动器的伸缩量,使反射面尽量贴近理想抛物面。将理想抛物面的顶点坐标,以及调节后反射面 300 米口径内的主索节点编号、位置坐标、各促动器的伸缩量等结果按照规定的格式(见附件 4)保存在“result.xlsx”文件中。
-
计算接收比:基于第2问的反射面调节方案,计算调节后馈源舱的接收比(馈源舱有效区域接收到的反射信号与300米口径内反射面的反射信号之比),并与基准反射球面的接收比进行比较。
题目附录
- 主动反射面共有主索节点 2226 个,节点间连接主索 6525 根,不考虑周边支承结构连接的部分反射面板,共有反射面板 4300 块。基准球面的球心在坐标原点,附件 1 给出了所有主索节点的坐标和编号,附件 2 给出了促动器下端点(地锚点)坐标、基准态时上端点(顶端)的坐标,以及促动器对应的主索节点编号,附件 3 给出了 4300 块反射面板对应的主索节点编号。
- 基准态下,所有主索节点均位于基准球面上。
- 每一块反射面板均为基准球面的一部分。反射面板上开有许多直径小于 5 毫米的小圆孔,用于透漏雨水。由于小孔的直径小于所观察的天体电磁波的波长,不影响对天体电磁波的反射,所以可以认为面板是无孔的。
- 电磁波信号及反射信号均视为直线传播。
- 主索节点调节后,相邻节点之间的距离可能会发生微小变化,变化幅度不超过 0.07%。
- 将主索节点坐标作为对应的反射面板顶点坐标。
- 通过促动器顶端的伸缩,可控制主索节点的移动变位,但连接主索节点与促动器顶端的下拉索的长度保持不变。促动器伸缩沿基准球面径向趋向球心方向为正向。假设基准状态下,促动器顶端径向伸缩量为 0,其径向伸缩范围为-0.6~+0.6 米。
- 天体 S 的方位可用方位角𝛼和仰角𝛽来表示(见图 5)。
问题分析
理想抛物面的确定
理想抛物面的方程由顶点坐标和焦点坐标决定。观测天体位于基准球面正上方( α = 0 ° \alpha = 0 \degree α=0° , β = 90 ° \beta= 90 \degree β=90°),其距离FAST很远,故可以认为天体的电磁波均平行主轴射入天眼,可将馈源舱设定为焦点。当待观测天体位于基准球面正上方时,抛物面顶点位于Z轴上,通过主索节点坐标可确定顶点坐标,进而得到抛物线方程,再经旋转得到理想抛物面方程。
反射面板调节模型的建立
问题 2 分为两部分,第一部分要求在天体 S 位于位于 α = 36.795 ° \alpha = 36.795 \degree α=36.795° , β = 78.169 ° \beta= 78.169 \degree β=78.169° 时,建立理想抛物面关系式,并求出顶点坐标;第二部分要求在第一部分的基础上,计算调节后反射面 300 米口径内的主索节点被调节后的位置坐标、各促动器的伸缩量等数据
第一部分
目标天体的位置发生变化,馈源仓的中心在一定范围内随之移动,始终保持“天体-球心-馈源舱”三点一线。理想抛物面的焦点坐标和顶点坐标发生变化,但是焦距始终保持不变,故抛物面整体形状保持不变。因此,对于不同观测角度的理想抛物面,可利用坐标系变换和布尔莎模型求解,类同问题一中的理想抛物面,即可得出所求理想抛物面关系式。
第二部分
通过促动器顶端沿基准球面径向的伸缩,可控制主索节点的移动变位,且促动器顶端径向伸缩范围为-0.6~+0.6 米。在确定理想抛物面后,为调节反射面尽量贴近理想抛物面,则需要计算基准态下主索节点距理想抛物面的径向距离,来确定促动器的伸缩量和主索节点的调节坐标。
需要注意的是,假如存在节点距离理想抛物面的径向距离超过0.6米,则此时促动器的伸缩量不足以满足要求,则需要调整理想抛物面。具体的调整方式为——不断调整顶点的位置,重新计算理想抛物面和径向距离,直至满足要求;或者,采用最优化的方法,遍历所有的顶点位置(顶点所对应促动器伸缩-0.6~+0.6 米),寻找所有促动器移动量最少的情况。
馈源舱接收比计算
通过微元对每个三角反射板进行调整,得到调整后的三角反射面板主索节点坐标,从而确定调整后三角板的位置。接着通过几何模型求取每个三角发射板的法向量,利用光的反射定理和几何公式计算待测天体电磁波经过三角面板反射之后的信号,通过投影统计反射信号经过馈源舱有效区域的量,进而计算调节后馈源舱的接收比,并与基准面的接收比进行比较。
根据反射面调节后主索节点的坐标,确定三角反射面板的位置和偏转角度,利用光的反射定律计算反射信号方向,通过平行投影法和微元法计算馈源舱接收比。
模型假设
- 计算过程中所有数据保留五位小数,假设导致的误差在合理范围内,对结果影响可忽略不计。
- 收集的数据真实有效。
符号说明
符号 | 含义 |
---|---|
S | 待观测天体 |
C | 基准面球心 |
P | 馈源舱 |
R | 基准球面半径 |
A | 抛物面顶点 |
B | 基准球面上随机点 |
E | 椭圆方程 |
α \alpha α | 方位角 |
β \beta β | 仰角 |
θ \theta θ | SB 连线与Z 轴夹角 |
p | 焦距 |
a | 半长轴长 |
b | 半短轴长 |
J | 主索点 |
Q | 与反射信号垂直的平面 |
O | 基准球面圆心 |
模型建立与求解
(一)基于光学性质的理想抛物面确定模型
- 抛物线方程的确定
- 顶点坐标:当待观测天体S位于基准球面正上方时,抛物面顶点位于 Z Z Z 轴上,设为 ( 0 , 0 , z ) (0,0,z) (0,0,z)。通过Excel计算附件1中节点 Z Z Z 坐标的平均值,得到顶点A的坐标为 ( 0 , 0 , − 300.4 ) (0,0,-300.4) (0,0,−300.4)。
- 焦点坐标:根据旋转抛物面的光学性质,将馈源舱P设定为焦点,半焦距 1 2 p = 0.466 R 即 p = 0.932 R = 279.6 \frac{1}{2}p=0.466R\quad \text{即} \quad p = 0.932R = 279.6 21p=0.466R即p=0.932R=279.6
则焦点的 Z Z Z 坐标为 − 300.4 + 0.466 R = − 160.6 -300.4+0.466R=-160.6 −300.4+0.466R=−160.6,焦点坐标为P ( 0 , 0 , − 160.6 ) (0, 0, -160.6) (0,0,−160.6)。 - 抛物线方程:将顶点与焦点坐标代入抛物线方程 x 2 = 2 p y x^{2}=2py x2=2py,可得理想抛物线方程为 y = 0.0017857 x 2 − 300.4 y = 0.0017857x^{2}-300.4 y=0.0017857x2−300.4
- 抛物面方程确定:理想抛物面由抛物线绕Z轴旋转180°得到,令 y = z y = z y=z, x 2 = x 2 + y 2 x^{2}=x^{2}+y^{2} x2=x2+y2,得到理想抛物面方程 z = 0.0017857 ( x 2 + y 2 ) − 300.4 z = 0.0017857(x^{2}+y^{2})-300.4 z=0.0017857(x2+y2)−300.4
- 理想抛物面的检验:考虑反射面板调节约束,促动器顶端径向伸缩范围为-0.6 ~ +0.6米。在基准球面上随机取点 B B B 与基准面圆心 O O O 进行连线,设定与 Z Z Z 轴的夹角为 θ \theta θ,工作态时反射面口径与基准球面半径构成等边三角形, θ \theta θ范围为 15 0 ∘ ∼ 18 0 ∘ 150^{\circ}\sim 180^{\circ} 150∘∼180∘。
改变 θ \theta θ 的范围,计算基准面和理想抛物面之间的偏移量,发现当 θ \theta θ在 15 5 ∘ ∼ 16 5 ∘ 155^{\circ} \sim 165^{\circ} 155∘∼165∘时,偏移量超过0.6米,不满足条件。
- 理想抛物面的调整:将理想抛物面顶点定义在 [ − 299.8 , 301.0 ] [-299.8,301.0] [−299.8,301.0]之间移动,通过MATLAB计算不同顶点位置的最大、最小偏移量,得到顶点的z轴坐标位于 [ − 300.9 , − 300.5 ] [-300.9,-300.5] [−300.9,−300.5]之间时满足条件。
选取平均中心点 ( 0 , 0 , − 300.7 ) (0,0,-300.7) (0,0,−300.7)作为顶点,得到优化后的抛物面方程 z = 0.00178207 ( x 2 + y 2 ) − 300.7 z = 0.00178207(x^{2}+y^{2})-300.7 z=0.00178207(x2+y2)−300.7
(二)基于布尔莎模型的反射面板调节
在天体 S 位于位于 α = 36.795 ° \alpha = 36.795 \degree α=36.795° , β = 78.169 ° \beta= 78.169 \degree β=78.169° 时,建立理想抛物面关系式,并求出顶点坐标;并计算调节后反射面 300 米口径内的主索节点被调节后的位置坐标、各促动器的伸缩量等数据
布尔莎七参数
布尔莎七参数(Bursa-Wolf 参数)是一种用于坐标转换的数学模型,主要用于将一个坐标系统中的点转换到另一个坐标系统中。
布尔莎七参数的定义
布尔莎七参数包括三个平移参数、三个旋转参数和一个缩放参数,具体如下:
-
三个平移参数(Translation Parameters)
用于调整两个坐标系的原点位置差异- T x T_x Tx:沿X轴的平移量
- T y T_y Ty:沿Y轴的平移量
- T z T_z Tz:沿Z轴的平移量
-
三个旋转参数(Rotation Parameters)
用于调整两个坐标系的轴向差异。这些旋转角度通常以弧度表示,且在实际应用中通常非常小- ε x \varepsilon_x εx:绕X轴的旋转角度
- ε y \varepsilon_y εy:绕Y轴的旋转角度
- ε z \varepsilon_z εz:绕Z轴的旋转角度
-
一个缩放参数(Scale Factor)
用于调整两个坐标系的尺度差异。缩放因子通常是一个接近于1的值(如1 + s),其中 s s s 表示尺度变化量- s s s:缩放因子
坐标转换公式
假设点 P P P 在源坐标系中的坐标为 ( X , Y , Z ) (X, Y, Z) (X,Y,Z),在目标坐标系中的坐标为 ( X ′ , Y ′ , Z ′ ) (X', Y', Z') (X′,Y′,Z′),则布尔莎七参数的转换公式为:
[ X ′ Y ′ Z ′ ] = R [ X Y Z ] + [ T x T y T z ] + [ s X s Y s Z ] \begin{bmatrix} X' \\ Y' \\ Z' \end{bmatrix}= R \begin{bmatrix} X \\ Y \\ Z \end{bmatrix} + \begin{bmatrix} T_x \\ T_y \\ T_z \end{bmatrix} + \begin{bmatrix} sX \\ sY \\ sZ \end{bmatrix} X′Y′Z′ =R XYZ + TxTyTz + sXsYsZ
其中:
- 第一个矩阵是旋转矩阵,用于调整旋转差异。
- 第二个向量是平移向量,用于调整原点位置差异。
- 第三个向量是缩放向量,用于调整尺度差异。
R R R是欧拉角旋转矩阵,具体可参考https://blog.csdn.net/m0_49842530/article/details/144094059
假设旋转顺序为 绕Z轴旋转 ε z \varepsilon_z εz → 绕Y轴旋转 ε y \varepsilon_y εy → 绕X轴旋转 ε x \varepsilon_x εx,则旋转矩阵可以表示为:
R = R x ( ε x ) ⋅ R y ( ε y ) ⋅ R z ( ε z ) R = R_x(\varepsilon_x) \cdot R_y(\varepsilon_y) \cdot R_z(\varepsilon_z) R=Rx(εx)⋅Ry(εy)⋅Rz(εz)
其中:
R x ( ε x ) = [ 1 0 0 0 cos ( ε x ) − sin ( ε x ) 0 sin ( ε x ) cos ( ε x ) ] R_x(\varepsilon_x) = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos(\varepsilon_x) & -\sin(\varepsilon_x) \\ 0 & \sin(\varepsilon_x) & \cos(\varepsilon_x) \end{bmatrix} Rx(εx)= 1000cos(εx)sin(εx)0−sin(εx)cos(εx)
R y ( ε y ) = [ cos ( ε y ) 0 sin ( ε y ) 0 1 0 − sin ( ε y ) 0 cos ( ε y ) ] R_y(\varepsilon_y) = \begin{bmatrix} \cos(\varepsilon_y) & 0 & \sin(\varepsilon_y) \\ 0 & 1 & 0 \\ -\sin(\varepsilon_y) & 0 & \cos(\varepsilon_y) \end{bmatrix} Ry(εy)= cos(εy)0−sin(εy)010sin(εy)0cos(εy)
R z ( ε z ) = [ cos ( ε z ) − sin ( ε z ) 0 sin ( ε z ) cos ( ε z ) 0 0 0 1 ] R_z(\varepsilon_z) = \begin{bmatrix} \cos(\varepsilon_z) & -\sin(\varepsilon_z) & 0 \\ \sin(\varepsilon_z) & \cos(\varepsilon_z) & 0 \\ 0 & 0 & 1 \end{bmatrix} Rz(εz)= cos(εz)sin(εz)0−sin(εz)cos(εz)0001
参数的求解
布尔莎七参数通常通过已知的控制点(即在两个坐标系中都有坐标的点)来求解。具体步骤如下:
- 选择多个已知的控制点(至少需要三个不共线的点)。
- 建立坐标转换方程组。
- 使用最小二乘法或其他优化方法求解七个参数。
参数的求解过程详细解释
-
选择多个已知的控制点
控制点是那些在两个坐标系中都有已知坐标的点。这些点作为转换的参考,用于建立坐标转换方程组。通常,选择的控制点数量越多,求解的精度越高,但至少需要三个不共线的点来求解布尔莎七参数。 -
建立坐标转换方程组
对于每个控制点 P i P_i Pi,假设其在源坐标系中的坐标为 ( X i , Y i , Z i ) (X_i, Y_i, Z_i) (Xi,Yi,Zi),在目标坐标系中的坐标为 ( X i ′ , Y i ′ , Z i ′ ) (X_i', Y_i', Z_i') (Xi′,Yi′,Zi′)。根据布尔莎七参数的坐标转换公式,可以写出以下方程组:
[ X i ′ Y i ′ Z i ′ ] = R [ X i Y i Z i ] + [ T x T y T z ] + [ s X i s Y i s Z i ] \begin{bmatrix} X_i' \\ Y_i' \\ Z_i' \end{bmatrix}= R \begin{bmatrix} X_i \\ Y_i \\ Z_i \end{bmatrix} + \begin{bmatrix} T_x \\ T_y \\ T_z \end{bmatrix} + \begin{bmatrix} sX_i \\ sY_i \\ sZ_i \end{bmatrix} Xi′Yi′Zi′ =R XiYiZi + TxTyTz + sXisYisZi
对于 n n n 个控制点,可以得到 3 n 3n 3n 个方程。这些方程可以写成矩阵形式:
A x = b \mathbf{A} \mathbf{x} = \mathbf{b} Ax=b
其中- A \mathbf{A} A 是系数矩阵
- x \mathbf{x} x 是参数向量(包含 T x , T y , T z , ε x , ε y , ε z , s T_x, T_y, T_z, \varepsilon_x, \varepsilon_y, \varepsilon_z, s Tx,Ty,Tz,εx,εy,εz,s)
- b \mathbf{b} b是观测向量(包含所有控制点的目标坐标与源坐标的差值)
-
使用最小二乘法或其他优化方法求解七个参数
由于方程组通常是超定的(即方程数量多于未知数数量),直接求解可能没有精确解。因此,通常使用最小二乘法来求解参数向量 x \mathbf{x} x,使得残差平方和最小。最小二乘解可以通过以下公式求得:
x = ( A T A ) − 1 A T b \mathbf{x} = (\mathbf{A}^T \mathbf{A})^{-1} \mathbf{A}^T \mathbf{b} x=(ATA)−1ATb
其中, A T \mathbf{A}^T AT 是 A \mathbf{A} A 的转置矩阵, ( A T A ) − 1 (\mathbf{A}^T \mathbf{A})^{-1} (ATA)−1 是 A T A \mathbf{A}^T \mathbf{A} ATA 的逆矩阵。 -
注意事项
- 旋转矩阵的非线性
- 旋转矩阵 R R R 是非线性的,因为它依赖于旋转角度 ε x , ε y , ε z \varepsilon_x, \varepsilon_y, \varepsilon_z εx,εy,εz的三角函数。因此,直接使用最小二乘法可能需要迭代求解。
- 通常需要将非线性问题转化为线性问题,或者使用非线性最小二乘法(如 Levenberg-Marquardt 算法)来求解。
- 初始值的重要性
对于大角度旋转,初始值的选择对求解结果有很大影响。如果初始值不合适,可能导致求解失败或收敛到错误的解。
- 旋转矩阵的非线性
理想抛物面求解模型的建立
-
理想抛物面求解模型的分析:由于目标天体平行电磁波入射时焦距不变,理想抛物面形状不变。在坐标系合理变换下,利用问题一所得方程,通过布尔莎模型将新坐标系中的方程转换为原坐标系中的方程。
-
基于坐标转换的新坐标系的建立:将圆心C设为原点,通过偏转原坐标系使抛物面对称轴落在(z)轴上,且保证抛物面顶点距圆心距离不变,这样问题一所得方程中的常数项无需改变。
考虑旋转顺序为 绕Z轴旋转 ε z \varepsilon_z εz → 绕Y轴旋转 ε y \varepsilon_y εy → 绕X轴旋转 ε x \varepsilon_x εx,则相对应的布尔莎模型的参数为- T x = 0 T_x = 0 Tx=0
- T y = 0 T_y = 0 Ty=0
- T z = 0 T_z = 0 Tz=0
- ε x = 0 \varepsilon_x = 0 εx=0
- ε y = β − 9 0 ∘ \varepsilon_y = \beta - 90^{\circ} εy=β−90∘
- ε z = − α \varepsilon_z = -\alpha εz=−α
- s = 0 s = 0 s=0
-
理想抛物面坐标的求解:以圆心C为原点建立坐标系,根据 α \alpha α、 β \beta β确定问题二中抛物面顶点坐标。通过MATLAB计算角度关系和转换矩阵,得到 α = 36.75 9 ∘ \alpha = 36.759^{\circ} α=36.759∘, β = 78.16 9 ∘ \beta = 78.169^{\circ} β=78.169∘时的理想抛物面表达式,并代入问题一方程得到在原坐标系中的方程。
假设问题一中的理想抛物线的坐标为 ( X , Y , Z ) (X,Y,Z) (X,Y,Z),问题二中的理想抛物线的坐标为 ( X ′ , Y ′ , Z ′ ) (X',Y',Z') (X′,Y′,Z′),我们尝试将问题二中的理想抛物线通过旋转与问题一中理想抛物线重合,我们根据旋转,获取 ( X , Y , Z ) (X,Y,Z) (X,Y,Z)与 ( X ′ , Y ′ , Z ′ ) (X',Y',Z') (X′,Y′,Z′)之间转换关系
[ X Y Z ] = [ 1 0 0 0 cos ( ε x ) − sin ( ε x ) 0 sin ( ε x ) cos ( ε x ) ] [ cos ( ε y ) 0 sin ( ε y ) 0 1 0 − sin ( ε y ) 0 cos ( ε y ) ] [ cos ( ε z ) − sin ( ε z ) 0 sin ( ε z ) cos ( ε z ) 0 0 0 1 ] [ X ′ Y ′ Z ′ ] + [ T x T y T z ] + [ s X ′ s Y ′ s Z ′ ] = [ 1 0 0 0 1 0 0 0 1 ] [ cos ( − 11.83 1 ∘ ) 0 sin ( − 11.83 1 ∘ ) 0 1 0 − sin ( − 11.83 1 ∘ ) 0 cos ( − 11.83 1 ∘ ) ] [ cos ( − 36.79 5 ∘ ) − sin ( − 36.79 5 ∘ ) 0 sin ( − 36.79 5 ∘ ) cos ( − 36.79 5 ∘ ) 0 0 0 1 ] [ X ′ Y ′ Z ′ ] = [ 0.784 0.587 − 0.198 − 0.599 0.800 0 0.158 0.119 0.980 ] [ X ′ Y ′ Z ′ ] \begin{bmatrix} X\\ Y\\ Z \end{bmatrix}= \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos(\varepsilon_x) & -\sin(\varepsilon_x) \\ 0 & \sin(\varepsilon_x) & \cos(\varepsilon_x) \end{bmatrix} \begin{bmatrix} \cos(\varepsilon_y) & 0 & \sin(\varepsilon_y) \\ 0 & 1 & 0 \\ -\sin(\varepsilon_y) & 0 & \cos(\varepsilon_y) \end{bmatrix} \begin{bmatrix} \cos(\varepsilon_z) & -\sin(\varepsilon_z) & 0 \\ \sin(\varepsilon_z) & \cos(\varepsilon_z) & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} X' \\ Y' \\ Z' \end{bmatrix}+ \begin{bmatrix} T_x \\ T_y \\ T_z \end{bmatrix}+ \begin{bmatrix} sX' \\ sY' \\ sZ' \end{bmatrix} \\ = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} \cos(-11.831^{\circ}) & 0 & \sin(-11.831^{\circ}) \\ 0 & 1 & 0 \\ -\sin(-11.831^{\circ}) & 0 & \cos(-11.831^{\circ}) \end{bmatrix} \begin{bmatrix} \cos(-36.795^{\circ}) & -\sin(-36.795^{\circ}) & 0 \\ \sin(-36.795^{\circ}) & \cos(-36.795^{\circ}) & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} X' \\ Y' \\ Z' \end{bmatrix} \\ = \begin{bmatrix} 0.784 & 0.587 & -0.198 \\ -0.599 & 0.800 & 0 \\ 0.158 & 0.119 & 0.980 \end{bmatrix} \begin{bmatrix} X' \\ Y' \\ Z' \end{bmatrix} XYZ = 1000cos(εx)sin(εx)0−sin(εx)cos(εx) cos(εy)0−sin(εy)010sin(εy)0cos(εy) cos(εz)sin(εz)0−sin(εz)cos(εz)0001 X′Y′Z′ + TxTyTz + sX′sY′sZ′ = 100010001 cos(−11.831∘)0−sin(−11.831∘)010sin(−11.831∘)0cos(−11.831∘) cos(−36.795∘)sin(−36.795∘)0−sin(−36.795∘)cos(−36.795∘)0001 X′Y′Z′ = 0.784−0.5990.1580.5870.8000.119−0.19800.980 X′Y′Z′
将所得 ( X , Y , Z ) (X, Y, Z) (X,Y,Z) 带入问题一所得方程: z = 0.00178207 ( x 2 + y 2 ) − 300.7 z = 0.00178207(x^2 + y^2) - 300.7 z=0.00178207(x2+y2)−300.7,得到 α = 36.75 9 ∘ \alpha = 36.759^\circ α=36.759∘, β = 78.16 9 ∘ \beta = 78.169^\circ β=78.169∘ 时的理想抛物面表达式:
z = 3.823 x + 2.859 y − 6674.7 × − 1.916 × 1 0 − 7 x 2 + 5.122 × 1 0 − 7 x y − 3.424 × 1 0 − 7 y 2 + 0.0008753 y + 1.048 + 6352.9 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} z=3.823x+2.859y−6674.7×−1.916×10−7x2+5.122×10−7xy−3.424×10−7y2+0.0008753y+1.048+6352.9
-
理想抛物面点坐标对应自变量的范围求解:通过对理想抛物面在XOY坐标系投影椭圆方程的求解,确定 x x x、 y y y的范围。
我们已知椭圆标准方程为: x 2 a 2 + y 2 b 2 = 1 \frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}} = 1 a2x2+b2y2=1 其中 2 a 2a 2a为长轴, 2 b 2b 2b为短轴。我们设置长短轴分别为 a = 300.7 2 c o s ( β ) = 147.15605 a = \frac{300.7}{2}cos(\beta)=147.15605 a=2300.7cos(β)=147.15605 和 b = 300.7 2 = 150.35 b = \frac{300.7}{2} = 150.35 b=2300.7=150.35将椭圆E的长轴与短轴代入可得其方程: x 2 147.1560 5 2 + y 2 150.3 5 2 = 1 \frac{x^{2}}{147.15605^{2}}+\frac{y^{2}}{150.35^{2}} = 1 147.156052x2+150.352y2=1
椭圆E绕原点逆时针旋转 α \alpha α得到椭圆 E 1 E_1 E1。根据椭圆旋转坐标变换公式,椭圆绕原点逆时针旋转 α \alpha α之后的坐标关系为
x = y ′ sin α + x ′ cos α x = y'\sin\alpha + x'\cos\alpha x=y′sinα+x′cosα
y = y ′ cos α − x ′ sin α y = y'\cos\alpha - x'\sin\alpha y=y′cosα−x′sinα
然后将 x x x、 y y y代入原方程可得: ( y ′ sin α + x ′ cos α ) 2 a 2 + ( y ′ cos α − x ′ sin α ) 2 b 2 = 1 \frac{(y'\sin\alpha + x'\cos\alpha)^{2}}{a^{2}}+\frac{(y'\cos\alpha - x'\sin\alpha)^{2}}{b^{2}} = 1 a2(y′sinα+x′cosα)2+b2(y′cosα−x′sinα)2=1
将 α = 36.79 5 ∘ \alpha = 36.795^{\circ} α=36.795∘代入得椭圆 E 1 E_1 E1方程为:
( 0.5989527 y + 0.80078364 x ) 2 147.1560 5 2 + ( 0.8007864 y − 0.5989537 x ) 2 150.3 5 2 = 1 \frac{(0.5989527y + 0.80078364x)^{2}}{147.15605^{2}}+\frac{(0.8007864y - 0.5989537x)^{2}}{150.35^{2}} = 1 147.156052(0.5989527y+0.80078364x)2+150.352(0.8007864y−0.5989537x)2=1
椭圆 E 1 E_1 E1沿与 y y y轴负方向夹角为 α \alpha α平移53.43288得到 E 2 E_2 E2,根据平面直角坐标中图像平移关系可得椭圆 E 2 E_2 E2方程:
z = ( ( 0.59 ( y + 42.788 ) + 0.80 ( x + 32.003 ) ) 2 ) 147.15 6 2 + ( ( 0.80 ( y + 42.788 ) − 0.59 ( x + 32.004 ) ) 2 ) 150.3 5 2 − 1 z = \frac{((0.59(y + 42.788)+ 0.80(x + 32.003))^{2})}{147.156^{2}}+\frac{((0.80(y + 42.788)- 0.59(x + 32.004))^{2})}{150.35^{2}}-1 z=147.1562((0.59(y+42.788)+0.80(x+32.003))2)+150.352((0.80(y+42.788)−0.59(x+32.004))2)−1
x x x、 y y y的范围在椭圆曲线所围区域。即理想抛物面 E 2 E2 E2方程中, x x x, y y y满足范围为:
( ( 0.59 ( y + 42.788 ) + 0.80 ( x + 32.003 ) ) 2 ) 147.15 6 2 + ( ( 0.80 ( y + 42.788 ) − 0.59 ( x + 32.004 ) ) 2 ) 150.3 5 2 ≤ 1 \frac{((0.59(y + 42.788)+ 0.80(x + 32.003))^{2})}{147.156^{2}}+\frac{((0.80(y + 42.788)- 0.59(x + 32.004))^{2})}{150.35^{2}}\leq1 147.1562((0.59(y+42.788)+0.80(x+32.003))2)+150.352((0.80(y+42.788)−0.59(x+32.004))2)≤1
-
检验所得关系式准确程度:将所得理想抛物面顶点坐标代入表达式,验证在不同坐标系中的坐标误差小于 0.0001 % 0.0001\% 0.0001%,证明关系式准确。
基于几何模型的相关调节数据的计算
-
利用几何模型计算主索节点调节后的坐标:通过联立直线方程与抛物面方程,计算主索节点调节后的坐标。
为理想抛物面与以 C C C 为顶点的圆面相交切面图,其中, P P P 为理想抛物面焦点, D D D 为理想抛物面顶点, G G G 为原始在圆面上未经调节的主索节点, F F F 即为主索节点调节后的理想坐标。
由于促动器的伸缩量体现在节点和球心 O O O的连线方向上,我们设置直线表达式直线CG的方程
x − x 0 k 1 = y − y 0 k 1 = z − z 0 k 1 \frac{x - x_0}{k_1}=\frac{y - y_0}{k_1}=\frac{z - z_0}{k_1} k1x−x0=k1y−y0=k1z−z0
联立直线CG与抛物线求F点坐标(主索节点在调节后的坐标)的联立方程组:
{ z = 3.823 x + 2.859 y − 6674.7 × − 1.916 × 1 0 − 7 x 2 + 5.122 × 1 0 − 7 x y − 3.424 × 1 0 − 7 y 2 + 0.0008753 y + 1.048 + 6352.9 x − x 0 k 1 = y − y 0 k 1 = z − z 0 k 1 \begin{cases} z = 3.823x + 2.859y - 6674.7\times\sqrt{-1.916\times10^{-7}x^{2}+5.122\times10^{-7}xy - 3.424\times10^{-7}y^{2}+0.0008753y + 1.048}+6352.9\\ \frac{x - x_0}{k_1}=\frac{y - y_0}{k_1}=\frac{z - z_0}{k_1} \end{cases} {z=3.823x+2.859y−6674.7×−1.916×10−7x2+5.122×10−7xy−3.424×10−7y2+0.0008753y+1.048+6352.9k1x−x0=k1y−y0=k1z−z0
求解公式,获取最终调节量 -
判断调节后主索节点是否在为得到理想抛物面需调节区域内:根据理想抛物面点坐标对应自变量的范围进行初步判断,再结合主索节点到圆心距离与调节后点到圆心距离的差值是否满足 ∣ C G − C F ∣ ≤ 0.6 |CG-CF|\leq0.6 ∣CG−CF∣≤0.6米的条件,确定最终需调节的主索节点。
(三)基于光反射定律和几何模型的馈源舱接收比的计算
- 主索节点的调节与三角反射面板转动的关系:反射面调节通过下拉索与促动器配合,根据主索节点调节后的坐标确定三角反射面板的准确位置和调节后的角度。
- 三角反射板法向量计算:已知三角板上三个主索点坐标,利用向量叉乘计算平面法向量。
- 馈源舱接收比的计算
- 光的反射定律:反射光线与入射光线、法线满足“三线共面,两线分居,两角相等”,光具有可逆性。
- 空间中关于直线的对称问题:通过建立方程组求解空间中直线关于平面法线对称的直线方向向量。
- 利用微元法、平行投影法和几何模型计算馈源舱接收比:利用平行投影法将三维问题转化为二维问题,对馈源舱接收信号的有效区域投影,将反射信号投射到同一平面。采用微元法在特定区域撒点,判断点是否在三角形与圆的重合部分,计算重合面积,进而得到馈源舱接收比。最终得到调整后抛物面的接收比为(1.55805%),基准球面的接收比为(0.00042594%)。
主索节点调节与三角反射面板转动的关系
反射面调节主要依靠下拉索和促动器协作,调节主索点来形成工作抛物面。在这个过程中,主索网上的三角反射面板节点调节数量会有所不同,调节3个主索点的反射面板能最大程度贴合理想抛物面,调节2个或1个主索点的会有偏移,0个主索点调节的面板则保持不动。
通过之前的计算,我们已经得到了理想抛物面的方程以及三角反射面板主索点调节后的坐标。基于这些坐标,我们就能进一步确定各个三角反射面板的准确位置和调节后的角度。这一步是后续计算的基础,它为我们确定反射信号的方向提供了关键信息。
三角反射板法向量计算
在三维空间直角坐标系中,已知三角板上三个主索点坐标 J 1 ( x 1 , y 1 , z 1 ) J_1(x_1, y_1, z_1) J1(x1,y1,z1)、 J 2 ( x 2 , y 2 , z 2 ) J_2(x_2, y_2, z_2) J2(x2,y2,z2)、 J 3 ( x 3 , y 3 , z 3 ) J_3(x_3, y_3, z_3) J3(x3,y3,z3),我们设平面法向量 n → = ( a , b , c ) \overrightarrow{n}=(a, b, c) n=(a,b,c)。根据向量叉乘的原理,法向量 n → = J 1 J 2 → × J 1 J 3 → \overrightarrow{n}=\overrightarrow{J_1J_2}×\overrightarrow{J_1J_3} n=J1J2×J1J3,通过行列式计算可得:
a = ( y 2 − y 1 ) ( z 3 − z 1 ) − ( y 3 − y 1 ) ( z 2 − z 1 ) a=(y_2 - y_1)(z_3 - z_1) - (y_3 - y_1)(z_2 - z_1) a=(y2−y1)(z3−z1)−(y3−y1)(z2−z1)
b = ( z 2 − z 1 ) ( x 3 − x 1 ) − ( z 3 − z 1 ) ( x 2 − x 1 ) b=(z_2 - z_1)(x_3 - x_1) - (z_3 - z_1)(x_2 - x_1) b=(z2−z1)(x3−x1)−(z3−z1)(x2−x1)
c = ( x 2 − x 1 ) ( y 3 − y 1 ) − ( x 3 − x 1 ) ( y 2 − y 1 ) c=(x_2 - x_1)(y_3 - y_1) - (x_3 - x_1)(y_2 - y_1) c=(x2−x1)(y3−y1)−(x3−x1)(y2−y1)
只要把三角反射板对应的主索节点坐标代入这些公式,就能求出三角反射板的法向量。这个法向量在后续计算反射信号方向时起着关键作用。
馈源舱接收比的计算
-
光的反射定律
光的反射定律是我们计算的重要依据,它的内容可以归纳为“三线共面,两线分居,两角相等”,也就是说反射光线与入射光线、法线在同一平面上,反射光线和入射光线分居在法线两侧,反射角等于入射角,并且光具有可逆性,光路上是相等的。 -
空间中关于直线的对称问题
在计算反射信号方向时,空间中两条直线关于某平面的法线对称。已知其中一条直线的向量(单位向量)为 ( m 0 , n 0 , p 0 ) (m_0, n_0, p_0) (m0,n0,p0),平面法向量为 ( a , b , c ) (a, b, c) (a,b,c),设对称直线的方向向量(单位向量)为 ( x , y , z ) (x, y, z) (x,y,z)。我们先计算出 ( m 0 , n 0 , p 0 ) (m_0, n_0, p_0) (m0,n0,p0)的单位向量 ( m , n , p ) (m, n, p) (m,n,p),然后通过解方程组:
{ m + x = t a n + y = t b p + z = t c x 2 + y 2 + z 2 = 1 \begin{cases} m + x = ta \\ n + y = tb \\ p + z = tc \\ x^2 + y^2 + z^2 = 1 \end{cases} ⎩ ⎨ ⎧m+x=tan+y=tbp+z=tcx2+y2+z2=1
就能求出对称直线的方向向量,也就是反射信号的方向向量。 -
利用平行投影法将三维问题转化为二维问题**
我们采用平行投影法模拟反射信号投射到空间中馈源舱的过程。馈源舱接收信号的有效区域是直径1米的中心圆盘,我们让投射线垂直于中心圆盘,这样就得到了中心圆盘的投影和投影平面。把入射电磁波经三角反射面板反射后得到的反射信号投射到这个平面上,此时反射信号投影与中心圆盘投影在同一平面上,投影三角形与投影圆重合区域的面积大小和圆面积之比,就是调节后的馈源舱的接收比。
-
利用微元法和几何模型求解投影三角形与投影圆重合区域面积
根据各三角反射面板的位置、反射信号的方向向量以及圆盘投影平面的法向量,我们建立平面直角坐标系,这样就能得到投影三角形顶点坐标和投影圆圆心坐标。由于圆与三角形重合的情况比较复杂,我们采用微元法来求解。
在 x = [ − 0.5 , 0.5 ] x=[-0.5, 0.5] x=[−0.5,0.5], y = [ − 0.5 , 0.5 ] y=[-0.5, 0.5] y=[−0.5,0.5]这个区域内平均投放10000个点,然后判断每个点是否位于三角形与圆的重合部分。
经过计算,
- 调整后的抛物面所有重合面积和为 1101.318742215 m 2 1101.318742215m^2 1101.318742215m2,接收比为 1.55805 % 1.55805\% 1.55805%;
- 调整前的基准球面所有重合面积和为 0.301082555556663 m 2 0.301082555556663m^2 0.301082555556663m2,接收比为 0.00042594 % 0.00042594\% 0.00042594%。
这表明调整后的抛物面能让馈源舱接收到更多的反射信号,效果更好。
可以改进的部分
- 可以将模型的坐标系由直角坐标系,转化为极坐标系;
- 起初为了简洁,将三角板设置成了平面,实际上应该是弧面(当年的评判标准中明确指出)
参考文献
[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.