巴马腕表批发销售联盟

WGS-84坐标与国家或地方坐标的转换在excel中的实现

工程测量技术2019-08-22 11:14:28

摘 要

随着GPS技术的发展,精度的提高,其以全天候,高精度及操作简单的特点被越来越广泛的运用。GPS平差后结果为大地坐标,而工程中我们常用的为国家坐标系或地方独立坐标系,所以需要进行坐标转换。本文简要介绍了WGS-84坐标系和西安80坐标系及北京54坐标系等常用坐标系。通过空间直角坐标系和大地坐标系间的关系公式,用EXCEL表实现了两者的相互转换。接下来又介绍了两个不同空间直角坐标系的关系,转换原理及模型,使用的是七参数布尔莎方法。用EXCEL表实现了不同空间直角坐标系间的互相转换,也就实现了WGS-84同其他坐标系间的互相转换,因为还涉及到换带计算,文中有添加了间接换带计算方法。

关键词:EXCEL表,WGS-84坐标系,坐标转换

WGS-84坐标与国家或地方坐标的相互转换在Excel中的实现

1973年12月美国国防部批准海陆空三军联合研制了一种新的军用卫星导航系统即GPS卫星全球定位系统。随着GPS技术的发展,精度的提高,其以全天候,高精度及操作简单的特点被越来越广泛的运用。它能方便快捷地测定出点位坐标,无论是操作上还是精度上,比全站仪等其他常规测量设备有明显的优越性。随着我国各地GPS差分台站的不断建立以及美国S A政策的取消,使得单机定位的精度大大提高,有的已经达到了亚米级精度,能够满足国土资源调查、土地利用更新、遥感监测、海域使用权清查等工作的应用。由此可以看出GPS定位技术有很大的应用潜力,近年来GPS接收机的小型化、小功耗以及卫星星座的不断改善与发展,技术也越来越成熟,实践证明,在缩短工期、降低成本和设计的灵活性方面,GPS测量较常规测量更为优越。

GPS平差后的解属WGS-84坐标系,而一些工程网要求使用原坐标系统,以充分利用原有的地形图、工程图、地下管道施工图等资料,所以进行坐标转换是必要的。如果原采用的是我国的8O坐标系或54坐标系(整体平差值),我们需要把WGS-84坐标转化为国家坐标。很多情况下为了减少投影变形等原因,采用的是地方独立坐标系,所以我们也要掌握WGS-84坐标向地方独立坐标转换的方法。还有很多时候我们为了计算等需要还可能反过来,把我们现有的国家坐标或地方坐标转换为WGS-84坐标。

坐标转换的方法有很多,其运算非常复杂,语言也很难理解,经常需要运用软件转换而且很多软件需要一个一个进行转换,但用EXCEL可以进行批量转换,极大的节省了时间。EXCEL表是大家熟知的一种比较普及的电子表格软件,其操作简单,可视化强,功能强大,不只能做简单的加减乘除运算而且能够处理复杂的数学运算,迭加,矩阵都可以在EXCEL表中实现,用它来进行坐标转换轻松自如。本文就是运用EXCEL表进行相关计算以及坐标转换的实现。

常用坐标系

进行坐标转换,首先要适当了解各个坐标系统。下面简要介绍了几种常用坐标系,以便更好的理解坐标系间的关系及转换。

坐标系

GPS卫星定位测量是用三维地心坐标系——WGS-84坐标系为依据来测定和表示总的空间位置,它既可以用地心空间坐标系(X,Y,Z)表示,也可以用托球大地坐标系(B,L,H)表示。WGS-84坐标系属于协议地球地心坐标系,最初是由美国国防部根据TRANSIT导航卫星系统的多普勒观测数据所建立的,1987年1月作为GPS卫星所发布的广播星历的坐标参照基准。

该坐标系的原点是地球的质心, Z轴指向BIH1984.0定义的协议地球极CTP方向,X轴指向BIH1984.0 零度子午面和CTP对应的赤道的交点,Y轴和Z、X轴构成右手坐标系,它的椭球如图1所示。

椭球的4个基本参数是:

长半轴 a=6378137m

地球引力常数(含大气层)

GM=3986005×108m3s-2

正常化二阶带球谐系数

2.0=-484.16685×10-6

地球自转角速度 ω=7292115×10-11rad/s

根据以上4个参数可以进一步求得:

地球扁率α =0.00335281066474

第一偏心率平方e2=0.0066943799013

第二偏心率平方e’2=0.00673949674227

赤道正常重力 γe=9.7803267714m/s2

极正常重力γp=9.8321863685m/s2

北京54坐标系

北京54坐标系是一个参心大地坐标系,原点是前苏联的普尔科沃,采用了克拉索夫斯基椭球参数,其中高程异常是以前苏联1955年大地水准面差距重新平差的结果为依据,按我国的天文水准路线换算过来的。中华人民共和国建立后,在全国范围内开展了正规的、全面的大地测量和测图工作,由此迫切建立了北京54坐标系,其与前苏联1942年坐标系进行联测,通过计算建立了我国大地坐标系,因此北京54坐标系可以认为是前苏联1942年坐标系的延伸。

西安80坐标系

西安80坐标系也是一个参心大地坐标系,原点为我国陕西省泾阳县永乐镇,采用国际大地测量和地球物理联合会1975年推荐的四个地球椭球基本参数。其椭球短轴平行于地球质心指向我国地极原点JYD1968.0方向,大地起始子午面平行于格林尼治平均天文台的子午面。西安80坐标系与北京54坐标系相比,椭球参数更加精确,椭球面也与我国大地水准面更加吻合,是在北京54坐标系的基础上建立起来的。

地方独立坐标系

在一些城市测量和工程测量中,若直接采用国家坐标系可能会由于远离中央子午线,造成变形严重,误差太大,还有些特殊的测量如桥梁大坝等工程使用国家坐标系也很不方便,因此,常常会建立适合本地区的独立坐标系。

空间直角坐标系和大地坐标系的互相转换

坐标转换中不可避免的会用到空间直角坐标和大地坐标之间的转换。比如有时我们GPS平差后的结果为大地坐标系,与其他坐标转换时需要先转换为空间直角坐标。当检核时,我们还需要反过来把空间直角坐标转换为大地坐标。所以先介绍一下两者之间的关系以及转换的实现。

由图2我们可以推出空间直角坐标与椭球大地坐标之间的公式,进而实现两者的互相转换。以下是公式:

表1 椭球参数


克拉索夫斯基椭球体1975年国际椭球体WGS-84椭球体
a6378245.00000000006378140.00000000006378137.0000000000
b6356863.01877304736356755.28815752876356752.3142
c6399698.90178271106399596.65198801056399593.6258
α1/298.31/298.2571/298.257223563
e20.0066934216229660.0066943849995880.0066943799013
e'20.0067385254146830.0067395018194730.00673949674227

每个参考椭球的形状和大小由五个基本参数来决定,长半轴a,短半轴b,椭圆的扁率α,椭圆的第一偏心率e,第二偏心率e′。但只需知道其中两个参数就(至少知道一个长度元素)可以算出其他元素。

我们以WGS-84椭球为例,进行空间直角坐标同椭球大地坐标的相互转换。

大地坐标系向空间直角坐标系的转换

表2 大地坐标向空间直角坐标转换的准备数据

表2中B,L,H是以度分秒的格式输入的,如38度48分7.38137秒写成38.480738137

B’,L’,H’列分别为把B,L,H转成度的形式。

表格H3中输入公式:

=INT(B3)+(INT(B3*100)-INT(B3)*100)/60+MOD(B3,INT(B3*100)/100)*1000/3600

表格I3中输入公式:

=INT(C3)+(INT(C3*100)-INT(C3)*100)/60+MOD(C3,INT(C3*100)/100)*1000/3600

表格J3中输入公式:

=INT(D3)+(INT(D3*100)-INT(D3)*100)/60+MOD(D3,INT(D3*100)/100)*1000/3600

表格M2为N的计算:

输入公式=K3/SQRT(1-L3*POWER(SIN(RADIANS(H3)),2))

根据公式1计算出X,Y,Z。

表格D3中输入公式=(M3+J3)*COS(RADIANS(H3))*COS(RADIANS(I3))计算出X

表格E3中输入公式=(M3+J3)*COS(RADIANS(H3))*SIN(RADIANS(I3))计算出Y

表格F3中输入公式=(M3*(1-L3)+J3)*SIN(RADIANS(H3))计算出Z

隐藏计算过程,使界面更直观简洁.

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

大地坐标向空间直角坐标转换完毕。

空间直角坐标系向大地坐标系的转换

步骤一:求出B的值

因此

式中 k=(3)

ti为前一次迭代值,第一次迭代令ti=t0

步骤二:求出L

步骤三:求出H

可以用EXCEL表完成相关计算及转换。一般迭代三四次就能算出B的值,在这里我们可以多迭代几次。以WGS-84坐标系椭球为例。先输入已知量X,Y,Z,c,e2,e’2,直接拖动鼠标即可。

表4 空间直角坐标转换成大地坐标已知量

根据式3,计算出t0,P,k,并拖到鼠标填充表格。

表5 空间直角坐标转换成大地坐标的已知计算量

G列输入公式

H列输入公式

I列输入公式

表6 迭加量

注:J2中输入=H2,即给t1赋值为t0。J3中输入公式=IF(ABS(J2-K2)<0.00000001,"",K2)。

IF(ABS(J2-K2)<0.00000001,"",K2)意思是如果J2和K2足够相近,那么也就求出了tanB的值,此空就没必要填了。如果不是还需要进行下一次迭代计算,此时令J3=K2,依次类推。

K2中填写公式=IF(J2="","",I2+G2*J2/SQRT(H2+J2*J2))。意思是如果J2为空,也就是已经求出tanB的值了,此时就不需要填写,否则继续迭代计算,最后算出B的值。

L和H的值计算比较简单,

运用公式求出L,H的值。

表7 所求值

注: M列公式,化为角度

注:N列公式

这种方法我们发现两大缺点,一是迭代时需要手动操作;二是我国地区所处经纬度都为正,上述显然考虑不周。

经过尝试我们现在选用一种更简单快捷的方法

输入已知坐标

并输入e2,e’2,c,P,k,t0,N,这些计算时需要的量

ti下面输入公式

=N3+I3*O3/SQRT(J3+O3*O3)

此时点击工具——选项——重新计算,勾选迭代计算。

确定后就求出了我们需要的值。接下来的计算就很简单。但是求L的值时我们要考虑当X大于零时,说明0

L下填入下面的公式:

=IF(B3>0,DEGREES(ATAN(C3/B3)),DEGREES(3.1415926+ATAN(C3/B3)))

为了界面美观简洁我们可以把不需要的列隐藏,即隐藏计算过程,最后界面如下:

表8 空间直角坐标向大地坐标转换方法2

我们可以根据已知数据进行精度评定。

表9空间直角坐标和大地坐标精度评定

J列输入公式=(E3B3)/B3

K列输入公式=(F3-C3)/C3

L列输入公式=(G3-D3)/D3

就此我们完成了大地坐标和空间直角坐标的转换。

坐标系和其他坐标系的转换

七参数布尔莎方法进行坐标转换

GPS的测量结果和我国的西安80或北京54系坐标相差很多,有时甚至100多米,由此可见我们必须进行坐标转换才行,以空间直角坐标转换为例。进行空间直角两坐标的转换,除对坐标原点实施三个平移参数外,当坐标轴间互不平行时还存在三个旋转角度参数,以及两个坐标系尺度不一样的一个尺度变换参数。我们选用布尔莎七参数公式进行推导,如果忽略旋转参数或尺度参数等,可以推出五参数,四参数,三参数转换公式。

图3为布尔莎模型

由布尔莎模型我们可以推出任意点在两坐标间的转换公式,如公式3所示

其中为平移参数,dK为尺度变换参数 ,为旋转后的坐标

根据公式3我们可以完成任意两个空间直角坐标系统相互转换。

当已知多个公共点(两个以上)时按最小二乘法求解转换参数。

(4)

式中i=1,2,3……,N,若设

,

Y=,

(5)

则公式3变为误差方程

设观测值等全观测,则PLΔx=E,则法方程

进而求出转换参数

单位权方差

协因数阵

若不是等权观测,则

单位权方差

协因数阵

为了方便计算及理解,在这里设为等权观测,然后把这些公式输入到相应的EXCEL表中,进行任意坐标转换的实现。

已知三对WGS-84和北京54的对应空间直角坐标坐标((-2142730.263,4492272.368,3975216.033)(-2142730.16,4492398.27, 3975268.731));((-2124722.344,4525873.151,3946900.096)(-2124722.276242,4525999.450226, 4525873.151));((-2161391.350,4508812.888)(-2161390.903938, 4508938.774386, 3946534.491616)

对(-2138714.416,4541513.853,3921454.813)进行坐标转换。

具体过程如下:

输入已知坐标,如表10所示

表10 已知坐标

在表格中输入公式求得LΔX

表格中的公式不用一一输入,只需拖拽即可。

然后根据前文中的公式在表格中输入B如表11所示。

表11 B矩阵

注:B中需要输入Xi等坐标的地方,需要输入公式,如K4中输入公式“=B2”。

表12 所求参数

注:根据公式6,选中P2:P8,输入公式

=MMULT(MMULT(MINVERSE(MMULT(TRANSPOSE(H2:N10),H2:N10)),TRANSPOSE(H2:N10)),G2:G10),按住CTRL+SHIFT+ENTER后,计算出各个参数。

在Q2:S4中输入,Xi,Yi,Zi为需要转换的坐标,输入时也要按公式输

入。选中U2:U4,输入公式=P2:P4+MMULT(T2:T4,P5)+MMULT(Q2:S4,P6:P8)+T2:T4,按shift+ctrl+enter,计算出最终坐标。

接下来根据公式7计算精度,进行精度评定。

表13 精度评定

精度评定步骤:V列按照公式5输入,选中V2:V8,输入=MMULT(H2:N10,P2:P8)-G2:G10,按CTRL+SHIFT+ENTER;W列输入已知坐标对数,X列按照公式7输入,选中X2,输入公式=MMULT(TRANSPOSE(V2:V8),V2:V8)/(3*W2-7)

转换完毕。

隐藏计算过程最终界面如下:

表14 不同空间直角坐标的转换

以上是WGS-84与北京54两个直角空间坐标系的转换。如果需要其他转换只需把已知坐标改成WGS-84坐标和西安80等你需要转换的坐标,就可以完成WGS-84和西安80和北京54等空间直角坐标的转换。

进行WGS-84与北京54、西安80坐标系间的转换一般只涉及X,Y坐标,此时我们可以设Z=0, εz=0,此时只有五个参数需要求,可简化计算。

换带转换计算

以上方法还有一个漏洞,就是已知坐标必须是在同一带中,所以当WGS-84坐标向国家,地方坐标转换时,如果已知坐标不在同一带,还需要考虑换带问题,把坐标都转到同一带下,一般选用有已知坐标比较多的带中。换带有两种方法,直接法和间接法,本文详细介绍一下应用高斯投影正反算公式间接进行换带计算的间接法。

应用高斯投影正反算公式间接进行换带计算的间接法

这种方法的实质是把椭球大地坐标作为过渡坐标。过程如下

(x,y)1→→→→→→(B,l) →→→→→→(x,y)2

具体通过一个实例说明一下。比如在中央子午线=123。的1带中,有一点的平面直角坐标x1=5728374.726m,y1=+210198.193m,现要求计算出在=129。的第2带的平面直角坐标。

先由x1,y1换算B1,L1。B1=51°38 ′43.9024″,L1=126°02′13.1362″,L=123°+126°02′13.1362″=249°02′13.1362″。 L2=249°02′13.1362″-129°=120°02′13.1362″,然后由B,L2换算x2=5728164.378m,y2=-205079.963m。

就此我们完成了此次坐标的相互转换。

结论及建议

根据上述论述,我们可以运用EXCEL表方便快捷的进行大地空间坐标与空间直角坐标的转换以及任意空间直角坐标系的转换。EXCEL表可视化强,操作简单,是进行坐标转换的很好选择。

进行不同空间直角坐标系的转换时运用的布尔莎模型即公式3有一个缺点:一般不同空间直角坐标系旋转角度很小,即设sinε=ε,cosε=1,如果角度很大误差也会相应变大。我们平常运用的北京54坐标或西安80坐标一般是空间直角坐标进行高斯投影后的坐标,高程也不是单纯的Z值。所以进行完WGS-84空间直角坐标与国家空间直角坐标转换后,要进行高斯投影。方法为国家空间直角坐标大地坐标按照国家椭球参数转换成大地坐标再进行高斯投影正算,求出平面坐标。这中间还涉及误差改正等问题,这里不再详细论述。




工程测量技术


搜索微信号: gcceliang 长按复制)


国内先进的测量技术交流圣地,钻研于最精确、最简单、最实用、最高效的测量技术。主要涉及隧道,桥梁,路基,水电站,矿山等工程。整合各类曲线程序,测量方法,测量器材。


不只是测量,每天推荐一篇值得品味的好文章





如何和查看历史文章:1、点击文章右上角的三个小点查看公众号,2、扫描上方二维码,3、点击标题下方蓝色的“工程测量技术”,4、搜索微信号:gcceliang 。都可以查看历史文章。




黑龙江市政公路测量及造价


搜索微信号:HLJSZGL (←长按复制)


时时发布路桥行业动态信息、求职招聘、技术交流、视频软件、资源共享、测量测绘、实战培训(测量、内页资料、造价)等。






建造师之家


微信号: jzshome (←长按复制)


每天推送考试动态、行业资讯、考试技巧、考生经验。定期更新一建、二建课件。定期发布各科实务专业口诀。每周推送一次证书挂靠资讯。 一建QQ群号:245867243 二建QQ群号:427983591


Copyright © 巴马腕表批发销售联盟@2017