椭球点空间直角坐标转换的研究

王仲锋,刘 凯,初凤婷

(1.长春工程学院 勘查与测绘工程学院，吉林 长春 130021;2.中水东北勘测设计研究有限责任公司,吉林 长春 130021)

目前，测量上用于地面点坐标转换的模型主要有平面四参数模型、布尔沙(Bursa)模型、三维七参数模型、二维七参数模型和多元回归模型等[1-2]。平面四参数模型是一种平面相似变换模型，当重合点所在坐标系的椭球长半轴和扁率差异不大、转换坐标的区域不大且不存在跨带问题时使用。

Bursa模型也是一种相似变换模型，需要知道重合点在源坐标系和目标坐标系中的大地高。但是目前使用的1954年北京坐标系、1980西安坐标系和2000国家大地坐标系中，只有2000国家大地坐标系中的点可以精确获得大地高。王解先研究认为，对于地面上100 km×100 km的范围，即使公共点中地方坐标的高程存在误差，在求得转换参数后，转换出来的点平面坐标仍基本不变，当多个公共点地方坐标高程有误差时，该结论仍然成立[3]。张锐研究认为：在应用“大地高不确定性主要表现在对高程的影响，而对平面坐标影响较小”这一结论时还是有一定的局限性[4]。张勤和王利研究认为，高程异常误差对坐标的影响随经度减小而增大；区域平均高程面的高度对平面坐标转换精度随高度增大而降低[5]。针对三维空间和高程未知的二维空间之间的七参数的求解问题，谢鸣宇、姚宜斌提出高程趋近法[6],甄佳宁、吴琼、康靖玉、杨国东提出基于曲面逼近的空间坐标转换方法[7]。

高程趋近法和曲面逼近法的思路基本相似，前者取公共点源坐标的大地高初值与目标坐标的初值一致，后者取源坐标与目标坐标的大地高初值均为0。高程趋近法和曲面逼近法均认为通过迭代可求解高精度参数，但在多次迭代后发现，每次迭代后所求参数并无明显变化。刘根友、朱耀仲、朱才连为解决三维空间和大地高未知(但高程已知)的二维空间之间的七参数的求解问题，提出用高程加高程异常代替大地高，带入Bursa模型，将高程异常作为未知参数和七参数一并求解的方法[8]。该方法忽略将高程异常作为未知参数和七参数一并求解带来的另一个问题，即病态方程问题，使得解算结果极不稳定、极不可靠，有时可能无解。

二维七参数模型是在大地微分公式(三维七参数模型)基础上，舍弃其中的大地高和大地高微分式而得来的一种用于不同坐标系下经纬度之间转换的模型[1-2]。该理论极其不严密，存在用较大的增量代替微分的问题(相关研究成果将专文发表)。在实践上，该公式也较为复杂，误差方程呈病态现象也很严重。

多元回归模型分为以平面坐标为变量的多元回归模型和以经纬度为变量的多元回归模型，前者适用于小区域不跨带的平面坐标转换，后者适用于省级及全国区域的坐标转换问题。从实践来看，多元回归坐标转换模型是高精度的转换模型，但并未较好地被接受。

本研究将通过推导，证明不同坐标系下的同名椭球点之三维空间直角坐标之间，利用“布尔沙”模型进行转换，并叙述转换的方法步骤及应用实例。

1 转换原理

由地面点的经度L、纬度B和大地高H可算得该点在三维空间直角坐标系中的坐标X，Y，Z，即:

(1)

当H=0时，即为地面点沿着椭球法线方向在椭球面上的投影点(简称椭球点)的空间直角坐标，简称椭球点空间直角坐标，记为标X0，Y0，Z0，即有：

(2)

其中：

(3)

e2=(a2-b2)/a2.

(4)

根据空间直角坐标转换的“布尔沙”模型，点在空间直角坐标系“1”中的空间直角坐标通过3个平移参数ΔX0，ΔY0，ΔZ0，3个旋转参数ωX，ωY，ωZ和1个尺度比参数m可转换为空间直角坐标系“2”中的空间直角坐标，即有：

(5)

或：

(6)

利用n个公共点的三维空间直角坐标公式(6)，列出3n个误差方程，即可利用最小二乘法解出7个坐标转换参数。但是，当点在空间直角坐标系“1”和“2”中的大地高均未知时，通常认为“布尔沙”模型不再适用，因此有学者推导出了所谓“二维七参数模型”，用于同名点不同坐标系之间经纬度转换。其中，“二维七参数模型”在理论上是不严密的，一是存在用较大的增量代替微分的问题，二是存在舍弃大地高的问题。

以下讨论当点在空间直角坐标系“1”和“2”中的大地高均未知时如何利用“布尔沙”模型求解坐标转换参数，实现同名点在不同坐标系之间经纬度的转换问题。在以下的推导中，假定同名点在“1”和“2”坐标系中的经度和纬度差极小(实际上均在10秒以下)。

由式(1)，(2)知：

(7)

记：

(8)

或

(9)

则：

(ΔX12)2+(ΔY12)2+(ΔZ12)2=(H2cosB2cosL2-H1cosB1cosL1)2+

(H2cosB2sinL2-H1cosB1sinL1)2+(H2sinB2-H1sinB1)2=

而：

(cosB2cosL2cosB1cosL1+cosB2sinL2cosB1sinL1+sinB2sinB1)=

cosB2cosB1(cosL2cosL1+sinL2sinL1)+sinB2sinB1=

cosB2cosB1cos(L2-L1)+sinB2sinB1≈cosB2cosB1+sinB2sinB1=

cos(B2-B1)≈1.

故有：

(ΔX12)2+(ΔY12)2+(ΔZ12)2=(H2-H1)2.

(10)

取H1=H2=0时，式(6)可改写为：

(11)

显然，利用n个公共点在“1”、“2”空间直角坐标系中于椭球面上的投影点的三维空间直角坐标和式(11)，列出3n个方程，如果不考虑观测误差的话，利用最小原理求解的7个坐标转换参数为满足[ΔXΔX+ΔYΔY+ΔZΔZ]=min=0的解。实际上，由于H1=H2=0时已经使得ΔX12=ΔY12=ΔZ12=0，故利用最小二乘原理求解的7个坐标转换参数实际上就是满足误差平方和等于最小的解。

由此可见，同名点在不同空间直角坐标系下沿法线方向投影到椭球面上的点的空间直角坐标之间的转换，同样可使用“布尔沙”模型。

当H1=H2≠0时，如果不考虑观测误差的话，利用最小二乘原理求解的7个坐标转换参数为满足[ΔXΔX+ΔYΔY+ΔZΔZ]≈0的解。

由式(2)知，椭球点的空间直角坐标，求出点的经度和纬度，即：

(12)

因此，同名点在不同空间直角坐标系下，沿法线方向投影到椭球面上的点的空间直角坐标之间的转换，实际上也是同名点在不同坐标系下的经纬度之间的转换。

2 转换参数计算步骤

1)选取n个公共点，要求点位均匀分布地覆盖测区，且点数不少于6个；

2)计算公共点在“1”，“2”空间直角坐标系中的椭球点坐标，即：

(i=1,2,…,n).

3)列误差方程：

(12)

(13)

解式(13)，便可得要解算的7个参数。

(14)

3 算 例

如图1所示，香港地区74个点(覆盖面积月1 000 km2)，其在HK80坐标系和WGS84坐标系下的椭球点三维空间直角坐标见表1(根据香港地政署测绘处网站下载资料计算与整理)，试求HK80坐标系下的椭球点三维空间直角坐标向和WGS84坐标系下的椭球点三维空间直角坐标转换的7个参数。

图1 香港地区74个控制点分布示意图

表1 香港地区74个控制点的椭球点之三维空间直角坐标(节选)

解：为了证明本文所述方法的可行性，取表中63-74号点作为建模点(如图1中“口”形符号所示)，1-62号点作为检验点(如图1中“◇”形符号所示)。

1)解算结果。利用中心化坐标求解，得3个平移参数-157.095 6 m、-267.509 8 m和-187.720 0 m；3个旋转参数分别为2.863 667 s、7.079 517 s、-6.658 663 s；尺度比为-0.604 489 ppm。

2)内符合精度。反应建模内符合精度的点位残差见图2。其中空间点位差最小者为70号点，差值为0.10 cm，最大者为64号点，差值为0.64 cm，平均为0.38 cm。

图2 内符合检验空间点位差散点

3)外符合精度。利用61个检验点进行外符合检验，相应的空间点位差分布情况见图3。其中空间点位差最小者为49号点，差值为0.07 cm，最大者为53号点，差值为1.26 cm，平均为0.46 cm。

图3 外符合检验空间点位差散点

4 结 论

本文推导不同坐标系下同名椭球点的三维坐标之间相互转换的原理，叙述转换的步骤和方法，并用实例验证了其可行性。通过研究分析可得如下结论：

1)用“布尔沙”模型进行椭球点的三维坐标相互转换理论严密、方法可行，适用于任何大小的区域，可代替“二维七参数模型”通过椭球点的三维空间直角坐标转换来间接实现不同坐标系下地面点的经纬度转换；

2)取同名点在“1”、“2”坐标系下的大地高H1=H2时(如H1代表北京54或西安80坐标系下的大地高(通常为未知值)，H2代表CGCS2000坐标系下的大地高(为已知值)，可取H54 or 80=H2 000)，其三维空间直角坐标也可用“布尔沙”模型进行相互转换，但如此转换带有一点“近似”的性质。