大地测量与地球动力学

算法	180阶	360阶	720阶	1 800阶	2 160阶
FFT	0.025	0.730	4.910	62.810	111.910
LS	0.110	0.980	14.040	277.690	423.190
WLS	0.310	2.737	24.200	716.190	1 231.390
AQ	0.026	0.340	3.860	40.510	84.690
FNM	0.026	0.360	3.710	43.890	90.730

算法	180阶	360阶	720阶	1 800阶	2 160阶
FFT	0.025	0.730	4.910	62.810	111.910
LS	0.110	0.980	14.040	277.690	423.190
WLS	0.310	2.737	24.200	716.190	1 231.390
AQ	0.026	0.340	3.860	40.510	84.690
FNM	0.026	0.360	3.710	43.890	90.730

算法	min	max	mean	STD
DTM2006.0	-7 923.6	5 542.2	-1 890.7	2 629.2
FFT	-7 948.1	5 592.2	-1 891.9	2 637.0
LS	-7 923.6	5 542.2	-1 890.7	2 629.2
WLS	-7 936.2	5 534.9	-1 890.7	2 629.3
AQ	-7 940.2	5 537.1	-1 891.9	2 637.0
FNM	-7 939.7	5 537.9	-1 891.1	2 629.4

算法	min	max	mean	STD
DTM2006.0	-7 923.6	5 542.2	-1 890.7	2 629.2
FFT	-7 948.1	5 592.2	-1 891.9	2 637.0
LS	-7 923.6	5 542.2	-1 890.7	2 629.2
WLS	-7 936.2	5 534.9	-1 890.7	2 629.3
AQ	-7 940.2	5 537.1	-1 891.9	2 637.0
FNM	-7 939.7	5 537.9	-1 891.1	2 629.4

多种球谐分析与球谐综合算法研究及其软件设计

唐晓雨 ¹ , 陈文进 ¹^,* , 龚枭 ² , 谭小龙 ¹

大地测量与地球动力学 | 地球物理解释 2025,45(11): 1121-1126, 1184

收起

大地测量与地球动力学 |地球物理解释 2025 , 45 (11) : 1121 -1126, 1184

多种球谐分析与球谐综合算法研究及其软件设计

全屏

[Author(id=1322506546580411372, tenantId=1146029695717560320, journalId=1281212831689347082, articleId=1322506542939755478, orderNo=0, firstName=null, middleName=null, lastName=null, nameCn=null, orcid=null, stid=null, country=null, authorPic=null, dead=0, email=6720230182@mail.jxust.edu.cn, emailSecond=null, emailThird=null, correspondingAuthor=0, authorType=1, ext={EN=AuthorExt(id=1322506546702046190, tenantId=1146029695717560320, journalId=1281212831689347082, articleId=1322506542939755478, authorId=1322506546580411372, language=EN, stringName=Xiaoyu TANG, firstName=Xiaoyu, middleName=null, lastName=TANG, prefix=null, suffix=null, authorComment=null, nameInitials=null, affiliation=null, department=null, xref=¹, address=1 School of Civil and Surveying and Mapping Engineering, Jiangxi University of Science and Technology, 86 Hongqi Road, Ganzhou 341000, China, bio={"content":"

TANG Xiaoyu, postgraduate, majors in geodesy and surveying engineering, E-mail: 6720230182@mail.jxust.edu.cn

"}, bioImg=null, bioContent=

TANG Xiaoyu, postgraduate, majors in geodesy and surveying engineering, E-mail: 6720230182@mail.jxust.edu.cn

, aboutCorrespAuthor=null), CN=AuthorExt(id=1322506546773349359, tenantId=1146029695717560320, journalId=1281212831689347082, articleId=1322506542939755478, authorId=1322506546580411372, language=CN, stringName=唐晓雨, firstName=null, middleName=null, lastName=null, prefix=null, suffix=null, authorComment=null, nameInitials=null, affiliation=null, department=null, xref=¹, address=1 江西理工大学土木与测绘工程学院, 江西省赣州市红旗大道86号, 341000, bio={"content":"

唐晓雨, 硕士生, 主要从事大地测量学与测量工程研究, E-mail: 6720230182@mail.jxust.edu.cn

"}, bioImg=null, bioContent=

唐晓雨, 硕士生, 主要从事大地测量学与测量工程研究, E-mail: 6720230182@mail.jxust.edu.cn

, aboutCorrespAuthor=null)}, companyList=[AuthorCompany(id=1322506545154348005, tenantId=1146029695717560320, journalId=1281212831689347082, articleId=1322506542939755478, xref=null, ext=[AuthorCompanyExt(id=1322506545162736614, tenantId=1146029695717560320, journalId=1281212831689347082, articleId=1322506542939755478, companyId=1322506545154348005, language=EN, country=null, province=null, city=null, postcode=null, companyName=null, departmentName=null, remark=1 School of Civil and Surveying and Mapping Engineering, Jiangxi University of Science and Technology, 86 Hongqi Road, Ganzhou 341000, China), AuthorCompanyExt(id=1322506545171125223, tenantId=1146029695717560320, journalId=1281212831689347082, articleId=1322506542939755478, companyId=1322506545154348005, language=CN, country=null, province=null, city=null, postcode=null, companyName=null, departmentName=null, remark=1 江西理工大学土木与测绘工程学院, 江西省赣州市红旗大道86号, 341000)])]), Author(id=1322506546911761393, tenantId=1146029695717560320, journalId=1281212831689347082, articleId=1322506542939755478, orderNo=1, firstName=null, middleName=null, lastName=null, nameCn=null, orcid=null, stid=null, country=null, authorPic=null, dead=0, email=cwjwhu@whu.edu.cn, emailSecond=null, emailThird=null, correspondingAuthor=1, authorType=1, ext={EN=AuthorExt(id=1322506546987258867, tenantId=1146029695717560320, journalId=1281212831689347082, articleId=1322506542939755478, authorId=1322506546911761393, language=EN, stringName=Wenjin CHEN, firstName=Wenjin, middleName=null, lastName=CHEN, prefix=null, suffix=null, authorComment=null, nameInitials=null, affiliation=null, department=null, xref=¹^,*, address=1 School of Civil and Surveying and Mapping Engineering, Jiangxi University of Science and Technology, 86 Hongqi Road, Ganzhou 341000, China, bio=null, bioImg=null, bioContent=null, aboutCorrespAuthor=null), CN=AuthorExt(id=1322506547071144948, tenantId=1146029695717560320, journalId=1281212831689347082, articleId=1322506542939755478, authorId=1322506546911761393, language=CN, stringName=陈文进, firstName=null, middleName=null, lastName=null, prefix=null, suffix=null, authorComment=null, nameInitials=null, affiliation=null, department=null, xref=¹^,*, address=1 江西理工大学土木与测绘工程学院, 江西省赣州市红旗大道86号, 341000, bio=null, bioImg=null, bioContent=null, aboutCorrespAuthor=null)}, companyList=[AuthorCompany(id=1322506545154348005, tenantId=1146029695717560320, journalId=1281212831689347082, articleId=1322506542939755478, xref=null, ext=[AuthorCompanyExt(id=1322506545162736614, tenantId=1146029695717560320, journalId=1281212831689347082, articleId=1322506542939755478, companyId=1322506545154348005, language=EN, country=null, province=null, city=null, postcode=null, companyName=null, departmentName=null, remark=1 School of Civil and Surveying and Mapping Engineering, Jiangxi University of Science and Technology, 86 Hongqi Road, Ganzhou 341000, China), AuthorCompanyExt(id=1322506545171125223, tenantId=1146029695717560320, journalId=1281212831689347082, articleId=1322506542939755478, companyId=1322506545154348005, language=CN, country=null, province=null, city=null, postcode=null, companyName=null, departmentName=null, remark=1 江西理工大学土木与测绘工程学院, 江西省赣州市红旗大道86号, 341000)])]), Author(id=1322506547142448118, tenantId=1146029695717560320, journalId=1281212831689347082, articleId=1322506542939755478, orderNo=2, firstName=null, middleName=null, lastName=null, nameCn=null, orcid=null, stid=null, country=null, authorPic=null, dead=0, email=null, emailSecond=null, emailThird=null, correspondingAuthor=0, authorType=1, ext={EN=AuthorExt(id=1322506547234722808, tenantId=1146029695717560320, journalId=1281212831689347082, articleId=1322506542939755478, authorId=1322506547142448118, language=EN, stringName=Xiao GONG, firstName=Xiao, middleName=null, lastName=GONG, prefix=null, suffix=null, authorComment=null, nameInitials=null, affiliation=null, department=null, xref=², address=2 Jiangxi Science and Technology Infrastructure Platform Center, 1416 Jinggangshan Road, Nanchang 330003, China, bio=null, bioImg=null, bioContent=null, aboutCorrespAuthor=null), CN=AuthorExt(id=1322506547310220281, tenantId=1146029695717560320, journalId=1281212831689347082, articleId=1322506542939755478, authorId=1322506547142448118, language=CN, stringName=龚枭, firstName=null, middleName=null, lastName=null, prefix=null, suffix=null, authorComment=null, nameInitials=null, affiliation=null, department=null, xref=², address=2 江西省科技基础条件平台中心, 南昌市井冈山大道1416号, 330003, bio=null, bioImg=null, bioContent=null, aboutCorrespAuthor=null)}, companyList=[AuthorCompany(id=1322506546475553768, tenantId=1146029695717560320, journalId=1281212831689347082, articleId=1322506542939755478, xref=null, ext=[AuthorCompanyExt(id=1322506546483942377, tenantId=1146029695717560320, journalId=1281212831689347082, articleId=1322506542939755478, companyId=1322506546475553768, language=EN, country=null, province=null, city=null, postcode=null, companyName=null, departmentName=null, remark=2 Jiangxi Science and Technology Infrastructure Platform Center, 1416 Jinggangshan Road, Nanchang 330003, China), AuthorCompanyExt(id=1322506546492330986, tenantId=1146029695717560320, journalId=1281212831689347082, articleId=1322506542939755478, companyId=1322506546475553768, language=CN, country=null, province=null, city=null, postcode=null, companyName=null, departmentName=null, remark=2 江西省科技基础条件平台中心, 南昌市井冈山大道1416号, 330003)])]), Author(id=1322506547402494971, tenantId=1146029695717560320, journalId=1281212831689347082, articleId=1322506542939755478, orderNo=3, firstName=null, middleName=null, lastName=null, nameCn=null, orcid=null, stid=null, country=null, authorPic=null, dead=0, email=null, emailSecond=null, emailThird=null, correspondingAuthor=0, authorType=1, ext={EN=AuthorExt(id=1322506547494769661, tenantId=1146029695717560320, journalId=1281212831689347082, articleId=1322506542939755478, authorId=1322506547402494971, language=EN, stringName=Xiaolong TAN, firstName=Xiaolong, middleName=null, lastName=TAN, prefix=null, suffix=null, authorComment=null, nameInitials=null, affiliation=null, department=null, xref=¹, address=1 School of Civil and Surveying and Mapping Engineering, Jiangxi University of Science and Technology, 86 Hongqi Road, Ganzhou 341000, China, bio=null, bioImg=null, bioContent=null, aboutCorrespAuthor=null), CN=AuthorExt(id=1322506547566072830, tenantId=1146029695717560320, journalId=1281212831689347082, articleId=1322506542939755478, authorId=1322506547402494971, language=CN, stringName=谭小龙, firstName=null, middleName=null, lastName=null, prefix=null, suffix=null, authorComment=null, nameInitials=null, affiliation=null, department=null, xref=¹, address=1 江西理工大学土木与测绘工程学院, 江西省赣州市红旗大道86号, 341000, bio=null, bioImg=null, bioContent=null, aboutCorrespAuthor=null)}, companyList=[AuthorCompany(id=1322506545154348005, tenantId=1146029695717560320, journalId=1281212831689347082, articleId=1322506542939755478, xref=null, ext=[AuthorCompanyExt(id=1322506545162736614, tenantId=1146029695717560320, journalId=1281212831689347082, articleId=1322506542939755478, companyId=1322506545154348005, language=EN, country=null, province=null, city=null, postcode=null, companyName=null, departmentName=null, remark=1 School of Civil and Surveying and Mapping Engineering, Jiangxi University of Science and Technology, 86 Hongqi Road, Ganzhou 341000, China), AuthorCompanyExt(id=1322506545171125223, tenantId=1146029695717560320, journalId=1281212831689347082, articleId=1322506542939755478, companyId=1322506545154348005, language=CN, country=null, province=null, city=null, postcode=null, companyName=null, departmentName=null, remark=1 江西理工大学土木与测绘工程学院, 江西省赣州市红旗大道86号, 341000)])])]

唐晓雨¹, 陈文进¹^,*, 龚枭², 谭小龙¹

作者信息

1 江西理工大学土木与测绘工程学院, 江西省赣州市红旗大道86号, 341000

2 江西省科技基础条件平台中心, 南昌市井冈山大道1416号, 330003

通讯作者:

陈文进, 博士, 副教授, 主要从事物理大地测量学研究, E-mail: cwjwhu@whu.edu.cn

CHEN Wenjin, PhD, associate professor, majors in physical geodesy, E-mail: cwjwhu@whu.edu.cn

作者简介:

唐晓雨, 硕士生, 主要从事大地测量学与测量工程研究, E-mail: 6720230182@mail.jxust.edu.cn

Various Spherical Harmonic Analysis and Spherical Harmonic Synthesis Algorithms and Their Software Design

Xiaoyu TANG¹, Wenjin CHEN¹^,*, Xiao GONG², Xiaolong TAN¹

Affiliations

1 School of Civil and Surveying and Mapping Engineering, Jiangxi University of Science and Technology, 86 Hongqi Road, Ganzhou 341000, China

2 Jiangxi Science and Technology Infrastructure Platform Center, 1416 Jinggangshan Road, Nanchang 330003, China

TANG Xiaoyu, postgraduate, majors in geodesy and surveying engineering, E-mail: 6720230182@mail.jxust.edu.cn

出版时间: 2025-11-15 doi: 10.14075/j.jgg.2024.11.540

文章导航

摘要

收起

详细推导5种球谐分析算法, 并采用科学计算语言MATLAB开发一款界面友好且易于操作的实用软件。为验证算法性能, 以全球地形数据DTM2006.0为例, 对快速傅里叶变换(FFT)算法、最小二乘(LS)算法、加权最小二乘(WLS)算法、近似正交(AQ)算法、第一黎曼算法(FNM)在180、360、720、1 800和2 160阶的计算效率和精度进行数值分析。结果表明, AQ算法在计算效率方面表现最优异, 而LS算法计算精度最高; 5种算法的残差都近似服从正态分布, 表明5种算法均无系统性偏差。

关键词

球谐分析 / 球谐综合 / 地球重力场 / MATLAB / 软件设计

Abstract

收起

We derive five spherical harmonic analysis algorithms in detail and develop a user-friendly and easy-to-operate practical software using the scientific computing language MATLAB. To validate the performance of five algorithms, we conduct a numerical analysis of the computational efficiency and accuracy of fast Fourier transform(FFT), least squares(LS), weighted least squares(WLS), approximate quadrature(AQ), and first Neumann method(FNM) algorithms at harmonic degrees of 180, 360, 720, 1 800 and 2 160, using the global topographic data DTM2006.0 as a case study. The results show that the AQ algorithm performs best in terms of computational efficiency, while the LS algorithm achieves the highest precision. The residuals of all five algorithms approximately follow a normal distribution, indicating that none of five algorithms exhibit systematic biases.

Key words

spherical harmonic analysis / spherical harmonic synthesis / Earth's gravity field / MATLAB / software design

引用本文

唐晓雨, 陈文进, 龚枭, 谭小龙. 多种球谐分析与球谐综合算法研究及其软件设计. 大地测量与地球动力学, 2025 , 45 (11) : 1121 -1126, 1184 . DOI: 10.14075/j.jgg.2024.11.540

Xiaoyu TANG, Wenjin CHEN, Xiao GONG, Xiaolong TAN. Various Spherical Harmonic Analysis and Spherical Harmonic Synthesis Algorithms and Their Software Design[J]. Journal of Geodesy and Geodynamics, 2025 , 45 (11) : 1121 -1126, 1184 . DOI: 10.14075/j.jgg.2024.11.540

正文

收起

球谐函数是一组定义在球坐标系中的正交基函数，广泛用于描述球面上连续变化的现象^[1]。在地球科学中，地球的形状和特性通常通过球面上的物理场进行建模，球谐函数因此成为分析与综合球面数据的重要工具。随着地球科学研究的深入，球谐函数计算对分辨率、精度和性能的需求不断提升，推动了球谐分析与球谐综合算法及工具的进步。具体而言，在算法研究方面，球谐分析的主要算法包括最小二乘(LS)算法、快速傅里叶变换(FFT)算法以及Gauss-Legendre积分算法。Colombo^[2]对这3种算法均有所提及，尤其对FFT算法进行了详细阐述。Sneeuw^[3]正式提出基于LS方法的球谐分析方法，并推导出近似正交(AQ)算法。随后，Sneeuw^[4]提出第一黎曼算法(FNM)和第二黎曼算法(SNM)，两者均可实现离散勒让德函数的精确正交，但在纬度点选择和积分方法上有所不同，其中SNM因采用Gauss-Legendre积分计算正交权重，又称Gauss-Legendre积分算法。对于球谐分析与球谐综合的FFT算法，其在计算效率方面具有显著优势，因此成为学术界的研究热点，如Gruber等^[5]和Huang等^[6]都围绕该算法开展深入研究。鉴于现有研究成果已十分丰富，本文不再赘述。除算法研究外，为满足大规模数据处理和复杂计算的需求，国内外学者相继开发出多种用于球谐分析与球谐综合的软件与工具。Hwang等^[7]基于FFT和复数代数方法开发两款Fortran程序，用于球谐分析与球谐综合，支持水文和大气重力变化分析及卫星重力场数据处理。Wittwer等^[8]研究开发出基于扩展算术范围的球谐分析和球谐综合软件，但计算效率较低，且误差随纬度和阶数增加而增大。Wieczorek等^[9]开发出开源软件SHTools，支持球谐变换、数据旋转和局部多锥谱分析，最高支持2 800阶。

尽管球谐分析与球谐综合在相关算法及软件的研究上已取得诸多成果，但目前尚缺乏一款开源且用户友好的软件工具。因此，本文基于多种球谐分析与球谐综合算法，采用MATLAB语言开发一款GUI软件(软件下载地址：https://github.com/violet66688/Spherical_Harmonic_Analysis_And_Synthesis.git)，旨在提供直观友好的界面，使用户能够便捷地选择算法和调整参数，以满足多样化需求，为地球科学等领域的研究与应用提供技术支持。

1 连续观测数据下的球谐分析与球谐综合基本理论和算法

收起

单位球(0≤θ≤π以及0≤λ≤2π)面上任意平方可积的解析函数$f(\theta, \lambda)$可以展开为面球谐函数的级数^[2]，即：

(1)

$\begin{gather*}f(\theta, \lambda)= \\\sum\limits_{l=0}^{\infty} \sum\limits_{m=0}^{l} \bar{P}_{l m}(\cos \theta)\left(\bar{C}_{l m} \cos m \lambda+\bar{S}_{l m} \sin m \lambda\right) \end{gather*}$

式中，l和m分别为截断阶数和截断次数，且$l \geqslant m; \bar{P}_{l m}(\cos \theta)$为完全归一化的缔合勒让德函数；$\bar{C}_{l m}、\bar{S}_{l m}$为完全归一化的球谐系数。

根据球谐函数的正交性，由式(1)可以得到球谐分析公式为：

(2)

$\left.\begin{array}{l}\bar{C}_{l m} \\\bar{S}_{l m}\end{array}\right\}=\frac{1}{4 \mathsf{π} } {\displaystyle{\iint\limits_{\sigma}}} f(\theta, \lambda) \bar{P}_{l m}(\cos \theta)\left\{\begin{array}{c}\cos m \lambda \\\sin m \lambda\end{array}\right\} \mathrm{d} \sigma$

式中，$\mathrm{d} \sigma=\sin \theta \mathrm{d} \theta \mathrm{d} \lambda$式(1)所表达的过程为球谐综合，式(2)所表达的过程为球谐分析。

1.1 基于FFT的球谐分析与球谐综合算法

为方便球谐分析，定义下式^[10]：

(3)

$\begin{gather*}f_{l m}=\bar{S}_{l m}-\mathrm{i} \bar{C}_{l m}= \\\frac{1}{4 \mathsf{π} } \displaystyle{\int\limits_{\theta=0 \lambda=0}^{\mathsf{π} }} \displaystyle{\int\limits_{0}^{2 \mathsf{π} }} f(\theta, \lambda)(\sin m \lambda-\mathrm{i} \cos m \lambda) \bar{P}_{l m}(\cos \theta) \mathrm{d} \sigma= \\\frac{1}{4 \mathsf{π} } \displaystyle{\int\limits_{\theta=0}^{\mathsf{π} }} F_{\lambda}[f(\theta, \lambda)] \bar{P}_{l m}(\cos \theta) \sin \theta \mathrm{d} \theta \end{gather*}$

式中，i表示复数算子。式(3)便是基于FFT的球谐分析算法，计算可采用对每条纬线圈数据进行FFT，即可表示为F_λ[f(θ, λ)]，而沿经线圈的计算可以采用矩阵相乘的方法。

同样，根据式(1)发现：

(4)

$\begin{gathered}f(k \Delta \theta, j \Delta \lambda)=\sum\limits_{l=0}^{\bar{l}} \sum\limits_{m=0}^l\left(\bar{C}_{l m} \cos m j \Delta \lambda+\right. \\\left.\bar{S}_{l m} \sin m j \Delta \lambda\right) \bar{P}_{l m}(\cos k \Delta \theta)= \\2 \operatorname{Re}\left\{\sum\limits_{m=0}^{\bar{l}} \frac{1}{2}\left(\bar{C}_m-\mathrm{i} \bar{S}_m\right) \mathrm{e}^{\frac{2 m j k}{\bar{l}}}\right\}= \\2 \operatorname{Re}\left\{F_{\lambda}^{-1}\left\{\frac{\bar{C}_{m}-\mathrm{i} \bar{S}_{m}}{2}\right\}\right\}\end{gathered}$

式中，$\bar{l}$ 为最大的球谐展开阶次；k表示第k个余纬采样点；j表示第j个经度采样点。式(4)便是基于FFT的球谐综合算法，计算可采用沿纬线圈进行FFT逆变换，即可表示为$F_{\lambda}^{-1}\left\{\frac{\bar{C}_{m}-\mathrm{i} \bar{S}_{m}}{2}\right\}$。

1.2 基于LS的球谐分析算法

将球谐分析公式中的积分进行离散化并分两步进行，其中经度和纬度信息分别独立处理。将式(1)分两步进行计算，在经度方向进行等角离散^[11]：

(5)

$\begin{array}{c}\left. {\begin{array}{*{20}{c}}{{A_m}\left( {{\theta _i}} \right)}\\{{B_m}\left( {{\theta _i}} \right)}\end{array}} \right\} = \\\left\{ {\begin{array}{*{20}{c}}{\frac{1}{{l\left( {1 + {\delta _{m0}} - {\delta _{ml}}} \right)}}}\\{\frac{1}{{l\left( {1 - {\delta _{m0}} + {\delta _{ml}}} \right)}}}\end{array}} \right\}\sum\limits_{k = 0}^{2l - 1} f (\theta ,\lambda )\left\{ {\begin{array}{*{20}{l}}{\cos m\lambda }\\{\sin m\lambda }\end{array}} \right.\end{array}$

式中，$\delta_{m 0}$ 为克罗内克函数：

(6)

$\delta_{m 0}= \begin{cases}1 & m=0 \\ 0 & m \neq 0\end{cases}$

在进行第2步时，选择不同纬度会出现相应的离散解，这些解可以分为近似解和精确解。前者忽略了正交性，或为每个平行的纬线圈分配一个与阶次无关的正交权重；后者可以通过设计权重使其满足特定条件来实现精确正交。

采用矩阵矢量符号重写离散分析和综合分析公式。引入矩阵$\boldsymbol{a}=\left[A_{m}\left(\theta_{1}\right) A_{m}\left(\theta_{2}\right) \cdots\right. \left.A_{m}\left(\theta_{N}\right)\right]^{\mathrm{T}}、\boldsymbol{b}=\left[B_{m}\left(\theta_{1}\right) B_{m}\left(\theta_{2}\right) \cdots B_{m}\left(\theta_{N}\right)\right]^{\mathrm{T}}、\boldsymbol{c}= \left[\begin{array}{llll}\bar{C}_{m m} & \bar{C}_{m+1, m} & \cdots & \bar{C}_{l m}\end{array}\right]^{\mathrm{T}}、\boldsymbol{s}=\left[\begin{array}{llll}\bar{S}_{m m} & \bar{S}_{m+1, m} & \cdots & \bar{S}_{l m}\end{array}\right]^{\mathrm{T}}$以及维度为$N \times(l-m+1)$的缔合勒让德函数矩阵：

(7)

$\boldsymbol{P}=\left[\begin{array}{ccccc}\bar{P}_{m m}\left(\cos \theta_1\right) & \bar{P}_{m+1, m}\left(\cos \theta_1\right) & \bar{P}_{m+2, m}\left(\cos \theta_1\right) & \cdots & \bar{P}_{l m}\left(\cos \theta_1\right) \\\bar{P}_{m m}\left(\cos \theta_2\right) & \bar{P}_{m+1, m}\left(\cos \theta_2\right) & \bar{P}_{m+2, m}\left(\cos \theta_2\right) & \cdots & \bar{P}_{l m}\left(\cos \theta_2\right) \\\bar{P}_{m m}\left(\cos \theta_3\right) & \bar{P}_{m+1, m}\left(\cos \theta_2\right) & \bar{P}_{m+2, m}\left(\cos \theta_3\right) & \cdots & \bar{P}_{l m}\left(\cos \theta_3\right) \\\vdots & \vdots & \vdots & & \vdots \\\bar{P}_{m m}\left(\cos \theta_N\right) & \bar{P}_{m+1, m}\left(\cos \theta_N\right) & \bar{P}_{m+2, m}\left(\cos \theta_N\right) & \cdots & \bar{P}_{l m}\left(\cos \theta_N\right)\end{array}\right]$

根据式(2)可知，离散形式的球谐分析可以表示为：

(8)

$\left\{\begin{array}{l}\boldsymbol{c}=\boldsymbol{P}^{-1} \boldsymbol{a} \\\boldsymbol{s}=\boldsymbol{P}^{-1} \boldsymbol{b}\end{array}\right.$

当采用LS方法进行求解时，可以得到如下等式：

(9)

$\left\{\begin{array}{l}\boldsymbol{c}=\left(\boldsymbol{P}^{\mathrm{T}} \boldsymbol{P}\right)^{-1} \boldsymbol{P}^{\mathrm{T}} \boldsymbol{a} \\\boldsymbol{s}=\left(\boldsymbol{P}^{\mathrm{T}} \boldsymbol{P}\right)^{-1} \boldsymbol{P}^{\mathrm{T}} \boldsymbol{b}\end{array}\right.$

1.3 基于WLS的球谐分析算法

在得到式(9)时，可以为每条平行的纬线圈赋予相应的权重S_i，其值与该纬线圈的余纬正弦值成正比：

(10)

$\begin{equation*}S_{i}=2 \sin \theta_{n} / \sum\limits_{n=1}^{N} \sin \theta_{n} \end{equation*}$

式中，θ_n为各纬线圈的度数；N为纬线圈个数。将权重S_i组成对角矩阵 S，其对角线元素即为各权重S_i。因此，式(9)可以改写为：

(11)

$\left\{\begin{array}{l}\boldsymbol{c}=\left(\boldsymbol{P}^{\mathrm{T}} \boldsymbol{S} \boldsymbol{P}\right)^{-1} \boldsymbol{P}^{\mathrm{T}} \boldsymbol{S} \boldsymbol{a} \\\boldsymbol{s}=\left(\boldsymbol{P}^{\mathrm{T}} \boldsymbol{S} \boldsymbol{P}\right)^{-1} \boldsymbol{P}^{\mathrm{T}} \boldsymbol{S} \boldsymbol{b}\end{array}\right.$

1.4 基于AQ的球谐分析算法

球谐函数的离散正交性使得上述加权最小二乘法中的法方程矩阵$\boldsymbol{P}^{\mathrm{T}} \boldsymbol{S P}$接近对角矩阵，对角元素可用标量2(2－δ_m0)统一表示，因此，该法方程矩阵可以近似等同于一个缩放后的单位矩阵：

(12)

$\begin{equation*}\boldsymbol{P}^{\mathrm{T}} \boldsymbol{S} \boldsymbol{P} \approx 2\left(2-\delta_{m 0}\right) \boldsymbol{I} \end{equation*}$

公式(11)可以改写为：

(13)

$\left\{\begin{array}{l}\boldsymbol{c}=\frac{1}{4}\left(1+\delta_{m 0}\right) \boldsymbol{P}^{\mathrm{T}} \boldsymbol{S} \boldsymbol{a} \\\boldsymbol{s}=\frac{1}{4}\left(1+\delta_{m 0}\right) \boldsymbol{P}^{\mathrm{T}} \boldsymbol{S} \boldsymbol{b}\end{array}\right.$

这样就根据WLS方法得到另一种计算方法，即AQ算法。

1.5 基于FNM的球谐分析算法

到目前为止的3种方法，均忽略了缔合勒让德函数在纬度方向上的正交性。通过定义赋予每个纬线圈的权重W_i，可以在离散情况下保持缔合勒让德函数的正交性^[12]，即：

(14)

$\begin{equation*}\boldsymbol{P}^{\mathrm{T}} \boldsymbol{W} \boldsymbol{P}=2\left(2-\delta_{m 0}\right) \boldsymbol{I} \end{equation*}$

式中，W为权重矩阵，其对角线由权重W_i组成，而权重W_i可表示为：

(15)

$\sum\limits_{i=1}^{N} W_{i} \cos ^{n} \theta_{i}=\left\{\begin{array}{cc}\frac{2}{n+1}, & n \text { 为偶数 } \\0, & n \text { 为奇数 }\end{array}\right.$

式中，n=0，1，…，N-1。通过式(15)反解出正交权重W_i，并以对角线形式保存在对角权重矩阵 W中，则WLS方法可表示为：

(16)

$\left\{\begin{array}{l}\boldsymbol{c}=\left(\boldsymbol{P}^{\mathrm{T}} \boldsymbol{W} \boldsymbol{P}\right)^{-1} \boldsymbol{P}^{\mathrm{T}} \boldsymbol{W} \boldsymbol{a} \equiv \frac{1}{4}\left(1+\delta_{m 0}\right) \boldsymbol{P}^{\mathrm{T}} \boldsymbol{W} \boldsymbol{a} \\\boldsymbol{s}=\left(\boldsymbol{P}^{\mathrm{T}} \boldsymbol{W} \boldsymbol{P}\right)^{-1} \boldsymbol{P}^{\mathrm{T}} \boldsymbol{W} \boldsymbol{b} \equiv \frac{1}{4}\left(1+\delta_{m 0}\right) \boldsymbol{P}^{\mathrm{T}} \boldsymbol{W} \boldsymbol{b}\end{array}\right.$

以上便是用FNM来确定正交权重后进行球谐分析的具体过程。

2 软件设计与功能

收起

软件采用MATLAB编程语言在Windows系统下开发完成，软件系统结构与软件界面如图 1所示。该软件集成了多种球谐分析和球谐综合的核心算法，能够实现球谐系数的精确计算、球面函数值的还原、数据的可视化处理以及相关的数值统计分析等功能。其操作流程图如图 2所示。

3 软件验证

收起

为验证多种球谐分析与球谐综合算法的有效性及软件的可靠性，探讨各算法的优劣势和特点，本文采用DTM2006.0全球数字地形模型^[13]进行对比分析。以下所有实验均是在配置为Windows 11操作系统、Intel Core i9-14900HX CPU和32 GB内存的个人计算机上完成。

3.1 数值结果与分析

采用本文开发的球谐分析与球谐综合软件，分别对DTM2006.0数据进行球谐分析至180、360、720、1 800和2 160阶，得到相应的球谐系数图像(图 3为2 160阶下取完对数log₁₀后的球谐系数图像)以及计算时间(表 1，单位s)。

由图 3可知，在各算法生成的球谐系数中，即使经过对数尺度放大，系数值之间依然未见显著差异，表明算法在系数计算结果上具有一致性。从表 1可以看出，在低阶计算中，FFT、AQ和FNM三种算法的计算速度相差不大；但随着阶数增大，AQ算法展现出最高的计算效率，FNM算法的计算速度略低于AQ算法，而其他算法的计算速度则明显较慢。这得益于编写AQ和FNM算法时采用了式(13)和式(16)，在高阶次球谐分析中相较其他算法减少了大规模的矩阵相乘运算，从而大幅减少了计算时间。

3.2 残差展示与统计

球谐分析与球谐综合互为逆过程，因此通过得到的数据可以相互验证和校准。为此，本文利用§3.1中通过多种球谐分析方法得到的球谐系数数据进行球谐综合计算，以验证综合数据的完整性和合理性；然后对球谐综合结果进行残差分析，展示其空间分布特征及统计特性。基于FFT和数值积分的球谐综合方法的计算结果在空间分布上高度相似，仅计算速度不同。为避免重复，本文未单独列出FFT方法的计算结果。

图 4为5种算法计算得到的球谐系数经球谐综合后的数据空间分布，可以看出，数据分布合理，且未观察到局部异常值，验证了球谐分析与球谐综合过程的有效性。表 2(单位m)则表明5种算法得到的球谐系数经球谐综合后还原的数值偏差不大。图 5为各算法计算得到的球谐系数经球谐综合还原后，与原始数值之间的残差分布。由图可知，FFT和LS算法的残差分布均匀，且符合预期；而WLS算法由于按纬线圈分配权重，其残差呈现明显的纬线圈带状分布。相比之下，AQ和FNM算法的残差主要集中在高纬度区域，且靠近两极时残差值逐渐增大。残差空间分布结果表明，残差分布合理且数据重构过程中未引入显著的系统性误差。图 6显示5种方法计算的残差值都集中在0附近，且近似符合正态分布，进一步证明了综合结果与原始数据之间的差异是无偏且随机的。以上结果均可表明所采用的算法有效，且软件实现正确。

4 结语

收起

1) 本文详细推导5种球谐分析算法，并基于DTM2006.0数据，通过对比不同算法球谐分析的计算效率和球谐综合后的残差分布，论证各算法的特点与有效性。结果表明，在这5种算法中，LS算法精度最高，而AQ算法计算速度最快。此外，各算法均能较好且完整地实现球面函数值的分解与还原，残差符合正态分布，无显著系统性偏差，验证了该软件在实际应用中的可靠性和准确性。

2) 为满足不同应用领域的需求，开发出集多种球谐分析与球谐综合算法的软件。用户可以根据实际需求，灵活选择最适合的算法，平衡计算速度和精度。同时，软件设计简洁、实用，具备友好的用户界面，操作直观易懂，用户即使无专业背景，也能迅速上手并高效完成任务，可极大提高工作效率。

基金

收起

国家自然科学基金(42264001)
江西省重点研发计划(20243BBG71036)
江西理工大学博士启动基金(205200100588)
稀有金属资源安全高效开采江西省重点实验室开放基金(2024XJKCCSZ05)

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

Ditmar

, Klees

, Kostenko

. Fast and Accurate Computation of Spherical Harmonic Coefficients from Satellite Gravity Gradiometry Data[J]. Journal of Geodesy, 2003, 76, 690- 705

https://doi.org/10.1007/s00190-002-0298-x

Colombo O L. Numerical Methods for Harmonic Analysis on the Sphere[R]. The Ohio State University, Columbus, 1981

Sneeuw N. Discrete Spherical Harmonic Analysis: Neumann's Approach[M]//Montag H, Reigber C. Geodesy and Physics of the Earth. Berlin, Heidelberg: Springer, 1993

Sneeuw

. Global Spherical Harmonic Analysis by Least-Squares and Numerical Quadrature Methods in Historical Perspective[J]. Geophysical Journal International, 1994, 118 (3): 707- 716

https://doi.org/10.1111/j.1365-246X.1994.tb03995.x

Gruber

, Novák

, Sebera

. FFT-Based High-Performance Spherical Harmonic Transformation[J]. Studia Geophysica et Geodaetica, 2011, 55 (3): 489- 500

https://doi.org/10.1007/s11200-011-0029-y

Huang

Q H

, Zhang

, Fang

. Two-Step Spherical Harmonics ESPRIT-Type Algorithms and Performance Analysis[J]. IEEE/ACM Transactions on Audio, Speech, and Language Processing, 2018, 26 (9): 1684- 1697

https://doi.org/10.1109/TASLP.2018.2836436

Hwang

, Kao

Y C

. Spherical Harmonic Analysis and Synthesis Using FFT: Application to Temporal Gravity Variation[J]. Computers and Geosciences, 2006, 32 (4): 442- 451

https://doi.org/10.1016/j.cageo.2005.07.006

Wittwer

, Klees

, Seitz

, et al. Ultra-High Degree Spherical Harmonic Analysis and Synthesis Using Extended-Range Arithmetic[J]. Journal of Geodesy, 2008, 82 (4): 223- 229

Wieczorek

M A

, Meschede

. SHTools: Tools for Working with Spherical Harmonics[J]. Geochemistry, Geophysics, Geosystems, 2018, 19 (8): 2574- 2592

https://doi.org/10.1029/2018GC007529

Elowitz

, Hill

, Duvall

T L

. A Test of a Modified Algorithm for Computing Spherical Harmonic Coefficients Using an FFT[J]. Journal of Computational Physics, 1989, 80 (2): 506- 511

https://doi.org/10.1016/0021-9991(89)90115-0

Tsoulis D. Spherical Harmonic Computations with Topographic/Isostatic Coefficients[R]. Institute for Astronomical and Physical Geodesy, München, 1999

Tobler

. Global Spatial Analysis[J]. Computers, Environment and Urban Systems, 2002, 26 (6): 493- 500

https://doi.org/10.1016/S0198-9715(02)00017-0

Pavlis N K, Factor J K, Holmes S A. Terrain-Related Gravimetric Quantities Computed for the Next EGM[C]. The 1st International Symposium of the International Gravity Field Service (IGFS), Istanbul, 2007

2025年第45卷第11期

引用本文

BibTeX

文章信息

doi: 10.14075/j.jgg.2024.11.540

接收时间：2024-11-27
首发时间：2026-10-29
出版时间：2025-11-15

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2024-11-27

基金

National Natural Science Foundation of China(42264001)

国家自然科学基金(42264001)

Key Research and Development Program of Jiangxi Province(20243BBG71036)

江西省重点研发计划(20243BBG71036)

Foundation for Doctors of Jiangxi University of Science and Technology(205200100588)

江西理工大学博士启动基金(205200100588)

Open Fund of Jiangxi Key Laboratory of Safe and Efficient Mining of Rare Metal Resources(2024XJKCCSZ05)

稀有金属资源安全高效开采江西省重点实验室开放基金(2024XJKCCSZ05)

作者信息

1 江西理工大学土木与测绘工程学院, 江西省赣州市红旗大道86号, 341000

2 江西省科技基础条件平台中心, 南昌市井冈山大道1416号, 330003

通讯作者:

陈文进, 博士, 副教授, 主要从事物理大地测量学研究, E-mail: cwjwhu@whu.edu.cn

CHEN Wenjin, PhD, associate professor, majors in physical geodesy, E-mail: cwjwhu@whu.edu.cn

参考文献

分享链接

https://castjournals.cast.org.cn/joweb/ddcl/CN/10.14075/j.jgg.2024.11.540

分享至

全文二维码

扫描看全文

引用本文

BibTeX

本文的引用情况

2种不同金属材料的力学参数

科 Family	属数 Number of genus	种数 Number of species	占总种数比例 Percentage of total species (%)	属 Genus	种数 Number of species	占总种数比例 Percentage of total species (%)
鹅膏菌科Amanitaceae	2	11	5.26	鹅膏菌属 Amanita	10	4.78
小菇科 Mycenaceae	2	12	5.74	丝盖伞属 Inocybe	5	2.39
多孔菌科 Polyporaceae	8	14	6.70	蜡蘑属 Laccaria	5	2.39
红菇科 Russulaceae	3	23	11.00	小皮伞属 Marasmius	6	2.87
				小菇属 Mycena	11	5.26
				光柄菇属 Pluteus	5	2.39
				红菇属 Russula	17	8.13
				栓菌属 Trametes	5	2.39

关闭全屏

BibTeX
EndNote
RefWorks
TxT