扰动重力梯度的球冠谐分析建模

王 燚，姜效典

中国海洋大学海洋地球科学学院，山东 青岛 266100

扰动重力梯度的球冠谐分析建模

王 燚，姜效典

中国海洋大学海洋地球科学学院，山东 青岛 266100

从球冠谐理论出发，详细推导了球冠坐标系下扰动重力梯度的无奇异性计算公式。基于Tikhonov正则化方法，利用GOCE卫星实际观测数据解算局部重力场球冠谐模型。数值计算表明，基于扰动重力梯度的球冠谐分析建模方法能够有效地恢复局部重力场中的短波信号，与GO_CONS_GCF_2_DIR_R5模型的差异在±0.3×10-5m/s2水平。

球冠谐和分析; 扰动重力; 重力梯度; Tikhonov正则化方法

随着CHAMP、GRACE和GOCE等3颗重力卫星的发射升空，利用卫星重力观测值恢复地球重力场的研究，已经成为当前地球物理学研究的热点问题之一。文献[1—2]的研究表明，重力梯度张量反映了重力位水准面的曲率和力线弯曲，对重力场中的短波信号变化敏感，所以重力梯度张量更能够反映重力场的精细结构，包含有大量的局部重力场的有效信息。

国内外众多学者就利用卫星重力观测值恢复地球重力场的方法做了大量的研究工作。文献[3]提出了广义轮胎调和分析方法，利用卫星重力梯度张量分量组合以较高的精度还原地球重力场。文献[4]提出了使用最小二乘法恢复地球重力场的方法，通过数值模拟，证明重力异常的梯度对轨道高度的变化不敏感。在文献[5]中研究了全张量重力梯度数据的全局和局部分量的广义球谐谱表示和轨道根数表示，给出了广义球谐函数与球谐函数之间的关系，从理论上得到了全张量重力梯度数据的描述方法和由全张量重力梯度网格数据恢复全球重力位谱系数的基本公式。然而这些研究工作在理论方面都是利用传统的球谐分析法。在表达全球重力场时，球谐分析法有明显的优势，但是对于局部区域重力场的精细结构和重力场高频信息的表述上，球谐分析法存在着不足和缺陷。例如球谐函数不能满足局部重力场的边界条件，在表述高分辨率重力场时，需要求解巨量的球谐系数等[6]。因此需要针对卫星重力梯度恢复局部重力场的理论方法进行系统的研究。

文献[7]在前人的研究基础上总结提出了球冠谐分析理论，对地球北极区域的磁异常做了数值分析。文献[8]利用航磁异常观测值，解算了中国大陆航磁异常的不同阶次球冠谐模型，结合大地构造理论对相关图件进行了分析解释。文献[9]研究了球冠谐函数在重力场中的表达式，并利用球冠谐理论建立了我国的局部重力场模型。文献[10]提出了一种基于点质量模型的多层球冠谐系数迭代构造算法，并对求解矩阵A的病态性和球冠半角与点质量埋藏深度和密度进行了详细的分析研究，通过数值模拟表明当迭代算法收敛时，球冠谐模型能够以高精度拟合重力观测值。文献[11]提出了使用球冠谐模型的重力异常径向导数和曲面导数建立曲面位场的方法。这些研究工作都取得了令人满意的结果。然而利用扰动重力梯度恢复局部球冠谐重力场的研究，到目前为止还不多见。由此，本文从球冠谐理论出发，推导了无δ奇异性球冠谐梯度计算公式，使用GOCE卫星数据产品中地球固定参考框架(terrestrial reference frame，TRF)下的重力梯度数据恢复了4个球冠谐重力场模型，最后在计算点上，与文献[12]中使用GOCE卫星梯度数据建立的GO_CONS_GCF_2_DIR_R5地球重力场模型的扰动重力数值进行了比对分析，检验该方法的正确性。

1 球冠谐和分析扰动位重力梯度模型

球冠谐和分析是一种在球冠坐标下的满足位理论边值条件的谱函数法，最早由加拿大人Haines在研究与地壳相关的局部区域磁异常时提出。通过利用非整阶勒让德函数代替整阶勒让德函数，余纬θ的定义域由[0,π]变换到[0，α](α是球冠半角)，同时将地球Z轴由北极点沿子午线旋转到球冠中心，以保证球冠上点的经度依然满足周期性的边界条件。采用分离变量法可以得到重力场扰动位在球冠坐标下的解[13]

(1)

(2)

(3)

系数Aj(m,l)的递推关系如下

(4)

定义

(5)

则由文献[7]推导的近似式可得

(6)

由上述推导过程可知

(7)

如果令

为球冠谐系数，则式(1)可以改写为

(8)

式中，G是万有引力常数；M是地球质量。因为δ与λc是球冠坐标下的余纬和经度(参见图2)，所以将球冠中心Q看作是新的北极点，定义球冠坐标系下的球冠局部指北坐标系(spherical cap local north-oriented reference frame，SCLNOF)的x轴指向球冠极点，z轴平行地心向径，方向向外，y轴与x和z轴构成右手坐标系，则与球谐函数的局部指北坐标系(local north-oriented reference frame，LNOF )下重力梯度张量计算一样，将式(8)代入局部指北坐标系梯度计算公式[16]，可以得到球冠局部指北坐标下的扰动重力和扰动重力梯度张量模型公式

(9)

(10)

式中

(11)

从式(9)、(10)中可以看出，在球冠极点附近，当sinδ趋于零时，存在δ奇异性，因此根据非整阶勒让德函数的计算公式，本文推导了如下的去δ奇异性球冠谐重力梯度张量计算公式。记

(12)

(13)

(14)

(15)

将式(15)代入式(10)，就可以得到非奇异的球冠局部指北坐标下扰动重力梯度张量的计算公式。根据球冠谐分析理论，整阶勒让德函数是非整阶勒让德函数的特例[17-18]，因此当球冠半角α=π/2时，式(15)中扰动位梯度张量的各个阶数l为整数，这时非整阶勒让德函数的计算结果与整阶勒让德函数的计算结果应当一致。下面取α=π/2，δ∈[0,π/2]，通过式(15)计算各个梯度中l=2,3,m=0,1,2,3对应的非整阶勒让德函数项，并与文献[19]中的扰动重力梯度张量公式中n=2,3,m=0,1,2,3整阶勒让德函数项计算结果进行比较，统计结果如表1所示，二者差值的均方根小于10-3数量级。从表1可见，当m≠0且m≠n时，二者误差较大，这是由于非整阶勒让德函数值使用近似式(6)计算，而当m=0或m=n时，式(6)不可用，必须用式(5)计算，因此误差范围不一致。图1是l=n=3,m=2时非整阶勒让德函数和整阶勒让德函数扰动重力梯度张量误差。

表1 整阶勒让德函数和非整阶勒让德函数的误差均方根Tab.1 RMS errors between Legendre function and Legendre function with non-integral degree

图1 扰动重力梯度误差(l=n=3、m=2)Fig.1 Error of disturbing gravity gradients(l=n=3、m=2)

2 不同坐标系下重力梯度张量之间的转换

如图2所示，GOCE卫星任务Level 2产品的EGG_TRF_2中重力梯度数据的轨道位置使用地心地固坐标系(θ,λ,r)表示，其中θ表示地心余纬，λ表示地心经度，r是卫星轨道距地心距离；重力梯度张量使用的是局部指北坐标系(LNOF)，因此使用GOCE卫星重力梯度恢复球冠谐重力场模型时，需要进行坐标变换。首先是将轨道位置坐标(θ,λ)转换到球冠坐标系[20]。定义计算点P在球坐标下的位置用(θ,λ)表示，球冠坐标系下的余纬和经度用δ和λc表示，设球冠中心点Q在球坐标系中为(θ0,λ0)，则根据球面三角公式可知，计算点P在球冠坐标系中的余纬δ和λc经度计算公式为

(16)

然后将重力梯度张量从局部指北坐标系(LNOF)转换到球冠指北坐标系(SCLNOF)中。定义LNOF坐标系到SCLNOF坐标系的转换矩阵。设β为球冠上点P到球冠中心点Q的方位角,如图2。显然SCLNOF坐标系就是将点P处的LNOF坐标系绕Z轴旋转一个β角度得到，所以LNOF坐标系到SCLNOF坐标系的梯度张量转换矩阵为

(17)

式中，方位角β的cosβ与sinβ的计算公式为

(18)

式中，δ为点P在球冠坐标下的余纬。

图2 球冠坐标系Fig.2 Spherical cap coordinate system

3 球冠谐梯度公式恢复局部重力场的数值模拟计算和分析

(19)

(20)

表2 扰动重力梯度Txz分量统计结果Tab.2 Disturbing gravity gradient Txz survey results E

表3 扰动重力统计结果Tab.3 Disturbing gravity survey results 10-5 m/s2

图3 球冠谐模型对应的L曲线Fig.3 L-curve according to spherical cap harmonic models

从表6可以看出，4个球冠谐模型的Laplace方程的精度都在±10-8E，由于GOCE卫星任务的重力梯度张量对角线分量的误差PSD在测量带宽范围内不大于±3.2×10-3E/Hz1/2(见文献[16])，因此文中计算的球冠谐模型完全满足卫星重力梯度测量的精度要求，可以用于恢复地球局部重力场模型。

表4 球冠谐模型扰动重力梯度的误差Tab.4 Disturbing gravity gradient error of the spherical cap harmonic models E

表5 球冠谐模型径向扰动重力的误差Tab.5 Disturbing gravity error of the spherical cap harmonic models 10-5 m/s2

表6 球冠谐模型计算的扰动重力梯度张量对角线分量所满足Laplace方程的精度Tab.6 The diagonal components of the disturbing gravity gradients computed by the spherical cap harmonic models that fulfill the precision of Laplace equation E

图4 球冠谐模型径向扰动重力误差Fig.4 Error of the spherical cap harmonic models

4 结 论

本文从球冠谐分析理论出发，详细推导了球冠坐标系下的扰动位重力梯度计算公式，将GOCE卫星任务重力梯度观测值延拓到平均轨道面上。网格化后，使用Tikhonov正则化算法反演恢复了4个球冠谐重力场模型，最后与使用GO_CONS_GCF_2_DIR_R5重力场模型计算的卫星平均轨道面扰动重力进行比较分析。计算结果表明，由于扰动重力梯度中各个分量对重力场频谱的敏感性是不一样的，因此单一径向扰动重力梯度分量不能恢复全部扰动重力场的信号，所以联合不同的重力梯度分量，利用其频谱的互补性恢复的局部重力场模型，能够包含更多的重力场信号。如果从扰动重力梯度数据中移除重力场的长波信号，这时无论是使用单一径向梯度分量，还是使用多个梯度分量联合的方式恢复的球冠谐模型，都能够很好地拟合扰动重力场中剩余的重力信号，与GO_CONS_GCF_2_DIR_R5模型的差异不大于±0.3x10-5m/s2。

[1] 吴晓平.利用卫星重力梯度与地面数据确定地球重力场的分析[J].测绘学院学报，1991(1)：19-32.

WU Xiaoping.Some Analyses of the Determination of the Earth’s Gravity Field from the Satellite Gravity Gradient and Terrestrial Data[J].Journal of the Institute of Surveying and Mapping，1991(1)：19-32.

[2] 孟嘉春，蔡喜楣.卫星重力梯度测量及其应用前景探讨[J].地球物理学报，1991，34(3)：369-376.

MENG Jiachun，CAI Ximei.Approach on Satellite Gravity Gradiometry and Its Vistas of Applications[J].Acta Geophysica Sinica，1991，34(3)：369-376.

[3] 吴星，张传定，赵东明.卫星重力梯度分量的广义轮胎调和分析[J].测绘学报，2009，38(2)：101-107.

WU Xing，ZHANG Chuanding，ZHAO Dongming.Generalized Torus Harmonic Analysis of Satellite Gravity Gradients Component[J].Acta Geodaetica et Cartographica Sinica，2009，38(2)：101-107.

[4] REED G B.Application of Kinematical Geodesy for Determining the Short Wave Length Components of the Gravity Field by Satellite Gradiometry[R].Columbus：Ohio State University，1973.

[5] 张传定，吴晓平，陆仲连.全张量重力梯度数据的谱表示方法[J].测绘学报，2000，29(4)：297-304.

ZHANG Chuanding，WU Xiaoping，LU Zhonglian.Spectral Representation of the Full Gravity Tensor[J].Acta Geodaetica et Cartographica Sinica，2000，29(4)：297-304.

[6] 王燚，姜效典，李德勇.多源重力数据球冠谐模型抗差融合法[J].测绘学报，2015，44(9)：952-957.DOI：10.11947/j.AGCS.2015.20140345.

WANG Yi，JIANG Xiaodian，LI Deyong.The Robust Fusion of Multi-source Gravity Data Based on the Spherical Cap Harmonic Model[J].Acta Geodaetica et Cartographica Sinica，2015，44(9)：952-957.DOI：10.11947/j.AGCS.2015.20140345.

[7] HAINES G V.Spherical Cap Harmonic Analysis[J].Journal of Geophysical Research：Solid Earth，1985，90(B3)：2583-2591.

[8] 申宁华，王光杰，王喜臣.球冠谐和分析算法及其在中国大陆航磁异常分析中的应用[J].物探化探计算技术，1998，20(1)：9-18.

SHEN Ninghua，WANG Guangjie，WANG Xichen.Spherical Cap Harmonic Analysis and Its Application to the Analysis of the Aeromagnetic Anomalies in China[J].Computing Techniques for Geophysical and Geochemical Exploration，1998，20(1)：9-18.

[9] LI Jiancheng，CHAO Dingbo,NING Jinsheng.Spherical Cap Harmonic Expansion for Local Gravity Field Representation[J].Manuscripta Geodaetica，1995，20(2)：265-277.

[10] 王燚，姜效典.地球重力场球冠谐模型的分层构造和分析[J].大地测量与地球动力学，2014，34(5)：30-34，39.

WANG Yi，JIANG Xiaodian.Structure and Analysis of Multilayered Spherical Cap Harmonic Model[J].Journal of Geodesy and Geodynamics，2014，34(5)：30-34，39.

[11] 吴招才，刘天佑，高金耀.局部重力场球冠谐分析中的导数计算及应用[J].海洋与湖沼，2006，37(6)：488-492.

WU Zhaocai，LIU Tianyou，GAO Jinyao.Derivative Calculation and Application in Spherical Cap Harmonic Analysis of Local Gravity Field[J].Oceanologia et Limnologia Sinica，2006，37(6)：488-492.

[12] BRUINSMA S L，FÖRSTE C，ABRIKOSOV O，et al.The New ESA Satellite-only Gravity Field Model via the Direct Approach[J].Geophysical Research Letters，2013，40(14)：3607-3612.

[13] 曹月玲，王解先.利用球冠谐分析拟合GPS水准高程[J].武汉大学学报(信息科学版)，2008，33(7)：740-743.

CAO Yueling，WANG Jiexian.Application of Spherical Cap Harmonic Analysis to Fit GPS Level Height[J].Geomatics and Information Science of Wuhan University，2008，33(7)：740-743.

[14] 彭富清，于锦海.球冠谐分析中非整阶Legendre函数的性质及其计算[J].测绘学报，2000，29(3)：204-208.

PENG Fuqing，YU Jinhai.The Characters and Computation of Legendre Function with Non-integral Degree in SCHA[J].Acta Geodaetica et Cartographica Sinica，2000，29(3)：204-208.

[15] HWANG C，CHEN S K.Fully Normalized Spherical Cap Harmonics：Application to the Analysis of Sea-level Data from TOPEX/Poseidon and ERS-1[J].Geophysical Journal International，1997，129(2)：450-460.

[16] 刘晓刚.GOCE卫星测量恢复地球重力场模型的理论与方法[D].郑州：信息工程大学，2011.

LIU Xiaogang.Theory and Methods of the Earth’s Gravity Field Model Recovery from GOCE Data[D].Zhengzhou：Information Engineering University，2011.

[17] HOBSON E W.The Theory of Spherical and Ellipsoidal Harmonics[M].Cambridge：Cambridge University Press，1931.

[18] 王竹溪，郭敦仁.特殊函数概论[M].北京：北京大学出版社，1989.

WANG Zhuxi，GUO Dunren.Introduction to Special Function[M].Beijing：Peking University Press，1989.

[19] 刘晓刚，吴晓平，赵东明，等.扰动重力梯度的非奇异表示[J].测绘学报，2010，39(5)：450-457.

LIU Xiaogang，WU Xiaoping，ZHAO Dongming，et al.Non-singular Expression of the Disturbing Gravity Gradients[J].Acta Geodaetica et Cartographica Sinica，2010，39(5)：450-457.

[20] 赵建虎，王胜平，刘辉，等.海洋局域地磁场球冠谐分析建模方法研究[J].测绘科学，2010，35(1)：50-52，9.

ZHAO Jianhu，WANG Shengping，LIU Hui，et al.Study on Establishing Local Geomagnetic Model Using Spherical Cap Harmonic Analysis[J].Science of Surveying and Mapping，2010，35(1)：50-52，9.

[21] 罗志才，钟波，宁津生，等.卫星重力梯度测量确定地球重力场的理论与方法[M].武汉：武汉大学出版社，2015.

LUO Zhicai，ZHONG Bo，NING Jinsheng，et al.Theory and Method for Determining the Earth’s Gravity Field from Satellite Gravity Gradiometry[M].Wuhan：Wuhan University Press，2015.

[22] 韦博成，鲁国斌，史建清.统计诊断引论[M].南京：东南大学出版社，1991.

WEI Bocheng，LU Guobin，SHI Jianqing.Introduction to Statistical Diagnostics[M].Nanjing：Publishing House of Dongnan University，1991.

[23] 黄谟涛，欧阳永忠，刘敏，等.融合海域多源重力数据的正则化点质量方法[J].武汉大学学报(信息科学版)，2015，40(2)：170-175.

HUANG Motao，OUYANG Yongzhong，LIU Min，et al.Regularization of Point-mass Model for Multi-source Gravity Data Fusion Processing[J].Geomatics and Information Science of Wuhan University，2015，40(2)：170-175.

[24] HANSEN P C.Analysis of Discrete Ill-posed Problems by Means of the L-curve[J].SIAM Review，1992，34(4)：561-580.

The Spherical Cap Harmonic Analysis Modeling Method Based on Disturbing Gravity Gradients

WANG Yi，JIANG Xiaodian

College of Marine Geo-science, Ocean University of China, Qingdao 266100，China

It is deduced that the non-singular computational formulae of the disturbing gravity gradients based on spherical cap harmonic analysis theory.On the basis of the Tikhonov regularization method，the spherical cap harmonic model of local gravity field is calculated using GOCE satellite real surveying data.It is shown that the short wavelength information of local gravity field can be recovered from the disturbing gravity gradient data using spherical cap harmonic analysis method，and the difference with the GO_CONS_GCF_2_DIR_R5 model is at ±0.3×10-5m/s2.

spherical cap harmonic analysis; disturbing gravity; gravity gradients;Tikhonov regulariza-tion method

The Oceanic Public Nonprofit Research Project(No.201305029-02)

WANG Yi(1974—)，male,PhD，lecturer，majors in geophysical exploration.

王燚，姜效典.扰动重力梯度的球冠谐分析建模[J].测绘学报，2017,46(11)：1802-1811.

10.11947/j.AGCS.2017.20160412.

WANG Yi，JIANG Xiaodian.The Spherical Cap Harmonic Analysis Modeling Method Based on Disturbing Gravity Gradients[J].Acta Geodaetica et Cartographica Sinica,2017,46(11):1802-1811.DOI:10.11947/j.AGCS.2017.20160412.

P223

A

1001-1595(2017)11-1802-10

海洋公益性行业科研专项经费(201305029-02)

(责任编辑：丛树平)

2016-08-18

修回日期：2017-07-26

王燚(1974—)，男，博士，讲师，研究方向为地球物理勘探方法与信息技术。

E-mail：wyfox009@163.com