坐標(biāo)轉(zhuǎn)換參數(shù)求取及坐標(biāo)轉(zhuǎn)換程序設(shè)計(jì)

畢業(yè)設(shè)計(jì)設(shè)計(jì)題目 坐標(biāo)轉(zhuǎn)換參數(shù)求取及坐標(biāo)轉(zhuǎn)換程序設(shè)計(jì) 學(xué)生姓名 張威 指導(dǎo)教師 杜繼亮 專 業(yè) 測(cè)繪工程 班 級(jí) 測(cè)繪12-2班 填寫日期 2016/4/6 礦業(yè)工程學(xué)院摘 要坐標(biāo)系統(tǒng)是測(cè)量工作中定位的基礎(chǔ)，坐標(biāo)系統(tǒng)有多種形式和基準(zhǔn)，由于各測(cè)量工作目的不同，所選用的坐標(biāo)基準(zhǔn)也會(huì)不同，根據(jù)不同的工作要求需要將不同的坐標(biāo)系下的坐標(biāo)進(jìn)行相互轉(zhuǎn)換。在這些坐標(biāo)轉(zhuǎn)換的過程中會(huì)用到很多坐標(biāo)轉(zhuǎn)換模型，但是坐標(biāo)系轉(zhuǎn)換模型過于復(fù)雜手算非常困難。本設(shè)計(jì)為了方便施工時(shí)遇到的坐標(biāo)轉(zhuǎn)換問題，設(shè)計(jì)利用Visual Basic 6.0編程語言編寫程序，用來實(shí)現(xiàn)坐標(biāo)系統(tǒng)之間的轉(zhuǎn)換以及轉(zhuǎn)換參數(shù)的求解，例如:大地坐標(biāo)與空間直角坐標(biāo)的相互轉(zhuǎn)換、高斯投影正反算、二維坐標(biāo)轉(zhuǎn)換與四參數(shù)計(jì)算、三維坐標(biāo)轉(zhuǎn)換與七參數(shù)轉(zhuǎn)換、同參考基準(zhǔn)下的坐標(biāo)換帶計(jì)算，以及坐標(biāo)數(shù)據(jù)的批量處理。關(guān)鍵字：坐標(biāo)系統(tǒng)，轉(zhuǎn)換模型，坐標(biāo)轉(zhuǎn)換，程序設(shè)計(jì)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 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 follows: 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 coordinate. 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 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 the transformation 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程序設(shè)計(jì)思路方法2第2章 坐標(biāo)理論和知識(shí)32.1地球橢球32.2基準(zhǔn)32.3測(cè)量常用坐標(biāo)系42.3.1大地坐標(biāo)系42.3.2空間直角坐標(biāo)系42.3.3平面坐標(biāo)系52.3.4地方獨(dú)立坐標(biāo)系52.4我國常用坐標(biāo)系62.4.1 1954年北京坐標(biāo)系62.4.2 1980國家大地坐標(biāo)系62.4.3 WGS-84世界大地坐標(biāo)系62.4.4 2000國家大地坐標(biāo)系7第3章 坐標(biāo)轉(zhuǎn)換程序設(shè)計(jì)83.1大地坐標(biāo)系與空間直角坐標(biāo)系的轉(zhuǎn)換93.1.1大地坐標(biāo)轉(zhuǎn)換成空間直角坐標(biāo)93.1.2空間直角坐標(biāo)轉(zhuǎn)成換大地坐標(biāo)103.2高斯平面坐標(biāo)與大地坐標(biāo)的轉(zhuǎn)換133.2.1 高斯正算133.2.2 高斯反算143.3 高斯投影換帶計(jì)算173.4 二維坐標(biāo)轉(zhuǎn)換與四參數(shù)計(jì)算203.5 三維坐標(biāo)轉(zhuǎn)換與七參數(shù)計(jì)算233.6 三、六度帶帶號(hào)與中央子午線計(jì)算27第4章 設(shè)計(jì)總結(jié)31致謝32參考文獻(xiàn)33IV第1章 緒論1.1研究背景和意義隨著大地測(cè)量學(xué)，攝影測(cè)量學(xué)的發(fā)展和電子計(jì)算機(jī)的普及，對(duì)各種坐標(biāo)系的研究變得越來越重要了1。隨著現(xiàn)在社會(huì)的快速發(fā)展，各種各樣的大型工程的建設(shè)，凡是工程施工就必定需要坐標(biāo)來定位才可以，而建造的地方又不同，又沒有一個(gè)滿足全部地方的坐標(biāo)系統(tǒng)，所以產(chǎn)生多種不同的坐標(biāo)系統(tǒng)，實(shí)際工作中測(cè)量人員必定會(huì)按照實(shí)際情況選擇最為合適的坐標(biāo)系統(tǒng)。而坐標(biāo)系統(tǒng)之間的轉(zhuǎn)換比較復(fù)雜，手算工作量巨大，因而各種坐標(biāo)轉(zhuǎn)換模型相繼出現(xiàn)，利用計(jì)算機(jī)強(qiáng)大的數(shù)據(jù)計(jì)算能力可以輕松應(yīng)對(duì)這些問題，提高工作效率。坐標(biāo)轉(zhuǎn)換的意義重大，不僅在我們熟知的工程領(lǐng)域中，在國防建設(shè)、航空航天科技、城市匯劃等眾多領(lǐng)域中都發(fā)揮著重要的作用，可以說對(duì)社會(huì)進(jìn)步有著必不可少的作用。1.2國內(nèi)外研究現(xiàn)狀自60年代以來，各國大地測(cè)量學(xué)者，經(jīng)過大量研究，提出了多種坐標(biāo)轉(zhuǎn)換模型及多種解算方法，北美1927基準(zhǔn)面（基于克拉克1966橢球體與北美1983基準(zhǔn)面（基于GRS1980橢球體）之間坐標(biāo)轉(zhuǎn)換是根據(jù)研究區(qū)內(nèi)一系列己知點(diǎn)的大地坐標(biāo)或網(wǎng)格坐標(biāo)改正量進(jìn)行插值進(jìn)行的坐標(biāo)系轉(zhuǎn)換；英國采用北向與東向的雙線性網(wǎng)格插值進(jìn)行坐標(biāo)轉(zhuǎn)換；挪威在海岸帶調(diào)查中，采用經(jīng)緯度多項(xiàng)式用于坐標(biāo)系轉(zhuǎn)換這種方法進(jìn)行新（ED87一歐洲1987基準(zhǔn)面）、舊（ED50一歐洲1950基準(zhǔn)面）坐標(biāo)系之間的轉(zhuǎn)換：歐洲石油勘探組織（EPSG）對(duì)新、舊坐標(biāo)系采用“雙線性插值” 進(jìn)行坐標(biāo)轉(zhuǎn)換2。國內(nèi)空間三維直角坐標(biāo)轉(zhuǎn)換中，一般采用布爾莎七參數(shù)模型。一般有7個(gè)轉(zhuǎn)換參數(shù)，即3個(gè)平移參數(shù)，3個(gè)旋轉(zhuǎn)參數(shù)和1個(gè)尺度參數(shù)。需要三個(gè)及已經(jīng)公共點(diǎn)時(shí)，才能利用平差的方法求出七參數(shù)。1.3研究的主要內(nèi)容本坐標(biāo)轉(zhuǎn)換程序可實(shí)現(xiàn)功能有:1、大地坐標(biāo)與空間直角坐標(biāo)的相互轉(zhuǎn)換， 2、高斯投影正反算，3、二維坐標(biāo)轉(zhuǎn)換與四參數(shù)計(jì)算，4、三維坐標(biāo)轉(zhuǎn)換與七參數(shù)轉(zhuǎn)換，5、同參考基準(zhǔn)下的坐標(biāo)換帶計(jì)算，以及坐標(biāo)數(shù)據(jù)的批量處理。1.4程序設(shè)計(jì)思路方法本程序名為萬能坐標(biāo)轉(zhuǎn)換器。設(shè)計(jì)前期收集相關(guān)資料，參考一些成熟的坐標(biāo)轉(zhuǎn)換軟件，確定程序應(yīng)有的功能以及界面設(shè)計(jì)。運(yùn)用VB編寫程序時(shí)，查閱相關(guān)書籍獲取理論知識(shí)以及轉(zhuǎn)換模型。完成程序后將已知正確數(shù)據(jù)帶入其中驗(yàn)證程序結(jié)果是否正確。若出現(xiàn)錯(cuò)誤則檢查每步代碼，直到程序完美運(yùn)行為止。第2章 基礎(chǔ)知識(shí)準(zhǔn)備2.1地球橢球地球橢球體又稱地球橢圓體或地球扁球體，代表地球大小和形狀的數(shù)學(xué)曲面，以長半徑和扁率表示，因它十分迫近于橢球體，故通常以參考橢球體表示地球橢球體的形狀和大小。通常所說地球的形狀和大小，實(shí)際上就是以參考橢球體的半長徑、半短徑和扁率來表示。1975年國際大地測(cè)量與地球物理聯(lián)合會(huì)推薦的數(shù)據(jù)為：半長徑6378140米，半短徑6356755米，扁率1298.257。在眾多橢球體中，WGS-84橢球體被認(rèn)為符合上述條件最好的橢球。2.2基準(zhǔn)所謂基準(zhǔn)是指為描述空間位置而定義的點(diǎn)、線和面。而大地測(cè)量基準(zhǔn)是指用以描述地球形狀的地球橢球參數(shù)，包含描述地球橢球幾何特征的長短半軸和物理特征的有關(guān)參數(shù)、地球在空間的定位及定向以及描述這些位置所采用的單位長度的定義。不同的坐標(biāo)系統(tǒng)會(huì)使用的基準(zhǔn)也不同，根據(jù)參考橢球所選原點(diǎn)位置不同，可以分為地心坐標(biāo)系和參心坐標(biāo)系。地心坐標(biāo)系是以地球的質(zhì)心為原點(diǎn)，有地心大地坐標(biāo)系和地心空間直角坐標(biāo)系兩種表述方法。地心空間直角坐標(biāo)系的定義為：以地球質(zhì)心為原點(diǎn)，X軸指向格林尼治子午面與地球赤道的交點(diǎn)，Z軸指向北極，Y軸過原點(diǎn)垂直于平面 XOZ，構(gòu)成右手空間直角坐標(biāo)系。地心大地坐標(biāo)系定義為：以地球的質(zhì)心作為原點(diǎn)，以地球自轉(zhuǎn)軸作為橢球的短軸，大地緯度B是過地面點(diǎn)的橢球法線與橢球赤道面之間的夾角，大地經(jīng)度L為過地面點(diǎn)的橢球子午面與格林尼治子午面之間的夾角，大地高度H為地面點(diǎn)沿橢球法線到橢球面的最短距離。WGS-84坐標(biāo)系，CGCS2000坐標(biāo)系，GLONASS是采用PZ-90坐標(biāo)，都是屬于地心坐標(biāo)系。參心坐標(biāo)系是選取一個(gè)參考橢球面作為基本的參考面，選一參考點(diǎn)作為大地測(cè)量的起算點(diǎn)，從而確定參考橢球在地球面的位置和方向。這時(shí)參考橢球的原點(diǎn)不會(huì)和地球質(zhì)心重合，所以稱為參心。北京54坐標(biāo)系、西安80坐標(biāo)系和新北京54坐標(biāo)系，都是參心坐標(biāo)系。它同樣具有參心大地坐標(biāo)系和參心空間直角坐標(biāo)系兩種表述方法，它們的定義與地心坐標(biāo)系的定義相似。2.3測(cè)量常用坐標(biāo)系2.3.1大地坐標(biāo)系空間一點(diǎn)的大地坐標(biāo)用大地經(jīng)度L、大地緯度B和大地高度H表示，地面上P地點(diǎn)的大地子午面NPS與起始大地子午面所構(gòu)成的二面角L稱P地點(diǎn)的大地經(jīng)度， P地點(diǎn)對(duì)于橢球的法線與赤道面的夾角B稱P地點(diǎn)的大地緯度。如圖2-1所示圖2-1大地坐標(biāo)系P地點(diǎn)沿法線到橢球面的距離H稱大地高，從橢球面起算，向外為正，向內(nèi)為負(fù)3。H = H正常 + （高程異常）H = H正 + N（大地水準(zhǔn)面差距） 2.3.2空間直角坐標(biāo)系空間直角坐標(biāo)系的坐標(biāo)原點(diǎn)與參考橢球的中心重合，Z軸正向指向參考橢球的北極，X軸正向指向起始子午面與赤道的交點(diǎn)，Y軸按右手系與X軸呈90夾角且位于赤道面上，某點(diǎn)在空間中的坐標(biāo)可用該點(diǎn)在此空間坐標(biāo)系的各個(gè)坐標(biāo)軸上的投影來表示4。如圖2-2所示：圖2-2空間直角坐標(biāo)系2.3.3平面坐標(biāo)系平面直角坐標(biāo)系是利用投影，將空間坐標(biāo)通過某種數(shù)學(xué)變換映射到平面上，這種變換稱為投影變換5。在我國一般采用的是高斯一克呂格投影，是目前測(cè)量上廣泛采用的正形投影，特點(diǎn)是沒有角度變形，在不同點(diǎn)上的長度比隨點(diǎn)位而異，但在同一點(diǎn)上各方向的長度比相同，也稱為高斯投影6。2.3.4地方獨(dú)立坐標(biāo)系在我國平面坐標(biāo)主要采用的是高斯投影，在該投影中，除中央子午線外，其它位置上的任何線段，投影后都會(huì)產(chǎn)生一定的長度變形，而且變形隨離開中央子午線的距離增加而增加7。因此一般采用分帶投影的辦法，來限制長度變形，我國規(guī)定了采用3度帶或6度帶進(jìn)行分帶投影。在城市、工礦等工程測(cè)量中，如果直接在國家分帶坐標(biāo)系中建立控制網(wǎng)，會(huì)使地面長度投影的變形較大，當(dāng)長度變形大于25 cm/km時(shí)，就難以滿足工程上的需要8。另一些特殊的測(cè)量，比如大橋施工測(cè)量，水利水壩測(cè)量，滑坡變形監(jiān)測(cè)等，采用國家坐標(biāo)系精度不能滿足工程要求，所以常常會(huì)建立適合本地區(qū)的地方獨(dú)立坐標(biāo)系9。2.4我國常用坐標(biāo)系2.4.1 1954年北京坐標(biāo)系1954年北京坐標(biāo)系，是采取先將我國一等鎖與前蘇聯(lián)遠(yuǎn)東一等鎖相聯(lián)接，然后以連接處呼瑪，吉拉林，東寧基線網(wǎng)擴(kuò)大邊端點(diǎn)的前蘇聯(lián)1942年普爾科沃坐標(biāo)系的坐標(biāo)為起算數(shù)據(jù)，平差我國東北及東部一等鎖10。橢球參數(shù)：長半軸a=6378245m，短半軸b=6356863.0188m，扁率=1/298.3，第一偏心率平方e2=0.00669342161454287，第二偏心率平方e2=0.00673852540614652。2.4.2 1980國家大地坐標(biāo)系1980國家大地坐標(biāo)系采用地球橢球基本參數(shù)為1975年國際大地測(cè)量與地球物理聯(lián)合會(huì)第十六屆大會(huì)推薦的數(shù)據(jù)，大地原點(diǎn)設(shè)在我國中部的陜西省涇陽縣永樂鎮(zhèn)，位于西安市西北方向約60公里。橢球參數(shù)：長半軸a=63781405m，短半軸b=6356755.2882m，扁率=1/298.257，第一偏心率平方e2 =0.00669438499959，第二偏心率平方e2=0.00673950181947。2.4.3 WGS-84世界大地坐標(biāo)系WGS-84坐標(biāo)系是一種國際上采用的地心坐標(biāo)系。坐標(biāo)原點(diǎn)為地球質(zhì)心，其地心空間直角坐標(biāo)系的Z軸指向BIH （國際時(shí)間服務(wù)機(jī)構(gòu)）1984.0定義的協(xié)議地球極（CTP）方向，X軸指向BIH 1984.0的零子午面和CTP赤道的交點(diǎn)，Y軸與Z軸、X軸垂直構(gòu)成右手坐標(biāo)系，稱為1984年世界大地坐標(biāo)系統(tǒng)。橢球參數(shù)：長半軸a=6378137，短半軸b=6356752.3142，扁率=1/298.2572236，第一偏心率平方e2=0.00669437999013 ，第二偏心率平方e2=0.006739496742227。2.4.4 2000國家大地坐標(biāo)系CGCS2000是（中國）2000國家大地坐標(biāo)系的縮寫，該坐標(biāo)系是通過中國GPS連續(xù)運(yùn)行基準(zhǔn)站、空間大地控制網(wǎng)以及天文大地網(wǎng)與空間地網(wǎng)聯(lián)合平差建立的地心大地坐標(biāo)系統(tǒng)11。Z軸指向BIH1984.0定義的協(xié)議極地方向（BIH國際時(shí)間局），X軸指向BIH1984.0定義的零子午面與協(xié)議赤道的交點(diǎn)，Y軸按右手坐標(biāo)系確定。橢球參數(shù)：長半軸a=6378137，短半軸b=6356752.31414，扁率=1/298.2572236，第一偏心率平方e2=0.00669437999013 ，第二偏心率平方e2=0.006739496742227。2.5坐標(biāo)轉(zhuǎn)換模型2.5.1大地坐標(biāo)系與空間直角坐標(biāo)系轉(zhuǎn)換模型將同一坐標(biāo)系下的大地坐標(biāo)（B，L，H）轉(zhuǎn)換成空間直角坐標(biāo)（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。 將同一坐標(biāo)系下的空間直角坐標(biāo)（X，Y，Z）轉(zhuǎn)換為大地坐標(biāo)（B，L，H）的公式為: L=arctan(YX) B=arctan(Z+Ne2sinBX2+Y2) H=ZsinB-N(1-e2)用公式進(jìn)行空間直角坐標(biāo)轉(zhuǎn)換大地坐標(biāo)時(shí)，需要采用迭代計(jì)算大地緯度B。具體計(jì)算時(shí)，可先根據(jù)下式求出大地緯度B的初值：B=arctan(ZX2+Y2)因?yàn)?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值之差小與允許誤差為止。當(dāng)?shù)玫酱蟮鼐暥菳后，代入公式即可求出大地高H。2.5.2高斯正反算轉(zhuǎn)換模型得到了點(diǎn)的大地坐標(biāo)(L，B)，就可以將其轉(zhuǎn)化為某投影帶的高斯平面坐標(biāo)，我們將橢球參數(shù)代入高斯投影正算公式得到更適用于電算的高斯坐標(biāo)計(jì)算的實(shí)用公式： 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-e2sin2B t=tanB =ecosB l=(L-L0) 其計(jì)算結(jié)果的精度可達(dá)0. 001m。只要得到了高斯平面坐標(biāo)(X, Y)后，便可通過高斯反算公式將其轉(zhuǎn)換成大地坐標(biāo)(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(1-e2sin2Bf)3 Nf=a1-e2sin2Bf f=tanBf 它們的計(jì)算精度,即平面坐標(biāo)可達(dá)0.0001。2.5.3坐標(biāo)轉(zhuǎn)換與參數(shù)計(jì)算轉(zhuǎn)換模型二維轉(zhuǎn)換模型：2個(gè)平移參數(shù)x y （原點(diǎn)不重合產(chǎn)生的）1個(gè)旋轉(zhuǎn)參數(shù)（坐標(biāo)軸不平行產(chǎn)生的）1個(gè)尺度參數(shù)m（兩坐標(biāo)系間的尺度不一致產(chǎn)生的）x2y2=1+mcossin-sincosx1y1+xy(x1,y1)為某點(diǎn)在A空間直角坐標(biāo)系中的坐標(biāo)。(x2,y2)為某點(diǎn)在B空間直角坐標(biāo)系中的坐標(biāo)。(x,y)為某點(diǎn)從A空間直角坐標(biāo)系轉(zhuǎn)換到B空間直角坐標(biāo)系的兩個(gè)平移參數(shù)。為從A空間直角坐標(biāo)系轉(zhuǎn)換到B空間直角坐標(biāo)系中標(biāo)系的兩個(gè)平移參數(shù)。為從A空間直角坐標(biāo)系轉(zhuǎn)換到B空間直角坐標(biāo)系中一個(gè)旋轉(zhuǎn)參數(shù)。m為從A空間直角坐標(biāo)系轉(zhuǎn)換到B空間直角坐標(biāo)系中的一個(gè)尺度參數(shù)。平面四參數(shù)求解步驟如下:利用公共點(diǎn)計(jì)算坐標(biāo)參數(shù)，但至少有兩個(gè)公共點(diǎn)，當(dāng)有i個(gè)公共點(diǎn)時(shí)，可利用最小二乘原理求解參數(shù)。將B直角坐標(biāo)系中的坐標(biāo)視為觀測(cè)值，設(shè)A直角坐標(biāo)系下的坐標(biāo)視為無誤差，列誤差方程為：vx1vy1vxivyiB=10100101y1-x1yi-xix1y1xiyixym-X1Y1XiYiB-X1Y1XiYiA寫成矩陣形式既：V=BX-L由于各點(diǎn)的坐標(biāo)可視為同精度獨(dú)立觀測(cè)值，因此P=I。把各點(diǎn)坐標(biāo)已知值帶入上述誤差方程，然后按下列公式求解出四參數(shù)：X=(BTB)-1(BTL)三維轉(zhuǎn)換模型：3個(gè)平移參數(shù)x y z（原點(diǎn)不重合產(chǎn)生的）3個(gè)旋轉(zhuǎn)參數(shù)x y z（坐標(biāo)軸不平行產(chǎn)生的）1個(gè)尺度參數(shù)m（兩坐標(biāo)系間的尺度不一致產(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)為某點(diǎn)在A空間直角坐標(biāo)系中的坐標(biāo)。(x2,y2z2)為某點(diǎn)在B空間直角坐標(biāo)系中的坐標(biāo)。(x,y,z)為某點(diǎn)從A空間直角坐標(biāo)系轉(zhuǎn)換到B空間直角坐標(biāo)系的三個(gè)平移參數(shù)。(x，y，z)為從A空間直角坐標(biāo)系轉(zhuǎn)換到B空間直角坐標(biāo)系中三個(gè)旋轉(zhuǎn)參數(shù)。m為從A空間直角坐標(biāo)系轉(zhuǎn)換到B空間直角坐標(biāo)系中的一個(gè)尺度參數(shù)。七參數(shù)求解步驟如下:利用公共點(diǎn)計(jì)算坐標(biāo)參數(shù)，但至少有三個(gè)公共點(diǎn)，當(dāng)有i個(gè)公共點(diǎn)時(shí)，可利用最小二乘原理求解參數(shù)。將B空間直角坐標(biāo)系中的坐標(biāo)視為觀測(cè)值，設(shè)A空間直角坐標(biāo)系下的坐標(biāo)視為無誤差，列誤差方程為：vx1vy1vz1vxivyivziB=1001000100100010010Z1-Y10Zi-Yi-Z10X1-Zi0Xi-Y1-X10-Yi-Xi0X1Y1Z1XiYiZixyzxyzm-X1Y1Z1XiYiZiB-X1Y1Z1XiYiZiA寫成矩陣形式既：V=BX-L由于各點(diǎn)的坐標(biāo)可視為同精度獨(dú)立觀測(cè)值，因此P=I。把各點(diǎn)坐標(biāo)已知值帶入上述誤差方程，然后按下列公式求解出七參數(shù)：X=(BTB)-1(BTL)2.5.4三、六度帶帶號(hào)與中央子午線計(jì)算模型求帶號(hào)、中央子午線公式為：3帶:帶號(hào)N=IntL-1.53+1 中央子午線L0=N*3 L為當(dāng)?shù)亟?jīng)度6帶:帶號(hào)N=IntL6+1 中央子午線L0=N*6-3 L為當(dāng)?shù)亟?jīng)度第3章 坐標(biāo)轉(zhuǎn)換程序設(shè)計(jì) 計(jì)劃程序有哪些功能 程序界面設(shè)計(jì) 搜集坐標(biāo)轉(zhuǎn)換公式 運(yùn)用VB編寫程序 檢查能否正常運(yùn)行 以及數(shù)據(jù)正確性 否 是完成程序圖3-1 編寫程序流程圖圖3-2程序主界面3.1大地坐標(biāo)系與空間直角坐標(biāo)系的轉(zhuǎn)換主要代碼如下：橢球參數(shù)設(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.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大地坐標(biāo)轉(zhuǎn)空間直角坐標(biāo)主要代碼 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(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空間直角坐標(biāo)轉(zhuǎn)大地坐標(biāo)主要代碼 x = Text4.Text: y = Text5.Text: z = Text6.Text L = HuToDu(Atn(y / x) 計(jì)算經(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) 用TanB作為未知量 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高斯平面坐標(biāo)與大地坐標(biāo)的轉(zhuǎn)換高斯平面坐標(biāo)(X, Y)與大地坐標(biāo)(B ,L)的相互關(guān)系式分為兩類：第一類稱高斯投影正算公式，由(B,L)求(X,Y)；第二類稱高斯投影反算公式，由(X,Y)求(B,L) 13。由于電子計(jì)算機(jī)和各種可編程序電子計(jì)算器在測(cè)量上廣泛應(yīng)用，因此在這里給出適合電算的高斯正反算公式。主要代碼如下：正算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 輸入中央子午線經(jīng)度yj = Text8.Text Y坐標(biāo)加常數(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) 極點(diǎn)子午線曲率半徑bb0 = 1 - (3 / 4) * e22 + (45 / 64) * 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 = 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(DuToHu(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) - 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(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) / (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 高斯投影換帶計(jì)算高斯投影雖然保證了角度沒有變形，但是長度變形比較嚴(yán)重，為了限制高斯投影的長度變形，必須以中央子午線進(jìn)行分帶，把投影范圍限制在中央子午線東、西側(cè)一定的范圍內(nèi)；但是這樣又使得統(tǒng)一的坐標(biāo)系分割成各帶的獨(dú)立坐標(biāo)系14。這樣就會(huì)產(chǎn)生新的問題，需要鄰帶換算來解決。主要代碼如下：l1 = Text3.Text : yy1 = Text4.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) 極點(diǎn)子午線曲率半徑bb0 = 1 - (3 / 4) * e22 + (45 / 64) * 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(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(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 * 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: yy2 = 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) 極點(diǎn)子午線曲率半徑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 4bbb8 = (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 * (nnn 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 二維坐標(biāo)轉(zhuǎn)換與四參數(shù)計(jì)算由于矩陣計(jì)算代碼繁瑣，由后臺(tái)調(diào)用Excel函數(shù)完成計(jì)算。其主要代碼如下：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) * Cos(DuToHu(jd) / k四參數(shù)計(jì)算：Dim xls As Excel.ApplicationDim book As Excel.WorkbookDim sheet As Excel.WorksheetSet xls = CreateObject(Excel.Application) 創(chuàng)建EXCEL對(duì)象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, 9).Value = x xls.Cells(m + 1, 6).Value = 0: xls.Cells(m + 1, 7).Value = 1: xls.