高斯投影正反算

1、高斯投影正反算學(xué)院：資源與環(huán)境工程工程學(xué)院專業(yè)：測繪工程學(xué)號：X51414012姓名：孫超一、高斯投影概述想象有一個橢圓柱面橫套在地球橢球體外面，并與某一條子午線相切，橢圓柱的中心軸通過橢球體的中心，然后用一定投影方法，將中央子午線兩側(cè)各一定經(jīng)差范圍內(nèi)的地區(qū)投影到橢圓柱面上，再將此柱面展開即成為投影面。高斯投影由于是正形投影，故保證了投影的角度不變性，圖形的相似性以及在某點(diǎn)各方向上長度比的同一性。由于采用了同樣法則的分帶投影，這即限制了長度變形，又保證了在不同投影帶中采用相同的簡便公式和數(shù)表進(jìn)行變形引起的各項改正的計算，并且?guī)c帶間的互相換算也能用相同的公式和方法進(jìn)行。高斯投影的這些優(yōu)點(diǎn)必將使

2、它得到廣泛的推廣和具有國際意義。二、高斯投影坐標(biāo)正算公式1.高斯投影必須滿足以下三個條件1）中央子午線投影后為直線2）中央子午線投影后長度不變3）投影具有正形性質(zhì)，即正形投影條件2.高斯正算公式推導(dǎo)1）由第一個條件可知，由于地球橢球體是一個旋轉(zhuǎn)橢球體，所以高斯投影必然有這樣一個性質(zhì)，即中央子午線東西兩側(cè)的投影必然對稱于中央子午線。2）由于高斯投影是換帶投影，在每帶內(nèi)經(jīng)差是不大的，是一個微小量，所以可以將X=X（，q），Y=Y（l，q）展開為經(jīng)差為l的冪級數(shù)，它可寫成如下的形式X=m0+m2l2+m4l4+Y=m1l+m3l2+m5l5+式中m0,m1,m2,是待定系數(shù)，他們都是緯度B的函數(shù)。3

3、）由第三個條件：yl=xq和xl=-yq，將上式分別對l和q求偏導(dǎo)可得到下式經(jīng)過計算可以得出三、高斯投影坐標(biāo)反算公式推導(dǎo)1.思路：級數(shù)展開，應(yīng)用高斯投影三個條件，待定系數(shù)法求解。2.投影公式在底點(diǎn)處展開展開為3.引入高斯投影條件之一：正形條件4.由于可得到，帶入上式可得到5.引入高斯投影條件之三：中央子午線投影后長度不變四、高斯投影的特點(diǎn)1.當(dāng)l等于常數(shù)時，隨著B的增加x的值增大，y的值減小，無論B值為正或為負(fù)，y值不變。這就是說，橢球面上除中央子午線外，其它子午線投影后，均向中央子午線彎曲，并向兩極收斂，同時還對稱于中央子午線和赤道。2.當(dāng)B等于常數(shù)時，隨著l的增加，x值和y值都增大。所以在

4、橢球面上對稱于赤道的緯圈，投影后仍成為對稱的曲線，同時與子午線的投影曲線互相垂直凹向兩極。3.據(jù)中央子午線越遠(yuǎn)的子午線，投影后彎曲越厲害，長度變形越大。五、MATLAB編程實現(xiàn)坐標(biāo)正反算1.編寫main函數(shù)function maindisp('歡迎使用高斯投影正反算及相鄰帶的坐標(biāo)換算程序');disp('1：高斯正算 2：高斯反算 3：換帶計算');K=0;while (K<1|K>3)K=input('請根據(jù)上列選擇計算類型 K=');switch Kcase 1GSZS;case 2GSFS;case 3HDJS;otherwis

5、edisp('K 值無效 （1-3）');enddisp('程序作者：亞里士多墩');disp('指導(dǎo)老師：亞里士多德');end2.編寫高斯正算GSZS函數(shù)function GSZS%GSZS 是將大地坐標(biāo)換算為高斯坐標(biāo)的子函數(shù)%此函數(shù)要調(diào)用 DHH 和 HHD 兩個子函數(shù)%此函數(shù)包含子午線收斂角的計算disp('你選擇的是高斯正算');B=input('輸入大地坐標(biāo) B=');L=input('輸入大地坐標(biāo) L=');L0=input('輸入所用中央子午線 L0=');B=DHH

6、(B);L=DHH(L);L0=DHH(L0); disp('1:克拉索夫斯基橢球 2:1975 年國際橢球 3:WGS-84 橢球');T=0;while (T<1|T>3)T=input('請根據(jù)上列選擇橢球模型 T=');switch T case 1 a=6378245.0000000000; b=6356863.0187730473; X=111134.861*(B*180/pi)-16036.480*sin(2*B)+16.828*sin(4*B)-0.022*sin(6*B); case 2 a=6378140.0000000000;

7、b=6356755.2881575287; X=111133.005*(B*180/pi)-16038.528*sin(2*B)+16.833*sin(4*B)-0.022*sin(6*B); case 3 a=6378137.0000000000; f=1/298.257223563; b=a*(1-f); X=6367449.1458*B-32009.8185*cos(B)*sin(B)-133.9975*cos(B)*(sin(B)3-0.6975*cos(B)*(sin(B)5; otherwise disp('T 值無效 （1-3）');endende=(sqrt(a

8、2-b2)/a;e1=(sqrt(a2-b2)/b;V=sqrt(1+(e12)*(cos(B)2);c=(a2)/b;M=c/(V3);N=c/V;t=tan(B);n=sqrt(e12)*(cos(B)2);l=L-L0;xp1=X;xp2=(N*sin(B)*cos(B)*l2)/2;xp3=(N*sin(B)*(cos(B)3)*(5-t2+9*n2+4*n4)*l4)/24;xp4=(N*sin(B)*(cos(B)5)*(61-58*t2+t4)*l6)/720;x=xp1+xp2+xp3+xp4;yp1=N*cos(B)*l;yp2=N*(cos(B)3*(1-t2+n2)*l3

9、/6;yp3=N*(cos(B)5*(5-18*t2+t4)*l5/120;y=yp1+yp2+yp3; r1=l*sin(B);r2=(1/3)*sin(B)*(cos(B)2*(l3)*(1+3*n2+2*n4);r3=(1/15)*sin(B)*(cos(B)2*(l5)*(2-t2);r=r1+r2+r3;format long gx yR=HHD(r)End3.編寫高斯反算公式GSFS函數(shù)function GSZS%GSZS 是將大地坐標(biāo)換算為高斯坐標(biāo)的子函數(shù)%此函數(shù)要調(diào)用 DHH 和 HHD 兩個子函數(shù)%此函數(shù)包含子午線收斂角的計算disp('你選擇的是高斯正算')

10、;B=input('輸入大地坐標(biāo) B=');L=input('輸入大地坐標(biāo) L=');L0=input('輸入所用中央子午線 L0=');B=DHH(B);L=DHH(L);L0=DHH(L0); disp('1:克拉索夫斯基橢球 2:1975 年國際橢球 3:WGS-84 橢球');T=0;while (T<1|T>3)T=input('請根據(jù)上列選擇橢球模型 T=');switch T case 1 a=6378245.0000000000; b=6356863.0187730473; X=1111

11、34.861*(B*180/pi)-16036.480*sin(2*B)+16.828*sin(4*B)-0.022*sin(6*B); case 2 a=6378140.0000000000; b=6356755.2881575287; X=111133.005*(B*180/pi)-16038.528*sin(2*B)+16.833*sin(4*B)-0.022*sin(6*B); case 3 a=6378137.0000000000; f=1/298.257223563; b=a*(1-f); X=6367449.1458*B-32009.8185*cos(B)*sin(B)-133.

12、9975*cos(B)*(sin(B)3-0.6975*cos(B)*(sin(B)5; otherwise disp('T 值無效 （1-3）');endende=(sqrt(a2-b2)/a;e1=(sqrt(a2-b2)/b;V=sqrt(1+(e12)*(cos(B)2);c=(a2)/b;M=c/(V3);N=c/V;t=tan(B);n=sqrt(e12)*(cos(B)2);l=L-L0;xp1=X;xp2=(N*sin(B)*cos(B)*l2)/2;xp3=(N*sin(B)*(cos(B)3)*(5-t2+9*n2+4*n4)*l4)/24;xp4=(N*sin(B)*(cos(B)5)*(61-58*t2+t4)*l6)/720;x=xp1+xp2+xp3+xp4;yp1=N*cos(B)*l;yp2=N*(cos(B)3*(1-t2+n2)*l3/6;yp3=N*(cos(B)5*(5-18*t2+t4)*l5/120;y=yp1+yp2+yp3; r1=l*sin(B);r2=(1/3)*sin(B)*(cos(B)2*(l3)*(1+3*n2+2*n4);r3=(1/15)*sin(B)*(cos(B)2*(l5)*(2-t2);r=r1+r