摘 要：子午线弧长计算的经典算法是对子午线曲率半径按照牛顿二项式定理进行展开，分项积分得到近似解析解。本文研究了五种常用的数值积分算法及其在子午线弧长计算中的应用，并用MATLAB软件予以实现。将数值积分结果与经典算法结果进行比较，结果表明：利用数值积分算法求解子午线弧长，简单易行，准确可靠。同时还证明了复合辛普森算法和龙贝格算法要优于其它几种算法。
关键词：子午线弧长 数值积分算法 MATLAB
引言
子午线弧长计算是椭球大地测量学的一项重要内容，是高斯投影坐标计算的基础。但因子午线弧长计算公式中的被积函数（即子午线曲率半径M）无法找到其原函数，故而不能直接积分，经典的解法是采用牛顿二项式定理进行展开，分项积分得到近似解析解[1]。
近些年，国内一些学者对该问题提出一些解算方法：文献[4]采用常微分方程的数值解法，得到了理想的结果；文献[5]基于第二类椭圆积分进行求解，结果理想但过程过于复杂，不便应用；文献[6]给出了任意精度的子午线弧长递归计算公式，可满足不同精度的弧长计算；文献[7]对子午线弧长的数值积分法进行了探讨，但得到的数值积分结果与经典算法结果相差几百米甚至上万米，这与数值积分算法的准确性相矛盾，不免让人有所疑惑。本文拟就数值积分算法求解子午线弧长进行研究，以验证此类算法究竟是否准确可靠。
一、子午线弧长计算的经典算法
大地测量学经典教材[1]给出了计算子午线弧长的基本公式：从赤道到大地纬度为B处的子午线弧长为 （1）
式中，M为子午线曲率半径，a为椭球长半轴，e为椭球第一偏心率。
对于（1）式而言，由于被积函数结构复杂，其原函数无法求出，故不能直接用牛顿-莱布尼茨积分进行计算。在经典教材中，采用近似解析法，即把被积函数M按照牛顿二项式定理展开为e的幂级数，并将正弦的幂函数展开为余弦的倍数函数，然后逐项进行积分，考虑到计算结果的目的和精度，得到以下两个实用计算公式：
M=m0+m2sin2B+m4sin4B+m6sin6B+m8sin8B （2）
（3）
（4） （5）
根据式（3）计算的结果已被诸多文献[4～6]证明了其准确性，故本文以经典算法的结果作为准确值，用以衡量各种数值积分算法的准确度。
3.数值积分算法
数值积分算法[2]为直接计算式（1）提供了可能性。由于数值积分是采用特定算法直接进行积分，这样就可以避开经典算法中的展开计算，且辅以MATLAB程序设计对算法进行实现，可大大减少计算工作量，得到高精度的计算结果。
数值积分算法有很多种，常用的主要有梯形算法、辛普森（Simpson）算法、复合梯形算法、复合辛普森算法、龙贝格（Romberg）算法、高斯-勒让德（Gauss-Legendre）算法以及蒙特卡罗（Monte Carlo）算法等，但梯形算法与辛普森算法已被证明计算结果精度有限[2]，故本文主要采用后五种算法分别进行子午线弧长的计算，并对各种算法的结果予以比较。
3.1 复合梯形算法
对于积分 而言，将积分区间[a，b]划分为n等分，步长 ，在每个子区间[xk，xk+1]（k=0，1，...，n-1）上分别采用梯形求积公式即可得到复合梯形算法计算公式：
（6）
复合梯形算法结果的准确度主要取决于步长h的大小，h越小，则计算结果越准确。
3.2 复合辛普森算法
类似复合梯形算法，若在每个子区间[xk，xk+1]上分别采用辛普森求积公式即可得到复合辛普森算法计算公式：
（7）
式中
3.3 龙贝格算法
对于复合梯形算法和复合辛普森算法而言，虽然计算过程简答，但收敛速度有限。龙贝格算法是在积分区间逐次分半的过程中，对用复合梯形算法和复合辛普森算法产生的近似值进行加权平均，不仅结果准确可靠，且迅速收敛。
由于龙贝格算法推导过程稍显复杂，这里直接给出最终的计算公式，详细推导过程请查阅文献[2]。
（8）
式中， 称为柯特斯（Cotes）公式， 。 Sn、Tn即为式（6）、（7）。
3.4 高斯-勒让德算法
n阶勒让德多项式定义为 （9）
高斯-勒让德求积公式为 （10）
其中xk为勒让德多项式在[-1，1]区间上的零点。权系数Ak为：
对于一般积分区间为[a，b]问题，可以做变换 ，得到
（11）
零点和权系数可以通过程序计算出来，也可以在很多计算手册中查找得到。
3.5 蒙特卡罗算法
蒙特卡洛算法类似于机械求积方法，采用随机数的方法产生结点，然后把各区间段上的函数值与区间段相乘后再相加。而实际上如果均匀分布的随机点比较密集，则随机数产生的区间近似均分整个积分区间段。具体方法如下：
对于积分 ，用函数产生0到1之间均匀分布的随机数t（1xn）向量，然后做区间转换xk=a+（b-a）tk，其中k=1，2...，n。如果随机数数目较大，则有 （12）
4、数值积分计算结果及分析
现以IUGG-1975椭球的子午线弧长为例进行计算，分别采用上述五种数值积分算法进行计算。计算工具选择数学计算软件MATLAB[3]，并根据各种算法原理及公式编写相应的子程序。除了龙贝格算法的程序代码为25行之外，其余四种算法的程序代码都在15行之内，非常简洁。同时在计算过程中还需要顾及到角度值与弧度值之间的转换，最后得到计算结果列于表1中，从表中可对比得到各种数值积分算法结果和经典算法结果间的差值（以下简称为差值）。
在计算过程中，对于复合梯形算法，步长h分别选择1o、1"和 ，1""可见当步长h=1o时，差值在2m以内；当h=1"时，差值不超过1mm；当h=1""时，差值为0mm，说明复合梯形算法要想提高计算结果精度，需要减小步长。当步长h=1""时，才能得到与经典算法相同的结果。
对于复合辛普森算法，通过计算发现，当步长h=1o时，差值为0mm，即可得到与经典算法相同的结果，这说明复合辛普森算法要优于复合梯形算法。
对于龙贝格算法，当限差△=10-3时（对应计算结果取到mm），差值为0mm，亦得到理想结果。
而高斯-勒让德算法当取到5阶时，差值在才1mm以内，而取到3阶或4阶的结果均与经典算法差值较大，故而不再列于表中。
蒙特卡罗算法当选择的n值为对应于h=1""换算的数值后，差值仍然为几十米，说明此种算法在子午线弧长计算中不适用，其计算结果亦未列出。
备注：①以上数值的单位为m；②椭球元素选用IUGG-1975椭球元素；③限于篇幅，考虑到各种数值积分算法结果的千米量级数值和经典算法结果数值一样，故在本表中略去了数值积分算法结果中大于千米的数据。
结束语
本文基于MATLAB软件，分别对复合梯形算法、复合辛普森算法、龙贝格算法、高斯-勒让德算法和蒙特卡洛算法在子午线弧长计算中的应用进行了探讨，通过与经典算法结果之间的对比，证明了数值积分算法的准确性和可靠性，说明文献[7]在计算过程中有失误之处。
基于计算结果及其分析比较，在利用数值积分进行子午线弧长计算时，复合辛普森算法和龙贝格算法要优于其它算法，建议在测量工作中实际计算时予以采用。
参考文献
[1].孔祥元，郭际明，刘宗泉.大地测量学基础[M].武汉：武汉大学出版社，2006
[2].李庆扬，王能超，毅大义.数值分析（第4版）[M].北京：清华大学出版社，2001