子午線弧長計算的數(shù)值積分算法及其比較

鄭紅曉 張紅方 雷偉偉

(1.河南省中緯測繪規(guī)劃信息工程有限公司，河南焦作 454000；2.河南理工大學測繪與國土信息工程學院，河南焦作 454000)

子午線弧長計算的數(shù)值積分算法及其比較

鄭紅曉1張紅方1雷偉偉2

(1.河南省中緯測繪規(guī)劃信息工程有限公司，河南焦作 454000；2.河南理工大學測繪與國土信息工程學院，河南焦作 454000)

子午線弧長計算的經典算法是對子午線曲率半徑按照牛頓二項式定理進行展開，分項積分得到近似解析解。研究了五種常用的數(shù)值積分算法及其在子午線弧長計算中的應用，并用Matlab軟件予以實現(xiàn)。將數(shù)值積分結果與經典算法結果進行比較，結果表明：利用數(shù)值積分算法求解子午線弧長，簡單易行，準確可靠。同時，還證明了復合辛普森算法和龍貝格算法要優(yōu)于其它幾種算法。

子午線弧長、數(shù)值積分算法、Matlab

子午線弧長計算是橢球大地測量學的一項重要內容，是高斯投影坐標計算的基礎。但因子午線弧長計算公式中的被積函數(shù)(即子午線曲率半徑M)無法找到其原函數(shù)，故而不能直接積分，經典的解法是采用牛頓二項式定理進行展開，分項積分得到近似解析解[1-2]。

近些年，國內一些學者對該問題提出一些解算方法[6-10]：文獻[6]采用常微分方程的數(shù)值解法，得到了理想的結果；文獻[7]基于第二類橢圓積分進行求解，結果理想但過程過于復雜，不便應用；文獻[8]給出了任意精度的子午線弧長遞歸計算公式，可滿足不同精度的弧長計算；文獻[9]利用幾何分析法進行了探討；文獻[10]對子午線弧長的數(shù)值積分法進行了研究，但得到的數(shù)值積分結果與經典算法結果相差幾百米甚至上萬米，這與數(shù)值積分算法的準確性相矛盾，不免讓人有所疑惑。本文擬就數(shù)值積分算法求解子午線弧長再次進行研究，以驗證此類算法究竟是否準確可靠。

1 子午線弧長計算的經典算法

大地測量學經典教材[1]給出了計算子午線弧長的基本公式：從赤道到大地緯度為B處的子午線弧長為

(1)

式(1)中，M為子午線曲率半徑，a為橢球長半軸，e為橢球第一偏心率。

對于(1)式而言，由于被積函數(shù)結構復雜，其原函數(shù)無法求出，故不能直接用牛頓-萊布尼茨積分進行計算。在經典教材中，采用近似解析法，即把被積函數(shù)M按照牛頓二項式定理展開為e的冪級數(shù)，并將正弦的冪函數(shù)展開為余弦的倍數(shù)函數(shù)，然后逐項進行積分?？紤]到計算結果的目的和精度，得到以下兩個實用計算公式

(2)

(3)

(4)

(5)

根據(jù)式(3)計算的結果已被諸多文獻[6-8]證明了其準確性，故本文以經典算法的結果作為準確值，用以衡量各種數(shù)值積分算法的準確度。

2 數(shù)值積分算法

數(shù)值積分算法[3]為直接計算式(1)提供了可能性。數(shù)值積分采用特定算法直接進行積分，可以避開經典算法中的展開計算，且輔以Matlab程序設計對算法進行實現(xiàn)，可大大減少計算工作量，得到高精度的計算結果。

數(shù)值積分算法有很多種，常用的主要有梯形算法、辛普森(Simpson)算法、復合梯形算法、復合辛普森算法、龍貝格(Romberg)算法、高斯-勒讓德(Gauss-Legendre)算法以及蒙特卡羅(Monte Carlo)算法等，但梯形算法與辛普森算法已被證明計算結果精度有限[3-4]，故本文主要采用后五種算法分別進行子午線弧長的計算，并對各種算法的結果予以比較。

2.1 復合梯形算法

(6)

復合梯形算法結果的準確度主要取決于步長h的大小，h越小，則計算結果越準確。

2.2 復合辛普森算法

類似復合梯形算法，若在每個子區(qū)間[xk，xk+1]上分別采用辛普森求積公式，即可得到復合辛普森算法計算公式

(7)

2.3 龍貝格算法

對于復合梯形算法和復合辛普森算法而言，雖然計算過程簡單，但收斂速度有限。龍貝格算法是在積分區(qū)間逐次分半的過程中，對用復合梯形算法和復合辛普森算法產生的近似值進行加權平均，不僅結果準確可靠，且迅速收斂。

龍貝格算法推導過程稍顯復雜，這里直接給出最終的計算公式，詳細推導過程請查閱文獻[3]。

(8)

2.4 高斯-勒讓德算法

n階勒讓德多項式定義為

其中xk為勒讓德多項式在[-1，1]區(qū)間上的零點。權系數(shù)Ak為

(9)

零點和權系數(shù)可以通過程序計算出來，也可以在很多計算手冊中查找得到。

2.5 蒙特卡羅算法

蒙特卡洛算法類似于機械求積方法，采用隨機數(shù)的方法產生結點，然后把各區(qū)間段上的函數(shù)值與區(qū)間段相乘后再相加。而實際上如果均勻分布的隨機點比較密集，則隨機數(shù)產生的區(qū)間近似均分整個積分區(qū)間段。具體方法如下：

(10)

3 數(shù)值積分計算結果及分析

現(xiàn)以IUGG—1975橢球的子午線弧長為例，分別采用上述五種數(shù)值積分算法進行計算。計算工具選擇數(shù)學計算軟件Matlab[4-5]，并根據(jù)各種算法原理及公式編寫相應的子程序。除了龍貝格算法的程序代碼為25行之外，其余四種算法的程序代碼都在15行之內，非常簡潔。同時在計算過程中還需要顧及到角度值與弧度值之間的轉換，最后得到計算結果列于表1，從表1中可對比得到各種數(shù)值積分算法結果和經典算法結果間的差值(以下簡稱為差值)。

在計算過程中，對于復合梯形算法，步長h分別選擇1°、1′和1″，可見當步長h=1°時，差值在2 m以內；當h=1′時，差值不超過1 mm；當h=1″時，差值為0 mm，說明復合梯形算法要想提高計算結果精度，需要減小步長。當步長h=1″時，才能得到與經典算法相同的結果。

表1 經典算法計算結果與各種數(shù)值積分算法計算結果

注：①以上數(shù)值的單位為m；②橢球元素選用IUGG-1975橢球元素；③限于篇幅，考慮到各種數(shù)值積分算法結果的千米量級數(shù)值和經典算法結果數(shù)值一樣，故在本表中略去了數(shù)值積分算法結果中大于千米的數(shù)據(jù)。

對于復合辛普森算法，通過計算發(fā)現(xiàn)，當步長h=1°時，差值為0 mm，即可得到與經典算法相同的結果，這說明復合辛普森算法要優(yōu)于復合梯形算法。

對于龍貝格算法，當限差Δ=10-3時(對應計算結果取到mm)，差值為0 mm，亦得到理想結果。

而高斯-勒讓德算法當取到5階時，差值在才1 mm以內，而取到3階或4階的結果均與經典算法差值較大，故而不再列于表中。

蒙特卡羅算法當選擇的n值為對應于h=1″換算的數(shù)值后，差值仍然為幾十米，說明此種算法在子午線弧長計算中不適用，其計算結果亦未列出。

4 結束語

基于Matlab軟件，分別對復合梯形算法、復合辛普森算法、龍貝格算法、高斯-勒讓德算法和蒙特卡洛算法在子午線弧長計算中的應用進行探討，通過與經典算法結果之間的對比，證明了數(shù)值積分算法的準確性和可靠性，說明文獻[10]在計算過程中有失誤之處。

基于計算結果及其分析比較，在利用數(shù)值積分進行子午線弧長計算時，復合辛普森算法和龍貝格算法要優(yōu)于其他算法，建議在測量工作中實際計算時予以采用。

[1] 孔祥元，郭際明，劉宗泉.大地測量學基礎[M].武漢：武漢大學出版社，2006

[2] 朱華統(tǒng).大地坐標系的建立[M].北京：測繪出版社，1986

[3] 李慶揚，王能超，毅大義.數(shù)值分析：第4版[M].北京：清華大學出版社，2001

[4] 宋葉志，賈東永.Matlab數(shù)值分析與應用[M].北京：機械工業(yè)出版社，2010

[5] 謝進，李大美.Matlab與計算方法實驗[M].武漢：武漢大學出版社，2009

[6] 王繼剛，崔旭升，張成.計算子午線弧長的常微分方程數(shù)值法[J].測繪通報，2012(9)：1-3

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

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

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

[10]劉修善.計算子午線弧長的數(shù)值積分法[J].測繪通報，2006(5)：4-6

TheCalculationandComparisonofMeridianArcLengthBasedonNumericalIntegralAlgorithm

ZHENG Hong-xiao1ZHANG Hong-fang1LEI Wei-wei2

2014-09-01

國家自然科學基金項目(41172199)；河南省科技創(chuàng)新人才項目(094100510023)

鄭紅曉(1981—)，女，2007年畢業(yè)于武漢大學地圖學與地理信息系統(tǒng)專業(yè)，理學碩士，工程師。

1672-7479(2014)06-0008-03

P22

： A