目录

经纬度坐标与平面坐标转换问题提出经纬度的定义如何计算经纬度变化经度变化计算纬度变化计算曲率半径推导纬度与直角坐标变换

经纬度坐标与平面坐标转换

问题提出

已知地球上的某点的经纬度信息（这里暂时不考虑高度信息），将这个点沿着某方向移动一段距离之后，其经纬度坐标是多少

最近在设计Gazebo模拟器的GPS传感器插件，遇到了这个问题，乍看一下感觉这个问题似乎很简单，但在做的过程中，发现有一些和以前认知不同的东西，在这里做个记录。

经纬度的定义

在讨论经纬度之前，首先需要对地球构建一个模型。依照我们的常识认为地球是一个椭球体，这个椭球可以看作是一个长轴为a，短轴为b的椭圆形绕着地轴旋转而成。在这个椭球体上我们定义：

经线：过旋转轴(地轴) 的平面与椭球面的截线赤道平面：垂直于地轴并通过地心的平面赤道：赤道平面与椭球面相交的交线纬线：过某一点与赤道面平行的平面,与椭球面的交线（截线）

经度的计算：

国际上公认通过英国格林尼治天文台的经线为本初子午线，对于地球上任意一点的子午圈截面与本初子午面之交角称之为经度。由本初子午线起，向东为正，称东经。由0度到+180度。 由本初子午线起，向西为负，称西经。由0度到-180度。

纬度的计算：

地球上任意一点纬线的法线与赤道面的交角。 从赤道起，向北为正，称“北纬”。纬度由0度到+90度； 从赤道起，向南为负，称“南纬”。纬度由0度到-90度。

在上图中，假设Q为本初子午线上的一点，β\betaβ为点H的经度，α\alphaα为点H的纬度。

如何计算经纬度变化

上面已经介绍了经纬度的定义，那么下一个问题就是我在地球椭球体上运动了一段距离，那么怎么计算经纬度的变化。

经度变化计算

过地球椭球体表面上的一点做与赤道面平行的平面，根据地球椭球体的定义，这个平面与椭球体的截面为一个标准圆。所以我们可以考虑在一个半径为RRR标准圆弧上从CCC点移动一段距离到DDD点，这段圆弧的长度为sss，其对应变化的角度α\alphaα是多少？

按照圆弧的定义:

α=sR\alpha = \frac{s}{R} α=Rs​

纬度变化计算

相比于经度的计算，纬度变化的计算就复杂很多。过地球椭球体表面上的一点且经过地轴的平面与地球椭球体的截面为椭圆形，所以这个问题可以简化成已知椭圆上的两个点DDD与EEE的弧长s，求点DDD与EEE分别对应的纬度差。

由于

α−β=(90+α)−(90+β)\alpha - \beta = (90+\alpha) - (90+\beta) α−β=(90+α)−(90+β)

我们可以发现点DDD与点EEE的纬度差可以转换为点DDD与点EEE在椭圆上的切线的转角。

记DE⌢\overset{\frown} {DE}DE⌢的长度为sss，如果sss相比地球半径短很多，我们可以将这段弧看作是标准圆弧。

这个时候问题就转换为求∣θ1−θ2∣|\theta_1 -\theta_2|∣θ1​−θ2​∣，经过简单的推导我们就可以得到

∣θ1−θ2∣=ε=sR|\theta_1 - \theta_2| = \varepsilon = \frac{s}{R} ∣θ1​−θ2​∣=ε=Rs​

曲率半径推导

那么问题就变成了如何得到这个RRR，这个时候就需要引出一个概念叫做曲率半径

对于曲线上的一点，该点的曲率半径是指最接近该点处曲线的圆弧的半径。曲率半径的值为曲率的倒数。曲率是针对曲线上某个点的切线方向角对弧长的转动率，可以通过微分来定义，表明曲线偏离直线的程度。

这里我们就简单的推导一下曲率半径的计算公式：

假设在平面中有一条曲线

y=f(x)y = f(x) y=f(x)

对于曲线上一点(x,y)(x,y)(x,y)，其切线的斜率为

y′=df(x)dx=tan⁡(θ)y' = \frac{df(x)}{dx} = \tan(\theta) y′=dxdf(x)​=tan(θ)

可以推出

θ(x)=arctan(y′)\theta(x) = arctan(y') θ(x)=arctan(y′)

θ(x)\theta(x)θ(x)也是xxx的函数，我们继续对xxx求导

dθdx=y′′1+y′2\frac{d\theta}{dx} = \frac{y''}{1+y'^2} dxdθ​=1+y′2y′′​

从上图可以看出

dsdx=(dx)2+(dy)2dx=1+(dy)2(dx)2=1+(y′)2\frac{ds}{dx} = \frac{\sqrt{(dx)^2+(dy)^2}}{dx} = \sqrt{1+\frac{(dy)^2}{(dx)^2}} = \sqrt{1+(y')^2} dxds​=dx(dx)2+(dy)2​​=1+(dx)2(dy)2​​=1+(y′)2​

按照曲率半径的微分定义

R=∣ds∣∣dθ∣=∣dsdx∣∣dθdx∣=1+(y′)2∣y′′∣1+y′2=(1+y′2)32∣y′′∣R = \cfrac{|ds|}{|d\theta|} =\cfrac{|\frac{ds}{dx}|}{|\frac{d\theta}{dx}|} = \cfrac{ \sqrt{1+(y')^2}}{ \cfrac{|y''|}{1+y'^2}} = \cfrac{(1+y'^2)^\frac{3}{2}}{|y''|} R=∣dθ∣∣ds∣​=∣dxdθ​∣∣dxds​∣​=1+y′2∣y′′∣​1+(y′)2​​=∣y′′∣(1+y′2)23​​

回到我们的椭圆方程：

x2a2+y2b2=1\frac{x^2}{a^2}+\frac{y^2}{b^2} =1 a2x2​+b2y2​=1

两边对xxx求导可以得到

2xa2+2yy′b2=0⇒y′=−b2xa2y\frac{2x}{a^2} + \frac{2yy'}{b^2} = 0 \Rightarrow y' = -\frac{b^2x}{a^2y} a22x​+b22yy′​=0⇒y′=−a2yb2x​

对xxx继续求导可以得到

y′′=−b2a2⋅y−xy′y2=−b2a2⋅y+x⋅b2xa2yy2=−b2a2⋅a2y2+b2x2a2yy2y'' = -\frac{b^2}{a^2}\cdot\frac{y-xy'}{y^2} = -\frac{b^2}{a^2} \cdot \cfrac{y+x\cdot\cfrac{b^2x}{a^2y}}{y^2} = -\cfrac{b^2}{a^2}\cdot\cfrac{\cfrac{a^2y^2+b^2x^2}{a^2y}}{y^2} y′′=−a2b2​⋅y2y−xy′​=−a2b2​⋅y2y+x⋅a2yb2x​​=−a2b2​⋅y2a2ya2y2+b2x2​​

根据椭圆公式有

b2x2+a2y2=a2b2b^2x^2+a^2y^2 = a^2b^2 b2x2+a2y2=a2b2

带入上式可以有

y′′=−b2a2⋅a2b2a2y3=−b4a2y3y'' = -\frac{b^2}{a^2}\cdot\frac{a^2b^2}{a^2y^3} = -\frac{b^4}{a^2y^3} y′′=−a2b2​⋅a2y3a2b2​=−a2y3b4​

将y′y'y′和y′′y''y′′带入曲率半径的计算公式

R=(1+y′2)32∣y′′∣=(1+b4x2a4y2)32⋅∣a2y3b4∣=(a4y2+b4x2)32∣a6y3∣⋅∣a2y3b4∣R = \cfrac{(1+y'^2)^\frac{3}{2}}{|y''|} = (1+\frac{b^4x^2}{a^4y^2})^\frac{3}{2}\cdot |\frac{a^2y^3}{b^4}| = \frac{(a^4y^2+b^4x^2)^\frac{3}{2}}{|a^6y^3|}\cdot|\frac{a^2y^3}{b^4}| R=∣y′′∣(1+y′2)23​​=(1+a4y2b4x2​)23​⋅∣b4a2y3​∣=∣a6y3∣(a4y2+b4x2)23​​⋅∣b4a2y3​∣

化简一下就能得到我们最终的公式

R=(a4y2+b4x2)32a4b4R= \frac{(a^4y^2+b^4x^2)^\frac{3}{2}}{a^4b^4} R=a4b4(a4y2+b4x2)23​​

纬度与直角坐标变换

回到求纬度变化的公式

∣θ1−θ2∣=ε=sR|\theta_1 - \theta_2| = \varepsilon = \frac{s}{R} ∣θ1​−θ2​∣=ε=Rs​

在上一节我们已经知道在椭圆上的一个点(x0,y0)(x_0,y_0)(x0​,y0​)的曲率半径计算公式，然而往往我们拿到的数据是经度和纬度信息，如何将经纬度信息转换成坐标信息呢？

假设点DDD的坐标为(x0,y0)(x_0,y_0)(x0​,y0​)，其纬度为α\alphaα，那么有

tan⁡(90ο⁡+α)=yD′=−b2x0a2y0=−1tan⁡α\tan{(90^{\operatorname{\omicron}} + \alpha)} = y'_D = -\frac{b^2x_0}{a^2y_0} = -\frac{1}{\tan\alpha} tan(90ο+α)=yD′​=−a2y0​b2x0​​=−tanα1​

可以推出

y0=b2a2x0⋅tan⁡αy_0 = \frac{b^2}{a^2}x_0\cdot\tan{\alpha} y0​=a2b2​x0​⋅tanα

又有

x02a2+y02b2=1\frac{x_0^2}{a^2} + \frac{y_0^2}{b^2} = 1 a2x02​​+b2y02​​=1

由上面两个式子可以联合求解得到

x0=acos⁡α1−e2sin⁡2α,y0=a(1−e2)sin⁡α1−e2sin⁡2αx_0 = \frac{a\cos{\alpha}}{\sqrt{1-e^2\sin^2\alpha}}, y_0 = \frac{a(1-e^2)\sin{\alpha}}{\sqrt{1-e^2\sin^2\alpha}} x0​=1−e2sin2α​acosα​,y0​=1−e2sin2α​a(1−e2)sinα​

其中eee表示椭圆的偏心率，计算公式为

e=a2−b2ae = \frac{\sqrt{a^2-b^2}}{a} e=aa2−b2​​

将x0x_0x0​，y0y_0y0​带入椭圆的曲率半径公式，就可以得到以纬度表示的曲率半径，这里就不做计算了。