兩種子午線弧長(zhǎng)計(jì)算方法的比較

1、兩種大地子午線弧長(zhǎng)計(jì)算方法的比較胡洋長(zhǎng)安大學(xué)西安2604070210摘要：在研究與大地橢球體有關(guān)的一些測(cè)量計(jì)算時(shí)，例如研究高斯投影計(jì)算和 弧度測(cè)量計(jì)算時(shí)，往往要用到大地子午線弧長(zhǎng)，而大地子午線弧長(zhǎng)的計(jì)算公式 涉及到橢圓積分，不能用普通方法求出被積函數(shù)。采用變步長(zhǎng)辛普森公式求它 的積分和按泰勒級(jí)數(shù)展開采用普通方法求出被積函數(shù)是兩種截然不同的計(jì)算方 法，本文就兩種不同的計(jì)算方法作出比較，用程序編程實(shí)現(xiàn)從而得到計(jì)算結(jié)果， 驗(yàn)證了兩種方法的正確性和可靠性以及按泰勒級(jí)數(shù)展開的方法求解大地子午線弧長(zhǎng)的優(yōu) 之處。關(guān)鍵字：大地子午線弧長(zhǎng)，變步長(zhǎng)辛普森公式，泰勒級(jí)數(shù)。引言：計(jì)算地球橢球子午線的弧長(zhǎng)，是大地測(cè)量、

2、天文測(cè)量、航空航天技術(shù)以及地理信息處理技術(shù) 中的一項(xiàng)基本內(nèi)容。計(jì)算子午線的弧長(zhǎng)涉及到橢圓積分，所以一般是將其展開成級(jí)數(shù)形式， 再用逐項(xiàng)積分的方法求出滿足一定精度要求的計(jì)算公式。本文在引出子午線弧長(zhǎng)的計(jì)算公式 之后，介紹了兩種計(jì)算方法的原理，并用C語(yǔ)言程序?qū)崿F(xiàn)了兩種算法的電算，繼而對(duì)兩種 算出的結(jié)果進(jìn)行了分析討論。一、子午線弧長(zhǎng)的計(jì)算公式我們知道，子午橢圓的一半，端點(diǎn)與極點(diǎn)相重合，而赤道又把子午線分成對(duì)稱的兩部分，因 此，推到從赤道開始到已知緯度B之間的子午線弧長(zhǎng)的計(jì)算公式即可。如下圖，今取子午線上的某微分弧段PP= d x，令P點(diǎn)緯度等于B, P點(diǎn)緯度為B+ dB, P 點(diǎn)的子午圈曲率半徑為M

3、，于是有dj MdB，因此，為了計(jì)算從赤道開始到任意緯度B的“ r Ba（1 一e2）平行圈之間的弧長(zhǎng)，必須求出下面的積分值X = j Md（式2）。又知M=（1 .）切（式3），故（式2）符合橢圓積分的定義（橢圓積分定義為可以表達(dá)為如下形式的任何函數(shù)x的積分：f （x） = j 0#t,P（t）dt，其中R是其兩個(gè)參數(shù)的有理函數(shù)，P是一個(gè)無(wú)重根的3或 4階多項(xiàng)式的平方根，而c是一個(gè)常數(shù)），因此（式1）為橢圓積分，無(wú)法用普通的積分方法 求出原函數(shù)。、兩種方法計(jì)算子午線弧長(zhǎng)的原理（一）變步長(zhǎng)辛普森公式求子午線弧長(zhǎng)的積分a (1 一 e 2)令M(B)=，計(jì)算步驟如下:V (1 - e2 sin2

4、 B)31、用梯形公式計(jì)算T = hM（0） + M（B）/2，其中n=1，h=B。且令Sn=Tnn2、3、1 _h n-1 一用變步長(zhǎng)梯形法則計(jì)算T2廣2 T2n + 2 M (B + h /2) k=0S = (4T -T )/3用辛普森公式計(jì)算2n 2 nn4、若IS - S l8，則令2 n n n，h/2 n h，轉(zhuǎn)2繼續(xù)進(jìn)行計(jì)算；否則結(jié)束即為所求的 2n n積分值。8是事先給定的積分精度。a (1 - e 2)，（二）泰勒級(jí)數(shù)展開后按普通方法積分令 M(B)=-V (1 - e2 sin2 B)31、按泰勒級(jí)數(shù)展開如下:M - m + m sin2 B + m sin4 B + m

5、 sin6 B + m sins B ；2、為便于積分將正弦的幕函數(shù)展開為余弦的倍數(shù)函數(shù)，由于sin 2 B =1-1 cos 2 B2 2 TOC o 1-5 h z .八311 fsin 4 B = 一一一 cos 2 B + cos 4 B 2sin 6 B = - cos 2 B + cos 4 B- cos 6 B16 321632.八 357 f 7 f 11sin 8 B = cos 2 B + cos 4 B - cos 6 B +cos 8B128 1632161283、將其代入上式，并經(jīng)過(guò)整理得到:M = a - a cos 2B + a cos 4B - a cos 6B

6、 + a cos8B 02468式中a = m + 2 + m + m +- m +.0028 4 16 6 128 8m157a =2 + m + m + m +.222 4 32 6 16 8137a = - m + m + m48 416 632 81284、將M = a - a cos 2B + a cos4B - a cos 6B + a cos8B代入到積分式中進(jìn)行積分， 02468整理后得到 X = a B 一 sin 2B + = sin 4B 一 與 sin 6B + 8 sin 8B 02468三、兩種方法計(jì)算子大地午線弧長(zhǎng)的數(shù)據(jù)和算例的討論筆者已將兩種方法的C語(yǔ)言程序?qū)懗?/p>

7、，限于篇幅本文未列出具體程序。該程序在任何PC機(jī) 上均可運(yùn)行。經(jīng)過(guò)用VC+6.0編寫的程序代碼調(diào)試，選用30度，45度，60度，90度三種 不同緯度的點(diǎn)來(lái)求相應(yīng)地點(diǎn)的大地子午線弧長(zhǎng)。采用的1954年北京坐標(biāo)系應(yīng)用的克拉索夫 斯基橢球。具體的計(jì)算結(jié)果如下表：用變步長(zhǎng)辛普森公式計(jì)算的結(jié)果X1按泰勒級(jí)數(shù)展開計(jì)算 的結(jié)果X2兩者之間的差值X1-X2緯度B3030對(duì)應(yīng)大地子午線 弧長(zhǎng)X3320172.406714743320172.40672017-0.00000543緯度B4545對(duì)應(yīng)大地子午線 弧長(zhǎng)X4985032.290475264985032.29047666-0.0000014緯度B6060對(duì)

8、應(yīng)大地子午線 弧長(zhǎng)X6654189.092215926654189.092214220.0000017緯度B9090對(duì)應(yīng)大地子午線 弧長(zhǎng)X10002137.497542910002137.49745320.000089從上表可以看出，兩種不同的數(shù)學(xué)方法計(jì)算大地子午線弧長(zhǎng)的精度相當(dāng)，差異只是在小數(shù)點(diǎn) 后邊的第五位。在實(shí)際的應(yīng)用過(guò)程當(dāng)中，我們只要求精確到0.1mm的數(shù)量級(jí)別就完全夠用 了。但是，我們從上表中也可以發(fā)現(xiàn)，兩種方法雖然在求得的結(jié)果上略有不同，但是發(fā)散增 長(zhǎng)的速度上差異還是比較大的。在中低緯度地區(qū)，由變步長(zhǎng)辛普森公式計(jì)算的結(jié)果要小于由 泰勒級(jí)數(shù)展開計(jì)算的結(jié)果；但到了高緯度地區(qū)，兩者計(jì)算的結(jié)果差變?yōu)檎龜?shù)。特別是在極點(diǎn) 處，兩者相差達(dá)到接近0.1mm。這說(shuō)明兩者之間的二階導(dǎo)數(shù)是步長(zhǎng)辛普森公式的方法大于泰勒級(jí)數(shù)展開的。兩者結(jié)果增長(zhǎng)的穩(wěn)定性上，按泰勒級(jí)數(shù)展開的方法求出的結(jié)果更加穩(wěn)定可 靠。雖然這兩種方法的精度略有不同，但是對(duì)