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

歐陽明達，張敏利，于 亮

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

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

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

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

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

勒讓德函數；遞推公式；球諧分析

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

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

1 締合勒讓德函數計算公式

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

(1)

1.1 跨階數遞推公式

(2)

其中

(3)

即可以寫成

(4)

(5)

當m≥2時，有公式

(6)

其中

(7)

(8)

(9)

(10)

1.2Belikov列推公式

Belikov列推法引入新的非正常化球諧函數

(11)

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

(12)

(13)

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

(14)

其中

(15)

2 結果精度分析與比較

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

(16)

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

圖1 當θ=3°,30°,60°,87°時，跨階次遞推法Tn結果統計圖

圖2 當θ=3°,30°,60°,87°時，Belikov遞推法Tn結果統計圖

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

(17)

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

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

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

表1 計算用時比較

3 結 論

本文介紹了適用于超高階締合勒讓德函數計算的跨階數遞推算法和Belikov遞推算法，數值計算結果表明，兩種方法計算結果的相對精度均優于10-12，跨階數遞推算法在計算速度上優于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] 張傳定，許厚澤，吳星.地球重力場調和分析中的“輪胎”問題[C]//大地測量與地球動力學進展.武漢：科學技術出版社，2004：302-314.

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

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

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

[9] 王建強，趙國強，朱廣彬.常用超高階次締合勒讓德函數計算方法對比分析[J].大地測量與地球動力學，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.

[責任編輯：劉文霞]

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列推和跨階次遞推方法計算超高階締合勒讓德函數[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