VB实现高斯正反算

1、VB实现高斯正反算专业班级：测绘0729班学生姓名：张智宁学生学号：0700001955指导教师：欧阳老师完成时间：2010年6月5日1前言在日益发展的社会中，能很好很快的完成各种任务才能取得更多的发展机会。用VB实现高斯正反算，就能很快实现高斯坐标系和大地坐标系的转换，为我们的计算和测量工作提供高效的服务。此次设计要求我们利用大学所学的VB编程和控制测量中的高斯平面坐标系为基础，用VB实现高斯正反算，做到能够完整写出VB实现高斯正反算的代码，而且还能生成运行文件进行高斯正反算。同时，也让我们体会到编程在测量中应用之广，在以后的工作中应该博学多思 。2控制测量的基本理论1.大地水准面大地水准面

2、高是大地水准面至地球椭球面的垂直距离。与静止海平面重合的重力等位面称为大地水准面。大地水准面是一个与地球内部密度分布有关的不规则曲面。实际观测结果表明，大地水准面非常接近一个旋转椭球面。选择适当的旋转椭球，使得大地水准面相对此椭球面的起伏的平方在椭球面上的积分最小，这个旋转椭球称为参考椭球。大地水准面与参考椭球面的垂直距离称为“大地水准面高”。在1980大地参考系中，参考椭球的赤道半径取为6 378.137米，扁率的倒数取为286 26，此时大地水准面的高度不超过110米。地面点的正高与该点的大地水准面高的和称为该点的大地高，大地高是地面点的大地坐标系中的一个坐标。大地水准面高在地面上的分布与

3、地球内部特别是地壳上地幔密度水平向分布不均匀有关，因此大地水准面高也是研究地球内部结构和地球内部动力学过程的一种。图1-13图1-14水准面静止的水面称为水准面，水准面是受地球重力影响而形成的，是一个处处与重力方向垂直的连续曲面，并且是一个重力场的等位面。设想一个静止的海水面扩展到陆地部分。这样，地球的表面就形成了一个较地球自然表面规则而光滑的曲面，这个曲面被称为水准面。图1-25高斯投影高斯-克吕格(Gauss-Kruger)投影简称“高斯投影”，又名等角横切椭圆柱投影”，地球椭球面和平面间正形投影的一种。德国数学家、物理学家、天文学家高斯（Carl Friedrichauss，1777一

4、1855）于十九世纪二十年代拟定，后经德国大地测量学家克吕格（Johannes Kruger，18571928）于 1912年对投影公式加以补充，故名。该投影按照投影带中央子午线投影为直线且长度不变和赤道投影为直线的条件，确定函数的形式，从而得到高斯一克吕格投影公式。投影后，除中央子午线和赤道为直线外， 其他子午线均为对称于中央子午线的曲线。设想用一个椭圆柱横切于椭球面上投影带的中央子午线，按上述投影条件，将中央子午线两侧一定经差范围内的椭球面正形投影于椭圆柱面。将椭圆柱面沿过南北极的母线剪开展平，即为高斯投影平面。取中央子午线与赤道交点的投影为原点，中央子午线的投影为纵坐标x轴，赤道的投影为

5、横坐标y轴，构成高斯克吕格平面直角坐标系。高斯-克吕格投影在长度和面积上变形很小，中央经线无变形，自中央经线向投影带边缘，变形逐渐增加，变形最大之处在投影带内赤道的两端。由于其投影精度高，变形小，而且计算简便（各投影带坐标一致，只要算出一个带的数据，其他各带都能应用），因此在大比例尺地形图中应用，可以满足军事上各种需要，能在图上进行精确的量测计算。高斯投影：它是一种横轴等角切圆柱投影。它把地球视为球体,假想一个平面卷成一个横圆柱面并把它套在球体外面,使横轴圆柱的轴心通过球的中心,球面上一根子午线与横轴圆柱面相切。这样,该子午线在圆柱面上的投影为一直线,赤道面与圆柱面的交线是一条与该子午线投影垂

6、直的直线。将横圆柱面展开成平面,由这两条正交直线就构成高斯-克吕格平面直角坐标系。为减少投影变形,高斯-克吕格投影分为3o带和6o带投影。 高斯克吕格投影是设想用一个椭圆柱横套在地球椭球的外面,并与设定的中央经线相切。 6高斯克吕格投影分带规定:该投影是国家基本比例尺地形图的数学基础,为控制变形,采用分带投影的方法,在比例尺 1:2.5万-1:50万图上采用6分带,对比例尺为 1:1万及大于1:1万的图采用3分带。 6分带法:从格林威治零度经线起,每6分为一个投影带,全球共分为60个投影带,东半球从东经0-6为第一带,中央经线为3,依此类推,投影带号为1-30。其投影代号n和中央经线经度L0的

7、计算公式为:L0=(6n-3);西半球投影带从180回算到0,编号为31-60,投影代号n和中央经线经度L0的计算公式为L0=360-(6n-3)。 3分带法:从东经130起,每3为一带,将全球划分为120个投影带,东经130-430,.17830-西经17830,.130-东经130。东半球有60个投影带,编号1-60,各带中央经线计算公式:L0=3n ,中央经线为3、6.180。西半球有60个投影带,编号1-60,各带中央经线计算公式:L0=360-3n ,中央经线为西经177、.3、0。我国规定将各带纵坐标轴西移500公里,即将所有y值加上500公里,坐标值前再加各带带号以18带为例,原

8、坐标值为y=243353.5m,西移后为y=743353.5,加带号通用坐标为y=18743353.5 。7VB简介vb是Visual Basic的简称，是由美国微软公司于1991年开发的一种可视化的、面向对象和采用事件驱动方式的结构化高级程序设计语言，可用于开发 Windows 环境下的各类应用程序。它简单易学、效率高，且功能强大可以与 Windows 专业开发工具SDK相媲美。在Visual Basic环境下，利用事件驱动的编程机制、新颖易用的可视化设计工具，使用Windows内部的广泛应用程序接口（API）函数，动态链接库（DLL）、对象的链接与嵌入（OLE）、开放式数据连接（ODBC）

9、等技术，可以高效、快速地开发Windows环境下功能强大、图形界面丰富的应用软件系统。8三 高斯正反算公式1，高斯投影正算已知椭球面上某点的大地坐标（L,B），该点在高斯投影平面上的直角坐标（x,y）， 即（L,B）=(X,Y)的坐标转换。 1.1 投影必须满足的条件（1）中央子午线投影后为直线。（2）中央子午线投影后长度不变。（3）投影具有正形性质，即正投影条件。1.2 正投影过程 在椭球上有对称于中央子午线的两点P1和P2，它们的大地坐标分别为（L,B）及（l,b），其中l为椭球面上p点的经度与中央子午线的经度L0之差：l=L-L0，P点在中央子午线之东，l为正，在西则为负，则投影后的平面

10、坐标一定为P1(X,Y)和P2(X,Y)。9高斯投影正算计算公式 2.高斯投影坐标反算公式2.1高斯投影反算 已知某点的高斯投影面上直角坐标（x,y），求该点在椭球面上的大地坐标（L,B），即（x,y）=(L,B)的坐标变换。2.2 投影变换必须满足的条件a. x坐标轴投影成中央子午线，是投影的对称轴。b. x轴上的长度投影保持不变。c.投影具有正形性质，即正形投影条件。2.3 根据x计算纵坐标在椭球面上的投影的底点纬度，接着按计算（-B）及经差，最后得到B= L=10高斯投影反算公式当要求转换精度至时 ，可以简化为下式：11VB设计界面以上就是VB实现高斯正反算得操作界面。12程序使用方法

11、这个程序使用很简单，首先确认中央子午线的，点击中央子午线下方的小黑三角形选择大地坐标或或高斯投影平面坐标所在的中央子午线的经度。在计算之前还要选择所依据的椭球面，椭球不一样，计算的参数不一样。若选择“北京54”，则点击“北京54”前面的小白圆圈，显黑点就说明已经选中，例如上面界面所显示的那样。在选定中央子午线和椭球后，要大地坐标转高斯投影平面坐标，则在大地坐标下输入B,L，按上方的转换按钮，则可以得到高斯投影平面坐标X,Y。反之，输入高斯投影平面坐标X,Y，按下方的转换按钮，则可以得到大地坐标B,L。13程序代码 这个程序使用很简单，首先确认中央子午线的，点击中央子午线下方的小黑三角形选择大地

12、坐标或或高斯投影平面坐标所在的中央子午线的经度。在计算之前还要选择所依据的椭球面，椭球不一样，计算的参数不一样。若选择“北京54”，则点击“北京54”前面的小白圆圈，显黑点就说明已经选中，例如上面界面所显示的那样。在选定中央子午线和椭球后，要大地坐标转高斯投影平面坐标，则在大地坐标下输入B,L，按上方的转换按钮，则可以得到高斯投影平面坐标X,Y。反之，输入高斯投影平面坐标X,Y，按下方的转换按钮，则可以得到大地坐标B,L。六程序代码 Public Const PI = 3.149Public a As Double, e2 As Double, e22 As DoublePublic B As

13、 Double, L As Double, L0 As DoublePublic Bf As DoublePublic x0 As Double, y0 As DoublePrivate ll As DoublePublic Number As IntegerPublic AA As Double, BB As Double, CC As Double, DD As Double, EE As Double, FF As Double, GG As DoublePublic Sub Gauss(mL0 As Double, mB As Double, mL As Double) Dim T A

14、s Double, Yita As Double mL0 = mGetrad(mL0): mB = mGetrad(mB): mL = mGetrad(mL) Yita = e22 * (Cos(mB) 2): T = Tan(mB) ll = mL - mL0 x0 = mGetX(0, mB) + mGetN(mB) * T * (Cos(mB) 2 * ll 2 * (0.5 + (5 - T 2 + 9 * Yita + 4 * Yita 2) * (Cos(mB) 2 * ll 2 / 24 + (61 - 58 * T 2 + T 4) * (Cos(mB) 4 * ll 4 /

15、720) y0 = mGetN(mB) * Cos(mB) * ll * (1 + (1 - T 2 + Yita) * (Cos(mB) 2 * ll 2 / 6 + (5 - 18 * T 2 + T 4 + 14 * Yita - 58 * T 2 * Yita) * (Cos(mB) 4 * ll 4 / 120)End SubPublic Sub mGauss(mL0 As Double, mB As Double, mY As Double)14Dim T As Double, Yita As Double Dim N As Double, M As Double Yita = e

16、22 * (Cos(mB) 2): T = Tan(mB) L0 = mGetrad(mL0) N = mGetN(mB) M = mGetM(mB) ll = mY / Cos(mB) / N * (1 - (1 + 2 * T 2 + Yita) * (mY / N) 2 / 6 + (5 + 28 * T 2 + 24 * T 4 + 6 * Yita + 8 * Yita * T 2) * (mY / N) 4 / 120) B = mB - T * mY * mY / N / M / 2 * (1 - (5 + 3 * T 2 + Yita - 9 * Yita * T 2) * (

17、mY / N) 2 / 12 + (61 + 90 * T 2 + 45 * T 4) * (mY / N) 4 / 360) L = L0 + llEnd SubPublic Sub mGetBf(mB As Double) If Abs(x0 - mGetX(0, mB) 0.001 Then Bf = mB Exit Sub Else mB = mB + (x0 - mGetX(0, mB) / mGetSS(mB) Number = Number + 1 Call mGetBf(mB) End IfEnd SubPrivate Function mGetX(B1 As Double,

18、B2 As Double)15 mGetX = mGetM(B1 + B2) / 2) * (B2 - B1) * (1 + e22 * Cos(B1 + B2) * (B1 - B2) 2 / 8)End FunctionPrivate Function mGetSS(mB As Double) Dim R As Double, S As Double, T As Double R = a * (1 - e2) S = AA - BB * Cos(2 * mB) - CC * (Sin(mB) 2 * (Cos(2 * mB) + 2 * (Cos(mB) 2) T = -DD * (Sin

19、(mB) 4 * (Cos(2 * mB) + 4 * (Cos(mB) 2) - EE * (Sin(mB) 6 * (Cos(2 * mB) + 6 * (Cos(mB) 2) - FF * (Sin(mB) 8 * (Cos(2 * mB) + 8 * (Cos(mB) 2) mGetSS = R * (S + T) mGetSS = a * (1 - e2) * (AA - BB * Cos(2 * mB) - CC * (Sin(mB) 2 * (Cos(2 * mB) + 2 * (Cos(mB) 2) - DD * (Sin(mB) 4 * (Cos(2 * mB) + 4 *

20、(Cos(mB) 2) - EE * (Sin(mB) 6 * (Cos(2 * mB) + 6 * (Cos(mB) 2) - FF * (Sin(mB) 8 * (Cos(2 * mB) + 8 * (Cos(mB) 2)End FunctionPrivate Function mGetrad(S As Double) S = mTransfer2(S) * PI / 180 mGetrad = SEnd FunctionPublic Function mTransfer1(S As Double) 将以度为单位的角度转换 16Dim du As Integer, fe As Intege

21、r, mi As Single du = Int(S) fe = Int(S - du) * 60) mi = (S - du - fe / 60) * 3600 mTransfer1 = du + fe / 100 + mi / 10000End FunctionPublic Function mTransfer2(R As Double) 将小数度分秒格式的角度转换为以度为单位的格式 Dim du As Integer, fe As Integer, mi As Single If R 0 Then R = Abs(R) du = Int(R) fe = Int(R - du) * 100

22、) mi = (R * 100 - Int(R * 100) * 100 mTransfer2 = -(du + fe / 60 + mi / 3600) Else du = Int(R) fe = Int(R - du) * 100) mi = (R * 100 - Int(R * 100) * 100 mTransfer2 = du + fe / 60 + mi / 3600 End If17End FunctionPublic Function mTransfer3(R As Double) 检验出单位秒是60时，进行进位 Dim du As Integer, fe As Integer, mi As Single If R 0 Then R = Abs(R) du = Int(R) fe = Int(R - du) * 100) mi = (R * 100 - Int(R * 100) * 100 If mi = 60 Then mi = 0: fe = fe + 1 If fe = 60 Then fe = 0