两种子午线弧长计算方法的比较

1、两种大地子午线弧长计算方法的比较胡洋长安大学西安2604070210摘要：在研究与大地椭球体有关的一些测量计算时，例如研究高斯投影计算和 弧度测量计算时，往往要用到大地子午线弧长，而大地子午线弧长的计算公式 涉及到椭圆积分，不能用普通方法求出被积函数。采用变步长辛普森公式求它 的积分和按泰勒级数展开采用普通方法求出被积函数是两种截然不同的计算方 法，本文就两种不同的计算方法作出比较，用程序编程实现从而得到计算结果， 验证了两种方法的正确性和可靠性以及按泰勒级数展开的方法求解大地子午线弧长的优 之处。关键字：大地子午线弧长，变步长辛普森公式，泰勒级数。引言：计算地球椭球子午线的弧长，是大地测量、

2、天文测量、航空航天技术以及地理信息处理技术 中的一项基本内容。计算子午线的弧长涉及到椭圆积分，所以一般是将其展开成级数形式， 再用逐项积分的方法求出满足一定精度要求的计算公式。本文在引出子午线弧长的计算公式 之后，介绍了两种计算方法的原理，并用C语言程序实现了两种算法的电算，继而对两种 算出的结果进行了分析讨论。一、子午线弧长的计算公式我们知道，子午椭圆的一半，端点与极点相重合，而赤道又把子午线分成对称的两部分，因 此，推到从赤道开始到已知纬度B之间的子午线弧长的计算公式即可。如下图，今取子午线上的某微分弧段PP= d x，令P点纬度等于B, P点纬度为B+ dB, P 点的子午圈曲率半径为M

3、，于是有dj MdB，因此，为了计算从赤道开始到任意纬度B的“ r Ba（1 一e2）平行圈之间的弧长，必须求出下面的积分值X = j Md（式2）。又知M=（1 .）切（式3），故（式2）符合椭圆积分的定义（椭圆积分定义为可以表达为如下形式的任何函数x的积分：f （x） = j 0#t,P（t）dt，其中R是其两个参数的有理函数，P是一个无重根的3或 4阶多项式的平方根，而c是一个常数），因此（式1）为椭圆积分，无法用普通的积分方法 求出原函数。、两种方法计算子午线弧长的原理（一）变步长辛普森公式求子午线弧长的积分a (1 一 e 2)令M(B)=，计算步骤如下:V (1 - e2 sin2

4、 B)31、用梯形公式计算T = hM（0） + M（B）/2，其中n=1，h=B。且令Sn=Tnn2、3、1 _h n-1 一用变步长梯形法则计算T2广2 T2n + 2 M (B + h /2) k=0S = (4T -T )/3用辛普森公式计算2n 2 nn4、若IS - S l8，则令2 n n n，h/2 n h，转2继续进行计算；否则结束即为所求的 2n n积分值。8是事先给定的积分精度。a (1 - e 2)，（二）泰勒级数展开后按普通方法积分令 M(B)=-V (1 - e2 sin2 B)31、按泰勒级数展开如下:M - m + m sin2 B + m sin4 B + m

5、 sin6 B + m sins B ；2、为便于积分将正弦的幕函数展开为余弦的倍数函数，由于sin 2 B =1-1 cos 2 B2 2 TOC o 1-5 h z .八311 fsin 4 B = 一一一 cos 2 B + cos 4 B 2sin 6 B = - cos 2 B + cos 4 B- cos 6 B16 321632.八 357 f 7 f 11sin 8 B = cos 2 B + cos 4 B - cos 6 B +cos 8B128 1632161283、将其代入上式，并经过整理得到:M = a - a cos 2B + a cos 4B - a cos 6B

6、 + a cos8B 02468式中a = m + 2 + m + m +- m +.0028 4 16 6 128 8m157a =2 + m + m + m +.222 4 32 6 16 8137a = - m + m + m48 416 632 81284、将M = a - a cos 2B + a cos4B - a cos 6B + a cos8B代入到积分式中进行积分， 02468整理后得到 X = a B 一 sin 2B + = sin 4B 一 与 sin 6B + 8 sin 8B 02468三、两种方法计算子大地午线弧长的数据和算例的讨论笔者已将两种方法的C语言程序写出

7、，限于篇幅本文未列出具体程序。该程序在任何PC机 上均可运行。经过用VC+6.0编写的程序代码调试，选用30度，45度，60度，90度三种 不同纬度的点来求相应地点的大地子午线弧长。采用的1954年北京坐标系应用的克拉索夫 斯基椭球。具体的计算结果如下表：用变步长辛普森公式计算的结果X1按泰勒级数展开计算 的结果X2两者之间的差值X1-X2纬度B3030对应大地子午线 弧长X3320172.406714743320172.40672017-0.00000543纬度B4545对应大地子午线 弧长X4985032.290475264985032.29047666-0.0000014纬度B6060对

8、应大地子午线 弧长X6654189.092215926654189.092214220.0000017纬度B9090对应大地子午线 弧长X10002137.497542910002137.49745320.000089从上表可以看出，两种不同的数学方法计算大地子午线弧长的精度相当，差异只是在小数点 后边的第五位。在实际的应用过程当中，我们只要求精确到0.1mm的数量级别就完全够用 了。但是，我们从上表中也可以发现，两种方法虽然在求得的结果上略有不同，但是发散增 长的速度上差异还是比较大的。在中低纬度地区，由变步长辛普森公式计算的结果要小于由 泰勒级数展开计算的结果；但到了高纬度地区，两者计算的结果差变为正数。特别是在极点 处，两者相差达到接近0.1mm。这说明两者之间的二阶导数是步长辛普森公式的方法大于泰勒级数展开的。两者结果增长的稳定性上，按泰勒级数展开的方法求出的结果更加稳定可 靠。虽然这两种方法的精度略有不同，但是对