赵 丹
(中铁第四勘察设计院集团有限公司,湖北武汉,430063)
Geodetic Coordinates Transformation Programming Based on VC++
Zhao Dan
摘要:介绍了工程应用中常见的大地坐标系统,并介绍了三参数、四参数、七参数坐标转换模型,以及大地坐标系与空间直角坐标系之间的转换模型,同时利用VC++开发工具设计了实用的坐标转换程序,实现了相应的坐标转换功能。
关键字: 坐标转换;WGS-84;BJ54;GZD80;VC++;
1 引言
常用的坐标转换一般包括各种空间直角坐标系与大地坐标系,地心空间直角坐标系与参心空间直角坐标系,以及不同参心空间直角坐标系之间的相互转换。当不同坐标系之间存在严密的数学转换模型时,可以采用相应的模型之间进行坐标转换。当不同坐标系之间不存在明确的函数关系时,常采用三参数、四参数、七参数坐标转换模型。本文介绍常用的坐标转换模型,同时利用VC++开发工具设计了实用的坐标转换程序,实现了相应的坐标转换功能。
2 坐标转换模型
2.1 七参数法
工程应用中常常需要将不同参考椭球下的空间直角坐标转换。对于既有平移、旋转,又有缩放的两个空间直角坐标系的坐标转换,相应地存在3个平移参数、3个旋转参数以及1个尺度参数,共计7个参数,称为七参数模型。相应的坐标转换模型为:
式中,为3个平移参数;为3个旋转参数,为尺度参数。为了求得这7个转换参数,至少需要3个公共点,当多余3个公共点时,按最小二乘法求得7个参数的最或然值。
2.2 三参数法
三参数法是七参数法的特殊情况,当测区所覆盖的范围不大,精度要求不高时,可以只考虑坐标轴的平移,而不考虑旋转和缩放,相应地只存在3个平移参数,即三参数模型。相应的坐标转换模型为:
2.3 四参数法
由于高斯投影带的划分不同,会造成常用的大地坐标系(如BJ54,GZD80,WGS-84等)在某些地区的投影不能满足工程要求的精度时,此时就需要建立地方坐标系。此时面临的问题是如何将常用的大地坐标系转换成地方坐标系。这里就要用到四参数模型进行转换,相应的坐标转换模型为:
2.2 空间直角坐标系与大地坐标系的相互转换[1]
同一地面点在空间直角坐标系中的坐标和在大地坐标系中的坐标可用如下两组公式进行转换。
从大地坐标转换到空间直角坐标的转换模型如下式:
利用点的空间直角坐标到大地坐标的转换模型如下式:
式中:e为参考椭球第一偏心率,可由长短半径按式算得。N为法线长度,可由式算得。大地纬度的计算比较复杂,通常采用迭代法。迭代时可取,用的初值计算和,按进行第二次迭代,直至最后两次值之差小于限差为止。
3 常用大地坐标系
3.1参心坐标系
参心大地坐标的应用十分广泛,它是经典大地测量的一种通用坐标系。根据地图投影理论,参心大地坐标系可以通过高斯投影计算转化为平面直角坐标系,为地形测量和工程测量提供控制基础。由于不同时期采用的地球椭球不同或其定位与定向不同,在新中国历史上先后建立了北京54坐标系(旧)(简记BJ54旧)、1980年国家大地坐标系(也称西安80坐标系,简记GZD80)和1954年北京坐标系(新)(简记BJ54新)等三种参心大地坐标系[2]。
1.1954年北京坐标系(BJ54(旧))
1954年总参谋部测绘局在有关方面的建议和支持下,采取中国和前苏联远东一等锁相连接、平差,这样传过来的坐标定名为1954年北京坐标系(旧)。几十年来,中国在该坐标系上完成了大量的测绘工作,实施了天文大地网局部平差,通过高斯-克吕格投影,得到点的平面坐标,测绘了各种比例尺的地形图。但是随着测绘新理论、新技术的不断发展,人们发现该坐标系存在诸多缺点。另外,该坐标系是按分区平差逐步提供大地点成果的,在分区的结合部产生了较大的不符值。但该坐标系确实在测绘生产中发挥了巨大的作用,至今仍在一些部门使用。1954 年北京坐标系(旧)的缺陷是: ①椭球参数有较大误差; ②参考椭球面与我国大地水准面存在着自西向东明显的系统性的倾斜; ③几何和物理大地测量应用的参考面不统一; ④定向不明确; ⑤当时使用的仪器、测量方法的落后,致使在大面积长距离传递中误差累计较大,且系统只是进行了局部的平差。因而不可避免地出现一些矛盾和不够合理的地方。后者优点是: ①椭球短轴平行于地球地轴; ②起始大地子午面平行于格林尼治天文台起始子午面; ③椭球面同似大地水准面在我国境内最为密合; ④系统经过了整体平差。所以向1980年国家大地坐标系过渡及1980年国家大地坐标系的全面使用就成了必然。
2.1980年国家大地坐标系(GDZ80)
为了解决1954年北京坐标系中存在的椭球参数不够精确、参考椭球与中国大地水准面密合不好等缺点,1978年4月于西安召开了全国天文大地网整体平差会议,决定建立新的国家大地坐标系——1980年国家大地坐标系(GDZ80),也称西安80坐标系。大地坐标系的原点设在陕西省泾阳县永乐镇。采用了国际大地测量与地球物理联合会(IUGG)1975年推荐的四个地球椭球基本参数,并根据这四个参数求解椭球扁率和其他参数。1980年国家大地坐标系的椭球短轴平行于地球质心指向我国地级原点JYD1968.0方向,大地起始子午面平行于格林尼治平均天文台的子午面。椭球定位参数以我国范围内高程异常值平方和最小为准则求解。
3.1954年北京坐标系(BJ54(新))
由于几十年来,我国数十万个国家控制点都是在1954年原北京坐标系内完成计算的,一切测量工程和测绘成果均无一例外地采用着这个系统,考虑到1980年国家大地坐标系有着它的先进性和严密性,于是就产生了1954年新北京坐标系。1954年北京坐标系(新)是由1980年国家大地坐标系转换得来的,简称BJ54(新)。BJ54(新)是在GDZ80的基础上,改变GDZ80相对应的IUGG1975椭球几何参数为克拉索夫斯基椭球参数,并将坐标原点(椭球中心)平移,使坐标轴保持平行而建立起来的。GZD80坐标系与BJ54(新)其实是一种椭球参数的转换,这种转换在不同的椭球之间的转换是不严密的,因此不存在一套转换参数可以全国通用的,在每个地方会不一样,因为它们是两个不同的椭球基准。在两个椭球间的坐标转换,一般而言比较严密的是用七参数法。
3.2 WGS-84坐标系
GPS所采用的坐标系为WGS-84地心坐标系,它是以地球的质心为原点的地心坐标系,X, Y 轴在地球赤道平面内, Z 轴与地球自转轴相重合。WGS-84大地坐标系是一个协议地球坐标参考系CTS(Conventional Terrestrial System) ,其几何定义是:原点位于地球质心, Z 轴指向B IHl984. 0定义的协议地极CTP(Conventional Terrestrial Pole)方向, X轴指向BIHl984.0的零子午面和CTP赤道的交点, Y轴与X, Z轴正交构成右手坐标系。对应于WGS-84大地坐标系有WGS-84椭球。WGS-84椭球采用的基本参数与克拉索夫斯基椭球的基本参数有所不同,表1列出了各椭球参数[4]
表1. 常用大地坐标系采用的椭球参数
克拉索夫斯基椭球 | IUGG1975 | WGS-84椭球 | |
f | 6378 245.000 000 000 0 (m) 6356 863.018 773 047 3 (m) 1/298.3 0.006 693 421 622 966 | 6378 140.000 000 000 0 (m) 6356 755.288 157 528 7(m) 1/298.257 221 010 0.06 694 384 999 588 | 6378 137.000 0 (m) 6356 752.314 2 (m) 1/298.257 223 563 0.006 694 379 901 3 |
4 坐标转换程序实现
首先根据坐标转换程序设计的需要,定义了点的坐标结构体以及参考椭球参数结构体:
//点的坐标结构体:
struct Point2D //点的平面坐标结构体
{
double X; //X坐标
double Y; //Y坐标
Point2D (){X = 0.0;Y = 0.0;} //类似于类的构造函数
};
struct Point3D //点的三维坐标结构体
{
double X; //X坐标
double Y; //Y坐标
double Z; //Z坐标
Point3D (){X = 0.0;Y = 0.0;Z = 0.0} //类似于类的构造函数
};
struct EllipsePara //参考椭球参数结构体
{
double a; //椭球的长半轴
double b; //椭球的短半轴
double e_square; //椭球第一偏心率的平方
double f; //椭球的扁率
EllipsePara(){a=0.0; b=0.0; e=0.0; f=0.0;}
};
然后定义了相应的坐标转换函数:
//选择参考椭球的类型,获取其参数(EllipseID参考椭球的类型)
BOOL GetEllipsePara(int EllipseID)
{
if(EllipseID==0) //IUGG1975
{
pEllipsePara->a= 6378140.0;
pEllipsePara->e_square= 0.006694384999588;
pEllipsePara->f= 1.0/298.25722101;
pEllipsePara->b= pEllipsePara->a*(1.0-pEllipsePara->f);
return true;
}
if(EllipseID==1) //WGS-84
{
pEllipsePara->a= 6378137.0;
pEllipsePara->e_square= 0.006694379901413;
pEllipsePara->f= 1.0/298.257223563;
pEllipsePara->b= pEllipsePara->a*(1.0-pEllipsePara->f);
return true;
}
if(EllipseID==2) //Krasovsky
{
pEllipsePara->a = 6378245.0;
pEllipsePara->e_square = 0.006693421622966;
pEllipsePara->f = 1.0/298.3;
pEllipsePara->b= 6356863.0187730473;
return true;
}
return false;
}
//将大地坐标转换为空间直角坐标(参考椭球的类型EllipseID, 点的大地坐标BLH, 点的空间直角坐标XYZ)
void BLHtoXYZ(int EllipseID, Point3D BLH, Point3D XYZ);
//将空间直角坐标转换为大地坐标(参考椭球的类型EllipseID, 点的空间直角坐标XYZ, 点的大地坐标BLH)
void XYZtoBLH(int EllipseID, Point3D XYZ, Point3D BLH);
//七参数法坐标转换函数(三个公共点在A坐标系中的坐标point1A,point2A,point3A, 三个公共点在B坐标系中的坐标point1B,point2B,point3B, 非公共点(即是待转换的点)在A,B坐标系中的坐标point4A,point4B)
void SevenParametersCoordinateTrans(Point3D point1A, Point3D point2A, Point3D point3A, Point3D point4A, Point3D point1B, Point3D point2B, Point3D point3B, Point3D point4B );
//三参数法坐标转换函数(公共点在A坐标系中的坐标point1A,公共点在B坐标系中的坐标point1B,非公共点(即是待转换的点)在A,B坐标系中的坐标point2A,point2B)
void ThreeParametersCoordinateTrans(Point3D point1A, Point3D point1B, Point3D point2A, Point3D point2B);
//四参数法坐标转换函数(两个公共点在A坐标系中的坐标point1A, point2A,公共点在B坐标系中的坐标point1B,point2B,非公共点(即是待转换的点)在A,B坐标系中的坐标point3A,point3B)
void FourParametersCoordinateTrans(Point2D point1A, Point2D point2A ,Point2D point1B, Point2D poin2B ,Point2D point3A, Point2D point3B)。
5 结语
本文介绍了目前工程应用中常见的大地坐标系统,并介绍了三参数、四参数、七参数坐标转换模型,以及大地坐标系与空间直角坐标系之间的转换模型,同时利用VC++开发工具设计了实用的坐标转换程序,实现了相应的坐标转换功能。利用VC++ 语言定义了相应的结构体和坐标转换函数, 实现了相应的坐标转换功能。
参考文献
[1] 孔祥元,郭际明,刘中泉.大地测量学基础[M]. 武汉:武汉大学出版社,2006.
[2] 朱华统.GPS坐标系统的变换[M]. 北京:测绘出版社,1994.
[3] 朱华统.常用大地坐标系及其变换[M]. 北京: 出版社,1990.
[4] 徐绍铨. GPS测量原理及应用[M]. 武汉:武汉测绘科技大学出版社,1998.