经纬度坐标与平面坐标转换
问题提出
已知地球上的某点的经纬度信息(这里暂时不考虑高度信息),将这个点沿着某方向移动一段距离之后,其经纬度坐标是多少
最近在设计Gazebo模拟器的GPS传感器插件,遇到了这个问题,乍看一下感觉这个问题似乎很简单,但在做的过程中,发现有一些和以前认知不同的东西,在这里做个记录。
经纬度的定义
在讨论经纬度之前,首先需要对地球构建一个模型。依照我们的常识认为地球是一个椭球体,这个椭球可以看作是一个长轴为a,短轴为b的椭圆形绕着地轴旋转而成。在这个椭球体上我们定义:
- 经线:过旋转轴(地轴) 的平面与椭球面的截线
- 赤道平面:垂直于地轴并通过地心的平面
- 赤道:赤道平面与椭球面相交的交线
- 纬线:过某一点与赤道面平行的平面,与椭球面的交线(截线)
经度的计算:
国际上公认通过英国格林尼治天文台的经线为本初子午线,对于地球上任意一点的子午圈截面与本初子午面之交角称之为经度。由本初子午线起,向东为正,称东经。由0度到+180度。 由本初子午线起,向西为负,称西经。由0度到-180度。
纬度的计算:
地球上任意一点纬线的法线与赤道面的交角。 从赤道起,向北为正,称“北纬”。纬度由0度到+90度; 从赤道起,向南为负,称“南纬”。纬度由0度到-90度。
在上图中,假设Q为本初子午线上的一点, β \beta β为点H的经度, α \alpha α为点H的纬度。
如何计算经纬度变化
上面已经介绍了经纬度的定义,那么下一个问题就是我在地球椭球体上运动了一段距离,那么怎么计算经纬度的变化。
经度变化计算
过地球椭球体表面上的一点做与赤道面平行的平面,根据地球椭球体的定义,这个平面与椭球体的截面为一个标准圆。所以我们可以考虑在一个半径为 R R R标准圆弧上从 C C C点移动一段距离到 D D D点,这段圆弧的长度为 s s s,其对应变化的角度 α \alpha α是多少?
按照圆弧的定义:
α = s R \alpha = \frac{s}{R} α=Rs
纬度变化计算
相比于经度的计算,纬度变化的计算就复杂很多。过地球椭球体表面上的一点且经过地轴的平面与地球椭球体的截面为椭圆形,所以这个问题可以简化成已知椭圆上的两个点 D D D与 E E E的弧长s,求点 D D D与 E E E分别对应的纬度差。
由于
α − β = ( 90 + α ) − ( 90 + β ) \alpha - \beta = (90+\alpha) - (90+\beta) α−β=(90+α)−(90+β)
我们可以发现点 D D D与点 E E E的纬度差可以转换为点 D D D与点 E E E在椭圆上的切线的转角。
记 D E ⌢ \overset{\frown} {DE} DE⌢的长度为 s s s,如果 s s s相比地球半径短很多,我们可以将这段弧看作是标准圆弧。
这个时候问题就转换为求 ∣ θ 1 − θ 2 ∣ |\theta_1 -\theta_2| ∣θ1−θ2∣,经过简单的推导我们就可以得到
∣ θ 1 − θ 2 ∣ = ε = s R |\theta_1 - \theta_2| = \varepsilon = \frac{s}{R} ∣θ1−θ2∣=ε=Rs
曲率半径推导
那么问题就变成了如何得到这个 R R R,这个时候就需要引出一个概念叫做曲率半径
对于曲线上的一点,该点的曲率半径是指最接近该点处曲线的圆弧的半径。曲率半径的值为曲率的倒数。曲率是针对曲线上某个点的切线方向角对弧长的转动率,可以通过微分来定义,表明曲线偏离直线的程度。
这里我们就简单的推导一下曲率半径的计算公式:
假设在平面中有一条曲线
y = f ( x ) y = f(x) y=f(x)
对于曲线上一点 ( x , y ) (x,y) (x,y),其切线的斜率为
y ′ = d f ( x ) d x = tan ( θ ) y' = \frac{df(x)}{dx} = \tan(\theta) y′=dxdf(x)=tan(θ)
可以推出
θ ( x ) = a r c t a n ( y ′ ) \theta(x) = arctan(y') θ(x)=arctan(y′)
θ ( x ) \theta(x) θ(x)也是 x x x的函数,我们继续对 x x x求导
d θ d x = y ′ ′ 1 + y ′ 2 \frac{d\theta}{dx} = \frac{y''}{1+y'^2} dxdθ=1+y′2y′′
从上图可以看出
d s d x = ( d x ) 2 + ( d y ) 2 d x = 1 + ( d y ) 2 ( d x ) 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