基于实数等角纬度的高斯投影变换
Gaussian Projection Transformation Based on Real Equiangular Latitude
DOI: 10.12677/gst.2024.124043, PDF, HTML, XML,    国家自然科学基金支持
作者: 刘大海*, 方春波#:深圳市地质局,广东 深圳;深圳地质建设工程公司,广东 深圳;陈永红, 洪声亮:深圳地质建设工程公司,广东 深圳;陈永冰:海军工程大学电气工程学院,湖北 武汉
关键词: 高斯投影正变换反变换等角纬度复数实数Excel VBAGaussian Projection Forward Transformation Inverse Transformation Isometric Latitude Complex Number Real Number Excel VBA
摘要: 传统的高斯投影变换采用大地纬度作为变换参数。李厚朴、边少锋等研究了基于等角纬度的高斯投影复变换,金立新给出了实数域的非迭代法的反算变换式。基于非迭代法的高斯投影变换,正反算变换为2套各不相同的变换式及其变换系数,表达冗长。本文从数值计算的实用角度出发,正反算变换均采用正算式,反算变换利用正算式迭代计算。算法简明,程序精短,易于在Excel VBA上实现。基于等角纬度的高斯投影变换,变换精度与径差l无关。仅就考虑变换精度而言,无需进行分带投影。
Abstract: The traditional Gaussian projection transformation uses geodetic latitude as the transformation parameter. Li Houpu, Bian Shaofeng, and others studied Gaussian projection complex transformations based on equiangular latitude, and Jin Lixin provided the inverse transformation formula for non-iterative methods in the real number field. The Gaussian projection transformation based on non-iterative method involves two sets of different transformation equations and their transformation coefficients, which are expressed in a lengthy manner. This article starts from the practical perspective of numerical calculation, and both forward and inverse transformations are calculated using forward equations, while inverse transformations are iteratively calculated using forward equations. The algorithm is concise, the program is short, and it is easy to implement on Excel VBA. The Gaussian projection transformation based on equiangular latitude is independent of the transformation accuracy and radial difference l. In terms of transformation accuracy alone, there is no need for band projection.
文章引用:刘大海, 方春波, 陈永红, 洪声亮, 陈永冰. 基于实数等角纬度的高斯投影变换[J]. 测绘科学技术, 2024, 12(4): 349-358. https://doi.org/10.12677/gst.2024.124043

1. 前言

高斯投影变换是一经典课题。高斯投影属于等角投影[1]-[3],根据复变函数理论,复变量 z w 之间的变换为[4]-[7]

复变量: { z=x+iy w=q+il

即正反变换为: { z=f( w )=x( q,l )+iy( q,l ) w= f 1 ( z )=q( x,y )+il( x,y )

在实数域,高斯投影正解(中央子午线弧长 X )是一关于大地纬度 B 的积分式[4]-[7]

X= 0 B MdB =a( 1 e 2 ) 0 B 1 1 e 2 sin ( B ) 2 3 dB

边少锋等(2001, 2004)将实数变量 B 拓展到复数变量 B c ,从而将高斯投影变换从实数域拓展到复数域,得到关于复数变量 B c 的积分[8]-[12]

z=x+iy= 0 B c MdB =a( 1 e 2 ) 0 B c 1 1 e 2 sin ( B ) 2 3 dB

上述积分式目前还没有找到闭合解。解决的方法是对被积分函数展开为三角函数级数式,然后积分得到三角函数级数近似解。

传统的高斯投影变换,以大地纬度 B 作为积分变量。李厚朴、边少锋等(2007, 2009, 2012, 2015~2017)研究了基于复变量等角纬度 φ c 的非迭代的正反解表达式[13]-[18]

金立新等(2017),在基于复变量等角纬度 φ c 高斯投影变换的基础上,将其解析到实数域,给出了实数域非迭代的正反解表达式及其变换系数[19]

实际上,非迭代的高斯投影变换,正解与反解为2套不同的级数表达式与系数,表达冗长。

本文在上述研究的基础上,从数值计算的实用角度出发,对实数域基于等角纬度 φ 的高斯投影变换,正反解变换均采用正解公式,反解变换利用正解公式迭代计算。算法简明,程序精短。

2. 等角纬度高斯投影变换路径

等角纬度的高斯投影变换,其变换路径为:

复数变换路径:

1) 正算路径: ( B,l ) w φ c z

2) 反算路径:

a. 迭代法: ( B,l ) w φ c z

b. 非迭代法: ( B,l ) w ψ c z

实数变换路径:

1) 正算路径: ( B,l ) ( q,l ) ( φ x , φ y ) ( x,y )

2) 反算路径:

a. 迭代法: ( B,l ) ( q,l ) ( φ x , φ y ) ( x,y )

b. 非迭代法: ( B,l ) ( q,l ) ( ψ x , ψ y ) ( x,y )

上述路径中,方框内的变换为主变换,其它为辅助变换。

其中:

( B,l ) 为实数纬度及径差

w=q+il 为复数等量纬度

φ c = φ x +i φ y 为复数等角纬度

ψ c = ψ x +i ψ y 为复数底点纬度(复数等距离纬度)

z=x+iy 为复数直角坐标

3. 复变量等角纬度的高斯投影变换

3.1. 复变量正算变换

复变量的正算变换,由以下步骤完成:

1) 正算辅助转换 ( B,l )w

{ q=arctanh( sinB )earctanh( esinB ) w=q+il (1)

2) 正算辅助转换 w φ c

φ c =arcsin[ tanh( w ) ] (2)

3) 正算主变换 φ c z

李厚朴等(2009, 2012, 2018, 2021)给出了基于等角纬度复变量的主变换式[12] [14]-[16]

z=a d 0 φ c +a j=1 n d 2j sin( 2j φ c ) (3)

其中系数 d

d 0 =1 1 4 e 2 3 64 e 4 5 256 e 6 175 16384 e 8 441 65536 e 10 43659 65536 e 12

d 2 = 1 8 e 2 1 96 e 4 9 1024 e 6 901 184320 e 8 16381 5898240 e 10 + 2538673 4587520 e 12

d 4 = 13 768 e 4 + 17 5120 e 6 311 737282 e 8 18931 20643840 e 10 1803171 9175040 e 12

d 6 = 61 15360 e 6 + 899 430080 e 8 + 18757 27525120 e 10 + 461137 20643840 e 12

d 8 = 49561 41287680 e 8 + 175087 165150720 e 10 869251 20643840 e 12

d 10 = 179101 41287680 e 10 25387 1290240 e 12

式中, e 为椭球偏心率。

3.2. 复变量反算变换(非迭代变换)

复变量的反算变换的迭代法,可顺序由式(3)、(2)、式(1)迭代求解,此不赘述。

对于非迭代法,李厚朴等(2009, 2012, 2018, 2021)给出的复变量的反算变换式为[12] [14]-[16]

1) 反算主变换 ψ c z

ψ c = z a d 0 (4)

2 )反算辅助转换 w ψ c

w=q+il=arctanh[ sin( ψ c ) ]+ j=1 5 ξ 2j1 sin[ ( 2j1 ) ψ c ] (5)

3) 反算辅助转换 (B,l)w

{ l=Im( w ) q=Re( w ) φ=arcsin[ tanh( q ) ] B=φ+ j=1 n b 2j sin( 2jφ ) (6)

式中,辅助变换由 ψ c 反算 w 的级数系数 ξ 为:

ξ 1 = 1 4 e 2 1 64 e 4 + 1 3072 e 6 + 33 16384 e 8 + 2363 1310720 e 10

ξ 3 = 1 96 e 4 13 3072 e 6 13 8192 e 8 1057 1966080 e 10

ξ 5 = 11 7680 e 6 29 24576 e 8 2897 3932160 e 10

ξ 7 = 25 86016 e 8 727 1966080 e 10

ξ 9 = 53 737280 e 10

式中,辅助变换由 φ 反算 B 的级数系数 b 为:

b 2 = 1 2 e 2 + 5 24 e 4 + 1 12 e 6 + 13 360 e 8 + 3 160 e 10

b 4 = 7 48 e 4 + 29 240 e 6 + 811 11520 e 8 + 81 2240 e 10

b 6 = 7 120 e 6 + 81 1120 e 8 + 3029 53760 e 10

b 8 = 4279 161280 e 8 + 883 20160 e 10

b 10 = 2087 161280 e 10

4. 实变量等角纬度的高斯投影变换

基于复变量等角纬度的高斯投影变换,计算精简,宜为首选方法。

目前,大众计算工具Excel还不支持复数计算,因此,基于实数等角纬度的高斯投影变换还有很大的实用价值。

4.1. 实变量等角纬度与等量纬度关系式

在实数域中,等角纬度与等量纬度的关系是等角纬度高斯投影变换的桥梁。

等角变换中,有以下复数关系式:

{ sin( φ c )=tanh( w ) cos( φ c )= 1 cosh( w ) tanh( w )= sinh( w ) cosh( w ) (7)

以式(7)为基础,可以分别得到式(8)及(9):

{ sin( φ x +i φ y )=sin( φ x )cosh( φ y )+icos( φ x )sinh( φ y ) cos( φ x +i φ y )=cos( φ x )cosh( φ y )isin( φ x )sinh( φ y ) (8)

{ tanh( q+il )= sinh( q )cosh( q ) cosh ( q ) 2 cos ( l ) 2 +sinh ( q ) 2 sin ( l ) 2 +i sin( l )cos( l ) cosh ( q ) 2 cos ( l ) 2 +sinh ( q ) 2 sin ( l ) 2 1 cosh( q+il ) = cosh( q )cos( l ) cosh ( q ) 2 cos ( l ) 2 +sinh ( q ) 2 sin ( l ) 2 i sinh( q )sin( l ) cosh ( q ) 2 cos ( l ) 2 +sinh ( q ) 2 sin ( l ) 2 (9)

由(8-1)与(9-1)的实部与虚部分别相等,得到:

{ sin( φ x )cosh( φ y )= sinh( q )cosh( q ) cosh ( q ) 2 cos ( l ) 2 +sinh ( q ) 2 sin ( l ) 2 cos( φ x )sinh( φ y )= sin( l )cos( l ) cosh ( q ) 2 cos ( l ) 2 +sinh ( q ) 2 sin ( l ) 2 (10)

由(8-2)与(9-2)的实部与虚部分别相等,得到:

{ cos( φ x )cosh( φ y )= cosh( q )cos( l ) cosh ( q ) 2 cos ( l ) 2 +sinh ( q ) 2 sin ( l ) 2 sin( φ x )sinh( φ y )= sinh( q )sin( l ) cosh ( q ) 2 cos ( l ) 2 +sinh ( q ) 2 sin ( l ) 2 (11)

则,由式(10)及(11)即可得到 ( φ x , φ y ) ( q,l ) 的转换解析式[2] [19] [20]

{ tan( φ x )= sinh( q ) cos( l ) tanh( φ y )= sin( l ) cosh( q ) (12)

式(12)是等角纬度 ( φ x , φ y ) 与等量纬度 ( q,l ) 的转换关系式。

4.2. 实变量正算变换(解析式)

1) 正算辅助转换 ( B,l )( q,l )

{ q=arctanh( sinB )earctanh( esinB ) l=L L 0 (13)

式中, L 0 为投影带中央子午线经度, L

2) 正算辅助转换 ( q,l )( φ x , φ y )

由关系式(12)得到:

{ φ x =arctan[ sinh( q ) cos( l ) ] φ y =arctanh[ sin( l ) cosh( q ) ] (14)

3) 正算主变换 ( φ x , φ y )( x,y )

由于: sin( φ x +i φ y )=sin φ x cosh φ y +icos φ x sinh φ y

将上式代入式(3),得到实数域等角纬度的高斯投影正算主变换式:

{ x=a d 0 φ x +a j=1 5 d 2j sin( 2j φ x )cosh( 2j φ y ) y=a d 0 φ y +a j=1 5 d 2j cos( 2j φ x )sinh( 2j φ y ) (15)

4.3. 实变量反算变换

4.3.1. 实变量反算迭代法

迭代算法,反算变换只需利用正算式及其变换系数,计算简明。

1) 反算主变换 ( φ x , φ y )( x,y )

由正算主变换式(15)可得到反算主变换式(16),迭代初值可设置为 ( φ x , φ y )=( 0,0 )

{ φ x = 1 a d 0 [ xa j=1 5 d 2j sin( 2j φ x ) cosh( 2j φ y ) ] φ y = 1 a d 0 [ ya j=1 5 d 2j cos( 2j φ x ) sinh( 2j φ y ) ] (16)

2) 反算辅助转换 ( q,l )( φ x , φ y )

由关系式(12),可得到方程组(17),迭代初值可设置为 ( q,l )=( 0,0 )

{ q=arcsinh[ tan( φ x )cos( l ) ] l=arcsin[ tanh( φ y )cosh( q ) ] (17)

3) 反算辅助转换 ( B,l )( q,l )

由关系式(13),可得到式(18),迭代初值可设置为 B=0

{ B=arcsin{ tanh[ q+earctanh( esinB ) ] } L= L 0 +l (18)

上述各迭代计算中,一般迭代5~6次,角度计算精度 < 1010 rad,直角坐标计算精度 < 1 mm,满足高斯投影变换精度要求。为保证更高计算精度及简化程序起见,本文迭代次数取固定值10次。

4.3.2. 实变量反算非迭代法

基于等角纬度的非迭代反算变换,转换式与正算完全不同、系数冗长[19]

1) 反算主变换 ( ψ x , ψ y )( x,y )

{ ψ x = x a d 0 ψ y = y a d 0 (19)

2) 反算辅助变换 ( q,l )( ψ x , ψ y )

{ q=arctanh( sin( ψ x ) cosh( ψ y ) )+ j=1 5 ξ 2j1 sin( ( 2j1 ) ψ x )cosh( ( 2j1 ) ψ y ) l=arctan( sinh( ψ y ) cos( ψ x ) )+ j=1 5 ξ 2j1 cos( ( 2j1 ) ψ x )sinh( ( 2j1 ) ψ y ) (20)

3) 反算辅助变换 ( B,l )( q,l )

( φ=arcsin( tanh( q ) ) B=φ+ j=1 5 b 2j sin( 2jφ ) (21)

5. 高斯变换与投影分带

传统的基于大地纬度的高斯投影变换式,展开为Taylor幂级数,初始展开点位于中央子午线,变换函数 ( x,y )=f( B,l ) 表达为径差 l 的幂级数,其变换精度与径差 l 的大小有关, l 越大精度越差,需要的级数项数越多。传统的基于大地纬度的高斯投影Taylor幂级数变换式,投影带宽限制在 < 7˚范围内,直角坐标的变换精度 < 1 mm [7]

基于等角纬度的高斯投影变换,变换式 ( x,y )=f( φ x , φ y ) 表达为等角纬度的三角函数, B q 的关系式(1),以及 ( q,l ) ( φ x , φ y ) 的关系式(2)均为解析式;正算变换式(3)与径差 l 无关。换言之,基于等角纬度的高斯投影变换,变换精度与径差 l 无关,投影带宽不限于7˚范围。仅就考虑变换精度而言,无需进行分带投影。由此,对于传统的3˚及6˚带的分带投影,可以进行跨带坐标变换,不限于传统相邻带坐标变换。

6. 变换实例

[实例1] 选自文献([5], p. 168)实例,点名1。坐标系为BJ54坐标系,投影带中央子午线 L 0 = 111 ,变换点大地坐标为: B= 31 04'41.6832" L= 111 47'24.8974"

解算:1) 由大地坐标 ( B,L ) 正算直角坐标 ( x,y )

2) 由直角坐标 ( x,y ) 反算大地坐标 ( B,L )

正算及迭代法反算结果详见表1。正算结果与原文献相同,本文计算精度高于原文献;反算结果 ( x,y ) 可精确还原到原大地坐标值 ( B,L )

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

等角纬度

φ x /rad

0.539495977396

0.539495977396

1

φ y /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= 45 ,径差 l= 3 ,中央子午线 L 0 = 0

b) B= 45 ,径差 l= 12 ,中央子午线 L 0 = 15

求上述大地坐标的直角坐标,以及由直角坐标反算溯源大地坐标。

不分带的正算结果详见表2,反算结果见表3。溯源反算的大地坐标与源坐标完全一致。实例2b的径差 l= 12 ,跨越了标准投影分带。

7. 结语

1) 传统的高斯投影变换采用大地纬度作为变换参数。本文利用李厚朴、边少锋、金立新等基于等角纬度作为变换参数的高斯投影变换的最新研究成果,在Excel VBA上实现了实数域的高斯投影变换。

2) 实数域等角纬度的高斯投影主变换式与复数域主变换式类似,同为三角函数级数表达式。

3) 根据复数域等角纬度 φ c 与等量纬度 w 的关系,导出了实数域 ( φ x , φ y ) ( q,l ) 的解析关系式。这一关系式是组构实数域等角纬度的高斯投影正反算变换的桥梁。

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次,角度计算精度 < 1010 rad,直角坐标变换精度 < 1 mm。迭代算法对计算效率的影响可忽略不计。

5) 基于等角纬度的高斯投影变换,变换精度与径差 l 无关,投影带宽不限于7˚范围。仅就考虑变换精度而言,无需进行分带投影。对于传统的3˚及6˚带的分带投影,可以进行跨带坐标变换,不限于传统相邻带坐标换带。

致 谢

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

基金项目

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

NOTES

*第一作者。

#通讯作者。

参考文献

[1] 杨启和. 等角投影数值变换的研究[J]. 测绘学报, 1982, 11(4): 268-282.
[2] 程阳. 复变函数与等角投影[J]. 测绘学报, 1985, 14(1): 51-60.
[3] 杨启和. 等角投影理论和方法[J]. 解放军测绘学院学报, 1994, 11(2): 133-139.
[4] 熊介. 椭球大地测量学[M]. 北京: 解放军出版社, 1988: 48-49.
[5] 陈健, 晁定波. 椭球大地测量学[M]. 北京: 测绘出版社, 1989: 18, 168.
[6] 杨启和. 地图投影变换原理与方法[M]. 北京: 测绘出版社, 1990: 56.
[7] 孔祥元, 郭际明, 刘宗泉. 大地测量学基础[M]. 第2版. 武汉: 武汉大学出版社, 1990: 115, 169.
[8] 边少锋, 张传定. Gauss投影的复变函数表示[J]. 测绘学院学报, 2001, 18(3): 157-159.
[9] 李厚朴, 边少锋. 高斯投影的复变函数表示[J]. 测绘学报, 2008, 37(1): 5-9.
[10] 边少锋, 许江宁. 计算机代数系统与大地测量数学分析[M]. 北京: 国防工业出版社, 2004: 107-110.
[11] 边少锋, 许江宁. 大地坐标系与大地基准[M]. 北京: 国防工业出版社, 2005: 99-101.
[12] 边少锋, 李厚朴. 大地测量计算机代数分析[M]. 北京: 科学出版社, 2018: 91-95.
[13] 李厚朴, 边少锋. 等量纬度展开式的新解法[J]. 海洋测绘, 2007, 27(4): 6-10.
[14] 李厚朴, 王瑞, 边少锋. 复变函数表示的高斯投影非迭代公式[J]. 海洋测绘, 2009, 29(6): 17-20.
[15] 李厚朴, 边少锋, 李海波. 常用等角投影及其解析变换的复变函数表示[J]. 测绘科学技术学报, 2012, 29(2): 109-112, 117.
[16] 边少锋, 李厚朴. 高斯投影的复变函数表示[M]. 北京: 科学出版社, 2021: 127-130.
[17] 李厚朴, 边少锋, 李海波. 利用复变函数高斯换带的方法[J]. 海军工程大学学报, 2016, 28(1): 15-19.
[18] 边少锋, 李厚朴, 李忠美. 地图投影计算机代数分析研究进展[J]. 测绘学报, 2017, 46(10): 1557-1569.
[19] 金立新, 许常文, 魏桂华. 高斯投影复变函数表示的实数解[J]. 海洋测绘, 2017, 37(2): 27-31.
[20] 刘勇, 李忠美, 刘佳奇, 钟业勋. 复数高斯投影若干数学性质分析[J]. 海洋测绘, 2023, 43(6): 77-82.