详细介绍：测绘工具箱（3）高斯投影下坐标的正算与反算

高斯投影下坐标的正算与反算

引言

大地坐标以经纬度形式描述地球椭球面上的点位，而高斯投影坐标则将椭球面投影到平面上，以直角坐标形式表示。这种转换广泛应用于地形测量​​、​​地图制图​​、​​工程测量​​以及​​GPS数据处理​​中。

一、核心参数

文中的计算以WGS84椭球为例。

  • ​​长半轴 a​​: 椭球赤道半径（WGS84为​​6378137.0米​​）
  • ​​扁率 f​​: 椭球的扁率（WGS84为​​1/298.257223563​​）
  • ​​第一偏心率平方 e²​​: e² = 2f - f²
  • ​​第二偏心率平方 e’²​​: e’² = e² / (1 - e²)
  • ​​​​中央子午线经度 (λ₀)​​由带号决定，也可自定义中央子午线
  • ​​比例因子 k0​​: ​​1
  • ​​假东偏移 (False Easting)​​: ​​500000米​​
  • ​​假北偏移 (False Northing)​​: 北半球为0，南半球为​​10000000米​​

二、高斯投影

高斯投影（Gauss-Krüger Projection）是一种​​横轴等角切椭圆柱投影​​，由德国数学家高斯提出并经克吕格完善，故又称高斯-克吕格投影。其基本特性包括：

  • ​​等角性质​​：投影后保持无限小图形的形状不变，即角度不变形。
  • ​​中央子午线为直线​​：且其长度投影后保持不变。
  • ​​投影变形​​：中央子午线无长度变形，离中央子午线越远，长度变形越大。

为了控制变形，地球表面通常按经差分带投影。常用分带方式包括：

  • ​​6°分带​​：从首子午线开始，每隔经差6°划分为一带，全球共60带。带号 N与中央子午线经度 L0 的关系为：.
    N = ⌊ L + 3 6 ⌋ , L 0 = 6 N − 3 N = \left\lfloor \frac{L + 3}{6} \right\rfloor, \quad L_{0} = 6N - 3 N=⌊6L+3​⌋,L0​=6N−3

  • 3°分带​​：从东经1.5°开始，每隔经差3°划分为一带，全球共120带。带号N与中央子午线经度 L0 的关系为：
    n = ⌊ L 3 + 0.5 ⌋ , L 0 = 3 n n = \left\lfloor\frac{L}{3} + 0.5\right\rfloor, \quad L_0 = 3n n=⌊3L​+0.5⌋,L0​=3n

在高斯平面直角坐标系中，​​X轴​​是中央子午线的投影，向北为正；​​Y轴​​是赤道的投影，向东为正。为了避免横坐标出现负值，通常会在Y坐标值上加上500000米，并在坐标前冠以带号。

三、高斯正算

高斯投影正算是指由大地坐标 (B,L)计算高斯平面直角坐标 (x,y)的过程。基本步骤如下

1. 计算经差

根据点的大地经度 L和所选用的分带方 ​​计算经差​​，需转换为弧度制。
l = L − L 0 l=L − L_{0} l=L−L0​

  • ​​辅助量计算​​：包括卯酉圈半径 N、子午线弧长 X等。

2. 计算卯酉圈半径 N

N = a 1 − e 2 sin ⁡ 2 B , η 2 = e ′ 2 cos ⁡ 2 B N=\frac{a}{\sqrt{1-e^{2}\sin^{2}B}},\quad\eta^{2}=e^{\prime 2}\cos^{2}B N=1−e2sin2B​a​,η2=e′2cos2B

3. 计算子午线弧长 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 - \frac{a_{2}}{2}\sin(2B) + \frac{a_{4}}{4}\sin(4B) - \frac{a_{6}}{6}\sin(6B) + \frac{a_{8}}{8}\sin(8B) X=a0​B−2a2​​sin(2B)+4a4​​sin(4B)−6a6​​sin(6B)+8a8​​sin(8B)

  • B为​​大地纬度​​，以弧度（rad）为单位。
  • a0 ,a2 ,a4 ,a6 ,a8 为与椭球参数有关的系数，通过中间量mi计算得出
    { 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 \begin{cases} m_{0} = a(1 - e^{2}) \\ m_{2} = \frac{3}{2}e^{2}m_{0} \\ m_{4} = \frac{5}{4}e^{2}m_{2} \\ m_{6} = \frac{7}{6}e^{2}m_{4} \\ m_{8} = \frac{9}{8}e^{2}m_{6} \end{cases} ⎩⎨⎧​m0​=a(1−e2)m2​=23​e2m0​m4​=45​e2m2​m6​=67​e2m4​m8​=89​e2m6​​
  • 然后计算ai
    { 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{cases} 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{cases} ⎩⎨⎧​a0​=m0​+21​m2​+83​m4​+165​m6​+12835​m8​+⋯a2​=21​m2​+21​m4​+3215​m6​+167​m8​a4​=81​m4​+163​m6​+327​m8​a6​=321​m6​+161​m8​a8​=1281​m8​​

4. 计算投影坐标

x = X + N 2 sin ⁡ B cos ⁡ B ⋅ l 2 + N 24 sin ⁡ B cos ⁡ 3 B ( 5 − tan ⁡ 2 B + 9 η 2 ) l 4 + N 720 sin ⁡ B cos ⁡ 5 B ( 61 − 58 tan ⁡ 2 B + tan ⁡ 4 B ) l 6 y = N cos ⁡ B ⋅ l + N 6 cos ⁡ 3 B ( 1 − tan ⁡ 2 B + η 2 ) l 3 + N 120 cos ⁡ 5 B ( 5 − 18 tan ⁡ 2 B + tan ⁡ 4 B + 14 η 2 − 58 tan ⁡ 2 B ⋅ η 2 ) l 5 \begin{aligned} x &= X + \frac{N}{2}\sin B \cos B \cdot l^{2} + \frac{N}{24}\sin B \cos^{3} B \left(5 - \tan^{2} B + 9\eta^{2}\right) l^{4} + \frac{N}{720}\sin B \cos^{5} B \left(61 - 58\tan^{2} B + \tan^{4} B\right) l^{6} \\ y &= N \cos B \cdot l + \frac{N}{6}\cos^{3} B \left(1 - \tan^{2} B + \eta^{2}\right) l^{3} + \frac{N}{120}\cos^{5} B \left(5 - 18\tan^{2} B + \tan^{4} B + 14\eta^{2} - 58\tan^{2} B \cdot \eta^{2}\right) l^{5} \end{aligned} xy​=X+2N​sinBcosB⋅l2+24N​sinBcos3B(5−tan2B+9η2)l4+720N​sinBcos5B(61−58tan2B+tan4B)l6=NcosB⋅l+6N​cos3B(1−tan2B+η2)l3+120N​cos5B(5−18tan2B+tan4B+14η2−58tan2B⋅η2)l5​

  • 辅助参数计算
    η 2 = e ′ 2 cos ⁡ 2 B \eta^{2}=e^{\prime 2}\cos^{2}B η2=e′2cos2B
    其中，e为第一偏心率，e′ 为第二偏心率。
// 转换为弧度
double B = bLCoordinate.Latitude * Math.PI / 180.0;
double L = bLCoordinate.Longitude * Math.PI / 180.0;
double L0 = projection.CentralMeridian * Math.PI / 180.0;
double Scale = projection.ScaleFactor;
double l = L - L0;
// 经差
// 计算辅助量
double N = ellipsoid.SemiMajorAxis / Math.Sqrt(1 - ellipsoid.EccentricitySquared * Math.Pow(Math.Sin(B), 2));
double t = Math.Tan(B);
//double eta2 = ellipsoid.EccentricitySquared / (1 - ellipsoid.EccentricitySquared) * Math.Pow(Math.Cos(B), 2);
double eta2 = ellipsoid.SecondEccentricitySquared * Math.Pow(Math.Cos(B), 2);
// 计算子午线弧长
double m0 = ellipsoid.SemiMajorAxis * (1 - ellipsoid.EccentricitySquared);
double m2 = 3.0 / 2.0 * ellipsoid.EccentricitySquared * m0;
double m4 = 5.0 / 4.0 * ellipsoid.EccentricitySquared * m2;
double m6 = 7.0 / 6.0 * ellipsoid.EccentricitySquared * m4;
double m8 = 9.0 / 8.0 * ellipsoid.EccentricitySquared * m6;
double a0 = m0 + m2 / 2.0 + 3.0 / 8.0 * m4 + 5.0 / 16.0 * m6 + 35.0 / 128.0 * m8;
double a2 = m2 / 2.0 + m4 / 2.0 + 15.0 / 32.0 * m6 + 7.0 / 16.0 * m8;
double a4 = m4 / 8.0 + 3.0 / 16.0 * m6 + 7.0 / 32.0 * m8;
double a6 = m6 / 32.0 + m8 / 16.0;
double a8 = m8 / 128.0;
double X = Scale * (a0 * B - a2 * Math.Sin(2 * B) / 2.0 + a4 * Math.Sin(4 * B) / 4.0
- a6 * Math.Sin(6 * B) / 6.0 + a8 * Math.Sin(8 * B) / 8.0);
// 高斯投影公式
X += N * t * (Math.Pow(l, 2) / 2 * Math.Pow(Math.Cos(B), 2) +
Math.Pow(l, 4) / 24 * Math.Pow(Math.Cos(B), 4) *(5 - Math.Pow(t, 2) + 9 * eta2));
double Y = Scale * N * (l * Math.Cos(B) +
Math.Pow(l, 3) / 6 * Math.Pow(Math.Cos(B), 3) * (1 - Math.Pow(t, 2) + eta2) +
Math.Pow(l, 5) / 120 * Math.Pow(Math.Cos(B), 5) * (5 - 18 * Math.Pow(t, 2) + Math.Pow(t, 4) + 14 * eta2 -58 * t * t * eta2));

在这里插入图片描述

四、高斯反算

1. 去除虚假位移

y = y − 500000 m \mathrm{y} = \mathrm{y} - 500000 m y=y−500000m
或​​剥离带号与偏移值​​：将通用横坐标转换为自然坐标
y 自然 = y 通用 − 500 ,  ⁣ 000 − 带号 × 1 0 6 y_{\text{自然}} = y_{\text{通用}} - 500,\!000 - \text{带号} \times 10^{6} y自然​=y通用​−500,000−带号×106

2. 迭代计算子午线弧长

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 - \frac{a_{2}}{2}\sin(2B) + \frac{a_{4}}{4}\sin(4B) - \frac{a_{6}}{6}\sin(6B) + \frac{a_{8}}{8}\sin(8B) X=a0​B−2a2​​sin(2B)+4a4​​sin(4B)−6a6​​sin(6B)+8a8​​sin(8B)

  • B为​​大地纬度​​，以弧度（rad）为单位。
  • a0 ,a2 ,a4 ,a6 ,a8 为与椭球参数有关的系数，通过中间量mi计算得出。参数计算方式与正算相同。

3. 计算辅助参数

η 2 = e ′ 2 cos ⁡ 2 B f \eta^{2}=e^{\prime 2}\cos^{2}B_{f} η2=e′2cos2Bf​
N f = a 1 − e 2 sin ⁡ 2 B f N_{f}=\frac{a}{\sqrt{1-e^{2}\sin^{2}B_{f}}} Nf​=1−e2sin2Bf​​a​
M f = a ( 1 − e 2 ) ( 1 − e 2 sin ⁡ 2 B f ) 1.5 M_{f} = \frac{a \left(1 - e^{2}\right)}{\left(1 - e^{2} \sin^{2}B_{f}\right)^{1.5}} Mf​=(1−e2sin2Bf​)1.5a(1−e2)​

4. 计算大地坐标

B = B f − tan ⁡ B f 2 M f N f y 2 + tan ⁡ B f 24 M f N f 3 ( 5 + 3 tan ⁡ 2 B f + η f 2 − 9 η f 2 tan ⁡ 2 B f ) y 4 l = y N f cos ⁡ B f − y 3 6 N f 3 cos ⁡ B f ( 1 + 2 tan ⁡ 2 B f + η f 2 ) + y 5 120 N f 5 cos ⁡ B f ( 5 + 28 tan ⁡ 2 B f + 24 tan ⁡ 4 B f ) \begin{aligned} B &= B_{f} - \frac{\tan B_{f}}{2M_{f}N_{f}}y^{2} + \frac{\tan B_{f}}{24M_{f}N_{f}^{3}}(5 + 3\tan^{2}B_{f} + \eta_{f}^{2} - 9\eta_{f}^{2}\tan^{2}B_{f})y^{4} \\ l &= \frac{y}{N_{f}\cos B_{f}} - \frac{y^{3}}{6N_{f}^{3}\cos B_{f}}(1 + 2\tan^{2}B_{f} + \eta_{f}^{2}) + \frac{y^{5}}{120N_{f}^{5}\cos B_{f}}(5 + 28\tan^{2}B_{f} + 24\tan^{4}B_{f}) \end{aligned} Bl​=Bf​−2Mf​Nf​tanBf​​y2+24Mf​Nf3​tanBf​​(5+3tan2Bf​+ηf2​−9ηf2​tan2Bf​)y4=Nf​cosBf​y​−6Nf3​cosBf​y3​(1+2tan2Bf​+ηf2​)+120Nf5​cosBf​y5​(5+28tan2Bf​+24tan4Bf​)​

// 移除东伪偏移并转换为以中央子午线为原点的坐标
double y = xy.YDirection - proj.YOffset;
double x = xy.XDirection;
// 1. 计算辅助参数
double a = ellipsoid.SemiMajorAxis;
double b = ellipsoid.SemiMinorAxis;
double e2 = ellipsoid.EccentricitySquared;
double e4 = e2 * e2;
double e6 = e4 * e2;
// 2. 计算Footprint纬度Bf（迭代求解）
// 计算子午线弧长幂次分解
double m0 = ellipsoid.SemiMajorAxis * (1 - ellipsoid.EccentricitySquared);
double m2 = 3.0 / 2.0 * ellipsoid.EccentricitySquared * m0;
double m4 = 5.0 / 4.0 * ellipsoid.EccentricitySquared * m2;
double m6 = 7.0 / 6.0 * ellipsoid.EccentricitySquared * m4;
double m8 = 9.0 / 8.0 * ellipsoid.EccentricitySquared * m6;
double a0 = m0 + m2 / 2.0 + 3.0 / 8.0 * m4 + 5.0 / 16.0 * m6 + 35.0 / 128.0 * m8;
double a2 = m2 / 2.0 + m4 / 2.0 + 15.0 / 32.0 * m6 + 7.0 / 16.0 * m8;
double a4 = m4 / 8.0 + 3.0 / 16.0 * m6 + 7.0 / 32.0 * m8;
double a6 = m6 / 32.0 + m8 / 16.0;
double a8 = m8 / 128.0;
double Bf = x / a0;
double FBf;
double B0 = Bf;
//迭代求数值为x坐标的子午线弧长对应的底点纬度
while ((Math.Abs(Bf - B0) >
0.0000001) || (B0 == Bf))
{
B0 = Bf;
FBf = -a2 / 2.0 * Math.Sin(2 * B0) + a4 / 4.0 * Math.Sin(4 * B0) - a6 / 6.0 * Math.Sin(6 * B0);
Bf = (x - FBf) / a0;
}
// 3. 计算辅助参数
double sinBf = Math.Sin(Bf);
double cosBf = Math.Cos(Bf);
double tanBf = Math.Tan(Bf);
double eta2 = e2 / (1 - e2) * cosBf * cosBf;
double Nf = ellipsoid.SemiMajorAxis / Math.Sqrt(1 - e2 * sinBf * sinBf);
double Mf = ellipsoid.SemiMajorAxis * (1 - e2) / Math.Pow(1 - e2 * sinBf * sinBf, 1.5);
// 4. 计算经差l和纬度B
double l = y / (Nf * cosBf) * (1 -
(1 + 2 * tanBf * tanBf + eta2) * y * y / (6 * Nf * Nf) +
(5 + 28 * tanBf * tanBf + 24 * Math.Pow(tanBf, 4)) * Math.Pow(y, 4) / (120 * Math.Pow(Nf, 4)));
double B = Bf - tanBf * y * y / (2 * Mf * Nf) * (1 -
(5 + 3 * tanBf * tanBf + eta2 - 9 * eta2 * tanBf * tanBf) * y * y / (12 * Nf * Nf));

在这里插入图片描述

posted @ 2025-09-11 19:59  yfceshi  阅读(178)  评论(0)    收藏  举报