高斯投影正反算.doc

高斯投影正反算学院：资源与环境工程工程学院专业：测绘工程学号：X51414012姓名：孙超一、高斯投影概述想象有一个椭圆柱面横套在地球椭球体外面，并与某一条子午线相切，椭圆柱的中心轴通过椭球体的中心，然后用一定投影方法，将中央子午线两侧各一定经差范围内的地区投影到椭圆柱面上，再将此柱面展开即成为投影面。高斯投影由于是正形投影，故保证了投影的角度不变性，图形的相似性以及在某点各方向上长度比的同一性。由于采用了同样法则的分带投影，这即限制了长度变形，又保证了在不同投影带中采用相同的简便公式和数表进行变形引起的各项改正的计算，并且带与带间的互相换算也能用相同的公式和方法进行。高斯投影的这些优点必将使它得到广泛的推广和具有国际意义。二、高斯投影坐标正算公式1.高斯投影必须满足以下三个条件1）中央子午线投影后为直线2）中央子午线投影后长度不变3）投影具有正形性质，即正形投影条件2.高斯正算公式推导1）由第一个条件可知，由于地球椭球体是一个旋转椭球体，所以高斯投影必然有这样一个性质，即中央子午线东西两侧的投影必然对称于中央子午线。2）由于高斯投影是换带投影，在每带内经差是不大的，是一个微小量，所以可以将X=X（，q），Y=Y（l，q）展开为经差为l的幂级数，它可写成如下的形式X=m0+m2l2+m4l4+Y=m1l+m3l2+m5l5+式中m0,m1,m2,是待定系数，他们都是纬度B的函数。3）由第三个条件：yl=xq和xl=-yq，将上式分别对l和q求偏导可得到下式经过计算可以得出三、高斯投影坐标反算公式推导1.思路：级数展开，应用高斯投影三个条件，待定系数法求解。2.投影公式在底点处展开展开为3.引入高斯投影条件之一：正形条件4.由于可得到，带入上式可得到5.引入高斯投影条件之三：中央子午线投影后长度不变四、高斯投影的特点1.当l等于常数时，随着B的增加x的值增大，y的值减小，无论B值为正或为负，y值不变。这就是说，椭球面上除中央子午线外，其它子午线投影后，均向中央子午线弯曲，并向两极收敛，同时还对称于中央子午线和赤道。2.当B等于常数时，随着l的增加，x值和y值都增大。所以在椭球面上对称于赤道的纬圈，投影后仍成为对称的曲线，同时与子午线的投影曲线互相垂直凹向两极。3.据中央子午线越远的子午线，投影后弯曲越厉害，长度变形越大。五、MATLAB编程实现坐标正反算1.编写main函数function maindisp(欢迎使用高斯投影正反算及相邻带的坐标换算程序);disp(1：高斯正算 2：高斯反算 3：换带计算);K=0;while (K3)K=input(请根据上列选择计算类型 K=);switch Kcase 1GSZS;case 2GSFS;case 3HDJS;otherwisedisp(K 值无效 （1-3）);enddisp(程序作者：亚里士多墩);disp(指导老师：亚里士多德);end2.编写高斯正算GSZS函数function GSZS%GSZS 是将大地坐标换算为高斯坐标的子函数%此函数要调用 DHH 和 HHD 两个子函数%此函数包含子午线收敛角的计算disp(你选择的是高斯正算);B=input(输入大地坐标 B=);L=input(输入大地坐标 L=);L0=input(输入所用中央子午线 L0=);B=DHH(B);L=DHH(L);L0=DHH(L0); disp(1:克拉索夫斯基椭球 2:1975 年国际椭球 3:WGS-84 椭球);T=0;while (T3)T=input(请根据上列选择椭球模型 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; 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(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+r2+r3;format long gx yR=HHD(r)End3.编写高斯反算公式GSFS函数function GSZS%GSZS 是将大地坐标换算为高斯坐标的子函数%此函数要调用 DHH 和 HHD 两个子函数%此函数包含子午线收敛角的计算disp(你选择的是高斯正算);B=input(输入大地坐标 B=);L=input(输入大地坐标 L=);L0=input(输入所用中央子午线 L0=);B=DHH(B);L=DHH(L);L0=DHH(L0); disp(1:克拉索夫斯基椭球 2:1975 年国际椭球 3:WGS-84 椭球);T=0;while (T3)T=input(请根据上列选择椭球模型 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; 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(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+r2+r3;format long gx yR=HHD(r)End六、代码测试1.L=1114724.8974，B=310441.6832，L0=1112.L=111.47532575，B=31.23484275，L0=1113.L=114.20，B=30.30，L0=1114.L=118.54152206，B=32.24576522，L0=117