空中飞行器无源定位数学模型

1、空中飞行器无源定位数学模型摘 要：本文根据题目要求，综合运用数据拟合、GM 1,1预测等数学模型，结合无源定位技术算出了飞行器的位置参数以及在可用卫星数较少时，飞行器在不同时刻 的位置，最后考虑到定位的精度与效率问题，综合分析了卫星的数量与空间几何分布等 因素的影响，给出了卫星优选方法。对于问题（1）,通过建立地心坐标系，将物理几何关系与题中所给数据联立成非线 性方程组，将问题转化成了用matlab求解方程组，在求解过程中将固定在飞行器上的两 个测向阵方向d1,d2定义为单位向量，考虑到测量精度，通过联立九颗卫星参数得出的 18个方程与两个单位向量的性质3个方程共21个方程求解得出飞行器的位置

2、参数p x,y,z 二-5202.27,6604.13,3129.29km ， 距地表高度 h = 2603.54 km ，d, =0707059, 0.70716, -0.00019， d2 =：-0.70715, 0.70706, 0.00087 。对于问题（2）,同样利用问题（1）联立方程组的方法，得出五个不同时刻飞行器 的位置参数，再利用最小二乘法将各时刻飞行器的位置在各个坐标轴上的参数值进行拟合 ， 得至U 拟合 方程； x = -4.19 106 - 3.62 103t ， y =4.60 106 - 2.77 10、，z=6.07 10 6 624.65 t单位：km。将不同时刻

3、的拟合值与计算所得值之间的空间距离作为原始非负时间序列，建立GM （1,1）模型，对距离进行灰色预测，得到其距离误差的 预测值，最后综合考虑得出，t =70s时飞行器的位置为以-3934.43, 4796.87, 6114.851单位:k m为圆心，以53.75km 为半径的球体内，将距离误差值53.75km作为其预测的可靠度。问题（3）考虑当定位卫星较多时，飞行器的最优定位方式。通过用第1问的数据进行仿真，比较不同卫星组合测得的飞行器位置参数，发现卫星数量越多，空间分布越 均匀，飞行器定位经度越高。而卫星数量的增加会降低定位的效率，综合考虑定位精度 和效率，得到选取5颗在空间均匀分布的卫星是

4、最优的定位方式。最后用最大误差限0.1 对最优卫星组合方式进行检验。将误差看成均匀正态分布，建立非线性最小二乘拟合模型，即可算出相对于误差限为0.1时的角误差E（x, y,z） =3.78。关键词：飞行器；无源定位；最小二乘法；GM 1,1预测1问题重述目标定位技术是导航与制导技术的重要基础。 在现有的导航与制导技术中，卫星定 位技术是精度最高的，也是较为理想的导航与制导技术。对于空中飞行器，在其飞行过程中很容易接收到太空卫星的信号。现在考虑通过测 量飞行器与地球同步卫星的方向角来实现空中飞行器的自定位。在球心坐标系下，空中 飞行器P的空间坐标记为x,y,z，不妨设它同时能接收到N颗同步卫星的

5、信号，其N颗 同步卫星Xi空间坐标分别记为Xi, y , Zi （i =1,2,N）。为了方便检测与同步卫星的方向 角，在空中飞行器上固定安装了两个相互垂直的测向阵列，它们的指向分别为 d! d1x,dy,dz和d2 d2x,d2y,d2z。地球同步卫星与空中飞行器P的位置关系示意图如图1所示，5 L分别表示空中飞行器P的测向阵列方向di,d2与地球同步卫星Xi i =1,2 ,N的夹角。现在请你们建立数学模型研究解决下面的问题:（1）通过测量空中飞行器测向阵列方向d1和d2与多颗地球同步卫星的夹角:i和，建立空中飞行器定位的数学模型； 对于附表1所给出的9颗同步卫星的数据，试确定 空中飞行器

6、P的位置参数。（2）在某些特殊情况下，空中飞行器能直接检测到的同步卫星数量较少，可以利用空中飞行器在匀速飞行过程中多次检测的结果来实现定位。针对这种情况，试建立空中飞行器定位的数学模型；对附表2中给出的3颗同步卫星的检测数据，确定空中飞行器 P在第70秒时的位置参数，并分析其可靠性。（3）当可用同步卫星数量较多时，为了提高定位精度和定位效率，需要对可用的 同步卫星进行一定的优选。试研究具体的优选策略，并通过仿真，分析在检测方向角误 差限为0.1时空中飞行器的定位方法和精度。2问题分析2.1问题（1）的分析该问题要求确定空中飞行器的位置参数。题目中给出了地球同步卫星的参数，包 括其所在的经度和其

7、与固定在飞行器上的垂直测向阵列方向d1x,d2x之间的夹角:，-i，于是，我们可以通过建立球心坐标系，把所要求的位置参数设成未知数，并与已知的数据 相结合，然后利用空间几何关系建立方程组求解：将飞行器上两个测向阵列方向d2x ,d2x设成单位向量，则与所求飞行器位置参数共9个未知数，利用向量的余弦定理，每一个卫星可分别与飞行器的两个垂直测向阵列方向建立两个方程， 利用单位向量本身上数学 定义建立3个方程，则任意3个卫星的定位数据就可以确定的方程组， 即可以求的所设 的未知数，即可得飞行器的位置参数。在题目给出的数据中共有 9组卫星，我们可以建 立21组方程，方程数大于未知数，联立此超定方程，利

8、用 matlab即可求得解值。2.2问题（2）的分析由问题（1），我们可以得知，由三个卫星就可以确定飞行器的位置，但精度存在问 题。问题（2）研究的就是在只有三颗卫星的情况下，通过多次检测结果对匀速飞行的 飞行器进行定位。题目中给了我们3颗卫星在5个时刻的飞星定位数据，通过这5组数据 我们可以算出飞行器在5个不同时刻的位置（存在一定的位置误差）。利用这5组不同时 刻的位置信息，再结合飞行器匀速飞行的条件，进行数据拟合，可以得到一条飞行器的 运动轨迹（此轨迹在短时间内有效），则可得飞行器在未来短时间内位置参数随时间的 变化关系，则可得卫星在t =70s时的位置参数。再利用此运动轨迹对前 5组数据

9、进行检 验，可以算得一组残差值，将此组值代入灰色预测模型，得到其残差预测值与拟合值相 加即为所要预测飞行器更为精确的位置参数。通过计算前五个残差值的的方差与飞行器 在相应阶段内飞行距离的比值，看其百分比，即可算得飞行器位置参数的可靠性。2.3问题（3）的分析问题（3）考虑的是当可供选择的卫星较多时，怎样优选卫星，才能提高确定定位 的精度和效率。对问题（1）数据进行分析，9个卫星进行定位时提供的位置参数较多， 综合求得的位置参数降低了个别偏差较大的结果的影响，所以我们认为如果不考虑定位效率的话，在一定范围内，卫星的数量较多为好，而且这样更能保持定位的稳定性。所 以在此题中，九颗卫星定位的飞行器位

10、置参数精度最高，可认为题目（1）求出的飞行器位置即为实际位置。但当卫星数量较多时，随着卫星数量的增加可能精度提高并不明 显，但定位效率会降低，所以在定位精度变化不大的情况下，应尽量选取较少的卫星。 于是选取不同卫星数的组合求得的飞行器位置参数，与实际值比较，可筛选出最优的定位卫星数量，然后分析各最优卫星数量组合的空间分布，比较定位误差的大小，可得出卫星空间分布对定位高度近代的影响。在考虑方向角误差限为0.1时，对测量角误差求和，算出最小的角误差，即为精度。3模型假设（1）假设地球为匀质标准球体，其质心与球心重合；（2）卫星环绕运动以地心为中心的圆；（3）飞行器在很短的时间内飞行曲面为平面；（4

11、）固定在飞行器上的两个垂直测向阵列为单位向量；（5）检测方向角误差正态分布；（6）题目给出的数据真实可靠（7）飞行器接收到的信号均为卫星直接发射，不考虑经地面或其它天体折射的电 磁波信号。4符号说明i表示第i i =1,2,3,4,56 7,8,9颗地球同步卫星；p代表飞行器P x, y,z表示飞行器的空间坐标；X表示地球同步卫星；Xi x,y,Zi表示第i颗地球同步卫星的空间坐标；d! ddydz表示第一个测向阵列的单位向量；d2 d2x,d2y,d2z表示第二个测向阵列的单位向量；:i表示空中飞行器的测向阵列d!与地球同步卫星Xi的夹角;表示空中飞行器的测向阵列d2与地球同步卫星Xi的夹角

12、; 表示第i颗卫星Xi的经度（正值为东经，负值为西经）；I表示卫星X到球心的距离；h为飞行器到球心的距离；5模型的建立与求解5.1模型1的建立与求解5.1.1球坐标系的建立和飞行器定位以地心为坐标原点，地球地心指向北极方向的轴线为z轴，地心指向本初子午线与 赤道交点方向的轴线X轴，垂直于Z轴和X轴的轴线为丫轴，建立如图1球心直角坐标 轴系2:P （飞行器）山Xi (卫星i)图1直角空间坐标轴系不妨假设d!指向飞行器飞行的前行方向，d2垂直于d!指向飞行器左方（从飞行器后方观察）。卫星到坐标原点的距离为I，卫星Xi的坐标可表示为X i （I cos K，l si n K，0）, 通过已知条件“取

13、地球半径为6367 Km，同步卫星高度为35800 Km ”，可以得知：I =6367 Km 35800 Km =42167 Km。设飞行器P的坐标为P（x,y,z），则由卫星指向飞行器的向量可表示为X i p =任-Icostyj - l sin飞，召）利用余弦定理及向量运算关系可联立得方程组（1）:X ZI(i =1,2,3,9);Xi p d1 cos( X i p, dj = cos、n1咲间cos( X i p, d2)二 cosXi p d2PXi d2(i =1,2,3,9);2 2 2d1xd2yd3z =1;d 2d 2d 2 =12x3y3zd1x*d2x d1y*d2y

14、d1z* d20将表1中9组卫星测得的数据代入方程组，利用 matlab求解得结果如下:N x, y, z =：-5202.27,6604.13,3129.29 单位：kmd1x= 0.7070590.70716-0.00019d2x= -0.707150.707062-0.00087根据题中dix，d2x方向的假设可知，飞行器此刻的航向即为dix的指向，此时飞行器距地表高度为:2603.54 km。5.2模型2的建立与求解5.2.1利用最小二乘法对数据进行拟合最小二乘法3进行数据拟合的基本思路：给定平面上的（X1，% ，（Xn，yn）,求f（X ）， 使在X1, X2，Xn处的函数值与实验数

15、据丫1,讨2，y的偏差的平方和为最小.禾I用同一时刻检测到地球同步卫星的相关数据，调用问题（1）所得函数方程组1 i =1,2,3，可以求出给出的t = 0，t =10, t = 20, t = 30, t = 40五个时刻飞行器的位置参 数如表1:单位：km表1:不同时刻飞行器的位置参数时刻X轴丫轴Z轴t=0-4185.01844606.52806070.9011t=10-4154.0318P 4621.58006052.6508t=20-4110.73954667.00396102.6247t=30-4088.1703:4687.81856126.4660t=40-4036.9748471

16、1.748960652.2607通过以上信息，以时间t为横坐标，分别以X轴的数据，丫轴的数据，Z轴的数据 为纵坐标，建立平面坐标轴，利用matlab得出五个不同时刻拟合值与实测值之间的图形 关系，见图2、图3、图4:6X的坐标值拟合方程：x - _4.19 10 3 3.62t单位：km-4.02x 10飞行器坐标拟合-1实测值-预测值1-4.04-4.06-4.0825356-4.14-4.16-4.18-4.2101520时间值-4.1标坐坐-4.123040图2飞行器坐标拟合图一Y的坐标值拟合方程：y =4.60103 - 2.77t单位：km6x 10飞行器坐标拟合时间图3飞行器坐标拟

17、合图83Z的坐标值拟合方程：z=6.07 10o.624Tt单位：kmx 10飞行器坐标拟合i-1实测值预测值356.136.126.116.16.096.086.076.066.0525图4飞行器坐标拟合图三20时间409#拟合值如表2所示:表2 :飞行器位置参数拟合值单位：km时间t ( s)X轴丫轴Z轴0-4187.376674603.599756071.0807110-4151.181814631.26780 16077.3273320-44414.986964658.935856083.5737330-4078.792104686.603896089.8202440-4042.597

18、254714.271946096.0667670-3934.219164797.070626114.85119由此可知，t =70s时，飞行器的位置参数为-3934.21916 ,4797.07062,6114.851195.2.2建立GM (1,1)模型对距离进行灰色预测建立GM (1,1)模型4的方法如下：设x(0)=X(1), x(0)(2)，x(0)( n)?为原始非负时间序列，x(t)为累加生成序列，即: (0)x (i) - 7 x (t), i =1,2,，nt AGM (1,1)的白化微分方程为：心 ax")dt式2中，a为待辨识参数，亦称发展系数；u为待辨识内生变量

19、，亦称为灰作用量,待辨识向量，a，按最小二乘法得a=(btB）BTy式中1 ( 1 ) ( 1 ) (x( )( 2 ) x ()(2(3 )1 ( 1 ) ( 1 ) -(x ( 2 ) x (3)20x ( 2x0( 3(x(12n -1) x (n)x0( n)10#于是可得到灰色预测的离散时间相应函数为:(t 1 x(0)(1)丄a丿u+ax（1）（t 1）为所得累加的预测值，将预测值还原即为:x (t 1x(1 )t( +"xt ( t，# 1，21,(3)(4 )#计算拟合位置与测得位置之间的距离作为原始时间序列，即xC )=(x(0 *1 ) x(0 尢2 ), xC

20、3 ), xf 丸4 ), xC *5 )=3.7641, 22.6628, 21.1204, 37.8462, 31.4505根据灰色预测模型得到t=0：70s时距离残差的预测值如表3所示: 表 3：距离残差预测值（单位： km ）时间t （ s）距离残差03.76411022.40232025.92053029.99134034.70135040.15106046.45667053.75243由此可知，t =70s时，飞行器的位置为以_3934.4256, 4796.8652, 6114.8512为圆心,以距离残差值预测值53.7524km为半径的球体内，位置可靠度即为53.7524 km

21、表4：灰色预测模型的检验:时间实际数据模拟数据残差相对误差(0)x (k)c(°)?(k)(0)c(0)呂(k) = x (k) (k)|x (k)|t =1022.662822.40230.26051.1492%t =2021.120420.64040.48002.2727%t =3037.846237.06070.78552.0755%t =4031.450431.77550.32491.0337%平均相对误差:1 51k (0.011492+0.022727+0.020755+0.010337)=1.6375%4 k丄45.3模型3的建立与求解5.3.1优选方法的确定考虑卫星

22、空间分布 对精度 的影响，对同等 数量的卫星，选取两 组分布差 别较 大的组合来计算定位精度。用问题（1）的方法联立方程组即可求 得飞行器的位置 参数，再比较与实际位置的距离，得出其精度。不同数量的卫星组合得到的飞行 器位置参数，见表5：表5 不同卫星组合 位置参数对照表卫星组合位置参数（单位:km）距离（km）(2, 4, 6, 8)(-5210.6254，6608.6498，3217.9274)83.8937(1,3,7, 8)(-5203.9396 ，6605.6841，3149.3620)15.1494(1,3,5, 7, 9)(-5198.7077 ，6604.1294，3122.9

23、055)11.9927(2, 4, 7, 8, 9)(-5202.8658 ，6605.9736，3166.5278)31.8727(1,3,4, 7, 8, 9)(-5198.2643 ，6605.1247，3124.3896)10.5407(2, 3,5, 7, 8, 9,)(-5201.3893 ，6605.7027，3157.3946)22.6420(1,2,3, 4, 6, 8,9)(-5200.2716 ，6605.2229，3139.0540)4.2707(1,2,3, 4, 5, 7,8)(-5202.5331，6605.5371，3147.9762)13.4573(1, 2,

24、 3,4,5, 6, 7, 8)(-5202.4130 ，6605.8972，3150.5806)16.0086(1,2,4, 5, 6, 7, 8, 9)(-5198.5059 ，6604.0057，3121.0125)13.900011由上表数据得 出相同数量卫星组合测得的距离误差平均值，见表 6:卫星数45678距离误差均值（km）49.521621.932716.59148.864014.9543表6 不同数量卫星组合距离误差均值表不同卫星数组合距离误差50454030252015105idata 3-144.555.566.577.58卫星数图5:不同卫星数组合得到的距离误差均值位单

25、差误离距从图我们可以 看出卫星组合1,2,3,468,9得到的距离值最小，表示由这组卫星所求得的飞 飞行器位置参 七颗卫星组合 到的飞行器位 位的成本是不 可认为，当可 量。行器位置与最优位置的距离最小 数精度最高，同时观察其他卫星 得到的另一个距离值13.4573km置参数精度较高，分析定位的效 断增加的，同时效率不断降低。 供选择的定位卫星数量较多时，所以可以认为此组合卫星测 组合所得到的距离值，可以看 亦较小，于是认为七颗卫星组 率，随着定位卫星数量的增加 结合定位的精度和效率综合考 选择5颗进行定位是最优的定得的 出由 合得，疋 虑， 位数在对定位卫星 位置进行优化，对相同数量不同经度

26、的两组卫星分析得知，当卫星在赤道上 空分布越均匀时，误差越小。见表7:12表7不同数量卫星组合经度差值表卫星数卫星组合相邻卫星经度差距离4颗卫星(2, 4, 6, 8)(36,11,25 )83.8937(1,3, 7, 8)(34,31,19 )P 15.14945颗卫星(1,3,5, 7, 9)(34,20,12,28 )11.9927(2, 4, 7, 8, 9)(36,17,19,9 )31.87276颗卫星(1,3,4, 7, 8, 9)(34,15,17,19,9 )10.5407(2, 3,5, 7, 8,9,)(21,20,12,19,9 )22.6427颗卫星(1,2, 3,

27、4, 6,8, 9)(13,21,15,11,25,9)4.2707(1,2, 3,4, 5,7, 8)(13,21,15,5,12,19)13.45738颗卫星(1,2,3, 4, 5, 6, 7, 8)(13,21,15,5,6,6,19,9)16.0086(1,2,4, 5, 6, 7, 8,9)(13,21,36,5,6,6,19,9)13.9说明：相邻卫星经度差(36,11,25 )中的36是指编号2、4两个卫星的经 度的差值，11是指编号4、6两个卫星的经度的差值，25是指编号& 8两个卫星的经度的差值，其余类推经过以上分析可知，卫星数量越多， 卫星在赤道上分布越均匀，所

28、测得的误 差值越小。经比 较，选取5颗夹角在1 5 25之间的卫星时，定位精 度和定位效率 最优。5.3.2检测误差角误差限为0.1度时飞行器的定位在多次测量中，设备的测量误差是正态分布的，而且均值为0。问题中给出的精度 0.10可以认为是测量误差的标准差。在这种理解下，用各自的误差限对测量误差进行无量纲化(也可以看成是一种加权法)处理是合理的，即求解的无约束优化问题更合理。由5.3.1的分析可知，当五颗卫星均匀分布时，定位精度和效率最优，故选取优化建模如下:M ini -6E (x, y, z)=i 413#其中COS片pXi a r ( iPXi d,PX i d2 PXi d272,3,

29、4, 5)-1,2, 3,4, 5)#由于目标函数是平方和的形式，因此这是一个非线性最小二乘拟合问题5 O建立优化模型，编写LIN GO程序，求得飞行器坐标为5198.51, 6604.01, 3121.01，这个解对应的 目标函数值大约为E x,y,z =3.78。优化建模这个误差数值较大，这是因为这里使用的 是相对于误差限的误差，分母二很小所以相对误差很大。其实，可以认为此时的目标 函数是5个标准正态分布的误差平方和，只要在 5以内都是合理的。6模型的评价及改进6.1模型的优点模型（1）建立球心坐标系使求解过程简单明了。模型（2）在（1）问题的分析基础上利用同一时刻给出的三颗卫星的参数算出

30、各时 刻对应的飞行器位置参数，然后利用最小二乘法拟合成曲线，推算出运动模型。并用GM 1,1模型预测拟合位置与实际位置距离，则其实际位置等于拟合值跟预测值的和， 将误差的变化考虑进模型，提高了定位的精度。模型3在问题的分析过程中，考虑题目所给的条件，将不同数量的卫星组合得到 的不同的飞行器位置，再与飞行器实际位置相比较，得出不同组合方式的测量精度。通 过多组的比较，选出最优组合，整个问题的求解过程中，大量运用了的图表，使运算过程与运算结果简单明了。6.2模型的不足模型（1）的求解过程中没有考虑到那些个别会引起测量偏差较大的卫星对结果的 影响，即没有对卫星进行优选。模型2仅仅利用三颗卫星提供的参

31、数得出飞行器在不同时刻的位置参数本身就存 在较大误差，解的过程中我们发现飞行器的位置参数在 Z轴上的数据残差过大。模型（3）直接将问题（1）中9颗卫星组合算出的飞行器位置作为实际位置有所欠 妥，同时没有充分考虑所有卫星组合，对结果产生一定的误差。6.3模型的改进方向模型（2）求解过程中，将的飞行器各轴上的位置参数分别进行数据拟合，最后得 到t =70s时飞行器的位置在一个球体范围内， 难以定位，可综合各轴位置参数进行拟合 得出飞行器的空间轨迹随时间得变化规律，可得预测到更为精确直观的飞行器位置参 数。模型3中，考虑到选取卫星组合的人为主观因素，可尽量多的选取卫星组合对飞行器进行定位，然后对比分

32、析出各定位精度大小与定位效率，而得出最优解。参考文献孙仲康 郭福成 冯道旺.单站无源定位跟踪技术国防大学出版社 2008.11第1版 徐天园.单站无源定位算法研究.杭州电子科技大学硕士学问论文.2009.12 Frank R.Giordano William P.Fox 等.机械工业出版社.2011.7第 1 版卓金武.Matlab在数学建模中的应用.北京航天航空大学出版社.2011.4第1版5百度文库.飞机的精确定位问题.2012.1115附表1:某时刻空中飞行器检测到地球同步卫星的相关数据卫星编号经度®i1E7621 .3369 .252E8935 .8254 .633E1106

33、0 .7329 .81自4E12579 .2411 .995E13085 .466 .956E13692 .915 .997E142100 .3211 .588E163125 .5836 .049E172135 .9346 .34注：取地球半径为6367 Km，同步卫星高度为35800 Km ； E表示东京。附表2：空中飞行器连续检测到地球同步卫星的相关数据检测时间（检测间隔10s）卫星编号经度*t =01E83 °38.45 &53.04 &2E136 °91.61 09.76 ©3E172132.3843.85t =101E89 °3

34、8.51 °52.99 °2E136 091.68 09.78 °3E172 °132.43 °43.89 °t =201E89 &38.56 &52.94 &2E136 °91.76 o9.79 °3E172 °132.47 °43.94 &t =301E89 °38.62 &52.89 &2E13691.839.813E172 °132.52 °43.99 °t =401E89 °38.67 &

35、#176;52.83 °2E136 °91.91 09.82 03E172 °132.57 °44.03 °注：取地球半径为6367 Km，同步卫星高度为35800 Km ； E表示东京。非线性方程求解程序：function F=f1(x0)a=x0(1);b=x0 (2);c=x0 (3);d=xO 4);e=x0(5);f=x0(6);x=x 0(7);y=xo(8)；z=x0(9);D=42167;e1=(a*(x-D*cos(76*pi/180)+b*(y-D*si n(76*pi/180)+c*z)/(xA2+yA2+zA2+DA2-

36、2*D *(x*cos(76*pi/180)+y*si n(76*pi/180)A0.5-cos(21.33*pi/180);e2=(a*(x-D*cos(89*pi/180)+b*(y-D*si n( 89*pi/180)+c*z)/(xA2+yA2+zA2+DA2-2*D *(x*cos(89*pi/180)+y*si n( 89*pi/180)A0.5-cos(35.82*pi/180);e3=(a*(x-D*cos(110*pi/180)+b*(y-D*si n(110*pi/180)+c*z)/(xA2+yA2+zA2+DA2-2 *D*(x*cos(110*pi/180)+y*si

37、 n(110*pi/180)A0.5-cos(60.73*pi/180);e4=(d*(x-D*cos(76*pi/180)+e*(y-D*s in (76*pi/180)+f*z)/(xA2+yA2+zA2+DA2-2*D *(x*cos(76*pi/180)+y*si n(76*pi/180)A0.5-cos(69.25*pi/180);e5=(d*(x-D*cos(89*pi/180)+e*(y-D*s in (89*pi/180)+f*z)/(xA2+yA2+zA2+DA2-2*D *(x*cos(89*pi/180)+y*si n( 89*pi/180)A0.5-cos(54.63*

38、pi/180);e6=(d*(x-D*cos(110*pi/180)+e*(y-D*si n(110*pi/180)+f*z)/(xA2+yA2+zA2+DA2-2 *D*(x*cos(110*pi/180)+y*si n(110*pi/180)A0.5-cos(29.81*pi/180);e7=(d*(x-D*cos(125*pi/180)+e*(y-D*si n(125*pi/180)+f*z)/(xA2+yA2+zA2+DA2-2*D*(x*cos(125*pi/180)+y*si n(125*pi/180)A0.5-cos(11.99*pi/180);e8=(a*(x-D*cos(12

39、5*pi/180)+b*(y-D*si n(125*pi/180)+c*z)/(xA2+yA2+zA2+DA2-2 *D*(x*cos(125*pi/180)+y*si n(125*pi/180)A0.5-cos(79.24*pi/180);e9=(a*(x-D*cos(130*pi/180)+b*(y-D*si n(130*pi/180)+c*z)/(xA2+yA2+zA2+DA2-2 *D*(x*cos(130*pi/180)+y*si n(130*pi/180)A0.5-cos(85.46*pi/180);e10=(d*(x-D*cos(130*pi/180)+e*(y-D*s in (

40、130*pi/180)+f*z)/(xA2+yA2+zA2+DA2-2*D*(x*cos(130*pi/180)+y*si n(130*pi/180)A0.5-cos(6.95*pi/180);e11=(a*(x-D*cos(136*pi/180)+b*(y-D*si n(136*pi/180)+c*z)/(xA2+yA2+zA2+DA2-2*D*(x*cos(136*pi/180)+y*si n(136*pi/180)A0.5-cos(92.91*pi/180);e12=(a*(x-D*cos(136*pi/180)+b*(y-D*si n(136*pi/180)+c*z)/(xA2+yA2

41、+zA2+DA2-2*D*(x*cos(136*pi/180)+y*si n(136*pi/180)A0.5-cos(5.99*pi/180);e13=(d*(x-D*cos(142*pi/180)+e*(y-D*si n(142*pi/180)+f*z)/(xA2+yA2+zA2+DA2-2*D*(x*cos(142*pi/180)+y*si n(142*pi/180)A0.5-cos(11.58*pi/180);e14=(a*(x-D*cos(142*pi/180)+b*(y-D*si n(142*pi/180)+c*z)/(xA2+yA2+zA2+DA2-2*D*(x*cos(142*p

42、i/180)+y*si n(142*pi/180)A0.5-cos(100.32*pi/180);e15=(a*(x-D*cos(163*pi/180)+b*(y-D*si n(163*pi/180)+c*z)/(xA2+yA2+zA2+DA2-2*D*(x*cos(163*pi/180)+y*si n(163*pi/180)A0.5-cos(125.58*pi/180);e16=(d*(x-D*cos(163*pi/180)+e*(y-D*s in (163*pi/180)+f*z)/(xA2+yA2+zA2+DA2-2*D*(x*cos(163*pi/180)+y*si n(163*pi/180)A0.5-cos(36.04*pi/180);e17=(a*(x-D*cos(172*pi/180)+b*(y-D*si n(1