子午線弧長(zhǎng)的計(jì)算方法及精度分析

劉學(xué)杰，楊麗坤

(1. 河南省中緯測(cè)繪規(guī)劃信息工程有限公司，河南 焦作 454000； 2. 鄭州工業(yè)貿(mào)易學(xué)校，河南 鄭州 450007)

子午線弧長(zhǎng)的計(jì)算方法及精度分析

劉學(xué)杰1，楊麗坤2

(1. 河南省中緯測(cè)繪規(guī)劃信息工程有限公司，河南 焦作 454000； 2. 鄭州工業(yè)貿(mào)易學(xué)校，河南 鄭州 450007)

計(jì)算子午線弧長(zhǎng)除了采用經(jīng)典的級(jí)數(shù)展開算法之外，還可通過(guò)數(shù)值積分與常微分方程數(shù)值解法進(jìn)行求解。為評(píng)價(jià)各種算法的精度，本文選取大地緯度自0°—90°、間隔距離為1°、1′、1″的3組樣本數(shù)據(jù)，分別基于傳統(tǒng)算法、數(shù)值積分算法和常微分方程數(shù)值算法3大類11種算法計(jì)算得到各組樣本所對(duì)應(yīng)的子午線弧長(zhǎng)結(jié)果，并從算法精度和運(yùn)算速度兩個(gè)方面對(duì)各種數(shù)值算法進(jìn)行了分析與評(píng)價(jià)。實(shí)例表明三階、四階Runge-Kutta算法不僅精度高，而且運(yùn)算效率是其他算法的2倍多，研究結(jié)果為計(jì)算子午線弧長(zhǎng)的提供了有效的算法模型。

子午線弧長(zhǎng)；數(shù)值積分；常微分方程；展開算法

計(jì)算子午線弧長(zhǎng)是橢球大地測(cè)量學(xué)的一項(xiàng)基本內(nèi)容，是高斯-克呂格投影坐標(biāo)正算的基礎(chǔ)問(wèn)題。子午線弧長(zhǎng)的計(jì)算公式本質(zhì)上是一個(gè)橢圓函數(shù)積分公式，由于被積函數(shù)的原函數(shù)無(wú)法用初等函數(shù)的形式予以表示，因此無(wú)法得到嚴(yán)格的解析結(jié)果。傳統(tǒng)做法是基于二項(xiàng)式定理將被積函數(shù)展開為級(jí)數(shù)形式，通過(guò)分項(xiàng)積分得到近似的解析結(jié)果。

近年來(lái)，國(guó)內(nèi)外相關(guān)學(xué)者對(duì)此問(wèn)題提出了一些新的算法，文獻(xiàn)[1—5]基于第二類橢圓積分，將子午線弧長(zhǎng)公式表達(dá)為有理函數(shù)和第二類橢圓積分之和等兩類形式不同但結(jié)果等價(jià)的公式，或是引入橢球的第三扁率和高斯超幾何函數(shù)，給出解算公式的簡(jiǎn)化形式，完善了子午線弧長(zhǎng)計(jì)算方法。文獻(xiàn)[6]給出了以空間直角坐標(biāo)為參數(shù)的子午線弧長(zhǎng)計(jì)算公式。文獻(xiàn)[7—8]給出了子午線弧長(zhǎng)計(jì)算公式的直接展開式。文獻(xiàn)[9]基于遞歸關(guān)系給出了任意精度的子午線弧長(zhǎng)計(jì)算公式，可滿足不同精度要求。文獻(xiàn)[10]基于數(shù)值積分算法對(duì)子午線弧長(zhǎng)進(jìn)行了計(jì)算，但結(jié)果表明數(shù)值積分算法與傳統(tǒng)級(jí)數(shù)展開算法存在較大差異。文獻(xiàn)[11]利用復(fù)合Simpson積分公式表明數(shù)值積分算法與傳統(tǒng)級(jí)數(shù)展開算法結(jié)果基本一致。

子午線弧長(zhǎng)計(jì)算公式直觀而言是一個(gè)橢圓函數(shù)積分問(wèn)題，但其實(shí)也是一個(gè)標(biāo)準(zhǔn)的一階常微分方程求解問(wèn)題?，F(xiàn)代數(shù)值分析對(duì)于這兩個(gè)數(shù)學(xué)問(wèn)題有著成熟的解算方法，前者通過(guò)數(shù)值積分予以解決，后者則基于常微分方程數(shù)值解法進(jìn)行解算，兩類方法的計(jì)算過(guò)程中都沒(méi)有涉及深?yuàn)W的數(shù)學(xué)理論知識(shí)，有完備的算法及公式，程序設(shè)計(jì)簡(jiǎn)單易行。本文利用這兩類方法對(duì)子午線弧長(zhǎng)進(jìn)行了求解，并分析比較了各種算法的精度及運(yùn)算速度，對(duì)算法質(zhì)量進(jìn)行了比較與評(píng)價(jià)。

1 計(jì)算子午線弧長(zhǎng)的三類算法

1.1 級(jí)數(shù)展開算法

根據(jù)橢球大地測(cè)量中曲線弧長(zhǎng)與曲率半徑的基本關(guān)系[12]，子午線上某點(diǎn)P從赤道到該點(diǎn)的子午線弧長(zhǎng)X與點(diǎn)P的大地緯度B之間滿足如下微分關(guān)系

(1)

式中，a、e、M分別為橢球的長(zhǎng)半徑、第一偏心率、子午線曲率半徑。則有

(2)

式(2)無(wú)法直接用定積分進(jìn)行求解，因?yàn)楸环e函數(shù)M的原函數(shù)無(wú)法用初等函數(shù)形式進(jìn)行表達(dá)。傳統(tǒng)做法是將被積函數(shù)M按照二項(xiàng)式定理展開為sin2B的冪級(jí)數(shù)，并將sinB偶次冪函數(shù)按照三角函數(shù)的倍角公式展開為余弦的倍數(shù)函數(shù)，轉(zhuǎn)換為可積函數(shù)，然后逐項(xiàng)進(jìn)行積分?？紤]到計(jì)算結(jié)果目的和精度要求，截?cái)嘀蟮玫揭韵聦?shí)用計(jì)算公式

(3)

(4)

(5)

對(duì)于我國(guó)常用的4個(gè)橢球(克拉索夫斯基橢球、IUGG-75橢球、WGS-84橢球、CGCS2000橢球)而言，式(3)最后一項(xiàng)的最大值均未超過(guò)0.03 mm，因此截?cái)嗾`差不超過(guò)0.1 mm，完全滿足測(cè)量工作的精度要求。

1.2 數(shù)值積分算法

子午線弧長(zhǎng)公式本質(zhì)上是一個(gè)橢圓函數(shù)的積分問(wèn)題，可采用數(shù)值積分直接進(jìn)行求解。數(shù)值積分有諸如復(fù)合梯形算法、復(fù)合Simpson算法、Romberg算法、Gauss-Legendre算法、Monte Carlo算法等幾種成熟的算法[13～14]，但計(jì)算精度與效率各不相同，需要進(jìn)行分析與比較。

(6)

復(fù)合Simpson算法公式為

(7)

Romberg算法公式為

(8)

Gauss-Legendre算法公式為

(9)

(10)

(11)

1.3 常微分方程數(shù)值解法

標(biāo)準(zhǔn)的一階常微分方程形式為

(12)

子午線弧長(zhǎng)計(jì)算公式本質(zhì)上即為一個(gè)標(biāo)準(zhǔn)的一階常微分方程，結(jié)合橢球?qū)嶋H，式(1)可寫為如下形式

(13)

對(duì)于式(12)的標(biāo)準(zhǔn)一階常微分方程，常用的數(shù)值計(jì)算方法主要有Euler算法、改進(jìn)的Euler算法，以及二階、三階、四階Runge-Kutta算法[13-14]。

(14)

改進(jìn)的Euler算法公式為

(15)

二階Runge-Kutta算法又稱為變形的Euler算法，其公式為

(16)

三階的Runge-Kutta算法公式為

(17)

四階的Runge-Kutta算法公式為

(18)

2 計(jì)算結(jié)果分析與評(píng)價(jià)

為分析對(duì)比上述3大類11種算法的計(jì)算精度與運(yùn)算效率，本文以2000國(guó)家大地坐標(biāo)系的橢球基準(zhǔn)CGCS2000橢球(a=6 378 137 m、f=1/298.257 222 101、e2=2f-f2)為例，選擇了3組大地緯度的樣本數(shù)據(jù)，分別計(jì)算得到了各組樣本相應(yīng)的結(jié)果，并統(tǒng)計(jì)了相應(yīng)的計(jì)算時(shí)間(見表1)。

表1 3組大地緯度樣本數(shù)據(jù)基本情況

基于Matlab平臺(tái)[15]，根據(jù)式(3)—式(18)分別對(duì)上述3大類11種算法予以編程實(shí)現(xiàn)。在對(duì)3組樣本數(shù)據(jù)計(jì)算過(guò)程中，當(dāng)數(shù)值算法程序涉及區(qū)間等分?jǐn)?shù)n與步長(zhǎng)h時(shí)，需要相應(yīng)樣本數(shù)據(jù)中的樣本數(shù)目和間隔距離保持一致。

根據(jù)各種算法計(jì)算得到3組樣本數(shù)據(jù)的子午線弧長(zhǎng)后，將同一樣本中各種數(shù)值算法結(jié)果與級(jí)數(shù)展開算法結(jié)果進(jìn)行求差，差值即可反映各種數(shù)值算法的精度情況，具體結(jié)果見表2。

采用Monte Carlo算法，3組樣本中差值的最大值分別達(dá)到4998、1670和301 m，表明該算法結(jié)果錯(cuò)誤，表2中不再列出具體數(shù)值。對(duì)于Gauss-Legendre算法，分別計(jì)算了二階、三階、四階、五階的計(jì)算結(jié)果，其中二階、三階、四階計(jì)算精度較差，因此表2中僅給出了五階算法的計(jì)算結(jié)果。對(duì)于Romberg算法，計(jì)算限差取10-5(對(duì)應(yīng)計(jì)算結(jié)果截至0.01 mm)。

各種算法的運(yùn)算速度也是衡量算法有效性的一個(gè)重要指標(biāo)，現(xiàn)將3組樣本中各種數(shù)值算法的計(jì)算時(shí)間列于表2中。其中第1—2組樣本中傳統(tǒng)算法的計(jì)算耗時(shí)均未超過(guò)1 s，第3組樣本中傳統(tǒng)算法的計(jì)算耗時(shí)為458 s。

從表2可以看出，在算法的計(jì)算精度方面，只有Romberg算法和三階、四階的Runge-Kutta算法的計(jì)算結(jié)果與傳統(tǒng)展開算法結(jié)果間的差值在0.1 mm以內(nèi)，能夠滿足測(cè)量工作的需要；而復(fù)合梯形算法、復(fù)合Simpson算法、Euler算法、改進(jìn)的Euler算法及二階的Runge-Kutta算法，隨著步長(zhǎng)h的縮小、區(qū)間等分點(diǎn)n的增多，計(jì)算精度隨之提高，但Euler算法直至步長(zhǎng)h=1″差值仍大于155 mm，由此可知這些算法精度不能滿足測(cè)量工作的需要。Gauss-Legendre算法結(jié)果精度與步長(zhǎng)、等分點(diǎn)數(shù)無(wú)關(guān)，計(jì)算精度僅隨階數(shù)的增大而提高。

表2 各種數(shù)值算法與傳統(tǒng)算法結(jié)果間差值及數(shù)值算法計(jì)算耗時(shí)統(tǒng)計(jì)

注：表2第一列中的梯形、Sim、Rom、GL5、Eu、Eu1、RK2、RK3、RK4分別代表復(fù)合梯形、復(fù)合Simpson、Romberg、五階Gauss-Legendre、Euler、改進(jìn)Euler、二階Runge-Kutta、三階Runge-Kutta、四階Runge-Kutta算法。

在算法的運(yùn)算速度方面，由于第1組和第2組樣本數(shù)目較少，各種算法的計(jì)算耗時(shí)基本都在1 s之內(nèi)，難以反映出算法的效率，但從第2組樣本的計(jì)算耗時(shí)可以看出復(fù)合梯形、復(fù)合Simpson算法的運(yùn)算速度要明顯慢于其他幾種算法。第3組樣本數(shù)目多達(dá)324 001個(gè)，各種數(shù)值算法的運(yùn)算速度可以非常清晰區(qū)分開來(lái)：Romberg算法、Gauss-Legendre算法的運(yùn)算速度與傳統(tǒng)算法大體相當(dāng)；傳統(tǒng)算法速度分別是復(fù)合梯形算法、復(fù)合Simpson算法的165倍和345倍，這兩種數(shù)值算法速度最慢，不適用于大數(shù)據(jù)的計(jì)算；而常微分方程的5種數(shù)值解法運(yùn)算速度大致相同，均為傳統(tǒng)算法速度的2倍多。

綜合3組樣本數(shù)據(jù)各種數(shù)值算法結(jié)果精度與運(yùn)算速度兩方面因素，三階、四階的Runge-Kutta算法最優(yōu)，不僅計(jì)算精度高，而且運(yùn)算速度是傳統(tǒng)算法的2倍多，優(yōu)于其他幾種數(shù)值積分算法。而復(fù)合梯形算法、復(fù)合Simpson算法雖然可以通過(guò)減小步長(zhǎng)及增大區(qū)間等分點(diǎn)數(shù)來(lái)提高計(jì)算結(jié)果的精度，但運(yùn)算速度卻急劇降低。

3 結(jié) 語(yǔ)

本文基于傳統(tǒng)的展開算法、數(shù)值積分和常微分方程數(shù)值解法分別對(duì)子午線弧長(zhǎng)進(jìn)行了計(jì)算，并從計(jì)算精度、運(yùn)算速度兩個(gè)方面對(duì)各種算法的有效性進(jìn)行了分析與評(píng)價(jià)。試驗(yàn)結(jié)果表明，采用常微分方程數(shù)值解法比傳統(tǒng)算法速度快、精度高；而數(shù)值積分雖然可以通過(guò)減小步長(zhǎng)來(lái)提高結(jié)果精度，但計(jì)算速度低，難以適用于大數(shù)據(jù)計(jì)算。本文研究結(jié)果為計(jì)算子午線弧長(zhǎng)提供了一條新的途徑和借鑒。

[2] DORRER E.From Elliptic Arc Length to Gauss-Kruger Coordinates by Analytical Continuation[C]∥Geodesy: The Challenge of the 3rd Millennium. Berlin: Springer, 2003: 293-298.

[3] KAWAS K.A General Formula for Calculating Meridian Arc Length and Its Application to Coordinate Conversion in the Gauss-Kruger Projection [J]. Bulletin of the Geospatial Information Authority of Japan, 2011, 59(2):1-13.

[4] 過(guò)家春，趙秀俠，徐麗，等．基于第二類橢圓積分的子午線弧長(zhǎng)公式變換及解算[J]．大地測(cè)量與地球動(dòng)力學(xué)，2011，31(4)：94-98．

[5] 過(guò)家春．子午線弧長(zhǎng)公式的簡(jiǎn)化及其泰勒級(jí)數(shù)解釋[J]．測(cè)繪學(xué)報(bào)，2014，43(2)：125-130．

[6] 牛卓立．以空間直角坐標(biāo)為參數(shù)的子午線弧長(zhǎng)計(jì)算公式[J]．測(cè)繪通報(bào)，2001(11)：14-15．

[7] 嚴(yán)伯鐸．橢球子午線弧長(zhǎng)的一種計(jì)算方法[J]．地礦測(cè)繪，2003，19(3)：7-10．

[8] 李厚樸，劉敏，孔海英，等．子午線弧長(zhǎng)和等面積緯度函數(shù)變換的直接展開式[J]．海洋測(cè)繪，2011，31(1)：17-19．

[9] 劉仁釗，伍吉倉(cāng)．任意精度的子午線弧長(zhǎng)遞歸計(jì)算[J]．大地測(cè)量與地球動(dòng)力學(xué)，2007，27(5)：59-62.

[10] 劉修善．計(jì)算子午線弧長(zhǎng)的數(shù)值積分法[J]．測(cè)繪通報(bào)，2006(5)：4-6．[11] 楊雙富．再議計(jì)算子午線弧長(zhǎng)的數(shù)值積分法[J]．城市勘測(cè)，2010(6)：140-142．

[12] 孔祥元，郭際明，劉宗泉．大地測(cè)量學(xué)基礎(chǔ)[M]．武漢：武漢大學(xué)出版社，2007：106-116．

[13] 李慶揚(yáng)，王能超，易大義．?dāng)?shù)值分析 [M]．4版．北京：清華大學(xué)出版社，2001：336-359．

[14] 易大義，沈云寶，李有法．計(jì)算方法[M]．杭州：浙江大學(xué)出版社，2000：243-266．

[15] 宋葉志，賈東永．MATLAB數(shù)值分析與應(yīng)用[M]．北京：機(jī)械工業(yè)出版社，2010：380-410．

Calculation Methods and Accuracy Analysis of Meridian Arc Length

LIU Xuejie1，YANG Likun2

(1. Zhongwei Surveying and Planning Information Engineering Co.,Ltd. of Henan Province, Jiaozuo 454000，China;2. Zhengzhou Industry and Trade School, Zhengzhou 450007，China)

There are several kinds of algorithms for calculating the meridian arc length except the classical expanded algorithm, such as numerical integration and ordinary differential equations numerical solution. In order to study the accuracy of each algorithm, this paper selected 3 sets of sample data within the geodetic latitude from 0° to 90°, whose intervals are 1°,1′,1″, respectively. Based on the traditional expanded algorithms, numerical integral algorithms and numerical solution of ordinary differential equations, the corresponding meridian arc length results were calculated and the quality of each numerical algorithm with regard to algorithm accuracy and computation speed were evaluated. The results show that the 3 and 4 order Runge-Kutta algorithm not only have high precision but the computing speed is twice more than other algorithms. This paper provides new, reliable algorithm with high speed for meridian arc length calculation. The results provide an effective algorithm for calculating the meridian arc length.

meridian arc length; numerical integral; ordinary differential equations; expanded algorithm

劉學(xué)杰，楊麗坤.子午線弧長(zhǎng)的計(jì)算方法及精度分析[J].測(cè)繪通報(bào)，2017(8)：106-109.

10.13474/j.cnki.11-2246.2017.0264.

2017-02-20

2016年國(guó)家重點(diǎn)研發(fā)計(jì)劃(2016YFC0803103)；河南省高校創(chuàng)新團(tuán)隊(duì)支持計(jì)劃(14IRTSTHN026)；河南省創(chuàng)新型科技創(chuàng)新團(tuán)隊(duì)支持計(jì)劃

劉學(xué)杰(1968—)，男，高級(jí)工程師，主要研究方向?yàn)闇y(cè)繪科學(xué)與技術(shù)。E-mail：zwchlxj@126.com

P226

A

0494-0911(2017)08-0106-04