基于第二类椭圆积分的子午线弧长公式变换及解算＊

过家春 赵秀侠 徐 丽 田劲松 高 飞

(1)安徽农业大学理学院,合肥 230036 2)安徽大学资源与环境工程学院,合肥 230039 3)合肥工业大学土木与水利工程学院,合肥 230009)

基于第二类椭圆积分的子午线弧长公式变换及解算＊

过家春1)赵秀侠2)徐 丽1)田劲松1)高 飞3)

(1)安徽农业大学理学院,合肥 230036 2)安徽大学资源与环境工程学院,合肥 230039 3)合肥工业大学土木与水利工程学院,合肥 230009)

通过推导,将子午线弧长公式变换为基于第二类椭圆积分的两种形式：“形式Ⅰ”将子午线弧长公式表达为一个有理函数和第二类椭圆积分之和,建立了以大地纬度B为自变量的子午线弧长公式与第二类椭圆积分之间的关系;“形式Ⅱ”给出了以归化纬度μ为自变量、直接利用第二类椭圆积分计算子午线弧长的公式。利用此两种形式的子午线弧长公式,在Matlab中编写程序,调用第二类椭圆积分函数 EllipticE(x,k)计算子午线弧长,精度和计算效率均优于经典算法。对 CGCS2000所采用的地球椭球子午线弧长的计算表明,此两种形式的子午线弧长公式建立了子午线弧长公式与第二类椭圆积分的关系,结构简洁,易于展开,一定程度上完善了子午线弧长理论,且便于手工计算及计算机程序实现。

子午线弧长;公式变换;椭圆积分;大地纬度;归化纬度

1 引言

子午线弧长是大地测量学中的一个基本量。计算从赤道开始到任意大地纬度 B的子午线弧长 S的基本公式为

式中,a为地球椭球长半轴;e为地球椭球第一偏心率;M为子午圈曲率半径,B为大地纬度。

显然,子午线弧长的计算涉及到勒让德第二类椭圆积分(简称为第二类椭圆积分),其原函数无法用初等函数的形式表达,不能直接求出。目前,世界各国对子午线弧长的经典计算方法是将 a(1-e2) (1-e2sinB)按二项式定理展开为级数形式,然后再逐项积分,得出一定精度的计算公式[1,2]：

式中,A0、B0、C0、D0、…为一组与 e有关的常系数,具体表达式可参阅文献[2]。

在此基础上,可以计算任意区间 [B1,B2]上的子午线弧长(本文称之为“经典算法”)。

为了使子午线弧长公式的表达式更为简洁,或达到精度更高、便于计算机实现等目的,不少学者对此展开了研究[3-9]。但从其研究成果来看,大都仍停留在以研究子午线弧长的计算为目的的表面问题上,尚未探求出子午线弧长公式与第二类椭圆积分标准形式之间的内在本质关系,这在一定程度上影响了子午线弧长公式理论上的完整性,也限制了椭圆积分研究成果在子午线弧长计算问题上的应用。

2 第二类椭圆积分及其展开式

2.1 第二类椭圆积分的标准形式

一般椭圆积分定义为[10-13]

其中,R(x,y)为 x和 y的有理函数,而

椭圆积分即求椭圆的弧长[10,11]。勒让德证明了一般椭圆积分能够化成 3种基本类型。其中,第二类椭圆积分的标准形式为

与此等价,作变量代换 x=sinφ,另外一种记法为

式中,k为椭圆模,且 0＜k＜1;φ称为幅度。

通常,称式(5)、(6)为第二类不完全椭圆积分。

第二类椭圆积分的几何意义即为椭圆的弧长。设一椭圆的参数方程为

图1 子午圈Fig.1 Meridian

2.2 第二类椭圆积分的基本关系式及其证明

讨论过程中用到的第二椭圆积分的关系式为：

证明如下：

因为

所以对于等式(8)的右边有

即等式(8)成立。

3 基于第二类椭圆积分的子午线弧长公式变换

由于旋转椭球上的子午圈为椭圆,所以子午线弧长的计算问题也即椭圆弧长问题,这就决定了子午线弧长公式必然可以变换为第二类椭圆积分形式。

3.1 基于第二类椭圆积分的子午线弧长公式变换Ⅰ

将式(11)中的 k改写为 e,φ改写为 B,得子午线弧长的第二类椭圆积分的表达形式为：

式(12)可进一步化简化为：

特别地,当B=90°时,

式 (12)和式 (13)将子午线弧长公式表达为一个有理函数和第二类椭圆积分之和,建立了以大地纬度为自变量的子午线弧长公式与第二类椭圆积分之间的关系,现简称为“形式Ⅰ”。

3.2 基于第二类椭圆积分的子午线弧长公式变换Ⅱ

因此有

化简后为

式(19)也可由第二类椭圆积分的几何意义直接得到。

特别地,当μ=90°时,式(19)化为式(14)。

式(19)给出了以归化纬度μ为自变量的直接利用第二类椭圆积分计算子午线弧长的公式,现简称为“形式Ⅱ”。

4 基于“形式Ⅰ”、“形式Ⅱ”的子午线弧长计算与分析

以 I UGG75椭球参数为例,笔者在Matlab软件中调用第二类椭圆积分函数 EllipticE(x,k)[14,15],分别编写了基于“形式Ⅰ”和“形式Ⅱ”的子午线弧长解算程序,实现了子午线的弧长计算功能。Matlab中数据显示精度可以任意设置,本文取至 10-8m。表 1为基于“形式Ⅰ”和“形式Ⅱ”的计算结果与“经典算法”的结果对照。

表1 基于不同计算方法的子午线弧长计算结果Tab.1 Comparison among the meridian arc lengths computed with different algorithm

表1显示,基于“形式Ⅰ”和“形式Ⅱ”的计算结果完全一致。由于其结果是在椭圆积分真值的基础上计算而得的,因此可视为相应纬度的子午线弧长真值。在Matlab、Mathematica、Maple等数学软件的特殊函数库、C++的 Boost库等,计算第二类椭圆积分的内部算法为卡尔松方法,计算效率高于按二项式定理的级数展开方法,图 2为在Matlab中分别利用本文算法、二项式定理展开方法同时计算子午线弧长的 CPU执行时间对比,结果表明,本文算法的计算效率提高到了 10倍左右。

程序还分别绘制了子午线弧长随大地纬度变化的曲线图(图 3)。

直观上,图 3中曲线接近于直线,这正反映了地球扁率很小、平均曲率半径很大的特点。

图2 不同子午线弧长计算方法的 CPU运行时间对比Fig.2 Comparison between the CPU ti mes for differentprocedures to compute the meridian arc length

图3 子午线弧长随大地纬度变化的曲线图Fig.3 Graph of the meridian arc length changingwith Geodetic Latitude

为方便读者应用,下面给出基于“形式Ⅰ”Matlab程序代码：

该程序中用于子午线弧长计算的代码仅一行,且精度没有损失。基于“形式Ⅱ”的程序结构与其类似,计算结果完全一致。

将程序中的椭球参数替换为中国 2000国家大地坐标系 (CGCS2000)所采用的地球椭球参数[15],即可计算得到 CGCS2000地球椭球的子午线弧长,结果如表2所示。

5 结论

“形式Ⅰ”和“形式Ⅱ”的两种子午线弧长公式将子午线弧长的计算转化为第二类椭圆积分的计算,简洁直观,便于手工计算及计算机程序实现。其中,按“形式Ⅱ”进行子午线弧长计算时,需分两步进行,即先按式(15)由大地纬度为B计算出归化纬度μ,然后再按式(19)进行计算,这对于计算机程序设计来说是容易实现的,但不如“形式Ⅰ”简洁。而且当B=μ=90°时,两种形式的公式直接将子午线弧长表达为第二类完全椭圆积分形式,这是“经典算法”所不具备的。

表2 CGCS2000椭球子午线弧长Tab.2 M eridian arc length of the CGCS2000 ellipsoid

由于许多数学软件、程序语言的函数库中自带第二类椭圆积分函数,如Matlab、Mathematica、Maple等数学软件的特殊函数库、C++的 Boost库等,在编写程序时按“形式Ⅰ”或“形式Ⅱ”的形式实现,直接调用其第二类椭圆积分函数即可,代码简短高效。

已有学者指出,子午线弧长的反解问题是目前仍然没得到完美解决的问题[6]。“形式Ⅰ”和“形式Ⅱ”两种形式的公式可以建立起子午线弧长公式与其他特殊函数的关系,如超几何函数、雅氏椭圆函数等[10-13],有望为子午线弧长的反解问题提供新的思路。尤其是“形式Ⅱ”式 (19)的第一项对于一定的地球椭球为一常数,将子午线弧长的反解问题直接转化为了第二类椭圆积分的求逆问题。

对于基于“形式Ⅰ”及“形式Ⅱ”的子午线弧长反解问题需要进一步研究。

1 Wolfgang Torge.Geodesy(3rd.)[M].Berlin：Walter De Gruyter,2001.

2 孔祥元,郭际明,刘宗泉.大地测量学基础[M].武汉：武汉大学出版社,2001.(Kong Xiangyuan,Guo Jiming and Liu Zongquan.Foundation of geodesy[M].Wuhan：Wuhan University Press,2001)

3 姜晨光,阎桂峰.椭球子午线弧长计算的新方法[J].测绘工程,1998,7(4)：38-42.(Jiang Chengguang and Yan Guifeng.A novel method for the calculation of meridian arc length of ellipsoid[J].Engineering of Surveying Mapping, 1998,7(4)：38-42)

4 牛卓立.以空间直角坐标系为参数的子午线弧长计算公式[J].测绘通报,2001,11：14-15.(Niu Zhuoli.Formulae for calculation ofmeridian arc length by the parametersof space rectangular coordinates[J].Bulletin of SurveyingMapping,2001,11：14-15)

5 刘修善.计算子午线弧长的数值积分法[J].测绘通报, 2006,5：4-6.(Liu Xiushan.Numerical integral method for calculating meridian arc length[J].Bulletin of Surveying Mapping,2006,5：4-6)

6 易维勇,边少峰,朱汉泉.子午线弧长的解析型幂级数确定[J].测绘学院学报,2000,17(3)：167-171.(Yi Weiyong,Bian Shaofeng and Zhu Hanquan.Determination of foot point latitude by analytic positive serires[J].Journal of Institute of Surveying Mapping,2000,17(3)：167-171)

7 边少峰,许江宁.计算机代数系统与大地测量数学分析[M].北京：国防工业出版社,2004.(Bian Shaofeng and Xu Jiangning.Computer algebra system and mathematical analysis in geodesy[M].Beijing：NationalDefence Industrial Press,2004)

8 刘仁钊,伍吉仓.任意精度的子午线弧长递归计算[J].大地测量与地球动力学,2007,(5)：59-62.(Liu Renzhao andWu Jicang.Recursive computation ofmeridian arc length with discretionary precision[J].Journalof Geodesy and Geodynamics,2007,(5)：59-62)

9 Klaus Hehl.C++and Java code for recursion formulas in mathematical geodesy[J].GPS Solution,2005,9：51-58.

10 王竹溪,郭敦仁.特殊函数概论[M].北京：北京大学出版社,2000.(Wang Zhuxi and Guo Dunren.Introduction to special function[M].Beijing：Peijing University Press, 2000)

11 陆源.特阿贝尔对椭圆函数论的贡献[D].内蒙古师范大学,2009.(Lu Yuan.A study on contribution of Abel to the theory of elliptic function[D].InnerMongolia Nor mal University,2009)

12 AbramowitzM and Stegun IA.Handbook of mathematical functions[M].New York：Dover Publications,1964.

13 刘式适,刘式达.特殊函数 [M].北京：气象出版社, 1988.(Liu Shishi and Liu Shida.Special function[M]. Beijing：ChinaMeteorological Press,1988)

14 TheMathWorks Inc.MATLAB 7.0(R14SP2).TheMath-Works Inc.,2005.

15 程鹏飞,等.2000国家大地坐标系椭球参数与 GRS80和WGS84的比较[J].测绘学报,2009,38(6)：189-194. (Cheng Pengfei,et al.Parameters of CGCS2000 ellipsoid and comparisonswith GRS80 and WGS84[J].Acta Geodaetica et Cartographica Sinica,2009,38(6)：189-194)

CALCULATING M ERI D IAN ARC LENGTH BY TRANSFORM ING ITS FORM ULAE INTO ELL IPTIC INTEGRAL OF SECOND KIND

Guo Jiachun1),Zhao Xiuxia2),Xu Li1),Tian Jinsong1)and Gao Fei3)

(1)School of Science,Anhui Agricultural University,Hefei 230036 2)School of Resources and Environm ental Engineering,Anhui University,Hefei 230039 3)School of Civil and Hydraulic Engineering,Hefei University of Technology,Hefei 230009)

A new idea that by transforming the meridian arc length for mula into two other forms was put forwand,which are expressed by the elliptic integrals of the second kind,theywere named as“FormⅠ”and“FormⅡ”.In“FormⅠ”,the meridian arc length formula is represented as the sum of a rational function and the elliptic integrals of the second kind by the geodetic latitude parameter,which establishs the function relations between the meridian arc length and the elliptic integrals of the second kind.Analogously,taking the reduced latitude as an independent variable,the“For mⅡ”formula give a si mpler form of themeridian arc length formula in ter msof the elliptic integrals of the second kind.On these bases of theoretical analysis,the computer program is compiled in MATLAB by calling the EllipticE(x,k)Function to calculate the meridian arc length. It was proved that this method improved greatly the accuracy and efficiency of previous calculation.Moreover,the meridian arc length of CGCS2000 was also calculated based on the principle that provided.The results indicate that the“FormⅠ”and“FormⅡ”for mula are simpler and more suitable for series expansions and the realization on computer than the classical form.Further more,it perfects the meridian theory.

meridian arc length;formula transformation;elliptic integrals;geodetic latitude;reduced latitude

1671-5942(2011)04-0094-05

2011-02-18

国家自然科学基金(40771117);国家农业信息化工程技术研究中心开放课题(KF2010W40-046)

过家春,男,1981年生,硕士,讲师,主要研究大地测量学.E-mail：guojiachun@ahau.edu.cn

赵秀侠,女,1981年生,博士,主要研究方向为地图学与地理信息系统.E-mail：xiuxia99@126.com

P226

A