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

第31卷第4期2011年8月

大地测量与地球动力学

J OURNAL OF GEODESY AND GEODYNAM I CS

V o.l31N o.4

A ug.,2011

文章编号:1671-5942(2011)04-0094-05

基于第二类椭圆积分的子午线弧长公式变换及解算*过家春1) 赵秀侠2) 徐 丽1) 田劲松1) 高 飞3)

1)安徽农业大学理学院,合肥 230036

2)安徽大学资源与环境工程学院,合肥 230039

3)合肥工业大学土木与水利工程学院,合肥 230009

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

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

中图分类号:P226 文献标识码:A

CALCULAT I NG M ERID I AN ARC LENGTH BY TRANSFORM I NG

ITS FOR MULAE I NTO ELL I PTIC I NTEGRAL OF S ECOND K I ND

Guo Jiachun1),Zhao X i u x ia2),Xu L i1),T ian Ji n song1)and Gao Fe i3)

1)School of Science,AnhuiAgricultural Universit y,H efei 230036

2)School of R esources and Environ m ental Engineering,Anhu i University,H efei 230039

3)School of C ivil and H ydraulic Engineeri n g,H efei Un i v ersit y of Technology,H efei 230009

Abstract A ne w idea that by transfor m i n g the m er i d i a n arc leng th for m u la into t w o other for m s w as put for-w and,w h i c h are expressed by t h e e lli p tic i n teg rals o f the second kind,they w ere na m ed as For m and For m .In For m ,the m eridian arc length for m ula is represented as the sum o f a rational functi o n and the e lli p tic i n tegra ls of the second k i n d by the geodetic latitude para m eter,wh ich establishs the function relations bet w een the m eridian arc length and the e lli p tic i n tegra ls of the second ki n d.Analogousl y,taking the reduced latit u de as an i n-dependent variable,the For m for m ula give a si m pler for m of the m er i d ian arc leng t h for m ula i n ter m s o f t h e e-l li p tic i n tegrals of t h e second k i n d.On these bases of theoretica l analysis,the co m puter progra m is co m p iled i n MATLAB by calli n g t h e E lli p tic E(x,k)Functi o n to calcu late the m er i d i a n arc length.It w as proved that t h is m ethod i m proved greatly the accuracy and effic iency of prev ious calcu lation.M oreover,the m eridian arc length o f CGCS2000w as a lso calcu lated based on the pri n ciple that prov ided.The resu lts i n d icate that the For m and

*收稿日期:2011-02-18

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

作者简介:过家春,男,1981年生,硕士,讲师,主要研究大地测量学.E-m ai:l guojiachun@ah https://www.doczj.com/doc/c95570087.html,

通讯作者:赵秀侠,女,1981年生,博士,主要研究方向为地图学与地理信息系统.E-m ai:l xi ux i a99@126.co m

第4期过家春等:基于第二类椭圆积分的子午线弧长公式变换及解算

For m f o r m ula are si m pler and m ore su itab le f o r series expansi o ns and t h e rea lizati o n on co m puter t h an the classica l for m .Further m ore ,it perfects the m eri d ian theo r y .

K ey w ords :m eri d ian arc length ;for m ula transfor m ati o n ;elli p tic i n tegra ls ;geodetic latitude ;reduced latitude

1 引言

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

S =

B

M d B =

B a (1-e 2

)

(1-e 2

si n B )

3

2

d B (1)

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

显然,子午线弧长的计算涉及到勒让德第二类椭圆积分(简称为第二类椭圆积分),其原函数无法用初等函数的形式表达,不能直接求出。目前,世界

各国对子午线弧长的经典计算方法是将a(1-e 2

)(1-e 2

si n B )-3

2按二项式定理展开为级数形式,然后

再逐项积分,得出一定精度的计算公式[1,2]

:S =a(1-e 2

){A 0B -B 02si n 2B +C 04si n 4B -D 06

si n 6B + }(2)

式中,A 0、B 0、C 0、D 0、 为一组与e 有关的常系数,具体表达式可参阅文献[2]。

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

为了使子午线弧长公式的表达式更为简洁,或达到精度更高、便于计算机实现等目的,不少学者对此展开了研究

[3-9]

。但从其研究成果来

看,大都仍停留在以研究子午线弧长的计算为目

的的表面问题上,尚未探求出子午线弧长公式与第二类椭圆积分标准形式之间的内在本质关系,这在一定程度上影响了子午线弧长公式理论上的完整性,也限制了椭圆积分研究成果在子午线弧长计算问题上的应用。

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

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

一般椭圆积分定义为

[10-13]

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

y 2=P (x )=ax 4+b x 3+cx 2

+dx +e

(4)

椭圆积分即求椭圆的弧长[10,11]

。勒让德证明

了一般椭圆积分能够化成3种基本类型。其中,第

二类椭圆积分的标准形式为

E ( ,k )=

1-k 2sin 2

d

(5)

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

为

E (x,k )=

k 1-k 2t

2

1-t

2d t (6)

式中,k 为椭圆模,且0

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

2

,x =si n =1时,称式(5)、(6)为第二类完全椭圆积分,记为E 、E (k )、E (

2

,

k )或E (1,k )。

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

x =a sin

y =b cos

(0

(7)

令k =a 2

-b 2

/a,容易证明,式(5)所表达的

第二类椭圆积分E ( ,k )的几何意义为:以式(7)为参数方程,当a =1,0

图1 子午圈F i g .1 M er i d i an

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

明

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

E ( ,k )=

k 2

sin cos

1-k 2

si n 2

+

(1-k 2

)

d (1-k 2

si n 2 )

3/2

(8)

证明如下:因为

95

大地测量与地球动力学31卷

d d k

k 2

si n cos 1-k 2

sin 2

=

k 2

cos 2

1-k 2si n 2

-k 2

si n 2

1-k 2si n 2

+k 4

si n 2

cos 2

(1-k 2sin 2 )3/2=k 2-2k 2sin 2 +k 4

si n 4

(1-k 2si n 2 )3/2

(9)

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

si n cos

1-k 2si n 2

+(1-k 2

)

d (1-k 2

sin 2 )

3/2=

k 2

-2k 2

si n 2

+k 4

sin 4

(1-k 2si n 2 )

3/2

d +(1-k 2

)

d

(1-k 2sin 2 )

3/2=

1-2k 2

si n 2

+k 4

si n 4

(1-k 2si n 2 )

3/2

d = 0 (1-k 2

si n 2

)

2

(1-k 2

si n 2

)3/2

d =

1-k 2

si n 2

d =E ( ,k )

(10)

即等式(8)成立。

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

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

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

式变换

令k =

a 2-

b 2

/a,其中,0

并两边同时乘以a,移项后可得a(1-k 2

)

d

(1-k 2

si n 2

)

1-k 2sin 2

=-

ak 2

si n cos

1-k 2

sin 2

+aE ( ,k )

(11)

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

S =

B

M d B =-

ae 2

sin B cos B 1-e 2

si n 2

B

+aE (B,e)(12)

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

S =

0B

M d B =-sin B cos B (a 2-b 2)

a 2cos 2B +

b 2sin 2

B

+aE (B,e)(13)

特别地,当B =90 时,

S 90 =aE (

2

,e)

(14)

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

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

式变换

如图1所示,P 点的大地纬度为B,归化纬度为

=

2

- 。其中,归化纬度 与大地纬度B 之间的关系为

tan =

1-e 2

tan B =

b a

tan B (15)

因此有

d B =a(1+tan 2

B )b(1+tan 2 )d =a (1+a 2

b 2tan 2

)b (1+tan 2

)d (16)对式(1)作变量代换 =arctan (

b

a

tan B )得S =

B

M d B = 0

B

a (1-e 2

)

(1-e 2sin 2B )

3/2d

B =

a (1-e 2

)(1-e 2a 2si n 2 a 2si n 2 +b 2cos 2

)

3/2

a (1+a 2

b

2tan 2

)

b (1+tan 2

)

d (17)

化简后为

S =

a 2si n 2 +

b 2cos 2

d

(18)

再由 =

2

- ,可得S =

a 2

sin 2

+b 2

cos 2

d

- /2 a 2

cos 2

+b 2

si n 2

d = 0 /2

a 2

cos 2

+b 2

sin 2

d - 0

a 2

cos 2

+b 2

si n 2

d =a 0

/2

1-e 2

si n 2

d -a 0

1-e 2

si n 2

d =

aE (

2

,e)-aE ( ,e)=aE (

2,e)-aE ( 2

- ,e)(19)

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

到。

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

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

96

第4期过家春等:基于第二类椭圆积分的子午线弧长公式变换及解算称为 形式 。

4 基于 形式 、 形式 的子午线

弧长计算与分析

以I U GG75椭球参数为例,笔者在M atlab软件

中调用第二类椭圆积分函数E lliptic E(x,k)[14,15],

分别编写了基于 形式 和 形式 的子午线弧

长解算程序,实现了子午线的弧长计算功能。M a-t

lab中数据显示精度可以任意设置,本文取至10-8

m。表1为基于 形式 和 形式 的计算结果

与 经典算法 的结果对照。

表1 基于不同计算方法的子午线弧长计算结果

Tab.1 Co m parison a m ong the m eridian arc lengths co m-

puted w ith d ifferen t algorith m

大地纬度B (归化纬度 )I UGG75椭球子午线弧长S(m)经典算法本文算法*

10 00 00

(09 58 01.7231 )

1105855.34791105855.34788751

20 00 00

(19 56 17.6474 )

2212367.28432212367.28427513

30 00 00

(29 55 00.2915 )

3320114.94503320114.94499593

40 00 00

(39 54 18.9975 )

4429531.09644429531.09638841

50 00 00

(49 54 18.7985 )

5540849.62905540849.62902309

60 00 00

(59 54 59.7878 )

6654075.93056654075.93045982

70 00 00

(69 56 17.0746 )

7768984.36447768984.36442706

80 00 00

(79 58 01.3492 )

8885144.0358*******.03581356

90 00 00

(90 00 00 )

10001970.421110001970.42122640注:本文算法指基于 形式 和 形式 调用第二类椭圆积分的算法,二者结果完全一致。

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

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

直观上,图3中曲线接近于直线,这正反映了地

球扁率很小、平均曲率半径很大的特点。

图2 不同子午线弧长计算方法的CPU运行时间对比

F i g.2 Co m pa rison bet ween t he CPU ti m es fo r different pro-

cedures t o com pute the m er i d i an arc leng t h

图3 子午线弧长随大地纬度变化的曲线图

F i g.3

G raph of the m er i dian arc l eng t h chang i ng w ith G eo-

deti c La titude

为方便读者应用,下面给出基于 形式 M a-t lab程序代码:

clear

a=6378140;%椭球参数赋值

b=6356755.2881575287;

e=sqrt((a^2-b^2)/a^2);

for i=1:10;

B(i)=(i-1)*p i/18;%以10 为步长

S(i)=a*m fun( E lliptic E ,sin(B(i)),e)-a*e^2 *si n(B(i))*cos(B(i))/(1-e^2*si n(B(i))^2)^ (1/2);%子午线弧长计算

end

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

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

5 结论

形式 和 形式 的两种子午线弧长公式将子午线弧长的计算转化为第二类椭圆积分的计算,简洁直观,便于手工计算及计算机程序实现。其中,按 形式 进行子午线弧长计算时,需分两步进行,即先按式(15)由大地纬度为B计算出归化纬

97

大地测量与地球动力学31卷

表2 CGCS2000椭球子午线弧长

Tab.2 M erid i an arc length of the CGCS2000e llipsoid

大地纬度B

(归化纬度 )

CGCS2000椭球子午线弧长(m)

10 00 00

(09 58 01.7231 )

1105854.83319845

20 00 00

(19 56 17.6474 )

2212366.25410298

30 00 00

(29 55 00.2915 )

3320113.39784502

40 00 00

(39 54 18.9975 )

4429529.03023659

50 00 00

(49 54 18.7985 )

5540847.04156097

60 00 00

(59 54 59.7878 )

6654072.81936744

70 00 00

(69 56 17.0746 )

7768980.72765552

80 00 00

(79 58 01.3492 )

8885139.87183676

90 00 00

(90 00 00 )

10001965.72923050

度 ,然后再按式(19)进行计算,这对于计算机程序

设计来说是容易实现的,但不如 形式 简洁。而

且当B= =90 时,两种形式的公式直接将子午线弧长表达为第二类完全椭圆积分形式,这是 经典算法 所不具备的。

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

已有学者指出,子午线弧长的反解问题是目前

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

对于基于 形式 及 形式 的子午线弧长

反解问题需要进一步研究。

参考文献

1 W o lfgang T o rge.G eodesy(3rd.)[M].Berli n:W a lter D e

G ruyte r,2001.

2 孔祥元,郭际明,刘宗泉.大地测量学基础[M].武汉:武

汉大学出版社,2001.(K ong X i angyuan,G uo Ji m i ng and L i u Zongquan.F oundati on o f geodesy[M].W uhan:W uhan U n i versity P ress,2001)

3 姜晨光,阎桂峰.椭球子午线弧长计算的新方法[J].测

绘工程,1998,7(4):38-42.(Jiang Chengguang and Y an

G u if eng.A novel m ethod for the ca lcu l a tion of m eridian arc

length of ellipsoid[J].Eng i neer i ng o f Survey i ng M appi ng, 1998,7(4):38-42)

4 牛卓立.以空间直角坐标系为参数的子午线弧长计算公

式[J].测绘通报,2001,11:14-15.(N i u Zhuo l.i F or mu-lae fo r calculati on o fm eridian arc length by t he para m eters of space rectangular coordi nates[J].Bull e tin o f Surveying M ap-p i ng,2001,11:14-15)

5 刘修善.计算子午线弧长的数值积分法[J].测绘通报,

2006,5:4-6.(L i u X i ushan.N u m er ica l i nteg ra l me t hod f o r ca lcu l a ting m eridian arc length[J].Bu lleti n of Survey i ng M app i ng,2006,5:4-6)

6 易维勇,边少峰,朱汉泉.子午线弧长的解析型幂级数确

定[J].测绘学院学报,2000,17(3):167-171.(Y i W e i yong,B i an Shao f eng and Zhu H anquan.D eter m i na ti on o f foo t po i nt latitude by analytic positi ve ser i res[J].Journa l o f Insti tute of Survey i ng M apping,2000,17(3):167-171)

7 边少峰,许江宁.计算机代数系统与大地测量数学分析

[M].北京:国防工业出版社,2004.(B ian Shao feng and Xu Jiangni ng.Co m puter a l gebra system and m at hema ti ca l a-na l ys i s i n geodesy[M].Beiji ng:N ationa lD efence Industr i a l P ress,2004)

8 刘仁钊,伍吉仓.任意精度的子午线弧长递归计算[J].

大地测量与地球动力学,2007,(5):59-62.(L iu R enzhao and W u J i cang.R ecursi ve co m puta tion of m eri d ian arc l eng th w it h disc retionary prec i sion[J].Journa l o fG eodesy and G eo-dyna m ics,2007,(5):59-62)

9 K laus H eh.l C++and Java code for recursion for m ulas in

m athem ati ca l g eodesy[J].GPS So l ution,2005,9:51-58.

10 王竹溪,郭敦仁.特殊函数概论[M].北京:北京大学出

版社,2000.(W ang Zhux i and Guo Dunren.Introduction to special functi on[M].Be ijing:Pe iji ng U niversity P ress, 2000)

11 陆源.特阿贝尔对椭圆函数论的贡献[D].内蒙古师范

大学,2009.(L u Y uan.A study on contr i bution of A be l to the theo ry o f e lli ptic functi on[D].Inner M ongo li a N or m a l U n i versity,2009)

12 A bra m ow itz M and Stegun I A.H andbook o f m athema ti ca l

functions[M].N e w Y ork:Dover Pub licati ons,1964.

13 刘式适,刘式达.特殊函数[M].北京:气象出版社,

1988.(L i u Sh i sh i and L i u Shida.Spec i a l f unction[M].

Be iji ng:Ch i na M eteoro log i ca l P ress,1988)

14 The M ath W o rks Inc.M ATLA B7.0(R14SP2).The M at h-

W o rks Inc.,2005.

15 程鹏飞,等.2000国家大地坐标系椭球参数与GRS80和

WG S84的比较[J].测绘学报,2009,38(6):189-194.

(Cheng Peng fe,i et a.l P ara m eters of CGCS2000ellipsoid and compar i sons w it h GRS80and WG S84[J].A cta G eo-daetica et Cartographica S i n i ca,2009,38(6):189-194)

98