過(guò)家春 趙秀俠 徐 麗 田勁松 高 飛

(1)安徽農(nóng)業(yè)大學(xué)理學(xué)院,合肥 230036 2)安徽大學(xué)資源與環(huán)境工程學(xué)院,合肥 230039 3)合肥工業(yè)大學(xué)土木與水利工程學(xué)院,合肥 230009)

基于第二類(lèi)橢圓積分的子午線弧長(zhǎng)公式變換及解算＊

過(guò)家春1)趙秀俠2)徐 麗1)田勁松1)高 飛3)

(1)安徽農(nóng)業(yè)大學(xué)理學(xué)院,合肥 230036 2)安徽大學(xué)資源與環(huán)境工程學(xué)院,合肥 230039 3)合肥工業(yè)大學(xué)土木與水利工程學(xué)院,合肥 230009)

通過(guò)推導(dǎo),將子午線弧長(zhǎng)公式變換為基于第二類(lèi)橢圓積分的兩種形式：“形式Ⅰ”將子午線弧長(zhǎng)公式表達(dá)為一個(gè)有理函數(shù)和第二類(lèi)橢圓積分之和,建立了以大地緯度B為自變量的子午線弧長(zhǎng)公式與第二類(lèi)橢圓積分之間的關(guān)系;“形式Ⅱ”給出了以歸化緯度μ為自變量、直接利用第二類(lèi)橢圓積分計(jì)算子午線弧長(zhǎng)的公式。利用此兩種形式的子午線弧長(zhǎng)公式,在Matlab中編寫(xiě)程序,調(diào)用第二類(lèi)橢圓積分函數(shù) EllipticE(x,k)計(jì)算子午線弧長(zhǎng),精度和計(jì)算效率均優(yōu)于經(jīng)典算法。對(duì) CGCS2000所采用的地球橢球子午線弧長(zhǎng)的計(jì)算表明,此兩種形式的子午線弧長(zhǎng)公式建立了子午線弧長(zhǎng)公式與第二類(lèi)橢圓積分的關(guān)系,結(jié)構(gòu)簡(jiǎn)潔,易于展開(kāi),一定程度上完善了子午線弧長(zhǎng)理論,且便于手工計(jì)算及計(jì)算機(jī)程序?qū)崿F(xiàn)。

子午線弧長(zhǎng);公式變換;橢圓積分;大地緯度;歸化緯度

1 引言

子午線弧長(zhǎng)是大地測(cè)量學(xué)中的一個(gè)基本量。計(jì)算從赤道開(kāi)始到任意大地緯度 B的子午線弧長(zhǎng) S的基本公式為

式中,a為地球橢球長(zhǎng)半軸;e為地球橢球第一偏心率;M為子午圈曲率半徑,B為大地緯度。

顯然,子午線弧長(zhǎng)的計(jì)算涉及到勒讓德第二類(lèi)橢圓積分(簡(jiǎn)稱(chēng)為第二類(lèi)橢圓積分),其原函數(shù)無(wú)法用初等函數(shù)的形式表達(dá),不能直接求出。目前,世界各國(guó)對(duì)子午線弧長(zhǎng)的經(jīng)典計(jì)算方法是將 a(1-e2) (1-e2sinB)按二項(xiàng)式定理展開(kāi)為級(jí)數(shù)形式,然后再逐項(xiàng)積分,得出一定精度的計(jì)算公式[1,2]：

式中,A0、B0、C0、D0、…為一組與 e有關(guān)的常系數(shù),具體表達(dá)式可參閱文獻(xiàn)[2]。

在此基礎(chǔ)上,可以計(jì)算任意區(qū)間 [B1,B2]上的子午線弧長(zhǎng)(本文稱(chēng)之為“經(jīng)典算法”)。

為了使子午線弧長(zhǎng)公式的表達(dá)式更為簡(jiǎn)潔,或達(dá)到精度更高、便于計(jì)算機(jī)實(shí)現(xiàn)等目的,不少學(xué)者對(duì)此展開(kāi)了研究[3-9]。但從其研究成果來(lái)看,大都仍停留在以研究子午線弧長(zhǎng)的計(jì)算為目的的表面問(wèn)題上,尚未探求出子午線弧長(zhǎng)公式與第二類(lèi)橢圓積分標(biāo)準(zhǔn)形式之間的內(nèi)在本質(zhì)關(guān)系,這在一定程度上影響了子午線弧長(zhǎng)公式理論上的完整性,也限制了橢圓積分研究成果在子午線弧長(zhǎng)計(jì)算問(wèn)題上的應(yīng)用。

2 第二類(lèi)橢圓積分及其展開(kāi)式

2.1 第二類(lèi)橢圓積分的標(biāo)準(zhǔn)形式

一般橢圓積分定義為[10-13]

其中,R(x,y)為 x和 y的有理函數(shù),而

橢圓積分即求橢圓的弧長(zhǎng)[10,11]。勒讓德證明了一般橢圓積分能夠化成 3種基本類(lèi)型。其中,第二類(lèi)橢圓積分的標(biāo)準(zhǔn)形式為

與此等價(jià),作變量代換 x=sinφ,另外一種記法為

式中,k為橢圓模,且 0＜k＜1;φ稱(chēng)為幅度。

通常,稱(chēng)式(5)、(6)為第二類(lèi)不完全橢圓積分。

第二類(lèi)橢圓積分的幾何意義即為橢圓的弧長(zhǎng)。設(shè)一橢圓的參數(shù)方程為

圖1 子午圈Fig.1 Meridian

2.2 第二類(lèi)橢圓積分的基本關(guān)系式及其證明

討論過(guò)程中用到的第二橢圓積分的關(guān)系式為：

證明如下：

因?yàn)?/p>

所以對(duì)于等式(8)的右邊有

即等式(8)成立。

3 基于第二類(lèi)橢圓積分的子午線弧長(zhǎng)公式變換

由于旋轉(zhuǎn)橢球上的子午圈為橢圓,所以子午線弧長(zhǎng)的計(jì)算問(wèn)題也即橢圓弧長(zhǎng)問(wèn)題,這就決定了子午線弧長(zhǎng)公式必然可以變換為第二類(lèi)橢圓積分形式。

3.1 基于第二類(lèi)橢圓積分的子午線弧長(zhǎng)公式變換Ⅰ

將式(11)中的 k改寫(xiě)為 e,φ改寫(xiě)為 B,得子午線弧長(zhǎng)的第二類(lèi)橢圓積分的表達(dá)形式為：

式(12)可進(jìn)一步化簡(jiǎn)化為：

特別地,當(dāng)B=90°時(shí),

式 (12)和式 (13)將子午線弧長(zhǎng)公式表達(dá)為一個(gè)有理函數(shù)和第二類(lèi)橢圓積分之和,建立了以大地緯度為自變量的子午線弧長(zhǎng)公式與第二類(lèi)橢圓積分之間的關(guān)系,現(xiàn)簡(jiǎn)稱(chēng)為“形式Ⅰ”。

3.2 基于第二類(lèi)橢圓積分的子午線弧長(zhǎng)公式變換Ⅱ

因此有

化簡(jiǎn)后為

式(19)也可由第二類(lèi)橢圓積分的幾何意義直接得到。

特別地,當(dāng)μ=90°時(shí),式(19)化為式(14)。

式(19)給出了以歸化緯度μ為自變量的直接利用第二類(lèi)橢圓積分計(jì)算子午線弧長(zhǎng)的公式,現(xiàn)簡(jiǎn)稱(chēng)為“形式Ⅱ”。

4 基于“形式Ⅰ”、“形式Ⅱ”的子午線弧長(zhǎng)計(jì)算與分析

以 I UGG75橢球參數(shù)為例,筆者在Matlab軟件中調(diào)用第二類(lèi)橢圓積分函數(shù) EllipticE(x,k)[14,15],分別編寫(xiě)了基于“形式Ⅰ”和“形式Ⅱ”的子午線弧長(zhǎng)解算程序,實(shí)現(xiàn)了子午線的弧長(zhǎng)計(jì)算功能。Matlab中數(shù)據(jù)顯示精度可以任意設(shè)置,本文取至 10-8m。表 1為基于“形式Ⅰ”和“形式Ⅱ”的計(jì)算結(jié)果與“經(jīng)典算法”的結(jié)果對(duì)照。

表1 基于不同計(jì)算方法的子午線弧長(zhǎng)計(jì)算結(jié)果Tab.1 Comparison among the meridian arc lengths computed with different algorithm

表1顯示,基于“形式Ⅰ”和“形式Ⅱ”的計(jì)算結(jié)果完全一致。由于其結(jié)果是在橢圓積分真值的基礎(chǔ)上計(jì)算而得的,因此可視為相應(yīng)緯度的子午線弧長(zhǎng)真值。在Matlab、Mathematica、Maple等數(shù)學(xué)軟件的特殊函數(shù)庫(kù)、C++的 Boost庫(kù)等,計(jì)算第二類(lèi)橢圓積分的內(nèi)部算法為卡爾松方法,計(jì)算效率高于按二項(xiàng)式定理的級(jí)數(shù)展開(kāi)方法,圖 2為在Matlab中分別利用本文算法、二項(xiàng)式定理展開(kāi)方法同時(shí)計(jì)算子午線弧長(zhǎng)的 CPU執(zhí)行時(shí)間對(duì)比,結(jié)果表明,本文算法的計(jì)算效率提高到了 10倍左右。

程序還分別繪制了子午線弧長(zhǎng)隨大地緯度變化的曲線圖(圖 3)。

直觀上,圖 3中曲線接近于直線,這正反映了地球扁率很小、平均曲率半徑很大的特點(diǎn)。

圖2 不同子午線弧長(zhǎng)計(jì)算方法的 CPU運(yùn)行時(shí)間對(duì)比Fig.2 Comparison between the CPU ti mes for differentprocedures to compute the meridian arc length

圖3 子午線弧長(zhǎng)隨大地緯度變化的曲線圖Fig.3 Graph of the meridian arc length changingwith Geodetic Latitude

為方便讀者應(yīng)用,下面給出基于“形式Ⅰ”Matlab程序代碼：

該程序中用于子午線弧長(zhǎng)計(jì)算的代碼僅一行,且精度沒(méi)有損失。基于“形式Ⅱ”的程序結(jié)構(gòu)與其類(lèi)似,計(jì)算結(jié)果完全一致。

將程序中的橢球參數(shù)替換為中國(guó) 2000國(guó)家大地坐標(biāo)系 (CGCS2000)所采用的地球橢球參數(shù)[15],即可計(jì)算得到 CGCS2000地球橢球的子午線弧長(zhǎng),結(jié)果如表2所示。

5 結(jié)論

“形式Ⅰ”和“形式Ⅱ”的兩種子午線弧長(zhǎng)公式將子午線弧長(zhǎng)的計(jì)算轉(zhuǎn)化為第二類(lèi)橢圓積分的計(jì)算,簡(jiǎn)潔直觀,便于手工計(jì)算及計(jì)算機(jī)程序?qū)崿F(xiàn)。其中,按“形式Ⅱ”進(jìn)行子午線弧長(zhǎng)計(jì)算時(shí),需分兩步進(jìn)行,即先按式(15)由大地緯度為B計(jì)算出歸化緯度μ,然后再按式(19)進(jìn)行計(jì)算,這對(duì)于計(jì)算機(jī)程序設(shè)計(jì)來(lái)說(shuō)是容易實(shí)現(xiàn)的,但不如“形式Ⅰ”簡(jiǎn)潔。而且當(dāng)B=μ=90°時(shí),兩種形式的公式直接將子午線弧長(zhǎng)表達(dá)為第二類(lèi)完全橢圓積分形式,這是“經(jīng)典算法”所不具備的。

表2 CGCS2000橢球子午線弧長(zhǎng)Tab.2 M eridian arc length of the CGCS2000 ellipsoid

由于許多數(shù)學(xué)軟件、程序語(yǔ)言的函數(shù)庫(kù)中自帶第二類(lèi)橢圓積分函數(shù),如Matlab、Mathematica、Maple等數(shù)學(xué)軟件的特殊函數(shù)庫(kù)、C++的 Boost庫(kù)等,在編寫(xiě)程序時(shí)按“形式Ⅰ”或“形式Ⅱ”的形式實(shí)現(xiàn),直接調(diào)用其第二類(lèi)橢圓積分函數(shù)即可,代碼簡(jiǎn)短高效。

已有學(xué)者指出,子午線弧長(zhǎng)的反解問(wèn)題是目前仍然沒(méi)得到完美解決的問(wèn)題[6]。“形式Ⅰ”和“形式Ⅱ”兩種形式的公式可以建立起子午線弧長(zhǎng)公式與其他特殊函數(shù)的關(guān)系,如超幾何函數(shù)、雅氏橢圓函數(shù)等[10-13],有望為子午線弧長(zhǎng)的反解問(wèn)題提供新的思路。尤其是“形式Ⅱ”式 (19)的第一項(xiàng)對(duì)于一定的地球橢球?yàn)橐怀?shù),將子午線弧長(zhǎng)的反解問(wèn)題直接轉(zhuǎn)化為了第二類(lèi)橢圓積分的求逆問(wèn)題。

對(duì)于基于“形式Ⅰ”及“形式Ⅱ”的子午線弧長(zhǎng)反解問(wèn)題需要進(jìn)一步研究。

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

2 孔祥元,郭際明,劉宗泉.大地測(cè)量學(xué)基礎(chǔ)[M].武漢：武漢大學(xué)出版社,2001.(Kong Xiangyuan,Guo Jiming and Liu Zongquan.Foundation of geodesy[M].Wuhan：Wuhan University Press,2001)

3 姜晨光,閻桂峰.橢球子午線弧長(zhǎng)計(jì)算的新方法[J].測(cè)繪工程,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 牛卓立.以空間直角坐標(biāo)系為參數(shù)的子午線弧長(zhǎng)計(jì)算公式[J].測(cè)繪通報(bào),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ì)算子午線弧長(zhǎng)的數(shù)值積分法[J].測(cè)繪通報(bào), 2006,5：4-6.(Liu Xiushan.Numerical integral method for calculating meridian arc length[J].Bulletin of Surveying Mapping,2006,5：4-6)

6 易維勇,邊少峰,朱漢泉.子午線弧長(zhǎng)的解析型冪級(jí)數(shù)確定[J].測(cè)繪學(xué)院學(xué)報(bào),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 邊少峰,許江寧.計(jì)算機(jī)代數(shù)系統(tǒng)與大地測(cè)量數(shù)學(xué)分析[M].北京：國(guó)防工業(yè)出版社,2004.(Bian Shaofeng and Xu Jiangning.Computer algebra system and mathematical analysis in geodesy[M].Beijing：NationalDefence Industrial Press,2004)

8 劉仁釗,伍吉倉(cāng).任意精度的子午線弧長(zhǎng)遞歸計(jì)算[J].大地測(cè)量與地球動(dòng)力學(xué),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 王竹溪,郭敦仁.特殊函數(shù)概論[M].北京：北京大學(xué)出版社,2000.(Wang Zhuxi and Guo Dunren.Introduction to special function[M].Beijing：Peijing University Press, 2000)

11 陸源.特阿貝爾對(duì)橢圓函數(shù)論的貢獻(xiàn)[D].內(nèi)蒙古師范大學(xué),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 劉式適,劉式達(dá).特殊函數(shù) [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國(guó)家大地坐標(biāo)系橢球參數(shù)與 GRS80和WGS84的比較[J].測(cè)繪學(xué)報(bào),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

國(guó)家自然科學(xué)基金(40771117);國(guó)家農(nóng)業(yè)信息化工程技術(shù)研究中心開(kāi)放課題(KF2010W40-046)

過(guò)家春,男,1981年生,碩士,講師,主要研究大地測(cè)量學(xué).E-mail：guojiachun@ahau.edu.cn

趙秀俠,女,1981年生,博士,主要研究方向?yàn)榈貓D學(xué)與地理信息系統(tǒng).E-mail：xiuxia99@126.com

P226

A