大地测量学基础实验报告.docx
- 文档编号:27977706
- 上传时间:2023-07-07
- 格式:DOCX
- 页数:35
- 大小:1.46MB
大地测量学基础实验报告.docx
《大地测量学基础实验报告.docx》由会员分享,可在线阅读,更多相关《大地测量学基础实验报告.docx(35页珍藏版)》请在冰豆网上搜索。
大地测量学基础实验报告
大地测量学基础
实习报告
学号:
姓名:
班级:
专业:
课程名称:
大地测量学基础
指导老师:
2014年5月
1.引言
《大地测量学基础》的程序设计实验是教学中的重要环节,是测绘工程专业的学生必修的内容。
编程计算在《大地测量学基础》的教学中不可或缺,是大地测量技术理论联系实际的重要环节,可以培养学生的动手能力,增强学生对大地测量数据处理在程序设计、算法阐述,程序实现等方面的认识。
它在培养严谨的治学态度、活跃的创新意识、理论联系实际和适应科技发展的综合应用能力等方面具有其他实践类课程不可替代的作用。
1.1.实习的基本任务
培养基本程序开发技能,提高自身的动手能力,加深对空间大地坐标与空间坐标转化、高斯投影正算与反算等内容的理解。
培养自身的科学思维和创新意识,使自己掌握实验研究的基本方法,提高自身的分析能力和创新能力。
提高自身的科学素养,培养自己理论联系实际和实事求是的科学作风,认真严谨的科学态度,积极主动的探索精神,遵守纪律,团结协作,爱护公共财产的优良品德。
1.2实验要求
(1)分析实习题目的要求。
(2)写出详细设计说明。
(3)编写程序代码,调试程序使其能正确运行。
(4)设计完成的软件要便于操作和使用。
(5)设计完成后提交实习报告。
2.实习内容
2.1实习环境
本程序可以在98/2000/XP下运行,可以用VC++6.0执行。
2.2实习主要内容介绍
2.2.1空间大地坐标转化为空间直角坐标
空间大地坐标以参考椭球的中心为坐标原点,椭球的短轴与参考椭球旋转轴重合,以过地面点的椭球法线与椭球赤道面的夹角为大地纬度B,以过地面点的椭球子午面与起始子午面之间的夹角为大地经度L,地面点沿椭球法线至椭球面的距离为大地高H,地面点的点位用(B,L,H)表示。
空间直角坐标以参心0为坐标原点,X轴与起始子午面和赤道的交线重合,Y轴在赤道面上与X轴垂直,构成右手直角坐标系0-XYZ,Z轴与参考椭球的短轴(旋转轴)相重合,地面点的点位用(X,Y,Z)表示。
在已知椭球的长半轴a、椭球的第一偏心率e,空间大地坐标B、L、H,可以求出空间直角坐标X、Y、Z。
2.2.1.1计算公式
对应计算公式如下:
公式中,N为椭球面卯酉圈的曲率半径,e为椭球的第一偏心率,a、b椭球的长短半径,f椭球扁率,W为第一辅助系数
或
2.2.2空间坐标转化为空间大地坐标
空间直角坐标以参心0为坐标原点,X轴与起始子午面和赤道的交线重合,Y轴在赤道面上与X轴垂直,构成右手直角坐标系0-XYZ,Z轴与参考椭球的短轴(旋转轴)相重合,地面点的点位用(X,Y,Z)表示。
空间大地坐标以参考椭球的中心为坐标原点,椭球的短轴与参考椭球旋转轴重合,以过地面点的椭球法线与椭球赤道面的夹角为大地纬度B,以过地面点的椭球子午面与起始子午面之间的夹角为大地经度L,地面点沿椭球法线至椭球面的距离为大地高H,地面点的点位用(B,L,H)表示。
在已知椭球的长半轴a、椭球的第一偏心率e,空间直角坐标X、Y、Z,可以求出空间大地坐标B、L、H。
2.2.2.1计算公式
对应计算公式如下:
公式中,N为椭球面卯酉圈的曲率半径,e为椭球的第一偏心率,a、b椭球的长短半径,f椭球扁率,W为第一辅助系数
或
2.2.3高斯投影正算
高斯投影正算,已知大地坐标(B,L)及中央子午线经度L0,计算高斯平面坐标(x,y)。
2.2.3.1计算公式
公式中,所用参数如下:
2.2.4高斯投影反算
高斯投影反算,已知高斯平面直角坐标(x,y)及指定中央子午线经度L0,计算大地坐标(B,L)。
2.2.4.1计算公式
公式中,所用参数如下:
3.1空间大地坐标转化为空间直角坐标
3.1.1基本程序
依据空间大地坐标转化为空间直角坐标的公式,计算X、Y、Z,其中椭球选择克拉索夫斯基椭球体,程序如下:
a=6378245;e2=0.006693421622966;
N=a/sqrt(1-e2*sin(radianB)*sin(radianB));
X=(N+H)*cos(radianB)*cos(radianL);
Y=(N+H)*cos(radianB)*sin(radianL);
Z=(N*(1-e2)+H)*sin(radianB);
printf("X=%lf\nY=%lf\nZ=%lf\n",X,Y,Z);
3.1.2程序改进
空间大地坐标转化为空间直角坐标的过程中椭球类型有4种:
克拉索夫斯基椭球体、1975年国际椭球体、WGS-84椭球体、2000中国大地坐标系,通过switch语句进行选择,程序如下:
printf("Yourchooseis\n1(克拉索夫斯基椭球体)\n2(1975年国际椭球体)\n3(WGS-84椭球体)\n4(2000中国大地坐标系):
\n");
scanf("%d",&tuoqiuti);
switch(tuoqiuti)
{
case1:
a=6378245;e2=0.006693421622966;break;
case2:
a=6378140;e2=0.006694384999588;break;
case3:
a=6378137;e2=0.006694379990130;break;
case4:
a=6378137;e2=0.006694380022900;break;
}
N=a/sqrt(1-e2*sin(radianB)*sin(radianB));
X=(N+H)*cos(radianB)*cos(radianL);
Y=(N+H)*cos(radianB)*sin(radianL);
Z=(N*(1-e2)+H)*sin(radianB);
printf("X=%lf\nY=%lf\nZ=%lf\n",X,Y,Z);
3.1.3程序完善
选择椭球类型时可以使用通过switch语句进行选择,程序如下:
#include
#include
#definePI3.14159265358979323846
voidmain()
{
intD1,M1,D2,M2,tuoqiuti,a;
doubleX,Y,Z,e2;
doubleMS1,S1,angleB,radianB,degree1;
doubleMS2,S2,angleL,radianL,degree2;
doubleN,H;
printf("PleaseinputangleB,angleL,H\n");
scanf("%lf%lf%lf",&angleB,&angleL,&H);
/***************************************
*********[纬度转化为弧度]***********
****************************************/
D1=(int)(angleB+0.0001);
MS1=angleB-D1;
M1=(int)((MS1)*100+0.0001);
S1=(MS1*100-M1)*100;
degree1=D1+M1/60.0+S1/3600.0;
radianB=degree1*PI/180.0;
/***************************************
*********[经度转化为弧度]***********
****************************************/
D2=(int)(angleL+0.0001);
MS2=angleL-D2;
M2=(int)(MS2+0.0001)*100;
S2=(MS2*100-M2)*100;
degree2=D2+M2/60.0+S2/3600.0;
radianL=degree2*PI/180.0;
printf("radianB=%lf\nradianL=%lf\n",radianB,radianL);
printf("Yourchooseis\n1(克拉索夫斯基椭球体)\n2(1975年国际椭球体)\n3(WGS-84椭球体)\n4(2000中国大地坐标系):
\n");
/***************************************
***********[P100表4-1]*************
****************************************/
scanf("%d",&tuoqiuti);
switch(tuoqiuti)
{
case1:
a=6378245;e2=0.006693421622966;break;
case2:
a=6378140;e2=0.006694384999588;break;
case3:
a=6378137;e2=0.006694379990130;break;
case4:
a=6378137;e2=0.006694380022900;break;
}
N=a/sqrt(1-e2*sin(radianB)*sin(radianB));
/****************************************
**********[P105公式4-30]************
****************************************/
X=(N+H)*cos(radianB)*cos(radianL);
Y=(N+H)*cos(radianB)*sin(radianL);
Z=(N*(1-e2)+H)*sin(radianB);
printf("X=%lf\nY=%lf\nZ=%lf\n",X,Y,Z);
}
3.1.4程序扩展
选择椭球类型时可以使用到判断、循环等语句,进行选择、重复运算。
3.2空间坐标转化为空间大地坐标
3.2.1基本程序
依据空间坐标转化为空间大地坐标的公式,计算B、L、H,程序如下:
L=atan(Y/X);
N=a/sqrt(1-e2*sin(B)*sin(B));
H=(sqrt(X*X+Y*Y)/cos(B))-N;
B=atan(Z/sqrt(X*X+Y*Y)+N*e2*sin(B)/sqrt(X*X+Y*Y));
printf("B=%lf\nL=%lf\nH=%lf\n",B,L,H);
3.2.2程序改进
依据空间坐标转化为空间大地坐标的公式,计算B、L、H,此程序关键在于迭代计算(如黑色加出部分)程序如下:
B=atan(Z/sqrt(X*X+Y*Y));
printf("B=%lf\n",B);
while(fabs(B-B0)>ERROR)
{
B0=B;
N=a/sqrt(1-e2*sin(B0)*sin(B0));
H=(sqrt(X*X+Y*Y)/cos(B0))-N;
B=atan(Z/sqrt(X*X+Y*Y)+N*e2*sin(B0)/sqrt(X*X+Y*Y));
}
Bf=B;
printf("Bf=%lf\n",Bf);
H=(sqrt(X*X+Y*Y)/cos(Bf))-N;
radianB=atan(Z/sqrt(X*X+Y*Y)+N*e2*sin(Bf)/sqrt(X*X+Y*Y));
degreeB=radianB*180/PI;
D1=(int)(degreeB+0.000001);
MS1=degreeB-D;
M1=(int)(MS1*60);
s1=(MS1*60-M1)*60;
angleB=D1+M1/100.0+s1/10000.0;
printf("angleB=%lf\nangleL=%lf\nH=%lf\n",angleB,angleL,H);
3.2.3程序完善
总述以上基本程序和迭代程序,可得如下程序:
##include
#include
#definePI3.14159265358979323846
#defineERROR1e-10//1e-10为1*10^(-10)
voidmain()
{
inttuoqiuti,a;
intD,M;
intD1,M1;
doubleX,Y,Z,e2;
doubleB,B0,Bf;
doubleH,N;
doubleradianL,degreeL,angleL,MS,s;
doubleradianB,degreeB,angleB,MS1,s1;
printf("PleaseinputX,Y,Z\n");
scanf("%lf%lf%lf",&X,&Y,&Z);
printf("Yourchooseis\n1(克拉索夫斯基椭球体)\n2(1975年国际椭球体)\n3(WGS-84椭球体)\n4(2000中国大地坐标系):
\n");
scanf("%d",&tuoqiuti);
/***************************************
***********[P100表4-1]*************
****************************************/
switch(tuoqiuti)
{
case1:
a=6378245;e2=0.006693421622966;break;
case2:
a=6378140;e2=0.006694384999588;break;
case3:
a=6378137;e2=0.006694379990130;break;
case4:
a=6378137;e2=0.006694380022900;break;
}
/****************************************
**********[P105公式4-31]************
****************************************/
radianL=atan(Y/X);
degreeL=radianL*180/PI;
D=(int)(degreeL+0.000001);
MS=degreeL-D;
M=(int)(MS*60);
s=(MS*60-M)*60;
angleL=D+M/100.0+s/10000.0;
printf("angleL=%lf\n",angleL);
B=atan(Z/sqrt(X*X+Y*Y));
printf("B=%lf\n",B);
/****************************************
******[迭代计算-->Bf公式4-32]********
****************************************/
while(fabs(B-B0)>ERROR)
{
B0=B;
N=a/sqrt(1-e2*sin(B0)*sin(B0));
H=(sqrt(X*X+Y*Y)/cos(B0))-N;
B=atan(Z/sqrt(X*X+Y*Y)+N*e2*sin(B0)/sqrt(X*X+Y*Y));
}
Bf=B;
printf("Bf=%lf\n",Bf);
/****************************************
********[P105公式4-32/4-35]*********
****************************************/
H=(sqrt(X*X+Y*Y)/cos(Bf))-N;
radianB=atan(Z/sqrt(X*X+Y*Y)+N*e2*sin(Bf)/sqrt(X*X+Y*Y));
degreeB=radianB*180/PI;
D1=(int)(degreeB+0.000001);
MS1=degreeB-D;
M1=(int)(MS1*60);
s1=(MS1*60-M1)*60;
angleB=D1+M1/100.0+s1/10000.0;
printf("angleB=%lf\nangleL=%lf\nH=%lf\n",angleB,angleL,H);
}
3.2.4程序扩展
此程序的迭代运算中,常常要使用到判断、循环等语句,进行选择、重复运算,可加以扩展。
3.3高斯投影正算
3.3.1基本程序
根据高斯投影正算公式,由大地坐标(B,L)及中央子午线经度L0,计算高斯平面坐标(x,y),选择椭球类型为WGS-84椭球体,程序如下:
N=a/sqrt(1-e2*sin(radianB)*sin(radianB));
t=tan(radianB);
n=e22*cos(radianB);x=X+N*sin(radianB)*cos(radianB)*(radianj*radianj)/2+N*sin(radianB)*cos(radianB)*cos(radianB)*cos(radianB)*(5-t*t+9*n*n)*(radianj*radianj*radianj*radianj)/24;y=N*cos(radianB)*(radianj)+N*cos(radianB)*cos(radianB)*cos(radianB)*(1-t*t+n*n)*(radianj*radianj*radianj)/6+N*cos(radianB)*cos(radianB)*cos(radianB)*cos(radianB)*cos(radianB)*(5-18*t*t+n*n*n*n)*(radianj*radianj*radianj*radianj*radianj)/120;
printf("x=%lf\ny=%lf\n",x,y);
3.3.2程序改进
高斯投影正算中,椭球类型有4种:
克拉索夫斯基椭球体、1975年国际椭球体、WGS-84椭球体、2000中国大地坐标系,通过switch语句进行选择,程序如下:
printf("Yourchooseis\n1(克拉索夫斯基椭球体)\n2(1975年国际椭球体)\n3(WGS-84椭球体)\n4(2000中国大地坐标系):
\n");
scanf("%d",&tuoqiuti);
switch(tuoqiuti)
{case1:
a=6378245;e2=0.006693421622966;e22=0.006738525414683;break;case2:
a=6378140;e2=0.006694384999588;e22=0.006739501819473;break;case3:
a=6378137;e2=0.006694379990130;e22=0.006739496742270;break;case4:
a=6378137;e2=0.006694380022900;e22=0.006739496775480;break;}
N=a/sqrt(1-e2*sin(radianB)*sin(radianB));
t=tan(radianB);
n=e22*cos(radianB);x=X+N*sin(radianB)*cos(radianB)*(radianj*radianj)/2+N*sin(radianB)*cos(radianB)*cos(radianB)*cos(radianB)*(5-t*t+9*n*n)*(radianj*radianj*radianj*radianj)/24;
y=N*cos(radianB)*(radianj)+N*cos(radianB)*cos(radianB)*cos(radianB)*(1-t*t+n*n)*(radianj*radianj*radianj)/6+N*cos(radianB)*cos(radianB)*cos(radianB)*cos(radianB)*cos(radianB)*(5-18*t*t+n*n*n*n)*(radianj*radianj*radianj*radianj*radianj)/120;
printf("x=%lf\ny=%lf\n",x,y);
3.3.3程序完善
整合以上程序可得如下结果:
#include
#include
#definePI3.14159265358979323846
voidmain()
{
intD1,M1,D2,M2,tuoqiuti,a,n1;
doubleX,e2,e22;
doubleMS1,S1,angleB,radianB,degree1;
doubleMS2,S2,angleL,radianL,degree2,l1,l2;
doubleradianj;
doubleN;
doublem0,m2,m4,m6,m8;
doublea0,a2,a4,a6,a8;
doublex,y;
doublet,n;
printf("PleaseinputangleB,angleL\n");
scanf("%lf%lf",&angleB,&angleL);
/***************************************
*********[纬度转化弧度]************
****************************************/
D1=(int)(angleB+0.0001);
MS1=angleB-D1;
M1=(int)((MS1)*100+0.0001);
S1=(MS1*100-M1)*100;
degree1=D1+M1/60.0+S1/3600.0;
radianB=degree1*PI/180.0;
/***************************************
*********[经度转化弧度]************
****************************************/
D2=(int)(angleL+0.0001);
MS2=angleL-D2;
M2=(int)(MS2+0.0001)*100;
S2=(MS2*100-M2)*100;
degree2=D2+M2/60.0+S2/3600.0;
radianL=degree2*PI/180.0;
/***************************************
***********
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 大地 测量学 基础 实验 报告
![提示](https://static.bdocx.com/images/bang_tan.gif)