劉學杰，楊麗坤

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

子午線弧長的計算方法及精度分析

劉學杰1，楊麗坤2

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

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

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

計算子午線弧長是橢球大地測量學的一項基本內容，是高斯-克呂格投影坐標正算的基礎問題。子午線弧長的計算公式本質上是一個橢圓函數(shù)積分公式，由于被積函數(shù)的原函數(shù)無法用初等函數(shù)的形式予以表示，因此無法得到嚴格的解析結果。傳統(tǒng)做法是基于二項式定理將被積函數(shù)展開為級數(shù)形式，通過分項積分得到近似的解析結果。

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

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

1 計算子午線弧長的三類算法

1.1 級數(shù)展開算法

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

(1)

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

(2)

式(2)無法直接用定積分進行求解，因為被積函數(shù)M的原函數(shù)無法用初等函數(shù)形式進行表達。傳統(tǒng)做法是將被積函數(shù)M按照二項式定理展開為sin2B的冪級數(shù)，并將sinB偶次冪函數(shù)按照三角函數(shù)的倍角公式展開為余弦的倍數(shù)函數(shù)，轉換為可積函數(shù)，然后逐項進行積分。考慮到計算結果目的和精度要求，截斷之后得到以下實用計算公式

(3)

(4)

(5)

對于我國常用的4個橢球(克拉索夫斯基橢球、IUGG-75橢球、WGS-84橢球、CGCS2000橢球)而言，式(3)最后一項的最大值均未超過0.03 mm，因此截斷誤差不超過0.1 mm，完全滿足測量工作的精度要求。

1.2 數(shù)值積分算法

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

(6)

復合Simpson算法公式為

(7)

Romberg算法公式為

(8)

Gauss-Legendre算法公式為

(9)

(10)

(11)

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

標準的一階常微分方程形式為

(12)

子午線弧長計算公式本質上即為一個標準的一階常微分方程，結合橢球實際，式(1)可寫為如下形式

(13)

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

(14)

改進的Euler算法公式為

(15)

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

(16)

三階的Runge-Kutta算法公式為

(17)

四階的Runge-Kutta算法公式為

(18)

2 計算結果分析與評價

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

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

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

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

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

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

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

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

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

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

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

3 結 語

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

[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] 過家春，趙秀俠，徐麗，等．基于第二類橢圓積分的子午線弧長公式變換及解算[J]．大地測量與地球動力學，2011，31(4)：94-98．

[5] 過家春．子午線弧長公式的簡化及其泰勒級數(shù)解釋[J]．測繪學報，2014，43(2)：125-130．

[6] 牛卓立．以空間直角坐標為參數(shù)的子午線弧長計算公式[J]．測繪通報，2001(11)：14-15．

[7] 嚴伯鐸．橢球子午線弧長的一種計算方法[J]．地礦測繪，2003，19(3)：7-10．

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

[9] 劉仁釗，伍吉倉．任意精度的子午線弧長遞歸計算[J]．大地測量與地球動力學，2007，27(5)：59-62.

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

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

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

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

[15] 宋葉志，賈東永．MATLAB數(shù)值分析與應用[M]．北京：機械工業(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

劉學杰，楊麗坤.子午線弧長的計算方法及精度分析[J].測繪通報，2017(8)：106-109.

10.13474/j.cnki.11-2246.2017.0264.

2017-02-20

2016年國家重點研發(fā)計劃(2016YFC0803103)；河南省高校創(chuàng)新團隊支持計劃(14IRTSTHN026)；河南省創(chuàng)新型科技創(chuàng)新團隊支持計劃

劉學杰(1968—)，男，高級工程師，主要研究方向為測繪科學與技術。E-mail：zwchlxj@126.com

P226

A

0494-0911(2017)08-0106-04