基于matlab的坐标转换 - 图文

河南理工大学本科毕业论文

本 科 毕 业 设 计（论文）

题目 基于MATLAB的坐标转换

院（系部）测绘与国土信息工程学院

专业名称 测绘工程 年级班级 10-5 学生姓名 刘威 指导教师 胡圣武

2014年 5月 25日

I

河南理工大学本科毕业论文

摘要

本文的目的是使用第四代编程软件MATLAB来实现大地测量学中坐标转换的程序编写。为实现国际测量数据的共享和我国不同时期下不同坐标系下数据成果的有效利用，坐标转换是一项具有现实意义的工作。而MATLAB软件界面友好，编程效率高，它以矩阵作为最基本的数据结构，具有强大而方便的矩阵转置、求逆、相乘等计算能力而且它能方便的进行大规模的数据处理，代码也十分简洁。这些优点使得它在工程计算中应用越来越广泛。因此使用MATLAB实现坐标转换乃至编制坐标转换系统具有工程应用上的意义。

本文先对MATLAB软件的发展历程、功能、特点、编程基本知识进行了介绍，然后叙述了大地测量学中坐标转换的基本概念，重点阐述了坐标转换的理论和方法及其数学模型。又对坐标转换中的空间直角坐标与大地坐标系之间的互换，高斯正反算以及七参数的空间直角转换模型进行了程序框图设计，并写出对应的MATLAB函数程序（附于附录中）。最后以珠穆朗玛峰峰顶坐标为实例进行计算并把计算结果与中海达坐标转换软件进行对比从而验证了文中空间直角坐标与大地坐标的计算的准确性，又对三个已知点坐标进行高斯正反算并将计算结果与中海达坐标转换软件进行了对比，验证了文中用MATLAB实现高斯正反算的计算的准确性和可靠性，并且以用文中提供的七参数计算程序，计算得出了已知WGS-84坐标和1954北京坐标系下三重合点坐标的从WGS-84坐标系向1954北京坐标系的七参数解。

关键字：MATLAB程序; 坐标转换; 程序框图; 七参数

II

河南理工大学本科毕业论文

ABSTRACT

The purpose of this article is to use the fourth generation programming software MATLAB to realize the program writing of coordinate transformation in geodesy. In order to achieve the sharing of international measurement data and effective using data results under different coordinate system during the different period of our country, coordinate transformation is a work full of realistic significance. The MATLAB is a Software with friendly interface, high efficiency, it takes matrix as the basic data structure, with a powerful and convenient matrix transpose, inverse, multiply computing capacity and it can be convenient for large-scale data processing, the code is also very concise. These advantages make it more and more widely applied in engineering calculation. So using MATLAB software to realize coordinate transformation or make coordinate conversion system is significant on engineering application.

In this article, first introduced the development, function, characteristic, basic knowledge of programming of MATLAB software, and then describes the basic concept of coordinate transformation in geodesy, in which expounds the theory and method of coordinate transformation and its mathematical model. The article displayed the program block diagram of the coordinate transformation between space rectangular coordinate and geodetic coordinate system, positive and negative gauss calculation and seven parameters of spatial orthogonal transformation model and writed the MATLAB function (attached to the appendix). As an actual example, the paper used Mt. Everest summit coordinate to calculate its geodetic coordinate and spatial orthogonal coordinate, and compared the result with Zhonghaida Software .In the positive and negative gauss calculation ,the paper calculated three known point and compared the result with Zhonghaida Software , proved the accurate and reliable of MATLAB. And the article gave out the seven parameters solution of WGS-84 coordinate system to 1954 Beijing coordinate system by three coincident point of them.

Keywords: MATLAB Procedures; Coordinate Transformation; Block Diagram

III

河南理工大学本科毕业论文

目录

第一章 引言 ................................................................................................................. 1 1.1绪论.......................................................................................................................... 1 1.2 研究现状 ................................................................................................................. 1 1.3 本文研究的主要内容 ............................................................................................. 2 第二章MATLAB软件的基本功能及特点 ................................................................ 4 2.1 MATLAB的发展历程及基本功能 ....................................................................... 4 2.2 MATLAB的优点 ................................................................................................... 5 2.3 MATLAB程序设计基础 ....................................................................................... 6 2.3.1 MATLAB程序设计的基本原则 ......................................................................... 7 2.3.2 MATLAB中的变量和常量 ................................................................................. 7 2.3.3矩阵运算基本操作 ............................................................................................... 7 第三章 大地测量坐标系统及其转换 ......................................................................... 9 3.1 地球、大地水准面及参考椭球面的基本定义及关系 ........................................ 9 3.2 大地测量基准 ....................................................................................................... 10 3.3 常用椭球坐标系及坐标形式 ............................................................................... 10 3.3.1地心坐标系 ......................................................................................................... 10 3.3.2参心坐标系 ......................................................................................................... 12 3.3.3常用坐标形式 ..................................................................................................... 14 3.4坐标系转换的模型及其公式表示 ....................................................................... 16 3.4.1同一参考椭球下不同坐标形式的转换 ............................................................. 17 3.4.2不同参考椭球下坐标的转换 ............................................................................. 20 第四章 坐标转换程序设计 ....................................................................................... 21 4.1坐标转换总思路设计 ........................................................................................... 21 4.2具体程序框图设计 ............................................................................................... 21 4.2.1同一椭球参数下大地坐标与空间直角坐标之间的转换框图设计 ................. 22 4.2.2同一椭球参数下大地坐标与高斯平面坐标的转换框图设计 ......................... 23 第五章 实例应用 ....................................................................................................... 24 5.1大地坐标与空间直角坐标间换算实例 ............................................................... 24 5.2高斯正反算实例 ................................................................................................... 27 5.3不同参考椭球基准下的坐标转换实例 ............................................................... 29

IV

河南理工大学本科毕业论文

第六章 结论与展望 ................................................................................................... 30 致谢.............................................................................................................................. 31 参考文献 ..................................................................................................................... 32 附录.............................................................................................................................. 33

V

河南理工大学本科毕业论文

1 引言

1.1绪论

随着全球经济和科技的飞速发展，越来越多的新技术新方法被应用于各个领域行业中，带给各领域新的发展契机和更为广阔的发展前景。作为与人们日常生活息息相关的测绘事业，也非常迫切地需要充分运用新的科技产品来更新测绘学科知识、提高测绘设备水平、改善测绘技能和方法，促进全球测绘一体化趋势的发展，从而更为有效地发挥测绘事业对社会、对国家、对人民大众服务的支撑作用。伴随大地测量领域的不断发展，坐标转换作为基础而重要的理论知识越来越多地被应用到规范和统一各类和各区域甚至全球测绘资料等相关方面。每个国家都会根据自身的发展情况和历史条件建立合适的测量坐标系统，这些测量坐标系基本都是基于地固坐标系统的参心坐标系或者地心坐标系，而且由于经济的发展或者习惯的变化，往往每个国家在不同的历史阶段也会采用不同的坐标系统，另外在实际的测量生产中，在某些未布设国家级大地控制网的地方或为了测量的方便建立不同于国家坐标系的地方独立坐标系。美国从1972年开始开发第二代卫星导航系统(即为现在通常所说的GPS)，它采用WGS-84地心坐标系统，凭借自身具有的全天候、连续实时、全球范围以及三维导航和定位等特点，在测绘领域很快被广泛采用，从而使测绘资料的全球一体化逐渐成为新的发展趋势。与此同时，坐标转换也成为一个越来越重要的课题，在测绘工作中发挥越来越重要的作用。

坐标转换本身是一个复杂的数值计算过程，倘若采用人工计算，必将增加计算的难度，费时费力且不能保证计算的精确度，因此对于坐标转换来讲，不但要研究更为严密的坐标转换方法，还要不断地借助技术更加优良的坐标转换工具来提高转换成果的质量以及简化计算的过程，从而使测绘资料得到更加有效地利用和统一。这不仅有利于我国的坐标系与国际通用坐标系相统一和共享，还有利于我国测绘事业的发展。 1.2 研究现状

由于坐标系的不同之处主要是坐标系的定位、定向以及尺度各异，因此若要保证坐标系间的转换精度，首先需要研究各种类型的坐标转换的理论方法。目前，测绘领域的坐标转换方法模型等研究已经相当成熟。不同参考椭球坐标转换的方法有直接参数法、多项式逼近法和相似变换法等。其中七参数法是坐标转换中的经典严密的常用坐标转换方法，我国通常采用的七参数数学模型为莫洛金斯基(Molodensky)模型、布尔莎(Bursa)模型和中国的武测模型等，这些转换模型已经广泛应用在各种工程测绘中。为确保坐标转换的质量，除了选取适当的转换模型之外，还需要不断地借助于一些新技术新工具，从而在提高坐

1

河南理工大学本科毕业论文

标转换精度的同时，也能够尽量简化程序编写过程。为此广大的测绘工作者利用计算机开发出许多内嵌于测量仪器或独立的坐标转换软件，并取得了很大成果。但由于采用的计算工具或方法的差异导致各种坐标转换产品转换的精度和操作的难易程度都不尽相同，甚至一些独立的小型坐标系统还无法保证正确的计算结果。而且，我国目前使用的坐标转换产品中很多都是国外软件，使用过程繁琐复杂且购买成本高。因此，需要开发一种性能更加良好的坐标转换工具来改善目前的情况。

目前，我国已正式启用最新的坐标系统—2000国家大地坐标(CGCS2000)，这就要求要把现阶段大量的1954年北京坐标系或1980年国家大地坐标系的成果转换到CGCS2000坐标系中，从而坐标转换再一次被提到了很重要的位置。今后，在常用坐标转换工具中也会出现各个坐标系向CGCS2000系坐标转换的模块。然而，我们目前运用坐标转换工具的编写大部分采用单一的VB、C++高级语言等，这些语言在编写程序时都存在着一些潜在的不足。例如：VB程序界面的编写操作非常简便、灵活，但在编写测量程序代码的过程中由于其内置的数学运算函数较少，导致程序编写很复杂，而且也给后续的程序调试工作带来很大的困难。因此，在刚刚启用CGCS2000坐标系的时候，一种编写程序简便、界面操作灵活的软件来编写坐标转换工具是非常必要的。 1.3本文研究的主要内容

本文的核心内容是用MATLAB数学计算软件解决坐标转问题。由于我国现有的大量测量成果大部分都涉及到1954年北京坐标系或1980年国家大地坐标系等参心坐标系与WGS-84坐标系的转换，其转换原理同CGCS2000坐标系和WGS-84坐标系间转换的原理基本一致，因此本文主要针对上述两种参心坐标系分别与WGS-84坐标系的坐标转换来研究。通过MATLAB数学运算软件在坐标转换方面的应用，认识到了MATLAB用于测绘行业的优势，并对这一技术应用到测绘领域的其它方面做了铺垫。论文的主要内容如下：

1．介绍了MATLAB的基本功能及特点，然后最后利用MATLAB软件编制了常用的坐标转换的一些模型函数（代码见附录）。

2．简要介绍了坐标转换方面的基本理论知识，包括地球椭球和参考椭球的定义、大地测量常用的坐标系；再详细介绍了不同参考椭球基准的坐标系转换和同一参考椭球基准常用坐标形式之间的换算过程和方法。

3．结合WGS-84坐标系和我国常用的1954年北京坐标系、1980年国家大地坐标系等，利用文中提出的MATLAB编制了坐标转换函数；最后以工程实例验证了MATLAB在解决坐标转换问题的有效性和可行性，并对本文进行总结和展望。

2

河南理工大学本科毕业论文

2 MATLAB软件的基本功能及特点

2.1 MATLAB的发展历程及基本功能

MATLAB的诞生源于对数值计算的需求。1980年，美国新墨西哥州大学计算机系主任Cleve Moler为了便于学生使用计算机计算，用Fortran语言编写出最初的MATLAB(即Matrix Laboratory的前三个字母)并受到广泛地欢迎和使用之后Moler继续对MATLAB进行开发和研究，并与一些数学家、软件专家合作成立Mathworks软件开发公司。MATLAB第一个商业版本在1984年问世；而后MATLAB通过不断地更新和添加不同的功能，如独特的符号运算、种类多样的绘图技术、多媒体应用以及与其他流行语言的接口功能等，逐渐发展壮大，8年之后，MATLAB4.0版本问世；经过很多功能的改进和增加，1999年该公司推出MATLAB5.0版本并随之推出了全新Simulink3.0版本，使MATLAB达到了新高度；之后又对操作界面做了很大的改观，并建立了程序发布窗口、变量管理窗口及历史信息窗口等，于2000年10月推出了更为便捷的MATLAB 6.0版本；2001年6月问世的MATLAB 6.1，其增加的虚拟显示工具箱能在三维实景下显示出更好的仿真结果；在界面设计不断完善的同时，MATLAB的核心数值算法和外部接口应用等功能得到了更多改善，并于2003年6月推出MATLAB Simulink 5.0版以及2004年9月推出MATLAB7.0，使MATLAB的功能更加庞大，之后每年均推出两个版本。MATLAB经过多年来的不断研究和完善，已经把最初只用于矩阵等数学计算发展为一种新型计算机编程高级语言，具有广泛的应用前景和推广价值。

MATLAB的功能主要包括MATLAB语言、MATLAB开发环境、MATLAB数学函数库、绘图系统和MATLAB应用程序接口(API)等几大部分。

1．MATLAB语言

MATLAB语言具有函数、程序流控制、输入输出、数据结构和面向对象编程等特征，是一种基于矩阵或者数组的高级计算机语言，尤其在大型编程方面其编写速度具备较大优势。

2．MATLAB开发环境

MATLAB函数和文件工具集就是MATLAB的开发环境。MATLAB的开发环境是一个集成的工作空间，方便人机交互式输入输出数据，而且具备对程序的编译和调试功能。它包括MATLAB桌面、执行命令窗口、历史命令窗口、M文件编辑调试器、MATLAB工作空间和在线帮助文档等。

3．数学函数库

MATLAB在解决数学计算中繁琐和困难的矩阵计算问题具备突出优势。它

3

河南理工大学本科毕业论文

以矩阵作为数据操作的基本单位，大大简化了矩阵计算，使其变得更加高效和方便。此外，MATLAB还提供了数量巨大种类丰富的数值计算函数库，可方便快速地解决多个领域的数值计算问题。

4．绘图技术

MATLAB的绘图技术主要是针对图形句柄的操作，操作种类分为高低两层。MATLAB可以十分方便地绘制各种图形，还可针对不同的需要来修饰和控制图形，具备强大的数值可视化共功能。

5．MATLAB应用程序接口(API)

MATLAB应用程序接口是一种函数库。该函数库能调用动态链接库(DLL)来完成MATLAB文件与C、Fortran等其他高级编程语言的数据交换，在MATLAB与其他应用程序间实现交互功能。

6．MATLAB工具箱(Toolbox)

工具箱(Toolbox)是MATLAB具备的一种特有的家族产品，用于解决不同领域的问题。它的实质是M文件和高级MATLAB语言函数库，可以使用户很容易地修改或者增加函数和代码。还可以使用户使用不同的工具箱来解决不同领域的问题。MATLAB工具箱一直在更新中，大致分为：应用类工具箱、控制类工具箱、信号处理类工具箱和其他常用工具箱。 2.2 MATLAB的优点

MATLAB与其他同类产品相比较，其在数值及符号运算、图形图像处理和可开发性等方面的功能更具优势。它具备简单易学、易操作、开放性、实用性强等优点，作为一种面向21世纪的科学计算语言，是目前科研和工程技术等众多领域的必备工具。

具体优点主要有： 1．编程易学且工作效率高

MATLAB不同于编译性语言，用户可以以一种更接近于书写计算公式的思维方式直接编写类似数学形式语言的程序。而且MATLAB程序文件还是一个扩展名为.m纯文本文件，若要对其进行编写和修改则可使用任何文字处理软件，同时也便于调试，有很强大的人机交互性。MATLAB是用C语言开发的，因此它的程序流控制语句同C语言极为相似，用户可以很容易就掌握它的使用方法。另外，MATLAB系统还有功能非常强大的帮助系统，能够以查询的方式帮助用户得到更多信息。MATLAB还有intro和demo等演示命令，为用户提供简单易懂的例子和演示。

2．平台独立性

MATLAB包含大量的平台独立功能，可以与Windows98/2000/NT系统和

4

河南理工大学本科毕业论文

UNIX的多种版本等很多操作系统兼容。同一个编写好的程序，在不同的平台上都可运行正常，同一个编写好的数据文件，在不同的平台上都可编译成功。这样极大地便于用户将MATLAB编写的程序和数据等移植到新平台使用，具有很强的平台独立性。

3．预定义函数

MATLAB自身含有一个用来解决基本工程问题的大型预定义函数库，函数库中的函数有成百上千个且均可直接使用。它可以避免输入例如数组下标、统计量等的诸多麻烦，使编程变成一件非常简单的事情。除此以外，针对一些重要领域的复杂问题，MATLAB还加载了多种专用工具箱来方便地解决。这些工具箱涉及的方面有信号处理、通信、控制系统、神经网络、图像处理和其它领域的诸多相关问题。

4．数学运算方面

除了拥有超强优势的数值计算功能，MATLAB还可凭借自己的符号运算功能帮助用户轻松地解决繁杂数学运算的分析问题，比如矩阵变换及运算、微积分运算、线性与非线性方程求解、插值求解等问题。另外，MATLAB的函数都具备算法先进的自适应能力，弥补了非可执行文件的MATLAB程序运行速度欠佳的不足，从而解决了计算机对算法的选择难题。

5．机制独立的画图

MATLAB作为一种可视化数据技术的卓越工具，有着其他语言所没有的许多画图或图像处理命令。MATLAB运行这些命令之后，相应的二维或者三维图形或图像就会显示出来。

6．MATLAB编译器

MATLAB编译器先将MATLAB程序编译成某种独立性程序，再以解释的方式运行该程序，很好地体现出MATLAB的平台独立性能。

7．可扩充性和可开发性

MATLAB具有可扩充和可开发功能，这些功能保障了MATLAB自身能不断地发展和完善。MATLAB相当于一个解释系统，它以一种解释执行的方式来运行函数程序，这种系统最大的特点就是MATLAB具备一个开放的系统，便于用户查看和添加程序代码。 2.3 MATLAB程序设计基础

2.3.1 MATLAB程序设计的基本原则

突破以往其他程序语言经常采用的循环思想，尽量用MATLAB矩阵式语言书写程序，使得程序简洁，执行效率高。在程序设计中尽量避免重复的脚本代码，多用MATLAB提供的函数。系统中的函数要比用一般代码编的函数执行效率高很

5

河南理工大学本科毕业论文

多。在编写比较大的程序时，应该对各个细节以函数或子过程方式处理，避免矩阵混淆。

在程序编制过程中，各个功能部分尽量封装在函数中，这样不但可以减少全局变量个数，而且对各个函数的修改要比对整个程序的修改方便得多。 2.3.2 MATLAB中的变量和常量

在MATLAB中，变量名可由字母A-Z、a-z、数字和下划线“_”组成，但第一个字符必须是字母。

注意：MATLAB是区分大小写字母的，如矩阵a和A是不一样的。 在变量使用之前，用户不需要指定一个变量的数据类型，也不必声明变量。MATLAB有许多不同的数据类型，这对决定变量的大小和形式是有价值的，特别适合于混合数据类型、矩阵、细胞矩阵、结构和对象。

变量有局部变量和全局变量两种。

局部变量（local）是存在于函数空间内部的中间变量，产生于该函数的运行过程中，其影响范围也仅限于函数本身。

全局变量（global）是在不同的工作空间以及基本工作空间中可以被共享的变量。必须用global逐个对具体变量加以专门定义，没有global定义的函数和基本空间，将无权享用全局变量。 2.3.3矩阵运算基本操作

MATLAB中最基本的数据结构是矩阵，矩阵的二维数据结构，能非常容易就成倍的存取数据元素，数据元素可以是数字、字符、逻辑真假或其他类型的MATLAB结构。MATLAB采用这种二维的矩阵存取单个的数，用1乘1的维数表示，也存储向量，用1乘n的维数表示向量的长度。下面介绍矩阵有关的操作。

1. 创建矩阵

在MATLAB中，建立矩阵最简单的方法，是利用矩阵构造操作符：方括号 []。在方括号中写入元素，元素之间用空格或逗号隔开，能建立矩阵的一行。2. 连接矩阵

C=[A B]是横向连接矩阵A和B，连接矩阵最简单的方法就是使用方括号[]。

要求A与B有相同的行数；C=[A;B]是纵向连接矩阵A、B，要求A和B有相同的列数。 3. 重置矩阵形状

获取矩阵的形状与大小信息，经常使用length、size和ndims 4个函数。

6

河南理工大学本科毕业论文

3 大地测量坐标系统及其转换

3.1 地球、大地水准面及参考椭球面的基本定义及关系

地球表面是一个由地球陆地和海洋所构成的变化异常、褶皱不平的似椭球面，故不便于作为测量基准。由于海洋约占全球面积的71％，故设想用海平面表述地球形体，其定义是：假定海洋的水平面在完全平衡和静止的状态，没有风浪潮汐等自然因素的影响，海洋的表面以及由它延长到大陆的下面，并处处保持与垂线方向相交成直角，具有这一特征的闭合形状的面，称为大地水准面。由它所包围的整个地球体，叫做大地体。由于地球形状复杂、质量分布不均，导致大地水准面实际上仍然是一个很不规则的闭合曲面，无法将其简单地表示为一个模型或者公式。实际应用中，一般选择某一平均海水面代替大地水准面，但是平均海水面相对地球重力等位面而言并不是等位面，它同样也是个褶皱不平的闭合曲面，在海洋学中被叫做海面地形。由此可知，不同的验潮站所推算的平均海水面必然各不相同，同样地，每个国家选定作为高程基准面的平均海水面也不尽相同，由此导致全球大地水准面不具备统一性，从而增加了大地测量工作的难度。

长期的理论研究和观测表明两极扁平的椭球体较为接近真实地球的几何形状和物理形态。我们把由一个通过南、北极的子午圈绕地球南北极轴旋转一周而成的椭球叫做旋转椭球，用来代表地球的旋转椭球叫做地球椭球，用来代表某一地区大地水准面且具有一定几何参数、定位及定向的地球椭球叫做参考椭球。参考椭球面被看作大地测量基准面，也被看作地球形状和地图投影研究的参考面。地球表面、大地水准面和椭球面的关系如图3．1所示。

图3-1 地球表面、大地水准面、似大地水准面之间的关系

3.2 大地测量基准

基准指的是为了唯一确定空间的某个位置或形态而选取或设定的参照物，大地测量基准指的是为了确定地球形状的参考椭球。该基准用参考椭球的长半轴、短半轴、参考椭球的定位和定向等参数表示。在大地测量中，不同的坐标

7

河南理工大学本科毕业论文

系具有不同的基准。例如地心坐标系和参心坐标系就是根据选取的原点不同选用了不同的参考椭球。基准面指的是选取最符合某一区域地球表面的参考椭球面。每个地区均有不同的参考椭球，因此也就有不同的基准面，比如我国常用的1954年北京坐标系、1980年国家大地坐标系椭球面就是两个不同的基准面。目前GPS卫星发布的广播星历的坐标参考基准是WGS-84世界大地坐标系，该坐标系基准面就是WGS-84椭球面。 3.3 常用椭球坐标系及坐标形式

坐标是用来表示一个点在相对参照系下的确定位置，而坐标系是表示用坐标来确定点位的方法，被看作测量参照系的数学依据。坐标系分为多种类型，比如大地极坐标系、大地坐标系和坐标轴相互正交的笛卡儿坐标系等。单独的用坐标系还不能将点的位置确定下来，需要再引入参考基准即参照系的概念，从而形成固定严密的坐标参照系来唯一地确定点的位置。

由于点位的表示形式会随着所选择的坐标系的改变而改变。甚至在同一种坐标系下，其坐标形式也有所不同。不管如何不同，各种坐标形式之间必然可进行某种相互转换，即坐标转换。坐标系统之间的坐标转换有不同参考椭球间的坐标系转换，比如参心坐标系间、地心坐标系间或者参心与地心坐标系间，也有同一参考椭球坐标系下的空间直角坐标与大地坐标间、大地坐标与高斯平面坐标间的坐标换算等。 3.3.1地心坐标系

地心坐标系是以地球质心为原点的椭球坐标系，一般有空间直角坐标和大地坐标两种常用坐标形式。地心空间直角坐标系的X轴指向格林尼治平均子午面与地球赤道面的交点，Z轴指向地球北极，Y轴过原点垂直于平面XOZ并与其他两轴构成右手空间直角坐标系。地心大地坐标系的大地纬度B是椭球面某点的椭球法线与赤道面的夹角，大地经度L是椭球面某点的子午面与格林尼治大地子午面之间的夹角，大地高H是某点沿椭球法线到椭球面的最短距离。

1. WGS-84世界大地坐标系

WGS-84坐标系是目前最为常用的地心坐标系，它最初由美国国防部根据 TRANSIT导航卫星系统的多普勒观测数据建立，从1987年1月开始作为GPS卫星所发布广播星历的坐标参照基准。它的三轴具体指向为：

X轴指向BIH1984.0的起始子午线和赤道面的交点； Z轴指向国际时间局(BIH)1984．O定义的协议地球极方向； Y轴和Z、X轴构成右手坐标系。 WGS-84椭球体的四个主要参数如下： 长半轴a=6378137 m；

8

河南理工大学本科毕业论文

地球引力常数(含大气层)：GM?3986005?108m3s?2

?6正常化二阶带球谐系数：C2.0??484.16685?10

地球自转角速度：??7292115?1011rad/s

根据以上4个参数可以进一步求得：地球扁率、第一偏心率、第二偏心率、赤道正常重力和极正常重力。

2. 2000国家大地坐标系

我国目前实际使用的两个大地坐标参照系—1954北京坐标系和1980西安坐标系，都属于参心系，它们都是采用传统地面测量技术建立起来的，并满足了当时实际应用的需求。但随着时代变迁和科学技术的发展，特别是空间技术的发展，一方面，越来越多的实际应用要求建立和采用地心系，另一方面，空间定位技术的发展，也是的建立地心系成为可能。为顺应这一趋势，我国提出了2000国家大地坐标系（CGCS 2000—China Geodetic Coordinate System 2000）。2000国家大地坐标系的定义如下：

? 原点:包括海洋和大气在内的整个地球的质心。

? 长度单位：国际单位制的米，与局部地心框架下的地心坐标时一致，通

过适当的相对论模型获得。

? 定向：初始定向与1984.0时的BIH（国际时间局）定向给定。 ? 定向的时变：定向的时变不产生相对于地壳的参与全球旋转。 ? CGCS 2000大地坐标系是右手地固直角坐标系。原点在地心，Z轴与IERS

参考极（IRP）方向一致，X轴为IERS参考子午面（IRM）与垂直于Z轴的赤道面的交线，Y轴与Z轴和X轴垂直并最终构成右手正交坐标系。

CGCS 2000的参考历元为2000.0。

参考椭球采用2000参考椭球，其相关常数定义为：

a?6378137m;

GM?3.986004418?1014m3?s?2; J2?1.082629832258?10?3;

??7.292115?10?5rad?s?1。 正常椭球与参考椭球一致。

CGCS 2000由一下三个层次的站网坐标和速度具体实现：

1) 第一层次为连续运行参考站。由它们构成CGCS 2000的基本骨架，其坐标精度为毫秒级，速度精度为1mm/年。

2) 第二层次为大地控制网。包括中国全部领土和领海内的高精度GPS网

9

河南理工大学本科毕业论文

点，其三维地心坐标精度为厘米级，速度精度为2-3mm/年。

3) 第三层次为天文大地网。包括经空间网与地面网联合平差的约5万个天文大地点，其大地经纬度误差不超过0.3m，大地高误差不超过0.5m。 3.3.2参心坐标系

参心坐标系的建立要素：(1)椭球定位(确定椭球中心的位置)；(2)椭球定向 (确定椭球短轴的指向)；(3)椭球大小和形状 (确定椭球的几何参数)；(4)建立大地原点。选定某一合适的点作为大地原点，并通过该点来进行精密天文大地测量和高程测量，对参考椭球进行定位和定向。这种坐标系的确定方法带有很大的地区适用性，建立的参考椭球中心一般不会重合于地球质心，因此被称之为参心坐标系。参心坐标系和地心坐标系的定义相似，也有空间直角坐标系和大地坐标系两种常用坐标形式。我国常用的参心坐标系有：1954年北京坐标系、1980年国家大地坐标系(1980西安坐标系)和新1954年北京坐标系、WGS-84坐标系以及2000国家大地坐标系，下面分别简要说明。

1．1954年北京坐标系

新中国建立后，我国的大地测量事业进入了全面发展时期，在全国范围内开展了正规的、全面的大地测量和测图工作，迫切地需要建立一个参心大地坐标系。受限于当时的条件，我国选取克拉索夫斯基椭球作为参考基准，选取前苏联境内普尔科沃天文台作为坐标起算点并与前苏联1942年坐标系联测，通过计算建立了我国大地坐标系，定名为1954年北京坐标系。高程异常参考前苏联1955年大地水准面差距的重新平差结果为依据，按我国的天文水准路线计算过来的。之后我国以1954年北京坐标系为基准建立了全国天文大地网并取得大量测绘成果。但是该坐标系存在的很多缺点也随着全球测绘理论及技术的不断研究和发展逐渐凸显：

1) 椭球参数有较大的误差。克拉索夫斯基椭球的长半轴比目前国际精确计算出的椭球约长109米；

2) 物理与几何大地测量各自的基准面不一致。我国处理重力数据采用的是赫尔默特1900-1909年正常重力公式，而与这个公式相应的赫尔默特椭球扁球不是旋转椭球，因而存在一定错误；

3) 我国大地水准面与克拉索夫斯基椭球面有明显自西向东的系统性倾斜，东部地区的大地水准面差距甚至可达+68米。这种现象导致大比例尺地图不能准确地表示地面实际情况，对观测元素的归算的要求也更为严格；

4) 定向不够明确。椭球短轴既不指向我国的地极原点JYD1968.0；也不指向国际上广泛采用的国际协议原点CIO(Conventional International Origin)：起始大地子午面与国际时间局BIH(Bureau International de I Heure)所定义的格林尼治

10

河南理工大学本科毕业论文

平均天文台子午面也稍有不同。

另外，因为1954年北京坐标系的大地点成果是由局部平差得来，所以肯定会产生一定误差和矛盾。

2．1980年国家大地坐标系

为了克服1954年北京坐标系的缺点建立了1980年国家大地坐标系(GDZ80)。该坐标系的特点是：

1) 用1975年国际大地测量与地球物理联合会(IUGG)第16届大会推荐的4个椭球基本参数。

地球椭球长半径a=6378140m；

地球重力场二阶带球谐系数J2?1.08263?10-3； 地心引力常数GM?3.986005?1014?1014m3/s2； 地球自转角速度??7.292115?10?5rad/s。

根据物理大地测量学中有关公式，可由上述4个参数算得。

地球椭球扁率??1/298.257；

2赤道正常重力值?0?9.78032m/s。

2) 心大地坐标系的建立基础是1954年北京坐标系；

3) 地原点位于我国陕西省西安市以北60km处的泾阳县永乐镇，简称为西安原点；

4) 向明确。椭球短轴与地球质心指向地极原点JYD1968.0的方向平行，起始大地子午面与我国起始天文子午面平行；

5) 点定位，使椭球面和似大地水准面最为密合于我国境内； 6) 取1956年黄海高程基准为大地高程基准。

该坐标系和1954年北京坐标系的不同之处除了它们各自的不同参考椭球和不同椭球定位、定向外，还有前者经过整体平差，而后者只做了局部平差。

3．新1954年北京坐标系

新1954年北京坐标系，简称BJ54新，原1954年北京坐标系又称为旧1954年北京坐标系BJ54旧。由于在全国的以GDZ80为基准的测绘成果建立之前，

BJ54旧的测绘成果仍将存在较长的时间，而BJ54旧与GDZ80两者之间差距较

大，给成果的使用带来不便，所以又建立了BJ54新作为过渡坐标系。经过过渡坐标系的转换，BJZBJ54新?ZGDZ80??Z0和BJ54旧的控制点的高斯平面坐标，其差值在全国80%地区内小于5m，局部地区最大达到12.9m，这种差值反映在1:5万以及更小的比例尺的地形图上的影响，图上位移绝大部分不超过0.1mm。这样采用BJ54新，对于小比例尺地形图可认为不受影响，在完全采用GDZ80测绘成果之后，1:5万以下的小比例尺地形图不必重新绘制。

11

河南理工大学本科毕业论文

北京54系的特点是：

1) 克拉索夫斯基参考椭球为基准；

2) 点定位，但椭球面与大地水准面不最佳拟合于我国境内； 3) 地起算数据与GDZ80不一致，但大地原点一致；

4) 向明确，坐标轴平行于GDZ80，椭球短轴平行于地球质心指向1968.0地极原点JYD1968.0的方向，起始大地子午面平行于我国起始天文子午面

?X??Y??Z?0；

5) 用1956年黄海高程系作为大地高程基准；

6) BJ54旧相比，所采用的椭球参数相同，其定位相近，但定向不同。BJ54旧的坐标是局部平差结果，而BJ54新是GDZ80整体平差结果的转换值，两者之间无全国统一的转换参数，只能进行局部转换。

4．地方独立坐标系

我国的平面坐标采用高斯投影。高斯投影之后，除中央子午线外的其他线段，投影后都会有长度变形，距中央子午线越远变形程度越大。我国采用6?带或3?带的国家分带投影来减小这种变形。然而在城市、工矿等相对较小区域的实际测量中，如果控制网根据国家分带坐标系来建立，就可能导致地面长度投影产生较大的变形。实践表明，长度变形若大于2.5 cm/km，就达不到某些工程测量精度要求。解决的办法是建立一个地方独立坐标系：将当地平均海拔高程面作为测量控制网高程基准，再选取适当的当地子午线作为中央子午线来缩小对其高斯投影之后的误差。该椭球的中心、扁率、轴向同国家参考椭球一致但长半轴不同，被叫做地方参考椭球。因为我国常用的坐标系都是参心坐标系，所以将地方独立坐标系也划分在参心坐标系范围内。 3.3.3常用坐标形式

不论是参心坐标系还是地心坐标系，都可以将坐标表示为空间直角坐标、大地坐标和经过平面投影之后的平面直角坐标。

1．空间直角坐标

以椭球中心为原点O，X轴为起始子午面与赤道面的交线，Z轴为椭球的短轴，Y轴与X轴正交于赤道面内，三轴相互垂直构成右手坐标系O—XYZ，如图3．2所示。在该坐标系中，A点的位置用(X，Y，Z)表示。任意选取空间直角坐标系的两轴就可确定一个坐标平面，共可确定出三个坐标平面，它们相互垂直，将空间分为八个象限。用OXY平面将空间分为上下两个部分。上面的空间，由X，Y，Z三轴的正半轴所组成的空间为第一象限，从第一象限开始，按逆时针方向依次称为第二、三、四象限，第一象限下面的象限为第五象限，同样按逆时针方向依次将下面的空间分为称为第五、六、七、八象限四个部分。

12

河南理工大学本科毕业论文

图 3.2 空间直角坐标系

2．大地坐标

大地坐标系中的点与空间直角坐标系中的点一一对应，如图3.3所示，只有坐标的表示形式不同。参考椭球面上点P的子午面与起始子午面所构成的二面角L，叫做P点的大地经度。它将起始子午面看作0?，向东为东经(0?～180?)；向西为西经(0?～?180?)。P点的法线与赤道面的夹角B叫做P点的大地纬度，它将赤道看作0?，向北为北纬(0?～90?)；向南为南纬(0?～?90?)。如果点不在椭球面上，还需增加一个参数：大地高H大地，它与正高H正及正常高H常的关系为：

H大地?H正?N H大地?H常??

式中：N表示大地水准面差距，?表示高程异常。在大地坐标系中，P点的位置用(B，L，H)表示。

图3.3 大地坐标系

3．平面直角坐标系

一般来讲，测量用的平面直角坐标系是指在二维平面内用一条垂直方向的X轴和一条水平方向的Y轴分割而成的坐标系，其中X轴正方向指向北而Y轴

13

河南理工大学本科毕业论文

正方向指向东。常见的大范围地形图及测量定位等往往都采用平面直角坐标系，该坐标系是通过椭球投影转换将一种基于参考椭球的大地坐标换算到平面直角坐标系上。由于地球椭球曲面无法准确地展成平面，因此在换算的过程中必定会产生长度、面积和角度等方面的变形。

我国测量生产中广泛采用高斯一克吕格正形投影(简称高斯投影)，它是一种等角(正形)投影，即角度无变形，同一位置点的任意方向的长度比相同，而不同位置点的长度比不同。高斯投影也是一种横轴切椭圆柱投影，该投影被看作是在地球椭球的外表面横套一个椭圆柱面，椭圆柱的中心通过椭球的中心，并与椭球相切于椭球的一条子午线(该子午线叫做中央子午线)，将椭球上的点投影到该椭圆柱后，将椭球柱展开，形成以中央子午线投影作为X轴(位于X轴上的点无长度变形)，以赤道在椭圆柱上的投影作为Y轴的高斯平面直角坐标系。

3.4坐标系转换的模型及其公式表示

目前，GPS卫星定位系统已经作为一种高效和便捷的测量工具被广泛使用，但是此系统的参考基准是WGS-84坐标系，与很多国家采用的基准不同。如果要使用GPS的测量成果，就必须把这些成果数据转换到各自建立的坐标系中；另外在同一国家坐标系基准下，不同的坐标形式也需要相互换算，由此就产生了同一参考椭球基准的不同坐标形式换算和不同参考椭球基准的坐标转换等问题。若要实现多种坐标系下的数据共享，就必须解决坐标系转换问题。 3.4.1同一参考椭球下不同坐标形式的转换

1. 空间直角坐标向大地坐标换算

??Y?L?arctan????X?????Z?N?H???? (3.1) ?B?arctan?222???X?YN1?e?H??ZX2?Y2?2H??N1?e或?N?sinBcosB???????

式中：参考椭球参数长轴为a，短轴为b，空间直角坐标为?X,Y,Z?，大地坐标为?B,L,H?， N为卯酉圈曲率半径且N?a1?esinB22 ，e为椭球的第一偏

a2?b2心率，且 e? 。

a

14

河南理工大学本科毕业论文

式中：B的求取需要通过迭代的方法，一般取迭代初值B?arctanZ?X2?Y2?，

再代入式3.1中进行计算直至最终B?的值与B的值小于某一允许值。

2.大地坐标向空间直角坐标换算

图3.4 大地坐标与空间直角坐标之间的换算

??X??N?H?cosBcosL?? (3.2) ?Y?(N?H)cosBsinL?2a2?Z?[N(1?e)?H]sinB?[N?H]sinB?b2?

式3.2中空间直角坐标为?X,Y,Z?，大地坐标为?B,L,H?，N为卯酉圈曲率半径且N?a1?esinB22 ，e为椭球的第一偏心率，且 e?a2?b2。 a3.大地坐标向高斯平面坐标的换算

图 3.5 高斯（横轴切椭圆柱）投影

15

河南理工大学本科毕业论文

图3.6 投影换带计算

1. 高斯正算公式

11x?X?t?N?cos2B?l2?t?N?5?t2?9?2?4?4cos4B?l4224(3.3)

1?t?N?61?58t2?t4?270?2?330?2t2?cos6B?l6720????y?N?cosB?l??11N?1?t2??2?cos3B?l3?N?6120 (3.4)

5?18t2?t4?14?2?58?2t2?cos5B?l5???高斯正算是由大地坐标?B,L?求取高斯平面上投影点坐标?x,y?，，在高斯正

??L?L0???e2cosB，算公式中l?，L0 表示中央子午线的经度，而t?tanB，

???求取x时关键在于子午线长X的计算，其数值计算公式列举如下：

X?a?1?e1??A0?B?A2?sin?2B??A4?sin?6B??A8?sin?8B?? （3.5）

2??3e45e1350e111025e1 式中： A0?1?1???4645121638424683e160e1525e117640e1?1??； A2???????2?46451216384??46815e1210e18820e1?1??； A4?????4?6451216384??6835e12520e1?1??； A6?????6?51216384??24681315e1； A8??8163848其中A0?A8为计算中涉及参数。 2.高斯反算公式

16

河南理工大学本科毕业论文

22Vftf??y???15?3tf??f2?9?f2tf2??B?Bf?2??N?12??f????y??N?f???161?90tf2?45tf4?360?4???y??N?f?y??N?f????6??（3.6） ????（3.7） ??1l?cosBf??y??N????f?1??1?2tf2??f2?6????y??N?f???15?28tf2?24tf4?6?f2?8?f2tf2?120?3??????5高斯反算公式是已知某点投影在高斯平面上的坐标?x,y?，要求该点在椭球上的大地坐标?B,L?，公式中tf?tan?Bf? ，Nf??1?e??1?e21a22?cos?Bf??，

Mf?a1?e12 ，?f?e2?cos?Bf? ，Vf?1?e22?cos2?Bf?。很显然想要求

出大地坐标，其关键在于底点纬度Bf的求解。以下为底点纬度数值计算公式，以及由它推导出来的?B,L?数值计算公式即电算公式；

电算时Bf的求解见式3.8：

?3.8）Bf???50221746?293622?2350?22cos2?cos2?cos2?10?10sin?cos???（

由Bf进一步计算?B,L?：

?B?Bf?1?b4?0.12Z2Z2Z2b2?????22 ?l?1?b3?b5ZZZ??? （3.9）

?L?L?l0????????

????????式（3.9）中， Z?y，

NfcosBf??x???；

6367558.4969b2??0.5?0.003369cos2Bf?sinBfcosBf； b3?0.333333??0.166667?0.001123cos2Bf?cos2Bf；

b4?0.25??0.16161?0.00562cos2Bf?cos2Bf； b5?0.2??0.1667?0.0088cos2Bf?cos2Bf；

而b1?b5均为求解参数。 3.4.2不同参考椭球下坐标的转换

不同坐标系统的转换本质上是不同基准间的转换，不同基准间的转换方法有很多，其中，最为常用的有布尔沙模型，又称为七参数转换法。

17

河南理工大学本科毕业论文

七参数转换法是：

设两空间直角坐标系间有七个转换参数：3个平移参数、3个旋转参数和1个尺度参数。

ZAZBXCBOA?X0OBXB?xXAYB?ZXCAC?YYA

图3.2 七参数转换

若：

(XA,YA,ZA)T为某点在空间直角坐标系A的坐标； (XB,YB,ZB)T为该点在空间直角坐标系B的坐标；

(?X0,?Y0,?Z0)T为空间直角坐标系A转换到空间直角坐标系B的平移

参数；

(?X,?Y,?Z)为空间直角坐标系A转换到空间直角坐标系B的旋转参

数；

m为空间直角坐标系A转换到空间直角坐标系B的尺度参数。

则由空间直角坐标系A到空间直角坐标系B的转换关系为：

?XB???X0??XA??Y????Y??(1?m)R???R???R????Y?

XYZ?A??B??0????ZB????ZA????Z0??

(3.10)

式中：

0?1?R??X???0cos?X?0?sin?X???sin?X? (3.11) cos?X??0 18

河南理工大学本科毕业论文

?cos?Y?R??Y???0?sin?Y??cos?Z?R??Z????sin?Z0?sin?Y??10? (3.12) 0cos?Y??sin?Zcos?Z0??0? (3.13) ??001??

一般?X、?Y和?Z均为小角度，可以认为：

cos??1sin???

则有：

??ZR(?)?R(?)?R(??1ZY)?R(?X)????1?Z??Y??X也可将转换公式表示为：

???X?B?????X0????XA????X? ?Y??Y?B????Y0???Y?A??K?? ?ZB?????Z0????Z?ZA??????m????ZAYAXA?

K??0?Z0?X?AAY?A? ??YAXA0ZA??19

??Y???X? (3.14) 1?? (3.15) (3.16)

河南理工大学本科毕业论文

4 坐标转换程序设计

4.1坐标转换总思路设计

如图4-1所示：本论文中采用先实现某椭球内部各形式坐标之间相互转换：包括空间直角坐标与大地坐标之间的转换、大地坐标与高斯平面直角坐标之间的转换。再通过空间直角坐标系之间旋转、平移、缩放的七参数转换实现不同参考椭球之间的坐标转换。

图 4.1 坐标转换总思路

4.2具体程序框图设计

在使用MATLAB编制坐标转换的程序时，要先进行程序框图设计，以使编程思路清晰，为方便读者理解编程的思路，下面列出各种坐标相互转换的程序框图。程序框图包括同一椭球参数下大地坐标与空间直角坐标之间的转换框图、大地坐标与高斯平面坐标的转换框图，以及不同椭球参数下两个空间直角坐标系之间的七参数转换框图。

4.2.1同一椭球参数下大地坐标与空间直角坐标之间的转换框图设计

1.大地坐标向空间直角坐标转换

20

河南理工大学本科毕业论文

图4.2大地坐标向空间直角坐标转换

2．空间直角坐标系之间的七参数转换框图设计

图4.3空间直角坐标系之间的七参数转换框图设计

3. 空间直角坐标向大地坐标转换

21

河南理工大学本科毕业论文

开始 1 输入X、Y、Z 输入椭球参数a、b 取B的迭代初值 计算椭球第一偏心率e B1?arctanZ?X2?Y2? 1 2 3 2 否 ?Z?Ne2sinB?按式B?arctan?? ??22?X?Y?B??B?0.00001? 迭代计算B值 是 3 B?B?输出B、L、H的值 ?Y?L?arctan???X?ZX2?Y22H??N1?e或?NsinBcosB??结束

图4.4空间直角坐标向大地坐标转换

22

河南理工大学本科毕业论文

4.2.2同一椭球参数下大地坐标与高斯平面坐标的转换框图设计 1. 高斯正算

图4.5 高斯正算

2. 高斯反算

图4.6 高斯反算

23

河南理工大学本科毕业论文

5 实例应用

本章从同一参考椭球下不同坐标形式之间的转换和不同参考椭球之间的坐标转换两个方面出发，结合具体工程算例，验证了上一章所述基于MATLAB编程实现坐标转换的可行性和准确度。 5.1大地坐标与空间直角坐标间换算实例

已知珠穆朗玛峰峰顶B?27?59?16.94241??，L?86?55?31.72137??，以及雪面的大地高H?8821.4016m，其坐标系为1954年北京坐标系。先用本文提供的m文件进行大地坐标向空间直角坐标换算，并与中海达坐标转换模块计算的相应值进行了比较，操作步骤如下：

1.在MATLAB系统中运行本文附录中提供的坐标转换的大地坐标转为空间直坐标的dd2kj.m，（输入时以北京54椭球的长轴a、扁率Bl以及珠穆朗玛峰峰顶大地坐标为参数）如图5.1所示；

2.得到大地坐标转为空间直角坐标的计算结果（图5.2）。

图5.1 本文程序运行

图5.2 MATLAB计算结果

3.用中海达软件中的坐标转换模块进行大地坐标向空间直角坐标的换算，界面如图5.3所示。

24

河南理工大学本科毕业论文

图5.3中海达软件大地坐标向空间直角坐标换算界面

将换算结果整理如表5.1所示：

表5.1与中海达软件从原大地坐标向空间直角坐标的换算结果 单位：m 本文换算值 本文换算与中海达软件换算差值 中海达软件换算值（X，Y，Z） （X，Y，Z）302726.854415229 5636102.39017883 2979527.61924339 302726.854415 5636102.390179 2979527.619243 0.000000229 0.00000017 0.00000039 再将计算出的空间直角坐标值用本系统和中海达软件的坐标转换模块分别换算为大地坐标并分别和原大地坐标对比，操作步骤如下；

1.在MATLAB系统中运行本文附录中提供的坐标转换的空间直角坐标转为大地坐标的kj2dd.m，（输入时以北京54椭球的长轴a、扁率Bl，以及上述换算得到的空间直角坐标为参数），得大地坐标值，如表5.2所示；

2.用中海达软件坐标转换模块进行空间直角向大地坐标的换算，界面如图5.4所示。

图5.4中海达软件大地坐标向空间直角坐标换算界面

将换算结果整理在表5.2：

25

河南理工大学本科毕业论文

表5.2 本文和中海达软件分别计算的大地坐标值与原大地坐标值对比 本文算得的大地中海达算得的大地本文与原大地中海达与原大坐标值 坐标值 坐标差值 地坐标差值 原大地坐标 （B，L，H）27?59?16.94241??27?59?16.94241?? 8821.4016m 27?59?16.942?? 86?55?31.721?? 8821.4016m 0 0 0 0.00041?? 0.00037?? 0 86?55?31.72137??86?55?31.72137?? 8821.4016m

由此例可知，本文将大地坐标换算为空间直角坐标的成果与中海达软件的坐

标转换模块计算出的空间直角坐标一致，计算成果准确；本文将空间直角坐标换算为大地坐标的成果与原大地坐标、中海达的计算成果均为一致，计算成果准确且小数点后保留位数较中海达坐标转换的模块多。 5.2 高斯正反算实例

某测区内现有点1、点2、点3三个点，已知它们在1954年北京坐标系下的

表5.3 点1、2、3的大地坐标 大地纬度B 22?15?58.98294?? 大地坐标，如表5.3所示。

点号 1 2 3 大地经度L 32?02?57.65221?? 30?23?46.65321?? 111?28?52.15387?? 118?54?15.22065?? 112?44?12.21227?? 利用本文和中海达软件的坐标转换模块对三个点的大地坐标值计算投影带为6?带的高斯平面坐标，再将二者的换算结果对比，操作步骤如下：

1. 在MATLAB下运行本文坐标转换的高斯正算的代码，得到中央子午线的经度和高斯平面直角坐标，见表5.4；

2. 用中海达的软件中的坐标转换模块进行高斯正算，界面如图5.5所示：

图5.5 中海达高斯正算操作界面

将换算结果整理如下表5.4

26

河南理工大学本科毕业论文

X1 Y1 X2 Y2 X3 Y3 表5.4 本文和中海达系的高斯正算成果及对比 本文的高斯平面坐标 中海达高斯平面坐标 2463420.565707 2463420.565707 549592.908438 549592.908438 3548973.8265572 3548973.8265572 679857.639639 679857.639639 3365384.741269 3365384.741269 666915.770959 666915.770959 差值 0 0 0 0 0 0 再将计算出的高斯平面坐标利用本文和中海达软件的坐标转换模块进行高斯反算并分别与原大地坐标对比，操作步骤如下：

1.在matlab下运行本文坐标转换的高斯反算的代码，得到中央子午线的经度和将高斯平面坐标转换为大地坐标，如表5.5所示；

2.用中海达软件的坐标转换模块进行高斯反算，界面如图5.6所示；

图5.6 中海达高斯反算界面

将换算结果整理如下表5.5： 表5.5 中海达软件与本文的高斯反算成果并分别与原大地坐标的对比 中海达软件计算的大地坐标 本文计算的大地坐原大地坐标（B，L）标 本文与原大中海达与原地坐标差值 大地坐标差值 0 0.00006?? 22?15?58.98294?? 22?15?58.98294?? 22?15?58.983?? 0 0.00003?? 111?28?52.15387?? 111?28?52.15387?? 111?28?52.154?? 0 -0.00021?? 32?02?57.65221?? 32?02?57.65221?? 32?02?57.652?? 0 0.00035?? 118?54?15.22065?? 118?54?15.22065?? 118?54?15.221?? 0 -0.00021?? 30?23?46.65321?? 30?23?46.65321?? 30?23?46.653?? 0 -0.00027?? 112?44?12.21227?? 112?44?12.21227?? 112?44?12.212?? 由本例可知，本文高斯正算成果与相应的中海达软件计算结果相比均一致计算成果准确；高斯反算成果与相应的中海达软件计算结果、原大地坐标相比均一致，计算成果准确且小数点后保留位数较中海达坐标转换模块多。

27

河南理工大学本科毕业论文

5.3 不同参考椭球基准下的坐标转换实例

由于某工程需要将WGS-84坐标系下的点位坐标转换到1954年北京坐标系

下点位坐标。已知三个公共点在两个椭球基准下的空间直角坐标分别如下：

表5.6 已知点WGS-84坐标系的坐标（单位：m） 点号 A1 A2 A3 XWGS?84 -2850017.4720 -2838514.0744 -2865534.8823 YWGS?84 4690744.5225 4704426.8235 4672522.4133 ZWGS?84 3237959.9725 3228266.3341 3250504.9271 表5.7 已知点WGS-84坐标系的坐标（单位：m） 点号 A1 A2 A3

XWGS?84 -2850023.9497 -2838520.5570 -2865541.3037 YWGS?84 4680915.4705 4694597.9568 4662693.0793 ZWGS?84 3237036.8089 3227343.6466 3249581.2094 本算例选择七参数空间转换模型来进行空间直角坐标系之间的转换。利用三个公共点通过七参数转换结果为：

表5.8 七参数结果 七参数 dX 146.9406 dY -9760.281 dZ -846.053260113258 m -3.3251144?X ?Y2.76240630512447E-2 ?Z-1.56541445702274E-2 -2.64749295656720817 4582117 7481091E-2 81441E-2 经MATLAB计算能得出七参数转换的结果，由于本人水平有限，在此不对七参数转换结果的精度进行分析，这个问题有待进一步的研究。

28

河南理工大学本科毕业论文

6 结论与展望

通过本文的研究，利用MATLAB软件进行应用性开发，把坐标转换的数学模型转化成为MATLAB中的函数子块，开发出了工程测量很关键的坐标转换程序。程序解决了工程测量技术人员在运算方面所能遇见的一些常见的坐标转换计算和数据分析处理问题。该程序能帮助工程技术人员、科研人员从繁锁的数据运算和数据分析处理中解脱出来，只要在MATLAB中调用函数，并输入函数的参数就可以实现坐标转换，避免了为解决坐标转换问题而专门花精力，大大减轻了用户的工作量，使用户有更多的精力用在其它方面。

由于时间紧再加个人能力的限制，本程序开发得还不够深入全面,未来需要进一步开展的工作，可以从下面几个方面来开展:

(1) 现在的程序需要在MATLAB环境下才能运行，下一步需要研究把坐标转换做成界面方式的，更为人性化，人机互动更为方便。

(2) 对数据的导入，本文中数据需要用户自己进行输入，但是工程应用中数据量往往是非常庞大的，下一步应开发出直接与GPS、全站仪等数据获取仪器交换数据的端口，从而方便的处理数据。同时也应开发出与Exel等数据处理软件交流的端口。

(3) 对运算结果的显示输出中，是用MATLAB软件输出,下一步要研究将其改进为GPS或是全站仪支持的数据格式进行文件输出，使用户直接进行工程测量工作。

(4) 在程序的功能方面，各坐标基准转换精度有待加强,下一步需要研究进一步解决，七参数计算的精度也有待验证。

解决了以上问题,才能真正满足目前施工一线测量人员的需求。这是我今后工作学习生活中应不断努力的，争取能够早日完善它。

29

河南理工大学本科毕业论文

致谢

在本文的写作过程中,特别要感谢胡圣武老师,他认真的阅读了拙作,对本人的论文提出了宝贵的建议,也为我解答了许多疑惑。同时也要感谢那些在论文写作过程中给予我无私帮助的同学们。

30

河南理工大学本科毕业论文

参考文献

[1] 刘亚静,毛善君,郭达志,等.基于VC++的坐标系统转换设计与实现[J].湖南科技大学学报（自然科学版），2006，9（3）：61-64. [2]宁津生.现代大地测量参考系统[J].测绘学报，2002，（5）：7-11. [3]王沫然.MATLAB与科学计算[M].北京:电子工业出版社,2004.2.

[4]姚连璧，周小平.基于MATLAB的控制网平差程序设计[M].上海:同济大学出版社,2006.4.

[5]边少锋，柴洪洲，金际航，等.大地坐标系与大地基准[M].北京:国防工业出版社,2005.5. [6]孔祥元，郭际明，刘宗泉，等.大地测量学基础[M].武汉:武汉大学出版社,2006.1. [7]李开友. 基于MATLAB的工程运算可视化系统的设计与实现:(硕士学位论文). 昆明理工大学,2006.

[8]董钧祥，杨德宏，等.测量坐标转换模型及其应用[J].昆明理工大学学报(理工版)，2006，4(3):1-4.

[9]刘大杰,施一民，等.全球定位系统（GPS）的原理与数据处理[M].上海:同济大学出版社,1996.

[10]孔祥元,梅是义，等.控制测量学(下)[M].武汉:武汉大学出版社,2002. [11]施一民.现代大地控制测量[M].上海:同济大学出版社,2003.6.

[12]陈清礼,胡家华.大地坐标转换程序[J].江汉石油学院学报，1998，3(1):45-49. [13]范新云.利用GPS测定地方坐标系转换的四参数法[J].海洋测绘，2005，7(4):35-37. [14]Bowring B. R Transformation from Spatial to Geographic Coordinates. Survey Review. No.181,1987.

[15]Paul,M. K. A Note on Computation of Geographic Coordinates. Bulletin Geodesique. No.108,1978.

[16]李岳.坐标系统的设计与实现：[D].武汉：中国地质大学，2010.

[17]徐生望，周建国，王超，等.坐标转换模型问题研究[J].西部探矿工，2009，（2）：162-165.

[18]黄海礼，官云兰，韩子太，等.基于MATLAB的坐标转换实践[J].测绘科学，2012，27（1）：20-22.

[19]郑航行.坐标转换研究：[D].郑州：解放军信息工程大学，2007.

[20]刘大杰，陶本藻，等.实用测量数据处理方法[M].北京：测绘出版社，2000. [21]姚吉利.三维坐标转换参数直接计算的严密公式[J].测绘通报，2006，（5）：7-10.

31

河南理工大学本科毕业论文

附录

1.大地坐标转空间直角的matlab计算代码 %大地坐标转换成空间直角坐标函数dd2kj.m

% -------------------------------------------------------------------- function XYZ=dd2kj(a,Bl,BB) format long g global XYZ;

B=BB(:,1);L=BB(:,2);H=BB(:,3); %将输入数据分配给变量 b=a*(1-Bl);

e=sqrt(a^2-b^2)/a; B=d2r(B); L=d2r(L);

N=a./sqrt(1-sin(B).^2.*e^2); X=(N+H).*cos(B).*cos(L); Y=(N+H).*cos(B).*sin(L); Z=(N.*(1-e^2)+H).*sin(B); XYZ=[X Y Z];

% ------------------------------------------------------------------------- % ------------------------------------------------------------------------

2.空间直角坐标转为大地坐标matlab计算代码 %空间直角坐标转换为大地坐标转换函数kj2dd.m

% -------------------------------------------------------------------- function BLH=kj2dd(a,Bl,XYZ) format long g global BLH BLH=[];

X=XYZ(1,1);Y=XYZ(1,2);Z=XYZ(1,3);%将输入数据分配给变量 b=a*(1-Bl);

e=sqrt(a^2-b^2)/a;

L=atan(Y/X);%求大地经度 if L<0

L=L+pi; end

B(1)=atan((Z/sqrt(X^2+Y^2))*(1+e^2)); for i=1:4%迭代求大地纬度

B(i+1)=atan((1/sqrt(X^2+Y^2))*(Z+(a*e^2*tan(B(i))/(sqrt(1+(1-e^2)*tan(B(i))))))); end

Bn=B(i+1);

H=sqrt(X^2+Y^2)*cos(Bn)+Z*sin(Bn)-a/sqrt(1-sin(Bn)^2*e^2)*(1-sin(Bn)^2*e^2);%求大地高

32

河南理工大学本科毕业论文

Bn=r2d(Bn); L=r2d(L);

BLH=[BLH;Bn L H]; end

% -------------------------------------------------------------------- % --------------------------------------------------------------------

3.高斯正反算matlab计算代码

正算：

%高斯正算[x,y]=BLtoxy(B,L,L0,ellipse) function [x,y]=BLtoxy(B,L,L0,ellipse) switch(ellipse)

case 0 %WGS-84椭球 a=6378137;

b=6356752.3142;

case 1 T椭球 a=6378245;

b=6356863.0187730473;

case 2 ?椭球 a=6378140;

b=6356755.2881575287; end

format long %定义数据类型为双精度型

l=(L-L0)*60*60/206265 ; %把角度转为弧度 B=pi*B/180;

%求取第一二偏心率以及卯酉圈、和各个参数的数值 e1=sqrt(a^2-b^2)/a; e2=sqrt(a^2-b^2)/b;

N=a/sqrt(1-e1^2*sin(B)^2); t=tan(B); n=e2*cos(B);

a0=1+3/4*e1^2+45/64*e1^4+350/512*e1^6+11025/16384*e1^8;

a2=-(3/4*e1^2+60/64*e1^4+60/64*e1^4+525/512*e1^6+17640/16384*e1^8)/2; a4=(15/64*e1^4+210/512*e1^6+8820/16384*e1^8)/4; a6=-(35/512*e1^6+2520/16384*e1^8)/6; a8=(315/16384*e1^8)/8;

X=a*(1-e1^2)*(a0*B+a2*sin(2*B)+a4*sin(4*B)+a6*sin(6*B)+a8*sin(8*B));

x=X+1/2*N*cos(B)^2*l^2+1/24*N*t*(5-t^2+9*n^2+4*n^4)*cos(B)^4*l^4+1/720*N*t*(61-58*t^2+t^4)*cos(B)^6*l^6

y=N*cos(B)*l+1/6*N*(1-t^2+n^2)*cos(B)^3*l^3+1/120*N*(5-18*t^2+t^4+14*n^2-58*n^2*t^2)*cos(B)^5*l^5; y=y+500000 反算：

%高斯反算[B,L]=xytoBL(x,y,L0,ellipse)

33

河南理工大学本科毕业论文

function [B,L]=xytoBL(x,y,L0,ellipse) format long

'ellipse值：选择WGS-84椭球输入0，选择54椭球输入1，选择80椭球输入2' switch(ellipse)

case 0 %WGS-84椭球 a=6378137;

b=6356752.3142;

case 1 T椭球 a=6378245;

b=6356863.0187730473;

case 2 ?椭球 a=6378140;

b=6356755.2881575287; end

e1=sqrt(a^2-b^2)/a; e2=sqrt(a^2-b^2)/b; y=y-500000;

A0=1+3/4*e1^2+45/64*e1^4+350/512*e1^6+11025/16384*e1^8+43659/65536*e1^10; B0=X/a/(1-e1^2)/A0;

k0=1/2*(3/4*e1^2+45/64*e1^4+350/512*e1^6+11025/16384*e1^8); k2=-1/3*(63/64*e1^4+1108/512*e1^6+58239/16384*e1^8); k4=1/3*(604/512*e1^6+68484/16384*e1^8); k6=-1/3*(26328/16384*e1^8); si=sin(B0)^2;

Bf=B0+sin(2*B0)*(k0+si*(k2+si*(k4+k6*si))); tf=tan(Bf);

ankef=e2*cos(Bf); c=a/sqrt(1+ankef^2); Nf=c/Vf;

Mf=a*(1-e1^2)/(1-e1^2*sin(Bf))^1.5; an=ankef^2;

B=Bf-1/2*Vf^2*tf*(y^2/Nf^2-1/12*(5+3*tf^2+an-9*an*tf^2)*y^4/Nf^4+1/360*(61+90*tf^2+45*tf^4)*y^6/Nf^6);

l=y/Nf/cos(Bf)-y^3*(1+2*tf^2+an)/6/Nf^3/cos(Bf)+y^5*(5+28*tf^2+24*tf^4+6*an+8*an*tf^2)/120/Nf^5/cos(Bf); B=B*180/pi; l=l*180/pi; L=L0+l;

4.空间直角与空间直角之间坐标转换matlab计算代码 %七参数计算函数burse.m

% -------------------------------------------------------------------- function y = burse(X) global dX

X1=X(:,1);Y1=X(:,2);Z1=X(:,3); X2=X(:,4);Y2=X(:,5);Z2=X(:,6); N=size(X1); n=N(1);

B=zeros(3*n,7);

34

河南理工大学本科毕业论文

for i=1:n

%V(3*i:3*(i+1),1)=[VX2(i) VY2(i) Vy3(i)];

B((3*i-2):3*i,1:7)=-[1 0 0 X1(i) 0 -Z1(i) Y1(i) 0 1 0 Y1(i) Z1(i) 0 -X1(i) 0 0 1 Z1(i) -Y1(i) X1(i) 0]; %dX=[dX0 dX0 dX0 a1 a2 a3 a4]';

L(3*i-2:3*i,1)=[X2(i) Y2(i) Z2(i)]'; end

P=eye(3*n);

dX=-inv(B'*P*B)*B'*P*L

% -------------------------------------------------------------------- % --------------------------------------------------------------------

35

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