(武漢大學(xué)測繪學(xué)院，湖北武漢 430000)

1 子午線正算經(jīng)典算法

參考橢球具有對稱性，若要求解從赤道開始到任意緯度B的子午線弧長，只需求出積分[1，3]

(1)

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

為了求出M原函數(shù)，根據(jù)牛頓二項式將M展開為冪級數(shù)，然后代入式(1)中進行積分，即可得到結(jié)果。根據(jù)牛頓二項式對其進行級數(shù)展開，展至8次項得[3]

M=m0+m2sin2B+m4sin4B+m6sin6B+m8sin8B

(2)

式中

(3)

再將正弦的冪函數(shù)展開為余弦的倍數(shù)函數(shù)

(4)

將上式代入式(1)，得

M=a0-a2cos2B+a4cos4B-a6cos6B+a8cos8B

(5)

式中

將式(5)代入式(1)進行積分得

(6)

根據(jù)式(6)很容易編寫出計算機程序。

2 常微分方程數(shù)值算法求解子午線弧長

(7)

子午線弧長可看做有初值的常微分方程(7)在B處的近似解。

2.1 歐拉迭代算法

對于一階帶有初值的常微分方程

(8)

在xn處，采用泰勒級數(shù)展開

y(xn+1)=y(xn+h)

略去余項，有

y(xn+1)=y(xn)+y'(xn)h

(9)

yn+1=yn+hf(xn，yn) (n=0，1，2，…)

(10)

式(10)即為歐拉公式。

2.2 歐拉-梯形迭代算法

從式(10)中不易求得yn+1，還需要在區(qū)間[xn，xn+1]上對微分方程進行積分

(11)

將式(11)右端用梯形求積公式，有

f(xn+1，y(xn+1))]

(12)

對式(12)等號右端，用近似值yn代替y(xn)，yn+1代替y(xn+1)，可得

(13)

式(13)稱為梯形公式，將(10)和式(13)合用，構(gòu)成如下表達式

k=0，1，2，…;n=0，1，2，…

(14)

2.3 歐拉預(yù)估-矯正算法

實際上，當h很小時，讓式(14)中的梯形公式只迭代一次就結(jié)束，精度也滿足要求，該式稱為歐拉預(yù)估-矯正公式

k=0，1，2，…;n=0，1，2，…

(15)

2.4 龍格-庫塔算法

龍格-庫塔算法推導(dǎo)較為復(fù)雜，這里直接給出龍格-庫塔算法常用的兩種形式。

(1)二階龍格-庫塔算法

(16)

(2)三階龍格-庫塔算法

(17)

3 子午線弧長反算經(jīng)典算法

(18)

然后開始迭代，每次都讓

(19)

直到|Bi+1-Bi|<ε停止迭代，此時Bi+1即為所求的大地緯度。

4 常微分方程與數(shù)值迭代算法

根據(jù)公式(7)，可將子午線弧長與緯度看作一個帶有初值的常微分方程，將數(shù)值迭代算法應(yīng)用在這個常微分方程上，即可解得大地緯度B。常用的數(shù)值迭代算法有牛頓迭代、割線法以及單點迭代法，每一種迭代算法都可以與常微分方程數(shù)值解法結(jié)合使用。這里使用牛頓迭代法來進行討論。

(20)

(21)

式(21)中的f(Bn)可由上述四種常微分數(shù)值解法求解(X已知)。因此，每次迭代都可以根據(jù)常微分方程數(shù)值解法求得每次迭代后的f(Bn)，然后進行牛頓迭代，進而求得大地緯度B。

5 程序設(shè)計與結(jié)果分析

根據(jù)上述算法，使用C#實現(xiàn)上述算法并設(shè)計了程序界面[8，9]，操作界面如圖1所示。在此基礎(chǔ)上實現(xiàn)高斯正反算及數(shù)據(jù)檢驗。

圖1 程序主界面

通過選擇不同的算法，可得到相應(yīng)算法下的結(jié)果，同時，程序會給出與經(jīng)典算法的差值，如圖2、圖3所示。

圖2 子午線弧長正算算法選擇

以1975國際橢球為例，分別采用上述所介紹的數(shù)值積分、常微分方程數(shù)值解法和數(shù)值迭代方法進行計算，所得結(jié)果見表1[2]。

表1 子午線弧長正算(常微分方程數(shù)值解法)

注：(1)所得子午線弧長單位均為m；(2)由于所得結(jié)果和經(jīng)典算法均在米級以下，因此表格中的幾種數(shù)值積分算法所得結(jié)果省去了大于km的數(shù)值。

表2 子午線弧長反算(常微分與數(shù)值迭代解算)

注：(1)所得子午線弧長單位均為m；(2)由于所得結(jié)果和經(jīng)典算法只是在(")上不同，最后三列省去了度分值；(3)歐拉迭代和歐拉預(yù)估-校正公式試步長為1/1 000，二階龍格庫塔算法步長為1/100，四階龍格庫塔算法步長為1/10；(4)牛頓迭代次數(shù)為5次。

由表1可知，在子午線弧長正算中，步長1/1 000情況下的歐拉公式結(jié)果與經(jīng)典算法相同，而龍格-庫塔算法在迭代次數(shù)方面優(yōu)于歐拉公式和經(jīng)典算法。

在表2中，常微分方程數(shù)值解法所得結(jié)果與經(jīng)典算法結(jié)果基本一致，最大相差0.006 7″(基本可以忽略此差值)，并且牛頓迭代法與歐拉迭代算法相結(jié)合，彌補了歐拉公式精度低且步長小的缺點。

6 結(jié)束語

首先驗證了數(shù)值積分[1]和數(shù)值迭代[2]算法在子午線正反算中的正確性，并在此基礎(chǔ)上使用歐拉迭代、歐拉預(yù)估-矯正、龍格庫塔三種常見的常微分數(shù)值解法對子午線弧長進行正反算，并與傳統(tǒng)的子午線正反算結(jié)果進行比較。

基于公式推導(dǎo)及計算結(jié)果，常微分數(shù)值解法結(jié)果和傳統(tǒng)算法結(jié)果基本一致，并且具有實現(xiàn)簡單，迭代次數(shù)少、速度快等優(yōu)點。

[1] 鄭紅曉，張紅方，雷偉偉.子午線弧長計算的數(shù)值積分算法及其比較[J].鐵道勘察，2014，40(6):8-10

[2] 鄭紅曉，張紅方，雷偉偉.計算底點緯度Bf的數(shù)值迭代算法及其比較[J].測繪與空間地理信息，2015，38(2):42-44

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

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

[5] 李信真，車剛明，歐陽潔，等.計算方法[M].西安：西北工業(yè)大學(xué)出版社，2010

[6] 利慶揚，王能超，毅大義，等.數(shù)值分析[M].北京:清華大學(xué)出版社，2001

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

[8] JonSkeet.深入理解C#[M].姚琪琳，譯.北京：人民郵電出版社，2014

[9] 里克特.CLR via C#[M].周靖，譯.北京：清華大學(xué)出版社，2010

[10] 易維勇，邊少鋒，朱漢泉.子午線弧長的解析型冪級數(shù)確定[J].測繪學(xué)院學(xué)報，2000(3):167-171

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

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