The paper presents a method of evaluating the terrain correction integral using the Fast Fourier Transform. The method requires height data on a regular grid and produces terrain corrections on all grid points. For a 1 km grid spacing, the accuracy is generally better than 1.5 mgal for typically rough areas and is rather insensitive to errors in the data. The required CPU time is proportional to NlogN, where N is the number of grid points. This paper also discusses the covariance function computation by Fourier techniques. The method is most suitable for application in the solution of geodetic boundary value problems, but it can also be used for other types of geodetic or geophysical problems involving terrain corrections.

The precise numerical evaluation of solutions of geodetic boundary value problems in rough topography depends on the proper use of gravity terrain reductions. In order to be able to efficiently handle the large amounts of gridded height data available today, formulas and computer software were developed based on the Fast Fourier Transform (FFT) algorithm. The theoretical problem with such techniques is that the kernel functions involved are often singular. Practically, this leads to a decrease in accuracy and efficiency due to the approximations adopted to avoid the singularity. Moreover, after linearizing the nonlinear integrals, only the linear and/or quadratic terms lend themselves to numerical evaluation. In this paper, the planar gravity terrain effect integrals are treated from the beginning without the usual approximations, except for the assumption of constant density which was kept for the sake of simplicity. New rigorous, simpler and more efficient formulas are developed, based on the successive application of Molodensky’s vertical derivative operator to powers of heights. Terms of any order can now be efficiently computed by FFT using the newly developed formulas, giving the effect of the terrain on gravity, gradiometry, deflections of the vertical and undulations of the geoid.

The multi-band spherical FFT approach is a generalization of the Strang van Hess spherical FFT geoid prediction method, allowing virtually error-free spherical FFT solutions through utilization of continously merged “stripes” of transforms. In the paper the errors of the planar and spherical FFT methods are investigated by comparing the geoids predicted from spherical harmonic models to geoids obtained by transforming the corresponding gravity anomalies. Test areas selected include northern areas extending to 85° N, representing all of the North Atlantic/Northern Europe region, and even for such a large area the multi-band FFT method yields very satisfactory results. As an example of practical applications the method has been used to predict gravimetric geoids in Greenland, using the spherical FFT routines for both the geoid prediction step and the computation of terrain effects. The ability to also transform height data directly from a geographical grid yields a significant ease in practical computations, especially in connection with major international geoid projects, where the handling of height grids in a multitude of national or UTM map projections might be very cumbersome.

This paper studies the Fast Hartley Transform (FHT) and uses it to compute discrete two dimensional convolutions. For real convolutions, especially when one of the two functions is even, the FHT, in addition to having all the advantages of the Fast Fourier Transform (FFT), can deal with twice the amount of data simultaneously and save one-third to one half of the computer time compared with the FFT. By using the FHT to compute geoid undulations and terrain corrections, numerical examples indicate that the RMS error introduced by FHT is 0.000 m and 0.000 mGal, respectively, as compared to the results by numerical integration.

The stochastic calibration of low-cost and consumer grade inertial sensors has recently become very important due to their wide-spread utilization in a multitude of mass-market applications like smartphone and drone navigation. The reason behind this is because if accurate stochastic modeling about the inertial sensor noise is obtained, then the estimation quality of the navigation solution may improve significantly. Generally, the mainstream methods for stochastic calibration consider only a single signal, collected under static conditions, to infer that knowledge. However, it has been observed that even though the stochastic model structure that characterizes each (static) calibration signal remains the same, its parameter values vary from one replicate to another. Even though techniques have been recently proposed to address this in a statistically efficient way, a very important factor has been neglected, namely the influence of outliers on the estimation process. In this paper, a robust multi-signal framework for the stochastic modeling of inertial sensor errors is proposed, which contains two layers of robustness: one that reduces the influence of outliers in each observed signal (data corruption) and one that safeguards the estimation process from the collection of calibration signal replicates with notably different stochastic behaviour compared to the majority (sample contamination). Furthermore, two estimators are defined from this framework, with each encompassing either one or both layers of robustness, and their efficiency in different data contamination scenarios is assessed in a simulation setting. Finally, real data collected from a consumer-grade MEMS-based device are used within a navigation simulator to evaluate the relationship between the quality of the stochastic models obtained by the two robust estimators in different data collection scenarios and the navigation solution stability.

The computation of gravity terrain corrections is implemented in this paper by the three-dimensional fast Fourier transform (3D FFT) method. By using density values on a 3D grid, a 3D grid of terrain corrections is produced from which the terrain corrections of the points on the Earth's surface are evaluated by interpolation. The technique gives directly the results at the geoid level, i.e., the indirect effect of the topographic reduction, and at a flight level, which will find a very important application in the upcoming measurement systems of airborne gravimetry and gradiometry.Numerical results by the 3D FFT are compared with those by the linear 2D FFT and numerical integration methods (NIM), in terms of computational accuracy and required time, with different interpolating methods. It is indicated by comparisons that more accurate numerical results can be achieved by the 3D FFT than by the linear 2D FFT when the grid spacing in the z-direction is small. Although it is possible to get comparable results with the 2D FFT method by evaluating more terms in the Taylor series expansion, there are two situations in which the 2D FFT method cannot be applied. These are (a) when the density varies in the z-direction and (b) when there are large terrain inclinations. In these cases, the 3D FFT method is the only efficient alternative to numerical integration. Besides the feasibility and the accuracy of the 3D FFT method, the paper also discusses some strategies for minimizing its required computer memory and CPU time.

The study of global-scale geophysical signals requires the modification of conventional spectral analysis and signal processing techniques from the real line to the sphere. These techniques often depend on the use of window functions (e.g., for localized spectral analysis and to improve the detection of periodic constituents). Normalized window functions are also utilized as averaging filters.In this work, we only focus on polynomial window functions. We present some families of polynomial windows that have been used in conventional signal processing, such as the B-spline, Singla-Singh, Kulkarni-type and generalized adaptive polynomial windows. We also demonstrate the possibility of approximating more sophisticated non-polynomial windows, such as the Kaiser, Lanczos and hyperbolic cosine windows, using their Taylor series expansion. The approach followed for their adaptation to the sphere results in isotropic (i.e., rotationally symmetric) window functions. We also examine their related filter kernels and provide expressions for their representation in the spatial domain.Recent advances on the evaluation of spherical harmonic coefficients of polynomial functions also enable us to assess the spectral characteristics of all window functions and filter kernels examined. We compare their main spectral characteristics, such as the main lobe width, first side lobe level and side lobe decay rate. Since all of these windows and filters have not been examined on the sphere before, the present work extends the current methods for localizing and filtering geophysical signals on the sphere.

Pose detection of objects is an important topic in object-level mapping and indoor localization. In the past, pose estimation methods were performed either with the help of artificial markers or natural features found on the object. However, due to the fact that the markers can only be utilized in controlled environment experiments, the application of marker-based approaches is very limited. Furthermore, methods that depend on the object's natural visual features require texture on the object and lack robustness to illumination and camera viewpoint variations. With the advent of Deep Learning (DL), the classical pose estimation methods have been outperformed. The DL-based pose estimation can detect deep features of the object and exhibits higher robustness to many distortions and variabilities caused by the changes in the illumination and viewpoint conditions. However, the massive training data set requirement is the main challenge with most DL-based methods. The training set is often a real set of images that have been manually labeled or annotated. In addition, such methods face problems related to the degradation of their predicted accuracy in the presence of uncertainties due to the symmetrical structure of many objects. To address the aforementioned issues, a novel and very fast method for generating synthetic data, as well as a contour-based technique for accurate pose estimation (that can handle pose ambiguities for a symmetrical object) are proposed in this paper. The tests that are conducted in multiple indoor scenarios demonstrate not only the effectiveness of the synthetic data generation but also exhibit, in many cases, the very high accuracy of the proposed pose estimation method.

Various aspects of gravity field modeling rely upon analytical mathematical functions for calculating spherical harmonic coefficients. Such functions allow quick and efficient evaluation of cumbersome convolution integrals defined on the sphere. In this work, we present a new analytical method for determining spherical harmonic coefficients of isotropic polynomial functions. This method in computationally flexible and efficient, since it makes use of recurrence relations. Also, its use is universal and could be extended to piecewise polynomials and polynomials with compact support. Our numerical investigation of the proposed method shows that certain recurrence relations lose accuracy as the order of implemented polynomials increases because of accumulation of numerical errors. Propagation of these errors could be mitigated by hybrid methods or using extended precision arithmetic. We demonstrate the relevance of our method in gravity field modeling and discuss two areas of application. The first one is the design of B-spline windows and filter kernels for the low-pass filtering of gravity field functionals (e.g., GRACE Follow-On monthly geopotential solutions). The second one is the calculation of spherical harmonic coefficients of isotropic polynomial covariance functions.

Modeling a Global Navigation Satellite System (GNSS) receiver clock instability is crucial to advancing GNSS-based navigation systems, especially regarding low-cost devices. This work proposes a way to obtain an insight to the drift that characterizes Local Oscillators (LOs) of GNSS receivers by using a raw GNSS baseband signal, namely the instantaneous code phases and then model its random behaviour using a state-of-the-art stochastic modeling framework, the Generalized Method of Wavelet Moments (GMWM). First, a GNSS Software-Defined Radio (SDR) is used to process three types of IntermediateFrequency (IF) GPS signals (i.e., L1 C/A, L5 data, and L5 pilot) and generate the respective code phase measurements, which are synthesized by using the LO from a GNSS frontend. Then, a pre-processor computes clock bias errors by subtracting the code phases from their references. Finally, the single-difference results of the pre-processor outputs, meaning the first order change increments of the previously computed clock bias errors, are inputted into the Robust GMWM (RGMWM) framework, which offers a stochastic modeling solution that is partially protected by the influence of outliers in the data at hand. The suggested methodology was tested by carrying out real-world experiments under open sky conditions, which validated the effectiveness of the RGMWM in stochastically modeling the LO instability using the instantaneous code phase of GNSS signals. In addition, a new type of sinusoidal noise has been identified in the input data computed from the GPS L1 C/A tracking and quantified via the RGMWM, while the L5 signals appeared to be unaffected by such noise type. This finding has the potential to contribute in the development of more advanced GNSS receivers by leveraging this more reliable modeling in either the baseband processors or navigators in the future.

In gravity field modeling, covariance functions are mainly associated with least squares collocation. Prior to the implementation of least squares collocation, the characteristics of the selected analytical covariance function need to be well understood. In this contribution, we study four polynomial covariance functions, i.e., the spherical, Askey, C^2 -Wendland and C^4 -Wendland models. All of them are defined on the sphere and correspond to isotropic, positive definite and compactly supported functions. We examine them in the spatial and spectral domains, and assess their characteristics, such as the correlation length, the curvature parameter, the spectral maximum and the spectral decay rate. We also provide analytical expressions and numerical estimates for these parameters.

The use of low-cost inertial sensors is nowadays wide-spread in many different mass-market applications, especially for navigation, for example in drones and smartphones. However, ensuring their performance is challenging since it is highly dependent on the available knowledge for the inertial sensor random error behavior, and which is hard to obtain accurately in practice. The main reason is that in many cases, the inertial sensor measurements collected during calibration contain outliers caused by either external (e.g., vibrations) or internal factors (e.g., ageing of the IMU). Therefore, it is essential that the modeling of that behavior is conducted by an estimator that has the capability to effectively reduce the influence of potential outliers from its estimation product (robustness), such that the estimated model parameters, when supplied to the chosen navigation algorithm, lead to optimal performances. The current state-of-the-art to handle the intricate nature of the stochastic errors, typically of low-cost inertial sensors, is the Generalized Method of Wavelet Moments (GMWM) and its Multi Signal extension (MS-GMWM). Although the GMWM possesses such a robustness feature thanks to an M-estimator for its fundamental quantity, the wavelet variance, it can be difficult to detect whether the analyzed data contain outliers or not, while this feature has not been extended to the multi signal approach. In this paper, it is demonstrated through simulations, that the utilization of the robust version of the GMWM in every scenario is a worthwhile trade-off between reduction of the outlier impact and reduction of the estimator’s efficiency. Furthermore, two new robust estimators in the context of the MS-GMWM are proposed and their performance is evaluated: the first one is able to decrease the influence of outliers in the analyzed data, while the second reinforces protection against calibration signal replicate(s) that present(s) significantly different stochastic error behavior compared to the others.

Gravity field and steady-state Ocean Circulation Explorer (GOCE) data are strongly affected by noise and long-wavelength errors outside the satellite measurement bandwidth (MBW). One of the main goals in utilizing GOCE data for gravity field modeling is the application of filtering techniques that can remove gross errors and reduce low-frequency errors and high-frequency noise while preserving the original signal. This paper aims to present and analyze three filtering strategies used to de-noise the GOCE Level 2 data from long-wavelength correlated errors and noise. These strategies are Finite Impulse Response (FIR), Infinite Impulse Response (IIR), and Wavelet Multi-resolution Analysis (WL), which have been applied to GOCE residual second order derivatives of the gravity potential. Several experiments were performed for each filtering scheme in order to identify the ideal filtering parameters. The outcomes indicate that all the suggested filtering strategies proved to be effective in removing low-frequency errors while preserving the signals in the GOCE MBW, with FIR filtering providing the overall best results.

The isotropic Gaussian filter has been used extensively in Gravity Recovery and Climate Experiment (GRACE) temporal gravity field solutions, and is still being applied to GRACE Follow-On products to remove high-frequency errors and improve the estimation of mass transport events on the Earth’s surface. For such applications, the only known rigorous method to calculate the spherical harmonic coefficients of an isotropic Gaussian filter is by the use of a second-order recurrence relation. As an alternative, an approximate expression is also used frequently. In this paper, we provide some additional expressions for the calculation of isotropic Gaussian filter kernels in the spherical harmonic domain. Specifically, we derive a new recurrence relation, a closed-form expression, expressions involving modified Bessel functions of the first kind, and a new approximate expression. We also examine and compare them from a computational viewpoint. The results of our numerical investigations indicate that the new recurrence relation and the closed-form expression are unstable in a way similar to the second-order recurrence relation that has been used so far. The expressions involving modified Bessel functions, and particularly the ones using exponentially scaled modified Bessel functions, provide a simple, elegant and stable way of calculating isotropic Gaussian filter coefficients, since routines for their stable evaluation are readily available in many programming languages. Alternatively, the new approximate expression can be used, which is also stable and offers better accuracy than previous approximations.

Changes in the density of the shallow crust has been previously related to co-seismic strain release during earthquakes, however, the influence of inter-seismic deformation on crustal density variations is poorly understood. Here we present gravity observations from the iGrav superconducting gravimeter in southern Vancouver Island, British Columbia, Canada which reveal a substantial gravity increase between July 2012 and April 2015. We identify a negative correlation between this gravity increase and crustal dilatation strain derived from horizontal GPS velocities. The overall increasing gravity trend is caused by the gravity increase during and immediately before and after episodic tremor and slip events, which is partially compensated by gravity decrease occurring between the events. We conclude that the observed gravity increase results from a density increase due to crustal compression and that this is mostly a result of inter-seismic strain accumulation during the subduction of the Juan de Fuca plate beneath the North American plate.

In this paper, we discuss some methods for the calculation of the Legendre polynomial difference P_n - 1(t) - P_n + 1(t) . Such differences frequently appear when calculating the spherical harmonic coefficients of truncated spatial filters. By examining the relation between the Legendre polynomials and other classical orthogonal polynomials (i.e., Gegenbauer and Jacobi polynomials), we provide analytical expressions for P_n - 1(t) - P_n + 1(t) in terms of the latter. We also derive a recurrence relation that avoids the direct evaluation of the Legendre polynomials. Analogous analytical expressions and recurrence relations are provided for the calculation of the P_n - 1(t) - P_n + 1(t) derivative and integral. We show that these recurrence relations can be derived from a more general recurrence relation; thus, the former are special cases of the latter. The expressions given in this study can be useful for the more efficient calculation of truncated kernels in the spherical harmonic domain. We demonstrate the applicability of these new expressions for the Pellinen (i.e., rectangular) filter kernel that is regularly used in regional gravity field modeling.

Abstract The Gravity Recovery and Climate Experiment (GRACE) mission has enabled mass changes and transports in the hydrosphere, cryosphere and oceans to be quantified with unprecedented resolution. However, while this legacy is currently being continued with the GRACE Follow-On (GRACE-FO) mission there is a gap of 11 months between the end of GRACE and the start of GRACE-FO which must be addressed. Here we bridge the gap by combining time-variable, low-resolution gravity models derived from European Space Agency’s Swarm satellites with the dominating spatial modes of mass variability obtained from GRACE. We show that the noise inherent in unconstrained Swarm gravity solutions is greatly reduced, that basin averages can have root mean square errors reduced to the order of $$\text {cm}$$ cm of equivalent water height, and that useful information can be retrieved for basins as small as $$1000 \times 1000\,\hbox {km}$$ 1000 × 1000 km . It is found that Swarm data contains sufficient information to inform the leading three global mass modes found in GRACE at the least. By comparing monthly reconstructed maps to GRACE data from December 2013 to June 2017, we suggest the uncertainty of these maps to be $$2{-}3\,\text {cm}$$ 2 - 3 cm of equivalent water height.

In this chapter, the author first defines the two main geodetic boundary value problems, namely Stokes' problem and Molodensky's problem. Then he summarizes theoretical solutions of the problems, which yield expressions for evaluating, followed by practical formulas for computing geoid undulations using spherical harmonic coefficients, gravity anomalies, and topographic data. Precise geoid determination over a large region is possible with a combination of spherical harmonic coefficients, gravity anomalies, and heights. The use of Stokes' Equation requires gravity anomalies all over the earth for the computation of a single geoid undulation. Obviously, this is impractical to say the least and thus, in practice, some modifications of the technique are necessary. Consequently, the use of the GPS technology with gravimetric geoids can provide vertical positioning for all sorts of mapping, positioning, and exploration projects in a very fast and cost effective manner.

The strength and direction of gravity and ground deformation anomaly vectors over normal and tight reservoirs after injection of CO2 depend on the reservoir and CO2 properties. Some of these properties, such as porosity, permeability, and size, define the reservoir type and therefore determine the existence or sign of the anomalies, indirectly. However, some other properties of reservoirs, such as reservoir depth, horizontal extension of CO2 plume, and CO2 mass amount and density, affect the strength of the measurements directly. Gravimetric and geodetic modelling of synthetic CO2 reservoirs can quantify each of the direct effects that represent the expected signals over different geological settings. We present gravity signal as the gravity increase due to CO2 mass attraction, and its combination with the gravity decrease caused by ground uplift. Our results indicate that the reservoir depth and horizontal extension, along with the injected mass, have significant influences on both ground deformation and gravity signals. The results also demonstrate that the horizontal extension of the CO2 distribution decreases the dependency on the depth. In addition, the observed gravity signal over CO2 reservoirs is dominated by the free-air effect from large ground deformation. Finally, the gravity effect over both normal and tight reservoirs and the ground surface deformation over tight reservoirs are highly dependent on the density change of the injected CO2 inside the reservoir as a result of depth dependant changes of pressure and temperature.