大地测量学课程设计报告

...---可修编.测量学课程设计设计题目：白塞尔主题解算学院：矿业学院专业：测绘工程班级：********学号：120801024*学生：武*指导教师：*2014年12月31日..--可修编.目录1.根本原理及思想12.白塞尔法主题正算步骤23.白塞尔法主题反算步骤44.同一平行圈弧长、子午线弧长与线比拟大小65.程序代码86.演算例如137.参考文献168.心得体会179.教师评语18..--可修编.白塞尔主题解算一：根本原理建立以椭球中心为中心，以任意长〔或单位长〕为半径的辅助球，按以下三个步骤计算。第一，按一定条件将椭球面元素投影到辅助球面上。第二，在球面上解算问题。第三，将求得的球面元素按投影关系换算到相应的椭球元素。关键：确定球面元素与椭球面元素的关系，即它们间的投影关系。二：白塞尔法解算主题的根本思想:以辅助球面为根底,将椭球面三角形转换为辅助球面的相应三角形,由三角形对应元素关系,将椭球面上的元素按照白塞尔投影条件投影到辅助球面上，然后在球面上进展主题解算，最后再将球面上的计算结果换算到椭球面上。三：在球面上进展主题解算球面上的主题正算：球面上的主题反算：球面三角元素间的相互关系：四：正反算步骤1.白塞尔法主题正算步骤、、、、()、S，计算、、()。〔1〕将椭球面元素投影到球面上由求：计算辅助量和,计算球面长度，将S化为式中系数分别为：上式右端含有，因此需要迭代计算。第一次迭代取近似值，第二次计算取以后计算用代换代入上式迭代计算，直到所要求的精度为止。一般取。解算球面三角形计算计算或计算〔3〕将球面元素换算到椭球面上由求或将球面经差化为椭球面经差l，求l式中式中的最大值为，故在计算时通常可以略去不计。象限的判定符号++--符号+-+--l 符号--++符号+-+-其中、为锐角。2.白塞尔法主题反算步骤、、、、、，计算()、()、S。〔1〕将椭球面元素投影到球面上由B求u，，l=,,采用逐次趋近法，由l计算在反算中，椭球面上经差l，球面经差上的对应经差未知，为了由l求，由下式可知还需计算、、，计算又还需量，故需要进展迭代计算。第一次趋近，取l；，或判断的象限p符号++--q符号+-+-判断象限+-l+仿照上述计算步骤迭代计算，直到为止。将球面元素换算到椭球面上或象限的判断与前面一致五：同一平行圈弧长、子午线弧长与线比拟大小子午线弧长计算公式：式中：平行圈弧长公式：〔〕不同纬度对应的一些弧长的数值B子午线弧长平行圈弧长1′1″1′1″0°110576m1842.94m30.716m111321m1855.36m30.923m15°1106561844.2630.7381075521792.5429.87630°1108631847.7130.795964881608.1326.80245°1111431852.3930.873788481314.1421.90260°1114231857.0430.95155801930.0215.575°1116251860.4231.00728902481.718.02890°1116961861.631.027000利用白塞尔主题反算求解线长S纬度为30°，经差为1°的平行圈弧长S=96488m，两点间线长为96487.595m经度为30°，纬度差为1°的子午线弧长X=110863m，两点间线长为110862.869m通过比拟可知，同一平行圈或同一子午线两点间线长度与对应的平行圈弧长或子午线弧长相等。六：程序代码#include<stdio.h>#include<math.h>doublehudu(double,double,double);/*度分秒转换为弧度*/doubledu(double);/*弧度转换为度*/doublefen(double);/*弧度转换为分*/doublemiao(double);/*弧度转换为秒*/#definePI3.1415926voidmain(void){intk; printf("请选择主题算法，假设执行正算，请输入1；假设执行反算，请输入2。\n"); scanf("%d",&k);/*主题正算*/ if(k==1) { doubleax,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; doublee2,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.2966; 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); elseif(sin(A1)>0&&tan(q)<0) q=PI-fabs(q); elseif(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); elseif(sin(A1)<0&&tan(A2)<0) A2=PI-fabs(A2); elseif(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); } /*主题反算*/ elseif(k==2) {doublefx,fy,fz,gx,gy,gz,hx,hy,hz,ix,iy,iz,jz,kz,B1,B2,L1,L2,S,A1,A2; intjx,jy,kx,ky; doublee2,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.2966; 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); elseif(p>0&&q<0) A1=PI-fabs(A1); elseif(p<0&&q<0) A1=PI+fabs(A1); elseA1=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); elseif(p<0&&q>0) A2=PI-fabs(A2); elseif(p>0&&q>0) A2=PI+fabs(A2); elseA2=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");}/*度分秒转换为弧度*/doublehudu(doublex,doubley,doublez){ doubleA0; A0=(x+y/60+z/3600)*PI/180;returnA0;}/*弧度转换为度*/++++++++++++++doubledu(doubleB0){ doublex0; x0=(int)(B0*180/PI); returnx0;}/*弧度转换为分*/doublefen(doubleC0){ double_y,y0;_y=(int)(C0*180/PI); y0=(fabs)((int)((C0*180/PI-_y)*60)); returny0;}/*弧度转换为秒*/doublemiao(doubleD0){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)); returnz0;}七：演算例如正算反算正算计算结果30°30′00.00″-0.505246786114°20′00.00″-0.610640770225°00′00.00″-0.707105781S10000000.000.3565378880.5308581490.003350562-0.6097894721.766987874E-060.7925632244.622109698E-10-0.692911606-1.1037772801.5714413545.270424445E-031.055977930-37°43′44.1351″1.39384099751°16′32.4976″1.57247898850°21′22.4896″-0.609789747反算计算结果0.530858149-0.656869486-0.3091512840.682919788-0.5266000700.400921954计算值趋近次数1234-1.100563429-1.103768055-1.103777253-1.103777280p-0.705956271-0.707102492-0.707105771-0.707105781q-0.708255369-0.707109082-0.707105789-0.707105780224°54′24.6737″224°59′59.0390″224°59′59.9973″225°00′00.0002″0.9999999960.9999985940.9999985840.9999985842.755182939E-04-1.677035942E-04-1.682644780E-04-1.682660878E-041.5705208091.5724733641.5724789731.5724789890.6087976010.6097869050.6097897390.609789747-0.692114035-0.692909315-0.692911599-0.6929116063.350558481E-033.350561868E-033.350561878E-033.350561878E-031.770381674E-61.766997603E-61.766987901E-61.766987874E-64.639919221E-104.622160703E-104.622109841E-104.622109695E-10**5.2638615645.2704055945.2704243925.270424446-1.103768055-1.103777253-1.103777280-1.1037772801.571441354-0.6097897471.055977929-0.5052467861.39384099650°21′22.4897″Y-4.367643449E-10S10000000.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].科学技术与工程,20