测量专业基本计算与编程

第二部分 基本计算与编程 3 一 编程语言简介及特点 3 1.1 Visual Basic 概述 3 1.2 VB的基本概念 3 1.3 数据类型 4 1.4 标准控件简介 4 1.5 过程、函数和方法 6 1.6 应用程序的设计 9 1.7 简单的编程实例 11

二 测量中几种常用的计算 12 2.1 角度化为弧度 12 2.2 弧度化为角度计算 13 2.3 坐标方位的计算 13

三 空间大地坐标与空间直角坐标之间的换算 3.1 由空间大地坐标计算空间直角坐标 14 3.2 空间直角坐标计算空间大地坐标 14 3.3 计算程序与范例 15 四 大地主题问题正反计算 17 4.1 高斯平均引数大地问题解算 17 4.2 白塞尔大地主题解算 24

4.3 嵌套系数法解算任意距离的大地主题 32 4.4 大地主题解算程序功能介绍 39 五 子午线弧长正反解算 39 5.1 子午线弧长正算 39 5.2 子午线弧长反算 41

六 高斯投影正反算和邻带换算 42 6.1 高斯投影正算 42 6.2 高斯投影反算 42 6.3 高斯投影邻带换算 43

6.4 高斯投影计算程序及其功能 43 七 平面直角坐标系换算 47 7.1 直接参数法 47

7.2 相似变换（赫尔墨特法） 48 7.3 正形变换法 48 7.4 多项式逼近法 49

7.5 算例与程序功能介绍 50

八 测量仪器与计算机间的数据通讯 63 8.1 MSComm通讯控件及其属性简介 63

8.2 使用MSComm控件设计通讯程序的步骤 64 8.3 数据传输的设置 65

8.4 计算机接收数据的通讯程序 65 九 平面控制网数据处理 67 9．1 平面控制网的概算 67

14 9．2 误差方程式与法方程的组成 77 9．3 法方程的求逆与平差计算 82 9．3 精度评定 88 9．4 图形的绘制 92

十 高程控制网数据处理 97 10．1高程控制网的概算 97 10．2 误差方程的组成 97 10．3 法方程式组成与计算 97 10．4 精度评定 97

第二部分 基本计算与编程

《大地测量学基础》教程包含着许多测量基本计算问题，如常见坐标系的计算，坐标系之间的转换计算，椭球大地计算、高斯投影计算、平面网与高程网平差数据处理等。对于上述相关计算必须要求测绘专业学生加以掌握，其目的是通过计算加深学生对相关基本理论与方法的理解与掌握，同时有益于学生能力的培养。

采用计算机编程进行数据处理计算，其计算语言较多，比如Basic语言、C语言、Tortran、Pascal、Delphi等，不同语言各自具有不同的特点，为了方便学生对语言的初步学习，这里对Visual Basic语言作一简单介绍，供学生参考。

一 编程语言简介及特点

1.1 Visual Basic 概述

Visual Basic（以下简称VB）的前身是QBASIC，语言基础是BASIC。自从微软推出VB后，VB便成为程序开发人员的首选工具。它是一套完全独立的 Windows 开发系统，是可视化、面向对象、采用事件驱动方式的结构化高级程序设计语言，利用其事件驱动的编程机制、新颖易用的可视化工具，并使用 Windows 内部应用程序接口（API）函数，采用动态连接库（DLL）动态数据交换（DDE）、对象的链接与嵌入（OLE）以及开放式数据库访问（ODBC） 等技术，可以高效、快速地建立 Windows环境下功能强大、图形界面丰富的应用软件系统。据统计，仅在数据库系统开发领域，VB就占了90％的份额。VB是基于对象的可视化程序开发工具，其优点在于快捷简易地建立Windows应用程序。按使用人员来分，VB有以下三个版本：

1．标准版 针对一般程序人员，适合普通应用软件的开发。

2．专业版 针对专业程序开发人员，在标准版的基础上提供了对数据库和Internet的支持。 3．企业版 适合于企业设计应用软件的程序开发。

对于一般的设计者而言，只要充分发挥自身的想象力，任何人均可在在较短的时间内，利用VB语言开发各种自己的实用程序。可视化编程的一个突出特点是具备集成开发环境，在相应的开发平台上集成了编辑器，编译连接工具，控制器箱辅助工具。在VB集成开发环境下包括一些主要性能：工具箱，工具栏，工具管理器窗口，属性窗口，窗体设计器，代码编辑窗口等，同时集成环境的设置也非常灵活，开发人员可以按照自己的编程习惯来加以配置。 1.2 VB的基本概念

对象，是指可以当作一个单元的代码和数据的组合，它可以是程序中的窗体或控件，也可以是整个程序。对象有自身的状态与方法。

属性，是指对象具有的性质，表示对象的状态。对象的属性设置，可以在程序设计时在属性

窗口中实现，也可以使用代码设置属性的值，其语法为： 对象名_属性＝新值

事件，是指发生在对象上的事情。Windows应用程序属于\事件驱动\模式。对象对事件的反应又称为\事件过程\。 事件过程的语法： Sub 对象名_事件（） 处理事件的代码 End Sub

事件驱动，只有当事件发生时，程序才会运行。在没有事件的时候，整个程序处于停滞状态。VB程序中流动的不是数据而是事件，如果说属性决定了对象的外观，方法决定了对象的行为，那么事件就决定了对象之间联系的手段。

方法，对象本身包含的函数和过程。方法决定了对象可以进行的动作，方法的内容（代码）是不可见的，当我们需要使用某个对象的方法是，只需要使用其格式：对象名.方法 如清除窗体Form1上的内容： Form1.cls.

以坐标（1920，1300）为圆心，以800为半径画园，其方法为： Form1.Circle(1920,1300), 800

过程，指事件发生时要执行的代码。

面向对象编程， 以对象为核心，支持对象的封装机制，多态机制和继承机制（VB不能真正支持继承机制，故从严格意义上讲VB不是真正意义上的面向对象编程）。 1.3 数据类型

VB中有丰富的数据类型，这里仅作简单介绍。 1．数字型：包括整型（Integer）、长整型（Long）、单精度（Single）、双精度（Double）、货币型（Currency）等。例如：Dim X As Integer 2. 字符型(String):字符型变量可以存储可变的字符串。

3. 布尔型(Boolean)：一个变量包含简单的Yes/No、Ture/Fals信息，则可以定义为布尔型变量。例如：Dim temp As Boolean.

4. 日期型(Date):专门用来表示时间的数值类型，可以有多种表达方式。

5．对象型(Object):对象变量存储的是对象的地址信息，它本身并不是一个变量，但定义为Object类型的变量可以通过赋值语句指向程序能识别的任何对象。例如： Dim Mydb As Object

Set Mydb=OpenDatabase(\

这样对Mydb进行访问时，实际上就是对Access数据库tempDB.mdb进行访问。

6: Variant类型：Variant变量类型可以存储所有的数据类型，VB会自动执行相应的转换。但Variant变量类型会耗用较多的资源，所以不提倡采用。

1.4 标准控件简介

VB对控件有三种广义的分类：

1、内部控件：内部控件就是在工具箱中默认出现的控件，如CommandButton、Frame控件等。内部控件总是出现在在工具箱中，不像ActiveX控件和可插入对象那样可以添加到工具箱中，或从工具箱中加以删除。 2、ActiveX控件：ActiveX 控件是扩展名为 .ocx 的独立文件，其中包括各种版本 Visual Basic提供的控件（如 DBGrid、DBCombo、DBList 控件等），另外还有许多第三方提供的 ActiveX 控件。

3、可插入对象：如一个 Excel、Word 工作表对象。因这些对象能添加到工具箱中，故可以把他们当作控件使用。其中一些对象还支持自动化（OLE 自动化），使用这些控件就可在 Visual Basic应用程序中编程控制另一应用程序的对象。 这里仅对其中一些常用控件略加介绍 CommandButton 按钮（命令按钮）：按钮是用户与应用程序交互的最简便的方法。用户点击按钮，就会调用 Click 事件过程。将代码写入Click 事件过程，执行想要执行的动作。 Label 控件（标签）：用来在窗体中显示文字，一本对没有Caption属性的控件（如TextBox）起标识作用。文本是只读的，用户不能修改。 TextBox 控件（文本框）：文本框控件可用来显示与输入文本。如果TextBox控件的Locked属性设为True, 则用户不能直接修改文本框的文本内容。如果TextBox控件的MultiLine属性设为True, 则可显示多行文本。

CheckBox 控件（复选框）:复选框像一个开关，表明一个特定的状态是选定（on）还是清除（off）。在应用程序中使用 CheckBox 为用户提供了 Yes/No 或 True/False 的选择。CheckBox 是独立工作的，用户可以同时选择任意多个复选框。CheckBox 的 Value 属性有三个值（0是没有选中即默认值，1为选中，2 为灰色状态）。 OptionButton 控件（选项按钮）：选项按钮表示给用户一组或更多的选择。它不同于 CheckBox控件，选择一个选项按钮就会清除该组中的其他按钮，一组中只有一个选项处于选项状态。OptionButton 的 Value 属性有两个选项 True/False。 Frame 控件（框架）：Frame 控件用来创建选项按钮组。如果将一组 OptionButton 选项按钮放在一个窗体上，那么所有选项按钮将构成一组。如果相创建其他选项按钮组，就必须将其中一组选项按钮放在 Frame 控件中。

ListBox 控件（列表框）和ComboBox 控件（组合框）：列表框与组合框是在有限的空间为用户提供大量选项的有效方法。ComboBox 兼有 TextBox 和 ListBox 两者的功能，该控件允许用户通过键入文本或选择列表中的项目来进行选择。 HScrollBar 和 VScrollBar 控件（垂直与水平滚动条）：滚动条相当于一个模糊的输入装置，当用户不需要精确设置数据只要大概范围的时候，使用滚动条控件进行输入是一种合适的选择。在早期VB版本中通常把滚动条作为输出设备，但目前建议采用滑板控件。 CommonDialog 控件（通用对话框）：通用对话框控件提供了一组标准的操作对话框，可以进行如文件的打开与保存，设置打印选项，选择颜色和字体操作等。通过运行 Windows 的帮助引擎控件还能显示帮助。

DBGrid 控件(数据绑定网格）: DBGrid (Data Bound Grid) 控件提供了一般Grid 控件显示信息，但与 Grid 控件不同，一旦绑定了某个数据表，不仅可以将数据显示出来，而且还可以修改每个字段的内容。

MSFLexGrid 控件：MSFLexGrid 控件显示和操作表格数据。对包含字符串和图片的表格提供了灵活的排序、插入数据和格式编排功能。当与Data控件绑定时，只显示只读数据。可将 MSFLexGrid控件与 TextBox 控件组合，使其控件具有单元编辑功能。 DriveListBox 控件（驱动器列表框）：驱动器列表框是下拉式列表框。在默认式在用户系统上显示当前驱动器。当该控件获得焦点时，用户可以输入任何有效的驱动器标识符，或者从下拉列表中选定新驱动器。 DirListBox 控件（目录列表框）：目录列表框可以显示当前驱动器的目录结构。如果把驱动器列表框的Dirve属性赋予目录列表框的Path 属性，则可以显示驱动器列表中选定的驱动器的目录：

Dir1.Path=Dirve1.Dirve FileListBox 控件（文件列表框）：文件列表框在运行时显示由Path属性指定的包含在

目录中的文件。例如可用下列语句显示当前目录中的所有文件： File.Path=Dir1.Path

Image 控件：Image 控件用来显示图形。它可显示如位图、图标、图元文件、增强型图元文件、JPEG 或GIF格式的图形文件。

Line 控件：直线控件用来在窗体、框架或图片框中创建简单的线段。

PictrueBox 控件：PictrueBox 控件用来显示图形。它可显示如位图、图标、图元文件、增强型图元文件、JPEG 或GIF格式的图形文件。但PictrueBox 控件包含了Image 控件不具有的功能，如作为其它控件的容器的功能和支持图形方法的功能。

Shape 控件：可用形状控件在窗体、框架或图片框中创建如矩形、正方形、椭圆形、园形、园角矩形或园角正方形形状。

Timer 控件：Timer 控件响应时间的流逝。它独立于用户，编程用来在一定的时间间隔执行某种操作。对于其它的后台处理该控件非常有用。 除上述标准控件以外，这里仅介绍三种 ActiveX 控件。

MSComm 控件：MSComm 控件为你的应用程序提供了串口通讯功能，允许通过串口发送和接收数据。每个 MSComm 控件都与一个串口对应，如果在应用程序中需要访问多个串口，就必需有多个MSComm 控件。

TabStrip 控件: 一个 TabStrip 控件就象一组文件夹上的标签，使用 TabStrip 控件用户可以在应用程序中在窗口或对话框中的同一区域定义多个数据页面。该控件包含一个或多个 Tab 对象，可以通过 TabStrip 控件的属性对话框来增加或删除Tab对象。

Toolbar 控件:一个Toolbar 控件包含了一个 Button 对象的集合，用来创建与应用程序关联的工具条。

对于其它控件的应用请参考有关VB专业书籍，这里就不一一加以介绍。

1.5 过程、函数和方法

Visual Basic程序由若干子程序构成，这些子程序称为过程、函数和方法。他们都在代码窗口中设计。

1、 过程（Procedure）

完成某种特定功能的一组程序代码称为过程。在Visual Basic程序中用关键字 Sub 和 End Sub表示过程的开始和结束。VB中共有两种类型的过程。 事件过程（Event Procedure） 当用户在窗体上设计图形界面时，针对每个对象均有多个事件与其相关连，每个对象与每个事件都可以构成一个事件驱动程序，也就是说当用户或系统在某一对象上触发某种事件时，就会引发去执行相应的事件驱动程序，完成特定的功能。事件过程是依附于每个对象上的，由特定事件引发的程序，是VB程序的主体。VB在运行时会自动通过事件过程名称来识别执行哪个事件驱动程序。 Sub 对象名_事件名称（）

处理事件的代码 End Sub

通用过程（General Procedure）

当多个事件过程都需要完成某种公共的功能，如完成一些公共的数据计算，或对某些变量进行共同的操作，那么用户可自己建立通用过程，编写公共代码模块，供其他程序调用。通用模块的声明如下：

Sub 过程名称（参数1，参数2，…） 程序语句代码

End Sub

2、 函数（Funtion）

Visual Basic中包含两大类函数；一类是VB本身提供的已被封装好的通用函数，它不需要用户去创建和声明及编程，只需要直接调用；另一类是用户自定义函数。 ? 常用函数：包括数学函数、字符串函数、日期函数、类型转换函数。 ① 数学函数

函数表达式 函数功能说明 Abs( 表达式) 绝对值函数。

Int(表达式) 取整函数, 将表达式的值转换为小于或等于该表达式的最大整数。 Fix(表达式) 取整函数, 将表达式的值的小数部分截去，直接取整数。 Atn(表达式) 反正切函数。 Sin(表达式) 正弦函数。 Cos(表达式) 余弦函数。 Tan(表达式) 正切函数。 Sqr(表达式) 平方根函数。

Exp(表达式) 以e为底的指数函数。

Sgn(表达式) 符号函数，表达式值为正数、零、负数，函数值分别为1、0、－1。 Rnd(表达式) 随机数生成函数，返回0与1之间的单精度的随机数。

②字符串函数

函数表达式 函数功能说明 函数表达式 功能

Len(字符串) 字符串长度函数 Str(表达式) 字符串函数 Instr

查找字符函数

Lcase$(字符串) 将字符串的字符转换成小写字母。 Ucase$(字符串) 将字符串的字符转换成大写字母。 Left$(字符串，表达式) 左截取字符函数 Right$(字符串，表达式) 右截取字符函数

Mid$(字符串，起始位置，长度) 从字符串的起始位置截取指定的长度的字符串。 LTrim$(字符串) 将字符串的左边空格去掉返回字符串。 RTrim$(字符串) 将字符串的右边空格去掉返回字符串。 Trim$(字符串) 将字符串的前后空格删除返回字符串。 InStr( ) 字符串匹配函数 StrComp() 字符串比较函数 Space(表达式) 返回由数值指定的一定数目的空格字符串。 ② 日期函数

函数表达式 函数功能说明

Year(日期表达式) 日期表达式为任意整数或字符型表达式，返回值为年份。 Month(日期表达式) 日期表达式为任意整数或字符型表达式，返回值为月份。

Weekday(日期表达式) 日期表达式为任意整数或字符型表达式，返回值为星期。 Day(日期表达式) 日期表达式为任意整数或字符型表达式，返回值为月份的天数。 Hour(日期表达式) 表达日期表达式对应天中的小时。 Minute(日期表达式) 表示指定时间的分钟数。 Second(日期表达式) 表示指定时间的秒数。 Now 返回计算系统的日期和时间。

Data$ 返回当前系统日期，以字符串格式：月－日－年或月/日/年表示。 Time$ 返回当前系统日期，以字符串格式：时：分：秒。

④转换函数及其它

函数表达式 函数功能说明 函数表达式 功能

Str$(数值表达式) 字符串函数，将数值表达式转换成字符串。 Val(字符串) 将文本字符串变为数值型数据。 Asc(字符) 返回字符的ASCII码值。

Chr$(整数) 将0－255之间的整数转换为对应的ASCII字符。

EOF(文件号) 读取文件进程函数，读到文件末尾，则其函数值为 True或1。 Error(错误代码) 错误响应函数，返回错误代码所对应的错误信息字符串。 Format(变量, 格式字符串) 输出格式函数。

RGB(R,G,B) 颜色定义函数， R、G、B为数值表达式。

Sell（字符串） 可执行程序运行函数，字符串：为执行程序的路径。 Tab（整数） 文件输出定位函数。

? 用户自定义函数（Funtion Procedure）

其用途与建立方法类似于通用过程，只是通用过程是单方向调用，只有参数传给过程，而没有参数值的返回；而用户自定义函数是双向的，调用时参数传入函数，函数执行完毕后返回其函数值。故用户自定义函数像变量一样有自己的类型，它决定了函数返回值的类型。其描述为

Funtion 函数名称（参数1，参数2，…）As 类型名称 程序语句代码 End Funtion

Funtion 函数名称（参数1，参数2，…）As 类型名称 程序语句代码1 Exit Funtion 程序语句代码2 End Funtion

3、 方法（Method）

面向对象的程序设计语言为程序设计人员提供了一种特殊的过程和函数称为方法。在VB中一些通用的函数与过程编好并封装起来，作为方法供用户直接使用。方法是针对特定对象执行一项任务的过程或函数。如在早期的 Basic 语言中，往屏幕上显示信息和像打印机打印该信息，其语句是有区别的，即用 PRINT 语句表示像屏幕打印，而 LPRINT才是向打印机打印某信息，不同的对象完成同一任务其命令语句是不同的。而 VB 语言将打印功能封装成一特殊的 Print 方法，向不同的对象上打印信息直接指明对象，调用同一方法即可完成。

方法调用的语法为：

对象名称.方法 1.6 应用程序的设计 1、设计应用程序的界面 设计一个窗体

窗体对象是VB应用程序设计的基本构造模块，是运行应用程序时与用户交互操作的实际窗口。窗体有自己的属性、事件和方法。设计窗体的第一步是设置它的属性，窗体的属性很多，它不仅控制着窗体的外观，还控制着窗体的位置、行为等其它方面。属性可以在设计程序时在\属性\窗体中设置，也可以在程序运行时由代码来实现。窗体的属性如：BorderStyle、Caption 、Height、 Left MaxButton、MinButton、Moveable、Name、ShowInTaskar、WindowState、Icon等。

向窗体上添加控件

根据自身的需要在窗体上加不同的控件，并使用代码控制来完成不同的任务。向窗体上添加控件要使用控件工具箱和窗体编辑器。实用工具箱向窗体添加控件有两种途径：其一在工具箱中的控件按钮上双击，则窗体的中央会出现一个相应的控件；其二在工具箱中的控件按钮上单击，则该按钮会凹下去，鼠标指针变为＋形状，然后在窗体的合适位置按下鼠标左键即可。要删除不合适的控件只要选中然后按下DEL键即可。 设置启动窗体

除了窗体的细节设计以外，还要考虑应用程序的开始与结束。每个应用程序都有自己的入口及开始执行的地方。这里可以使用两种方法来加以实现。其一设置启动窗体，从\工程\菜单中选择\工程属性\命令，在显示对话框中选择\通用\选项卡，在\启动对象\列表中选取新启动的窗体，单击\确认\按钮即可。其二采用不使用启用窗体开始运行程序的方法，可在标准模块中创建一个名为 Main 的过程，如 Sub Main( ) 过程代码 End If Main 过程必须在一个标准模块内，不能在窗体模块内。要将 Sub Main 过程设为启动对象，可从\工程\菜单中选择\工程属性\命令，在显示对话框中选择\通用\选项卡，在\启动对象\列表中选定\，单击\确认\按钮即可。 使用函数生成的对话框 在应用程序中，可能会需要显示一些暂时性的简短的错误或警告信息，可以引起用户的注意。用户可以设计一个窗体来完成这个任务，但最简便的方法是使用 MsgBox 函数来完成则更直接、更为方便。MsgBox 函数可以用来在对话框中显示消息，并等待用户单击按钮，然后返回一个整形的值，让程序了解用户单击的是哪个按钮。MsgBox 函数的语法为： MsgBox( prompt [,buttons] [,title] [,helpfile,context])

另外可采用 InputBox 对话框实现一些简单的数据或信息的输入，并返回包含文本框的内容的字符串，InputBox 对话框的语法为：

InputBox（prompt [,title] [，default] [,xpos] [,ypos][,helpfile,context])

2、编写程序代码 赋值语句

VB的程序代码由语句、常数和声明部分组成。其中赋值语句使用频率最高，其语法为: 对象属性或变量＝表达式 程序的书写规则

注释，注释语句可用来说明编写的某段代码或声明某个变量的目的，以方便以后阅读这些源代码。要添加注释，使用\符号作为注释文字的开头。

断行，如果一个很长，打印和阅读恨不方便，可采用续行符\－\（一空格后紧跟一下划线）将长语句分成多行。

将多语句写在一行，VB通常是一行一条语句，如果在一行中写下多条语句，可使用\：\作为分隔符号。 变量

变量的命名：必须以字母开头，不能在变量名中出现 . 、空格或嵌入！，#，$，%，&等符号，变量名的长度不得超过255个字符。

变量的声明：变量的声明语句为 Dim 变量名 As 类型。

变量的作用范围：如果同一窗体的所有过程分享同一变量，就应该把它定义为模块级变量，其方法是在窗体模块的声明段中定义该变量。单击窗体模块代码窗口的对象列表框，从中选择\通用\选项即可。在窗体模块的声明段声明变量，在除了使用 Dim 关键字外，还可以使用Public 和Private 关键字。用Public 关键字声明模块级变量，变量在整个应用程序中由效，称为公共变量或全局变量，其它模块中的过程也可以使用这个变量。用Private 关键字声明模块级变量，本窗体中的过程可以访问它，但其它模块中的过程不能使用这个变量。与模块级变量相对，在过程中声明的变量被称为局部变量，局部变量只能在过程执行期间有效，其它代码不能使用。如果过程结束以后还需保持过程中变量的值，可使用 Static 关键字声明变量为静态变量。

不同作用范围的3种变量的声明方式 作用范围 局部变量 模块级变量 公共变量 声明方式 Dim, Static Dim, Private Public

变量声明的位置 过程中 模块的声明段中 模块的声明段中 是否被本模块中其它过程访问 × √ √ 能否被其它模块访问 × × √

常数

在应用程序之中，往往要用到一些不变的量即常量，如pi=3.11415926.有在VB中，声明常数的语法为

[Public|Private] Const 常数名[ As 类型]=表达式 运算符号

算术运算：加 +、减 －、乘 *、浮点数除法 /、整数除法 \\、幂运算 ^、求余数 MOD。 比较运算：大于 >、小于 <、大于或等于 >= 、小于或等于 <=、等于 = 、不等于 <>。 连接运算：+ 、& 。 逻辑运算：逻辑非 Not 、逻辑与 And、逻辑或 Or、逻辑异或 Xor、逻辑等 Eqv、蕰含 Imp。

流程控制语句

Visual Basic支持的条件判定结构有三种： If…Then 结构：If 条件 Then 语句 或 If 条件 Then 若干语句

End if If…Then…Else 结构：

If 条件1 Then 若干语句1

[Else If 条件2 Then 若干语句2 ]

Else

若干语句3 End if Select Case 结构：

Select Case 条件 Case 表达式1 若干语句1 Case 表达式2 若干语句2 ： ： Case else

若干语句n End Select

Visual Basic支持的条件循环控制： D0 While 循环条件 若干语句 [Exit Do]

Loop D0

若干语句

Loop While 循环条件

For 计算器 = 初值 To 终止值 Step 增量 若干语句

Next 计算器 1.7 简单的编程实例

运行VB6，新建一个标准的EXE工程，从工具箱中双击 CommandButton 控件，在主窗体 Form1上设置两个命令按钮 Command1 与 Command2，双击Text文本控件，在主窗体 Form1上设置三个文本控件，其 Name 属性依次为 Text1、Text2、Text3，最后拖动一个标签 Label 控件放在主窗体的上中间位置。双击窗体 Form1，打开代码窗口，输入以下代码 Private Sub Command1_click() ' 单击按钮事件 Text3.text=Text1.text+Text2.text End sub

Private Sub Command2_click()' 单击按钮事件 End ' 程序运行结束退出 End sub

Private Sub Form_load( ) ' 窗体启动调入内存 Label1.Caption=\大地测量学基础实验\Command1.Caption=\计算\Command2.Caption=\退出\

Text1.text=\：Text2.text=\：Text3.text=\给文本框赋初值 Text1.SetFocus' Text1文本框获得焦点 End sub

然后运行，在Text1文本框与Text2文本框中输入数据，单击计算按钮，在Text3文本框中显示计算结果。

一个较好的应用程序必须要有好的用户界面，在大多应用程序中大多包含着许多通用的东西，如工具栏、状态栏、工具提示、上下文菜单、帮助以及选项卡对话框等。VB 具有把所有这些添加到应用程序中的能力，因此应该很好的学习并加以掌握运用。

二 测量中几种常用的计算 角度化为弧度、弧度化为角度以及坐标方位角的计算在测量数据处理中是经常遇到的计算问题。也是最基本的计算问题。

2.1 角度化为弧度

将度、分、秒形式的角度angle化为弧度，采用函数功能来实现，其函数值为返回的弧度值。其数据格式为: 参数angle的整数为表示度，小数点后两位表示分, 小数点后第三位开始表示秒， 如180度34分54.23秒, 即 angle=180.345423。 Function Radian(ByVal angle As Double) As Double Dim mm As Double, a As Double a = Abs(xx)

a = a + 0.0000001 dd% = Int(a)

ii% = Int((a - dd%) * 100) mm = (a - dd%) * 100 - ii% mm = dd% + ii% / 60 + mm / 36 Radian = pi * mm / 180

Radian = Sgn(angle) * Radian End Function

2.2 弧度化为角度计算

将弧度值radian 化为度、分、秒的角度形式，采用函数功能实现，计算返回角度值。其数据格式为: 函数值qdms的整数为表示度，小数点后两位表示分, 小数点后第三位开始表示秒， 如160度14分23.03秒, 即qdms=160.142303。 Function qdms(ByVal radian As Double) As Double Dim a1, xx , second As Double a1 = Abs(radian) a1 = a1 * 180 / pi degree% = Int(a1)

xx = (a1 - degree%) * 60 minute% = Int(xx)

second = (xx - minute%) * 60

qdms = degree% + minute% / 100 + second / 10000 qdms = Sgn(radian) * qdms End Function

2.3 坐标方位的计算

已知点号为 和 两点的平面坐标 、 、 、 ，求 和 两点的坐标方位角, 先计算 ，在按如下计算公式如下：

Function qiua(ByVal px As Double, ByVal py As Double) As Double Dim PZ As Double If px = 0 Then PZ = pi / 2#

If py < 0 Then PZ = pi * 1.5 Else

PZ = Atn(py / px)

If px < 0 Then PZ = PZ + pi

If py < 0 And px > 0 Then PZ = PZ + pi * 2 End If qiua = PZ End Function

三 空间大地坐标与空间直角坐标之间的换算 3.1 由空间大地坐标计算空间直角坐标

已知椭球 、 ，大地坐标直角坐标X、Y、Z，求大地坐标B、L、H，计算公式如下 式中 。

3.2 空间直角坐标计算空间大地坐标

已知椭球 、 ，大地坐标直角坐标X、Y、Z，求大地坐标B、L、H。 1、 直接解法计算公式（参见熊介编著：《椭球大地测量学》） 计算辅助量 计算大地经度 计算大地纬度

计算大地高

2、 迭代解法计算公式汇编 计算辅助量

， ， 计算大地纬度

计算大地纬度

令 ， ，

迭代循环计算 直到满足 ，以保证 的计算精度至 计算大地高

3.3 计算程序与范例

空间大地坐标与空间直角坐标之间的计算程序界面如下图所示。在窗体XYZ-BLH上设置3个Frame控件，在 Frame1 内放置3个选择按钮 option1(0)、option1(1)、option(2)，控制椭球参数；在Frame2 内放置2个选择按钮 option2(0) 与 option2(1)，控制计算方法；在Frame3 内放置1个MSFlexGrid 控件，4个 CommandButton 命令按钮，其 Caption 属性分别为文件、添加、删除、保存，其功能主要用来编辑数据与保存数据；在窗体 XYZ-BLH上 设置3个 CommandButton 命令按钮，其 Caption 属性分别为计算、显示、退出，其功能为计算、用 Word 显示计算结果、退出。其计算界面直观、方便。 计算范例1：

已知数据 ， ， 。

计算结果 ， ， 。 计算范例2：

已知数据 ， ， 。 计算结果 ， ， 。

图3－1 空间大地坐标与空间直角坐标换算程序界面

上述计算采用直接算法，计算子程序如下： Sub BLH_XYZ(b, l, H, x, y, Z, ByVal K As Integer) EPS = e2 / (1# - e2) BB = a / Sqr(1# + EPS) Select Case K Case 1

p = Sqr(x ^ 2 + y ^ 2) UO = Atn(Z * a / p / BB)

b = Atn((Z + EPS * BB * Sin(UO) ^ 3) / (p - e2 * a * Cos(UO) ^ 3)) p1 = Sqr(1# + (1# - e2) * Tan(b) ^ 2) b = Atn(Z / p + a * e2 * Tan(b) / p / p1) l = Atn(y / x)

If b < 0# Then b = b + 2# * pi If l < 0# Then l = l + 2# * pi U = Atn(BB / a * Tan(b))

H = Sqr((p - a * Cos(U)) ^ 2 + (Z - BB * Sin(U)) ^ 2) If (p - a * Cos(U)) < 0# Then H = -H Case 2

RN = a / Sqr(1# - e2 * Sin(b) * Sin(b))

x = (RN + H) * Cos(b) * Cos(l) y = (RN + H) * Cos(b) * Sin(l) Z = (RN * (1# - e2) + H) * Sin(b) End Select End Sub

四 大地主题问题正反计算

参考椭球面是大地测量计算的基准面。大地坐标是椭球面上的基本坐标系，根据大地测量的观测成果（如距离与方向），从大地原点出发，逐点计算在椭球面上的大地坐标；或根据两点的大地坐标，计算它们之间的大地线长度和大地方位角，这类计算称为大地问题解算，又叫大地主题解算。大地问题解算的用途是多方面的，随着现代空间技术和航空航天、航海等领域的发展，大地问题解算（尤其是大地反算）有着更为重要的作用。鉴于各种用途与要求不同，产生了不同的大地解算方法与公式。

椭球面上大地坐标的解算比平面坐标的解算要复杂得多，正是由于这种复杂性导致了大地问题解算公式的多样化，其解算方法多大几十种。若按解算的距离来分类，一般分为短距离（400公里以内），中距离（400～1000公里）和长距离（1000～2000公里）；若按解法分类可分为直接解法和间接解法；若按解算精度来分类，又可分为精密公式与近似公式等。这里介绍三种方法即高斯平均引数法、白塞尔方法与嵌套系数法的计算过程、步骤与计算程序。

4.1 高斯平均引数大地问题解算

1、高斯平均引数正算计算公式（ km） 1）计算辅助量公式

2) 计算 、 、 的初值 3）计算 、 、

4）再次计算 、 、

5）重复计算（3）,直到计算满足

如按弧度计算可取 ，按角度计算可取 。 6）计算 、 、 的最后值

如果计算距离小于70km时，上述计算公式中的 可以略去，得到如下简化计算方法： 式中

编程计算步骤如下： 1）计算辅助量

2) 计算 、 、 的初值

3）计算两点的平均值 4）迭代计算 、 、 5）计算辅助量

6）计算 、 、 的终值

7）重复计算（3）-(6)步,直到计算满足

如按弧度计算可取 ，按角度计算可取 。 8）、计算 、 、 的最后值

2、高斯平均引数反算公式（ km） 1）、 2)、 式中各系数

注：这里对教材公式中相应系数进行了修正与改进。 3）、

4)、

5)、 或 ，

3、算例与计算程序 1）、正算范例 算例 例1( 30km ) 例2(50 km) 例3(80 km) 例4( 400 km) 参考椭球 克拉索夫斯基 克拉索夫斯基 克拉索夫斯基已 知 数 据

克拉索夫斯基

m m m m 参 考 值

注：例1取自朱华统：《椭球大地计算》P97算例；例2取自陈键、晁定波主编《椭球大地测量学》P90面的数据；例3、例4取自赵文光著《椭球大地测量学》。

大地主题正算

==================

计算方法： 高斯平均引数法 ( S< 200公里 ) 参考椭球：克拉索夫斯基椭球

N01 NO2 B1/B2 L1/L2 A12/A21 S12 1 2 30.29582043 120.05402184 247.2750428 28230.9350 30.24058365 119.49233864 67.19353743

3 4 40.02356784 130.10122676 328.12367500 48741.7580 40.24574365 129.52031587 148.00533348

5 6 40.02356784 130.10122676 1.49430000 80000.0000 40.45479037 130.12011050 181.50535462

7 8 40.02356784 115.10000000 36.12010300 414306.5380 *43.00558744 118.10030113 218.11268087 2)、反算范例 算例 例1( 30km ) 例2(50 km) 例3(80 km) 例4( 410 km) 参考椭球 克拉索夫斯基 贝塞尔椭球 克拉索夫斯基 克拉索夫斯基 已 知 数 据 m 参 考 值

注：例1取自朱华统：《椭球大地计算》P100算例；例2自陈键、晁定波主编《椭球大地测量学》P98面的数据；例3、例4取自赵文光著《椭球大地测量学》。

大地主题反算

================== 计算方法： 高斯平均引数法 ( S< 200公里 ) 参考椭球：克拉索夫斯基椭球

N01 NO2 B1/B2 L1/L2 A12/A21 S12 1 2 30.29582043 120.05402184 247.27504227 28230.9379 30.24058354 119.49233853 67.19353680 30.24058354 119.49233853 67.19354031

*3 4 53.50028809 10.12041772 25.16319782 47652.5967 54.13152891 10.30472430 205.31408796

5 6 40.02356784 130.10122627 1.49432983 80000.0044 40.45479027 130.12011040 181.50538462

*7 8 40.02356784 115.10000000 36.12012391 414306.5531 43.00558784 118.10030000 218.11265892

注：例2参考椭球为贝塞尔椭球；例4计算结果未能满足计算精度。 3)、高斯平均引数正反算计算程序

Sub DDZT_GS(B1 As Double, L1 As Double, A1 As Double, S12 As Double, B2 As Double, L2 As Double, A2 As Double, k As Integer) eps = e2 / (1# - e2) Select Case k Case 1

t1 = Sin(B1): t2 = Cos(B1)

NB = a / Sqr(1# - e2 * t1 * t1)

MB = a * (1 - e2) / Sqr((1# - e2 * t1 * t1) ^ 3) dB = S12 * Cos(A1) / MB

dL = S12 * Sin(A1) / NB / Cos(B1) dA = dL * Sin(B1) Do

dB0 = dB: dL0 = dL: dA0 = dA

Bm = B1 + dB0 / 2: LM = L1 + dL0 / 2: Am = A1 + dA0 / 2 t1 = Sin(Am)：t2 = Cos(Am)：t = Tan(Bm) q = eps * Cos(Bm) * Cos(Bm) n = a / Sqr(1# - e2 * Sin(Bm) * Sin(Bm)) MB = a * (1 - e2) / Sqr((1# - e2 * Sin(Bm) * Sin(Bm)) ^ 3) db1 = t1 * t1 * (2 + 3 * t * t + 2 * q) db2 = 3 * t2 * t2 * q * (t * t - 1 - (1 + 4 * t * t) * q)

dB = S12 * t2 * (1 + S12 * S12 * (db1 + db2) / n / n / 24) / MB xt = t * t xt1 = t1 * t1 xt2 = t2 * t2 dL1 = S12 * S12 * (xt * xt1 - xt2 * (1 + (1 - 9 * xt) * q)) / n / n / 24 dL = S12 * t1 * (1 + dL1) / n / Cos(Bm) da1 = xt2 * (2 + (7 + 9 * xt + 5 * q) * q) da2 = xt1 * (2 + xt + 2 * q)

dA = S12 * t * t1 * (1 + S12 * S12 * (da1 + da2) / NB / NB / 24) / NB XdB = Abs(dB - dB0) * P0 '常数P0=206265 XdL = Abs(dL - dL0) * P0 XdA = Abs(dA - dA0) * P0

Loop While XdB > 0.0001 And XdL > 0.0001 And XdA > 0.0001 B2 = B1 + dB: L2 = L1 + dL: A2 = A1 + dA If A1 > pi Then A2 = A2 - pi If A1 < pi Then A2 = A2 + pi Case 2

Bm = (B1 + B2) / 2: dB = B2 - B1: dL = L2 - L1 CosBm = Cos(Bm)

NB = a / Sqr(1# - e2 * Sin(Bm) * Sin(Bm)) t = Tan(Bm)

q = eps * CosBm * CosBm

v2 = 1 + q: v4 = v2 * v2: v6 = v4 * v2 r10 = NB * CosBm

r12 = NB * CosBm * (1 + q - 9 * q * t * t + q * q) / 24 / v4 r30 = -NB * CosBm ^ 3 * t ^ 2 / 24 s01 = NB / v2

s21 = -NB * CosBm ^ 2 * (2 + 3 * t * t + 2 * q) / 24 / v2 s03 = NB * q * (1 - t * t) / 8 / v6 't10 = CosBm * t

't12 = CosBm * t * (3 + 2 * q - 2 * q * q) / 24 't30 = CosBm ^ 3 * t * (1 + q) / 12 t10 = CosBm * t

t12 = CosBm * t * (2 + 7 * q + 9 * t * t * q) / 24 / v4 t30 = CosBm ^ 3 * t * (2 + t * t + 2 * q) / 24

SsinAm = r10 * dL + r12 * dL * dB ^ 2 + r30 * dL ^ 3 ScosAm = s01 * dB + s21 * dB * dL ^ 2 + s03 * dB ^ 3 dA = t10 * dL + t12 * dL * dB ^ 2 + t30 * dL ^ 3 '计算大地正反方位角

Am = Atn(SsinAm / ScosAm) If Am < 0 Then Am = Am + pi If SsinAm < 0 Then Am = Am + pi A1 = Am - dA / 2 A2 = Am + dA / 2

If A1 > pi Then A2 = A2 - pi If A1 < pi Then A2 = A2 + pi S12 = SsinAm / Sin(Am) End Select End Sub

4.2 白塞尔大地主题解算

高斯平均大地问题解算，是将大地经差、大地纬差和方位角展开为大地线的幂级数，在一定的精度条件下，距离愈长，公式的结构就愈复杂，收敛就愈慢。因此这种公式不适合于解长距离大地问题。

白塞尔大地问题解算公式是首先建立以椭球中心为中心，以任意长（一般为单位长）为半径的辅助球，按以下三个步骤解算：

1．按一定条件将椭球面元素投影到球面上； 2．在球面上解算大地问题；

3．将求得的球面元素按投影关系换算为相应的球面元素。 白塞尔大地问题解算公式的三个投影条件为：

1．使投影后球面上点的球面纬度等于椭球面上对应点的归化纬度； 2．椭球面上两点间的大地线投影到辅助球面上为大园弧； 3．大地方位角投影后数值保持不变。 1、白塞尔大地主题正算 已知 ，计算 。

1） 将椭球面元素投影到球面上 （1） 由 求 ： （2） 计算辅助量 和 ，

（3）计算球面长度，将 化为 式中分别系数为

上式右端含有代求量 ，因此需要迭代计算。第一次迭代取近似值 ，第二次计算取

以后计算用 代换 代如上式迭代计算，直到所要求的精度为止。一般取 。 2）解算球面三角形 （1）计算 （2）计算 或

（3）计算

或

3）将球面元素换算到椭球面上 （1）由 求 或

（2）将球面经差 化为椭球面经差 ，求 式中

式中 的最大值为 ，故在计算时通常可以略去不计。 （3）象限的判定（参照《大地测量基础》教材 P102）

SinA1符号 + + － － 符号

+ － + － | |

SinA1符号 － － + + TanA2符号 + － + －

其中 ， 为锐角。

2、 白塞尔大地主题反算 已知 ，计算 ， 。

1）将椭球面元素投影到球面上 （1）由 求

（2）采用逐次趋近方法，由 计算

在反算中，已知椭球面上经差 ，球面经差上的对应经差 未知，为了由 求 ，由下式可知还需计算 、 、 ，计算 又需要 量，故需要进行迭代计算。 第一次趋近，取 ； ， 或 在此要对 象限进行判断

P符号 + + － － q符号 + － － + =

在此要对 象限进行判断

符号 + － =

仿照上述计算步骤迭代计算，直到 为止。 2）将球面元素换算到椭球面上

或

在此要对 象限进行判断

3、 白塞尔大地主题算例与计算程序 1）正算算例 算例 例1( 80km ) 例2(400km) 例3(8000 km) 例4( 15000 km) 参考椭球 克拉索夫斯基 克拉索夫斯基 克拉索夫斯基 已 知 数 据 m m m m 参 考 值

克拉索夫斯基

注：例1、例2取自赵文光著《椭球大地测量学》P88和P93例，例3取自陈键、晁定波主编《椭球大地测量学》P111面的数据。例4取自周江华：《测绘通报》，2002，5（6）：108－111，编程计算结果如下：

大地主题正算

================== 计算方法： 贝塞尔算法( S< 20000 公里 ) 参考椭球：克拉索夫斯基椭球

N01 NO2 B1/B2 L1/L2 A12/A21 S12 1 2 40.02356784 130.10122676 1.49430000 80000.0000 40.45479037 130.12011050 181.50535462

3 4 40.02356784 115.10000000 36.12010300 414306.5380 43.00558793 118.10030013 218.11268010

5 6 68.58000000 33.05000000 339.4956385 7999606.380 37.44599767 -122.26000057 9.01078092

7 8 35.00002200 90.00001100 100.00003300 15000000.200 -30.29209655 215.59043332 290.32533894 2） 反算算例 算例 例1(80 km) 例2( 410 km) 例4( 8000 km) 例5( 410 km) 参考椭球 克拉索夫斯基 克拉索夫斯基 克拉索夫斯基 克拉索夫斯基 已 知 数 据

参 考 值

注：例1、例2取自赵文光著《椭球大地测量学》P88和P93例，例3取自陈键、晁定波主编《椭球大地测量学》P113面的数据。例4取自周江华：《测绘通报》，2002，5（6）：108－111，编程计算结果如下：

大地主题反算

================== 计算方法： 贝塞尔算法( S< 20000 公里 )

参考椭球：克拉索夫斯基椭球

N01 NO2 B1/B2 L1/L2 A12/A21 S12 1 2 40.02356784 130.10122627 1.49432980 80000.0045 40.45479027 130.12011040 181.50538463

3 4 40.02356784 115.10000000 36.12010267 414306.5362 43.00558784 118.10030000 218.11267964

5 6 68.58000000 33.05000000 339.49563867 7999606.3891 37.44599755 -122.26000057 9.01078088

7 8 35.00002200 90.00001100 100.00003279 15000000.3380 -30.29209640 215.59043380 290.32533887

3） 贝塞尔算法计算程序

Sub DDZT_Bessel (B1 As Double, L1 As Double, A1 As Double, S12 As Double, _ B2 As Double,L2 As Double, A2 As Double, k As Integer) eps = e2 / (1# - e2) b = a / Sqr(1# + eps) Select Case k Case 1

U1 = Atn(Sqr(1 - e2) * Tan(B1))

sinA0 = Cos(U1) * Sin(A1): cosA0 = Sqr(1 - sinA0 * sinA0) sigma1 = Atn(Tan(U1) / Cos(A1)) xk2 = eps * cosA0 * cosA0 xk4 = xk2 *xk2 xk6 = xk4 * xk2

alpha = (1 - xk2 / 4 + 7 * xk4 / 64 - 15 * xk6 / 256) / b beta = xk2 / 4 - xk4 / 8 + 37 * xk6 / 512 gamma = xk4 / 128 - xk6 / 128 sigma = alpha * S12 Do

sigma0 = sigma

sigma = alpha * S12 + beta * Sin(sigma0) * Cos(2 * sigma1 + sigma0) sigma = sigma + gamma * Sin(2 * sigma0) * Cos(4 * sigma1 + 2 * sigma0) Dsigma = Abs(sigma - sigma0) * P0 '常数P0=206265 Loop While Dsigma > 0.0001 '计算反方位角A2

sinA2 = Cos(U1) * Sin(A1)

cosA2 = Cos(U1) * Cos(sigma) * Cos(A1) - Sin(U1) * Sin(sigma) tanA2 = sinA2 / cosA2: A2 = Abs(Atn(sinA2 / cosA2)) SinA1 = Sin(A1)

If SinA1 < 0 And tanA2 > 0 Then A2 = A2 If SinA1 < 0 And tanA2 < 0 Then A2 = pi - A2 If SinA1 > 0 And tanA2 > 0 Then A2 = pi + A2 If SinA1 > 0 And tanA2 < 0 Then A2 = 2 * pi - A2 '计算大地纬度B2

sinU2 = Sin(U1) * Cos(sigma) + Cos(U1) * Cos(A1) * Sin(sigma) B2 = Atn(sinU2 / Sqr(1 - e2) / Sqr(1 - sinU2 * sinU2)) '计算大地经度L2

sinl = Sin(A1) * Sin(sigma)

cosl = Cos(U1) * Cos(sigma) - Sin(U1) * Sin(sigma) * Cos(A1) tanlambda = sinl / cosl: lambda = Abs(Atn(sinl / cosl)) If tanlambda > 0 And SinA1 > 0 Then lambda = lambda If tanlambda < 0 And SinA1 > 0 Then lambda = pi - lambda If tanlambda < 0 And SinA1 < 0 Then lambda = -lambda If tanlambda > 0 And SinA1 < 0 Then lambda = lambda - pi e4 = e2 * e2 e6 = e4 * e2

xk2 = e2 * cosA0 * cosA0 xk4 = xk2 * xk2 xk6 = xk4 * xk2

alpha1 = (e2 / 2 + e4 / 8 + e6 / 16) - e2 * (1 + e2) * xk2 / 16 + 3 * xk4 * e2 / 128 beta1 = e2 * (1 + e2) * xk2 / 16 - e2 * xk4 / 32 gamma1 = e2 * xk4 / 256

xx = alpha1 * sigma + beta1 * Sin(sigma) * Cos(2 * sigma1 + sigma) xx = xx + gamma1 * Sin(2 * sigma) * Cos(4 * sigma1 + 2 * sigma) l = lambda - sinA0 * xx L2 = L1 + l Case 2

U1 = Atn(Sqr(1 - e2) * Tan(B1)): U2 = Atn(Sqr(1 - e2) * Tan(B2)) DL = L2 - L1

sa1 = Sin(U1) * Sin(U2): sa2 = Cos(U1) * Cos(U2) cb1 = Cos(U1) * Sin(U2): cb2 = Sin(U1) * Cos(U2) lambda = DL Do

lambda0 = lambda

p = Cos(U2) * Sin(lambda0): q = cb1 - cb2 * Cos(lambda0) '计算方位角A1 A1 = Abs(Atn(p / q))

If p > 0 And q > 0 Then A1 = A1 If p > 0 And q < 0 Then A1 = pi - A1 If p < 0 And q < 0 Then A1 = pi + A1 If p < 0 And q > 0 Then A1 = 2 * pi - A1 '计算sigma

Ssigma = p * Sin(A1) + q * Cos(A1) Csigma = sa1 + sa2 * Cos(lambda0) sigma = Abs(Atn(Ssigma / Csigma)) If Csigma > 0 Then sigma = sigma If Csigma < 0 Then sigma = pi - sigma '计算A0 与 sigma1

sinA0 = Cos(U1) * Sin(A1)

sigma1 = Atn(Tan(U1) / Cos(A1)) '计算椭球面经差lambda

cosA0 = Sqr(1 - sinA0 * sinA0) e4 = e2 * e2 e6 = e4 * e2

xk2 = e2 * cosA0 * cosA0 xk4 = xk2 * xk2 xk6 = xk4 * xk2

alpha1 = (e2 / 2 + e4 / 8 + e6 / 16) - e2 * (1 + e2) * xk2 / 16 + _ 3 * xk4 * e2 / 128

beta1 = e2 * (1 + e2) * xk2 / 16 - e2 * xk4 / 32 gamma1 = e2 * xk4 / 256

xx = alpha1 * sigma + beta1 * Sin(sigma) * Cos(2 * sigma1 + sigma) xx = xx + gamma1 * Sin(2 * sigma) * Cos(4 * sigma1 + 2 * sigma) lambda = DL + sinA0 * xx

Dlambda = Abs(lambda - lambda0) * P0 'P0=206265 Loop While Dlambda > 0.0001 '计算椭球面距离S12

cosA0 = Sqr(1 - sinA0 * sinA0) xk2 = eps * cosA0 * cosA0 xk4 = xk2 * xk2 xk6 = xk2 * xk4

alpha = (1 - xk2 / 4 + 7 * xk4 / 64 - 15 * xk6 / 256) / b beta = xk2 / 4 - xk4 / 8 + 37 * xk6 / 512 gamma = xk4 / 128 - xk6 / 128

xs12 = gamma * Sin(2 * sigma) * Cos(4 * sigma1 + 2 * sigma)

S12 = (sigma - beta * Sin(sigma) * Cos(2 * sigma1 + sigma) - xs12) / alpha '计算椭球面反方位角A21 sinA2 = Cos(U1) * Sin(A1)

cosA2 = Cos(U1) * Cos(sigma) * Cos(A1) - Sin(U1) * Sin(sigma) tanA2 = sinA2 / cosA2: A2 = Abs(Atn(sinA2 / cosA2)) SinA1 = Sin(A1)

If SinA1 < 0 And tanA2 > 0 Then A2 = A2 If SinA1 < 0 And tanA2 < 0 Then A2 = pi - A2 If SinA1 > 0 And tanA2 > 0 Then A2 = pi + A2 If SinA1 > 0 And tanA2 < 0 Then A2 = 2 * pi - A2 End Select End Sub

4.3 嵌套系数法解算任意距离的大地主题

目前大地主题解算方法虽然比较完善，但大多数解算方法（特别是反算）仍然存在一些缺陷，主要表现在应用范围受到一定的限制。受其计算精度的影响，它们几乎都只适合一定范围而不能解任意大地距离的主题问题。另外产生奇异性，当两大地点位于同一子午圈或同

在赤道时算法不能正常计算。电脑的广泛应用为改进大地主题的解算方法提供了新的条件，在电算中进行迭代计算，简单易行，这里介绍一种嵌套系数法精密解算任意距离的大地主题问题，仅作为教材内容的补充。 1、符号意义

， 参考椭球的长、短半轴； 参考椭球的扁率； 参考椭球的二偏心率；

、 椭球面上两点的大地纬度； 、 椭球面上两点的大地经度； 两点大地经差；

椭球面上两点的大地线长度； 、 两点的大地线方位角； 、 两点的归化纬度；

大地线S 的最高纬度； 两点的球面经差； 两点的球面角距；

大地线自赤道交点至1点间的球面角距； 大地线自赤道交点至2点间的球面角距；

大地线自赤道交点至1、2点间的球面角距之和 ； 球面经差 的改正项； 球面经差 的改正项； 、 计算嵌套系数的算子；

、 计算椭球面距离改正项的嵌套系数；

计算椭球面经差改正项的嵌套系数。 2、嵌套系数法大地主题正算 已知： ，计算 。 1）计算准备

2）嵌套系数计算 ，

3）球面角距 改正项 的迭代计算； （首次取 ） 或 ；

采取迭代算法，直到 的变化小于限差为止，一般取限差为 弧度。 4）计算 ； 或

，在此需要判断象限； ； ；

，在此需要判断象限； 5） 检验计算 。

3、嵌套系数法大地主题反算 已知： ，计算 。 1）计算准备

2）球面经差 的改正数 的计算 （首次取 ） ，并判断象限 ；

采取迭代算法，直到 的变化小于限差为止，一般取限差为 弧度。 3）球面距离 改正项 的迭代计算

4）计算 、 、 ；

，并判断象限； ，并判断象限； 5）检核计算 。

4、嵌套系数法大地主题算例与计算程序 1）正算算例

大地主题正算

================== 计算方法：嵌套系数法( 适合任意距离 ) 参考椭球：克拉索夫斯基椭球

N01 NO2 B1/B2 L1/L2 A12/A21 S12 1 2 30.29582043 120.05402184 247.27504300 28230.9350 30.24058365 119.49233864 67.19353763

3 4 47.46526740 35.49363300 44.12136600 44797.2820 48.04096664 36.14450516 224.30535480

5 6 40.02356784 130.10122676 1.49430000 80000.0000 40.45479037 130.12011050 181.50535462

7 8 40.02356784 115.10000000 36.12010300 414306.5380 43.00558793 118.10030013 218.11268010

9 10 68.58000000 33.05000000 339.49563900 7999606.3800 37.44599774 -122.26000117 9.01078071

11 12 35.00002200 90.00001100 100.00003300 15000000.20 -30.29209655 215.59043332 290.32533894

2）反算算例

大地主题反算

================== 计算方法： 嵌套系数法( 适合任意距离 ) 参考椭球：克拉索夫斯基椭球

N01 NO2 B1/B2 L1/L2 A12/A21 S12 1 2 30.29582043 120.05402184 247.27504216 28230.9375 30.24058354 119.49233853 67.19353669

3 4 47.46526470 35.49363300 44.12102618 44796.5558 48.04096384 36.14450004 224.30501114

5 6 40.02356784 130.10122627 1.49432980 80000.0045 40.45479027 130.12011040 181.50538463

7 8 40.02356784 115.10000000 36.12010267 414306.5361 43.00558784 118.10030000 218.11267964

9 10 68.58000000 33.05000000 339.49563867 7999606.3881 37.44599755 -122.26000057 9.01078088

11 12 35.00002200 90.00001100 100.00003279 15000000.3385 -30.29209640 215.59043380 290.32533886

注：上述正反算数据例1取自朱华统：《椭球大地计算》P100算例；例2自陈键、晁定波主编《椭球大地测量学》P98面的数据；例3、例4取自赵文光著《椭球大地测量学》P88和P93例，例5取自陈键、晁定波主编《椭球大地测量学》P111和P113面的数据。例54取自周江华：《测绘通报》，2002，5（6）：108－111。

3）嵌套系数法正反算计算程序

Sub DDZT_Nest(B1 As Double, L1 As Double, A1 As Double, S12 As Double, _ B2 As Double, L2 As Double, A2 As Double, k) eps = e2 / (1# - e2) '计算椭球二偏心率 alpha = 1 - Sqr(1 - e2) b = a / Sqr(1# + eps) Select Case k Case 1 '正算

u1 = Atn(Sqr(1 - e2) * Tan(B1)) sigma1 = Atn(Tan(u1) / Cos(A1))

cosUn = Cos(u1) * Sin(A1): sinUn2 = 1 - cosUn * cosUn

t = eps * sinUn2 / 4 V = alpha * sinUn2 / 4

xk1 = 1 + t * (1 - t * (3 - 5 * t + 11 * t * t) / 4) xk2 = t * (1 - t * (2 - t * (37 - 94 * t) / 8))

xk3 = V * (1 + alpha + alpha * alpha - V * (3 + 7 * alpha - 13 * V)) Dsigma = 0 Do

Dsigma0 = Dsigma

sigma = S12 / xk1 / b + Dsigma0 Msigma = 2 * sigma1 + sigma

xx = Cos(sigma) * Cos(2 * Msigma) + xk2 * (1 + 2 * Cos(2 * sigma)) _ * Cos(3 * Msigma) / 6

Dsigma = xk2 * Sin(sigma) * (Cos(Msigma) + xk2 * xx / 4) Loop While Abs(Dsigma - Dsigma0) * P0 > 0.00001

sinB2 = Sin(u1) * Cos(sigma) + Cos(u1) * Sin(sigma) * Cos(A1) sigma2 = sigma1 + sigma

cosB2 = Sqr(1 - e2) * Sqr(1 - sinUn2 * Sin(sigma2) * Sin(sigma2)) B2 = Atn(sinB2 / cosB2)

'sinU2 = Sin(U1) * Cos(sigma) + Cos(U1) * Cos(A1) * Sin(sigma) 'B2 = Atn(sinU2 / Sqr(1 - e2) / Sqr(1 - sinU2 * sinU2)) sinw = Sin(sigma) * Sin(A1)

cosw = Cos(u1) * Cos(sigma) - Sin(u1) * Sin(sigma) * Cos(A1) tanw = sinw / cosw: w = Abs(Atn(sinw / cosw)) sinA1 = Sin(A1)

If tanw > 0 And sinA1 > 0 Then w = w If tanw < 0 And sinA1 > 0 Then w = pi - w If tanw < 0 And sinA1 < 0 Then w = -w If tanw > 0 And sinA1 < 0 Then w = w - pi

dw0 = sigma + xk3 * Sin(sigma) * (Cos(Msigma) + xk3 * Cos(sigma) * Cos(2 * Msigma))

dw = (1 - xk3) * alpha * cosUn * dw0 L2 = L1 + w - dw

SinA2 = Cos(u1) * Sin(A1)

cosA2 = Cos(u1) * Cos(sigma) * Cos(A1) - Sin(u1) * Sin(sigma) tanA2 = SinA2 / cosA2: A2 = Abs(Atn(SinA2 / cosA2)) If sinA1 < 0 And tanA2 > 0 Then A2 = A2 If sinA1 < 0 And tanA2 < 0 Then A2 = pi - A2 If sinA1 > 0 And tanA2 > 0 Then A2 = pi + A2 If sinA1 > 0 And tanA2 < 0 Then A2 = 2 * pi - A2 H = Cos(u1) * sinA1 - Cos(u2) * Sin(A2) '检核计算 Case 2 '反算

u1 = Atn(Sqr(1 - e2) * Tan(B1)): u2 = Atn(Sqr(1 - e2) * Tan(B2)) dL = L2 - L1

sa1 = Sin(u1) * Sin(u2): sa2 = Cos(u1) * Cos(u2)

cb1 = Cos(u1) * Sin(u2): cb2 = Sin(u1) * Cos(u2) '计算椭球面经差w的改正数dw dw = 0 Do

dw0 = dw w = dL + dw0

xx1 = Cos(u2) * Sin(w): xx2 = cb1 - cb2 * Cos(w) Ssigma = Sqr(xx1 ^ 2 + xx2 ^ 2) csigma = sa1 + sa2 * Cos(w)

sigma = Abs(Atn(Ssigma / csigma)) If csigma > 0 Then sigma = sigma If csigma < 0 Then sigma = pi - sigma

cosUn = sa2 * Sin(w) / Sin(sigma): sinUn = Sqr(1 - cosUn * cosUn) Un = Atn(sinUn / cosUn)

sinUn2 = 1 - Cos(Un) * Cos(Un)

CMsigma = Cos(sigma) - 2 * sa1 / sinUn2 SMsigma = Sqr(1 - CMsigma * CMsigma) Msigma = Atn(SMsigma / CMsigma)

If CMsigma < 0 Then Msigma = pi + Msigma V = alpha * sinUn2 / 4

xk3 = V * (1 + alpha + alpha * alpha - V * (3 + 7 * alpha - 13 * V))

xdw = sigma + xk3 * Sin(sigma) * (Cos(Msigma) + xk3 * Cos(sigma) * Cos(2 * Msigma))

dw = (1 - xk3) * alpha * cosUn * xdw dww = Abs((dw - dw0)) * P0 Loop While dww > 0.0001 t = eps * sinUn2 / 4

xk1 = 1 + t * (1 - t * (3 - 5 * t + 11 * t * t) / 4) xk2 = t * (1 - t * (2 - t * (37 - 94 * t) / 8)) xx0 = Cos(sigma) * Cos(2 * Msigma)

xx0 = xk2 * (xx0 + xk2 * (1 + 2 * Cos(2 * sigma)) * Cos(3 * Msigma) / 6) / 4 Dsigma = xk2 * Sin(sigma) * (Cos(Msigma) + xx0) '计算椭球面距离S12

S12 = xk1 * b * (sigma - Dsigma) ' 计算计算椭球面方位角A1

p = Cos(u2) * Sin(w): q = cb1 - cb2 * Cos(w) A1 = Abs(Atn(p / q))

If p > 0 And q > 0 Then A1 = A1 If p > 0 And q < 0 Then A1 = pi - A1 If p < 0 And q < 0 Then A1 = pi + A1 If p < 0 And q > 0 Then A1 = 2 * pi - A1 '计算椭球面反方位角A2 p = Cos(u1) * Sin(w)

q = Cos(u1) * Cos(w) * Sin(u2) - Cos(u2) * Sin(u1)

A2 = Abs(Atn(p / q))

If p > 0 And q > 0 Then A2 = A2 If p > 0 And q < 0 Then A2 = pi - A2 If p < 0 And q < 0 Then A2 = pi + A2 If p < 0 And q > 0 Then A2 = 2 * pi - A2 If A1 > pi Then A2 = A2 - pi If A1 < pi Then A2 = A2 + pi

H = Cos(u1) * Sin(A1) - Cos(u2) * Sin(A2) '检核计算 End Select End Sub

4.4 大地主题解算程序功能介绍

大地主题解算的计算程序界面如下图5-1所示。在大地主题解算窗体上设置4个Frame控件，在Frame1内放置3个选择按钮option1(0)、option1(1)、option1(2)、option1(3)控制椭球参数；在Frame2内放置3个选择按钮option2(0)、option2(1)与option2(2)控制计算方法；在Frame3内放置2个选择按钮option3(0)与option3(1)控制计算方式；在Frame4内放置1个MSFlexGrid控件，4个CommandButton命令按钮，其Caption 属性分别为文件、添加、删除、保存，其功能主要用来编辑数据与保存数据；在大地主题解算窗体上设置3个CommandButton命令按钮，其Caption 属性分别为计算、显示、退出，其功能为计算。其计算界面直观、方便。计算的起算数据可通过文本编辑器事先编辑好通过磁盘文件打开，如果是正算则其文件名为BLAS.txt, 反算文件名为BLBL.txt,或者选择直接输入方式输入原始数据；其原始数据可存盘也可以不存盘不影响计算。计算结果直接存盘, 正算其文件名BLAS_BLA2.txt,反算结果文件名为BLBL_A1A2S.txt. 其计算结果单击显示按钮通过Word文档显示计算结果，然后单击退出按钮退回到主菜单界面。

图4－1大地主题解算的计算程序界面

五 子午线弧长正反解算 5.1 子午线弧长正算 1、数学模型

● 计算方法1 其中

为了便于计算, 常常将子午线弧长倍角计算公式转化为幂级数函数式,注意到式中 以后的各项最大不超过0.03mm, 小于投影计算所需精度0.5mm, 计算时只取到 项. 化算过程参见陈健、晁定波主编《椭球大地测量学》）. 其中

● 计算方法2（参见孔祥元等编著《大地测量学基础》） ， ， ， ，

5.2 子午线弧长反算

利用子午线弧长反算大地纬度，高斯投影坐标反算公式中要用到此项计算，反解公式可采迭代解法和直接解法。 1、 迭代解法

取初值

令 迭代计算

直到迭代计算满足 为止，以保证 的计算精度至 。在高斯投影坐标反算中由子午线弧长反算出的大地纬度称之为底点纬度。 2、 直接解法 其中

式中 、 、 由下式计算

上述算法参见陈健、晁定波主编《椭球大地测量学》。

六 高斯投影正反算和邻带换算 6.1 高斯投影正算

高斯投影正算：已知大地坐标 及指定中央子午线经度 ，计算高斯平面坐标 。计算公式如下：

其中 ， ， ， 。 的计算见子午线弧长正算公式。

6.2 高斯投影反算

高斯投影反算已知高斯平面直角坐标 及指定中央子午线经度 ，计算大地坐标 。计算公式如下：

式中 ， ，

底点纬度 的计算参见子午线弧长反算公式。

6.3 高斯投影邻带换算

为了限制投影变形，高斯投影采用分带投影的方法，但也带来了投影不连续的缺点。两相邻带的公共边缘子午线在两带平面上的描写形，弯曲方向相反，这样使得位于边缘子午线附近，分别属于两带的地形图不能拼接起来，为了解决这种困难，在相邻带之间设立重叠部分，在重叠部分的三角点要计算两带的坐标，这就产生坐标邻带换算的问题。此外当三角网跨越两个投影带，而平差计算是在一个带内进行时，也会遇到坐标换算问题。另外、在工程测量中常常使用3度带或1.5度或任意度带，而国家三角点成果是6度带，有时需要将6度坐标换算成3度、1.5度或任意度带坐标，或与此相反。 坐标邻带换算采用高斯投影正反算公式进行，坐标换算的实质是把椭球面上的大地坐标作为过渡坐标。已知 ，新与旧带轴子午线的经度 、 ，计算过程如下： 1) 由高斯投影反算公式计算： ； 2) ；

3) 由高斯投影正算公式计算：

6.4 高斯投影计算程序及其功能

高斯投影计算程序界面如下图6-1所示。在高斯投影计算窗体上设置3个Frame控件，在Frame1内放置3个选择按钮option1(0)、option1(1)、option1(2)控制椭球参数；在Frame2内放置3个选择按钮option2(0)、option2(1)与option2(2)控制计算方式；在Frame3内放置1个MSFlexGrid控件，3个CommandButton命令按钮，其Caption 属性分别为文件、添加、删除，其功能主要用来编辑数据；在高斯投影计算窗体上设置3个CommandButton命令按钮，其Caption 属性分别为计算、显示、退出。其计算界面直观、方便。计算的起算数据可通过文本编辑器事先编辑好通过磁盘文件打开，如果是正算则其文件名为BL.txt, 反算文件名为xy.txt,换代计算文件名为x1y1.txt，或者选择直接输入方式输入原始数据；计算结果直接存盘, 正算其文件名BL_xy.txt，反算结果文件名为xy_BL.txt，换代计算结果文件名文xyxy.txt. 其计算结果单击显示按钮通过Word文档显示计算结果，然后单击退出按钮退回到主菜单界面。

图6-1 高斯投影计算程序界面

● 高斯投影正反算计算子程序

Sub XY_BL(b As Double, l As Double, x As Double, y As Double, l0 As Double, K As Integer) Dim b11(1 To 4), r11(1 To 4), d11(1 To 4) Dim c As Double Dim m As Double

x2 = e2

X4 = x2 * x2 X6 = x2 * X4 X8 = X4 * X4 X10 = x2 * X8

B10 = 1# - Sqr(1# - x2) c = a / Sqr(1# - x2) EA = a * (1# - x2)

aa = 1# + 3# * x2 / 4# + 45# * X4 / 64# + 175# * X6 / 256# BB = 3# * x2 / 4# + 15# * X4 / 16# + 525# * X6 / 512# CC = 15# * X4 / 64# + 105# * X6 / 256# dd = 35# * X6 / 512#

aa = aa + 11025# * X8 / 16384# + 43659# * X10 / 65536# BB = BB + 2205# * X8 / 2048# + 72765# * X10 / 65536# CC = CC + 2205# * X8 / 4096# + 10395# * X10 / 16384# dd = dd + 315# * X8 / 2048# + 31185# * X10 / 13072# EE = 315# * X8 / 16384# + 3465# * X10 / 65536# FF = 693# * X10 / 131072# A1 = aa * EA

a2 = -BB * EA / 2# a3 = CC * EA / 4# a4 = -dd * EA / 6# A5 = EE * EA / 8# A6 = -FF * EA / 10#

R1 = 2# * a2 + 4# * a3 + 6# * a4 R2 = -8# * a3 - 32# * a4 R3 = 32# * a4 b11(1) = -a2 / A1 r11(1) = -a3 / A1 d11(1) = -a4 / A1 For i = 1 To 3 BIH = b11(1) + b11(1) * r11(i) - 2# * r11(1) * b11(i) b11(i + 1) = BIH - 3# * b11(i) * b11(i) * b11(1) / 2# r11(i + 1) = r11(1) + b11(1) * b11(i) DIH = d11(1) + b11(1) * r11(i) d11(i + 1) = DIH + b11(1) * b11(i) * b11(i) / 2# + 2# * r11(1) * b11(i) Next i

D0 = 1# / A1

D1 = 2# * b11(4) + 4# * r11(4) + 6# * d11(4) D2 = -8# * r11(4) - 32# * d11(4) D3 = 32# * d11(4) E1 = x2 / (1# - x2) Select Case K Case 1

BF = D0 * x t2 = Cos(BF): T3 = Sin(BF) BF = BF + t2 * T3 * (D1 + T3 * T3 * (D2 + D3 * T3 * T3)) T1 = Tan(BF): t2 = Cos(BF) N1 = E1 * t2 * t2 q = 1# + N1 N3 = c / Sqr(q) m1 = y / N3: m = m1 * m1 G = m / 30# * (61# + (90# + 45# * T1 * T1) * T1 * T1) NI10 = 3# - 9# * N1 b = BF - m * q * T1 / 2# * (1# - m / 12# * (5# + N1 + NI10 * T1 * T1 - G)) NI8 = 5# + (6# + 8# * T1 * T1) * N1 + (28# + 24# * T1 * T1) * T1 * T1 l = m1 / t2 * (1# - m / 6# * (q + 2# * T1 * T1 - m / 20# * NI8)) l = l + l0 Case 2

T1 = Tan(b): t2 = Cos(b): T3 = Sin(b) l1 = l - l0

N1 = E1 * t2 * t2

q = 1# + N1: N2 = c / Sqr(q)

x0 = A1 * b + t2 * T3 * (R1 + T3 * T3 * (R2 + R3 * T3 * T3)) m1 = t2 * l1 m = m1 * m1

NI4 = (4# * N1 + 5#) * q + m / 30# * (61# + 270 * N1 + (T1 * T1 _ - 58# - 330 * N1) * T1 * T1) x = x0 + m * N2 * T1 / 2# * (1# + m / 12# * (NI4 - T1 * T1))

NI5 = 5# + (T1 * T1 - 18#) * T1 * T1 - (58# * T1 * T1 - 14#) * N1 y = N2 * m1 * (1# + m / 6# * (q - T1 * T1 + m / 20# * NI5)) End Select End Sub

● 高斯投影换带计算子程序

Sub XY_BL(x1 As Double, y1 As Double, x2 As Double, y2 As Double, L01 As Double, _LO2 As Integer) DIM b, l As Double

Call XY_BL(b, l, x1, y1, LO1, 1) Call XY_BL(b, l, x2, y2, LO2, 2) End Sub

● 计算范例

1、已知参考椭球为克拉索夫斯基椭球 1）、 ， ， ； 2）、 ， ， ； 3）、 ， ， ； 4）、 ， ， ；

试计算 。 参考答案 1）、 , ； 2）、 , ； 3）、 , ； 4）、 , ；

2、设 点坐标 m ， m，求 点在 带的坐标？ 参考答案： 1）、计算 带轴子午线的经度，由 可知， 20， ； 2）、高斯投影反算： ， ，正算公式检核； 3）、计算 带轴子午线的经度， ，则 4）、 ； 5）、高斯投影正算： ， ，反算公式检核。

七 平面直角坐标系换算 众所周知，坐标系之间的差异主要取决于坐标系的定位与定向，椭球参数以及坐标系的尺度定义。新旧坐标系的换算是一个常见的测绘使用计算问题，一般旧坐标系是先前建立控制点坐标系，而新旧坐标系是后来建立的控制点的坐标系。新旧坐标系的换算方法与有严密方法和近似方法。从原理上讲，严密方法是将旧网的全部观测资料重新归算到新坐标系中，重新平差计算出个点的新坐标。而近似方法是在旧网原始观测资料不足或其它工程急需的情况下，常采用的一种方法。采用近似方法进行新旧坐标换算必须有足够的新旧坐标控制点，根据重合点的差值，按一定的规律修正旧网各点的坐标值，使旧网与新网达到最佳吻合。下面介绍几种常见的坐标换算的数学模型。

7.1 直接参数法

直接参数法就是利用两套坐标系两个已知公共点的坐标 、 、 、 求出坐标转换平移参数、尺度因子、旋转参数。 数学模型：

平移参数：

尺度因子： 旋转参数： 则其它点 的坐标转换公式：

上述方法是直接根据两公共坐标求转换参数，然后根据转换参数求坐标增量的转换值，最后求出转换点在新坐标系下的坐标。 7.2 相似变换（赫尔墨特法）

从理论上讲，就相似变换而言，二维坐标转换可以采用4参数模型，三维坐标可以采用7参数模型。二维坐标转换4参数模型具体表达如下

设 ， 表示新坐标的转换值， 表示新坐标的固定值， 表示旧坐标，即：

其中 平移参数， 为尺度比参数， 为旋转参数。公共点新坐标系转换坐标与固定之差为： 令 ， 。则上式可以写为

选取 个公共点，建立误差方程 ， 其法方程为

求出转换参数，按下式求任意点在新坐标系下的坐标

相似变换特点是不变更旧网的几何形状，将旧网整体平移，旋转尺度缩放配合到新坐标系中，其缺点在公共点有间隙存在，而且间隙可能还比较大，为了克服上述缺点，可以采用以下方法进一步逐步逼近。

7.3 正形变换法

经相似变换以后，坐标残差 、 被认为局部系统误差部分和偶然误差部分，其局部系统误差可以通过对网实施局部变形加以消除，这种局部变形必须满足正形条件， 根据第一次相似变换，然后根据转换值 与固定值 的差异（其残差）按正形变换理论进行第二次变换，其方法如下：

设 表示坐标系的系统部分， 表示偶然部分，则在公共点上第一次变换后

以残差 作为观测值，作为等权观测量，按间接平差列出误差方程式：

在上式求解中，共有10个未知数，至少需要5个公共点，经过上述处理，任一点正形变换的坐标公式为：

7.4 多项式逼近法

. 多项式逼近法在于选取多项式逼近待求的新旧坐标系统变换函数，由多项式逼近任意连续函数时，从理论上讲，选择适当的多项式阶数和系数，可以逼近到任意的程度，并保证点与点之间一一对应的可逆连续变换的特性。

设 、 表示旧坐标系的坐标， 、 表示新坐标系的坐标，多项式变换的一般公式为：

式中 为变换中心附近的一个原点坐标， 、 是待定系数，按 (2-7-23) 与（2-7-24）式求解待定系数至少六个公共点，在六个公共点上经变换不产生间隙。然后用所求的多项式系数变换其它点的坐标。若变换的区域较大，公共点较多，则可选择更高阶的多项式进行变换。若有 个公共点，可列出误差方程，用间接平差求解待定系数，误差方程如下：

式中 ，若观测量的坐标差 ，彼此独立等精度，则可分别求 的两个独立的法方程式：

求 的法方程阵系数完全一样，只是常数项中 换成 即可。最小二乘多项式拟合公共点经变换后坐标不再与固定坐标重合，而是保持公共点上的间隙平方和为最小。 7.5 算例与程序功能介绍

某城市GPS定位所到的北京54坐标、该市的独立坐标系坐标，采用\直接参数法\相似变换法\多项式拟合法\编程计算. 1、 已知数据与计算结果

表7-1 某城市北京54系与独立坐标系部分公共点坐标 点号 北京54坐标系坐标 独立坐标系坐标 X(m) Y(m) x(m) y(m)

1 2496657.820 508202.617 21931.951 117563.229 2 2495670.370 548381.370 20258.393 157719.614 3 2501866.993 486785.114 27506.217 96237.631 4 2517661.770 491937.908 43210.824 101659.490 5 2515160.425 486957.034 40794.893 96636.575 6 2520692.825 489367.657 46285.362 99141.362 7 2505545.925 510742.498 30775.459 120254.569 8 2516906.787 511422.913 42123.144 121128.932 9 2524169.605 522970.739 49187.733 132799.223 10 2490588.890 505566.376 15908.879 114823.693

表7-2 某城市北京54坐标→独立坐标转换结果 点号 直接参数法 相似变换法 多项式拟合法 X（M） Y（M） X（M） Y（M） X（M） Y（M）

1 * 21931.951 117563.229 * 21931.952 117563.229 *21931.951 2 * 20258.393 157719.614 * 20258.393 157719.614 *20258.393 3 27506.216 96237.631 * 27506.218 096237.630 *27506.217 4 43210.822 101659.490 * 43210.823 101659.490 *43210.824 5 40794.891 96636.576 * 40794.893 096636.576 *40794.893 6 46285.361 99141.360 46285.362 099141.361 *46285.362 7 30775.458 120254.569 30775.459 120254.569 30775.462 8 42123.143 121128.931 42123.144 121128.932 42123.149 9 49187.732 132799.223 49187.733 132799.224 49187.741 10 15908.879 114823.693 15908.879 114823.693 15908.875

新旧坐标转换计算成果 (相似变换法) ――――――――-----------------――

1.旧坐标系统坐标 ----------------------------------------------------------------------------- No Name x(m) y(m) 1 1 2496657.8200 508202.6170 2 2 2495670.3700 548381.3700 3 3 2501866.9930 486785.1140

117563.229 157719.614 96237.631 101659.49 96636.575 99141.362 120254.567 121128.932 132799.231 114823.697 4 4 2517661.7700 491937.9080 5 5 2515160.4250 486957.0340 6 6 2520692.8250 489367.6570 7 7 2505545.9250 510742.4980 8 8 2516906.7870 511422.9130 9 9 2524169.6050 522970.7390 10 10 2490588.8900 505566.3760 ------------------------------------------------------------------------------ 2.新坐标系统坐标

------------------------------------------------------------------------------ No Name x(m) y(m) 1 1 21931.9510 117563.2290 2 2 20258.3930 157719.6140 3 3 27506.2170 96237.6310 4 4 43210.8240 101659.4900 5 5 40794.8930 96636.5750 6 6 46285.3620 99141.3620 --------------------------------------------------------------------------

3.旧坐标--> 新坐标 ----------------------------------------------------------------------------

No Name x/y(m) x/y(m) v(mm) (Before) (After)

----------------------------------------------------------------------------

1 1 21931.9510 21931.95174 .74 117563.2290 117563.2290 .02 2 2 20258.3930 20258.39273 -.27 157719.6140 157719.6143 .30 3 3 27506.2170 27506.21759 .59 96237.6310 96237.6306 -.43 4 4 43210.8240 43210.82298 -1.02 101659.4900 101659.4902 .23 5 5 40794.8930 40794.89276 -.24 96636.5750 96636.5759 .92 6 6 46285.3620 46285.36221 .21 99141.3620 99141.3610 -1.04 7 7 30775.4589 120254.5691 8 8 42123.1438 121128.9319 9 9 49187.7329 132799.2238 10 10 15908.8793 114823.6931 ---------------------------------------------------------------------------

4.转换参数

--------------------------------------------------------------------------- X轴平移参数= -2465703.974 (m) Y轴平移参数= -433212.156 (m) 旋转参数= -0.58431(度.分秒) 尺度比参数= 8.92338e-6

------------------------------------------------------------------------------ 2、 程序界面及功能

坐标转换计算程序界面如下图7-1、图7-2所示。在坐标转换计算程序的窗体上设置1个TapStrip控件，在TapStrip控件的属性页的选项卡标添加的三个标题，分别为\坐标转换方法\、\旧坐标系统坐标\、\新坐标系统坐标\；设置3个Frame控件，在Frame1内放置3个选择按钮option1(0)、option1(1)、option1(2)控制坐标转换的计算方法；在Frame2内放置1个MSFlexGrid1控件，一个text1控件，3个CommandButton命令按钮，其Caption 属性分别为文件、添加、删除，用来编辑旧坐标数据；在Frame3内放置1个MSFlexGrid2控件，一个text2控件，3个CommandButton命令按钮，其Caption 属性分别为文件、添加、删除，主要用来编辑新坐标数据；在坐标转换计算程序的窗体上设置4个CommandButton命令按钮，其Caption 属性分别为保存、计算、显示、退出。起算数据新旧坐标可通过文本编辑器事先编辑好通过磁盘文件打开，其文件名分别为old.txt 与newxy.txt。或者选择直接输入方式输入新旧坐标原始数据；先存盘然后计算，计算结果直接存盘, 文件名为zhuanhuan.txt。计算结果单击显示按钮通过Word文档显示计算结果，然后单击退出按钮退回到主菜单界面。

3、 程序的主要功能按钮程序清单： Private Sub computer_Click()

If xyxyForm.option(0).Value = True Then Call computer1 If xyxyForm.option(1).Value = True Then Call computer2 If xyxyForm.option(2).Value = True Then Call computer3 view.Enabled = True: view.SetFocus End Sub

Private Sub computer1() '按直接法计算

Dim dx As Double, dy As Double, dx0 As Double, dy0 As Double Dim ddx As Double, ddy As Double Dim dis As Double, ddis As Double Dim xy1 As Double, xy2 As Double Dim a As Double, aa As Double On Error GoTo Line1

Open projectdir + \oldnum = 0

Do While Not EOF(1) oldnum = oldnum + 1 Input #1, nno$, xy1, xy2

oldname(oldnum) = nno$: oldx(oldnum) = xy1: oldy(oldnum) = xy2 Loop

Close #1

Open projectdir + \newnum = 0

Do While Not EOF(1)

newnum = newnum + 1 Input #1, nno$, xy1, xy2

newname(newnum) = nno$: newx(newnum) = xy1: newy(newnum) = xy2 Loop Close #1

'----------------寻找两公共点-------------- KK = 0

For J1 = 1 To 2

For j2 = 1 To oldnum

If KK = 0 And newname(J1) = oldname(j2) Then k1 = J1: i1 = j2: KK = KK + 1

nno$ = oldname(J1): oldname(J1) = oldname(j2): oldname(j2) = nno$ xy3 = oldx(J1): oldx(J1) = oldx(j2): oldx(j2) = xy3 xy3 = oldy(J1): oldy(J1) = oldy(j2): oldy(j2) = xy3 Exit For

ElseIf KK = 1 And newname(J1) = oldname(j2) Then k2 = J1: i2 = j2: KK = KK + 1

nno$ = oldname(J1): oldname(J1) = oldname(j2): oldname(j2) = nno$ xy3 = oldx(J1): oldx(J1) = oldx(j2): oldx(j2) = xy3 xy3 = oldy(J1): oldy(J1) = oldy(j2): oldy(j2) = xy3 Exit For End If Next j2

If KK = 2 Then Exit For Next J1

'计算转换参数

dx = oldx(2) - oldx(1): dy = oldy(2) - oldy(1)

ddx = newx(2) - newx(1): ddy = newy(2) - newy(1)

dis = Sqr(dx * dx + dy * dy): ddis = Sqr(ddx * ddx + ddy * ddy) a = qiua(dx, dy): aa = qiua(ddx, ddy)

q = a - aa: m = (ddis - dis) / ddis: m1 = m + 1 '计算转换坐标 For j = 1 To oldnum

dx = oldx(j) - oldx(1): dy = oldy(j) - oldy(1) ddx = Cos(q) * dx * m1 + Sin(q) * dy * m1 ddy = -Sin(q) * dx * m1 + Cos(q) * dy * m1 zhx(j) = newx(1) + ddx: zhy(j) = newy(1) + ddy Next j

Call save_result

Value = MsgBox(\坐标转换计算已完成!\坐标转换(\直接法)\

Exit Sub Line1:

msg = \错误：\

Value = MsgBox(msg, 32, \错误信息窗口\End Sub

Private Sub computer2() '按相似变换法计算 Dim IJ As Integer

Dim xy1 As Double, xy2 As Double, xy3 As Double ReDim c(100), w(10) On Error GoTo line1 '读取旧坐标文件数据

If projectdir = \Open projectdir + \ oldnum = 0

Do While Not EOF(1) oldnum = oldnum + 1 Input #1, nno$, xy1, xy2

oldname(oldnum) = nno$: oldx(oldnum) = xy1: oldy(oldnum) = xy2 Loop Close #1

'读取新坐标文件数据

Open projectdir + \ newnum = 0

Do While Not EOF(1) newnum = newnum + 1 Input #1, nno$, xy1, xy2

newname(newnum) = nno$: newx(newnum) = xy1: newy(newnum) = xy2 Loop Close #1

'根据新坐标调整旧坐标的顺序并确定公共点 For J1 = 1 To newnum For j2 = 1 To oldnum

If J1 <> j2 And newname(J1) = oldname(j2) Then

nno$ = oldname(J1): oldname(J1) = oldname(j2): oldname(j2) = nno$ xy3 = oldx(J1): oldx(J1) = oldx(j2): oldx(j2) = xy3 xy3 = oldy(J1): oldy(J1) = oldy(j2): oldy(j2) = xy3 Exit For End If Next j2 Next J1

'组成法方程 For i = 1 To 4

w(i) = 0 For j = i To 4

IJ = j * (j - 1) / 2 + i c(IJ) = 0 Next j: Next i

For i = 1 To newnum c(1) = c(1) + 1 c(2) = 0

c(3) = c(3) + 1

c(4) = c(4) + oldx(i) c(5) = c(5) + oldy(i)

c(6) = c(6) + oldx(i) * oldx(i) + oldy(i) * oldy(i) c(7) = c(5) c(8) = -c(4) c(9) = 0 c(10) = c(6)

w(1) = w(1) + newx(i) w(2) = w(2) + newy(i)

w(3) = w(3) + oldx(i) * newx(i) + oldy(i) * newy(i) w(4) = w(4) + oldy(i) * newx(i) - oldx(i) * newy(i) Next i

'解法方程系数矩阵求逆 nx = 4

ReDim xx(1 To nx)

Dim H(1 To 2000) As Double, h1 As Double, h2 As Double For k = nx To 1 Step -1 h1 = c(1): i1 = 1 For i = 2 To nx

i2 = i1: i1 = i1 + i h2 = c(i2 + 1) H(i) = h2 / h1

If i <= k Then H(i) = -H(i) J1 = i2 + 2 For j = J1 To i1

c(j - i) = c(j) + h2 * H(j - i2) Next j Next i

i2 = i2 - 1: c(i1) = 1 / h1 For i = 2 To nx c(i2 + i) = H(i) Next i Next k

'求转换未知数 For i = 1 To nx

xx(i) = 0

For j = 1 To nx If j > i Then

KK0 = j * (j - 1) / 2 + i Else

KK0 = i * (i - 1) / 2 + j End If

xx(i) = xx(i) + c(KK0) * w(j) Next j Next i

m = Sqr(xx(3) * xx(3) + xx(4) * xx(4)) - 1 q = Atn(xx(4) / xx(3)): q = qdms(q) For i = 1 To oldnum

zhx(i) = xx(1) + xx(3) * oldx(i) + xx(4) * oldy(i) zhy(i) = xx(2) + xx(3) * oldy(i) - xx(4) * oldx(i) Next i

Call save_result

Open projectdir + \

Print #1, \

Print #1, \ 4.转换参数 Print #1, \Print #1, \轴平移参数= \Print #1, \轴平移参数= \Print #1, \旋转参数= \度.分秒)\Print #1, \尺度比参数= \

Print #1, \Close #1

Value = MsgBox(\坐标转换计算已完成!\坐标转换(\相似变换法）\Exit Sub Line1:

msg = \错误：\

Value = MsgBox(msg, 32, \错误信息窗口\End Sub

Private Sub computer3() ' 多项式逼近法

Dim KK As Integer, IJ As Integer

Dim a(15) As Double, l(1 To 8) As Double, m(1 To 8) As Double Dim xy1 As Double, xy2 As Double, xy3 As Double Dim ll0 As Double, ll1 As Double ReDim c(200), w(15) Dim w1(15) As Double

Dim px As Double, py As Double On Error GoTo Line1

\ Open projectdir + \oldnum = 0

Do While Not EOF(1) oldnum = oldnum + 1 Input #1, nno$, xy1, xy2

oldname(oldnum) = nno$: oldx(oldnum) = xy1: oldy(oldnum) = xy2 Loop Close #1

Open projectdir + \newnum = 0

Do While Not EOF(1)

newnum = newnum + 1 Input #1, nno$, xy1, xy2

newname(newnum) = nno$: newx(newnum) = xy1: newy(newnum) = xy2 Loop Close #1

'---根据新坐标调整旧坐标的顺序并确定公共点--- For J1 = 1 To newnum For j2 = 1 To oldnum

If J1 <> j2 And newname(J1) = oldname(j2) Then

nno$ = oldname(J1): oldname(J1) = oldname(j2): oldname(j2) = nno$ xy3 = oldx(J1): oldx(J1) = oldx(j2): oldx(j2) = xy3 xy3 = oldy(J1): oldy(J1) = oldy(j2): oldy(j2) = xy3 Exit For End If Next j2 Next J1

'--多项式逼近组成法方程-- nx = 6

For i = 1 To nx

w(i) = 0#: w1(i) = 0# For j = i To nx

IJ = j * (j - 1) / 2 + i c(IJ) = 0# Next j Next i

'--求公共点重心坐标-- px = 0: py = 0

For i = 1 To newnum px = px + oldx(i) py = py + oldy(i) Next i

px = px / newnum: py = py / newnum '--求误差方程系数并组成法方程--

For i = 1 To newnum

a(1) = 1#：a(2) = oldx(i) - px： a(3) = oldy(i) - py

a(4) = a(2) * a(2)：a(5) = a(3) * a(3)：a(6) = a(2) * a(3) ll0 = oldx(i) - newx(i) ll1 = oldy(i) - newy(i) For J1 = 1 To nx

w(J1) = w(J1) + a(J1) * ll0 w1(J1) = w1(J1) + a(J1) * ll1 For j2 = J1 To nx

IJ = j2 * (j2 - 1) / 2 + J1 c(IJ) = c(IJ) + a(J1) * a(j2) Next j2 Next J1 Next i

'-法方程系数矩阵求逆- ReDim xx(1 To nx)

Dim H(1 To 200) As Double, h1 As Double, h2 As Double For k = nx To 1 Step -1 h1 = c(1): i1 = 1 For i = 2 To nx

i2 = i1: i1 = i1 + i h2 = c(i2 + 1) H(i) = h2 / h1

If i <= k Then H(i) = -H(i) J1 = i2 + 2 For j = J1 To i1

c(j - i) = c(j) + h2 * H(j - i2) Next j Next i

i2 = i2 - 1: c(i1) = 1 / h1 For i = 2 To nx c(i2 + i) = H(i) Next i Next k

'-求法方程未知数 For i = 1 To nx xx(i) = 0

For j = 1 To nx If j > i Then

IJ = j * (j - 1) / 2 + i Else

IJ = i * (i - 1) / 2 + j End If

xx(i) = xx(i) - c(IJ) * w(j)

Next j Next i

'---------------------------------- For i = 1 To nx l(i) = xx(i) Next i

'--------------------------------- For i = 1 To nx xx(i) = 0

For j = 1 To nx If j > i Then

IJ = j * (j - 1) / 2 + i Else

IJ = i * (i - 1) / 2 + j End If

xx(i) = xx(i) - c(IJ) * w1(j) Next j Next i

For i = 1 To nx m(i) = xx(i) Next i

' 求转换坐标

For i = 1 To oldnum

a(1) = 1#: a(2) = oldx(i) - px：a(3) = oldy(i) - py

a(4) = a(2) * a(2)：a(5) = a(3) * a(3)：a(6) = a(2) * a(3) zhx(i) = oldx(i) + l(1) + l(2) * a(2) + l(3) * a(3) zhx(i) = zhx(i) + l(4) * a(4) + l(5) * a(5)

zhy(i) = oldy(i) + m(1) + m(2) * a(2) + m(3) * a(3) zhy(i) = zhy(i) + m(4) * a(4) + m(5) * a(5) Next i

Call save_result

Value = MsgBox(\坐标转换计算已完成!\坐标转换(\多项式拟合法）\Exit Sub Line1:

msg = \错误：\

Value = MsgBox(msg, 32, \错误信息窗口\End Sub

Private Sub save_result()

Open projectdir + \

Print #1, \ 新旧坐标转换计算成果 \Print #1, \

Print #1, \ 1.旧坐标系统坐标 \Print #1, \

本文来源：https://www.bwwdw.com/article/kcgp.html