高斯投影正反算

高斯投影正反算

高斯投影正反算

学院：资源与环境工程工程学院

专业：测绘工程 学号：X51414012 姓名：孙超

一、高斯投影概述

想象有一个椭圆柱面横套在地球椭球体外面，并与某一条子午线相切，椭圆柱的中心轴通过椭球体的中心，然后用一定投影方法，将中央子午线两侧各一定经差范围内的地区投影到椭圆柱面上，再将此柱面展开即成为投影面。高斯投影由于是正形投影，故保证了投影的角度不变性，图形的相似性以及在某点各方向上长度比的同一性。由于采用了同样法则的分带投影，这即限制了长度变形，又保证了在不同投影带中采用相同的简便公式和数表进行变形引起的各项改正的计算，并且带与带间的互相换算也能用相同的公式和方法进行。高斯投影的这些优点必将使它得到广泛的推广和具有国际意义。

二、高斯投影坐标正算公式

1.高斯投影必须满足以下三个条件 1）中央子午线投影后为直线 2）中央子午线投影后长度不变

3）投影具有正形性质，即正形投影条件 2.高斯正算公式推导

1）由第一个条件可知，由于地球椭球体是一个旋转椭球体，所以高斯投影必然有这样一个性质，即中央子午线东西两侧的投影必然对称于中央子午线。 2）由于高斯投影是换带投影，在每带内经差ｌ是不大的，𝜌是一个微小量，所以可以将

ｌX=X（ｌ，q），Y=Y（l，q）

展开为经差为l的幂级数，它可写成如下的形式

X=m0+m2l+m4l+…

24

Y=m1l+m3l+m5l+…

式中m0,m1,m2,…是待定系数，他们都是纬度B的函数。 3）由第三个条件：∂l=∂q和∂l=-∂q，将上式分别对l和q求偏导

xm0m1lm2l2m3l3m4l4.....yn0n1ln2ln3ln4l......23425

∂y∂x∂x∂y

可得到下式

dm01dm11dm21dm3n,n,n,n1dq22dq33dq44dq,L mdn0,m1dn1,m1dn2,m1dn3,L1234dq2dq3dq4dq经过计算可以得出

NNsinBcosBl2sinBcos3B(5t29244)l4224N sinBcos5B(6158t2t4)l6720NyNcosBlcos3B(1t22)l3 6N cos5B(518t2t4142582t2)l5120xX三、高斯投影坐标反算公式推导

1.思路：级数展开，应用高斯投影三个条件，待定系数法求解。 2.投影公式在底点处展开

qf1(x,y)lf2(x,y)展开为

m1ym2y2m3y3m4y4.....qm0n1yn2yn3yn4y......ln0234

3.引入高斯投影条件之一：正形条件

0，带入上式可得到 4.由于可得到n0m2y2m4y4.....qm0yn3yn5y......ln135

5.引入高斯投影条件之三：中央子午线投影后长度不变

qfm0secBfndm01dxNftfsecBf1dn1m222dx2NfsecBf22n(12t3ff)36NftsecBf22f(56tm4ff)424NfsecBf22n(528t65ff)5120Nf

四、高斯投影的特点

1.当l等于常数时，随着B的增加x的值增大，y的值减小，无论B值为正或为负，y值不变。这就是说，椭球面上除中央子午线外，其它子午线投影后，均向中央子午线弯曲，并向两极收敛，同时还对称于中央子午线和赤道。 2.当B等于常数时，随着l的增加，x值和y值都增大。所以在椭球面上对称于赤道的纬圈，投影后仍成为对称的曲线，同时与子午线的投影曲线互相垂直凹向两极。

3.据中央子午线越远的子午线，投影后弯曲越厉害，长度变形越大。

五、MATLAB编程实现坐标正反算

1.编写main函数 function main

disp('欢迎使用高斯投影正反算及相邻带的坐标换算程序'); disp('1：高斯正算 2：高斯反算 3：换带计算'); K=0;

while (K<1||K>3)

K=input('请根据上列选择计算类型 K='); switch K case 1 GSZS; case 2 GSFS; case 3 HDJS; otherwise

disp('K 值无效 （1-3）'); end

disp('程序作者：亚里士多墩'); disp('指导老师：亚里士多德'); end

2.编写高斯正算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 (T<1||T>3)

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）'); end end

e=(sqrt(a^2-b^2))/a; e1=(sqrt(a^2-b^2))/b; V=sqrt(1+(e1^2)*(cos(B))^2); c=(a^2)/b; M=c/(V^3); N=c/V; t=tan(B);

n=sqrt((e1^2)*(cos(B))^2); l=L-L0; xp1=X;

xp2=(N*sin(B)*cos(B)*l^2)/2;

xp3=(N*sin(B)*((cos(B))^3)*(5-t^2+9*n^2+4*n^4)*l^4)/24; xp4=(N*sin(B)*((cos(B))^5)*(61-58*t^2+t^4)*l^6)/720; x=xp1+xp2+xp3+xp4; yp1=N*cos(B)*l;

yp2=N*(cos(B))^3*(1-t^2+n^2)*l^3/6; yp3=N*(cos(B))^5*(5-18*t^2+t^4)*l^5/120; y=yp1+yp2+yp3; r1=l*sin(B);

r2=(1/3)*sin(B)*(cos(B))^2*(l^3)*(1+3*n^2+2*n^4); r3=(1/15)*sin(B)*(cos(B))^2*(l^5)*(2-t^2); r=r1+r2+r3; format long g x

y R=HHD(r) End

3.编写高斯反算公式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 (T<1||T>3)

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）'); end end

e=(sqrt(a^2-b^2))/a; e1=(sqrt(a^2-b^2))/b; V=sqrt(1+(e1^2)*(cos(B))^2); c=(a^2)/b; M=c/(V^3); N=c/V; t=tan(B);

n=sqrt((e1^2)*(cos(B))^2); l=L-L0; xp1=X;

xp2=(N*sin(B)*cos(B)*l^2)/2;

xp3=(N*sin(B)*((cos(B))^3)*(5-t^2+9*n^2+4*n^4)*l^4)/24; xp4=(N*sin(B)*((cos(B))^5)*(61-58*t^2+t^4)*l^6)/720; x=xp1+xp2+xp3+xp4;

yp1=N*cos(B)*l;

yp2=N*(cos(B))^3*(1-t^2+n^2)*l^3/6; yp3=N*(cos(B))^5*(5-18*t^2+t^4)*l^5/120; y=yp1+yp2+yp3; r1=l*sin(B);

r2=(1/3)*sin(B)*(cos(B))^2*(l^3)*(1+3*n^2+2*n^4); r3=(1/15)*sin(B)*(cos(B))^2*(l^5)*(2-t^2); r=r1+r2+r3; format long g x y R=HHD(r) End

六、代码测试

1.L=111°47'24〞.8974，B=31°04'41〞.6832，L0=111°

2.L=111.47532575，B=31.23484275，L0=111

3.L=114.20，B=30.30，L0=111

4.L=118.54152206，B=32.24576522，L0=117