由高斯-克吕格投影平面直角坐标反解地理坐标的方法
钟业勋 魏文展
（广西师范学院 南宁市明秀东路19号 530001）
【摘 要】 根据高斯-克吕格投影平面直角坐标x、y的公式，设计了适用于CASIO fx-4800P计算器的程序。提出了利用正算程序采用有理逼近法反解地理坐标 、 的步骤和方法。经试验表明，本法可反解出精确的、值且收敛较快。
【关键词】 高斯-克吕格投影 平面直角坐标 地理坐标 计算程序 反算
我国现行的大于1：50万比例尺的各种地形图，都采用高斯——克吕格投影[1]。根据地面点的纬度和该点对投影带的中央经线的经差，可用坐标公式求解该点的平面直角坐标X、Y。由于用X、Y反解、的方程建立困难，笔者经过大量试验，探讨了根据正解公式，用工程测绘中常用的CASIO fx-4800P计算器建立程序，通过计算逐渐逼近的反解、的数值解法。
1 由X、Y坐标正算公式反解地理坐标的原理
根据给定的坐标X、Y，可从已知的坐标表中确定所在的球面上的区间，即具有一定纬差和经差的球面梯形内。在给定点附近，可以任意给定一个确知其纬度为 ，经差为 的点，由，可按正算程序算出 、 ，由于，并非求解的真值，故X与的差的绝对值和Y与的差的绝对值必然大于零，即
 （1）
在高斯——克吕格投影中，同一纬差 对应的X△坐标增量X 为纬度的函数，其趋势是随纬度的增加而增大。点位对中央经线的距离不同，同一经差 所对应的Y坐标增量△Y 随着点位离中央经线愈远而递减。这一规律已在表1中呈现。表1以纬度间隔为5°，经差间隔为30′计算了各经纬线交点的坐标值，并以0.1′的纬差和经差，分别算出相应的坐标增△X、△Y。根据X、Y所在区间的△X、△Y值，可对第一次给定的近似值 和进行改正：
 （2）
由 ， ，可以正算出X1、Y1。设 为改正的次数，求解 的通式为：
 （3）
当 时，即为（2）式。每次改正，将使 更加接近，最终会出现 的结果。而此时的 对应的 即为所求。确定X、Y所在的区间及其相应的纬度每增加0.1分的纵坐标增量△X和纬差每增加0.1分的横坐标增量△Y可参考表1。
高斯——克吕格投影平面直角坐标x、y及坐标增量△x、△y
单位：m 表1
x
y |
0 °00′ |
0 °30′ |
1 °00′ |
1 °30′ |
X |
△x |
x |
△x |
x |
△x |
x |
△x |
Y |
△y |
y |
△y |
y |
△y |
y |
△y |
85°00′ |
9443514.6
0.0 |
186.1
16.2 |
9443535.7
4867.3 |
186.3
14.6 |
9443599.2
9734.5 |
186.4
13.0 |
9443705.0
14601.8 |
186.5
11.4 |
80°00′ |
8885144.0
0.0 |
186.1
32.3 |
8885185.7
9696.8 |
186.4
30.7 |
8885310.7
19393.7 |
186.6
29.1 |
8885519.0
29090.9 |
186.9
27.5 |
75°00′ |
8326941.5
0.0 |
186.0
48.2 |
8327002.4
14451.0 |
186.4
46.6 |
8327185.1
28902.0 |
186.8
45.0 |
8327489.6
43353.0 |
187.2
43.5 |
70°00′ |
7768984.4
0.0 |
185.9
63.6 |
7769062.4
19093.3 |
186.5
62.1 |
7769297.5
38186.6 |
187.0
60.6 |
7769689.0
57280.0 |
187.4
59.0 |
65°00′ |
7211342.5
0.0 |
185.8
78.6 |
7211435.8
23587.2 |
186.4
77.1 |
7211715.6
47171.1 |
187.0
75.4 |
7212181.9
70748.3 |
187.7
73.4 |
60°00′ |
6654075.9
0.0 |
185.7
93.0 |
6654181.4
27900.1 |
186.3
91.6 |
6654497.6
55800.7 |
187.1
90.2 |
6655024.8
83702.4 |
187.8
88.8 |
55°00′ |
6097233.2
0.0 |
185.5
106.7 |
6097347.5
31997.3 |
186.3
105.3 |
6097690.6
63996.0 |
187.1
104.0 |
6098262.6
95997.4 |
187.8
102.7 |
50°00′ |
5540849.6
0.0 |
185.4
119.5 |
5540969.4
35848.2 |
186.2
118.2 |
5541328.9
71697.8 |
187.0
117.1 |
5541928.2
107550.6 |
187.7
115.8 |
45°00′ |
4984946.7
0.0 |
185.2
131.4 |
4985068.3
39423.4 |
186.1
130.3 |
4985433.3
78846.9 |
186.8
129.1 |
4986041.5
118270.3 |
187.7
128.0 |
40°00′ |
4429531.1
0.0 |
185.1
142.3 |
4429650.8
42697.4 |
185.9
141.3 |
4430010.1
85397.4 |
186.7
140.2 |
4430609.1
128102.6 |
187.4
139.2 |
35°00′ |
3874594.7
0.0 |
184.9
152.1 |
3874708.9
45644.8 |
185.7
151.2 |
3875051.7
91293.5 |
186.4
150.3 |
3875623.0
136950.2 |
187.2
149.4 |
30°00′ |
3320114.9
0.0 |
184.8
160.8 |
3320220.2
48243.6 |
185.5
160.0 |
3320536.0
96490.0 |
186.2
159.2 |
3321062.4
144742.0 |
186.9
158.4 |
25°00′ |
2766055.5
0.0 |
184.6
168.2 |
2766148.5
50469.6 |
185.3
167.4 |
2766427.7
100906.2 |
185.9
165.9 |
2766892.9
151277.1 |
186.5
162.9 |
20°00′ |
2212367.3
0.0 |
184.5
174.4 |
2212445.4
52323.7 |
185.0
173.8 |
2212679.7
104647.9 |
185.5
173.3 |
2213070.2
156973.3 |
186.1
172.8 |
15°00′ |
1658990.4
0.0 |
184.4
179.2 |
1659051.1
53775.3 |
184.8
178.8 |
1659233.3
107550.6 |
185.3
178.4 |
1659537.0
161325.9 |
185.7
178.0 |
10°00′ |
1105855.3
0.0 |
184.4
182.7 |
1105896.9
54824.2 |
184.6
182.5 |
1106021.5
109675.5 |
185.0
182.2 |
1106229.4
164580.9 |
185.2
181.9 |
5°00′ |
552885.7
0.0 |
184.3
184.8 |
552906.8
55449.8 |
184.4
184.7 |
552970.1
110901.8 |
184.6
184.5 |
553075.5
166358.3 |
184.8
184.5 |
0°00′ |
0.0
0.0 |
184.3
185.5 |
0.0
55660.5 |
184.3
185.5 |
0.0
111325.2 |
184.3
185.6 |
0.0
166998.5 |
184.4
185.6 |
续表1
x
y |
2 °00′ |
2 °30′ |
3 °00′ |
x |
△x |
x |
△x |
x |
△x |
y |
△y |
y |
△y |
y |
△y |
85°00′ |
9443853.1
19469.1 |
186.6
9.7 |
9444043.5
24336.3 |
186.7
8.2 |
9444276.2
29203.6 |
186.0
6.5 |
80°00′ |
8885810.7
38788.6 |
187.1
25.9 |
8886185.8
48486.8 |
187.3
24.3 |
8886644.2
58185.8 |
187.5
22.6 |
75°00′ |
8327916.0
57804.0 |
187.0
41.9 |
8328464.2
72255.0 |
187.9
40.4 |
8329134.2
86706.1 |
188.3
38.7 |
70°00′ |
7770237.0
76373.4 |
187.9
57.5 |
7770941.7
95467.0 |
188.4
56.0 |
7771803.0
114560.7 |
188.8
54.4 |
65°00′ |
7212834.8
94315.4 |
188.2
70.9 |
7213674.0
117869.2 |
188.8
67.6 |
7214699.6
141406.2 |
189.3
63.6 |
60°00′ |
6655763.0
111605.7 |
188.4
87.4 |
6656712.1
139511.2 |
189.1
85.9 |
6657872.4
167419.2 |
189.8
84.6 |
55°00′ |
6099063.3
128002.9 |
188.6
101.3 |
6100093.1
160013.8 |
189.3
100.0 |
6101352.1
192031.6 |
190.0
98.7 |
50°00′ |
5542767.2
143408.0 |
188.6
114.6 |
5543846.2
179271.6 |
193.4
113.3 |
5545165.4
215142.9 |
190.1
112.0 |
45°00′ |
4986893.3
157693.8 |
188.4
126.7 |
4987988.5
197117.3 |
189.3
125.5 |
4989327.6
236540.8 |
190.0
124.2 |
40°00′ |
4431447.7
170815.6 |
188.3
138.2 |
4432526.4
213539.1 |
189.1
137.2 |
4433845.3
256275.7 |
189.9
136.1 |
35°00′ |
3876423.2
182618.9 |
187.9
148.4 |
3877452.3
228303.4 |
188.8
147.6 |
3878710.7
274007.9 |
189.6
146.6 |
30°00′ |
3321799.6
193002.2 |
187.6
157.6 |
3322747.8
241273.5 |
188.3
156.8 |
3323907.2
289558.7 |
189.1
156.0 |
25°00′ |
2767543.9
201549.1 |
187.1
158.1 |
2768380.4
251689.3 |
187.8
150.6 |
2769402.2
301664.6 |
188.4
140.0 |
20°00′ |
2213617.1
209300.5 |
186.7
172.2 |
2214320.6
261630.0 |
187.2
171.5 |
2215180.7
313962.4 |
187.9
170.9 |
15°00′ |
1659962.4
215101.3 |
186.1
177.5 |
1660509.5
268876.8 |
186.6
177.0 |
1661178.5
322652.3 |
187.1
176.5 |
10°00′ |
1106520.6
219567.6 |
185.6
181.5 |
1106895.5
274662.9 |
185.9
180.9 |
1107354.2
329894.3 |
186.3
179.9 |
5°00′ |
553223.3
221821.7 |
184.9
184.3 |
553413.3
277294.0 |
185.1
184.2 |
553645.6
332777.7 |
185.4
184.1 |
0°00′ |
0.0
222684.6 |
184.4
185.7 |
0.0
278387.8 |
184.5
185.7 |
0.0
334112.4 |
184.5
185.7 |
注：纬度增加0.1分的x的增量为△x；经差增加0.1分的y的增量为△y。
2 由高斯——克吕格投影直角坐标反解地理坐标的方法
2.1 高斯——克吕格投影X、Y坐标公式[2] [3]
  （1）
 （2）
（1）式中S为自赤道起到某一纬度的经线长度（IAG-75椭球）：
  （3）
（1）、（2）式中的N为地球椭球卯酉圈曲率半径：
 （4）
为所求点的纬度，以度计；为所求点对中央经线的经差，以弧度计。当以度为单位时，X、Y公式中的应用下式变换后再代入；
 （5）
（4）式中的 为椭球长半径，e为第一偏心率。
（1）、（2）式中的 为辅助函数， 。
对于1980西安坐标系，采用IAG-75椭球，其参数为：=6378140米，e =0.081819221，e′=0.082094469
2.2 用CASIO fx4800P计算器计算高斯——克吕格投影直角坐标的程序
程序中的符号与公式中的符号对应如下：
 ， ，K-已知点的X坐标，C-纬度增加0.1分的X坐标增量，Q-已在点的Y坐标，D-经差增加0.1分的Y坐标增量，V-解出的纬度值，W-解出的经差值。
2.2.1 程序
在文件名[NSXYAB]下输入下列程序：
Fix5：lbl1：{A，B}：I=sinA：J=cosA：T=tanA 2：G=（0.082094469J）2 ：
：S=111133.0046A-16038.528sin2A+16.833sin4A- 0.022sin6A+0.00003sin8A：L=Bπ 180 ◢
lbl2：
◢
lbl3：  ◢
lbl4：{K，C，Q，D}∶Abs（K-X）>0=>V=（A+（K-X）÷（600C∶A=V∶Abs（Q-Y）>0=>W=（B+（Q-Y）÷
（600D∶B=W∶Goto1∶=>Fixm{V，W}◢
2.2.2 建立程序的操作过程
步 骤 |
键 操 作 |
显 示 画 面 |
1 |
（NEW）NSXYAB |
Filename?
[NSXYAB]
|
2 |
|
PGM：NSXYAB
1. COMP 2.BASE-N
3.SD 4.LR
5.Save formula
|
3 |
（COMP） （输入程式2.2.1） |
▂ |
3 算 例
已知K=5300000米，Q=120096米，求该点的纬度和对中央经线的经差。已知的坐标X、Y即K和Q。
解：根据已知的X、Y，可从表1得知其区间为45°00′～50°00′，1°30′～2°00′。为检验解算精度和解算的正确性，选取三组△X、△Y值，即在程序中采用三组C、D值：I：C=187.7，D=128，Ⅱ：C=189.3，D=115.8，Ⅲ：C=188.5，D=121.9，第Ⅲ组的C、D为Ⅰ、Ⅱ组的算术平均值。
初始值 （即输入的初始A、B值）也选取三组：I：45°00′，1°30′，Ⅱ：46°00′，1°34′，Ⅲ：48°00′，1°39′以K、Q与各组A、B、C、D值组合，输入计算器，其结果如表2。
表2
Ⅰ |
C=187.7 D=128 |
ι· |
A（度） |
B（度） |
X（米） |
Y（米） |
1 |
45.00000 |
1.50000 |
4986041.538 |
118270.3298 |
2 |
47.78777 |
1.52377 |
5295957.105 |
114192.8786 |
4 |
47.82367 |
1.60064 |
5300064.738 |
119870.8526 |
5 |
47.82309 |
1.60357 |
530005.377 |
120091.8042 |
6 |
47.82304 |
1.60362 |
5300000.154 |
120096.0085 |
7 |
47.82304 |
1.60362 |
5300000.002 |
120096.0034 |
Ⅱ |
C=189.3 D=115.8 |
ι· |
A（度） |
B（度） |
X（米） |
Y（米） |
1 |
46.00000 |
1.56667 |
5097256.772 |
120082.0878 |
2 |
47.78503 |
1.55020 |
5295691.549 |
116180.2048 |
| |