1、徐州師范大學本科畢業設計論文1摘要：首先簡要論述了空間參照系統、常用坐標系等基礎知識。然后研究了坐標系轉換的算法，其中主要研究了坐標系之間的坐標轉換方法，如大地坐標與空間直角坐標的轉換、不同坐標下的空間直角坐標轉換、高斯投影坐標變換算法和簡化算法。另外介紹了 WGS 一 84 坐標變換為 BJ 一 54 坐標的四參數轉換模型。最后通過 Visual Baisc 語言編程構建坐標轉換實驗系統，實現 Windows 窗口界面化的輸入和顯示方式，快速實現了坐標系之間的一系列轉換。對于求底點緯度 Bf和空間直角坐標轉換為大地坐標求 B 時，采用迭代法進行求解比較方便，迭代次數較少，精度也很高。在求解四

2、參數的過程中，按間接平差最小二乘原理得到精確的轉換參數。 經過計算實例的驗證，這些程序可以用于實際作業,具有一定的實用價值。關鍵字： 坐標系統；空間直角坐標系；大地坐標系；高斯平面坐標系；1954北京坐標系；1980 西安坐標系；WGS-84 世界大地坐標系；2000 國家大地坐標系；坐標系統轉換2Abstract： First of all, briefly addressed the spatial reference system used, such as basic knowledge of coordinate system. Coordinate system and then

3、study the conversion algorithm, which studied the coordinates between the coordinate system conversion methods, such as earth and space coordinates of the rectangular coordinate conversion, coordinate different Cartesian coordinate space conversion, Gaussian projection coordinate transformation algo

4、rithm and simplify the algorithm . Also introduced WGS coordinate transformation of a 84 a 54 coordinates for the BJs four-parameter transformation model. Finally, through the programming language Visual Baisc coordinate transformation to build the experimental system, the Windows window interface o

5、f input and display, and rapid realization of the series between the coordinate system conversion. For the end point for Bf and latitude coordinates into a Cartesian space geodetic coordinates for B, the iterative method used to solve the more convenient, less number of iterations, the accuracy is a

6、lso high. Solving the four parameters in the process, according to the principle of indirect least squares adjustment of the conversion parameters to be accurate. After verification example, these procedures can be used in practical operation, has a certain practical value.Keyword:Coordinate system;

7、 space Cartesian coordinate system; the earth coordinate system; Gauss plane coordinate system; coordinate system 1954 Beijing; Xian in 1980 coordinate system; WGS-84 World Geodetic Coordinate System; 2000 National Geodetic coordinate system; coordinate system conversion3前言：為了確定地面點的空間位置，需要建立坐標系，一個點在

8、空間的位置需要三個坐標量來表示。由于衛星在大地測量中的迅速發展，我們使用大地坐標系、空間直角坐標系和 WGS-84 坐標系，而一般的測量工作中，為了方便工程建設規劃、設計，我們通常采用高斯平面直角坐標系。自從建國以來，我國分別根據當時的技術要求建立的 1954 北京坐標系和 1980 西安坐標系，測繪了各種比例尺的地形圖，在各種生產建設中發揮了重要的重用。由于技術條件的限制，54 坐標系采用的是克拉索夫斯基橢球體參數，該橢球在計算和定位的過程中，沒有采用中國的數據，該系統在我國范圍內符合得不好，不能滿足高精度定位以及地球科學、空間科學和戰略武器發展的需要；繼而在 70-80 年代，我們的測繪工

9、作者經過大量的實地測繪，建立了 1980 西安坐標系，它采用的是 1975 年 IUGG 第十六屆大會推薦的參考橢球參數。由于以全球衛星導航定位系統為主的現代空間定位技術快速發展，要求我們采用地心坐標系，能和 WGS-84 坐標系一樣，在我國航天、民航、海事、海洋、交通、地震、水利、建設、規劃、地質調查、國土資源管理等部門得到廣泛的應用，所以我們建立了 2000 國家大地坐標系。這一系列坐標系的建立與使用，固然對我國的國民經濟建設、社會發展產生了巨大的社會效益，但是由于資料的保存不當，以及實際參數的修改，很多已有的地形圖轉換到現有的坐標系中需要進行大量的工作，進行坐標系統的轉換。比如當礦山進行

10、新的工程建設時,如果舊坐標系統符合礦山測量的要求,可將測量成果統一至舊坐標系統上,即進行新、舊坐標系統的換算;反之,如果因開采范圍的擴大或舊坐標系統原本不正確使舊坐標系統不符合要求時,則應重新建立新的統一坐標系統。本文主要介紹的是幾種常用的坐標系統，進而討論他們之間的轉換方式，建立他們之間的數學轉換模型，利用 VB 編程實現坐標系統轉換模型的建立與設計。4目錄：摘要： 關鍵字： 前言： 第一章 坐標系統的簡介 .11.1 坐標系統概述.11.2 天球坐標系.11.3 地球坐標系.21.3.1 地心坐標系.31.3.2 參心坐標系.51.3.3 高斯平面直角坐標系.6第二章 坐標系統的轉換 .7

11、2.1 大地坐標系統（BLH）和空間直角坐標系統（XYZ）的轉換.72.2 大地坐標與高斯平面坐標之間的轉換.92.3 新舊坐標系之間的轉換.11第三章 坐標系統轉換程序的設計與調試 .133.1 流程圖.133.2 程序界面及使用方法.143.2.1 大地坐標和空間直角坐標的轉換.143.2.2 高斯坐標正算.153.2.3 高斯坐標反算.163.2.4 四參數轉換.17第四章 程序代碼 .194.1 大地坐標和空間直角坐標的轉換.194.2 高斯坐標轉換的正算.224.3 高斯坐標轉換的反算.234.4 四參數轉換.26結束語： .36致謝： .37參考文獻： .38徐州師范大學本科畢業設

12、計論文1第一章 坐標系統的簡介1.1 坐標系統概述 測繪工作必然需要一個完整的坐標體系，假如沒有一個完整獨立的坐標系，我們就無法把空間物體的位置、形狀和相互之間的關系表示出來。比如礦山的開采和礦井的挖掘與鋪設就離不開坐標系統的建立；對于人造衛星和天體的研究同樣需要坐標參考系統來支持。坐標參考系統分為天球坐標系統和地球坐標系統。天球坐標系統用于研究天體和人造衛星的定位與運動。地球坐標系用于研究地球上物體的定位與運動，是以旋轉橢球為參照體建立的坐標系統，它又分為大地坐標系和空間直角坐標系兩種形式。當然，由于分類的不同還有其他形式的坐標系統，本章主要介紹地球坐標系，根據坐標系原點所處的位置不同，地球

13、坐標系可以分為地心坐標系和參心坐標系。地心坐標系與參心坐標系的區別是就是地心坐標系是以地球的質心為原點建立的坐標參考系統，而參心坐標系是以參考橢球的中心為原點建立的坐標參考系統。比如 WGS-84 坐標系統和 2000 國家大地坐標系就屬于地心坐標系，我國采用的 1954 北京坐標系和 1980 國家大地坐標系就屬于參心坐標系。以上所說的關系可以簡單的用下圖表示：坐標系統天球坐標系地球坐標系地心坐標系參心坐標系圖 1-1 坐標系統的分類1.2 天球坐標系所謂天球，是指以地球質心M為中心，半徑r為任意長度的一個假想的球體。好似晴朗的夜晚仰望天空，滿天的星斗就分布在這個半圓球面上。2EsrPsEn

14、PnO天球赤道黃道在天文學中，通常均把天體投影到天球的球面上，并利用球面坐標系統來表達或研究天體的位置及天體之間的關系。為了建立球面坐標系統，必須確定球面上的一些參考點、線、面和圈。在全球定位系統中，為描述衛星的位置也將涉及到這些概念。對此，現按圖1-2介紹如下： 圖1-2 天球天軸與天極-地球自轉軸的延伸直線為天軸；天軸與天球的交點Pn和Ps稱為天極，其中Pn稱為北天極，Ps為南天極。天球赤道面與天球遲到-通過地球質心M與天軸垂直的平面，稱為天球赤道面。這時天球赤道面與地球赤道面相重。該赤道面與天球相交的大圓稱為天球赤道。天球子午面與子午圈-包含天軸并通過地球上任一點的平面，稱為天球子午面。

15、而天球子午面與天球相交的大圓稱為天球子午圈。時圈-通過天軸的平面與天球相交的半個大圓。黃道-地球公轉的軌道面與天球相交的大圓，即當地球繞太陽公轉時，地球上的觀測者所見到的太陽在天球上運動的軌跡。黃道面與赤道面的夾角稱為黃赤夾角，約為23.5。黃極-通過天球中心，且垂直與黃道面的直線與天球的交點，其中靠近北極的交點n稱為北黃極，靠近南天極的交點s為南黃極。春分點-當太陽在黃道上從天球南半球向北半球運行時，黃道與天球赤道的交點。在天文學和衛星大地測量學中，春分點和天球赤道面，是建立參考系的重要基準點和基準面。 （參考文獻【2】 ）1.3 地球坐標系地球坐標系用于研究地球上物體的定位與運動，是以旋轉

16、橢球為參照體建立的坐標系統，它又分為大地坐標系和空間直角坐標系兩種形式。大地坐標系如圖 1-3 所示，P 點的子午面 NPS 與起始子午面 NGS 所構 3成的二面角 L，叫做 P 點的大地經度，由起始子午面起算，向東為正，稱東經（0-180）向西為負，稱西經（0-180）,P 點的法線 Pn與赤道面的夾角 B，稱做 P 點的大地維度，由赤道面起算，向北為正，稱北緯（0-90） ，向南為負，稱南緯（0-90） 。在該坐標系中，P 點的位置用（L,B）表示。如果點不在橢球面上，表示點的位置除（L,B）外，還要附加另一參數大地高 H，它是從觀測點沿橢球的法線方向到橢球面的距離?？臻g直角坐標系如圖

17、1-4 所示，空間任意點的坐標用（X,Y,Z）表示，坐標原點位在總地球質心或者參考橢球中心，Z 軸與地球平均自轉軸相重合，亦即指向某一時刻的平均北極點，X 軸指向平均自轉軸與平均格林尼治天文臺所決定的子午面與赤道面的交點 Ge，而 Y 軸與 XOZ 平面垂直，且指向東為正。 N E P OB nNS G 圖 1-3大地坐標系Z X O Z P X Y Ge P2 Y圖 1-4 空間直角坐標系 根據坐標系原點所處的位置不同，地球坐標系可以分為地心坐標系和參心坐標系。1.3.1 地心坐標系地心坐標系按其應用有兩種表達形式，如圖（1）地心空間直角坐標系地心空間直角坐標系的定義是：原點O與地球質心重合

18、；Z軸指向地球北極，X軸指向格林尼治平子午面與地球赤道的交點E，Y軸垂直與XOY平面構成右手坐標系。（2）地心大地坐標系E 4BIH 定義的零子午圈（1984.0）地球質心Z X Y O CTP 赤道地心大地坐標系的定義是：地球橢球的中心與地球質心重合，橢球的短軸與地球自轉軸相合，大地緯度B為過地面點的橢球法線與橢球赤道面的夾角，大地經度L為過地面點的橢球子午面與格林尼治平大地子午面之間的夾角，大地高H為地面點沿橢球法線至橢球面的距離。Z X Y O L B H E 圖 1-5 地心坐標系 圖 1-6 WGS-84 世界大地坐標系下面介紹 2 種常用的地心坐標系：1、WGS-84 世界大地坐標

19、系WGS-84 世界大地坐標系是為 GPS 全球定位系統使用而建立的坐標系統，通過遍布世界的衛星觀測站觀測到的坐標建立。WGS-84 世界大地坐標系是一個協議地球參考系。該坐標系的原點是地球的質心，Z 軸指向 BIH1984.0定義的協議地球極 CTP 方向，X 軸指向BIH1984.0零度子午面和 CTP 對應的赤道的交點，Y 軸和 Z,X 軸構成右手坐標系。它的地球橢球基本參數:地球橢球長半徑 a=6378137m,地球橢球扁率a=0.003352810664742、2000 國家大地坐標系2000 中國大地坐標系符合 ITRS(國際地球參考系)的如下定義:1)原點在包括海洋和大氣的整個地

20、球的質量中心;2)長度單位為米(SI)。這一尺度同地心局部框架的 TCG(地心坐標時)時間坐標一致;3)定向在 1984. 0 時與 BIH(國際時間局)的定向一致;4)定向隨時間的演變由整個地球的水平構造運動無凈旋轉條件保證。以上定義對應一個直角坐標系,它的原點和軸定義如下(圖 1-6):5IERS 參考極IERS 參考子午面參考橢球地球質心Z X Y O 圖 1-6 2000 國家大地坐標系1)原點:地球的質量中心;2)Z 軸:指向 IERS 參考極方向;3)X 軸: IERS 參考子午面與通過原點且同 Z 軸正交的赤道面的交線;4)Y 軸:完成右手地心地固直角坐標系。CGCS2000 的

21、參考橢球的幾何中心與坐標系的原點重合,旋轉軸與坐標系的 Z 軸一致。參考橢球既是幾何應用的參考面,又是地球表面上及空間正常重力場的參考面。CGCS2000 參考橢球的定義常數是:長半軸 a= 6 378 137. 0 m扁率 f= 1/298. 257 222 101短半軸 b=6 356 752.3141 m1.3.2 參心坐標系參心坐標系按其應用又分為參心大地坐標系和參心空間直角坐標系兩種。 參心大地坐標系的應用十分廣泛，它是經典大地測量的一重通用坐標系。根據地圖投影理論，參心大地坐標系可以通過高斯投影計算轉化為平面直角坐標系，為地形測量和工程測量提供控制基礎。由于不同時期采用的地球橢球不

22、同或其定位與定向不同，我國歷史上出現的參心大地坐標系，主要有BJZ54（原） 、GDZ80和BJZ54等三種參心大地坐標系。 參心空間直角坐標系是用三維坐標 x、y、z 表示點位的，它可按一定的數學公式與參心大地直角坐標系相互換算。通常在由 GPS 定位結果（地心空間直角坐標系）計算參心大地直角坐標系時，作為中間一種過度換算的坐標系。 （參考文獻【1】 ）6 （1）1954 北京坐標系可以認為是前蘇聯 1942 年坐標系的延伸。他的原點不在北京，而在前蘇聯的普爾科沃。相應的橢球為克拉所夫斯基橢球。克拉所夫斯基橢球的參數為地球長半軸 a=6378245m 扁率 f=1：298.3（2）1980

23、年國家大地坐標系的大地原點定在我國中部，具體選址是陜西省涇陽縣永樂鎮；它采用國際大地測量和地球物理聯合會 1975 年推薦的四個地球橢球基本參數（a，J2，GM，w） ，并根據這四個參數求解橢球扁率和其他參數；這個坐標系的橢球短軸平行于地球質心指向我國地極原點JYD1968.0方向，大地起始子午面平行于格林尼治平均天文臺的子午面；橢球定位參數以我國范圍內高程異常值平方和等于最小為條件求解。1980 年國家大地坐標系的特點是：1. 地球橢球長半徑 a=6378140m 地球橢球扁率 a=1：298.2572. 參心大地坐標系是在 1954 年北京坐標系基礎上建立起來的。3. 橢球面同似大地水準面

24、在我國境內最為密合，是多點定位。4. 定向明確。橢球短軸平行于地球質心指向地極原點 JYD1968.0 的方向，起始大地子午面平行與我國起始子午面。5. 大地原點地處我國中部，位于西安市以北 60 千米處的涇陽縣永樂鎮。6. 大地高程采用 1956 年黃海高程系。1.3.3 高斯平面直角坐標系高斯平面直角坐標系采用的是等角投影，用一個橢圓柱橫套在地球的外面，并且是橢圓與地球的中央子午線相切，橢圓柱的中心必須和橢球的中心一致，然后采用等角投影，在中央子午線兩側按一定得級差把地球表面投影到橢圓柱上，再把橢圓柱展開就是高斯投影面。 中央子午線ONS 赤道SxyON7 圖 1-7 高斯投影 圖 1-8

25、 高斯平面坐標系在投影面上，中央子午線和赤道的投影都是直線，并且以中央子午線和赤道的交點 O 作為坐標原點，以中央子午線的投影為縱坐標抽，以赤道的投影為橫坐標軸，這樣便形成了高斯平面直角坐標系。在我國坐標 x 坐標都是正的，y 坐標的最大值約為 330km。為了避免出現負的橫坐標，可在橫坐標上加上 50000m。在坐標前面冠以帶號，這種坐標稱為國家統一坐標。8Z X Y O L B H T第二章 坐標系統的轉換2.1 大地坐標系統（BLH）和空間直角坐標系統（XYZ）的轉換對于任一地面點 T 在地球坐標系中的坐標，可表示為（X,Y,Z）或（B,L,H） 。如圖 1-9 可以清楚的看清楚大地坐標

26、系和空間直角坐標系的關系。大地坐標系轉換為空間直角坐標系的公式為： 圖 1-9 大地坐標系與空間直角坐標系LBcoscosHNX （2-1）LBHNYsincosBHeNZsin12BN22sine1a式中，N 為卯酉圈曲率半徑；e 為橢球第一偏心率。 不難看出，由大地坐標換算空間直角坐標的正解問題，可以由上訴公式直接加以解算。當用空間直角坐標系轉換為大地坐標系時，通常用下式來進行計算：2/1222sin1sin1arctanBeZBaetgB (2-2)XYLarctan9NBRHcoscos 式中， 2/122arctanYXZ2/1222ZYXRa,e 分別為所取橢球的長半徑和偏心率。在

27、這個公式里面求 B 采用迭代法比較好，按下述方法求 B 的初值 B0，可以減少迭代次數。BB0 (2-3)/(sin1RZ)2cos21 (sin2AAB)sin12/(222eRaeA2.2 大地坐標與高斯平面坐標之間的轉換 我國地形圖的測量控制和城市、礦山等區域性的測量控制一般都是采用的高斯平面直角坐標系，而現代大地測量很多時候是測得的 GPS 數據。所以我們需要把 GPS 坐標準換到高斯平面直角坐標系上，下面就從高斯投影坐標正算和反算公式來說明這個問題。首先介紹一下高斯投影坐標正算公式，也就是把大地坐標轉換成高斯平面直角坐標。以下是高斯投影坐標正算公式：6425644223422)586

28、1(cossin720N495cossin24cossin2xlttBBltBBNBlBNX 52224255322335814185cos1201cos6cosNyltttBNltBNBl式中：X 為自赤道量起的子午線弧長， （2-4）10BeN221sin1acosBetanBt，對于 X 的計算，我們一般采用下列實用公式進行計算，6/ )6sin(4/ )4sin(2/ )2sin()1 (43212BABABABAeaX25617564454316421eeeA 5125251615436422eeeA2561056415643eeA5123564eA 其中 a 為橢球長半徑，b 為橢

29、球短半徑，e 為橢球第一偏心率，X 為由赤道到緯度的子午線弧長，緯度的單位為弧度。為了對高斯投影坐標正算進行檢核，又產生了高斯坐標反算公式，也就是把高斯平面直角坐標轉換成大地坐標。以下是高斯坐標轉換反算公式：6425422223245tt9061720tt9t 3524tcos2tyNMyNMyBNMBBfffffffffffffffff522242532238624285cos120121cos61cos1ytttBNytBNyBNlfffffffffffff式中：為底點緯度，fBBeNtff221sin1acosBetanB，（2-5）（2-6）11，2/3222sin11ffBeeaM對

30、于其中底點緯度的計算，參考文獻【6】的內容，采用迭代法求fB解，不僅精度能夠滿足要求，而且免去了高斯投影坐標反算公式的推導過程的繁瑣。具體計算過程見程序代碼。2.3 新舊坐標系之間的轉換根據以上要求我們需要把以前的舊的坐標系轉換到 80 坐標系中，如第一章所介紹的，我國最先使用的是 1954 北京坐標系，繼而使用 1980 國家大地坐標系，如今逐漸轉而使用 2000 國家大地坐標系，當然我們現在用的GPS 是建立在 WGS-84 坐標系上的。 通常我們熟知的是布爾莎七參數模型，利用至少 3 個已知點的坐標來計算出轉換參數。由于實際情況的需要和精度要求的不同，我們也會采用四參數轉換來實現新舊坐標

31、系之間的轉換。四參數是同一個橢球內不同坐標系之間進行轉換的參數。需要特別注意的是參予計算的控制點原則上至少要用兩個或兩個以上的點，控制點等級的高低和分布直接決定了四參數的控制范圍。四參數的四個基本項分別是：X 平移、Y 平移、旋轉角和比例參數。四參數法轉換相比七參數法轉換有一定得缺陷，但是四參數法轉換比較容易實現，它需要的條件相對較低。它的基本思想可以概括如下：（1） 把有公共點的坐標分別投影到高斯平面上，得到它們在高斯平面對應的坐標。（2） 先求出平面內一組精度不太高的轉換參數，把所有點在某一（如WGS84）平面坐標轉換為另一（如 BJ54）坐標系下的近似平面坐標。（3） 從已知的公共點上求

32、出近似坐標與其準確坐標值之間的差值。（4） 利用最小二乘法求出需要的四參數，最終得到精確的轉換后的坐標。其計算方法如下（以 WGS84 和 BJ54 為例）：1.有一組公共點 P1,P2,PiPn，對于公共點 Pi 有：12 (2.7) 084848484i845484541001mYXXYYXYYXXi其中， 為一組精度不太高的初始坐標轉換參數。TmYXX00,取其中任意 2 個公共點 Pi 和 Pj，分別代入方程式（2.7）中，可以求得一組轉換參數。TmYXX00,2.用求出的精度不太高的轉換參數，將 WGS84 高斯平面坐標84WGSYX，轉換成近似的 54 北京高斯平面坐標，則坐標轉換

33、異54BJYX，常為： (2.8) YYXX5454yx3.利用最小二乘原理，求出平面坐標轉換的 4 個參數：，即，其中：TmYXX00,PLBPBBXTT1，TmYXX00,i1001iXYYXBiiYYXXL，P 取單位陣。TBBBBBn,i, 2, 1TLLiLLLn, 2, 14.利用求出的轉換參數將需要的點的坐標轉換為 BJ54 高斯平面坐標， YXYXmYXYX001（2.9） 13第三章 坐標系統轉換程序的設計與調試3.1 流程圖本節主要介紹如何把理論的公式的推導轉換為 VB 的流程圖，使程序的編譯與識別簡潔明了。分析需要轉換的坐標系，整理需要使用的轉換公式確定轉換需要使用的參數

34、（如橢球長半軸 a，偏心率 e）用 VB 語言對公式的計算進行電算化數據的輸入，并進行計算，得到結果結果的輸出，并進行適當的檢核，確保結果的正確程序的應用與說明圖 3-1 流程圖這里以四參數坐標轉換為例：1.如 2.3 介紹，首先確定四參數坐標轉換的轉換公式為（2.9） ，為了求解公式（2.9）里的轉換參數，又需要用到公式（2.7） ，同時為了保證精度，還需要用最小二乘法進行平差。2.有了轉換的公式和基礎，接下來就是確定在哪 2 個橢球之間進行轉換，從而就知道了橢球的長半徑與偏心率。3.公式的電算化，用 VB 語言把公式進行簡化，能通過計算機計算。（假設已知 4 組公共點）4.數據的輸入，通過

35、 VB 的可視化界面，很清楚的把需要的數據輸入程序。5.數據的輸出，同樣，VB 的界面可以直接顯示我們求出去的結果。143.2 程序界面及使用方法3.2.1 大地坐標和空間直角坐標的轉換1、.程序界面2、 說明在運行的對話框中，需要輸入橢球的長半徑 a 和橢球的第一偏心率，這 2 個數據可以根據資料查詢所得。當進行大地坐標轉換為直角坐標的時候，需要輸入 B，L 他們的單位為度。分秒，由于在三角函數中使用的是弧度，我們就需要把 B，L 轉換為弧度；然后就是依據公式進行數據的計算，在接下來的文本框把計算結果輸出。當進行直角坐標轉換為大地坐標的時候，基本過程也是數據讀取、計算、15輸出，在輸出的時候

36、需要把弧度轉換為角度表示方法。3、 運行結果表 3.1 結果驗證以 WGS 84 橢球參數為例：長半徑 a=6378137m 第一偏心率 e= 0.0818191902997090大地坐標轉換為空間直角坐標B50 度 32 分 28秒X115026.591809885L88 度 22 分 38秒Y4060183.69759044H376mZ4901550.48730854空間直角坐標轉換為大地坐標X120000B50 度 16 分 3.23 秒Y4100000L88 度 19 分 24.70 秒Z4902000H26180.7320015852m結論：此運算結果經過檢核與比較，XYZ 的轉換精

37、度在 0.0001m 以上，角度的轉換精度在 0.001 秒以上。3.2.2 高斯坐標正算1、 程序界面162、說明在運行的對話框中，需要輸入橢球的長半徑 a 和橢球的第一偏心率，這 2 個數據可以根據資料查詢所得。 然后輸入需要進行計算的數據，分別以度分秒輸入對應的文本框，點擊轉換，即可完成計算，在后面的文本框輸出數據。3、運行結果表 3.2 結果驗證以克氏橢球為例：長半徑 a=6378245，第一偏心率e=0.081813334016931494319BLXY轉換前32 度 24 分57.6522 秒118 度 54分 15.2206秒轉換后3589644.285190420679136.

38、43821917結論：此運算結果參考文獻【6】 ，進行比較，結果是正確的。3.2.3 高斯坐標反算1、程序界面2、運行結果表 3.3 結果驗證以克氏橢球為例：長半徑 a=6378245，第一偏心率e=0.081813334016931494319XYBL轉換前3589644.285190420679136.4382196轉換后32 度 24 分57.652180864秒118 度 54 分15.220600291 秒結論：此運算結果參照參考文獻【6】 ，進行比較，結果是正確的。 183.2.4 四參數轉換1、程序界面：2、 說明：根據界面的表示，假設已有 4 個公共點，首先拿 2 個公共點在第

39、一步進行運算，求出一組近似的轉換參數；第二步就是用最小二乘法對轉換參數進行平差，求出精確的轉換參數；第三步就是轉換需要轉換的點。19第四章 程序代碼4.1 大地坐標和空間直角坐標的轉換Private Sub Command1_Click()Dim B1 As Double, L1 As Double, H1 As DoubleDim X1 As Double, Y1 As Double, Z1 As DoubleDim B11 As Double, B12 As Double, B13 As DoubleDim L11 As Double, L12 As Double, B13 As Doub

40、leDim N As Double, e As Double, a As DoubleB11 = Text1.Text: B12 = Text19.Text: B13 = Text18.TextL11 = Text2.Text: L12 = Text17.Text: L13 = Text16.Text20H1 = Text3.Texta = Text13.Texte = Text14.TextB1 = B11 + B12 / 60 + B13 / 3600 把度分秒轉換為度來表示L1 = L11 + L12 / 60 + L13 / 3600B1 = B1 * 3.14159265358979

41、 / 180 把角度轉換為弧度L1 = L1 * 3.14159265358979 / 180N = a / Sqr(1 - e * e * Sin(B1) * Sin(B1)X1 = (N + H1) * Cos(B1) * Cos(L1)Y1 = (N + H1) * Cos(B1) * Sin(L1)Z1 = (N * (1 - e * e) + H1) * Sin(B1)Text4.Text = Str(X1)Text5.Text = Str(Y1)Text6.Text = Str(Z1)End SubPrivate Sub Command2_Click()Dim B2 As Doub

42、le, L2 As Double, H2 As DoubleDim X2 As Double, Y2 As Double, Z2 As DoubleDim B21 As Double, B22 As Double, B23 As DoubleDim L21 As Double, L22 As Double, L23 As DoubleDim N As Double, e As Double, a As Double, j1 As DoubleDim R As Double, j2 As Double, C As Double, B0 As Double, dB As DoubleX2 = Te

43、xt9.TextY2 = Text10.TextZ2 = Text7.Texta = Text13.Text21e = Text14.Textj1 = Atn(Z2 / Sqr(X2 * X2 + Y2 * Y2) 計算大地經度 BR = Sqr(X2 * X2 + Y2 * Y2 + Z2 * Z2)j2 = Atn(Z2 / R) / Sqr(-(Z2 / R) * (Z2 / R) + 1)C = a * e * e / (2 * R * Sqr(1 - e * e * Sin(j2) * Sin(j2)dB = C * Sin(j2) * Sin(j2) * (1 + 2 * C *

44、Cos(2 * j2)B0 = j2 + dB Do B2 = B0 B0 = Atn(Z2 / Sqr(X2 * X2 + Y2 * Y2) * (1 + a * e * e * Sin(B2) / (Z2 * Sqr(1 - e * e * Sin(B2) * Sin(B2) Loop Until Abs(B0 - B2) 0, 1, 0) 得到帶號L1 = L1 - 6 * ZoneNum + 3 把精度轉換到當前帶23B1 = B1 * 3.141592653 / 180 把度轉換成弧度L1 = L1 * 3.141592653 / 180t = Tan(B1)yi = e1 * Co

45、s(B1)N = a / Sqr(1 - e2 * Sin(B1) * Sin(B1)m = Cos(B1) * L1X = a * (1 - e2) * (1 + 3 * e2 / 4 + 45 * e2 * e2 / 64 + 175 * e2 * e2 * e2 / 256) * B1 _ - (3 * e2 / 4 + 15 * e2 * e2 / 16 + 525 * e2 * e2 * e2 / 512) * Sin(2 * B1) / 2 _ + (15 * e2 * e2 / 64 + 105 * e2 * e2 * e2 / 256) * Sin(4 * B1) / 4 _

46、- (35 * e2 * e2 * e2 / 512) * Sin(6 * B1) / 6)X1 = X + N * t * (1 / 2 + (5 - t * t + 9 * yi * yi + 4 * yi * yi * yi * yi) / 24 + (61 - 58 * t * t + t * t * t * t) / 720 * m * m) * m * m) * m * m)Y1 = N * (1 + (1 - t * t + yi * yi) / 6 + (5 - 18 * t * t + t * t * t * t + 14 * yi * yi - 58 * yi * yi *

47、 t * t) / 120 * m * m) * m * m) * m)Y1 = Y1 + 500000Y1 = 1000000 * ZoneNum + Y1Text3.Text = Str(X1)Text4.Text = Str(Y1)End Sub4.3 高斯坐標轉換的反算Private Sub Command2_Click()Dim X2 As Double, Y2 As Double, B2 As Double, L2 As DoubleDim a As Double, e As Double, e1 As Double, yi As Double, t As DoubleDim B2

48、1 As Double, B22 As Double, B23 As Double24Dim L21 As Double, L22 As Double, L23 As DoubleDim iZone As IntegerDim Bf As DoubleDim sBf As DoubleDim Bfi As DoubleDim x0 As Double Dim V As DoubleDim s As DoubleDim yn As DoubleX2 = Text5.TextY2 = Text6.Texta = Text9.Texte = Text10.Texte1 = Sqr(e * e / (

49、1 - e * e)iZone = Int(Y2 / 1000000)Y2 = Y2 - iZone * 1000000 帶號設置Y2 = Y2 - 500000 減 500 公里常量x0 = 111134.8611Bf = X2 / x0While (Abs(Bfi - Bf) 0.00001) Bfi = Bf sBf = Sin(Bfi * 3.141592653 / 180) Bf = (X2 + (32005.7799 * sBf + 133.9238 * sBf * sBf * sBf + _ 0.6973 * sBf * sBf * sBf * sBf * sBf + _ 0.0

50、039 * sBf * sBf * sBf * sBf * sBf * sBf * sBf) * _ Cos(Bf * 3.141592 / 180) / x0WendBf = Bf * 3.141592653 / 180t = Tan(Bf)s = Sin(Bf)25yi = e1 * Cos(Bf)V = Sqr(1 + e1 * e1 * Cos(Bf) * Cos(Bf)yn = Y2 * Sqr(1 - e * e * s * s) / aB2 = Bf * 180 / 3.141592653 - V * V * t * 180 / 3.141592 / 2 * (yn * yn -

51、 _ (5 + 3 * t * t + yi * yi - 9 * yi * yi * t * t) * yn * yn * yn * yn / 12 + _ (61 + 90 * t * t + 45 * t * t * t * t) * yn * yn * yn * yn * yn * yn / 360) L2 = 180 / 3.141592653 / Cos(Bf) * (yn - (1 + 2 * t * t + yi * yi) * yn * yn * yn / 6 + _ (5 + 28 * t * t + 24 * t * t * t * t + 6 * yi * yi + 8

52、 * yi * yi * t * t) _ * yn * yn * yn * yn * yn / 120) L2 = L2 + iZone * 6 - 3B21 = Int(B2) 把以度表示轉換為度分秒表示B22 = Int(B2 - B21) * 60)B23 = (B2 - B21) * 60 - B22) * 60L21 = Int(L2)L22 = Int(L2 - L21) * 60)L23 = (L2 - L21) * 60 - L22) * 60Text7.Text = Str(B21): Text18.Text = Str(B22): Text16.Text = Str(B2

53、3)Text8.Text = Str(L21): Text17.Text = Str(L22): Text15.Text = Str(L23)End SubPrivate Sub Option1_Click()If Option1.Value = True Then Frame2.Visible = False: Frame1.Visible = TrueEnd SubPrivate Sub Option2_Click()If Option2.Value = True Then Frame1.Visible = False: Frame2.Visible = True26End Sub4.4

54、四參數轉換Option Base 1Option ExplicitFunction MCinv(N As Integer, S1() As Double, S2() As Double) As Boolean 矩陣求逆 Dim x As Integer, y As Integer Dim SS(7, 7) As Double For x = 1 To N For y = 1 To N SS(x, y) = S1(x, y) Next y Next x ReDim nIs(N) As Integer, nJs(N) As Integer Dim i As Integer, j As Intege

55、r, K As Integer Dim D As Double, p As Double, s As Double, t As Double, q As Double, B As Double For K = 1 To N D = 0# For i = K To N For j = K To N p = S1(i, j) * S1(i, j) + S2(i, j) * S2(i, j) If (p D) Then D = p nIs(K) = i nJs(K) = j End If Next j Next i27 If (D + 1# = 1#) Then MCinv = False Exit

56、 Function End If If (nIs(K) K) Then For j = 1 To N t = S1(K, j) S1(K, j) = S1(nIs(K), j) S1(nIs(K), j) = t t = S2(K, j) S2(K, j) = S2(nIs(K), j) S2(nIs(K), j) = t Next j End If If (nJs(K) K) Then For i = 1 To N t = S1(i, K) S1(i, K) = S1(i, nJs(K) S1(i, nJs(K) = t t = S2(i, K) S2(i, K) = S2(i, nJs(K

57、) S2(i, nJs(K) = t Next i End If S1(K, K) = S1(K, K) / D S2(K, K) = -S2(K, K) / D For j = 1 To N If (j K) Then28 p = S1(K, j) * S1(K, K) q = S2(K, j) * S2(K, K) s = (S1(K, j) + S2(K, j) * (S1(K, K) + S2(K, K) S1(K, j) = p - q S2(K, j) = s - p - q End If Next j For i = 1 To N If (i K) Then For j = 1

58、To N If (j K) Then p = S1(K, j) * S1(i, K) q = S2(K, j) * S2(i, K) s = (S1(K, j) + S2(K, j) * (S1(i, K) + S2(i, K) t = p - q B = s - p - q S1(i, j) = S1(i, j) - t S2(i, j) = S2(i, j) - B End If Next j End If Next i For i = 1 To N If (i K) Then p = S1(i, K) * S1(K, K) q = S2(i, K) * S2(K, K) s = (S1(

59、i, K) + S2(i, K) * (S1(K, K) + S2(K, K) S1(i, K) = q - p S2(i, K) = p + q - s End If29 Next i Next K For K = N To 1 Step -1 If (nJs(K) K) Then For j = 1 To N t = S1(K, j) S1(K, j) = S1(nJs(K), j) S1(nJs(K), j) = t t = S2(K, j) S2(K, j) = S2(nJs(K), j) S2(nJs(K), j) = t Next j End If If (nIs(K) K) Th

60、en For i = 1 To N t = S1(i, K) S1(i, K) = S1(i, nIs(K) S1(i, nIs(K) = t t = S2(i, K) S2(i, K) = S2(i, nIs(K) S2(i, nIs(K) = t Next i End If Next K MCinv = True For x = 1 To N For y = 1 To N S2(x, y) = S1(x, y) Next y Next x For x = 1 To N30 For y = 1 To N S1(x, y) = SS(x, y) Next y Next xEnd Functio