鄭紅曉 張紅方 雷偉偉

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

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

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

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

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

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

子午線弧長(zhǎng)計(jì)算是橢球大地測(cè)量學(xué)的一項(xiàng)重要內(nèi)容，是高斯投影坐標(biāo)計(jì)算的基礎(chǔ)。但因子午線弧長(zhǎng)計(jì)算公式中的被積函數(shù)(即子午線曲率半徑M)無(wú)法找到其原函數(shù)，故而不能直接積分，經(jīng)典的解法是采用牛頓二項(xiàng)式定理進(jìn)行展開(kāi)，分項(xiàng)積分得到近似解析解[1-2]。

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

1 子午線弧長(zhǎng)計(jì)算的經(jīng)典算法

大地測(cè)量學(xué)經(jīng)典教材[1]給出了計(jì)算子午線弧長(zhǎng)的基本公式：從赤道到大地緯度為B處的子午線弧長(zhǎng)為

(1)

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

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

(2)

(3)

(4)

(5)

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

2 數(shù)值積分算法

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

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

2.1 復(fù)合梯形算法

(6)

復(fù)合梯形算法結(jié)果的準(zhǔn)確度主要取決于步長(zhǎng)h的大小，h越小，則計(jì)算結(jié)果越準(zhǔn)確。

2.2 復(fù)合辛普森算法

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

(7)

2.3 龍貝格算法

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

龍貝格算法推導(dǎo)過(guò)程稍顯復(fù)雜，這里直接給出最終的計(jì)算公式，詳細(xì)推導(dǎo)過(guò)程請(qǐng)查閱文獻(xiàn)[3]。

(8)

2.4 高斯-勒讓德算法

n階勒讓德多項(xiàng)式定義為

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

(9)

零點(diǎn)和權(quán)系數(shù)可以通過(guò)程序計(jì)算出來(lái)，也可以在很多計(jì)算手冊(cè)中查找得到。

2.5 蒙特卡羅算法

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

(10)

3 數(shù)值積分計(jì)算結(jié)果及分析

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

在計(jì)算過(guò)程中，對(duì)于復(fù)合梯形算法，步長(zhǎng)h分別選擇1°、1′和1″，可見(jiàn)當(dāng)步長(zhǎng)h=1°時(shí)，差值在2 m以內(nèi)；當(dāng)h=1′時(shí)，差值不超過(guò)1 mm；當(dāng)h=1″時(shí)，差值為0 mm，說(shuō)明復(fù)合梯形算法要想提高計(jì)算結(jié)果精度，需要減小步長(zhǎng)。當(dāng)步長(zhǎng)h=1″時(shí)，才能得到與經(jīng)典算法相同的結(jié)果。

表1 經(jīng)典算法計(jì)算結(jié)果與各種數(shù)值積分算法計(jì)算結(jié)果

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

對(duì)于復(fù)合辛普森算法，通過(guò)計(jì)算發(fā)現(xiàn)，當(dāng)步長(zhǎng)h=1°時(shí)，差值為0 mm，即可得到與經(jīng)典算法相同的結(jié)果，這說(shuō)明復(fù)合辛普森算法要優(yōu)于復(fù)合梯形算法。

對(duì)于龍貝格算法，當(dāng)限差Δ=10-3時(shí)(對(duì)應(yīng)計(jì)算結(jié)果取到mm)，差值為0 mm，亦得到理想結(jié)果。

而高斯-勒讓德算法當(dāng)取到5階時(shí)，差值在才1 mm以內(nèi)，而取到3階或4階的結(jié)果均與經(jīng)典算法差值較大，故而不再列于表中。

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

4 結(jié)束語(yǔ)

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

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

[1] 孔祥元，郭際明，劉宗泉.大地測(cè)量學(xué)基礎(chǔ)[M].武漢：武漢大學(xué)出版社，2006

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

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

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

[5] 謝進(jìn)，李大美.Matlab與計(jì)算方法實(shí)驗(yàn)[M].武漢：武漢大學(xué)出版社，2009

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

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

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

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

[10]劉修善.計(jì)算子午線弧長(zhǎng)的數(shù)值積分法[J].測(cè)繪通報(bào)，2006(5)：4-6

TheCalculationandComparisonofMeridianArcLengthBasedonNumericalIntegralAlgorithm

ZHENG Hong-xiao1ZHANG Hong-fang1LEI Wei-wei2

2014-09-01

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

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

1672-7479(2014)06-0008-03

P22

： A