毕业设计
设计题目 坐标转换参数求取及坐标转换程序设计
学生姓名 张威
指导教师 杜继亮
专 业 测绘工程
班 级 测绘12-2班
填写日期 2016/4/6
矿业工程学院
坐标系统是测量工作中定位的基础，坐标系统有多种形式和基准，由于各测量工作目的不同，所选用的坐标基准也会不同，根据不同的工作要求需要将不同的坐标系下的坐标进行相互转换。在这些坐标转换的过程中会用到很多坐标转换模型，但是坐标系转换模型过于复杂手算非常困难。本设计为了方便施工时遇到的坐标转换问题，设计利用Visual Basic 6.0编程语言编写程序，用来实现坐标系统之间的转换以及转换参数的求解，例如:大地坐标与空间直角坐标的相互转换、高斯投影正反算、二维坐标转换与四参数计算、三维坐标转换与七参数转换、同参考基准下的坐标换带计算，以及坐标数据的批量处理。
关键字：坐标系统，转换模型，坐标转换，程序设计
The 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 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 follows: 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 coordinate. 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 in 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 the transformation between the coordinate system and transformation parameters.
Key words : coordinate systems transformation model coordinate transform programming
1.1研究背景和意义
随着大地测量学，摄影测量学的发展和电子计算机的普及，对各种坐标系的研究变得越来越重要了[1]。随着现在社会的快速发展，各种各样的大型工程的建设，凡是工程施工就必定需要坐标来定位才可以，而建造的地方又不同，又没有一个满足全部地方的坐标系统，所以产生多种不同的坐标系统，实际工作中测量人员必定会按照实际情况选择最为合适的坐标系统。而坐标系统之间的转换比较复杂，手算工作量巨大，因而各种坐标转换模型相继出现，利用计算机强大的数据计算能力可以轻松应对这些问题，提高工作效率。坐标转换的意义重大，不仅在我们熟知的工程领域中，在国防建设、航空航天科技、城市汇划等众多领域中都发挥着重要的作用，可以说对社会进步有着必不可少的作用。
1.2国内外研究现状
自60年代以来，各国大地测量学者，经过大量研究，提出了多种坐标转换模型及多种解算方法，北美1927基准面（基于克拉克1966椭球体与北美1983基准面（基于GRS1980椭球体）之间坐标转换是根据研究区内一系列己知点的大地坐标或网格坐标改正量进行插值进行的坐标系转换；英国采用北向与东向的双线性网格插值进行坐标转换；挪威在海岸带调查中，采用经纬度多项式用于坐标系转换这种方法进行新（ED87一欧洲1987基准面）、旧（ED50一欧洲1950基准面）坐标系之间的转换：欧洲石油勘探组织（EPSG）对新、旧坐标系采用“双线性插值” 进行坐标转换[2]。
国内空间三维直角坐标转换中，一般采用布尔莎七参数模型。一般有7个转换参数，即3个平移参数，3个旋转参数和1个尺度参数。需要三个及已经公共点时，才能利用平差的方法求出七参数。
1.3研究的主要内容
本坐标转换程序可实现功能有:1、大地坐标与空间直角坐标的相互转换， 2、高斯投影正反算，3、二维坐标转换与四参数计算，4、三维坐标转换与七参数转换，5、同参考基准下的坐标换带计算，以及坐标数据的批量处理。
1.4程序设计思路方法
本程序名为万能坐标转换器。设计前期收集相关资料，参考一些成熟的坐标转换软件，确定程序应有的功能以及界面设计。运用VB编写程序时，查阅相关书籍获取理论知识以及转换模型。完成程序后将已知正确数据带入其中验证程序结果是否正确。若出现错误则检查每步代码，直到程序完美运行为止。
2.1地球椭球
地球椭球体又称地球椭圆体或地球扁球体，代表地球大小和形状的数学曲面，以长半径和扁率表示，因它十分迫近于椭球体，故通常以参考椭球体表示地球椭球体的形状和大小。通常所说地球的形状和大小，实际上就是以参考椭球体的半长径、半短径和扁率来表示。1975年国际大地测量与地球物理联合会推荐的数据为：半长径6378140米，半短径6356755米，扁率1∶298.257。在众多椭球体中，WGS-84椭球体被认为符合上述条件最好的椭球。
2.2基准
所谓基准是指为描述空间位置而定义的点、线和面。而大地测量基准是指用以描述地球形状的地球椭球参数，包含描述地球椭球几何特征的长短半轴和物理特征的有关参数、地球在空间的定位及定向以及描述这些位置所采用的单位长度的定义。不同的坐标系统会使用的基准也不同，根据参考椭球所选原点位置不同，可以分为地心坐标系和参心坐标系。
地心坐标系是以地球的质心为原点，有地心大地坐标系和地心空间直角坐标系两种表述方法。地心空间直角坐标系的定义为：以地球质心为原点，X轴指向格林尼治子午面与地球赤道的交点，Z轴指向北极，Y轴过原点垂直于平面 XOZ，构成右手空间直角坐标系。地心大地坐标系定义为：以地球的质心作为原点，以地球自转轴作为椭球的短轴，大地纬度B是过地面点的椭球法线与椭球赤道面之间的夹角，大地经度L为过地面点的椭球子午面与格林尼治子午面之间的夹角，大地高度H为地面点沿椭球法线到椭球面的最短距离。WGS-84坐标系，CGCS2000坐标系，GLONASS是采用PZ-90坐标，都是属于地心坐标系。
参心坐标系是选取一个参考椭球面作为基本的参考面，选一参考点作为大地测量的起算点，从而确定参考椭球在地球面的位置和方向。这时参考椭球的原点不会和地球质心重合，所以称为参心。北京54坐标系、西安80坐标系和新北京54坐标系，都是参心坐标系。它同样具有参心大地坐标系和参心空间直角坐标系两种表述方法，它们的定义与地心坐标系的定义相似。
2.3测量常用坐标系
2.3.1大地坐标系
空间一点的大地坐标用大地经度L、大地纬度B和大地高度H表示，地面上P地点的大地子午面NPS与起始大地子午面所构成的二面角L称P地点的大地经度， P地点对于椭球的法线与赤道面的夹角B称P地点的大地纬度。如图2-1所示
图2-1大地坐标系
P地点沿法线到椭球面的距离H称大地高，从椭球面起算，向外为正，向内为负[3]。
H = H正常 + ζ（高程异常）
H = H正 + N（大地水准面差距）
2.3.2空间直角坐标系
空间直角坐标系的坐标原点与参考椭球的中心重合，Z轴正向指向参考椭球的北极，X轴正向指向起始子午面与赤道的交点，Y轴按右手系与X轴呈90°夹角且位于赤道面上，某点在空间中的坐标可用该点在此空间坐标系的各个坐标轴上的投影来表示[4]。如图2-2所示：
图2-2空间直角坐标系
2.3.3平面坐标系
平面直角坐标系是利用投影，将空间坐标通过某种数学变换映射到平面上，这种变换称为投影变换[5]。在我国一般采用的是高斯一克吕格投影，是目前测量上广泛采用的正形投影，特点是没有角度变形，在不同点上的长度比随点位而异，但在同一点上各方向的长度比相同，也称为高斯投影[6]。
2.3.4地方独立坐标系
在我国平面坐标主要采用的是高斯投影，在该投影中，除中央子午线外，其它位置上的任何线段，投影后都会产生一定的长度变形，而且变形随离开中央子午线的距离增加而增加[7]。因此一般采用分带投影的办法，来限制长度变形，我国规定了采用3度带或6度带进行分带投影。在城市、工矿等工程测量中，如果直接在国家分带坐标系中建立控制网，会使地面长度投影的变形较大，当长度变形大于2．5 cm/km时，就难以满足工程上的需要[8]。另一些特殊的测量，比如大桥施工测量，水利水坝测量，滑坡变形监测等，采用国家坐标系精度不能满足工程要求，所以常常会建立适合本地区的地方独立坐标系[9]。
2.4我国常用坐标系
2.4.1 1954年北京坐标系
1954年北京坐标系，是采取先将我国一等锁与前苏联远东一等锁相联接，然后以连接处呼玛，吉拉林，东宁基线网扩大边端点的前苏联1942年普尔科沃坐标系的坐标为起算数据，平差我国东北及东部一等锁[10]。椭球参数：长半轴a=6378245m，短半轴b=6356863.0188m，扁率α=1/298.3，第一偏心率平方e2=0.00669342161454287，第二偏心率平方e’2=0.00673852540614652。
2.4.2 1980国家大地坐标系
1980国家大地坐标系采用地球椭球基本参数为1975年国际大地测量与地球物理联合会第十六届大会推荐的数据，大地原点设在我国中部的陕西省泾阳县永乐镇，位于西安市西北方向约60公里。椭球参数：长半轴a=6378140±5m，短半轴b=6356755.2882m，扁率α=1/298.257，第一偏心率平方e2 =0.00669438499959，第二偏心率平方e’2=0.00673950181947。
2.4.3 WGS-84世界大地坐标系
WGS-84坐标系是一种国际上采用的地心坐标系。坐标原点为地球质心，其地心空间直角坐标系的Z轴指向BIH （国际时间服务机构）1984.0定义的协议地球极（CTP）方向，X轴指向BIH 1984.0的零子午面和CTP赤道的交点，Y轴与Z轴、X轴垂直构成右手坐标系，称为1984年世界大地坐标系统。椭球参数：长半轴a=6378137，短半轴b=6356752.3142，扁率α=1/298.2572236，第一偏心率平方e2=0.00669437999013 ，第二偏心率平方e’2=0.006739496742227。
2.4.4 2000国家大地坐标系
CGCS2000是（中国）2000国家大地坐标系的缩写，该坐标系是通过中国GPS连续运行基准站、空间大地控制网以及天文大地网与空间地网联合平差建立的地心大地坐标系统[11]。Z轴指向BIH1984.0定义的协议极地方向（BIH国际时间局），X轴指向BIH1984.0定义的零子午面与协议赤道的交点，Y轴按右手坐标系确定。椭球参数：长半轴a=6378137，短半轴b=6356752.31414，扁率α=1/298.2572236，第一偏心率平方e2=0.00669437999013 ，第二偏心率平方e’2=0.006739496742227。
2.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。具体计算时，可先根据下式求出大地纬度B的初值：
因为 ,
带入B=arctan(
)
得:
令:
则式子可写成
然后利用该初值B代入公式右端tanB中，将等式左边的结果再次代入右端tanB，直到最后两次B值之差小与允许误差为止。当得到大地纬度B后，代入公式即可求出大地高H。
2.5.2高斯正反算转换模型
得到了点的大地坐标(L，B)，就可以将其转化为某投影带的高斯平面坐标，我们将椭球参数代入高斯投影正算公式得到更适用于电算的高斯坐标计算的实用公式：
其计算结果的精度可达0. 001m。
只要得到了高斯平面坐标(X, Y)后，便可通过高斯反算公式将其转换成大地坐标(B, L)，高斯投影反算公式为:
它们的计算精度,即平面坐标可达。
2.5.3坐标转换与参数计算转换模型
二维转换模型：
2个平移参数（原点不重合产生的）
1个旋转参数（坐标轴不平行产生的）
1个尺度参数（两坐标系间的尺度不一致产生的）
为某点在A空间直角坐标系中的坐标。
为某点在B空间直角坐标系中的坐标。
为某点从A空间直角坐标系转换到B空间直角坐标系的两个平移参数。β为从A空间直角坐标系转换到B空间直角坐标系中标系的两个平移参数。β为从A空间直角坐标系转换到B空间直角坐标系中一个旋转参数。m为从A空间直角坐标系转换到B空间直角坐标系中的一个尺度参数。
平面四参数求解步骤如下:
利用公共点计算坐标参数，但至少有两个公共点，当有i个公共点时，可利用最小二乘原理求解参数。
将B直角坐标系中的坐标视为观测值，设A直角坐标系下的坐标视为无误差，列误差方程为：
写成矩阵形式既：
由于各点的坐标可视为同精度独立观测值，因此P=I。
把各点坐标已知值带入上述误差方程，然后按下列公式求解出四参数：
三维转换模型：
3个平移参数（原点不重合产生的）
3个旋转参数（坐标轴不平行产生的）
1个尺度参数（两坐标系间的尺度不一致产生的）
一般为微小转角，模型可简化为：
为某点在A空间直角坐标系中的坐标。
为某点在B空间直角坐标系中的坐标。
为某点从A空间直角坐标系转换到B空间直角坐标系的三个平移参数。
为从A空间直角坐标系转换到B空间直角坐标系中三个旋转参数。m为从A空间直角坐标系转换到B空间直角坐标系中的一个尺度参数。
七参数求解步骤如下:
利用公共点计算坐标参数，但至少有三个公共点，当有i个公共点时，可利用最小二乘原理求解参数。
将B空间直角坐标系中的坐标视为观测值，设A空间直角坐标系下的坐标视为无误差，列误差方程为：
写成矩阵形式既：
由于各点的坐标可视为同精度独立观测值，因此P=I。
把各点坐标已知值带入上述误差方程，然后按下列公式求解出七参数：
2.5.4三、六度带带号与中央子午线计算模型
求带号、中央子午线公式为：
3°带:带号 中央子午线
L为当地经度
6°带:带号 中央子午线
L为当地经度
计划程序有哪些功能
程序界面设计
搜集坐标转换公式
运用VB编写程序
检查能否正常运行
以及数据正确性
否
是
完成程序
图3-1 编写程序流程图
图3-2程序主界面
3.1大地坐标系与空间直角坐标系的转换
主要代码如下：
椭球参数设置
If Combo1.Text = "WGS-84" Then
a1 = 6378137: b1 = 6356752.3142 'a1为长半轴，b1为短半轴
Text7.Text = a1: Text8.Text = b1
ElseIf Combo1.Text = "CGCS2000" Then
a1 = 6378137: b1 = 6356752.31414
Text7.Text = a1: Text8.Text = b1
ElseIf Combo1.Text = "北京54" Then
a1 = 6378245: b1 = 6356863.0188
Text7.Text = a1: Text8.Text = b1
ElseIf Combo1.Text = "西安80" Then
a1 = 6378140: b1 = 6356755.2882
Text7.Text = a1: Text8.Text = b1
ElseIf Combo1.Text = "自定义" Then
Text7.Enabled = True: Text8.Enabled = True
a1 = Text7.Text: b1 = Text8.Text
End If
大地坐标转空间直角坐标主要代码
B = Text1.Text: L = Text2.Text: H = Text3.Text
e12 = (a1 ^ 2 - b1 ^ 2) / a1 ^ 2
N = 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: Text5.Text = y: Text6.Text = z
空间直角坐标转大地坐标主要代码
x = Text4.Text: y = Text5.Text: z = Text6.Text
L = HuToDu(Atn(y / x)) '计算经度
If y / x < 0 Then
L = HuToDu(Atn(y / x) + 3.14159265358979)
End If
e12 = (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) / Sqr(x ^ 2 + y ^ 2)
k = 1 + (a1 ^ 2 - b1 ^ 2) / b1 ^ 2
TB1 = t + (p * TB) / Sqr(k + TB ^ 2) '用TanB作为未知量
Do Until Abs(TB1 - TB) < 0.000000000001 '迭代法
TB = TB1
TB1 = t + (p * TB) / Sqr(k + TB ^ 2)
Loop
B = Atn(TB)
N = a1 / Sqr(1 - e12 * (Sin(B) ^ 2)) '卯酉圈曲率半径
H = z / Sin(B) - N * (1 - e12) '大地高
Text1.Text = HuToDu(B): Text2.Text = L: Text3.Text = H
3.2高斯平面坐标与大地坐标的转换
高斯平面坐标(X, Y)与大地坐标(B ,L)的相互关系式分为两类：第一类称高斯投影正算公式，由(B,L)求(X,Y)；第二类称高斯投影反算公式，由(X,Y)求(B,L) [13]。由于电子计算机和各种可编程序电子计算器在测量上广泛应用，因此在这里给出适合电算的高斯正反算公式。
主要代码如下：
正算
Dim B, L, x, y As Double
Dim N, a1, b1, c, e12, e22, e1, e2, ll, l0, nn, t As Double
Dim XX, bb0, bb2, bb4, bb6, bb8 As Double
Dim yj As Double
Dim mf, nf, nnf, tf, bf, bf1 As Double
B = Text1.Text: L = Text2.Text
l0 = 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 = 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 ^ 4
bb2 = bb0 - 1
bb4 = (15 / 32) * e22 ^ 2 - (175 / 384) * e22 ^ 3 + (3675 / 8192) * e22 ^ 4
bb6 = -(35 / 96) * e22 ^ 3 + (735 / 2048) * e22 ^ 4
bb8 = (315 / 1024) * e22 ^ 4
ll = (Du(L) - l0) / 57.295779513
N = a1 / Sqr(1 - e12 * (Sin(DuToHu(B)) ^ 2)) '卯酉圈曲率半径
t = Sin(DuToHu(B)) / Cos(DuToHu(B))
nn = e2 * Cos(DuToHu(B))
XX = 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 + (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 + nn ^ 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.Text
y = y - yj
bf = x / (c * bb0) '迭代
Do
bf1 = ((x / c) - (bb2 * Cos(bf) + bb4 * (Cos(bf) ^ 3) + bb6 * (Cos(bf) ^ 5) + bb8 * (Cos(bf) ^ 7)) * Sin(bf)) / bb0
Loop Until Abs(bf - bf1) < 0.000000001
mf = (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 * (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 高斯投影换带计算
高斯投影虽然保证了角度没有变形，但是长度变形比较严重，为了限制高斯投影的长度变形，必须以中央子午线进行分带，把投影范围限制在中央子午线东、西侧一定的范围内；但是这样又使得统一的坐标系分割成各带的独立坐标系[14]。这样就会产生新的问题，需要邻带换算来解决。
主要代码如下：
l1 = Text3.Text : yy1 = Text4.Text '加常数和中央子午线
x = Text9.Text : y = Text10.Text
反算
y = y - yy1
e12 = (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(1 - e12) '极点子午线曲率半径
bb0 = 1 - (3 / 4) * e22 + (45 / 64) * e22 ^ 2 - (175 / 256) * e22 ^ 3 + (11025 / 16384) * e22 ^ 4
bb2 = bb0 - 1
bb4 = (15 / 32) * e22 ^ 2 - (175 / 384) * e22 ^ 3 + (3675 / 8192) * e22 ^ 4
bb6 = -(35 / 96) * e22 ^ 3 + (735 / 2048) * e22 ^ 4
bb8 = (315 / 1024) * e22 ^ 4
bf1 = x / (c * bb0) '迭代 单位弧度
Do
bf = bf1
bf1 = ((x / c) - (bb2 * Cos(bf) + bb4 * (Cos(bf) ^ 3) + bb6 * (Cos(bf) ^ 5) + bb8 * (Cos(bf) ^ 7)) * Sin(bf)) / bb0
Loop Until Abs(bf - bf1) < 0.000000001
mf = (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)) / (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 * (nnf ^ 2) * (tf ^ 2)) * (y ^ 5) / (120 * (nf ^ 5) * Cos(bf))
L = L + DuToHu(l1)
正算
l2 = 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 = 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 ^ 4
bbb2 = bbb0 - 1
bbb4 = (15 / 32) * ee22 ^ 2 - (175 / 384) * ee22 ^ 3 + (3675 / 8192) * ee22 ^ 4
bbb6 = -(35 / 96) * ee22 ^ 3 + (735 / 2048) * ee22 ^ 4
bbb8 = (315 / 1024) * ee22 ^ 4
L = HuToDu(L)
ll = (Du(L) - l2) / 57.295779513
XX = 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 * (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 + yy2
3.4 二维坐标转换与四参数计算
由于矩阵计算代码繁琐，由后台调用Excel函数完成计算。
其主要代码如下：
Dim x#, y#, x1#, y1#, dx#, dy#, jd#, k#
正算：
x1 = (Cos(DuToHu(jd)) * x - Sin(DuToHu(jd)) * y) * k + dx
y1 = (Sin(DuToHu(jd)) * x + Cos(DuToHu(jd)) * y) * k + dy
反算：
x = ((x1 - dx) * Cos(DuToHu(jd)) + (y1 - dy) * Sin(DuToHu(jd))) / k
y = ((x1 - dx) * (-Sin(DuToHu(jd))) + (y1 - dy) * Cos(DuToHu(jd))) / k
四参数计算：
Dim xls As Excel.Application
Dim book As Excel.Workbook
Dim sheet As Excel.Worksheet
Set xls = CreateObject("Excel.Application") '创建EXCEL对象
Set book = xls.Workbooks.Open(Text22.Text)'打开已经存在的工件簿文件
Set sheet = book.Worksheets(1)
sheet.Activate
a = sheet.UsedRange.Rows.Count
m = 1
For i = 2 To a
x = xls.Cells(i, 2).Value: y = xls.Cells(i, 3).Value
XX = xls.Cells(i, 4).Value: yy = xls.Cells(i, 5).Value
xls.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 = x
xls.Cells(m + 1, 6).Value = 0: xls.Cells(m + 1, 7).Value = 1: xls.Cells(m + 1, 8).Value = -x: xls.Cells(m + 1, 9).Value = y
m = m + 2
Next i
Range(Cells(1,11),Cells(4,2*(a+1)+10))=Application.WorksheetFunction.Transpose(Range(Cells(1, 6), Cells(2 * (a - 1), 9))) '参数B转置
Range(Cells(5,11),Cells(8,14))=Application.WorksheetFunction.MMult(Range(Cells(1, 11), Cells(4, 2 * (a - 1) + 10)), (Range(Cells(1, 6), Cells(2 * (a - 1), 9)))) 'B转置*B
Range(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(Cells(1, 11), Cells(4, 2 * (a - 1) + 10)), (Range(Cells(1, 10), Cells(2 * (a - 1), 10)))) 'B转置*L
Range(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 = True
xls.Quit
3.5 三维坐标转换与七参数计算
由于矩阵计算代码繁琐，由后台调用Excel函数完成计算。
其主要代码如下：
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 + dx
y1 = (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(jdx)) * Sin(DuToHu(jdy)) * Sin(DuToHu(jdz))) + z * Sin(DuToHu(jdx)) * Cos(DuToHu(jdy))) * k + dy
z1 = (x*(Sin(DuToHu(jdx))*Sin(DuToHu(jdz))+Cos(DuToHu(jdx))* Sin(DuToHu(jdy)) * Cos(DuToHu(jdz))) + y * (-Sin(DuToHu(jdx)) * Cos(DuToHu(jdz)) + Cos(DuToHu(jdx)) * Sin(DuToHu(jdy)) * Sin(DuToHu(jdz))) + z * Cos(DuToHu(jdx)) * Cos(DuToHu(jdy))) * k + dz
反算：
x = ((x1 - dx) * Cos(DuToHu(jdy)) * Cos(DuToHu(jdz)) + (y1 - dy) * (-Cos(DuToHu(jdy)) * Sin(DuToHu(jdz)) + Sin(DuToHu(jdx)) * Sin(DuToHu(jdy)) * Cos(DuToHu(jdz))) + (z1 - dz) * (Sin(DuToHu(jdx)) * Sin(DuToHu(jdz)) + Cos(DuToHu(jdx)) * Sin(DuToHu(jdy)) * Cos(DuToHu(jdz)))) / k
y = ((x1 - dx) * Cos(DuToHu(jdy)) * Sin(DuToHu(jdz)) + (y1 - dy) * (Cos(DuToHu(jdx)) * Cos(DuToHu(jdz)) + Sin(DuToHu(jdx)) * Sin(DuToHu(jdy)) * Sin(DuToHu(jdz))) + (z1 - dz) * (-Sin(DuToHu(jdx)) * Cos(DuToHu(jdz)) + Cos(DuToHu(jdx)) * Sin(DuToHu(jdy)) * Sin(DuToHu(jdz)))) / k
z = ((x1 - dx) * (-Sin(DuToHu(jdy))) + (y1 - dy) * Sin(DuToHu(jdx)) * Cos(DuToHu(jdy)) + (z1 - dz) * Cos(DuToHu(jdx)) * Cos(DuToHu(jdy))) / k
七参数计算：
m = 1
For i = 2 To a
x = xls.Cells(i, 2).Value: y = xls.Cells(i, 3).Value: z = xls.Cells(i,4).Value
XX = xls.Cells(i, 5).Value: yy = xls.Cells(i, 6).Value: zz = xls.Cells(i,7).Value
xls.Cells(m, 15).Value = XX - x: xls.Cells(m + 1, 15).Value = yy - y: xls.Cells(m + 2, 15).Value = zz - z '参数L
xls.Cells(m, 8).Value = 1: xls.Cells(m, 9).Value = 0: xls.Cells(m,10).Value = 0: xls.Cells(m, 11).Value = 0: xls.Cells(m, 12).Value =-z: xls.Cells(m, 13).Value = y: xls.Cells(m, 14).Value = x
xls.Cells(m + 1, 8).Value = 0: xls.Cells(m + 1, 9).Value = 1: xls.Cells(m + 1, 10).Value = 0: xls.Cells(m + 1, 11).Value = z: xls.Cells(m + 1, 12).Value = 0: xls.Cells(m + 1, 13).Value = -x: xls.Cells(m + 1, 14).Value = y
xls.Cells(m + 2, 8).Value = 0: xls.Cells(m + 2, 9).Value = 0: xls.Cells(m + 2, 10).Value = 1: xls.Cells(m + 2, 11).Value = -y: xls.Cells(m + 2, 12).Value = x: xls.Cells(m + 2, 13).Value = 0: xls.Cells(m + 2, 14).Value = z
m = m + 3
Next i
Range(Cells(1,16),Cells(7,3*(a-1)+15))= Application.WorksheetFunction.Transpose(Range(Cells(1, 8), Cells(3 * (a - 1), 14))) '参数B转置
Range(Cells(8,16),Cells(14,22))=Application.WorksheetFunction.MMult(Range(Cells(1, 16), Cells(7, 3 * (a - 1) + 15)), (Range(Cells(1, 8), Cells(3 * (a - 1), 14)))) 'B转置*B
Range(Cells(15,16),Cells(21,22)=Application.WorksheetFunction.MInverse(Range(Cells(8, 16), Cells(14, 22))) '矩阵求逆
Range(Cells(22,16),Cells(28,16))=Application.WorksheetFunction.MMult(Range(Cells(1, 16), Cells(7, 3 * (a - 1) + 15)), (Range(Cells(1, 15), Cells(3 * (a - 1), 15)))) 'B转置*L
Range(Cells(22,18),Cells(28,18))=Application.WorksheetFunction.MMult(Range(Cells(15, 16), Cells(21, 22)), (Range(Cells(22, 16), Cells(28, 16)))) '七参数
3.6 三、六度带带号与中央子午线计算
在采用分带的投影坐标系统中，我们最常用的是高斯-克吕格投影，它是由德国数学家、物理学家、天文学家高斯（Carl Friedrich Gauss，1777—1855）于十九世纪二十年代拟定，后经德国大地测量学家克吕格（Johannes Kruger，1857～1928）于1912年对投影公式加以补充，所以因此而得名。它是横轴墨卡托投影的一个变种，高斯-克吕格只是它通俗的名称，比较专业的名称叫做横轴等角切椭圆柱投影。我国各种大、中比例尺地形图采用了不同的高斯-克吕格投影带。其中大于1∶1万的地形图采用3°带；1∶2.5万至1∶50万的地形图采用6°带[15]。
其主要代码如下：
Dim dfm, i, d, f, m As Single
i = Val(Text1.Text)
d = Int(i)
f = Int((i - d) * 100)
m = Int((((i - d) * 100) - Int((i - d) * 100)) * 100)
dfm = d + f / 60 + m / 3600
Text2.Text = Int(dfm / 6) + 1
Text3.Text = (Int(dfm / 6) + 1) * 6 - 3 & "°"
Text4.Text = Int((dfm - 1.5) / 3) + 1
Text5.Text = (Int((dfm - 1.5) / 3) + 1) * 3 & "°"
一些主要函数代码如下
声明弧度转度.分秒函数：
Public Function HuToDu(ByVal hu As Double) As Double
Dim Du#, fen#, miao#
hu = hu * 180 / 3.1415926535898
Du = Fix(hu)
hu = (hu - Du) * 60
fen = Fix(hu)
hu = (hu - fen) * 60
miao = hu + 0.00001
miao = Format(miao, "0.00")
If miao = 60 Then
fen = fen + 1
miao = 0
If fen = 60 Then
Du = Du + 1
fen = 0
End If
End If
HuToDu = Du + fen / 100 + miao / 10000
End Function
声明度.分秒转弧度函数：
Public Function DuToHu(ByVal dfm As Double) As Double
Dim d#, f#, m#, angle#
d = Fix(dfm)
dfm = (dfm - d + 0.0000000001) * 100
f = Fix(dfm)
m = (dfm - f) * 100
angle = Format(d + f / 60 + m / 3600, "0.00000000000")
DuToHu = angle * 3.1415926535898 / 180
End Function
声明度.分秒转度函数：
Public Function Du(ByVal dfm As Double) As Double
Dim d#, f#, m#, angle#
d = Fix(dfm)
dfm = (dfm - d + 0.00000001) * 100
f = Fix(dfm)
m = (dfm - f) * 100
Du = Format(d + f / 60 + m / 3600, "0.000000000000")
End Function
声明秒转度.分秒函数：
Public Function MtoDFM(ByVal s As Double) As Double
Dim d#, f#, m#
d = s / 3600
f = Fix(d - Fix(d)) * 60
m = ((d - Fix(d)) * 60 - f) * 60
d = Fix(s / 3600)
MtoDFM = d + f / 100 + m / 10000
End Function
批量处理主要代码：
On Error GoTo hh
CD1.ShowOpen
Dim xls As Excel.Application
Dim book As Excel.Workbook
Dim sheet As Excel.Worksheet
Set xls = CreateObject("Excel.Application") '创建EXCEL对象
Set book = xls.Workbooks.Open(CD1.FileName) '打开已经存在的test.xls工件簿文件
Set sheet = book.Worksheets(1)
sheet.Activate
xls.Visible = False
a = sheet.UsedRange.Rows.Count '获取已经使用多少行
……
主要转换代码……
……
CD1.Filter = "(*.xlsx)|*.xlsx"
CD1.ShowSave
sheet.SaveAs (CD1.FileName)
MsgBox "保存位置:" & CD1.FileName, , "恭喜您保存成功"
xls.Quit
hh:
本章将我编写的软件计算结果与中海达HGO数据处理软件中坐标转换计算结果相对比，已验证程序的计算精度。
1.验证大地坐标与空间直角坐标相互转换
图4-1 BLH转换XYZ界面
图4-2 中海达软件验证界面
2.验证高斯正反算
图4-3 高斯正反算界面
图4-4 中海达软件验证界面
3.验证高斯换带计算
图4-5 换带计算界面
图4-6 中海达软件验证界面
4.验证二维坐标转换
图4-7 二维坐标转换界面
图4-8 中海达软件验证界面
5.验证三维坐标转换
图4-9 二维坐标转换界面
图4-10 中海达软件验证界面
6.验证利用公共点结算参数
公共点为大地测量学课本上的，经验证结果正确，再次不继续演示。
本坐标转换程序的设计主要是为了方便在施工过程中控制点在不同坐标系之间转换的问题。在实际应用中工程必然会涉及到测量坐标转换问题，所以开发这样的程序是十分有必要的。由于本人的知识水平和计算机编程语言掌握的成度有限，在编写程序时查阅了大量参考书籍，并且参考了一些成熟软件的界面设计。因为初次编写软件，程序只能达到最基本的要求，一定还存在许多不足之处，恳请批评指正。
我要特别感谢杜继亮老师在我大学的最后学习阶段——毕业设计阶段对我设计的细心指导，从选题和大方向上的指导到修改、定稿都有老师足迹，由衷的感谢老师为我完成这篇论文提供了巨大的帮助。老师认真的工作态度诚信宽厚的为人处世态度都给我留下了难以磨灭的印象也为我今后的工作树立了优秀的榜样。
同时感谢所有任课老师和所有同学在这四年来给自己的指导和帮助，是他们教会了我专业知识教会了我如何学习教会了我如何做人。感谢同学们在四年的大学中对我学习与生活中的帮助和关心。在这次毕业设计的过程中把坐标转换方面的理论知识更好的理解一遍，自学能力在这里得提升，VB编程能力也有了很大的提高，对以后工作中遇到的问题也可以用程序来解决。
掩卷而思，唯有感激，唯有自省，唯有不懈地追寻！
[1]蔡习文.大地测量坐标转换算法研究及系统开发[J].矿山测量，2009，(3) :73 -75
[2]陈航.不同坐标系的转换及其程序设计[D].江苏:江苏师范大学，2009:1-2
[3] 张宁红．坐标转换在雷达实时显示软件中的应用[J].现代雷达，2002，24(5):30 -32
[4]张国舫，赵春年.一种气象雷达车方位角仰角数据精度修正的设计[J].第五届长三角科技论坛——长三角气象科技创新论坛论文集，2008，()： -
[5]潘开志.手持式GPS接收机应用中坐标系统的转换[J].重庆林业科技，2010，(1):29 -32
[6]王清华，鄂栋臣.南极地区常用地图投影及其应用[J].极地研究，2002，14(3)：226 -233
[7]徐生望.基准转换与城市平面控制网若干问题研究[D].中南大学，2008
[8]张瑛.城乡规划中的控制基准建设[J].城市勘测，2007，(3):3 -5
[9]赵继宏，余荣宣.矿区测量中选择适宜坐标系统的方法研讨[J].西部探矿工程，2012，(8) :97 -99
[10]许家琨.常用大地坐标系的分析比较[J].海洋测绘，2005，25(6) :71 -74
[11]史俊莉，牛鹏涛. CGCS2000坐标系的推广与应用[J].价值工程，2014，33(26)： 229-230
[12]李海燕，王晓峰.浅谈坐标系的建立和转换[J].新课程：教研版，2010，(3)：93 -94
[13]刘大地. p2p网络拓扑发现[D].北京工业大学，2007
[14]王虎，尹水清，胡琨.大连地方投影带选取问题的探讨[J].矿山测量，2007，(3):51 -52
[15]王永成，杨保山，侯来福.浅谈地形图分幅与编号的计算方法[J].西部探矿工程，2010，22(6):136 -138
本文来源：https://www.2haoxitong.net/k/doc/fd23c79ec4da50e2524de518964bcf84b9d52da7.html
文档为doc格式