WGS84坐标转换到BJ54坐标的方法的实验研究
第一章绪论
随着GPS 定位精度的不断提高,GPS 技术在测量中的应用也越来越广泛。但是由于GPS 卫星星历表示于WGS84 坐标系中,算得的GPS 定位结果只能表示在WGS84 全球坐标系中。WGS84坐标系是一种质心坐标系统，其坐标原点位于地球的质心上。而我国的国土测量成果和在进行工程施工时大都采用BJ54坐标系,它是一种参心坐标系，它以克拉索夫斯基椭球为参考椭球, 并采用高斯—克吕格投影（等角横切圆锥投影）方式进行投影, 如何实现WGS84地心空间直角坐标系与BJ54在平面直角坐标系的之间的转换，一直是各个部门关心的热点。
在进行WGS84坐标系和BJ54坐标系转换时有两种转换思想和模型，即平面转换模型和空间转换模型。在平面转换模型中，首先要假定两种坐标系的中心和坐标轴的方向一致，所以只适合小范围内国土测量和工程设计使用，平面转换模型原理简单,数值稳定可靠。要进行大范围的GPS测量，应该使用空间转换模型。按实际情况又分为7 参数转换和3 参数转换两种。鉴于54 坐标点的大地高通常不能精确得知,对这两种转换方法得到的平面坐标的精度进行了比较,得出大地高精度主要表现为对高程的影响,对平面坐标影响较小的结论。此外,还讨论了7 参数与3 参数模型对转换结果的影响。
第二章坐标系统简介
第一节坐标系统
2.1坐标系统
所谓的坐标系指的是描述空间位置的表达形式，即采用什么方法表示空间位置。人们为了描述空间位置，采用了多种方法，从而产生了不同的坐标系。在各种测量中，经常使用的坐标系有3种类型。
2.1.1 空间直角坐标系
空间直角坐标系的原点位于参考椭球的中心，Z轴指向参考椭球的北极，X轴指向起始子午面与赤道的交点，Y轴位于赤道的平面上，且按右手系于X轴成90°夹角。（见图2—1）
图2—1空间直角坐标系
2.1.2 空间大地坐标系
空间大地坐标系是采用大地经度和大地高来描述空间位置的。纬度（B）是空间的点与参考椭球面的法线与赤道面的夹角，经度（L）是空间中的点与参考椭球的自转轴所在的面与参考椭球的起始子午面的夹角，大地高（H）是空间点沿参考椭球的法线的方向到参考椭球面的距离。（见图2—2）。
图2—2 空间大地坐标系
2.1.3 平面直角坐标系
平面直角坐标系是利用投影变换，将空间坐标（主要指空间大地坐标）通过某种数学变换映射到平面上，这种变换有称为投影变换。在我国经常采用的投影方式是高斯—克吕格投影。
第二节 WGS84坐标系与BJ54坐标系
2.2.1 WGS84坐标[1]
WGS84(World Geodetic System,1984年)是目前国际通用的地心坐标系，其坐标的几何意义是：原点在地球的质心，Z轴指向BIH 1984.0定义的协议地球极（CTP）方向，X轴指向BIH 1984.0的零子午面CTP赤道的交点。Y轴与Z,X轴构成右手系坐标。WGS84椭球采用IUGG第17届大会大地测量常数的推荐值。采用的四个基本参数是：
长半轴 a=6378137m
地球引力常数（含大气层） GM=3986005×108m3s-2
正常化二阶带球谐系数 C2.0=-484.16685×10-6
地球自转角速度 w=7292115×10-11rad/s
根据以上四个参数可以进一步求得：
地球扁率a=0.00335281066474
第一偏心率平方 e2=0.0066943799013
第二偏心率平方e′2=0.00673949674227
赤道正常重力 r e=9.7803267714m/s2
极正常重力 r p=9.8321863685 m/s2
自1987年1月10之后，GPS卫星星历均采用WGS84坐标系统。因此GPS网的测站坐标及测站之间的坐标差均属于WGS84系统。为了求得GPS测站点在地面坐标系（属于参心坐标）中的坐标，就必须进行坐标系的转换。
2.2.2 BJ54坐标[1]
1954年北京坐标系可以认为是前苏联1942年坐标系的延伸。他的坐标原点不在北京，而在前苏联的普尔科沃。1954年北京坐标系建立以来，我国依据这个坐标系建成了全国天文大地网，完成了大量的测绘任务。但是随着测绘新理论，新技术的不断发展，人们发现该坐标系存在缺点。但我国的国土测量成果和在进行工程施工时大都仍采用BJ54坐标系。北京54坐标有以下特点：
(1)属于参心大地坐标系；
(2)采用克拉索夫斯基椭球的两个几何参数；
(3)大地原点在原苏联的普尔科沃；
(4)高程基准为1956年青岛验湖站求出的黄海平均海水面；
(5)高程异常以原苏联1955年大地水准面重新平差结果为起算数据，按我国天文水准路线推算而得。
第三章 WGS84转换到BJ54的方法
第一节 坐标转换的基本思想
在进行空间转换时，首先必须假定WGS84坐标测定点中有3个已知的BJ54坐标系的平面坐标。根据这3个公共点的坐标首先进行7个参数的测定，然后进行WGS84坐标和BJ54坐标的转换。
在WGS84坐标与BJ54坐标的转换过程中，主要先求出坐标转换参数。无论使用哪一种方法，只有知道了转换参数，才能进行坐标转换。各个国家和地区都根据本地区的精度要求解算出本地区的转换参数，所求出的转换参数因坐标系不同而异。为保证变换精度，在一个国家或地区，只使用一个变换参数，必然带来极大的变换误差。本文仅针对某一地区进行实地测量，使用两种方法确定变换参数，再进行坐标变换，并比较转换精度。
可以通过地面点在两个坐标系的坐标求解转换参数。本文中WGS84坐标，BJ54坐标系平差联网测得，都为大地坐标。
坐标变换可以分为7参数换算模式和3参数换算模式。由于WGS84坐标系为地心坐标系,BJ54坐标系为参心坐标系，它们所对应的空间直角坐标系是不同的（见图3—1）。
图3—1 ＷＧＳ84与ＢＪ５４的关系
因此坐标转换参数不但包括平移参数，还包括旋转参数和尺度变化参数。设
[]T
D i D i
D i
D i Z Y X =X 和[]T
G i G i
G i
G i Z Y X =X 分别为BJ54坐标和WGS84坐标
的坐标向量。
WGS84坐标与BJ54坐标的坐标变换，可用下列步骤实现： （1）将两个坐标系的坐标都变换为直角坐标；
（2）按所采用的变换方法(7参数或3参数）求解转换参数； （3）根据所求参数进行坐标变换； （4）将直角坐标变为大地坐标。 直角坐标和大地坐标的关系如图4所示。
图3—2 参考椭球
其中直角坐标与大地坐标的相互转换，可以采用以下公式进行
???
??
?
??
???--=+-++==+-=+=+=)1(sin /)]})1(()/[()(arctan{)/arctan(sin ))1((sin cos )(cos cos )(2
22/1222e N B Z H H e N Y X H N Z B X Y L B H e N Z L
B H N Y L B H N X （3—1）
式中，L 为大地纬度，B 为大地经度，H 为大地高程，
2/122)sin 1/(B e a N -= （3—2） N 为该点的卯酉圈曲率半径；
2222/)(a b a e -= （3—3） a,e 分别为该大地坐标系对应椭球的长半轴和第一偏心率。
第二节 坐标转换的数学模型
GPS 测量得到的是WGS284 中的地心空间直角坐标,而工程施工中通常使用地方独立坐标系,要求得到地方平面坐标。如何实现两者的转换,下面介绍从GPS 定位结果至平面坐标转换模型。坐标转换模型主要有两种，平面转换模型和空间转换模型。
3.2.1平面转换模型
假设北京54 椭球的中心和坐标轴方向与WGS284 椭球相一致,可通过平面转换模型,将GPS 定位得到的大地经纬度和大地高( B 84 L 84 h 84)T ，通过以下过程转换成平面坐标( X g Yg)T ;
(1) 由WGS84的椭球参数, 即椭球长半径和扁率,将( B 84 L 84 h 84)T 换算至
空间直角坐标( X Y Z)T 的公式[4]
为
??
?
??
+-=+=+=B h e N Z L B h N Y L B h N X sin ])1([sin cos )(cos cos )(2 (3—4) (2) 由BJ54椭球参数,将( X Y Z)T 换算至大地坐标形式( B 54 L 54 h 54)T 的公式为:
??
?
??
-+=++==N B Y X h Y X B Ne Z B X Y L sec )(])/()sin arctan[()
/arctan(2/1222/1222 （3—5）
(3) 根据工程需要,确定中央子午线、投影面高程及北向、东向平移量,将( B 54 L 54 h 54)T 投影为Gauss 坐标( X ′g Y ’g)T 高斯投影公式为：
??
??
?
?
??????+-++-?++-?+?='???+-++-??+++-?+?+='522242322362224
25442232)5814185(cos 120/)1(cos 6/cos )330720586(cos sin 720/)495(cos sin 24/cos sin 2/l t t t B N l t B N l B N Y l t t t l B B N l t B B N Bl B N X X g g
ηηηηηηηη （3—6） 以上步骤是在假定54 椭球与WGS84 椭球的中心与坐标轴相同的前提下进行的, 但实际应用中还应考虑旋转平移缩放的问题。若GPS 测定的大量点中,已知部分点的平面坐标为( X g Yg)T , 则可写出这些点的平面坐标( X ′g Y ′g)T 与已知坐标之间的关系: ????
?
?
?ψ++????
?
?=???
?? ?
?
'
'
y x y
x
y
x g g
g g
R r )()1(00
（3—7） 其中: ( X 0 Y 0)T 为坐标平移量；r 为缩放尺度；
???
?
ψψ ??ψ-ψ=ψ)cos()sin()sin()cos()(R （3—8）
为旋转矩阵, Ψ 为旋转角。
为求出(3—8) 式中的平移、缩放尺度和旋转参数,至少需要已知两个平面点,如多于两个点,可按最小二乘求解。
对所有的GPS 测定点经过以上3 个步骤及公式(3—8) 的计算,即可求得当地平面坐标。水准高程可以由大地高h54扣除大地水准面差距求得。大地水准面差距可以根据大地水准面模型(如EGM96) 或水准重合点拟合求得。
平面转换模型原理简单,数值稳定可靠,可用于RTK 手簿软件。但由于(3--8) 式是一个线性变换公式,而Gauss 投影变形是非线性的, 它的一次项与y2g 成正比,因此平面转换模型只适合范围较小的工程使用。对于大范围的GPS 测量,应使用空间转换模型。
3.2.2空间转换模型
若GPS 测定的点中部分点的平面坐标已知,对这些已知的平面坐标T g g
Y X )(
进行Gauss 投影反算,可得到大地坐标)(5454L B ,再加上大地高h 54 ,由54 椭球
参数按( 1) 式转换成空间坐标, 以T Z Y X )(5454
54表示。GPS 直接测定点的空
间坐标以T Z Y X )(8484
84
表示,则两者的转换关系为:
?????
?
??++?????? ??=?????? ??Z Y X R R R D D D Z Y X k z y x 8484843
21545454)()()()1(γβα （3—9）
式中: T z y
x D D D )( 是空间转换坐标平移量; k 为缩放尺度参数;α、β、γ为
旋转参数。
当已知的平面点多于3 个时,由(5) 式可以反求出这7 个转换参数。由GPS 测定点的大地坐标
()T Z L B 84
84
84
,可以通过以下步骤转换成平面坐标T Y X
g g
)(:
( 1 ) 将GPS 测定的T h L B ][848484由WGS84 椭球参数按（3-4）式转换成空间坐标形式T Z Y X )(848484;
( 2 ) 按(3-9) 式将T Z Y X ][848484换算为T Z Y X ][545454 ;
(3) 根据54 椭球的椭球参数, 按(3-5) 式将T Z Y X ][545454 换算为大地坐标
T h L B ][545454 ;
(4) 由Gauss 投影正算公式(3) 求得平面坐标T
Y X g g )(。
第三节 坐标转换的简化模型
空间转换模型适用于大范围GPS 测量,但在实际施工过程中,根据施工精度的要求,又分为3 种情况:在空间转换模型(3-9) 式中,需求解7 个参数,故称为7 参数转换模型;若其中的缩放比例不变,不需求尺度参数,则称为6 参数转换;若尺度参数和旋转参数均不求,则称为3 参数转换。
对于7 参数模型的求解,至少需要3 个公共点;6 参数模型的求解也至少需要3 个公共点,因为尽管两个公共点有6 个坐标分量,按(3-9) 式可以列出6个观测方程,但这6 个坐标分量中,由于两点间的距离是固定的,只有5 个是独立的;3 参数模型可在只有1 个公共点的情况下求解。对于54 椭球,点的大地高往往不能精确得到,这对实际空间转换带来的影响就值得关注,尤其是对平面坐标的精度影响。下面分别介绍7参数和3参数模型。
3.3.1七参数转换模型[13]
?????
?
??????????????---+??????????++?????????????=??????????Z Y X Gi Gi Gi x
y
x z y z
Gi Gi Gi Di Di Di Z Y X Z Y X k Z Y X 000)1(εεεεεε (3—10) 式中，[△X,△Y,△Z]为平移转换参数，εx, εy,εz 为旋转参数，k 为 尺度因子，上式中共7个参数，称为7参数变换法，是一种比较精确的转换模型。 令),,,,,,(z y x k Z Y X R εεε???=
????
?
?????---=0
1
0000
1000
011Gi
Gi
Gi
Gi Gi Gi
Gi Gi Gi X Y Z X Z Y Y Z X C （3—11） 则（3—8）式可简写为：
R C X X i G i D i += （3—12）
现测定3个点，i=1，2，3
则可以根据最小二乘法解R ，方法如下： 令
????
??????=321C C C A ????
?
?????---=332211
C D C D C D X X X X X X b （3—13） 建立方程组
AR=b （3--14）
可求得R 时，纯量（)()(b AR b AR T --取最小值，把变量的纯量函数 取极值的定义化为方程：
0=-b A AR A T T （3—15）
解方程
b A A A R T T 1)(-= （3—16）
最后解得R ，求得7个转换参数。
3.3.2三参数转换模型[14]
采用3参数变换方法只考虑3个平移参数，方法如下：
???
???
?????Z ?Y ?X +??????????=??????????Gi Gi Gi Di Di Di Z Y X Z Y X （3—17）
当重合点较多时，利用一点求出转换参数后可以再用其他点进行修正。
第四章坐标转换的实例计算
第一节坐标转换的研究的区域
上一章介绍了7参数和3参数的换算模型，下面讨论在具体实例中的应用。
首先看了解一下所研究的模型区域，下面是所研究区域的有关数据和控制网图表4—1和表4—2是实测的WGS84和BJ54的坐标成果：
表4－1 某区域GPS控制网WGS84坐标成果表
表4-2 某区域GPS控制网BJ54-6度带坐标成果表
图4-1 控制点分布图
把以上两表中的已知坐标都转换成空间直角坐标系下的坐标，首先取点HN15，把BJ54的坐标（X54,Y54,h54）=（4372952.906，274820.104，1381.393m ）转换成（B54,L54,h54）=(39.460304737,111.3833745763,1381.393)在应用公式（3—4）转换成(X54,Y54,Z54)=（6379626.377，352.079824，290.642543）同样把点HN15的WGS84坐标（84B ,84L ，84H ）=（?273835643.39，?230303119.108，1344.432m ）应用公式（3—4）转换为（84,84,84Z Y X ）=（6379481.421，254.551743，284.907019）；
同样把其他的点分别转换成(X54,Y54,Z54)，（84,84,84Z Y X ）成果如下：
表4-3 ＢＪ５４平面坐标转换为空间坐标成果
表4-4 ＷＧＳ８４大地坐标转换为空间坐标成果
第二节 七参数转换
由上一章公式（3—10）可知WGS84坐标转换到BJ54坐标七参数的数学转换模型如下：
?????
?
??????????????---+??????????++?????????????=??????????Z Y X 848484848484545454000)1(Z Y X Z Y X k Z Y X x
y
x z y z
εεεεεε (4—1) 式中，[△X,△Y,△Z]为平移转换参数，εx, εy,εz 为旋转参数，k 为 尺度因子，上式中共7个参数，称为7参数变换法，是一种比较精确的转换模型。 令),,,,,,(z y x k Z Y X R εεε???=
????
??????---=01
0000
1000
0184
84
84
848484
848484X Y Z X Z Y Y Z X C （4—2） 则（4—1）式可简写为：
R C X X i G i D i += （4—3）
现测定3个点，i=1，2，3(分别代表点HN15,HN16,HN,18) 则可以根据最小二乘法解R ，方法如下： 令
??????????=321C C C A ????
?
?????---=Gi Di Gi Di Gi Di
X X X X X X b （4—4） 建立方程组 AR=b （4—5） 可求得R 时，纯量（)()(b AR b AR T
--取最小值，把变量的纯量函数取极值的定义化为方程：
0=-b A AR A T T （4—6）
解方程
b A A A R T T 1)(-= （4—7）
最后解得R ， 可以求得7个转换参数。 方案一：
取点HN15,点HN16,HN18作为解算点，点HN21,HN22，HN25，ＨＮ２６作
为检核点。把三个点的空间直角坐标系下的坐标代入以上公式，可以求得WGS84坐标转换到BJ54坐标的7个参数为：
表4-5 方案一求得的七参数
有了坐标转换的7个参数，就可以算出其他点转换后的坐标，以及与BJ54坐标的差值结果如下表：
表4-6 方案一ＷＧＳ８４转换为ＢＪ５４的成果和坐标残差
方案二：
取点HN15,HN16,HN18,HN21作为解算点，点HN21,HN22，HN25，ＨＮ２６作为检核点。把三个点的空间直角坐标系下的坐标代入以上公式，可以求得WGS84坐标转换到BJ54坐标的7个参数为：
表4-7 方案二求得的七参数
表4-8 方案二ＷＧＳ８４转换为ＢＪ５４的成果和坐标残差
方案三：
取点HN15,点HN16,HN18,HN21,HN22作为解算点，点HN21,HN22，HN25，ＨＮ２６作为检核点。把三个点的空间直角坐标系下的坐标代入以上公式，可以求得WGS84坐标转换到BJ54坐标的7个参数为：
表4-8 方案三求得的七参数
表4-9 方案三ＷＧＳ８４转换为ＢＪ５４的成果和坐标残差
方案四：
取点HN15,HN16,HN18,HN21，HN22,HN25作为解算点，点HN21,HN22，HN25，ＨＮ２６作为检核点。把三个点的空间直角坐标系下的坐标代入以上公式，可以求得WGS84坐标转换到BJ54坐标的7个参数为：
表4-10 方案四求得的七参数
表4-11 方案四ＷＧＳ８４转换为ＢＪ５４的成果和坐标残差
方案五：
取点HN15,HN16,HN18,HN21，HN22，HN25，HN26作为解算点，点HN21，HN22，HN25，ＨＮ２６作为检核点。把三个点的空间直角坐标系下的坐标代入以上公式，可以求得WGS84坐标转换到BJ54坐标的7个参数为：
表4-12 方案五求得的七参数
表4-13 方案五ＷＧＳ８４转换为ＢＪ５４的成果和坐标残差
第三节 三参数转换
4.1.2三参数模型的计算
三参数算法取最简单的数学模型，仅考虑坐标原点的平移参数，为了增进三参数法在我国地域内的换算精度，把平移参数的取值作些调整，适当补足简化数学模型时的精度丢失，称之为综合平移参数，其转换数学模型为:
??????
????+??????????=??????????DZ DY DX Z Y X Z Y X 848484545454 （4—8） 式中：DX 、DY 、DZ 为综合平移参数改正。
显然 8454X X DX -=；8454Y Y DY -=；8454Z Z DZ -=。 （4—9） 取四个点（HN16,HN18,HN22,HN26）由公式(4—9)可计算出四个位置的DX ，
DY ，DZ 把计算的结果直接填入表
表4-14 转换坐标与原坐标的差值
这里4组的DX ，DY ，DZ ，它们是非常接近。并把这4组的DX ，DY ，DZ 取平均 值做为最后的综合平移参数改正。
得DX=145.9836118; DY=97.7340275; DZ=7.3125735; 单位为米。 所以具体的三参数模型为: