基于四階Runge-Kutta方法的大地主題正算*

朱寧寧，張趙興，梁曉莉，魏祖帥

(河南理工大學(xué)礦山空間信息技術(shù)國(guó)家測(cè)繪地理信息局重點(diǎn)實(shí)驗(yàn)室，河南焦作454003)

橢球面上點(diǎn)的大地經(jīng)度、大地緯度，兩點(diǎn)間的大地線長(zhǎng)度及其正、反大地方位角，統(tǒng)稱為大地元素。通過(guò)已知的某些大地元素推求另一些大地元素，這樣的計(jì)算稱為大地主題解算。大地主題解算是大地主題正算和反算的統(tǒng)稱，是橢球面大地測(cè)量學(xué)研究的核心問(wèn)題之一［1］。

大地主題解算常用的一種方法是運(yùn)用勒讓德級(jí)數(shù)將它們展開(kāi)為大地線長(zhǎng)度的升冪級(jí)數(shù)，再逐項(xiàng)計(jì)算以達(dá)到主題解算的目的。但是這種級(jí)數(shù)式收斂緩慢，項(xiàng)目繁多，計(jì)算起來(lái)非常復(fù)雜，雖然經(jīng)過(guò)不同學(xué)者的改進(jìn)，可以達(dá)到精度很高的解算值，但該方法涉及的數(shù)學(xué)技巧多，要求掌握較多的橢球大地測(cè)量和空間解析幾何知識(shí)［1-2］。計(jì)算機(jī)的發(fā)展和數(shù)值計(jì)算算法的提高，為一些原本封閉的方程和偏微分方程等在傳統(tǒng)手工方法下難以解算的問(wèn)題，提供了數(shù)值方法的解算途徑［3-4］。這些問(wèn)題的數(shù)值解可以人為控制精度，并以很快的速度進(jìn)行解算。對(duì)于數(shù)值計(jì)算方法，目前文獻(xiàn)主要討論的是數(shù)值積分法，而在此類問(wèn)題中，應(yīng)用微分方程的數(shù)值解法進(jìn)行求解在目前的文獻(xiàn)中尚不多見(jiàn)。本文提出一種用四階Runge-Kutta方法求解大地線微分方程的數(shù)值解法，該方法通過(guò)橢球面上各大地元素之間的嚴(yán)密函數(shù)關(guān)系，導(dǎo)出了橢球面上的大地線微分方程組計(jì)算模型，并用四階Runge-Kutta方法求解經(jīng)過(guò)組合后的微分方程組，完成了大地主題正算。

1 四階Runge-Kutta方法的分析及模型

1.1 龍格-庫(kù)塔(Runge-Kutta)法解大地線微分方程分析

龍格-庫(kù)塔(Runge-Kutta)法解大地線微分方程的實(shí)質(zhì)是用若干點(diǎn)函數(shù)值的線性組合來(lái)代替泰勒級(jí)數(shù)展開(kāi)式中的導(dǎo)數(shù)計(jì)算。四階Runge-Kutta方法屬于數(shù)值積分算法，它既不采用勒讓德級(jí)數(shù)也不用借用輔助面。使用這種方法可以達(dá)到較高的解算精度，可用于解算150 km以內(nèi)的大地主題正解。該算法簡(jiǎn)單、易于編程，并且可以達(dá)到和經(jīng)典的高斯平均引數(shù)法相當(dāng)?shù)挠?jì)算量。

1.2 大地線微分方程

橢球參數(shù)計(jì)算示意圖如圖1所示，設(shè)P為大地線上任意一點(diǎn)，其經(jīng)度為L(zhǎng)，緯度為B，大地線方位角為A。當(dāng)大地線增加dS到點(diǎn)P1時(shí)，則上述各量相應(yīng)變化dL、dB及dA。大地線的微分方程可表達(dá)為dL、dB、dA與dS的關(guān)系式。由圖1可知，dS在子午圈上的分量P2P1=M·dB，在平行圈上的分量PP2=N·cosB·dL。三角形PP2P1是個(gè)一微分直角三角形，它結(jié)合球面直角三角形P1P3N，可求得大地線微分方程組中dL、dB、dA與dS的關(guān)系式:

將式(1)化為dL、dB、dA與dS的關(guān)系式:

式(2)稱為大地線微分方程，是本文所用四階龍格 －庫(kù)塔(Runge-Kutta)法的基本函數(shù)模型。

圖1 橢球參數(shù)計(jì)算示意圖

1.3 四階Runge-Kutta模型

Runge-Kutta方法實(shí)質(zhì)上是間接地使用Taylor級(jí)數(shù)法的一種方法。考慮差商，根據(jù)微分中值定理，存在0＜θ＜1，使得

于是，利用所給方程y'=f(x，y)得到

設(shè)K*=f(xn+θ·h，y(xn+θ·h))，稱K*為區(qū)間［xn，xn+1］上的平均斜率。由此可見(jiàn)，精確求解微分方程數(shù)值解的關(guān)鍵是找到一個(gè)適當(dāng)?shù)钠骄甭?。采用四階Runge-Kutta方法確定的平均斜率，可使求解值具有四階精度，即局部截?cái)嗾`差是O(h5)。具體做法是在區(qū)間［xn，xn+1］上用四個(gè)點(diǎn)處的斜率加權(quán)作為平均斜率的近似值。四階Runge-Kutta具有多種格式，本文采用經(jīng)典格式模型［5］。

1.4 四階Runge-Kutta法求解大地線微分方程

大地線微分方程中分別包含dL、dB、dA與dS的關(guān)系式，各個(gè)變量之間具有相關(guān)性，經(jīng)典格式只包含單一方程，為將其用于解算大地元素，需將經(jīng)典格式改寫(xiě)為方程組形式，得到求解模型。

大地主題正解微分方程組可表示為:

改進(jìn)的四階Runge-Kutta經(jīng)典形式為:

大地主題正算的初始條件為 B0、L0、A0、S，則有

式中:h為等分后步長(zhǎng);si為選取的大地線的第i個(gè)等分節(jié)點(diǎn)到起始點(diǎn)的大地線長(zhǎng)度。

由此可以看出，用改進(jìn)的經(jīng)典四階Runge-Kutta法進(jìn)行大地主題的正算，公式簡(jiǎn)潔、規(guī)律明顯、易于理解，而且易于編程實(shí)現(xiàn)。

2 算例分析

本文使用C++編程，實(shí)現(xiàn)上述算法，在克拉索夫斯基橢球上進(jìn)行計(jì)算，橢球的參數(shù)為長(zhǎng)半軸a=6 378 245.0 m，橢球第一偏心率為e，e2=0.006 693 421 622 966，將經(jīng)典算法的解算值作為參考值，比較本文算法與經(jīng)典算法的差別。表1的四個(gè)例子中大地線S等分的個(gè)數(shù)分別選為10、10、100、100.

表1 大地主題正解解算結(jié)果 °'″

續(xù)表

表1算例包含了短距離(例1)，中距離(例2)，長(zhǎng)距離(例3、4)等情況，最大的偏差在例3中B2-B2’的解算中，有亞秒級(jí)的差別，其余的差別大多在毫秒級(jí)別上，分析偏差的原因，主要有兩項(xiàng)，一是距離的影響，二是等分節(jié)點(diǎn)個(gè)數(shù)的選擇。當(dāng)距離一定時(shí)，合理的選擇節(jié)點(diǎn)可以很大程度上提高精度。因此，若需提高精度，簡(jiǎn)單的方法是僅增加等分節(jié)點(diǎn)的個(gè)數(shù)，不需要改變程序的其它部分。這與經(jīng)典的基于級(jí)數(shù)展開(kāi)的算法有很大不同，經(jīng)典算法如果要提高精度，需要增加多項(xiàng)式展開(kāi)的次數(shù)，這會(huì)導(dǎo)致求解公式變得復(fù)雜，程序也要作較大改動(dòng)。

圖2 解算結(jié)果比較圖

3 結(jié)論

本文提出的大地主題正解的數(shù)值解法，采用嚴(yán)密公式進(jìn)行計(jì)算，而且公式表達(dá)簡(jiǎn)單，特別適合計(jì)算機(jī)計(jì)算。本文對(duì)正解的微分方程采用單步的龍格－庫(kù)塔方法進(jìn)行解算。理論上，如果采用多步法，如Adams方法，那么對(duì)于作為光滑曲線的大地線來(lái)說(shuō)將可以得到任意精度的解算結(jié)果。

［1］ 孔祥元，郭際明，劉宗全.大地測(cè)量學(xué)基礎(chǔ)［M］.武漢:武漢大學(xué)出版社，2002.

［2］ 周振宇，郭廣禮，賈新果.大地主題解算方法綜述［J］.測(cè)繪科學(xué)，2007(4):190－191.

［3］ 王建強(qiáng)，胡明慶.貝賽爾大地主題解算分析［J］.測(cè)繪科學(xué)，2012(1):30－31.

［4］ 范業(yè)明，王解先，劉慧芹.大地主題的數(shù)值解法［J］.工程勘察，2007，38(1):6－9.

［5］ 李慶揚(yáng)，王能超，易大義.數(shù)值分析［M］.武漢:華中科技大學(xué)出版社，2006.