缔合Legendre函数的快速插值计算及其在六边形网格模型重力异常计算中的应用

doi:10.11947/j.AGCS.2024.20230110

摘要/Abstract

摘要：

鉴于局部区域非等纬度分布点的超高阶次球谐综合计算效率较低，本文深入研究了缔合Legendre函数的插值算法，结合模型重力异常求解的展开式，提出了超高阶模型空间重力异常插值快速计算方法，为了验证本文方法在东西狭长分布的点集更具有应用优势，采用球谐旋转变换技术进一步提升了计算效率。试验结果表明，利用2160阶次EGM2008模型计算日本地区同一高度的30 303个非等纬度分布的六边形网格点处模型重力异常，插值快速计算方法相比逐点计算，在误差不超±0.005 mGal水平下，计算耗时从3 669.41 s缩减到98.05 s，同时经过球谐旋转变换，将南北狭长的日本区域点集旋转为东西狭长分布，使得上述相应计算内容的耗时进一步由98.05 s缩减到19.06 s，计算效率相比最初方法提速比达到近200倍，有效解决了非等纬度分布点的超高阶模型重力异常快速解算的效率难题，且验证了该方法在东西狭长分布的情况下具有更高的提速比。

关键词: 地球重力场模型, 球谐综合, 勒让德函数, 六边形网格, 球谐旋转变换, 插值方法

Abstract:

In view of the low efficiency of the calculation of ultra-high-degree spherical harmonic synthesis for non-equal latitude distributed points of regional area, this article has conducted a deep research on the interpolation algorithm of associated Legendre functions. Combined with the harmonic expansion of the model gravity anomaly, a fast calculation method using interpolation was proposed for ultra-high-degree model free-air gravity anomaly. In order to verify the advantages of this new method over point sets with narrow distribution from east to west, the efficiency of the new method was further improved by adopting the spherical harmonic rotation(SHR) technology. The experimental results showed that the proposed method using interpolation reduced the calculation time from 3 669.41 s to 98.05 s compared with point-by-point approaches, with errors not exceeding ±0.005 mGal, when we used EGM2008 model up to 2160 degree and order to achieve the model free-air gravity anomalies of 30 303 hexagonal grid points at the same height level in Japan with non-equal latitude distribution. Meanwhile, by applying SHR and rotating these points distributed in the north-south direction to the east-west direction, the time consumption for the corresponding solution was further reduced from 98.05 s to 19.06 s, with a speed-up ratio of nearly 200 times to the original method. The method proposed in this article effectively solves the problem of low efficiency in solving ultra-high-degree model free-air gravity anomalies for regional non-equal latitude distributed points, and has a higher computing speed-up effect in the case of east-west elongated distribution.

Key words: Earth gravitational model, spherical harmonic synthesis, Legendre function, hexagonal grid, spherical harmonic transformation under rotation, interpolation method

中图分类号:

P223

李新星, 范昊鹏, 万宏发, 范雕, 冯进凯. 缔合Legendre函数的快速插值计算及其在六边形网格模型重力异常计算中的应用[J]. 测绘学报, 2024, 53(9): 1737-1747.

Xinxing LI, Haopeng FAN, Hongfa WAN, Diao FAN, Jinkai FENG. A fast method for interpolation of the associated Legendre functions and its application to the calculation of local hexagonal grid point gravity anomalies using an ultra-high-degree gravity field model[J]. Acta Geodaetica et Cartographica Sinica, 2024, 53(9): 1737-1747.

图/表 13

图1

图2

表1

图3

图4

表2

图5

表3

图6

图7

图8

表4

参考文献 35

[1]	LEMOINE F G, SMITH D E, KUNZ L, et al. The development of the NASA GSFC and NIMA joint geopotential model[C]//Proceedings of 1997 International Association of Geodesy Symposia. Berlin: Springer, 1997: 461-469.
[2]	PAVLIS N K, HOLMES S A, KENYON S C, et al. The development and evaluation of the earth gravitational model 2008 (EGM2008)[J]. Journal of Geophysical Research (Solid Earth), 2012, 117(B4):B04406.
[3]	ZINGERLE P, PAIL R, GRUBER T, et al. The combined global gravity field model XGM2019e[J]. Journal of Geodesy, 2020, 94(7):66.
[4]	李建成. 最新中国陆地数字高程基准模型：重力似大地水准面CNGG2011[J]. 测绘学报, 2012, 41(5):651-660,669.
	LI Jiancheng. The recent Chinese terrestrial digital height datum model: gravimetric quasi-geoid CNGG2011 [J]. Acta Geodaeticaet Cartographica Sinica, 2012, 41(5):651-660,669.
[5]	宁津生, 罗志才, 杨沾吉, 等. 深圳市1 km高分辨率厘米级高精度大地水准面的确定[J]. 测绘学报, 2003, 32(2):102-107.
	NING Jinsheng, LUO Zhicai, YANG Zhanji, et al. Determination of Shenzhen geoid with 1 km resolution and centimeter accuracy[J]. Acta Geodaetica et Cartographica Sinica, 2003, 32(2):102-107.
[6]	HIRT C, KUHN M. Evaluation of high-degree series expansions of the topographic potential to higher-order powers[J]. Journal of Geophysical Research: Solid Earth, 2012, 117(B12):2012JB009492.
[7]	HOLMES S A, FEATHERSTONE W E. A unified approach to the Clenshaw summation and the recursive computation of very high degree and order normalised associated Legendre functions[J]. Journal of Geodesy, 2002, 76(5):279-299.
[8]	XING Zhibin, LI Shanshan, TIAN Miao, et al. Numerical experiments on column-wise recurrence formula to compute fully normalized associated Legendre functions of ultra-high degree and order[J]. Journal of Geodesy, 2019, 94(1):2.
[9]	吴星, 张传定, 王凯. 卫星重力梯度边值问题的点质量调和分析[J]. 测绘学报, 2011, 40(2):213-219.
	WU Xing, ZHANG Chuanding, WANG Kai. Point-mass harmonic analysis of satellite gradiometry boundary value problem[J]. Acta Geodaetica et Cartographica Sinica, 2011, 40(2):213-219.
[10]	于锦海, 曾艳艳, 朱永超, 等. 超高阶次Legendre函数的跨阶数递推算法[J]. 地球物理学报, 2015, 58(3):748-755.
	YU Jinhai, ZENG Yanyan, ZHU Yongchao, et al. A recursion arithmetic formula for Legendre functions of ultra-high degree and order on every other degree[J]. Chinese Journal of Geophysics, 2015, 58(3):748-755.
[11]	吴星, 刘雁雨. 多种超高阶次缔合勒让德函数计算方法的比较[J]. 测绘科学技术学报, 2006, 23(3):188-191.
	WU Xing, LIU Yanyu. Comparison of computing methods of the ultra-high degree and order[J]. Journal of Geomatics Science and Technology, 2006, 23(3):188-191.
[12]	CLENSHAW C W. A note on the summation of Chebyshev series[J]. Mathematics of Computation, 1955, 9(51):118-120.
[13]	JEKELI C, LEE J K, KWON J H. On the computation and approximation of ultra-high-degree spherical harmonic series[J]. Journal of Geodesy, 2007, 81(9):603-615.
[14]	赵德军, 吴晓平. 基于Clenshaw求和法的重力场元计算[J]. 海洋测绘, 2004, 24(6):13-15, 26.
	ZHAO Dejun, WU Xiaoping. Computations of gravity quantities based on clenshaw summation method[J]. Hydrographic Surveying and Charting, 2004, 24(6):13-15, 26.
[15]	COLOMBO O L. Numerical methods for harmonic analysis on the sphere [R]. Columbus: The Ohio State University, 1981.
[16]	SNEEUW N, BUN R. Global spherical harmonic computation by two-dimensional Fourier methods[J]. Journal of Geodesy, 1996, 70(4):224-232.
[17]	RIZOS C. An efficient computer technique for the evaluation of geopotential from spherical harmonics [J]. Australian Journal of Geodesy and Photogrammetry, 1979, 31.
[18]	李新星, 李建成, 刘晓刚, 等. 球谐旋转变换结合非全次Legendre方法的局部六边形网格重力场球谐综合[J]. 地球物理学报, 2021, 64(11):3933-3947.
	LI Xinxing, LI Jiancheng, LIU Xiaogang, et al. Spherical harmonic synthesis of local hexagonal grid point gravity anomalies with non-full-order Legendre method combined with spherical harmonic rotation transformation[J]. Chinese Journal of Geophysics, 2021, 64(11):3933-3947.
[19]	MOAZEZI S, ZOMORRODIAN H, SIAHKOOHI H R, et al. Fast ultrahigh-degree global spherical harmonic synthesis on nonequispaced grid points at irregular surfaces[J]. Journal of Geodesy, 2016, 90(9):853-870.
[20]	KUNIS S, POTTS D. Fast spherical Fourier algorithms[J]. Journal of Computational and Applied Mathematics, 2003, 161(1):75-98.
[21]	BUCHA B, JANÁK J. A Matlab-based graphical user interface program for computing functionals of the geopotential up to ultra-high degrees and orders: efficient computation at irregular surfaces[J]. Computers & Geosciences, 2014, 66:219-227.
[22]	邢志斌, 李姗姗, 曲政豪, 等. 局部区域模型重力异常快速算法[J]. 测绘科学技术学报, 2016, 33(6):566-571.
	XING Zhibing, LI Shanshan, QU Zhenghao, et al. Fast method for model gravity anomalies calculation in local areas[J]. Journal of Geomatics Science and Technology, 2016, 33(6):566-571.
[23]	范昊鹏, 李姗姗. 局部区域模型重力异常快速算法研究[J]. 大地测量与地球动力学, 2013, 33(6):28-30, 35.
	FAN Haopeng, LI Shanshan. Study on a fast algorithm for model gravity anomalies in local areas[J]. Journal of Geodesy and Geodynamics, 2013, 33(6):28-30, 35.
[24]	JIANG Tao, XU Xinyu, CHU Yonghai, et al. Review of the research progress on static earth gravity field and vertical datum in China during 2019-2023[J]. Journal of Geodesy and Geoinformation Science, 2023, 6(3):76-86.
[25]	DANG Y. Preface to the special issue on the 2019—2023 China national report on geodesy[J]. Journal of Geodesy and Geoinformation Science, 2023, 6(3):3-4.
[26]	SEIF M R, SHARIFI M A, ESHAGH M. Polynomial approximation for fast generation of associated Legendre functions[J]. Acta Geodaetica et Geophysica, 2018, 53(2):275-293.
[27]	LI Xinxing, LI Jiancheng, TONG Xiaochong, et al. The employment of quasi-hexagonal grids in spherical harmonic analysis and synthesis for the earth's gravity field[J]. Journal of Geodesy, 2022, 96(11):89.
[28]	FUKUSHIMA T. Numerical computation of spherical harmonics of arbitrary degree and order by extending exponent of floating point numbers[J]. Journal of Geodesy, 2012, 86(4):271-285.
[29]	刘晓刚. GOCE卫星测量恢复地球重力场模型的理论与方法[D]. 郑州: 信息工程大学, 2011.
	LIU Xiaogang. Theory and method of recovering earth gravity field model by GOCE satellite measurement[D]. Zhengzhou: Information Engineering University, 2011.
[30]	HOFMANN-WELLENHOF B, MORITZ H. Physical geodesy [M]. New York: Springer, 2002.
[31]	贲进, 童晓冲, 张永生, 等. 球面等积六边形离散网格的生成算法及变形分析[J]. 地理与地理信息科学, 2006, 22(1):7-11.
	BEN Jin, TONG Xiaochong, ZHANG Yongsheng, et al. A generation algorithm of spherical equal-area hexagonal discrete gird and analysis of its deformation[J]. Geography and Geo-Information Science, 2006, 22(1):7-11.
[32]	STOUGH T, CRESSIE N, KANG E L, et al. Spatial analysis and visualization of global data on multi-resolution hexagonal grids[J]. Japanese Journal of Statistics and Data Science, 2020, 3(1):107-128.
[33]	许厚泽, 蒋福珍. 关于重力异常球函数展式的变换[J]. 测绘学报, 1964, 7(4):252-260.
	XU Houze, JIANG Fuzhen. On the transformation of the gravity anomaly expansion of spherical function[J]. Acta Geodaetica et Cartographica Sinica, 1964, 7(4):252-260.
[34]	金文华, 何涛, 刘晓平, 等. 基于有序简单多边形的平面点集凸包快速求取算法[J]. 计算机学报, 1998, 21(6):533-539.
	JIN Wenhua, HE Tao, LIU Xiaoping, et al. A fast convex hull algorithm of planar point set based on sorted simple polygon[J]. Chinese Journal of Computers, 1998, 21(6):533-539.
[35]	SHAMOS M. Computational geometry[D]. New Haven: Yale University, 1978.

最高阶次	行推方法	列推方法	跨阶次递推	Belikov递推	插值方法
2160	23.8	17.5	18.5	6.6	3.5
5400	135.3	115.0	115.2	68.7	21.6

重力异常点数	CPU耗时/s	重力异常值/mGal
重力异常点数	CPU耗时/s	最大值	最小值	平均值	标准差
30 303	3 669.41	177.48	－52.84	62.41	35.81

步长Δθ	节点数	耗时/s	模型重力异常计算误差/mGal
步长Δθ	节点数	耗时/s	最大值	最小值	平均值	标准差
10′	117	47.94	7.71×10^－1	－5.42×10^－1	1.65×10^－3	7.69×10^－2
5′	233	98.05	3.34×10^－3	－2.90×10^－3	－3.06×10^－6	3.14×10^－4
2.5′	465	180.85	8.45×10^－6	－7.32×10^－6	9.87×10^－10	7.41×10^－7
1′	1163	439.56	2.20×10^－9	－1.88×10^－9	6.84×10^－13	2.03×10^－10
30″	2327	884.37	1.27×10^－10	－1.32×10^－10	1.79×10^－14	9.28×10^－12

算法：旋转卡壳法
输入：凸多边形P={p¹,…,pⁿ}
输出：最短宽度、最短宽度的方向
	initial width[n]
	for i=1,2,…,n do
		j=i+1
		while Area(pⁱ,pⁱ⁺¹,p^j)<Area(pⁱ,pⁱ⁺¹,p^j⁺¹)do
			j=j+1
			if j>n then
				j=1
			end if
		end while
		width[i]=Width And Direction Compute(pⁱ,pⁱ⁺¹,p^j)
	end for
	return Width, Direction=Minimum (width)

步长Δθ	节点数	耗时/s	模型重力异常计算误差/mGal
步长Δθ	节点数	耗时/s	最大值	最小值	平均值	标准差
10′	21	11.13	3.57×10^－1	－1.62×10^－1	1.96×10^－3	1.88×10^－2
5′	41	19.06	3.61×10^－3	－4.10×10^－3	－6.53×10^－6	3.41×10^－4
2.5′	81	35.62	1.01×10^－5	－1.11×10^－5	－9.77×10^－10	8.20×10^－7
1′	199	85.74	2.85×10^－9	－2.44×10^－9	1.50×10^－12	2.26×10^－10
30″	395	167.13	1.33×10^－10	－1.40×10^－10	8.62×10^－13	1.08×10^－11