一、坐标系转换

1.1 测站坐标系rae（极坐标系）和空间直角坐标系xyz（东北天坐标系）互转

1.1.1、原理

1.1.2、公式

1.1.3、python代码

1.2 空间直角坐标系xyz（东北天坐标系）和地心非惯性系直角坐标系XYZ互转

1.2.1、原理

1.2.2、公式

1.2.3、python代码

1.3 经纬高LBH和地心非惯性系直角坐标系XYZ互转

1.3.1、原理

1.3.2、公式

1.3.3、python代码

1.4 测站坐标系rae（极坐标系）和经纬高LBH互转

1.4.1、python代码

二、参考链接

一、坐标系转换

1.1 测站坐标系rae（极坐标系）和空间直角坐标系xyz（东北天坐标系）互转

1.1.1、原理

这个原理过于简单，此处不进行详述。要注意的是我这个地方a（方位）是从正北方向为0°，顺时针为正，e（俯仰）水平面朝上为0~90°，水平面朝下为0~-90°。x轴指向东，y轴指向北，z轴指向天。

1.1.2、公式

a）rae转xyz

$x = r * cos(e) * sin(a)$

$y = r * cos(e) * cos(a)$
$z = r * sin(e)$

b）xyz转rae

$r = \sqrt{x^{2} + y^{2} + z^{2}}$
$a = arctan(x/y)$
$e = arctan(z / \sqrt{x^{2} + y^{2}})$

1.1.3、python代码

import math
import numpy as npdef rae2xyz(r: float, a: float, e: float):"""rae转xyz（测站xyz）:param r: 距离，单位m:param a: 方位，单位（°），从正北顺时针计算0~360°:param e: 俯仰，（°）水平面朝上为0~90°，朝下为0~-90°:return: x,y,z(m)，北为y轴，东为x轴。水平面朝上为z轴"""# 转为弧度a = a * math.pi/180e = e * math.pi/180x = r * math.cos(e) * math.sin(a)y = r * math.cos(e) * math.cos(a)z = r * math.sin(e)return x, y, zdef xyz2rae(x: float, y: float, z: float):"""xyz转rae（测站xyz）:param x: 东为x轴，单位m:param y: 北为y轴，单位m:param z: 水平面朝上为z轴，单位m:return: r(m)、a、e。方位，单位（°），从正北顺时针计算0~360°,俯仰，（°）水平面朝上为0~90°，朝下为0~-90°。距离单位（m）"""r = math.sqrt(math.pow(x, 2) + math.pow(y, 2) + math.pow(z, 2))a = math.atan2(x, y) / math.pi * 180e = math.atan2(z, math.sqrt(math.pow(x, 2) + math.pow(y, 2))) / math.pi * 180return r, a, e

1.2 空间直角坐标系xyz（东北天坐标系）和地心非惯性系直角坐标系XYZ互转

1.2.1、原理

空间直角坐标系的转换主要涉及到平移和旋转。
        （1）旋转
        旋转主要使用欧拉角公式，任何两个同原点的空间直角坐标系都可以通过三次旋转进行转换，分别是绕X轴、绕Y轴、绕Z轴旋转。
        绕X轴逆时针的转换为：

$\begin{bmatrix} x'\\ y'\\ z' \end{bmatrix} = \begin{bmatrix} 1 &0 &0 \\ 0&cos\alpha & sin\alpha\\ 0& -sin\alpha & cons\alpha \end{bmatrix} * \begin{bmatrix} x\\ y\\ z \end{bmatrix}$ ，转换矩阵： $C_{x} = \begin{bmatrix} 1 &0 &0 \\ 0&cos\alpha & sin\alpha\\ 0& -sin\alpha & cons\alpha \end{bmatrix}$

绕y轴逆时针的转换为：

$\begin{bmatrix} x'\\ y'\\ z' \end{bmatrix} = \begin{bmatrix} cos\beta &0 &-sin\beta \\ 0&1 & 0\\ sin\beta& 0 &cos\beta \end{bmatrix} * \begin{bmatrix} x\\ y\\ z \end{bmatrix}$ ，转换矩阵为： $C_{y} = \begin{bmatrix} cos\beta &0 &-sin\beta \\ 0&1 & 0\\ sin\beta& 0 &cos\beta \end{bmatrix}$

绕z轴逆时针的转换为：

$\begin{bmatrix} x'\\ y'\\ z' \end{bmatrix} = \begin{bmatrix} cos\gamma &sin\gamma & 0\\ -sin\gamma &cos\gamma & 0\\ 0&0 &1 \end{bmatrix} * \begin{bmatrix} x\\ y\\ z \end{bmatrix}$ ，转换矩阵为： $C_{z} = \begin{bmatrix} cos\gamma &sin\gamma & 0\\ -sin\gamma &cos\gamma & 0\\ 0&0 &1 \end{bmatrix}$

实际在转换时，根据绕x、y、z轴旋转的先后顺序调用计算即可。设初始点为P，那么经过旋转后我们得到了P'点。

（2）平移

经过上一步我们已经完成了旋转，再经过一次平移就能够完成空间中任意两个直角坐标系的转换。

经过第一个步骤旋转之后我们得到了P'点，要得到P'点在x''y''z''坐标系下的坐标，则计算出O''P'向量值即可。平移公式如下：

$\underset{O''P'}{\rightarrow} = \underset{O''O'}{\rightarrow} + \underset{O'P'}{\rightarrow}$

1.2.2、公式

a）xyz转XYZ

sX、sY、sZ代表测站在XYZ坐标系下的位置，这是个已知值。

$\begin{bmatrix} X\\ Y\\ Z \end{bmatrix} = \begin{bmatrix} sX\\ sY\\ sZ \end{bmatrix} + \begin{bmatrix} -sin(L)&-cos(L)*sin(B) &cos(L)*cos(B) \\ cos(L) &-sin(L)*sin(B) &sin(L)*cos(B) \\ 0&cos(B) &sin(B) \end{bmatrix}* \begin{bmatrix} x\\ y\\ z \end{bmatrix}$

b）XYZ转xyz

$\begin{bmatrix} x\\ y\\ z \end{bmatrix} = \begin{bmatrix}-sin(L) &cos(L) & 0\\ -sin(B)*cos(L)&-sin(B)*sin(L) &cos(B) \\ cos(B)*cos(L)&cos(B)*sin(L) &sin(B) \end{bmatrix}*\begin{bmatrix} X\\ Y\\ Z \end{bmatrix} - \begin{bmatrix} sX\\ sY\\ sZ \end{bmatrix}$

1.2.3、python代码

import math
import numpy as npdef xyz2XYZ(x:float, y:float, z:float, L:float, B:float, H:float, coorSys:str):"""测站xyz（东北天坐标系转地心非惯性直角坐标系），测站LBH是在什么坐标系的，转出来的XYZ就是相应坐标系下的。比如我测站是在wgs84系的，那么我转出来的XYZ就是wgs84系下的XYZ，因为不同坐标系椭球模型本身的差异，会导致转出来的XYZ有相应的差异:param x: 东为x轴，单位m:param y: 北为y轴，单位m:param z: 水平面朝上为z轴，单位m:param L: 测站经度,单位°:param B: 测站纬度（地理纬度），单位°:param H: 测站高度，单位m:param coorSys: 测站LBH所属坐标系:return: 地心空间直角坐标系XYZ，单位m"""sX, sY, sZ = LBH2XYZ(L, B, H, coorSys)  # 获取测站位置XYZ# 转为弧度L = L/180 * math.piB = B/180 * math.piMs2e = np.mat([[-math.sin(L), -math.cos(L)*math.sin(B), math.cos(L)*math.cos(B)],[math.cos(L),  -math.sin(L)*math.sin(B), math.sin(L)*math.cos(B)],[0, math.cos(B), math.sin(B)]])eP = Ms2e * np.mat([[x], [y], [z]]) + np.mat([[sX], [sY], [sZ]])X = eP[0, 0]Y = eP[1, 0]Z = eP[2, 0]return X, Y, Zdef XYZ2xyz(X:float, Y:float, Z:float, L:float, B:float, H:float, coorSys:str):"""地心空间直角坐标系XYZ转测站xyz(东北天坐标系，东为x轴，北为y轴，水平面朝上为z轴)，这个XYZ是什么坐标系下的XYZ，测站位置就应当用相应的LBH，否则转换就不正确。:param X: X轴刻度值，单位m。指向0°经线与协议赤道焦点:param Y: Y轴刻度值，单位m。Y轴与X、Z成右手正交，由Z到X大拇指指向方向为Y:param Z: Z轴刻度值，单位m。指向协议北极:param L: 测站经度,单位°:param B: 测站纬度（地理纬度），单位°:param H: 测站高度，单位m:param coorSys: 测站LBH所属坐标系:return: 测站xyz，单位m"""sX, sY, sZ = LBH2XYZ(L, B, H, coorSys)  # 获取测站位置XYZ# 转为弧度L = L/180 * math.piB = B/180 * math.piMe2s = np.mat([[-math.sin(L), math.cos(L), 0],[-math.sin(B)*math.cos(L), -math.sin(B)*math.sin(L), math.cos(B)],[math.cos(B)*math.cos(L), math.cos(B)*math.sin(L), math.sin(B)]])sP = Me2s * (np.mat([[X], [Y], [Z]]) - np.mat([[sX], [sY], [sZ]]))x = sP[0, 0]y = sP[1, 0]z = sP[2, 0]return x, y, z

1.3 经纬高LBH和地心非惯性系直角坐标系XYZ互转

1.3.1、原理

如上图所示，已知P'点的L、B、H求X、Y、Z。画出P点所在子午圈平面图如下：

P点子午圈平面图

如上图所示，P(v, z)，半长轴a，半短轴b，有

$\frac{v^{2}}{a^{2}} + \frac{z^{2}}{b^{2}} = 1$ ①椭圆公式

对①式进行微分得：

$\frac{dz}{dv} = -\frac{b^{2}v}{a^{2}z} = -\frac{1}{tan(B)}$ ②切线斜率与纬度的关系

由②式得：

$z = \frac{b^{2}}{a^{2}}tan(B)v$ ③

又有偏心率公式：

$e = \frac{\sqrt{a^{2} - b^{2}}}{a}$ ④偏心率公式

将④式代入③式得：

$z= v(1 - e^{2})tan(B)$ ⑤

将⑤式代入①式得（化简过程省略）：

$v = \sqrt{\frac{a^{2}cos^{2}(B)}{1-e^{2}sin^{2}(B)}} = \frac{a}{\sqrt{1-e^{2}sin^{2}(B)}}cos(B) = MP*cos(B)$ ⑥

令:

$N = \frac{a}{\sqrt{1-e^{2}sin^{2}(B)}}$ ⑦

可知MP = N，综上所述，可知P点在空间直角坐标系下的坐标（x,y,z）：

$x = v * cos(L) = N * cos(B) * cos(L)$
$y = v * sin(L) = N * cos(B) * sin(L)$
$z= v(1 - e^{2})tan(B) = Ncos(B)(1 - e^{2})tan(B) = N(1 - e^{2})sin(B)$

以上过程我们求得了椭圆上一点P的x、y、z值，但我们实际是要求椭圆外一点P'的X、Y、Z值。由“P点子午圈平面图”我们不难看出，P'点的X、Y、Z值也就是OP'向量的值，又：

$\underset{OP'}{\rightarrow} = \underset{OP}{\rightarrow} + \underset{PP'}{\rightarrow}$ ⑧

有：

$\underset{OP}{\rightarrow} = \begin{bmatrix} Ncos(B)cos(L)\\ Ncos(B)sin(L)\\ N(1-e^{2})sin(B) \end{bmatrix}$
$\underset{PP'}{\rightarrow} = H * \begin{bmatrix} cos(B)cos(L)\\ cos(B)sin(L)\\ sin(B) \end{bmatrix}$

代入⑧式，得：

$\underset{OP'}{\rightarrow} = \begin{bmatrix} Ncos(B)cos(L)\\ Ncos(B)sin(L)\\ N(1-e^{2})sin(B) \end{bmatrix} + H * \begin{bmatrix} cos(B)cos(L)\\ cos(B)sin(L)\\ sin(B) \end{bmatrix} = \begin{bmatrix} (N+H)cos(B)cos(L)\\ (N+H)cos(B)sin(L)\\ (N(1-e^{2}) + H)sin(B) \end{bmatrix}$

由XYZ转LBH过程如下：

$L = arctan(Y/X)$

做辅助线画图如下：

所以，

$tan(B) = \frac{NP'}{MN} = \frac{Z + OM}{\sqrt{X^{2} + Y^{2}}}$ ⑨

又由⑤式知，

$z= N(1 - e^{2})sin(B) = (N - Ne^{2})sin(B) = (MP -MO' )sin(B)$

所以，

$MO' = Ne^{2}$

$OM = MO'*sin(B) = Ne^{2}sin(B)$ ⑩

所以将⑩式代入⑨式，得：

$B=arctan((Z + N*e^{2}sin(B)) / \sqrt{X^{2} + Y^{2}})$ //迭代求解

$H = O'P' - (MP - MO')=Z / sin(B) - N*(1 - e^{2})$

1.3.2、公式

a）LBH转XYZ：

a为半长轴，e为偏心率

$N = a / \sqrt{1 - e^{2}*sin^{2}B}$
$X = (N + H)*cosB*cosL$
$Y = (N + H)*cosB*sinL$
$Z = [N*(1 - e^{2}) + H]*sinB$

b）XYZ转LBH：

a为半长轴，e为偏心率

$N = a / \sqrt{1 - e^{2}*sin^{2}B}$
$L = arctan(Y/X)$
$B=arctan((Z + N*e^{2}sin(B)) / \sqrt{X^{2} + Y^{2}})$ //此处为超越方程，需要迭代求解
$H = Z / sin(B) - N*(1 - e^{2})$

1.3.3、python代码

import math
import numpy as npdef LBH2XYZ(L:float, B:float, H:float, coorSys:str):"""经纬高转地心XYZ（X指向0°经线与协议赤道焦点，Z指向协议北极，Y轴与X、Z成右手正交，由Z到X大拇指指向方向为Y）:param L: 经度,单位°:param B: 纬度（地理纬度），单位°:param H: 高度，单位m:param coorSys 坐标系选择，参数可以为"wgs84"、"cgcs2000"，两种坐标系所采取的椭圆半长轴均为6378137m，扁率不同，wgs84的扁率为1/298.257223564，半短轴为6356752.314245251，偏心率e为0.08181919084248535cgcs2000的扁率为1/298.257222101，半短轴为6356752.314140356，偏心率e为0.08181919104281517:return: 地心XYZ 单位m"""L = L/180 * math.pi   # 先转为弧度B = B/180 * math.pi   # 先转为弧度a = 6378137   # 半长轴if coorSys == "wgs84":e = 0.08181919084248535elif coorSys == "cgcs2000":e = 0.08181919104281517N = a / math.sqrt(1 - math.pow(e, 2) * math.pow(math.sin(B), 2))X = (N + H) * math.cos(B) * math.cos(L)Y = (N + H) * math.cos(B) * math.sin(L)Z = (N * (1 - math.pow(e, 2)) + H) * math.sin(B)return X, Y, Zdef XYZ2LBH(X:float, Y:float, Z:float, coorSys:str):"""地心非惯性坐标系（地固系）XYZ转换为L（经度）B（纬度）H（高程）:param X: X轴刻度值，单位m。指向0°经线与协议赤道焦点:param Y: Y轴刻度值，单位m。Y轴与X、Z成右手正交，由Z到X大拇指指向方向为Y:param Z: Z轴刻度值，单位m。指向协议北极:param coorSys: 坐标系选择，参数可以为"wgs84"、"cgcs2000"，:return: LBH  单位°、m"""a = 6378137  # 半长轴if coorSys == "wgs84":e = 0.08181919084248535elif coorSys == "cgcs2000":e = 0.08181919104281517L = math.atan2(Y, X) / math.pi * 180tB = 0N = a / math.sqrt(1 - math.pow(e, 2) * math.pow(math.sin(tB), 2))B = math.atan2((Z + N * math.pow(e, 2) * math.sin(tB)), math.sqrt(X*X + Y*Y))while math.fabs(B - tB) > 1e-16:  # 超越方程迭代计算tB = BN = a / math.sqrt(1 - math.pow(e, 2) * math.pow(math.sin(tB), 2))B = math.atan2((Z + N * math.pow(e, 2) * math.sin(tB)), math.sqrt(X * X + Y * Y))N = a / math.sqrt(1 - math.pow(e, 2) * math.pow(math.sin(B), 2))if B == 0:  # 防止赤道上目标刚好为纬度0°时，无法计算的问题H = math.sqrt(X * X + Y * Y) / math.cos(B) - Nelse:H = Z / math.sin(B) - N * (1 - e * e)B = B / math.pi * 180return L, B, H

1.4 测站坐标系rae（极坐标系）和经纬高LBH互转

1.4.1、python代码

结合上文描述的几个函数使用。

def rae2LBH(r: float, a: float, e: float, sL:float, sB:float, sH:float, coorSys:str):"""rae转LBH:param r: 距离，单位m:param a: 方位，单位（°），从正北顺时针计算0~360°:param e: 俯仰，（°）水平面朝上为0~90°，朝下为0~-90°:param sL: 测站经度,单位°:param sB: 测站纬度（地理纬度），单位°:param sH: 测站高度，单位m:param coorSys: 坐标系选择，参数可以为"wgs84"、"cgcs2000"，测站所属什么坐标系，这个地方就应当使用什么坐标系:return: L(经度)、B(纬度)、H（高度） 单位°、m"""x, y, z = rae2xyz(r, a, e)X, Y, Z = xyz2XYZ(x, y, z, sL, sB, sH, coorSys)L, B, H = XYZ2LBH(X, Y, Z, coorSys)return L, B, Hdef LBH2rae(L: float, B: float, H: float, sL:float, sB:float, sH:float, coorSys:str):"""LBH转rae:param L: 经度,单位°:param B: 纬度（地理纬度），单位°:param H: 高度，单位m:param sL: 测站经度,单位°:param sB: 测站纬度（地理纬度），单位°:param sH: 测站高度，单位m:param coorSys: 坐标系选择，参数可以为"wgs84"、"cgcs2000"，测站所属什么坐标系，这个地方就应当使用什么坐标系:return: r(m)、a、e。方位，单位（°），从正北顺时针计算0~360°,俯仰，（°）水平面朝上为0~90°，朝下为0~-90°。距离单位（m）"""X, Y, Z = LBH2XYZ(L, B, H, coorSys)x, y, z = XYZ2xyz(X, Y, Z, sL, sB, sH, coorSys)r, a, e = xyz2rae(x, y, z)return r, a, e

二、参考链接

此处可参考链接：

（1）第四章地球椭球及其数学详解

（2）第四章椭球数学变换16节

（3）大地测量学基础[1](7)(控制)

测站坐标系、地心非惯性系、经纬高互转相关推荐

经纬高坐标系-ECEF坐标系-ENU坐标系
无人机搭载的RTK获得的经纬高坐标要转换为东北天坐标,才能用于局部的导航和定位.为了这个目的,查阅资料,越查越懵逼,竟然这么多的坐标系,略懂之后,将学到的信息记录如下,很多跟我的目的:"RT ...
经纬高坐标系转到东北天坐标系
经纬高坐标系转到东北天坐标系基本思路:首先把经纬高(大地坐标系.lla.llh)转到直角坐标系(地心地固直角坐标系(ECEF).xyz),然后再转为局部坐标系下(东北天坐标系.以第一点作为东北天坐标 ...
地心直角坐标系转经纬高
/**@func TheGeoCoordinateSysRotatesToLonAndLat*@brief 地心直角坐标系转经纬高*@input*@output*@return 经纬高*@author ...
经纬高坐标转东北天坐标
目录 1 问题描述 2 解决方案 2.1 经纬高转ECEF 2.2 ECEF转东北天 2.3 代入求解 2.3.1 东向化简 2.3.2 北向化简 2.3.3 天向化简 2.4 总结 1 问题描述已 ...
LLA（经纬高）坐标转换成ENU（东北天）坐标的详细推导
这是一篇经纬高(LLA)坐标转东北天坐标(ENU)的详细推导,并给出近似转换的过程和结果参考资料: https://blog.csdn.net/qq_34213260/article/details ...
已知经纬高求其在无人机图像中的像素点坐标
背景:已知无人机拍摄时刻的经纬高(WGS84坐标系)以及物点的经纬高,求解该物点在图像中的像素点位置三种像空间坐标系定义: ①德国汉诺威大学(Hannover)定义的像平面坐标系BLUH:原点位于像 ...
OpenCV图像坐标系与行列宽高的关系
这篇文章挺好 OpenCV图像坐标系与行列宽高的关系图片坐标系,与从小到大见到的xy坐标系,x轴方向相同,只是y轴方向相反.
站心直角坐标系转经纬高
序:网上很多算法,基本都是错的,还得去找论文才行.原理就不讲了,我会贴出论文引用. 1.转换过程 (1)确定站心直角坐标系的原点(X0,Y0,Z0)=(0,0,0)的经纬高(B0,L0,H0)作为基准 ...
经纬高、ECEF、ENU坐标系相互转换与实现
ENU与ecef与WGS84相互转换(c++和python) - CodeAntenna 原理: 坐标系之间的转换关系(ECEF.LLA.ENU)_可即的博客-CSDN博客_lla坐标系

测站坐标系、地心非惯性系、经纬高互转

一、坐标系转换

1.1 测站坐标系rae（极坐标系）和空间直角坐标系xyz（东北天坐标系）互转

1.1.1、原理

1.1.2、公式

1.1.3、python代码

1.2 空间直角坐标系xyz（东北天坐标系）和地心非惯性系直角坐标系XYZ互转

1.2.1、原理

1.2.2、公式

1.2.3、python代码

1.3 经纬高LBH和地心非惯性系直角坐标系XYZ互转

1.3.1、原理

1.3.2、公式

1.3.3、python代码

1.4 测站坐标系rae（极坐标系）和经纬高LBH互转

1.4.1、python代码

二、参考链接

测站坐标系、地心非惯性系、经纬高互转相关推荐

最新文章

热门文章