GNSS中常见坐标转换

GNSS中常用的4种坐标转换：大地坐标与地心地固坐标之间的相互转换、地心地固坐标与站心地平坐标的相互转换

本文内容包括对上述每种坐标转换程序的设计思路、预期功能、算例及结果分析的阐述，以及笔者对编程过程中一些常见问题和注意事项的总结。

代码详见：

4种坐标转换MATLAB程序资源-CSDN文库https://download.csdn.net/download/m0_58307078/88583888

1 实验内容及目的

1. 了解常用坐标系统的基本知识；
2. 掌握地心地固坐标系下坐标与大地坐标系下坐标之间相互转换的方法，编写程序进行实现；
3. 了解并掌握地心坐标与站心坐标间的相互转换。

2 数据来源及编程测试环境

1．算例数据来源：PPP精密单点定位原理实验平台

2．编程环境：MATLAB 2016a

3．测试环境：MATLAB 2016a及PPP精密单点定位原理实验平台

3 实验内容及成果

3.1 大地坐标与地心地固坐标之间的相互转换

3.1.1 大地坐标转地心地固坐标——BLH2XYZ(pos)

1．程序设计

预期功能：将输入参数pos（B,L,H ）转换为X,Y,Z 并输出。

设计思路：按B,L,H 与X,Y,Z 的换算关系编写函数，程序计算流程及步骤如图3-1所示：

图3-1 大地坐标转地心地固坐标计算流程示意

2．遇到的问题及解决过程

问题描述：由于椭球参数长半轴a和偏心率e选取差异，导致计算结果出现差异，超出PPP精密单点定位原理实验平台正确答案的范围。

解决过程：查阅资料时发现，椭球长半轴a的取值有两种说法：1） a=6378137；2）a=6378140。使用第二种取值计算时发现许多次结果验证错误，只有偶然的几次是正确的，于是尝试第一种取值，即a=6378137，经验证，计算结果正确。

注意事项：使用sind()、cosd()函数可直接输入角度进行计算，无需将角度转化为弧度，使用方便。

具体代码实现如下图所示：

图3-2 大地坐标转地心地固坐标代码截图

3．算例结果

大地坐标转地心地固坐标的算例结果如下表所示：

表3-1 大地坐标转地心地固坐标算例结果

B	L	H	X	Y	Z
32.55165258	-76.05517747	497.85	1296948.55	-5223200.13	3412420.66
-5.70390942	-99.08322398	884.85	-1002096.49	-6268048.78	-629773.15
-14.65804010	17.00186163	86.67	5902212.93	1804697.28	-1603545.47
43.65184484	62.08014161	252.91	2164396.15	4084409.15	4380358.41
-46.55232398	-58.19086276	933.62	2316352.24	-3734563.76	-4608360.04

3.1.2地心地固坐标转大地坐标——XYZ2BLH(coord)

1．程序设计

预期功能：将输入参数coord（X,Y,Z ）转换为B,L,H 并输出。

设计思路：按X,Y,Z 与B,L,H 的换算关系编写函数，值得注意的是：大地纬度的计算需要通过迭代实现，而大地高的计算又依赖于大地纬度B，因此可以先计算大地经度，再通过设置循环迭代计算大地纬度B和卯酉圈半径N。

具体计算流程及步骤如图3-3所示：

图3-3 地心地固坐标转大地坐标计算流程示意

2．遇到的问题及解决过程

问题描述：对于B、H、N，可以参考的计算公式有多种形式（常见的几种见图3-4），不同的计算形式的公式计算出的结果存在差异，影响结果的正确性。

解决过程：参考不同来源的资料，尝试不同形式的公式并验证计算结果，选择合适计算公式进行编程。

注意事项：由于大地经度L在[-180°,180°]，利用反三角函数计算时，使用atan2d()函数，既可以返回正确的L（根据输入参数的符号自动判断象限），又可以直接输入角度计算（无需转换为弧度），使用方便。

图3-4 两种不同形式的计算公式截图

3．算例结果

地心地固坐标转大地坐标的一些算例结果如下表所示：

表3-2 地心地固坐标转大地坐标的算例结果

X	Y	Z	B	L	H
1745851.74	5324845.06	3037595.70	28.62208494	71.84728536	823.31
1296949.16	-5223202.58	3412422.26	-32.41248422	18.53999448	949.61
-5923515.14	-1641755.91	-1700020.69	-15.55894265	-164.50880458	964.63
5434792.50	1661775.91	2885854.23	27.07662761	17.00186166	165.77
337907.31	3810513.32	5087370.76	53.24320097	84.93240142	729.63

3.2 地心地固坐标与站心地平坐标的相互转换

3.2.1 地心坐标转站心坐标——XYZ2NEU(A,B)

1．程序设计

预期功能：通过输入站心A的空间直角坐标X,Y,Z 和待转换位置B的空间直角坐标X,Y,Z ，将B点的X,Y,Z 转换为N,E,U 并输出。

设计思路：按X,Y,Z 与N,E,U 的换算关系编写函数，程序计算流程及步骤如图3-5所示：

图3-5 地心坐标转站心坐标计算流程示意

2．遇到的问题及解决过程

问题描述：计算旋转矩阵时忘记把角度转换为弧度导致计算错误。

解决过程：由于计算结果和正确结果相差极大，起初怀疑计算公式有误，查阅资料后发现公式正确；查验中间步骤（step2）时突然想到step3中没有进行单位换算，发现错误后将旋转矩阵计算用的sin()、cos()函数改为sind()、cosd()函数（函数输入为角度）后经验证，计算结果正确。

具体代码实现如下图所示：

图3-6 地心坐标转站心坐标代码截图

3．算例结果

地心坐标转站心坐标的算例结果及计算误差如表3-3所示：

表3-3 地心坐标转站心坐标算例结果

站心A(BJFS)的空间直角坐标	坐标类型	BJSH坐标	JIXN坐标
X0=-2148744.2580 Y0=4426641.2470 Z0=4044655.8790	X	-2154109.4234	-2259012.3602
	Y	4373150.5330	4333892.0191
	Z	4099357.1061	4084475.2137
	N	71328.1837	53173.9316
	E	28185.1135	139700.9321
	U	-394.1085	-1798.7857

3.2.2 站心坐标转地心坐标——NEU2XYZ(A,B)

1．程序设计

预期功能：通过输入站心A的空间直角坐标X,Y,Z

和待转换位置B的站心地平坐标N,E,U ，将B点的N,E,U 转换为X,Y,Z 并输出。

设计思路：按X,Y,Z 与N,E,U 的换算关系编写函数，具体计算流程及步骤如图3-7所示：

图3-7 站心坐标转地心坐标计算流程示意

注意事项：该变换与地心坐标转站心坐标相逆，计算时需要对旋转矩阵进行转置操作。具体代码实现如下图所示：

图3-8 站心坐标转地心坐标代码截图

2．算例结果

站心坐标转空间直角坐标的一些算例结果及计算误差如表3-4所示：

表3-4 站心坐标转地心坐标的算例结果

站心A(BJFS)的空间直角坐标	坐标类型	BJSH坐标	JIXN坐标
X0=-2148744.2580 Y0=4426641.2470 Z0=4044655.8790	N	71328.1837	53173.9317
	E	28185.1135	139700.9321
	U	-394.1048	-1798.7830
	X	-2154109.4246	-2259012.3611
	Y	4373150.5356	4333892.0209
	Z	4099357.1085	4084475.2156

3.2.3 地心坐标与站心坐标相互转换算例误差

下表为地心地固坐标与站心地平坐标相互转化算例的误差示意：

表3-5 站心坐标转地心坐标的算例结果

类型	BJSH坐标	BJSH计算结果	BJSH误差	JIXN坐标	JIXN计算结果	JIXN误差
N	71328.1837	71328.18367	0.0000	53173.9317	53173.9316	-0.0001
E	28185.1135	28185.1135	0.0000	139700.9321	139700.9321	0.0000
U	-394.1048	-394.1085285	-0.0037	-1798.783	-1798.785745	-0.0027
X	-2154109.423	-2154109.425	-0.0012	-2259012.36	-2259012.361	-0.0009
Y	4373150.533	4373150.536	0.0026	4333892.019	4333892.021	0.0018
Z	4099357.106	4099357.108	0.0024	4084475.214	4084475.216	0.0019

不难看出，各个坐标计算误差极小，说明计算结果正确。

代码详见：4种坐标转换MATLAB程序资源-CSDN文库https://download.csdn.net/download/m0_58307078/88583888