由于采用的是 Japan rectangular coordinate system，所以在其他区域并不准确，于是借此机会重新推导一下高斯投影正算。
先明确一下相关参数及其关系：
长半径：
a
a
a
短半径：
b
b
b
扁率：
f
=
a
−
b
a
f={a-b \over a}
f=aa−b
第一偏心率：
e
=
a
2
−
b
2
a
e={\sqrt{a^2-b^2} \over a}
e=aa2−b2
第二偏心率：
e
′
=
a
2
−
b
2
b
e'={\sqrt{a^2-b^2} \over b}
e′=ba2−b2
第一辅助函数：
W
=
1
−
e
2
s
i
n
2
B
W=\sqrt{1-e^2sin^2B}
W=1−e2sin2B
第二辅助函数：
V
=
1
+
e
′
2
c
o
s
2
B
V=\sqrt{1+e'^2cos^2B}
V=1+e′2cos2B
卯酉圈曲率半径：
N
=
a
W
=
c
V
N={a \over W} = {c \over V}
N=Wa=Vc
下面结合GPS定位的开源代码：
/*WGS84 Parameters*/
//长半径
AW = 6378137.0;
//扁率
FW = 1.0 / 298.257222101;
//椭球第一偏心率 e
Pe = static_cast<double>(std::sqrt(2.0 * FW - std::pow(FW, 2)));
//第二偏心率 e'
Pet = static_cast<double>(std::sqrt(std::pow(Pe, 2) / (1.0 - std::pow(Pe, 2))));
在这里Pe的形式似乎与上面定义的第一偏心率有所不同，但是经过结合扁率f的化简就可以发现，实质上是一样的，在此就不再详细推导。
接下来是计算子午圈弧长:
X
=
a
(
1
−
e
2
)
[
A
′
B
2
−
B
1
ρ
−
B
′
(
s
i
n
2
B
2
−
s
i
n
2
B
1
)
+
C
′
(
s
i
n
4
B
2
−
s
i
n
4
B
1
)
X=a(1-e^2)[A'{B_2-B_1 \over \rho} - B'(sin2B_2-sin2B_1)+C'(sin4B_2-sin4B_1)
X=a(1−e2)[A′ρB2−B1−B′(sin2B2−sin2B1)+C′(sin4B2−sin4B1)
−
D
′
(
s
i
n
6
B
2
−
s
i
n
6
B
1
)
+
E
′
(
s
i
n
8
B
2
−
s
i
n
8
B
1
)
−
F
′
(
s
i
n
10
B
2
−
s
i
n
B
1
)
+
.
.
.
]
-D'(sin6B_2-sin6B_1)+E'(sin8B_2-sin8B_1)-F'(sin10B_2-sinB_1)+...]
−D′(sin6B2−sin6B1)+E′(sin8B2−sin8B1)−F′(sin10B2−sinB1)+...]
其中
ρ
=
3600
∗
180
/
P
I
\rho = 3600*180/ PI
ρ=3600∗180/PI
A
′
=
1
+
3
4
e
2
+
45
64
e
4
+
175
256
e
6
+
11025
16384
e
8
+
43659
65536
e
10
A'=1+{3\over 4}e^2+{45\over 64}e^4+{175\over 256}e^6+{11025\over 16384}e^8+{43659\over 65536}e^{10}
A′=1+43e2+6445e4+256175e6+1638411025e8+6553643659e10
B
′
=
3
8
e
2
+
15
32
e
4
+
525
1024
e
6
+
2205
4096
e
8
+
72765
131072
e
10
B'={3\over 8}e^2+{15\over 32}e^4+{525\over 1024}e^6+{2205\over 4096}e^8+{72765\over 131072}e^{10}
B′=83e2+3215e4+1024525e6+40962205e8+13107272765e10
C
′
=
15
256
e
4
+
105
1024
e
6
+
2205
16384
e
8
+
10395
65536
e
10
C'={15\over 256}e^4+{105\over 1024}e^6+{2205\over 16384}e^8+{10395\over 65536}e^{10}
C′=25615e4+1024105e6+163842205e8+6553610395e10
D
′
=
35
3072
e
6
+
105
4096
e
8
+
10395
262144
e
10
D'={35\over 3072}e^6+{105\over 4096}e^8+{10395\over 262144}e^{10}
D′=307235e6+4096105e8+26214410395e10
E
′
=
315
131072
e
8
+
3465
524288
e
10
E'={315\over 131072}e^8+{3465\over 524288}e^{10}
E′=131072315e8+5242883465e10
F
′
=
639
1310720
e
10
F'={639\over 1310720}e^{10}
F′=1310720639e10
原理相同，精度可能会有所提升。到此计算的到子午圈弧长，其中日本区域划分为19个，因此可在plane=20的情况下，将附近某个点经纬度设置成原点，即 B 1 B_1 B1，其中 B 2 B_2 B2为接收到的经纬度信息。
其中在高斯投影正算中使用卯酉圈曲率半径，但是由于根据不同经纬度查所代表的实际距离并不相同，因此在我理解下，需要改变这个值来使精度更高，对于纬度来说全球相差不多，但是经度相差较高，具体可以通过网站查找你所在区域每1度所代表的实际距离。
参考文章：Autoware 中 GPS 定位问题 - 简书
Length Of A Degree Of LatitudeAnd Longitude Calculator
通过网上查询资料，得到该区域每1经度所代表的实际距离：95285.86。
由于在源码中，经纬度都换算成了弧度制，在这我们也将距离换成经度的每1弧度所代表的实际距离：
ORGIN_POINT_LON_LENGTH = 95285.86/180*M_PI;
同理纬度也是，然后替换原有代码，直接用经纬度差乘上该地区实际距离：
m_x = static_cast<double>((m_lon - lon_org)*ORGIN_POINT_LON_LENGTH);
m_y = static_cast<double>((m_lat - lat_org)*ORGIN_POINT_LAT_LENGTH);
最后将原点于初始位置经纬度设置大致相同，建图时需要将此时的原点记录下，因为地图建立在map下，导航时需要将此事的GPS POSE转换到和建图时同一坐标系下，也可以添加World坐标系，通过World->Map的TF和World->Base_link的TF将其关联。
扫一扫
2万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
