matlab大地测量高斯投影正反算程序设计实验_精品文档.pdf

matlab大地测量高斯投影正反算程序设计实验_精品文档.pdf

《matlab大地测量高斯投影正反算程序设计实验_精品文档.pdf》由会员分享，可在线阅读，更多相关《matlab大地测量高斯投影正反算程序设计实验_精品文档.pdf（10页珍藏版）》请在冰豆网上搜索。

高斯投影正反算高斯投影正反算一、实验目的一、实验目的编写高斯投影正反算的程序，并对格式化文件数据进行计算，验证程序。

二、实验内容：

二、实验内容：

1、高斯投影正算公式高斯投影正算公式已知大地坐标（B,L）及中央子午线经度L0，计算高斯平面坐标（x，y），公式如下：

6222425442232）3302705861）（cos）sin（720）495）（cos）sin（24）cos（）sin（2ltttBBNltBBNlBBNXx52224253223）5814185）（cos120）1）（cos6）cos（ltttBNltBNlBNy其中：

B为纬度，0LLl，单位为弧度，BeaN22sin1，为卯酉圈曲率半径，Bttan,Be222cos，bbae22为第二偏心率，a为旋转椭球长半轴，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椭球参数：

长半轴a=6378137扁率f=1/298.257223563b=a（1-f）c=a*a/b;e2=（a*a-b*b）/（a*a）;-1-e2_=（a*a-b*b）/（b*b）;GaussianMapDirect函数函数如下如下：

functionx,y,L0=GaussianMapDirect（B,L,X）%WGS84椭球参数：

a=6378137;%长半轴f=1/298.257223563;%扁率b=a*（1-f）;%短半轴e=（sqrt（a2-b2）/a;%第一偏心率e_=（sqrt（a2-b2）/b;%第二偏心率%当地中央子午线决定于当地的直角坐标系统，首先确定您的直角坐标系统是3度带还是6度带投影，然后再根据如下公式推算。

Q=input（请选择投影带：

n）;ifQ=6%6度带：

M=round（L+3）./6）;%带号M=round（L+3）/6，即对（L+3）/6的值四舍五入取整数，L为当地经度；L0=6.*M-3;%则中央子午线经度L0=6M-3else%3度带：

M=round（L./3）;%带号M=round（L/3），即对（L/3）的值四舍五入取整数，L为当地经度；L0=3.*M;%则中央子午线经度L0=3Mendl=（L-L0）.*pi/180;N=a./（sqrt（1-（e2）.*（sin（B）.2）;t=tan（B）;p=e_.*cos（B）;%p表示yita-2-%计算高斯平面坐标（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）;endlatitude2meridian函数如下：

函数如下：

functionx=latitude2meridian（B,a,e）%由纬度B求对应的子午线弧长x，计算公式m0=a*（1-e2）;m2=（3*（e2）*m0）/2;m4=（5*（e2）*m2）/4;m6=（7*（e2）*m4）/6;m8=（9*（e2）*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度分秒转度函数如下：

度分秒转度函数如下：

functiondegree=dms2degree（jiaodu）%度分秒（dd.mmss）度-3-degree=fix（jiaodu）;mimute=fix（jiaodu-degree）*100）;second=（jiaodu-degree-mimute/100）*10000;degree=degree+mimute/60+second/3600;end度转度分秒函数如下：

度转度分秒函数如下：

functiondms=degree3dms（du）%度度分秒（dd.mmss）degree=fix（du）;min=fix（du-degree）*60）;second=（du-degree）*60-min）*60）;dms=degree+min/100+second/10000;end度分秒度分秒转转弧度弧度dms2rad函数如下函数如下：

functionrad=dms2rad（n）deg=fix（n）;%度取所给数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;end2、高斯反算公式高斯反算公式已知高斯平面坐标（x，y）及指定中央子午线经度L0，计算大地坐标（B,L）：

64254222232）459061（720）935（242yttNMtyttNMtyNMtBBffffffffffffffff522242532230）8624285（cos1201）21（cos61cos1ytttBNytBNyBNLLfffffffffffff-4-式中，ffBeaN22sin1，3222）sin1（）1（ffBeeaM，ffBe222cos，ffBttan，fB为根据子午线弧长x反算的底点纬度。

创建以GaussianMapInverse命名的函数，函数输入输出格式为B,L=GaussianMapInverse（x,y,RefEllipsoid）或者B,L=GaussianMapInverse（x,y,a,f）GaussianMapInverse函数如下：

函数如下：

functionL,B=GaussianMapInverse（x,y,L0）%WGS84椭球参数：

a=6378137;%长半轴f=1/298.257223563;%扁率b=a*（1-f）;%短半轴e=（sqrt（a2-b2）/a;%第一偏心率e_=（sqrt（a2-b2）/b;%第二偏心率%根据中央子午线弧长x反算底点纬度BfBf=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）./（24.*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）;-5-end子午线弧长反算函数子午线弧长反算函数meridian2latitude如下：

如下：

functionB=meridian2latitude（x,a,e）m0=a*（1-e2）;m2=（3*（e2）*m0）/2;m4=（5*（e2）*m2）/4;m6=（7*（e2）*m4）/6;m8=（9*（e2）*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的初始值while1F=-（a2.*sin（2.*B0）./2+（a4.*sin（4.*B0）./4-（a6.*sin（6.*B0）./6+（a8.*sin（8.*B0）./8;B=（x-F）./a0;ifabs（B0-B）10.-6break;endabs（B0-B）B0=B;endend弧度转度分秒弧度转度分秒rad2dms函数如下函数如下：

functiondms=rad2dms（azimuth）%弧度转度分秒-6-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函数如下函数如下：

functionrad=dms2rad（n）deg=fix（n）;%度取所给数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;end3、实例计算验证实例计算验证首先读取文件data1.txt中的数据，计算其在相应六度带高斯投影带内的高斯平面直角坐标，并存贮在文件data2.txt中，data1.txt格式为：

经度（dd.mmss）纬度（dd.mmss）大地高data2.txt格式为：

x（m）y（m）中央子午线经度（dd.mmss）test61data文件程序如下文件程序如下:

clear;clc;filename,pathname=uigetfile（*.*）;%文件查找窗口file=fullfile（pathname,filename）;%合并路径文件名-7-A=importdata（file）;%将生成的文件导入工作空间，变量名为A%RefEllipsoid为椭球参数%RefEllipsoid=a,b,c,f,e2,e2_;a=6378137;%WGS84椭球参数：

长半轴f=1/298.257223563;%扁率b=a*（1-f）;%WGS84椭球参数：

短半轴e=（sqrt（a2-b2）/a;X=latitude2meridian（dms2rad（A.data（:

2）,a,e）;%X为子午线弧长，有纬度B算出x,y,L0=GaussianMapDirect（dms2rad（A.data（:

2）,dms2degree