1. 前言
高斯投影变换是一经典课题。高斯投影属于等角投影[1]-[3],根据复变函数理论,复变量
与
之间的变换为[4]-[7]:
复变量:
即正反变换为:
在实数域,高斯投影正解(中央子午线弧长
)是一关于大地纬度
的积分式[4]-[7]:
边少锋等(2001, 2004)将实数变量
拓展到复数变量
,从而将高斯投影变换从实数域拓展到复数域,得到关于复数变量
的积分[8]-[12]:
上述积分式目前还没有找到闭合解。解决的方法是对被积分函数展开为三角函数级数式,然后积分得到三角函数级数近似解。
传统的高斯投影变换,以大地纬度
作为积分变量。李厚朴、边少锋等(2007, 2009, 2012, 2015~2017)研究了基于复变量等角纬度
的非迭代的正反解表达式[13]-[18]。
金立新等(2017),在基于复变量等角纬度
高斯投影变换的基础上,将其解析到实数域,给出了实数域非迭代的正反解表达式及其变换系数[19]。
实际上,非迭代的高斯投影变换,正解与反解为2套不同的级数表达式与系数,表达冗长。
本文在上述研究的基础上,从数值计算的实用角度出发,对实数域基于等角纬度
的高斯投影变换,正反解变换均采用正解公式,反解变换利用正解公式迭代计算。算法简明,程序精短。
2. 等角纬度高斯投影变换路径
等角纬度的高斯投影变换,其变换路径为:
复数变换路径:
1) 正算路径:
2) 反算路径:
a. 迭代法:
b. 非迭代法:
实数变换路径:
1) 正算路径:
2) 反算路径:
a. 迭代法:
b. 非迭代法:
上述路径中,方框内的变换为主变换,其它为辅助变换。
其中:
为实数纬度及径差
为复数等量纬度
为复数等角纬度
为复数底点纬度(复数等距离纬度)
为复数直角坐标
3. 复变量等角纬度的高斯投影变换
3.1. 复变量正算变换
复变量的正算变换,由以下步骤完成:
1) 正算辅助转换
(1)
2) 正算辅助转换
(2)
3) 正算主变换
李厚朴等(2009, 2012, 2018, 2021)给出了基于等角纬度复变量的主变换式[12] [14]-[16]:
(3)
其中系数
:
式中,
为椭球偏心率。
3.2. 复变量反算变换(非迭代变换)
复变量的反算变换的迭代法,可顺序由式(3)、(2)、式(1)迭代求解,此不赘述。
对于非迭代法,李厚朴等(2009, 2012, 2018, 2021)给出的复变量的反算变换式为[12] [14]-[16]:
1) 反算主变换
(4)
2 )反算辅助转换
(5)
3) 反算辅助转换
(6)
式中,辅助变换由
反算
的级数系数
为:
式中,辅助变换由
反算
的级数系数
为:
4. 实变量等角纬度的高斯投影变换
基于复变量等角纬度的高斯投影变换,计算精简,宜为首选方法。
目前,大众计算工具Excel还不支持复数计算,因此,基于实数等角纬度的高斯投影变换还有很大的实用价值。
4.1. 实变量等角纬度与等量纬度关系式
在实数域中,等角纬度与等量纬度的关系是等角纬度高斯投影变换的桥梁。
等角变换中,有以下复数关系式:
(7)
以式(7)为基础,可以分别得到式(8)及(9):
(8)
(9)
由(8-1)与(9-1)的实部与虚部分别相等,得到:
(10)
由(8-2)与(9-2)的实部与虚部分别相等,得到:
(11)
则,由式(10)及(11)即可得到
与
的转换解析式[2] [19] [20]:
(12)
式(12)是等角纬度
与等量纬度
的转换关系式。
4.2. 实变量正算变换(解析式)
1) 正算辅助转换
(13)
式中,
为投影带中央子午线经度,
2) 正算辅助转换
由关系式(12)得到:
(14)
3) 正算主变换
由于:
将上式代入式(3),得到实数域等角纬度的高斯投影正算主变换式:
(15)
4.3. 实变量反算变换
4.3.1. 实变量反算迭代法
迭代算法,反算变换只需利用正算式及其变换系数,计算简明。
1) 反算主变换
由正算主变换式(15)可得到反算主变换式(16),迭代初值可设置为
:
(16)
2) 反算辅助转换
由关系式(12),可得到方程组(17),迭代初值可设置为
:
(17)
3) 反算辅助转换
由关系式(13),可得到式(18),迭代初值可设置为
:
(18)
上述各迭代计算中,一般迭代5~6次,角度计算精度 < 10−10 rad,直角坐标计算精度 < 1 mm,满足高斯投影变换精度要求。为保证更高计算精度及简化程序起见,本文迭代次数取固定值10次。
4.3.2. 实变量反算非迭代法
基于等角纬度的非迭代反算变换,转换式与正算完全不同、系数冗长[19]。
1) 反算主变换
(19)
2) 反算辅助变换
(20)
3) 反算辅助变换
(21)
5. 高斯变换与投影分带
传统的基于大地纬度的高斯投影变换式,展开为Taylor幂级数,初始展开点位于中央子午线,变换函数
表达为径差
的幂级数,其变换精度与径差
的大小有关,
越大精度越差,需要的级数项数越多。传统的基于大地纬度的高斯投影Taylor幂级数变换式,投影带宽限制在 < 7˚范围内,直角坐标的变换精度 < 1 mm [7]。
基于等角纬度的高斯投影变换,变换式
表达为等角纬度的三角函数,
与
的关系式(1),以及
与
的关系式(2)均为解析式;正算变换式(3)与径差
无关。换言之,基于等角纬度的高斯投影变换,变换精度与径差
无关,投影带宽不限于7˚范围。仅就考虑变换精度而言,无需进行分带投影。由此,对于传统的3˚及6˚带的分带投影,可以进行跨带坐标变换,不限于传统相邻带坐标变换。
6. 变换实例
[实例1] 选自文献([5], p. 168)实例,点名1。坐标系为BJ54坐标系,投影带中央子午线
,变换点大地坐标为:
,
。
解算:1) 由大地坐标
正算直角坐标
;
2) 由直角坐标
反算大地坐标
。
正算及迭代法反算结果详见表1。正算结果与原文献相同,本文计算精度高于原文献;反算结果
可精确还原到原大地坐标值
。
Table 1. Gaussian projection transform of real variables with equal angular latitude
表1. 实变量等角纬度高斯投影变换表
正算 |
参量 |
原文正算 |
本文正算 |
本文反算 |
反算 |
序次 |
中央子午线 |
L0/deg |
111 |
111 |
111 |
序次 |
0 |
纬度 |
B |
31˚04'41.6832" |
31˚04'41.6832" |
31˚04'41.6832" |
3 |
|
经度 |
L |
111˚47'24.8974" |
111˚47'24.8974" |
111˚47'24.8974" |
|
1 |
等量纬度 |
q/rad |
|
0.567699332168 |
0.567699332168 |
2 |
|
经差 |
l/rad |
|
0.013792451809 |
0.013792451809 |
|
2 |
等角纬度 |
/rad |
|
0.539495977396 |
0.539495977396 |
1 |
|
|
/rad |
|
0.011833952487 |
0.011833952487 |
|
3 |
直角坐标 |
x/m |
3439978.971 |
3439978.9701 |
3439978.9701 |
0 |
|
|
y/m |
757412.873 |
75412.8724 |
75412.8724 |
|
|
|
y5/m |
|
575412.8724 |
575412.8724 |
|
[实例2] 选自文献([16], p. 139)例10.6。坐标系为CGCS2000坐标系,变换点大地坐标为:
a)
,径差
,中央子午线
;
b)
,径差
,中央子午线
;
求上述大地坐标的直角坐标,以及由直角坐标反算溯源大地坐标。
不分带的正算结果详见表2,反算结果见表3。溯源反算的大地坐标与源坐标完全一致。实例2b的径差
,跨越了标准投影分带。
7. 结语
1) 传统的高斯投影变换采用大地纬度作为变换参数。本文利用李厚朴、边少锋、金立新等基于等角纬度作为变换参数的高斯投影变换的最新研究成果,在Excel VBA上实现了实数域的高斯投影变换。
2) 实数域等角纬度的高斯投影主变换式与复数域主变换式类似,同为三角函数级数表达式。
3) 根据复数域等角纬度
与等量纬度
的关系,导出了实数域
与
的解析关系式。这一关系式是组构实数域等角纬度的高斯投影正反算变换的桥梁。
Table 2. Gaussian projection non band transform calculation table
表2. 高斯投影不分带变换正算表
项目 |
参数名 |
参数符号 |
单位 |
正算变换值 |
原文献值 |
实例2a |
纬度 |
B |
˚ |
45 |
例10.6 |
|
经度 |
L |
˚ |
3 |
|
|
中央子午线 |
L0 |
˚ |
0 |
|
|
径差 |
l |
˚ |
3 |
|
正算直角坐标 |
直角坐标 |
x |
m |
4989325.234673 |
4.989325235*1E6 |
|
直角坐标 |
y |
m |
236540.642360 |
236540.6424 |
实例2b |
纬度 |
B |
˚ |
45 |
|
|
经度 |
L |
˚ |
3 |
|
|
中央子午线 |
L0 |
˚ |
15 |
|
|
径差 |
l |
˚ |
−12 |
|
正算直角坐标 |
直角坐标 |
x |
m |
5055522.235134 |
|
|
直角坐标 |
y |
m |
−946127.113917 |
|
Table 3. Gaussian projection non band transform inverse calculation table
表3. 高斯投影不分带变换反算表
项目 |
参数名 |
参数符号 |
单位 |
反算变换值 |
原文献值 |
实例2a |
直角坐标 |
x |
m |
4989325.234673 |
例10.6 |
|
直角坐标 |
y |
m |
236540.642360 |
|
|
中央子午线 |
L0 |
° |
0 |
|
|
径差 |
l |
° |
3 |
|
反算大地坐标 |
纬度 |
B |
° |
45 |
45 |
|
经度 |
L |
° |
3 |
3 |
实例2b |
直角坐标 |
x |
m |
5055522.235134 |
|
|
直角坐标 |
y |
m |
−946127.113917 |
|
|
中央子午线 |
L0 |
° |
15 |
|
|
径差 |
l |
° |
−12 |
|
反算大地坐标 |
纬度 |
B |
° |
45 |
|
|
经度 |
L |
° |
3 |
|
4) 基于等角纬度的高斯投影变换,非迭代法的正算及反算为2套各不相同的级数式及级数系数,公式及系数表达冗长。本文从数值计算的实用角度出发,反算变换利用正算式迭代计算。算法简明,程序精短。迭代5~6次,角度计算精度 < 10−10 rad,直角坐标变换精度 < 1 mm。迭代算法对计算效率的影响可忽略不计。
5) 基于等角纬度的高斯投影变换,变换精度与径差
无关,投影带宽不限于7˚范围。仅就考虑变换精度而言,无需进行分带投影。对于传统的3˚及6˚带的分带投影,可以进行跨带坐标变换,不限于传统相邻带坐标换带。
致 谢
海军工程大学边少锋教授对课题研究给予了大力支持,在此表示衷心感谢!
基金项目
国家自然科学基金项目(42074010)。
NOTES
*第一作者。
#通讯作者。