白塞尔大地主题解算
- 格式:doc
- 大小:122.00 KB
- 文档页数:10
大地测量上机实习报告
一、实习时间:第9周星期三
二、实习地点:系办三楼机房
三、实习内容:白塞尔大地主题解算
四、实习目的:(1)尝试编写白塞尔主题解算的程序
(2)通过实例检验解算的结果
五、实习步骤:
1、先了解书本上关于白塞尔大地主题解算的步骤;
2、双击打开白塞尔大地主题解算的程序,如下:
3、在“选择椭球”中选择“克氏椭球”,“选择正反算”中选择“正算”;
4、然后在输入框中输入相应的数据,如下:
5、接下来单击菜单栏中的“数据正确”,记得正算结果,如下所示:
6、重复第3步,选择反算,并输入相应数据:
7、点击“数据正确”,即得反算结果:
六、实习总结:通过白塞尔大地主题结算方法的上机操作,懂得了白塞尔法结算大地主题是将椭球面上的大地元素按照白塞尔投影条件投影到辅助球面上,继而在球面上进行大地主题解算,再将球面上的计算结果换算到椭球面上的原理,并熟悉白塞尔大地主题解算的程序。
白塞尔大地主题解算的基本思想
首先,白塞尔大地主题解算的基本思想之一是建立地壳变形的力学模型。
地壳变形是地球表面的一项重要现象,是由于地质作用和地球内部力
学过程的结果。
地壳变形的力学模型是研究地壳变形的重要方法。
常见的
地壳变形的力学模型有弹性模型、弹塑性模型和拟弹性模型等。
这些模型
可以描述地壳的变形过程,通过对地壳的变形过程进行建模,可以更好地
理解地壳变形的机制和动力学过程。
其次,白塞尔大地主题解算的思想之一是研究现今地球表面的动力学
过程。
地球表面的动力学过程包括板块运动、地震活动、火山喷发等。
这
些过程在地球的长期演化中起着重要的作用。
通过对现今地球表面动力学
过程的研究,可以揭示地球内部的结构和动力学机制,进而更好地理解地
壳变形的成因和发展变化的规律。
第三,白塞尔大地主题解算的基本思想之一是研究地壳变形的成因。
地壳变形的成因是地壳运动的基本原因,也是地球科学研究的一个重要问题。
地壳变形的成因包括构造运动、地壳应力状态的改变、地震活动等。
通过研究地壳变形的成因,可以更好地了解地壳变形的机制和规律,进而
为地震预测和地壳运动的控制提供科学依据。
总之,白塞尔大地主题解算的基本思想是通过建立地壳变形的力学模型,研究现今地球表面的动力学过程,探究地壳变形的成因,揭示地壳变
形与地球动力学过程之间的相互关系,以推动地壳运动的理论和应用研究。
这一思想对于研究地壳运动的机制和规律,了解地壳变形的成因和动力学
过程,提高地震预测和地壳运动控制的水平具有重要意义。
1.垂线偏差:地面一点上的重力向量g和相应椭球面上的法线向量n之间的夹角定义为该点的垂线偏差。
2.参考椭球:具有确定参数(长半径a和扁率α),经过局部定位和定向,同某一地区大地水准面最佳拟合的地球椭球,叫参考椭球。
3.大地线:椭球面上两点间的最短程曲线叫做大地线。
4.力高:水准面在纬度45度处的正常高。
5.大地主题解算:已知某些大地元素推求另一些大地元素的计算工作叫大地主题解算。
6.大地主题正算:已知P1点的大地坐标(L1,B1),P1至P2的大地线长S及其大地方位角,计算P2点的大地坐标(L2,B2)和大地线S在P2点的反方位角A21,这类问题叫做大地主题正算。
7.大地基准:是指能够最佳拟合地球形状的地球椭球的参数及椭球定位和定向8.高斯投影:横轴椭圆柱等角投影(假象有一个椭圆柱横套在地球椭球体外,并与某一条子午线相切,椭球柱的中心轴通过椭球体中心,然后用一定投影方法,将中央子午线两侧各一定范围内的地区投影到椭圆柱上,再将此柱面展开成投影面)。
9.大地测量学:是指在一定的时间与空间参考系中,测量和描绘地球形状及其重力场并监测其变化,为人类活动提供关于地球的空间信息的一门科学。
10.理论闭合差:由水准面不平行而引起的水准环线闭合差,称为理论闭合差。
11.地心坐标系:地心坐标系是在大地体内建立的O-XYZ坐标系。
原点O设在大地体的质量中心,用相互垂直的X,Y,Z三个轴来表示,X轴与首子午面与赤道面的交线重合,向东为正。
Z轴与地球旋转轴重合,向北为正。
Y 轴与XOZ平面垂直构成右手系。
12.高斯投影正、反算公式进行换带计算的步骤。
这种方法的实质是把椭球面上的大地坐标作为过度坐标。
首先把某投影带内有关点的平面坐标(x,y)1利用高斯投影反算公式换算成椭球面上的大地坐标(B,l),进而得到L=L0+l,然后再由大地坐标(B,l),利用投影正算公式换算成相邻带的平面坐标(x,y)2在计算时,要根据第2带的中央子午线来计算经差l,亦即此时l=L-L0。
大地测量实验报告实验名称:白塞尔大地主题解算(正算和反算)实验目的:1.通过编写白塞尔大地主题电算程序进一步掌握白塞尔法解算大地主题的基本思想。
2.熟练掌握将椭球面上的大地元素按照白塞尔投影条件投影到辅助球面上,继而在球面上进行大地主题解算,最后再将球面上的计算结果换算到椭球面上的基本方法和步骤。
3.学会掌握计算机编程的基本能力。
实验环境:Microsoft Visual C++注意事项:1.在编写程序的过程当中要注意代码的前后统一和重复。
2.注意数值类型的转换和度分秒的换算。
实验步骤:正算:1.计算起点的规划纬度2.计算辅助函数值3.计算系数A,B,C及d,e.4计算球面长度5.计算经度差改正数6.计算终点大地坐标及大地方位角。
反算:1.进行计算辅助函数值2.用逐次趋近法同时计算起点大地方位角、球面长度及经差。
3.计算系数A,B,C及大地线长度S.4.计算反方位角及确定符号。
程序源代码:正算:#include<stdio.h>#include<math.h>#define ee 0.006694384999588#define I 3.141592653double F(double,double,double);void main(void){double A1,B1,L1,S,A2,B2,L2; double x1,x2,x3,y1,y2,y3,z1,z2,z3; double W1,sinu1,sinu2,cosu1,sinA0;doublecota1,cos2a1,sin2a1,cosA0A0;double A,B,C,d,e,a0,a1,m;double n,a,Q,R;printf("请输入数据B1= "); scanf("%lf %lf %lf",&x1,&x2,&x3);B1=F(x1,x2,x3);printf("请输入数据L1= "); scanf("%lf %lf %lf",&y1,&y2,&y3);L1=F(y1,y2,y3);printf("请输入A1= ");scanf("%lf %lf %lf",&z1,&z2,&z3);A1=F(z1,z2,z3);printf("请输入S= ");scanf("%lf",&S);printf("B1=%f\n",B1);printf("L1=%f\n",L1);printf("A1=%f\n",A1);printf("S=%f\n",S);/*计算起点的规划纬度*/W1=sqrt(1-ee*sin(B1)*sin(B1));sinu1=sin(B1)*sqrt(1-ee)/W1;cosu1=cos(B1)/W1;printf("W1=%f\n",W1);printf("sinu1=%f\n",sinu1);printf("cosu1=%f\n",cosu1);/*计算辅助函数值*/sinA0=cosu1*sin(A1);cota1=cosu1*cos(A1)/sinu1;sin2a1=2*cota1/(cota1*cota1+1);cos2a1=(cota1*cota1-1)/(cota1*cota1+1);printf("sinA0=%f\n",sinA0);printf("cota1=%f\n",cota1);printf("sin2a1=%f\n",sin2a1);printf("cos2a1=%f\n",cos2a1);/*计算系数ABC及de*/cosA0A0=1-sinA0*sinA0;A=6356755.288+(10710.341-(13.534*cosA0A0))*cosA0A0;B=(5355.171-9.023*cosA0A0)*cosA0A0;C=(2.256*(cosA0A0))*cosA0A0+0.006;d=691.46768-(0.58143-0.00144*cosA0A0)*cosA0A0;e=(0.2907-cosA0A0*0.0010)*cosA0A0;printf("cosA0A0=%f\n",cosA0A0);printf("A=%f\n",A);printf("B=%f\n",B);printf("C=%f\n",C);printf("d=%f\n",d);printf("e=%f\n",e);/*计算球面长度*/a0=(S-(B+C*cos2a1)*sin2a1)/A;m=sin2a1*cos(2*a0)+cos2a1*sin(2*a0);n=(cos2a1)*(cos(2*a0))-(sin2a1)*(sin(2*a0));a=a0+((B+5*C*n))*m/A;printf("a0=%f\n",a0);printf("m=%f\n",m);printf("n=%f\n",n);printf("a=%f\n",a);/*计算经度差改正数*/Q=(d*a+(e*(m-sin2a1))/3600/180*I)*sinA0;printf("Q=%f\n",Q);/*计算终点大地坐标及大地方位角*/sinu2=sinu1*cos(a)+cosu1*cos(A1)*sin(a);B2=180*atan(sinu2/((sqrt(1-ee))*(sq rt(1-sinu2*sinu2))))/I;R=180*atan(sin(A1)*sin(a)/(cosu1*co s(a)-sinu1*sin(a)*cos(A1)))/I;printf("sinu2=%f\n",sinu2);printf("B2=%f\n",B2);printf("R=%f\n",R*180/I);/*确定R的值*/if(sin(A1)>0 && tan(R)>0)R=abs(R);else if(sin(A1)>0 && tan(R)<0)R=I-abs(R);else if(sin(A1)<0 && tan(R)<0)R=-abs(R);elseR=abs(R)-I;/*确定L2A2的值*/L2=(L1*180/I+R-(Q/206265*180/I));A2=atan(cosu1*sin(A1)/(cosu1*cos(a) *cos(A1)-sinu1*sin(a)));if(sin(A1)<0&&tan(A2)>0)A2=(fabs(A2))*180/I;else if(sin(A1)<0&&tan(A2)<0)A2=(I-fabs(A2))*180/I;else if(sin(A1)>0&&tan(A2)>0)A2=(I+fabs(A2))*180/I;elseA2=(2*I-fabs(A2))*180/I;printf("A2=%3f\n B2=%3f\nL2=%3f\n",A2,B2,L2);}double F(double a2,double b2,doublec2){double d2;d2=(double)(a2+1.0*b2/60+1.0*c2/3600);d2=(d2/180)*I;return (d2);}注:A1,B1,L1,S分别为大地线起点的大地方位角、纬度、经度、大地线长;B2,L2,A2为大地线终点纬度、经度及方位角。
白塞尔大地主题解算方向:学号:姓名:一.基本思路:基本思想:将椭球面上的大地元素按照白塞尔投影条件投影到辅助球面上,继而在球面上进行大地主题解算,最后在将球面上的计算结果换算到椭球面上。
其关键问题是找出椭球面上的大地元素与球面上相应元素之间的关系式,同时解决在球面上进行大地主题解算的方法。
正算流程:1.计算起点的归化纬度2.计算辅助函数值,解球面三角形3.按公式计算相关系数A,B,C 以及α,β4.计算球面长度5.计算纬度差改正数6.计算终点大地坐标及大地方位角011122S B C A{sin (cos )}σσσ=-+10101022222sin ()sin sin cos cos σσσσσσ+=+10101022222cos ()cos cos sin sin σσσσσσ+=-001101522B C A[cos ()]sin ()σσσσσσ=++++010122L A sin [(sin ()sin )]λδασβσσσ-==++-2111u u u A sin sin cos cos cos sin σσ=+2222222222222222222222111111111e u B u B W W u e B u u B B arctan e u e u sin sin cos cos tan tan sin sin tan cos -cos ⎧-==⎪⎪⎨⎪=-⎪⎩⎡⎤⎢⎥==--⎢⎥-⎣⎦1111A arctan u u A sin sin []cos cos sin sin cos σλσσ=-21L L λδ=+-112111u A A arctan u A u cos sin cos cos cos sin sin σσ⎡⎤=⎢⎥-⎣⎦反算流程:1.辅助计算2.用逐次趋近法同时计算起点大地方位角、球面长度及经差,第一次趋近时,取δ=0。
计算下式,重复上述计算过程2.3.计算大地线长度S4.计算反方位角二.已知数据序号B1(DD.MMSS)L1 (DD.MMSS)A12(DD.MMSS)S12(m)1 41.01356874 130.10122676 1.4943 8000L λδ=+211212u pA u u u u qsin cos tan cos sin sin cos cos λλ==-2121p p u q b b A arctanqsin cos cos λλ==-=11p A q A sin cos tan cos σσ+=11p A q A sin sin cos σ=+12a a cos cos σλ=+arctan sin cos σσσ⎛⎫=⎪⎝⎭011A u A sin cos sin =111u A tan tan sec σ=21+σσσ=02122L A sin [(sin sin )]λδασβσσ-==+-L λδ=+11222222S A B C B C sin (cos )sin (cos )σσσσσ=++-+1212u A b b cos sin arctan cos λλ⎡⎤=⎢⎥-⎣⎦三.源代码:#include <stdio.h>#include <math.h>#define e 0.081813334016931499 //克拉索夫斯基椭球体第一偏心率void main(){int k,B10,B11,L10,L11,A10,A11,B20,B21,L20,L21,A20,A21;double B12,L12,A12,B22,L22,A22;double B1,L1,A1,S,B2,L2,A2,L,pi;double A,B,C,afa,beta;double a1,a2,b1,b2,p,q,x,y;doubleW1,W2,sinu1,sinu2,cosu1,cosu2,sinA0,cotsigma1,sin2sigma1,cos2sigma1,sigma0,sin2,cos2,sigm a,sins,coss,delta0,delta,lamda;pi=4*atan(1);printf("白塞尔大地主题正算请输入1\n白塞尔大地主题反算请输入2\n");scanf("%d",&k);if(k==1){printf("请输入大地线起点纬度B经度L,大地方位角A及大地线长度S:\n");scanf("%d%d%lf%d%d%lf%d%d%lf%lf",&B10,&B11,&B12,&L10,&L11,&L12,&A10,&A 11,&A12,&S);B1=(B10+(float)B11/60+B12/3600)*pi/180;L1=(L10+(float)L11/60+L12/3600)*pi/180;A1=(A10+(float)A11/60+A12/3600)*pi/180;W1=sqrt(1-e*e*sin(B1)*sin(B1)); //计算起点规划纬度sinu1=sin(B1)*sqrt(1-e*e)/W1; //计算起点规划纬度cosu1=cos(B1)/W1; //计算起点规划纬度sinA0=cosu1*sin(A1); //计算辅助函数值cotsigma1=cosu1*cos(A1)/sinu1; //计算辅助函数值sin2sigma1=2*cotsigma1/(cotsigma1*cotsigma1+1); //计算辅助函数值cos2sigma1=(cotsigma1*cotsigma1-1)/(cotsigma1*cotsigma1+1); //计算辅助函数值A=6356863.020+(10708.949-13.474*(1-sinA0*sinA0))*(1-sinA0*sinA0);B=(5354.469-8.798*(1-sinA0*sinA0))*(1-sinA0*sinA0);C=(2.238*(1-sinA0*sinA0))*(1-sinA0*sinA0)+0.006;afa=691.46768-(0.58143-0.00144*(1-sinA0*sinA0))*(1-sinA0*sinA0);beta=(0.2907-1.0E-3*(1-sinA0*sinA0))*(1-sinA0*sinA0);sigma0=(S-(B+C*cos2sigma1)*sin2sigma1)/A;sin2=sin2sigma1*cos(2*sigma0)+cos2sigma1*sin(2*sigma0);cos2=cos2sigma1*cos(2*sigma0)-sin2sigma1*sin(2*sigma0);sigma=sigma0+(B+5*C*cos2)*sin2/A;delta=(afa*sigma+beta*(sin2-sin2sigma1))*sinA0; //计算经度差改正数delta=delta/3600*pi/180;sinu2=sinu1*cos(sigma)+cosu1*cos(A1)*sin(sigma);B2=atan(sinu2/(sqrt(1-e*e)*sqrt(1-sinu2*sinu2)));lamda=atan(sin(A1)*sin(sigma)/(cosu1*cos(sigma)-sinu1*sin(sigma)*cos(A1))); if(sin(A1)>0){if(tan(lamda)>0)lamda=fabs(lamda);elselamda=pi-fabs(lamda);}else{if(tan(lamda)>0)lamda=fabs(lamda)-pi;elselamda=-1*fabs(lamda);}L2=L1+lamda-delta;A2=atan(cosu1*sin(A1)/(cosu1*cos(sigma)*cos(A1)-sinu1*sin(sigma)));if(sin(A1)>0){if(tan(A2)>0)A2=pi+fabs(A2);elseA2=2*pi-fabs(A2);}else{if(tan(A2)>0)A2=fabs(A2);elseA2=pi-fabs(A2);}B2=B2*180*3600/pi;L2=L2*180*3600/pi;A2=A2*180*3600/pi;B20=(int)B2/3600;B21=(int)B2/60-B20*60;B22=B2-B20*3600-B21*60;L20=(int)L2/3600;L21=(int)L2/60-L20*60;L22=L2-L20*3600-L21*60;A20=(int)A2/3600;A21=(int)A2/60-A20*60;A22=A2-A20*3600-A21*60;printf("正算得到的终点大地经度和大地纬度及A2:\n%d %d %lf\n%d %d %lf\n%d %d %lf\n",B20,B21,B22,L20,L21,L22,A20,A21,A22);}else{printf("请输入大地线起点和终点的坐标BL\n");scanf("%d%d%lf%d%d%lf%d%d%lf%d%d%lf",&B10,&B11,&B12,&L10,&L11,&L12,&B 20,&B21,&B22,&L20,&L21,&L22);B1=(B10+(double)B11/60+B12/3600)*pi/180;L1=(L10+(double)L11/60+L12/3600)*pi/180;B2=(B20+(double)B21/60+B22/3600)*pi/180;L2=(L20+(double)L21/60+L22/3600)*pi/180;W1=sqrt(1-e*e*sin(B1)*sin(B1));W2=sqrt(1-e*e*sin(B2)*sin(B2));sinu1=sin(B1)*sqrt(1-e*e)/W1;sinu2=sin(B2)*sqrt(1-e*e)/W2;cosu1=cos(B1)/W1;cosu2=cos(B2)/W2;L=L2-L1;a1=sinu1*sinu2;a2=cosu1*cosu2;b1=cosu1*sinu2;b2=sinu1*cosu2;delta0=0;lamda=L+delta0;p=cosu2*sin(lamda);q=b1-b2*cos(lamda);A1=atan(p/q);if(p>0){if(q>0)A1=fabs(A1);elseA1=pi-fabs(A1);}else{if(q>0)A1=2*pi-fabs(A1);elseA1=pi+fabs(A1);}sins=p*sin(A1)+q*cos(A1); //计算sigma的正弦值coss=a1+a2*cos(lamda); //计算sigma的余弦值sigma=atan(sins/coss);if(coss>0)sigma=fabs(sigma);elsesigma=pi-fabs(sigma);sinA0=cosu1*sin(A1);x=2*a1-(1-sinA0*sinA0)*cos(sigma);afa=(33523299-(28189-70*(1-sinA0*sinA0))*(1-sinA0*sinA0))*1.0e-10; beta=(28189-94*(1-sinA0*sinA0))*1.0e-10;delta=(afa*sigma-beta*x*sin(sigma))*sinA0;lamda=L+delta;while(fabs(delta-delta0)>4.8e-10){delta0=delta;p=cosu2*sin(lamda);q=b1-b2*cos(lamda);A1=atan(p/q);if(p>0){if(q>0)A1=fabs(A1);elseA1=pi-fabs(A1);}else{if(q>0)A1=2*pi-fabs(A1);elseA1=pi+fabs(A1);}sins=p*sin(A1)+q*cos(A1); //计算sigma的正弦值coss=a1+a2*cos(lamda); //计算sigma的余弦值sigma=atan(sins/coss);if(coss>0)sigma=fabs(sigma);elsesigma=pi-fabs(sigma);sinA0=cosu1*sin(A1);x=2*a1-(1-sinA0*sinA0)*cos(sigma);afa=(33523299-(28189-70*(1-sinA0*sinA0))*(1-sinA0*sinA0))*1.0e-10;beta=(28189-94*(1-sinA0*sinA0))*1.0e-10;delta=(afa*sigma-beta*x*sin(sigma))*sinA0;lamda=L+delta;}A=6356863.020+(10708.949-13.474*(1-sinA0*sinA0))*(1-sinA0*sinA0);B=10708.938-17.956*(1-sinA0*sinA0);C=4.487;y=((1-sinA0*sinA0)*(1-sinA0*sinA0)-2*x*x)*cos(sigma);S=A*sigma+(B*x+C*y)*sin(sigma);A2=atan((cosu1*sin(lamda))/(b1*cos(lamda)-b2));if(sin(A1)>0){if(tan(A2)>0)A2=pi+fabs(A2);elseA2=2*pi-fabs(A2);}else{if(tan(A2)>0)A2=fabs(A2);elseA2=pi-fabs(A2);}A1=A1*3600*180/pi;A2=A2*3600*180/pi;A10=(int)A1/3600;A11=(int)A1/60-A10*60;A12=A1-A10*3600-A11*60;A20=(int)A2/3600;A21=(int)A2/60-A20*60;A22=A2-A20*3600-A21*60;printf("反算得到的方位角A1A2及大地线长S:\n%d %d %lf\n%d %d %lf\n%lf\n",A10,A11,A12,A20,A21,A22,S);}}四.程序执行结果:。
《大地测量学基础》习题与思考题一 绪论1.试述您对大地测量学的理解?2.大地测量的定义、作用与基本内容是什么?3.简述大地测量学的发展概况?大地测量学各发展阶段的主要特点有哪些?4.简述全球定位系统(GPS )、激光测卫(SLR )、 甚长基线干涉测量(VIBL )、 惯性测量系统(INS )的基本概念? 二 坐标系统与时间系统1.简述是开普勒三大行星定律? 2.什么是岁差与章动?什么是极移? 3.什么是国际协议原点 CIO?4.时间的计量包含哪两大元素?作为计量时间的方法应该具备什么条件? 5.恒星时、 世界时、 历书时与协调时是如何定义的?其关系如何? 6.什么是大地测量基准?7.什么是天球?天轴、天极、天球赤道、天球赤道面与天球子午面是如何定义的 ? 8.什么是时圈 、黄道与春分点?什么是天球坐标系的基准点与基准面? 9.如何理解大地测量坐标参考框架?10.什么是椭球的定位与定向?椭球的定向一般应该满足那些条件? 11.什么是参考椭球?什么是总地球椭球?12.什么是惯性坐标系?什么协议天球坐标系 、瞬时平天球坐标系、 瞬时真天球坐标系?13.试写出协议天球坐标系与瞬时平天球坐标系之间,瞬时平天球坐标系与瞬时真天球坐标系的转换数学关系式。
14.什么是地固坐标系、地心地固坐标系与参心地固坐标系?15.什么协议地球坐标系与瞬时地球坐标系?如何表达两者之间的关系?16.如何建立协议地球坐标系与协议天球坐标系之间的转换关系,写出其详细的数学关系式。
17.简述一点定与多点定位的基本原理。
18.什么是大地原点?大地起算数据是如何描述的?19.简述1954年北京坐标系、1980年国家大地坐标系、 新北京54坐标系的特点以及它们之间存在相互关系。
20.什么是国际地球自传服务(IERS )、国际地球参考系统(ITRS) 、国际地球参考框架(ITRF)? ITRS 的建立包含了那些大地测量技术,请加以简要说明?21. 站心坐标系如何定义的?试导出站心坐标系与地心坐标系之间的关系?22.试写出不同平面直角坐标换算、不同空间直角坐标换算的关系式?试写出上述两种坐标转换的误差方程式? 23.什么是广义大地坐标微分方程(或广义椭球变换微分方程)?该式有何作用? 三 地球重力场及地球形状的基本理论1.简述地球大气中平流层、对流层与电离层的概念。
存档日期:存档编号:江苏师范大学科文学院本科生毕业设计(论文)论文题目:贝塞尔大地主题正反算及程序设计*名:**系别:环境与测绘系专业:测绘工程年级、学号:08测绘、************师:***江苏师范大学科文学院教务部印制摘要在大地测量计算过程中,大地主题解算计算繁琐复杂,手工计算易于出错,而且费时费力。
随着计算机技术的高速发展,计算机计算的速度快、准确度高、计算机语言的丰富、编程可视化等优点为我们将复杂烦琐的计算过程简单、简洁、高效化带来了契机。
为了便于工程计算,本课题着眼于研究借助计算机及其编程语言MATLAB来实现大地主题解算问题。
大地主题解算方法,主要有高斯平均引数法、勒让德级数法、贝塞尔法。
前两种方法受到大地线长度的制约,随着大地线两端点的距离加大,其解算精度明显降低。
而贝塞尔法具有不受大地线长度制约的优点,解算精度最大不超过5毫米,是大地主题解算方法中解算精度最高的一种。
因此,本文就以贝塞尔法为研究对象,开发贝塞尔大地主题解算小程序。
关键词:贝塞尔大地主题正反算,程序设计A b s t r a c tIn Geodetic computation process, the solution of geodetic problem computational complexity of manual calculation, error prone, and took the time and trouble. With the rapid development of computer technology, computational speed, high accuracy, computer language, the advantages of rich programming visualization for we will complex complicated calculating process is simple, concise, efficient change brings opportunity. For the convenience of engineering calculation, this paper focus on the research of have the aid of computer and programming language MATLAB to realize the geodetic problem solving.Solution of geodetic problem method, mainly Gauss average argument method, Legendre series expansion method, Bessel method. The former two methods by geodesic length constraints, along with the line ends point distance increase, the calculation precision significantly reduced. Bessel law is not affected by the advantages of geodesic length restriction, calculation accuracy of less than 5mm, is the theme of the earth solution method of calculating precision is highest kind. Therefore, this article on Bessel law as the object of study, the development of Bessel solution of geodetic problem of small procedures.Key words:Direct and inverse solution of geodetic problem,The designing of program目录摘要 .................................................... I Abstract (II)1. 椭球面和球面上对应元素间的关系 (1)1.1 贝塞尔法解算大地问题的基本思想 (1)1.2 对应元素关系式 (1)2. 在球面上进行大地主题解算 (5)2.1 球面上大地主题正解方法 (6)2.2 球面上大地主题反解方法 (6)3. 贝塞尔微分方程的积分 (8)3.1 用于大地主题反算时的大地线长度公式 (8)3.2 用于大地主题正算时的大地线长度公式 (10)3.3 椭球面大地线端点经差与球面经差的关系式 (11)3.4 反解时,大地线长度和球面长度关系式的简化 (13)4. 贝塞尔大地主题正解算步骤 (15)4.1 计算起点的归化纬度 (15)4.2 计算辅助函数值 (15)4.3 计算系数 (15)4.4 计算球面长度 (15)4.5 计算经差改正数 (16)4.6 计算终点大地坐标及大地方位角 (16)5. 贝塞尔大地主题正解算MATLAB 程序设计 (17)5.1 正算流程 (17)5.2 界面设计及功能模块编写 (18)6. 贝赛尔大地主题反解算步骤 (23)6.1 辅助计算 (23)6.2 用逐次趋近法同时计算起点大地方位角、球面长度及经差δλ+=l :236.3 计算系数及大地线长度S (24)6.4 计算反方位角 (24)7. 贝塞尔大地主题反解算MATLAB 程序设计 (25)7.1 反算流程 (25)7.2 界面设计及功能模块编写 (26)8 总结 (29)参考文献 (31)致谢 (29)1. 椭球面和球面上对应元素间的关系1.1 贝塞尔法解算大地问题的基本思想基本思想:将椭球面上的大地元素按照贝塞尔投影条件投影到辅助球面上,继而在球面上进行大地问题解算,最后再将球面上的计算结果换算到椭球面上。
大地测量学知识点第一篇:大地测量学知识点1.大地坐标系:地面点在参考椭圆的位置用大地经度和纬度表示,若地面的点不在椭球面上,它沿法线到椭球面的距离称为大地高2.空间大地直角坐标系:是大地坐标系相应的三维大地直角坐标系3.地心坐标系:定义大地坐标系时,如果选择的旋转椭球为总地球椭球,椭球中心就是地质中心,再定义坐标轴的指向,此时建立的大地坐标系叫做地心坐标系大地方位角:p点的子午面与过p点法线及Q点的平面所成的角度正高系统:地面上一点沿铅垂线到大地水准面的距离正常高系统:一点沿铅垂线到似水准面的距离国家水准网布设的原则:从高级到低级,从整体到局部,分为四个等级布设,逐级控制,逐级加密4.理论闭合差:在闭合的环形水准路线中,由于水准面不平行所产生的闭合差5.大地高系统:地面一点沿法线到椭球面的距离6.平面控制网的测量方法三角测量:在地面上按一定的要求选定一系列的点,他们与周围的邻近点通视,并构成相互联接的三角网状图形,称为三角网,网中各点称为三角点,在各点上可以进行水平角测量,精确观测各三角内角,另外至少精确测量一条三角形边长度D和方位角,作为网的起始边长和起始方位角,推算边长,方位角进而推算各点坐标三边测量:根据三角形的余弦公式,便可求出三角形内角,进而推算出各边的方位角和各点坐标7.国家高程基准的参考面有平均海水面,大地水准面,似大地水准面,参考椭球面1956年黄海高程系统1985年国家高程基准8.角度观测误差分析视准轴误差:视准轴不垂直于水平轴产生水平轴误差:水平轴不垂直于垂直轴产生这2个的消除误差方法为取盘左盘右读数取平均值垂直轴倾斜误差:垂直轴本身偏离铅垂线的位置,即不竖直解决的方法:观测时,气泡不得偏离一格,测回之间重新整理仪器,观测目标的垂直角大于3度,按气泡偏离的格数计算垂直轴倾斜改正9.方向观测法是在一测回内将测站上所有要观测的方向先置盘左位置,逐一照准进行观测,再盘右的位置依次观测,取盘左盘右的平均值作为各方向的观测值。
大地主题解算(正算)代码与白塞尔大地主题解算大地主题解算(正算)代码:根据经纬度和方向角以及距离计算另外一点坐标新建模块->拷贝下面的大地主题(正算)代码,调用方法示例:起点经度:116.235(度)终点纬度:37.435(度)方向角:50(度)长度:500(米)终点经纬度("经度,纬度")=Computation(37.435,116.235,50,500)Const Pi = 3.1415926535898Private a, b, c, alpha, e, e2, W, V As Double'a 长轴半径'b 短轴'c 极曲率半径'alpha 扁率'e 第一偏心率'e2 第二偏心率'W 第一基本纬度函数'V 第二基本纬度函数Private B1, L1, B2, L2 As Double'B1 点1的纬度'L1 点1的经度'B2 点1的纬度'L2 点2的经度Private S As Double '''''大地线长度Private A1, A2 As Double'A1 点1到点2的方位角'A2 点2到点1的方位角Function Computation(STARTLAT, STARTLONG, ANGLE1, DISTANCE As Double) As StringB1 = STARTLATL1 = STARTLONGA1 = ANGLE1S = DISTANCEa = 6378245b = 6356752.3142c = a ^ 2 / balpha = (a - b) / ae = Sqr(a ^ 2 - b ^ 2) / ae2 = Sqr(a ^ 2 - b ^ 2) / bB1 = rad(B1)L1 = rad(L1)A1 = rad(A1)W = Sqr(1 - e ^ 2 * (Sin(B1) ^ 2))V = W * (a / b)Dim W1 As DoubleE1 = e ''''第一偏心率'// 计算起点的归化纬度W1 = W ''Sqr(1 - e1 * e1 * Sin(B1 ) * Sin(B1 )) sinu1 = Sin(B1) * Sqr(1 - E1 * E1) / W1cosu1 = Cos(B1) / W1'// 计算辅助函数值sinA0 = cosu1 * Sin(A1)cotq1 = cosu1 * Cos(A1)sin2q1 = 2 * cotq1 / (cotq1 ^ 2 + 1)cos2q1 = (cotq1 ^ 2 - 1) / (cotq1 ^ 2 + 1)'// 计算系数AA,BB,CC及AAlpha, BBeta的值。
白塞尔大地主题解算。
方向:学号:姓名:\一.基本思路:;基本思想:将椭球面上的大地元素按照白塞尔投影条件投影到辅助球面上,继而在球面上进行大地主题解算,最后在将球面上的计算结果换算到椭球面上。
其关键问题是找出椭球面上的大地元素与球面上相应元素之间的关系式,同时解决在球面上进行大地主题解算的方法。
正算流程:1.计算起点的归化纬度2.计算辅助函数值,解球面三角形3.按公式计算相关系数A,B,C 以及α,β4.计算球面长度#5.计算纬度差改正数6.计算终点大地坐标及大地方位角、011122S B C A{sin (cos )}σσσ=-+10101022222sin ()sin sin cos cos σσσσσσ+=+10101022222cos ()cos cos sin sin σσσσσσ+=-001101522B C A[cos ()]sin ()σσσσσσ=++++010122L A sin [(sin ()sin )]λδασβσσσ-==++-2111u u u A sin sin cos cos cos sin σσ=+22222222222222221111u B u B W u B u B B arctan e u sin cos cos tan sin tan cos ⎧==⎪⎪⎨⎪=⎪⎩⎡⎤==--1111A arctan u u A sin sin []cos cos sin sin cos σλσσ=-21L L λδ=+-反算流程:1.辅助计算2.用逐次趋近法同时计算起点大地方位角、球面长度及经差,第一次趋近时,取δ=0。
!计算下式,重复上述计算过程2.3.计算大地线长度S;4.计算反方位角二.已知数据L λδ=+211212u pA u u u u qsin cos tan cos sin sin cos cos λλ==-2121p p u q b b A arctanqsin cos cos λλ==-=11p A q A sin cos tan cos σσ+=11p A q A sin sin cos σ=+12a a cos cos σλ=+arctan sin cos σσσ⎛⎫=⎪⎝⎭011A u A sin cos sin =111u A tan tan sec σ=21+σσσ=02122L A sin [(sin sin )]λδασβσσ-==+-L λδ=+11222222S A B C B C sin (cos )sin (cos )σσσσσ=++-+…三.源代码:#include <>#include <>#define e //克拉索夫斯基椭球体第一偏心率void main(){int k,B10,B11,L10,L11,A10,A11,B20,B21,L20,L21,A20,A21;double B12,L12,A12,B22,L22,A22;double B1,L1,A1,S,B2,L2,A2,L,pi;【double A,B,C,afa,beta;double a1,a2,b1,b2,p,q,x,y;doubleW1,W2,sinu1,sinu2,cosu1,cosu2,sinA0,cotsigma1,sin2sigma1,cos2sigma1,sigma0,sin2,cos2,sigma ,sins,coss,delta0,delta,lamda;pi=4*atan(1);printf("白塞尔大地主题正算请输入1\n白塞尔大地主题反算请输入2\n");scanf("%d",&k);if(k==1){printf("请输入大地线起点纬度B经度L,大地方位角A及大地线长度S:\n");scanf("%d%d%lf%d%d%lf%d%d%lf%lf",&B10,&B11,&B12,&L10,&L11,&L12,&A10,&A11,&A1 2,&S);`B1=(B10+(float)B11/60+B12/3600)*pi/180;L1=(L10+(float)L11/60+L12/3600)*pi/180;A1=(A10+(float)A11/60+A12/3600)*pi/180;W1=sqrt(1-e*e*sin(B1)*sin(B1)); //计算起点规划纬度sinu1=sin(B1)*sqrt(1-e*e)/W1; //计算起点规划纬度cosu1=cos(B1)/W1; //计算起点规划纬度sinA0=cosu1*sin(A1); //计算辅助函数值cotsigma1=cosu1*cos(A1)/sinu1; //计算辅助函数值sin2sigma1=2*cotsigma1/(cotsigma1*cotsigma1+1); //计算辅助函数值cos2sigma1=(cotsigma1*cotsigma1-1)/(cotsigma1*cotsigma1+1); //计算辅助函数值;A=+ B= C=*(1-sinA0*sinA0))*(1-sinA0*sinA0)+;afa= beta= sigma0=(S-(B+C*cos2sigma1)*sin2sigma1)/A;sin2=sin2sigma1*cos(2*sigma0)+cos2sigma1*sin(2*sigma0);cos2=cos2sigma1*cos(2*sigma0)-sin2sigma1*sin(2*sigma0);sigma=sigma0+(B+5*C*cos2)*sin2/A;delta=(afa*sigma+beta*(sin2-sin2sigma1))*sinA0; //计算经度差改正数delta=delta/3600*pi/180;sinu2=sinu1*cos(sigma)+cosu1*cos(A1)*sin(sigma);B2=atan(sinu2/(sqrt(1-e*e)*sqrt(1-sinu2*sinu2)));lamda=atan(sin(A1)*sin(sigma)/(cosu1*cos(sigma)-sinu1*sin(sigma)*cos(A1))); $if(sin(A1)>0){if(tan(lamda)>0)lamda=fabs(lamda);elselamda=pi-fabs(lamda);}else{if(tan(lamda)>0)~lamda=fabs(lamda)-pi;elselamda=-1*fabs(lamda);}L2=L1+lamda-delta;A2=atan(cosu1*sin(A1)/(cosu1*cos(sigma)*cos(A1)-sinu1*sin(sigma)));if(sin(A1)>0){if(tan(A2)>0)A2=pi+fabs(A2);【elseA2=2*pi-fabs(A2);}else{if(tan(A2)>0)A2=fabs(A2);elseA2=pi-fabs(A2);}~B2=B2*180*3600/pi;L2=L2*180*3600/pi;A2=A2*180*3600/pi;B20=(int)B2/3600;B21=(int)B2/60-B20*60;B22=B2-B20*3600-B21*60;L20=(int)L2/3600;L21=(int)L2/60-L20*60;L22=L2-L20*3600-L21*60;A20=(int)A2/3600;《A21=(int)A2/60-A20*60;A22=A2-A20*3600-A21*60;printf("正算得到的终点大地经度和大地纬度及A2:\n%d %d %lf\n%d %d %lf\n%d %d %lf\n",B20,B21,B22,L20,L21,L22,A20,A21,A22);}else{printf("请输入大地线起点和终点的坐标BL\n");scanf("%d%d%lf%d%d%lf%d%d%lf%d%d%lf",&B10,&B11,&B12,&L10,&L11,&L12,&B20,&B2 1,&B22,&L20,&L21,&L22);B1=(B10+(double)B11/60+B12/3600)*pi/180;L1=(L10+(double)L11/60+L12/3600)*pi/180;、B2=(B20+(double)B21/60+B22/3600)*pi/180;L2=(L20+(double)L21/60+L22/3600)*pi/180;W1=sqrt(1-e*e*sin(B1)*sin(B1));W2=sqrt(1-e*e*sin(B2)*sin(B2));sinu1=sin(B1)*sqrt(1-e*e)/W1;sinu2=sin(B2)*sqrt(1-e*e)/W2;cosu1=cos(B1)/W1;cosu2=cos(B2)/W2;L=L2-L1;a1=sinu1*sinu2;!a2=cosu1*cosu2;b1=cosu1*sinu2;b2=sinu1*cosu2;delta0=0;lamda=L+delta0;p=cosu2*sin(lamda);q=b1-b2*cos(lamda);A1=atan(p/q);if(p>0){"if(q>0)A1=fabs(A1);elseA1=pi-fabs(A1);}else{if(q>0)A1=2*pi-fabs(A1);else)A1=pi+fabs(A1);}sins=p*sin(A1)+q*cos(A1); //计算sigma的正弦值coss=a1+a2*cos(lamda); //计算sigma的余弦值sigma=atan(sins/coss);if(coss>0)sigma=fabs(sigma);elsesigma=pi-fabs(sigma);sinA0=cosu1*sin(A1);>x=2*a1-(1-sinA0*sinA0)*cos(sigma);afa=(-(28189-70*(1-sinA0*sinA0))*(1-sinA0*sinA0))*;beta=(28189-94*(1-sinA0*sinA0))*;delta=(afa*sigma-beta*x*sin(sigma))*sinA0;lamda=L+delta;while(fabs(delta-delta0)>{delta0=delta;p=cosu2*sin(lamda);q=b1-b2*cos(lamda);A1=atan(p/q);if(p>0){if(q>0)A1=fabs(A1);elseA1=pi-fabs(A1);}else{》if(q>0)A1=2*pi-fabs(A1);elseA1=pi+fabs(A1);}sins=p*sin(A1)+q*cos(A1); //计算sigma的正弦值coss=a1+a2*cos(lamda); //计算sigma的余弦值sigma=atan(sins/coss);if(coss>0)sigma=fabs(sigma);%elsesigma=pi-fabs(sigma);sinA0=cosu1*sin(A1);x=2*a1-(1-sinA0*sinA0)*cos(sigma);afa=(-(28189-70*(1-sinA0*sinA0))*(1-sinA0*sinA0))*;beta=(28189-94*(1-sinA0*sinA0))*;delta=(afa*sigma-beta*x*sin(sigma))*sinA0;lamda=L+delta;}A=+ B= C=;~y=((1-sinA0*sinA0)*(1-sinA0*sinA0)-2*x*x)*cos(sigma);S=A*sigma+(B*x+C*y)*sin(sigma);A2=atan((cosu1*sin(lamda))/(b1*cos(lamda)-b2));if(sin(A1)>0){if(tan(A2)>0)A2=pi+fabs(A2);elseA2=2*pi-fabs(A2);}~else{if(tan(A2)>0)A2=fabs(A2);elseA2=pi-fabs(A2);}A1=A1*3600*180/pi;A2=A2*3600*180/pi;A10=(int)A1/3600;A11=(int)A1/60-A10*60;A12=A1-A10*3600-A11*60;A20=(int)A2/3600;A21=(int)A2/60-A20*60;A22=A2-A20*3600-A21*60;printf("反算得到的方位角A1A2及大地线长S:\n%d %d %lf\n%d %d %lf\n%lf\n",A10,A11,A12,A20,A21,A22,S);}}四.程序执行结果:&。