地球物理学报  2013, Vol. 56 Issue (7): 2245-2256   PDF    
引用本文
刘金钊, 柳林涛, 梁星辉 等. 基于实测重力异常和地形数据确定重力梯度的研究. 地球物理学报, 2013, 56(7): 2245-2256, doi: 10.6038/cjg20130712.  
Liu J Z, Liu L T, Liang X H, et al. Based on measured gravity anomaly and terrain data to determine the gravity gradients. Chinese J. Geophys. (in Chinese), 2013, 56(7): 2245-2256, doi: 10.6038/cjg20130712.  
基于实测重力异常和地形数据确定重力梯度的研究
刘金钊1,2 , 柳林涛1 , 梁星辉1 , 叶周润1,2 , 苏晓庆1,2     
1. 中国科学院测量与地球物理研究所大地测量与地球动力学国家重点实验室, 武汉 430077;
2. 中国科学院大学, 北京 100049
收稿日期 2012-07-04, 2013-06-18 收修定稿
基金项目: 科技部国家重大科学仪器设备开发专项(2011YQ120045)和国家自然科学基金(41074050)资助
作者简介: 刘金钊, 男, 1986年生, 博士研究生, 从事航空重力梯度方向的研究.E-mail:whiggliujinzhao@126.com
通讯作者: 柳林涛, 男, 1967年生, 研究员, 博士生导师, 主要从事时频分析, 航空/海洋重力测量.E-mail:llt@asch.whigg.ac.cn
摘要: 本文利用澳大利亚北领地West Arnhem Land地区实测重力异常数据并联合DEM (9")和SRTM3(3")地形高程数据, 使用移去-恢复技术和Stokes积分方法计算了该地区两条剖面的重力梯度及其功率谱密度, 使用FFT方法解算了整个地区的重力梯度值, 结果证明了联合重力异常数据和高分辨率地形高程数据能有效地提高重力梯度的解算精度; 功率谱密度的计算结果与国外成熟的重力梯度功率谱密度模型相吻合, 表明高于0.3 Hz频率范围的功率谱密度可看做噪声, 为重力梯度数据处理中噪声的辨别和剔除提供了借鉴, 另外对重力梯度辅助导航基准图的构建以及重力梯度测量系统的标定提供了有益的探索.
关键词: 重力梯度      移去-恢复法      SRTM3      Stokes积分      功率谱密度     
Based on measured gravity anomaly and terrain data to determine the gravity gradients
LIU Jin-Zhao1,2, LIU Lin-Tao1, LIANG Xing-Hui1, YE Zhou-Run1,2, SU Xiao-Qing1,2     
1. State Key Laboratory of Geodesy and Earth's Dynamics, Institute of Geodesy and Geophysics, CAS, Wuhan 430077, China;
2. University of Chinese Academy of Sciences, Beijing 100049, China
Abstract: Based on gravity anomaly data and terrain elevation data (DEM (9") and SRTM3(3")), this paper calculates the gravity gradients and its power spectral density of two profiles in West Arnhem Land of Australia's northern territory. We implement Remove-Recovery and Stokes integral method to obtain the outcomes. We also implement FFT methods to calculate the gravity gradients covering the entire area. The results show that combination of gravity anomaly data with high spatial resolution terrain elevation data can effectively improve the accuracy of gravity gradients. Calculated power spectral density matches well with the power spectral density model, it shows that power spectral density greater than 0.3 Hz can be regarded as noise. The results provide a reference for data processing in identifying and eliminating the noise in gradients. In addition, it provides beneficial exploration for reference map of gravity gradient assisted navigation and for gravity gradient measurement system calibration..
Key words: Gravity gradients      Remove-Recovery      SRTM3      Stokes integral      PSD     
1 引言

地球表面的重力梯度值非常微小，大约只有几十到几千个E(1E=10-9 s-2)，相对于重力异常，重力梯度场能够更好地反映出地球表面浅层地质体的密度变化和结构特征，对矿藏的探测以及重力梯度匹配辅助导航都有着十分重要的意义；另外，随着卫星重力梯度和航空重力梯度测量系统的试验和发展，对星载或航空重力梯度测量系统的有效性进行检核，就需要选取一块实际重力梯度值已知的标定场，通过对上述测量系统在该标定场内进行试验和验证，剔除系统中不合理的误差因素，获得仅由地质结构变化引起的重力梯度值.相对来说，目前国内航空重力梯度测量难以大范围展开，在现有水平条件下，利用重力异常和高分辨率地形高程数据解算重力梯度场就成为一种行之有效的方法，具有十分重要的现实意义[1-3].

航空重力梯度的解算方法很多：由于重力梯度对地球表层地形变化十分敏感，所以很多学者提出了单独利用地形数据来解算重力梯度，如文献[4-7]中利用地形数据解算重力梯度的研究，但是由于没有相应密度异常数据，单独利用地形数据解算重力梯度的方法通常假定地层密度值为某一常数(一般代以地层平均密度值2.67 g/cm3)，这种方法难以符合地质结构的实际情况；另外，利用重力异常解算重力梯度也是一种常见的方法，许多学者对此也做过大量的研究[8]；在现有条件下，高分辨率的重力异常数据难以获取，一定程度上限制了重力梯度的解算精度，而重力异常在一定程度上恰恰反映了地层结构的密度差异，所以综合利用重力异常数据和高分辨率地形数据解算重力梯度能够更好地反映地质结构重力梯度的实际变化，这方面文献[9-11]作出了相应的研究.

本文利用澳大利亚北领地某次重力测量公布的重力异常数据及高分辨率地形高程数据[12-13]，利用Stokes积分法和FFT法计算了相应的重力梯度及其功率谱密度，与国外成熟的重力梯度功率谱密度模型进行了比较，结果证明了这种方法解算重力梯度的合理性和有效性，对重力梯度标定场的确定和重力梯度辅助导航基准图的构建都有一定的现实意义[14-17].

2 重力梯度解算的原理和方法

利用EGM2008模型和地形数据分别解算重力梯度时，我们发现相对高分辨率地形数据，EGM2008模型分辨率过低，尺度过大，解算得到的重力梯度与地形数据解算得到的重力梯度相比，数值较小，所以在后面的计算中我们将由EGM2008解算得到的重力梯度看成偏差.利用EGM2008计算重力梯度的公式为[18-20]：

(1)

其中Γjk为重力梯度分量，p为计算点，GM为地球引力常数和地球质量的乘积，R为平均地球半径；

(2)

Cnm，Snm为球谐系数，Pnm为完全规格化缔合Legendre函数，rp、θp、λp为计算点的地心距、余纬和经度.

球谐坐标(r，λ，θ)与局部坐标(x1，x2，x3)二阶导数的关系为

(3)

利用地形数据解算重力梯度公式如下：

(4)

式中：

(5)

为计算点和积分点的空间距离.

①模型一

物理大地测量中，Stokes积分作为经典重力位边值问题的解，将大地水准面上的重力异常和扰动位联系起来了，利用扰动位的二阶微分可以得到重力梯度，一般来说，自由空间重力异常表征了大地水准面以上地形和地下介质的重力场影响，能够完整地反映实际地区的重力变化，是大地测量中应用最为广泛的一种重力异常数据[14-21].

重力梯度总体来说可分为三个部分，由EGM2008模型计算的重力梯度可作为总重力梯度场的长波分量，局部地区重力梯度主要来源于地形变化，高精度的地形高程数据解算的重力梯度可以看做总重力梯度场的短波分量，而中长波部分的重力梯度可由重力异常求得，同时重力异常数据也能够弥补高分辨率地形高程数据无法反映地质体密度变化的缺陷.

模型一通过利用移去-恢复技术，来联合地形数据和重力异常数据解算重力梯度：首先将自由空间重力异常中地形影响部分的重力异常扣除(这部分重力异常可根据高分辨率地形高程数据求得)，由剩余重力异常求得重力梯度后，再恢复地形影响产生的重力梯度，得到最后总的重力梯度值[9].模型见公式(6)，离散表达式见附录.

(6)

式中R是地球平均半径，Δgfree-air和Δgterrain分别表示自由空间重力异常和地形影响重力异常，S(r，ψ)为Stokes函数，σ为积分区域，Γjk|terrain为地形影响部分的重力梯度.

(7)

(x1′，x2′，x3′)为积分点，(xp，yp，hp)为计算点，而

(8)

(9)

式中

(10)

②模型二

我们同时获取了自由空间重力异常和完全布格重力异常数据，完全布格重力异常是在自由空间重力异常的基础上进行了剩余地形改正、布格片改正和曲率改正、消除了大地水准面以上地形质量的影响得到的[21]：

(11)

式中δgFA为自由空间重力异常，δgBA为完全布格重力异常，δgBC为布格片改正，δgBB为曲率改正，δgTC为剩余地形改正.

我们联合地形高程数据和完全布格重力异常重新计算了自由空间重力异常，将公式中地形改正部分用地形数据进行计算，并根据快速傅里叶算法，利用计算得到的自由空间重力异常解算了整个区域的重力梯度值，并和实际自由空间重力异常解算的重力梯度进行了比较.

模型见公式(12)：

(12)

ϑ[DEM]为数字高程模型数据重新计算剩余地形改正.根据文献[22]中介绍的FFT算法建立各重力梯度分量的傅里叶变换与重力异常傅里叶变换的关系：

(13)

式中，分别为傅里叶变换及其逆变换，

(14)

(15)

kx，ky，kz分别是x，y，z方向的波数.

3 实测数据实验

2003年8月27日-9月15日，Fugro Airborne Surveys Pty Ltd和Canadian Micro Gravity Pty Ltd两家公司利用俄罗斯GT-1A航空重力仪，在澳大利亚北领地WestArnhem Land(见图 1)进行了为期三周的航空重力测量，飞机飞行高度为655 m，总测线长度约为4725.5 km，随后，Geoscience Australia网站公布了滤波后的重力异常数据[12-13].

图 1 航空重力测量测区范围 Fig. 1 Surveying area of airborne gravimetry

下载得到的数据包括时长为60 s，80 s和107 s(s为时间秒)滤波处理后的自由空间重力异常和完全布格重力异常，同时，我们还下载了该区域的DEM(9″)数据，以及相同测区的SRTM3(3″)数据(″为空间秒).

我们选取的计算范围为实际测区中的一块格网区域，重力异常数据和地形数据见图 2及表 1.

图 2 重力异常数据和地形数据 (a)自由空间重力异常；(b)完全布格重力异常；(c)DEM(9″)数据；(d)SRTM3(3″)数据. Fig. 2 Gravity anomaly data and terrain elevation data (a) Free air gravity anomaly; (b) Complete Bouguer gravity anomaly; (c) DEM (9″) data; (d) SRTM3(3″) data.
表 1 数据信息统计 Table 1 Data information statistics

SRTM3数据与前面的重力异常数据参考基准不同，需要进行一定的延拓，我们这里为简单起见，直接进行了计算.

我们在测区模拟了两条剖面，剖面1的数据分辨率为0.001°，长度约23 km，剖面2的数据分辨率为0.001°，长度约42 km，见图 3，并计算了相应两条剖面的重力梯度值.根据DEM(9″)数据和SRTM3(3″)数据，我们首先利用棱柱法和直接积分法计算了该地区地形质量在剖面1上的重力梯度，以及EGM2008模型计算得到的剖面重力梯度，结果见图 4.

图 3 两剖面在测区中的位置 Fig. 3 Location of two profiles in surveying area
图 4 地形数据和EGM2008模型计算的剖面1，剖面2重力梯度值 (a)DEM(9″)计算结果；(b)SRTM(3″)计算结果. Fig. 4 Gravity gradients of profile 1 and 2 derived form terrain data and EGM2008 model (a) Results derived from DEM (9″) data; (b) Results derived from SRTM (3″) data.

从图 4中可以看出，因为分辨率的关系，由EGM2008模型计算的重力梯度值相对地形数据计算的重力梯度非常小，局部地区重力梯度主要来源于地形质量影响.

利用棱柱法和直接积分法对两种分辨率地形数据计算结果的差值列于表 2.从表中可以看出，利用空间分辨率为3″(~90 m)的地形数据计算得到的差值更小，结果更精确，为了作图简洁，我们后面仅给出利用SRTM3(3″)地形数据和重力异常数据融合的计算结果.

表 2 棱柱法与直接积分法计算重力梯度之差(单位：E) Table 2 Differences of gravity gradients calculated by prism and direct integral method

利用模型一，我们计算了剖面1上各部分的重力梯度值，见图 5，其中，灰色虚线表示地形影响得到的重力梯度；黑色实线表示剩余重力异常得到的重力梯度；黑色虚线表示EGM2008计算的重力梯度.综合利用重力异常和地形数据计算得到的总重力梯度见图 6.

图 5 各部分重力梯度值 Fig. 5 Gravity gradient of eachpart
图 6 两条剖面总重力梯度值 (a)剖面1总重力梯度值；(b)剖面2总重力梯度值；黑色虚线为自由空间重力异常计算的梯度值，灰色实线为综合地形数据和重力异常计算的梯度值. Fig. 6 Total gravity gradients of two profiles (a) Total gravity gradients of profile 1;(b) Total gravity gradients of profile 2;Black dashed line represents gradients only derived from free air gravity anomaly and gary line represents gradients derived from gravity anomaly and terrain data comprehensively.

图 6中，黑色虚线为直接利用自由空间重力异常解算得到的重力梯度，灰色实线为综合利用重力异常和高分辨率地形数据计算的重力梯度，由于重力异常数据分辨率比地形数据要低，解算得出的重力梯度曲线变化也较为平缓，而高分辨地形数据的融入，增强了剖面重力梯度的细节纹理.另外利用模型二，根据快速傅里叶算法，我们解算了整个区域的重力梯度，结果见图 7，其中(a，c，e，g，i，k)为直接利用已知自由空间重力异常数据解算得到的结果，(b，d，f，h，j，l)为综合利用完全布格重力异常数据和地形数据解算得到的结果.

图 7 整个区域的重力梯度分量 (a，c，e，g，i，k)直接利用自由空间重力异常解算的重力梯度值；(b，d，f，h，j，l)综合完全布格重力异常和地形数据解算的重力梯度值. Fig. 7 Gravity gradient components of entire area (a, c, e, g, i, k) Gravity gradients derived from free air gravity anomaly by FFT method; (b, d, f, h, j, l) Gravity gradients derived from terrain data and free air gravity anomaly comprehensively by FFT method.

从图 7可以看出高分辨率地形数据和重力异常数据的融合，增强了梯度数据解算的纹理细节部分.对比地形高程图，以重力梯度分量Γzz(图 7k和7l)为例，红色圆圈内为相应地形图上比较明显的地形特征，利用自由空间重力异常解算的梯度分量Γzz没有很好地反映该地形特征，而综合解算的Γzz分量却明显地反映了该地形结构变化，此外分量Γxx，Γxz，Γyy也从不同方面确定了测区内的地形结构变化，说明了重力梯度比重力异常数据在地质结构的表现上更加丰富，这在利用梯度数据探测地下矿藏应用方面表现的更为明显.

4 功率谱密度解算

针对图 6解算得到的剖面重力梯度，我们计算了相应的功率谱密度(Power Spectral Density，PSD).设飞机飞行的高度为655 m，速度为250 km/h，空间频率和时间频率的转换关系根据公式[23-25]

(16)

功率谱密度解算根据周期图法[26-29]，公式见(17)，(18).

(17)

(18)

式中f表示时间频率，单位是Hz；v表示速度，φ表示时间域功率谱密度，S表示空间域功率谱密度；X(ωp)是f(xk)的一维离散傅里叶变换，K是总采样点数，xk=k·Δx，Δx是数据间隔，D是数据距离长度，ωp=p/D是频率，p=0，1，…，K-1；φ(ωp)是一维功率谱密度.

功率谱密度解算结果见图 8(a、b)，其中黑色虚线是根据周期图法解算得到的剖面重力梯度功率谱密度，灰色实线为根据模型计算得到的功率谱密度[30-32].

图 8 剖面数值解算PSD与模型PSD比较 (a)剖面1综合解算PSD与模型PSD比较；(b)剖面2综合解算PSD与模型PSD比较. Fig. 8 Compare between synthetically calculated PSD and model PSD (a) Synthetically calculated PSD and model PSD of profile 1;(b) Synthetically calculated PSD and model PSD of profile 2.

同时我们也计算了剖面2直接由自由空间重力异常解算的重力梯度的功率谱密度值，与模型功率谱的比较结果见图 9.

图 9 剖面2直接解算PSD与模型PSD比较 Fig. 9 Directly calculated PSD and model PSD of profile 2

我们将联合解算的重力梯度功率谱密度和由自由空间重力异常直接解算的重力梯度功率谱密度分别与功率谱密度模型求差，相应统计量见表 3.

表 3 重力梯度分量功率谱密度差值统计量(单位:) Table 3 Statistics of PSD differences of gravity gradient components

综合图 8-9，我们发现利用高分辨地形数据和重力异常数据联合解算的重力梯度功率谱密度与已有的重力梯度功率谱模型相吻合，从图中我们可以看出重力梯度的功率谱密度呈衰减趋势，大部分能量集中在0.3 Hz频率范围以内，这为我们在梯度数据处理中，剔除不合理的高频噪声提供了依据.因为地球表面重力梯度场比重力场衰减的更加迅速，所以航空重力梯度测量中对飞行载体的高度也有一定选择，另外，飞行载体的飞行速度一定程度上也会对数据分辨率产生影响.由表 3，可以看到除了分量PΓ22的标准差(4.8946＞3.9194)，其他重力梯度分量的功率谱密度之差的标准差在融入高精度地形数据后都得到提高，所以我们认为，本文介绍的高精度地形数据和重力异常综合解算重力梯度，提高了重力梯度的解算精度.

5 结论

本文利用重力异常数据和高分辨SRTM3地形高程数据联合解算了两条剖面及整块测区的重力梯度，结果显示了高分辨率地形数据的融入有效地增强了重力梯度的细节纹理.重力梯度是扰动位的二次微分，比重力异常对地质结构的变化更加敏感，同时，重力梯度在空间上有9个分量(5个独立分量)，比重力异常对地质结构的反映上也更为丰富，但是随着载体飞行高度的增加，重力梯度衰减的也更为迅速，所以重力梯度主要反映局部地区和浅层地表的地质结构特征.由联合解算的剖面重力梯度功率谱密度和C.Jekeli的功率谱密度模型比较发现，二者吻合较好，间接说明了本文介绍的重力梯度解算方法的可行性和有效性，对重力梯度实验场的确立和重力梯度辅助匹配导航提供了有益的探索；此外利用功率谱分析方法可以看出，在重力梯度功率谱密度中，重力梯度能量集中在0.3 Hz以内，频率大于该阈值的各种因素可视为噪声，进行剔除，这对提高重力梯度测量的信噪比，以及后期重力梯度的数据处理提供了借鉴.最后，根据重力梯度功率谱密度差值统计量，我们发现重力异常数据在融入高分辨率地形高程数据后，6个重力梯度分量中有5个分量的标准差得到了提高，所以我们认为本文的方法对提高重力梯度的解算精度具有一定的参考价值.

附录

其中

参考文献
[1] Lee J B. FALCON gravity gradiometer technology. Exploration Geophysics , 2001, 32(4): 247-250. DOI:10.1071/EG01247
[2] While J, Jackson A, Smit D, et al. Spectral analysis of gravity gradiometry profiles. Geophysics , 2006, 71(1): J11-J22. DOI:10.1190/1.2169848
[3] Moody M V. A superconducting gravity gradiometer for measurements from a moving vehicle. Review of Scientific Instruments , 2011, 82(9): 094501. DOI:10.1063/1.3632114
[4] 武凛, 马杰, 周瑶, 等. 重力场匹配导航的全张量重力梯度基准图模拟. 系统仿真学报 , 2009, 21(22): 7037–7041. Wu L, Ma J, Zhou Y, et al. Modelling full-tensor gravity gradient maps for gravity matching navigation. Journal of System Simulation (in Chinese) , 2009, 21(22): 7037-7041.
[5] 边少锋, 张赤军. 地形起伏对重力垂直梯度影响的计算. 物探化探计算技术 , 1999, 21(2): 133–140. Bian S F, Zhanf C J. Computation of topographic effects on vertical gravity gradient. Computing Techniques for Geophysical and Geochemical Exploration (in Chinese) , 1999, 21(2): 133-140.
[6] 张赤军. 用地形数据确定重力异常垂直梯度. 科学通报 , 1999, 44(6): 656–661. Zhang C J. Using topographic data to determine the vertical gradient of gravity anomaly. Chinese Science Bulletin (in Chinese) , 1999, 44(6): 656-661.
[7] Jekeli C, Zhu L Z. Comparison of methods to model the gravitational gradients from topographic data bases. Geophysical Journal International , 2006, 166(3): 999-1014. DOI:10.1111/gji.2006.166.issue-3
[8] 武凛, 胡维, 马杰, 等. 基于重力异常分析的重力梯度图制备方法. 华中科技大学学报(自然科学版) , 2009, 37(11): 57–60. Wu L, Hu W, Ma J, et al. Modeling full-tensor gravity gradient maps using gravity anomaly analysis. Journal of Huazhong University of Science & Technology (Natural Science Edition) (in Chinese) , 2009, 37(11): 57-60.
[9] Zhu L Z, Jekeli C. Gravity gradient modeling using gravity and DEM. Journal of Geodesy , 2009, 83(6): 557-567. DOI:10.1007/s00190-008-0273-2
[10] Zhu L, Jekeli C. Combining gravity and topographic data for local gradient modeling. Dynamic Planet , 2007, 130: 288-295. DOI:10.1007/978-3-540-49350-1
[11] Zhu L Z. Gradient Modeling with Gravity and DEM. Ohio: Ohio State University, 2007. www.geology.osu.edu/~jekeli.1/OSUReports/reports/report_483.pdf.
[12] Australia G. Geophysical Archive Data Delivery System (GADDS)[EB/OL].[2010-05-10]. http://www.ga.gov.au/gadds.
[13] Gabell A R, Tuckett H. Project 200380 West Arnhem Land GT-1A Airborne Gravity Survey. Canberra: Geoscience Australia and Northern Territory Geological Survey, 2004 .
[14] Moritz H. Kinematical Geodesy II, 1971. Report No. 92. Reports of the Department of Geodetic Science. www.dtic.mil/dtic/tr/fulltext/u2/666052.pdf.
[15] Forsberg R. Gravity field terrain effect computations by FFT. Bull. Geod. , 1985, 59(4): 342-260. DOI:10.1007/BF02521068
[16] Dransfield M H. Airborne Gravity Gradiometry[Ph. D. thesis]. Australia: University of Western Australia, 1994.
[17] Bell R. Gravity gradiometry. Scientific American , 1998, 278: 74-79.
[18] Petrovskaya M S, Vershkov A N. Development of the Second-order derivatives of the Earth's potential in the local north-oriented reference frame in orthogonal series of modified spherical harmonics. Journal of Geodesy , 2008, 82(12): 929-944. DOI:10.1007/s00190-008-0223-z
[19] Pertrovikaya M S, Vershkov A N. Non-singular expressions for the gravity gradients in the local north-oriented and orbital reference frames. Journal of Geodesy , 2006, 80(3): 117-127. DOI:10.1007/s00190-006-0031-2
[20] 郭春喜, 宁津生, 陈俊勇, 等. 珠峰地区似大地水准面精化与珠峰顶正高的确定. 地球物理学报 , 2008, 51(1): 101–107. Guo C X, Ning J S, Chen J Y, et al. Improvement of regional quasi-geoid in Qomolangma (Mt. Everest) and determination of orthometric elevation. Chinese J. Geophys (in Chinese) , 2008, 51(1): 101-107.
[21] 钱东, 刘繁明, 李艳, 等. 导航用重力梯度基准图构建方法的比较研究. 测绘学报 , 2011, 40(6): 736–744. Qian D, Liu F M, Li Y, et al. Comparison of gravity gradient reference map composition for navigation. Acta Geodaetica et Cartographica Sinica (in Chinese) , 2011, 40(6): 736-744.
[22] Mickus K L, Hinojosa J H. The complete gravity gradient tensor derived from the vertical component of gravity: a Fourier transform technique. Journal of Applied Geophysics , 2001, 46(3): 159-174. DOI:10.1016/S0926-9851(01)00031-3
[23] Jekeli C. Correlation modeling of the gravity field in classical geodesy.//Freedn W, Nashed M Z, Sonar T eds. Handbook of Geomathematics. Heidelberg: Springer Verlag, 2010: 833-864.
[24] Jekeli C. Airborne gradiometry error analysis. Surveys in Geophysics , 2006, 27(2): 257-275. DOI:10.1007/s10712-005-3826-4
[25] Jekeli C. Statistical analysis of moving-base gravimetry and gravity gradiometry. Report no. 466, Geodetic Science. Columbus: Ohio State University.http://www.geology.osu.edu/~jekeli.1/OSUReports/reports/report_466.pdf.
[26] Heck B. On Helmert's methods of condensation. Journal of Geodesy , 2003, 77(3-4): 155-170. DOI:10.1007/s00190-003-0318-5
[27] Flury J. Short-wavelength spectral properties of the gravity field from a range of regional data sets. Journal of Geodesy , 2006, 79(10-11): 624-640. DOI:10.1007/s00190-005-0011-y
[28] 柳林涛, 许厚泽. 航空重力测量数据的小波滤波处理. 地球物理学报 , 2004, 47(3): 490–494. Liu L T, Xu H Z. Wavelets in airborne gravimetry. Chinese J. Geophys (in Chinese) , 2004, 47(3): 490-494. DOI:10.1002/cjg2.v47.3
[29] 万晓云, 于锦海, 曾艳艳. GOCE引力梯度的频谱分析及滤波. 地球物理学报 , 2012, 55(9): 2909–2916. Wan X Y, Yu J H, Zeng Y Y. Frequency analysis and filtering processing of gravity gradients data from GOCE. Chinese J. Geophys (in Chinese) , 2012, 55(9): 2909-2916.
[30] Fukushima T. Numerical computation of spherical harmonics of arbitrary degree and order by extending exponent of floating point numbers: II first-, second-, and third-order derivatives. Journal of Geodesy , 2012, 86(11): 1019-1028. DOI:10.1007/s00190-012-0561-8
[31] Stelkens-Kobsch T H. The Airborne Gravimeter CHEK AN-A at the Institute of Flight Guidance (IFF). www.upf.pf/ICET/cd_rom_bgi/Cours/Heyen/PAPER.pdf.
[32] Vasco D W, Taylor C. Inversion of airborne gravity gradient data, southwestern Oklahoma. Geophysics , 1991, 56(1): 90-101. DOI:10.1190/1.1442961