Open AccessArticle

High Accurate Mathematical Tools to Estimate the Gravity Direction Using Two Non-Orthogonal Inclinometers

Daniele Mortari

and

Anthony Gardner

Aerospace Engineering, Texas A&M University, 746C H.R. Bright Building, 3141 TAMU, College Station, TX 77843-3141, USA

Author to whom correspondence should be addressed.

Sensors 2021, 21(17), 5727; https://doi.org/10.3390/s21175727

Submission received: 23 July 2021 / Revised: 18 August 2021 / Accepted: 23 August 2021 / Published: 25 August 2021

(This article belongs to the Special Issue Attitude Estimation Based on Data Processing of Sensors)

Download

Browse Figures

Figure 1
Geometry of the inclinometer system and associated reference frame. "> Figure 2
Geometrical explanation justifying the result provided by Equation (<a href="#FD6-sensors-21-05727" class="html-disp-formula">6</a>). "> Figure 3
Monte Carlo tests: histograms of the angles and the gravity components. "> Figure 4
Monte Carlo validation for the <math display="inline"><semantics> <mrow> <mi>P</mi> <mo>(</mo> <mn>1</mn> <mo>,</mo> <mn>3</mn> <mo>)</mo> </mrow> </semantics></math>, <math display="inline"><semantics> <mrow> <mi>P</mi> <mo>(</mo> <mn>2</mn> <mo>,</mo> <mn>3</mn> <mo>)</mo> </mrow> </semantics></math> and <math display="inline"><semantics> <mrow> <mi>P</mi> <mo>(</mo> <mn>3</mn> <mo>,</mo> <mn>3</mn> <mo>)</mo> </mrow> </semantics></math> terms. "> Figure 5
Extended Kalman Filter test results. "> Figure 6
Geocentric (<math display="inline"><semantics> <mi>θ</mi> </semantics></math>) and geodetic (<math display="inline"><semantics> <mi>ϕ</mi> </semantics></math>) latitude angles. "> Figure 7
Gravity description with respect to geoid and reference ellipsoid surfaces. "> Figure 8
East–West and North–South geoid deviation (Himalaya region). "> Figure 9
Residual <math display="inline"><semantics> <msub> <mi>L</mi> <mn>2</mn> </msub> </semantics></math> norms. "> Figure 10
Least-squares condition number. ">

Review Reports Versions Notes

Abstract

This study provides two mathematical tools to best estimate the gravity direction when using a pair of non-orthogonal inclinometers whose measurements are affected by zero-mean Gaussian errors. These tools consist of: (1) the analytical derivation of the gravity direction expectation and its covariance matrix, and (2) a continuous description of the geoid model correction as a linear combination of a set of orthogonal surfaces. The accuracy of the statistical quantities is validated by extensive Monte Carlo tests and the application in an Extended Kalman Filter (EKF) has been included. The continuous geoid description is needed as the geoid represents the true gravity direction. These tools can be implemented in any problem requiring high-precision estimates of the local gravity direction.

Keywords:

inclinometer data processing; covariance analysis; geoid correction

1. Introduction

The direct measurement of gravity using inclinometers has a variety of practical applications. Tilt sensors, such as accelerometers or inclinometers, are the hardware of choice for measuring a planet’s gravitational field. While the gravity is perturbed on a global scale by other massive bodies (e.g., Moon, Sun), local perturbations to the gravity vector’s pointing direction depend on nearby mass distributions and anomalies, or a lack thereof. A detailed review of various kinds of solid tilt sensors with comparisons between the most typical tilt sensing techniques can be found in [1].

From a sensory perspective, inclinometers are commonly used to measure inclination angle magnitudes and structural deformations. Measured values are typically provided as a percentage or angular deflection with respect to a level reference surface, whose nominal orientation is perpendicular to the local gravity vector. Other classic inclinometer use cases include measuring the slope gradient during activities such as tunneling, de-watering and excavation, as well as those that require monitoring the integrity of the ground around a structure.

In land surveying and mapping, an inclinometer can provide a rapid measurement of the slope of a geographic feature [2] or be used for analyzing the magnitude, rate, direction, depth, and type of landslides [3]. Continuous monitoring of landslides is important for early warning purposes. In general, inclinometers can be used for continuous monitoring of any movement, providing valuable information for landslide autonomous geometry reconstruction [4]. In the oil and gas industries, inclinometers are used to measure the strike and dip of geologic formations [5] in order to detect oil and mineral deposits. In forestry, tree height measurements can be made by an inclinometer using standardized methods [6]. Major artillery guns may have an associated inclinometer used to facilitate the aiming of shells over long distances [7,8]. Permanently installed inclinometers are employed at major earthworks such as dams [9,10] to monitor the long-term stability of the structure. Inclinometers can also be used in robotics [11] and to measure bridge deflections [12] for bridge safety evaluation. Specifically, the authors of [13] developed and implemented a measuring system to assess railway bridges using measurements of the structural response to passing trains. Furthermore, low-cost tilt sensors are now commonplace in commercially available smartphones in order to make relative orientation estimates of the device.

Inclinometers are also used in more complex systems, such as in problems of estimating the position and orientation of a calibrated camera observing a set of n 3D points, known as Perspective-n-Point (PnP) problems [14], and in “Stellar Positioning Systems” [15,16,17], where the gravity direction along with celestial references can be used to estimate the geographical position in scenarios where real-time GPS information is unavailable (e.g., on Mars). Such systems can also be used on Earth as a backup system during instances of GPS jamming or spoofing.

Today, new types of inclinometers are emerging, such as the two-dimensional highly sensitive fiber-optic inclinometer proposed by [18], which can provide long-term continuous monitoring of inclinations or real-time feedback control of tilt angles, especially in harsh environments with violent temperature variations. Microelectromechanical (MEMS) inclination sensors are also being developed to enhance the properties of design, fabrication, and signal measurement precision [19].

Prior to being deployed in the field, inclinometers are efficiently calibrated in a laboratory setting. Because of vibrations and subsequent thermal fluctuations, the orientation of the inclinometer assembly degrades over time, requiring re-calibration methods to properly account for new measurement errors. Reference [20] proposes a new inclinometer assembly error calibration and horizontal image correction method utilizing plumb lines, where the horizontal correction is derived by a homography matrix.

The monitoring of possible road subsidence is proposed using a combination of theodolite measurements and inclinometers [21] to detect ground motion near critical infrastructure during excavation in the vicinity. The proposed method, validated by Monte Carlo simulation, is formulated as a covariance weighted Hermite approximation problem. On the contrary, this paper addresses the problem of the best estimation of the gravity direction by non-approximated (analytical) derivations for a pair of inclinometers only. The proposed method can then be adopted in almost all the scenarios involving two inclinometers.

This study considers the classic system of two non-orthogonal inclinometers affected by zero-mean Gaussian error (hypothesis) to measure the gravity direction, and it provides the non-approximated mathematical tools to best estimate the gravity direction for any static and/or dynamic scenario (objective). In particular, we provide:

The analytical tools to best estimate the gravity direction. These tools consist of: (a) an exact estimate of the gravity direction when inclinometer measurements are affected by zero-mean Gaussian noise, and (b) the associated gravity covariance matrix. These tools can be used in static and filtered dynamic scenarios to accurately estimate the actual gravity direction and to model the local subtle gravity variations.
A continuous mathematical model to describe the true gravity direction deviation with respect to the latest (1984) revision of the approximated World Geodetic System (WGS-84) ellipsoidal model. Two deflection correction models are proposed for the north–south and east–west deflection variations. These models are expressed using combinations of a set of N linearly independent orthogonal surfaces that are derived from Chebyshev orthogonal polynomials of the first kind.

These tools and models find relevance in a variety of gravity measurement tasks and associated use, especially in filtered systems, such as the one developed for geo-technical monitoring in [22].

System Definition

The system under consideration in this study is shown in Figure 1, which describes the axes of the two non-orthogonal inclinometers,

[I_{x}, I_{y}]

, and an orthonormal inclinometer reference frame,

{x, y, z}

. The

I_{x}

-axis of the inclinometer frame is assumed to be coincident with the

x

-axis reference frame, and the plane defined by the two inclinometers axes is also coincident with the

[x, y]

-plane. The relative orientation of the two inclinometer axes deviates from the nominal/optimal orthogonal displacement by an angle

ε

. This angle, which takes into account the small deviation from orthogonality due to mounting errors, can be estimated during lab calibration and, possibly, re-calibrated later to account for variations caused by vibrations and/or thermal expansions.

The gravity direction in the

[x, y, z]

frame is

g = \{\begin{matrix} g_{x} \\ g_{y} \\ g_{z} \end{matrix}\} = \{\begin{matrix} cos ϑ_{x} \\ \frac{cos ϑ_{y} - cos ϑ_{x} sin ε}{cos ε} \\ - \sqrt{1 - g_{x}^{2} - g_{y}^{2}} \end{matrix}\},

(1)

where the third component enforces the unit-norm constraint. The inclinometer measurement angles,

ϑ_{x}

and

ϑ_{y}

(see Figure 1), allow the estimation of the gravity direction in the reference frame.

2. Gravity Direction Estimation

This section contains the covariance analysis of the gravity direction estimated by two identical inclinometers, as shown in Figure 1. The analysis is performed for the generic case of a non-orthogonal mounting, where the inclinometer axes differ from orthogonality by the angle

ε

Let the two measured angles,

ϑ_{x}

and

ϑ_{y}

, be affected by Gaussian noise

ϑ_{x} \sim N (μ_{x}, σ^{2}) and ϑ_{y} \sim N (μ_{y}, σ^{2}) .

where

μ_{x}

and

μ_{y}

are the mean values and

σ^{2}

is the measurement variance, identical for the two inclinometers. The mean values,

μ_{x}

and

μ_{y}

, constitute the two main components of the EKF state vector. The quantities

ε

and

σ^{2}

are initially estimated by lab calibration methods. However, vibrations and/or thermal expansions will likely change their values over time. In this case, these parameters should be re-calibrated and their values re-estimated by adding them as new variables to the EKF state vector (see the Discussion section).

The angles,

ϑ_{x}

and

ϑ_{y}

, can be represented by a zero-mean Gaussian angle,

δ

, as

\{\begin{matrix} ϑ_{x} = μ_{x} + δ \\ ϑ_{y} = μ_{y} + δ \end{matrix}, where δ \sim N (0, σ^{2}) .

To estimate the gravity direction, Equation (1) requires the evaluation of

E {cos ϑ_{x}}

and

E {cos ϑ_{y}}

. Dropping the axis-specific subscript, we can write,

E {cos ϑ} = E {cos (μ + δ)} = cos μ \cdot E {cos δ} - sin μ \cdot E {sin δ} .

(2)

The term

sin δ

is an odd function. Therefore, its expectation is zero. Hence, Equation (2) simplifies to

E {cos ϑ} = cos μ \cdot E {cos δ} .

To compute

E {cos δ}

, let us expand it by Maclaurin series [23]:

E {cos δ} = E \{\sum_{k = 0}^{\infty} \frac{{(- 1)}^{k}}{(2 k)!} δ^{2 k}\} = \sum_{k = 0}^{\infty} \frac{{(- 1)}^{k}}{(2 k)!} E \{δ^{2 k}\} .

(3)

Reference [24] provides the following identity:

E \{δ^{2 k}\} = 1 \cdot 3 \cdot 5 \cdot \dots \cdot (2 k - 1) σ^{2 k} = \frac{(2 k)!}{2^{k} k!} σ^{2 k} .

(4)

Substituting Equation (4) into Equation (3), the exact expectation of the cosine of the zero-mean Gaussian angle,

δ

, is

E {cos δ} = \sum_{k = 0}^{\infty} \frac{{(- 1)}^{k}}{2^{k} k!} σ^{2 k} = \sum_{k = 0}^{\infty} \frac{{(- 1)}^{k}}{k!} {(\frac{σ^{2}}{2})}^{k} = e^{- σ^{2} / 2} .

Therefore, we obtain

\{\begin{matrix} E {cos ϑ_{x}} = e^{- σ^{2} / 2} cos μ_{x} \\ E {cos ϑ_{y}} = e^{- σ^{2} / 2} cos μ_{y} \end{matrix} .

(5)

This result allows us to find the expected values of the first two components of the gravity direction:

\begin{matrix} {\hat{g}}_{x} & = E {g_{x}} = e^{- σ^{2} / 2} cos μ_{x} \end{matrix}

(6)

\begin{matrix} {\hat{g}}_{y} & = E {g_{y}} = e^{- σ^{2} / 2} (\frac{cos μ_{y}}{cos ε} - cos μ_{x} tan ε) . \end{matrix}

(7)

Figure 2 provides a geometrical interpretation of why the expectation of the cosine is lower than the cosine of the mean when the measured angle is affected by zero-mean error. If the inclinometer is not biased, then the error of the inclinometer can be described by a small cone whose axis is coincident with the inclinometer axis. This implies that inclinometer measurements affected by errors will fall within this cone. If the true angle with respect to the gravity direction is

ϑ

, then all measurements providing the true angle belong to the surface of a cone whose axis is the gravity direction and aperture is the angle

ϑ

. This cone intersects the measurement error cone by splitting it into two spherical areas that, as shown in Figure 2, are identified as “A” and “B.” All measurement angles falling in the “A” area provide an angle greater than the true value, while all the measurement angles falling into the “B” area provide an angle lower than the true value. Since the “A” area is greater than the “B” area, a zero-mean inclinometer will provide an angle greater than the true value. This means that a zero-mean inclinometer provides biased angles. This phenomenon is generated by our three-dimensional spherical world, and only vanishes under two situations: (a) when the measurement angle is perfectly

ϑ = 90^{\circ}

, and (b) when the inclinometer standard deviation (dictating the aperture of the measurement error cone) is

σ = 0^{\circ}

(idealized inclination sensor).

This bias effect is actually very small, even for poorly accurate commercial inclinometers. However, since the inclinometer axis direction is normally placed perpendicular to the gravity direction, this effect is even smaller, because the inclinometer works around the nominal orthogonality displacement of

ϑ = 90^{\circ}

3. Gravity Direction Covariance Matrix

The gravity direction covariance matrix is

P = E \{(g - \hat{g}) {(g - \hat{g})}^{T}\} = [\begin{matrix} E {g_{x}^{2}} - {\hat{g}}_{x}^{2} & E {g_{x} g_{y}} - {\hat{g}}_{x} {\hat{g}}_{y} & E {g_{x} g_{z}} - {\hat{g}}_{x} {\hat{g}}_{z} \\ (sym) & E {g_{y}^{2}} - {\hat{g}}_{y}^{2} & E {g_{y} g_{z}} - {\hat{g}}_{y} {\hat{g}}_{z} \\ (sym) & (sym) & E {g_{z}^{2}} - {\hat{g}}_{z}^{2} \end{matrix}] .

To compute the

P (1, 1)

P (1, 2)

, and

P (2, 2)

terms of this matrix, the following estimations are required:

\{\begin{matrix} E {g_{x}^{2}} = E \{{cos}^{2} ϑ_{x}\} \\ E {g_{x} g_{y}} = E \{cos ϑ_{x} cos ϑ_{y}\} \frac{1}{cos ε} - E \{{cos}^{2} ϑ_{x}\} tan ε \\ E {g_{y}^{2}} = \frac{E \{{cos}^{2} ϑ_{y}\} - 2 E \{cos ϑ_{x} cos ϑ_{y}\} sin ε + E \{{cos}^{2} ϑ_{x}\} {sin}^{2} ε}{{cos}^{2} ε} \end{matrix}

(8)

The term

E {cos ϑ_{x} cos ϑ_{y}}

can be expressed as

E {cos ϑ_{x} cos ϑ_{y}} = E {cos ϑ_{x}} E {cos ϑ_{y}} = cos μ_{x} cos μ_{y} e^{- σ^{2}},

(9)

because

cos ϑ_{x}

and

cos ϑ_{y}

are statistically independent, while the computation of the

E {{cos}^{2} ϑ}

terms requires the estimation of

E {{cos}^{2} δ}

and

E {{sin}^{2} δ}

. We have

E {{cos}^{2} ϑ} = E {{(cos μ cos δ - sin μ sin δ)}^{2}} = {cos}^{2} μ E {{cos}^{2} δ} + {sin}^{2} μ E {{sin}^{2} δ} .

The

E {{cos}^{2} δ}

and

E {{sin}^{2} δ}

expectations can also be computed by Maclaurin series:

\begin{matrix} E {{cos}^{2} δ} & = & E \{\frac{1}{2} + \frac{1}{2} \sum_{k = 0}^{\infty} \frac{{(- 1)}^{k} 2^{2 k}}{(2 k)!} δ^{2 k}\} = \frac{1}{2} + \frac{1}{2} \sum_{k = 0}^{\infty} \frac{{(- 1)}^{k} 2^{2 k}}{(2 k)!} E \{δ^{2 k}\} \\ = & \frac{1}{2} + \frac{1}{2} \sum_{k = 0}^{\infty} \frac{{(- 1)}^{k} 2^{2 k}}{(2 k)!} \cdot \frac{(2 k)!}{2^{k} k!} \cdot σ^{2 k} = \frac{1}{2} + \frac{1}{2} \sum_{k = 0}^{\infty} \frac{{(- 1)}^{k} 2^{k}}{k!} \cdot σ^{2 k} \\ = & \frac{1}{2} (1 + e^{- 2 σ^{2}}) \end{matrix}

and, similarly,

E {{sin}^{2} δ} = E \{\frac{1}{2} - \frac{1}{2} \sum_{k = 0}^{\infty} \frac{{(- 1)}^{k} 2^{2 k}}{(2 k)!} δ^{2 k}\} = \dots = \frac{1}{2} (1 - e^{- 2 σ^{2}}) .

Therefore,

E {{cos}^{2} ϑ} = \frac{{cos}^{2} μ}{2} (1 + e^{- 2 σ^{2}}) + \frac{{sin}^{2} μ}{2} (1 - e^{- 2 σ^{2}}) = \frac{2 {cos}^{2} μ - 1}{2} e^{- 2 σ^{2}} + \frac{1}{2},

and, specifically,

\{\begin{matrix} E {{cos}^{2} ϑ_{x}} = \frac{2 {cos}^{2} μ_{x} - 1}{2} e^{- 2 σ^{2}} + \frac{1}{2} \\ E {{cos}^{2} ϑ_{y}} = \frac{2 {cos}^{2} μ_{y} - 1}{2} e^{- 2 σ^{2}} + \frac{1}{2} \end{matrix} .

(10)

Equations (9) and (10) allow us to obtain the following expectation expressions:

\begin{matrix} E {g_{x}^{2}} & = \frac{2 {cos}^{2} μ_{x} - 1}{2} e^{- 2 σ^{2}} + \frac{1}{2} \end{matrix}

(11)

\begin{matrix} E {g_{x} g_{y}} & = \frac{cos μ_{x} cos μ_{y} e^{- σ^{2}}}{cos ε} - [\frac{2 {cos}^{2} μ_{x} - 1}{2} e^{- 2 σ^{2}} + \frac{1}{2}] tan ε \end{matrix}

(12)

\begin{matrix} E {g_{y}^{2}} & = \frac{1}{{cos}^{2} ε} [\frac{2 {cos}^{2} μ_{y} - 1}{2} e^{- 2 σ^{2}} + \frac{1}{2}] - cos μ_{x} cos μ_{y} e^{- σ^{2}} \frac{2 sin ε}{{cos}^{2} ε} \\ + [\frac{2 {cos}^{2} μ_{x} - 1}{2} e^{- 2 σ^{2}} + \frac{1}{2}] {tan}^{2} ε \end{matrix}

(13)

Using Equations (6) and (11), the term,

P (1, 1) = E {g_{x}^{2}} - {\hat{g}}_{x}^{2}

, of the covariance matrix can be written as

P (1, 1) = \frac{2 {cos}^{2} μ_{x} - 1}{2} e^{- 2 σ^{2}} + \frac{1}{2} - {(e^{- σ^{2} / 2} cos μ_{x})}^{2}

This expression can be simplified into the following more compact form:

P (1, 1) = (\frac{1 - e^{- 2 σ^{2}}}{2}) - (e^{- σ^{2}} - e^{- 2 σ^{2}}) {cos}^{2} μ_{x}

(14)

To derive the analytical expression of the term,

P (2, 2) = E {g_{y}^{2}} - {\hat{g}}_{y}^{2}

, we make use of Equations (7) and (13):

\begin{matrix} P (2, 2) & = \frac{1}{{cos}^{2} ε} [\frac{2 {cos}^{2} μ_{y} - 1}{2} e^{- 2 σ^{2}} + \frac{1}{2}] - cos μ_{x} cos μ_{y} e^{- σ^{2}} \frac{2 sin ε}{{cos}^{2} ε} \\ + [\frac{2 {cos}^{2} μ_{x} - 1}{2} e^{- 2 σ^{2}} + \frac{1}{2}] {tan}^{2} ε - e^{- σ^{2}} {(\frac{cos μ_{y}}{cos ε} - cos μ_{x} tan ε)}^{2} \end{matrix}

After some simple manipulations and simplifications, the previous equation can be written in the more compact form

P (2, 2) = \frac{1 + {sin}^{2} ε}{{cos}^{2} ε} (\frac{1 - e^{- 2 σ^{2}}}{2}) + \frac{{cos}^{2} μ_{x} {sin}^{2} ε - {cos}^{2} μ_{y}}{{cos}^{2} ε} (e^{- σ^{2}} - e^{- 2 σ^{2}})

(15)

Note that, as a sanity check, if the two inclinometers were perfectly orthogonal, meaning if

ε = 0

, then Equation (15) becomes formally identical to Equation (14).

As for the term

P (1, 2) = E {g_{x} g_{y}} - {\hat{g}}_{x} {\hat{g}}_{y}

, we have

\begin{matrix} P (1, 2) & = \frac{cos μ_{x} cos μ_{y} e^{- σ^{2}}}{cos ε} - [\frac{2 {cos}^{2} μ_{x} - 1}{2} e^{- 2 σ^{2}} + \frac{1}{2}] tan ε \\ - e^{- σ^{2}} (\frac{cos μ_{x} cos μ_{y}}{cos ε} - {cos}^{2} μ_{x} tan ε) \end{matrix}

which can be simplified to

P (1, 2) = [(e^{- σ^{2}} - e^{- 2 σ^{2}}) {cos}^{2} μ_{x} - (\frac{1 - e^{- 2 σ^{2}}}{2})] tan ε,

(16)

This equation highlights an unexpected result: the covariance term,

P (1, 2)

, is not a function of

μ_{y}

. This result, which has been numerically verified by extensive Monte Carlo tests, remains with no explanation. Note also that, if the two inclinometers were perfectly orthogonal (

ε = 0

), then we would obtain

P (1, 2) = 0

Computation of the P(1,3), P(2,3) and P(3,3) Terms via Covariance Law

The analytical expressions for the three missing terms of the covariance matrix,

P (1, 3)

P (2, 3)

and

P (3, 3)

, were not found. However, a highly accurate estimation of these terms can be obtained using the covariance law [25,26], sometimes called the error propagation law.

The covariance law connects the two covariance matrices that are expressed in terms of two different (but equivalent) parameterizations by the Jacobian,

J

, of the transformation between parameterizations. In this specific case, the covariance law connects the

3 \times 3

covariance matrix of the gravity unit-vector, P, with the

2 \times 2

covariance matrix written in terms of the cosine of the measurements, C. This relationship consists of the identity

P \approx J C J^{T}

(17)

where the analytical expressions of the covariance matrix terms,

C = [\begin{matrix} E \{{cos}^{2} ϑ_{x}\} - E {\{cos ϑ_{x}\}}^{2} & E \{cos ϑ_{x} cos ϑ_{y}\} - E \{cos ϑ_{x}\} E \{cos ϑ_{y}\} \\ (sym) & E \{{cos}^{2} ϑ_{y}\} - E {\{cos ϑ_{y}\}}^{2} \end{matrix}]

can all be computed using Equations (5), (9) and (10), while the Jacobian of the transformation between the two parameterizations is

J = [\begin{matrix} \frac{\partial g_{x}}{\partial cos ϑ_{x}} & \frac{\partial g_{x}}{\partial cos ϑ_{y}} \\ \frac{\partial g_{y}}{\partial cos ϑ_{x}} & \frac{\partial g_{y}}{\partial cos ϑ_{y}} \\ \frac{\partial g_{z}}{\partial cos ϑ_{x}} & \frac{\partial g_{z}}{\partial cos ϑ_{y}} \end{matrix}] .

The first four elements of the Jacobian can be computed from Equation (1). Their expressions are

\{\begin{matrix} \frac{\partial g_{x}}{\partial cos ϑ_{x}} = 1, \\ \frac{\partial g_{x}}{\partial cos ϑ_{y}} = 0 \end{matrix} and \{\begin{matrix} \frac{\partial g_{y}}{\partial cos ϑ_{x}} = - tan ε, \\ \frac{\partial g_{y}}{\partial cos ϑ_{y}} = - \frac{1}{cos ε} \end{matrix}

while the last two elements can be computed using the chain rule as

\{\begin{matrix} \frac{\partial g_{z}}{\partial cos ϑ_{x}} & = \frac{\partial g_{z}}{\partial g_{x}} \cdot \frac{\partial g_{x}}{\partial cos ϑ_{x}} + \frac{\partial g_{z}}{\partial g_{y}} \cdot \frac{\partial g_{y}}{\partial cos ϑ_{x}} = \frac{\partial g_{z}}{\partial g_{x}} - \frac{\partial g_{z}}{\partial g_{y}} tan ε \\ \frac{\partial g_{z}}{\partial cos ϑ_{y}} & = \frac{\partial g_{z}}{\partial g_{x}} \cdot \frac{\partial g_{x}}{\partial cos ϑ_{y}} + \frac{\partial g_{z}}{\partial g_{y}} \cdot \frac{\partial g_{y}}{\partial cos ϑ_{y}} = - \frac{\partial g_{z}}{\partial g_{y}} \cdot \frac{1}{cos ε} \end{matrix}

where

\begin{matrix} \frac{\partial g_{z}}{\partial g_{x}} & = \frac{g_{x}}{\sqrt{1 - g_{x}^{2} - g_{y}^{2}}} = \frac{cos ϑ_{x} cos ε}{\sqrt{{cos}^{2} ε - {cos}^{2} ϑ_{x} - {cos}^{2} ϑ_{y} + 2 cos ϑ_{x} cos ϑ_{y} sin ε}} \end{matrix}

(18)

\begin{matrix} \frac{\partial g_{z}}{\partial g_{y}} & = \frac{g_{y}}{\sqrt{1 - g_{x}^{2} - g_{y}^{2}}} = \frac{cos ϑ_{y} - cos ϑ_{x} sin ε}{\sqrt{{cos}^{2} ε - {cos}^{2} ϑ_{x} - {cos}^{2} ϑ_{y} + 2 cos ϑ_{x} cos ϑ_{y} sin ε}} \end{matrix}

(19)

Finally, the expression of the Jacobian of the transformation is

J = [\begin{matrix} 1 & 0 \\ - tan ε & - \frac{1}{cos ε} \\ \frac{\partial g_{z}}{\partial g_{x}} - \frac{\partial g_{z}}{\partial g_{y}} tan ε & - \frac{\partial g_{z}}{\partial g_{y}} \cdot \frac{1}{cos ε} \end{matrix}]

where the expressions of the two partials,

\frac{\partial g_{z}}{\partial g_{x}}

and

\frac{\partial g_{z}}{\partial g_{y}}

, are provided by Equation (18) and Equation (19), respectively.

The covariance law can actually be used to estimate all the elements of the covariance matrix. However, this law has been used here only to estimate the three missing elements of the covariance matrix,

P (1, 3)

P (2, 3)

and

P (3, 3)

. The main reason is that the covariance law is based on a linear Taylor expansion of the transformation and is therefore an approximated method, even if it provides a highly accurate estimation of the covariance matrix. On the contrary, the terms

P (1, 1)

P (2, 2)

and

P (1, 2)

, provided by the Equations (14)–(16), are analytically exact.

It is easy to derive them by setting

Q = C (1, 2) J (3, 1) + C (2, 2) J (3, 2)

. In doing so, the expressions of the last three terms of the covariance matrix can be written as

\{\begin{matrix} P (1, 3) = C (1, 1) J (3, 1) + C (1, 2) J (3, 2) \\ P (2, 3) = J (2, 1) P (1, 3) + J (2, 2) Q \\ P (3, 3) = J (3, 1) P (1, 3) + J (3, 2) Q \end{matrix}

(20)

4. Covariance Analysis Numerical Validation

The numerical validation of the covariance matrix was performed by Monte Carlo tests using the R2016b version of MATLAB software. All of the tests were carried out using two identical inclinometers affected by zero-mean Gaussian error with standard deviation

1 σ = 0.1

deg., and with axes differing from orthogonality by the angle

ε = 5

deg. A set of

10^{7}

Monte Carlo tests were performed on 1000 gravity directions, uniformly generated in a cone with 30 deg. (

π / 6

rad.) aperture around the nominal gravity direction

g = {0, 0, - 1}^{T}

. The uniformly distributed gravity directions were generated as

g = \{\begin{matrix} sin ϕ cos λ \\ sin ϕ sin λ \\ cos ϕ \end{matrix}\} where \{\begin{matrix} λ \sim U [0, 2 π] \\ cos ϕ \sim U [- 1, cos (π - π / 6)] \end{matrix}

where

U [a, b]

indicates a uniform distribution in the

[a, b]

range. The gravity direction components,

(g_{x}, g_{y}, g_{z})

, and the two angles,

(ϑ_{x}, ϑ_{y})

, generated by this distribution have the histograms reported in Figure 3.

Since the three terms

P (1, 1)

P (1, 2)

and

P (2, 2)

are mathematically correct (meaning, not approximated), the Monte Carlo numerical validation tests were performed simply to quantify the accuracy of the three elements of the covariance matrix,

P (1, 3)

P (2, 3)

and

P (3, 3)

, whose expressions were obtained by the covariance law.

The validation was performed by comparing the three elements of the covariance matrix,

T (1, 3)

T (2, 3)

and

T (3, 3)

, numerically estimated by the Monte Carlo tests, and their approximated analytical expressions,

P (1, 3)

P (2, 3)

and

P (3, 3)

, obtained using Equation (20).

The top three plots of Figure 4 show the histograms of

T (1, 3)

T (2, 3)

and

T (3, 3)

, while the bottom three plots show the histograms of the differences with the three terms estimated using Equation (20). The results of these Monte Carlo tests highlight the accurate estimation of the approximated terms,

P (1, 3)

P (2, 3)

and

P (3, 3)

, computed via covariance law. The values of the terms

T (1, 3)

and

T (2, 3)

, numerically estimated by the Monte Carlo tests, fall within the

[- 10^{- 6}, + 10^{- 6}]

range. The differences with the terms

P (1, 3)

P (2, 3)

, estimated via the covariance law, are in the

[- 10^{- 9}, + 10^{- 9}]

range, a difference better than three orders of magnitude with respect to the numerically estimated values. The same accuracy error level is experienced for the

P (3, 3)

term with respect to the numerically estimated

T (3, 3)

term, whose values fall within the

[0, + 10^{- 6}]

range.

4.1. Extended Kalman Filter (EKF) Validation

The variable to estimate is represented by the gravity direction, i.e., by a three-component vector. However, since the third component of the gravity vector,

g_{z}

, is not independent (actually, it is derived from

g_{x}

and

g_{y}

), then a reduced two-state EKF estimator must be developed and used. In order to quantify the accuracy gain obtained using the closed-form expression of the covariance matrix, a static scenario was selected to validate the EKF estimator. This choice does not depend on the dynamics under consideration or the accuracy of the adopted dynamical model. The EKF summary is provided as follows:

State model:

$x_{k} = {\{\begin{matrix} g_{x} \\ g_{y} \end{matrix}\}}_{k} = x_{k - 1} + w_{k - 1}$
Observation model:

$y_{k} = {\{\begin{matrix} {\tilde{ϑ}}_{x} \\ {\tilde{ϑ}}_{y} \end{matrix}\}}_{k} = h (x_{k}) + v_{k} = {cos}^{- 1} [e^{σ^{2} / 2} {\{\begin{matrix} g_{k} \\ g_{k} cos ε + g_{k} sin ε \end{matrix}\}}_{k}] + v_{k}$
Initialization: Given initial state vector, $x_{0}^{-}$ , and initial covariance matrix, $P_{0}^{-}$ , and assuming no process noise, $Q_{k} = E {w_{k} w_{k}^{T}} = 0_{2 \times 2}$ , the measurement’s covariance is

$R_{k} = E {v_{k} v_{k}^{T}} = σ^{2} I_{2 \times 2}$

(21)

The initial state, which is given by Equation (1), can be computed using the first measurements,

{\tilde{ϑ}}_{x}

and

{\tilde{ϑ}}_{y}

x_{0}^{-} = \{\begin{matrix} cos {\tilde{ϑ}}_{x} \\ \frac{cos {\tilde{ϑ}}_{y} - cos {\tilde{ϑ}}_{x} sin ε}{cos ε} \end{matrix}\} = \{\begin{matrix} x_{0} \\ y_{0} \end{matrix}\}

(22)

while the state covariance matrix is

P_{k}^{-} = [\begin{matrix} P_{k} (1, 1) & P_{k} (1, 2) \\ P_{k} (1, 2) & P_{k} (2, 2) \end{matrix}]

(23)

where

P_{k} (1, 1)

P_{k} (1, 2)

and

P_{k} (2, 1)

, are provided by Equation (14), Equation (16) and Equation (15), respectively. Specifically, they can be written in terms of the state vector,

x_{k}^{-}

\begin{matrix} P_{k} (1, 1) & = & (e^{- 2 σ^{2}} - e^{- σ^{2}}) x_{k}^{2} + \frac{1}{2} (1 - e^{- 2 σ^{2}}) \\ P_{k} (1, 2) & = & [(e^{- σ^{2}} - e^{- 2 σ^{2}}) x_{k}^{2} - \frac{1}{2} (1 - e^{- 2 σ^{2}})] tan ε \\ P_{k} (2, 2) & = & (y_{k}^{2} + 2 x_{k} y_{k} tan ε + 2 x_{k}^{2} {tan}^{2} ε) (e^{- 2 σ^{2}} - e^{- σ^{2}}) + \\ + (\frac{2 - {cos}^{2} ε}{{cos}^{2} ε}) \frac{1}{2} (1 - e^{- 2 σ^{2}}) \end{matrix}

Predictor: Since the problem under analysis is static,

$\begin{matrix} x_{k}^{+} & = & x_{k - 1}^{-} \end{matrix}$

(24)

$\begin{matrix} P_{k}^{+} & = & P_{k - 1}^{-} \end{matrix}$

(25)
Corrector:

$\begin{matrix} K_{k} & = & P_{k}^{+} J_{h}^{T} (x_{k}) {(J_{h} (x_{k}) P_{k}^{+} J_{h}^{T} (x_{k}) + R_{k})}^{- 1} \end{matrix}$

(26)

$\begin{matrix} x_{k}^{-} & = & x_{k}^{+} + K_{k} (y_{k} - h (x_{k}^{+})) \end{matrix}$

(27)

$\begin{matrix} P_{k}^{-} & = & [I - K_{k} J_{h} (x_{k})] P_{k}^{+} \end{matrix}$

(28)

where the Jacobian, $J_{h} (x_{k})$ , of the observation is given by

$J_{h} (x_{k}) = \frac{\partial h (x_{k})}{\partial x_{k}} = [\begin{matrix} - \frac{1}{\sqrt{1 - x_{k}^{2} e^{σ^{2}}}} & 0 \\ - \frac{sin ε}{d_{k}} & - \frac{cos ε}{d_{k}} \end{matrix}] e^{σ^{2} / 2},$

(29)

where

$d_{k} = \sqrt{1 - e^{σ^{2}} {(x_{k} sin ε + y_{k} cos ε)}^{2}} .$

(30)

4.2. Numerical Tests

The true direction of the gravity in the inclinometer reference frame has been selected as

g_{true} = - \{\begin{matrix} sin α cos β \\ sin α sin β \\ cos α \end{matrix}\} where \{\begin{matrix} θ = 30 \deg \\ β = 135 \deg \end{matrix}

The inclinometer measurements are corrupted by a zero-mean Gaussian noise with a standard deviation of

σ = 0.1

deg., which is a common value for commercial inclinometers. In addition, the deviation from orthogonality was selected as

ε = 2

deg.

Figure 5 shows the results obtained using 200 samples. In particular, the top-left figure shows for the

g_{x}

component: (1) the noisy measurements (black markers), (2) the true value (solid black line), (3) the filtered estimates (solid blue line) and (4) the

\pm 3 σ

error bounds (dotted red lines). The top-right figure shows the same quantities for the

g_{x}

component.

The bottom left and right plots show the absolute errors for the

g_{x}

and

g_{y}

components, respectively. These error plots highlight, thanks to the analytically exact expressions of the covariance matrix terms, a very fast convergence of the EKF estimator.

5. Continuous Gravity Description

The irregular Earth shape (and mass distribution) has presented significant challenges to surveyers and cartographers throughout history. The early map makers considered the Earth to be a perfect sphere prior to Isaac Newton positing an oblate spheroid shape model. The axial–symmetric ellipsoid description of Earth, which has been formalized as the World Geodetic System (WGS-84) [27], is reliably used as the standard shape model for most of today’s GPS applications. The WGS-84 is defined by estimates of the Earth’s equatorial (

R_{E}

) and polar (

R_{P}

) radii, and a flattening factor (f)

\{\begin{matrix} R_{E} = 6, 378, 137 m \\ R_{P} = 6, 356, 752.3 m \end{matrix} and f = \frac{R_{E} - R_{P}}{R_{E}} .

Describing a point on the Earth’s surface can be done by Cartesian coordinates (

x, y, z

) or by spherical coordinates, using geocentric (

θ

) or geodetic (

ϕ

) latitude and longitude (

λ

). Figure 6 shows the two definitions of latitude for an ellipsoid-shaped model.

The relationship between these two latitude definitions is [28]:

tan θ = {(1 - f)}^{2} tan ϕ .

(31)

The greatest angular difference between geocentric and geodetic latitudes is, for the WGS-84 model, approximately 700 arcseconds, a difference that would introduce over 21 km of position error on Earth. However, even if the WGS-84 ellipsoidal model describes the gravity direction better than the sphere model, the actual gravity direction is dictated by the geoid model [29].

The geoid was first described by Gauss as the shape that the ocean’s surface would take if the only acceleration forces acting on it were Earth’s gravitational field and the rotational dynamics of the planet. In reality, this is still a simplification of the surface always orthogonal to the gravity direcion. The geoid is particularly useful for describing the effect that mass anomalies (such as mountain ranges or caverns) have on the local pointing direction of gravity, which can deviate from the reference ellipsoid by tens, hundreds or even thousands of arcseconds. Figure 7 shows the main parameters defining the relationship between the geoid surface and the WGS-84 ellipsoid. P is a point on the Earth’s surface that is located at altitude h normal to the ellipsoid and altitude H normal to the geoid. This point on the geoid surface is located at altitude N with respect to the WGS-84 ellipsoid. The angle between the true gravity direction (normal to the geoid and provided by the inclinometers) and the fictitious gravity direction (normal to WGS-84 reference model) is the geoid correction.

It is also worth pointing out that, unlike the WGS-84 axial–symmetric ellipsoid model, the geoid cannot be described by a simple equation because it is a function of the Earth’s mass distribution. However, the geometry of the geoid surface can be derived with respect to the WGS-84 model using the derivatives of the Earth’s gravitational potential. The Earth’s gravitational potential can be described in terms of spherical harmonics:

V (r, ϕ, λ) = \frac{G M}{r} \{1 + \sum_{n = 0}^{\infty} {(\frac{R_{E}}{r})}^{n} \sum_{m = 0}^{n} P_{n}^{m} (sin ϕ) [C_{n m} cos (m λ) + S_{n m} sin (m λ)]\},

(32)

where

G M = 3.986013 \times 10^{14}

^{3}

^{2}

is the Earth gravitational constant,

(r, ϕ, λ)

are the spherical coordinates at which the evaluation is made,

P_{n}^{m} (sin ϕ)

are the associated Legendre polynomials of degree n, order m and argument

sin ϕ

, and

C_{n m}

and

S_{n m}

are spherical harmonic coefficients.

Taking normalized partial derivatives of Equation (32) with respect to the geodetic latitude and longitude, we obtain a mathematical definition for the north–south

(δ N)

and east–west

(δ E)

deflections of the gravity direction with respect to the WGS-84 ellipsoidal model:

δ N = - \frac{1}{γ r} \cdot \frac{\partial V}{\partial ϕ} and δ E = - \frac{1}{γ r cos ϕ} \cdot \frac{\partial V}{\partial λ}

(33)

where

γ

is a scalar representing the theoretical normal gravity, whose value and clear explanation are provided in Reference [30].

5.1. Continuous Description of Geoid

The analytical complexities associated with the infinite series of Equations (32) and (33) make them unable to evaluate the gravity direction without introducing truncation error. From satellite measurement data, a geoid model can be constructed in order to generate

δ E

and

δ N

deflection reference databases for a specified grid of points on the Earth’s surface, defined by geodetic latitude

(ϕ_{k})

and longitude

(λ_{k})

. To use the geoid database, interpolation between discontinuous points is needed. The accuracy of this interpolation is obtained at the expense of carrying a larger database (fine grid geoid description). In this study, a simple continuous geoid model is derived in order to use geoid corrections in specific locations and/or in dynamic situations.

The continuous geoid model is provided by expressing the

δ E

and

δ N

deflections as linear combinations of a set of N orthogonal surfaces:

δ E (ϕ, λ) = \sum_{k = 1}^{N} α_{k} S_{k} (x (ϕ), y (λ)) and δ N (ϕ, λ) = \sum_{k = 1}^{N} β_{k} S_{k} (x (ϕ), y (λ)),

(34)

where the

α_{k}

and

β_{k}

are the unknown coefficients and the

S_{k} (x (ϕ), y (λ))

are a set of orthogonal surfaces satisfying

\int_{- 1}^{+ 1} \int_{- 1}^{+ 1} S_{i} (x, y) S_{i} (x, y) d x d y = 0, if i \neq j .

These orthogonal surfaces are built using orthogonal polynomials, such as Chebyshev polynomials,

T_{k} (x)

, which are orthogonal in the

x \in [- 1, + 1]

range. This means that, when dealing with a specific geographical range of interest,

ϕ \in [ϕ_{min}, ϕ_{max}]

and

λ \in [λ_{min}, λ_{max}]

, a mapping between the geographical coordinates

[ϕ, λ]

and the Chebyshev variables

[x, y]

must be defined. The most simple mapping is the linear

x = \frac{ϕ - ϕ_{min}}{ϕ_{max} - ϕ_{min}} and y = \frac{λ - λ_{min}}{λ_{max} - λ_{min}} .

Using Chebyshev orthogonal polynomials, the orthogonal surfaces can be defined as

S_{k} (x, y) = T_{i} (x (ϕ)) T_{j} (y (λ)) where i, j = 0, 1, 2, \dots

where any two orthogonal surfaces cannot be associated with the same

[i, j]

pair of Chebyshev polynomials.

Chebyshev orthogonal polynomials of the first kind can conveniently be computed in a recursive way. Starting with

T_{0} (x) = 1

and

T_{1} (x) = x

, the subsequent terms are computed by the recursive relationship,

T_{k} (x) = 2 x T_{k - 1} (x) - T_{k - 2} (x), k = 2, 3, \dots .

(35)

The problem of fitting the geoid data grids,

δ E (ϕ_{i}, λ_{j})

and

δ N (ϕ_{i}, λ_{j})

, consists of two over-determined linear systems:

δ E (ϕ_{i}, λ_{j}) = \sum_{k = 1}^{N} α_{k} S_{k} (x (ϕ_{i}), y (λ_{j})) and δ N (ϕ_{i}, λ_{j}) = \sum_{k = 1}^{N} β_{k} S_{k} (x (ϕ_{i}), y (λ_{j})) .

where the total number of the geoid database points is

N_{p} > N

, for the problem to admit a solution. The previous equation describes two linear systems that can be written in matrix form:

S_{E} α = ξ and S_{N} β = η,

where the two coefficient vectors,

α

and

β

, can be computed by least squares:

α = {(S_{E}^{T} S_{E})}^{- 1} S_{E}^{T} ξ and β = {(S_{N}^{T} S_{N})}^{- 1} S_{N}^{T} η .

This continuous description of the geoid is applied in the next subsection to a geographical region where the geoid gravity direction deviations are particularly significant.

5.2. Example of Continuous Geoid Description: Himalaya Region

The results of a numerical example of this procedure are shown in Figure 8 providing north–south and east–west deviations, in arcsec, of the gravity direction with respect to the direction normal to the WGS-84 ellipsoid surface. The region of interest selected is in the coordinates range

[ϕ_{min}, ϕ_{max}] = [28 . 083^{\circ}, 28 . 875^{\circ}]

and

[λ_{min}, λ_{max}] = [83 . 667^{\circ}, 84 . 458^{\circ}]

The associated geoid database consists of a

10 \times 10

grid (

N_{p} = 100

points), spanning a region of 2000 km

\times 2000

km. The number of orthogonal surfaces used for the continuous geoid model was

N = 55

. This is also the number of unknowns (size of

α

and

β

vectors). The data points of the geoid grid are shown with black markers while the continuous models are indicated by the surfaces. The

L_{2}

norms of the least-squares residuals are shown in Figure 9. The condition number of the matrix to invert (in the least-squares process), which is provided in Figure 10 as a function of the number of orthogonal surfaces, does not exceed 200.

By selecting orthogonal polynomials (instead of monomials), the condition number of the least-squares matrix to invert is reduced and, consequently, the numerical process becomes more robust and provides a more reliable numerical estimation of the unknown coefficient vectors,

α

and

β

6. Discussion

This study introduced the mathematical tools necessary for extracting a highly accurate estimation of the local gravity (plumb-line) direction from a pair of non-orthogonal inclinometers. The statistical analyses performed herein can be used to quantify the performance of static tilt sensors, since the measured direction of gravity is known to be normal to the reference geoid surface. The next section provides an example use case of these tools for high-precision continuous gravity direction estimates under static operating conditions.

The static Stellar Positioning System [15,16,17] introduces a framework for estimating the geographical position of a vehicle on the surface of a planetary body using a pair of inclinometers along with an atomic clock and a night-sky camera to observe celestial directions (e.g., stars, visible planets). This system, which can be used in GPS-denied scenarios or in locations where the GPS signal is not available (i.e., on Moon or Mars), assumes that a reliably accurate estimate of the local gravity direction is available, which in turn would require a geoid correction. The resulting mathematical estimation problem is nonlinear and, therefore, can be solved by iterative nonlinear least squares. The mathematical procedure requires computing the Jacobian of the gravity direction and, consequently, the first derivatives of Equation (35). The first derivatives can be obtained using the recursive equation,

\frac{d T_{k} (x)}{d x} = 2 (T_{k - 1} (x) + x \frac{d T_{k - 1} (x)}{d x}) - \frac{d T_{k - 2} (x)}{d x} .

More accurate models will require the computation of the Hessian or even higher-order terms. The subsequent derivatives needed in these cases can be computed recursively:

\frac{d^{n} T_{k}}{d x^{n}} = 2 (n \frac{d^{n - 1} T_{k - 1}}{d x^{n - 1}} + x \frac{d^{n} T_{k - 1}}{d x^{n}}) - \frac{d^{n} T_{k - 2}}{d x^{n}}, n = 2, 3, \dots .

Future Research Directions

The authors are planning to validate the proposed gravity direction estimation method using two non-orthogonal inclinometers by:

Including the mounting angle deviation from orthogonality, $ε$ , as an element of the EKF state vector. This would keep the system calibrated from variations of $ε$ caused by vibrations and/or thermal expansions.
Integrating the proposed system with an Inertial Measurement Unit (IMU) to extend the use of the static “Stellar Positioning System” to dynamical scenarios. Similar to the previous point, this will extend the EKF state vector to account for the variations in the transformation matrix between the inclinometer and IMU reference frames.
Including in the EKF state estimator, for the static use “Stellar Positioning System,” the transformation matrix between the inclinometer and the night camera reference frames.

As for the continuous geoid model, the authors are planning to complete the analysis on the geoid north–south and east–west deflection estimates by:

Investigating the convergence radius as a function of the geographical coordinates.
Investigating the accuracy of the convergence as a function of the geoid database grid resolution.
Including, in the current nonlinear position estimation algorithm, a divergence identification algorithm that re-initializes the estimation process using a new random guess selected within an admissible distance bound.

7. Conclusions

This study introduces a new approach to accurately estimate the gravity direction using a pair of non-orthogonal inclinometers affected by zero-mean Gaussian errors. The proposed approach estimates the gravity direction and its

3 \times 3

covariance matrix, with the main three terms derived analytically. The accuracy of these statistical parameters is numerically validated by Monte Carlo simulation. The proposed method is then validated using an Extended Kalman Filter, which highlights a fast convergence to a very accurate solution.

Since the gravity direction is normal to the geoid surface, a continuous geoid model for the north–south and east–west corrections is provided as a linear combination of a set of orthogonal surfaces that are derived using Chebyshev polynomials of the first kind. An example application is provided for a specific region (centered in the Himalayas) to show the continuous description of the gravity direction deviations with respect to the WGS-84 ellipsoid model.

This covariance analysis, along with the continuous description of the geoid, can be applied to any static or dynamic scenarios involving the measurement of gravity using a two-axis inclinometer.

Author Contributions

Conceptualization, D.M.; methodology, D.M.; software, D.M. and A.G.; validation, A.G.; formal analysis, D.M.; writing—original draft preparation, D.M.; writing—review and editing, D.M. and A.G.; supervision, D.M. All authors have read and agreed to the published version of the manuscript.

Funding

This research was in part supported by a NASA grant#: 80NSSC19P1369.

Acknowledgments

The author acknowledges the initial technical investigation by Hunter Johnston and the support in coding provided by Carl Leake.

Conflicts of Interest

The author declares no conflict of interest.

References

Łuczak, S.; Ekwińska, M. Electric-Contact Tilt Sensors: A Review. Sensors 2021, 21, 1097. [Google Scholar] [CrossRef] [PubMed]
Sun, H.; Deng, X. Accuracy analysis of inclinometer in slope monitoring. Hydropower Autom. Dam. Monit. 2015, 1. Available online: https://en.cnki.com.cn/Article_en/CJFDTotal-DBGC201501008.htm (accessed on 23 August 2021).
Stark, T.D.; Choi, H. Slope inclinometers for landslides. Landslides 2008, 5, 339–350. [Google Scholar] [CrossRef]
Ruzza, G.; Guerriero, L.; Revellino, P.; Guadagno, F.M. A Multi-Module Fixed Inclinometer for Continuous Monitoring of Landslides: Design, Development, and Laboratory Testing. Sensors 2020, 20, 3318. [Google Scholar] [CrossRef] [PubMed]
Di Maio, C.; Vassallo, R.; Vallario, M.; Pascale, S.; Sdao, F. Structure and kinematics of a landslide in a complex clayey formation of the Italian Southern Apennines. Eng. Geol. 2010, 116, 311–322. [Google Scholar] [CrossRef]
Ziemer, R.R. Measurement of Soil Creep by Inclinometer; Engineering Technology Report ETR-7100-4; US Department of Agriculture, Forest Service: Washington, DC, USA, 1977; 10p.
Duan, C.; Mao, L.; Wang, C. Application of Sensor Technology to Detection of Zero-spot and Zero-line of Artilleries. J. Test Meas. Technol. 2012, 5. [Google Scholar]
Mao, L.; Ling, L. Compensation Method of Cannon Tilted Disturbance Torque. Ordnance Ind. Autom. 2011, 3. [Google Scholar]
Clarkson, L.; Williams, D.; Seppälä, J. Real-time monitoring of tailings dams. Georisk Assess. Manag. Risk Eng. Syst. Geohazards 2021, 15, 113–127. [Google Scholar] [CrossRef]
Scott, M.; Lo, R.; Thavaraj, T. Use of instrumentation to safeguard stability of a tailings Dam. In Proceedings of the 7th FMGM 2007: Field Measurements in Geomechanics, Boston, MA, USA, 24–27 September 2007; pp. 1–13. [Google Scholar]
Besnard, S.; Khalil, W. Calibration of parallel robots using two inclinometers. In Proceedings of the 1999 IEEE International Conference on Robotics and Automation (Cat. No. 99CH36288C), Detroit, MI, USA, 10–15 May 1999; Volume 3, pp. 1758–1763. [Google Scholar]
Hou, X.; Yang, X.; Huang, Q. Using inclinometers to measure bridge deflection. J. Bridge Eng. 2005, 10, 564–569. [Google Scholar] [CrossRef]
Olaszek, P.; Wyczałek, I.; Sala, D.; Kokot, M.; Świercz, A. Monitoring of the Static and Dynamic Displacements of Railway Bridges with the Use of Inertial Sensors. Sensors 2020, 20, 2767. [Google Scholar] [CrossRef] [PubMed]
D’Alfonso, L.; Garone, E.; Muraca, P.; Pugliese, P. On the use of the inclinometers in the PnP Problem. In Proceedings of the 2013 European Control Conference (ECC), Zurich, Switzerland, 17–19 July 2013; pp. 4112–4117. [Google Scholar]
Parish, J.J.; Parish, A.S.; Swanzy, M.; Woodbury, D.P.; Mortari, D.; Junkins, J.L. The Stellar Positioning System (Part I): An Autonomous Position Determination Solution. Navigation 2010, 57, 1–12. [Google Scholar] [CrossRef]
Woodbury, D.P.; Parish, J.J.; Parish, A.S.; Swanzy, M.; Denton, R.; Mortari, D.; Junkins, J.L. The Stellar Positioning System (Part II): Improving Accuracy During Implementation. Navigation 2010, 57, 13–24. [Google Scholar] [CrossRef]
Mortari, D. Lander State Determination via Stellar Positioning System; Technical Report Project: 80NSSC19P1369; Texas A&M University: College Station, TX, USA, 2020. [Google Scholar]
Chang, H.Y.; Chang, Y.C.; Liu, W.F. A Highly Sensitive Two-Dimensional Inclinometer Based on Two Etched Chirped-Fiber-Grating Arrays. Sensors 2017, 17, 2922. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Xu, H.; Zhao, Y.; Zhang, K.; Jiang, K. A Capacitive MEMS Inclinometer Sensor with Wide Dynamic Range and Improved Sensitivity. Sensors 2020, 20, 3711. [Google Scholar] [CrossRef] [PubMed]
Kong, X.; Chen, Q.; Wang, J.; Gu, G.; Wang, P.; Qian, W.; Ren, K.; Miao, X. Inclinometer Assembly Error Calibration and Horizontal Image Correction in Photoelectric Measurement Systems. Sensors 2018, 18, 248. [Google Scholar] [CrossRef] [PubMed] [Green Version]
O’Leary, P.; Harker, M. Inverse boundary value problems with uncertain boundary values and their relevance to inclinometer measurements. In Proceedings of the 2014 IEEE International Instrumentation and Measurement Technology Conference (I2MTC) Proceedings, Montevideo, Uruguay, 12–15 May 2014; pp. 165–169. [Google Scholar] [CrossRef]
Li, C.; Azzam, R.; Fernández-Steeger, T.M. Kalman filters in geotechnical monitoring of ground subsidence using data from MEMS sensors. Sensors 2016, 16, 1109. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Mortari, D.; Majji, M. Multiplicative Measurement Model. J. Astronaut. Sci. 2009, 57, 47–60. [Google Scholar] [CrossRef]
Stark, H.; Woods, J.W. Probability and Random Processes with Applications to Signal Processing; Prentice-Hall: Upper Saddle-River, NJ, USA, 2002. [Google Scholar]
Schleppe, J.B. Development of a Real-Time Attitude System Using a Quaternion Parameterization and Non-Dedicated GPS Receivers. Master’s Thesis, Department of Geomatics Engineering, University of Calgary, Calgary, AB, Canada, July 1996; pp. 69–71. [Google Scholar]
Vanícek, P.; Krakiwsky, E. Geodesy: The Concepts; Elsevier Science: Amsterdam, The Netherlands, 1986; p. 197. [Google Scholar]
International Civil Aviation Organization. World Geodetic System—1984 (WGS-84) Manual. Available online: https://gis.icao.int/eganp/webpdf/REF08-Doc9674.pdf (accessed on 23 August 2021).
Snyder, J.P. Map Projections—A Working Manual; U.S. Government Publishing Office: Washington, DC, USA, 1987.
Lerch, F.; Wagner, C.; Smith, D.; Sandson, M.; Brownd, J.; Richardson, J. Gravitational Field Models for the Earth (GEM 1&2); Technical Report X-553-72-146; Goddard Space Flight Center: Greenbelt, MD, USA, 1972. [Google Scholar]
NIMA/Geodesy and Geophysics Department. Department of Defense World Geodetic System 1984: Its Definition and Relationships with Local Geodetic Systems; Technical Report 8350.2; National Imagery and Mapping Agency: Springfield, VA, USA, 2000. [Google Scholar]

Figure 1. Geometry of the inclinometer system and associated reference frame.

Figure 2. Geometrical explanation justifying the result provided by Equation (6).

Figure 3. Monte Carlo tests: histograms of the angles and the gravity components.

Figure 4. Monte Carlo validation for the

P (1, 3)

P (2, 3)

and

P (3, 3)

terms.

Figure 4. Monte Carlo validation for the

P (1, 3)

P (2, 3)

and

P (3, 3)

terms.

Figure 5. Extended Kalman Filter test results.

Figure 6. Geocentric (

θ

) and geodetic (

ϕ

) latitude angles.

Figure 6. Geocentric (

θ

) and geodetic (

ϕ

) latitude angles.

Figure 7. Gravity description with respect to geoid and reference ellipsoid surfaces.

Figure 8. East–West and North–South geoid deviation (Himalaya region).

Figure 9. Residual

L_{2}

norms.

Figure 9. Residual

L_{2}

norms.

Figure 10. Least-squares condition number.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Mortari, D.; Gardner, A. High Accurate Mathematical Tools to Estimate the Gravity Direction Using Two Non-Orthogonal Inclinometers. Sensors 2021, 21, 5727. https://doi.org/10.3390/s21175727

AMA Style

Mortari D, Gardner A. High Accurate Mathematical Tools to Estimate the Gravity Direction Using Two Non-Orthogonal Inclinometers. Sensors. 2021; 21(17):5727. https://doi.org/10.3390/s21175727

Chicago/Turabian Style

Mortari, Daniele, and Anthony Gardner. 2021. "High Accurate Mathematical Tools to Estimate the Gravity Direction Using Two Non-Orthogonal Inclinometers" Sensors 21, no. 17: 5727. https://doi.org/10.3390/s21175727

APA Style

Mortari, D., & Gardner, A. (2021). High Accurate Mathematical Tools to Estimate the Gravity Direction Using Two Non-Orthogonal Inclinometers. Sensors, 21(17), 5727. https://doi.org/10.3390/s21175727

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

High Accurate Mathematical Tools to Estimate the Gravity Direction Using Two Non-Orthogonal Inclinometers

Abstract

1. Introduction

System Definition

2. Gravity Direction Estimation

3. Gravity Direction Covariance Matrix

Computation of the P(1,3), P(2,3) and P(3,3) Terms via Covariance Law

4. Covariance Analysis Numerical Validation

4.1. Extended Kalman Filter (EKF) Validation

4.2. Numerical Tests

5. Continuous Gravity Description

5.1. Continuous Description of Geoid

5.2. Example of Continuous Geoid Description: Himalaya Region

6. Discussion

Future Research Directions

7. Conclusions

Author Contributions

Funding

Acknowledgments

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI