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

2、空間直角坐標的相互轉換、高斯投影正反算、二維坐標轉換與四參數計算、三維坐標轉換與七參數轉換、同參考基準下的坐標換帶計算，以及坐標數據的批量處理。關鍵字：坐標系統，轉換模型，坐標轉換，程序設計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

3、. So It 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 fo

4、llows: 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 c

5、oordinate. 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 use

6、d in the 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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

21、標系與空間直角坐標系轉換模型將同一坐標系下的大地坐標（B，L，H）轉換成空間直角坐標（X，Y，Z）的轉換公式為： 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）轉換為大地坐標（B，L，H）的公式為: L=arctan(YX) B=arctan(Z+Ne2sinBX2+Y2) H=ZsinB-N(1-e2)用公式進行空間直角坐標轉換大地坐標時，需要采用迭代計算大地緯度B。具體計算時，可先根據下式求出大地緯度B的初值：

22、B=arctan(ZX2+Y2)因為 N=c1+e'2cos2B , 1cos2B=1+tan2B 帶入B=arctan(Z+Ne2sinBX2+Y2) 得:tanB=ZX2+Y2+ce2tanBX2+Y2*1+e'2+tan2B令:t=ZX2+Y2，p=ce2X2+Y2，k=1+e'2則式子可寫成:tanB'=t+p*tanBk+tan2B然后利用該初值B代入公式右端tanB中，將等式左邊的結果再次代入右端tanB，直到最后兩次B值之差小與允許誤差為止。當得到大地緯度B后，代入公式即可求出大地高H。2.5.2高斯正反算轉換模型得到了點的大地坐標(L，B)，就可

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-34e'2+4564e'4-175256e'6+1102516384e'8 2=0-1 4=4564e'4-

24、175256e'6+1102516384e'8 6=-175256e'6+1102516384e'8 8=3151024e'8 N=a1-e2sin2B t=tanB =e'cosB l=(L-L0)'''' 其計算結果的精度可達0. 001m。只要得到了高斯平面坐標(X, Y)后，便可通過高斯反算公式將其轉換成大地坐標(B, L)，高斯投影反算公式為: B=Bf-tf2MfNfy2+tf24MfNf35+3tf2+f2-9f2tf2y4-tf720MfNf5(61+90tf2+45tf4)y6 L=1NfcosB

25、fy-16Nf3cosBf1+2tf2+f2y3+1120Nf5cosBf(5+28tf2+24tf4)y5 Bf=Xc-2cosB+4cos3B+6cos5B+8cos7BsinB0 tf=tanBf Mf=a(1-e2)2(1-e2sin2Bf)3 Nf=a1-e2sin2Bf f=tanBf 它們的計算精度,即平面坐標可達0.0001''。2.5.3坐標轉換與參數計算轉換模型二維轉換模型：2個平移參數x y （原點不重合產生的）1個旋轉參數（坐標軸不平行產生的）1個尺度參數m（兩坐標系間的尺度不一致產生的）x2y2=1+mcossin-sincosx1y1+xy(x1,y

26、1)為某點在A空間直角坐標系中的坐標。(x2,y2)為某點在B空間直角坐標系中的坐標。(x,y)為某點從A空間直角坐標系轉換到B空間直角坐標系的兩個平移參數。為從A空間直角坐標系轉換到B空間直角坐標系中標系的兩個平移參數。為從A空間直角坐標系轉換到B空間直角坐標系中一個旋轉參數。m為從A空間直角坐標系轉換到B空間直角坐標系中的一個尺度參數。平面四參數求解步驟如下:利用公共點計算坐標參數，但至少有兩個公共點，當有i個公共點時，可利用最小二乘原理求解參數。將B直角坐標系中的坐標視為觀測值，設A直角坐標系下的坐標視為無誤差，列誤差方程為：vx1vy1vxivyiB=10100101y1-x1yi-x

27、ix1y1xiyixym-X1Y1XiYiB-X1Y1XiYiA寫成矩陣形式既：V=BX-L由于各點的坐標可視為同精度獨立觀測值，因此P=I。把各點坐標已知值帶入上述誤差方程，然后按下列公式求解出四參數：X=(BTB)-1(BTL)三維轉換模型：3個平移參數x y z（原點不重合產生的）3個旋轉參數x y z（坐標軸不平行產生的）1個尺度參數m（兩坐標系間的尺度不一致產生的） x2y2z2=1+mR(x)R(y)R(z)x1y1z1+xyzRx=1 000cosxsinx0-sinxcosxRy=cosy0-siny010sinycosy0Rz=coszsinz0-sinzcosz0001一般

28、x，x，x為微小轉角，模型可簡化為：x2y2z2=1+m1z-y-z1xy-x1x1y1z1+xyz(x1,y1,z1)為某點在A空間直角坐標系中的坐標。(x2,y2z2)為某點在B空間直角坐標系中的坐標。(x,y,z)為某點從A空間直角坐標系轉換到B空間直角坐標系的三個平移參數。(x，y，z)為從A空間直角坐標系轉換到B空間直角坐標系中三個旋轉參數。m為從A空間直角坐標系轉換到B空間直角坐標系中的一個尺度參數。七參數求解步驟如下:利用公共點計算坐標參數，但至少有三個公共點，當有i個公共點時，可利用最小二乘原理求解參數。將B空間直角坐標系中的坐標視為觀測值，設A空間直角坐標系下的坐標視為無誤差

29、，列誤差方程為：vx1vy1vz1vxivyivziB=1001000100100010010Z1-Y10Zi-Yi-Z10X1-Zi0Xi-Y1-X10-Yi-Xi0X1Y1Z1XiYiZixyzxyzm-X1Y1Z1XiYiZiB-X1Y1Z1XiYiZiA寫成矩陣形式既：V=BX-L由于各點的坐標可視為同精度獨立觀測值，因此P=I。把各點坐標已知值帶入上述誤差方程，然后按下列公式求解出七參數：X=(BTB)-1(BTL)2.5.4三、六度帶帶號與中央子午線計算模型求帶號、中央子午線公式為：3°帶:帶號N=IntL-1.53+1 中央子午線L0=N*3 L為當地經度6°

30、帶:帶號N=IntL6+1 中央子午線L0=N*6-3 L為當地經度第3章 坐標轉換程序設計 計劃程序有哪些功能 程序界面設計 搜集坐標轉換公式 運用VB編寫程序 檢查能否正常運行 以及數據正確性 否 是完成程序圖3-1 編寫程序流程圖圖3-2程序主界面3.1大地坐標系與空間直角坐標系的轉換主要代碼如下：橢球參數設置If Combo1.Text = "WGS-84" Then a1 = 6378137: b1 = 6356752.3142 'a1為長半軸，b1為短半軸 Text7.Text = a1: Text8.Text = b1ElseIf Combo1.Tex

31、t = "CGCS2000" Then a1 = 6378137: b1 = 6356752.31414 Text7.Text = a1: Text8.Text = b1ElseIf Combo1.Text = "北京54" Then a1 = 6378245: b1 = 6356863.0188 Text7.Text = a1: Text8.Text = b1ElseIf Combo1.Text = "西安80" Then a1 = 6378140: b1 = 6356755.2882 Text7.Text = a1: Text8.

32、Text = b1ElseIf Combo1.Text = "自定義" Then Text7.Enabled = True: Text8.Enabled = True a1 = Text7.Text: b1 = Text8.TextEnd If大地坐標轉空間直角坐標主要代碼 B = Text1.Text: L = Text2.Text: H = Text3.Text e12 = (a1 2 - b1 2) / a1 2 N = a1 / Sqr(1 - e12 * (Sin(DuToHu(B) 2) x = (N + H) * Cos(DuToHu(B) * Cos(DuT

33、oHu(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空間直角坐標轉大地坐標主要代碼 x = Text4.Text: y = Text5.Text: z = Text6.Text L = HuToDu(Atn(y / x) '計算經度 If y / x < 0 Then L = HuToDu(Atn(y / x) + 3.14159265358979) End I

34、f 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) '用TanB作為未知量 Do Until Abs(TB1 - TB) < 0.000000000001 '迭代法 TB = TB1 TB1 = t

35、 + (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高斯平面坐標與大地坐標的轉換高斯平面坐標(X, Y)與大地坐標(B ,L)的相互關系式分為兩類：第一類稱高斯投影正算公式，由(B,L)求(X,Y)；第二類稱高斯投影反算公式，由(X,Y)求(B,L) 13。由于電子計算機和

36、各種可編程序電子計算器在測量上廣泛應用，因此在這里給出適合電算的高斯正反算公式。主要代碼如下：正算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 = Text2.Textl0 = Text5.Text '輸入中央子午線經度yj = Text8.Te

37、xt 'Y坐標加常數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 / 64) * e22 2 - (1

38、75 / 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 = Sin(DuToHu

39、(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(DuToHu(B)

40、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) - 58

41、 * (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)

42、/ Cos(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) *

43、 (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)3.3 高斯投影換帶計算高斯投影雖然保證了角度沒有變形，但是長度變形比較嚴重，為了限制高斯投影的長度變形，必須以中央子午線進行分帶，把投影范圍限制在中央子午線東、西側一定的范圍內；但是這樣又使得統一的坐標系分割成各帶的獨立坐標系14。這樣就會產生新的問題，需要鄰帶換算來解決。主要代碼如下：l1 = Text3.Text : yy1 =

44、 Text4.Text '加常數和中央子午線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)

45、 * 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 4bf1 = x / (c * bb0) '迭代 單位弧度Do bf = bf1 bf1 = (x / c) - (bb2 * Cos(bf) + bb4 * (Cos(b

46、f) 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(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

47、 * (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) / (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

48、) * Cos(bf)L = L + DuToHu(l1)正算l2 = Text7.Text: yy2 = Text8.Text '加常數和中央子午線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)

49、 / Cos(B)nnn = ee2 * Cos(B)cc = a1 / Sqr(1 - ee12) '極點子午線曲率半徑bbb0 = 1 - (3 / 4) * 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 4b

50、bb8 = (315 / 1024) * ee22 4 L = HuToDu(L) ll = (Du(L) - l2) / 57.295779513 XX = cc * (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 * (n

51、nn 2) + 4 * (nnn 4) + (nnn / 720) * Sin(B) * (Cos(B) 5) * (ll 6) * (61 - 58 * (tt 2) + (tt 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 二維坐標轉換與四參數計算由于矩陣計

52、算代碼繁瑣，由后臺調用Excel函數完成計算。其主要代碼如下：Dim x#, y#, x1#, y1#, dx#, dy#, jd#, k#正算： x1 = (Cos(DuToHu(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) *

53、Cos(DuToHu(jd) / k四參數計算：Dim xls As Excel.ApplicationDim book As Excel.WorkbookDim sheet As Excel.WorksheetSet xls = CreateObject("Excel.Application") '創建EXCEL對象Set book = xls.Workbooks.Open(Text22.Text)'打開已經存在的工件簿文件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(