高斯投影正反算表达式及代码
一、带号及分带中央经线计算
提示：这里分成3度带和6度带的情况：
1.1 3度带带号计算
z
o
n
e
=
f
l
o
o
r
(
(
L
+
1.5
)
/
3
)
zone = floor((L + 1.5) / 3)
zone=floor((L+1.5)/3)
代码：
int zone = (int)floor((lbh.Longitude() + 1.5) / 3);
1.2 3度带中央经线计算
c
e
n
t
r
a
l
_
m
e
r
i
d
i
a
n
=
z
o
n
e
∗
3
central\_meridian = zone * 3
central_meridian=zone∗3
代码：
int zone = (int)floor((lbh.Longitude() + 1.5) / 3);
double central_meridian = zone * 3.0;
1.3 6度带带号计算
z
o
n
e
=
f
l
o
o
r
(
L
/
6
)
+
1
zone = floor(L / 6) + 1
zone=floor(L/6)+1
代码：
int zone = (int)floor(lbh.Longitude().GetRadixDegree() / 6) + 1;
1.4 6度带中央经线计算
c
e
n
t
r
a
l
_
m
e
r
i
d
i
a
n
=
z
o
n
e
∗
6
−
3
central\_meridian = zone * 6 - 3
central_meridian=zone∗6−3
代码：
int zone = (int)floor(lbh.Longitude().GetRadixDegree() / 6) + 1;
double central_meridian = zone * 6.0 - 3;
二、高斯投影正算表达式及代码
提示：参考自《大地测量学基础》第二版：p.169高斯投影坐标正算公式 p.115-116子午线弧长计算公式 p.110-111子午圈曲率半径计算公式，《CH∕T 2014-2016 大地测量控制点坐标转换技术规范》中高斯正算表达式存在异议的地方，未采用。本文采用两者结合的方式。
2.1 高斯投影正算表达式
x = X + N t cos 2 B ℓ 2 ρ 2 [ 0.5 + 1 24 ( 5 − t 2 + 9 η 2 + 4 η 2 ) cos 2 B ℓ 2 ρ 2 + 1 720 ( 61 − 58 t 2 + t 4 ) cos 4 B ℓ 4 ρ 4 ] y = N cos B ℓ ρ [ 1 + 1 6 ( 1 − t 2 + η 2 ) cos 2 B ℓ 2 ρ 2 + 1 120 ( 5 − 18 t 2 + t 4 + 14 η 2 − 58 η 2 t 2 ) cos 4 B ℓ 4 ρ 4 ] \begin{equation} \begin{split} x =& X+Nt\cos^{2}B\frac{{\ell}^2}{{\rho}^2}[0.5 +\frac{1}{24}(5-t^2+9{\eta}^2+4{\eta}^2)\cos^2B\frac{\ell^2}{\rho^2}+ \\ & \frac{1}{720}(61-58t^2+t^4)\cos^4B\frac{\ell^4}{\rho^4}] \\ y =& N{\cos}B\frac{\ell}{\rho}[1+\frac{1}{6}(1-t^2+\eta^2)\cos^2B\frac{\ell^2}{\rho^2} + \\ & \frac{1}{120}(5-18t^2+t^4+14\eta^2-58\eta^2t^2)\cos^4B\frac{\ell^4}{\rho^4}] \end{split} \end{equation} x=y=X+Ntcos2Bρ2ℓ2[0.5+241(5−t2+9η2+4η2)cos2Bρ2ℓ2+7201(61−58t2+t4)cos4Bρ4ℓ4]NcosBρℓ[1+61(1−t2+η2)cos2Bρ2ℓ2+1201(5−18t2+t4+14η2−58η2t2)cos4Bρ4ℓ4]
式中：
式中：
式中：
L
,
B
,
H
——
点位经度、纬度、大地高，经纬度单位为弧度
(
r
a
d
)
，
大地高单位为米
(
m
)
；
ρ
——
角度与弧度间的转换量，单位为角秒
(
′
′
)
，
ρ
=
3600
×
180
π
；
ℓ
——
经度差，单位为角秒
(
′
′
)
，
ℓ
=
L
输入位置
−
L
中央经线
；
e
2
——
第一偏心率的平方，
f
为椭球扁率，
a
、
b
为椭球长、短半轴，
e
2
=
2
f
−
f
2
=
1
−
b
2
a
2
=
e
′
2
1
+
e
′
2
，
f
=
a
−
b
a
，
b
=
a
1
−
e
2
，
c
=
a
2
b
；
e
′
——
第二偏心率，
e
′
=
a
2
−
b
2
b
=
e
2
1
−
e
2
，
η
2
=
e
′
2
cos
2
B
，
t
=
tan
B
；
N
——
卯酉圈曲率半径，
N
=
a
1
−
e
2
sin
2
B
=
c
1
+
e
′
2
cos
2
B
；
X
——
子午线弧长，设有子午线上两点
p
1
和
p
2
，
p
1
在赤道上，
p
2
纬度为
B
′
，
B
′
单位为弧度
(
r
a
d
)
，
p
1
、
p
2
间的子午线弧长
X
的计算公式为：
X
=
a
0
B
′
−
a
2
2
sin
2
B
′
+
a
4
4
sin
4
B
′
−
a
6
6
sin
6
B
′
+
a
8
8
sin
8
B
′
，
或
X
=
a
0
B
′
−
sin
B
′
cos
B
′
[
(
a
2
−
a
4
+
a
6
)
+
(
2
a
4
−
16
3
a
6
)
sin
2
B
′
+
16
3
a
6
sin
4
B
′
]
，式子中：
{
m
0
=
a
(
1
−
e
2
)
m
2
=
3
2
e
2
m
0
m
4
=
5
4
e
2
m
2
m
6
=
7
6
e
2
m
4
m
8
=
9
8
e
2
m
6
{
a
0
=
m
0
+
1
2
m
2
+
3
8
m
4
+
5
16
m
6
+
35
128
m
8
+
⋯
a
2
=
1
2
m
2
+
1
2
m
4
+
15
32
m
6
+
7
16
m
8
a
4
=
1
8
m
4
+
3
16
m
6
+
7
32
m
8
a
6
=
1
32
m
6
+
1
16
m
8
a
8
=
1
128
m
8
\begin{align*} L,B,H \quad \text{\textemdash}\text{\textemdash} & 点位经度、纬度、大地高，经纬度单位为弧度(rad)，\\ &大地高单位为米(m)； \\ \rho \quad \text{\textemdash}\text{\textemdash} & 角度与弧度间的转换量，单位为角秒('')，\rho=3600\times\frac{180}{\pi}； \\ \ell \quad \text{\textemdash}\text{\textemdash} & 经度差，单位为角秒('')，\ell=L_{输入位置}-L_{中央经线}； \\ e^2 \quad \text{\textemdash}\text{\textemdash} & 第一偏心率的平方，f为椭球扁率，a、b为椭球长、短半轴，\\ & e^2=2f-f^2=1-\frac{b^2}{a^2}=\frac{{e^{\rq}}^2}{1+{e^{\rq}}^2}，f=\frac{a-b}{a}，b=a\sqrt{1-e^2}，c=\frac{a^2}{b}；\\ e^{\rq} \quad \text{\textemdash}\text{\textemdash} & 第二偏心率，e^{\rq}=\frac{\sqrt{a^2-b^2}}{b}=\frac{e^2}{1-e^2}，\eta^2={e^{\rq}}^2\cos^2B，t=\tan{B}； \\ N \quad \text{\textemdash}\text{\textemdash} & 卯酉圈曲率半径，N=\frac{a}{\sqrt{1-e^2\sin^2B}}=\frac{c}{\sqrt{1+{e^{\rq}}^2\cos^2B}}；\\ X \quad \text{\textemdash}\text{\textemdash} & 子午线弧长，设有子午线上两点p_1和p_2，p_1在赤道上，p_2纬度为B^{\rq}， \\ & B^{\rq}单位为弧度(rad)，p_1、p_2间的子午线弧长X的计算公式为：\\ & X = a_0B^{\rq}-\frac{a_2}{2}\sin{2B^{\rq}}+\frac{a_4}{4}\sin{4B^{\rq}}-\frac{a_6}{6}\sin{6B^{\rq}}+\frac{a_8}{8}\sin{8B^{\rq}}，\\ & 或X=a_0B^{\rq}-\sin{B^{\rq}}\cos{B^{\rq}}[(a_2-a_4+a_6)+(2a_4-\frac{16}{3}a_6)\sin^2B^{\rq}+ \\ & \frac{16}{3}a_6\sin^4B^{\rq}]，式子中：\\ & \begin{dcases} m_0 =& a(1-e^2) \\ m_2 =& \frac{3}{2}e^2m_0 \\ m_4 =& \frac{5}{4}e^2m_2 \\ m_6 =& \frac{7}{6}e^2m_4 \\ m_8 =& \frac{9}{8}e^2m_6 \end{dcases} \\ & \begin{dcases} a_0 =& m_0 + \frac{1}{2}m_2+\frac{3}{8}m_4+\frac{5}{16}m_6+\frac{35}{128}m_8+\cdots \\ a_2 =& \frac{1}{2}m_2+\frac{1}{2}m_4+\frac{15}{32}m_6+\frac{7}{16}m_8 \\ a_4 =& \frac{1}{8}m_4+\frac{3}{16}m_6+\frac{7}{32}m_8 \\ a_6 =& \frac{1}{32}m_6+\frac{1}{16}m_8 \\ a_8 =& \frac{1}{128}m_8 \end{dcases} \end{align*}
L,B,H——ρ——ℓ——e2——e′——N——X——点位经度、纬度、大地高，经纬度单位为弧度(rad)，大地高单位为米(m)；角度与弧度间的转换量，单位为角秒(′′)，ρ=3600×π180；经度差，单位为角秒(′′)，ℓ=L输入位置−L中央经线；第一偏心率的平方，f为椭球扁率，a、b为椭球长、短半轴，e2=2f−f2=1−a2b2=1+e′2e′2，f=aa−b，b=a1−e2，c=ba2；第二偏心率，e′=ba2−b2=1−e2e2，η2=e′2cos2B，t=tanB；卯酉圈曲率半径，N=1−e2sin2Ba=1+e′2cos2Bc；子午线弧长，设有子午线上两点p1和p2，p1在赤道上，p2纬度为B′，B′单位为弧度(rad)，p1、p2间的子午线弧长X的计算公式为：X=a0B′−2a2sin2B′+4a4sin4B′−6a6sin6B′+8a8sin8B′，或X=a0B′−sinB′cosB′[(a2−a4+a6)+(2a4−316a6)sin2B′+316a6sin4B′]，式子中：⎩
⎨
⎧m0=m2=m4=m6=m8=a(1−e2)23e2m045e2m267e2m489e2m6⎩
⎨
⎧a0=a2=a4=a6=a8=m0+21m2+83m4+165m6+12835m8+⋯21m2+21m4+3215m6+167m881m4+163m6+327m8321m6+161m81281m8
katex数学公式太长如何分成多行：
使用split命令。split命令可以将一个长公式分割成多行，并在适当的位置自动添加换行符：\\
对齐需要使用对齐符号：&
示例代码如下：
\begin{equation}
\begin{split}
a &= b + c \\
&= d + e + f \\
&= g + h + i + j
\end{split}
\end{equation}
2.2《CH∕T 2014-2016 大地测量控制点坐标转换技术规范》异议点
《CH∕T 2014-2016 大地测量控制点坐标转换技术规范》中给出的高斯正算表达式有几点异议需要注意：
- (1) 表达式计算y分量多了一个N：
- (2) 计算子午线弧长关系式中，arcB表示B的弧度值：
2.3 代码实现
C++实现代码如下：
void ForwardWith3Degree(std::vector<mono::Coordinate3d>& targetProjs, const std::vector<mono::SLonLatHeight3d>& originLBHs, const CSpheroid* spheriod)
{
size_t pnum = originLBHs.size();
if (pnum == 0)
return;
//const double q = 3600.0 * (180.0 / mono::M_PI); // 1弧度 = q角秒
//const double q_pow2 = q * q;
const double q = 180.0 / mono::M_PI; // 1弧度 = q角度
const double f = spheriod->GetFlatting(); // 椭球扁率
const double a = spheriod->GetSemimajorAxis(); // 椭球长半轴
const double b = spheriod->GetSemiminorAxis(); // 椭球短半轴
const double e_pow2 = 2 * f - f * f; // (a * a - b * b) / a * a; 2 * f - f * f; // 第一篇心率的平方 e^2
const double e_sec_pow2 = e_pow2 / (1 - e_pow2); // e^2 / (1 - e^2); (a / b) * (a / b) - 1; (a * a - b * b) / b * b 此表达式存在大数值的平方运算 导致计算结果异常; // 第二偏心率的平方 e'^2
#ifdef GAUSS_KRUGER_CHT_2014_2016
const double e_pow2_pow2 = e_pow2 * e_pow2; // e^2 * e^2
const double e_pow2_pow3 = e_pow2_pow2 * e_pow2; // e^2 * e^2 * e^2
const double e_pow2_pow4 = e_pow2_pow3 * e_pow2; // e^2 * e^2 * e^2 * e^2
const double e_pow2_pow5 = e_pow2_pow4 * e_pow2; // e^2 * e^2 * e^2 * e^2 * e^2
const double e_pow2_pow6 = e_pow2_pow5 * e_pow2; // e^2 * e^2 * e^2 * e^2 * e^2 * e^2
// 幂级数展开中间变量
const double A_exp = 1 + (3.0 / 4) * e_pow2 + (45.0 / 64) * e_pow2_pow2 +
(175.0 / 256) * e_pow2_pow3 + (11025.0 / 16384) * e_pow2_pow4 +
(43659.0 / 65536) * e_pow2_pow5 +
(693693.0 / 1048576) * e_pow2_pow6;
const double B_exp = (3.0 / 8) * e_pow2 + (15.0 / 32) * e_pow2_pow2 +
(525.0 / 1024) * e_pow2_pow3 + (2205.0 / 4096) * e_pow2_pow4 +
(72765.0 / 131072) * e_pow2_pow5 +
(297297.0 / 524288) * e_pow2_pow6;
const double C_exp = (15.0 / 256) * e_pow2_pow2 + (105.0 / 1024) * e_pow2_pow3 +
(2205.0 / 16384) * e_pow2_pow4 +
(10395.0 / 65536) * e_pow2_pow5 +
(1486485.0 / 8388608) * e_pow2_pow6;
const double D_exp = (35.0 / 3072) * e_pow2_pow3 +
(105.0 / 4096) * e_pow2_pow4 +
(10395.0 / 262144) * e_pow2_pow5 +
(55055.0 / 1048576) * e_pow2_pow6;
const double E_exp = (315.0 / 131072) * e_pow2_pow4 +
(3465.0 / 524288) * e_pow2_pow5 +
(99099.0 / 8388608) * e_pow2_pow6;
const double F_exp = (693.0 / 1310720) * e_pow2_pow5 +
(9009.0 / 5242880) * e_pow2_pow6;
const double G_exp = (1001.0 / 8388608) * e_pow2_pow6;
// 子午线弧长
const double X_exp = a * (1 - e_pow2);
#else
double m0 = a * (1 - e_pow2);
double m2 = (3.0 / 2.0) * e_pow2 * m0;
double m4 = (5.0 / 4.0) * e_pow2 * m2;
double m6 = (7.0 / 6.0) * e_pow2 * m4;
double m8 = (9.0 / 8.0) * e_pow2 * m6;
double a0 = m0 + (1.0 / 2) * m2 + (3.0 / 8) * m4 + (5.0 / 16) * m6 + (35.0 / 128) * m8;
double a2 = (1.0 / 2) * m2 + (1.0 / 2) * m4 + (15.0 / 32) * m6 + (7.0 / 16) * m8;
double a4 = (1.0 / 8) * m4 + (3.0 / 16) * m6 + (7.0 / 32) * m8;
double a6 = (1.0 / 32) * m6 + (1.0 / 16) * m8;
double a8 = (1.0 / 128) * m8;
#endif // GAUSS_KRUGER_CHT_2014_2016
double central_meridian = 0;
double l_delta = 0;
double l_delta_pow2 = 0;
double B = 0;
double X = 0;
double N = 0;
double t = 0;
double t_pow2 = 0;
//double l_q_comm = 0;
//double l_q_comm_pow2 = 0;
double cos_B_pow2 = 0;
double n_pow2 = 0;
double x = 0;
double y = 0;
std::vector<mono::Coordinate3d> projs;
projs.resize(pnum);
// 遍历坐标点
for (size_t pi = 0; pi < pnum; pi++) {
const auto& lbh = originLBHs[pi];
// 计算可能所处的投影带中央子午线
central_meridian = b3Degree ? GetCentralMeridianWith3Degree(lbh) : GetCentralMeridianWith6Degree(lbh);
//l_delta = TORADAIN(lbh.Longitude().GetRadixDegree()) - TORADAIN(central_meridian);
//l_delta = (lbh.Longitude().GetRadixDegree() - central_meridian) * 3600; // 角秒
l_delta = (lbh.Longitude().GetRadixDegree() - central_meridian) / q; // 角度 --> 弧度
B = TORADAIN(lbh.Latitude().GetRadixDegree());
// 子午线弧长
#ifdef GAUSS_KRUGER_CHT_2014_2016
// 《CH∕T 2014-2016 大地测量控制点坐标转换技术规范》算法
X = X_exp * (A_exp * B - B_exp * sin(2 * B) + C_exp * sin(4 * B) - D_exp * sin(6 * B) + E_exp * sin(8 * B) - F_exp * sin(10 * B) + G_exp * sin(12 * B));
#else
// 《大地测量学基础》p.169高斯投影坐标正算公式 p.115-116子午线弧长计算公式 p.110-111子午圈曲率半径计算公式
// https://zhuanlan.zhihu.com/p/649655085
// https://blog.csdn.net/weixin_47156401/article/details/124363834
// url所给的算法
X = a0 * B - a2 / 2.0 * sin(2 * B) + a4 / 4 * sin(4 * B) - a6 / 6 * sin(6 * B) + a8 / 8 * sin(8 * B);
#endif
l_delta_pow2 = l_delta * l_delta;
N = a / sqrt(1 - e_pow2 * sin(B) * sin(B));
t = tan(B);
t_pow2 = t * t;
//l_q_comm = l_delta / q;
//l_q_comm_pow2 = l_delta * l_delta / q_pow2;
cos_B_pow2 = cos(B) * cos(B);
n_pow2 = e_sec_pow2 * cos_B_pow2;
// x = X + N * t * cos_B_pow2 * l_q_comm_pow2 * (
// 0.5 +
// 1.0 / 24 * (5.0 - t_pow2 + 9 * n_pow2 + 4 * n_pow2 * n_pow2) * cos_B_pow2 * l_q_comm_pow2 +
// 1.0 / 720 * (61.0 - 58 * t_pow2 + t_pow2 * t_pow2) * cos_B_pow2 * cos_B_pow2 * l_q_comm_pow2 * l_q_comm_pow2);
// y = N * cos(B) * l_q_comm * (
// 1 +
// 1.0 / 6 * (1 - t_pow2 + n_pow2) * cos_B_pow2 * l_q_comm_pow2 +
// 1.0 / 120 * (5 - 18 * t_pow2 + t_pow2 * t_pow2 + 14 * n_pow2 - 58 * n_pow2 * t_pow2) * cos_B_pow2 * cos_B_pow2 * l_q_comm_pow2 * l_q_comm_pow2);
x = X + N * t * cos_B_pow2 * l_delta_pow2 * (
0.5 +
1.0 / 24 * (5.0 - t_pow2 + 9 * n_pow2 + 4 * n_pow2 * n_pow2) * cos_B_pow2 * l_delta_pow2 +
1.0 / 720 * (61.0 - 58 * t_pow2 + t_pow2 * t_pow2) * cos_B_pow2 * cos_B_pow2 * l_delta_pow2 * l_delta_pow2);
y = N * cos(B) * l_delta * (
1 +
1.0 / 6 * (1 - t_pow2 + n_pow2) * cos_B_pow2 * l_delta_pow2 +
1.0 / 120 * (5 - 18 * t_pow2 + t_pow2 * t_pow2 + 14 * n_pow2 - 58 * n_pow2 * t_pow2) * cos_B_pow2 * cos_B_pow2 * l_delta_pow2 * l_delta_pow2);
projs[pi][0] = y + 500000; // 向东偏移500km
projs[pi][1] = x;
projs[pi][2] = lbh.Height();
}
targetProjs = std::move(projs);
}
三、高斯投影反算表达式及代码
提示：参考自《大地测量学基础》第二版：p.171高斯投影坐标反算公式 p.115-116子午线弧长计算公式 p.110-111子午圈曲率半径计算公式，未采用《CH∕T 2014-2016 大地测量控制点坐标转换技术规范》。本文采用两者结合的方式。
2.1 高斯投影反算表达式
B = B f − t f 2 M f y ( y N f ) [ 1 − 1 12 ( 5 + 3 t f 2 + η f 2 − 9 η f 2 t f 2 ) ( y N f ) 2 + 1 360 ( 61 + 90 t f 2 + 45 t f 4 ) ( y N f ) 4 ] ℓ = 1 cos B f ( y N f ) [ 1 − 1 6 ( 1 + 2 t f 2 + η f 2 ) ( y N f ) 2 + 1 120 ( 5 + 28 t f 2 + 24 t f 4 + 6 η f 2 + 8 η f 2 t f 2 ) ( y N f ) 4 ] y = e n h . e ( ) − 500000 B = B ∗ 180 π L = L 0 + ℓ ∗ 180 π \begin{equation} \begin{split} B =& B_f - \frac{t_f }{2M_f}y(\frac{y}{N_f})[1 - \frac{1}{12}(5+3{t_f}^2+{\eta_f}^2-9{\eta_f}^2{t_f}^2)(\frac{y}{N_f})^2+ \\ & \frac{1}{360}(61+90{t_f}^2+45{t_f}^4)(\frac{y}{N_f})^4] \\ \ell = & \frac{1}{\cos{B_f}}(\frac{y}{N_f})[1-\frac{1}{6}(1+2{t_f}^2+{\eta_f}^2)(\frac{y}{N_f})^2+ \\ & \frac{1}{120}(5+28{t_f}^2+24{t_f}^4+6{\eta_f}^2+8{\eta_f}^2{t_f}^2)(\frac{y}{N_f})^4] \\ y =& enh.e() - 500000 \\ B =& B * \frac{180}{\pi} \\ L =& L_0 + \ell*\frac{180}{\pi} \end{split} \end{equation} B=ℓ=y=B=L=Bf−2Mftfy(Nfy)[1−121(5+3tf2+ηf2−9ηf2tf2)(Nfy)2+3601(61+90tf2+45tf4)(Nfy)4]cosBf1(Nfy)[1−61(1+2tf2+ηf2)(Nfy)2+1201(5+28tf2+24tf4+6ηf2+8ηf2tf2)(Nfy)4]enh.e()−500000B∗π180L0+ℓ∗π180
式中：
式中：
式中：
e
n
h
——
点位投影平面坐标，
e
n
h
.
e
(
)
为东向坐标，
e
n
h
.
n
(
)
为北向坐标，
e
n
h
.
h
(
)
为大地高单位为米
(
m
)
；
L
,
B
,
H
——
点位经度、纬度、大地高，经纬度单位为弧度
(
r
a
d
)
，
大地高单位为米
(
m
)
；
L
0
——
高斯投影分带的中央经线，单位为弧度
(
r
a
d
)
；
ℓ
——
经度差，单位为弧度
(
r
a
d
)
；
B
f
——
利用子午线弧长反算得到的大地纬度，单位为弧度
(
r
a
d
)
，
B
f
i
为基于迭代法利用子午线弧长反算得到的大地纬度，
迭代初始值：
B
f
1
=
X
a
0
，
迭代值：
B
f
i
+
1
=
X
−
F
(
B
f
i
)
a
0
；
e
2
——
第一偏心率的平方，
f
为椭球扁率，
a
、
b
为椭球长、短半轴，
e
2
=
2
f
−
f
2
=
1
−
b
2
a
2
=
e
′
2
1
+
e
′
2
，
f
=
a
−
b
a
，
b
=
a
1
−
e
2
，
c
=
a
2
b
；
e
′
——
第二偏心率，
e
′
=
a
2
−
b
2
b
=
e
2
1
−
e
2
，
η
f
2
=
e
′
2
cos
2
B
f
i
，
t
f
=
tan
B
f
i
；
M
——
子午圈曲率半径，
M
=
a
(
1
−
e
2
)
(
1
−
e
2
sin
2
B
f
i
)
3
=
c
(
1
+
e
′
2
cos
2
B
f
i
)
3
；
N
——
卯酉圈曲率半径，
N
=
a
1
−
e
2
sin
2
B
f
i
=
c
1
+
e
′
2
cos
2
B
f
i
；
X
——
子午线弧长，设有子午线上两点
p
1
和
p
2
，
p
1
在赤道上，
p
2
纬度为
B
f
i
，
B
f
i
单位为弧度
(
r
a
d
)
，
p
1
、
p
2
间的子午线弧长
X
的计算公式为：
X
=
a
0
B
f
i
−
a
2
2
sin
2
B
f
i
+
a
4
4
sin
4
B
f
i
−
a
6
6
sin
6
B
f
i
+
a
8
8
sin
8
B
f
i
，
或
X
=
a
0
B
f
i
−
sin
B
f
i
cos
B
f
i
[
(
a
2
−
a
4
+
a
6
)
+
(
2
a
4
−
16
3
a
6
)
sin
2
B
f
i
+
16
3
a
6
sin
4
B
f
i
]
，
F
(
B
f
i
)
=
−
a
2
2
sin
2
B
f
i
+
a
4
4
sin
4
B
f
i
−
a
6
6
sin
6
B
f
i
+
a
8
8
sin
8
B
f
i
式子中：
{
m
0
=
a
(
1
−
e
2
)
m
2
=
3
2
e
2
m
0
m
4
=
5
4
e
2
m
2
m
6
=
7
6
e
2
m
4
m
8
=
9
8
e
2
m
6
{
a
0
=
m
0
+
1
2
m
2
+
3
8
m
4
+
5
16
m
6
+
35
128
m
8
+
⋯
a
2
=
1
2
m
2
+
1
2
m
4
+
15
32
m
6
+
7
16
m
8
a
4
=
1
8
m
4
+
3
16
m
6
+
7
32
m
8
a
6
=
1
32
m
6
+
1
16
m
8
a
8
=
1
128
m
8
\begin{align*} enh \quad \text{\textemdash}\text{\textemdash} & 点位投影平面坐标，enh.e()为东向坐标，enh.n()为北向坐标， \\ & enh.h()为大地高单位为米(m)； \\ L,B,H \quad \text{\textemdash}\text{\textemdash} & 点位经度、纬度、大地高，经纬度单位为弧度(rad)，\\ &大地高单位为米(m)； \\ L_0 \quad \text{\textemdash}\text{\textemdash} & 高斯投影分带的中央经线，单位为弧度(rad)； \\ \ell \quad \text{\textemdash}\text{\textemdash} & 经度差，单位为弧度(rad)； \\ B_f \quad \text{\textemdash}\text{\textemdash} & 利用子午线弧长反算得到的大地纬度，单位为弧度(rad)，\\ & B^i_f为基于迭代法利用子午线弧长反算得到的大地纬度， \\ & 迭代初始值：B^1_f = \frac{X}{a_0}， \\ & 迭代值：B^{i+1}_f = \frac{X-F(B^i_f)}{a_0}； \\ e^2 \quad \text{\textemdash}\text{\textemdash} & 第一偏心率的平方，f为椭球扁率，a、b为椭球长、短半轴，\\ & e^2=2f-f^2=1-\frac{b^2}{a^2}=\frac{{e^{\rq}}^2}{1+{e^{\rq}}^2}，f=\frac{a-b}{a}，b=a\sqrt{1-e^2}，c=\frac{a^2}{b}；\\ e^{\rq} \quad \text{\textemdash}\text{\textemdash} & 第二偏心率，e^{\rq}=\frac{\sqrt{a^2-b^2}}{b}=\frac{e^2}{1-e^2}，{\eta_f}^2={e^{\rq}}^2\cos^2B^i_f，t_f=\tan{B^i_f}； \\ M \quad \text{\textemdash}\text{\textemdash} & 子午圈曲率半径，M=\frac{a(1-e^2)}{(\sqrt{1-e^2\sin^2B^i_f})^3}=\frac{c}{(\sqrt{1+{e^{\rq}}^2\cos^2B^i_f})^3}；\\ N \quad \text{\textemdash}\text{\textemdash} & 卯酉圈曲率半径，N=\frac{a}{\sqrt{1-e^2\sin^2B^i_f}}=\frac{c}{\sqrt{1+{e^{\rq}}^2\cos^2B^i_f}}；\\ X \quad \text{\textemdash}\text{\textemdash} & 子午线弧长，设有子午线上两点p_1和p_2，p_1在赤道上，p_2纬度为B^i_f， \\ & B^i_f单位为弧度(rad)，p_1、p_2间的子午线弧长X的计算公式为：\\ & X = a_0B^i_f-\frac{a_2}{2}\sin{2B^i_f}+\frac{a_4}{4}\sin{4B^i_f}-\frac{a_6}{6}\sin{6B^i_f}+\frac{a_8}{8}\sin{8B^i_f}，\\ & 或X=a_0B^i_f-\sin{B^i_f}\cos{B^i_f}[(a_2-a_4+a_6)+(2a_4-\frac{16}{3}a_6)\sin^2B^i_f+ \\ & \frac{16}{3}a_6\sin^4B^i_f]，\\ & F(B^i_f) = -\frac{a_2}{2}\sin{2B^i_f}+\frac{a_4}{4}\sin{4B^i_f}-\frac{a_6}{6}\sin{6B^i_f}+\frac{a_8}{8}\sin{8B^i_f} \\ & 式子中：\\ & \begin{dcases} m_0 =& a(1-e^2) \\ m_2 =& \frac{3}{2}e^2m_0 \\ m_4 =& \frac{5}{4}e^2m_2 \\ m_6 =& \frac{7}{6}e^2m_4 \\ m_8 =& \frac{9}{8}e^2m_6 \end{dcases} \\ & \begin{dcases} a_0 =& m_0 + \frac{1}{2}m_2+\frac{3}{8}m_4+\frac{5}{16}m_6+\frac{35}{128}m_8+\cdots \\ a_2 =& \frac{1}{2}m_2+\frac{1}{2}m_4+\frac{15}{32}m_6+\frac{7}{16}m_8 \\ a_4 =& \frac{1}{8}m_4+\frac{3}{16}m_6+\frac{7}{32}m_8 \\ a_6 =& \frac{1}{32}m_6+\frac{1}{16}m_8 \\ a_8 =& \frac{1}{128}m_8 \end{dcases} \end{align*}
enh——L,B,H——L0——ℓ——Bf——e2——e′——M——N——X——点位投影平面坐标，enh.e()为东向坐标，enh.n()为北向坐标，enh.h()为大地高单位为米(m)；点位经度、纬度、大地高，经纬度单位为弧度(rad)，大地高单位为米(m)；高斯投影分带的中央经线，单位为弧度(rad)；经度差，单位为弧度(rad)；利用子午线弧长反算得到的大地纬度，单位为弧度(rad)，Bfi为基于迭代法利用子午线弧长反算得到的大地纬度，迭代初始值：Bf1=a0X，迭代值：Bfi+1=a0X−F(Bfi)；第一偏心率的平方，f为椭球扁率，a、b为椭球长、短半轴，e2=2f−f2=1−a2b2=1+e′2e′2，f=aa−b，b=a1−e2，c=ba2；第二偏心率，e′=ba2−b2=1−e2e2，ηf2=e′2cos2Bfi，tf=tanBfi；子午圈曲率半径，M=(1−e2sin2Bfi)3a(1−e2)=(1+e′2cos2Bfi)3c；卯酉圈曲率半径，N=1−e2sin2Bfia=1+e′2cos2Bfic；子午线弧长，设有子午线上两点p1和p2，p1在赤道上，p2纬度为Bfi，Bfi单位为弧度(rad)，p1、p2间的子午线弧长X的计算公式为：X=a0Bfi−2a2sin2Bfi+4a4sin4Bfi−6a6sin6Bfi+8a8sin8Bfi，或X=a0Bfi−sinBficosBfi[(a2−a4+a6)+(2a4−316a6)sin2Bfi+316a6sin4Bfi]，F(Bfi)=−2a2sin2Bfi+4a4sin4Bfi−6a6sin6Bfi+8a8sin8Bfi式子中：⎩
⎨
⎧m0=m2=m4=m6=m8=a(1−e2)23e2m045e2m267e2m489e2m6⎩
⎨
⎧a0=a2=a4=a6=a8=m0+21m2+83m4+165m6+12835m8+⋯21m2+21m4+3215m6+167m881m4+163m6+327m8321m6+161m81281m8
2.2 子午线弧长求大地纬度迭代计算式精度测算
- 1角秒 = 1/3600角度 = 1/3600 * pi/180弧度 = pi/648000弧度 = pi/648000 * 3600弧秒
- 1纬度大约等于111km ==> 1角秒大约等于30.833333米
- 为了确保厘米级精度，按照1角秒 / 100000000，此时对应距离大约等于30.833333 / 100000000 = 0.00000030833333米
- B2-B1 <= 0.00000030833333 表示B2 == B1
- 则精度计算阈值 = pi/648000 * 3600 * 1/100000000
2.3 代码实现
/**
* @brief 给定一个 > 0.0 and < 1.0 的小数 求出将输入的小数变为一位整数位的小数 需要移动的小数位数
*
* @param v > 0.0 and < 1.0 的小数
*
* @return 将输入的小数变为一位整数位的小数 需要移动的小数位数
*/
static int GetNearIntegerOffset(double v) {
int offset = 0;
int integer = 0;
double cv = 0;
double lv = 0;
lv = v;
while (true)
{
++offset;
cv = lv * 10;
integer = (int)floor(cv);
if (integer == 0) {
lv = cv;
continue;
}
break;
}
return offset;
}
void Inverse(std::vector<mono::SLonLatHeight3d>& targetLBHs, const std::vector<mono::Coordinate3d>& originProjs, const CSpheroid* spheriod, const CAngularAngle& centralMeridian)
{
size_t pnum = originProjs.size();
if (pnum == 0)
return;
// 《大地测量学基础》p.171高斯投影坐标反算公式 p.118-119由子午线弧长求大地纬度公式 p.115-116子午线弧长计算公式 p.110-111子午圈曲率半径计算公式
const double q = 3600.0 * (180.0 / mono::M_PI); // 1弧度 = q角秒
const double q_pow2 = q * q;
const double f = spheriod->GetFlatting(); // 椭球扁率
const double a = spheriod->GetSemimajorAxis(); // 椭球长半轴
const double b = spheriod->GetSemiminorAxis(); // 椭球短半轴
const double e_pow2 = 2 * f - f * f; // (a * a - b * b) / a * a; 2 * f - f * f; // 第一篇心率的平方 e^2
const double e_sec_pow2 = (a / b) * (a / b) - 1; // (a * a - b * b) / b * b 此表达式存在大数值的平方运算 导致计算结果异常; // 第二偏心率的平方 e'^2
// 计算子午线弧长公式系数
double m0 = a * (1 - e_pow2);
double m2 = (3.0 / 2.0) * e_pow2 * m0;
double m4 = (5.0 / 4.0) * e_pow2 * m2;
double m6 = (7.0 / 6.0) * e_pow2 * m4;
double m8 = (9.0 / 8.0) * e_pow2 * m6;
double a0 = m0 + (1.0 / 2) * m2 + (3.0 / 8) * m4 + (5.0 / 16) * m6 + (35.0 / 128) * m8;
double a2 = (1.0 / 2) * m2 + (1.0 / 2) * m4 + (15.0 / 32) * m6 + (7.0 / 16) * m8;
double a4 = (1.0 / 8) * m4 + (3.0 / 16) * m6 + (7.0 / 32) * m8;
double a6 = (1.0 / 32) * m6 + (1.0 / 16) * m8;
double a8 = (1.0 / 128) * m8;
// 迭代因子
double B_f = 1.0;
double B_f_1 = 1.0;
double B_f_cur = 1.0;
double B_f_nxt = 1.0;
double B_f_delta = 0.0;
double F_B_f = 1.0;
double t_f = 1.0;
double t_f_pow2 = 1.0;
double eta_f_pow2 = 1.0;
double W_f = 1.0;
double M_f = 1.0;
double N_f = 1.0;
double B = 0;
double l = 0;
double L = 0;
std::vector<mono::SLonLatHeight3d> lbhs;
lbhs.resize(pnum);
// 精度计算：
// 1角秒 = 1/3600角度 = 1/3600 * pi/180弧度 = pi/648000弧度 = pi/648000 * 3600弧秒
// 1纬度大约等于111km ==> 1角秒大约等于30.833333米
// 为了确保厘米级精度，按照1角秒 / 100000000，此时对应距离大约等于30.833333 / 100000000 = 0.00000030833333米
// p2-p1 <= 0.00000030833333 表示p2 == p1
// 则精度计算阈值 = pi/648000 * 3600 * 1/100000000
const double tmp = mono::M_PI / 648000 * 3600 * 1e-8;
const int offset = GetNearIntegerOffset(tmp);
// const int tmp_int = (int)floor(tmp * pow(10, offset));
// const double precision = tmp_int / pow(10.0, offset);
const double precision = 1.0 / pow(10.0, offset);
// 遍历坐标点
double e_pos, n_pos;
for (size_t pi = 0; pi < pnum; pi++) {
const auto& enh = originProjs[pi];
e_pos = enh.x() - 500000; // 减去500km的偏移量
n_pos = enh.y();
B_f_1 = n_pos / a0;
B_f_nxt = B_f_1;
do {
B_f_cur = B_f_nxt;
F_B_f = -a2 / 2 * sin(2 * B_f_cur) + a4 / 4 * sin(4 * B_f_cur) - a6 / 6 * sin(6 * B_f_cur) + a8 / 8 * sin(8 * B_f_cur);
B_f_nxt = (n_pos - F_B_f) / a0;
B_f_delta = abs(B_f_nxt * 3600 - B_f_cur * 3600); // 弧度要转换弧秒 因为精度计算按照弧秒计算
} while (B_f_delta > precision);
B_f = B_f_nxt;
t_f = tan(B_f);
t_f_pow2 = t_f * t_f;
eta_f_pow2 = e_sec_pow2 * cos(B_f) * cos(B_f);
W_f = sqrt(1 - e_pow2 * sin(B_f) * sin(B_f));
M_f = a * (1 - e_pow2) / pow(W_f, 3);
N_f = a / W_f;
B = B_f - t_f / (2 * M_f) * e_pos * (e_pos / N_f) * (
1 -
t_f / 12 * (5 + 3 * t_f_pow2 + eta_f_pow2 - 9 * eta_f_pow2 * t_f_pow2) * pow(e_pos / N_f, 2) +
t_f / 360 * (61 + 90 * t_f_pow2 + 45 * t_f_pow2 * t_f_pow2) * pow(e_pos / N_f, 4));
l = 1 / cos(B_f) * (e_pos / N_f) * (
1 -
1.0 / 6 * (1 + 2 * t_f_pow2 + eta_f_pow2) * pow(e_pos / N_f, 2) +
1.0 / 120 * (5 + 28 * t_f_pow2 + 24 * t_f_pow2 * t_f_pow2 + 6 * eta_f_pow2 + 8 * eta_f_pow2 * t_f_pow2) * pow(e_pos / N_f, 4));
// 转换为角度制
B = B * 180 / mono::M_PI;
L = centralMeridian.GetRadixDegree() + l * 180 / mono::M_PI;
lbhs[pi].Set(L, B, enh.z());
}
targetLBHs = std::move(lbhs);
}
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
