坐標轉(zhuǎn)換參數(shù)求取及坐標轉(zhuǎn)換程序設計

1、畢業(yè)設計設計題目 坐標轉(zhuǎn)換參數(shù)求取及坐標轉(zhuǎn)換程序設計 學生姓名 張威 指導教師 杜繼亮 專 業(yè) 測繪工程 班 級 測繪12-2班 填寫日期 2016/4/6 礦業(yè)工程學院摘 要坐標系統(tǒng)是測量工作中定位的基礎，坐標系統(tǒng)有多種形式和基準，由于各測量工作目的不同，所選用的坐標基準也會不同，根據(jù)不同的工作要求需要將不同的坐標系下的坐標進行相互轉(zhuǎn)換。在這些坐標轉(zhuǎn)換的過程中會用到很多坐標轉(zhuǎn)換模型，但是坐標系轉(zhuǎn)換模型過于復雜手算非常困難。本設計為了方便施工時遇到的坐標轉(zhuǎn)換問題，設計利用Visual Basic 6.0編程語言編寫程序，用來實現(xiàn)坐標系統(tǒng)之間的轉(zhuǎn)換以及轉(zhuǎn)換參數(shù)的求解，例如:大地坐標與空間直角坐標

2、的相互轉(zhuǎn)換、高斯投影正反算、二維坐標轉(zhuǎn)換與四參數(shù)計算、三維坐標轉(zhuǎn)換與七參數(shù)轉(zhuǎn)換、同參考基準下的坐標換帶計算，以及坐標數(shù)據(jù)的批量處理。關鍵字：坐標系統(tǒng)，轉(zhuǎn)換模型，坐標轉(zhuǎn)換，程序設計AbstractThe base of coordinate system in surveying work. there are many forms and benchmarks in the coordinate system. However, in general engineering, the control point and coordinate. System are the same. So I

3、t is necessary to transform the control point. coordinate during the construction process. Due to different purposes of each measurement and the selected. different coordinate references, there will be many different coordinate systems. Coordinate systems used in the measurement work are as follows:

4、 WGS-84 World Geodetic System, China Geodetic Coordinate System 2000, National Geodetic Coordinate System 1980, Beijing coordinate system 1954 and Local Coordinate System. There are space rectangular coordinate, geodetic coordinate and plane coordinate in the way of the reference in the same coordin

5、ate. According to the requirements of different tasks, we need to convert coordinates under the different coordinate systems. On condition that the coordinates of the reference standard can be obtained. the normal construction work can be done. A lot of coordinate transformation models are used in t

6、he process of the coordinate transformation. But the coordinate transformation model is very complex and difficult. Nowadays the conversion formula is suitable for the computerization whose language is easy to learn. So in the design I make use of Visual Basic 6 programming language to realize the t

7、ransformation between the coordinate system and transformation parameters.Key words : coordinate systems transformation model coordinate transform programming目錄摘 要IAbstractII第1章 緒論11.1研究背景和意義11.2國內(nèi)外研究現(xiàn)狀11.3研究的主要內(nèi)容21.4程序設計思路方法2第2章 坐標理論和知識32.1地球橢球32.2基準32.3測量常用坐標系42.3.1大地坐標系42.3.2空間直角坐標系42.3.3平面坐標系52.

8、3.4地方獨立坐標系52.4我國常用坐標系62.4.1 1954年北京坐標系62.4.2 1980國家大地坐標系62.4.3 WGS-84世界大地坐標系62.4.4 2000國家大地坐標系7第3章 坐標轉(zhuǎn)換程序設計83.1大地坐標系與空間直角坐標系的轉(zhuǎn)換93.1.1大地坐標轉(zhuǎn)換成空間直角坐標93.1.2空間直角坐標轉(zhuǎn)成換大地坐標103.2高斯平面坐標與大地坐標的轉(zhuǎn)換133.2.1 高斯正算133.2.2 高斯反算143.3 高斯投影換帶計算173.4 二維坐標轉(zhuǎn)換與四參數(shù)計算203.5 三維坐標轉(zhuǎn)換與七參數(shù)計算233.6 三、六度帶帶號與中央子午線計算27第4章 設計總結(jié)31致謝32參考文獻3

9、3第1章 緒論1.1研究背景和意義隨著大地測量學，攝影測量學的發(fā)展和電子計算機的普及，對各種坐標系的研究變得越來越重要了1。隨著現(xiàn)在社會的快速發(fā)展，各種各樣的大型工程的建設，凡是工程施工就必定需要坐標來定位才可以，而建造的地方又不同，又沒有一個滿足全部地方的坐標系統(tǒng)，所以產(chǎn)生多種不同的坐標系統(tǒng)，實際工作中測量人員必定會按照實際情況選擇最為合適的坐標系統(tǒng)。而坐標系統(tǒng)之間的轉(zhuǎn)換比較復雜，手算工作量巨大，因而各種坐標轉(zhuǎn)換模型相繼出現(xiàn)，利用計算機強大的數(shù)據(jù)計算能力可以輕松應對這些問題，提高工作效率。坐標轉(zhuǎn)換的意義重大，不僅在我們熟知的工程領域中，在國防建設、航空航天科技、城市匯劃等眾多領域中都發(fā)揮著重

10、要的作用，可以說對社會進步有著必不可少的作用。1.2國內(nèi)外研究現(xiàn)狀自60年代以來，各國大地測量學者，經(jīng)過大量研究，提出了多種坐標轉(zhuǎn)換模型及多種解算方法，北美1927基準面（基于克拉克1966橢球體與北美1983基準面（基于GRS1980橢球體）之間坐標轉(zhuǎn)換是根據(jù)研究區(qū)內(nèi)一系列己知點的大地坐標或網(wǎng)格坐標改正量進行插值進行的坐標系轉(zhuǎn)換；英國采用北向與東向的雙線性網(wǎng)格插值進行坐標轉(zhuǎn)換；挪威在海岸帶調(diào)查中，采用經(jīng)緯度多項式用于坐標系轉(zhuǎn)換這種方法進行新（ED87一歐洲1987基準面）、舊（ED50一歐洲1950基準面）坐標系之間的轉(zhuǎn)換：歐洲石油勘探組織（EPSG）對新、舊坐標系采用“雙線性插值” 進行坐

11、標轉(zhuǎn)換2。國內(nèi)空間三維直角坐標轉(zhuǎn)換中，一般采用布爾莎七參數(shù)模型。一般有7個轉(zhuǎn)換參數(shù)，即3個平移參數(shù)，3個旋轉(zhuǎn)參數(shù)和1個尺度參數(shù)。需要三個及已經(jīng)公共點時，才能利用平差的方法求出七參數(shù)。1.3研究的主要內(nèi)容本坐標轉(zhuǎn)換程序可實現(xiàn)功能有:1、大地坐標與空間直角坐標的相互轉(zhuǎn)換， 2、高斯投影正反算，3、二維坐標轉(zhuǎn)換與四參數(shù)計算，4、三維坐標轉(zhuǎn)換與七參數(shù)轉(zhuǎn)換，5、同參考基準下的坐標換帶計算，以及坐標數(shù)據(jù)的批量處理。1.4程序設計思路方法本程序名為萬能坐標轉(zhuǎn)換器。設計前期收集相關資料，參考一些成熟的坐標轉(zhuǎn)換軟件，確定程序應有的功能以及界面設計。運用VB編寫程序時，查閱相關書籍獲取理論知識以及轉(zhuǎn)換模型。完成程

12、序后將已知正確數(shù)據(jù)帶入其中驗證程序結(jié)果是否正確。若出現(xiàn)錯誤則檢查每步代碼，直到程序完美運行為止。第2章 基礎知識準備2.1地球橢球地球橢球體又稱地球橢圓體或地球扁球體，代表地球大小和形狀的數(shù)學曲面，以長半徑和扁率表示，因它十分迫近于橢球體，故通常以參考橢球體表示地球橢球體的形狀和大小。通常所說地球的形狀和大小，實際上就是以參考橢球體的半長徑、半短徑和扁率來表示。1975年國際大地測量與地球物理聯(lián)合會推薦的數(shù)據(jù)為：半長徑6378140米，半短徑6356755米，扁率1298.257。在眾多橢球體中，WGS-84橢球體被認為符合上述條件最好的橢球。2.2基準所謂基準是指為描述空間位置而定義的點、線

13、和面。而大地測量基準是指用以描述地球形狀的地球橢球參數(shù)，包含描述地球橢球幾何特征的長短半軸和物理特征的有關參數(shù)、地球在空間的定位及定向以及描述這些位置所采用的單位長度的定義。不同的坐標系統(tǒng)會使用的基準也不同，根據(jù)參考橢球所選原點位置不同，可以分為地心坐標系和參心坐標系。地心坐標系是以地球的質(zhì)心為原點，有地心大地坐標系和地心空間直角坐標系兩種表述方法。地心空間直角坐標系的定義為：以地球質(zhì)心為原點，X軸指向格林尼治子午面與地球赤道的交點，Z軸指向北極，Y軸過原點垂直于平面 XOZ，構(gòu)成右手空間直角坐標系。地心大地坐標系定義為：以地球的質(zhì)心作為原點，以地球自轉(zhuǎn)軸作為橢球的短軸，大地緯度B是過地面點的

14、橢球法線與橢球赤道面之間的夾角，大地經(jīng)度L為過地面點的橢球子午面與格林尼治子午面之間的夾角，大地高度H為地面點沿橢球法線到橢球面的最短距離。WGS-84坐標系，CGCS2000坐標系，GLONASS是采用PZ-90坐標，都是屬于地心坐標系。參心坐標系是選取一個參考橢球面作為基本的參考面，選一參考點作為大地測量的起算點，從而確定參考橢球在地球面的位置和方向。這時參考橢球的原點不會和地球質(zhì)心重合，所以稱為參心。北京54坐標系、西安80坐標系和新北京54坐標系，都是參心坐標系。它同樣具有參心大地坐標系和參心空間直角坐標系兩種表述方法，它們的定義與地心坐標系的定義相似。2.3測量常用坐標系2.3.1大

15、地坐標系空間一點的大地坐標用大地經(jīng)度L、大地緯度B和大地高度H表示，地面上P地點的大地子午面NPS與起始大地子午面所構(gòu)成的二面角L稱P地點的大地經(jīng)度， P地點對于橢球的法線與赤道面的夾角B稱P地點的大地緯度。如圖2-1所示圖2-1大地坐標系P地點沿法線到橢球面的距離H稱大地高，從橢球面起算，向外為正，向內(nèi)為負3。H = H正常 + （高程異常）H = H正 + N（大地水準面差距） 2.3.2空間直角坐標系空間直角坐標系的坐標原點與參考橢球的中心重合，Z軸正向指向參考橢球的北極，X軸正向指向起始子午面與赤道的交點，Y軸按右手系與X軸呈90夾角且位于赤道面上，某點在空間中的坐標可用該點在此空間坐

16、標系的各個坐標軸上的投影來表示4。如圖2-2所示：圖2-2空間直角坐標系2.3.3平面坐標系平面直角坐標系是利用投影，將空間坐標通過某種數(shù)學變換映射到平面上，這種變換稱為投影變換5。在我國一般采用的是高斯一克呂格投影，是目前測量上廣泛采用的正形投影，特點是沒有角度變形，在不同點上的長度比隨點位而異，但在同一點上各方向的長度比相同，也稱為高斯投影6。2.3.4地方獨立坐標系在我國平面坐標主要采用的是高斯投影，在該投影中，除中央子午線外，其它位置上的任何線段，投影后都會產(chǎn)生一定的長度變形，而且變形隨離開中央子午線的距離增加而增加7。因此一般采用分帶投影的辦法，來限制長度變形，我國規(guī)定了采用3度帶或

17、6度帶進行分帶投影。在城市、工礦等工程測量中，如果直接在國家分帶坐標系中建立控制網(wǎng)，會使地面長度投影的變形較大，當長度變形大于25 cm/km時，就難以滿足工程上的需要8。另一些特殊的測量，比如大橋施工測量，水利水壩測量，滑坡變形監(jiān)測等，采用國家坐標系精度不能滿足工程要求，所以常常會建立適合本地區(qū)的地方獨立坐標系9。2.4我國常用坐標系2.4.1 1954年北京坐標系1954年北京坐標系，是采取先將我國一等鎖與前蘇聯(lián)遠東一等鎖相聯(lián)接，然后以連接處呼瑪，吉拉林，東寧基線網(wǎng)擴大邊端點的前蘇聯(lián)1942年普爾科沃坐標系的坐標為起算數(shù)據(jù)，平差我國東北及東部一等鎖10。橢球參數(shù)：長半軸a=6378245m

18、，短半軸b=6356863.0188m，扁率=1/298.3，第一偏心率平方e2=0.00669342161454287，第二偏心率平方e2=0.00673852540614652。2.4.2 1980國家大地坐標系1980國家大地坐標系采用地球橢球基本參數(shù)為1975年國際大地測量與地球物理聯(lián)合會第十六屆大會推薦的數(shù)據(jù)，大地原點設在我國中部的陜西省涇陽縣永樂鎮(zhèn)，位于西安市西北方向約60公里。橢球參數(shù)：長半軸a=63781405m，短半軸b=6356755.2882m，扁率=1/298.257，第一偏心率平方e2 =0.00669438499959，第二偏心率平方e2=0.00673950181

19、947。2.4.3 WGS-84世界大地坐標系WGS-84坐標系是一種國際上采用的地心坐標系。坐標原點為地球質(zhì)心，其地心空間直角坐標系的Z軸指向BIH （國際時間服務機構(gòu)）1984.0定義的協(xié)議地球極（CTP）方向，X軸指向BIH 1984.0的零子午面和CTP赤道的交點，Y軸與Z軸、X軸垂直構(gòu)成右手坐標系，稱為1984年世界大地坐標系統(tǒng)。橢球參數(shù)：長半軸a=6378137，短半軸b=6356752.3142，扁率=1/298.2572236，第一偏心率平方e2=0.00669437999013 ，第二偏心率平方e2=0.006739496742227。2.4.4 2000國家大地坐標系CGC

20、S2000是（中國）2000國家大地坐標系的縮寫，該坐標系是通過中國GPS連續(xù)運行基準站、空間大地控制網(wǎng)以及天文大地網(wǎng)與空間地網(wǎng)聯(lián)合平差建立的地心大地坐標系統(tǒng)11。Z軸指向BIH1984.0定義的協(xié)議極地方向（BIH國際時間局），X軸指向BIH1984.0定義的零子午面與協(xié)議赤道的交點，Y軸按右手坐標系確定。橢球參數(shù)：長半軸a=6378137，短半軸b=6356752.31414，扁率=1/298.2572236，第一偏心率平方e2=0.00669437999013 ，第二偏心率平方e2=0.006739496742227。2.5坐標轉(zhuǎn)換模型2.5.1大地坐標系與空間直角坐標系轉(zhuǎn)換模型將同一坐

21、標系下的大地坐標（B，L，H）轉(zhuǎn)換成空間直角坐標（X，Y，Z）的轉(zhuǎn)換公式為： X=N+H*cosB*cosL Y=N+H*cosB*sinL Z=N*1-e2+H*sinB N=a1-e2sin2B式中：a為參考橢球長半軸，e為第一偏心率，N為卯酉圈的半徑12。 將同一坐標系下的空間直角坐標（X，Y，Z）轉(zhuǎn)換為大地坐標（B，L，H）的公式為: L=arctan(YX) B=arctan(Z+Ne2sinBX2+Y2) H=ZsinB-N(1-e2)用公式進行空間直角坐標轉(zhuǎn)換大地坐標時，需要采用迭代計算大地緯度B。具體計算時，可先根據(jù)下式求出大地緯度B的初值：B=arctan(ZX2+Y2)因為

22、 N=c1+e2cos2B , 1cos2B=1+tan2B 帶入B=arctan(Z+Ne2sinBX2+Y2) 得:tanB=ZX2+Y2+ce2tanBX2+Y2*1+e2+tan2B令:t=ZX2+Y2，p=ce2X2+Y2，k=1+e2則式子可寫成:tanB=t+p*tanBk+tan2B然后利用該初值B代入公式右端tanB中，將等式左邊的結(jié)果再次代入右端tanB，直到最后兩次B值之差小與允許誤差為止。當?shù)玫酱蟮鼐暥菳后，代入公式即可求出大地高H。2.5.2高斯正反算轉(zhuǎn)換模型得到了點的大地坐標(L，B)，就可以將其轉(zhuǎn)化為某投影帶的高斯平面坐標，我們將橢球參數(shù)代入高斯投影正算公式得到更

23、適用于電算的高斯坐標計算的實用公式： x=X+N2t*cos2B*l2+N24t5-t2+92+44cos4B*l4+ N720t61-58t2+t4cos6B*l6 y=N*cosB*l+N61-t2+2cos3B*l3+N120(5-18t2+t4+ 142-582t2)cos5B*l5 X=c0B+2cosB+4cos3B+6cos5B+8cos7BsinB 0=1-34e2+4564e4-175256e6+1102516384e8 2=0-1 4=4564e4-175256e6+1102516384e8 6=-175256e6+1102516384e8 8=3151024e8 N=a1

24、-e2sin2B t=tanB =ecosB l=(L-L0) 其計算結(jié)果的精度可達0. 001m。只要得到了高斯平面坐標(X, Y)后，便可通過高斯反算公式將其轉(zhuǎn)換成大地坐標(B, L)，高斯投影反算公式為: B=Bf-tf2MfNfy2+tf24MfNf35+3tf2+f2-9f2tf2y4-tf720MfNf5(61+90tf2+45tf4)y6 L=1NfcosBfy-16Nf3cosBf1+2tf2+f2y3+1120Nf5cosBf(5+28tf2+24tf4)y5 Bf=Xc-2cosB+4cos3B+6cos5B+8cos7BsinB0 tf=tanBf Mf=a(1-e2)2

25、(1-e2sin2Bf)3 Nf=a1-e2sin2Bf f=tanBf 它們的計算精度,即平面坐標可達0.0001。2.5.3坐標轉(zhuǎn)換與參數(shù)計算轉(zhuǎn)換模型二維轉(zhuǎn)換模型：2個平移參數(shù)x y （原點不重合產(chǎn)生的）1個旋轉(zhuǎn)參數(shù)（坐標軸不平行產(chǎn)生的）1個尺度參數(shù)m（兩坐標系間的尺度不一致產(chǎn)生的）x2y2=1+mcossin-sincosx1y1+xy(x1,y1)為某點在A空間直角坐標系中的坐標。(x2,y2)為某點在B空間直角坐標系中的坐標。(x,y)為某點從A空間直角坐標系轉(zhuǎn)換到B空間直角坐標系的兩個平移參數(shù)。為從A空間直角坐標系轉(zhuǎn)換到B空間直角坐標系中標系的兩個平移參數(shù)。為從A空間直角坐標系轉(zhuǎn)換

26、到B空間直角坐標系中一個旋轉(zhuǎn)參數(shù)。m為從A空間直角坐標系轉(zhuǎn)換到B空間直角坐標系中的一個尺度參數(shù)。平面四參數(shù)求解步驟如下:利用公共點計算坐標參數(shù)，但至少有兩個公共點，當有i個公共點時，可利用最小二乘原理求解參數(shù)。將B直角坐標系中的坐標視為觀測值，設A直角坐標系下的坐標視為無誤差，列誤差方程為：vx1vy1vxivyiB=10100101y1-x1yi-xix1y1xiyixym-X1Y1XiYiB-X1Y1XiYiA寫成矩陣形式既：V=BX-L由于各點的坐標可視為同精度獨立觀測值，因此P=I。把各點坐標已知值帶入上述誤差方程，然后按下列公式求解出四參數(shù)：X=(BTB)-1(BTL)三維轉(zhuǎn)換模型：

27、3個平移參數(shù)x y z（原點不重合產(chǎn)生的）3個旋轉(zhuǎn)參數(shù)x y z（坐標軸不平行產(chǎn)生的）1個尺度參數(shù)m（兩坐標系間的尺度不一致產(chǎn)生的） x2y2z2=1+mR(x)R(y)R(z)x1y1z1+xyzRx=1 000cosxsinx0-sinxcosxRy=cosy0-siny010sinycosy0Rz=coszsinz0-sinzcosz0001一般x，x，x為微小轉(zhuǎn)角，模型可簡化為：x2y2z2=1+m1z-y-z1xy-x1x1y1z1+xyz(x1,y1,z1)為某點在A空間直角坐標系中的坐標。(x2,y2z2)為某點在B空間直角坐標系中的坐標。(x,y,z)為某點從A空間直角坐標系轉(zhuǎn)

28、換到B空間直角坐標系的三個平移參數(shù)。(x，y，z)為從A空間直角坐標系轉(zhuǎn)換到B空間直角坐標系中三個旋轉(zhuǎn)參數(shù)。m為從A空間直角坐標系轉(zhuǎn)換到B空間直角坐標系中的一個尺度參數(shù)。七參數(shù)求解步驟如下:利用公共點計算坐標參數(shù)，但至少有三個公共點，當有i個公共點時，可利用最小二乘原理求解參數(shù)。將B空間直角坐標系中的坐標視為觀測值，設A空間直角坐標系下的坐標視為無誤差，列誤差方程為：vx1vy1vz1vxivyivziB=1001000100100010010Z1-Y10Zi-Yi-Z10X1-Zi0Xi-Y1-X10-Yi-Xi0X1Y1Z1XiYiZixyzxyzm-X1Y1Z1XiYiZiB-X1Y1Z

29、1XiYiZiA寫成矩陣形式既：V=BX-L由于各點的坐標可視為同精度獨立觀測值，因此P=I。把各點坐標已知值帶入上述誤差方程，然后按下列公式求解出七參數(shù)：X=(BTB)-1(BTL)2.5.4三、六度帶帶號與中央子午線計算模型求帶號、中央子午線公式為：3帶:帶號N=IntL-1.53+1 中央子午線L0=N*3 L為當?shù)亟?jīng)度6帶:帶號N=IntL6+1 中央子午線L0=N*6-3 L為當?shù)亟?jīng)度第3章 坐標轉(zhuǎn)換程序設計 計劃程序有哪些功能 程序界面設計 搜集坐標轉(zhuǎn)換公式 運用VB編寫程序 檢查能否正常運行 以及數(shù)據(jù)正確性 否 是完成程序圖3-1 編寫程序流程圖圖3-2程序主界面3.1大地坐標系

30、與空間直角坐標系的轉(zhuǎn)換主要代碼如下：橢球參數(shù)設置If Combo1.Text = WGS-84 Then a1 = 6378137: b1 = 6356752.3142 a1為長半軸，b1為短半軸 Text7.Text = a1: Text8.Text = b1ElseIf Combo1.Text = CGCS2000 Then a1 = 6378137: b1 = 6356752.31414 Text7.Text = a1: Text8.Text = b1ElseIf Combo1.Text = 北京54 Then a1 = 6378245: b1 = 6356863.0188 Text7.

31、Text = a1: Text8.Text = b1ElseIf Combo1.Text = 西安80 Then a1 = 6378140: b1 = 6356755.2882 Text7.Text = a1: Text8.Text = b1ElseIf Combo1.Text = 自定義 Then Text7.Enabled = True: Text8.Enabled = True a1 = Text7.Text: b1 = Text8.TextEnd If大地坐標轉(zhuǎn)空間直角坐標主要代碼 B = Text1.Text: L = Text2.Text: H = Text3.Text e12 =

32、 (a1 2 - b1 2) / a1 2 N = a1 / Sqr(1 - e12 * (Sin(DuToHu(B) 2) x = (N + H) * Cos(DuToHu(B) * Cos(DuToHu(L) y = (N + H) * Cos(DuToHu(B) * Sin(DuToHu(L) z = (N * (1 - e12) + H) * Sin(DuToHu(B) Text4.Text = x: Text5.Text = y: Text6.Text = z空間直角坐標轉(zhuǎn)大地坐標主要代碼 x = Text4.Text: y = Text5.Text: z = Text6.Text

33、L = HuToDu(Atn(y / x) 計算經(jīng)度 If y / x 0 Then L = HuToDu(Atn(y / x) + 3.14159265358979) End If e12 = (a1 2 - b1 2) / a1 2 第一偏心率平方 TB = z / Sqr(x 2 + y 2) 初始緯度B正切值 c1 = (a1 2) / b1 t = z / Sqr(x 2 + y 2) p = (c1 * e12) / Sqr(x 2 + y 2) k = 1 + (a1 2 - b1 2) / b1 2 TB1 = t + (p * TB) / Sqr(k + TB 2) 用Tan

34、B作為未知量 Do Until Abs(TB1 - TB) 0.000000000001 迭代法 TB = TB1 TB1 = t + (p * TB) / Sqr(k + TB 2) Loop B = Atn(TB) N = a1 / Sqr(1 - e12 * (Sin(B) 2) 卯酉圈曲率半徑 H = z / Sin(B) - N * (1 - e12) 大地高 Text1.Text = HuToDu(B): Text2.Text = L: Text3.Text = H3.2高斯平面坐標與大地坐標的轉(zhuǎn)換高斯平面坐標(X, Y)與大地坐標(B ,L)的相互關系式分為兩類：第一類稱高斯投影

35、正算公式，由(B,L)求(X,Y)；第二類稱高斯投影反算公式，由(X,Y)求(B,L) 13。由于電子計算機和各種可編程序電子計算器在測量上廣泛應用，因此在這里給出適合電算的高斯正反算公式。主要代碼如下：正算Dim B, L, x, y As DoubleDim N, a1, b1, c, e12, e22, e1, e2, ll, l0, nn, t As DoubleDim XX, bb0, bb2, bb4, bb6, bb8 As DoubleDim yj As DoubleDim mf, nf, nnf, tf, bf, bf1 As DoubleB = Text1.Text: L

36、= Text2.Textl0 = Text5.Text 輸入中央子午線經(jīng)度yj = Text8.Text Y坐標加常數(shù)a1 = Text6.Text: b1 = Text7.Text 橢球長短半軸e12 = (a1 2 - b1 2) / a1 2 第一偏心率平方e22 = (a1 2 - b1 2) / b1 2 第二偏心率平方e1 = Sqr(a1 2 - b1 2) / a1 2) 第一偏心率e2 = Sqr(a1 2 - b1 2) / b1 2) 第二偏心率c = a1 / Sqr(1 - e12) 極點子午線曲率半徑bb0 = 1 - (3 / 4) * e22 + (45 / 6

37、4) * e22 2 - (175 / 256) * e22 3 + (11025 / 16384) * e22 4bb2 = bb0 - 1bb4 = (15 / 32) * e22 2 - (175 / 384) * e22 3 + (3675 / 8192) * e22 4bb6 = -(35 / 96) * e22 3 + (735 / 2048) * e22 4bb8 = (315 / 1024) * e22 4 ll = (Du(L) - l0) / 57.295779513 N = a1 / Sqr(1 - e12 * (Sin(DuToHu(B) 2) 卯酉圈曲率半徑 t =

38、Sin(DuToHu(B) / Cos(DuToHu(B) nn = e2 * Cos(DuToHu(B) XX = c * (bb0 * DuToHu(B) + (bb2 * Cos(DuToHu(B) + bb4 * (Cos(DuToHu(B) 3) + bb6 * (Cos(DuToHu(B) 5) + bb8 * (Cos(DuToHu(B) 7) * Sin(DuToHu(B) 子午線弧長 x = XX + (N / 2) * Sin(DuToHu(B) * Cos(DuToHu(B) * (ll 2) + (N / 24) * Sin(DuToHu(B) * (Cos(DuToH

39、u(B) 3) * (ll 4) * (5 - t 2 + 9 * (nn 2) + 4 * (nn 4) + (N / 720) * Sin(DuToHu(B) *(Cos(DuToHu(B) 5) * (ll 6) * (61 - 58 * (t 2) + (t 4) y = N * Cos(DuToHu(B) * ll + (N / 6) * (Cos(DuToHu(B) 3) * (ll 3) * (1 - t 2 + nn 2) + (N / 120) * (Cos(DuToHu(B) 5) * (ll 5) * (5 - 18 * (t 2) + t 4 + 14 * (nn 2)

40、 - 58 * (t 2) * (nn 2)反算x = Text3.Text: y = Text4.Texty = y - yjbf = x / (c * bb0) 迭代Do bf1 = (x / c) - (bb2 * Cos(bf) + bb4 * (Cos(bf) 3) + bb6 * (Cos(bf) 5) + bb8 * (Cos(bf) 7) * Sin(bf) / bb0Loop Until Abs(bf - bf1) 0.000000001mf = (a1 * (1 - e12) / (Sqr(1 - e12 * (Sin(bf) 2) 3)tf = Sin(bf) / Cos

41、(bf)nnf = e2 * Cos(bf)nf = a1 / Sqr(1 - e12 * (Sin(bf) 2)B = bf - (tf * (y 2) / (2 * mf * nf) + (tf * (y 4) * (5 + 3 * (tf 2) + (nnf 2) - 9 * (nnf 2) * (tf 2) / (24 * mf * (nf 3) - (tf * (y 6) * (61 + 90 * (tf 2) + 45 * (tf 4) / (720 * mf * (nf 5)L = y / (nf * Cos(bf) - (1 + 2 * tf 2 + nnf 2) * (y 3

42、) / (6 * nf 3 * Cos(bf) + (5 + 28 * (tf 2) + 24 * (tf 4) + 6 * (nnf 2) + 8 * (nnf 2) * (tf 2) * (y 5) / (120 * (nf 5) * Cos(bf)3.3 高斯投影換帶計算高斯投影雖然保證了角度沒有變形，但是長度變形比較嚴重，為了限制高斯投影的長度變形，必須以中央子午線進行分帶，把投影范圍限制在中央子午線東、西側(cè)一定的范圍內(nèi)；但是這樣又使得統(tǒng)一的坐標系分割成各帶的獨立坐標系14。這樣就會產(chǎn)生新的問題，需要鄰帶換算來解決。主要代碼如下：l1 = Text3.Text : yy1 = Text

43、4.Text 加常數(shù)和中央子午線x = Text9.Text : y = Text10.Text反算y = y - yy1e12 = (a1 2 - b1 2) / a1 2 第一偏心率平方e22 = (a1 2 - b1 2) / b1 2 第二偏心率平方e1 = Sqr(a1 2 - b1 2) / a1 2) 第一偏心率e2 = Sqr(a1 2 - b1 2) / b1 2) 第二偏心率c = a1 / Sqr(1 - e12) 極點子午線曲率半徑bb0 = 1 - (3 / 4) * e22 + (45 / 64) * e22 2 - (175 / 256) * e22 3 + (1

44、1025 / 16384) * e22 4bb2 = bb0 - 1bb4 = (15 / 32) * e22 2 - (175 / 384) * e22 3 + (3675 / 8192) * e22 4bb6 = -(35 / 96) * e22 3 + (735 / 2048) * e22 4bb8 = (315 / 1024) * e22 4bf1 = x / (c * bb0) 迭代 單位弧度Do bf = bf1 bf1 = (x / c) - (bb2 * Cos(bf) + bb4 * (Cos(bf) 3) + bb6 * (Cos(bf) 5) + bb8 * (Cos(b

45、f) 7) * Sin(bf) / bb0Loop Until Abs(bf - bf1) 0.000000001mf = (a1 * (1 - e12) / (Sqr(1 - e12 * (Sin(bf) 2) 3)tf = Sin(bf) / Cos(bf)nnf = e2 * Cos(bf)nf = a1 / Sqr(1 - e12 * (Sin(bf) 2)B L單位弧度B = bf - (tf * (y 2) / (2 * mf * nf) + (tf * (y 4) * (5 + 3 * (tf 2) + (nnf 2) - 9 * (nnf 2) * (tf 2) / (24 *

46、 mf * (nf 3) - (tf * (y 6) * (61 + 90 * (tf 2) + 45 * (tf 4) / (720 * mf * (nf 5)L = y / (nf * Cos(bf) - (1 + 2 * tf 2 + nnf 2) * (y 3) / (6 * nf 3 * Cos(bf) + (5 + 28 * (tf 2) + 24 * (tf 4) + 6 * (nnf 2) + 8 * (nnf 2) * (tf 2) * (y 5) / (120 * (nf 5) * Cos(bf)L = L + DuToHu(l1)正算l2 = Text7.Text: yy

47、2 = Text8.Text 加常數(shù)和中央子午線ee12 = (a1 2 - b1 2) / a1 2 第一偏心率平方ee22 = (a1 2 - b1 2) / b1 2 第二偏心率平方ee1 = Sqr(a1 2 - b1 2) / a1 2) 第一偏心率ee2 = Sqr(a1 2 - b1 2) / b1 2) 第二偏心率nn = a1 / Sqr(1 - ee12 * (Sin(B) 2) 卯酉圈曲率半徑tt = Sin(B) / Cos(B)nnn = ee2 * Cos(B)cc = a1 / Sqr(1 - ee12) 極點子午線曲率半徑bbb0 = 1 - (3 / 4) *

48、 ee22 + (45 / 64) * ee22 2 - (175 / 256) * ee22 3 + (11025 / 16384) * ee22 4bbb2 = bbb0 - 1bbb4 = (15 / 32) * ee22 2 - (175 / 384) * ee22 3 + (3675 / 8192) * ee22 4bbb6 = -(35 / 96) * ee22 3 + (735 / 2048) * ee22 4bbb8 = (315 / 1024) * ee22 4 L = HuToDu(L) ll = (Du(L) - l2) / 57.295779513 XX = cc *

49、(bbb0 * B + (bbb2 * Cos(B) + bbb4 * (Cos(B) 3) + bbb6 * (Cos(B) 5) + bbb8 * (Cos(B) 7) * Sin(B) 子午線弧長 x = XX + (nn / 2) * Sin(B) * Cos(B) * (ll 2) + (nn / 24) * Sin(B) * (Cos(B) 3) * (ll 4) * (5 - tt 2 + 9 * (nnn 2) + 4 * (nnn 4) + (nnn / 720) * Sin(B) * (Cos(B) 5) * (ll 6) * (61 - 58 * (tt 2) + (tt

50、 4)y = nn * Cos(B) * ll + (nn / 6) * (Cos(B) 3) * (ll 3) * (1 - tt 2 + nnn 2) + (nn / 120) * (Cos(B) 5) * (ll 5) * (5 - 18 * (tt 2) + tt 4 + 14 * (nnn 2) - 58 * (tt 2) * (nnn 2) y = y + yy23.4 二維坐標轉(zhuǎn)換與四參數(shù)計算由于矩陣計算代碼繁瑣，由后臺調(diào)用Excel函數(shù)完成計算。其主要代碼如下：Dim x#, y#, x1#, y1#, dx#, dy#, jd#, k#正算： x1 = (Cos(DuToHu

51、(jd) * x - Sin(DuToHu(jd) * y) * k + dx y1 = (Sin(DuToHu(jd) * x + Cos(DuToHu(jd) * y) * k + dy反算： x = (x1 - dx) * Cos(DuToHu(jd) + (y1 - dy) * Sin(DuToHu(jd) / k y = (x1 - dx) * (-Sin(DuToHu(jd) + (y1 - dy) * Cos(DuToHu(jd) / k四參數(shù)計算：Dim xls As Excel.ApplicationDim book As Excel.WorkbookDim sheet As

52、 Excel.WorksheetSet xls = CreateObject(Excel.Application) 創(chuàng)建EXCEL對象Set book = xls.Workbooks.Open(Text22.Text)打開已經(jīng)存在的工件簿文件Set sheet = book.Worksheets(1)sheet.Activatea = sheet.UsedRange.Rows.Countm = 1For i = 2 To a x = xls.Cells(i, 2).Value: y = xls.Cells(i, 3).Value XX = xls.Cells(i, 4).Value: yy = xls.Cells(i, 5).Value xls.Cells(m, 10).Value = XX - x: xls.Cells(m + 1, 10).Value = yy - y 參L xls.Cells(m, 6).Value = 1: xls.Cells(m, 7).Value = 0: xls.Cells(m, 8).Value = y: xls.Cells(m