范文一：对弧长的曲线积分

对弧长的曲线积分

一、概念的引进

假设xoy 面内有一段曲线弧L 具有质量, 在L 上任一点(x , y ) 处的线密度为ρ(x , y ) , 且ρ(x , y ) 在L 上连续, A 与B 分别是弧L 的端点, 现计算弧L 的质量m 。

在L 上任意地插入n +1个分点

A =M 0, M 1, , M i -1, M i , , M n -1, M n =B

将L 分划成n 个小弧段。对于第 i 个小弧段弧M i -1M i , 由于线密度函数

ρ(x , y ) 在L 上连续, 当该小弧段的长度充分小时, 它的质量近似地等于

ρ(ξi , ηi ) ?s i

n

?(ξi , ηi ) 弧M i -1M i ,

?s i 表示弧M i -1M i 的长度

于是, 整个曲线弧L 的质量近似值为

m ≈∑ρ(ξi , ηi ) ??s i

i =1

用λ表示这n 个小弧段长度的最大者, 即

λ=max {?s i }

1≤i ≤n

为了得到质量m 的精确值, 只需对上述和式取极限, 令λ→0,

即

m =lim ∑ρ(ξi , ηi ) ??s i

λ→0i =1

n

(1)

撇开上例的物理意义, 我们引入对弧长的曲线积分的概念。 【定义】设L 为xoy 面内的一条光滑曲线弧, 函数任意地插入n +1点,

f (x , y ) 在L 上有界, 在L 内

A =M 0, M 1, , M i -1, M i , , M n -1, M n =B

它把L 分成n 个小弧段, 设第i 个小段弧M i -1M i 的长度为?s i , (ξi , ηi ) 为

弧M i -1M i 上任取的一点, 记

n

λ=max {?s i }

1≤i ≤n

作和式 i =1

∑f (ξi , ηi ) ??s i

如果极限 λ→0i =1

lim ∑f (ξi , ηi ) ??s i

存在,

n

这个极限值就叫做函数

f (x , y ) 在曲线弧L 上对弧长的曲线积分, 记作

L

?f (x , y ) ds

。

亦即 L 其中:

lim ∑f (ξi , ηi ) ??s i ?f (x , y ) ds =λ→0

i =1

n

f (x , y ) 叫做被积函数, L 叫做积分弧段。

注记:

1、L

?f (x , y ) ds

中的被积函数

f (x , y ) 的定义域为L 上的一切点。

f (x , y , z ) 在Γ上有界, 则

2、上述定义可类似地推广到空间曲线的情形,

设Γ是空间的一条光滑曲线, 函数

Γ

lim ∑f (ξi , ηi , ζi ) ??s i ?f (x , y , z ) ds =λ→0

i =1

n

。

3、若L 为一条封闭曲线, 一般将L

?f (x , y ) ds

记为 L

f (x , y ) ds

二、对弧长的曲线积分的性质

利用对弧长的曲线积分定义, 我们可以证明下述性质 1、L

?[f (x , y ) ±g (x , y )]ds =?f (x , y ) ds ±?g (x , y ) ds

L

L

2、若k 为常数, L

?k ?f (x , y ) ds =k ??f (x , y ) ds

L

3、L

?ds =L 的长度

4、若在L 上,

f (x , y ) ≤g (x , y ) , 则 L

?f (x , y ) ds ≤?g (x , y ) ds

L

L 1

L 2

1+L 2, 则 L 5、若L =L

?f (x , y ) ds =?f (x , y ) ds +?f (x , y ) ds

上述性质均不加以证明, 有兴趣的同学可以查阅有关书籍。

三、对弧长曲线积分的计算法

假设曲线L 由参数方程

x =?(t ) , y =φ(t ) (α≤t ≤β)

给出, 且函数?(t )

, φ(t ) 在[α, β]上具有一阶连续导数；函数f (x , y ) 在L

上连续；当参数t 由α变至β时, 依点A 至点B 的方向描出曲线L 。

在L 上取一系列的点

A =M 0, M 1, , M i -1, M i , M n -1, M n =B

设它们对应于一列单调增加的参数值

α=t 0<><>

依定义

L

lim ∑f (ξi , ηi ) ?s i ?f (x , y ) ds =λ→0

i =1

n

这里的(ξi 则 ξi

, ηi ) ∈弧M i -1M i , 并设点(ξi , ηi ) 对应于参数值τi

t i -1≤τi ≤t i

=?(τi ) , ηi =φ(τi )

t i

由弧长计算公式与定积分中值定理有

?s i =?

?'(t ) ]2+[φ'(t ) ]2dt

t i -1

=

从而

22) +φ(τ) ??t i ''?'(τ'][]i i

(t i -1≤τ'i ≤t i , ?t i =t i -t i -1)

2

2

L

?f (x , y ) ds =lim ∑f [?(τi ), φ(τi ?'(τ'i ) ]+[φ'(τ'i ) ]??t i

λ→0i =1

n

(2)

由于函数

?'(t ) ]2+[φ'(t ) ]2在[α, β]上连续, 在λ→0时, 小区间

[t i -1, t i ]的长度?t i →0(i =1, 2, , n ) 。 那么在[t i -1, t i ]上,

2222

?(τ) +φ(τ) '') +φ(τ) ''][]?'(τ'][]i i i i 与

只相差一个?t i 的高阶无穷小, 因此, 我们可以把(2)式右端的τ'换成τi , 有

22

lim ∑f [?(τi ), φ(τi ?'(τi ) ]+[φ'(τi ) ]??t i ?f (x , y ) ds =λ→0

i =1n

L

22

''f [?(t ) , φ(t ?(t ) +φ(t ) 而右端和式的极限, 就是函数

在区间

[α, β]上的定积分。由于函数是连续的, 故此定积分存在, 因此, 上式左端的曲线

积分亦存在, 且有

β

?f (x , y ) ds =?f [?(t ), φ(t ?'(t ) 2

+φ'(t ) 2

dt

L

α

(3)

强调指出, (3)式中的定积分下限α一定要小于上限β, 理由是 (2)式中的?t i 由表达式

?s i =

?'(τ'i ) ]2+[φ'(τ'i

) ]2

??t i 给出, 因小弧段的长度?s i >0, 从而

?t i >0?t i -t i -1>0?t i >t i -1(i =1, 2, , n )

因此

α=t 0<><>

利用(3)式, 可导出如下几种对弧长的曲线积分计算公式 1、曲线L 由方程

y =φ(x ) (a ≤x ≤b )

给出时,

b

?f (x , y ) ds =?f [x , φ(x +[φ'(x ) ]2

dx

L

a

2、曲线L 由方程

x =?(y ) (c ≤y ≤d )

给出时,

L

?f (x , y ) ds =?f [?(y ), y +[?'(y ) ]dy

c

d

2

3、空间曲线Γ由参数方程

?x =?(t ) ?

?y =φ(t ) (α≤t ≤β) ?z =ω(t ) ?

给出时,

?f (x , y , z ) ds =?f [?(t ), φ(t ), ω(t )]

Γ

β

?'(t ) 2

+φ'(t ) +ω'(t ) dt

2

2

α

【例1】计算L

22

x +y ds 22

x +y =ax (a >0) , 其中L 为圆周

【解法一】L 可化为参数方程

a a

(x -) 2+y 2=() 2

22

a ?

?x =2(1+cos θ) ?

?y =a sin θ??2

(0≤θ≤2π)

[x ']2+[y ']

2

a ?a ??a ?

=?-sin θ?+?cos θ?=

2 ?2??2?

22

2

a θ22

x +y =(1+cos θ) =a cos

22

22π2πa θθθ222

x +y ds =?cos d θ=?a cos d () L 020222

π

=?a cos t dt =2a ?cos tdt =2a 2

2

2

π

2

1是在x 轴上方的一支, 则方程应为 【解法二】曲线L 关于x 轴对称, 设L

y =ax -x 2

(0≤x ≤a )

而被积函数

f (x , y ) =x 2+y 2在L 上关于x 轴偶对称, 故

L

22

x +y ds

=2?x 2+y 2ds

L 1

=2?L 1

=2?

a

ax ?

?

?a -2x ?

+?dx 2?

?2ax -x ?

2

a

=2?

4(ax -x 2) +(a -2x ) 2

dx 4(ax -x )

a 2

=2?ax dx

4x (a -x ) 0

a

dx

=a a ?

0a -x

a

=a -2[

]0

a

=2a 2

【例2】计算半径为R , 中心角为2α的圆弧L 对于它的对称轴的转动惯量I (设线密度为ρ

=1) 。

解:建立如图所示的坐标系

I =?y 2ds

则

L

?x =R cos θL :?

?y =R sin θ而

于是

(-α≤θ≤α)

I =?R 2sin 2θ(-R sin θ) 2+(R cos θ) 2d θ

-α

α

3

2

3

=?R sin θd θ=2R ?sin 2θd θ

-α

αα

=2R 3?

α1-cos 2θ

2

θ

1?1?=2R 3?θ-sin 2θ?

4?2?0 1

=R (α-sin 2α)

2

3

α

范文二：【doc】基于第二类椭圆积分的子午线弧长公式变换及解算

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

换及解算

第31卷第4期

2011年8月

大地测量与地球动力学

JOURNALOFGE0DESYANDGEODYNAMICS Vo1.3lNo.4

Aug.,2011

文章编号:1671-5942(2011)04-0094-05 基于第二类椭圆积分的子午线弧长公式变换及解算

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

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

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

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

摘要通过推导,将子午线弧长公式变换为基于第二类椭圆积分的两种形式:"形式I"将子午线弧长公式表

达为一个有理函数和第二类椭圆积分之和,建立了以大地纬度B为自变量的子午线弧长公式与第二类椭圆积分之

间的关系;"形式?"给出了以归化纬度为自变量,直接利用第二类椭圆积分计算子午线弧长的公式.利用此两

种形式的子午线弧长公式,在Matlab中编写程序,调用第二类椭圆积分函数EllipticE(,k)计算子午线弧长,精度

和计算效率均优于经典算法.对CGCS2000所采用的地球椭球子午线弧长的计算表明,此两种形式的子午线弧长

公式建立了子午线弧长公式与第二类椭圆积分的关系,结构简洁,易于展开,一定程度上完善了子午线弧长理论,

且便于手工计算及计算机程序实现.

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

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

CALCULATINGMERIDIANARCLENGTHBYTRANSFoRMING

ITSFoRMULAEINToELLIPTICINTEGRALoFSECoNDKIND

GuoJiachun,ZhaoXiuxia,XuLi?

,

TianJinsongandGaoFei

Schoolo厂

Schoolof

Schoolof

AbstractAnewideathatbytransformingthemeridianarclengthformulaintotwootherformswasputfor—

wand.whichareexpressedbytheellipticintegralsofthesecondkind,theywerenamedas"FormI''and"Form

?

".In"FormI",themeridianarelengthformulaisrepresentedasthesumofarationalfunctionandtheelliptic

integralsofthesecondkindbythegeodeticlatitudeparameter,whichestablishsthefunctionrelationsbetweenthe

meridianarclengthandtheellipticintegralsofthesecondkind.Analogously,takingthereducedlatitudeasanin

dependentvariable.the"Form?

"formulagiveasimplerformofthemeridianarelengthformulaintermsoftheel—

lipticintegralsofthesecondkind.Onthesebasesoftheoreticalanalysis,thecomputerprogramiscompiledin

MATLABbycallingtheEllipticE(,k)Functiontocalculatethemeridianarelength.Itwasprovedthatthis

methodimprovedgreatlytheaccuracyandefficiencyofpreviouscalculation.Moreover,themeridianarelengthof

CGCS2000wasalsocalculatedbasedontheprinciplethatprovided.Theresultsindicatethatt

he"FormI"and

收稿日期:2011-02—18

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

课题(KF2010W40—046)

作者简介:过家春,男,1981年生,硕士,讲师,主要研究大地测量学.E—

mail:guojiachun@ahau.edu.cn

通讯作者:赵秀侠,女,1981年生,博士,主要研究方向为地图学与地理信息系统.E—

mail:xiuxia99@126.con 96大地测量与地球动力学

d/生\

dI厂

k2cos2k2sin2~

.

sinc0s

,

.

(1-k2sin(D)一

霉(9)(1一

k2sin)

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

+(?一)J=—

Jfd+J...........................—J—.,...,..................一r?,,, o(1一k2sin).

(1一k2)=

簪

J(二d=E(cp,)(10)

即等式(8)成立.

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

公式变换

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

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

令:/a,其中,0<b<0,代人式(8) 并两边同时乘以a,移项后可得

口(一)—二赢

一

+aE(~p,)(11

~/1一sin(D

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

=一

ae2s

iinBcosB+aE(日,e)(12)

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

S:fMdB:,竺+.E(B,)

0?0COSB+bsinB

(13)

特别地,当B=90.时,

S=.E(,e)(14)

=3盯-一0.其中,归化纬度与大地纬度B之间的 ta:tanB:tanB(15) 等舞:6,

对式(1)作变量代换=arctan(AtanB)得 s=J=dB=JcdB=

easin~).(1一

+

6"

S:f(18)

再由=一,可得

S:f

.f,/二d一.f.~/—1-e2—sin20d0:./0JO 0E(,e)一aE(,e)=

o(詈,e)一.E(詈一/z,e)(19)

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

称为"形式?".

4基于"形式I","形式?"的子午线

弧长计算与分析

以IUGG75椭球参数为例,笔者在Matlab软件 中调用第二类椭圆积分函数EllipticE(,)?'", 分别编写了基于"形式I"和"形式?"的子午线弧 长解算程序,实现了子午线的弧长计算功能.Mat— lab中数据显示精度可以任意设置,本文取至10 m.表1为基于"形式I"和"形式?"的计算结果 与"经典算法"的结果对照.

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

putedwithdifferentalgorithm

大地纬度B

(归化纬度)

IUGG75椭球子午线弧长S(m)

经典箅法本文算法

注:本文算法指基于"彤式I"和"形式?"调用第二类椭圆积分的算

法,二者结果完全一致.

表l显示,基于"形式I"和"形式?"的计算结 果完全一致.由于其结果是在椭圆积分真值的基础 上计算而得的,因此可视为相应纬度的子午线弧长 真值.在Matlab,Mathematica,Maple等数学软件的 特殊函数库,C++的Boost库等,计算第二类椭圆 积分的内部算法为卡尔松方法,计算效率高于按二 项式定理的级数展开方法,图2为在Matlab中分别 利用本文算法,二项式定理展开方法同时计算子午

线弧长的CPU执行时间对比,结果表明,本文算法的计算效率提高到了lO倍左右.

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

直观上,图3中曲线接近于直线,这正反映了地 球扁率很小,平均曲率半径很大的特点. 图2不同子午线弧长计算方法的CPU运行时间对比 Fig.2ComparisonbetweentheCPUtimesfordifferentpro— g

邑

暮

cedurestocomputethemeridianarclength

/

/

//jl

l/l

图3子午线弧长随大地纬度变化的曲线图 Fig.3GraphofthemeridianarclengthchangingwithGeo— deticLatitude

为方便读者应用,下面给出基于"形式I"Mat— lab程序代码:

clear

a=6378140;%椭球参数赋值

b=6356755.2881575287; e=sqrt((a^2一b^2)/a^2);

fori=1:10:

B(i)=(i一1)pi/18;%以10.为步长

S(i)=amfun('EllipticE',sin(B(i)),e)一ae^2 :l:sin(B(i))COS(B(i))/(1一esin(B(i))^2) (1/2);%子午线弧长计算

end

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

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

5结论

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

表2CGCS2000椭球子午线弧长

Tab.2MeridianarclengthoftheCGCS2000ellipsoid

()cGCs2000椭球子午线弧长(m)

10.00O0

f09.5801.7231)

20.00O0"

f19.5617.6474)

30.0000"

f29.5500.2915)

40.0O00

f39.5418.9975") 50.0000"

f49.5418.7985") 60.0000

(59.5459.7878)

70.0000

f69.5617.0746)

80.0000

(79.5801.3492") 90.0000"

(90.00001

l105854.83319845 2212366.25410298 3320113.39784502 4429529.03023659 5540847.04156097 6654072.81936744 7768980.72765552 8885139.87183676 10001965.72923050 度肛,然后再按式(19)进行计算,这对于计算机程序 设计来说是容易实现的,但不如"形式I"简洁.而 且当B=/x=90.时,两种形式的公式直接将子午线 弧长表达为第二类完全椭圆积分形式,这是"经典 算法"所不具备的.

由于许多数学软件,程序语言的函数库中自带

第二类椭圆积分函数,如Matlab,Mathematica,Maple 等数学软件的特殊函数库,c++的Boost库等,在

编写程序时按"形式I"或"形式?"的形式实现,直

接调用其第二类椭圆积分函数即可,代码简短高效.

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

仍然没得到完美解决的问题."形式I"和"形式

?"两种形式的公式可以建立起子午线弧长公式与

其他特殊函数的关系,如超几何函数,雅氏椭圆函数

等"],有望为子午线弧长的反解问题提供新的思

路.尤其是"形式?"式(19)的第一项对于一定的

地球椭球为一常数,将子午线弧长的反解问题直接

转化为了第二类椭圆积分的求逆问题.

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

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

参考文献

1WolfgangTorge.Geodesy(3rd.)[M].Berlin:WalterDe Gruyter,2001.

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

汉大学出版社,2001.(KongXiangyuan,GuoJimingand LiuZongquan.Foundationofgeodesy[M].Wuhan:Wuhan UniversityPress,2001)

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

绘工程,1998,7(4):38—42.(JiangChengguangandYan Guifeng.Anovelmethodforthecalculationofmeridianarc lengthofellipsoid[J].EngineeringofSurveyingMapping, 1998,7(4):38—42)

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

式[J].测绘通报,2001,11:14—15.(NiuZhuoli.Formu—

laeforcalculationofmeridianarclengthbytheparametersof

spacerectangularcoordinates[J].BulletinofSurveyingMap—

ping,2001,11:14—15)

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

2006,5:4—6.(LiuXiushan.Numericalintegralmethodfor calculatingmeridianarclength[J].BulletinofSurveying Mapping,2006,5:4—6)

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

定[J].测绘学院,2000,17(3):167—171.(Yi

Weiyong,BianShaofengandZhuHanquan.Determination offootpointlatitudebyanalyticpositiveserires[J].Journal ofInstituteofSurveyingMapping,2000,17(3):167—

171)

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

[M].北京:国防工业出版社,2004.(BianShaofengand

XuJiangning.Computeralgebrasystemandmathenlaticala—

nalysisingeodesy[M].Beijing:NationalDefenceIndustrial Press,2004)

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

大地测量与地球动力学,2007,(5):59—62.(LiuRenzhao

andWuJicang.Reeursivecomputationofmeridianarclength withdiscretionaryprecisionJ].JournalofGeodesyandGeo—

dynamics,2007,(5):59—62)

9KlausHeh1.C++andJavacodeforrecursionformulasin mathematicalgeodesy[J].GPSSolution,2005,9:51—58.

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

版社,2000.(WangZhuxiandGuoDunren.Introductionto specialfunction[M].Beijing:PeijingUniversityPress, 2000)

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

大学,2009.(LuYuan.AstudyoncontributionofAbelto thetheoryofellipticfunction[D].InnerMongoliaNormal University,2009)

12AbramowitzMandStegunIA.Handbookofmathematical functions[M].NewYork:DoverPublications,1964. 13刘式适,刘式达.特殊函数[M].北京:气象出版社,

1988.(LiuShishiandLiuShida.Specialfunction[M]. Beijing:ChinaMeteorologicalPress,1988) 14TheMathWorksInc.MATLAB7.0(R14SP2).TheMath—

WorksInc.,2005.

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

WGS84的比较[J].测绘,2009,38(6):189—194.

(ChengPengfei,eta1.ParametersofCGCS2000ellipsoid andcomparisonswithGRS80andWGS84[J].ActaGeo—

daeticaetCartographicaSiniea,2009,38(6):189—194)

范文三：基于第二类椭圆积分的子午线弧长公式变换及解算（可编辑）

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

. .

第 卷第 期 大地测量与地球动力 学 年 月 .,

文章编号:

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

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

安徽农业大学理学院,合肥 、

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

合肥工业大学土木与水利工程学院,合肥 /

摘 要 通过推导,将子午线弧长公式变换为基于第二类椭圆积分的两种形式:“形式 ”将子午线弧长公式表

达为一个有理函数和第二类椭圆积分之和,建立了以大地纬度 为自变量的子午线弧长公式与第二类椭圆积分之

间的关系;“形式?”给出了以归化纬度 为自变量、直接利用第二类椭圆积分计算子午线弧长的公式。利用此两

种形式的子午线弧长公式,在 中编写程序,调用第二类椭圆积分函数 , 计算子午线弧长,精度

和计算效率均优于经典算法。对 所采用的地球椭球子午线弧长的计算表明,此两种形式的子午线弧长

公式建立了子午线弧长公式与第二类椭圆积分的关系,结构简洁,易于展开,

一定程度上完善了子午线弧长理论,

且便于手工计算及计算机程序实现。

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

中图分类号: 文献标识码:, , ?

, 厂 ? ., “ ’’ “?”. “ ”,, .,. “?”?. , ,. . , .“”收稿日期:?

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

作者简介:过家春,男,年生,硕士,讲师,主要研究大地测量学. ? :. .

通讯作者:赵秀侠,女,年生,博士,主要研究方向为地图学与地理信息系统. ? : . 大地测量与地球动力学

/生 \厂.

。 一

、

霉一所以对于等式的右边有

盯 一

?一 ? 。 其中,归化纬度 与大地纬度 之间的

...........................? ?.,...,..................一?,,、 一 。

::一

等舞 : ,

簪

对式作变量代换得

二

即等式成立。

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

。 一

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

式。

. 基于第二类椭圆积分的子午线弧长公 再由一 ,可得

式变换

:

令 : / ,其中, ,代人式并两边同时乘以 ,移项后可得

口一 ?二 赢

一 ,

/ 一将式 中的 改写为 , 改写为 ,得子午 。 、/ 二 一。。/? ? :

线弧长的第二类椭圆积分的表达形式为: ./ 一 日, , 一,

式可进一步化简化为:

詈, 一。 詈一/ ::~ 竺 。 ,?

特别地,当 。时, 。 , 第 期 过家春等:基于第一类椭圆积分的子午线弧

长公式变换及解算称为“形式?”。

基于“形式 、“形式?”的子午线

弧长计算与分析

以 椭球参数为例,笔者在软件

中调用第二类椭圆积分函数 , ? ’” ,

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

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

长解算程序,实现了子午线的弧长计算功能。 ?. ?中数据显示精度可以任

意设置,本文取至 。表 为基于“形式 ”和“形式?”的计算结果与“经典算法”

的结果对照。

邑

/

/

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

/ /

/ 椭球子午线弧长

大地纬度归化纬度

经典箅法 本文算法

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

为方便读者应用,下面给出基于“形式 ” ?程序代码: ;%椭球参数赋值 . ;

一 /;: : 一/ ;%以 。为步长 ‘ ’,, 一::/ 一注:本文算法指基于“彤式 ”和“形式?”调用第二类椭圆积分的算

/ ;%子午线弧长计算

法,二者结果完全一致。

表 显示,基于“形式 ”和“形式?”的计算结 该程序中用于子午线弧长计算的代码仅一行,

果完全一致。由于其结果是在椭圆积分真值的基础 且精度没有损失。基于“形式?”的程序结构与其

上计算而得的,因此可视为相应纬度的子午线弧长

类似,计算结果完全一致。

真值。在、 、 等数学软件的 将程序中的椭球参数替换为中国 国家大

特殊函数库、的库等,计算第二类椭圆

地坐标系 所采用的地球椭球参数 ,

积分的内部算法为卡尔松方法,计算效率高于按二 即可计算得到 地球椭球的子午线弧长,

项式定理的级数展开方法,图 为在中分别

结果如表 所示。

利用本文算法、二项式定理展开方法同时计算子午

结论

线弧长的执行时间对比,结果表明,本文算法

的计算效率提高到了 倍左右。 “形式 ”和“形式?”的两种子午线弧长公式

将子午线弧长的计算转化为第二类椭圆积分的计 程序还分别绘制了子午线弧长随大地纬度变化 算,简洁直观,便于手工计算及计算机程序实现。其 的曲线图 图 。

中,按“形式?”进行子午线弧长计算时,需分两步 直观上,图 中曲线接近于直线,这正反映了地 进行,即先按式由大地纬度为 计算出归化纬 球扁率很小、平均曲率半径很大的特点。

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

第 3 1卷 第 4大 地 测 量 与 地 球 测 力 l. 31 No. Vo期4学

2 0 1 1 年 8A ug. , JOURNAL O F GEOD ESY AND GEOD YNAM

( ) 文章编号 : 1671 25942 2011 04

20094 205

3

基于第二编编编编分的子午编弧编公式编编及

))) ))1 2 1 1 3 测 测 徐高测家春 测秀侠 田测松

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

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

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

通测推测 将子午测弧测公式测测测基于第二测测测测分的两测形式 “形式 ?”子午测弧将测公式表摘 要 ,:

,B 达个测一有理函数和第二测测测测分之和 建立了以大地测度 测自测量的子午测弧测公式与第二测测测测分之

μ测的测系 ;“形式 ?”测出了以测化测度 测自测量 、直接利用第二测测测测分测算子午测弧测的公式 。利用此 测形式两的

( )子午测弧测公式 ,在 M a tlab中测程序写 ,测用第二测测测测分函数 E llip tic E x, k 测算子午测弧测 ,精度 和测算效率均测于

测典算法 。测 CGCS2000所采用的地球测球子午测弧测的测算表明 ,此两测形式的子午测弧测 公式建立了子午测弧测公式

与第二测测测测分的测系 ,测测构测 ,易于展测 ,一定程度上完善了子午测弧测理测 , 且便于手工测算及测算机程序测测 。

;;;;编编编 子午测弧测 公式测测 测测测分 大地测度 测化测度

: : 中编分编号 文献编编编

P226A

CAL C UL A T ING M ER ID IA N A RC L ENGTH BY TRA NSFO

RM ING

ITS FO RM UL A E INTO ELL IPT IC INTEGRAL O F SECO ND K

IND

1 ) 2 ) 1 ) 1 ) Guo J iachun, Zhao X iuxia, Xu L i, Tian J in songand Gao 3 )Fe i

) 1 S chool of S cience, A nhu i A g ricu ltu ra l U n iversity, H efei 230036

A new idea tha t by tran sfo rm ing the m e rid ian a rc length fo rm u la in to two o the r fo rm s A b stra c

wa s p u t fo r2twand, wh ich a re exp re ssed by the e llip tic in tegra ls of the second k ind, they we re nam ed a s“Fo rm

?”and“Fo rm

?”. In“Fo rm ?”, the m e rid ian a rc length fo rm u la is rep re sen ted a s the sum of a ra tiona l func tion and the e llip tic in tegra ls of the second k ind by the geode tic la titude p a ram e te r, wh ich e stab lish s the func

tion re la tion s be tween the m e rid ian a rc length and the e llip tic in tegra ls of the second k ind. A na logou sly,

tak ing the reduced la titude a s an in2 dep enden t va riab le, the“Fo rm ?”fo rm u la give a simp le r fo rm of the m e rid ian a rc length fo rm u la in te rm s of the e l2 lip tic in tegra ls of the second k ind. O n the se ba se s

( ) of theo re tica l ana lysis, the comp u te r p rogram is comp iled in MA TLAB by ca lling the E llip tic E x, k Func tion to ca lcu la te the m e rid ian a rc length.It wa s p roved tha t th

is m e thod imp roved grea tly the accu racy and effic iency of p reviou s ca lcu la tion. Mo reove r, the m e rid ian a

3 收稿日期 : 2011 202 218

( ) ( ):40771117 ; KF2010W 40 - 046 基金编目 国家自然科学基金 国家测测信息化工程技测研究中心测放测测

:,, 1981 ,,,. E - m a il: guo jiachun@ ahau. edu. cn作者编介 测家春 男 年生 测士 测测 主要研究大地测量学

:,, 1981 ,,. E - m a il: xiuxia99 @ 126. 通编作者 测秀侠 女 年生 博士 主要研究方向测地测学与地理信息系测

com

“Fo rm ?”fo rm u la a re simp le r and mo re su itab le fo r se rie s exp an sion s and the rea liza tion on comp

u te r than the

c la ssica l fo rm. Fu rthe rmo re, it p e rfec ts the m e rid ian theo ry.

Key word s:m e rid ian a rc length; fo rm u la tran sfo rm a tion; e llip tic in tegra ls; geode tic la titude; reduced la

φ22φ() ( =1 - sin E , k k 5 0φφ1 引言d?)φ, 作测量代测 x = sin, 另外一测测法与此等价 子午测弧测是大地测量学中的一个基本量 。测

测算从赤道测 始到 任 意大 地测 度 B 的子 午测 S2 2k弧测1 - kt( ) ( E x, k =d t6 2?02 BB的基本公式测1 - t)( )a 1 - e( dB1 3 02 φ?0式中 , k 测测测模 且 0 < k="">< 1="" ;="" 称测幅度="" 。,s="M" db=")" ()1="" -="" esinb="" )="" (="" ),="" 5="" 6="" 通常="" 。="" 二测不完全测测测分="" 、测第称式式中="" ,="" a="" 测地球测球测半测="" ;="" e测地球测球第一偏心π="" ;="" m="" ,="" b="" 率="" 子午圈曲率半径="" 大地测度="" 。测测(="" φ="" φ,="," x="sin=" 1="" ,="" 5="" 特测地="" 式="" 当测称2测然="" ,="" 子午测弧测的测算涉及到勒测德第二测="" 测测测分="" )="" 、π="" (="" )(="" )="" (6="" ,="" ee="" k="" e="" ,第二测完全测测测分="" 测="" 测测、、()="" 测测称第二测测测测分="" ,="" 其原函无数法2用初等函的形数式表达="" ,="" 不能直接求出="" 。目前="" ,="" )(="" )="" k="" 或="" e="" 1,="" k="" 2="" (="" )世界="" 各国测子午测弧测的测典测算方法是将="" a="" 1="" -="" e。第二测测测测分的何意几测测即测测的弧测="" 。2="" -="" 2="" (="" )="" 1="" -="" esinb="" 定理展测测测形数式="" ,="" 后然按二测式测一测测的方参数程测[="" ,="" ]="" 12="" 再逐测测分="" ,="" :精度的测算公式="" 得出一定x="" a="(" (0="">< b="">< a="" 7="" b="" cd0="" 0="" 0="" θsin2="" ))(="" )="" +s="a" 1="" -="" e{="" a="" b="" -="" sin2b="" +="" sin4b="" -="" 02462="" 2y="bco" (="" )="" a令="" k="-" b/="" a,="" 易测明="" ,="" 5="" 表达="" 容式所bsin6(="" )}2="" 的φ()(="" )="" 第二测测测测分="" e="" ,="" k="" 的几何意测测="" :以式="" 7="" 测式中="" ,="" a、b="" 、c、d、="" 测一测e,="" 测的常系数="" 有具0="" 0="" 0="" 0="" 参数方程="" ,="" 当="" a="1," 0="">< b="">< 1="" 测表示的测测“自其短="" 与体达="" 表式可测参文献="" [="" 2="" ]。()="" 测测点到测测上任一点="" p="" 的弧测="" ”测="" 1="" 。可上,="" [="" b="" ,="" b="" ]在此基测上="" 以测算任意区测="" 的1="" 2="">

() 子午测弧测 本文称之测“测典算法 ”。 测 了 使

子 午 测 弧 测 公 式 的 表 达 式 更 测 测 测 ,

或 达 到 精 度 更 高 、便 于 测 算 机 测 测 等

, 目 的 少不[ 3 - 9 ] 学 者 测 此 展 测 了 研 究 。但 从 其 研 究

成 果 看来 , 大 都 仍 停 留 在 以 研 究 子 午

测 弧 测 的 测 算 测 目 的 的 表 面 测 测 上 , 尚 未

探 求 出 子 午 测 弧 测 公 式 第与 二 测 测 测 测

分 测 准 形 式 之 测 的 内 在 本 测 测 系 , 测 在

一 定 程 度 上 影 响 了 子 午 测 弧 测 公 式 理

测 上 的 完 整 性 , 也 限 制 了 测 测 测 分 研 究 测 1 子午圈成 果 在 子 午 测 弧 测 测 算 测 测 上 的 测 用 。F ig. 1 M e rid ian2 第二测测测测分及其展测式

第二编编编编分的基本编系式及2. 2 2. 1 第二编编编编分的编准形式其编[ ]10 - 13 一般测测测分定测测 明

( R x, y 2( 测测测程中用到的第二测测测分的测系式测 :3 φk sinco ?φ(=+E , k )2 2φsφ1 - ksin)) dxφ φ d2 ( 4 ( 8 ( )1 - k2 23 / [ , ] ?( )1011 0其中 , R x, y x y 理函数 , 测和的有而( φ) 1 - ksin))测测测分即 求 测 测 的 弧 测 勒 测 德 测 。2测明如下 :明

因测了一般测测测分能测化成 3 测基本测型 。其中 , 第 二测

测测测分的测准形式测

2φ( )( ) k sinco 12 和式 13 将子午测弧测公式表测达一个有式 d=2 2φsdkφ理函数和第二测测测测分之和 , 建立了以大地测度 测自测1 - ksin

2 22 242 2φφφ k sin 量的子午测弧测公式与第二测测测测分之测 的测系 , 测测称测 k c o s k sin co +=- 23 / 222 22 φs( ) 1 - k sinφ“形式 ?”。φ 1 - ksin1 - k2φsin23. 2 基于第二编编编编分的子午编弧编k( 9 2 24 / 42 23 2 φ) (kk- 2 1 - sinksinφ+ sinφ)公 式编编 ?

( )8 所以测于等式 的右测有1, P B , 如测 所示 点的大地测度测 测化测度测φ 2 φφksinco d2 π=( )+ 1 - k223 / 22 2?μ μθ0, B =- 。其中 测度 与大地测度 的测化之测φ( s1 - k sin φ1 - ksin2φ)φ 2 2 24 φ k- 2 ksin+ k测系测φ d42 23 / ?0φ( φ) sin1 - ksin+ 22μ ( tan15 φe B B1 - tan=tanbφd2 a)==( )1 - k2 23 / ?0( φ) 1 - ksin因此有2φ 2 24 2φ 1 - 2 ksin+ k2aφ d42 2( )a 1 +tanμ0? 2 φ( ) sin1 - ksin2φ( =a 1 + tanB bμ ( d16 μB d=d=2 2φ2μ3 / 2( φ) ( ) ( 1 - ksinb 1 + tanb 1 + 2)22 23 / 2 μ)tan?0( ) 1 - ksinφb ( )μ ( )1 = a rc tan tanB 测式 φ作测量代测得φ d=22aφφφ (1 - sin d= k E , ( 10 ?02 BB))( a 1 - eS =M dB =dB =2 2 )k ( )即等式 8 成立 。0?0?( ) 1 - esinB

23 / 22 ( a 1 +tan 3 基于第二测测测测分的子午测弧测2 μ2( )a 1 - eμ)bμd22 22?0公式测测( b 1 + tan sine a μ3 / 2( 1 -22 2 2) μ)μ a sin + b co 由于旋测测球上的子午圈测测测 , 所以子午测 弧测的μs( 17 测算测测也测即测弧测测测 , 测就决定了子 午测弧测公式必然)化测后测可以测测测 第二 测测 测 测分 形 式 。μ22 2 2μ ( asin + b co S =18 3. 1 基于第二编编编编分的子午编弧编0?)μμsdπ 公θμ=- , 再由 可2式编编 ?μ得μ2 2 2 22 = sin+ co S ab( a令 k - b/ a, 中 , 0 < b="">< a="" ,="" 入式="" 其代?0μ=")8" θ并两测同测乘以="" a,="" 移测后可得μ2sd2="" 2="" 2="" 2θa="" co="" s+="" b="" sin="" -φφπd2="" θθ="(" )d="a" 1="" -="" k2="" 2="" 2?π="" 02/="" (="" φ1="" -="" k1="" -="" ksin2="" 2="" 2="" 2?θa="" co="" s+="" b="" sin="" 22?0φ)φφsin="" ak="" sinco="" sθθd-(="" -11="" ae="" +="" θθ2="" 2φ1="" -="" ksin2="" 22="" )co="" s+="" ab?0()φ,="" k="" θφ(="" )将式="" 11="" 中的="" k="" 改写测="" e,="" 改写测="" b="" ,="" 得子θπ="" 2θ2sind="2222午=a1" -="" e="" sin="" d="" -="" a="" 1="" -="" e="" sin="" dθθθθ0?0?2b:测弧测的第二测测测测分的表形式测达="" sinco="" ae="" b="" π()(="" +="" ae="" b="" ,="" e12="" 0sθb()="" ()="" ae="" ,="" e-="" ae="" ,="" e2="" 2="" s="M" db="-2)esin1" -="">

ππ( )式 12 可测一步化测化测 :μ() ( (aE , e- aE - , 19 2 B22( sinB co sB )- b))e()+ aE B , e202 22 2( )a式 19 也可由第二测测测测分的几何意测直接得S =M dB = -a co s B + b sin ?B到 。 (13

)μ ( )( , , 特测地 = 90 测 19 14 当?式化测式特测地 , 当 B = 90 ?测

,) 。π( )S= aE , e( 14 ? 90 ( )μ式 2 19 测出了以测化测度 测自测量的直接 利)

用第二测测测测分测算子午测弧测的公式 , 测测

称测“形式 ?”。

4 基于“形式 ?”、“形式 ?”的子午测

弧测测算与分析

以 IU GG75 测球参数测例 , 笔者在 M a tlab 测件[ 14 , 15 ] ( 中测用第二测测测测分函 数 E llip tic E x,) k ,分测测写了基于“形式 ?”和“形式 ?”的子午测

测 2不同子午测弧测测算方法的 CPU 运行测测测比弧

Comp a rison be tween the CPU tim e s fo r d iffe ren F ig. , M a t测解算程序 2测了子午测的弧测测算功能 。测- 8lab中数据测示精度可以任意测置 , 本文取至2t p ro210 m。表 1测基于“形式 ?”和“形式 ?”的测算测cedu re s to comp u te the m e rid ian a rc length果

与“测典算法 ”的测果测照 。

表 1 基于不同编算方法的子午编弧编编算

编果

Ta b. 1 C om pa r ison am on g the m er id ian a rc len g th () I U G G 7 5 测球子午测弧测 S m 大地测度 B3()μ 测化测度 测典算法本文算法

10?00′01 105 855. 347 91 105 855. 347 887 0″测 3 子午测弧测随大地测度测化的曲测测51( 095801. 723 ?′F ig. 3 Grap h of the m e rid ian a rc length changing w ith 2 212 367. 284 32 212 367. 284 275 )1″Geo21320?00′0de tic L a titude0″3 320 114. 945 03 320 114. 944 995 ( 测方便测者测用 ,下面测出基于“形式 ?”M a t2195617. 647 93?′)4″4 429 531. 096 44 429 531. 096 388 lab: 程序代测30?00′041c lea 0″5 540 849. 629 05 540 849. 629 023 ( 295500. 291 ?′r09)5″a = 6378140; % 球测测参数测6 654 075. 930 56 654 075. 930 459 40?00′0820″b = 6356755. 2881575287;

( 395418. 997 ?′7 768 984. 364 47 768 984. 364 427 ( ( ) e = sq rt a^2 - b^2 / )506″) a^2 ;50?00′08 885 144. 035 88 885 144. 035 813 0″fo r i = 1: 56( 495418. 798 ?′10;10 001 970. 421 110 001 970. 421 226 )5″40( ) ( ) B i= i - 1 3p i /18; % 10 测步以?60?00′0注 :本文算法指基于“形式 ?”和“形式 ?”测用第二测测测测分的测算

法 ,二者测果完全一致 。( ) ((( ) ) ) S i= am fun E llip ticE, sin B i, e- a33 ‘’

表 1测示 ,基于“形式 ?”和“形式 ?”的测e^2算测(( ) ) (( ) ) ( (( ) ) 3 sin B i3 co s B i/ 1 - e^2 3 sin B i果完全一致 。由于其测果是在测测测分真测的基测 上测) ( ) ^2 ^ 1 /2 ; % 子午测弧测测算算而得的 ,因此可测测相测测度的子午测弧测end真测 。在 M a tlabM a them a ticaM ap le 测、、等数学,测程序中用于子午测弧测测算的代测测一行 件的且精度没有测失 。基于“形式 ?”的程序测构与特殊函数测 、C + +的 Boo st测等 , 测算第二测测测 测其 测似 ,测算测果完全一致 。分的部算内法测测卡松方法 ,测算效率高于按二 测式定

5 测测理的测数展测方法 ,测 2 测在 M a tlab 中分测

利用本文算法 、二测式定理展测方法同测测算子午 测“形式 ?”和“形式 ?”的两测子午测弧测公

式 子午将测弧测的测算测化测第 二测 测 测测 分的 测 算 弧测的 CPU 测行测测测比 ,测果表明 ,本文算法 的测算

,测测直测 ,便于手工测算及测算机程序测测 。其 中 ,按效率提高到了 10倍左右 。“形式 ?”测行子午测弧测测算测 , 需分两步 测行 ,程序测分测测制了子午测弧测大随地测度测化 的曲测测 ( )即先按式 15 由大地测度测 B测算出测化测() 测 3 。

直测上 ,测 3中曲测接近于直测 ,测正反映了地 球

( 测大 学 出 版 社 , 2001. Kong X iangyuan, Guo J im

ing and表 2 C GC S2000编球子午编弧

编L iu Zongquan. Founda tion of geode sy [ M ]. W uhan: W uhanTa b. 2 M er id ian a rc len g th of the C GC S2000 e ll 3大地测度 BU n ive rsity P re ss, )(CGCS2000测球子午测弧测 m)()2001μ 测化测度

姜晨光 ,测桂峰 . 测球子午测弧测测算的新方法 [ J ]. 测 测工10 ?00 1 105 854. 833 198 45( ( ) ′00″程 , 1998, 7 4 : 38 - 42. J iang Chengguang and Yan ( 095801. 7231 ?′Gu ifeng. A nove l m e thod fo r the ca lcu la tion of m e rid 42 212 366. 254 102 98)″ian a rc length of e llip so id [ J ]. Enginee ring of Su rveying

20 ?00 M app ing,′00″3 320 113. 397 845 02( ) 1998 , 7 4: 38 - ( 195617. 6474 ?′)42

)″牛卓立 . 以空测直角坐测系测的子参数午测弧测测算公 式 [ J 4 429 529. 030 236 5930 ?00 5( ]. 测测通测 , 2001 , 11: 14 - 15. N iu Zhuo li. Fo rm u2 lae ′00″fo r ca lcu la tion of m e rid ian a rc length by the p a ram e te rs 5 540 847. 041 560 97( 295500. 2915 ?′of sp ace rec tangu la r coo rd ina te s[ J ]. B u lle tin of Su )″rveying M ap 26 654 072. 819 367 4440 ?00 6p ing, 2001, 11: 14 - 15 ′00″)( 395418. 9975 ?′7 768 980. 727 655 52. [ J ]. ,刘修善 测算子午测弧测的数测测分法 测测通测 )″(2006 , 5: 4 - 6. L iu X iu shan. N um e rica l in tegra l m e 50 ?00 8 885 139. 871 8367 6thod fo r ca lcu la ting m e rid ian a rc length [ J ]. B u lle tin ′00″)of Su rveying M app ing, 2006, 5: 4 - 6( 495418. 7985 ?′易测勇 ,测少峰 ,朱测泉 . 子午测弧测的解析型测测数确10 001 965. 729 230 507)″( ) ( [ J ]. , 2000, 17 3 : 167 - 171. 定 测 测 学 院 学 测 60 ?00 Yiμ( )度 , 然后再按式 19 测行测算 , 测测于测算机程序W e iyong, B ian Shaofeng and Zhu H anquan. D e te rm 测测来测是容易测测的 , 但不如“形式 ?”测测 。而 且ina tion of foo t po in t la titude by ana lytic po sitive se rire s[ μ 当 B == 90 ?测 ,两测形式的公式直接将子午测 J ]. Jou rna l of In stitu te of Su rveying M app ing, 2000 , 8弧测表达测第二测完全测测测分形式 , 测是“测典 算法 ( ) 17 3 : 167 -”所不具测的 。)171由于测多数学测件 、程序测言的函数测中自测 第二测测少峰 ,测江宁 . 测算机代系测数与大地测量分数学析 [M ]. 测测测分函数 ,如 M a tlab、M a them a tica、M ap le 等(北 京 : 国 防 工 测 出 版 社 , 2004. B ian Shaofeng 数学测件的特殊函数测 、C + +的 Boo st测等 , 在 测and Xu J iangn ing. Comp u te r a lgeb ra system and m a 9写程序测按“形式 ?”或“形式 ?”的形式测测 ,直 them a tica l a2 na lysis in geode sy[M ]. B e ijing: N a tiona l 接测用其第二测测测测分函数即可 ,代测测短高效 。)D efence Indu stria l P re ss, 2004 10已有学者指出 ,子午测弧测的反解测测是目前 仍刘仁测 ,伍吉测 . 任意精度的子午测弧测测测测算 [ J ]. 大地测量[ 6 ] 然没得到完美解的测测决 。“形式 ?”和“形式 ( ) (与地球测力学 , 2007, 5 : 59 - 62. L iu R enzhao and W u ?”测两形式的公式可以建立起子午测弧测公式 其与J icang. R ecu rsive comp u ta tion of m e rid ian a rc length

w ith d isc re tiona ry p rec ision [ J ]. Jou rna l of Geode sy 他特殊函的数测系 ,如超几数何函 、雅氏测测函数11and Geo2[ 10 - 13 ] 等 ,有望测子午测弧测的反解测测提供新的思 路 ( ) 5: 59 - dynam ic s, 2007 , ( ) 尤其是“形式 ?”式 19 的第一测测于一定的 )62

Klau s H eh l. C + + and J ava code fo r recu rsion fo rm u 地球测球测一常数 ,将子午测弧测的反解测测直接 测化测12la s in m a them a tica l geode sy[ J ]. GPS So lu tion, 2005, 了第二测测测测分的求逆测测 。9: 51 - 58. 王竹溪 ,郭敦仁 . 特殊函数概测 [M ]. 北京 : 13测于基于“形式 ?”及“形式 ?”的子午测弧测(北京大学出 版社 , 2000. W ang Zhuxi and Guo D un ren. 反解测测需要测一步究研 。In troduc tion to sp ec ia l func tion [ M ]. B e ijing: Pe ijing

U n ive rsity P re ss,14)2000

测源 . 特阿测测测测测函测的测数献 [ D ]. 内蒙古测范 大学 , 15(参考文献2009. L u Yuan. A study on con tribu tion of A be l to the

theo ry of e llip tic func tion [ D ]. Inne r Mongo lia No rm a ( [ M ].B e rlin: W a lte r 1 W o lfgang To rge. Geode sy 3 )l U n ive rsity, 2009D e) Gruyterd. r, 2001.A b ramow itz M and Stegun I A. H andbook of m a them 2 孔祥元 ,郭测明 ,刘宗泉 . 大地测量基学测 [M ]. 武测 : 武a tica l func tion s[M ]. N ew Yo rk: Dove r Pub lica tion s,

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

第 31卷 第 4期 2011年 8月

大 地 测 量 与 地 球 动 力 学

J OURNAL OF GEODESY AND GEODYNAM I CS

V o. l 31N o . 4

A ug . , 2011

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

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

*

过家春 1)

赵秀侠 2)

徐 丽 1)

田劲松 1)

高 飞

3)

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

2) 安徽大学资源与环境工程学院 , 合肥 2300393) 合肥工业大学土木与水利工程学院 , 合肥 摘 要 通过推导 , 将子午线弧长公式变换为基 于第二 类椭圆 积分的 两种形 式 : 形式 将子午 线弧长 公式表

达为一个有理函数和第 二类椭圆积分之和 , 建立了以大地纬度 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 Jiachun 1)

, Zhao X i u x ia 2)

, Xu L i 1)

, T ian Ji n song

1)

and Gao Fe i

3)

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 2300393)School of C ivil and H ydraulic Engineeri n g, H efei Un i v ersit y of Technology , H efei 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); 国家农业信息化工程技术研究中心开放课题 (KF2010W 40-046) 作者简介 :过家春 , 男 , 1981年生 , 硕士 , 讲师 , 主要研究大地测量学 . E-m ai:l guojiachun @ah au . edu . cn :女 , , 博士 , . E -m ai:l 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 = 0B M d B =0B 2

(1-e 2si n B ) 32

d B (1)

式中 , a 为地球椭球长半轴 ; e 为地球椭球第一偏心

率 ; M 为子午圈曲率半径 , B 为大地纬度。

显然 , 子午线弧长的计算涉及到勒让德第二类

椭圆积分 (简称为第二类椭圆积分 ), 其原函数无法

用初等函数的形式表达 , 不能直接求出。目前 , 世界

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

(1-e 2si n B ) -32按二项式定理展开为级数形式 , 然后

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

S =a(1-e 2) {A0B - B 0

2

si n 2B +

C 0

4

si n 4B -

D 0 6 si n 6B +

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

在此基础上 , 可 以计算任意区间 [B1, 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 sin d (5) 与此等价 , 作变量代换 x =sin , 另外一种记法 为

E (x,k ) =

k k t

1-t

d t (6) 式中 , k 为椭圆模 , 且 0<1; 称为幅度。="" 通常="" ,="" 称式="" (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

图 1 子午圈

F i g . 1 M er i d i an

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

讨论过程中用到的第二椭圆积分的关系式为 : E ( , k ) =

2

1-k si n

+

(1-k 2)

(1-k si n )

(8) 证明如下 :

因为

95

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

d d k

2

1-k sin

=

22

1-k si n

-

22

1-k si n

+

42 cos 2 (1-k 2sin 2 ) 3/2 =

k 2-2k 2sin 2 +k 4si n 4

(1-k si n )

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

2

1-k si n +(1-k 2)

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

0 222 +k 44 (1-k si n )

d +

(1-k 2) 0 d (1-k 2sin 2 ) 3/2 =

0 2244

(1-k si n )

d =

0 22 ) 2 (1-k 2si n 2 ) 3/2 d =

1-k si n d =E ( , k ) (10)

即等式 (8) 成立。

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

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

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

令 k =a -b /a,其中 , 0

a(1-k 2) 0

(1-k 2si n 2 ) sin =

-

2

1-k sin

+aE ( , k ) (11) 将式 (11) 中的 k 改写为 e , 改写为 B, 得子午

线弧长的第二类椭圆积分的表达形式为 :

S = 0B M d B =-21-e si n B +aE (B,e) (12) 式 (12) 可进一步化简化为 :

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

+aE (B,e) (13) 特别地 , 当 B =90 时 ,

S 90 =aE (

, e) (14)

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

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

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

=

2

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

tan =1-e tan B =

a

tan B (15) 因此有

d B =

2

b(1+tan 2 )

d =

a (1+

2

2

tan 2 ) b (1+tan 2 )

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

b

a

tan B ) 得

S = 0B M d B =0B 2

(1-e 2sin 2B ) 3/2

d B =

a (1-e 2)

(1-

222

a si n +b cos

) 3/2

a (1+

2

b

tan 2 ) b (1+tan 2 )

d

(17) 化简后为

S =

a si n cos d (18)

再由 =

2

- , 可得

S =

a sin +b cos d

- /2 a cos si n d =

0 /2a cos sin d -

a cos si n d =

a 0 /21-e si n d -a 0 1-e si n 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. 03588885144. 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 B 7. 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

转载请注明出处范文大全网 » 对弧长的曲线积分