過家春 趙秀俠 徐 麗 田勁松 高 飛

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

基于第二類橢圓積分的子午線弧長公式變換及解算＊

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

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

通過推導,將子午線弧長公式變換為基于第二類橢圓積分的兩種形式：“形式Ⅰ”將子午線弧長公式表達為一個有理函數和第二類橢圓積分之和,建立了以大地緯度B為自變量的子午線弧長公式與第二類橢圓積分之間的關系;“形式Ⅱ”給出了以歸化緯度μ為自變量、直接利用第二類橢圓積分計算子午線弧長的公式。利用此兩種形式的子午線弧長公式,在Matlab中編寫程序,調用第二類橢圓積分函數 EllipticE(x,k)計算子午線弧長,精度和計算效率均優(yōu)于經典算法。對 CGCS2000所采用的地球橢球子午線弧長的計算表明,此兩種形式的子午線弧長公式建立了子午線弧長公式與第二類橢圓積分的關系,結構簡潔,易于展開,一定程度上完善了子午線弧長理論,且便于手工計算及計算機程序實現。

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

1 引言

子午線弧長是大地測量學中的一個基本量。計算從赤道開始到任意大地緯度 B的子午線弧長 S的基本公式為

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

顯然,子午線弧長的計算涉及到勒讓德第二類橢圓積分(簡稱為第二類橢圓積分),其原函數無法用初等函數的形式表達,不能直接求出。目前,世界各國對子午線弧長的經典計算方法是將 a(1-e2) (1-e2sinB)按二項式定理展開為級數形式,然后再逐項積分,得出一定精度的計算公式[1,2]：

式中,A0、B0、C0、D0、…為一組與 e有關的常系數,具體表達式可參閱文獻[2]。

在此基礎上,可以計算任意區(qū)間 [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執(zhí)行時間對比,結果表明,本文算法的計算效率提高到了 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程序代碼：

該程序中用于子午線弧長計算的代碼僅一行,且精度沒有損失?；凇靶问舰颉钡某绦蚪Y構與其類似,計算結果完全一致。

將程序中的橢球參數替換為中國 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]?！靶问舰瘛焙汀靶问舰颉眱煞N形式的公式可以建立起子午線弧長公式與其他特殊函數的關系,如超幾何函數、雅氏橢圓函數等[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 邊少峰,許江寧.計算機代數系統(tǒng)與大地測量數學分析[M].北京：國防工業(yè)出版社,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);國家農業(yè)信息化工程技術研究中心開放課題(KF2010W40-046)

過家春,男,1981年生,碩士,講師,主要研究大地測量學.E-mail：guojiachun@ahau.edu.cn

趙秀俠,女,1981年生,博士,主要研究方向為地圖學與地理信息系統(tǒng).E-mail：xiuxia99@126.com

P226

A