大地测量学课程设计.doc

大地测量学课程设计

设计题目：白塞尔大地主题解算

学 院： 矿业学院
专 业： 测绘工程
班 级： ********
学 号： 120801024*
学生姓名： 武 *
指导教师： 张 *

2014年12月31日

目录

1.基本原理及思想 1
2.白塞尔法大地主题正算步骤 2
3.白塞尔法大地主题反算步骤 4
4.同一平行圈弧长、子午线弧长与大地线比较大小 6
5.程序代码 8
6.演算示例 13
7.参考文献 16
8.心得体会 17
9.教师评语 18

3
白塞尔大地主题解算
一：基本原理
建立以椭球中心为中心，以任意长（或单位长）为半径的辅助球，按以下三个步骤计算。
第一，   按一定条件将椭球面元素投影到辅助球面上。
第二，   在球面上解算大地问题。
第三，   将求得的球面元素按投影关系换算到相应的椭球元素。

关键：确定球面元素与椭球面元素的关系，即它们间的投影关系。
二：白塞尔法解算大地主题的基本思想:
以辅助球面为基础,将椭球面三角形转换为辅助球面的相应三角形,由三角形对应元素关系,将椭球面上的大地元素按照白塞尔投影条件投影到辅助球面上，然后在球面上进行大地主题解算，最后再将球面上的计算结果换算到椭球面上。

三：在球面上进行大地主题解算
球面上的大地主题正算：

球面上的大地主题反算：

球面三角元素间的相互关系：

四：正反算步骤
1.白塞尔法大地主题正算步骤
已知 、、、、()、S，计算、、()。
（1）将椭球面元素投影到球面上
由求：
‚计算辅助量和
,
ƒ计算球面长度，将S化为

式中系数分别为：

上式右端含有，因此需要迭代计算。第一次迭代取近似值，第二次计算取

以后计算用代换代入上式迭代计算，直到所要求的精度为止。一般取。
（2） 解算球面三角形
计算

‚计算
或
ƒ计算

（3）将球面元素换算到椭球面上
由求
或
‚将球面经差化为椭球面经差l，求
l
式中

式中的最大值为，故在计算时通常可以略去不计。
ƒ象限的判定
符号
+
+
-
-
符号
+
-
+
-

-

l
符号
-
-
+
+
符号
+
-
+
-

其中、为锐角。
2.白塞尔法大地主题反算步骤
已知 、、、、、，计算()、()、S。
（1）将椭球面元素投影到球面上
由B求u
，，l=
,
,
‚采用逐次趋近法，由l计算
在反算中，已知椭球面上经差l，球面经差上的对应经差未知，为了由l求，由下式可知还需计算、、，计算又还需量，故需要进行迭代计算。
第一次趋近，取l；
，
或
判断的象限
p符号
+
+
-
-
q符号
+
-
+
-

判断象限

+
-

l+
仿照上述计算步骤迭代计算，直到为止。
（2） 将球面元素换算到椭球面上

或

象限的判断与前面一致
五：同一平行圈弧长、子午线弧长与大地线比较大小
子午线弧长计算公式：
式中：
平行圈弧长公式：（）
不同纬度对应的一些弧长的数值
B
子午线弧长
平行圈弧长

1′
1″

1′
1″
0°
110 576m
1 842.94m
30.716m
111 321m
1 855.36m
30.923m
15°
110 656
1 844.26
30.738
107 552
1 792.54
29.876
30°
110 863
1 847.71
30.795
96 488
1 608.13
26.802
45°
111 143
1 852.39
30.873
78 848
1 314.14
21.902
60°
111 423
1 857.04
30.951
55 801
930.02
15.5
75°
111 625
1 860.42
31.007
28 902
481.71
8.028
90°
111 696
1 861.6
31.027
0
0
0
利用白塞尔大地主题反算求解大地线长S
纬度为30°，经差为1°的平行圈弧长S=96 488m，两点间大地线长为96 487.595m

经度为30°，纬度差为1°的子午线弧长X=110 863m，两点间大地线长为110 862.869m

通过比较可知，同一平行圈或同一子午线两点间大地线长度与对应的平行圈弧长或子午线弧长相等。
六：程序代码
#include<stdio.h>
#include<math.h>
double hudu(double,double,double); /*度分秒转换为弧度*/
double du(double); /*弧度转换为度*/
double fen(double); /*弧度转换为分*/
double miao(double); /*弧度转换为秒*/
#define PI 3.1415926
void main (void)
{
int k;
printf("请选择大地主题算法，若执行正算，请输入1；若执行反算，请输入2。\n");
scanf("%d",&k);
/*大地主题正算*/
if(k==1)
{
double ax,ay,az,bx,by,bz,cx,cy,cz,S,dz,ez,fz,B1,B2,L1,L2,A1,A2;
int dx,dy,ex,ey,fx,fy;
double e2,W1,sinu1,cosu1,sinA0,coto1,sin2o1,cos2o1,sin2o,cos2o,A,B,C,r,t,o0,o,g,sinu2,q;
/*输入度分秒数据*/
printf("请输入大地线起点纬度度分秒\n");
scanf("%lf%lf%lf",&ax,&ay,&az);
printf("请输入大地线起点经度度分秒\n");
scanf("%lf%lf%lf",&bx,&by,&bz);
printf("请输入大地方位角度分秒\n");
scanf("%lf%lf%lf",&cx,&cy,&cz);
printf("请输入大地线长度\n");
scanf("%lf",&S);
/*调用函数*/
B1=hudu(ax,ay,az);
L1=hudu(bx,by,bz);
A1=hudu(cx,cy,cz);
/*白塞尔大地主题解算*/
e2=0.006693421622966;
W1=sqrt(1-e2*sin(B1)*sin(B1));
sinu1=sin(B1)*(sqrt(1-e2))/W1;
cosu1=cos(B1)/W1;
sinA0=cosu1*sin(A1);
coto1=cosu1*cos(A1)/sinu1;
sin2o1=2*coto1/(coto1*coto1+1);
cos2o1=(coto1*coto1-1)/(coto1*coto1+1);
A=6356863.020+(10718.949-13.474*(1-sinA0*sinA0))*(1-sinA0*sinA0);
B=(5354.469-8.978*(1-sinA0*sinA0))*(1-sinA0*sinA0);
C=(2.238*(1-sinA0*sinA0))*(1-sinA0*sinA0)+0.006;
r=691.46768-(0.58143-0.00144*(1-sinA0*sinA0))*(1-sinA0*sinA0);
t=(0.2907-0.0010*(1-sinA0*sinA0))*(1-sinA0*sinA0);
o0=(S-(B+C*cos2o1)*sin2o1)/A;
sin2o=sin2o1*cos(2*o0)+cos2o1*sin(2*o0);
cos2o=cos2o1*cos(2*o0)-sin2o1*sin(2*o0);
o=o0+(B+5*C*cos2o)*sin2o/A;
g=(r*o+t*(sin2o-sin2o1))*sinA0;
/*求B2*/
sinu2=sinu1*cos(o)+cosu1*cos(A1)*sin(o);
B2=atan(sinu2/(sqrt(1-e2)*sqrt(1-sinu2*sinu2)));
/*求L2*/
q=atan(sin(A1)*sin(o)/(cosu1*cos(o)-sinu1*sin(o)*cos(A1)));
/*判断q*/
if(sin(A1)>0 && tan(q)>0)
q=fabs(q);
else if(sin(A1)>0 && tan(q)<0)
q=PI-fabs(q);
else if(sin(A1)<0 && tan(q)<0)
q=-fabs(q);
else
q=fabs(q)-PI;
L2=L1+q-g/3600/180*PI;
/*求A2*/
A2=atan(cosu1*sin(A1)/(cosu1*cos(o)*cos(A1)-sinu1*sin(o)));
/*判断A2*/
if(sin(A1)<0 && tan(A2)>0)
A2=fabs(A2);
else if(sin(A1)<0 && tan(A2)<0)
A2=PI-fabs(A2);
else if(sin(A1)>0 && tan(A2)>0)
A2=PI+fabs(A2);
else
A2=2*PI-fabs(A2);
/*调用函数*/
dx=(int)(du(B2));
dy=(int)(fen(B2));
dz=miao(B2);
ex=(int)(du(L2));
ey=(int)(fen(L2));
ez=miao(L2);
fx=(int)(du(A2));
fy=(int)(fen(A2));
fz=miao(A2);
printf("大地线终点纬度度分秒分别为：\n%d\n%d\n%lf\n",dx,dy,dz);
printf("大地线终点经度度分秒分别为：\n%d\n%d\n%lf\n",ex,ey,ez);
printf("终点大地方位角度分秒分别为：\n%d\n%d\n%lf\n",fx,fy,fz);
}
/*大地主题反算*/
else if(k==2)
{
double fx,fy,fz,gx,gy,gz,hx,hy,hz,ix,iy,iz,jz,kz,B1,B2,L1,L2,S,A1,A2;
int jx,jy,kx,ky;
double e2,W1,W2,sinu1,sinu2,cosu1,cosu2,L,a1,a2,b1,b2,g,g2,g0=0.0,r,p,q,sino,coso,o,sinA0,x,t1,t2,A,B,C,y;
/*输入度分秒数据*/
printf("请输入大地线起点纬度度分秒\n");
scanf("%lf%lf%lf",&fx,&fy,&fz);
printf("请输入大地线起点经度度分秒\n");
scanf("%lf%lf%lf",&gx,&gy,&gz);
printf("请输入大地线终点纬度度分秒\n");
scanf("%lf%lf%lf",&hx,&hy,&hz);
printf("请输入大地线终点经度度分秒\n");
scanf("%lf%lf%lf",&ix,&iy,&iz);
/*调用函数*/
B1=hudu(fx,fy,fz);
L1=hudu(gx,gy,gz);
B2=hudu(hx,hy,hz);
L2=hudu(ix,iy,iz);
/*白塞尔大地主题解算*/
e2=0.006693421622966;
W1=sqrt(1-e2*sin(B1)*sin(B1));
W2=sqrt(1-e2*sin(B2)*sin(B2));
sinu1=sin(B1)*sqrt(1-e2)/W1;
sinu2=sin(B2)*sqrt(1-e2)/W2;
cosu1=cos(B1)/W1;
cosu2=cos(B2)/W2;
L=L2-L1;
a1=sinu1*sinu2;
a2=cosu1*cosu2;
b1=cosu1*sinu2;
b2=sinu1*cosu2;
/*逐次趋近法求解A1*/
g=0;
r=L;
while(1)
{
p=cosu2*sin(r);
q=b1-b2*cos(r);
A1=atan(p/q);
/*判断A1*/
if(p>0 && q>0)
A1=fabs(A1);
else if(p>0 && q<0)
A1=PI-fabs(A1);
else if(p<0 && q<0)
A1=PI+fabs(A1);
else
A1=2*PI-fabs(A1);
sino=p*sin(A1)+q*cos(A1);
coso=a1+a2*cos(r);
o=atan(sino/coso);
/*判断o*/
if(coso>0)
o=fabs(o);
else
o=PI-fabs(o);
sinA0=cosu1*sin(A1);
x=2*a1-(1-sinA0*sinA0)*coso;
t1=(33523299-(28189-70*(1-sinA0*sinA0))*(1-sinA0*sinA0))*1e-10;
t2=(28189-94*(1-sinA0*sinA0))*1e-10;
g2=(t1*o-t2*x*sino)*sinA0;
if(fabs(g2-g0)<=1e-100)
break;
else
{
r=L+g2;
g0=g2;
}
}
/*求解S*/
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)*coso;
S=A*o+(B*x+C*y)*sino;
/*求解A2*/
A2=atan(cosu1*sin(r)/(b1*cos(r)-b2));
/*判断A2*/
if(p<0 && q<0)
A2=fabs(A2);
else if(p<0 && q>0)
A2=PI-fabs(A2);
else if(p>0 && q>0)
A2=PI+fabs(A2);
else
A2=2*PI-fabs(A2);
/*调用函数*/
jx=(int)(du(A1));
jy=(int)(fen(A1));
jz=miao(A1);
kx=(int)(du(A2));
ky=(int)(fen(A2));
kz=miao(A2);
printf("起点大地方位角度分秒分别为：\n%d\n%d\n%lf\n",jx,jy,jz);
printf("终点大地方位角度分秒分别为：\n%d\n%d\n%lf\n",kx,ky,kz);
printf("大地线长度为：%lf\n",S);
}
/*数据错误*/
else
printf("数据错误，请重新输入\n");
}
/*度分秒转换为弧度*/
double hudu(double x,double y,double z)
{
double A0;
A0=(x+y/60+z/3600)*PI/180;
return A0;
}
/*弧度转换为度*/++++++++++++++
double du(double B0)
{
double x0;
x0=(int)(B0*180/PI);
return x0;
}
/*弧度转换为分*/
double fen(double C0)
{
double _y,y0;
_y=(int)(C0*180/PI);
y0=(fabs)((int)((C0*180/PI-_y)*60));
return y0;
}
/*弧度转换为秒*/
double miao(double D0)
{
double _z1,_z2,z0;
_z1=(int)(D0*180/PI);
_z2=(int)((D0*180/PI-_z1)*60);
z0=(fabs)((double)(((D0*180/PI-_z1)*60-_z2)*60));
return z0;
}

七：演算示例
正算

反算

正算计算结果

30°30′00.00″

-0.505 246 786

114°20′00.00″

-0.610 640 770

225°00′00.00″

-0.707 105 781
S
10000000.00

0.356 537 888

0.530858149

0.003 350 562

-0.609789472

1.766 987 874E-06

0.792563224

4.622 109 698E-10

-0.692911606

-1.103 777 280

1.571 441 354

5.270 424 445E-03

1.055 977 930

-37°43′44.1351″

1.393 840 997

51°16′32.4976″

1.572 478 988

50°21′22.4896″

-0.609 789 747

反算计算结果

0.530 858 149

-0.656 869 486

-0.309 151 284

0.682 919 788

-0.526 600 070

0.400 921 954
计算值
趋近次数
1
2
3
4

-1.100 563 429
-1.103 768 055
-1.103 777 253
-1.103 777 280
p
-0.705 956 271
-0.707 102 492
-0.707 105 771
-0.707 105 781
q
-0.708 255 369
-0.707 109 082
-0.707 105 789
-0.707 105 780

224°54′24.6737″
224°59′59.0390″
224°59′59.9973″
225°00′00.0002″

0.999 999 996
0.999 998 594
0.999 998 584
0.999 998 584

2.755 182 939E-04
-1.677 035 942E-04
-1.682 644 780E-04
-1.682 660 878E-04

1.570 520 809
1.572 473 364
1.572 478 973
1.572 478 989

0.608 797 601
0.609 786 905
0.609 789 739
0.609 789 747

-0.692 114 035
-0.692 909 315
-0.692 911 599
-0.692 911 606

3.350 558 481E-03
3.350 561 868E-03
3.350 561 878E-03
3.350 561 878E-03

1.770 381 674E-6
1.766 997 603E-6
1.766 987 901E-6
1.766 987 874E-6

4.639 919 221E-10
4.622 160 703E-10
4.622 109 841E-10
4.622 109 695E-10
xx
5.263 861 564
5.270 405 594
5.270 424 392
5.270 424 446

-1.103 768 055
-1.103 777 253
-1.103 777 280
-1.103 777 280

1.571 441 354

-0.609789747

1.055 977 929

-0.505246786

1.393 840 996

50°21′22.4897″
Y
-4.367 643 449E-10
S
10000000.003

八：参考文献
[1]史国友,周晓明,贾传荧. 贝塞尔大地主题正解的改进算法[J]. 大连海事大学学报,2008,01:15-19.
[2]郭际明,丁士俊,苏新洲,刘宗泉.大地测量学基础实践教程.武汉:武汉大学出版社,2009
[3]裴连磊. 用C语言实现大地主题解算[J]. 价值工程,2013,20:235-236.
[4]周振宇,郭广礼,贾新果. 大地主题解算方法综述[J]. 测绘科学,2007,04:190-191+174+200.
[5]丁士俊,杨艳梅,史俊波,程新明. 大地主题解算几种不同算法在计算中应注意的问题[J]. 黑龙江工程学院学报(自然科学版),2013,03:1-5.
[6]王建强,胡明庆. 贝赛尔大地主题解算分析[J]. 测绘科学,2012,01:30-31.
[7]徐晓晗,谢云开,李亚军. 大地主题解算实用算法[J]. 科学技术与工程,2012,09:2062-2068.
[8]许厚泽. 关于正反大地主题解算方法的综合研究[J]. 测量制图学报,1958,04:274-288.
.

心得体会

学习课程设计是培养学生综合运用所学知识,发现,提出,分析和解决实际问题,锻炼实践能力的重要环节,是对学生实际工作能力的具体训练和考察过程.随着科学技术发展的日新日异，测绘技术也不断发展，这次课程设计让我进一步了解和掌握大地线是椭球面上最短程曲线的特性以及大地主题解算在测绘学科中的地位和在实际应用中的意义，同时，我也学会利用计算机解决问题的基本技能得到训练和提高，为我走向工作岗位打下一定实践基础。
这次课程设计最难得地方就是程序设计，虽然之前学过C语言，但是只是懂得一些基础，而且之前学习了的c语言也忘记的差不多了，所以我又重新开始看书，复习了一下c语言才开始做课程设计在这次课程设计中，我发现教材上有一些印刷错误，虽然问题不大，但对于还不熟悉这门课的我们来说，在学习过程中难免会造成疑惑，多走弯路。此外，教材上的算法虽然简单易懂，但是它在计算系数时是用椭球参数直接代入后化简的公式，一个程序只能依据一个参考椭球，不够灵活。我在此基础上稍微做了修改，可以以不同参考椭球为依据，使程序更具普遍性，让课程设计更完善。
回想起起此次课程设计，至今我仍感慨颇多，从理论到实践，在这段日子里，可以说得是苦多于甜，但是可以学到很多很多的东西，同时不仅可以巩固了以前所学过的知识，而且学到了很多在书本上所没有学到过的知识。通过这次课程设计使我懂得了理论与实际相结合是很重要的，只有理论知识是远远不够的，只有把所学的理论知识与实践相结合起来，从理论中得出结论，才能真正学以致用，从而提高自己的实际动手能力和独立思考的能力。在设计的过程中遇到问题，可以说得是困难重重，但可喜的是最终都得到了解决。
虽然课程设计已经结束了，但是在这次的课程设计中不仅检验了我所学习的知识，也培养了我如何去把握一件事情，如何去做一件事情，又如何完成一件事情。学会了合作，学会了运筹帷幄。课程设计是我们专业课程知识综合应用的实践训练，着是我们迈向社会，从事职业工作前一个必不少的过程．”千里之行始于足下”，通过这次课程设计，我深深体会到这句千古名言的真正含义．我今天认真的进行课程设计，学会脚踏实地迈开这一步，就是为明天能稳健的在社会大潮中走的更远。我会继续努力专业知识，将自己所学学以致用。

十：教师评语