坐标转换参数求取及坐标转换程序设计教材

1、毕业设计设计题目 坐标转换参数求取及坐标转换程序设计学生姓名张威指导教师杜继亮专 业测绘工程班 级 测绘 12-2 班 填写日期2016/4/6矿业工程学院摘要坐标系统是测量工作中定位的基础，坐标系统有多种形式和基准，由于 各测量工作目的不同，所选用的坐标基准也会不同，根据不同的工作要求需 要将不同的坐标系下的坐标进行相互转换。在这些坐标转换的过程中会用到 很多坐标转换模型，但是坐标系转换模型过于复杂手算非常困难。本设计为 了方便施工时遇到的坐标转换问题，设计利用 Visual Basic 6.0 编程语言编写 程序，用来实现坐标系统之间的转换以及转换参数的求解，例如:大地坐标与空间直角坐标的

2、相互转换、高斯投影正反算、二维坐标转换与四参数计算、 三维坐标转换与七参数转换、同参考基准下的坐标换带计算，以及坐标数据 的批量处理。关键字： 坐标系统，转换模型，坐标转换，程序设计AbstractThe base of coordinate system in surveying work. there are many forms and benchmarks in the coordinate system. However, in general engineering, the control point and coordinate. System are the same. So

3、 It is necessary to transform the control point. coordinate during the construction process. Due to different purposes of each measurement and the selected. different coordinate references, there will be many different coordinate systems. Coordinate systems used in the measurement work are as follow

4、s: WGS-84 World Geodetic System, China Geodetic Coordinate System 2000, National Geodetic Coordinate System 1980, Beijing coordinate system 1954 and Local Coordinate System. There are space rectangular coordinate, geodetic coordinate and plane coordinate in the way of the reference in the same coord

5、inate. According to the requirements of different tasks, we need to convert coordinates under the different coordinate systems . On condition that the coordinates of the reference standard can be obtained. the normal construction work can be done. A lot of coordinate transformation models are used i

6、n the process of the coordinate transformation. But the coordinate transformation model is very complex and difficult. Nowadays the conversion formula is suitable for the computerization whose language is easy to learn. So in the design I make use of Visual Basic 6 programming language to realize th

7、e transformation between the coordinate system and transformation parameters.Key words : coordinate systems transformation model coordinate transform programmingII目录摘 要 IAbstract II第 1 章 绪论 11.1 研究背景和意义 11.2 国内外研究现状 11.3 研究的主要内容 21.4 程序设计思路方法 2第 2 章 坐标理论和知识 . 32.1 地球椭球 32.2 基准 32.3 测量常用坐标系 42.3.1 大地

8、坐标系 42.3.2 空间直角坐标系 42.3.3 平面坐标系 52.3.4 地方独立坐标系 52.4 我国常用坐标系 62.4.1 1954 年北京坐标系 62.4.2 1980 国家大地坐标系 62.4.3 WGS-84 世界大地坐标系 62.4.4 2000 国家大地坐标系 7第 3 章 坐标转换程序设计 . 143.1 大地坐标系与空间直角坐标系的转换 153.1.1 大地坐标转换成空间直角坐标 错误！未定义书签。3.1.2 空间直角坐标转成换大地坐标 错误！未定义书签。3.2 高斯平面坐标与大地坐标的转换 17III3.2.1 高斯正算 错误！未定义书签。3.2.2 高斯反算 错误！

9、未定义书签。3.3 高斯投影换带计算 183.4 二维坐标转换与四参数计算 203.5 三维坐标转换与七参数计算 223.6 三、六度带带号与中央子午线计算 23第 4 章 设计总结 27致谢 33参考文献 34IV第 1章绪论1.1 研究背景和意义随着大地测量学，摄影测量学的发展和电子计算机的普及，对各种坐标 系的研究变得越来越重要了 1。随着现在社会的快速发展，各种各样的大型 工程的建设，凡是工程施工就必定需要坐标来定位才可以，而建造的地方又 不同，又没有一个满足全部地方的坐标系统， 所以产生多种不同的坐标系统， 实际工作中测量人员必定会按照实际情况选择最为合适的坐标系统。而坐标 系统之间

10、的转换比较复杂，手算工作量巨大，因而各种坐标转换模型相继出 现，利用计算机强大的数据计算能力可以轻松应对这些问题， 提高工作效率。 坐标转换的意义重大，不仅在我们熟知的工程领域中，在国防建设、航空航 天科技、城市汇划等众多领域中都发挥着重要的作用，可以说对社会进步有 着必不可少的作用。1.2 国内外研究现状自 60 年代以来，各国大地测量学者，经过大量研究，提出了多种坐标 转换模型及多种解算方法，北美 1927 基准面（基于克拉克 1966 椭球体与北 美 1983 基准面（基于 GRS1980椭球体）之间坐标转换是根据研究区内一系 列己知点的大地坐标或网格坐标改正量进行插值进行的坐标系转换；

11、英国采 用北向与东向的双线性网格插值进行坐标转换；挪威在海岸带调查中，采用 经纬度多项式用于坐标系转换这种方法进行新 （ED87 一欧洲 1987 基准面）、 旧（ ED50 一欧洲 1950 基准面）坐标系之间的转换：欧洲石油勘探组织（EPSG） 对新、旧坐标系采用“双线性插值” 进行坐标转换 2。国内空间三维直角坐标转换中，一般采用布尔莎七参数模型。一般有 7 个转换参数，即 3 个平移参数， 3 个旋转参数和 1 个尺度参数。需要三个及 已经公共点时，才能利用平差的方法求出七参数。1.3 研究的主要内容本坐标转换程序可实现功能有 :1、大地坐标与空间直角坐标的相互转换， 2、高斯投影正反

12、算， 3、二维坐标转换与四参数计算， 4、三维坐标转换与 七参数转换， 5、同参考基准下的坐标换带计算，以及坐标数据的批量处理。1.4 程序设计思路方法本程序名为万能坐标转换器。设计前期收集相关资料，参考一些成熟的 坐标转换软件，确定程序应有的功能以及界面设计。运用 VB 编写程序时， 查阅相关书籍获取理论知识以及转换模型。完成程序后将已知正确数据带入 其中验证程序结果是否正确。若出现错误则检查每步代码，直到程序完美运 行为止。第 2 章 基础知识准备2.1 地球椭球地球椭球体又称地球椭圆体或地球扁球体，代表地球大小和形状的数学 曲面，以长半径和扁率表示，因它十分迫近于椭球体，故通常以参考椭球

13、体 表示地球椭球体的形状和大小。通常所说地球的形状和大小，实际上就是以 参考椭球体的半长径、半短径和扁率来表示。 1975 年国际大地测量与地球物 理联合会推荐的数据为：半长径 6378140米，半短径 6356755 米，扁率 1 298.257。在众多椭球体中， WGS-84 椭球体被认为符合上述条件最好的椭球。2.2 基准所谓基准是指为描述空间位置而定义的点、线和面。而大地测量基准是 指用以描述地球形状的地球椭球参数，包含描述地球椭球几何特征的长短半 轴和物理特征的有关参数、地球在空间的定位及定向以及描述这些位置所采 用的单位长度的定义。不同的坐标系统会使用的基准也不同，根据参考椭球 所

14、选原点位置不同，可以分为地心坐标系和参心坐标系。地心坐标系是以地球的质心为原点，有地心大地坐标系和地心空间直角 坐标系两种表述方法。地心空间直角坐标系的定义为：以地球质心为原点， X 轴指向格林尼治子午面与地球赤道的交点， Z 轴指向北极， Y 轴过原点垂 直于平面 XOZ ，构成右手空间直角坐标系。地心大地坐标系定义为：以地 球的质心作为原点，以地球自转轴作为椭球的短轴，大地纬度 B 是过地面点 的椭球法线与椭球赤道面之间的夹角，大地经度 L 为过地面点的椭球子午面 与格林尼治子午面之间的夹角， 大地高度 H 为地面点沿椭球法线到椭球面的 最短距离。 WGS-84 坐标系， CGCS2000

15、 坐标系， GLONASS 是采用 PZ-90 坐标，都是属于地心坐标系。参心坐标系是选取一个参考椭球面作为基本的参考面，选一参考点作为 大地测量的起算点，从而确定参考椭球在地球面的位置和方向。这时参考椭 球的原点不会和地球质心重合，所以称为参心。北京 54 坐标系、西安 80 坐 标系和新北京 54 坐标系，都是参心坐标系。它同样具有参心大地坐标系和 参心空间直角坐标系两种表述方法，它们的定义与地心坐标系的定义相似。2.3 测量常用坐标系2.3.1 大地坐标系空间一点的大地坐标用大地经度 L、大地纬度 B 和大地高度 H 表示，地 面上 P 地点的大地子午面 NPS 与起始 大地子午面 所构

16、成 的二面角 L 称 P 地 点的大地经度， 如图 2-1 所示P 地点对于椭球的法线与赤道面的夹角 B 称 P 地点的大地纬度图 2-1 大地坐标系P 地点沿法线到椭球面的距离 H 称大地高，从椭球面起算，向外为正， 向内为负 3 。H = H 正常 + （高程异常）H = H 正 + N（大地水准面差距）2.3.2 空间直角坐标系空间直角坐标系的坐标原点与参考椭球的中心重合， Z 轴正向指向参考 椭球的北极， X 轴正向指向起始子午面与赤道的交点， Y 轴按右手系与 X 轴 呈 90夹角且位于赤道面上，某点在空间中的坐标可用该点在此空间坐标系2.3.3 平面坐标系平面直角坐标系是利用投影，

17、将空间坐标通过某种数学变换映射到平面 上，这种变换称为投影变换 5 。在我国一般采用的是高斯一克吕格投影，是 目前测量上广泛采用的正形投影，特点是没有角度变形，在不同点上的长度 比随点位而异，但在同一点上各方向的长度比相同，也称为高斯投影 62.3.4 地方独立坐标系在我国平面坐标主要采用的是高斯投影， 在该投影中， 除中央子午线外， 其它位置上的任何线段，投影后都会产生一定的长度变形，而且变形随离开 中央子午线的距离增加而增加 7 。因此一般采用分带投影的办法，来限制长 度变形，我国规定了采用 3 度带或 6 度带进行分带投影。在城市、工矿等工 程测量中，如果直接在国家分带坐标系中建立控制网

18、，会使地面长度投影的 变形较大，当长度变形大于 25 cm/km 时，就难以满足工程上的需要 8。另 一些特殊的测量，比如大桥施工测量，水利水坝测量，滑坡变形监测等，采 用国家坐标系精度不能满足工程要求，所以常常会建立适合本地区的地方独 立坐标系 9 。2.4 我国常用坐标系2.4.1 1954 年北京坐标系1954 年北京坐标系，是采取先将我国一等锁与前苏联远东一等锁相联 接，然后以连接处呼玛，吉拉林，东宁基线网扩大边端点的前苏联 1942 年 普尔科沃坐标系的坐标为起算数据，平差我国东北及东部一等锁 10 。椭球参 数：长半轴 a=6378245m，短半轴 b=6356863.0188m

19、，扁率 =1/298.3，第一 偏 心 率 平 方 e2=0.00669342161454287 ， 第 二 偏 心 率 平 方 e2=0.0067385254061465。22.4.2 1980 国家大地坐标系1980 国家大地坐标系采用地球椭球基本参数为 1975 年国际大地测量与 地球物理联合会第十六届大会推荐的数据，大地原点设在我国中部的陕西省 泾阳县永乐镇，位 于西 安市 西北 方向 约 60 公里 。椭球参数：长半轴 a=63781405m，短半轴 b=6356755.2882m，扁率 =1/298.257 ，第一偏心率 平方 e2 =0.00669438499959，第二偏心率平

20、方 e2=0.0067395018194。72.4.3 WGS-84 世界大地坐标系WGS-84 坐标系是一种国际上采用的地心坐标系。 坐标原点为地球质心，其地心空间直角坐标系的 Z 轴指向 BIH （国际时间服务机构） 1984.0 定义6的协议地球极（ CTP）方向， X 轴指向 BIH 1984.0 的零子午面和 CTP 赤道的 交点， Y 轴与 Z 轴、 X 轴垂直构成右手坐标系，称为 1984 年世界大地坐标 系统 。椭球参数： 长半轴 a=6378137 ，短半轴 b=6356752.3142 ，扁率 =1/298.2572236，第一偏心率平方 e2=0.006694379990

21、13 ，第二偏心率平方 e2=0.00673949674222。72.4.4 2000 国家大地坐标系CGCS2000 是（中国） 2000 国家大地坐标系的缩写，该坐标系是通过中 国 GPS 连续运行基准站、 空间大地控制网以及天文大地网与空间地网联合平 差建立的地心大地坐标系统 11 。Z 轴指向 BIH1984.0 定义的协议极地方向 （ BIH 国际时间局）， X 轴指向 BIH1984.0 定义的零子午面与协议赤道的交 点，Y 轴 按右 手坐标系确定。椭 球参 数：长半轴 a=6378137， 短半轴 b=6356752.31414 ， 扁 率 =1/298.2572236 ， 第

22、一 偏 心 率 平 方 e2=0.00669437999013 ，第二偏心率平方 e 2=0.00673949674222。72.5 坐标转换模型2.5.1 大地坐标系与空间直角坐标系转换模型将同一坐标系下的大地坐标（ B， L ，H）转换成空间直角坐标（ X， Y， Z）的转换公式为：式中： a为参考椭球长半轴， e为第一偏心率， N 为卯酉圈的半径 12。 将同一坐标系下的空间直角坐标（ X， Y，Z）转换为大地坐标（ B， L ， H ）的公式为 :L=arctan( )B=arctan( )H=用公式进行空间直角坐标转换大地坐标时，需要采用迭代计算大地纬度B 。具体计算时，可先根据下式

23、求出大地纬度B 的初值：因为 , 带入 B=arctan( )得:令:则式子可写成然后利用该初值 B 代入公式右端 tanB 中，将等式左边的结果再次代入 右端 tanB，直到最后两次 B 值之差小与允许误差为止。当得到大地纬度 B 后，代入公式即可求出大地高 H 。2.5.2 高斯正反算转换模型得到了点的大地坐标 （L ，B），就可以将其转化为某投影带的高斯平面坐标，我们将椭球参数代入高斯投影正算公式得到更适用于电算的高斯坐标计 算的实用公式：其计算结果的精度可达 0. 001m 。只要得到了高斯平面坐标 (X, Y) 后，便可通过高斯反算公式将其转换成 大地坐标 (B, L) ，高斯投影反

24、算公式为 :9它们的计算精度 ,即平面坐标可达2.5.3 坐标转换与参数计算转换模型二维转换模型：2 个平移参数 （原点不重合产生的）1 个旋转参数 （坐标轴不平行产生的）1 个尺度参数 （两坐标系间的尺度不一致产生的）为某点在 A 空间直角坐标系中的坐标。 为某点在 B 空间直 角坐标系中的坐标。 为某点从 A 空间直角坐标系转换到 B空间直角坐 标系的两个平移参数。 为从 A 空间直角坐标系转换到 B空间直角坐标系中 标系的两个平移参数。 为从 A 空间直角坐标系转换到 B空间直角坐标系中10一个旋转参数。 m 为从 A 空间直角坐标系转换到 B 空间直角坐标系中的一个 尺度参数。平面四参

25、数求解步骤如下 : 利用公共点计算坐标参数，但至少有两个公共点，当有 i 个公共点时， 可利用最小二乘原理求解参数。将 B直角坐标系中的坐标视为观测值，设 A 直角坐标系下的坐标视为无 误差，列误差方程为：写成矩阵形式既：由于各点的坐标可视为同精度独立观测值，因此P=I。把各点坐标已知值带入上述误差方程，然后按下列公式求解出四参数：三维转换模型：3 个平移参数（原点不重合产生的）3 个旋转参数（坐标轴不平行产生的）1 个尺度参数 （两坐标系间的尺度不一致产生的）11般 为微小转角，模型可简化为：为某点在 A 空间直角坐标系中的坐标。为某点在 B 空间直角坐标系中的坐标。 为某点从 A 空间直角

26、坐标系转换到 B 空 间直角坐标系的三个平移参数。 为从 A 空间直角坐标系转换到 B 空间直角坐标系中三个旋转参数。 m 为从 A 空间直角坐标系转换到 B 空间直 角坐标系中的一个尺度参数。七参数求解步骤如下 : 利用公共点计算坐标参数，但至少有三个公共点，当有 i 个公共点时， 可利用最小二乘原理求解参数。将 B 空间直角坐标系中的坐标视为观测值， 设 A 空间直角坐标系下的坐 标视为无误差，列误差方程为：12写成矩阵形式既：由于各点的坐标可视为同精度独立观测值，因此P=I。把各点坐标已知值带入上述误差方程，然后按下列公式求解出七参数：2.5.4 三、六度带带号与中央子午线计算模型中央子

27、午线中央子午线L 为当地经度L 为当地经度求带号、中央子午线公式为：3带:带号6带:带号13第 3 章 坐标转换程序设计图 3-1 编写程序流程图14图 3-2 程序主界面3.1 大地坐标系与空间直角坐标系的转换主要代码如下：椭球参数设置If Combo1.Text = WGS-84 Thena1 = 6378137: b1 = 6356752.3142 a1为长半轴， b1 为短半轴 Text7.Text = a1: Text8.Text = b1ElseIf Combo1.Text = CGCS2000 Thena1 = 6378137: b1 = 6356752.31414Text7.T

28、ext = a1: Text8.Text = b1ElseIf Combo1.Text = 北京 54 Thena1 = 6378245: b1 = 6356863.0188Text7.Text = a1: Text8.Text = b1ElseIf Combo1.Text = 西安 80 Then15a1 = 6378140: b1 = 6356755.2882Text7.Text = a1: Text8.Text = b1ElseIf Combo1.Text = 自定义 ThenText7.Enabled = True: Text8.Enabled = Truea1 = Text7.Tex

29、t: b1 = Text8.TextEnd If 大地坐标转空间直角坐标主要代码B = Text1.Text: L = Text2.Text: H = Text3.Texte12 = (a1 2 - b1 2) / a1 2N = a1 / Sqr(1 - e12 * (Sin(DuToHu(B) 2)x = (N + H) * Cos(DuToHu(B) * Cos(DuToHu(L) y = (N + H) * Cos(DuToHu(B) * Sin(DuToHu(L) z = (N * (1 - e12) + H) * Sin(DuToHu(B) Text4.Text = x: Text

30、5.Text = y: Text6.Text = z 空间直角坐标转大地坐标主要代码x = Text4.Text: y = Text5.Text: z = Text6.TextL = HuToDu(Atn(y / x) 计算经度If y / x 0 ThenL = HuToDu(Atn(y / x) + 3.14159265358979)End Ife12 = (a1 2 - b1 2) / a1 2 第一偏心率平方 TB = z / Sqr(x 2 + y 2) 初始纬度 B 正切值 c1 = (a1 2) / b1 t = z / Sqr(x 2 + y 2) p = (c1 * e12)

31、 / Sqr(x 2 + y 2) k = 1 + (a1 2 - b1 2) / b1 2TB1 = t + (p * TB) / Sqr(k + TB 2) 用 TanB 作为未知量 Do Until Abs(TB1 - TB) 0.000000000001 迭代法 TB = TB1TB1 = t + (p * TB) / Sqr(k + TB 2)LoopB = Atn(TB)N = a1 / Sqr(1 - e12 * (Sin(B) 2) 卯酉圈曲率半径H = z / Sin(B) - N * (1 - e12) 大地高Text1.Text = HuToDu(B): Text2.Te

32、xt = L: Text3.Text = H163.2 高斯平面坐标与大地坐标的转换高斯平面坐标 (X, Y) 与大地坐标 (B ,L)的相互关系式分为两类：第一类称 高斯投影正算公式，由 (B,L)求(X,Y) ；第二类称高斯投影反算公式，由 (X,Y) 求(B,L) 13 。由于电子计算机和各种可编程序电子计算器在测量上广泛应用， 因此在这里给出适合电算的高斯正反算公式。主要代码如下：正算Dim B, L, x, y As DoubleDim N, a1, b1, c, e12, e22, e1, e2, ll, l0, nn, t As DoubleDim XX, bb0, bb2, b

33、b4, bb6, bb8 As DoubleDim yj As DoubleDim mf, nf, nnf, tf, bf, bf1 As DoubleB = Text1.Text: L = Text2.Textl0= Text5.Text 输入中央子午线经度yj = Text8.Text Y 坐标加常数a1 = Text6.Text: b1 = Text7.Text 椭球长短半轴e12 = (a1 2 - b1 2) / a1 2 第一偏心率平方e22 = (a1 2 - b1 2) / b1 2 第二偏心率平方e1 = Sqr(a1 2 - b1 2) / a1 2) 第一偏心率e2 =

34、Sqr(a1 2 - b1 2) / b1 2) 第二偏心率c = a1 / Sqr(1 - e12) 极点子午线曲率半径bb0 = 1 - (3 / 4) * e22 + (45 / 64) * e22 2 - (175 / 256) * e22 3 + (11025 / 16384) * e22 4bb2 = bb0 - 1bb4 = (15 / 32) * e22 2 - (175 / 384) * e22 3 + (3675 / 8192) * e22 4bb6 = -(35 / 96) * e22 3 + (735 / 2048) * e22 4bb8 = (315 / 1024)

35、* e22 4ll= (Du(L) - l0) / 57.295779513N = a1 / Sqr(1 - e12 * (Sin(DuToHu(B) 2) 卯酉圈曲率半径t = Sin(DuToHu(B) / Cos(DuToHu(B)nn = e2 * Cos(DuToHu(B)17XX = c * (bb0 * DuToHu(B) + (bb2 * Cos(DuToHu(B) + bb4 (Cos(DuToHu(B) 3) + bb6 * (Cos(DuToHu(B) 5) + bb8 * (Cos(DuToHu(B) 7) * Sin(DuToHu(B) 子午线弧长x = XX + (

36、N / 2) * Sin(DuToHu(B) * Cos(DuToHu(B) * (ll 2) + (N / 24)* Sin(DuToHu(B) * (Cos(DuToHu(B) 3) * (ll 4) * (5 - t 2 + 9 * (nn 2) +4 * (nn 4) + (N / 720) * Sin(DuToHu(B) *(Cos(DuToHu(B) 5) * (ll 6) *(61 - 58 * (t 2) + (t 4)y = N * Cos(DuToHu(B) * ll + (N / 6) * (Cos(DuToHu(B) 3) * (ll 3) * (1 - t 2 + n

37、n 2) + (N / 120) * (Cos(DuToHu(B) 5) * (ll 5) * (5 - 18 * (t 2) + t 4 + 14 * (nn 2) - 58 * (t 2) * (nn 2)反算x = Text3.Text: y = Text4.Texty = y - yjbf = x / (c * bb0) 迭代Dobf1 = (x / c) - (bb2 * Cos(bf) + bb4 * (Cos(bf) 3) + bb6 * (Cos(bf) 5) + bb8 * (Cos(bf) 7) * Sin(bf) / bb0Loop Until Abs(bf - bf1)

38、 0.000000001mf = (a1 * (1 - e12) / (Sqr(1 - e12 * (Sin(bf) 2) 3)tf = Sin(bf) / Cos(bf)nnf = e2 * Cos(bf)nf = a1 / Sqr(1 - e12 * (Sin(bf) 2)B = bf - (tf * (y 2) / (2 * mf * nf) + (tf * (y 4) * (5 + 3 * (tf 2) + (nnf 2) - 9 * (nnf 2) * (tf 2) / (24 * mf * (nf 3) - (tf * (y 6) * (61 + 90 * (tf 2) + 45

39、* (tf 4) / (720 * mf * (nf 5)L = y / (nf * Cos(bf) - (1 + 2 * tf 2 + nnf 2) * (y 3) / (6 * nf 3 * Cos(bf) + (5 + 28 * (tf 2) + 24 * (tf 4) + 6 * (nnf 2) + 8 * (nnf 2) * (tf 2) * (y 5) / (120 * (nf 5) * Cos(bf)3.3 高斯投影换带计算高斯投影虽然保证了角度没有变形，但是长度变形比较严重，为了限制 高斯投影的长度变形，必须以中央子午线进行分带，把投影范围限制在中央 子午线东、西侧一定的范围内

40、；但是这样又使得统一的坐标系分割成各带的 独立坐标系 14 。这样就会产生新的问题，需要邻带换算来解决。18主要代码如下：l1 = Text3.Text : yy1 = Text4.Text 加常数和中央子午线x = Text9.Text : y = Text10.Text反算y = y - yy1e12 = (a1 2 - b1 2) / a1 2 第一偏心率平方e22 = (a1 2 - b1 2) / b1 2 第二偏心率平方e1 = Sqr(a1 2 - b1 2) / a1 2) 第一偏心率e2 = Sqr(a1 2 - b1 2) / b1 2) 第二偏心率c = a1 / Sqr

41、(1 - e12) 极点子午线曲率半径bb0 = 1 - (3 / 4) * e22 + (45 / 64) * e22 2 - (175 / 256) * e22 3 + (11025 / 16384) * e22 4bb2 = bb0 - 1bb4 = (15 / 32) * e22 2 - (175 / 384) * e22 3 + (3675 / 8192) * e22 4bb6 = -(35 / 96) * e22 3 + (735 / 2048) * e22 4bb8 = (315 / 1024) * e22 4bf1 = x / (c * bb0) 迭代 单位弧度Dobf = b

42、f1bf1 = (x / c) - (bb2 * Cos(bf) + bb4 * (Cos(bf) 3) + bb6 * (Cos(bf) 5) + bb8 * (Cos(bf) 7) * Sin(bf) / bb0Loop Until Abs(bf - bf1) 0.000000001mf = (a1 * (1 - e12) / (Sqr(1 - e12 * (Sin(bf) 2) 3)tf = Sin(bf) / Cos(bf)nnf = e2 * Cos(bf)nf = a1 / Sqr(1 - e12 * (Sin(bf) 2)B L 单位弧度B = bf - (tf * (y 2)

43、/ (2 * mf * nf) + (tf * (y 4) * (5 + 3 * (tf 2) + (nnf 2) - 9 * (nnf 2) * (tf 2) / (24 * mf * (nf 3) - (tf * (y 6) * (61 + 90 * (tf 2) + 45 * (tf 4) / (720 * mf * (nf 5)L = y / (nf * Cos(bf) - (1 + 2 * tf 2 + nnf 2) * (y 3) / (6 * nf 3 Cos(bf) + (5 + 28 * (tf 2) + 24 * (tf 4) + 6 * (nnf 2) + 8 * (nn

44、f 2) * (tf 2) * (y 5) / (120 * (nf 5) * Cos(bf)L = L + DuToHu(l1)正算19l2= Text7.Text: yy2 = Text8.Text 加常数和中央子午线ee12 = (a1 2 - b1 2) / a1 2 第一偏心率平方ee22 = (a1 2 - b1 2) / b1 2 第二偏心率平方ee1 = Sqr(a1 2 - b1 2) / a1 2) 第一偏心率ee2 = Sqr(a1 2 - b1 2) / b1 2) 第二偏心率nn = a1 / Sqr(1 - ee12 * (Sin(B) 2) 卯酉圈曲率半径tt =

45、 Sin(B) / Cos(B)nnn = ee2 * Cos(B)cc = a1 / Sqr(1 - ee12) 极点子午线曲率半径bbb0 = 1 - (3 / 4) * ee22 + (45 / 64) * ee22 2- (175 / 256) * ee22 3 +(11025 / 16384) * ee22 4bbb2 = bbb0 - 1bbb4 = (15 / 32) * ee22 2 - (175 / 384) * ee22 3 + (3675 / 8192) * ee22 4bbb6 = -(35 / 96) * ee22 3 + (735 / 2048) * ee22 4b

46、bb8 = (315 / 1024) * ee22 4L = HuToDu(L)ll = (Du(L) - l2) / 57.295779513XX = cc * (bbb0 * B + (bbb2 * Cos(B) + bbb4 * (Cos(B) 3) + bbb6 * (Cos(B) 5) + bbb8 * (Cos(B) 7) * Sin(B) 子午线弧长x = XX + (nn / 2) * Sin(B) * Cos(B) * (ll 2) + (nn / 24) * Sin(B) (Cos(B) 3) * (ll 4) * (5 - tt 2 + 9 * (nnn 2) + 4 *

47、 (nnn 4) + (nnn / 720) * Sin(B) * (Cos(B) 5) * (ll 6) * (61 - 58 * (tt 2) + (tt 4)y = nn * Cos(B) * ll + (nn / 6) * (Cos(B) 3) * (ll 3) * (1 - tt 2 + nnn 2) + (nn / 120) * (Cos(B) 5) * (ll 5) * (5 - 18 * (tt 2) + tt 4 + 14 * (nnn 2) - 58 * (tt 2) * (nnn 2)y = y + yy23.4 二维坐标转换与四参数计算由于矩阵计算代码繁琐，由后台调用

48、Excel 函数完成计算。 其主要代码如下：Dim x#, y#, x1#, y1#, dx#, dy#, jd#, k#正算：x1 = (Cos(DuToHu(jd) * x - Sin(DuToHu(jd) * y) * k + dxy1 = (Sin(DuToHu(jd) * x + Cos(DuToHu(jd) * y) * k + dy20反算：x = (x1 - dx) * Cos(DuToHu(jd) + (y1 - dy) * Sin(DuToHu(jd) / ky = (x1 - dx) * (-Sin(DuToHu(jd) + (y1 - dy) * Cos(DuToHu(

49、jd) / k 四参数计算：Dim xls As Excel.ApplicationDim book As Excel.WorkbookDim sheet As Excel.WorksheetSet xls = CreateObject(Excel.Application) 创建 EXCEL 对象Set book = xls.Workbooks.Open(Text22.Text) 打开已经存在的工件簿文件 Set sheet = book.Worksheets(1) sheet.Activate a = sheet.UsedRange.Rows.Countm = 1For i = 2 To a

50、x = xls.Cells(i, 2).Value: y = xls.Cells(i, 3).ValueXX = xls.Cells(i, 4).Value: yy = xls.Cells(i, 5).Valuexls.Cells(m, 10).Value = XX - x: xls.Cells(m + 1, 10).Value = yy - y 参 L xls.Cells(m, 6).Value = 1: xls.Cells(m, 7).Value = 0: xls.Cells(m, 8).Value = y: xls.Cells(m, 9).Value = xxls.Cells(m + 1

51、, 6).Value = 0: xls.Cells(m + 1, 7).Value = 1: xls.Cells(m +1, 8).Value = -x: xls.Cells(m + 1, 9).Value = ym = m + 2 Next i Range(Cells(1,11),Cells(4,2*(a+1)+10)=Application.WorksheetFunction.Tr anspose(Range(Cells(1, 6), Cells(2 * (a - 1), 9) 参数 B 转置Range(Cells(5,11),Cells(8,14)=Application.Workshe

52、etFunction.MMult(Ran ge(Cells(1, 11), Cells(4, 2 * (a - 1) + 10), (Range(Cells(1, 6), Cells(2 * (a - 1), 9) B 转置 *BRange(Cells(9,11),Cells(12,14)=Application.WorksheetFunction.MInverse( Range(Cells(5, 11), Cells(8, 14) 矩阵求逆Range(Cells(13,11),Cells(16,11)=Application.WorksheetFunction.MMult( Range(Ce

53、lls(1, 11), Cells(4, 2 * (a - 1) + 10), (Range(Cells(1, 10), Cells(2 * (a - 1), 10) B 转置 *LRange(Cells(13,13),Cells(16,13)=Application.WorksheetFunction.MMult( Range(Cells(9, 11), Cells(12, 14), (Range(Cells(13, 11), Cells(16, 11) 四参数 book.Saved = True21xls.Quit3.5 三维坐标转换与七参数计算由于矩阵计算代码繁琐，由后台调用 Excel

54、 函数完成计算 其主要代码如下：Dim x#, y#, z#, x1#, y1#, z1#, dx#, dy#, dz#, jdx#, jdy#, jdz#, k# 正算：x1 = (x * Cos(DuToHu(jdy) * Cos(DuToHu(jdz) + y * Cos(DuToHu(jdy) * Sin(DuToHu(jdz) + z * (-Sin(DuToHu(jdy) * k + dxy1 = (x * (-Cos(DuToHu(jdy) * Sin(DuToHu(jdz) + Sin(DuToHu(jdx) * Sin(DuToHu(jdy)*Cos(DuToHu(jdz)+ y * (Cos(DuToHu(jdx) * Cos(DuToHu(jdz) + Sin(DuToHu(j