GPS从入门到放弃（三） --- GPS坐标系

GPS从入门到放弃（三） — GPS坐标系

定位就需要坐标，坐标当然是相对坐标系而言的，我们描述一个物体的位置，首先就需要建立坐标系。

按大类来分，坐标系可以分为惯性坐标系和非惯性坐标系。惯性坐标系是在空间静止或者做匀速直线运动的坐标系，其他都是非惯性坐标系。

GPS涉及到的坐标系大体有五个，在说这五大坐标系之前，我们需要先了解一些基本概念。

基本概念

地极：（Polar）地球自转轴与地球表面的两个交点，北边的叫北极，南边的叫南极。赤道面：（Equator Plane）通过地心并于地球自转轴垂直的平面。赤道：（Equator）赤道面与地球表面相交的大圆。天球：（Celestial Sphere）天文学概念，指一个以地心为中心，半径为任意长的假想球体。其目的是将天体沿观测者视线投影到球面上，以便于研究天体及其相互关系。黄道：（Ecliptic）太阳中心在天球上视运动的轨迹。即地球绕太阳公转的轨道平面与地球表面相交的大圆。黄赤交角：（Ecliptic Obliquity）黄道面与赤道面的夹角，约23.5°。春分点：（Vernal Equinox）黄道与赤道有两个交点，其中太阳投影沿黄道从南向北通过赤道的点，称为春分点，另一点为秋分点。岁差：（Axial Precession ）地球自转轴长期进动，引起春分点沿黄道西移，致使回归年短于恒星年的现象。周期约为25800年。主要有日月岁差和行星岁差。章动：（Nutation ）月球在白道上运行，白道与黄道相交成5°9′的角，月球围绕地球公转导致地球在公转轨道上左右摇摆，以18.6年为周期，这种现象称为章动。极移：（Polar Wandering ）地球自转轴相对于地球并不固定，这种运动称地极移动，简称极移。子午面：（Meridian Plane）包含地球自转轴的平面。本初子午面：（Prime Meridian Plane）通过英国伦敦格林尼治天文台与地球自转轴构成的平面，是地球上计算经度的起始经线。

有了这些基本概念后，就可以进一步了解五大坐标系了。

地心惯性坐标系（ECI: Earth Centered Inertial） 地心惯性坐标系是太阳系内的一个惯性坐标系，不随地球而转动，也不受地球、太阳运行的章动和岁差的影响。坐标原点位于地心；X轴位于赤道平面内，指向某一特定年(历元时刻)的太阳春分点位置；Z轴指向那一年地球北极的平均位置处；Y轴位于赤道平面内，与X轴垂直，且与X、Z轴构成右手直角坐标系。由于采用的历元时间不同，可以有各种不同的地心惯性坐标系，目前国际上通用的地心惯性坐标系是J2000历元坐标系，它是以公元2000年的春分点为基准的历元坐标系。 地心地固直角坐标系（ECEF: Earth Centered Earth Fixed） 地固坐标系固定在地球上而随地球一起在空间做公转和自转运动，因此地球上任一固定点在地球坐标系的坐标就不会由于地球旋转而变化。坐标原点位于地心；X轴指向参考子午面与地球赤道的交点；Z轴与地球自转轴重合并指向地球北极；Y轴位于赤道平面内，与X轴垂直，且与X、Z轴构成右手直角坐标系。因为有极移，所以采用了协议地极，以1900年到1905年间的地极实际位置的平均值作为基准点。 大地坐标系：也叫经纬高坐标系（LLA: Longitude Latitude Altitude） 也是地固坐标系。坐标原点位于地心。基于基准椭球体（基准椭球体是定义的与地球几何最吻合的椭球体）。大地纬度 ϕ\phiϕ 是过该点的基准椭球面法线与赤道面的夹角。纬度值在-90°到+90°之间。北半球为正，南半球为负。大地经度 λ\lambdaλ 是过该点的子午面与本初子午面之间的夹角。经度值在-180°到+180°之间。大地高度 hhh 是过该点到基准椭球面的法线距离，基准椭球面以内为负，以外为正。 站心坐标系：也叫东北天坐标系（ENU: East North Up） 是以观测站为原点的坐标系，主要用于了解以观察者为中心的其他物体运动规律。三个坐标轴分别指向相互垂直的东向、北向和天向，因而又称东北天坐标系。可用于计算卫星在用户处的观测向量、仰角和方位角。 WGS-84: World Geodetic System-1984 Coordinate System 是一个地心地固直角坐标系。坐标原点为地心，Z轴指向国际时间服务机构（BIH）1984年定义的协议地球极（CTP: Conventional Terrestrial Pole）方向，X轴指向本初子午面和CTP赤道的交点，Y轴与Z轴、X轴垂直构成右手坐标系。GPS广播星历是以WGS-84坐标系为基准的。

最后列举一下坐标之间的转换关系。

LLA => ECEF

x=(N+h)cos⁡ϕcos⁡λy=(N+h)cos⁡ϕsin⁡λz=[N(1−e2)+h]sin⁡ϕx = (N+h)\cos\phi\cos\lambda \\ y = (N+h)\cos\phi\sin\lambda \\ z = \left[N(1-e^2)+h\right]\sin\phi x=(N+h)cosϕcosλy=(N+h)cosϕsinλz=[N(1−e2)+h]sinϕ

其中eee是椭球偏心率，NNN是基准椭球体的卯酉圈曲率半径。若基准椭球体长半径为aaa，短半径为bbb，则

e=a2−b2a2N=a1−e2sin⁡2ϕe=\sqrt{\frac{a^2-b^2}{a^2}} \\ N=\frac{a}{\sqrt{1-e^2\sin^2\phi}} e=a2a2−b2​​N=1−e2sin2ϕ​a​ECEF => LLA

λ=arctan⁡(yx)h=x2+y2cos⁡ϕ−Nϕ=arctan⁡[zx2+y2(1−e2NN+h)−1]\lambda = \arctan\left(\frac{y}{x}\right) \\ h = \frac{\sqrt{x^2+y^2}}{\cos\phi} - N \\ \phi = \arctan\left[\frac{z}{\sqrt{x^2+y^2}}\left(1-e^2\frac{N}{N+h}\right)^{-1}\right] λ=arctan(xy​)h=cosϕx2+y2​​−Nϕ=arctan[x2+y2​z​(1−e2N+hN​)−1]ECEF => ENU 和 ENU => ECEF

用户到卫星的观测向量 [ΔxΔyΔz]T[\Delta x \ \Delta y \ \Delta z]^T[ΔxΔyΔz]T 是卫星位置与用户位置的差，其变换到ENU坐标系可通过两次旋转得到。

第一次旋转将ECEF绕Z轴旋转 λ+90°\lambda + 90°λ+90°，第二次选择绕新的X轴旋转 90°−ϕ90°-\phi90°−ϕ。于是得到坐标变换矩阵

S1=[cos⁡(λ+π2)sin⁡(λ+π2)0−sin⁡(λ+π2)cos⁡(λ+π2)0001]=[−sin⁡λcos⁡λ0−cos⁡λ−sin⁡λ0001]\boldsymbol{S_1} = \left[ \begin{array}{ccc} \cos(\lambda+\frac{\pi}{2}) & \sin(\lambda+\frac{\pi}{2}) & 0\\ -\sin(\lambda+\frac{\pi}{2}) & \cos(\lambda+\frac{\pi}{2}) & 0\\ 0 & 0 & 1 \end{array} \right] = \left[ \begin{array}{ccc} -\sin\lambda & \cos\lambda & 0\\ -\cos\lambda & -\sin\lambda & 0\\ 0 & 0 & 1 \end{array} \right] S1​=⎣⎡​cos(λ+2π​)−sin(λ+2π​)0​sin(λ+2π​)cos(λ+2π​)0​001​⎦⎤​=⎣⎡​−sinλ−cosλ0​cosλ−sinλ0​001​⎦⎤​

S2=[1000cos⁡(π2−ϕ)sin⁡(π2−ϕ)0−sin⁡(π2−ϕ)cos⁡(π2−ϕ)]=[1000sin⁡ϕcos⁡ϕ0−cos⁡ϕsin⁡ϕ]\boldsymbol{S_2} = \left[ \begin{array}{ccc} 1 & 0 & 0\\ 0 & \cos(\frac{\pi}{2}-\phi) & \sin(\frac{\pi}{2}-\phi)\\ 0 & -\sin(\frac{\pi}{2}-\phi) & \cos(\frac{\pi}{2}-\phi) \end{array} \right] = \left[ \begin{array}{ccc} 1 & 0 & 0\\ 0 & \sin\phi & \cos\phi\\ 0 & -\cos\phi & \sin\phi \end{array} \right] S2​=⎣⎡​100​0cos(2π​−ϕ)−sin(2π​−ϕ)​0sin(2π​−ϕ)cos(2π​−ϕ)​⎦⎤​=⎣⎡​100​0sinϕ−cosϕ​0cosϕsinϕ​⎦⎤​

S=S2⋅S1=[−sin⁡λcos⁡λ0−sin⁡ϕcos⁡λ−sin⁡ϕsin⁡λcos⁡ϕcos⁡ϕcos⁡λcos⁡ϕsin⁡λsin⁡ϕ]\boldsymbol{S} = \boldsymbol{S_2} \cdot \boldsymbol{S_1} = \left[ \begin{array}{ccc} -\sin\lambda & \cos\lambda & 0\\ -\sin\phi\cos\lambda & -\sin\phi\sin\lambda & \cos\phi\\ \cos\phi\cos\lambda & \cos\phi\sin\lambda & \sin\phi \end{array} \right] S=S2​⋅S1​=⎣⎡​−sinλ−sinϕcosλcosϕcosλ​cosλ−sinϕsinλcosϕsinλ​0cosϕsinϕ​⎦⎤​

而用户到卫星的观测向量在ENU坐标系中表示为：

[ΔeΔnΔu]=S⋅[ΔxΔyΔz]\left[ \begin{array}{ccc} \Delta e\\ \Delta n\\ \Delta u \end{array} \right] = \boldsymbol{S} \cdot \left[ \begin{array}{ccc} \Delta x\\ \Delta y\\ \Delta z \end{array} \right] ⎣⎡​ΔeΔnΔu​⎦⎤​=S⋅⎣⎡​ΔxΔyΔz​⎦⎤​

由上式可得：

[ΔxΔyΔz]=S−1⋅[ΔeΔnΔu]\left[ \begin{array}{ccc} \Delta x\\ \Delta y\\ \Delta z \end{array} \right] = \boldsymbol{S}^{-1} \cdot \left[ \begin{array}{ccc} \Delta e\\ \Delta n\\ \Delta u \end{array} \right] ⎣⎡​ΔxΔyΔz​⎦⎤​=S−1⋅⎣⎡​ΔeΔnΔu​⎦⎤​

有了用户到卫星的观测向量，我们就可以计算该卫星相对用户的方位角和仰角了。仰角 θ\thetaθ是高出水平面的角度：

θ=arcsin⁡(Δu(Δe)2+(Δn)2+(Δu)2)\theta = \arcsin\left(\frac{\Delta u}{\sqrt{(\Delta e)^2+(\Delta n)^2+(\Delta u)^2}}\right) θ=arcsin((Δe)2+(Δn)2+(Δu)2​Δu​)

方位角α\alphaα是观测向量在水平面内的投影指的方位，以北向为0度顺时针为正，即：

α=arctan⁡(ΔeΔn)\alpha = \arctan\left(\frac{\Delta e}{\Delta n}\right) α=arctan(ΔnΔe​)