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

1空中飞行器无源定位数学模型摘 要：本文根据题目要求，综合运用数据拟合、 预测等数学模型，结合1,GM无源定位技术 1算出了飞行器的位置参数以及在可用卫星数较少时，飞行器在不同时刻的位置，最后考虑到定位的精度与效率问题，综合分析了卫星的数量与空间几何分布等因素的影响，给出了卫星优选方法。对于问题（1）,通过建立地心坐标系，将物理几何关系与题中所给数据联立成非线性方程组，将问题转化成了用 求解方程组，在求解过程中将固定在飞行器上matlb的两个测向阵方向 定义为单位向量，考虑到测量精度，通过联立九颗卫星参数得12,d出的 18 个方程与两个单位向量的性质 3 个方程共 21 个方程求解得出飞行器的位置参数 ，距地表高度 ，,-50.7,64.,19.2pxyzk2603.54hkm， 。1.79,1d 0.715,.,.87d对于问题（2），同样利用问题（1）联立方程组的方法，得出五个不同时刻飞行器的位置参数，再利用最小二乘法将各时刻飞行器的位置在各个坐标轴上的参数值进行拟合，得到拟合方程； ， ，634.91.2tx63y4.01 2.0t。将不同时刻的拟合值与计算所得值之间的空间距6.07124.5ztkm单 位 :离作为原始非负时间序列，建立 模型，对距离进行灰色预测，得到其距离误差(,)GM的预测值，最后综合考虑得出， 时飞行器的位置为以70ts为圆心，以 为半径的球体内，将距离394.,76.8,14.5单 位 :k53.7km误差值 作为其预测的可靠度。5km问题（3）考虑当定位卫星较多时，飞行器的最优定位方式。通过用第 问的数1据进行仿真，比较不同卫星组合测得的飞行器位置参数，发现卫星数量越多，空间分布越均匀，飞行器定位经度越高。而卫星数量的增加会降低定位的效率，综合考虑定位精度和效率，得到选取 颗在空间均匀分布的卫星是最优的定位方式。最后用最大误5差限 对最优卫星组合方式进行检验。将误差看成均匀正态分布，建立非线性最小0.1二乘拟合模型，即可算出相对于误差限为 时的角误差 。0.1(,)3.78Exyz关键词：飞行器；无源定位；最小二乘法； 预测,GM21 问题重述目标定位技术是导航与制导技术的重要基础。在现有的导航与制导技术中，卫星定位技术是精度最高的，也是较为理想的导航与制导技术。对于空中飞行器，在其飞行过程中很容易接收到太空卫星的信号。现在考虑通过测量飞行器与地球同步卫星的方向角来实现空中飞行器的自定位。在球心坐标系下，空中飞行器 的空间坐标记为 ，不妨设它同时能接收到 颗同步卫星的信号，P,xyzN其 颗同步卫星 空间坐标分别记为 。为了方便检测与同步卫NiX,(1,2)ii星的方向角，在空中飞行器上固定安装了两个相互垂直的测向阵列，它们的指向分别为 和 。地球同步卫星与空中飞行器 的位置关系示意图11,xyzd22,xyzd P如图1所示， 分别表示空中飞行器P的测向阵列方向 与地球同步卫星,i 12,d的夹角。现在请你们建立数学模型研究解决下面的问题： ,2iXN（1）通过测量空中飞行器测向阵列方向 和 与多颗地球同步卫星的夹角 和12 i，建立空中飞行器定位的数学模型；对于附表1所给出的9颗同步卫星的数据，试确i定空中飞行器P的位置参数。 （2）在某些特殊情况下，空中飞行器能直接检测到的同步卫星数量较少，可以利用空中飞行器在匀速飞行过程中多次检测的结果来实现定位。针对这种情况，试建立空中飞行器定位的数学模型；对附表2中给出的3颗同步卫星的检测数据，确定空中飞行器P在第70秒时的位置参数，并分析其可靠性。 （3）当可用同步卫星数量较多时，为了提高定位精度和定位效率，需要对可用的同步卫星进行一定的优选。试研究具体的优选策略，并通过仿真，分析在检测方向角误差限为 时空中飞行器的定位方法和精度。 0.12 问题分析2.1 问题（1）的分析该问题要求确定空中飞行器的位置参数。题目中给出了地球同步卫星的参数，包括其所在的经度和其与固定在飞行器上的垂直测向阵列方向 之间的夹角 ，12,xd,i于是，我们可以通过建立球心坐标系，把所要求的位置参数设成未知数，并与已知的数据相结合，然后利用空间几何关系建立方程组求解：将飞行器上两个测向阵列方向设成单位向量，则与所求飞行器位置参数共 9 个未知数，利用向量的余弦定理，2,xd3每一个卫星可分别与飞行器的两个垂直测向阵列方向建立两个方程，利用单位向量本身上数学定义建立 3 个方程，则任意 3 个卫星的定位数据就可以确定的方程组，即可以求的所设的未知数，即可得飞行器的位置参数。在题目给出的数据中共有 9 组卫星，我们可以建立 组方程，方程数大于未知数，联立此超定方程，利用 即可求得21 matlb解值。2.2 问题（2）的分析由问题（1） ，我们可以得知，由三个卫星就可以确定飞行器的位置，但精度存在问题。问题（2）研究的就是在只有三颗卫星的情况下，通过多次检测结果对匀速飞行的飞行器进行定位。题目中给了我们 颗卫星在 个时刻的飞星定位数据，通过这 组355数据我们可以算出飞行器在 个不同时刻的位置（存在一定的位置误差） 。利用这 组5不同时刻的位置信息，再结合飞行器匀速飞行的条件，进行数据拟合，可以得到一条飞行器的运动轨迹（此轨迹在短时间内有效） ，则可得飞行器在未来短时间内位置参数随时间的变化关系，则可得卫星在 时的位置参数。再利用此运动轨迹对前 组70ts数据进行检验，可以算得一组残差值，将此组值代入灰色预测模型，得到其残差预测值与拟合值相加即为所要预测飞行器更为精确的位置参数。通过计算前五个残差值的的方差与飞行器在相应阶段内飞行距离的比值，看其百分比，即可算得飞行器位置参数的可靠性。2.3 问题（3）的分析问题（3）考虑的是当可供选择的卫星较多时，怎样优选卫星，才能提高确定定位的精度和效率。对问题（1）数据进行分析，9 个卫星进行定位时提供的位置参数较多，综合求得的位置参数降低了个别偏差较大的结果的影响，所以我们认为如果不考虑定位效率的话，在一定范围内，卫星的数量较多为好，而且这样更能保持定位的稳定性。所以在此题中，九颗卫星定位的飞行器位置参数精度最高，可认为题目（1）求出的飞行器位置即为实际位置。但当卫星数量较多时，随着卫星数量的增加可能精度提高并不明显，但定位效率会降低，所以在定位精度变化不大的情况下，应尽量选取较少的卫星。于是选取不同卫星数的组合求得的飞行器位置参数，与实际值比较，可筛选出最优的定位卫星数量，然后分析各最优卫星数量组合的空间分布，比较定位误差的大小，可得出卫星空间分布对定位高度近代的影响。在考虑方向角误差限为 时，对测0.1量角误差求和，算出最小的角误差，即为精度。3 模型假设（1）假设地球为匀质标准球体，其质心与球心重合；（2）卫星环绕运动以地心为中心的圆；（3）飞行器在很短的时间内飞行曲面为平面；（4）固定在飞行器上的两个垂直测向阵列为单位向量；（5）检测方向角误差正态分布；（6）题目给出的数据真实可靠（7）飞行器接收到的信号均为卫星直接发射，不考虑经地面或其它天体折射的电4磁波信号。4 符号说明表示第 颗地球同步卫星；i1,2345,6789i代表飞行器p表示飞行器的空间坐标；,Pxyz表示地球同步卫星；X表示第 颗地球同步卫星的空间坐标；,iii表示第一个测向阵列的单位向量；11xyzd表示第二个测向阵列的单位向量；22,z表示空中飞行器的测向阵列 与地球同步卫星 的夹角；i1diX表示空中飞行器的测向阵列 与地球同步卫星 的夹角；i2i表示第 颗卫星 的经度（正值为东经，负值为西经） ；iiix表示卫星 到球心的距离；lX为飞行器到球心的距离；h5 模型的建立与求解5.1 模型 1 的建立与求解5.1.1 球坐标系的建立和飞行器定位以地心为坐标原点，地球地心指向北极方向的轴线为 轴，地心指向本初子午线Z与赤道交点方向的轴线 轴，垂直于 轴和 轴的轴线为 轴，建立如图 1 球心直角XZXY坐标轴系 2：5图 1 直角空间坐标轴系不妨假设 指向飞行器飞行的前行方向， 垂直于 指向飞行器左方（从飞行器1d 2d1后方观察） 。卫星到坐标原点的距离为 ，卫星 的坐标可表示为 ,liXiX(cos,in0)ll通过已知条件“取地球半径为 ，同步卫星高度为 ”，可以得知：637Km3580Km。设飞行器 的坐标为 ，则由卫星指向飞637580421lKmP(,)xyz行器的向量可表示为 利用余弦定理及向量运算关系可联(cos,in,)iii iXpxlylz立得方程组(1)：112212321211cos(,)cos(,3,9);(,)(,);*0ii iixyzxyzXpdpdiPXidd 1将表 1 中 9 组卫星测得的数据代入方程组，利用 求解得结果如下：matlb,-0.7,64.3,29.Nzk单 位 ：61=0.759.0716.9xd2 2087根据题中 方向的假设可知，飞行器此刻的航向即为 的指向 ，此时飞行器距12,xd 1xd地表高度为： 。603.54km5.2 模型 2 的建立与求解5.2.1 利用最小二乘法对数据进行拟合最小二乘法 3进行数据拟合的基本思路：给定平面上的 ，求1,nxy，使在 处的函数值与实验数据 的偏差的平方和为最小. fx12,nx 12,ny利用同一时刻检测到地球同步卫星的相关数据，调用问题（1）所得函数方程组，可以求出给出的 五个时刻飞行器的位置参,3i0,304ttt数如表 1： km单 位 ：表1：不同时刻飞行器的位置参数时刻 X轴 Y轴 Z轴t=0 -4185.0184 4606.5280 6070.9011t=10 -4154.0318 4621.5800 6052.6508t=20 -4110.7395 4667.0039 6102.6247t=30 -4088.1703 4687.8185 6126.4660t=40 -4036.9748 4711.7489 60652.2607通过以上信息，以时间 为横坐标，分别以 轴的数据， 轴的数据， 轴的数据tXYZ为纵坐标，建立平面坐标轴，利用 得出五个不同时刻拟合值与实测值之间的图matlb形关系，见图 2、图 3、图 4：7的坐标值拟合方程：X34.190.62txkm单 位 ：的坐标值拟合方程：Y3y4.601 2.7tkm单 位 ：图 2 飞行器坐标拟合图一 0 5 10 15 20 25 30 35 40-4.2-4.18-4.16-4.14-4.12-4.1-4.08-4.06-4.04-4.02x 106 不不X不不不 不不不不不不不 不不不不不不图 3 飞行器坐标拟合图二 0 5 10 15 20 25 30 35 404.64.624.644.664.684.74.72x 106不不Y不不不不不不不不不不不不不不不不8的坐标值拟合方程：Z36.071.6247ztkm单 位 ：拟合值如表 2 所示：表 2：飞行器位置参数拟合值 km单 位 ：时间 （s ）tX 轴 Y 轴 Z 轴0 -4187.37667 4603.59975 6071.0807110 -4151.18181 4631.26780 6077.3273320 -44414.98696 4658.93585 6083.5737330 -4078.79210 4686.60389 6089.8202440 -4042.59725 4714.27194 6096.0667670 -3934.21916 4797.07062 6114.85119由此可知， 时，飞行器的位置参数为 70ts 394.216 ,79.02,614.8595.2.2 建立 模型对距离进行灰色预测(1,)GM建立 模型 4的方法如下：(,)设 为原始非负时间序列， 为累加生成序列，即：0(0)(0)1,2,xxn (1)xt图 4 飞行器坐标拟合图三0 5 10 15 20 25 30 35 406.056.066.076.086.096.16.116.126.13x 106不不Z不不不不不不不不不不不不不不不不9(1)(0)1,2,itxin 2的白化微分方程为：(1,)GM(1)(1)dxtatu3式 中， 为待辨识参数，亦称发展系数； 为待辨识内生变量，亦称为灰作用量，2a设待辨识向量， ，按最小二乘法得 式中u1()TaBy(1)(1)()()(1)(1)232xBxn 00(2)3()xyxn于是可得到灰色预测的 离散时间相应函数为：(1)(0)1atuxtxe (3)为所得累加的预测值，将预测值还原即为：(1)xt(0)(1)(1),2,)tttn (4)计算拟合位置与测得位置之间的距离作为原始时间序列，即 000001,2,3,4,5xxx3.764.8.17.86231.根据灰色预测模型得到 时距离残差的预测值如表 3 所示：ts:表 3：距离残差预测值（单位： ）km时间 （s ）t距离残差0 3.764110 22.402320 25.920530 29.991340 34.701350 40.151060 46.456670 53.7524310由此可知， 时，飞行器的位置为以 为圆70ts 394.256,7.82,614.5心，以距离残差值预测值 为半径的球体内，位置可靠度即为 。53.24km37km表 4：灰色预测模型的检验：平均相对误差：51(0.492+.70.25+.1037).65%4k :5.3 模型 3 的建立与求解5.3.1 优选方法的确定考 虑 卫 星 空 间 分 布 对 精 度 的 影 响 ， 对 同 等 数 量 的 卫 星 ， 选 取 两 组 分 布 差 别 较大 的 组 合 来 计 算 定 位 精 度 。 用 问 题 （ 1） 的 方 法 联 立 方 程 组 即 可 求 得 飞 行 器 的 位置 参 数 ， 再 比 较 与 实 际 位 置 的 距 离 ， 得 出 其 精 度 。 不 同 数 量 的 卫 星 组 合 得 到 的 飞 行器 位 置 参 数 ， 见 表 5：表 5 不 同 卫 星 组 合 位 置 参 数 对 照 表卫 星 组 合 位 置 参 数 (单 位 :km) 距 离 (km)(2,468)（ -5210.6254， 6608.6498， 3217.9274） 83.8937137(-5203.9396， 6605.6841， 3149.3620) 15.1494,59(-5198.7077， 6604.1294， 3122.9055) 11.9927()(-5202.8658， 6605.9736， 3166.5278) 31.8727,48,(-5198.2643， 6605.1247， 3124.3896) 10.5407237(-5201.3893， 6605.7027， 3157.3946) 22.6420(1,6,)(-5200.2716， 6605.2229， 3139.0540) 4.27075(-5202.5331， 6605.5371， 3147.9762) 13.4573,4,8(-5202.4130， 6605.8972， 3150.5806) 16.0086(279)(-5198.5059， 6604.0057， 3121.0125) 13.9000实际数据 模拟数据 残差 相对误差时间 (0)xk)(0kx)()()(00kx )(0kxk1t22.6628 22.4023 0.2605 1.1492%221.1204 20.6404 0.4800 2.2727%3037.8462 37.0607 0.7855 2.0755%4t31.4504 31.7755 0.3249 1.0337%11由 上 表 数 据 得 出 相 同 数 量 卫 星 组 合 测 得 的 距 离 误 差 平 均 值 ， 见 表 6：表 6 不 同 数 量 卫 星 组 合 距 离 误 差 均 值 表卫 星 数 4 5 6 7 8距 离 误 差 均 值 km49.5216 21.9327 16.5914 8.8640 14.95434 4.5 5 5.5 6 6.5 7 7.5 85101520253035404550不不不不不不不不不不不km不不不不不不不不不不不不data 3从 图 我 们 可 以 看 出 卫 星 组 合 得 到 的 距 离 值 最 小 ， 表 示 由 这 组 卫1,2346,89星 所 求 得 的 飞 行 器 位 置 与 最 优 位 置 的 距 离 最 小 ， 所 以 可 以 认 为 此 组 合 卫 星 测 得 的飞 行 器 位 置 参 数 精 度 最 高 ， 同 时 观 察 其 他 卫 星 组 合 所 得 到 的 距 离 值 ， 可 以 看 出 由 七颗 卫 星 组 合 得 到 的 另 一 个 距 离 值 13.4573km 亦 较 小 ， 于 是 认 为 七 颗 卫 星 组 合 得 到的 飞 行 器 位 置 参 数 精 度 较 高 ， 分 析 定 位 的 效 率 ， 随 着 定 位 卫 星 数 量 的 增 加 ， 定 位 的成 本 是 不 断 增 加 的 ， 同 时 效 率 不 断 降 低 。 结 合 定 位 的 精 度 和 效 率 综 合 考 虑 ， 可 认 为 ，当 可 供 选 择 的 定 位 卫 星 数 量 较 多 时 ， 选 择 5 颗 进 行 定 位 是 最 优 的 定 位 数 量 。在 对 定 位 卫 星 位 置 进 行 优 化 ， 对 相 同 数 量 不 同 经 度 的 两 组 卫 星 分 析 得 知 ， 当卫 星 在 赤 道 上 空 分 布 越 均 匀 时 ， 误 差 越 小 。 见 表 7：图 5：不同卫星数组合得到的距离误差均值12表 7 不 同 数 量 卫 星 组 合 经度差值表卫星数 卫星组合 相邻卫星经度差 距离(2,468)（36,11,25） 83.89374 颗卫星 13（34,31,19） 15.1494,59（34,20,12,28） 11.99275 颗卫星 (7)（36,17,19,9） 31.8727,48,（34,15,17,19,9） 10.54076 颗卫星 23（21,20,12,19,9） 22.642(1,6,9)（13,21,15,11,25,9） 4.27077 颗卫星 57（13,21,15,5,12,19） 13.4573,4,8（13,21,15,5,6,6,19,9） 16.00868 颗卫星 (2)（13,21,36,5,6,6,19,9） 13.9说 明 ： 相邻卫星经度差（36,11,25）中的 36 是指编号 2、4 两个卫星的经度的差值，11 是指编号 4、6 两个卫星的经度的差值，25 是指编号6、8 两个卫星的经度的差值，其余类推。经 过 以 上 分 析 可 知 ， 卫 星 数 量 越 多 ， 卫 星 在 赤 道 上 分 布 越 均 匀 ， 所 测 得 的 误 差值 越 小 。 经 比 较 ， 选 取 5 颗 夹 角 在 之 间 的 卫 星 时 ， 定 位 精 度 和 定 位 效 率 最152优 。5.3.2 检测误差角误差限为 0.1 度时飞行器的定位在多次测量中，设备的测量误差是正态分布的，而且均值为 。问题中给出的精0度 可以认为是测量误差的标准差。在这种理解下，用各自的误差限对测量误差进行0.1无量纲化( 也可以看成是一种加权法)处理是合理的，即求解的无约束优化问题更合理。由 5.3.1 的分析可知，当五颗卫星均匀分布时，定位精度和效率最优，故选取优化建模如下：Min2251(,)iiiiiiExyz(5)其中 1122cos(,345)(,)iiiiPXd由于目标函数是平方和的形式，因此这是一个非线性最小二乘拟合问题 5。建立优化模型，编写 程序，求得飞行器坐标为 ，这个解对应LINGO5198.,604.1,32.13的目标函数值大约为 。优化建模这个误差数值较大，这是因为这里使,3.78Exyz用的是相对于误差限 的误差，分母 很小所以相对误差很大。其实，可以认为此时的目标函数是 5 个标准正态分布的误差平方和，只要在 5 以内都是合理的。6 模型的评价及改进6.1 模型的优点模型（1）建立球心坐标系使求解过程简单明了。模型（2）在（1）问题的分析基础上利用同一时刻给出的三颗卫星的参数算出各时刻对应的飞行器位置参数，然后利用最小二乘法拟合成曲线，推算出运动模型。并用 模型预测拟合位置与实际位置距离，则其实际位置等于拟合值跟预测值的,GM和，将误差的变化考虑进模型，提高了定位的精度。模型 在问题的分析过程中，考虑题目所给的条件，将不同数量的卫星组合得到3的不同的飞行器位置，再与飞行器实际位置相比较，得出不同组合方式的测量精度。通过多组的比较，选出最优组合，整个问题的求解过程中，大量运用了的图表，使运算过程与运算结果简单明了。6.2 模型的不足模型（1）的求解过程中没有考虑到那些个别会引起测量偏差较大的卫星对结果的影响，即没有对卫星进行优选。模型 仅仅利用三颗卫星提供的参数得出飞行器在不同时刻的位置参数本身就存2在较大误差，解的过程中我们发现飞行器的位置参数在 轴上的数据残差过大。Z模型（3）直接将问题（1）中 9 颗卫星组合算出的飞行器位置作为实际位置有所欠妥，同时没有充分考虑所有卫星组合，对结果产生一定的误差。6.3 模型的改进方向模型（2）求解过程中，将的飞行器各轴上的位置参数分别进行数据拟合，最后得到 时飞行器的位置在一个球体范围内，难以定位，可综合各轴位置参数进行拟70ts合得出飞行器的空间轨迹随时间得变化规律，可得预测到更为精确直观的飞行器位置参数。模型 3 中，考虑到选取卫星组合的人为主观因素，可尽量多的选取卫星组合对飞行器进行定位，然后对比分析出各定位精度大小与定位效率，而得出最优解。14参考文献1 孙仲康 郭福成 冯道旺. 单站无源定位跟踪技术. 国防大学出版社. 2008.11 第 1 版2 徐天园. 单站无源定位算法研究. 杭州电子科技大学硕士学问论文. 2009.123 Frank R.Giordano William P.Fox 等. 机械工业出版社. 2011.7 第 1 版4 卓金武. Matlab 在数学建模中的应用. 北京航天航空大学出版社. 2011.4 第 1 版5 百度文库.飞机的精确定位问题. 2012.1115附 录附表1：某时刻空中飞行器检测到地球同步卫星的相关数据卫星编号 经度 iiE76213 692528958 43310607 8142594 53 6956E621 71403 18858 3049729 6注：取地球半径为 ，同步卫星高度为 ； 表示东京。63KmKmE附表 2：空中飞行器连续检测到地球同步卫星的相关数据检测时间（检测间隔10s）卫星编号 经度 ii1 83E8.453.042 169169760t3 72851 22 3.8.1t3 1E4391 89.56.42 617720t3 721 382892 1E913t3 5431 8.67.2 36198240t3 1720注：取地球半径为 ，同步卫星高度为 ； 表示东京。67Km580KmE16非线性方程求解程序：function F=f1(x0)a=x0(1);b=x0(2);c=x0(3);d=x0(4);e=x0(5);f=x0(6);x=x0(7);y=x0(8);z=x0(9);D=42167;e1=(a*(x-D*cos(76*pi/180)+b*(y-D*sin(76*pi/180)+c*z)/(x2+y2+z2+D2-2*D*(x*cos(76*pi/180)+y*sin(76*pi/180)0.5-cos(21.33*pi/180);e2=(a*(x-D*cos(89*pi/180)+b*(y-D*sin(89*pi/180)+c*z)/(x2+y2+z2+D2-2*D*(x*cos(89*pi/180)+y*sin(89*pi/180)0.5-cos(35.82*pi/180);e3=(a*(x-D*cos(110*pi/180)+b*(y-D*sin(110*pi/180)+c*z)/(x2+y2+z2+D2-2*D*(x*cos(110*pi/180)+y*sin(110*pi/180)0.5-cos(60.73*pi/180);e4=(d*(x-D*cos(76*pi/180)+e*(y-D*sin(76*pi/180)+f*z)/(x2+y2+z2+D2-2*D*(x*cos(76*pi/180)+y*sin(76*pi/180)0.5-cos(69.25*pi/180);e5=(d*(x-D*cos(89*pi/180)+e*(y-D*sin(89*pi/180)+f*z)/(x2+y2+z2+D2-2*D*(x*cos(89*pi/180)+y*sin(89*pi/180)0.5-cos(54.63*pi/180);e6=(d*(x-D*cos(110*pi/180)+e*(y-D*sin(110*pi/180)+f*z)/(x2+y2+z2+D2-2*D*(x*cos(110*pi/180)+y*sin(110*pi/180)0.5-cos(29.81*pi/180);e7=(d*(x-D*cos(125*pi/180)+e*(y-D*sin(125*pi/180)+f*z)/(x2+y2+z2+D2-2*D*(x*cos(125*pi/180)+y*sin(125*pi/180)0.5-cos(11.99*pi/180);e8=(a*(x-D*cos(125*pi/180)+b*(y-D*sin(125*pi/180)+c*z)/(x2+y2+z2+D2-2*D*(x*cos(125*pi/180)+y*sin(125*pi/180)0.5-cos(79.24*pi/180);e9=(a*(x-D*cos(130*pi/180)+b*(y-D*sin(130*pi/180)+c*z)/(x2+y2+z2+D2-2*D*(x*cos(130*pi/180)+y*sin(130*pi/180)0.5-cos(85.46*pi/180);e10=(d*(x-D*cos(130*pi/180)+e*(y-D*sin(130*pi/180)+f*z)/(x2+y2+z2+D2-2*D*(x*cos(130*pi/180)+y*sin(130*pi/180)0.5-cos(6.95*pi/180);e11=(a*(x-D*cos(136*pi/180)+b*(y-D*sin(136*pi/180)+c*z)/(x2+y2+z2+D2-2*D*(x*cos(136*pi/180)+y*sin(136*pi/180)0.5-cos(92.91*pi/180);e12=(a*(x-D*cos(136*pi/180)+b*(y-D*sin(136*pi/180)+c*z)/(x2+y2+z2+D2-2*D*(x*cos(136*pi/180)+y*sin(136*pi/180)0.5-cos(5.99*pi/180);e13=(d*(x-D*cos(142*pi/180)+e*(y-D*sin(142*pi/180)+f*z)/(x2+y2+z2+D2-2*D*(x*cos(142*pi/180)+y*sin(142*pi/180)0.5-cos(11.58*pi/180);e14=(a*(x-D*cos(142*pi/180)+b*(y-D*sin(142*pi/180)+c*z)/(x2+y2+z2+D2-2*D*(x*cos(142*pi/180)+y*sin(142*pi/180)0.5-cos(100.32*pi/180);e15=(a*(x-D*cos(163*pi/180)+b*(y-D*sin(163*pi/180)+c*z)/(x2+y2+z2+D2-2*D*(x*cos(163*pi/180)+y*sin(163*pi/180)0.5-cos(125.58*pi/180);e16=(d*(x-D*cos(163*pi/180)+e*(y-D*sin(163*pi/180)+f*z)/(x2+y2+z2+D2-172*D*(x*cos(16