歐陽明達，張敏利，于 亮

(1.地理信息工程國家重點實驗室，陜西 西安 710054；2.西安測繪總站，陜西 西安 710054)

采用Belikov列推和跨階次遞推方法計算超高階締合勒讓德函數(shù)

歐陽明達1，2，張敏利2，于 亮2

(1.地理信息工程國家重點實驗室，陜西 西安 710054；2.西安測繪總站，陜西 西安 710054)

超高階球諧重力場模型的精確構(gòu)制與快速計算取決于締合勒讓德函數(shù)的計算方法。在前人研究的基礎(chǔ)上，文中對適合超高階締合勒讓德函數(shù)計算的Belikov列推和跨階次遞推方法進行介紹，為驗證精度，通過兩種途徑對計算結(jié)果進行檢驗，并比較其計算速度。結(jié)果表明，采用兩種算法得到的每個勒讓德函數(shù)的絕對精度均優(yōu)于10-12，在低階，跨階次遞推方法的計算用時大約是Belikov列推法的2倍，隨著階數(shù)的升高，跨階次遞推算法表現(xiàn)出明顯的速度優(yōu)勢。

勒讓德函數(shù)；遞推公式；球諧分析

隨著地球重力場模型的不斷精化，超高階次締合勒讓德函數(shù)的計算已經(jīng)成為了地球重力場和相關(guān)領(lǐng)域中的重要研究課題[1-8]。目前，遞推計算方法包括：標準前向列推法、標準前向行推法、跨階次遞推法和Belikov列推法。王建強等對4種方法的計算程序做了改進，從計算速度、計算精度兩方面分析比較了階次高至3 060階的各種方法的優(yōu)劣[9]；劉纘武等研究發(fā)現(xiàn)，超過2 000階次的締合勒讓德函數(shù)在接近兩極時達到極大的數(shù)量級，導(dǎo)致現(xiàn)有遞推方法在計算締合勒讓德函數(shù)值及其導(dǎo)數(shù)值時出現(xiàn)溢出，他提出通過插入壓縮因子，修改遞推算法，并結(jié)合使用Horner求和技術(shù)計算了球諧級數(shù)的部分和[7]; 于錦海等計算了超高階次勒讓德函數(shù)，對計算結(jié)果的精度進行了檢核，使得可以在雙精度數(shù)的范圍內(nèi)對任意階次的勒讓德函數(shù)進行計算[8]。

吳星等詳細介紹了現(xiàn)有多種締合勒讓德函數(shù)的遞推計算方法，通過數(shù)值試驗證明，Belikov列推法和跨階次遞推法是計算超高階次締合勒讓德函數(shù)較優(yōu)的方法[6]。在此研究的基礎(chǔ)上，本文采用兩種途徑對超高階締合勒讓德函數(shù)計算方法的結(jié)果精度進行了檢驗，并比較了其計算速度。

1 締合勒讓德函數(shù)計算公式

勒讓德函數(shù)，通常以符號pn(x)表示，n為階數(shù)，它滿足于勒讓德微分方程

(1)

1.1 跨階數(shù)遞推公式

(2)

其中

(3)

即可以寫成

(4)

(5)

當(dāng)m≥2時，有公式

(6)

其中

(7)

(8)

(9)

(10)

1.2Belikov列推公式

Belikov列推法引入新的非正?；蛑C函數(shù)

(11)

pnm(cosθ)是非正常化的勒讓德函數(shù)，可以得到其遞推關(guān)系為

(12)

(13)

由式(11)和式(13)可得

(14)

其中

(15)

2 結(jié)果精度分析與比較

引入檢核公式[10]，對其精度進行評估。

(16)

式(16)能夠用于評價某個特定階次勒讓德函數(shù)的計算精度，最高階數(shù)取N=3 000，考慮到θ值在接近于兩極附近時會出現(xiàn)溢出，令θ=3°,30°,60°,87°，圖1和圖2分別給出了采用跨階次遞推法和Belikov遞推法得到的Tn結(jié)果的統(tǒng)計圖?？芍?，兩種遞推方法的計算結(jié)果精度均優(yōu)于10-13，跨階次遞推法的相對精度值震蕩曲線較為劇烈，角度越小，振幅越大，在θ<30°左右的范圍，階數(shù)在N=500階內(nèi)，相對精度不斷提升，當(dāng)N>500階，相對精度呈下降趨勢，當(dāng)θ<30°時，階數(shù)越高，相對精度越低。隨著階數(shù)的升高，Belikov遞推法的相對精度不斷降低，且震蕩曲線較為平滑。

圖1 當(dāng)θ=3°,30°,60°,87°時，跨階次遞推法Tn結(jié)果統(tǒng)計圖

圖2 當(dāng)θ=3°,30°,60°,87°時，Belikov遞推法Tn結(jié)果統(tǒng)計圖

式(16)并不能用來評價同一階數(shù)內(nèi)勒讓德函數(shù)值平方和的精度，為此，引入另一個檢核公式[11]

(17)

圖3 當(dāng)N=90,360,1 000,1 600,2 160時，跨階次遞推法η(θ)的對數(shù)值曲線

圖4 當(dāng)N=90,360,1 000,1 600,2 160時，Belikov遞推法η(θ)的對數(shù)值曲線

利用締合勒讓德函數(shù)計算地球重力場模型要求在保證精度的情況下，盡可能實現(xiàn)快速計算。表1給出了在同一臺計算機上采用跨階數(shù)遞推法和Belikov遞推法計算不同階數(shù)締合勒讓德函數(shù)從θ=1°至θ=89°的時間，可以看出，在低階，跨階數(shù)遞推算法的計算時間大概是Belikov遞推算法的2～3倍，隨著階數(shù)升高，Belikov遞推算法的耗時性逐漸體現(xiàn)，原因在于本文采用Matlab平臺進行矩陣運算，隨著階數(shù)升高，矩陣之間的運算越加耗時。

表1 計算用時比較

3 結(jié) 論

本文介紹了適用于超高階締合勒讓德函數(shù)計算的跨階數(shù)遞推算法和Belikov遞推算法，數(shù)值計算結(jié)果表明，兩種方法計算結(jié)果的相對精度均優(yōu)于10-12，跨階數(shù)遞推算法在計算速度上優(yōu)于Belikov遞推算法，實用計算時建議使用跨階次遞推算法。

[1] COLMOBO O L.Numerical methods for harmonic analysis on the sphere[J]. OSU Rep 310.Columbus,Ohio:The Ohio state university，1981.

[2] WENZEL G. Ultra-high degree geopotential models GPM98A, B, and C to degree 1800.Paper to the joint meeting of the International Gravity Commission and International Geoid Commission[J],1998，7-12 September,Trieste.

[3] JEKELI C,LEE J K,KWON J H.On the computation and approximation of ultra-high-degree spherical harmonic series[J].J Geod.,2007，81(9):603-615.

[4] FUKUSHIMA T.Numerical computation of spherical harmonics of arbitrary degree and order by extending exponent of floating point numbers[J].J Geod.，2012，86(4)：271-285.

[5] 張傳定，許厚澤，吳星.地球重力場調(diào)和分析中的“輪胎”問題[C]//大地測量與地球動力學(xué)進展.武漢：科學(xué)技術(shù)出版社，2004：302-314.

[6] 吳星，劉雁雨.多種超高階締合勒讓德函數(shù)計算方法的比較[J].測繪科學(xué)技術(shù)學(xué)報,2006，23(3)：188-191.

[7] 劉纘武，劉世晗，黃歐.超高階次勒讓德函數(shù)遞推計算中的壓縮因子和Horner求和技術(shù)[J].測繪學(xué)報，2011，40(4)：454-458.

[8] 于錦海，曾艷艷，朱永超，等.超高階次勒讓德函數(shù)的跨階次遞推算法[J].地球物理學(xué)報，2015，58(3)：748-755.

[9] 王建強，趙國強，朱廣彬.常用超高階次締合勒讓德函數(shù)計算方法對比分析[J].大地測量與地球動力學(xué)，2009，4(2)：126-130.

[10] PAUL M R.Recurrence relations for integral of associatedlegendre functions[R]. Bulletin Geodesique,1978，53:177-190.

[11] HOLMES S A,FEATHERSTONE W E.A Unified approach to the clenshaw summation and the recursive computation of very high degree and order normalized associated legendre functions[J].J Geod.,2002，76(5):279-299.

[責(zé)任編輯：劉文霞]

Calculating the ultra-high-order associated legendre functions by Belikov column method and recursive method between every other order and degreeOU

YANG Mingda1，2，ZHANG Minli2,YU Liang2

(1.Stake Key Laboratory of Geo-information Engineering,Xi’an 710054,China; 2.Technical Division of Surveying and Mapping, Xi’an 710054,China)

Precision construction and rapid calculation of ultra-high-order spherical harmonic gravity field model,depend on the calculation method of the associated legendre functions. On the basis of previous studies,the suitable Belikov column method and recursion method between every other order and degree for ultra-high-order legendre function are introduced. The accuracy of calculations verified results in two ways after are their calculation speeds compared.The result shows that:every associated legendre function calculated by this two algorithms is obtained with absolute accuracy better than 10-12. In low-order, the recursion method between every other order and degree takes time as twice as Belikov column method extrapolation. As the order increasing, the recursion method between every other order and degree shows a significant speed advantage.

Legendre function; recursion function; spherical harmonic analysis

著錄：歐陽明達，張敏利，于亮.采用Belikov列推和跨階次遞推方法計算超高階締合勒讓德函數(shù)[J].測繪工程，2017，26(7)：12-15，21.

10.19349/j.cnki.issn1006-7949.2017.07.003

2016-06-24

地理信息工程國家重點實驗室開放研究基金資助項目(SKLGIE2015-M-1-2；SKLGIE2016-M-2-2)

歐陽明達(1986-)，男，助理工程師.

P223

A

1006-7949(2017)07-0012-04