由高斯-克呂格投影平面直角坐標反解地理坐標的方法

1、由高斯-克呂格投影平面直角坐標反解地理坐標的方法鐘業(yè)勛 魏文展（廣西師范學(xué)院 南寧市明秀東路19號 530001） 【摘要】 根據(jù)高斯-克呂格投影平面直角坐標x、y的公式，設(shè)計了適用于CASIO fx-4800P計算器的程序。提出了利用正算程序采用有理逼近法反解地理坐標、的步驟和方法。經(jīng)試驗表明，本法可反解出精確的、值且收斂較快。測繪信息網(wǎng)【關(guān)鍵詞】 高斯-克呂格投影 平面直角坐標 地理坐標 計算程序 反算我國現(xiàn)行的大于1：50萬比例尺的各種地形圖，都采用高斯克呂格投影1。根據(jù)地面點的緯度和該點對投影帶的中央經(jīng)線的經(jīng)差，可用坐標公式求解該點的平面直角坐標X、Y。由于用X、Y反解、的方程建立困難

2、，筆者經(jīng)過大量試驗，探討了根據(jù)正解公式，用工程測繪中常用的CASIO fx-4800P計算器建立程序，通過計算逐漸逼近的反解、的數(shù)值解法。測繪信息網(wǎng)1 由X、Y坐標正算公式反解地理坐標的原理根據(jù)給定的坐標X、Y，可從已知的坐標表中確定所在的球面上的區(qū)間，即具有一定緯差和經(jīng)差的球面梯形內(nèi)。在給定點附近，可以任意給定一個確知其緯度為，經(jīng)差為的點，由，可按正算程序算出、，由于，并非求解的真值，故X與的差的絕對值和Y與的差的絕對值必然大于零，即（1）在高斯克呂格投影中，同一緯差對應(yīng)的X坐標增量X 為緯度的函數(shù)，其趨勢是隨緯度的增加而增大。點位對中央經(jīng)線的距離不同，同一經(jīng)差所對應(yīng)的Y坐標增量Y 隨著點位

3、離中央經(jīng)線愈遠而遞減。這一規(guī)律已在表1中呈現(xiàn)。表1以緯度間隔為5°，經(jīng)差間隔為30計算了各經(jīng)緯線交點的坐標值，并以0.1的緯差和經(jīng)差，分別算出相應(yīng)的坐標增X、Y。根據(jù)X、Y所在區(qū)間的X、Y值，可對第一次給定的近似值和進行改正：測繪信息網(wǎng)（2）由，可以正算出X1、Y1。設(shè)為改正的次數(shù)，求解的通式為：（3）當(dāng)時，即為（2）式。每次改正，將使更加接近，最終會出現(xiàn)的結(jié)果。而此時的對應(yīng)的即為所求。確定X、Y所在的區(qū)間及其相應(yīng)的緯度每增加0.1分的縱坐標增量X和緯差每增加0.1分的橫坐標增量Y可參考表1。高斯克呂格投影平面直角坐標x、y及坐標增量x、y單位：m表1x y 0 °00 0

4、 °30 1 °00 1 °30 X x x x x x x x Y y y y y y y y 85°009443514.60.0 186.116.2 9443535.74867.3 186.314.6 9443599.29734.5 186.413.0 9443705.014601.8 186.511.4 80°008885144.00.0 186.132.3 8885185.79696.8 186.430.7 8885310.719393.7 186.629.1 8885519.029090.9 186.927.5 75°008

5、326941.50.0 186.048.2 8327002.414451.0 186.446.6 8327185.128902.0 186.845.0 8327489.643353.0 187.243.5 70°007768984.40.0 185.963.6 7769062.419093.3 186.562.1 7769297.538186.6 187.060.6 7769689.057280.0 187.459.0 65°007211342.50.0185.878.6 7211435.823587.2 186.477.1 7211715.647171.1 187.075

6、.4 7212181.970748.3 187.773.4 60°006654075.90.0 185.793.0 6654181.427900.1 186.391.6 6654497.655800.7 187.190.2 6655024.883702.4 187.888.8 55°006097233.20.0 185.5106.7 6097347.531997.3 186.3105.3 6097690.663996.0 187.1104.0 6098262.695997.4 187.8102.7 50°005540849.60.0 185.4119.5 5540

7、969.435848.2 186.2118.2 5541328.971697.8 187.0117.1 5541928.2107550.6 187.7115.8 45°004984946.70.0 185.2131.4 4985068.339423.4 186.1130.3 4985433.378846.9 186.8129.1 4986041.5118270.3 187.7128.0 40°004429531.10.0 185.1142.3 4429650.842697.4 185.9141.3 4430010.185397.4 186.7140.2 4430609.11

8、28102.6 187.4139.2 35°003874594.70.0 184.9152.1 3874708.945644.8 185.7151.2 3875051.791293.5 186.4150.3 3875623.0136950.2 187.2149.4 30°003320114.90.0 184.8160.8 3320220.248243.6 185.5160.0 3320536.096490.0 186.2159.2 3321062.4144742.0 186.9158.4 25°002766055.50.0 184.6168.2 2766148.5

9、50469.6 185.3167.4 2766427.7100906.2 185.9165.9 2766892.9151277.1 186.5162.9 20°002212367.30.0 184.5174.4 2212445.452323.7 185.0173.8 2212679.7104647.9 185.5173.3 2213070.2156973.3 186.1172.8 15°001658990.40.0 184.4179.2 1659051.153775.3 184.8178.8 1659233.3107550.6 185.3178.4 1659537.0161

10、325.9 185.7178.0 10°001105855.30.0 184.4182.7 1105896.954824.2 184.6182.5 1106021.5109675.5 185.0182.2 1106229.4164580.9 185.2181.9 5°00552885.70.0 184.3184.8 552906.855449.8 184.4184.7 552970.1110901.8 184.6184.5 553075.5166358.3 184.8184.5 0°000.00.0 184.3185.5 0.055660.5 184.3185.5

11、 0.0111325.2 184.3185.6 0.0166998.5 184.4185.6  續(xù)表1x y 2 °00 2 °30 3 °00 x x x x x x y y y y y y 85°009443853.119469.1 186.69.7 9444043.524336.3 186.78.2 9444276.229203.6 186.06.5 80°008885810.738788.6 187.125.9 8886185.848486.8 187.324.3 8886644.258185.8 187.522.6 75&#

12、176;008327916.057804.0 187.041.9 8328464.272255.0 187.940.4 8329134.286706.1 188.338.7 70°007770237.076373.4 187.957.5 7770941.795467.0 188.456.0 7771803.0114560.7 188.854.4 65°007212834.894315.4 188.270.9 7213674.0117869.2 188.867.6 7214699.6141406.2 189.363.6 60°006655763.0111605.7

13、188.487.4 6656712.1139511.2 189.185.9 6657872.4167419.2 189.884.6 55°006099063.3128002.9 188.6101.3 6100093.1160013.8 189.3100.0 6101352.1192031.6 190.098.7 50°005542767.2143408.0 188.6114.6 5543846.2179271.6 193.4113.3 5545165.4215142.9 190.1112.0 45°004986893.3157693.8 188.4126.7 49

14、87988.5197117.3 189.3125.5 4989327.6236540.8 190.0124.2 40°004431447.7170815.6 188.3138.2 4432526.4213539.1 189.1137.2 4433845.3256275.7 189.9136.1 35°003876423.2182618.9 187.9148.4 3877452.3228303.4 188.8147.6 3878710.7274007.9 189.6146.6 30°003321799.6193002.2 187.6157.6 3322747.824

15、1273.5 188.3156.8 3323907.2289558.7 189.1156.0 25°002767543.9201549.1 187.1158.1 2768380.4251689.3 187.8150.6 2769402.2301664.6 188.4140.0 20°002213617.1209300.5 186.7172.2 2214320.6261630.0 187.2171.5 2215180.7313962.4 187.9170.9 15°001659962.4215101.3 186.1177.5 1660509.5268876.8 18

16、6.6177.0 1661178.5322652.3 187.1176.5 10°001106520.6219567.6 185.6181.5 1106895.5274662.9 185.9180.9 1107354.2329894.3 186.3179.9 5°00553223.3221821.7 184.9184.3 553413.3277294.0 185.1184.2 553645.6332777.7 185.4184.1 0°000.0222684.6 184.4185.7 0.0278387.8 184.5185.7 0.0334112.4 184.5

17、185.7 注：緯度增加0.1分的x的增量為x；經(jīng)差增加0.1分的y的增量為y。2 由高斯克呂格投影直角坐標反解地理坐標的方法2.1高斯克呂格投影X、Y坐標公式2 3 測繪信息網(wǎng)（1）（2）（1）式中S為自赤道起到某一緯度的經(jīng)線長度（IAG-75橢球）：（3）（1）、（2）式中的N為地球橢球卯酉圈曲率半徑：（4）為所求點的緯度，以度計；為所求點對中央經(jīng)線的經(jīng)差，以弧度計。當(dāng)以度為單位時，X、Y公式中的應(yīng)用下式變換后再代入；（5）（4）式中的為橢球長半徑，e為第一偏心率。測繪信息網(wǎng)（1）、（2）式中的為輔助函數(shù)，。對于1980西安坐標系，采用IAG-75橢球，其參數(shù)為：=6378140米，e =

18、0.081819221，e=0.0820944692.2 用CASIO fx4800P計算器計算高斯克呂格投影直角坐標的程序程序中的符號與公式中的符號對應(yīng)如下： ，K-已知點的X坐標，C-緯度增加0.1分的X坐標增量，Q-已在點的Y坐標，D-經(jīng)差增加0.1分的Y坐標增量，V-解出的緯度值，W-解出的經(jīng)差值。2.2.1 程序測繪信息網(wǎng)在文件名NSXYAB下輸入下列程序：Fix5：lbl1：A，B：I=sinA：J=cosA：T=tanA2：G=（0.082094469J）2 ：S=111133.0046A-16038.528sin2A+16.833sin4A- 0.022sin6A+0.0000

19、3sin8A：L=B180 lbl2： lbl3： 測繪信息網(wǎng)lbl4：K，C，Q，DAbs（K-X）>0=>V=（A+（K-X）÷（600CA=VAbs（Q-Y）>0=>W=（B+（Q-Y）÷ （600DB=WGoto1=>FixmV，W2.2.2 建立程序的操作過程步 驟鍵 操 作顯 示 畫 面1（NEW）NSXYABFilename?NSXYAB 2PGM：NSXYAB1. COMP 2.BASE-N3.SD 4.LR5.Save formula3（COMP） （輸入程式2.2.1）3 算 例已知K=5300000米，Q=120096米，

20、求該點的緯度和對中央經(jīng)線的經(jīng)差。已知的坐標X、Y即K和Q。解：根據(jù)已知的X、Y，可從表1得知其區(qū)間為45°0050°00，1°302°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為、組的算術(shù)平均值。初始值（即輸入的初始A、B值）也選取三組：I：45°00，1°30，：46°00，1°34，：48°00，1°39以K、Q與各組A、B、C、D值組合，輸

21、入計算器，其結(jié)果如表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

22、 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 3 47.82296 1.60656 5299995.240 120316.2434 4 47.82300 1.60339 5299994.970 120078.6857 5 47.82305 1.60364 5300000.281 120097.2498 6 47.82304 1.60362 5299999.97

23、8 120095.9081 7 47.82304 1.60362 5300000.002 120096.0067 C=188.5 D=121.9 · A（度） B（度） X（米） Y（米） 1 48.00000 1.65000 5319747.884 123145.0736 2 47.82539 1.60831 5300268.757 120441.8469 3 47.82302 1.60358 5299997.199 120093.2011 4 47.82304 1.60362 5300000.012 120096.0098 5 47.82304 1.60362 5300000.0

24、00 120096.0000 比較三組的計算結(jié)果，知所求值為A=47.82304°，B=1.60362°。選定的初始值A(chǔ)0、B0，不會影響計算結(jié)果，但以愈接近真值的初始值愈好，收斂會更快。A、B被鎖定后，即為所求值，若用此值再次輸入，得到的Xi、Yi會更貼近已知的K、Q值，這是由于第5位以后的小數(shù)在改變而獲得更精確的結(jié)果。4 結(jié) 論測繪信息網(wǎng)在測繪工作中，控制點的平面直角坐標可以從成果表中獲得精確的數(shù)值，但該點的地理坐標卻不能直接獲得。通過某種比例尺地圖進行量測，由于受量算條件的影響，也不易得到較精確的經(jīng)緯度數(shù)值。本文提出了利用正解平面直角坐標X、Y的公式，采用逐漸逼近的解

25、析解法，可以解決由平面直角坐標反解地理坐標的問題并可獲得足夠的精確度，完全能滿足1：500測圖的需要，因而本方法具有實踐意義和應(yīng)用價值。參考文獻：測繪信息網(wǎng)1 尹貢白、王家耀、田德森 、黃采芝，地圖概論. 北京：測繪出版社，1991年2 胡毓鉅，龔劍文，地圖投影（第二版），北京：測繪出版社，1992年3 方炳炎編. 地圖投影計算用表. 北京：測繪出版社，1992年作者簡介：鐘業(yè)勛，教授，現(xiàn)主要從事地圖學(xué)理論研究。代表成果：地圖模糊矩陣模型與制圖術(shù)語表述數(shù)學(xué)化研究；地圖易讀性度量研究；地圖數(shù)學(xué)定義的研究；地圖圖像系統(tǒng)為布爾代數(shù)系的證明等。已發(fā)表論文60 余篇。E-mail：weiwzThe Me

26、thod of Inverse Calculation GeographicCoordinate from Rectangular Coordinate inThe plane of Gauss-Krüger ProjectionZHONG Yexun WEI Wenzhan（Guangxi Teachers Collage， 19 East Mingxiu Road， Nanning， China， 530001）Abstract：According to the formula of rectangular coordinate in the plane of Gauss-Krüger projection， t