子午線弧長公式的簡化及其泰勒級數(shù)解釋

過家春

1.安徽農(nóng)業(yè)大學 理學院，安徽 合肥230036；2.江西省數(shù)字國土重點實驗室，江西 撫州344000

1 引 言

子午線弧長計算是經(jīng)典大地測量問題之一［1－4］，圍繞這一問題的計算和應(yīng)用，近年來各國學者提出了許多新的方法和見解［5－19］。因子午線弧長問題涉及橢圓積分，不能直接求出，其經(jīng)典算法是按二項式定理展開的級數(shù)展開，國內(nèi)學者常采用的形式為［20］

式中

為提高式（2）的收斂速度，文獻［8—9］以第二偏心率e′來代換第一偏心率e，式（2）由正項級數(shù)轉(zhuǎn)為交錯級數(shù)。而國際上則多以橢球的第三扁率n（the third flattening）來代換e［1，3－4，14］。其中文獻［3，14］指出在子午線弧長的計算中，以e為參數(shù)的收斂性不及以n為參數(shù)的收斂性好。另外，文獻［10—14］則論證了子午線弧長與第二類橢圓積分的關(guān)系，文獻［12］提出調(diào)用第二類橢圓積分函數(shù)庫以達到計算任意精度子午線弧長的建議。以上研究成果豐富了子午線弧長的理論與應(yīng)用，但迄今尚未有文獻給出按二項式定理計算子午線弧長的誤差估計理論。本文通過引入高斯超幾何函數(shù)，對子午線弧長公式進一步簡化，并給出其泰勒級數(shù)解釋，分析其精度。

2 子午線弧長公式的簡化

引入?yún)?shù)橢球的第3扁率

可得e2＝4n／（1＋n）2，并由此可得

將其代入式（2），可得

根據(jù)文獻［12—13］，子午線弧長公式與高斯超幾何函數(shù)之間的關(guān)系為

式中，F(xiàn)為高斯超幾何函數(shù)的縮寫

對式（1）提取系數(shù)A′，由高斯超幾何函數(shù)F與A′的關(guān)系式（6），式（1）可化為

式中，B″、C″、…、G″等系數(shù)仍以n的冪級數(shù)形式給出

綜合式（7）—式（9），各系數(shù)取至n6，可得

至此，將子午線弧長計算公式（1）及式（2）化為公式（10），相對簡化了其系數(shù)結(jié)構(gòu)。

3 子午線弧長公式的泰勒級數(shù)解釋及其誤差估計

3.1 公式的泰勒級數(shù)解釋

由式（1）及式（2）可知，若視子午線弧長S為e的函數(shù)S＝f（e）（視e為變量），顯然S對e在以e0＝0為中心的鄰域內(nèi)無窮可導(dǎo)，因此式（1）在e0＝0處可展開為泰勒級數(shù)

該泰勒級數(shù)通過手工推導(dǎo)展開是較為困難的。筆者在數(shù)學軟件Mathematica 8.0中實現(xiàn)其級數(shù)展開，按下式輸入命令并運行

得到結(jié)果與式（1）、式（2）完全一致。式中Series［］命令為泰勒級數(shù)展開命令，其他各命令功能參見 Mathematica手冊［24］。

同理，將式（3）代入子午線弧長公式并化簡可得

在此基礎(chǔ)上，視子午線弧長S為n的函數(shù)（視n為變量），則式（13）在n0＝0處可展開為泰勒級數(shù)，按式（12）形式在 Mathematica 8.0中輸入相應(yīng)命令，并提取F，可得到與式（10）完全一致的結(jié)果。

以上分析表明按二項式定理展開與按泰勒級數(shù)展開求解子午線弧長是統(tǒng)一的，而有了子午線弧長公式的泰勒級數(shù)解釋，即可按泰勒級數(shù)的拉格朗日型余項來估計其誤差。

3.2 誤差估計

為討論方便，設(shè)F＝F（n）及

可得S＝aFY，則子午線弧長S的解算誤差由F的誤差和Y的誤差兩方面引起。記F和Y展開至n6的拉格朗日型余項分別為RF6（n）、RY6（n），則有

式中，ξ、ζ分別滿足0＜ξ＜n，0＜ζ＜n。F（7）ξ（）、Y（7）ζ（）的表達式為

顯然，應(yīng)用拉格朗日型余項（15）估計誤差的關(guān)鍵是確定F（7）（ξ）、Y（7）（ζ）的上限MF、MY。本文取MF、MY為的最大值，下面求之。

因為S＝aFY，所以F、Y均收斂，其中F為n的冪級數(shù)；對于［0，π／2］內(nèi)任意的大地緯度值，Y也為n的冪級數(shù)，由“收斂的冪級數(shù)在其收斂半徑內(nèi)逐項微分所得級數(shù)仍收斂［25］”的性質(zhì)可知F（7）（ξ）、Y（7）（ζ）也均收斂。對于F（7）（ξ），在 區(qū)間 ［0，n］上，恒有F（7）（ξ）＜0，F(xiàn)（8）（ξ）＞0，所以F（7）（ξ）為單調(diào)遞增函數(shù)，故

對于Y（7）（ζ），由式（16）可知其為周期性函數(shù)，理論上其最小正周期隨緯度倍角的增加而趨于無窮小。但事實上，經(jīng)驗證，隨著式（16）的展開，其緯度倍角的正弦的系數(shù)將逐漸趨于0，sin 14B以后各項的周期性影響很小，Y（7）的最小正周期趨于π／7。應(yīng)用多元函數(shù)極值分析，可求得其最大值725／32，具體過程從略。

再由誤差傳播定律可得

式（17）即為公式（10）的誤差估計。

同理，對子午線弧長公式（1）也可作泰勒級數(shù)解釋，其誤差估計為

式中

與F（7）ξ（）情況類似，容易證明Me為收斂的冪級數(shù)。類似的，可得展開至其他各階的誤差估計式，結(jié)構(gòu)與上述公式類似，此不贅述。

比較式（17）、式（18）兩誤差估計式，可得

對地球橢球而言，ΔSn／ΔSe≈1／1.12×104?？梢姼倪M后的公式（10）精度提高顯著。其主要原因在于：以n替換e，并提取高斯超幾何函數(shù)后，原公式（1）中各系數(shù)由以e為參數(shù)的正項級數(shù)轉(zhuǎn)換為公式（10）中以n為參數(shù)的交錯級數(shù)，且n僅約為e2的1／4，及至n6項，僅約e12的1／4014，加速了各項系數(shù)的收斂速度，從而有效提高了子午線弧長的計算精度。

下面以具體實例進行驗證分析。

4 算例驗證分析

為驗證公式泰勒級數(shù)解釋下估計誤差的正確性，以 WGS－84橢球參數(shù)為例，按文獻［12］給出的子午線弧長與第二類橢圓積分關(guān)系，在Mathematica 8.0中調(diào)用第二類橢圓積分函數(shù)，得到子午線弧長的任意精度解，進而得到分別按公式（1）展開至e6sin 6B、e8sin 8B、e10sin 10B、e12sin 12B項及公式（10）展開至n3sin 6B、n4sin 8B、n5sin 10B、n6sin 12B項的子午線弧長解算誤差，誤差對比曲線如圖1所示，最大估計誤差與實際最大絕對誤差對比列于表1。

圖1 子午線弧長解算誤差曲線圖Fig.1 Error curves of the solution of the meridian arc length

表1 最大估計誤差與實際最大絕對誤差驗證對比Tab.1 Verification and comparison of the maximum estimation errors and the worst absolute actual errors m

驗算結(jié)果表明，按公式（10）分別展開至n3sin 6B、n4sin 8B、n5sin 10B、n6sin 12B項比按公式（1）分別展開至e6sin 6B、e8sin 8B、e10sin 10B、e12sin 12B項的精度提高了3～4個數(shù)量級，精度提高顯著，且驗算誤差與估計誤差吻合，驗證了按泰勒級數(shù)解釋所給誤差估計理論的正確性。

5 結(jié) 語

本文通過引入新的參數(shù)，得到了子午線弧長的簡化形式，加快了收斂速度，精度提高顯著，體現(xiàn)了高斯超幾何函數(shù)的引入對子午線弧長在公式簡化及精度提高方面的意義，而本文給出的子午線弧長公式的泰勒級數(shù)解釋則給子午線弧長公式的誤差分析提供了理論依據(jù)。

［1］ TORGE W.Geodesy［M］.3rd ed.Berlin：Walter De Gruyter，2001：91－98.

［2］ KONG Xiangyuan，GUO Jiming，LIU Zongquan.Foundation of Geodesy［M］.Wuhan：Wuhan University Press，2001：64－73.（孔祥元，郭際明，劉宗泉.大地測量學基礎(chǔ)［M］.武漢：武漢大學出版社，2001：64－73.）

［3］ HELMERT F R.Die Mathematischen und Physikalischen Theorien der H?heren Geod?sie（English Translation Version）［M］.Leipzig：Printing and Publishing House of B.G.Teubner，1880：46－48.

［4］ DEAKIN R E，Hunter M N.Geometric Geodesy：Part A［R］.Melbourne：RMIT University，2010：60－77.

［5］ DING Jiabo.The Transforming the Zones of Gauss Projection from Latitudes of Low Points［J］.Acta Geodaetica et Cartographica Sinica，1993，22（3）：212－217.（丁佳波.利用底點緯度進行高斯投影換代計算［J］.測繪學報，1993，22（3）：212－217.）

［6］ LI Houpu，BIAN Shaofeng.The Expressions of Gauss Projection by Complex Numbers［J］.Acta Geodaetica et Cartographica Sinica，2008，37（1）：5－9.（李厚樸，邊少峰.高斯投影的復(fù)變函數(shù)表示［J］.測繪學報，2008，37（1）：5－9.）

［7］ YANG Yuanxi.The Respective Roles and Contributions of Various Observation in Integrated Geodesy［J］.Acta Geodaetica et Cartographica Sinica，1989，18（3）：232－238.（楊元喜.整體大地測量中各類觀測值的分工與貢獻［J］.測繪學報，1989，18（3）：232－238.）

［8］ CHENG Pengfei，WEN Hanjiang，CHENG Yingyan，et al.Parameters of the CGCS 2000Ellipsoid and Comparisons with GRS 80and WGS 84［J］.Acta Geodaetica et Cartographica Sinica，2009，38（3）：189－194.（程鵬飛，成英燕，文漢江，等.2000國家大地坐標系橢球參數(shù)與GRS 80和 WGS 84的比較［J］.測繪學報，2009，38（3）：189－194.）

［9］ LIU Zhengcai.Simplification of Formula of Meridian Arc Length ＆Program of Gauss Projection［J］.Engineering of Surveying and Mapping，2001，10（1）：55－56.（劉正才.子午線弧長公式的簡化及通用高斯投影計算程序介紹［J］.測繪工程，2001，10（1）：55－56.）

［10］ DORRER E.From Elliptic Arc Length to Gauss－Krüger Coordinates by Analytical Continuation［C］∥Geodesy：The Challenge of the 3rd Millennium.Berlin：Springer，2003：293－298.

［11］ BERMEJO－SOLERA M，OTERO J.Simple and Highly Accurate Formulas for the Computation of Transverse Mercator Coordinates from Longitude and Isometric Latitude［J］.Journal of Geodesy，2009，83：1－12.

［12］ GUO Jiachun，ZHAO Xiuxia，XU Li，etal.Calculating Meridian Arc Length by Transforming Its Formula into Elliptic Integral of Second Kind［J］.Journal of Geodesy and Geodynamics，2011，31（4）：94－98.（過家春，趙秀俠，徐麗，等.基于第二類橢圓積分的子午線弧長公式變換及解算［J］.大地測量與地球動力學，2011，31（4）：94－98.）

［13］ GUO Jiachun.New Method for Inverse Solution of Meridian Based on Elliptic Integral of the Second Kind［J］.Journal of Geodesy and Geodynamics，2012，32（3）：116－120.（過家春.基于第二類橢圓積分的子午線弧長反解新方法［J］.大地測量與地球動力學，2012，32（3）：116－120.）

［14］ KAWASE K.A General Formula for Calculating Meridian Arc Length and Its Application to Coordinate Conversion in the Gauss－Krüger Projection ［J］.Bulletin of the Geospatial Information Authority of Japan，2011，59：1－13.

［15］ YI Weiyong，BIAN Shaofeng，ZHU Hanquan.Determination of Foot Point Latitude by Analytic Positive Series［J］.Journal of Institute of Surveying and Mapping，2000，17（3）：167－171.（易維勇，邊少峰，朱漢泉.子午線弧長的解析型冪 級 數(shù) 確 定 ［J］.測 繪 學 院 學 報，2000，17（3）：167－171.）

［16］ LIU Renzhao，WU Jicang.Recursive Computation of Meridian Arc Length with Discretionary Precision［J］.Journal of Geodesy and Geodynamics，2007，27（5）：59－62.（劉仁釗，伍吉倉.任意精度的子午線弧長遞歸計算［J］.大地測量與地球動力學，2007，27（5）：59－62.）

［17］ BIAN S F，CHEN Y B.Solving an Inverse Problem of a Meridian Arc in Terms of Computer Algebra System［J］.Journal of Surveying Engineering，2006，132（1）：7－10.

［18］ NIU Zhuoli.Formulae for Calculation of Meridian Arc Length by the Parameters of Space Rectangular Coordinates［J］.Bulletin of Surveying and Mapping，2001（11）：14－15.（牛卓立.以空間直角坐標系為參數(shù)的子午線弧長計算公式［J］.測繪通報，2001（11）：14－15.）

［19］ LIU Xiushan.Numerical Integral Method Calculating Meridian Arc Length［J］.Bulletin of Surveying and Mapping，2006（5）：4－6.（劉修善.計算子午線弧長的數(shù)值積分法［J］.測繪通報，2006（5）：4－6.）

［20］ CHENG Pengfei，CHENG Yingyan，WEN Hanjiang，et al.Practical Manual on CGCS2000 ［M］. Beijing：Surveying and Mapping Press，2008：147－148.（程鵬飛，成英燕，文漢江，等.2000國家大地坐標系實用寶典［M］.北京：測繪出版社，2008：147－148.）

［21］ LIU Shishi，LIU Shida.Special Function［M］.Beijing：China Meteorological Press，1988：656－745.（劉式適，劉式達.特殊函數(shù)［M］.北京：氣象出版社，1988：656－745.）

［22］ ZHU Huatong.Review on the Methods for Calculating Latitude of Low Points［J］.Bulletin of Surveying and Mapping，1978（5）：10－14.（朱華統(tǒng).底點緯度計算方法評述［J］.測繪通報，1978（5）：10－14.）

［23］ SUN Qun，YANG Qihe.The Research on the Computation of the Foot－point Latitude and the Inverse Solution of Isometric Latitude and Function［J］.Journal of Institute of Surveying ang Mapping，1985（2）：64－75.（孫群，楊啟和.底點緯度解算以及等量緯度和面積函數(shù)反解問題的探討［J］.測繪學院學報，1985（2）：64－75.）

［24］ WOLFRAM S.The Mathematica Book［M］.5th ed.Champaign：Wolfram Media Inc，2003.

［25］ BRONSHTEIN I N，SEMENDIAEY KA，HIRSCH K A.Handbook of Mathematics［M］.4th ed.New York：Van Nostrand Reinhold，2003：412－416.