Efficient propagation of error covariance matrices of gravitational models: application to GRACE and GOCE

We have applied efficient methods for computing variances and covariances of functions of a global gravity field model expanded in spherical harmonics, using the full variance–covariance matrix of the coefficients. Examples are given with recent models derived from GRACE (up to degree and order 150), and with simulated GOCE derived solutions (up to degree and order 200).     

Comparison of undulation difference accuracies using gravity anomalies and gravity disturbances

Errors are considered in the outer zone contribution to oceanic undulation differences as obtained from a set of potential coefficients complete to degree 180. It is assumed that the gravity data of the inner zone (a spherical cap), consisting of either gravity anomalies or gravity disturbances, has negligible error. This implies that error estimates of the total undulation difference are analyzed. If the potential coefficients are derived from a global field of 1°×1° mean anomalies accurate to ε_Δg=10 mgal, then for a cap radius of 10°, the undulation difference error (for separations between 100 km and 2000 km) ranges from 13 cm to 55 cm in the gravity anomaly case and from 6 cm to 36 cm in the gravity disturbance case. If ε_Δg is reduced to 1 mgal, these errors in both cases are less than 10 cm. In the absence of a spherical cap, both cases yield identical error estimates: about 68 cm if ε_Δg=1 mgal (for most separations) and ranging from 93 cm to 160 cm if ε_Δg=10 mgal. Introducing a perfect 30-degree reference field, the latter errors are reduced to about 110 cm for most separations.     

The Meissl scheme for the geodetic ellipsoid

We present a variant of the Meissl scheme to relate surface spherical harmonic coefficients of the disturbing potential of the Earth’s gravity field on the surface of the geodetic ellipsoid to surface spherical harmonic coefficients of its first- and second-order normal derivatives on the same or any other ellipsoid. It extends the original (spherical) Meissl scheme, which only holds for harmonic coefficients computed from geodetic data on a sphere. In our scheme, a vector of solid spherical harmonic coefficients of one quantity is transformed into spherical harmonic coefficients of another quantity by pre-multiplication with a transformation matrix. This matrix is diagonal for transformations between spheres, but block-diagonal for transformations involving the ellipsoid. The computation of the transformation matrix involves an inversion if the original coefficients are defined on the ellipsoid. This inversion can be performed accurately and efficiently (i.e., without regularisation) for transformation among different gravity field quantities on the same ellipsoid, due to diagonal dominance of the matrices. However, transformations from the ellipsoid to another surface can only be performed accurately and efficiently for coefficients up to degree and order 520 due to numerical instabilities in the inversion.     

A comparison of gravimetric undulations computed by the modified Molodensky truncation method and the method of least squares spectral combination by optimal integral kernels

Gravimetric geoid undulations have been computed by the modified Molodensky truncation method (the Meissl procedure) and by the method of least squares spectral combination by optimal kernels (the Wenzel procedure). These undulations have been compared in two manners. One comparison used Doppler-derived undulation-at 65 stations in the United States as references. A second comparison used Geos-3 derived undulations in 30°×30° areas in the Indian Ocean and Tonga Trench as references. the mean difference of undulation-computed by the Wenzel procedure was 0.6 m smaller than that of the Meissl procedure when compared to the Doppler derived undulations. The standard deviations of the differences of both procedures appeared to be not significantly different. There are no significant changes in the mean differences of both procedures when compared to Geos-3 derived undulations. The standard deviations of the differences computed by the Wenzel procedure were of the order of 0.2 m smaller than those computed by the Meissl procedure.     

基于卡尔曼滤波的地形反演方法及其仿真研究

推导基于卡尔曼滤波的海底地形反演算法,算法通过利用重力梯度仪测量重力梯度异常,提高地形反演的精度和分辨率.结合地形坡度计算理论和矩形棱柱法,建立基于卡尔曼滤波的系统方程和量测方程.对算法进行仿真,分析地形均值误差、随机误差和梯度仪测量精度对反演精度的影响,最后的区域反演仿真结果表明:反演后的地形在精度和分辨率上都有所提...     

Performance of three types of Stokes's kernel in the combined solution for the geoid

When regional gravity data are used to compute a gravimetric geoid in conjunction with a geopotential model, it is sometimes implied that the terrestrial gravity data correct any erroneous wavelengths present in the geopotential model. This assertion is investigated. The propagation of errors from the low-frequency terrestrial gravity field into the geoid is derived for the spherical Stokes integral, the spheroidal Stokes integral and the Molodensky-modified spheroidal Stokes integral. It is shown that error-free terrestrial gravity data, if used in a spherical cap of limited extent, cannot completely correct the geopotential model. Using a standard norm, it is shown that the spheroidal and Molodensky-modified integration kernels offer a preferable approach. This is because they can filter out a large amount of the low-frequency errors expected to exist in terrestrial gravity anomalies and thus rely more on the low-frequency geopotential model, which currently offers the best source of this information. Received: 11 August 1997 / Accepted: 18 August 1998     

Spherical harmonic modelling to ultra-high degree of Bouguer and isostatic anomalies

The availability of high-resolution global digital elevation data sets has raised a growing interest in the feasibility of obtaining their spherical harmonic representation at matching resolution, and from there in the modelling of induced gravity perturbations. We have therefore estimated spherical Bouguer and Airy isostatic anomalies whose spherical harmonic models are derived from the Earth’s topography harmonic expansion. These spherical anomalies differ from the classical planar ones and may be used in the context of new applications. We succeeded in meeting a number of challenges to build spherical harmonic models with no theoretical limitation on the resolution. A specific algorithm was developed to enable the computation of associated Legendre functions to any degree and order. It was successfully tested up to degree 32,400. All analyses and syntheses were performed, in 64 bits arithmetic and with semi-empirical control of the significant terms to prevent from calculus underflows and overflows, according to IEEE limitations, also in preserving the speed of a specific regular grid processing scheme. Finally, the continuation from the reference ellipsoid’s surface to the Earth’s surface was performed by high-order Taylor expansion with all grids of required partial derivatives being computed in parallel. The main application was the production of a 1′ × 1′ equiangular global Bouguer anomaly grid which was computed by spherical harmonic analysis of the Earth’s topography–bathymetry ETOPO1 data set up to degree and order 10,800, taking into account the precise boundaries and densities of major lakes and inner seas, with their own altitude, polar caps with bedrock information, and land areas below sea level. The harmonic coefficients for each entity were derived by analyzing the corresponding ETOPO1 part, and free surface data when required, at one arc minute resolution. The following approximations were made: the land, ocean and ice cap gravity spherical harmonic coefficients were computed up to the third degree of the altitude, and the harmonics of the other, smaller parts up to the second degree. Their sum constitutes what we call ETOPG1, the Earth’s TOPography derived Gravity model at 1′ resolution (half-wavelength). The EGM2008 gravity field model and ETOPG1 were then used to rigorously compute 1′ × 1′ point values of surface gravity anomalies and disturbances, respectively, worldwide, at the real Earth’s surface, i.e. at the lower limit of the atmosphere. The disturbance grid is the most interesting product of this study and can be used in various contexts. The surface gravity anomaly grid is an accurate product associated with EGM2008 and ETOPO1, but its gravity information contents are those of EGM2008. Our method was validated by comparison with a direct numerical integration approach applied to a test area in Morocco–South of Spain (Kuhn, private communication 2011) and the agreement was satisfactory. Finally isostatic corrections according to the Airy model, but in spherical geometry, with harmonic coefficients derived from the sets of the ETOPO1 different parts, were computed with a uniform depth of compensation of 30?km. The new world Bouguer and isostatic gravity maps and grids here produced will be made available through the Commission for the Geological Map of the World. Since gravity values are those of the EGM2008 model, geophysical interpretation from these products should not be done for spatial scales below 5 arc minutes (half-wavelength).     

On the existence of solutions for the method of Bjerhammar in the continuous case

The method of Bjerhammar is studied in the continuous case for a sphere. By varying the kernel function, different types of unknowns (u^*) are obtained at the internal sphere (the Bjerhammar sphere). It is shown that a necessary condition for the existence of u^* is that the degree variances (σ _n ² ) of the observations are of an order less than n⁻². According to Kaula’s rule this condition is not satisfied for the earth’s gravity anomaly field (σ _n ² =n⁻¹) but well for the geopotential (σ _n ² =n⁻³).     

远区重力数据对扰动位径向二阶梯度的影响

提出了利用地面重力异常数据计算地面扰动位径向二阶梯度时，远距离带区影响的误差公式，结果以级数形式给出，以其数学模型为基础，利用全球重力位模型GPM98，计算确定了不同平均异常网格大小情况下截断半径取一系列值时的远区影响，在将卫星重力梯度和航空重力梯度观测数据与地面重力数据进行联合处理时，文中提出了远区影响的数量关系，对于确定合理的截断半径大小具有一定的参考价值。     

A numerical study of the divergence of spherical harmonic series of the gravity and height anomalies at the earth's surface

The problem of the divergence of the geopotential spherical harmonic series at the earth's surface is investigated from a numerical, rather than a theoretical, approach. A representative model of the earth's potential is devised on the basis of a density layer, which, in the spherical approximation, generates a gravity field whose harmonic constituents decay according to an accepted degree variance model. This field, expanded to degree 300, and a topographic surface specified to a corresponding resolution of 67 km are used to compute the differences between truncated inner and outer series of the gravity and height anomalies at the surface of the earth model. Up to degree 300, these differences attain RMS values from 0.33 μgal to 86 μgal for the gravity anomaly and from 0.32 μm to 410 μm for the height anomaly, in areas ranging respectively from near the equator to the vicinity of the pole. In addition to these values, there is an expected truncation effect, caused by the neglect of higher degree components of the inner series, of about 30 mgal and 36 cm, respectively. The field is then subjected to a Gaussian filter which effectively cuts off information at degree 300 (at the 5% level). The RMS error to degree 300 is thereby reduced by factors of 10 to 20, with a concomitant reduction in the truncation effect to about 0.3 mgal and 0.7 cm.     

Short-wavelength Spectral Properties of the Gravity Field from a Range of Regional Data Sets

The GRACE (gravity recovery and climate experiment) and GOCE (gravity field and steady-state ocean circulation explorer) dedicated gravity satellite missions are expected to deliver the long-wavelength scales of the Earth’s gravity field with extreme precision. For many applications in Earth sciences, future research activities will have to focus on a similar precision on shorter scales not recovered by satellite missions. Here, we investigate the signal power of gravity anomalies at such short scales. We derive an average degree variance and power spectral density model for topography-reduced gravity anomalies (residual terrain model anomalies and de-trended refined Bouguer anomalies), which is valid for wavelengths between 0.7 and 100  km. The model is based on the analysis of gravity anomalies from 13 test regions in various geographical areas and geophysical settings, using various power spectrum computation approaches. The power of the derived average topography-reduced model is considerably lower than the Tscherning–Rapp free air anomaly model. The signal power of the individual test regions deviates from the obtained average model by less than a factor of 4 in terms of square-root power spectral amplitudes. Despite the topographic reduction, the highest signal power is found in mountainous areas and the lowest signal power in flat terrain. For the derived average power spectral model, a validation procedure is developed based on least-squares prediction tests. The validation shows that the model leads to a good prediction quality and realistic error measures. Therefore, for least-squares prediction, the model could replace the use of autocovariance functions derived from local or regional data.     

小波分析在重力界面反演中的应用

基于Parker-Oldenburg反演方法，分析了低通滤波器B（K）的不同阈值对位场分离的影响。分别采用小波变换和低通滤波对青藏高原地区布格重力异常进行了位场分解和莫霍面反演，并将反演结果与地震资料作比较。结果表明，在重力位场分解中小波多尺度分解可以代替低通滤波器B（K），避免了闻值的选取。     

The gravity field model IGGT_R1 based on the second invariant of the GOCE gravitational gradient tensor

Based on tensor theory, three invariants of the gravitational gradient tensor (IGGT) are independent of the gradiometer reference frame (GRF). Compared to traditional methods for calculation of gravity field models based on the gravity field and steady-state ocean circulation explorer (GOCE) data, which are affected by errors in the attitude indicator, using IGGT and least squares method avoids the problem of inaccurate rotation matrices. The IGGT approach as studied in this paper is a quadratic function of the gravity field model’s spherical harmonic coefficients. The linearized observation equations for the least squares method are obtained using a Taylor expansion, and the weighting equation is derived using the law of error propagation. We also investigate the linearization errors using existing gravity field models and find that this error can be ignored since the used a-priori model EIGEN-5C is sufficiently accurate. One problem when using this approach is that it needs all six independent gravitational gradients (GGs), but the components \(V_{xy}\) and \(V_{yz}\) of GOCE are worse due to the non-sensitive axes of the GOCE gradiometer. Therefore, we use synthetic GGs for both inaccurate gravitational gradient components derived from the a-priori gravity field model EIGEN-5C. Another problem is that the GOCE GGs are measured in a band-limited manner. Therefore, a forward and backward finite impulse response band-pass filter is applied to the data, which can also eliminate filter caused phase change. The spherical cap regularization approach (SCRA) and the Kaula rule are then applied to solve the polar gap problem caused by GOCE’s inclination of \(96.7^{\circ }\). With the techniques described above, a degree/order 240 gravity field model called IGGT_R1 is computed. Since the synthetic components of \(V_{xy}\) and \(V_{yz}\) are not band-pass filtered, the signals outside the measurement bandwidth are replaced by the a-priori model EIGEN-5C. Therefore, this model is practically a combined gravity field model which contains GOCE GGs signals and long wavelength signals from the a-priori model EIGEN-5C. Finally, IGGT_R1’s accuracy is evaluated by comparison with other gravity field models in terms of difference degree amplitudes, the geostrophic velocity in the Agulhas current area, gravity anomaly differences as well as by comparison to GNSS/leveling data.     

Generation and study of satellite gravity over Gujarat,India and their possible correlation with earthquake occurrences

High-resolution satellite gravity data have been generated and utilized to infer subsurface geological structures in the area of devastating earthquake that struck the Bhuj region in Gujarat on 26 January 2001. Latitudinal gravity profiles have been generated in the Bhuj, Anjar and IBF regions across the epicentres (23.5° N, 69.8° E/M_w 7.0 in 2001; 23.2° N, 70° E/M_w 7.0 in 1956; 24.2° N, 69.2° E/M_w 7.8 in 1819). Substantial differences in gravity anomaly patterns as high as 37 mGal could be observed existing near the epicentre regions. These gravitational differences might have caused due to the plate tectonic processes and due to the changes in densities of different lithospheric zones/sedimentary layers. Temporal variations of the satellite-derived gravity and their probable relations with already occurred major earthquakes in this region have been studied. Hence we conclude that drastic changes in gravity anomalies can be considered as a precursor for occurrences of substantially large earthquakes.     

基于逆Vening-Meinesz公式的测高重力中央区效应精密计算

为提高利用逆Vening-Meinesz公式反演测高重力中央区效应的精度,视中央区为矩形域,将垂线偏差分量表示成双二次多项式插值形式,引入非奇异变换,推导出了重力异常的计算公式。以低纬度区域2′×2′的垂线偏差实际数据为背景场进行了计算,结果表明,当中央区包含4个网格时,传统公式与推导出的重力异常计算公式误差的最大值大于1 mGal。推导出的公式可为高精度测高重力中央区效应的计算提供理论依据。     

Determination of geoidal heights and deflections of the vertical for the hellenic area using heterogeneous data

Mean gravity anomalies, deflections of the vertical, and a geopotential model complete to degree and order180 are combined in order to determine geoidal heights in the area bounded by [34°≦ϕ≤42°, 18°≦λ≦28°]. Moreover, employing point gravity anomalies simultaneously with the above data, an attempt is made to predict deflections of the vertical in the same area. The method used in the computations is least squares collocation. Using empirical covariance functions for the data, the suitable errors for the different sources of observations, and the optimum cap radius around each point of evaluation, an accuracy better than±0.60m for geoidal heights and±1″.5 for deflections of the vertical is obtained taking into account existing systematic effects. This accuracy refers to the comparison between observed and predicted values.     

基于FFT的精确高程异常估计问题

本文讨论了FFT估计高程异常及Romberg算法解求截断系数的有关问题;并利用雪林-拉普阶方差估计FFT法的误差值,指出ψ_0=1°.0的球帽是最佳选择。     

中国西太平洋海域1′×1′垂线偏差模型及精度评估

联合多种测高数据和重力异常数据,设计了观测点距离和测高精度融合的定权方法,采用最小二乘方法和Vening-Meinesz公式,分别构建了西太平洋海域(0°N—40°N,105°E—145°E)1′×1′网格化垂线偏差数字模型。选取两个不同特征区域,将垂线偏差的两个数字模型和EGM2008模型三者进行相互比较分析。结果表明:卯酉分量η的均方根差大于子午分量ξ的均方根差,海底地形复杂的南海特征区域的垂线偏差均方根差大于西太平洋中部的均方根差,构建的两个垂线偏差模型总体均方根差优于1.6″。     

Detailed gravimetric geoid for the North Atlantic

A detailed gravimetric geoid in the North Atlantic Ocean, named DGGNA-77, has been computed, based on a satellite and gravimetry derived earth potential model (consisting in spherical harmonic coefficients up to degree and order 30) and mean free air surface gravity anomalies (35180 1°×1° mean values and 245000 4′×4′ mean values). The long wavelength undulations were computed from the spherical harmonics of the reference potential model and the details were obtained by integrating the residual gravity anomalies through the Stokes formula: from 0 to 5° with the 4′×4′ data, and from 5° to 20° with the 1°×1° data. For computer time reasons the final grid was computed with half a degree spacing only. This grid extends from the Gulf of Mexico to the European and African coasts. Comparisons have been made with Geos 3 altimetry derived geoid heights and with the 5′×5′ gravimetric geoid derived byMarsh andChang [8] in the northwestern part of the Atlantic Ocean, which show a good agreement in most places apart from some tilts which porbably come from the satellite orbit recovery.