大地主题解算深度干货超精

合集下载

白赛尔大地主题解算

dS

M
cos A c
dL

dS

sin A N cos B

V c
sec B sin B
dA

dS

tan B sin N
AV c
tan B sin
A
二阶导数：
d2B dS 2

B
( dB dS
)
dB dS

A
( dB dS
)
dA dS

V4 c2
t (3 2
10
Fundation of Geodesy
4.7.3 高斯平均引数正算公式
高斯平均引数正算公式推导的基本思想：
首先把勒让德级数在 P１点展开改在大地线长度中点M展开，以使级数公式项数减少，收敛快，精度高；其次，考虑到求定中点 M 的复杂性，将 M 点用大地线两端点平均纬度及平均方位角相对应的 m 点来代替，并借助迭代计算便可顺利地实现大地主题正解。

f (BM , AM )
F (Bm
BM
Bm
,
Am

AM

Am )
+
dB ( dS )M

f f (Bm , Am ) ( Bm )(BM
f Bm ) ( Am )( AM
Am )
+
dB
dB
dB ( dS )M

( )
f (Bm , Am ) (
dS Bm

S sin Am 24 N m2
[S 2tm2
sin 2
Am

S 2 cos2 Am (1m2 9m2tm2 )]

用C语言实现大地主题解算

６
４
６
４
（争寺＋
‘ ）ｃｏ５４Ａ ’
裴连磊ＰＥＩＬｉａｎ — ｌｅｉ
（新疆地矿局测绘大队，乌鲁木齐８３００１７）（ＸｉｎｊｉａｎｇＧｅｏｌｏｇｙａｎｄＭｉｎｅｒａｌＢｕｒｅａｕ，ＳｕｒｖｅｙｉｎｇａｎｄＭａｐｐｉｎｇＢｒｉｇａｄｅ，Ｕｒｕｍｑｉ８３００１７，Ｃｈｉｎａ）
１－５计算经度差改正数
本文的程序实现是采用的白塞尔大地主题解算的方法，根据白塞尔大地主题解算的方法，分析得出适合计算机编程的大地主题正反算的具体实现步骤。１白塞尔法大地主题正算思想已知量：大地线起点的纬度Ｂ，经度Ｌ，大地方位角Ａ１及大地线长度Ｓ。求解量：大地线终点的纬度Ｂ：，经度Ｌ２
ｃ番一番＋．．
ｃ＿ｂ（薷４一６一）
２
＝
ｃ０ｓｕ２＝
Ｌａ１＝ｓｉｎｕｌｓｉｎｕ２Ｉａ２＿ｃｏｓｕｌｃ０ｓｕ２ｂ＝
ＣＯＳＵｌｓｉｎｕ２，ｂ２＝ｓｉｎｕｌＣＯＳＵ２
之值。即下式：
求解量：大地线长度Ｓ及起、终点处的大地方位角Ａ
及Ａ＇。
Ａ：ｂ（１一鲁十击ｋ６＿－
２．１辅助计算ｗ：

大地主题解算深度干货超精

大地主题解算深度干货超精在我们生活的世界上，大地是我们所处的基本环境。

它是被压迫、被掠夺的，也是我们努力保护和改善的对象。

深度解算是一种精确的测量和计算方法，可以帮助我们更好地了解和利用大地资源。

本文将介绍大地主题解算和其带来的深度干货。

一、什么是大地主题解算？大地主题解算是一种基于卫星观测数据和大地测量学原理的技术，用于精确计算大地物体的三维位置和形态变化。

通过对不同时刻的卫星图像进行分析和计算，可以得到地面物体的精确坐标、高程等信息。

大地主题解算的核心思想是通过计算大地物体的三维位置和形态变化，来揭示地球表面的变化过程。

这种技术可以用于监测地表的沉降、地壳的变形、建筑物的结构和变化等。

它在测绘、地质灾害监测、城市规划等领域具有广泛的应用前景。

二、大地主题解算的应用领域1. 地质灾害监测地质灾害是世界各地普遍存在的问题，如地震、山体滑坡、地裂缝等。

通过大地主题解算技术，可以实时监测和预测地质灾害的发生和演变趋势。

这对于提前采取相应的防灾措施、保护人民的生命财产安全具有重要意义。

2. 城市规划随着城市化进程的不断加快，城市规划也变得愈发重要。

大地主题解算可以为城市规划提供精确的地面信息，包括土地利用、道路交通、建筑物布局等。

这有助于提高城市规划的科学性和有效性，减少资源浪费和环境污染。

3. 地表沉降监测地表沉降是地下水开采、矿井开采等人类活动造成的一个重要问题。

通过大地主题解算技术，可以实时监测地表的沉降情况，并对地下水开采和矿井开采等活动进行合理调整和管理。

这有助于减少地表沉降对城市建设和生态环境的不利影响。

4. 环境保护大地主题解算可以为环境保护提供准确的数据支持。

例如，通过对森林面积、湿地面积等进行监测和计算，可以及时发现和预防环境变化和破坏。

这对于保护生态环境和维护生态平衡具有重要意义。

三、大地主题解算的优势和挑战大地主题解算作为一种先进的测量和计算技术，具有很多优势。

首先，它可以实时获取和分析大地物体的位置和形态变化，具有高精度和高时效性。

分段累加法大地主题解算及高斯投影

大地测量编程实习报告班号：XXXX 学号：XXXXXXXXXX 姓名：XXX大地主题解算（分段累加法正算）结果截图：主要代码：privatevoid button1_Click(object sender, EventArgs e){double bx, by, bz, B1, lx, ly, lz, L1, ax, ay, az, A1, S;double dB, dL;double e2 = 0.006693421622966, a = 6378245.0000000000;double M, N, W, C;double B2 = 0, L2 = 0, A2 = 0;int Bx, By, Lx, Ly, Ax, Ay;double Bz, Lz, Az;bx = Convert.ToDouble(du1.Text); by = Convert.ToDouble(fen1.Text); bz = Convert.ToDouble(miao1.Text);lx = Convert.ToDouble(du2.Text); ly = Convert.ToDouble(fen2.Text); lz = Convert.ToDouble(miao2.Text); ax = Convert.ToDouble(du3.Text); ay = Convert.ToDouble(fen3.Text); az = Convert.ToDouble(miao3.Text); S = Convert.ToDouble(m.Text);double PI = Math.PI;B1 = (bx + by / 60 + bz / 3600) * PI / 180;L1 = (lx + ly / 60 + lz / 3600) * PI / 180;A1 = (ax + ay / 60 + az / 3600) * PI / 180;W = Math.Sqrt(1 - e2 * (Math.Sin(B1) * Math.Sin(B1)));M = a * (1 - e2) / (W * W * W);N = a / W;C = N * Math.Cos(B1) * Math.Sin(A1);int n ,i; double s;n = (int)(S / 4000 + 1);s = S / n;for (i = 1; i <= n; i++){dB = Math.Cos(A1) * s / M;dL = Math.Sin(A1) * s / N / Math.Cos(B1);B2 = B1 + dB;L2 = L1 + dL;W = Math.Sqrt(1 - e2 * (Math.Sin(B2) * Math.Sin(B2)));M = a * (1 - e2) / (W * W * W);N = a / W;A2 = Math.Asin(C / N / Math.Cos(B2));B1 = B2; L1 = L2; A1 = A2;}Bx =(int)( B2 / PI * 180);By = (int)((B2 / PI * 180 - Bx) * 60);Bz = Math.Abs(((B2 / PI * 180 - Bx) * 60 - By) * 60);By = Math.Abs(By);Lx = (int)(L2 / PI * 180);Ly = (int)((L2 / PI * 180 - Lx) * 60);Lz = Math.Abs(((L2 / PI * 180 - Lx) * 60 - Ly) * 60);Ly = Math.Abs(Ly);Ax = (int)(A2 / PI * 180);Ay = (int)((A2 / PI * 180 - Ax) * 60);Az = Math.Abs(((A2 / PI * 180 - Ax) * 60 - Ay) * 60);Ay = Math.Abs(Ay);du4.Text = Bx.ToString(); fen4.Text = By.ToString(); miao4.Text = Bz.ToString("0.0000"); du5.Text = Lx.ToString(); fen5.Text = Ly.ToString(); miao5.Text = Lz.ToString("0.0000"); du6.Text = Ax.ToString(); fen6.Text = Ay.ToString(); miao6.Text = Az.ToString("0.0000"); }高斯投影正、反算结果截图：正算：反算：主要代码：正算：privatevoid button1_Click(object sender, EventArgs e){double a = 0 , e12 = 0, e2 = 0;double PI = Math.PI;string tt;try{tt = boBox1.Items[boBox1.SelectedIndex].ToString();}catch{MessageBox.Show("请选择坐标系！"); return;}if (pareTo("1954年北京坐标系") == 0){a = 6378245;e2 = 0.006693421622966;e12 = 0.006738525414683;}if (pareTo("1980年国家大地坐标系") == 0){a = 6378140;e2 = 0.006694384999588;e12 = 0.006739501819473;}if (pareTo("WGS-84") == 0){a = 6378137;e2 = 0.00669437999013;e12 = 0.00673949674227;}if (pareTo("CGCS2000") == 0){a = 6378137;e2 = 0.00669438002290;e12 = 0.00673949677548;}double m0, m2, m4, m6, m8, a0, a2, a4, a6, a8, N, B, L, X, Bb, cosB ,sinB , ρ, η2, t, l, d, x, y;B = Convert.ToDouble(du1.Text) + (Convert.ToDouble(fen1.Text) / 60) + (Convert.ToDouble(miao1.Text) / 3600); L = Convert.ToDouble(du2.Text) + (Convert.ToDouble(fen2.Text) / 60) +(Convert.ToDouble(miao2.Text) / 3600);Bb = B * PI / 180;m0 = a * (1 - e2);m2 = 3 * e2 * m0 / 2;m4 = 5 * e2 * m2 / 4;m6 = 7 * e2 * m4 / 6;m8 = 9 * e2 * m6 / 8;a0 = m0 + (m2 / 2) + ((3 * m4) / 8) + ((5 * m6) / 16) + ((35 * m8) / 128);a2 = m2 / 2 + m4 / 2 + 15 * m6 / 32 + 7 * m8 / 16;a4 = m4 / 8 + 3 * m6 / 16 + 7 * m8 / 32;a6 = m6 / 32 + m8 / 16;a8 = m8 / 128;ρ = 180 * 3600 / PI;cosB = Math.Cos(Bb);sinB = Math.Sin(Bb);η2 = e12 * cosB * cosB;t = Math.Tan(Bb);d = Math.Floor(L / 6) + 1;l = Math.Abs(L - (6 * d - 3)) * 3600;N = a / Math.Sqrt(1 - e2 * sinB * sinB);X = a0 * Bb - sinB * cosB * ((a2 - a4 + a6) + (2 * a4 - 16 * a6 / 3) * sinB * sinB + 16 * a6 * sinB * sinB * sinB * sinB / 3);x = X + N * sinB * cosB * l * l / (2 * ρ * ρ) + N * sinB * cosB * cosB * cosB * (5 - t * t + 9 * η2) * l * l * l * l / (24 * ρ * ρ * ρ * ρ) + N * sinB * cosB * cosB * cosB * cosB * cosB * (61 - 58 * t * t + t * t * t * t) * l * l * l * l * l * l / (720 * ρ * ρ * ρ * ρ * ρ * ρ);y = N * cosB * l / ρ + N * cosB * cosB * cosB * (1 - t * t + η2) * l * l * l / (6 * ρ * ρ* ρ)+ N * cosB * cosB * cosB * cosB * cosB * (5 - 18 * t * t + t * t * t * t) * l * l * l * l * l/ (120 * ρ * ρ * ρ * ρ * ρ);x1.Text = x.ToString("0.0000")+"m";y1.Text = y.ToString("0.0000")+"m";}反算：privatevoid button2_Click(object sender, EventArgs e){double a = 0, e12 = 0, e2 = 0;double PI = Math.PI;string tt;try{tt = boBox2.Items[boBox2.SelectedIndex].ToString();}catch{MessageBox.Show("请选择坐标系！"); return;}if (pareTo("1954年北京坐标系") == 0){a = 6378245;e2 = 0.006693421622966;e12 = 0.006738525414683;}if (pareTo("1980年国家大地坐标系") == 0){a = 6378140;e2 = 0.006694384999588;e12 = 0.006739501819473;}if (pareTo("WGS-84") == 0){a = 6378137;e2 = 0.00669437999013;e12 = 0.00673949674227;}if (pareTo("CGCS2000") == 0){a = 6378137;e2 = 0.00669438002290;e12 = 0.00673949677548;}double x, y, m0, m2, m4, m6, m8, a0, a2, a4, a6, a8, Bf, Bf1, B , Mf , W, Nf, tf, ηf2, n2, n4, n6, n8 , l;x = Convert.ToDouble(x2.Text);y = Convert.ToDouble(y2.Text);m0 = a * (1 - e2);m2 = 3 * e2 * m0 / 2;m4 = 5 * e2 * m2 / 4;m6 = 7 * e2 * m4 / 6;m8 = 9 * e2 * m6 / 8;a0 = m0 + (m2 / 2) + ((3 * m4) / 8) + ((5 * m6) / 16) + ((35 * m8) / 128);a2 = m2 / 2 + m4 / 2 + 15 * m6 / 32 + 7 * m8 / 16;a4 = m4 / 8 + 3 * m6 / 16 + 7 * m8 / 32;a6 = m6 / 32 + m8 / 16;a8 = m8 / 128;Bf1 = 0;Bf = x / a0;while (Math.Abs(Bf - Bf1) > 1E-10){Bf1 = Bf;double sinB = Math.Sin(Bf1);double cosB = Math.Cos(Bf1);double sin2B = sinB * cosB * 2;double sin4B = sin2B * (1 - 2 * sinB * sinB) * 2;double sin6B = sin2B * Math.Sqrt(1 - sin4B * sin4B) + sin4B * Math.Sqrt(1 - sin2B * sin2B);Bf = (x - (-a2 / 2.0 * sin2B + a4 / 4.0 * sin4B - a6 / 6.0 * sin6B)) / a0;}double sinBf = Math.Sin(Bf);double cosBf = Math.Cos(Bf);ηf2 = e12 * cosBf * cosBf;tf = Math.Tan(Bf);W = Math.Sqrt(1 - e2 * sinBf * sinBf);Mf = a * (1 - e2) / (W * W * W);Nf = a / W;B = Bf - tf * y * y / (2 * Mf * Nf) + tf * (5 + 3 * tf * tf + ηf2 - 9 * ηf2 * tf * tf) * y * y * y * y / (24 * Mf * Nf * Nf * Nf) - tf * (61 + 90 * tf * tf + 45 * tf * tf * tf * tf) * y * y * y * y * y * y / (720 * Mf * Nf * Nf * Nf * Nf * Nf);l = y / (Nf * cosBf) - (1 + 2 * tf * tf + ηf2) * y * y * y / (6 * Nf * Nf * Nf * cosBf) + (5 + 28 * tf * tf + 24 * tf * tf * tf * tf + 6 * ηf2 + 8 * ηf2 * tf * tf) * y * y * y * y * y / (120 * Nf * Nf * Nf * Nf * Nf * cosBf);int Bx, By, lx, ly;double Bz, lz;Bx = (int)(B / PI * 180);By = (int)((B / PI * 180 - Bx) * 60);Bz = Math.Abs(((B / PI * 180 - Bx) * 60 - By) * 60);By = Math.Abs(By);lx = (int)(l / PI * 180);ly = (int)((l / PI * 180 - lx) * 60);lz = Math.Abs(((l / PI * 180 - lx) * 60 - ly) * 60);ly = Math.Abs(ly);du3.Text = Bx.ToString();fen3.Text = By.ToString();miao3.Text = Bz.ToString("0.0000");du4.Text = lx.ToString(); fen4.Text = ly.ToString();miao4.Text = lz.ToString("0.0000");}。

大地主题解算

大地主题解算一、实验目的：1.提高运用计算机语言编程开发的能力；2.加深对大地主题解算计算公式及辅助参数的理解并掌握计算步骤；3.通过编程语言实现大地主题解算。

二、工具：Windows XP Mode 环境下的Microsoft Visual C++ 6.0三、注意事项：1.计算所需变量多，容易混淆；2.正反算函数的编写；3.函数调用；4.弧度与角度之间的转化。

四、实验要求：1.提交报告，实验总结，编写代码；2.独立编程，调试运行；3.上交成果：编写思想，编写过程，问题分析，源代码，计算结果；五、编程过程实现：1.对白塞尔法大地主题解算有一定的了解，并参考教材P148-P150；2.由于参数较多，而在C语言环境下很多符号无法定义，需要符合要求的定义符号替代书本上那些无法直接在C语言环境下定义的符号来达到实现实验的目的；3.程序中采用弧度与度分秒之间转换的函数定义与调用，减轻一定的实验麻烦；4.在C语言环境下，数学函数fabs代替abs起绝对值作用，atan代替arctan起反函数作用；5.程序中尤其注意弧度与角度之间转换，在C语言环境下电脑默认为弧度。

六、源程序代码：#include<stdio.h>#include<math.h>double hudu(double,double,double); /*度分秒转换为弧度*/double du(double); /*弧度转换为度*/double fen(double); /*弧度转换为分*/double miao(double); /*弧度转换为秒*/#define PI 3.1415926void main (void){int k;printf("请选择执行正算或者反算，若执行正算，请输入1；若执行反算，请输入2。

\n");scanf("%d",&k);/*正算*/if(k==1){double bz,lz,az,S,bz2,lz2,az2,B1,L1,A1,B2,L2,A2,bx,by,lx,ly,ax,ay;int bx2,by2,lx2,ly2,ax2,ay2;doublee2,W1,sinu1,cosu1,sinA0,coto1,sin2o1,cos2o1,sin2o,cos2o,A,B,C,r,t,o0,o,g,si nu2,q;/*以度分秒顺序输入数据*/printf("请输入大地线起点纬度度分秒\n");scanf("%lf%lf%lf",&bx,&by,&bz);printf("请输入大地线起点经度度分秒\n");scanf("%lf%lf%lf",&lx,&ly,&lz);printf("请输入大地方位角度分秒\n");scanf("%lf%lf%lf",&ax,&ay,&az);printf("请输入大地线长度\n");scanf("%lf",&S);/*调用函数*/B1=hudu(bx,by,bz);L1=hudu(lx,ly,lz);A1=hudu(ax,ay,az);/*白塞尔大地主题解算*/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);elseq=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);elseA2=2*PI-fabs(A2);/*调用函数*/bx2=(int)(du(B2));by2=(int)(fen(B2));bz2=miao(B2);lx2=(int)(du(L2));ly2=(int)(fen(L2));lz2=miao(L2);ax2=(int)(du(A2));ay2=(int)(fen(A2));az2=miao(A2);printf("大地线终点纬度度分秒分别为：\n%d\n%d\n%lf\n",bx2,by2,bz2);printf("大地线终点经度度分秒分别为：\n%d\n%d\n%lf\n",lx2,ly2,lz2);printf("终点大地方位角度分秒分别为：\n%d\n%d\n%lf\n",ax2,ay2,az2);}/*反算*/else if(k==2){doublebz,lz,bz2,lz2,az,az2,B1,L1,B2,L2,S,A1,A2,bx,by,lx,ly,bx2,by2,lx2,ly2;int ax,ay,ax2,ay2;doublee2,W1,W2,sinu1,sinu2,cosu1,cosu2,L,a1,a2,b1,b2,g,g2,g0,r,p,q,sino,coso,o,si nA0,x,t1,t2,A,B,C,y;/*以度分秒顺序输入数据*/printf("请输入大地线起点纬度度分秒\n");scanf("%lf%lf%lf",&bx,&by,&bz);printf("请输入大地线起点经度度分秒\n"); scanf("%lf%lf%lf",&lx,&ly,&lz);printf("请输入大地线终点纬度度分秒\n"); scanf("%lf%lf%lf",&bx2,&by2,&bz2); printf("请输入大地线终点经度度分秒\n"); scanf("%lf%lf%lf",&lx2,&ly2,&lz2);/*调用函数*/B1=hudu(bx,by,bz);L1=hudu(lx,ly,lz);B2=hudu(bx2,by2,bz2);L2=hudu(lx2,ly2,lz2);/*白塞尔大地主题解算*/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*/g0=-662.904266/3600*PI/180; 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);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);elseo=PI-fabs(o);sinA0=cosu1*sin(A1);x=2*a1-(1-sinA0*sinA0)*coso;t1=(33523299-(28189-70*(1-sinA0*sinA0))*(1-sinA0*sinA0))*0.0000000001;t2=(28189-94*(1-sinA0*sinA0))*0.0000000001;g2=(t1*o-t2*x*sino)*sinA0;/*检验循环次数*/printf("\ng2=%lf\ng0=%lf\n",g2,g0);if(g2<=g0)break;elser=L+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);elseA2=2*PI-fabs(A2);/*调用函数*/ax=(int)(du(A1));ay=(int)(fen(A1));az=miao(A1);ax2=(int)(du(A2));ay2=(int)(fen(A2));az2=miao(A2);printf("起点大地方位角度分秒分别为：\n%d\n%d\n%lf\n",ax,ay,az);printf("终点大地方位角度分秒分别为：\n%d\n%d\n%lf\n",ax2,ay2,az2);printf("大地线长度为：%lf\n",S);}/*数据错误*/elseprintf("数据错误，请重新输入\n");}/*度分秒转换为弧度*/double hudu(double a0,double b0,double c0){double A0;A0=(a0+b0/60+c0/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;}大地主题解算正算：大地主题解算反算：以上检验数据来自书本P151，P152。

大地主题解算方法综述

第 32卷第 4期 2007年 7月
测绘科学 Sc ience o f Survey ing and M app ing
V o l132 N o14 Ju l1
大地主题解算方法综述
周振宇, 郭广礼, 贾新果
(中国矿业大学环测学院, 江苏徐州 221008)
【摘要】大地主题解算是大地测量中的重要问题, 由于椭球的复杂性, 随之产生的解算方法也是多种多样。本
和方位角差。
像这类方法还有很多, 就不再一一赘述, 其基本原理
大致相同。
21 2 以白塞尔大地投影为基础
可以近似认为地球椭球的形状与圆球区别不大。在椭
球面上解算大地主题问题时可借助于球面三角学公式简单
而严密的进行。因此, 如将椭球面上的大地线长度投影到
在于: 解算精度与距离有关, 距离越长, 收敛越慢, 因此
只适用于较短的距离。当公式为四次项时 , 在中纬度地区,
公式可用于 200km 以下的大地主题解算, 经纬度计算精度
可达 01 0001", 方位角计算精度达 01001", 计算边长可精
确至厘米。是短距离解算的理想公式。当距离小于 70km
椭球的过渡。
白塞尔 ( Besse l) 首先提出并解决了投影条件, 使这一解
法得以实现。
这类公式的特点是计算公式展开 e2 或的幂级数, 解算
精度与距离长短无关。因此它既适用于短距离解算, 也适用于长距离解算 [ 5] 。其主要缺点在于: 由 S 求 R、由 L 求 K, 或相反的运算, 需要进行迭代。同时还要预先计算辅助

大地主题解算-C#

大地主题解算-正算-C#大地主题解算-正算-程序using System;using System.Collections.Generic;using ponentModel;using System.Data;using System.Drawing;using System.Linq;using System.Text;using System.Windows.Forms;namespace WindowsFormsApplication2{public partial class Form1 : Form{public Form1(){InitializeComponent();}private void button1_Click(object sender, EventArgs e) {int n,du,fen,miao;double B1, L1, A1, S, dS, dB1, dL1, M1, N1, W1;double E = 0.006694384999588, a = 6378140, C;B1 =( Convert.ToDouble(textBox1.Text) + Convert.ToDouble(textBox9.Text)/60 + Convert.ToDouble(textBox10.Text)/3600 )* Math.PI / 180;L1 = (Convert.ToDouble(textBox2.Text) + Convert.ToDouble(textBox11.Text)/60 + Convert.ToDouble(textBox12.Text)/3600) * Math.PI / 180; ;A1 = (Convert.ToDouble(textBox3.Text) + Convert.ToDouble(textBox13.Text)/60 + Convert.ToDouble(textBox14.Text)/3600) * Math.PI / 180; ;S = Convert.ToDouble(textBox5.Text);n = Convert.ToInt16(textBox4.Text);double[] L = new double[100000];double[] B = new double[100000];double[] dL = new double[100000];double[] dB = new double[100000];double[] A = new double[100000];double[] W = new double[100000];double[] N = new double[100000];double[] M = new double[100000];W1 = Math.Pow(1 - E * Math.Pow(Math.Sin(B1), 2), 0.5);N1 = a / W1;M1 = a * (1 - E) / (W1 * W1 * W1);C = N1 * Math.Sin(A1) * Math.Cos(B1);dS = S / n;dB1 = dS * Math.Cos(A1) / M1;dL1 = dS * Math.Sin(A1) / N1 / Math.Cos(B1);B[0] = B1;L[0] = L1;dB[0] = dB1;dL[0] = dL1;A[0] = A1;M[0] = M1;N[0] = N1;W[0] = W1;int i;for (i = 1; i < n; i++){dB[i - 1] = dS * Math.Cos(A[i - 1]) / M[i - 1];B[i] = B[i - 1] + dB[i - 1];W[i] = Math.Pow(1 - E * Math.Pow(Math.Sin(B[i]), 2), 0.5);N[i] = a / W[i];M[i] = a * (1 - E) / (W[i] * W[i] * W[i]);dL[i - 1] = dS * Math.Sin(A[i - 1]) / (N[i - 1] * Math.Cos(B[i - 1])); L[i] = L[i - 1] + dL[i - 1];A[i] = Math.Asin(C / (Math.Cos(B[i]) * N[i]));}du = Convert.ToInt16(Math.Floor(B[i - 1] * 180 / Math.PI));//dufen = Convert.ToInt16(Math.Floor((B[i - 1] * 180 / Math.PI - Math.Floor(B[i - 1] * 180 / Math.PI)) * 60)); //fenmiao = Convert.ToInt16(Convert.ToInt16(((B[i - 1] * 180 / Math.PI - Math.Floor(B[i - 1] * 180 / Math.PI)) * 60 - Math.Floor((B[i - 1] * 180 / Math.PI - Math.Floor(B[i - 1] * 180 / Math.PI)) * 60)) * 60));//秒?textBox6.Text = Convert.ToString(du) + "度è" + Convert.ToString(fen) + "分?" + Convert.ToString(miao) + "秒?";du = Convert.ToInt16(Math.Floor(L[i - 1] * 180 / Math.PI));//dufen = Convert.ToInt16(Math.Floor((L[i - 1] * 180 / Math.PI - Math.Floor(L[i - 1] * 180 / Math.PI)) * 60)); //fenmiao = Convert.ToInt16(Convert.ToInt16(((L[i - 1] * 180 / Math.PI - Math.Floor(L[i - 1] * 180 / Math.PI)) * 60 - Math.Floor((L[i - 1] * 180 / Math.PI - Math.Floor(L[i - 1] * 180 / Math.PI)) * 60)) * 60));//秒?textBox7.Text = Convert.ToString(du) + "度è" + Convert.ToString(fen) + "分?" + Convert.ToString(miao) + "秒?";du = Convert.ToInt16(Math.Floor((180 + A[i - 1] * 180 / Math.PI)));//dufen = Convert.ToInt16(Math.Floor(((180 + A[i - 1] * 180 / Math.PI) - Math.Floor((180 + A[i - 1] * 180 / Math.PI))) * 60)); //fenmiao = Convert.ToInt16(Convert.ToInt16((((180 + A[i - 1] * 180 / Math.PI) - Math.Floor((180 + A[i - 1] * 180 / Math.PI))) * 60 - Math.Floor(((180 + A[i - 1] * 180 / Math.PI) - Math.Floor((180 + A[i - 1] * 180 / Math.PI))) * 60)) * 60));//秒?textBox8.Text = Convert.ToString(du) + "度è" + Convert.ToString(fen) + "分?" + Convert.ToString(miao) + "秒?";}private void Form1_Load(object sender, EventArgs e){} }}。

白塞尔大地主题解算(正算和反算)

大地测量实验报告实验名称：白塞尔大地主题解算（正算和反算）实验目的：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为大地线终点纬度、经度及方位角。

椭球面上大地问题的解算概要

高斯平均引数公式
若取大地线中点展开，得：
2 2 3 3 dB S d B S d B S B2 BM 2 3 dS M 2 dS M 8 dS M 48
N
2 2 3 3 dB S d B S d B S B1 BM 2 3 dS M 2 dS M 8 dS M 48
2 dB cos A Vm cos Am 由大地线的微分公式： dS M Nm
由大地线的微分公式，得其一阶导数为：
dB cos A dS M dL sin A dS N cos B dA tan B sin A dS N
二阶和三阶导数采用复合函数求导法计算：
d 2 B dB dB dB dA 2 dS B dS dS A dS dS d 3 B d 2 B dB d 2 B dA 3 2 2 dS B dS dS A dS dS
同理可求出四阶以上的导数和L、A的高阶导数，代入展开式即可。
三、高斯平均引数正解公式
（一）基本思想
首先把勒让德级数在P1点展开改为在大地线长度中点M展开，以使级数公式项数减少、收敛快、精度高；其次，考虑到求定中点M的复杂性，将M点用大地线两端的平均纬度及平均方位角相对应的 m点来代替，并借助迭代计算，便可顺利的实现大地问题的正解。
d2A S2 Am AM dS 2 8 M
1 Am A1 A2 2
1 Lm L1 L2 , 2
展开成级数，得： 2
dB dB dB dB Bm BM Am AM A dS m dS M dS m B dS m

大地测量学第四章-5大地测量主题解算

2R
2
ab sin C
A A0 , 3
B B0 , 3
bc sin A ac sin B ab sin C 2 2 2
C C0 3
化算平面角需要球面角超，而球面角超的计算又需要平面角，因此直接用球面角计算球面角超就带有误差。当边长不大于90km时，这种误差小于0.0005″，故可直接用球面角代替平面角计算球面角超ε
dB d 3 B B2 B1 B ( ) M S ( 3 ) M S 3 5次项 dS 24 dS
三、高斯平均引数正算公式
(1)建立级数展开式: 同理可得:
dL d 3 L L2 L1 L ( ) M S ( 3 ) M S 3 5次项 dS 24 dS
Bm BM ,
Am AM
三、高斯平均引数正算公式
(3)求以Bm、Am为依据的导数: 经整理得：
2 Vm S2 2 2 B S cos Am {1 [sin 2 Am ( 2 3t m 2 m ) 2 Nm 24 N m 2 2 2 2 3 m cos 2 Am ( 1 m 9 m t m )]} 5次
第四章 Ⅴ大地测量主题解算
——大地主题解算思路 ——勒让德级数式 ——高斯平均引数正算公式 ——高斯平均引数反算公式
上一讲应掌握的内容 1、垂线偏差改正垂线偏差对水平方向的影响 "u ( "sin A1 "cos A1) tan 1 2、标高差改正 e2 由照准点高度而引起的改正 "h H 2 cos2 B2 sin 2 A1
以大地线在大地坐标系中的微分方程为基础主要特点：解算精度与距离有关，距离越长，收敛越慢，因此只适用于较短的距离。典型解法：高斯平均引数法

白塞尔大地主题解算

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

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

正算流程：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);}}四．程序执行结果：。

大地测量主题解算

7、计算(B2)3、(A21)3
一般情况下主项趋近3次，改正项趋近2次就可满足要求。
四、高斯平均引数反算公式
高斯平均引数反算公式可以依正算公式导出:
S
sin
Am

L" "
Nm
cos Bm

S sin Am
24
N
2 m
[S 2tm2
sin 2
Am

S 2 cos2 Am (1 m2 9m2tm2 )]
t
2 m

2m2 )]}
5次
B2 B1 B, L2 L1 L, A21 A12 A 180
Am

1 2
( A12

A21 )

A12

1 2
A
Bm

1 2
(B2

B1)

B1

1 2
B
注意：
从公式可知，欲求ΔＬ，ΔＢ及ΔＡ，必先有Bm及 Am。但由于Ｂ2和Ａ21未知，故精确值尚不知，为
2f
2
上一讲应掌握的内容
四、整周数N值解算的一般原理
有可变频率法和固定频率法两种
五、全站仪中测距新技术
• 使用高频测距技术
· 温控与动态频率校正技术
• 无棱镜测距技术
· 目标自动识别技术
六、测距的误差分析和精度表达式
D c0 C 4fn
mD2
( c0 )2 m2
4 fn 2

Nm Vm2
高斯平均引数反算公式（续）
S sin Am r01L" r21B"2L" r03L"3 S cos Am s10B" s12B"L"2 s30B"3

大地测量学资源课-第4.7节

dB dS

A

dA dS
dA dS

V2 c2
sin
Acos
A(1
2t
2
2 )
d 2L dS 2

B

dL dS

dB dS

A

dL dS

dA dS

2V 2 c2
t
sec Bsin
Acos
A
三阶导数 d3B V5 cosA[sin2 A(13t2 2 9t22)32 cosA2(1t2 2 52t2)]
M
S3 8

B1 BM
dB dS M
S 2

1 2

d2B dS2
M
S2 4

1 6

d3B dS3
M
S3 8

上述两式相减可得：
(B2

B1 ) ''

B ''

''

dB dS
M
S

''
24

d 3B dS 3
B 1 B 0 , L1 L 0 , A12 A 0
B2
B1

B

dnB dS n
1
Sn n!
dB S dS 1

d 2B dS 2
1
S2 2!

d 3B dS 3
1
S3 3!

L2
L1

L

新型大地坐标系中的大地主题解算

（3）
s L1 ⎞ sin T1 ⎛ sin 2T1 ⎛ 2 4C 0 s L1 n1 ⎞ 2 ⎟ ⎜ ⎟ T2 = T1 − ⎜ ⎜ C 0 − R 2 ⎟ n S + 4n 2 ⎜ 2C 0 − R 2 + R 2 ⎟ S 1 0 ⎠ 1 0 0 ⎠ ⎝ ⎝ ⎡ 2C 3 sin 3T C sin 3T 3C sin 2T cos T1 ⎤ S 3 −⎢ 0 3 1 + 0 2 2 1 + 0 ⎥ n1 n1 R0 2 R02 n12 ⎣ ⎦ 6
（4）
2 大地坐标正反解的算法和公式
2.1 平均引数正算公式
以大地线中点 PM （即令 S P1PM = S PM P2 ）取代大地线起点 P 1 作为展开点，显然能提高数值计算的精度，相应的展开式则为
[3 ]
1 ⎛ d 3sL ⎞ 3 ⎛ ds L ⎞ ⎟ S ∆s L = s L2 − s L1 = ⎜ ⎟ S+ ⎜ 3 ⎟ 24 ⎜ dS ⎝ dS ⎠ M ⎝ ⎠M 1 ⎛ d 3sB ⎛ ds ⎞ ∆s B = s B2 − s B1 = ⎜ B ⎟ S + ⎜ 3 24 ⎜ ⎝ dS ⎠ M ⎝ dS ⎞ 3 ⎟ ⎟ S ⎠M
s Lm = 1 1 s L1 + s L2 及平均方位角 Tm = (T1 + T2 ) 作为平均引数来逼近 s LM 及 TM ，其间的 2 2
[3 ]
(
)
差异为二阶小量，舍去 5 次项后，即为
s LM − s Lm
S2 =− 8 S2 8
⎛ d 2 sL ⎜ ⎜ dS 2 ⎝ ⎛ d 2 sB ⎜ ⎜ dS 2 ⎝
⎞ S2 ⎟ ⎟ =− 8 ⎠M ⎞ S2 ⎟ = − ⎟ 8 ⎠M

大地主题解算(深度干货-超精)

大地主题解算(正算)代码与白塞尔大地主题解算大地主题解算(正算)代码：根据经纬度和方向角以及距离计算另外一点坐标新建模块->拷贝下面的大地主题(正算)代码，调用方法示例：起点经度：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的值。

大地主题解算方法综述

大地主题解算方法综述
大地主题解算方法综述
大地主题解算是大地测量中的重要问题,由于椭球的复杂性,随之产生的解算方法也是多种多样.本文分析了大地主题解算的一般方法及其公式的适用范围,指出了它们的特点和不足,讨论了现今常用解算方法存在的问题和进一步的研究方向.
作者：周振宇郭广礼贾新果 ZHOU Zhen-yu GUO Guang-li JIA Xin-guo 作者单位：中国矿业大学,环测学院,江苏徐州,221008 刊名：测绘科学ISTIC PKU英文刊名：SCIENCE OF SURVEYING AND MAPPING 年，卷(期)：2007 32(4) 分类号：P22 关键词：大地主题解算大地主题正算大地主题反算高斯平均引数法白塞尔公式。

分段累加法大地主题解算及高斯投影

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

大地主题解算(正算)代码与白塞尔大地主题解算大地主题解算(正算)代码：根据经纬度和方向角以及距离计算另外一点坐标新建模块->拷贝下面的大地主题(正算)代码，调用方法示例：起点经度：116.235(度)终点纬度：37.435(度)方向角：50（度）长度：500（米）终点经纬度（"经度,纬度"）=Computation(37.435,116.235,50,500)Const Pi = 3.14Private 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的值。

cos2A0 = 1 - sinA0 ^ 2e2 = Sqr(a ^ 2 - b ^ 2) / bk2 = e2 * e2 * cos2A0Dim aa, BB, CC, EE22, AAlpha, BBeta As Doubleaa = b * (1 + k2 / 4 - 3 * k2 * k2 / 64 + 5 * k2 * k2 * k2 / 256)BB = b * (k2 / 8 - k2 * k2 / 32 + 15 * k2 * k2 * k2 / 1024)CC = b * (k2 * k2 / 128 - 3 * k2 * k2 * k2 / 512)e2 = E1 * E1AAlpha = (e2 / 2 + e2 * e2 / 8 + e2 * e2 * e2 / 16) - (e2 * e2 / 16 + e2 * e2 * e2 / 16) * cos2A0 + (3 * e2 * e2 * e2 / 128) * cos2A0 * cos2A0BBeta = (e2 * e2 / 32 + e2 * e2 * e2 / 32) * cos2A0 - (e2 * e2 * e2 / 64) * cos2A0 * cos2A0'// 计算球面长度q0 = (S - (BB + CC * cos2q1) * sin2q1) / aasin2q1q0 = sin2q1 * Cos(2 * q0) + cos2q1 * Sin(2 * q0)cos2q1q0 = cos2q1 * Cos(2 * q0) - sin2q1 * Sin(2 * q0)q = q0 + (BB + 5 * CC * cos2q1q0) * sin2q1q0 / aa'// 计算经度差改正数theta = (AAlpha * q + beta * (sin2q1q0 - sin2q1)) * sinA0'// 计算终点大地坐标及大地方位角sinu2 = sinu1 * Cos(q) + cosu1 * Cos(A1) * Sin(q)B2 = Atn(sinu2 / (Sqr(1 - E1 * E1) * Sqr(1 - sinu2 * sinu2))) * 180 / Pilamuda = Atn(Sin(A1) * Sin(q) / (cosu1 * Cos(q) - sinu1 * Sin(q) * Cos(A1))) * 180 / PiIf (Sin(A1) > 0) ThenIf (Sin(A1) * Sin(q) / (cosu1 * Cos(q) - sinu1 * Sin(q) * Cos(A1)) > 0) Thenlamuda = Abs(lamuda)Elselamuda = 180 -Abs(lamuda)End IfElseIf (Sin(A1) * Sin(q) / (cosu1 * Cos(q) - sinu1 * Sin(q) * Cos(A1)) > 0) Thenlamuda = Abs(lamuda) - 180Elselamuda = -Abs(lamuda)End IfEnd IfL2 = L1 * 180 / Pi + lamuda - theta * 180 / PiA2 = Atn(cosu1 * Sin(A1) / (cosu1 * Cos(q) * Cos(A1) - sinu1 * Sin(q))) * 180 / PiIf (Sin(A1) > 0) ThenIf (cosu1 * Sin(A1) / (cosu1 * Cos(q) * Cos(A1) - sinu1 * Sin(q)) > 0) ThenA2 = 180 + Abs(A2)ElseA2 = 360 - Abs(A2)End IfElseIf (cosu1 * Sin(A1) / (cosu1 * Cos(q) * Cos(A1) - sinu1 * Sin(q)) > 0) ThenA2 = Abs(A2)ElseA2 = 180 - Abs(A2)End IfEnd IfComputation = L2 & "," & B2End FunctionPrivate Function rad(ByVal angle_d As Double) As Doublerad = angle_d * Pi / 180End Function白塞尔大地主题解算一、正算代码：#include<stdio.h>#include<math.h>#define ee 0.2966#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;double cota1,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=%.9f\n",B1);printf("L1=%.9f\n",L1);printf("A1=%.9f\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=%.9f\n",W1);printf("sinu1=%.9f\n",sinu1);printf("cosu1=%.9f\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=%.9f\n",sinA0);printf("cota1=%.9f\n",cota1);printf("sin2a1=%.9f\n",sin2a1);printf("cos2a1=%.9f\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=%.9f\n",cosA0A0);printf("A=%.3f\n",A);printf("B=%.6f\n",B);printf("C=%.9f\n",C);printf("d=%.7f\n",d);printf("e=%.9f\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=%.9f\n",a0);printf("m=%.9f\n",m);printf("n=%.9f\n",n);printf("a=%.9f\n",a);/*计算经度差改正数*/Q=(d*a+(e*(m-sin2a1))/3600/180*I)*sinA0;printf("Q=%.7f\n",Q);/*计算终点大地坐标及大地方位角*/sinu2=sinu1*cos(a)+cosu1*cos(A1)*sin(a);B2=180*atan(sinu2/((sqrt(1-ee))*(sqrt(1-sinu2*sinu2))))/I;R=180*atan(sin(A1)*sin(a)/(cosu1*cos(a)-sinu1*sin(a)*cos(A1)))/I;printf("sinu2=%.9f\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\n L2=%3f\n",A2,B2,L2);}double F(double a2,double b2,double c2){double d2;d2=(double)(a2+1.0*b2/60+1.0*c2/3600);d2=(d2/180)*I;return (d2);}运行结果：二、反算代码：#include<stdio.h>#include<math.h>#define ee 0.2966#define I 3.141592653double F(double,double,double);void main(void){double B1,B2,L1,L2,A1,A2,S,Y;double W1,W2,L,Q,R,A,B,C;double x,y,z,p,q;double x1,x2,x3,y1,y2,y3,z1,z2,z3,w1,w2,w3;double a1,a2,b1,b2,m,n;double sinp,cosp,sinu1,sinu2,cosu1,cosu2,sinA0,cosA0;Q=0;q=0;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("请输入起始数B2= ");scanf("%lf %lf %lf",&z1,&z2,&z3);B2=F(z1,z2,z3);printf("请输入起始数据L2= ");scanf("%lf %lf %lf",&w1,&w2,&w3);L2=F(w1,w2,w3);printf("B1=%f\n",B1);printf("L1=%.9f\n",L1);printf("B2=%.9f\n",B2);printf("L2=%.9f\n",L2);/*辅助计算*/W1=sqrt(1-ee*sin(B1)*sin(B1));W2=sqrt(1-ee*sin(B2)*sin(B2));sinu1=sin(B1)*sqrt(1-ee)/W1;sinu2=sin(B2)*(sqrt(1-ee))/W2;cosu1=cos(B1)/W1;cosu2=cos(B2)/W2;L=L2-L1;a1=sinu1*sinu2;a2=cosu1*cosu2;b1=cosu1*sinu2;b2=sinu1*cosu2;printf("W1=%.9f\n",W1);printf("W2=%.9f\n",W2);printf("sinu1=%.9f\n",sinu1);printf("sinu2=%.9f\n",sinu2);printf("cosu1=%.9f\n",cosu2);printf("L=%.9f\n",L);printf("a1=%.9f\n",a1);printf("a2=%.9f\n",a2);printf("b1=%.9f\n",b1);printf("b2=%.9f\n",b2);/*用逐次趋近法同时计算起点大地方位角、球面长度及经差R*/ do{R=L+Q;x=cosu2*sin(R);y=b1-b2*cos(R);A1=atan(x/y);if(x>0&&y>0)A1=fabs(A1);else if(x>0&&y<0)A1=I-fabs(A1);else if(x<0&&y<0)A1=I+fabs(A1);elseA1=2*I-fabs(A1);sinp=x*sin(A1)+y*cos(A1);cosp=a1+a2*cos(R);p=atan(sinp/cosp);if(cosp>0)p=fabs(p);sinA0=cosu1*sin(A1);cosA0=sqrt(1-sinA0*sinA0);p=I-fabs(p);z=2*a1-cosA0*cosA0*cos(p);m=(33523299-(28189-70*cosA0*cosA0))*(1e-10);n=(28189-94*cosA0*cosA0)*(1e-10);Q=q;q=(m*p-m*z*sin(p))*sinA0;}while(fabs(206265*(q-Q))>(1e-4));printf("Q=%f\n",Q);/*计算系数ABC及大地线长度*/A=6356863.020+(10708.949-13.474*cosA0*cosA0)*cosA0*cosA0;B=10708.938-17.956*cosA0*cosA0;C=4.487;Y=(cosA0*cosA0*cosA0*cosA0-2*z*z)*cos(p);S=A*p+(B*z+C*Y)*sinp;/*计算反方位角*/A2=atan((cosu1*sin(R))/(b1*cos(R)-b2));A2=A2*180/I;A1=A1*180/I;if(A1>0)A2=abs(A2);elseA2=-abs(A2);printf("A1=%3f\n A2=%3f\n S=%3f\n",A1,A2,S);}double F(double a,double b,double c){double d;if(a<0)d=(double)(a-1.0*b/60-1.0*c/3600);elsed=(double)(a+1.0*b/60+1.0*c/3600);d=(d/180)*I;return d; }运行结果：。