高斯投影变形解析分析

Analytic Analysis of Gaussian Projection De-formation

摘要: 作为一种常用的等角投影方式，高斯投影被广泛用在大地测量学的各个应用领域中。传统的高斯投影方法将地球视为椭球体，其数学公式主要为经差的幂级数展开式，公式冗长且计算量大，计算精度不足，也不能很清晰地反映高斯投影的本质和投影变形规律。而将地球视为球体时，高斯投影与横墨卡托投影是等价的，因此高斯投影公式可表示为形式紧凑的闭合形式。针对传统方法中存在的问题，本文通过球面公式计算特殊点处的投影变形和一条6度带条带内的投影变形，以此分析高斯投影的投影变形规律。

Abstract: As a commonly used isometric projection, Gauss projection is widely used in various applications of geodesy. Regarding the earth as an ellipsoid, the traditional Gauss projection mathematical formu-la’s computation is too complex to clearly reflect the nature of Gauss’s projection and its distortion. When the Earth is considered a sphere, the Gauss projection is equivalent to the transverse Merca-tor projection, so the Gauss projection formula can be expressed as a compact closed form. In order to solve the problems in the traditional method, the projection deformation at a special point and the projection deformation in a 6-degree strip are calculated by the spherical formula in this paper, and the projection deformation rule of Gauss projection is analyzed.

1. 引言

大地测量学中经常使用的高斯投影的全名是高斯–克吕格投影，又叫等角横切椭圆柱投影。高斯投影一般满足以下三个条件 [1] ：1) 正形投影，即等角投影；2) 中央经线和赤道投影后卫互相垂直的直线，且为投影的对称轴；3) 中央子午线投影后长度保持不变。

在过去的投影方法中，一般将地球视为一个椭球体，将高斯投影用椭球面来表示 [2] [3] 。传统高斯投影数学公式虽然有容易理解和直观的优点，但是表达式相当冗长，计算繁琐，且被表示为经差的幂级数形式，适用范围有限 [4] 。因此，根据高斯投影的等角性质，国内外学者提出了高斯投影复变函数的表示形式 [5] - [15] 。相比于传统公式，复变函数表示形式的高斯投影数学公式精度更高，且具有“不分带”的优势，但引入了椭球偏心率e，表达式必须展开成椭球偏心率的幂函数，使得表达式也非常复杂。此外，由于该表示形式中包含复数变量，导致所得结果不够直观，难以直接对其进行代数分析，不能直观反应高斯投影的本质和高斯投影主要要素(例如长度比、子午线收敛角)的变化规律。

考虑到椭球面高斯投影方法存在的不足，近年来有多名学者对球面高斯投影方法展开了一定的研究 [16] [17] [18] 。该方法将地球视为一个球体，在此情况下高斯投影等价于横轴墨卡托投影 [19] ，因此有关高斯投影的一些重要公式可以表示为形式紧凑的闭合形式，方便数学计算。此外，考虑到常用的椭球偏心率非常小，因此在球面情况下更能有效地反映高斯投影的变化规律和数学分析性质 [20] 。本文以先前专家学者们对球面高斯投影的研究为基础，详细阐述球面高斯投影的数学分析及其应用，然后计算球面情况下特殊点处高斯投影变形情况，最后将球面方法和传统数学方法所得高斯投影变形结果进行对比，体现球面高斯投影的可靠性以及独特的优势。

2. 传统高斯投影表示式

2.1. 传统高斯投影正反解表示式

传统高斯投影公式通常表现为经差的幂级数形式。高斯投影正反解方法主要包括赫里斯托夫法和待定系数法，其中赫里斯托夫法高斯投影正解公式计算精度更高，当大地经差 $l\le \text{3}.\text{5}˚$ 时，该方法精度为0.001米。其正解公式如下：

$\{\begin{array}{l}x=X+Nt{\cos }^{2}B\frac{l^2}{{\rho }^2}\left[0.5+\frac{1}{24}\left(5-t^2+9{\eta }^2+4{\eta }^4\right){\cos }^{2}B\frac{l^2}{{\rho }^2}+\frac{1}{720}\left(61-58t^2+t^4\right){\cos }^{4}B\frac{l^4}{{\rho }^4}\right]\\ y=N\cos B\frac{l}{\rho }\left[1+\frac{1}{6}\left(1-t^2+{\eta }^2\right){\cos }^{2}B\frac{l^2}{{\rho }^2}+\frac{1}{120}\left(5-18t^2+t^4+14{\eta }^2-58{\eta }^2{t}^2\right)N{\cos }^{4}B\frac{l^4}{{\rho }^4}\right]\end{array}$ (1)

式中，x和y分别为高斯投影平面坐标，X为子午线弧长，B为大地纬度，l为经差，N为卯酉圈曲率半径， $\rho$ 为球心向径， $\eta ={e^{\prime }}^2{\cos }^{2}B$ ( $e^{\prime }$ 为椭球第二偏心率)， $t=\tan B$ 。

其反解公式如下：

$\{\begin{array}{l}B=B_f-\frac{\rho t_f}{2M_f}y\left(\frac{y}{N_f}\right)\left[1-\frac{1}{12}\left(5+3t_f^2+{\eta }_f^2-9{\eta }_f^2{t}_f^2\right){\left(\frac{y}{N_f}\right)}^2+\frac{1}{360}\left(61+90t_f^2+45t_f^4\right){\left(\frac{y}{N_f}\right)}^4\right]\\ l=\frac{\rho }{\cos B_f}\left(\frac{y}{N_f}\right)\left[1-\frac{1}{6}\left(1+2t_f^2+{\eta }_f^2\right){\left(\frac{y}{N_f}\right)}^2+\frac{1}{120}\left(5+28t_f^2+24t_f^4+6{\eta }_f^2+8{\eta }_f^2{t}_f^2\right){\left(\frac{y}{N_f}\right)}^4\right]\end{array}$ (2)

式中，B为大地纬度，l为经差， $B_f$ 为底点纬度。

2.2. 传统高斯投影长度比和子午线收敛角表示式

除了中央子午线外，其他位置上的任何线段在投影后均会产生长度变形，且离中央子午线越远，变形程度越大。设高斯投影平面上有一点，称该点领域内某线段投影后与投影前比值为长度比，一般用m表示。称该点的真子午线与位于此点所在的投影带的中央子午线之间的夹角，即在高斯平面上的真子午线与坐标纵线的夹角为子午线收敛角，一般用 $\gamma$ 表示。从子午线投影曲线量至纵轴，顺时针方向为正，反之为负。

根据投影长度变形定义，高斯投影长度比和子午线收敛角公式可表示为经差l的幂级数形式，取至 $l^4$ 项的长度比公式为

$m=1+\frac{1}{2}{l}^2{\cos }^{2}B\left(1+{\eta }^2\right)+\frac{1}{24}{l}^4{\cos }^{4}B\left(5-4{\tan }^{2}B\right)$ (3)

其中，B为大地纬度，l为经差， $\eta ={e^{\prime }}^2{\cos }^{2}B$ 。为了便于计算，子午线收敛角公式取至 $l^3$ 项，即

$\gamma =l\sin \phi +\frac{1}{3}{l}^3\sin \phi {\cos }^2\phi \left(1+3{\eta }^2\right)$ (4)

上式中的 $\gamma$ 的单位为弧度，且中央子午线以东 $\gamma$ 为负值，反之为正。

由上述数学公式可见，传统的高斯投影数学公式表达式相当冗长，计算起来十分繁琐。此外，传统公式普遍表示为经差的幂级数形式，该幂级数表达式使得高斯投影必须限制在一条狭窄的投影带内，精度还存在局限性。对于以上问题，本文引入球面高斯投影的概念，将在球面情况下，对高斯投影进行数学性质分析。

3. 球面高斯投影公式

3.1. 球面高斯投影正反解表示式

如图1所示，球面高斯投影是一种双重投影，主要步骤如下：1) 先将地球椭球面投影到一个等距离球面上，保持经度 不变，球面纬度 $\phi$ 是大地纬度 的函数；2) 椭球面投影到球面后，再通过横轴圆柱投影到平面上，保证中央经线和赤道的投影互相垂直。由于高斯投影需满足中央子午线投影后长度保持不变，因此椭球面投影到球面上时也必须要保持中央经线投影后长度不变的特征，可以得到该球面就必然是等距离纬度球面，即球面纬度 $\phi$ 就是等距离纬度 $\psi$ 。

由前人研究可知 [21] ，等距离纬度 $\psi$ 关于大地纬度B的正解展开式为

$\psi =B+{\alpha }_2\sin 2B+{\alpha }_4\sin 4B+{\alpha }_6\sin 6B+{\alpha }_8\sin 8B+{\alpha }_{10}\sin 10B$ (5)

系数为

$\{\begin{array}{l}{\alpha }_2=-\frac{3}{8}{e}^2-\frac{3}{16}{e}^4-\frac{111}{1024}{e}^6-\frac{141}{2048}{e}^8-\frac{1533}{32768}{e}^{10}\\ {\alpha }_4=\frac{15}{256}{e}^4+\frac{15}{256}{e}^6+\frac{405}{8192}{e}^8+\frac{165}{4096}{e}^{10}\\ {\alpha }_6=-\frac{35}{3072}{e}^6-\frac{35}{2048}{e}^8-\frac{4935}{262144}{e}^{10}\\ {\alpha }_8=\frac{315}{131072}{e}^8+\frac{315}{65536}{e}^{10}\\ {\alpha }_{10}=-\frac{693}{1310720}{e}^{10}\end{array}$ (6)

Figure 1. The principle of spherical Gauss projection

图1. 球面高斯投影方法

基于幂级数展开法的等距离纬度反解展开式为

$B=\psi +a_2\sin 2\psi +a_4\sin 4\psi +a_6\sin 6\psi +a_8\sin 8\psi +a_{10}\sin 10\psi$ (7)

系数为

$\{\begin{array}{l}{a}_2=\frac{3}{8}{e}^2+\frac{3}{16}{e}^4+\frac{213}{2048}{e}^6+\frac{255}{4096}{e}^8+\frac{20861}{524288}{e}^{10}\\ a_4=\frac{21}{256}{e}^4+\frac{21}{256}{e}^6+\frac{533}{8192}{e}^8+\frac{197}{4096}{e}^{10}\\ a_6=\frac{151}{6144}{e}^6+\frac{151}{4096}{e}^8+\frac{5019}{131072}{e}^{10}\\ a_8=\frac{1097}{131072}{e}^8+\frac{1097}{65536}{e}^{10}\\ a_{10}=\frac{8011}{2621440}{e}^{10}\end{array}$ (8)

和国内普遍使用高斯投影不同，国外大比例尺地形图普遍使用横轴墨卡托投影。将墨卡托投影的圆柱面旋转90˚后，使该圆柱面与球体的中央子午线相切，以此法得到的投影可视为横轴墨卡托投影。根据前人的研究，横轴墨卡托投影中将地球视为球体，取球面上一点坐标为 $\left(\phi ,\lambda \right)$ ， $\phi ,\lambda \in \left[-\frac{\pi }{2},\frac{\pi }{2}\right]$ ，对应的平面坐标为(x, y)，等距离球体的半径为 $R_{\psi }$ ，球面横轴墨卡托投影公式可以表示为：

$\{\begin{array}{l}x=R_{\psi }\arctan \left(\tan \phi \sec \lambda \right)\\ y=R_{\psi }\text{arctanh}\left(\cos \phi \sin \lambda \right)\end{array}$ (9)

其中，由前人研究可知 [14] ，等距离球体半径 $R_{\psi }$ 为

$R_{\psi }=\frac{2X\left(\frac{\pi }{2}\right)}{\pi }=a\left(1-\underset{n=1}{\overset{\infty }{{\displaystyle \sum }}}\frac{1}{2n-1}{\left(\frac{\left(2n-1\right)!!}{2^{n}n!}\right)}^2{e}^{2n}\right)$ (10)

式中，a为椭球长半轴，e是椭球第一偏心率，X是子午线弧长，公式为

$X={\displaystyle {\int }_0^{B}M\text{d}B}={\displaystyle {\int }_0^B\frac{a\left(1-e^2\right)}{{\left(1-e^2{\sin }^{2}B\right)}^{\frac{3}{2}}}\text{d}B}$ (11)

式中M为子午圈曲率半径。

从地球椭球面上投影到等距离球面上时，经度 $\lambda$ 不变，纬度 $\phi$ 是大地纬度B与椭球第一偏心率e的函数，即 $\phi =X/R_{\psi }$ 。由李忠美等学者的研究可知 [19] ，球面高斯投影和横轴墨卡托投影是等价的，因此球面高斯投影可以表示为闭合公式的形式。结合公式(9)以及等距离纬度关于大地纬度的正解展开式(5)，即等距离球面高斯投影正解公式。结合等距离纬度反解公式(7)，可得球面高斯投影的反解公式为：

$\{\begin{array}{l}\phi =\arcsin \left(\sin \left(x/R_{\psi }\right)\text{sech}\left(y/R_{\psi }\right)\right)\\ \lambda =\arctan \left(\sec \left(x/R_{\psi }\right)\sinh \left(y/R_{\psi }\right)\right)\end{array}$ (12)

该反解公式可以根据任意一点的平面坐标推导其球面坐标。由以上公式可知，等距离球面高斯投影可直接利用球面坐标计算，简化了高斯投影中的数学公式计算。

3.2. 长度比和子午线收敛角公式

由定义可知，横轴墨卡托投影都是等角投影，即从任意一点出发，各方向的长度比都是相等的，因此任一方向长度比即该点处的长度比。本文以经线方向为例( $\text{d}l=0$ )计算球面上任意点处的长度比。设ds为球面上的微分线段， $\text{d}{s}^{\prime }$ 为其所在平面的投影，r为球半径，可得

$\{\begin{array}{l}\text{d}s=r\text{d}\phi \\ \text{d}{s}^{\prime }=\sqrt{{\left(\text{d}x\right)}^2+{\left(\text{d}y\right)}^2}\\ \text{d}x=\frac{\partial x}{\partial \phi }\text{d}\phi +\frac{\partial x}{\partial \lambda }\text{d}\lambda \\ \text{d}y=\frac{\partial y}{\partial \phi }\text{d}\phi +\frac{\partial y}{\partial \lambda }\text{d}\lambda \end{array}$ (13)

该点的长度比m为

$m=\frac{\text{d}{s}^{\prime }}{\text{d}s}=2\sqrt{\frac{1}{3+\cos 2\lambda -2\cos 2\phi {\sin }^2\lambda }}=\frac{1}{\sqrt{1-{\cos }^2\phi {\sin }^2\lambda }}$ (14)

图2和图3分别为长度变形随纬度、经差的变化曲线，由图和以上公式可以得出结论：1) 球面高斯投影在其投影中央经线上不存在长度变形，而在其余的经线上，长度变形随着纬度的增加而减小；2) 除在极点外的同一条纬线上，长度变形随着经度的增加而增大，而中央经线和极区附近的长度变形较小；3) 球面高斯投影在经度较大(如 $\lambda =75˚$ )且纬度较小(如 $\phi =15˚$ )处长度变形明显增大。

根据定义，子午线收敛角可表示为

$\gamma =\arctan \left(\frac{\partial x}{\partial \lambda }/\frac{\partial y}{\partial \lambda }\right)$ (15)

Figure 2. Curves of length deformation ratio with latitude

图2. 长度变形比随纬度变化曲线

Figure 3. Curves of length deformation ratio with longitude

图3. 长度变形比随经差变化曲线

将球面高斯投影正解公式代入并化简，可得球面下的高斯投影子午线收敛角公式为

$\gamma =\arctan \left(\sin \phi \tan \lambda \right)$ (16)

图4和图5分别为子午线收敛角分别随纬度和经差的变化图。由图可得以下结论：1) 当经差一定时，子午线收敛角随着纬度的增加而增大；2) 当纬度一定时，子午线收敛角随着经差的增加而增大。

同理，考虑到高斯投影关于赤道以及中央子午线的对称性，可得以下结论：1) 在同一子午线上，子午线收敛角的绝对值随着点逐渐靠近极点而递增；2) 在同一纬线圈上，子午线收敛角的绝对值随着点逐渐远离中央子午线而递增。

4. 高斯投影变形解析分析

4.1. 特殊点处投影长度变形计算

根据第3节中所得的球面高斯投影长度变形和子午线收敛角闭合公式，可计算任意特殊点处的长度比和子午线收敛角数值。本文分别选取并计算纬度 $\phi$ 为0˚、15˚、30˚、45˚、60˚、75˚，经差 $\lambda$ 为0˚、15˚、30˚、45˚、60˚、75˚时对应的长度比，所得结果如表1所示。

Figure 4. Curves of meridian convergence with latitude

图4. 子午线收敛角随纬度变化曲线

Figure 5. Curves of meridian convergence with longitude

图5. 子午线收敛角随经差变化曲线

Table 1. Length ratios at special points

表1. 特殊点处长度比

由表1可知，在中央子午线上( $\lambda =0˚$ )的长度比始终为1，球面高斯投影不存在长度变形。由该表可直接获得部分重要特殊点处的投影变形比，比如在赤道上( $\phi =0˚$ )，当经度 $\lambda =45˚$ 时，该点处长度比为 $\sqrt{2}$ ，而经度 $\lambda =60˚$ 时，该点处长度比为2。通过分析各点处的长度比值也可以对应3.2节中总结的规律：1) 除了中央子午线处，当经度一定时，纬度越高，长度比越小；2) 当纬度一定时，经度越大，长度比越大。

4.2. 特殊点处子午线收敛角正切计算

本文分别选取并计算纬度 $\phi$ 为0˚、15˚、30˚、45˚、60˚、75˚，经差 $\lambda$ 为0˚、15˚、30˚、45˚、60˚、75˚时对应的子午线收敛角的正切值，所得结果如表2所示。

Table 2. Tangent values of meridian convergence at special points

表2. 特殊点处子午线收敛角正切值

由表2可知，在中央子午线上( $\lambda =0˚$ )和在赤道上( $\phi =0˚$ )的任意一点的子午线收敛角的正切值始终为0。由该表可直接获得部分重要特殊点处的子午线收敛角正切值，比如在经度 $\lambda =45˚$ 、纬度 $\phi =30˚$ ，或在经度 $\lambda =30˚$ 、纬度 $\phi =60˚$ 处，子午线收敛角的正切值为0.5。通过分析各点处的子午线收敛角正切值可知：1) 除了中央子午线处，当经度一定时，纬度越高，子午线收敛角正切值越大；2) 在赤道以外，当纬度一定时，经度越大，子午线收敛角正切值越大。以上结论与3.2节中所得结论一致。

综上，由以上特殊点处的长度比和子午线收敛角的值可以证明3.2节中规律分析的正确性。此外，由两表可知，在纬度和经差较大时，球面高斯投影仍然有对应的确定性表示，可以有效避免传统公式受限于分带的问题。

4.3. 一个6度带内的误差计算

为验证球面高斯投影公式的正确性，本文将球面高斯投影(球面横轴墨卡托投影)的数学公式同传统的幂级数公式进行比较，主要计算在一个6度带条带内，球面高斯投影长度比和子午线收敛角闭合公式与传统公式在球面情况下的差值表达式。

4.3.1. 高斯投影长度比公式的差值计算

取球面上一点为 $\left(\phi ,\lambda \right)$ ，球面高斯投影长度比公式所得结果记为 $m_1$ ，传统幂级数公式计算结果记为 $m_2$ ，即

$\{\begin{array}{l}{m}_1=\frac{1}{\sqrt{1-{\cos }^2\phi {\sin }^2\lambda }}\\ m_2=1+\frac{1}{2}{\lambda }^2{\cos }^2\phi +\frac{1}{24}{\lambda }^4{\cos }^4\phi \left(5-4{\tan }^2\phi \right)\end{array}$ (17)

可得球面高斯投影长度比公式和传统公式的差值 $\Delta m=m_1-m_2$ 。本文分别选取并计算球面纬度 $\phi$ 为0˚、15˚、30˚、45˚、60˚、75˚，经差 $\lambda$ 为0.5˚、1˚、1.5˚、2˚、2.5˚、3˚时对应的球面高斯投影长度比公式和传统公式的差值，所得结果如表3所示。

Table 3. The differences between the closed formula and the traditional power series expansion of length ratio

表3. 球面形式的高斯投影长度比闭合公式与传统幂级数展开式之差

由表3和图6、图7可知，在经差 $\lambda$ 在0˚到3˚范围内时，球面高斯投影长度比闭合公式和传统幂级数展开式所求结果的绝对差值均小于10⁻⁸，而当经差 $\lambda$ 值接近或等于3˚时，两者的绝对差值随着纬度 $\phi$ 的减小和迅速增加。

Figure 6. Absolute differences between the closed formula and the traditional power series expansion of length ratio

图6. 球面长度比公式和传统长度比公式的绝对差值

4.3.2. 高斯投影子午线收敛角公式的差值计算

同4.3.1，球面高斯投影长度比公式所得结果记为 ${\gamma }_1$ ，传统幂级数公式计算结果记为 ${\gamma }_2$ ，即

$\{\begin{array}{l}{\gamma }_1=\arctan \left(\sin \phi \tan \lambda \right)\\ {\gamma }_2=\lambda \sin \phi +\frac{1}{3}{\lambda }^3\sin \phi {\cos }^2\phi \end{array}$ (18)

可得球面高斯投影子午线收敛角公式和传统高斯投影公式的差值 $\Delta \gamma ={\gamma }_1-{\gamma }_2$ 。分别选取并计算球面纬度 $\phi$ 为0˚、15˚、30˚、45˚、60˚、75˚，经差 $\lambda$ 为0.5˚、1˚、1.5˚、2˚、2.5˚、3˚时对应的球面高斯投影子午线收敛角公式和传统高斯投影公式的差值，幂级数形式取到经差的三次方项，所得结果如表4所示，单位为弧度。

Figure 7. The contour map of absolute differences between the closed formula and the traditional power series of length ratio

图7. 球面长度比公式和传统长度比公式的绝对差值的等值线图

Table 4. The differences between the closed formula and the traditional power series expansion of meridian convergence

表4. 球面形式的高斯投影子午线收敛角闭合公式与传统幂级数展开式之差

由表4和图8、图9可知，在经差 $\lambda$ 在0˚到3˚范围内时，球面高斯投影子午线收敛角闭合公式和传统幂级数展开式所求结果的绝对差值均小于10⁻⁷ rad，而当经差 $\lambda$ 值大约在2˚和3˚之间、纬度 $\phi$ 大约在0˚和50˚之间时，两者的绝对差值随着经差的增加而迅速增加。

Figure 8. Absolute differences between the closed formula and the traditional power series expansion of meridian convergence

图8. 球面子午线收敛角公式和传统公式的绝对差值

Figure 9. The contour map of absolute differences between the closed formula and the traditional power series of meridian convergence

图9. 球面子午线收敛角公式和传统公式的绝对差值的等值线图

综上所述，球面高斯投影长度比和子午线收敛角闭合公式精度较高且具有高可靠性，完全可以在误差允许的范围内代替经典的高斯投影幂级数展开式，可以有效简化计算量，并进行可视化分析。

5. 结论

本文分析了球面高斯投影的原理及重要公式，通过计算球面情况下高斯投影在特殊点处的长度比和子午线收敛角正切值，并将球面公式和传统公式进行对比，可得以下结论：

1) 传统高斯投影公式主要表示为经差的幂级数形式，计算量大且受经差的限制，而球面高斯投影闭合公式基本上不受经差的限制，并且能进行一定的可视化分析，可以在较大程度上反映高斯投影的变化规律；

2) 高斯投影是等角投影，投影到平面上会产生长度变形，本文计算了特殊点处高斯投影的长度比和子午线收敛角的正切值，一定程度上反映了高斯投影变形的变化规律；

3) 在一个6˚带条带内，球面闭合形式的长度比和子午线收敛角公式可以在误差允许的范围内替代经典的椭球高斯投影级数展开式。

基金项目

国家自然科学基金项目42074010。

NOTES

^*通讯作者。