下载文档原格式
高斯平均引数大地主题解算程序设计田桂娥;谢露;马广涛【摘要】介绍了大地主题解算基本方法,以Visual Studio 2010作为开发平台,采用C#语言编写,设计实现了一套基于高斯平均引数的大地主题解算实用程序,指出该程序功能强大,实现了单点和批量大地主题解算,能方便的将计算结果保存在文本,且精度满足大地测量及相应工程的需求.【期刊名称】《山西建筑》【年(卷),期】2016(042)027【总页数】2页(P192-193)【关键词】大地主题解算;高斯平均引数法;Visual Studio 2010;批量解算【作者】田桂娥;谢露;马广涛【作者单位】华北理工大学,河北唐山063000;中铁十一局集团第一工程有限公司,湖北襄阳441104;河北省制图院,河北石家庄050000【正文语种】中文【中图分类】TU198在天文大地测量中，为了获得点的大地坐标，需要在椭球面上进行控制点间的坐标解算。
椭球面上两点间的大地坐标(大地经度、大地纬度)、大地线、大地方位角称为大地元素，已知一些大地元素，推求另一些大地元素，通常称为大地主题解算。
大地主题解算包含大地主题正算和大地主题反算两种，已知椭球面上一点P1的大地坐标(L1，B1)，P1到P2点的大地线长度S及其大地方位角A12，计算P2点的大地坐标(L2，B2)和大地线长度S在P2点的反方位角A21，称为大地主题正解；反之，已知P1和P2点的大地坐标(L1，B1)和(L2，B2)，计算P1至P2的大地线长度S及其正反方位角A12和A21，称为大地主题反解。
由于椭球计算的复杂性，带来大地主题解算的复杂性，有的需要进行迭代计算逐步趋近，给人工计算带来极大困难。
随着计算机技术的飞速发展，计算机在大地主题解算上的应用也得到了快速的发展，迭代计算已经不再是难题，而且，可以根据精度的需要而自行确定迭代次数，极大的提高了计算效率。
同时，随着大地主题解算在空间技术领域的广泛运用，大地主题解算已经成为一项重要的研究工作。
该投影按照投影带中央子午线投影为直线且长度不变和赤道投影为直线的条件，确定函数的形式，从而得到高斯一克吕格投影公式。
投影后，除中央子午线和赤道为直线外，其他子午线均为对称于中央子午线的曲线。
设想用一个椭圆柱横切于椭球面上投影带的中央子午线，按上述投影条件，将中央子午线两侧一定经差范围内的椭球面正形投影于椭圆柱面。
将椭圆柱面沿过南北极的母线剪开展平，即为高斯投影平面。
取中央子午线与赤道交点的投影为原点，中央子午线的投影为纵坐标x轴，赤道的投影为横坐标y轴，构成高斯克吕格平面直角坐标系。
高斯-克吕格投影在长度和面积上变形很小，中央经线无变形，自中央经线向投影带边缘，变形逐渐增加，变形最大之处在投影带内赤道的两端。
由于其投影精度高，变形小，而且计算简便（各投影带坐标一致，只要算出一个带的数据，其他各带都能应用），因此在大比例尺地形图中应用，可以满足军事上各种需要，能在图上进行精确的量测计算。
高斯-克吕格投影分带 按一定经差将地球椭球面划分成若干投影带,这是高斯投影中限制长度变形的最有效方法。
分带时既要控制长度变形使其不大于测图误差，又要使带数不致过多以减少换带计算工作，据此原则将地球椭球面沿子午线划分成经差相等的瓜瓣形地带,以便分带投影。
通常按经差6度或3度分为六度带或三度带。
六度带自0度子午线起每隔经差6度自西向东分带，带号依次编为第1、2…60带。
三度带是在六度带的基础上分成的，它的中央子午线与六度带的中央子午线和分带子午线重合，即自 1.5度子午线起每隔经差3度自西向东分带，带号依次编为三度带第1、2…120带。
我国的经度范围西起73°东至135°，可分成六度带十一个，各带中央经线依次为75°、81°、87°、……、117°、123°、129°、135°，或三度带二十二个。
六度带可用于中小比例尺（如1：250000）测图，三度带可用于大比例尺（如1：10000）测图，城建坐标多采用三度带的高斯投影。
《大地测量学基础》习题与思考题一 绪论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．简述地球大气中平流层、对流层与电离层的概念。
该投影按照投影带中央子午线投影为直线且长度不变和赤道投影为直线的条件，确定函数的形式，从而得到高斯一克吕格投影公式。
投影后，除中央子午线和赤道为直线外，其他子午线均为对称于中央子午线的曲线。
设想用一个椭圆柱横切于椭球面上投影带的中央子午线，按上述投影条件，将中央子午线两侧一定经差范围内的椭球面正形投影于椭圆柱面。
将椭圆柱面沿过南北极的母线剪开展平，即为高斯投影平面。
取中央子午线与赤道交点的投影为原点，中央子午线的投影为纵坐标x轴，赤道的投影为横坐标y轴，构成高斯克吕格平面直角坐标系。
高斯-克吕格投影在长度和面积上变形很小，中央经线无变形，自中央经线向投影带边缘，变形逐渐增加，变形最大之处在投影带内赤道的两端。
由于其投影精度高，变形小，而且计算简便（各投影带坐标一致，只要算出一个带的数据，其他各带都能应用），因此在大比例尺地形图中应用，可以满足军事上各种需要，能在图上进行精确的量测计算。
高斯-克吕格投影分带 按一定经差将地球椭球面划分成若干投影带,这是高斯投影中限制长度变形的最有效方法。
分带时既要控制长度变形使其不大于测图误差，又要使带数不致过多以减少换带计算工作，据此原则将地球椭球面沿子午线划分成经差相等的瓜瓣形地带,以便分带投影。
通常按经差6度或3度分为六度带或三度带。
六度带自0度子午线起每隔经差6度自西向东分带，带号依次编为第1、2…60带。
三度带是在六度带的基础上分成的，它的中央子午线与六度带的中央子午线和分带子午线重合，即自 1.5度子午线起每隔经差3度自西向东分带，带号依次编为三度带第1、2…120带。
我国的经度范围西起73°东至135°，可分成六度带十一个，各带中央经线依次为75°、81°、87°、……、117°、123°、129°、135°，或三度带二十二个。
六度带可用于中小比例尺（如1：250000）测图，三度带可用于大比例尺（如1：10000）测图，城建坐标多采用三度带的高斯投影。
大地主题的数值解法范业明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 引言一直以来由于计算工具的限制,大地主题解算一般都采用级数展开形式,如短距离的高斯平均引数法,长距离的贝塞尔-赫尔默特方法等等。
大地测量编程实习报告班号：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");}。
