matlab大地测量高斯投影正反算程序设计实验

e2_ = (a*a-b*b)/(b*b);
GaussianMapDirect 函数如下： function [x,y,L0] = GaussianMapDirect(B, L, X) %WGS84 椭球参数： a = 6378137 ;%长半轴 f = 1/298.257223563; %扁率 b=a*(1-f) ; %短半轴 e=(sqrt(a^2-b^2))/a;%第一偏心率 e_=(sqrt(a^2-b^2))/b;%第二偏心率 %当地中央子午线决定于当地的直角坐标系统，首先确定您的直角坐标系统是 3 度带还是 6 度带投影，然后再根据如下公式推算。 Q=input('请选择投影带： \n'); if Q==6 %6 度带： M=round((L+3)./6); %带号 M=round[(L+3)/6]，即对(L+3)/6 的值四舍五入
rad=(deg+min/60.00+sec/3600)*pi/180.0; end 2、高斯反算公式 已知高斯平面坐标（ x， y）及指定中央子午线经度 L0，计算大地坐标（ B,L ） ：
B Bf tf
tf 2M f N f
5 f
y2
tf 24M f N 3 f
2 2 2 4 (5 3t 2 f f 9t f f ) y
720M f N
4 6 (61 90t 2 f 45t f ) y
L L0
1 1 2 3 y (1 2t 2 f f ) y 3 N f cos B f 6 N f cos B f
1 4 2 2 2 5 (5 28t 2 f 24t f 6 f 8t f f ) y 120 N cos B f
其中： B 为纬度， l L L0 ，单位为弧度， N a
2
1 e 2 sin 2 B
，为卯酉圈曲率半
' 径， t tan B , 2 e ' cos 2 B ， e
a2 b2 为第二偏心率，a 为旋转椭球长半轴， b
b 为短半轴，X 为子午线弧长。
根据上式创建以 GaussianMapDirect 命名的函数，函数输入输出格式为 [x,y,L0] = GaussianMapDirect(B, L, RefEllipsoid) 或者[x,y,L0] = GaussianMapDirect(B, L, a, f) RefEllipsoid 为椭球参数 RefEllipsoid = [a, b, c, f, e2, e2_]; WGS84 椭球参数：长半轴 扁率 b=a(1-f) c = a*a/b; e2 = (a*a-b*b)/(a*a); a = 6378137 f = 1/298.257223563
-1-
%计算高斯平面坐标（x，y） x=X+(N.*(sin(B)).*(cos(B)).*(l.^2))./2+(N.*(sin(B)).*((cos(B)).^3).*(5-((t).^2)+9.* (p.^2)+4.*(p.^4)).*(l.^4))./24 ... +(N.*sin(B).*(cos(B)).^5.*(61-58.*(t.^2)+(t.^4)+270.*(p.^2)-330.*(p.^2).*(t.^2)).*(l .^6))./720; y=N.*cos(B).*l+(N.*(cos(B)).^3.*(1-t.^2+p.^2).*(l.^3))./6+(N.*((cos(B)).^5).*(5-18 .*(t.^2)+ ... t.^4+14.*(p.^2)-58.*(p.^2).*(t.^2)).*(l.^5))./120; L0=degree3dms(L0); end latitude2meridian 函数如下： function x = latitude2meridian(B,a,e) %由纬度 B 求对应的子午线弧长 x，计算公式 m0=a*(1-e^2); m2=(3*(e^2)*m0)/2; m4=(5*(e^2)*m2)/4; m6=(7*(e^2)*m4)/6; m8=(9*(e^2)*m6)/8; a0=m0+m2/2+3*m4/8+5*m6/16+35*m8/128; a2=m2/2+m4/2+15*m6/32+7*m8/16; a4=m4/8+3*m6/16+7*m8/32; a6=m6/32+m8/16; a8=m8/128; x=a0*B-(a2*sin(2*B))/2+(a4*sin(4*B))/4-(a6*sin(6*B))/6+(a8*sin(8*B))/8; end 度分秒转度函数如下： function degree = dms2degree(jiaodu) %度分秒(dd.mmss) －>度
-4-
end 子午线弧长反算函数 meridian2latitude 如下： function B = meridian2latitude(x,a,e) m0=a*(1-e^2); m2=(3*(e^2)*m0)/2; m4=(5*(e^2)*m2)/4; m6=(7*(e^2)*m4)/6; m8=(9*(e^2)*m6)/8; a0=m0+m2/2+3*m4/8+5*m6/16+35*m8/128; a2=m2/2+m4/2+15*m6/32+7*m8/16; a4=m4/8+3*m6/16+7*m8/32; a6=m6/32+m8/16; a8=m8/128; %%纬度 B 的计算 B0=x./a0; %B 的初始值 while 1
F=-(a2.*sin(2.*B0))./2+(a4.*sin(4.*B0))./4-(a6.*sin(6.*B0))./6+(a8.*sin(8.*B0))./8; B=(x-F)./a0; if abs(B0-B)<10.^-6 break; end abs(B0-B) B0=B; end end 弧度转度分秒 rad2dms 函数如下： function dms=rad2dms(azimuth)%弧度转度分秒
-5-
dgree_tem=azimuth*180/pi; dgree=fix(dgree_tem); min_tem=((dgree_tem-dgree)*60); min=fix(min_tem); second=((min_tem-min)*60); dms=dgree+min/100+second/10000; end 度分秒转弧度 dms2rad 函数如下： function rad=dms2rad(n) deg=fix(n); min_tem=(n-deg)*100; min=fix(min_tem); sec=(min_tem-min)*100; 位的数字 rad=(deg+min/60.00+sec/3600)*pi/180.0; end % 度取所给数 n 的整数部分 %去掉整数部分后小数点后移两位 % 分取整 %秒是小数点再向后移两
取整数， L 为当地经度； L0=6.*M-3 ; else %3 度带： M=round(L./3); 为当地经度； L0=3.*M; end l=(L-L0).*pi/180; N=a./(sqrt(1-(e^2).*(sin(B).^2))); t=tan(B); p=e_.*cos(B); %p 表示 yita %则中央子午线经度 L0=3 × M %带号 M=round(L/3)，即对(L/3)的值四舍五入取整数，L %则中央子午线经度 L0=6 × M-3
高斯投影正反算
一、实验目的 编写高斯投影正反算的程序，并对格式化文件数据进行计算，验证程序。 二、实验内容： 1、高斯投影正算公式 已知大地坐标（B,L）及中央子午线经度 L0，计算高斯平面坐标（ x，y） ，公式 如下：
xX
N N sin( B) cos( B)l 2 sin( B) cos 3 ( B)(5 t 2 9 2 4 4 )l 4 2 24
N sin( B) cos 5 ( B)(61 58t 2 t 4 270 2 330 2 t 2 )l 6 720 N N y N cos( B)l cos 3 ( B)(1 t 2 2 )l 3 cos 5 ( B)(5 18t 2 t 4 14 2 58 2 t 2 )l 5 6 120
-2-
degree = fix(jiaodu); mimute = fix((jiaodu-degree)*100); second = (jiaodu-degree-mimute/100)*10000; degree = degree+mimute/60+second/3600; end 度转度分秒函数如下： function dms=degree3dms(du) degree=fix(du); min=fix((du-degree)*60); second=(((du-degree)*60-min)*60); dms=degree+min/100+second/10000; end 度分秒转弧度 dms2rad 函数如下： function rad=dms2rad(n) deg=fix(n); min_tem=(n-deg)*100; min=fix(min_tem); sec=(min_tem-min)*100; %度取所给数 n 的整数部分 %去掉整数部分后小数点后移两位 %分取整 %秒是小数点再向后移两位的数字 %度－>度分秒 (dd.mmss)
5 f
-3-
式中， Nf a
1 e sin B f
2 2
MfΒιβλιοθήκη Baidu，
2 a(1 e )
(1 e sin B f )
2 2
3
f ，
2
e' cos 2 B f ，
2
t f tan B f ， B f 为根据子午线弧长 x 反算的底点纬度。
创建以 GaussianMapInverse 命名的函数，函数输入输出格式为 [B,L] = GaussianMapInverse(x, y, RefEllipsoid) 或者[B,L] = GaussianMapInverse (x, y, a, f) GaussianMapInverse 函数如下： function [L,B] = GaussianMapInverse(x,y,L0) %WGS84 椭球参数： a = 6378137 ;%长半轴 f = 1/298.257223563; %扁率 b=a*(1-f) ; %短半轴 e=(sqrt(a^2-b^2))/a;%第一偏心率 e_=(sqrt(a^2-b^2))/b;%第二偏心率 %根据中央子午线弧长 x 反算底点纬度 Bf Bf = meridian2latitude(x,a,e); Nf=a./sqrt(1-(e.^2).*(sin(Bf).^2)); Mf=a.*(1-e.^2)./sqrt((1-(e.^2).*(sin(Bf).^2)).^3); pf=e_.*cos(Bf); tf=tan(Bf); %已知高斯平面坐标（ x，y）及指定中央子午线经度 L0 ，计算大地坐标（B,L） B=Bf-tf.*(y.^2)./(2.*Mf.*Nf)+tf.*(5+3.*(tf.^2)+pf.^2-9.*(tf.^2).*(pf.^2)).*(y.^4)./(2 4.*Mf.*(Nf.^3))... +tf.*(61+90.*(tf.^2)+45.*(tf.^4)).*(y.^6)./(720.*Mf.*(Nf.^5)); L=L0+y./(Nf.*cos(Bf))-(1+2.*(tf.^2)+pf.^2).*y.^3./(6.*(Nf.^3).*cos(Bf))... +(5+28.*tf.^2+24.*tf.^4+6.*pf.^2+8.*(tf.^2).*pf.^2).*y.^5./(120.*Nf.^5.*cos(Bf)); B=rad2dms(B); L=rad2dms(L);