VB實(shí)現(xiàn)高斯正反算

1、VB實(shí)現(xiàn)高斯正反算專業(yè)班級(jí)：測(cè)繪0729班學(xué)生姓名：張智寧學(xué)生學(xué)號(hào)：0700001955指導(dǎo)教師：歐陽老師完成時(shí)間：2010年6月5日1前言在日益發(fā)展的社會(huì)中，能很好很快的完成各種任務(wù)才能取得更多的發(fā)展機(jī)會(huì)。用VB實(shí)現(xiàn)高斯正反算，就能很快實(shí)現(xiàn)高斯坐標(biāo)系和大地坐標(biāo)系的轉(zhuǎn)換，為我們的計(jì)算和測(cè)量工作提供高效的服務(wù)。此次設(shè)計(jì)要求我們利用大學(xué)所學(xué)的VB編程和控制測(cè)量中的高斯平面坐標(biāo)系為基礎(chǔ)，用VB實(shí)現(xiàn)高斯正反算，做到能夠完整寫出VB實(shí)現(xiàn)高斯正反算的代碼，而且還能生成運(yùn)行文件進(jìn)行高斯正反算。同時(shí)，也讓我們體會(huì)到編程在測(cè)量中應(yīng)用之廣，在以后的工作中應(yīng)該博學(xué)多思 。2控制測(cè)量的基本理論1.大地水準(zhǔn)面大地水準(zhǔn)面

2、高是大地水準(zhǔn)面至地球橢球面的垂直距離。與靜止海平面重合的重力等位面稱為大地水準(zhǔn)面。大地水準(zhǔn)面是一個(gè)與地球內(nèi)部密度分布有關(guān)的不規(guī)則曲面。實(shí)際觀測(cè)結(jié)果表明，大地水準(zhǔn)面非常接近一個(gè)旋轉(zhuǎn)橢球面。選擇適當(dāng)?shù)男D(zhuǎn)橢球，使得大地水準(zhǔn)面相對(duì)此橢球面的起伏的平方在橢球面上的積分最小，這個(gè)旋轉(zhuǎn)橢球稱為參考橢球。大地水準(zhǔn)面與參考橢球面的垂直距離稱為“大地水準(zhǔn)面高”。在1980大地參考系中，參考橢球的赤道半徑取為6 378.137米，扁率的倒數(shù)取為286 26，此時(shí)大地水準(zhǔn)面的高度不超過110米。地面點(diǎn)的正高與該點(diǎn)的大地水準(zhǔn)面高的和稱為該點(diǎn)的大地高，大地高是地面點(diǎn)的大地坐標(biāo)系中的一個(gè)坐標(biāo)。大地水準(zhǔn)面高在地面上的分布與

3、地球內(nèi)部特別是地殼上地幔密度水平向分布不均勻有關(guān)，因此大地水準(zhǔn)面高也是研究地球內(nèi)部結(jié)構(gòu)和地球內(nèi)部動(dòng)力學(xué)過程的一種。圖1-13圖1-14水準(zhǔn)面靜止的水面稱為水準(zhǔn)面，水準(zhǔn)面是受地球重力影響而形成的，是一個(gè)處處與重力方向垂直的連續(xù)曲面，并且是一個(gè)重力場(chǎng)的等位面。設(shè)想一個(gè)靜止的海水面擴(kuò)展到陸地部分。這樣，地球的表面就形成了一個(gè)較地球自然表面規(guī)則而光滑的曲面，這個(gè)曲面被稱為水準(zhǔn)面。圖1-25高斯投影高斯-克呂格(Gauss-Kruger)投影簡(jiǎn)稱“高斯投影”，又名等角橫切橢圓柱投影”，地球橢球面和平面間正形投影的一種。德國(guó)數(shù)學(xué)家、物理學(xué)家、天文學(xué)家高斯（Carl Friedrichauss，1777一

4、1855）于十九世紀(jì)二十年代擬定，后經(jīng)德國(guó)大地測(cè)量學(xué)家克呂格（Johannes Kruger，18571928）于 1912年對(duì)投影公式加以補(bǔ)充，故名。該投影按照投影帶中央子午線投影為直線且長(zhǎng)度不變和赤道投影為直線的條件，確定函數(shù)的形式，從而得到高斯一克呂格投影公式。投影后，除中央子午線和赤道為直線外， 其他子午線均為對(duì)稱于中央子午線的曲線。設(shè)想用一個(gè)橢圓柱橫切于橢球面上投影帶的中央子午線，按上述投影條件，將中央子午線兩側(cè)一定經(jīng)差范圍內(nèi)的橢球面正形投影于橢圓柱面。將橢圓柱面沿過南北極的母線剪開展平，即為高斯投影平面。取中央子午線與赤道交點(diǎn)的投影為原點(diǎn)，中央子午線的投影為縱坐標(biāo)x軸，赤道的投影為

5、橫坐標(biāo)y軸，構(gòu)成高斯克呂格平面直角坐標(biāo)系。高斯-克呂格投影在長(zhǎng)度和面積上變形很小，中央經(jīng)線無變形，自中央經(jīng)線向投影帶邊緣，變形逐漸增加，變形最大之處在投影帶內(nèi)赤道的兩端。由于其投影精度高，變形小，而且計(jì)算簡(jiǎn)便（各投影帶坐標(biāo)一致，只要算出一個(gè)帶的數(shù)據(jù)，其他各帶都能應(yīng)用），因此在大比例尺地形圖中應(yīng)用，可以滿足軍事上各種需要，能在圖上進(jìn)行精確的量測(cè)計(jì)算。高斯投影：它是一種橫軸等角切圓柱投影。它把地球視為球體,假想一個(gè)平面卷成一個(gè)橫圓柱面并把它套在球體外面,使橫軸圓柱的軸心通過球的中心,球面上一根子午線與橫軸圓柱面相切。這樣,該子午線在圓柱面上的投影為一直線,赤道面與圓柱面的交線是一條與該子午線投影垂

6、直的直線。將橫圓柱面展開成平面,由這兩條正交直線就構(gòu)成高斯-克呂格平面直角坐標(biāo)系。為減少投影變形,高斯-克呂格投影分為3o帶和6o帶投影。 高斯克呂格投影是設(shè)想用一個(gè)橢圓柱橫套在地球橢球的外面,并與設(shè)定的中央經(jīng)線相切。 6高斯克呂格投影分帶規(guī)定:該投影是國(guó)家基本比例尺地形圖的數(shù)學(xué)基礎(chǔ),為控制變形,采用分帶投影的方法,在比例尺 1:2.5萬-1:50萬圖上采用6分帶,對(duì)比例尺為 1:1萬及大于1:1萬的圖采用3分帶。 6分帶法:從格林威治零度經(jīng)線起,每6分為一個(gè)投影帶,全球共分為60個(gè)投影帶,東半球從東經(jīng)0-6為第一帶,中央經(jīng)線為3,依此類推,投影帶號(hào)為1-30。其投影代號(hào)n和中央經(jīng)線經(jīng)度L0的

7、計(jì)算公式為:L0=(6n-3);西半球投影帶從180回算到0,編號(hào)為31-60,投影代號(hào)n和中央經(jīng)線經(jīng)度L0的計(jì)算公式為L(zhǎng)0=360-(6n-3)。 3分帶法:從東經(jīng)130起,每3為一帶,將全球劃分為120個(gè)投影帶,東經(jīng)130-430,.17830-西經(jīng)17830,.130-東經(jīng)130。東半球有60個(gè)投影帶,編號(hào)1-60,各帶中央經(jīng)線計(jì)算公式:L0=3n ,中央經(jīng)線為3、6.180。西半球有60個(gè)投影帶,編號(hào)1-60,各帶中央經(jīng)線計(jì)算公式:L0=360-3n ,中央經(jīng)線為西經(jīng)177、.3、0。我國(guó)規(guī)定將各帶縱坐標(biāo)軸西移500公里,即將所有y值加上500公里,坐標(biāo)值前再加各帶帶號(hào)以18帶為例,原

8、坐標(biāo)值為y=243353.5m,西移后為y=743353.5,加帶號(hào)通用坐標(biāo)為y=18743353.5 。7VB簡(jiǎn)介vb是Visual Basic的簡(jiǎn)稱，是由美國(guó)微軟公司于1991年開發(fā)的一種可視化的、面向?qū)ο蠛筒捎檬录?qū)動(dòng)方式的結(jié)構(gòu)化高級(jí)程序設(shè)計(jì)語言，可用于開發(fā) Windows 環(huán)境下的各類應(yīng)用程序。它簡(jiǎn)單易學(xué)、效率高，且功能強(qiáng)大可以與 Windows 專業(yè)開發(fā)工具SDK相媲美。在Visual Basic環(huán)境下，利用事件驅(qū)動(dòng)的編程機(jī)制、新穎易用的可視化設(shè)計(jì)工具，使用Windows內(nèi)部的廣泛應(yīng)用程序接口（API）函數(shù)，動(dòng)態(tài)鏈接庫（DLL）、對(duì)象的鏈接與嵌入（OLE）、開放式數(shù)據(jù)連接（ODBC）

9、等技術(shù)，可以高效、快速地開發(fā)Windows環(huán)境下功能強(qiáng)大、圖形界面豐富的應(yīng)用軟件系統(tǒng)。8三 高斯正反算公式1，高斯投影正算已知橢球面上某點(diǎn)的大地坐標(biāo)（L,B），該點(diǎn)在高斯投影平面上的直角坐標(biāo)（x,y）， 即（L,B）=(X,Y)的坐標(biāo)轉(zhuǎn)換。 1.1 投影必須滿足的條件（1）中央子午線投影后為直線。（2）中央子午線投影后長(zhǎng)度不變。（3）投影具有正形性質(zhì)，即正投影條件。1.2 正投影過程 在橢球上有對(duì)稱于中央子午線的兩點(diǎn)P1和P2，它們的大地坐標(biāo)分別為（L,B）及（l,b），其中l(wèi)為橢球面上p點(diǎn)的經(jīng)度與中央子午線的經(jīng)度L0之差：l=L-L0，P點(diǎn)在中央子午線之東，l為正，在西則為負(fù)，則投影后的平面

10、坐標(biāo)一定為P1(X,Y)和P2(X,Y)。9高斯投影正算計(jì)算公式 2.高斯投影坐標(biāo)反算公式2.1高斯投影反算 已知某點(diǎn)的高斯投影面上直角坐標(biāo)（x,y），求該點(diǎn)在橢球面上的大地坐標(biāo)（L,B），即（x,y）=(L,B)的坐標(biāo)變換。2.2 投影變換必須滿足的條件a. x坐標(biāo)軸投影成中央子午線，是投影的對(duì)稱軸。b. x軸上的長(zhǎng)度投影保持不變。c.投影具有正形性質(zhì)，即正形投影條件。2.3 根據(jù)x計(jì)算縱坐標(biāo)在橢球面上的投影的底點(diǎn)緯度，接著按計(jì)算（-B）及經(jīng)差，最后得到B= L=10高斯投影反算公式當(dāng)要求轉(zhuǎn)換精度至?xí)r ，可以簡(jiǎn)化為下式：11VB設(shè)計(jì)界面以上就是VB實(shí)現(xiàn)高斯正反算得操作界面。12程序使用方法

11、這個(gè)程序使用很簡(jiǎn)單，首先確認(rèn)中央子午線的，點(diǎn)擊中央子午線下方的小黑三角形選擇大地坐標(biāo)或或高斯投影平面坐標(biāo)所在的中央子午線的經(jīng)度。在計(jì)算之前還要選擇所依據(jù)的橢球面，橢球不一樣，計(jì)算的參數(shù)不一樣。若選擇“北京54”，則點(diǎn)擊“北京54”前面的小白圓圈，顯黑點(diǎn)就說明已經(jīng)選中，例如上面界面所顯示的那樣。在選定中央子午線和橢球后，要大地坐標(biāo)轉(zhuǎn)高斯投影平面坐標(biāo)，則在大地坐標(biāo)下輸入B,L，按上方的轉(zhuǎn)換按鈕，則可以得到高斯投影平面坐標(biāo)X,Y。反之，輸入高斯投影平面坐標(biāo)X,Y，按下方的轉(zhuǎn)換按鈕，則可以得到大地坐標(biāo)B,L。13程序代碼 這個(gè)程序使用很簡(jiǎn)單，首先確認(rèn)中央子午線的，點(diǎn)擊中央子午線下方的小黑三角形選擇大地

12、坐標(biāo)或或高斯投影平面坐標(biāo)所在的中央子午線的經(jīng)度。在計(jì)算之前還要選擇所依據(jù)的橢球面，橢球不一樣，計(jì)算的參數(shù)不一樣。若選擇“北京54”，則點(diǎn)擊“北京54”前面的小白圓圈，顯黑點(diǎn)就說明已經(jīng)選中，例如上面界面所顯示的那樣。在選定中央子午線和橢球后，要大地坐標(biāo)轉(zhuǎn)高斯投影平面坐標(biāo)，則在大地坐標(biāo)下輸入B,L，按上方的轉(zhuǎn)換按鈕，則可以得到高斯投影平面坐標(biāo)X,Y。反之，輸入高斯投影平面坐標(biāo)X,Y，按下方的轉(zhuǎn)換按鈕，則可以得到大地坐標(biāo)B,L。六程序代碼 Public Const PI = 3.149Public a As Double, e2 As Double, e22 As DoublePublic B As

13、 Double, L As Double, L0 As DoublePublic Bf As DoublePublic x0 As Double, y0 As DoublePrivate ll As DoublePublic Number As IntegerPublic AA As Double, BB As Double, CC As Double, DD As Double, EE As Double, FF As Double, GG As DoublePublic Sub Gauss(mL0 As Double, mB As Double, mL As Double) Dim T A

14、s Double, Yita As Double mL0 = mGetrad(mL0): mB = mGetrad(mB): mL = mGetrad(mL) Yita = e22 * (Cos(mB) 2): T = Tan(mB) ll = mL - mL0 x0 = mGetX(0, mB) + mGetN(mB) * T * (Cos(mB) 2 * ll 2 * (0.5 + (5 - T 2 + 9 * Yita + 4 * Yita 2) * (Cos(mB) 2 * ll 2 / 24 + (61 - 58 * T 2 + T 4) * (Cos(mB) 4 * ll 4 /

15、720) y0 = mGetN(mB) * Cos(mB) * ll * (1 + (1 - T 2 + Yita) * (Cos(mB) 2 * ll 2 / 6 + (5 - 18 * T 2 + T 4 + 14 * Yita - 58 * T 2 * Yita) * (Cos(mB) 4 * ll 4 / 120)End SubPublic Sub mGauss(mL0 As Double, mB As Double, mY As Double)14Dim T As Double, Yita As Double Dim N As Double, M As Double Yita = e

16、22 * (Cos(mB) 2): T = Tan(mB) L0 = mGetrad(mL0) N = mGetN(mB) M = mGetM(mB) ll = mY / Cos(mB) / N * (1 - (1 + 2 * T 2 + Yita) * (mY / N) 2 / 6 + (5 + 28 * T 2 + 24 * T 4 + 6 * Yita + 8 * Yita * T 2) * (mY / N) 4 / 120) B = mB - T * mY * mY / N / M / 2 * (1 - (5 + 3 * T 2 + Yita - 9 * Yita * T 2) * (

17、mY / N) 2 / 12 + (61 + 90 * T 2 + 45 * T 4) * (mY / N) 4 / 360) L = L0 + llEnd SubPublic Sub mGetBf(mB As Double) If Abs(x0 - mGetX(0, mB) 0.001 Then Bf = mB Exit Sub Else mB = mB + (x0 - mGetX(0, mB) / mGetSS(mB) Number = Number + 1 Call mGetBf(mB) End IfEnd SubPrivate Function mGetX(B1 As Double,

18、B2 As Double)15 mGetX = mGetM(B1 + B2) / 2) * (B2 - B1) * (1 + e22 * Cos(B1 + B2) * (B1 - B2) 2 / 8)End FunctionPrivate Function mGetSS(mB As Double) Dim R As Double, S As Double, T As Double R = a * (1 - e2) S = AA - BB * Cos(2 * mB) - CC * (Sin(mB) 2 * (Cos(2 * mB) + 2 * (Cos(mB) 2) T = -DD * (Sin

19、(mB) 4 * (Cos(2 * mB) + 4 * (Cos(mB) 2) - EE * (Sin(mB) 6 * (Cos(2 * mB) + 6 * (Cos(mB) 2) - FF * (Sin(mB) 8 * (Cos(2 * mB) + 8 * (Cos(mB) 2) mGetSS = R * (S + T) mGetSS = a * (1 - e2) * (AA - BB * Cos(2 * mB) - CC * (Sin(mB) 2 * (Cos(2 * mB) + 2 * (Cos(mB) 2) - DD * (Sin(mB) 4 * (Cos(2 * mB) + 4 *

20、(Cos(mB) 2) - EE * (Sin(mB) 6 * (Cos(2 * mB) + 6 * (Cos(mB) 2) - FF * (Sin(mB) 8 * (Cos(2 * mB) + 8 * (Cos(mB) 2)End FunctionPrivate Function mGetrad(S As Double) S = mTransfer2(S) * PI / 180 mGetrad = SEnd FunctionPublic Function mTransfer1(S As Double) 將以度為單位的角度轉(zhuǎn)換 16Dim du As Integer, fe As Intege

21、r, mi As Single du = Int(S) fe = Int(S - du) * 60) mi = (S - du - fe / 60) * 3600 mTransfer1 = du + fe / 100 + mi / 10000End FunctionPublic Function mTransfer2(R As Double) 將小數(shù)度分秒格式的角度轉(zhuǎn)換為以度為單位的格式 Dim du As Integer, fe As Integer, mi As Single If R 0 Then R = Abs(R) du = Int(R) fe = Int(R - du) * 100

22、) mi = (R * 100 - Int(R * 100) * 100 mTransfer2 = -(du + fe / 60 + mi / 3600) Else du = Int(R) fe = Int(R - du) * 100) mi = (R * 100 - Int(R * 100) * 100 mTransfer2 = du + fe / 60 + mi / 3600 End If17End FunctionPublic Function mTransfer3(R As Double) 檢驗(yàn)出單位秒是60時(shí)，進(jìn)行進(jìn)位 Dim du As Integer, fe As Integer, mi As Single If R 0 Then R = Abs(R) du = Int(R) fe = Int(R - du) * 100) mi = (R * 100 - Int(R * 100) * 100 If mi = 60 Then mi = 0: fe = fe + 1 If fe = 60 Then fe = 0