割圆锥投影的正反算变换——Excel VBA实数计算与Mathcad复数变换

Forward and Inverse Transformation of Secant Conic Projection—Excel VBA Real Number Calculation and Mathcad Complex Number Transformation

摘要: 圆锥投影图广泛用于中小比例尺地图，如国际航空图、全国地图、中小比例尺分省地图、地质图、专用地图等。利用计算工具Excel VBA及Mathcad分别开发了适用于BJ54、XIAN80、WGS84、CGCS2000坐标系割圆锥投影的实数及复数正反算变换程序，为野外地质调查进行快捷的批量坐标转换提供了普适工具。

Abstract: Conic projection maps are widely used in small and medium-sized maps, such as international aerial maps, national maps, provincial maps at small and medium-sized scales, geological maps, and specialized maps. Real and complex forward and inverse transformation programs suitable for cutting cone projection in BJ54, XIAN80, WGS84, and CGCS2000 coordinate systems were developed using calculation tools Excel VBA and Mathcad, respectively, providing a universal tool for fast batch coordinate transformation in field geological surveys.

1. 前言

圆锥投影图，因其制图简单、等角变换、在中纬度地区变形小等特点，从而广泛用于中小比例尺地图，如国际航空图、全国地图、中小比例尺的分省地图、地质图、专用地图等。

随着我国工程建设的发展，对投影变形的控制要求较高，尤其是地铁及高铁，其长度变形精度要求小于1 m/100km。采用高斯投影时，为满足变形控制要求，投影分带需<57 km，长线路工程的坐标转换非常繁琐。由于割圆锥投影的南北基准纬线长度不变形，投影区域变形小，因此，割圆锥投影变换对东西向条带状的投影变换有其独特优势。

国家行业标准地质图编绘规范规定，1:50万、1:100万地质图，需采用割圆锥投影[1]-[3]。

我国圆锥投影变换研究始于上世纪60年代。杨启和(1965, 1982, 1983)系统研究了等角投影变换[4]-[6]，讨论了圆锥投影的一般特性；黄伟(1979)介绍了国际上1:100万地图及我国1978年前后1:100万地图采用圆锥投影编制应用的基本情况[7]，肖学勤(1985)研究了正轴双标纬线等角圆锥投影[8]，开发了Fortran 77及Casiofx-180F正算程序，用于1:50万、1:100万、1:200万的土地资源及草地资源调查；庄培德(1995)开发了PC-1500袖珍机的等角割圆锥投影正反解程序[9]；钟业勋、李占元(2002)针对广西的地理位置及分布形状，研究了适用于广西壮族自治区的正轴等角割圆锥投影，计算了XIAN80坐标系的投影参数、长度变形及面积变形[10]；韩尤杰、尹继法等(2007)基于VC++开发了圆锥投影的图形界面的正反解程序(1940克拉索夫斯基椭球、1967国际椭球、1975国际椭球)，实现了(等角、等距离、等面积的)地理坐标与平面坐标之间的转换[11]；戴勤奋(2006)、达恒(中国)有限公司(2007)编制了图形界面的等角投影正反变换程序[12] [13]；李厚朴、边少锋等(2012, 2015, 2024)研究了复数域的等角投影及等角圆锥投影的参数优化[14]-[16]；赵淑湘(2016)研究了甘肃省的割圆锥投影(CGCS2000)，计算了投影参数，给出了161组坐标变换值[17]；丁士俊、李鹏鹏等(2022)研究了等角圆锥投影反解的不同算法，利用Mathematica导出了等量纬度反解大地纬度的直接算法[18]；焦晨晨、李松林等(2023)研究了等角圆锥投影基准纬线，利用Mathenatica导出了计算基准纬线的非迭代算法[19] [20]；割圆锥投影的1:50万广东省地质图于1999年出版。

本文在前人研究的基础上，从方便快捷、批量转换的实用角度出发，利用Excel VBA及Mathcad计算工具，开发了割圆锥投影正反算变换程序，适用于BJ54、XIAN80、WGS84、CGCS2000坐标系的坐标变换。

2. 割圆锥投影

圆锥投影，又称正轴等角圆锥投影，适用于中小比例尺地质图，如1:50万、1:100万地质图。

切圆锥投影，是用一个圆锥面正轴切于球面，将圆锥面作为投影面，然后沿圆锥面的一条母线切开展成平面的一种投影[21]-[23]，切纬线B₀长度不变形，见图1(a)。

割圆锥投影，又称兰勃特投影，是德国数学家J. H. Lembert于1772年拟定。该投影是用一个正圆锥割于球面两条标准纬线B₁及B₂ [21]-[23]，割纬线B₁、B₂长度不变形，见图1(b)。

投影坐标系，是以中央经线L₀为x轴，以基准纬线B_s为y轴，x轴与y轴的交点为投影原点O (0,0)，见图2(a)；球面与平面的投影线段微分关系见图2(b)。

(a) 切圆锥 (b) 割圆锥

Figure 1. Schematic diagram of isometric cone projection

图1. 等角圆锥投影示意图

(a) 坐标系关系图 (b) 投影线段微分关系图

Figure 2. Relationship between conic projection ellipsoid and plane

图2. 圆锥投影椭球面与平面关系图

3. 割圆锥投影的实数变换

3.1. 实数正算

正算，即由大地坐标(B, L)转换为直角坐标(x, y)。

数学上，圆锥投影平面直角坐标(x, y)的计算通式为[21]-[23]：

$\{\begin{array}{l}x={\rho }_s-\rho \cdot \cos \delta \\ y=\rho \cdot \sin \delta \end{array}$ (1)

其中：

$\{\begin{array}{l}\rho =C\cdot \exp \left(-{\alpha }_c\cdot q\right)\\ \delta ={\alpha }_c\cdot l\end{array}$ (2)

$\{\begin{array}{l}m=\frac{A^{\prime }{D}^{\prime }}{AD}=-\frac{1}{M}\cdot \frac{d\rho }{dB}\\ n=\frac{A^{\prime }{E}^{\prime }}{AE}=-\frac{\rho }{r}\cdot \\ {\alpha }_c=\frac{d\delta }{dl}\end{array}$ (3)

$\{\begin{array}{l}\delta ={\displaystyle {\int }_{L_0}^L{\alpha }_c\cdot dl}={\alpha }_c\cdot l\\ l=L-L_0\end{array}$ (4)

$\{\begin{array}{l}e=\sqrt{1-{\left(1-f\right)}^2}\\ N=\frac{a}{\sqrt{1-{\left(e\cdot \sin B\right)}^2}}\\ r=N\cdot \cos B\\ U=\tan \left(\frac{\pi }{4}+\frac{B}{2}\right)\cdot {\left(\frac{1-e\cdot \sin B}{1+e\cdot \sin B}\right)}^{\frac{e}{2}}\end{array}$ (5)

a	为椭球长半轴
f	为椭球扁率
e	为椭球偏心率
M	子午圈曲率半径
N	卯酉圈曲率半径
r	纬圈半径
U	等角表象函数
L0	中央经线经度
L	经线经度
l	椭球面经线与中央经线夹角
δ	平面经线与中央经线夹角
ρ	平面经线长
ρs	中央经线顶点到基准纬线的径线长
m	经线变形系数
n	纬线变形系数

对于割圆锥投影，B₁、B₂两条纬线长度不变形，即：

$\{\begin{array}{l}{n}_1=1\\ n_2=1\end{array}$ (6)

从而有：

$\{\begin{array}{l}{\alpha }_c=\frac{d\delta }{dl}=-\frac{\ln r{}_2-\ln r_1}{\ln U_2^{{\alpha }_c}-\ln U_1^{{\alpha }_c}}\\ C=\frac{r_1\cdot U_1^{{\alpha }_c}}{{\alpha }_c}=\frac{r_2\cdot U_2^{{\alpha }_c}}{{\alpha }_c}\\ \rho =\frac{C}{U^{{\alpha }_c}}\end{array}$ (7)

式中：

C	赤道投影半径
αc	平面角与球面角的投影变换常数

确定α_c及C后，就可由式(1)将坐标(B, L)转换为直角坐标(x, y)。

3.2. 实数反算

反算，即由直角坐标(x, y)转换为大地坐标(B, L)。

由式(7)确定α_c、C、ρ_s参量后，就可由等量纬度q确定大地坐标(B, L) [14] [15] [18]：

$\{\begin{array}{l}q=\frac{1}{{\alpha }_c}\cdot \frac{\ln \left(C\right)}{\ln \left(\sqrt{{\left({\rho }_s-x\right)}^2+y^2}\right)}\\ l=\frac{1}{{\alpha }_c}\cdot \arctan \left(\frac{y}{{\rho }_s-x}\right)\\ \phi =\arcsin \left(\tanh \left(q\right)\right)\\ B=\phi +{\displaystyle \sum _{t=1}^5{a}_t\cdot \sin \left(2\cdot t\cdot \phi \right)}\\ L=L_0+l\end{array}$ (8)

等角变换系数a可由下式确定(丁士俊，2022) [18]：

$\{\begin{array}{l}{a}_2=2\cdot e_3-\frac{2}{3}\cdot e_3{}^2-2\cdot e_3{}^3+\frac{116}{45}\cdot e_3{}^4+\frac{26}{45}\cdot e_3{}^5\\ a_4=\frac{7}{3}\cdot e_3{}^2-\frac{8}{5}\cdot e_3{}^3-\frac{227}{45}\cdot e_3{}^4+\frac{2704}{315}\cdot e_3{}^5\\ a_6=\frac{56}{15}\cdot e_3{}^3-\frac{136}{35}\cdot e_3{}^4-\frac{1262}{105}\cdot e_3{}^5\\ a_8=\frac{4279}{630}\cdot e_3{}^4-\frac{332}{35}\cdot e_3{}^5\\ a_{10}=\frac{4174}{315}\cdot e_3{}^5\end{array}$ (9)

其中第3偏心率：

$e_3=\frac{1-\sqrt{1-e^2}}{1+\sqrt{1-e^2}}$ (10)

4. 割圆锥投影的复数变换

4.1. 复数正算

$\{\begin{array}{l}q=\mathrm{arctanh}\left(\sin B-e\cdot \mathrm{arctanh}\left(e\cdot \sin B\right)\right)\\ l=L-L_0\\ w=q+i\cdot l\\ z={\rho }_s-C\cdot \exp \left(-{\alpha }_c\cdot w\right)\end{array}$ (11)

4.2. 复数反算

$\{\begin{array}{l}w=\frac{1}{{\alpha }_c}\cdot \frac{\ln \left(C\right)}{\ln \left({\rho }_s-z\right)}\\ q=\mathrm{Re}\left(w\right)\\ l=\mathrm{Im}\left(w\right)\\ \phi =\arcsin \left(tanh\left(q\right)\right)\\ B=\phi +{\displaystyle \sum _{t=1}^5{a}_{2t}\cdot \sin \left(2t\cdot \phi \right)}\\ L=L_0+l\end{array}$ (12)

5. 割圆锥投影变换程序开发

5.1. Excel VBA实域开发

Excel目前不支持复数运算，Excel平台的割圆锥投影变换只能在实数域中开发。

割圆锥投影正算变换宏为scp( )，反算变换宏为Inv_scp( )。

坐标转换程序使用方法：

1) 人工选择坐标系(三角下拉菜单)：BJ54、XIAN80、WGS84、CGCS2000；

2) 人工输入投影参数：L₀、B_s、B₁、B₂；

3) 人工输入源坐标：(B, L)或(x, y)；

4) 人工输入坐标点序号；

5) 运行宏，一键得到目标坐标。

5.2. Mathcad复域开发

Mathcad割圆锥投影复域变换的核心代码非常简洁：

1) 正算核心代码：

$\begin{array}{l}q:=\mathrm{atanh}\left(\sin \left(B\right)\right)-e\cdot \mathrm{atanh}\left(e\cdot \sin \left(B\right)\right)\\ l:=L-L_0\\ w:=q+l\cdot 1i\\ z:={\rho }_s-C\cdot \exp \left(-{\alpha }_c\cdot w\right)\end{array}$

2) 反算核心代码：

$w:=\frac{1}{{\alpha }_c}\cdot \frac{\ln \left(C\right)}{\ln \left({\rho }_s-z\right)}$

$\begin{array}{l}q:=\mathrm{Re}\left(w\right)\\ l:=\mathrm{Im}\left(w\right)\\ \phi :=\mathrm{asin}\left(tanh\left(q\right)\right)\\ B:=\phi +{\displaystyle \sum _{t=1}^5{a}_{2t}\cdot \sin \left(2t\cdot \phi \right)}\\ L:=L_0+l\end{array}$

6. 计算实例

文献[17]刊载了161组(B, L)源数据的CGCS2000坐标系的正算结果。割圆锥投影中央经度L₀ =101˚，基准纬度B_s = 32˚，第1标准纬线B₁ = 34˚，第2标准纬线B₂ = 41˚。

本文选取其中8组(B, L)源数据，用Excel VBA及Mathcad 15进行了投影变换验算，两者计算结果一致。表1为正算结果，本文换算精度高于原文献；表2为反算结果，可精确还原到原文献经纬度坐标。

Table 1. Conical projection calculation table

表1. 割圆锥投影正算表

坐标点	源坐标		原文献正算坐标		本文正算坐标
序号	B/˚	L/˚	x/m	y/m	x/m	y/m
1	43	101	1,220,470	0	1220468.242	0.000
2	43	110	1,255,650	734,810	1255650.125	734810.617
3	40	101	887,040	0	887043.311	0.000
4	40	110	923,750	766,670	923750.373	766665.594
5	35	101	332,960	0	332963.006	0.000
6	35	110	372,200	819,600	372204.587	819601.692
7	32	101	0	0	0.000	0.000
8	32	110	40,760	851,410	40764.648	851412.537

Table 2. Reverse calculation table for conical projection

表2. 割圆锥投影反算表

坐标点	源坐标		反算坐标		原文献坐标
序号	x/m	y/m	B/˚	L/˚	B/˚	L/˚
1	1220468.242	0.000	43.000000	101.000000	43	101
2	1255650.125	734810.617	43.000000	110.000000	43	110
3	887043.311	0.000	40.000000	101.000000	40	101
4	923750.373	766665.594	40.000000	110.000000	40	110
5	332963.006	0.000	35.000000	101.000000	35	101
6	372204.587	819601.692	35.000000	110.000000	35	110
7	0.000	0.000	32.000000	101.000000	32	101
8	40764.648	851412.537	32.000000	110.000000	32	110

7. 结语

1) 割圆锥投影的南北标准纬线长度不变形，投影区域变形小。因此，割圆锥投影变换对东西向条带状的投影变换有其独特优势，尤其适用于东西向长线路的地铁及高铁工程的地图投影。

2) 本文在割圆锥投影变换中，等角变换系数采用了丁士俊(2022)以第3偏心率表达的最新研究成果，其系数比以用第1偏心率表达的公式简洁。

3) 本文利用割圆锥投影理论，分别在实域的Excel VBA平台及复域的Mathcad平台开发了割圆锥投影的正反算变换工具。利用公开发表的相关论文数据进行了计算验证对比，本文转换工具计算精度更高，结果更准确可靠。

4) Mathcad是一款所见即所得的计算软件，其代码表达式与数学表达式十分接近。相比于实数域，复数域的Mathcad割圆锥变换工具，其代码简洁，可读性好。

5) 本文Excel VBA变换工具，虽未进行图形界面设计，但使用十分简单，可实现坐标的一键批量转换。在Excel表中输入相应的原始数据，直接运行正变换或反变换宏，就可得到目标坐标，简洁明快。

6) Excel变换表可直接作为工作用表。作为一款批量坐标转换的高效普适工具，具有较高的使用价值，特别适用于野外工作者进行坐标的快捷转换。

致 谢

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

基金项目

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

NOTES

^*第一作者。

^#通讯作者。