异面高斯投影变换与Excel VBA实现——以实数等角纬度为变换参数

Gaussian Projection Transformation of Different Surfaces and Its Implementation in Excel VBA—Using Real Isometric Latitude as Transformation Parameters

摘要: 随着我国工程建设的发展，对地球投影变换精度提出了更高的要求，尤其是高铁工程建设，其投影变换精度要求< 1 m/100 km。在传统标准椭球高斯投影变换的基础上，增加椭球膨胀计算模块，组构了异面高斯投影变换：反算、正算及换带。正反算变换均利用正算式计算，正算采用非迭代法，反算采用迭代法。参照大地纬度B表达的高斯变换式，本文给出了同一表达格式的基于等角纬度φ表达的高斯投影变换式。从实用角度出发，利用Excel VBA计算平台，以实数等角纬度为变换参数，采用大地坐标增量法及空间直角坐标法作为膨胀变换方法，开发了异面高斯投影变换宏，一键实现坐标的批量变换。椭球膨胀坐标偏位$\left(D_x,D_y\right)$ 主要由椭球长半轴增量$da$ 所贡献，且与$da$ 呈线性增长。

Abstract: With the development of engineering construction in our country, higher requirements have been put forward for the accuracy of earth projection transformation, especially for high-speed railway engineering construction, which requires a projection transformation accuracy of < 1 m/100 km. On the basis of the traditional standard ellipsoid Gaussian projection transformation, an ellipsoid dilation calculation module is added to construct different surface Gaussian projection transformations: inverse calculation, forward calculation, and strip transformation. Both forward and backward transformations are calculated using forward equations, with non iterative methods used for forward calculations and iterative methods used for backward calculations. Referring to the Gaussian transform expression of geodetic latitude B, this article presents a Gaussian projection transform expression based on equiangular latitude φ expression using the same expression format. From a practical perspective, using the Excel VBA computing platform, taking real isometric latitude as the transformation parameter, and adopting the geodetic coordinate increment method and spatial Cartesian coordinate method as the expansion transformation method, different surface Gaussian projection transformation macros were developed to achieve batch transformation of coordinates with one click. The offset $\left(D_x,D_y\right)$ of ellipsoidal expansion coordinates is mainly contributed by the increment $da$ of the ellipsoidal semi axis, and it increases linearly with $da$ .

1. 引言

随着我国工程建设的发展，对地图投影变换的精度提出了更高的要求，尤其是铁路工程建设，其投影变换精度要求< 1 m/100 km。由此，高斯投影需要建立城市独立坐标系或工程独立坐标系。

独立坐标系或工程坐标系，一般是在标准椭球高斯投影的基础上，选定中央子午线，必要时进行高程面抵偿以最大限度地拟合投影区内的平均高程。

高程面抵偿即是椭球膨胀。进入本世纪，许多学者对膨胀变换建立独立坐标系、坐标变换及其应用进行了研究。邓兴升、汤仲安等(2010)为建立区域独立坐标系初步研究了椭球膨胀法的高斯投影正反解算法[1]；李亚辉、石德斌等(2010)为解决高速铁路的线路设计研究了高速铁路坐标系的搭接问题[2]；梅熙、王国祥(2012)探讨了高速铁路坐标的转换方法[3]；赵淑湘(2015)针对国家3˚带坐标系投影长度超限的城市，阐述了椭球膨胀法建立独立坐标系的方法和步骤[4]；李亚磊(2016)针对东西向高铁的平面控制投影方式进行了详细研究[5]；舒鹏瑞、蒋璐阳等(2016)基于GPS-RTK技术及椭球膨胀法建立地方坐标系[6]；蔡婧、郭际明等(2019)研究了利用椭球变换建立独立坐标系，研究了椭球膨胀法、变形法、平移法变换[7]；陈春(2019)，李文军、李军等(2020)针对跨度大、高差大的长距离航道工程应用椭球膨胀法建立航道工程独立坐标系[8] [9]；于志威(2020)，姜韶、张坤先等(2023)应用椭球膨胀法建立独立坐标系，编制了C#语言转换程序[10] [11]；陈再辉(2022)利用空间直角坐标对独立坐标系的椭球变换与坐标换算进行了研究[12]；关东润(2022)，孟静娟、赵辉等(2024)为解决投影变形超限研究了高程抵偿坐标系的方法及独立坐标系的建立[13] [14]。

高斯投影变换，李厚朴、边少锋等(2007, 2009, 2012)研究了基于复数等角纬度为参数的高斯投影变换[15]-[17]；李立新(2017)导出了基于实数等角纬度的反变换非迭代变换式[18]；刘大海(2024)基于实数等角纬度的高斯投影反变换采用迭代法计算[19]。

本文在前人研究的基础上，从方便快捷、批量转换的实用角度出发，以实数等角纬度作为高斯变换参数，利用Excel VBA开发了适用于BJ54、XIAN80、WGS84、CGCS2000坐标系的椭球膨胀异面高斯投影变换，正反算变换均利用正算式计算，正算采用非迭代法，反算采用迭代法，计算简明。进一步，分析了椭球膨胀偏位，x及y的主要偏位均与椭球长半轴增量$da$ 成线性增长。

2. 异面高斯投影变换框图

传统的高斯投影变换是同一标准椭球面上的高斯投影变换。异面高斯投影变换是在不同椭球面上进行的高斯投影变换，即在传统同面高斯投影变换的基础上，加入椭球膨胀变换，组构了异面高斯投影变换，见图1。

异面高斯反算：$E_1$ 面高斯反算 + 膨胀变换；

异面高斯正算：膨胀变换 + $E_2$ 面高斯正算；

异面高斯换带：$E_1$ 面高斯反算 + 膨胀变换 + $E_2$ 面高斯正算。

Figure 1. Block diagram of Gaussian projection transformation with different surfaces

图1. 异面高斯投影变换框图

3. 异面高斯投影变换算法

3.1. 椭球膨胀变换

椭球膨胀变换，是异面高斯投影变换的桥梁，它是由$E_1$ 椭球向$E_2$ 椭球的膨胀变换。

设$E_0$ 为标准椭球，椭球长半径为$a_0$ 。椭球膨胀时，设定：球心位置、椭球定向、椭球扁率、偏心率不变，不平移，不旋转，仅椭球长半轴膨胀变化$da$ 。椭球$E_1$ 的投影面大地高为$H_{E1}$ ，椭球长半径为$a_1$ ，椭球E₂的投影面大地高为$H_{E2}$ ，椭球长半径为$a_2$ 。变换点P分别对椭球$E_0$ 、$E_1$ 、$E_2$ 作法线，其P点的相应大地纬度为$B_0$ 、$B_1$ 、$B_2$ ，见图2。

Figure 2. Schematic diagram of ellipsoidal expansion

图2. 椭球膨胀示意图

椭球的几个相关参量：

$da=a_2-a_1=H_{E2}-H_{E1}$

$e=\sqrt{1-{\left(1-f\right)}^2}$

$e_2=\frac{e}{1-f}$

$W=\sqrt{1-{\left(e\cdot \sin B\right)}^2}$

$N=\frac{a}{\sqrt{1-{\left(e\cdot \sin B\right)}^2}}$

$M=a\cdot \left(1-e^2\right)\cdot {\sqrt{1-{\left(e\cdot \sin B\right)}^2}}^{-3}$

式中：

$H_{E1}$ ：$E_1$ 椭球投影面大地高；

$H_{E2}$ ：$E_2$ 椭球投影面大地高；

$da$ ：椭球长半轴膨胀量；

$a$ ：椭球长半径(变量)；

$a_0$ ：$E_0$ 参考椭球长半径；

$a_1$ ：$E_1$ 椭球长半径；

$a_2$ ：$E_2$ 椭球长半径；

$f$ ：椭球扁率；

$e$ ：椭球偏心率；

$e_2$ ：椭球第2偏心率；

$W$ ：椭球辅助参数；

$M$ ：子午圈曲率半径；

$N$ ：卯酉圈曲率半径。

3.1.1. 大地坐标增量法

大地坐标增量法，是根据椭球长半轴增量$da$ ，计算由$E_1$ 投影面膨胀到$E_2$ 投影面的大地坐标增量，从而确定转换点P在$E_2$ 椭球上的大地坐标[5]：

$\{\begin{array}{l}dL=0\\ dB=\frac{e^2\sin B\cos B}{W\left(M+H\right)}\\ dH=-W\cdot da\end{array}$ (1)

$\{\begin{array}{l}{B}_2=B_1+dB\\ L_2=L_1+dL\\ H_2=H_1+dH\end{array}$ (2)

3.1.2. 空间直角坐标法

空间直角坐标法，是根据投影面$E_1$ 上的大地坐标$\left(B_1,L_1,H_1\right)$ ，转换为空间坐标$\left(X,Y,Z\right)$ ，再由空间坐标转换为P点在$E_2$ 上的大地坐标$\left(B_2,L_2,H_2\right)$ 。

1) 空间法正算：由$E_1$ 椭球的大地坐标$\left(B_1,L_1,H_1\right)$ 转换为空间直角坐标$\left(X,Y,Z\right)$ 。

$E_1$ 椭球参量：

$\{\begin{array}{l}a=a_1\\ B=B_1\\ L=L_1\\ H=H_1\end{array}$ (3)

计算空间直角坐标：

$\{\begin{array}{l}X=\left(N+B\right)\cos B\cos L\\ Y=\left(N+B\right)\cos B\sin L\\ Z=\left[N\left(1-e^2\right)+H\right]\sin B\end{array}$ (4)

2) 空间法反算：由空间直角坐标$\left(X,Y,Z\right)$ 转换为$E_2$ 椭球的大地坐标$\left(B_2,L_2,H_2\right)$ 。一般用迭代法求算。

$E_2$ 椭球参量初值：

$\{\begin{array}{l}a=a_2\\ B=B_1\\ L=L_1\end{array}$ (5)

迭代法求解$\left({\text{N}}_2{\text{,B}}_2{\text{,H}}_2\right)$ ：

$\{\begin{array}{l}N=\frac{a}{\sqrt{1-{\left(e\sin B\right)}^2}}=\frac{a}{W}\\ B=\arctan \left[\frac{Z+N\cdot e^2\sin B}{\sqrt{X^2+Y^2}}\right]\\ H=\frac{\sqrt{X^2+Y^2}}{\cos B}-N\end{array}$ (6)

3.2. 高斯投影变换

高斯投影变换，为同一投影面($E_1$ 面，或$E_2$ 面)的高斯投影变换。

正反算变换，李厚朴、边少锋等(2009, 2012, 2018, 2021)给出了复数等角纬度的非迭代式[15]-[17] [20] [21]；金立新(2017)给出了以实数等角纬度的非迭代变换式的最新研究成果[18]。非迭代变换，需要不同的变换式及级数系数，表达冗长。

本文采用金立新(2017)以实数等角纬度为变换参数的高斯投影正算式[18]。正算变换采用金氏正算式，反算变换采用金氏正算式迭代计算[19]。以Excel VBA为计算工具，以椭球膨胀为模型(增量法及空间法)，实现异面高斯投影的批量变换。

特别地，基于等角纬度的高斯投影变换，不限于7˚投影带宽。仅从转换精度考虑，无需进行分带投影[19]。

3.2.1. E₁面高斯反算(迭代算法)

本文高斯反算，采用金氏正算式迭代计算，仅需1套变换式及级数系数，计算简洁。

1) 反算主变换：$\left({\phi }_x,{\phi }_y\right)\leftarrow \left(x,y\right)$

$E_1$ 椭球参量：

$\{\begin{array}{l}{L}_0=L_{01}\\ a=a_1\end{array}$ (7)

式中，$L_{01}$ 为$E_1$ 投影带中央子午线经度，$a_1$ 为$E_1$ 椭球的长半轴。

反算主变换，等角纬度$\phi$ 可直接由正算式的变形公式进行迭代计算，迭代初值可设置为$\left({\phi }_x,{\phi }_y\right)=\left(0,0\right)$ ：

$\{\begin{array}{l}{\phi }_x=\frac{1}{d_0}\cdot \left[\frac{x}{a}-{\displaystyle \sum _{j=1}^5{d}_{2\cdot j}\cdot \sin \left(2\cdot j\cdot {\phi }_x\right)}\cdot \cosh \left(2\cdot j\cdot {\phi }_y\right)\right]\\ {\phi }_y=\frac{1}{d_0}\cdot \left[\frac{y}{a}-{\displaystyle \sum _{j=1}^5{d}_{2\cdot j}\cdot \cos \left(2\cdot j\cdot {\phi }_x\right)}\cdot \sinh \left(2\cdot j\cdot {\phi }_y\right)\right]\end{array}$ (8)

其中，级数系数$d$ [20] [21]：

$d_0=1-\frac{1}{4}\cdot e^2-\frac{3}{64}\cdot e^4-\frac{5}{256}\cdot e^6-\frac{175}{16384}\cdot e^8-\frac{441}{65536}\cdot e^{10}-\frac{43659}{65536}\cdot e^{12}$

$d_2=\frac{1}{8}\cdot e^2-\frac{1}{96}\cdot e^4-\frac{9}{1024}\cdot e^6-\frac{901}{184320}\cdot e^8-\frac{16381}{5898240}\cdot e^{10}+\frac{2538673}{4587520}\cdot e^{12}$

$d_4=\frac{13}{768}\cdot e^4+\frac{17}{5120}\cdot e^6-\frac{311}{737282}\cdot e^8-\frac{18931}{20643840}\cdot e^{10}-\frac{1803171}{9175040}\cdot e^{12}$

$d_6=\frac{61}{15360}\cdot e^6+\frac{899}{430080}\cdot e^8+\frac{18757}{27525120}\cdot e^{10}+\frac{461137}{20643840}\cdot e^{12}$

$d_8=\frac{49561}{41287680}\cdot e^8+\frac{175087}{165150720}\cdot e^{10}-\frac{869251}{20643840}\cdot e^{12}$

$d_{10}=\frac{179101}{41287680}\cdot e^{10}-\frac{25387}{1290240}\cdot e^{12}$

式中，$e$ 为椭球偏心率。

为统一表达格式，参照以大地纬度B表达的高斯投影变换式，式(8)可改写为式(8a)：

$\{\begin{array}{l}{\phi }_x=\frac{1}{{\kappa }_0}\cdot \left[\frac{x}{a\cdot \left(1-e^2\right)}-{\displaystyle \sum _{j=1}^5{\kappa }_{2\cdot j}\cdot \sin \left(2\cdot j\cdot {\phi }_x\right)}\cdot \cosh \left(2\cdot j\cdot {\phi }_y\right)\right]\\ {\phi }_y=\frac{1}{{\kappa }_0}\cdot \left[\frac{y}{a\cdot \left(1-e^2\right)}-{\displaystyle \sum _{j=1}^5{\kappa }_{2\cdot j}\cdot \cos \left(2\cdot j\cdot {\phi }_x\right)}\cdot \sinh \left(2\cdot j\cdot {\phi }_y\right)\right]\end{array}$ (8a)

式(8a)，用复数表达更为简洁：

${\phi }_c=\frac{1}{{\kappa }_0}\cdot \left[\frac{z}{a\cdot \left(1-e^2\right)}-{\displaystyle \sum _{j=1}^5{\kappa }_{2\cdot j}\cdot \sin \left(2\cdot j\cdot {\phi }_c\right)}\right]$ (8a-1)

复数坐标：$\text{z}=x+i\cdot y$ ；

复数等角纬度：${\phi }_c={\phi }_x+i\cdot {\phi }_y$ 。

${\kappa }_0=1+\frac{3}{4}\cdot e^2+\frac{45}{64}\cdot e^4+\frac{175}{256}\cdot e^6+\frac{11025}{16384}\cdot e^8+\frac{43659}{65536}\cdot e^{10}$

${\kappa }_2=\frac{1}{8}\cdot e^2+\frac{11}{96}\cdot e^4+\frac{325}{3072}\cdot e^6+\frac{18599}{184320}\cdot e^8+\frac{192929}{1966080}\cdot e^{10}$

${\kappa }_4=\frac{13}{768}\cdot e^4+\frac{311}{15360}\cdot e^6+\frac{14617}{737280}\cdot e^8+\frac{26023}{1376256}\cdot e^{10}$

${\kappa }_6=\frac{61}{15360}\cdot e^6+\frac{869}{143360}\cdot e^8+\frac{5195}{768432}\cdot e^{10}$

${\kappa }_8=\frac{49561}{41287680}\cdot e^8+\frac{53333}{23592960}\cdot e^{10}$

${\kappa }_{10}=\frac{34729}{82575360}\cdot e^{10}$

2) 反算辅助转换：$\left(q,l\right)\leftarrow \left({\phi }_x,{\phi }_y\right)$

$\left(q,l\right)$ 值由下式迭代计算，迭代初值可设置为$\left(q,l\right)=\left(0,0\right)$ ：

$\{\begin{array}{l}q=\mathrm{arc}\sinh \left[\tan \left({\phi }_x\right)\cdot \cos \left(l\right)\right]\\ l=\arcsin \left[\tanh \left({\phi }_y\right)\cdot \cosh \left(q\right)\right]\end{array}$ (9)

3) 反算辅助转换：$\left(B,l\right)\leftarrow \left(q,l\right)$

大地纬度$B$ 由下式迭代计算，迭代初值可设置为$B=0$ ：

$\{\begin{array}{l}B=\arcsin \left\{\tanh \left[q+e\cdot \mathrm{arc}\tanh \left(e\cdot \sin B\right)\right]\right\}\\ L=L_0+l\end{array}$ (10)

上述各迭代计算中，一般迭代6次，角度计算精度 < 10⁻¹⁰ rad，直角坐标计算精度 < 1 mm。为确保计算精度，本文迭代次数取为10次，角度计算精度 < 10⁻¹³ rad，直角坐标计算精度 < 0.01 mm。

3.2.2. E₂面高斯正算

E₂椭球参量：$\{\begin{array}{l}{L}_0=L_{02}\\ a=a_2\\ B=B_2\\ l=L-L_0\end{array}$ (11)

1) 正算辅助转换：$\left(B,l\right)\to \left(q,l\right)$

$\{\begin{array}{l}q=\mathrm{arc}\tanh \left(\sin B\right)-e\cdot \mathrm{arc}\tanh \left(e\cdot \sin B\right)\\ l=L-L_0\end{array}$ (12)

式中，$L_0$ 为投影代中央子午线经度，$L$ 为转换点经度。

2) 正算辅助转换：$\left(q,l\right)\to \left({\phi }_x,{\phi }_y\right)$

$\{\begin{array}{l}{\phi }_x=\mathrm{arc}\tan \left[\frac{\sinh \left(q\right)}{\cos \left(l\right)}\right]\\ {\phi }_y=\mathrm{arc}\tanh \left[\frac{\sin \left(l\right)}{\cosh \left(q\right)}\right]\end{array}$ (13)

3) 正算主变换：$\left({\phi }_x,{\phi }_y\right)\to \left(x,y\right)$

$\{\begin{array}{l}x=a\cdot \left[d_0\cdot {\phi }_x+{\displaystyle \sum _{j=1}^5{d}_{2\cdot j}}\cdot \sin \left(2\cdot j\cdot {\phi }_x\right)\cdot \cosh \left(2\cdot j\cdot {\phi }_y\right)\right]\\ y=a\cdot \left[d_0\cdot {\phi }_y+{\displaystyle \sum _{j=1}^5{d}_{2\cdot j}}\cdot \cos \left(2\cdot j\cdot {\phi }_x\right)\cdot \sinh \left(2\cdot j\cdot {\phi }_y\right)\right]\end{array}$ (14)

为统一表达格式，参照以大地纬度B表达的高斯投影变换式，式(14)可改写为式(14a)：

$\{\begin{array}{l}x=a\cdot \left(1-e^2\right)\cdot \left[{\kappa }_0\cdot {\phi }_x+{\displaystyle \sum _{j=1}^5{\kappa }_{2\cdot j}}\cdot \sin \left(2\cdot j\cdot {\phi }_x\right)\cdot \cosh \left(2\cdot j\cdot {\phi }_y\right)\right]\\ y=a\cdot \left(1-e^2\right)\cdot \left[{\kappa }_0\cdot {\phi }_y+{\displaystyle \sum _{j=1}^5{\kappa }_{2\cdot j}}\cdot \cos \left(2\cdot j\cdot {\phi }_x\right)\cdot \sinh \left(2\cdot j\cdot {\phi }_y\right)\right]\end{array}$ (14a)

式(14a)，用复数表达更为简洁：

$z=x+i\cdot y=a\cdot \left(1-e^2\right)\cdot \left[{\kappa }_0\cdot {\phi }_c+{\displaystyle \sum _{j=1}^5{\kappa }_{2\cdot j}}\cdot \sin \left(2\cdot j\cdot {\phi }_c\right)\right]$ (14a-1)

3.2.3. 异面高斯换带

异面高斯换带，由$L_{01}$ 投影带的$E_1$ 面高斯反算、膨胀变换、及$L_{02}$ 投影带的$E_2$ 面高斯正算等3个计算模块构成，顺次变换就可得到异面高斯换带结果，此不赘述。

特别地，由于$B$ 与$q$ 的关系式(12)及$\left(q,l\right)$ 与$\left({\phi }_x,{\phi }_y\right)$ 的关系式(13)均为解析式，正算主变换式(14)与径差$l$ 无关。换言之，基于等角纬度的高斯投影变换精度不受投影带宽的影响，仅从转换精度的角度考虑，无需进行分带投影[19]。

4. 异面高斯投影变换的Excel VBA实现与实例

在Excel工作簿中设置了异面反算、异面正算、异面换带3个变换工作表及相应的变换宏。

(1) 椭球坐标系选择项：BJ54；XIAN80；WGS84；CGCS2000。

(2) 通用栏目：数据序号；点号编码；投影带中央经线$\left(L_{01},L_{02}\right)$ 、膨胀前后的投影面高$\left(H_{E1},H_{E2}\right)$ 。

(3) 在相关工作页内，选定椭球坐标系，在Excel表中输入源数据后，执行相应变换宏就可得到相应的转换结果，方便快捷。

4.1. 异面高斯投影换带变换实例

实例选自文献[10]。坐标系为BJ54，坐标点6组，国家标准3˚投影带中央子午线$L_{01}=111\mathrm{deg}$ ，投影面高$H_{E1}=0\text{m}$ ；独立坐标系中央子午线为$L_{02}=111.5\mathrm{deg}$ ，投影面高$H_{E2}=350\text{m}$ 。求解转换BJ54坐标为独立坐标(异面高斯换带)。y15、y25为y + 500000坐标值。

异面坐标换带，空间法与增量法换带结果相同，Excel VBA与Mathcad计算结果一致。坐标转换精度 < 0.01 mm，详见表1。本文坐标转换精度高于原文献，与原文献坐标误差 < 0.4 mm，见表2。

Table 1. Results of Gaussian transform on different faces (Spatial Method)

表1. 异面高斯换带结果表(空间法)

点号	中央经线/deg		投影面高/m		E1源坐标/m			E2换带坐标(VBA)/m			E2换带坐标(Mathcad)/m
编码	L01	L02	HE1	HE2	x1	y15	H1	x2	y25	H2	x2	y25	H2
A1	111	111.5	0	350	4389438.66500	556534.69300	0	4389485.42233	513611.93742	−349.52301	4389485.42233	513611.93742	−349.52301
A2	111	111.5	0	350	4394542.07700	565709.05200	0	4394537.72261	522814.69672	−349.52210	4394537.72261	522814.69672	−349.52210
A3	111	111.5	0	350	4386452.52800	553598.35300	0	4386515.61582	510658.96703	−349.52355	4386515.61582	510658.96703	−349.52355
A4	111	111.5	0	350	4388258.03900	556563.85400	0	4388304.63119	513634.52656	−349.52323	4388304.63119	513634.52656	−349.52323
A5	111	111.5	0	350	4386990.65600	561670.38600	0	4387008.82986	518734.00256	−349.52347	4387008.82986	518734.00256	−349.52347
A6	111	111.5	0	350	4390460.75500	563643.23500	0	4390467.93321	520726.15880	−349.52284	4390467.93321	520726.15880	−349.52284

Table 2. Gaussian strip switching errors on different surfaces

表2. 异面高斯换带误差表

点号	原文换带结果/m		E2本文换带结果/m		换带误差（本文–原文）/m
编码	x2	y25	x2	y25	$\Delta x$	$\Delta y$
A1	4389485.4220	513611.9374	4389485.42233	513611.93742	0.00033	0.00002
A2	4394537.7230	522814.6967	4394537.72261	522814.69672	−0.00039	0.00002
A3	4386515.6160	510658.9670	4386515.61582	510658.96703	−0.00018	0.00003
A4	4388304.6310	513634.5266	4388304.63119	513634.52656	0.00019	−0.00004
A5	4387008.8300	518734.0026	4387008.82986	518734.00256	−0.00014	−0.00004
A6	4390467.9330	520726.1588	4390467.93321	520726.15880	0.00021	0.00000

4.2. 椭球膨胀变换结果比较

大地坐标增量法，计算简单；空间直角坐标法计算稍繁，但逻辑清晰明了。

两种膨胀变换方法，角度计算精度 < 10⁻¹² deg，详见表3、表4。

Table 3. Results of geodetic coordinate transformation on different surfaces (Spatial Method)

表3. 异面大地坐标变换结果(空间法)

点号	E1中间结果			空间坐标			E2中间结果
编码	B1/deg	L1/deg	H1/m	X/m	Y/m	Z/m	B2/deg	L2/deg	H2/m
A1	39.636360862454	111.658552570291	0	−1815346.34964	4571387.05577	4047045.81403	39.636371242115	111.658552570291	−349.52301
A2	39.681666024432	111.765920250404	0	−1822719.10963	4564995.71637	4050918.37015	39.681676407135	111.765920250404	−349.52210

续表

A3	39.609656301174	111.624108514755	0	−1813294.69308	4574235.34300	4044762.00971	39.609666679028	111.624108514755	−349.52355
A4	39.625726194086	111.658791375557	0	−1815643.37008	4572079.45519	4046136.42949	39.625736573028	111.658791375557	−349.52323
A5	39.613959753352	111.718144259821	0	−1820686.90501	4570970.12638	4045130.10467	39.613970131498	111.718144259821	−349.52347
A6	39.645066723023	111.741449593920	0	−1821730.02327	4568182.97988	4047790.16084	39.645077103270	111.741449593920	−349.52284

Table 4. Results of geodetic coordinate transformation on different surfaces (Incremental Method)

表4. 异面大地坐标变换结果(增量法)

点号	E1中间结果			膨胀增量			E2中间结果
编码	B1/deg	L1/deg	H1/m	dB/deg	dL/deg	dH/m	B2/deg	L2/deg	H2/m
A1	39.636360862454	111.658552570291	0	1.037965974E−05	0.000000000E+00	−349.52301	39.636371242114	111.658552570291	−349.52301
A2	39.681666024432	111.765920250404	0	1.038270233E−05	0.000000000E+00	−349.52210	39.681676407134	111.765920250404	−349.52210
A3	39.609656301174	111.624108514755	0	1.037785414E−05	0.000000000E+00	−349.52355	39.609666679028	111.624108514755	−349.52355
A4	39.625726194086	111.658791375557	0	1.037894177E−05	0.000000000E+00	−349.52323	39.625736573027	111.658791375557	−349.52323
A5	39.613959753352	111.718144259821	0	1.037814572E−05	0.000000000E+00	−349.52347	39.613970131498	111.718144259821	−349.52347
A6	39.645066723023	111.741449593920	0	1.038024642E−05	0.000000000E+00	−349.52284	39.645077103269	111.741449593920	−349.52284

4.3. 椭球膨胀偏位分析

E₁→E₂椭球膨胀异面换带坐标与E₁→E₁椭球同面换带坐标相比，存在坐标偏位，见表5。

Table 5. Ellipsoid expansion transformation deviation table

表5. 椭球膨胀变换偏位表

点号	E1→E2异面换带坐标			E1→E1同面换带坐标			椭球膨胀偏位 E2异面–E1同面
编码	x2/m	y25/m	H2/m	x2’/m	y25’/m	H2’/m	Dx/m	Dy/m
A1	4389485.42213	513611.93742	−349.52301	4389243.41415	513611.19255	0.00000	242.00798	0.74487
A2	4394537.72241	522814.69672	−349.52210	4394295.43686	522813.44827	0.00000	242.28556	1.24845
A3	4386515.61562	510658.96703	−349.52355	4386273.77080	510658.38375	0.00000	241.84482	0.58328
A4	4388304.63100	513634.52656	−349.52323	4388062.68789	513633.78045	0.00000	241.94311	0.74610
A5	4387008.82967	518734.00256	−349.52347	4386766.95775	518732.97741	0.00000	241.87192	1.02516
A6	4390467.93302	520726.15880	−349.52284	4390225.87106	520725.02463	0.00000	242.06196	1.13417

由表5可见，$Dx$ 偏位大，$Dy$ 偏位小。投影抵偿面仅350 m，$Dx$ 偏位就已达到近242 m。为此，有必要对偏位值进行进一步分析。

不失一般性，将椭球近似为圆球，则有：

$\{\begin{array}{l}x=B\cdot a\\ dx=B\cdot da+a\cdot dB\\ dx1=B\cdot da\\ dx2=a\cdot dB\end{array}$ (14)

$\{\begin{array}{c}y=l\cdot r=l\cdot a\cdot \cos B\\ dy=l\cdot \cos B\cdot da-l\cdot a\cdot \sin B\cdot dB\\ dy1=l\cdot \cos B\cdot da\\ dy2=-l\cdot a\cdot \sin B\cdot dB\end{array}$ (15)

由式(14)、(15)，得到圆球膨胀变换偏位表，见表6。

Table 6. Deviation table of spherical expansion transformation

表6. 圆球膨胀变换偏位表

点号	E1椭球大地坐标/deg		膨胀变换/deg		圆球膨胀偏位/m
编码	B1	L1	dB	L = L2-L02	dx1	dx2	dx	dy1	dy2	dy
A1	39.636360860714	111.658552570275	1.037965973982E−05	0.158552570	242.12475	0.00363	242.12838	0.74588	−0.11687	0.62901
A2	39.681666022692	111.765920250386	1.038270233318E−05	0.265920250	242.40150	0.00363	242.40514	1.25016	−0.19625	1.05390
A3	39.609656299434	111.624108514740	1.037785413764E−05	0.124108515	241.96162	0.00363	241.96525	0.58407	−0.09141	0.49266
A4	39.625726192346	111.658791375541	1.037894177000E−05	0.158791376	242.05979	0.00363	242.06342	0.74712	−0.11701	0.63011
A5	39.613959751613	111.718144259804	1.037814572195E−05	0.218144260	241.98791	0.00363	241.99154	1.02655	−0.16069	0.86586
A6	39.645066721283	111.741449593902	1.038024642484E−05	0.241449594	242.17793	0.00363	242.18156	1.13571	−0.17801	0.95770

对比表5、表6，椭球膨胀变换偏位$\left(Dx,Dy\right)$ 与圆球膨胀变换偏位$\left(dx,dy\right)$ 基本接近，从而证实了异面高斯投影变换的正确性。

从式(14)、(15)可见，圆球膨胀偏位$\left(dx,dy\right)$ 主要由椭球长半轴增量$da$ 所贡献，偏位均与$da$ 呈线性增长，$dx$ 偏位大，其增长系数为$B$ ；$dy$ 偏位小，其比例系数为$l\cdot \cos B$ 。

5. 结语

1) 参照基于大地纬度B表达的高斯投影变换式，本文给出了同一格式基于等角纬度φ表达的高斯投影变换式及其幂级数系数。

2) 异面高斯投影变换，正反算变换均利用正算式，反算采用正算式迭代计算。Excel VBA变换宏，可一键实现坐标的批量变换，具有一定的工程应用价值。

3) 椭球膨胀变换的大地坐标增量法及空间直角坐标法变换结果完全相同。

4) 膨胀偏位分析表明，椭球膨胀偏位$\left(Dx,Dy\right)$ ，主要由椭球长半轴增量$da$ 所贡献，且与$da$ 呈线性增长。$dx$ 偏位大，其增长系数为$B$ ；$dy$ 偏位小，其比例系数为$l\cdot \cos B$ 。

5) 基于等角纬度的高斯投影变换，投影带不限定于传统的7°范围。仅从转换精度的角度考虑，无需进行分带投影。

致 谢

海军工程大学边少锋教授对课题研究给予了大力支持，在此表示衷心感谢！

基金项目

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

NOTES

^*第一作者。

^#通讯作者。