大地主题解算方法综述

白塞尔大地主题解算的基本思想
首先，白塞尔大地主题解算的基本思想之一是建立地壳变形的力学模型。

地壳变形是地球表面的一项重要现象，是由于地质作用和地球内部力
学过程的结果。

地壳变形的力学模型是研究地壳变形的重要方法。

常见的
地壳变形的力学模型有弹性模型、弹塑性模型和拟弹性模型等。

这些模型
可以描述地壳的变形过程，通过对地壳的变形过程进行建模，可以更好地
理解地壳变形的机制和动力学过程。

其次，白塞尔大地主题解算的思想之一是研究现今地球表面的动力学
过程。

地球表面的动力学过程包括板块运动、地震活动、火山喷发等。

这
些过程在地球的长期演化中起着重要的作用。

通过对现今地球表面动力学
过程的研究，可以揭示地球内部的结构和动力学机制，进而更好地理解地
壳变形的成因和发展变化的规律。

第三，白塞尔大地主题解算的基本思想之一是研究地壳变形的成因。

地壳变形的成因是地壳运动的基本原因，也是地球科学研究的一个重要问题。

地壳变形的成因包括构造运动、地壳应力状态的改变、地震活动等。

通过研究地壳变形的成因，可以更好地了解地壳变形的机制和规律，进而
为地震预测和地壳运动的控制提供科学依据。

总之，白塞尔大地主题解算的基本思想是通过建立地壳变形的力学模型，研究现今地球表面的动力学过程，探究地壳变形的成因，揭示地壳变
形与地球动力学过程之间的相互关系，以推动地壳运动的理论和应用研究。

这一思想对于研究地壳运动的机制和规律，了解地壳变形的成因和动力学
过程，提高地震预测和地壳运动控制的水平具有重要意义。

大地主题的数值解法范业明1,王解先1,2,刘慧芹1(1 同济大学测量与国土信息工程系,上海 200092;2 现代工程测量国家测绘局重点实验室,上海 200092)摘要:本文基于椭球面上大地线的微分方程,将法截弧方位角与大地线方位角之间的关系作为初始条件,通过用数值方法求解大地线的微分方程,进行大地主题的正反解。

并以实际数据验证了其正确性与可行性。

本法均采用封闭公式计算,精度高,公式简单,特别适用于计算机解算。

关键词:大地线;微分方程;数值解法中图分类号:P226+1文献标识码:BAbstract :Based on differential equation of geodesic on the surface of ellipsoid,the problem of direct and inverse solution of geodetic is solved through numerical method.The relationship of normal arc and geodesic is deduced and introduced as the initial condition for the differential equation.Numerical experiments are given,and the validity and feasibility of the method proposed in this paper have been proved.Besides,all the formulas are close,simple and easy to be realized by computer.Key words :geodesic;differential equation;numerical method 收稿日期:2006 02 15;修订日期:2006 05 10作者简介:范业明(1982-),男(汉族),辽宁沈阳人,硕士研究生.0 引言一直以来由于计算工具的限制,大地主题解算一般都采用级数展开形式,如短距离的高斯平均引数法,长距离的贝塞尔-赫尔默特方法等等。

大地测量实验报告实验名称：白塞尔大地主题解算（正算和反算）实验目的：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、最小二乘法：最小二乘法是大地测量主题正算问题的另一种数值方法，它能够根据观测数据求出最优解。

它不仅可以计算出两个观测点之间的差距，而且还能够计算出观测点到大地磁北极点的差距。

白塞尔大地主题解算方向：学号：姓名：一．基本思路：基本思想：将椭球面上的大地元素按照白塞尔投影条件投影到辅助球面上，继而在球面上进行大地主题解算，最后在将球面上的计算结果换算到椭球面上。

其关键问题是找出椭球面上的大地元素与球面上相应元素之间的关系式，同时解决在球面上进行大地主题解算的方法。

正算流程：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.用逐次趋近法同时计算起点大地方位角、球面长度及经差，第一次趋近时，取δ＝０。

计算下式,重复上述计算过程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);}}四．程序执行结果：。