# **International Association of Geodesy Symposia**

# Jeffrey T. Freymueller Laura Sánchez  *Editors*

# 5th Symposium on Terrestrial Gravimetry: Static and Mobile Measurements (TG-SMM 2019)

Proceedings of the Symposium in Saint Petersburg, Russia, October 1 – 4, 2019

# International Association of Geodesy Symposia

Jeffrey T. Freymueller, Series Editor Laura Sánchez, Series Assistant Editor

#### **Series Editor**

Jeffrey T. Freymueller Endowed Chair for Geology of the Solid Earth Department of Earth and Environmental Sciences Michigan State University East Lansing, MI, USA

#### **Assistant Editor**

Laura Sánchez Deutsches Geodätisches Forschungsinstitut Technische Universität München Munich, Germany

# International Association of Geodesy Symposia

Jeffrey T. Freymueller, Series Editor Laura Sánchez, Series Assistant Editor

Symposium 112: Geodesy and Physics of the Earth: Geodetic Contributions to Geodynamics Symposium 113: Gravity and Geoid Symposium 114: Geodetic Theory Today Symposium 115: GPS Trends in Precise Terrestrial, Airborne, and Spaceborne Applications Symposium 116: Global Gravity Field and Its Temporal Variations Symposium 117: Gravity, Geoid and Marine Geodesy Symposium 118: Advances in Positioning and Reference Frames Symposium 119: Geodesy on the Move Symposium 120: Towards an Integrated Global Geodetic Observation System (IGGOS) Symposium 121: Geodesy Beyond 2000: The Challenges of the First Decade Symposium 122: IV Hotine-Marussi Symposium on Mathematical Geodesy Symposium 123: Gravity, Geoid and Geodynamics 2000 Symposium 124: Vertical Reference Systems Symposium 125: Vistas for Geodesy in the New Millennium Symposium 126: Satellite Altimetry for Geodesy, Geophysics and Oceanography Symposium 127: V Hotine-Marussi Symposium on Mathematical Geodesy Symposium 128: A Window on the Future of Geodesy Symposium 129: Gravity, Geoid and Space Missions Symposium 130: Dynamic Planet - Monitoring and Understanding . . . Symposium 131: Geodetic Deformation Monitoring: From Geophysical to Engineering Roles Symposium 132: VI Hotine-Marussi Symposium on Theoretical and Computational Geodesy Symposium 133: Observing our Changing Earth Symposium 134: Geodetic Reference Frames Symposium 135: Gravity, Geoid and Earth Observation Symposium 136: Geodesy for Planet Earth Symposium 137: VII Hotine-Marussi Symposium on Mathematical Geodesy Symposium 138: Reference Frames for Applications in Geosciences Symposium 139: Earth on the Edge: Science for a Sustainable Planet Symposium 140: The 1st International Workshop on the Quality of Geodetic Observation and Monitoring Systems (QuGOMS'11) Symposium 141: Gravity, Geoid and Height Systems (GGHS2012) Symposium 142: VIII Hotine-Marussi Symposium on Mathematical Geodesy Symposium 143: Scientific Assembly of the International Association of Geodesy, 150 Years Symposium 144: 3rd International Gravity Field Service (IGFS) Symposium 145: International Symposium on Geodesy for Earthquake and Natural Hazards (GENAH) Symposium 146: Reference Frames for Applications in Geosciences (REFAG2014) Symposium 147: Earth and Environmental Sciences for Future Generations Symposium 148: Gravity, Geoid and Height Systems 2016 (GGHS2016) Symposium 149: Advancing Geodesy in a Changing World Symposium 150: Fiducial Reference Measurements for Altimetry Symposium 151: IX Hotine-Marussi Symposium on Mathematical Geodesy Symposium 152: Beyond 100: The Next Century in Geodesy Symposium 153: Terrestrial Gravimetry: Static and Mobile Measurements (TG-SMM 2019)

# 5th Symposium on Terrestrial Gravimetry: Static and Mobile Measurements (TG-SMM 2019)

Proceedings of the Symposium in Saint Petersburg, Russia, October 1 – 4, 2019

Edited by

Jeffrey T. Freymueller, Laura Sánchez

Jeffrey T. Freymueller Endowed Chair for Geology of the Solid Earth Department of Earth and Environmental Sciences Michigan State University East Lansing, MI USA

#### *Assistant Editor*

Laura Sánchez Deutsches Geodätisches Forschungsinstitut Technische Universität München Munich Germany

*Series Editor Associate Editors*

Oleg Stepanov State Research Center of the Russian Federation Concern CSRI Elektropribor, JSC St. Petersburg Russia

Anton Krasnov State Research Center of the Russian Federation Concern CSRI Elektropribor, JSC St. Petersburg Russia

Olga Yashnikova State Research Center of the Russian Federation Concern CSRI Elektropribor, JSC St. Petersburg Russia

Dmitry Taranovskiy State Research Center of the Russian Federation Concern CSRI Elektropribor, JSC St. Petersburg Russia

ISSN 0939-9585 ISSN 2197-9359 (electronic) International Association of Geodesy Symposia ISBN 978-3-031-25901-2 ISBN 978-3-031-25902-9 (eBook) https://doi.org/10.1007/978-3-031-25902-9

© The Editor(s) (if applicable) and The Author(s) 2023. This book is an open access publication.

**Open Access** This book is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this book are included in the book's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the book's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

The use of general descriptive names, registered names, trademarks, service marks, etc. in this publication does not imply, even in the absence of a specific statement, that such names are exempt from the relevant protective laws and regulations and therefore free for general use.

The publisher, the authors, and the editors are safe to assume that the advice and information in this book are believed to be true and accurate at the date of publication. Neither the publisher nor the authors or the editors give a warranty, expressed or implied, with respect to the material contained herein or for any errors or omissions that may have been made. The publisher remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

This Springer imprint is published by the registered company Springer Nature Switzerland AG The registered company address is: Gewerbestrasse 11, 6330 Cham, Switzerland

# **Preface**

For the fifth time, the International Association of Geodesy (IAG) held the Symposium on Terrestrial Gravimetry: Static and Mobile Measurements (TG-SMM2019) in St. Petersburg, Russia, October 1–4, 2019. The Symposium was hosted by the *State Research Center of the Russian Federation Concern CSRI Elektropribor, JSC* and was attended by 75 participants from 15 different countries. 32 oral and 20 poster contributions were presented in four sessions:


This volume contains 20 selected papers from the all sessions of the symposium. All published papers were peer-reviewed, and we warmly recognize the contributions and support of the Volume Editors and Reviewers (see the list in later pages).

# **Contents**

#### **Part I Terrestrial, Shipboard and Airborne Gravimetry**


#### **Part II Absolute Gravimetry**



**Part I**

**Terrestrial, Shipboard and Airborne Gravimetry**

# **Measurement of Absolute Gravity and Deflections of the Vertical at Sea**

A. V. Sokolov, A. A. Krasnov, N. V. Kuz'mina, and Yu. F. Stus'

#### **Abstract**

Methods to measure absolute gravity and deflections of the vertical on a moving base are presented. The breadboard of integrated gravimetric system is described. The first results of experimental studies confirmed the possibility of high-precision measurement of the absolute values of gravity and deflection of the vertical at sea.

#### **Keywords**

Absolute gravity - Deflections of the vertical - Gravity measurements - Integrated gravimetric system

#### **1 Introduction**

Knowledge of absolute values of gravity and deflection of the vertical (DOV) is essential for solving a number of problems in geodesy, high-precision inertial navigation and fundamental position, navigation and time support. In water areas, the absolute values of gravity are calculated as a sum of the gravity a priori value measured at an onshore reference station, and of the observed gravity measured with relative gravimeters from marine and air vehicles. The errors in determining the drift and scale of relative gravimeters reduce the accuracy of determining the absolute values of gravity, and the necessity to regularly reference the results of offshore measurements to the onshore absolute value imposes serious operational limitations on the gravity survey procedure. In view of the fact that the geological exploration tasks require only the knowledge of character of gravity anomalies in the survey area, most of modern geophysical works are carried out without precise referencing of measurements to

A. V. Sokolov · A. A. Krasnov (-) · N. V. Kuz'mina Concern CSRI Elektropribor, JSC, ITMO University, St. Petersburg, Russia e-mail: anton-krasnov@mail.ru

Y. F. Stus'

the reference stations, which considerably increases the error of absolute values of gravity measured in water areas.

There is a known technique for the relative gravimeters' measurements re-calculation to the absolute level, using the global geopotential models (Zheleznyak et al. 2015). However, limited spatial resolution of the global models, as well as their significant errors in the areas with high gravity anomalies can make it impossible to precisely determine the absolute values of gravity using the above technique. At present, there are no commercial devices for measuring the absolute value of gravity from moving vehicles. A number of companies are currently studying and developing such devices (Bidel et al. 2018; Baumann et al. 2012).

Integrating the gravity anomalies according to Vening-Meinesz formulae has been the main method of highprecision determining of DOV in water areas for as long as 100 years (Vening-Meinesz 1928). However, this method is extremely labor-consuming, since it requires accumulation of background gravimetric data for an area that is much larger than the point specified for determining the DOV. Moreover, the error of the absolute value of gravity discussed above is a methodological error of the DOV calculating according to Vening-Meinesz formulae. Also it should be noted that the existing methods of space geodesy do not provide the DOV precisely enough, especially in regions with large gravity anomalies (Koneshov et al. 2013).

A design concept of an integrated gravimetric system for measuring the absolute gravity on a moving base, and

Institute of Automation and Electrometry, Siberian Branch of the Russian Academy of Sciences, Novosibirsk, Russia

J. T. Freymueller, L. Sánchez (eds.), *5th Symposium on Terrestrial Gravimetry: Static and Mobile Measurements (TG-SMM 2019)*, International Association of Geodesy Symposia 153, https://doi.org/10.1007/1345\_2022\_140

the first results of its bench tests were discussed at the Symposium TGSMM-2016 (Peshekhonov et al. 2016). Such system is useful to improve accuracy of gravity surveys that were carried out without reference measurements. Further development of the concept resulted in the integrated system construction and incorporation in the equipment for the DOV direct measurement by astrogeodetic method at sea which is essential in inertial navigation (Chatfield 1997).

A breadboard of gyrostabilized zenith telescope was made, and its operation methodology was developed. This paper presents the principles of the integrated system construction, including an absolute gravimeter, a zenith telescope, a system for their gyroscopic stabilization, and receiving equipment of global navigation satellite systems (GNSS). The integrated system is intended for measuring the absolute values of gravity with 1 mGal RMSE and DOV with 1 arcsec RMSE at sea. The system operation methodology is described, and the results of field testing of the system breadboard are discussed.

#### **2 Principles of Measuring the Absolute Values of Gravity and DOV at Sea**

Land gravity instruments for high-precision determination of the absolute values of gravity are commercially available and are based on measuring the time and length intervals of a test body fall in vacuum (Vitushkin 2015). Offshore applications of such devices are limited by the effect of inertial accelerations, caused by pitch/roll and orbital motion, on the measuring system of the gravimeter. For this reason, in case of moving vehicles, measurements are taken using relative gravimeters with a gyroscopic system for the sensitive element stabilization in the horizon plane. At the same time, the influence of the vertical component of inertial acceleration is compensated by low-frequency filtering methods involving the external navigation data (Stepanov and Koshaev 2010).

Obviously, absolute measurement of gravity at sea requires the sensitive axis of the absolute gravimeter to be stabilized in the direction of the local vertical with an accuracy of about 15 arcsec, in order to remove the effect of horizontal inertial accelerations. However, due to the effect of vertical inertial acceleration on the absolute gravimeter, the vast majority of its measurements are unreliable, and the methods of frequency filtering cannot be used in absence of a model of the test body motion in the field of inertial accelerations. The solution to this problem was found by integration of the initial data of the absolute and relative gravimeters. The idea of data integration is as follows. Based on the relative gravimeter data, current inertial vertical accelerations are determined, and the absolute gravimeter measurements at which the accelerations were minimal are selected. The resulting measurements of the absolute value of gravity are used for compensating for the errors in the relative gravimeter, caused by drift and nonlinearity of the scale, as well as the measurements referencing to the absolute level (Peshekhonov et al. 2016).

Field digital zenith cameras have been developed and widely used for measuring another important geodetic parameter DOV on the ground (Hirt and Bürki 2002; Tian et al. 2014; Halicioglu et al. 2012; Gerstbach and Pichler 2003). They help to implement the astrogeodetic method of DOV absolute measurement based on determining the direction to the stars located at the zenith, with known equatorial coordinates (right ascension ' and declination •). At that, the equivalence of astronomical coordinates (latitude ®, longitude œ) of the observation point and the equatorial coordinates of the stars is used. This equivalence is due to the validity of the following relations (Torge 2001):

$$\begin{array}{l} \varphi = \delta; \\ \lambda = \alpha - \theta, \end{array}$$

where ™ is Greenwich apparent sidereal time.The Helmert, or astronomic, DOV are calculated in accordance with the basic expressions (Jekeli 1999):

$$\begin{array}{l} \xi = \phi - B; \\ \eta = (\lambda - L)\cos\phi, \end{array}$$

where Ÿ is DOV projection on the meridian plane; ˜ is DOV projection on the prime vertical plane; *B, L* are geodetic coordinates of the observation point (latitude and longitude).As with the absolute gravimeter measuring absolute values of gravity at sea, the zenith telescope obviously needs to be stabilized for measuring DOV on a moving vehicle. However, to ensure the DOV measurement accuracy at the level of few seconds of arc, the required precision of stabilization should be comparable. Even on a slightly oscillating base, it is almost impossible to achieve such precision of stabilization by gyroscopic means only. Therefore, the only solution is to stabilize the telescope's sight axis along the normal to the reference ellipsoid rather than in the direction of the local vertical. The above normal is drawn according to the data from two sources: the zenith telescope itself, and the GNSS receiver (Peshekhonov et al. 1995).

A schematic of a gyrostabilized zenith telescope is shown in Fig. 1. Stabilization motors are controlled by the signals of mismatch between the geodetic coordinates and the astronomical coordinates of the crossing point of the zenith telescope's sight axis with the celestial sphere.

A block diagram of the algorithm for determining the DOV components on a moving base is presented in Fig. 2. The objective of near-zenith stars observation is to record a sequence of frames containing images of the stars, using a TV camera, with the frame record time being simultaneously **Fig. 1** Schematic of a gyrostabilized zenith telescope

fixed. In each frame, the coordinates of energy centers of the stellar images are determined and then used, along with the star catalog data, for stars identification (Mantsvetov et al. 2006).

As a result of identification, a data set is formed, where the coordinates of stars on the images are matched to the equatorial coordinates of stars from the catalogue. Based on this set, the parameters of transformation between the frame of the TV camera photodetector and the standard frame are calculated; these parameters are used for transforming the coordinates of the photodetector central point into equatorial coordinates and then, taking into account Greenwich apparent sidereal time, into astronomical coordinates. Moreover, the transformation parameters are useful in determining the angle of the photodetector frame rotation relative to the standard frame (azimuth of the TV camera row), which is used for feedback in the azimuth stabilization loop. After that, deviation of the astronomical coordinates of the crossing point of the zenith telescope's axis of sight with the celestial sphere relative to the geodetic vertical is determined. The deviation values are taken into account in the readings of accelerometers, and also used in the stabilization loops. Then linear accelerations are subtracted from the accelerometer readings using GNSS receiver data.

In order to compensate for the tilt of the sight axis relative to the rotation axis of the optronic device, and to prevent the effect of accelerometer bias, observations are carried out in two diametrically opposed positions, with the zenith telescope being turned to 180 degrees. The final values of the DOV components are determined by averaging the accelerometer data obtained in two positions of the zenith telescope.

#### **3 Integrated System Operation Methodology**

Combined operation of an absolute gravimeter and a relative one does not suggest any serious changes in the methodology of marine gravity surveys. The gravimeters are placed on a common base onboard a vessel, as close to the vessel's metacenter as possible (Fig. 3). The navigation data, including geographic coordinates, time stamps, and the data from the absolute and relative gravimeters are received to a common laptop which is also used for further post-processing of the data.

The survey starts from taking reference measurements at port, to determine the absolute value of gravity according to the methodology described above, as well as the initial value of the drift of the relative gravimeter. The gravity increments in the survey profiles are measured with the relative gravimeter. If the sea is weak (sea state 0,1,2), it is possible to determine the absolute value of gravity; to do so, the vessel should stay at the specified point of the water area for one to two hours. The survey results are processed in offline mode. Some additional offshore reference points would provide more detailed estimate of the relative gravimeter's drift and reduce the effect of the scale factor error during operations at a significant distance from the reference gravity station.

The zenith telescope is installed on the gyrostabilizer instead of an absolute ballistic gravimeter. In view of the fact that DOV measurements are taken at night and only in case weak sea and clear sky, this circumstance does not impose any significant limitations on the methodology of the integrated system operation during gravity measurements.

**Fig. 2** Block diagram of the algorithm for determining the DOV components on a moving base

#### **4 Integrated System Breadboard**

To verify the proposed methods for measuring the absolute values of the gravitational field parameters, a breadboard of the integrated system was assembled and tested. The breadboard composition was based on commercially available equipment.

For continuous measurement of gravity increments, the breadboard included a mobile gravimeter Chekan-AM (Shelf-E model) designed by Concern CSRI Elektropribor, JSC (Krasnov et al. 2014; Evstifeev et al. 2014). This gravimeter is based on a gravity sensor with double quartz elastic system, installed in a two-axis gyrostabilizer. Rootmean-square (RMS) error of gravity increment measurement with the gravimeter Chekan-AM under the action of dynamic disturbances does not exceed 0.4 mGal (Sokolov et al. 2016).

In the breadboard, the absolute value of gravity is measured with a ballistic gravimeter GABL-PM designed for field operation by the Institute of Automation and Electrometry of the Siberian Branch of the RAS (Bunin et al. 2010). RMS error of measurement of the gravity with the gravimeter does not exceed 5 Gal. The dimensions of the gravimeter GABL-PM is 42 47 93 cm and its weight is 59 kg, which made it possible to install it in the gyrostabilizer commercially manufactured by Concern CSRI Elektropribor,

**Fig. 3** A schematic of integrated system

JSC. The operating principle and the sensitive elements of this gyrostabilizer are similar to the stabilization system of the gravimeter Chekan-AM, while its weight and size are several times greater, so it could be used for the absolute gravimeter stabilization.

The main elements of the zenith telescope are a catadioptric lens with a focal distance of 2 m, and a specially designed TV camera -µ-62 with intrinsic function of automated determination of the coordinates of energy centers on the images of point objects. Such a system of a lens and a camera can determine the position of point objects on the TV camera plane with an accuracy of 1/20 pixel, taking into account the influence of various external factors, which corresponds to 0.1 arcsec. To determine DOV, electronic inclination sensors Zerotronic (accelerometers) with an error of 0.2 arcsec are installed on the zenith telescope base. The error in determination DOV components for such system is less than 0.500. The zenith telescope has an intrinsic azimuth rotation drive, and due to its configuration it can be installed in the gyrostabilizer instead of the absolute gravimeter.

To determine the geographic coordinates of location and to receive the time stamps, the breadboard included a GNSS receiver Javad of geodetic accuracy grade.

#### **5 Results of the Bench Tests**

Bench tests of the gravimetric system breadboard were carried out at the Concern CSRI Elektropribor. The test program included tests on the breadboard system accuracy on a fixed base as well as on dynamic benches of vertical displacements and a sea wave simulator.

The tests of the breadboard instrumental accuracy were carried out on a fixed base during 3 days. The methodology of continuous measurement of the absolute gravity consisted in a series of initial measurements by GABL-PM, continuous measurements of the gravity increments by Chekan-AM, and a series of final measurements by GABL-PM, which were used to refine the drift of the Chekan-AM gravimeter. The measurement results are shown in Fig. 4.

The RMS error of the initial and final measurements taken with the GABL-PM gravimeter was ¢gabs <sup>D</sup> 0.02 mGal. After taking into account the drift of the Chekan-AM gravimeter using the GABL-PM gravimeter readings and after the introduction of correction for the lunar-solar tide effects, calculated theoretically, the RMS error of the Chekan-AM gravimeter measurement was ¢grel <sup>D</sup> 0.03 mGal. Since the measurements of the gravimeters GABL-PM and Chekan-AM are independent, the instrumental accuracy of continuous measurement of the absolute gravity obtained with the mock-up integrated system can be estimated from the formula:

$$
\sigma\_{\text{g\_{int}}} = \sqrt{\sigma\_{\text{g\_{abs}}}^2 + \sigma\_{\text{g\_{rel}}}^2} = 0.04mGal\dots
$$

During three nights we put the zenith telescope on the gyrostabilizer and provided measurements outside the test room. The results of DOV measurements are shown in Fig. 5.

Standard deviation of the DOV measurements did not exceed 0.500. So installation of the zenith telescope on the gyrostabilizer do not affect on the accuracy of the system on the fixed base.

Another aim of the bench tests was to estimate the accuracy of absolute gravity measurements on a weakly oscillating base. We didn't provide dynamic tests of the gyrostabilized zenith telescope because it is impossible to observe stars at the test room.

As can be seen from Fig. 6, the spread in the GABL-PM readings is 200 mGal even with minimum heaving with an amplitude of 0.2 m and a period of 200 s, which additionally confirms the necessity of using both types of gravimeters for occasional determination of absolute gravity aboard a vessel.

The data from the relative gravimeter were used to determine vertical accelerations *Wz*, which made it possible to choose from a set of discrete gravity measurements of the absolute gravimeter, the measurements that were performed under the best conditions. The criterion for choosing *g* measurements was the limiting value of vertical accelerations of less than 3 mGal. It took about 2 h to obtain a set of 100 reliable measurements of *g* at heaving with an amplitude of 0.2 m and a period of 200 s, and the standard deviation of the GABL-PM gravimeter readings was 0.91 mGal. However, as can be seen from Fig. 7, the difference in the mean values of the measurements of the GABL-PM gravimeter before, during and after heaving by modulus was 0.11 mGal,

**Fig. 4** The results of gravity measurements on a fixed base

**Fig. 5** The results of DOV measurements on a fixed base

**Fig. 6** Measurements of the mock-up gravimetric system on a weakly oscillating base

**Fig. 7** Measurements results obtained with the absolute gravimeter on a fixed base and a weakly oscillating base



**Fig. 8** System breadboard aboard the vessel

which was the first reliable demonstration of the feasibility of measuring the absolute gravity on a weakly oscillating base.

The final stage of the bench tests was the breadboard testing on the sea-wave simulator, which allows specifying three-component angular motions of a mobile platform with the equipment under test fixed to it. The method for estimating the breadboard accuracy at the sea-wave simulator was identical to that described above for testing at the bench of vertical displacements. The results of two series of mock-up tests on dynamic benches are summarized in Table 1.

From the results of the bench tests it can be seen that imitation of smooth sea waves does not significantly affect the measurements of the absolute gravity with the breadboard gravimetric system. The standard deviation of the heaving measurements was less than 1 mGal, and the difference in the mean values in statics and dynamics did not exceed 0.1 mGal. Based on the positive results of the bench tests, we proceeded to the next stage – the sea trials of the breadboard gravimetric system.

#### **6 Results of the Sea Trials**

The system breadboard was tested at Priozersk seaport and in the Ladoga Lake from aboard a vessel with displacement of 120 tons (Fig. 8). The scope of tests included the measurements of the absolute value of gravity and DOV on moored vessel and during the vessel positioning at some points of the water area.

Before installing the absolute gravimeter and the zenith telescope on board, gravity and DOV were measured at the pier on a massive concrete base. These measurements were further used as reference. To account for the difference in altitudes between the location of the gravity measurements on a concrete base and the place where the absolute gravimeter was installed on the vessel, we performed leveling. The difference in the altitudes was 2.1 m. That was taken into account in the further data processing.

The sea tests started from measuring the absolute value of gravity on the vessel moored at the pier. Within 1.5 h, a required data set of 100 reliable measurements taken at the inertial disturbances less than 3 mGal was formed. The difference between the gravity measurements at the pier and aboard the vessel was less than 0.1 mGal (Fig. 9). This result confirmed the possibility of offshore surveys without referencing to onshore gravity stations.

The second stage of the sea tests was verification of gravity measurement accuracy when positioning the vessel at a point of the Ladoga Lake water area. The measurements were taken for 2 h at the sea state 2. The inertial accelerations reached 1.5 Gal, and the vessel pitching/rolling was up to 5ı. In spite of significant accelerations (Fig. 10),

**Fig. 9** Breadboard measurements at the pier and aboard the moored vessel

**Fig. 10** Gravity measurements at water area point

25 reliable measurements were obtained. The difference of gravity values relative to the pier, calculated by the readings of the absolute and relative gravimeters, was C4.20 and C5.22 mGal, respectively.

The obtained results confirmed the possibility of measuring the absolute value of gravity in real sea conditions at the sea state up to 2 points, using a system comprising the absolute and relative gravimeters.

DOV measurements with the gyrostabilized zenith telescope were carried out in the water area of the Ladoga Lake, at two points located 100 m away from each other, on four nights. The results of the measurements were compared to the known values of DOV according to the global model EGM-2008, since they can be used as a reference in a region with small gravity anomalies and smooth topography (Fig. 11).

The standard deviations of DOV measurements were 0.65<sup>00</sup> (Ÿ-component) and 0.85<sup>00</sup> (˜-component). The differences with EGM-2008 data were 0.96<sup>00</sup> (Ÿ-component) and 0.36 (˜-component). So the results of the field tests demonstrated the possibility of rapid high-precision measurement of the absolute values of DOV in marine conditions, using the astrogeodetic method.

#### **7 Conclusions**

The first results of experimental studies confirmed the possibility of high-precision measurement of the absolute values of gravity and deflection of the vertical at sea. The use of an integrated system will improve the quality of geodetic data in regions of the World Ocean with large gravity anomalies.

**Fig. 11** The results of DOV measurements at water area point

**Acknowledgement** This work was supported by the Russian Science Foundation, project no. 18-19-00627.

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/ licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Multi-Scenario Evaluation of the Direct Method in Strapdown Airborne and Shipborne Gravimetry**

Felix Johann, David Becker, Matthias Becker, and E. Sinem Ince

#### **Abstract**

In recent years, it was shown that the quality of strapdown airborne gravimetry using a navigation-grade strapdown inertial measurement unit (IMU) could be on par with "classical" airborne gravimeters as the 2-axis stabilized LaCoste and Romberg S-type gravimeter. Basically, two processing approaches exist in strapdown gravimetry. Applying the indirect method (also referred to as "inertial navigation approach" or "one-step approach"), all observations – raw GNSS observations or position solutions, IMU specific force and angular rate measurements – are combined in a single Kalman Filter. Alternatively, applying the direct method (also referred to as "accelerometry approach" or "cascaded approach"), GNSS position solutions are numerically differentiated twice to get the vehicle's kinematic acceleration, which is then directly removed from the IMU specific force measurement in order to obtain gravity. In the scope of this paper, test runs for the application of strapdown airborne and shipborne gravimetry are evaluated using an iMAR iNAV-RQH-1003 IMU. Results of the direct and the indirect methods are compared to each other. Additionally, a short introduction to the processing scheme of the Chekan-AM gravimeter data is given and differences between Chekan-AM and strapdown results of the shipborne campaigns are analysed. Using the same data set, the cross-over residuals suggest a similar accuracy of 0.39 mGal for the Chekan-AM and 0.41 mGal for the adjusted strapdown results (direct method).

#### **Keywords**

Airborne gravimetry - Gravity disturbance - IMU - Shipborne gravimetry -Strapdown

F. Johann (-) · M. Becker

e-mail: johann@psg.tu-darmstadt.de; becker@psg.tu-darmstadt.de D. Becker

School of Earth and Planetary Sciences, Curtin University, Perth, WA, Australia

#### e-mail: david.becker@curtin.edu.au

#### E. S. Ince

#### **1 Introduction**

Airborne and shipborne gravimetry enables gravity determination with a medium accuracy, compared to satellite and terrestrial gravimetry. Especially in remote areas, where terrestrial measurements are cost-intensive and time-consuming or even impossible, it is widely used to enhance the accuracy and resolution of satellite-based gravity data. In Polar Regions without satellite gravity data, it is the only practicable gravity determination approach.

While shipborne gravimetry using gas-pressured, pendulum or spring gravimeters had already been employed in the first half of the twentieth century (Marson 2012; Nabighian

J. T. Freymueller, L. Sánchez (eds.), *5th Symposium on Terrestrial Gravimetry: Static and Mobile Measurements (TG-SMM 2019)*, International Association of Geodesy Symposia 153, https://doi.org/10.1007/1345\_2020\_127

Physical and Satellite Geodesy, Technische Universität Darmstadt, Darmstadt, Germany

Global Geomonitoring and Gravity Field, Helmholtz-Centre Potsdam, GFZ German Research Center for Geosciences, Potsdam, Germany e-mail: elmas.sinem.ince@gfz-potsdam.de

et al. 2005), first airborne gravity test flight results using a horizontally stabilized spring gravimeter have not been published before 1960 (Nettleton et al. 1960). In the 1990s, a major improvement in the obtained accuracy followed the advent of the Global Positioning System (GPS), followed by further Global Navigation Satellite Systems (GNSS). The application of inertial measurement units (IMUs) as gravimeters could be shown to be capable of generating results on an accuracy level similar to that of the "classical" approach, where horizontally stabilized spring gravimeters are used (Glennie et al. 2000; Becker et al. 2016).

Main advantages of the IMU usage in the so-called "strapdown" approach are a lower instrument weight, less space and power requirements as well as lowered purchase and maintenance costs. Furthermore, the measurement principle enables vector (3-D) gravimetry and the system is less impaired by turbulences and flight manoeuvres. Strong sensor drifts affecting strapdown gravimetry can be reduced using calibration methods that can lower their impact significantly. It was demonstrated that the application of a simple warm-up calibration of the IMU's vertical accelerometer can eliminate most of the sensor drift (Becker et al. 2015b; Becker 2016).

Even though the vehicle dynamics differ in airborne and shipborne gravimetry, processing strategies virtually agree.

This paper gives an overview on strapdown gravimetry and recent campaigns with participation of the chair of Physical and Satellite Geodesy (PSGD) of the Technical University of Darmstadt. An introduction in the processing approaches in strapdown gravimetry is given in Sect. 2. To enable an evaluation of the presented strapdown result's exterior accuracy, the data processing in "classical" kinematic gravimetry is briefly presented in Sect. 3 using the example of the Chekan-AM gravimeter of the German Research Centre for Geosciences (GFZ Potsdam). The airborne and shipborne campaigns evaluated in the scope of this paper are discussed in Sects. 4 (campaign introduction), 5 (strapdown processing results) and 6 (comparison to Chekan-AM results), followed by the concluding Sect. 7.

#### **2 Processing Approaches in Strapdown Gravimetry**

There exist two fundamentally different processing approaches in (kinematic) strapdown gravimetry: the indirect and the direct method. This paper shortly introduces some basics of both approaches; for a more detailed view, the reader is referred to Becker (2016) and Johann et al. (2019).

The "indirect" method as it was introduced by Jekeli (2001) has also been called "inertial navigation approach" (Ayres-Sampaio et al. 2015), "traditional way" (Kwon and Jekeli 2001) or "one-step approach" (Becker 2016). Here, all observations from an IMU and a GNSS receiver are integrated in a Kalman filter with a state vector including position, velocity, attitude, sensor biases and gravity. In the prediction step, the previous epoch's state vector is propagated using the IMU observations (accelerations and angular velocities). If GNSS observations are available as pseudoranges and carrier phases (tightly coupled integration) or position and velocity solutions (loosely coupled integration), the Kalman filter will form a weighted average between the predicted and the actual observations. The weighting depends on the accuracy ratio of the predicted states and GNSS observations as reflected by the stochastic model.

"Accelerometry approach" (Ayres-Sampaio et al. 2015; Kwon and Jekeli 2001), or "cascaded approach" (Becker 2016) have been alternative terms for the "direct" method. In contrast to the indirect method, where gravity determination is done in the position domain, in the direct method, vertical gravity is directly obtained in the acceleration domain. The gravity disturbance

$$\boldsymbol{\mathfrak{\delta g}}^{n} = \vec{r}^{n} - \mathbf{C}\_{b}^{n} f^{b} + \boldsymbol{\mathfrak{\delta g}}\_{\mathrm{eot}}^{n} - \mathbf{y}^{n} \tag{1}$$

in the navigation frame being the difference of gravity and normal gravity <sup>n</sup> is directly obtained by subtracting the specific force C <sup>n</sup> <sup>b</sup>f <sup>b</sup> from the kinematic acceleration r<sup>R</sup> <sup>n</sup> (Wei and Schwarz 1998). All these terms are given in the navigation frame n whose origin is defined to be at the IMU's centre of observations with its orthogonal axes pointing towards North, East and Down. The specific force f <sup>b</sup> as measured by the IMU accelerometers in the body frame b is transformed to the navigation frame by the multiplication with the rotation matrix C <sup>n</sup> <sup>b</sup> . Alternatively, quaternions can be used to represent the attitude. The GNSS position solutions are numerically differentiated twice in order to obtain the kinematic acceleration. The vehicle's attitude used as input for the rotation matrix is computed in a GNSS/IMU integration algorithm. The addition of the velocity-dependent Eötvös correction ıg<sup>n</sup> eot eliminates the impacts of the Coriolis acceleration and the centrifugal acceleration caused by the movement of the vehicle relative to the Earth. The obtained gravity disturbance needs to be low-pass filtered in order to eliminate high-frequency noise.

On the one hand, the indirect method implements a wellknown optimal estimation procedure, the Kalman filter, and is commonly deemed to be the more rigorous approach (Becker 2016). On the other hand, if the direct method is applied, the GNSS/IMU integration can be performed using available commercial navigation software (e.g. the NovAtel Waypoint InertialExplorer 8.60 used in the scope of this paper), significantly reducing the burden of software development.

In order to get the best out of the strapdown data, several corrections need to be applied. A detailed description can be found in Johann et al. (2019).


In the scope of this paper, the precision of the results is examined using a cross-over analysis. For this, the survey plan includes intersecting lines. When the intersections are at nearly the same height and one assumes that there is no significant temporal gravity variation, the root mean square (RMS) of the cross-over residuals at these intersections is a measure for the precision of the results. Assuming an equal accuracy on both lines, the RMS error (RMSE), which is given by the RMS divided by <sup>p</sup>2, indicates the uncertainty of the obtained gravity disturbance values at the lines (Forsberg and Olesen 2010; Becker 2016).

Since the gravity disturbance results have already been fixed to the reference values at the airports/harbours, the linear trend is easily removed from the results. If a sufficient number of cross-over points is available, the impact of further drifts during a measurement flight/cruise can be reduced by a cross-over adjustment (also called "levelling"). During the least-squares adjustment, a bias per line is estimated. If the number of cross-over points is low, the resulting precision value is susceptible to being too optimistic, which can be avoided by the use of correction factors (Becker 2016; Johann et al. 2019).

#### **3 Chekan-AM Data Processing at GFZ**

To enable a comparison with another measurement system, processing methods and some results using a Chekan-AM gravimeter are presented in this paper in addition to the strapdown analysis.

Since 2011, the GFZ Potsdam has been involved in various shipborne and airborne gravimetry campaigns with its mobile gravimeter Chekan-AM whose working principle is based on angle variation measurements of two quartz sensors that are positioned in a viscous liquid. The performance of this equipment has been verified in GFZ's marine test campaign in Lake Müritz, and other dedicated campaigns such as Lake Constance and the GEOHALO mission (Petrovic et al. 2015). The same equipment is further reliably used in the shipborne gravimetry campaigns in the Baltic Sea within the subactivity 2.1 of the "Finalising Surveys for the Baltic Motorways of the Sea" (FAMOS) project (Swedish Maritime Organisation 2019).

The recovery of the measurements in terms of accelerations is based on the mathematical model and calibration constants provided by the manufacturer Electropribor, St. Petersburg. It is known that the survey vehicle and the measurement conditions have a strong impact on the gravimeter measurements. In principle, the processing scheme presented in Krasnov et al. (2011) is followed in GFZ's data processing routine but the processing needs to be adopted based on the quality of the raw measurements and measurement conditions. Obviously, in mobile gravimeter measurements, the raw recordings need to be converted into acceleration units. Like in strapdown gravimetry, measurements need to be corrected for the Eötvös effect and, most importantly, low-pass filtered to eliminate the high frequency noise components and retrieve meaningful results. Occasionally, different low pass filter lengths can be applied to different quality of datasets. Based on their experiences, GFZ found that, in shipborne gravimetry, a filter length of 400 s provides satisfying results in terms of measurement accuracy w.r.t. global models and other campaign results computed in a cross-over analysis and is a good tradeoff between the spatial resolution and final data accuracy (Lu et al. 2019).

One of the most challenging tasks of relative gravimetry is the drift estimation of the gravimeter, which needs to be computed as good as possible and taken into account in the data processing. The gravimeter recordings at the harbours, where reference gravity values are available, are taken before and after completing relevant tracks and are used to compute the drift behaviour of the instrument sensors. After correcting the drift, since the measured values are not absolute gravity but gravity differences measured along the measurement track they need to be tied to absolute gravity values at reference stations that are ideally co-located to the harbour tie points.

It is worth mentioning that the Chekan-AM measurements even after applying a low pass filter are not of very high quality during the turns of the ship. Therefore, these periods lasting about from few to tens of minutes as well as disturbed measurements due to harsh sea conditions or instrumental errors (e.g. unexpected drift behaviour, sensor aging) are manually detected and removed from the final delivered products. Good quality measurements (e.g. performed under optimum measurement conditions during dedicated gravimetry campaigns) may not require detailed investigations.

#### **4 Investigated Airborne and Shipborne Strapdown Gravimetry Campaigns**

Since 2013, PSGD co-organized or supported various strapdown gravimetry campaigns. An IMU of the type iMAR iNAV-RQH-1003 (Fig. 1, bottom right) (IMAR Navigation 2012) was used as gravimeter. The IMU was connected to GNSS antennas for the purpose of time synchronization. The following campaign results are presented in this paper:

– "MY2014": In August 2014, 12 flights were undertaken above the South China Sea northwest of Borneo, Malaysia, using a medium-size aircraft, type Beechcraft King Air 350 (Fig. 2a), with autopilot. For details,

**Fig. 1** Chekan-AM and iMAR iNAV-RQH-1003 with iTempStab-AddOn aboard DENEB 2018

including a scatter plot of the results and an elaborate evaluation, see (Johann et al. 2019).


**Fig. 2** Vehicles: (**a**) Beechcraft King Air 350 (MY2014); (**b**) Grob G109B (ODW2017); (**c**) Cessna 206 "Stationair 6" (ODW2018); (**d**) German research vessel DENEB (BTS2017/18)



aBTS2017: Three cruises are remaining if one of the four cruises was rejected assuming an outlier

bBTS2017: The heading-dependent correction has not been applied in the indirect method

cThe cross-over adjustment of the BTS campaigns was performed stint-wise (complete run from harbour to harbour), while the other campaigns were adjusted line-wise. Usually, lower RMSE values are obtained via line-wise adjustment

housing, the iMAR iTempStab-AddOn, which stabilises the IMU's temperature using two Peltier elements.

Dual frequency GNSS receivers tracking GPS and GLONASS were installed in all campaigns. Positions were calculated in the Precise Point Positioning (PPP) mode with the NovAtel Waypoint InertialExplorer 8.60 using precise satellite orbit and clock products by the Center for Orbit Determination in Europe (CODE). The campaign statistics are summarised in Table 1. The turbulence metric and the obtained precision values (RMSE) will be discussed in Sect. 5.

#### **5 Strapdown Results**

The obtained vertical gravity disturbance without crossover adjustment for both Odenwald and both Baltic Sea campaigns are illustrated in Figs. 3, 4, 5, 6. For the Malaysia

**Fig. 3** Vertical gravity disturbance [mGal] ODW2017 (Map data: Google, GeoBasis-DE/BKG, Europa Technologies)

campaign, equivalent figures including preliminary results of the horizontal components can be found in (Johann et al. 2019). RMSE obtained in cross-over analyses with and

**Fig. 4** Vertical gravity disturbance [mGal] ODW2018 (Map data: Google, GeoBasis-DE/BKG, Europa Technologies)

**Fig. 5** Vertical gravity disturbance [mGal] BTS2017 (Map data: Google, GeoBasis-DE/BKG, Landsat/Copernicus)

**Fig. 6** Vertical gravity disturbance [mGal] BTS2018 (Map data: Google, GeoBasis-DE/BKG, Landsat/Copernicus)

without cross-over adjustment for all four campaigns are noted in Table 1, where available results of the indirect method of strapdown gravimetry are given in parentheses. For the Malaysia 2014 campaign, the results of the direct and the indirect method are on a par with an RMSE of about 1.3 mGal without and 0.65 mGal with line-wise crossover adjustment. For the Odenwald 2017 campaign (Fig. 3), results are considerably worse with an RMSE of 3.1 mGal without and 1.1 mGal with adjustment. The indirect method performs slightly better (3.0 and 0.9 mGal). Without adjustment, the results for Odenwald 2018 (Fig. 4) are even worse (4.0 mGal), maybe caused by temperature calibration problems due to the low environment temperature. In contrast, after adjustment, the RMSE is at the precision level of the Malaysia campaign (0.64 mGal).

For the Baltic Sea 2017 campaign (Fig. 5), without adjustment, an RMSE of 1.28 mGal is obtained applying the direct method and 1.18 mGal applying the indirect method. It should be noted that the empirical heading-dependent correction was not used in the indirect method. Hence, further improvements might be possible. After a cruisewise cross-over adjustment, RMSE improved to 1.00 mGal (direct method) and 0.81 mGal (indirect method). A linewise adjustment was not possible, since isolated nets would occur causing a datum defect. One of the four cruises has been identified as an outlier. The rejection resulted in an improved accuracy of 0.81 mGal without adjustment (direct method). A cruise-wise adjustment did not enhance the results after the outlier rejection. In the Baltic Sea 2018 campaign (Fig. 6) with the installed temperature stabilisation housing, significantly better results with an RMSE of 0.71 mGal before and 0.55 mGal after the cruise-wise adjustment have been noted. Again, a line-wise adjustment was not applicable. To test the repeatability of the results, both Baltic Sea campaigns have also been examined together (without the outlier cruise in 2017) leading to RMSE of 0.87 mGal without and 0.66 mGal after adjustment, medium accuracies compared to the campaign's single results.

The accuracy of airborne gravimetry results is suspected to depend on the flight turbulence level. A turbulence metric that directly indicates the reaction of the aircraft turbulences is the RMS-g. It is defined as the moving standard deviation of the vertical acceleration during the flight (Becker 2016). Other turbulence metrics like the Eddy Dissipation rate, which rate the atmospheric turbulence state, eliminate some aircraft-dependent effects acting on the IMU. In the scope of this paper, the RMS-g is computed applying a 50 s moving window on 1 Hz GNSS vertical acceleration data. For the presented campaigns, a trend between cross-over RMSE and RMS-g might be indicated by the results given in Table 1 if the non-adjusted RMSE of ODW2018 was declared as an outlier caused by temperature calibration issues. Note that accuracy differences depend on multiple campaign parameters like the local gravity field, the low-pass filter length and GNSS data quality.

#### **6 Comparison to Chekan-AM Results**

Using the cross-over analysis algorithm implemented in the direct method of strapdown gravimetry, an RMSE of 0.33 mGal and 0.39 mGal resulted for the Chekan-AM data of the Baltic Sea campaigns 2017 and 2018, respectively. To enable a comparison, line segments affected by harsh sea conditions have also been removed from the strapdown results resulting

**Fig. 7** Comparison of estimated vertical gravity disturbances for BTS2018, cruise 2: (**a**) blue: strapdown (direct), green: Chekan-AM; (**b**) Difference (RMS: 0.9 mGal)

in a slightly improved RMSE for BTS2017 (improvement of about 0.1 mGal compared to Table 1). For BTS2018, the strapdown RMSE decreased to 0.52 mGal before and 0.41 mGal after cruise-wise adjustment. Regarding the absolute differences to the non-adjusted strapdown results (direct method), the RMS are 1.24 mGal in 2017 (without the outlier cruise) and 0.97 mGal in 2018. This indicates remaining systematic differences, possibly due to different sensor drift behavior. The comparison of the obtained vertical gravity disturbances of cruise 2, BTS2018, are shown in Fig. 7.

#### **7 Conclusions**

After applying several corrections, the precision of the presented strapdown gravimetry campaigns is between 0.7 and 4.0 mGal without cross-over adjustment. The best RMSE for an airborne campaign (1.26 mGal) was reached at Malaysia 2014, which was conducted with a medium-size aircraft with autopilot under relatively stable flight conditions.

After a line-wise adjustment, the airborne RMSE obtained at Malaysia 2014 and Odenwald 2018 were at the same level (about 0.6 mGal). The corresponding precision of the Odenwald 2017 campaign is significantly worse (1.1 mGal), probably suffering from GNSS acquisition problems.

The successful combination of both Baltic Sea campaigns indicates a good repeatability of the direct method of strapdown gravimetry. In the implementations created at PSGD, the accuracy of this method is on a par with indirectly obtained results for most evaluated campaigns. Many additional airborne and shipborne surveys processed with both methods would be required to enable statistically firm statements about the accuracy ratio of both approaches under various conditions.

The first test of PSGD's new thermal housing at the Baltic Sea 2018 campaign indicates that its installation improves the long-term stability significantly. Comparing line segments without harsh sea conditions, the precision of the cruise-wise adjusted strapdown gravity disturbance results of 0.41 mGal is competitive with the Chekan-AM results (0.39 mGal). Remaining systematic differences between the strapdown and Chekan results are a subject of research.

Advantages of the strapdown approach like the possibility of vector gravimetry let this approach be an auspicious alternative to the classical approach using horizontally stabilized spring gravimeters.

**Acknowledgements** The financial support by the Deutsche Forschungsgemeinschaft (DFG) as a part of the project "Innovative calibration methods for strapdown airborne vector gravimetry aboard HALO" within the framework of SPP 1294 "Atmospheric and Earth System Research with the "High Altitude and Long Range Research Aircraft" (HALO)" is gratefully acknowledged. The Malaysia airborne gravity campaign was financed by the Department of Survey and Mapping, Malaysia (JUPEM) and was conducted in cooperation with the Technical University of Denmark (DTU Space). The Odenwald 2017 campaign was conducted by PSGD in cooperation with DTU Space, the Technical University of Dresden (TU Dresden) and the Technical University of Munich (TU Munich). The Odenwald 2018 campaign was conducted by iMAR Navigation in cooperation with DTU Space. The DENEB 2017 and 2018 campaigns have been carried out by the German Federal Agency for Cartography and Geodesy (BKG), the German Research Centre for Geosciences (GFZ Potsdam) and the German Federal Maritime and Hydrographic Agency (BSH) within subactivity 2.1 of the project "Finalising Surveys for the Baltic Motorways of the Sea" (FAMOS), co-financed by the European Union. Joachim Schwabe, BKG, adapted the GFZ Chekan-AM results of the 2017 Baltic Sea campaign. The Malaysia airborne gravity campaign was financed by the Department of Survey and Mapping, Malaysia (JUPEM) and was conducted in cooperation with the Technical University of Denmark (DTU Space).

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons. org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Improving Gravity Estimation Accuracy for the GT-2A Airborne Gravimeter Using Spline-Based Gravity Models**

Vadim S. Vyazmin, Yuri V. Bolotin, and Anton O. Smirnov

#### **Abstract**

In this article, we propose a technique for processing airborne gravimetry measurements at repeat drape lines. To model gravity along a line, we use cubic B-splines. The unknown parameters of the model (the spline coefficients) are to be estimated from the measurements. In order to take into account dependance of gravity on position, the same set of the spline coefficients is assumed for gravity modeling at each repeat line. Gravity estimation is performed with the Kalman filter. The spline coefficients are included in the state vector of the filter as unknown constants. The developed estimation algorithm was tested using measurements of the GT-2A gravimeter collected at repeat drape lines. The gravity estimates obtained with the use of B-splines are more accurate compared with those obtained with the standard technique, in which gravity is modeled as a stochastic process in the time domain.

#### **Keywords**

Airborne gravimetry - Gravity disturbance - GT-2A gravimeter - Kalman filter - Repeat drape lines -Splines

#### **1 Introduction**

Airborne scalar gravimetry is widely used to obtain high accuracy information about the Earth gravity in hard-toreach regions (shelves, mountains, polar areas, etc.). A typical platformed airborne gravimetry system consists of a gravity sensor (a single-axis accelerometer) mounted on a gyro-stabilized platform and two (or more) global navigation satellite system (GNSS) receivers (one on board the aircraft, others on the ground). Postprocessing of airborne gravimetry measurements includes three main stages such as determination of GNSS positioning and velocity solutions, estimation of the platform tilt angles via the inertial navigation system (INS) and GNSS integration, and estimation of gravity using the Kalman filtering technique or low-pass filtering. Both

V. S. Vyazmin (-) · Y. V. Bolotin · A. O. Smirnov filtering techniques are based on representing gravity in the time domain (e.g., (Lu et al. 2017) or modeling gravity as a stochastic process (Bolotin and Golovan 2013; Kwon and Jekeli 2002; Stepanov et al. 2015)).

In this research, we focus on the last stage of postprocessing (i.e., gravity estimation) given airborne measurements at repeat drape lines. Using a stable platform gravimetry system (e.g., the GT-2A (Berzhitzky et al. 2002)) in a drape flight is challenging, since such a system is more suitable for benign dynamics (Studinger et al. 2008). For the GT-2A, underestimation of the scale factor of the gravity sensor can be observed during postprocessing of drape lines. For this reason, specific calibration lines are flown, at which the scale factor is reliably estimated. However, accuracy of gravity estimation at repeat drape lines is still lower than in the case of a flight at a constant height. To our opinion, the standard approach to gravity estimation based on modeling gravity in the time domain is not well suitable for processing repeat lines, as gravity modeled along a line, strictly speaking, differs from gravity modeled along a repeat line. Actual gravity is, of course, the same at repeat lines assuming

Navigation and Control Laboratory, Department of Applied Mechanics and Control, Lomonosov Moscow State University, Moscow, Russia e-mail: v.vyazmin@navlab.ru

that they were flown at the same altitude. Therefore, in the case of a repeat line flight, it is preferable to model gravity in the spatial domain rather than in the time domain.

In this work, we propose an approach based on spatial modeling of gravity to process repeat drape lines flown over the same ground track. We parameterize the gravity disturbance along a line using B-splines of order three. To take into account dependance of gravity on position, we use the same set of unknown parameters of the model (the spline coefficients) for gravity modeling at each repeat line. The model is easily described in the state space assuming that the spline coefficients are constant over the time during the aircraft's flight. The standard stochastic estimation problem is posed, and solved via Kalman filtering. The state vector for the filter includes the spline coefficients and the parameters of the gravimeter systematic errors (the calibration parameters).

The developed approach was applied to processing airborne measurements collected with the GT-2A system during a repeat drape line flight in South Africa in 2009. The obtained gravity estimates are compared with those provided by the standard approach based on modeling gravity in the time domain, and with the gravity values derived from the global model EGM2008. An increase in accuracy in gravity estimation using the new approach is shown.

#### **2 Gravity Estimation Methodology**

*Basic Equation and Gravimeter Measurements* The gravity sensor of the GT-2A is mounted on a gyrostabilized platform, which makes the instrument sensitivity axis to be vertical. The gravity disturbance ıg is determined from the fundamental equation of airborne scalar gravimetry, which is written as:

$$
\dot{V}\_3 = -\delta \mathbf{g} - \mathbf{g}\_0 + \mathbf{g}\_{Evv} + f\_3 \tag{1}
$$

where V3 is the vertical velocity of the gravimeter proof mass in the ENU coordinate frame, g0 is the normal gravity at the flight height, gEtv is the Eötvös correction term, f3 is the vertical projection of the specific force.

Raw measurements of the gravimeter f <sup>0</sup> <sup>3</sup> can be described by the equation

$$(1 + k\_3)f\_3' + \mathfrak{r}f\_3' = f\_3 - \delta f\_{hor} + \delta f\_{illt} + \delta f\_{drift} + q\_f \dots$$

Here k3 is the scale factor error, is the time delay of gravimeter measurements (due to physical delay and averaging of gravity sensor measurements performed in the internal software of the gravimeter), ıfhor is the horizontal acceleration influence on f <sup>0</sup> <sup>3</sup> due to imperfect placement of the gravimeter on the platform, i.e., more explicitly, ıfhor D k1f <sup>0</sup> <sup>1</sup> C k2f <sup>0</sup> <sup>2</sup> , where f <sup>0</sup> <sup>1</sup> , f <sup>0</sup> <sup>2</sup> are the platform accelerometer measurements, k1, k2 are unknown factors that can differ from one flight to another. Further, ıftilt is the horizontal acceleration due to the platform tilt error (it is determined from INS-GNSS integration), ıfdr if t is the sum of the gravity sensor bias and drift. Assuming linear drift, the term ıfdr if t can be computed given pre- and post-flight measurements and terrestrial gravity value. The term qf is the noise.

The frequency of the GT-2A raw measurements is 18.75 Hz, and the time delay is about 2 s.

As the terms ıftilt and ıfdr if t can be determined and compensated for using standard procedures, we will omit them further on.

*State-Space Modeling* Rewrite Eq. (1) replacing unknown variables on the right-hand side of the equation by measurements and omitting the gravity disturbance term:

$$
\dot{V}'\_3 = -\mathbf{g}\_0 + \mathbf{g}\_{Etv} + f'\_3 \tag{2}
$$

where g0 and gEtv are known from the GNSS position and velocity solutions (assuming that the lever arm effect is compensated for).

Denote by V3 the vertical velocity error defined as (Bolotin and Golovan 2013)

$$
\Delta V\_3 = V\_3 - V\_3' - \mathfrak{r} \, f\_3'. \tag{3}
$$

Then we obtain the error dynamics equation of airborne scalar gravimetry by subtracting Eq. (2) from Eq. (1) and using the above gravimeter measurement model

$$
\Delta \dot{V}\_3 = -\delta \mathbf{g} + k\_1 f\_1' + k\_2 f\_1' + k\_3 f\_3' - q\_f \tag{4}
$$

where k1; k2 are the residual platform misalignment errors in the sum with errors of placement of the gravity sensor on the platform, and k3 is the scale factor error of the gravity sensor.

Representing the GNSS velocity solution as V GNSS <sup>3</sup> D V3 C evel , where evel is the GNSS velocity random error, and introducing the velocity observation as <sup>y</sup> <sup>D</sup> <sup>V</sup> GNSS <sup>3</sup> V <sup>0</sup> <sup>3</sup> , we obtain the observation equation

$$
\mathbf{y} = \Delta V\_3 + \mathbf{r} f\_3' + e\_{vel} \,. \tag{5}
$$

Equations (4) and (5) contain six unknown variables:

$$\Delta V\_3, \delta \mathbf{g}, \mathbf{r}, k\_1, k\_2, k\_3.$$

The variables qf and evel in Eqs. (4) and (5) are assumed to be zero-mean white noises with known variances. We adopt the following models for the gravimeter systematic error parameters (calibration parameters):

$$
\dot{\pi} = q\_{\pi}, \quad \dot{k}\_1 = q\_{k\_1}, \quad \dot{k}\_2 = q\_{k\_2}, \quad \dot{k}\_3 = q\_{k\_3} \tag{6}
$$

where q- ; qk1 ; qk2 ; qk3 are zero-mean white noises with known (sufficiently small) variances.

By introducing an equation for ıg (a gravity model), we obtain the system of equations in the state-space form. The gravity estimation algorithm includes the following steps:


#### **3 Modeling Gravity Along Repeat Lines Using B-Splines**

To arrive at estimation algorithm, we introduce a spatial gravity model assuming a repeat line flight. We use Bsplines of order three for gravity modeling, which form a basis in the space of cubic splines in R<sup>1</sup> and each function of which has finite support (Fig. 1). Recall that the B-spline, by definition, is a fourfold convolution of a rectangular function with itself (De Boor 2001). We assume that the repeat lines are flown along the same ground track at the same altitude and have equal length <sup>L</sup>. Let <sup>x</sup> <sup>D</sup> x.t / <sup>2</sup> <sup>R</sup><sup>1</sup> denote the distance travelled by the aircraft (more precisely, the gravity sensor proof mass) along the first line. Then x varies from 0 to L. For simplicity, assume that the subsequent repeat lines coincide with the first line. Thus, x varies from 0 to L or from L to 0 at the repeat lines, depending on the direction of aircraft's flight.

We represent the gravity disturbance at a point of a line as follows

$$\delta \mathbf{g}(\mathbf{x}) = \sum\_{i=0}^{N-1} c\_i \mathcal{B}\_i(\mathbf{x}) \tag{7}$$

**Fig. 1** Cubic B-splines

where Bi.x/ is the cubic B-spline, ci is the i-th spline coefficient to be estimated, N is the total number of the spline coefficients. Note that the same spline coefficients ci , i D 0; : : : ; N 1, are used for describing gravity at each repeat line. Note also that the B-spline Bi.x/ in fact is the function of .x <sup>x</sup><sup>S</sup> <sup>i</sup> /=x, where x<sup>S</sup> <sup>0</sup> ;:::;x<sup>S</sup> N-<sup>1</sup> are the knots located at the first line, x is the length of the knot step (measured along the line). The value of x is to be chosen.

Assuming that the spline coefficients are constant over the time during the flight, we can easily transform the gravity model Eq. (7) into the state-space form by writing the obvious equations

$$
\dot{c}\_i = 0, \quad i = 0, \dots, N - 1. \tag{8}
$$

Therefore, we obtain the complete state-space system described by Eqs. (4)–(6) and (8) given the gravity model (7). Using Kalman filtering, the optimal (in the minimum mean squared error sense) estimates of the spline coefficients, ci , i D 0; : : : ; N 1, the velocity error V3 and the calibration parameters, -; k1; k2; k3, are obtained. The total dimension of the state vector is N C 5. At the initial iteration of the filter, the a priori variance of the state vector, including the variance of the spline coefficients, should be chosen.

*Transfer Function of the Optimal Filter* The Fourier transform of a cubic B-spline can be written as (De Boor 2001):

$$\hat{B}\_0(\alpha) = \left(\frac{e^{j\alpha T} - 1}{j\alpha T}\right)^4$$

where T D x=V is the time step of the knots, V is the average speed of the aircraft at the lines, <sup>j</sup> <sup>D</sup> p<sup>1</sup>. By varying the time step of the knots T , one can obtain B-splines with different frequency bandwidths. Since then, it is desirable to determine how the spatial resolution of the gravity estimate provided by the Kalman filter depends on T .

For this, let us derive the expression for the transfer function of the filter. Consider the state-space equations Eqs. (4)–(8) in the discrete-time form denoting by t the measurement time step. Neglect for simplicity the gravimeter calibration parameters and the gravity sensor noise qf in the equations. Let the measurement time step t and the spline knot step T satisfy the relationship T D mt for some integer m 2 and assume no a priori information on the variances of the spline coefficients (i.e., the variance of each coefficient tends to infinity). As above, assume that the GNSS velocity error that appears in Eq. (5) is a white noise. Then it can be shown that the optimal filter transfer function which maps the acceleration observations to the

**Fig. 2** Transfer functions of the optimal filter in the spline-based approach and the standard approach based on the Gauss-Markov model of order 2 for gravity, and the Butterworth filter (two-pass) of order 2

gravity disturbance estimate can be written as (Unser et al. 1993)

1

$$H(z) = \overline{B}\_0^2(z) \left( \frac{\Delta t}{T} \sum\_{k=0}^{m-1} \overline{B}\_0^2(ze^{j2\pi k/m}) \right)^{-1}$$

where B0.*z*/ is the Z-transform of the B-spline of order 4.

The obtained transfer function is close to the Butterworth filter (two-pass) of order 2 (see Fig. 2). Figure 2 also shows the transfer function of the optimal filter in the standard approach to gravity estimation, which is based on stochastic modeling of gravity (e.g., using a Gauss-Markov process of order 2 with infinite correlation time, that is, a process whose second time derivative is a white noise). In this case, the transfer function of the filter is close to the Butterworth filter (two-pass) of order 4.

Thus, with the use of the above expression for the transfer function of the optimal filter (for the spline-based approach), the relationship between the filter cutoff wavelength g and the length of the spline knot step x (or the time step of the knots T ) can be derived as

$$
\lambda\_{\mathfrak{g}} \approx \mathfrak{Z} \Delta \mathfrak{x}, \quad \Delta \mathfrak{x} = VT. \tag{9}
$$

#### **4 Results**

*Survey Flight Overview* The developed approach based on spatial gravity modeling using B-splines was applied to airborne gravimetry data collected within a flight over the Vrederfort Dome impact crater (120 km southwest of Johannesburg, South Africa) on March 21, 2009. The GT-2A gravimetry system used in this survey was

**Fig. 3** Survey area and the ground track of the repeat line flight

installed in a BN-2T Islander aircraft operated by Xcalibur Airborne Geophysics (Olson 2010). Two dual-frequency GPS receivers operating at 10 Hz were used during the flight (the rover on board the aircraft and the base placed on the ground).

The flight consisted of six repeat drape lines (with the ID labels R01–R06) and two calibration lines (with the ID labels R07–R08). The latter were flown in order to obtain more accurate estimate of the gravimeter scale factor error k3. All the lines were flown in the north-south and southnorth direction over the same ground track (Fig. 3). The drape lines were flown at the heights from 1,400 to 1,600 m above the WGS-84 ellipsoid and accurately followed the terrain (Fig. 4). The calibration lines were flown at the height of 3,000 m (Fig. 5).

The length of each line is 85 km. Average flight speed is 57 m/s with standard deviation (STD) 1.3 m/s. The flight duration is 3 h 50 min.

**Fig. 4** (**a**) GPS height at the repeat drape lines; (**b**) Elevation from SRTM-3

**Fig. 5** GPS height at the calibration lines

*Postprocessing Strategy* Postprocessing was performed using the GTNAV/GTGRAV suite of programs developed by the Navigation and Control Laboratory of the MSU (Bolotin and Golovan 2013) and included the following main steps:


At the last step, gravity was modeled in the time domain as a Gauss-Markov process of order 2 with infinite correlation time.

*Accuracy of Estimating the Scale Factor Error* k3 The maximum vertical acceleration at the repeat drape lines (with the Eötvös effect and the normal gravity subtracted) equals 0.1 m/s<sup>2</sup> (Fig. 6). Hence, for gravity estimation with the accuracy of 0.5 mGal the scale factor error k3 should be estimated with the accuracy of 0.00005 or better.

**Table 1** Estimates of the gravimeter calibration parameters: time delay - (s), misalignment errors k1, k2 (arcmin), and scale factor error k3


**Table 2** Standard deviations of the estimate errors for - (s), k1, k2 (arcmin), and k3


Table 1 shows the calibration parameter estimates obtained with the developed spline-based approach from processing the drape lines and the calibration lines. The corresponding standard deviations of the estimate errors provided by the Kalman filter are shown in Table 2. Equal estimation accuracy of 0.00001 for k3 is achieved at the drape lines and at the calibration lines. This accuracy is sufficient for the gravity estimation error due to k3 to be at the level of 0.5 mGal. Moreover, since the two estimates of the scale factor error k3 are close enough to each other (the

**Fig. 6** Vertical acceleration smoothed with the cutoff frequency of 1/100 Hz

**Fig. 7** The along-track gravity disturbance estimates obtained with the two approaches for the cutoff wavelength of (**a**) 100 s and (**b**) 70 s, and gravity from EGM2008; (**c**) GPS height (averaged over six repeat drape lines)

difference equals 0.00005), it is sufficient to process only the repeat drape lines to obtain a reliable estimate of k3, and there is no necessity to use the calibration lines. It should be noted here that the standard approach, which is used in the GTGRAV software and based on modeling gravity in the time domain, provides a reliable estimate of k3 only from processing the calibration lines.

#### **5 Discussion**

The gravity disturbance was estimated with the new approach based on spatial modeling of gravity using B-splines and with the standard approach based on modeling gravity in

the time domain. Figure 7 shows the along-track gravity disturbance estimates obtained for the filter cutoff wavelength (inverse cutoff frequency) of 100 and 70 s, which are equivalent to the spatial resolution g=2 of 3 and 2 km, respectively. From Eq. (9), the corresponding values of the length of the spline knot step x are 2 and 1.4 km. These values resulted in N D 46 and N D 64 spline coefficients in the gravity model Eq. (7), respectively.

Figure 8 shows how the estimates of the spline coefficients evolve during processing of the repeat drape lines. The estimate changes occur due to the measurement update step of the Kalman filter. Since B-splines have finite support in the time domain, the estimate changes are local. As the spline **Fig. 8** Evolution of the estimates of two spline coefficients over the time at six repeat drape lines. The vertical coloured stripes indicate the time intervals of the aircraft turns

coefficients depend on position, each new repeat line slightly refines the coefficient estimates.

Mean value of the difference of the gravity estimates obtained with the two approaches (100 s cutoff wavelength) is 0.05 mGal, the STD is 1.13 mGal. It can be seen from Fig. 7 that the gravity estimates provided by the new algorithm are more detailed and better correspond to the height profile in the interval from 60 to 70 km.

In the absence of ground gravity data, we compare the gravity estimates with the global gravity model EGM2008, which is sufficiently accurate in the area of interest, as dense gravity datasets were available at the model construction stage. Table 3 summarizes the statistics for comparing the gravity estimates and the gravity values derived from EGM2008. It can be seen that both approaches (spatial modeling and modeling in the time domain) show the same result for the cutoff wavelength of 180 s. For 320 s (equivalent to the nominal spatial resolution of EGM2008, which equals 9 km), the new approach showed greater discrepancy than the standard one. The reason is likely in short-wavelength errors present in the gravity estimate provided by the new approach, as the optimal filter has less steep slope than in the case of the standard approach (see Fig. 2).

At last, the gravity estimates were computed by processing a different number of repeat drape lines (from 1 to 6) to analyze how each new repeat line affects the gravity estimate. The estimates were obtained with the two approaches for the cutoff wavelength of 100 s and compared to EGM2008. The results are shown in Fig. 9. The dependance on the number of the processed lines is different for the two approaches. For the new approach, the discrepancy with EGM2008 grows with the number of processed lines, which can be explained by the more and more detailed along-track estimate of gravity. For the standard approach, on the contrary, the along-track gravity estimate obtained by averaging over the lines becomes less detailed with the number of processed lines, and tends to gravity from EGM2008.

**Table 3** Standard deviations of the difference of the along-track gravity estimates obtained with the two approaches (spline-based modeling and Gauss-Markov modeling) and gravity values derived from EGM2008 (mGal)


**Fig. 9** Standard deviations of the difference of the gravity estimates obtained for the cutoff wavelength of 100 s and gravity from EGM200

#### **6 Conclusions**

In this research, an approach to gravity estimation at repeat drape lines was developed using B-splines for spatial modeling of gravity. This approach was compared with the standard one based on modeling gravity in the time domain by processing the GT-2A gravimeter data. The spatial gravity modeling showed an increase in accuracy of estimating the gravity sensor scale factor significantly. The achieved accuracy was sufficient for gravity estimation along the repeat drape lines with the accuracy of 0.5 mGal. Moreover, the calibration lines, which are necessary for accurate estimation of the scale factor when using the standard approach, are not required in the case of the developed approach. The reliable estimate of the scale factor was obtained from processing the repeat drape lines.

The gravity estimate provided by the new approach is in good agreement with EGM2008 (for the filter cutoff wavelength close to the minimum wavelength represented by the global model). For the shorter cutoff wavelengths, the new approach provides more detailed gravity estimates than the standard one.

**Acknowledgements** This work was performed with financial support from the Russian Foundation for Basic Research (grant no. 19-01- 00179). We are very thankful to two anonymous reviewers for their valuable comments.

#### **References**

Berzhitzky VN, Bolotin YuV, Golovan AA, Ilyin VN, Parusnikov NA et al (2002) GT-1A inertial gravimeter system: results of flight tests. Lomonosov Moscow State University, Faculty of Mechanics and Mathematics, Moscow, Russia


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons. org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Gravimetric Studies in the Sea of Japan**

Maksim Georgievich Valitov, Ruslan Grigoryevich Kulinich, Zoya Nikolaevna Proshkina, and Tatyana Nikolaevna Kolpashchikova

#### **Abstract**

Marine gravimetric studies were carry out in the waters of the Sea of Japan and the Tatar Strait. Maps of the gravitational field of the earth were completed, structural density modeling was executed, the depth of the Mohoroviciˇ c discontinuity was determined. ´

**Keywords**

Earth's crust - Marine gravimetric studies - Mohoroviciˇ c discontinuity ´ - Sea of Japan - Structural density modeling

#### **1 Introduction**

The Sea of Japan, washing the southeastern coast of Russia, cuts almost across the cut the southern tip of the Sikhote-Alin and Laoelin-Grodekov fold systems. Here, at a short distance, a radical restructuring of the earth's crust takes place: the transition from a mature continent to a young oceanic crust with the disappearance or significant processing of the upper sialic shell rich in ore and non-metallic minerals. According to modern concepts, this, like the formation of the Sea of Japan as a whole, is the result of Meso-Cenozoic destruction of the margin of the Asian continent and rifting, which was replaced by spreading with the formation of a young oceanic crust in the eastern part of the deepsea Japanese Sea basin. As a result, a region was formed within which two radically different types of the earth's crust closely "coexist".

The deep structure of the earth's crust of this region is still an actual object of study. The decisive role in this is assigned to geophysical methods, in particular, seismic sounding and gravimetry. The capabilities of gravimetry to study the problem under consideration are based on the fact that one of the parameters distinguishing the continental crust from the oceanic crust is their significantly different thickness (the depth of the Mohoroviciˇ c discontinuity or Moho (M)). This ´ parameter with acceptable reliability is determined using gravimetry in combination with the reference data of seismic sensing, which makes it possible to determine the relief of the base of the earth's crust in a given area. On the other hand, with a "fixed" structural framework, density modeling of the geological environment is possible, which makes it possible to study changes in the material composition of the crust masses. The combination of both situations makes it possible to realize structural-density (structural-material) modeling and, thus, to study the deep structure of the earth's crust at the junction of its heterogeneous types.

#### **2 Marine Gravimetric Studies**

The first measurements of gravity in the Sea of Japan were made by Soviet researchers in 1937. Professor of Moscow State University L.V. Sorokin in a submarine made measurements of gravity at 74 points (Stroev et al. 2007).

Until the mid-1980s of the last century, the main contribution to the study of the anomalous gravitational field of the Sea of Japan was made by the P.K. Sternberg's State Astronomical Institute, who is part of Moscow State University (GAISH Moscow State University), in particular, his staff V.L. Panteleev, A.D. Gainanov, P.A. Stroyev and others. Unfortunately, the measurements were performed on

M.G. Valitov (-) · R.G. Kulinich · Z.N. Proshkina · T.N. Kolpashchikova

V.I. Il'ichev Pacific Oceanological Institute Far Eastern Branch Russian Academy of Science, Vladivostok, Russia e-mail: valitov@poi.dvo.ru

J. T. Freymueller, L. Sánchez (eds.), *5th Symposium on Terrestrial Gravimetry: Static and Mobile Measurements (TG-SMM 2019)*, International Association of Geodesy Symposia 153, https://doi.org/10.1007/1345\_2020\_128

disconnected polygons and profiles. Profile measurements were taken over the Pervenets Rise, over the southern part of the Bogorov Ridge, detailed areal observations within the Yamato Rise, at three polygons: near the Korea peninsula, and near the mainland slope of northern Primorye (Belousov et al. 1973; Stroev and Kovylin 1973).

At the same time, the entire eastern half of the Sea of Japan was covered by gravimetric surveys by the Hydrographic Service of the Ministry of Defense of Japan and Japanese universities (Tomoda 1973; Tomoda and Fujimoto 1981).

Since the mid-1980s, the V.I. Il'ichev Pacific Oceanological Institute Far Eastern Branch Russian Academy of Science (POI FEB RAS) has concentrated geophysical work in the northwestern part of the Sea of Japan (Fig. 1), within the economic zone of the Russia (formerly USSR). In addition to onboard gravimetry, seismic profiling and magnetometry were included in the set of these works (Prokudin et al. 2018).

The first measurements of gravity at this stage (1987– 1989), simultaneously with continuous seismic profiling (NSP) and hydro-magnetic surveying, were performed using on-board gravimeters GAG-ZhZ and GMN-K with analog recording of recordings on recorders. At that time, the navigational reference was very weak, positioning in the open sea was carried out no more than once in 4 h, in coastal areas, data from coastal radio navigation stations were used. The coordinates of the observation points were calculated by interpolation and navigational calculations, location errors introduced serious errors in the final accuracy of the survey. The measurement error was ˙ (2.4–3.0 mGal).

In 1990 (the 6th cruise of the R/V "Professor Gagarinsky"), survey was carried out on the border of the economic zones of Russia and North Korea, as well as on the Yamato

**Fig. 1** Scheme of the gravimetric study of the Sea of Japan. OP54–54 cruise R/V "Academician Oparin", LV81–81 cruise R/V "Academician M. Lavrentiev", OP55–55 cruise R/V Academician Oparin, LV85–85 cruise R/V "Academic M. Lavrentiev", GA – cruises made on vessels of the Far Eastern Branch of the Russian Academy of Sciences in the period from 1982 to 2010

Rise. The survey was carried out by six GMN-K instruments with analog registration (error ˙ 2.4 mGal). The direction of the profiles is meridional on the Yamato Rise and northwest in the coastal part of the Sea of Japan, the distance between the profiles is 10 miles.

The survey was continued in 1991 and 1994 in the Central basin of the Sea of Japan (cruise No. 11 and 14 of the R/V "Professor Gagarinsky"). The observations were carried out by five GMN-K instruments with analog and then with digital recording. The error of observations in 1991 was ˙2.1 mGal, and in 1994 ˙ 1.6 mGal. The direction of the profiles is meridional, the distance between the profiles on the 11th cruise was 10 miles, and on the 14th cruise – 5 miles. Characteristically, since 1994, GPS receivers began to be used for navigational support for marine surveys of the POI FEB RAS, which significantly increased the accuracy of the survey.

The Pervenets Rise is studied in detail in the 18th cruise of the R/V "Professor Gagarinsky" (1996). The survey was carried out by six GMN-K instruments with digital recording (error ˙ 1.4 mGal). The direction of the profiles is northeast, the distance between the profiles is 5 miles.

In 1996 (the 27th cruise of the R/V "Akademik M.A. Lavrentyev"), a single profile survey was taken of the northern end of the Central Basin and the Bogorov Ridge. The survey was carried out by six GMN-K instruments with digital recording (error ˙ 1.4 mGal).

The study of the northern part of the Central Basin continued on the 21st cruise of the R/V "Professor Gagarinsky" (1997). The survey was carried out by five GMN-K instruments with digital registration (error ˙ 2.0 mGal). The direction of the profiles is northwest, the distance between the profiles is 5 miles.

In 1998, at the R/V "Professor Gagarinsky" (24th cruise), a regional gravimetric survey was carried out in the Central Basin using a series of meridional and latitudinal profiles. The distance between the profiles was 25 km. The observation error was ˙1.6 mGal.

In 2003, POI FEB RAS (cruise 36) of the R/V Professor Gagarinsky carried out a detailed geophysical (including gravimetric) survey of Peter the Great Bay. The average distance between the profiles did not exceed 2 km. The observation error was ˙1.4 mGal. The location of the observations was determined using GPS satellite system.

In 2010, the R/V "Akademik M.A. Lavrentyev" (cruise 52) of the POI FEB RAS carried out a gravimetric survey (˙2.0 mGal error) at the northern closing of the Central Basin, in the vicinity of the Vityaz Rise.

This completed the first stage of areal geophysical work in the Central Basin of the Sea of Japan, within the economic zone of the Russian Federation. The root-mean-square error, calculated from the intersection of various gravimetric surveys, did not exceed 5 mGal. Measurement processing was carried out according to standard methods. At the first stage, straight sections of the profiles were chosen, them the vessel moved with uniform speed and constant heading. The Eötvös correction, correction for the cross-dropping effect, for the inertia of the measuring systems of gravimeters, and zero drift were introduced. From the observed values, the normal gravitational field of the Earth, calculated according to the international formula of 1967, was subtracted. The anomalous field thus obtained was analyzed for intersections, after which gravimetric data catalogs were created. The construction of maps was carried out in the SURFER Golden Software program. The main operation at the initial stage of mapping was the creation of a regular grid of gravitational field values with a square cell (grid) based on an irregular pseudo-rectangular grid of observed data. When building the grid in the specified program, the so-called "search radius" was used, i.e. the distance within which data is searched to interpolate a regular grid value. In our calculations, the optimal search radius was set from 5 to 25 km depending on the degree of exploration of the water area, as a result of which three nearby profiles participated in the formation of the interpolation value of gravity at each grid node. The map constructed in this way is shown in Fig. 2a.

Starting from 2017, the POI FEB RAS, at a new qualitative level, resumed geophysical exploration of the Sea of Japan. Marine gravimetric survey was carried out with Chekan-AM mobile gravimeters (Shelf-E modification) serial number 47. Studies were concentrated in the northern part of the Sea of Japan: the closure area of the deep-water Central Basin and the Tatar Strait. To date (October 2019), four expeditions have already been carried out in this area under the project "Integrated geological-geophysical, gasgeochemical and oceanographic studies in the Sea of Japan and the Tatar Strait."

So in 2017 at the R/V "Akademik Oparin" a geophysical survey was carried out in the southern part of the Tatar Strait from the latitude of the port of Sovetskaya Gavan in the north to the latitude of about. Moneron in the south. 112 geophysical profiles with a total length of 3,860 miles were completed here. The survey error did not exceed ˙0.49 mGal.

In 2018, two expeditions were completed. The research of the first expedition, to the R/V "Akademik M.A. Lavrentyev", was concentrated in the junction zone of the Central Basin of the Sea of Japan with the Tatar Strait trough valley. The second, at the R/V "Akademik Oparin", passed north, from the latitude of the Soviet harbor to Cape Surkum. The survey error, in the first case, did not exceed ˙0.38 mGal, in the second ˙0.4 mGal.

In May–June 2019, expeditionary studies were carried out in the remaining, not studied by us, part of the Tatar Strait, south of Moneron Island. During the expedition more than 113 tacks were completed, with a total length of more than 6,170 km. The survey error did not exceed ˙0.38 mGal.

**Fig. 2** Map of gravitational anomalies in free air (conditional level), made according to ship on-board measurements for the period of expeditionary research from 1982 to 2010 ( **a**) and from 2017 to 2019 ( **b**). The numbers indicate: 1 – Peter the Great Seamount; 2 – Bogorov Ridge; Rises: 3 – Pervenets, 4 – Vityaz, 5 – Yamato, 6 – Alpatov, 7 – Lavrentiev

#### **3 Results**

The works listed above made it possible to close white spots in the scientific geophysical study of the entire northwestern part of the Sea of Japan adjacent to the continent (Primorye, northern part of the Korea Peninsula). In recent years, the area of articulation of the Tatar Strait with the northern closure of the Sea of Japan and most of the Tatar Strait itself has been studied. This allowed, firstly, to build a series of conditional geophysical maps on a scale of 1:1,000,000 (including gravimetric) and, secondly, to study the structural and genetic relationship of the Central Japanese basin with the margins of the adjacent continent. One of the directions of such studies was the gravitational (density) modeling of the earth's crust at the junction of the indicated basin with the shelf and coastal geological structures of southern and south-western Primorye.

Modeling is one of the effective methods for studying the deep structure, physical properties and processes occurring in the bowels of the tectonosphere, based on the use of available data and ideas about the environment. The authors, in particular, in the framework of these studies quite widely used the so-called gravitational (density, structural-density) modeling. A gravitational model is usually understood as a density model, the calculated gravitational field from which, up to a constant component, is adequate to the initial field of gravity (Bryansky 1991). "By calculating models of the density distribution in the geological environment that create an effect similar to the observed gravitational field, we can obtain information about the structure of the geological environment" (Krasovsky 1981).

Structural-density modeling (Fig. 3) was carried out according to a series of profiles crossing the transition zone from the edge of the continent to the central deepsea basin of the Sea of Japan, which quite clearly showed the characteristic features of the transformation of the continental crust into the oceanic crust.

Transformation of the crust begins on the continent at a considerable distance from the upper edge of the continental slope. In different sections, this distance is not the same.

The transformation of the continental crust is both structural and substantial nature. The main factor determining this transformation is the reduction in the thickness of the "basalt" layer from the continent towards the basin of the Sea of Japan and its replacement with a mantle substrate. The supra-basaltic crust (crystalline basement) in most cases practically does not change its thickness to the continental slope. However, the substantial composition of this foundation in different parts of the studied area does not remain the same.

The final transformation of the continental crust occurs within a narrow strip of the outer shelf – the bottom of the continental slope. Here its sialic part disappears completely, the thickness of the lower crust ("basalt" substrate) is significantly reduced, which is replaced by mantle masses from below. In all the sections presented in the slope part, an increase in the thickness of the "transition" layer is observed, which indicates the processes of effusion of effusive rocks that took place here, with the formation, in some cases, clearly expressed in the relief of volcanic structures (Seamount Peter the Great).

At the base of the slope and the slope of the deep-sea basin, a narrow zone of deformations of all layers of the oceanic crust is recorded on all profiles, which can be a morphological expression of a powerful tectonic structure (tectonic seam) at the junction of dissimilar types of the earth's crust.

Another area of research was the depth's determination of the Mohoroviciˇ c discontinuity (Moho surface) bedding ´ within the Sea of Japan. It is well known that on the Mohoroviciˇ c discontinuity there is a sharp increase in the ´ speed of seismic waves and, accordingly, the density of the underlying rocks. Similar changes in density are observed during the transition from the water column to the upper part of the sedimentary layer and during the transition from the lower part of the sedimentary layer to consolidated rocks of the underlying layers (acoustic foundation). The observed gravitational field reflects the integral effect of the above boundaries and can be used to determine the depth of their occurrence by statistical methods. To determine the depth of the Mohoroviciˇ c discontinuity within the study area, we used ´ the statistical dependence of the magnitude of the averaged gravitational anomalies, sedimentary cover thickness and sea depth on the depth of the Mohoroviciˇ c discontinuity, ´ proposed in (Su Datsuan 1982):

$$H = 33.49 - 0.063 \Delta \text{g}\_{free} - 0.00482 H\_{sea} - 0.0017 h\_{sed}$$

where H is the depth to the Mohoroviciˇ c discontinuity, km; ´ 4*gfree* – anomaly in free air, mGal; *Hse<sup>Ã</sup>* – depth of the seabed, m; *hsed* – thickness of the sedimentary layer, m. The depths of the seabed (*HseÃ*) were determined according to the sonar data, *hsed* were determined according on the NSP data.

As is known, there is a statistical relationship between the thickness (the depth of Mohoroviciˇ c discontinuity) and ´ the type of the earth's crust (Belousov and Pavlenkova 1985). Moho depths in the indicated area are shown in the Figure 4. As can be seen from the figure, the relief of the base of the crust of the Central Basin of the Sea of Japan and its neast framing is of considerable complexity, which indicates the heterogeneity of different parts of this region. According to this feature, the studied area can be divided, first of all, into two large sections, the border between which runs approximately along the meridian of the underwater Pervenets Rise (132ı300). To the east of this boundary is the deepest part of the basin with a maximum protrusion of

**Fig. 3** Structural-density models of the earth' crust in the transition zone from the mainland to the Sea of Japan. 1 – water column; 2 – sedimentary layer; 3 – volcanic sedimentary layer; 4 – "basalt" layer; 5 – crystalline basement, 6 – gabbroids, 7 – pre-Mesozoic sedimentary

formations, 8 – granites, 9 – faults, 10 – a) observed gravitational field, b) calculated gravitational field, c) magnetic field. Numbers – densities (g/cm3) of rocks. The inset shows the position of the profiles. RF – Russian Federation, PRC – People's Republic of China

**Fig. 4** The map of the Moho discontinuity isodepths (kilometers below the sea level) and the types of crust beneath the north–western part Sea of Japan based on the marine gravimetric data. 1 – continental crust;

2 – reduced continental crust; 3 – oceanic crust; 4 – presumable oceanic crust with the thickened sedimentary cover (suboceanic?)

the Mohoroviciˇ c discontinuity, the depth of which from the ´ day surface varies from 14 km in the west to 12 km or less in the east. Without a water layer, the crust thickness here is 10.5–8.5 km, respectively, which is in good agreement with the results of seismic studies (Hirata et al. 1992). The Earth crust here is of the oceanic type. This area has a wedge-shaped east-north-east strike, complicated by local thickening of the crust up to 16–20 km within the underwater elevations.

The western half of the basin differs significantly from the previous one by the great depths of the Mohoroviciˇ c discon- ´ tinuity (16–18 km) and a different pattern of its relief. The central place here is occupied by a vast, laid out section of the western end of the basin of almost isometric shape, within which the depth of Moho is 15–16 km. Two "apophyses " with a relatively high position of Moho, separated by an East Korean Rise, depart from this section in the south-west and south directions. The thickness of the consolidated crust here is 13–15 km, which allows us to attribute it to the suboceanic type.

The thickness of the crust within large submarine hills reaches 26 km, which is characteristic of the Earth's crust of a subcontinental type. The transition from the deepsea basin to the continent is accompanied by the intense lowering of Moho to depths of 30 km and more, which corresponds to the minimum thickness of the continentaltype crust.

#### **4 Conclusion**

For more than half a century, Russian and foreign researchers have been studying the geological structure and natural fields of the Earth in the waters of the Sea of Japan. Scientists of the Pacific Oceanological Institute of the Far Eastern Branch of the Russian Academy of Sciences made the most complete contribution to this process, in the waters of the economic zone of Russia (formerly the USSR). Gravimetric studies, as one of the deepest geophysical methods, made it possible to determine the thickness and types of the earth's crust, study its internal substantial-density structure, separate blocks and trace tectonic faults. The integration of gravimetry with other geophysical methods, especially with deep seismic sounding, significantly improved the accuracy of these studies and reduce the ambiguity of the result. Thus, we can say that the first stage of the study of the gravitational field of the Central Basin of the Sea of Japan has been completed. In the future, it is necessary to focus on a detailed study of individual morphostructures, the implementation of near sea bottom gravimetric observations and the gravity field's gradients calculation, which, in combination with other geophysical methods, will provide new information on the structure and substantial composition of the earth's crust in the region.

The fundamental study and interpretation of the gravitational field in the waters of the Tatar Strait, in the northern part of the Sea of Japan, has just begun. The four expeditions carried out, unfortunately, did not allow us to characterize the gravitational field of this morphological structure with the required detail, especially in the area of transition from the Central Basin of the Sea of Japan to the southern part of the Tatar Strait. The authors hope that research in this area will be continued, including in cooperation with Japanese colleagues.

**Acknowledgments** The expeditionary work of 2017–2019 was carried out with the financial support of the Ministry of Education and Science of the Russian Federation, within the framework of the State Assignment of the POI FEB RAS: Subject 0271-2019-0002 (AAAA-A17-117030110032-3), Subject 0271-2019-0006 (AAAA-A17-117030110035-4), as well as the Priority Program of the Far East Branch of the Russian Academy of Sciences "Far East" (grant 18-1-010).

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/ licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **About Identification of Instrument Error Parameters for a Gravity Gradiometer**

A. A. Golovan, E. V. Gorushkina, and I. A. Papusha

#### **Abstract**

The article presents the description of two algorithms used for processing of the raw data of a gravity gradiometer. These algorithms are intended for estimation of some instrument errors. The first algorithm is applicable for the instrument operation in its stationary mode, the second proposes the use of a special test bench. Rotary gravity gradiometer of the accelerometric type was taken as a prototype for relevant mathematical models. Nowadays this type of gradiometer is brought to the stage of practical implementation and serial production.

#### **Keywords**

Gravity gradiometer - Calibration - Accelerometers -Instrument errors

#### **1 Introduction**

The gravity gradiometer, that is an analogue of the model developed by Bell Aerospace Laboratory (USA) and currently produced by Lockheed Martin (USA), is considered (Dransfield et al. 2010; Murphy 2010; Hammond and Murphy 2003).

The instrument design is based on the slowly rotating disk with diameter equals 0.2 m and angular rate 0.25 Hz. The parallel pairs of accelerometers are set on the diametrically opposite sides of the mentioned rotating disk and at the same distance from the center of the disk. The accelerometer sensitivity axes are tangent-oriented (Fig. 1). Further, for certainty, we assume that disk has the vertical

A. A. Golovan · E. V. Gorushkina

Lomonosov Moscow State University, Moscow, Russia e-mail: aagolovan@yandex.ru

I. A. Papusha (-) Lomonosov Moscow State University, Moscow, Russia

Trapeznikov Institute of Control Sciences, Russian Academy of Sciences, Moscow, Russia e-mail: ipapusha@yandex.ru

axis of rotation. This construction of a gradiometer is called GGI (gravity gradiometer instrument).

Gradiometer designed for airborne gravimetry. Gradiometer should be installed on a gyro stabilized platform.

The linear combination of accelerometer readings generates the GGI output signal. The rotation of the disk induces forced harmonic oscillations, and as a result the output signal becomes modulated at a frequency equal to twice angular rate:

$$\begin{aligned} W &= (f\_1 + f\_2) - (f\_3 + f\_4) = \\ &= l[(\varGamma\_{22} - \varGamma\_{11} + \varPi\_2^2 - \varPi\_1^2)\sin 2\varDelta t + \\ &+ (\varGamma\_{12} + \varPi\_1 \varPi\_2)\cos 2\varDelta t]. \end{aligned} \qquad (1)$$

Here

– W is the GGI output;


© The Author(s) 2021

J. T. Freymueller, L. Sánchez (eds.), *5th Symposium on Terrestrial Gravimetry: Static and Mobile Measurements (TG-SMM 2019)*, International Association of Geodesy Symposia 153, https://doi.org/10.1007/1345\_2021\_130

**Fig. 1** GGI model

The problem of estimation of tensor component using W measurements solved in post-processing. In this case one have to use current navigation parameters of high accuracy, in particular !.

The basic method used for increasing the accuracy of the GGI operation is the compensation of its instrument errors. The most stable parameters of these errors can be estimated during calibration that proposes the usage of a specialized test bench. There were discussed several options of the mentioned calibration procedure in recent publications. The first variant proposes the usage of a centrifuge (Yu abd Cai 2018). The second variant is based on the special motion of the disturbing mass (Deng et al. 2018).

The pre-start calibration mode of GGI is also possible, that consists of the regular operation at a fixed point for a short time.

It should be noted that the gradiometer error model includes the combinations of the instrumental errors of the paired accelerometers. So, it is important to do estimation of these combinations exclusively, and not to do estimation of errors of single accelerometers. This is why the calibration procedure should be performed for the assembled GGI.

#### **2 Model of Instrumental Errors of the Gradiometer**

The accepted model of instrumental errors is as follows. Let f be the measured value, f <sup>0</sup> be the result of the measurement, then f <sup>0</sup> <sup>D</sup> <sup>f</sup> <sup>C</sup> f , where f is the instrumental error. We will use the following designations (Golovan et al. 2018):


The error model of the i-th accelerometer is as follows

$$f\_i' = (1 + k\_i)f\_i + \Delta f\_{i0} + \Delta f\_{is} + \chi\_i f\_{i\perp} \tag{2}$$

where the component ifi? reflects the cross-coupling effect of the specific force acting on the proof mass of the accelerometer orthogonally to its longitudinal axis due to a skew in the disk's plane.

Let W <sup>0</sup> <sup>D</sup> .f <sup>0</sup> 1 <sup>C</sup> <sup>f</sup> <sup>0</sup> 2 / .f <sup>0</sup> 3 <sup>C</sup> <sup>f</sup> <sup>0</sup> 4 / be the measurement of the W (1). Measurement error W <sup>D</sup> W <sup>0</sup> <sup>W</sup> contains the frequency components 0; 1; 2:

$$
\Delta W = \Delta\_0 + \Delta\_1 + \Delta\_2. \tag{3}
$$

Taking into account the mentioned definitions, the error model of the gradiometer output takes the form (Golovan et al. 2018):

$$
\Delta\_0 = \left(\Delta f\_{10} + \Delta f\_{20} - \Delta f\_{30} - \Delta f\_{40}\right) +
$$

$$
+ \left(\Delta f\_{1s} + \Delta f\_{2s} - \Delta f\_{3s} - \Delta f\_{4s}\right) - \tag{4}
$$

$$
$$

$$
+ 2\frac{\Delta l\_1 + \Delta l\_2 - \Delta l\_3 - \Delta l\_4}{l}\Big|,
$$

$$
\Delta\_1 = \left[(k\_1 - k\_2 + \chi\_3 - \chi\_4)w\_2 + \right.
$$

$$
+ (k\_3 - k\_4 - \chi\_1 + \chi\_2)w\_1\right] \cos \Omega t - \tag{5}
$$

$$
$$

$$
$$

where *<sup>w</sup>*1,*w*2 are the horizontal components of a gyro platform's specific force, !<sup>P</sup> 3 is the vertical component of the angular acceleration,

$$\begin{aligned} \Delta\_2 &= 2l \left[ \delta\_1 (\varGamma\_{23} + \varalpha\_2 \omega\_3) + \delta\_2 (\varGamma\_{13} + \varalpha\_1 \omega\_3) - \\ &- 2\delta\_3 (\varGamma\_{12} + \varalpha\_1 \omega\_2) \right] \sin 2\varOmega t + \\ &+ 2l \left[ 2\delta\_3 (\varGamma\_{22} - \varGamma\_{11} + \varalpha\_2^2 - \varalpha\_1^2) + \delta\_1 (\varGamma\_{13} + \varalpha\_1 \omega\_3) - \right. \\ &- \delta\_2 (\varGamma\_{23} + \varalpha\_2 \omega\_3) \Big] \cos 2\varOmega t + \\ &+ \left( k\_1 + k\_2 - k\_3 - k\_4 + \varchi\_1 + \varchi\_2 - \varchi\_3 - \chi\_4 + \\ &+ 2 \frac{\Delta l\_1 + \varDelta l\_2 - \varDelta l\_3 - \varDelta l\_4}{l} \right) \times \end{aligned} (6)$$

$$\times \left[ \frac{l}{4} (\varGamma\_{22} - \varGamma\_{11} + \varPi\_2^2 - \varPi\_1^2) \sin 2\varOmega t + 1 \right.$$

$$+ \frac{l}{2} (\varGamma\_{12} + \varPi\_1 \varPi\_2) \cos 2\varOmega t \Big].$$

#### **3 Calibration Algorithm**

The following combinations of the GGI instrument error are of interest:

$$\begin{aligned} \gamma\_1 &= \delta\_1, & \gamma\_2 &= \delta\_2, & \gamma\_3 &= \delta\_3, \\ \gamma\_4 &= \Delta f\_{10} + \Delta f\_{20} - \Delta f\_{30} - \Delta f\_{40}, \\ \gamma\_5 &= k\_1 + k\_2 - k\_3 - k\_4 + \chi\_1 + \chi\_2 - \chi\_3 - \chi\_4 + \\ & + 2\frac{\Delta l\_1 + \Delta l\_2 - \Delta l\_3 - \Delta l\_4}{l}. \end{aligned} \quad (7)$$

It is obvious that the constructive instrument errors – i (that are the mutual misalignments of the sensitive axes of accelerometers) and li (that is the inaccuracy of disk's radius value) remain unchanged in GGI operation, so it is necessary to include them in the relevant calibration problem. At the same time ıi (that is the GGI installation error), fi 0 and ki (that are the biases and scale factor errors) may vary from run to run. So these parameters should be estimated during GGI operation.

#### **3.1 Case of Fixed Point**

The following reference information is used for calibration of the gradiometer at a fixed point:

– the vector of earth rotation with respect to the geographical reference frame

$$\begin{aligned} \omega\_1 &= 0 \\ \omega\_2 &= \mu \cos \varphi \\ \omega\_3 &= \mu \sin \varphi \end{aligned} \tag{8}$$



$$z = W' - W = \Delta W = \Delta[(f\_1 + f\_2) - (f\_3 + f\_4)] = 0$$

$$= \chi\_4 + (\Delta f\_{1s} + \Delta f\_{2s} - \Delta f\_{3s} - \Delta f\_{4s}) - $$

$$\begin{aligned} &+2l\left[\chi\_1(\varGamma\_{23}+u^2\sin\varphi\cos\varphi)+\varvarphi\_2\varGamma\_{13}-2\varvarphi\_3\varGamma\_{12}+\\ &+\frac{l}{4}\chi\_5(\varGamma\_{22}-\varGamma\_{11}+u^2\cos^2\varphi)\right]\sin 2\varOmega t+\\ &+2l\left[2\varvarphi\_3(\varGamma\_{22}-\varGamma\_{11}+u^2\cos^2\varphi)+\varvarphi\_1\varGamma\_{13}-\\ &-\chi\_2(\varGamma\_{23}+u^2\sin\varphi\cos\varphi)+\frac{l}{2}\varvarphi\_5\varGamma\_{12}\right]\cos 2\varOmega t.\end{aligned} (9)$$

Denote by x1, x2, x3 the amplitudes of the harmonics in the *z* measurement:

$$x\_1 = \gamma\_4 \tag{10}$$

$$\chi\_2 = 2l \left[ \chi\_1 (\varGamma\_{23} + \mu^2 \sin \varphi \cos \varphi) + \varvarphi\_2 \varGamma\_{13} - 2\varvarphi\_3 \varGamma\_{12} \right] + $$

$$+ \frac{l}{4} \chi\_5 (\varGamma\_{22} - \varGamma\_{11} + \mu^2 \cos^2 \varphi) \tag{11}$$

$$\tau$$

$$\chi\_3 = 2l \left[ 2\chi\_3(\varGamma\_{22} - \varGamma\_{11} + \varmu^2 \cos^2 \varphi) + \vargamma\_1 \varGamma\_{13} - \varGamma\_{23} \right]$$

$$-\vargamma\_2(\varGamma\_{23} + \varmu^2 \sin \varphi \cos \varphi) \Big] + \frac{l}{2} \vargamma\_5 \varGamma\_{12}.\tag{12}$$

Let assume that instrumental errors are constant during the experiment:

$$
\Delta \dot{f\_{i0}} = 0, \ i = 1, \ldots, 4$$

$$
\dot{\delta}\_i = 0, \ i = 1, \ldots, 3$$

$$
\dot{k}\_i = 0, \ i = 1, \ldots, 4 \qquad \qquad \qquad (13)$$

$$
\dot{\chi}\_i = 0, \ i = 1, \ldots, 4$$

$$
\Delta \dot{l}\_i = 0, \ i = 1, \ldots, 4.$$

Then the estimation problem for the state vector x <sup>D</sup> .x1; x2; x3/<sup>T</sup> of linear dynamic system can be set as

$$\begin{aligned} \dot{\mathbf{x}} &= \mathbf{0}, \\ \mathbf{z} &= \mathbf{x}\_1 + \mathbf{x}\_2 \sin 2\Omega t + \mathbf{x}\_3 \cos 2\Omega t + r, \qquad (14) \end{aligned}$$

where r is the white noise of a given intensity that is composed by the linear combination of noises of single accelerometers

$$r = \Delta f\_{1s} + \Delta f\_{2s} - \Delta f\_{3s} - \Delta f\_{4s} \tag{15}$$

The problem posed can be solved by application of the relevant Kalman filter.

Since the set of functions 1;sin 2˝t; cos 2˝t are linearly independent, the constant factors in these functions (14) become observable.

So, the estimates of combinations (7) can be obtained and above algorithm can be used in as a pre-start calibration algorithm.

#### **3.2 The Usage of a Stand of Linear Movement**

The estimates for combinations k1 <sup>C</sup>k2 k3k4 <sup>C</sup>1 <sup>C</sup>2 3 4 can be derived using test bench of linear movement. Calibration procedure proposes that the gradiometer moves in the horizontal plane with known linear accelerations along a given direction. The algorithm does not require reference values of the gravitational gradient tensor components in this case.

One can form measurements in the following way

*<sup>z</sup>*1 <sup>D</sup> <sup>f</sup> <sup>0</sup> 1 <sup>f</sup> <sup>0</sup> 2 <sup>D</sup> f10 f20 <sup>C</sup> f1s f2s Œ.k1 <sup>C</sup> k2/*w*1 <sup>C</sup> .1 <sup>C</sup> 2/*w*2sin ˝t<sup>C</sup> <sup>C</sup> Œ.k1 <sup>C</sup> k2/*w*2 <sup>C</sup> .1 <sup>C</sup> 2/*w*1 cos ˝t (16)

*<sup>z</sup>*2 <sup>D</sup> <sup>f</sup> <sup>0</sup> 3 <sup>f</sup> <sup>0</sup> 4 <sup>D</sup> f30 f40 <sup>C</sup> f3s f4s Œ.k3 <sup>C</sup> k4/*w*2 .3 <sup>C</sup> 4/*w*1sin ˝t Œ.k3 <sup>C</sup> k4/*w*1 <sup>C</sup> .3 <sup>C</sup> 4/*w*2 cos ˝t: (17)

When one selects the orientation of the stands rail in the North-South direction (the linear acceleration component *<sup>w</sup>*1 <sup>D</sup> <sup>0</sup>), then the measurements *<sup>z</sup>*1;*z*2 take the form:

$$\begin{aligned} z\_1 &= \Delta f\_{10} - \Delta f\_{20} + \Delta f\_{1s} - \Delta f\_{2s} - \\ &- (\chi\_1 + \chi\_2) w\_2 \sin \Omega t + (k\_1 + k\_2) w\_2 \cos \Omega t, \\ z\_2 &= \Delta f\_{30} - \Delta f\_{40} + \Delta f\_{3s} - \Delta f\_{4s} - \\ &- (k\_3 + k\_4) w\_2 \sin \Omega t - (\chi\_3 + \chi\_4) w\_2 \cos \Omega t. \end{aligned} \tag{18}$$

Let introduce:

$$\begin{aligned} \mathbf{x}\_1 &= \Delta f\_{10} - \Delta f\_{20} \\ \mathbf{x}\_2 &= k\_1 + k\_2 \\ \mathbf{x}\_3 &= \chi\_1 + \chi\_2 \end{aligned} \qquad \begin{aligned} \mathbf{x}\_3 &= \Delta f\_{30} - \Delta f\_{40} \\ \mathbf{x}\_4 &= k\_3 + k\_4 \\ \mathbf{x}\_5 &= \chi\_3 + \chi\_4. \end{aligned} \quad (19)$$

Then

$$\begin{aligned} z\_1 &= x\_1 + x\_2 w\_2 \cos \Omega t - x\_3 w\_2 \sin \Omega t + r\_1 \\ z\_2 &= x\_4 - x\_5 w\_2 \sin \Omega t - x\_6 w\_2 \cos \Omega t + r\_2. \end{aligned} \quad (20)$$

Thus, the problem is reduced to two independent estimation problems

$$\begin{aligned} \dot{\dot{x}}\_1 &= 0 & \dot{x}\_4 &= 0\\ \dot{x}\_2 &= 0 & \text{and} & \dot{x}\_5 &= 0\\ \dot{x}\_3 &= 0 & \dot{x}\_6 &= 0 \end{aligned} \tag{21}$$

with the help of aiding measurement *<sup>z</sup>*1;*z*2.

Since the functions 1;sin ˝t; cos ˝t are linearly independent, then the constant factors in these functions are observed. Therefore, the following combination k1Ck2; k3<sup>C</sup> k4; 1 <sup>C</sup> 2; 3 <sup>C</sup> 4 can be estimated, e.g. applying Kalman Filter.

The similar algorithm can be implemented for the estimation of k1 <sup>C</sup> k2 and k3 <sup>C</sup> k4 combinations in GGI regular operation mode (O'Keefe et al. 1997).

The 1; 2; 3; 4; 5 components (7) in the GGI error model (4, 6) can be compensated on the basis of corresponding estimates derived by Kalman filters.

#### **4 Conclusion**

The algorithms were elaborated for the estimation of instrument errors of GGI. The implementation of these algorithms is fairly easy, they can be included into calibration procedure.

**Acknowledgements** This research was supported by RFBR (project number 19-01-00179).

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons. org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Numerical Model of Moving-Base Rotating Accelerometer Gravity Gradiometer**

Mingbiao Yu and Tijing Cai

#### **Abstract**

In development of a rotating accelerometer gravity gradiometer (RAGG), it is difficult for us to distinguish the measurement signals and error components in the RAGG output without a verified and correctness RAGG analytical model. In addition, many key techniques, such as RAGG analytical model validation, online error compensation, post motion error compensation, are difficult to be verified. RAGG numerical model can provide validation platform for solving all the above problems, which can speed up the development of the RAGG. In this study, based on the principle and configuration of the RAGG, we synthetically consider almost all the error factors, such as circuit gain mismatch, installation errors, accelerometer scale-factor imbalance, and accelerometer second-order error coefficients, construct a parameters adjustable RAGG numerical model. In multi-frequency gravitational gradient simulation experiment, we use the RAGG numerical model simulating the situation that a test mass rotates about the RAGG with time-varying angular velocity to generate multi-frequency gravitational gradient excitations; the experiment results are consistent with the theoretical ones; the RAGG numerical model can recur some phenomenons of a actual RAGG.

#### **Keywords**

Rotating accelerometer gravity gradiometer - Numerical model - Center gravitational gradients -Virtual rotating accelerometer gravity gradiometer

#### **1 Introduction**

Airborne gravity gradiometry is an advanced technology for surveying a gravity field; it acquires gravity field information with high efficiency and high spatial resolution. Compared with gravity information, the gravity gradient tensor provides more information on the field source such as orientation, depth, and shape (Tang and Hu 2018; Yan and Ma 2015). There are many different types of gravity gradiometer: rotating accelerometer gravity gradiometers, superconducting gravity gradiometers, cold atomic interferometer gravity gradiometers (Paik 2007; Difrancesco 2007; Moody 2011; Liu 2014; Hao 2013). Rotating accelerometer gravity gradiometer (RAGG) was developed by Ernest Metzger of Bell-Aerospace in the 1980s (Heard 1988). But, over the past decade, rotating accelerometer gravity gradiometers are the only commercial moving base gravity gradiometer in the word; all other types are either in fight testing or in a laboratory setting. Companies that operate commercial rotating accelerometer gravity gradiometer systems are: Bell Geospace (Air-FTG), ARKeX (FTGeX), GEDEX (HD-AGG), and FUGRO (AGG-Falcon) (Rogers 2009; Metzger 1977). In this study, we synthetically consider almost all the unperfect factors, for example accelerometer installation errors, accelerometer scale-factor imbalance, etc., and construct a numerical model. The RAGG numerical model

M. Yu · T. Cai (-)

Instrument Science and Engineering College, Southeast University, Nanjing, China e-mail: caitij@seu.edu.cn

J. T. Freymueller, L. Sánchez (eds.), *5th Symposium on Terrestrial Gravimetry: Static and Mobile Measurements (TG-SMM 2019)*, International Association of Geodesy Symposia 153, https://doi.org/10.1007/1345\_2020\_114

is a virtual RAGG with a comprehensive set of precisely adjustable parameters; based on it, many key techniques, for example RAGG automatic calibration, post error compensation, and self-gradient modeling, can be verified.

#### **2 RAGG Numerical Model**

A high-precision numerical model of the RAGG is established, as shown in Fig. 1. In the numerical model, each accelerometer has six mounting error parameters: radial distance, initial phase angle, altitude angle, and misalignment error angles. Among them, the radial distance, initial phase angle, and altitude angle determine the mounting position of the accelerometer; the misalignment error angles determine the orientation deviation between the accelerometer measurement frame and the accelerometer nominal frame of the actual mounting position. The detailed definition of the mounting error parameters can be found in literature (Yu and Cai 2019). Moreover, each accelerometer has nine other output model parameters: zero bias (Kj 0), linear scale factors (Kj 1), second-order error coefficients (Kj 2, Kj 4, Kj 5, Kj 6, Kj 7, Kj 8), and current to voltage gain (kjV =I ). We use a test mass to produce gravitational gradients to excite the RAGG. The specific force of the accelerometer Aj in the RAGG numerical model is given by:

$$\begin{cases} \boldsymbol{f}\_{j} = \boldsymbol{f}\_{cmm} + \dot{\boldsymbol{\omega}}\_{im} \times \boldsymbol{r}\_{o\_{m}A\_{j}} + \boldsymbol{\omega}\_{im} \times \left(\boldsymbol{\omega}\_{im} \times \boldsymbol{r}\_{o\_{m}A\_{j}}\right) \\ -Gm\boldsymbol{A}\_{j}\boldsymbol{\mathbf{S}} \sqrt{\left|\boldsymbol{A}\_{j}\boldsymbol{\mathbf{S}}\right|^{3}} \,. \end{cases} (1)$$

Where f cmm is the specific force of the RAGG; !P im is the angular acceleration of the RAGG with respect to the inertial frame; !im is the angular velocity of the RAGG with respect to the inertial frame; G is the gravitational constant; romAj is the position vector of accelerometer Aj in the RAGG measurement frame; m represents the weight of the test mass; and AjS is the position vector from accelerometer Aj to the test mass. If the test mass is not a point mass, the gravitational acceleration that the RAGG accelerometers undergo produced by the test mass can be calculated using finite element analysis. In addition, the test mass can be in motion with respect to the RAGG; in this case, AjS is time varying. The specific forces of accelerometer Aj in the accelerometer nominal frame of the actual mounting position (fjx, fjy, fj *<sup>z</sup>*) can be calculated from:

$$\begin{array}{l} f\_{jx} = \mathbf{f}\_{\,j} \cdot \mathbf{r}\_{jx} \,, \\ f\_{jy} = \mathbf{f}\_{\,j} \cdot \mathbf{r}\_{jy} \,, \\ f\_{jz} = \mathbf{f}\_{\,j} \cdot \mathbf{r}\_{jz} \,. \end{array} \tag{2}$$

**Fig. 1** Principle and program flow of the RAGG numerical model. (**a**) Principle of the RAGG numerical model. (**b**) Program flow of the RAGG numerical model

Where jx, jy, and <sup>j</sup> *<sup>z</sup>* are unit vectors of the accelerometer nominal frame of the actual mounting position in the directions of the x-, y-, and z-axes. The specific forces in the accelerometer measurement frame are:

$$
\begin{bmatrix} f\_{ji} \\ f\_{jo} \\ f\_{jp} \end{bmatrix} = \mathbf{C} \begin{bmatrix} f\_{jx} \\ f\_{jy} \\ f\_{jz} \end{bmatrix} \,, \tag{3}
$$

where C is the transformation matrix from the accelerometer nominal frame of the actual mounting position to the accelerometer measurement frame. To make the numerical model approximate the actual RAGG, we add accelerometer noise to the accelerometer model:

$$\begin{array}{c} \frac{V\_i}{k\_{jV/I}K\_{j1}} = f\_{jnoise} + f\_{ji} + K\_{j0} + K\_{j2}f\_{ji}^{\ \ 2} + K\_{j5}f\_{jo}^{\ \ 2} \\ + K\_{j7}f\_{jp}^{\ \ 2} + K\_{j6}f\_{ji}f\_{jo} + K\_{j4}f\_{ji}f\_{jp} + K\_{j8}f\_{jo}f\_{jp} \ . \end{array} \tag{4}$$

The accelerometer noise fj noise is simulated by a power spectral density model (Jekeli 2006):

$$
\Phi(f)\_{noise} = \alpha f^{-b} + \alpha\_T \,, \tag{5}
$$

where ˛ and b represent the amplitude and low-frequency growth of the red noise, and !T denotes the amplitude of the white noise.

Figure 1b is the program flow of the RAGG numerical model. Firstly, the RAGG simulation parameters are set up, including test masses parameters, RAGG rotating disk parameters, accelerometer mounting parameters, accelerometer model parameters, RAGG motion parameters, etc. Then substituting the parameters into the formula (1)–(3) calculates the specific force in accelerometer measurement frame at time t. According to the formula (4), calculating the output voltage of the RAGG accelerometer, the RAGG output before demodulation at time t is calculated by: Go*u*t.t / D V1.t / C V2.t / V3.t / V4.t /. The above process is repeated until time t is equal to the simulation duration time. Finally, the RAGG output data is input to the quadrature amplitude modulation (QAM) demodulator to extract gravitational gradient.

Let xx, xy, <sup>x</sup>*<sup>z</sup>*, yy, and <sup>y</sup>*<sup>z</sup>* represent the five independent gravitational gradient elements at the origin of the RAGG measurement frame. When mass is far enough away from the RAGG, the gravitational acceleration measured by the RAGG accelerometers is a first-order approximation of the gravitational acceleration and gravitational gradient tensor at the center of the rotating disc; in this case, the inline channel measurement and the cross channel measurement of the RAGG approximate xx xy and xy; otherwise, the inline channel measurement and cross channel measurement of the RAGG are the sum of xx xy, xy, and high-order gravitational gradient tensor elements. To distinguish between xx xy, xy and the measurements of the RAGG, xx xy, xy is called as center gravitational gradients; xx xy is the inline channel of the center gravitational gradients; xy is the cross channel of the center gravitational gradients. In RAGG analytical model, gravitational acceleration that the RAGG accelerometer undergoes is a first-order taylor approximation of the gravitational acceleration and gravitational gradient tensor at the center of the rotating disc; but in the numerical model, the gravitational accelerations are calculated using Newton's law of gravitation instead of a linear approximation. Therefore the numerical model are close to the actual RAGG.

#### **3 Multi-frequency Gravitational Gradient Simulation Experiment**

In multi-frequency gravitational gradient simulation experiment, a test mass rotates about the RAGG with time-varying angular velocity producing multi-frequency gravitational gradient excitations. Based on the angular velocity of the test mass and its initial coordinate in the RAGG measurement frame, we can obtain the coordinates of the test mass in the RAGG measurement frame at any time. We then calculate the gravitational gradient tensor at the origin of the RAGG measurement frame and then calculate the center gravitational gradients. The RAGG numerical model simulates a perfect RAGG without accelerometer mounting errors, accelerometer scale-factor imbalances, accelerometer second-order error coefficients and accelerometer noise, so we set the accelerometer mounting errors, the accelerometer second-order error coefficients, accelerometer noise parameters to zero. The linear scale factor of the four accelerometers is kj 1 D 10 mA/g, the currentto-voltage gain is kjV =I D 10<sup>9</sup> ohm, the nominal mounting radius R is 0.1 m, and the rotation frequency

**Fig. 2** A test mass rotating about the RAGG with time-varying angular velocity

**Fig. 3** Accelerometer voltage outputs in RAGG numerical model

**Fig. 4** Output voltage of the numerical model before demodulation

of the RAGG disc is 0.25 Hz. Figure 2 shows a point mass of 486 kg with an initial position in the RAGG measurement frame of (1:1; 0:5; 0) and rotating about the RAGG with time-varying angular speed !.t / D 500 C 400 sin.0:0628t /ı=h.

Figure 3 shows the voltage outputs of the four accelerometers in the RAGG numerical model excited by the rotating point mass. Figure 4 shows the output voltage before demodulation of the numerical model. Figures 5 and 6 show the demodulated gravitational gradient comparison among the RAGG model and the center gravitational gradients; we can see that the inline channel and the cross channel of the RAGG numerical model are overlapped with that of the center gravitational gradients and only one curve is displayed; to distinguish the outputs of numerical model and the center gravitational gradients, their differences are calculated and shown in Figs. 7 and 8; the differences are in the order of 10 -<sup>1</sup> E, and they are caused by the high-order gravitational gradient tensor.

**Fig. 5** Demodulated gravitational gradient comparison between numerical model and center gravitational gradients: Inline channel

**Fig. 6** Demodulated gravitational gradient comparison between numerical model and center gravitational gradients: Cross-channel

**Fig. 7** The difference between center gravitational gradients and numerical model: Inline channel

**Fig. 8** The difference between center gravitational gradients and numerical model: Cross-channel

#### **4 Conclusion**

Based on the measurement principle and configuration of the RAGG, we considered the factors of circuit gain mismatch, installation error, accelerometer scale-factor imbalance, and accelerometer second-order error coefficients, then developed a high-precision numerical model. The parameters of the RAGG prototypes are unknowable and uncontrollable, but the parameters of the RAGG numerical model are adjustable and knowable; to some extent, the RAGG numerical model are more suitable for verifying some technique solutions of the RAGG.

**Acknowledgements** This work is supported by National Key R&D Program of China under Grant No. 2017YFC0601601, 2016YFC0303006 and International Special Projects for Scientific and Technological Cooperation under Grant No. 2014DFR80750.

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons. org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **New Algorithm for Gravity Vector Estimation from Airborne Data Using Spherical Scaling Functions**

Vadim S. Vyazmin

#### **Abstract**

The paper presents an approach to determination of the gravity disturbance vector from airborne gravimetry measurements at an aircraft's flight path. A navigation-grade inertial navigation system (INS) and the carrier-phase differential mode of GNSS are assumed. To improve observability of the gravity horizontal components, which are observed in combination with the INS systematic errors, we use a spatial model of gravity. We parameterize the disturbing potential in the observation area using the spherical scaling functions. The unknown coefficients of the parameterization and the INS systematic errors are estimated simultaneously using the Kalman filter. Due to ill-conditioning of the estimation problem, the information form of the Kalman filter and regularization are used. The numerical results obtained from simulated data processing shows that the approach based on spatial modeling is capable to improve accuracy of the gravity horizontal component determination comparing to a typical modeling of gravity in the time domain.

#### **Keywords**

Airborne gravimetry - Gravity disturbance vector - Kalman filter - Spherical scaling functions

#### **1 Introduction**

Airborne gravimetry is capable to provide Earth's gravity observations of high accuracy and spatial resolution. A typical airborne gravimetry system includes an inertial navigation system (INS) (strapdown or a gyro-stabilized platform) and global navigation satellite system (GNSS) receivers operating in the carrier-phase differential mode. Airborne gravimetry is well established for determining the vertical component of the gravity disturbance vector (scalar gravimetry based on platformed INS) and is capable to provide all three gravity vector components simultaneously from airborne measurements (vector gravimetry) (Kwon and Jekeli 2001; Becker et al. 2016). The common method used in vector gravimetry is based on INS/GNSS integration via Kalman filtering. A high accuracy (navigation grade) strapdown INS is assumed. The main difficulty is in determining the gravity horizontal components (deflections of the vertical) as these are observed in combination with the INS systematic errors (the inertial sensor biases, the attitude errors, etc.), which have a similar spectral content. An additional information (e.g., an a priori gravity model, external gravity data) is required for separation.

The common approach to airborne vector gravimetry is based on modeling gravity in the time domain during Kalman filtering, e.g., assuming each gravity component to be a stationary stochastic process (Becker et al. 2016; Jensen et al. 2019). As an alternative, a gravity model may not be used explicitly, and the gravity components can be determined from the Kalman filter residuals (Kwon and Jekeli 2001). However, remains of the systematic errors may be present in the gravity estimates of both approaches, and additional corrections should be applied.

V. S. Vyazmin (-)

Navigation and Control Laboratory, Department of Applied Mechanics and Control, Lomonosov Moscow State University, Moscow, Russia e-mail: v.vyazmin@navlab.ru

J. T. Freymueller, L. Sánchez (eds.), *5th Symposium on Terrestrial Gravimetry: Static and Mobile Measurements (TG-SMM 2019)*, International Association of Geodesy Symposia 153, https://doi.org/10.1007/1345\_2020\_113

In this study, we develop an approach that takes into account spatial correlation of gravity at adjacent survey lines to improve observability of the gravity horizontal components. The approach is based on spatial modeling of the disturbing potential in the observation area using the spherical scaling functions (SSF) of a certain resolution level, which are harmonic and spatially localized (decrease with distance from their origins) (Freeden and Michel 2004). The SSFs and other spherical radial basis functions are widely used in local and regional gravity field modeling and are mainly applied for processing satellite gravity data and combining gravity data, e.g., Lieb et al. (2016) and Klees et al. (2008). We use the SSFs to airborne vector gravimetry problem solving. The developed estimation algorithm is based on simultaneous estimation of the INS systematic errors and the scaling coefficients (the coefficients of the parameterization of the disturbing potential) via Kalman filtering given measurements at the aircraft's flight path. The information form of the Kalman filter (Kailath et al. 2000) and regularization are applied due to the ill-conditioning of the estimation problem. The approach was partially presented in Bolotin and Vyazmin (2016).

The aim of this work is to investigate the performance of the developed approach by processing simulated airborne data and by comparing with the common approach based on modeling the gravity vector components in the time domain. The numerical results obtained with the developed approach are optimistic and show that an increase in the accuracy of the gravity horizontal component determination can be achieved.

#### **2 Airborne Vector Gravimetry Equations**

The gravity disturbance vector determination is based on INS/GNSS integration. The basic equation is the equation of motion of the INS proof mass M expressed in the navigation (local-level) frame M x as

$$
\dot{\boldsymbol{\sigma}}\_{\boldsymbol{x}} = (\hat{\boldsymbol{\sigma}}\_{\boldsymbol{\chi}\boldsymbol{\xi}}^{\boldsymbol{\chi}} + \hat{\boldsymbol{\sigma}}\_{\boldsymbol{\eta}\boldsymbol{\xi}}^{\boldsymbol{\chi}}) \boldsymbol{\nu}\_{\boldsymbol{x}} + L\_{z\boldsymbol{x}}^{T} \boldsymbol{f}\_{z} + \boldsymbol{\mathcal{y}}\_{\boldsymbol{x}} + \boldsymbol{\delta} \mathbf{g}\_{\boldsymbol{x}}.\tag{1}
$$

Here <sup>v</sup><sup>x</sup> is the <sup>3</sup> <sup>1</sup> vector of the aircraft velocity relative to the Earth expressed in M x frame coordinates, <sup>ı</sup>g<sup>x</sup> is the gravity disturbance vector, <sup>x</sup> is the reference (normal) gravity vector, <sup>f</sup> *<sup>z</sup>* is the specific force expressed in coordinates of the body frame (denoted by M*z*). Further, L*<sup>z</sup>*<sup>x</sup> is the transformation matrix from the navigation frame (M x) to the body frame which satisfies the kinematic equation (Farrell 2008), !<sup>x</sup> x is the (absolute) angular velocity of the navigation frame with respect to the Earth-centered inertial frame, expressed in M x, !<sup>x</sup> - is the Earth absolute angular velocity vector expressed in M x, !<sup>O</sup> <sup>x</sup> is a skew-symmetric matrix formed by the components of the vector !<sup>x</sup> in such a way that !O x <sup>v</sup><sup>x</sup> D !<sup>x</sup> - <sup>v</sup>x, where means the cross product.

From Eq. (1), we arrive at the INS mechanization equations (Farrell 2008) by replacing true values in Eq. (1) with measurements from the INS sensors (measurements of the body-frame angular velocity and the specific force) and with GNSS position and velocity, and omitting the gravity disturbance vector term. Here, we assume GNSS positioning and velocity solutions to be obtained using the carrierphase differential mode. In this case, inaccuracies in GNSS positions can be neglected (Vyazmin and Bolotin 2019). Denote by v<sup>0</sup> <sup>x</sup> the velocity solution obtained from solving the INS mechanization equations and by L<sup>0</sup> *<sup>z</sup>*<sup>x</sup> the transformation matrix from the navigation frame to the computed body frame.

By subtracting Eq. (1) from the INS mechanization equations, the INS error dynamics equations can be obtained as follows (Farrell 2008; Vyazmin and Bolotin 2019)

$$
\delta\dot{\boldsymbol{\sigma}}\_{\boldsymbol{x}} = (\hat{\boldsymbol{\sigma}}\_{\boldsymbol{x}\boldsymbol{\xi}}^{\boldsymbol{x}} + \hat{\boldsymbol{\sigma}}\_{\eta\boldsymbol{\xi}}^{\boldsymbol{x}})\delta\boldsymbol{\sigma}\_{\boldsymbol{x}} + \hat{\boldsymbol{\mathcal{Y}}}\_{\boldsymbol{x}}\boldsymbol{\alpha}\_{\boldsymbol{x}} - \delta\mathbf{g}\_{\boldsymbol{x}} + \tag{2}
$$

$$+L\_{\boldsymbol{x}\boldsymbol{\varepsilon}}^{\prime}\Delta\boldsymbol{f} - \hat{\boldsymbol{v}}\_{\boldsymbol{x}}^{\prime}L\_{\boldsymbol{x}\boldsymbol{\varepsilon}}^{\prime}\Delta\boldsymbol{\omega},$$

$$\dot{\boldsymbol{\alpha}}\_{\boldsymbol{x}} = \hat{\boldsymbol{\omega}}\_{\boldsymbol{x}\boldsymbol{\xi}}^{\boldsymbol{x}}\boldsymbol{\alpha}\_{\boldsymbol{x}} - L\_{\boldsymbol{x}\boldsymbol{\varepsilon}}^{\prime}\Delta\boldsymbol{\omega} \tag{3}$$

where ˛<sup>x</sup> is the attitude error, <sup>ı</sup>v<sup>x</sup> is the velocity error defined as the sum of the total velocity error v<sup>x</sup> <sup>D</sup> v0 <sup>x</sup> <sup>v</sup><sup>x</sup> and the term <sup>v</sup><sup>O</sup> <sup>x</sup>˛<sup>x</sup>. Note that the equations do not contain the INS positioning solution as we assume that the INS proof mass position is perfectly determined from GNSS.

The terms !, f are the gyroscope and accelerometer errors, respectively, for which the following models are assumed in the study:

$$
\Delta \mathfrak{o} = \mathfrak{b}\_{\mathfrak{o}} + \mathfrak{q}\_{\mathfrak{o}}, \quad \Delta f = \mathfrak{b}\_f + \mathfrak{q}\_f,\tag{4}
$$

$$
\dot{\mathbf{b}}\_{\boldsymbol{\omega}} = \mathbf{s}\_{\boldsymbol{\omega}}, \quad \dot{\mathbf{b}}\_{f} = \mathbf{s}\_{f}, \tag{5}
$$

where <sup>b</sup>!, <sup>b</sup><sup>f</sup> are the gyroscope and accelerometer biases, respectively, the vectors <sup>q</sup>!, <sup>q</sup><sup>f</sup> , "!, "<sup>f</sup> are the random errors. The error models in Eq. (5) represent the sensor bias drifts.

Denote by vGNSS <sup>x</sup> the GNSS velocity of the aircraft relative to the Earth and by <sup>y</sup><sup>x</sup> the observation computed as <sup>y</sup><sup>x</sup> <sup>D</sup> <sup>v</sup><sup>0</sup> <sup>x</sup> <sup>v</sup>GNSS <sup>x</sup> . The observation equation is written as (Vyazmin and Bolotin 2019)

$$\mathbf{y}\_{\times} = \delta \mathbf{v}\_{\times} - \hat{\mathbf{v}}\_{\times}^{\prime} \mathbf{a}\_{\times} - \mathbf{e}\_{v} \tag{6}$$

where <sup>e</sup><sup>v</sup> is the GNSS velocity error.

#### **3 Local Gravity Model and Estimation Algorithm**

In order to use a Kalman filter framework, we introduce a model of the gravity disturbance vector. We use a local spatial model by parameterizing the disturbing potential at the points at the aircraft's flight path by the SSFs (Freeden and Michel 2004; Klees et al. 2008)

$$T(\mathbf{r}) = \sum\_{i=1}^{N} a\_i \Phi\_j(\mathbf{z}\_i, \mathbf{r}) \tag{7}$$

where ai , i D 1; : : : ; N, are the unknown parameters (the scaling coefficients) to be estimated from measurements, *z*<sup>i</sup> is a knot of a regular grid on a reference sphere ˝R of radius <sup>R</sup>, r is a geocentric position vector, <sup>j</sup>rj <sup>R</sup>, j is the resolution level, which is defined by the spatial resolution of measurements. The number of the scaling coefficients N depends on the density of the regular grid and on the size of the modeling area (see Fig. 1) as the SSFs are spatially localized (decrease with distance from r).

The SSF is harmonic outside the sphere ˝R and expressed via the Legendre polynomials Pn as (Freeden and Michel 2004)

$$\Phi\_j(\mathbf{z}\_i, r) = \frac{1}{4\pi} \sum\_{n=2}^{\infty} (2n+1) \left(\frac{R}{r}\right)^{n+1} b\_{j,n} \, P\_n(\mathbf{z}\_i^{\diamondsuit} \hat{\mathbf{r}}) \quad (8)$$

where <sup>r</sup> D jrj, *<sup>z</sup>*V<sup>i</sup> , <sup>r</sup><sup>V</sup> denote the unit vectors, <sup>b</sup><sup>n</sup> <sup>j</sup> is the Legendre coefficient, which determines the spectral behaviour of the SSF. We use the Abel–Poisson SSF, which is defined by bj ;n D b<sup>n</sup> <sup>j</sup> , bj D e-1=2j . The Abel– Poisson SSF has a closed-form expression (Freeden and Michel 2004), so there is no necessity to calculate the sum in Eq. (8). Figure 2 shows the spectral content of the Abel–Poisson SSFs of different resolution levels j . A greater value of j corresponds to a larger bandwidth of the SSF in the spherical harmonic domain and therefore to a higher spatial resolution of the gravity model Eq. (7).

The gravity disturbance vector can be represented via the scaling coefficients by differentiating Eq. (7) given the formula (Torge 2001)

$$
\delta \mathbf{g}(\mathbf{r}) = \text{grad } T(\mathbf{r}),
$$

as follows

$$\delta \mathbf{g}(\mathbf{r}) = \sum\_{i=1}^{N} a\_i \operatorname{grad} \Phi\_j(\mathbf{z}\_i, \mathbf{r}). \tag{9}$$

**Fig. 1** The ground track of the aircraft's flight path and the grid of the knots of the scaling coefficients

**Fig. 2** Legendre coefficients b<sup>n</sup> <sup>j</sup> of the Abel–Poisson scaling function for different resolution levels j

Assuming that each scaling coefficient ai is constant over the time during the aircraft's flight, the following equations hold:

$$\dot{a}\_i = 0, \quad i = 1, \ldots, N. \tag{10}$$

Thus, we obtain the state-space model, which is completely described by Eqs. (2)–(6) and Eq. (10) given the gravity spatial model Eq. (9).

*The Kalman Filter in the Information Form* Let us assume that the random errors <sup>q</sup>!, <sup>q</sup><sup>f</sup> , "!, "<sup>f</sup> , <sup>e</sup><sup>v</sup> in Eqs. (4)–(6) are zero-mean white noises with known variances. Transforming the state-space model Eqs. (2)–(6), (9)–(10) to the discretetime form, we can further use a Kalman filtering framework. The system's state vector <sup>x</sup><sup>k</sup> at the time moment tk, <sup>k</sup> <sup>D</sup> 0; : : : ; K, where K is the total number of observations, includes N C 12 unknown variables

$$
\delta \boldsymbol{\sigma}\_{\boldsymbol{x}}, \boldsymbol{a}\_{\boldsymbol{x}}, \boldsymbol{b}\_{\boldsymbol{a}}, \boldsymbol{b}\_{\boldsymbol{f}}, \boldsymbol{a}\_{1}, \dots, \boldsymbol{a}\_{N}.
$$

Using Kalman filtering, an optimal (in the sense of the minimum mean squared error) estimate of the state vector is obtained.

Since the state vector <sup>x</sup><sup>k</sup> at each time moment tk includes the scaling coefficients for the entire observation area, the estimation problem is ill-conditioned. For this reason, we use the Kalman filter in the information form (Kailath et al. 2000; Park and Kailath 1995). The filter provides the information matrix Qk D P -1 <sup>k</sup> and the estimate of the information vector *u*<sup>k</sup> D P -1 <sup>k</sup> <sup>x</sup><sup>k</sup> at each time moment tk given the initial values: Q0 D 0, *u*<sup>0</sup> D 0. Here Pk denotes the covariance matrix of the state vector estimate error. The estimates of the scaling coefficients are obtained at the last time moment tK from the state vector estimate <sup>x</sup><sup>O</sup> <sup>K</sup> , which is obtained by solving the equation QKx<sup>O</sup> <sup>K</sup> D O*u*<sup>K</sup> . Here, the information vector estimate *u*O <sup>K</sup> and the information matrix QK are provided by the filter at the last time moment tK . As the information matrix QK is ill-conditioned (due to extension of the observation area to reduce edge effects as shown in Fig. 1, and due to the inverse problem of determining the disturbing potential from the gravity disturbance vector), regularization is required. Using Tikhonov regularization with a unit matrix I , we finally obtain the state vector estimate and the covariance matrix of the estimate error, respectively, as

$$
\hat{\mathfrak{x}}\_K = P\_K \hat{\mathfrak{u}}\_K, \quad P\_K = \left(\mathcal{Q}\_K + \lambda I\right)^{-1} \tag{11}
$$

where >0 is the regularization parameter.

#### **4 Results**

*Data Simulation* The developed approach was applied to processing simulated airborne data. Error-free measurements of the inertial sensors and error-free GNSS data (positions and velocities) were simulated using the software from Bogdanov and Golovan (2017). The simulated aircraft's flight path contained 4 lines flown in the south-north direction (Fig. 1) at a constant height of 1,600 m above the reference ellipsoid. The initial geodetic latitude is 0.3ı. The line spacing is about 9.5 km, the length of each line is 100 km. The ground speed at each line is constant and equals 100 m/s. The attitude angles are constant at each line (the pitch and roll angles equal 0ı). The gravity disturbance vector at the aircraft's flight path is synthesized using EGM2008 up to d/o 1,500 (equivalent to minimum wavelength of about 26 km).

The inertial sensor errors with characteristics shown in Table 1 were simulated and added to the error-free data. The gyroscope errors were obtained from measurements of the calibrated FOG gyroscopes by Optolink recorded during a static test at constant temperature. A small constant bias

**Table 1** Statistics for the simulated inertial sensor errors


of 0.001ı/h was added to the simulated gyroscope measurements.

The accelerometer measurement error was simulated as a sum of a zero-mean white noise and a bias drift modeled as a random walk process (see Table 1). The constant bias of 2 mGal was added to the horizontal accelerometer measurements. The bias of the vertical accelerometer was assumed to be determined and removed by comparing preand post-flight static measurements with a terrestrial gravity value.

The simulated INS measurement frequency is 100 Hz. The GNSS velocity error <sup>e</sup><sup>v</sup> is modeled as the time derivative of a zero-mean white noise at 20 Hz, the standard deviation (STD) of <sup>e</sup><sup>v</sup> is 2–3 cm/s.

*Data Processing* The data processing procedure consisted of four steps. First, the INS mechanization equations were solved using data from the entire flight. The gravity vector <sup>x</sup> was associated with the reference gravity field represented as the sum of the normal gravity field and a longwavelength part of the anomalous field provided by a global gravity model (EGM2008 up to d/o 100 was used). The residual gravity disturbance vector (i.e., the gravity disturbance vector with a long-wavelength part subtracted) to be estimated varied from 15 to 15 mGal for the horizontal components and from 20 to 10 mGal for the vertical component.

Second, the spatial model of the residual disturbing potential was designed using Eq. (7). The resolution level of the Abel–Poisson SSFs was chosen equal to 8. The 10 km 4.5 km latitude-longitude regular grid of the scaling coefficient knots was used (Fig. 1). The total number of the scaling coefficients N equals 154. The dimension of the state vector in the Kalman filter is 166.

Third, the scaling coefficients of the gravity model and the INS systematic error parameters were estimated using the Kalman filter in the information form and regularization Eq. (11) applied at the last time moment tK. The value of the regularization parameter was set equal to s110-<sup>14</sup>, where s1 is the maximum singular value of the information matrix QK, which was found to be suitable.

Finally, the scaling coefficient estimates obtained at the last time moment tK were used to build the along-path estimates of the three components of the residual gravity disturbance vector using Eq. (9), and the long-wavelength part of the gravity disturbance vector was readded to the along-path estimates (EGM2008 up to d/o 100). The final

estimated using the spatial gravity modeling (SSF) and the modeling in the time domain (Fourier), and the true values. The vertical colored stripes mark the aircraft's turns

**Fig. 3** Gravity disturbances

**Fig. 4** Errors in the gravity disturbances estimated using the spatial gravity modeling (SSF) and the modeling in the time domain (Fourier). The vertical colored stripes mark the aircraft's turns

along-path estimates of the gravity disturbance vector components are shown in Fig. 3. The estimate errors are shown in Fig. 4.

*Numerical Results* The statistics for the errors of the alongline gravity estimates are given in Table 2. The standard deviation varies from 0.6 to 1.9 mGal for the east component and from 0.8 to 1.3 mGal for the north component. The average STDs are 1.2 and 1.0 mGal for the east and north component, respectively. The mean values are less than 1.7 mGal for the east and 0.7 mGal for the north component. A slightly lower accuracy of the east component, in our opinion, is due to the south-north direction of the lines and the small total number of the lines. It should be noted here that in the case of a smaller number of processed lines, the estimation accuracy for the east component is significantly worse than for the north component.

**Table 2** Statistics for the errors of the gravity disturbances estimated using the spatial modeling of gravity (mGal)


The accuracy of the along-line estimates of the vertical component is better than 0.5 mGal (1), the mean value is not greater than 1.5 mGal.


**Table 3** Statistics for the errors of the gravity disturbances estimated with the use of gravity modeling in the time domain (mGal)

*Comparison with Gravity Modeling in the Time Domain* The gravity disturbance vector was also estimated using the standard approach based on modeling gravity in the time domain. Each component of the gravity disturbance vector was represented using the Fourier expansion of degree 15 (equivalent to 25 km wavelength). The total number of the unknown coefficients (amplitudes of the trigonometric functions) equaled 93. The estimation results are shown in Fig. 3 and Table 3. The accuracy of the horizontal component estimation is significantly lower than in the case of the spatial gravity modeling. Large biases and drifts can be observed in the gravity estimates. The average STDs are 3.1 and 2.8 mGal for the east and north component, respectively.

For the vertical component, the estimation accuracy (1) is 0.3 mGal that is better than in the case of the developed approach. The better accuracy is probably due to a better approximation of the band-limited gravity signal using the Fourier expansion than using the non-bandlimited SSF. The mean values of the vertical component estimate errors are approximately the same for the two approaches.

#### **5 Discussion and Conclusions**

An approach based on spatial modeling of gravity using the spherical scaling functions was developed for the airborne vector gravimetry problem solving. The approach assumes the use of a navigation-grade INS and the carrier-phase differential mode of GNSS. From the results of simulated data processing, the new approach provides significantly more accurate estimates of the gravity horizontal components in comparison with the standard approach based on modeling gravity in the time domain. The results, however, may be too optimistic as simplified models were used for simulating the GNSS errors and the inertial sensor errors (only the bias drift and noise were simulated). Further, a long-wavelength part of the gravity vector was assumed precisely known from a lowresolution global gravity model. Despite this, the developed approach, in our opinion, can be put into practice in the near future.

There are, however, a number of issues, the study of which is of interest and probably will allow to improve accuracy of the gravity estimates. Namely, it should be investigated in more detail how the accuracy of gravity estimation depends on inaccuracies in the inertial sensor calibration results (e.g., on the bias level or a long-term drift in the horizontal accelerometers). At the same time, possible ways to improve the gravity estimation accuracy should be also studied, e.g., using a greater number of parallel survey lines, cross-lines, or several overlapping flights. In addition, more careful design of the spatial gravity model (e.g., selecting the type of the SSF, regularization parameter value, etc.) may also improve the resulting estimates.

**Acknowledgements** This work was performed with financial support from the Russian Foundation for Basic Research (grant no. 19-01- 00179). The author is thankful to two anonymous reviewers for their very valuable comments and suggestions.

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons. org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Results of Astro-Measurements of the Deflection of Vertical Using the New Observation Technique**

### Murat Murzabekov, Vyacheslav Fateev, Aleksei Pruglo, and Sergei Ravdin

#### **Abstract**

The initial information for the development of high-degree models of the Earth's gravitational field (EGF) are the results of satellite and ground-based measurements. At the same time, satellite measurements carry information on the long-wave structure of the EGF. Information on the short-wave structure of the EGF can be obtained only on the basis of ground-based measurements. Having organized the determination of deflection of vertical (DOV) with a resolution of several kilometers, the local structure of the EGF can be restored with the highest possible resolution. This can be done using digital zenith camera systems (DZCS). They are automated and allow to determine the components of the DOV at the point of placement in real time. The article presents the developed measurement technique with a DZCS and the results of its tests at various geographical points in the field. The proposed technique, unlike the existing traditional technique, allows to evaluate and take into account the calibration coefficients of the DZCS in each series of observations. In addition, the new proposed technique does not impose requirements on the accuracy of rotation of the telescope around the axis in the horizontal plane and the rigidity of the base of the DZCS. The test results of the new technique showed that the standard deviation of measurements is about 0.100–0.300.

### **Keywords**

Calibration coefficients - Deflection of vertical - Digital zenith camera - Observation technique

#### **1 Introduction**

Components of the DOV can be found in several ways:

1. If the astronomical ˆ, ƒ and geodetic *B, L* coordinates of the location are known, the DOV components -*,* are calculated as (Torge 2001):

$$\xi = \Phi - B, \qquad \eta = (\Lambda - L) \cdot \cos B. \tag{1}$$

This method is implemented in the traditional technique used in existing DZCS (Albayrak et al. 2019; Hirt et al. 2010; Tian et al. 2014; Somieski 2008).

2. If the values of the components *gx*, *gy* of the gravity vector **g** are known, the components of the DOV are found by the following formulas (Brovar 1983):

$$\xi = -\arcsin\left(\frac{\mathbf{g}\_x}{|\mathbf{g}|}\right), \qquad \eta = -\arcsin\left(\frac{\mathbf{g}\_y}{|\mathbf{g}|}\right). \qquad (2)$$

A minus sign means that the DOV in the direction of increasing coordinates are considered negative.

Normalize the gravity vector **g** in the form:

$$\frac{\mathbf{g}}{|\mathbf{g}|} = \hat{\mathbf{g}} = \boldsymbol{\gamma} \cdot \hat{\mathbf{g}} = \begin{pmatrix} \hat{\mathbf{g}}\_x & \hat{\mathbf{g}}\_y & \sqrt{1 - \hat{\mathbf{g}}\_x^2 - \hat{\mathbf{g}}\_y^2} \end{pmatrix}^T,\qquad(3))$$

J. T. Freymueller, L. Sánchez (eds.), *5th Symposium on Terrestrial Gravimetry: Static and Mobile Measurements (TG-SMM 2019)*, International Association of Geodesy Symposia 153, https://doi.org/10.1007/1345\_2021\_136

M. Murzabekov (-) · V. Fateev · A. Pruglo · S. Ravdin VNIIFTRI, Mendeleevo, Moscow Region, Russia

e-mail: murzabekov@vniiftri.ru; fateev@vniiftri.ru; my@dinfo.ru; wurf@yandex.ru

where

$$
\hat{\mathbf{g}}\_{\times} = \frac{\mathbf{g}\_{\times}}{|\mathbf{g}|}, \ \hat{\mathbf{g}}\_{\times} = \frac{\mathbf{g}\_{\times}}{|\mathbf{g}|}.
$$

With this in mind, formula (2) will take the form:

$$\xi = -\arcsin\left(\hat{\mathbf{g}}\_{\times}\right), \qquad \eta = -\arcsin\left(\hat{\mathbf{g}}\_{\mathbf{y}}\right). \qquad (4)$$

Thus, the task of determining the values of the DOV components is reduced to calculating the components of the normalized vector **g**O.

#### **2 New Proposed Technique**

The main components of DZCS are a telescope, a CCD camera, an inclinometer, a GNSS receiver and auxiliary equipment. Measurement data with DZCS in each position of the telescope are: a frame of the starry sky, geodetic coordinates of the location, exposure time of the frame of the starry sky, the current tilt of the telescope from the readings of an inclinometer and a temperature value that can be determined, for example, from the built-in temperature sensor in the inclinometer.

The unknown parameters of the DZCS are:

1. Orientation angles ®, , (Euler angles) between the inclinometer coordinate systems (CS) and the CCD camera CS. They can be determined using rotation matrices around the axes (Zharov 2006):

$$R = R\_z(\boldsymbol{\varphi}) \cdot R\_y\left(\boldsymbol{\theta}\right) \cdot R\_z\left(\boldsymbol{\psi}\right). \tag{5}$$

2. The scale factors *mx* and *my* and angle " between the axes of the inclinometer. They can be calculated using a matrix of the form:

$$M = \begin{pmatrix} m\_x \cdot \sin \varepsilon \ m\_x \cdot \sin \varepsilon & 0\\ 0 & m\_y & 0\\ 0 & 0 & 1 \end{pmatrix} \tag{6}$$

3. Temperature coefficients *kx* and *ky* of the inclinometer axes. They can be determined if the temperature change *µ* during the observation is known:

$$n'\_{\underline{x}} = n\_{\underline{x}}^{meas} + k\_{\underline{x}} \cdot \Delta T,\tag{7}$$

$$
\hat{n'\_y} = \hat{n'\_y}^{meas} + \hat{k\_y} \cdot \Delta T,\tag{8}
$$

$$
\sigma\_y' = n\_y^{meas} + k\_y \cdot \Delta T,\tag{8}
$$

$$
\Delta T = T\_{end} - T\_0,\tag{9}
$$

where *kx*, *ky* – temperature coefficients; nmeas <sup>x</sup> ; nmeas <sup>y</sup> – measured inclinometer readings without taking into account displacement due to temperature; *T*0, *Tend* – temperature in the first and last stationary position of the telescope in a single series; n<sup>0</sup> x; n<sup>0</sup> <sup>y</sup> – corrected inclinometer readings taking into account changes due to temperature.


Based on the foregoing, the parameters of the DZCS are:

	- (a) CCD sensor orientation matrix in local CS, *A*;
	- (b) measured inclinometer readings on two axes, nmeas <sup>x</sup> ; nmeas <sup>y</sup> I
	- (c) temperature change during a single series, *T*.
	- (a) orientation angles between the inclinometer CS and CCD camera, ®, , ;
	- (b) scale factors and angle between the axes of the inclinometer, *mx*, *my*, ";
	- (c) temperature coefficients of the inclinometer axes, *kx*, *ky*;
	- (d) the components of the normalized gravity vector **g**O, gOx; gO<sup>y</sup> .

Thus, the vector of unknowns of DZCS is as follows:

$$X = \left(\varphi, \theta, \psi, m\_{\times}, m\_{\times}, \varepsilon, k\_{\times}, k\_{\times}, \hat{\mathbf{g}}\_{\times}, \hat{\mathbf{g}}\_{\times}\right). \tag{10}$$

#### **3 Model of the Proposed Technique**

Multiply the normalized vector **g**O by the matrix *A*, and then by the matrix *R*. This will allow to obtain the components of **g**O in the inclinometer CS. In order take into account the inclinometer parameters, multiply by the matrix *M* and to take into account the shifts caused by changes in temperature, add term *k* -*T*. Based on this, the model of the new technique can be represented as follows:

$$\begin{cases} \begin{array}{c} \left(\boldsymbol{n}^{mod}\right)\_{1} = M \ \boldsymbol{\lambda} \ \boldsymbol{R} \times A\_{1} \ \boldsymbol{\lambda} \ \hat{\mathbf{g}} + k \ \boldsymbol{\lambda} \ \boldsymbol{\Delta T},\\ \left(\boldsymbol{n}^{mod}\right)\_{2} = M \ \boldsymbol{\lambda} \ \boldsymbol{R} \ \boldsymbol{\lambda} \ \boldsymbol{A}\_{2} \ \boldsymbol{\lambda} \ \hat{\mathbf{g}} + k \ \boldsymbol{\lambda} \ \boldsymbol{\Delta T},\\ \boldsymbol{\lambda} \ \boldsymbol{\lambda} \ \boldsymbol{\lambda} \end{array} \end{cases} ,$$
 
$$\begin{array}{c} \boldsymbol{\lambda} = \boldsymbol{\lambda}, \boldsymbol{\lambda} = \boldsymbol{\lambda}, \boldsymbol{\lambda} = \boldsymbol{\lambda}, \boldsymbol{\lambda} = \boldsymbol{\lambda}, \boldsymbol{\lambda} \ \boldsymbol{\xi} + k \ \boldsymbol{\lambda} \ \boldsymbol{\Delta T},\\ \left(\boldsymbol{n}^{mod}\right)\_{N} = M \ \boldsymbol{\lambda} \ \boldsymbol{\lambda} \ \boldsymbol{R} \ \boldsymbol{\lambda} \ \boldsymbol{A}\_{N} \ \boldsymbol{\lambda} \ \hat{\mathbf{g}} + k \ \boldsymbol{\lambda} \ \boldsymbol{\Delta T}, \end{array} ,\tag{11}$$

where *N* – number of measurements (number of stationary positions of the telescope, N 10); nmod <sup>D</sup> nmod <sup>x</sup> nmod y <sup>T</sup> – modeled inclinometer measurements along the axes *OX* and *OY*; *M* – matrix for estimating inclinometer parameters (6); *R* – rotation matrix from the CS of the CCD sensor to the CS of the inclinometer (5); *A* – CCD sensor orientation matrix in the local CS; **g**O – normalized gravity vector in local CS (3); k D kx ky T – temperature coefficients of inclinometer axes; *µ* – temperature change during a single series. To evaluate all unknown model parameters, after processing the measurement data in a single series, the objective function is formed:

$$\begin{split} E &= \sum\_{i=1}^{N} \left| \left( \left( \left( n\_x^{mod} \left( X, A\_i, \Delta T \right) \right)\_i - \left( n\_x^{meas} \right)\_i \right)^2 \right. \\ &\left. + \left( \left( n\_y^{mod} \left( X, A\_i, \Delta T \right) \right)\_i - \left( n\_y^{meas} \right)\_i \right)^2 \right) \right| \end{split} \tag{12}$$

where *N* – number of measurements; nmeas x i ; nmeas y i – measured inclinometer readings in the *i*-th stationary position of the telescope, recounted in the projection of the normalized gravity vector (sines of the respective inclinometer reading angles). All unknown parameters are estimated by minimizing objective function (12). Nonlinear optimization starts with initial values:

$$X = \left(0, 0, 0, 1, 1, \frac{\pi}{2}, 0, 0, 0, 0\right).$$

Least squares optimization is performed using the Marquardt-Levenberg method with the numerical calculation of derivatives. At the same time, all parameters of the model are evaluated simultaneously, i.e. DZCS "autocalibration" occurs in each series. After estimation gO<sup>x</sup> ; gOy, DOV components calculated with (4).

Thus, the new observation technique with the DZCS involves obtaining frames of the starry sky, the values of the inclinometer in a single series of measurements at different tilts of the telescope and temperature change during a single series. Measurements in each series can be performed in arbitrary directions of the optical axis of the telescope and at arbitrary angles in the horizontal plane and differ from series to series. The main requirement for applying the new technique is the rigidity of the telescope - CCD camera - inclinometer system. The accuracy of the rotary device and the rigidity of the horizontal plane of the DZCS's location do not affect the accuracy of the final values.

The algorithm of the new technique is presented in Fig. 1. The algorithm of the new technique includes three stages:


One the of ways to implement the new proposed technique is presented in Fig. 2.

In this case, the observation process with the new technique involves rotating the telescope, CCD camera and inclinometer around the vertical axis with a fixed number of steps in the horizontal plane twice: with the "initial" zenith angle "<sup>1</sup> (first observation cycle) and with set zenith angle "<sup>2</sup> (second observation cycle). The "initial" zenith angle is understood as the state of the initial alignment of the DZCS in the horizontal plane according to the readings of the inclinometer.

The advantages of the proposed technique, compared with the existing traditional technique, are:


**Fig. 1** The algorithm of the new observation technique

**Fig. 2** An example of the implementation of a new measurement technique

observations on any solid foundation (dirt, asphalt roads and sites). This is especially important when measuring in the field.

#### **4 Research of a New Technique**

#### **4.1 Research of Change of the Calibration Coefficients**

Tests of the new measurement technique in the field were performed on a DZCS (Murzabekov 2017) at five various geographical points during 16 stellar observation nights. The average number of stars observed in the frame is about 100. At each point, at least six series of measurements were performed. According to the research results, it was found that the number of stationary positions of the telescope equal to 24 (12 positions for each rotation) is sufficient. Their further increase does not lead to an increase in accuracy.

In the traditional measurement technique, calibration coefficients determined before the start of measurements are used as constant values for observations. Studies have been conducted to evaluate changes in calibration coefficients between series and the effect of these changes on current values of DOV.

An example of the values of the calibration coefficients *mx* and *my* (inclinometer scale factors) for each series and their

change between the series during the observation at one point are presented in Fig. 3.

As can be seen from Fig. 3, there are a changes in the coefficients *mx* and *my*. Consider how the change in the coefficients *mx* and *my* between the series affects the DOV values for each series and the average values for all series. Figure 4 shows the changes in the calculated values of DOV depending on the series number (curves 1 and 2) and the change in calibration coefficients *mx* and *my* (curves 3 and 4).

As can be seen from Fig. 4, there is a clear correlation between and and the change in the coefficients *mx* and *my*. In this case, the differences of single series can reach up to 0.400, the average -D 0.1000, and D –0.0100.

Thus, an uncontrolled change in calibration coefficients leads to a shift in the values of DOV, i.e. to the appearance of an additional calculation error in DOV. This confirms the need to clarify the values of the calibration coefficients in each series during observations at each point.

#### **4.2 Research of the Influence of the Choice of Methods for Processing Observational Data on the Accuracy of DOV**

In the process of research of a new technique we reviewed:


The impact of the choice of processing method on the accuracy of the DOV is estimated. The total impact does not exceed 0.0300.

#### **4.3 Measurement Model with an DZCS**

A measurement model with an DZCS can be represented as follows:

$$\begin{cases} \begin{aligned} h\left(\xi, \mathcal{B}, \Phi, n\_{\chi}^{meas}, \chi p, \text{yp}\right) = 0, \\ h\left(\eta, L, \Lambda, n\_{\text{y}}^{meas}, t\_{UTC}, \chi p, \text{yp}\right) = 0, \end{aligned} \end{cases} \tag{13}$$

where *tUTC* – exposure time of the frame of the starry sky in the Universal Coordination Time (UTC); *xp, yp* – current polar motion parameters. In accordance with the measurement model (13), the formulas for calculating the measurement error of the DOV components with an DZCS are written in the following form:

$$\begin{cases} \begin{aligned} m\_{\xi} &= \pm K \cdot \sqrt{c\_{1}^{2} \cdot m\_{B}^{2} + c\_{2}^{2} \cdot m\_{\Phi}^{2} + c\_{3}^{2} \cdot m\_{n\_{x}^{\max}}^{2} + c\_{4}^{2} \cdot m\_{xp}^{2} + c\_{5}^{2} \cdot m\_{pp}^{2}}, \\ m\_{\eta} &= \pm K \cdot \sqrt{c\_{6}^{2} \cdot m\_{L}^{2} + c\_{7}^{2} \cdot m\_{\Lambda}^{2} + c\_{8}^{2} \cdot m\_{n\_{\eta}^{\max}}^{2} + c\_{9}^{2} \cdot m\_{\mathrm{trrc}}^{2} + c\_{10}^{2} \cdot m\_{xp}^{2} + c\_{11}^{2} \cdot m\_{\mathrm{pp}}^{2}, \end{aligned}} \end{cases} \tag{14}$$

**Fig. 4** Difference and , calculated with and without evaluation of calibration coefficients and change of *mx* and *my*

where *K* D 1.1 (with a confidence level of P D 0.95); *Ô<sup>i</sup>* – sensitivity coefficients for each component of the error (*i* D 1 ::: 11); *mj* – *j*-th component of the error (*j* D 1...9).

In more detail, the values of each error, sensitivity coefficients and calculation examples for several points are given in work (Murzabekov et al. 2021).

For example, for a point on the territory of FSUE VNI-IFTRI, the errors of the DOV components according to formulas (14) are: *m*- D 0.3600, *m* D 0.2400. Differences in errors are due to the dependence of the sensitivity coefficients for the DOV component in longitude on the latitude of the measurement point.

#### **4.4 Test Results**

The test results of the new technique are presented in Fig. 5.

As can be seen from Fig. 5, the value of the standard deviation for determining components of DOV is in the range 0.100–0.300.

**Fig. 5** Standard deviation for determining DOV at five geographical points during 16 observing starry nights

#### **5 Summary and Conclusions**

Thus, a new technique for performing observations with DZSC was developed. It provides the ability to "auto-calibrate" the parameters of the device during the measurement session in each series. In addition, the proposed measurement technique does not require the installation of a special hard measuring concrete base and precise rotation of the telescope around a vertical axis.

According to the results of testing the device in field conditions, the standard deviation was obtained in the range of 0.100–0.300, which is at the level of world analogues.

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/ licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Observations with gPhone Gravimeter in Moscow**

## Ekaterina Chistiakova

#### **Abstract**

Long-term recordings of temporal gravity variations were made with the gPhoneX # 117 spring gravimeter at the Center of Geodesy, Cartography and SDI (TsNIIGAiK). The purpose of the observations is to obtain reliable parameters of the Earth tides and non-tidal component for a site located within the territory of Moscow. At a later stage, these parameters will be used to process absolute gravity observations. One year of recordings were chosen (April 1, 2016–March 03, 2017). Preliminary data processing is fully automated and runs in a program written in Python. In the course of the research, the parameters of the instrument drift were determined by the method of piecewise linear approximation of the observations. The estimated parameters of the Earth tides were obtained according to the algorithms of the ETERNA 3.4 program (Wenzel, Bulletin d'Informations des Marées Terrestres 124:9425–9439, 1996). The non-tidal gravity changes (the difference between the experimental and local model data) were compared with the model of temporal variations of the gravity field caused by atmospheric loading. Theoretical values were provided by EOST Loading Service.

The studies presented here, are unique for Moscow site and they are one of the first attempts to use high-precision recordings of spring gravimeters while computing observed tidal and non-tidal parameters.

#### **Keywords**

Atmosphere loading estimates - Earth tides - Gravity variations - Moscow - Non-tidal gravity variations

#### **1 Introduction**

High-precision observations of gravity play a principle role in the modern study of the figure of the Earth and its external gravity field. It is important to know not only the absolute value of gravity at gravity network stations, but also the nature of the gravity variations over time. Due to the complex interior structure of our planet, gravity varies according to a complicated rule. Variations include phenomena such as the Earth tides, ocean and atmospheric tides, pole tides, changes in groundwater levels, non-tidal loading deformations, etc.

An algorithm for identifying, accounting, and forecasting all of these factors is developed in each case separately. There are some studies for the Earth tides modeling taking place in Russia (Spiridonov et al. 2015) but there are very few of them based on gPhone gravimeter observations (Antonov et al. 2016). In order to obtain and predict the local parameters of the Earth tides and analyze loading estimates at the TsNI-IGAiK site (® D 55.8550ı, œ D 37.5160ı), regular observations of temporary gravity variations have been made using a modern gPhone gravimeter. Such studies have never been taken place at Moscow site. It is very important to research gravity variation parameters at TsNIIGAiK site thus it is the main Fundamental Site of Russian State Gravity Network.

E. Chistiakova (-)

Center of Geodesy, Cartography and SDI (TsNIIGAiK), Moscow, Russia

J. T. Freymueller, L. Sánchez (eds.), *5th Symposium on Terrestrial Gravimetry: Static and Mobile Measurements (TG-SMM 2019)*, International Association of Geodesy Symposia 153, https://doi.org/10.1007/1345\_2022\_142

#### **2 Data Pre-Processing**

The gPhoneX # 117 spring gravimeter is installed on a pillar on the ground floor of the Center for Geodesy, Cartography and SDI (Russia, Moscow, Onegskaya St., 26). Continuous gravity measurements are being taken from November 2014 to the present. The interruptions of the recordings have been taken place due to power outages and calibration of the gravimeter. For this experiment, a continuous 350-days data series was selected (from April 01, 2016 to March 03, 2017). Frequency of measurements is 5 Hz. The output data used for the analysis are: averaged values of gravity variations, ambient temperature and pressure, levels positions with 1 s resolution.

Data pre-processing was automated in the Python 3 language and it consisted of eliminating spikes, noise and changing the sampling rate to 1 measurement per hour. Spikes with deviation greater than 5 -Gal were replaced by linearly interpolated values. Noises were cut by applying the Savitzky Golay filter (Savitzky and Golay 1964). Discretization was made taking into account the Nyquist-Shannon sampling theorem (Kotelnikov 1933).

The polar motion effect on gravity was calculated by the Formula, specified in the IAGBN Absolute Observations Data Processing Standards (1992). The Earth rotation parameters are freely available on the website of the International Earth Rotation Service (ftp://hpier.obspm.fr/iers/eop/ eopc04/).

The gravity sensor of the gPhone system consists of a low drift, metal, Zero Length Spring system (Microg LaCoste 2012), whose characteristics are determined by the mechanics of elastic material. Thus the gPhone spring gravimeter is characterized by a stable quasi-linear drift of small magnitude (in comparison with quartz gravimeters). The process of fitting instrumental drift is made as described by Kang et al. (2011). First, according to the characteristics of varying trends of gravity residuals time-series (gravity residuals after corrections of the solid Earth tide, ocean tidal loading and pole tide), the entire time-series are set into segments, then segmental polynomial fitting is applied to each segment, and finally combined together to compensate entire period for a long term gravimeter drift (Kang et al. 2011).

#### **3 Tidal Analysis**

Let's compare the differences in the experimental drift corrected data series with several models of the Earth tides. These differences will be called hereinafter a local gravity variation component. Theoretical variations were obtained for three models of the Earth tides: the harmonic expansion of Tamura (1987) (hereinafter—Tamura) for the solid Earth (recommended for data processing by the manufacturer of the gravimeter (Microg LaCoste 2012)), tide model WDD (Dehant et al. 1999) for an elastic Earth model (obtained in the program for processing force variations gravity TSoft (Van Camp and Vauterin 2005), as well as the model calculated in the Russian program ATLANTIDA3.1\_2014 (Spiridonov et al. 2015) using models of the Earth's structure IASP91 (Kennett and Engdahl 1991) and oceans FES2012 (Carrere et al. 2012) (hereinafter—Atlantida).

On average over the year, the standard and absolute deviations of the local gravity variation component are 22.0 and 122.9 -Gal for the WDD model, 14.8 and 82.0 -Gal for the Tamura model, 7.6 -Gal and 43.7 -Gal for the Atlantida model, respectively. Thus, the Atlantida model describes the tidal gravity variations at the TsNIIGAiK site in the best way. It is mainly connected with the fact that Atlantida model includes the ocean tides loading effect.

For analysis, we divide the annual interval into 14-days periods. During the year, the local gravity variation component has the same character, but different absolute and standard deviations. As an example, Fig. 1 shows a plot of the local gravity variation component for the three models in one of the two-week periods.

Absolute and standard deviations show different values in different seasons. Table 1 shows that the largest deviations are observed in the summer and winter months.

Thus, theoretical model of the Earth tides Atlantida can be used for some issues at TSNIIGAiK site. However, local parameters of the Earth tides are necessary to conduct highprecision studies. The calculations were performed in the program ETERNA renewed by Schueller (Schueller 2019). The parameters were obtained as corrections to the Tamura model. Characteristics of the local parameters are a topic for another article. Further we will discuss the differences between the experimental data and the local model. These differences will be called non-tidal gravity variations.

Here is the comparison of the local gravity variation component and non-tidal component. As an example, Fig. 2 shows a plot in one of the two-week periods.

It is possible to conclude that the non-tidal component does not exceed 2 -Gal in the standard deviation and 15 -Gal in the absolute deviations. The use of the local model of the Earth tides reduced unaccounted gravity variations by 3 times in comparison with using the theoretical model of Atlantida.

#### **4 Atmospheric Loading Effects**

Let's now consider the non-tidal gravity variations. The main factor that makes them up is the so-called atmospheric loading estimates.

**Fig. 1** Local gravity variation component: differences between observed drift corrected gravity variations and calculated gravity variations using the Earth tides models Tamura, WDD, Atlantida


**Table 1** Absolute and standard deviations of the local gravity variation component

Temporal variations in the distribution of atmospheric density cause changes in the gravitational attraction of air masses during local observations. Also, the load on the atmospheric masses deforms the Earth's crust and the sea surface. It is known that local measurements of gravity are currently corrected for atmospheric pressure variations using a standard correlation coefficient of 0.3 microGal/mbar, which corresponds to the Resolution IAG No. 9 of 1983. However, the current accuracy of relative gravity measurements makes it possible to identify the real dependences of temporary gravity variations not only on barometric pressure, but also on the state of the atmosphere as a whole (temperature, wind, humidity, etc.). For this purpose, the University of Strasbourg has founded the EOST Loading Service, where one can find models of the Earth's gravity field variations, tilts and displacements of the Earth's surface caused by atmospheric, hydrological and induced ocean non-tidal loads. There are stations that constantly monitor and analyze the influence of the state of the atmosphere on gravity variations to maintain this service.

The EOST Loading Service provides simulated time series of gravity variations caused by non-tidal atmospheric loading effects according to the data on the state of the atmosphere in the Moscow region [http://loading. u-strasbg.fr/GGP/atmos]. We have chosen the ERA5 model (Copernicus Climate Change Service 2017) for

**Fig. 2** Local gravity variation component and non-tidal component

**Fig. 3** Unaccounted gravity variations after correcting measurements with drift, Polar motion, local tide and atmosphere loading effects (ERA5 model)

processing as recommended in some studies (Olauson 2018).

Figure 3 shows non-tidal gravity variations at the TSNI-IGAiK point and gravity variations caused by atmospheric loading deformations computed using the ERA5 model. We

can see that the ERA5 model closely approximates the non-tidal gravity variations. Thus, the unaccounted gravity variations are reduced to half of the non-tidal variations after discounting the atmospheric loading effect, as it is shown in Table 2.

**Table 2** Absolute and standard deviations before and after using ERA5 model


#### **5 Conclusions**

In this paper, we have shown that the Atlantida model is best suited for modeling tidal gravity variations at the site located in Moscow. However, when researching high-precision gravity issues, theoretical modeling is not recommended to be used. Using local tidal parameters delivered from the analysis of 1-year gravity data allowed us to improve the results up to 3 times. The non-tidal signal component was less than 15 -Gal peak to peak.

Finally, in gravity monitoring studies, it is recommended to apply atmospheric loading corrections, like that provided by the EOST Loading Service (ERA5 model). Doing so reduces by 50% the unaccounted variations of the local tidal model.

#### **References**

Antonov AY, Valitov MG, Ducarm B, Ardyukov DG, Timofeev AV, Kulinich RG, Kolpashikova TN, Proshkina ZN, Sizikov IS, Nosov DA, Naumov SB (2016) Tidal effects by gravity observation, models and liquid core effect. Vestnik SGUGiT 2(34):34–47. [in Russian]


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/ licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

**Part II**

**Absolute Gravimetry**

# **On Uncertainties of Laser Interferometric Absolute Ballistic Gravimeters Due to Magnetic Effects in the Free-Fall Gravity Measurements**

L. F. Vitushkin, D. D. Gidaspov, F. F. Karpeshin, O. A. Orlov, I. S. Khasiev, and V. I. Sheremet

#### **Abstract**

At present, the uncertainty in measurements of free fall acceleration (FFA) by the designed absolute ballistic gravimeter (ABG), developed in VNIIM, attains microgal values (1 Gal D 1 cm/s2). At such a high level of accuracy, the effects of the interaction of a test body (TB) falling in the vacuum chamber of the ABG in the geomagnetic field could comprise sources of systematic errors. Furthermore, individual units and systems of the ABG itself also can be sources of the magnetic field (MF). We report about more rigorous calculations of the possible effects than those performed in the past. A consecutive selfconsistent method of calculation of the desired correction to the FFA has been developed. It results in the differential equation, including not only the field itself, but also its vertical first and second derivatives and variation of the velocity along the trajectory. The desired correction is obtained by its numerical solution on the basis of the parameters of the field generated by the magnet of the stepper motor. The derived correction proved to be at the level of a tenth part of microgal.

#### **Keywords**

Absolute ballistic gravimeters - Free fall acceleration - Geomagnetic field - Magnetic induction -Systematic error

#### **1 Introduction**

At present, the uncertainty in measurements of free fall acceleration (FFA) by absolute ballistic gravimeters (ABG) attains microgal values (1 Gal D 1 cm/s2) (e.g., (Vitushkin et al. 2020)). The general operation principles of laser-interference ABG and the principles of absolute measurements of FFA with their help are presented, for example, in (Vitushkin 2014; Absolutnye gravimetry 2017; Niebauer et al. 1995). At such a high level of accuracy, it is necessary to take into account a number of effects that influence the result of FFA

measurements. Such effects, being the sources of systematic errors, include, in particular, the effects of the interaction of a test body (TB) falling in the vacuum chamber of the ABG in the geomagnetic field. Furthermore, individual units and systems of the ABG itself also can be sources of the magnetic field (MF). Those units are the ion (magnetic discharge) pump, electric motor used in some ABG designs, and others. Estimates of a possible effect of inhomogeneous MF on ABG readings have been obtained in the past, e.g. in (Niebauer et al. 1995; Murata 1978, 1980) and references therein.

In (Niebauer et al. 1995), the effect of the interaction of a conducting TB with an inhomogeneous MF has been estimated in the case of the FG5 type ABG. The effect is due to the emerging Foucault currents. The measured fields of the ion pump magnet, servo motor and Faraday optical isolators happened to be less than the geomagnetic field, estimated at 50 -T, and do not create significant effects on FFA measurements with this type of ABG. It was proved experimentally that there was no any effect of this kind.

**V. I. Sheremet** was deceased at the time of publication.

L. F. Vitushkin · D. D. Gidaspov · F. F. Karpeshin (-) · O. A. Orlov · I. S. Khasiev

D.I.Mendeleyev Institute for Metrology (VNIIM), St. Petersburg, Russia

J. T. Freymueller, L. Sánchez (eds.), *5th Symposium on Terrestrial Gravimetry: Static and Mobile Measurements (TG-SMM 2019)*, International Association of Geodesy Symposia 153, https://doi.org/10.1007/1345\_2021\_133

To this end, MFs of the order of 100 -T were applied to the system with a set of Helmholtz coils. This experiment showed absence of effects at a level of 1 -Gal. It should be noted that in FG5-type gravimeters, the free-fall path of the TB is about 21 cm and the fall velocity attains 2 m/s.

In (Murata 1978, 1980), the motion of a TB charged due to the supposed contact phenomena in a geomagnetic field under the action of the Lorentz force was considered. The estimates also show a negligible effect of the geomagnetic field on the motion of the TB. As a consequence, gravimeters are usually designed without adjusting for the magnetic effects of eddy currents. In the present work, we report about more rigorous calculations, we have been undertaken for the design of the ABG Grot. The correction is calculated on the basis of the parameters of the MF generated by the magnet of the stepper motor. The derived correction proved to be at the level of a tenth part of microgal. To a great extent, it can be taken into account by calibrating and comparing the gravimeters. However, in some cases, the spread can significantly exceed the indicated value. This may be the case if magnetized parts are used. Therefore, it is important to understand the nature of the correction. Without this understanding, the readings of gravimeters of different types can lead to incorrect results.

Bearing this in mind, in the first place the principles of the ABG operation are reminded in the next section. The method of assessment of the correction is described in Sects. 3 and 4. Numerical estimations of the integral correction to apply to the determined value of FFA are made in Sect. 5. In the last Sect. 6, the results are summarized, conclusions are drawn, and further prospects are considered.

#### **2 Operation Principle of the ABG Grot**

Let us evaluate the uncertainty caused by induction currents and their interaction with the geomagnetic field, as well as with the MF produced by individual parts of the gravimeter. The calculation procedure for all types of MF sources is the same. The ABG scheme is usual. The TB falls from a small height of 15 cm. Its motion is monitored by the laser-interference method. The system generates a set of instantaneous values of time *ti (zi), i* D *1, 2, ..., n*, into which the TB passes the points *zi,* fixed in height and spaced halfwave apart from each other. Up to the vertical gradient , the FFA at the point *z* can be written as follows:

$$\mathbf{g}\left(\mathbf{z}\right) = \mathbf{g}\_0 + \boldsymbol{\chi}\left(\boldsymbol{z} - \boldsymbol{z}\_0\right),\tag{1}$$

where *g*<sup>0</sup> is the FFA value at the point *z*0. Choosing point *z*<sup>0</sup> as that where the TB is at the initial moment of time *t0* D 0, and designating the velocity of the TB at this point v*0*, we arrive at the following set of coupled equations, relating the successive times *ti* with the points *zi* on the trajectory:

$$z\_i = z\_0 + v\_0 t\_i + \frac{\mathbf{g}\_0 t\_i^{\;2}}{2} + \nu \left(\frac{v\_0}{6} t\_i^{\;3} + \frac{\mathbf{g}\_0}{24} t\_i^{\;4}\right). \qquad (2)$$

The *z0* and v*<sup>0</sup>* values at the initial time *t0* D 0 form a set of initial conditions, which unambiguously defines the trajectory.

#### **3 Physical Premises for the Uncertainties Arising in the Presence of MF**

Let us consider a simple model. A conductive loop drops down city of v*(z)* along the *z* axis in an inhomogeneous MF with induction **B**, remaining in the horizontal plane. Since the MF is inhomogeneous, the magnetic flux through the contour changes during the motion (Fig. 1). As a result, an electromotive force appears in the contour, which in turn causes circular electric current. As is known, the closed frame with current has the magnetic moment **m**:

$$\mathbf{m} = \frac{S^2}{R} \left( \mathbf{v} \frac{\partial B\_z}{\partial z} \right) \equiv A \left( \mathbf{v} \frac{\partial B\_z}{\partial z} \right) . \tag{3}$$

**Fig. 1** Scheme of the fall of the conductive loop in an inhomogeneous MF with induction **B**. The magnetic flux passing through the loop is conventionally depicted by the number of lines of force. When the loop falls, the number of field lines changes. As a consequence, this induces an electric current along the loop and the conjugate magnetic moment of the circuit **m**

In Eq. (3), *S* is the area of the loop, *R* – its electrical resistance. Parameter *A* characterizes the physical properties of the TB: its area, electrical conductivity, and the shape in general case. In turn, the produced magnetic moment interacts with the external MF. If the latter is inhomogeneous, then the loop is drawn into the region with the maximum gradient of the MF. Therefore, the force acting on it arises, which can be expressed as follows (Landau and Lifshitz 1960):

$$\mathbf{F} = -\nabla \left( \mathbf{m} \mathbf{B} \right) \equiv -\nabla U\_{\text{pot}},\tag{4}$$

where *U*pot is the potential energy of the interaction of the induced magnetic moment with the MF. Dividing the force Eq. (4) by the mass of the loop *m*0, one obtains the desired correction to FFA:

$$
\Delta \mathbf{g} = -\frac{A}{m\_0} \frac{\partial}{\partial z} \left( \upsilon B\_z \frac{\partial B\_z}{\partial z} \right) . \tag{5}
$$

#### **4 Equation for the Motion of the TB**

The fall of the TB is described by the second Newton's law:

$$\frac{d^2z}{dt^2} = \mathbf{g}(z) + a\_{e.m.}(z) \tag{6}$$

In Eq. (6), *g(z)* is the desired FFA (1). *ae.m. (z) g(z)* is the correction Eq. (5) taking into account the geomagnetic and other electromagnetic fields. The correction for the geomagnetic field is small. The correction for the MF of the ionic pump can be leveled by the optimal orientation of its magnet. Let us examine the correction for the MF created by the permanent magnet of the stepper motor.

Performing differentiation in Eq. (5), one can expand Eq. (4) to read:

$$\begin{split} \Delta \mathbf{g} &= -\frac{A}{m\_0} \left( B\_z \frac{\partial B\_r}{\partial z} \frac{\partial v}{\partial z} + v B\_z \frac{\partial^2 B\_r}{\partial z^2} + v \left( \frac{\partial B\_r}{\partial z} \right)^2 \right) \\ &\equiv F\_1 + F\_2 + F\_3. \end{split} \tag{7}$$

As one can see from Eq. (7), the desired correction depends not only on the height *z*, but also on the velocity of the TB at a given point, and moreover, on its own acceleration in the point, as

$$d\,v/dz = \frac{1}{v(z)}\frac{d^2z}{dt^2}.$$

Therefore, <sup>d</sup>2*<sup>z</sup>* dt <sup>2</sup> is included in Eq. (6) twice: both on the left and on the right. Thus, the force induced by the MF depends on the acceleration of the TB at the given moment and at the given point in the trajectory, which is itself determined by this force. For this reason, Eq. (6) must be solved in a self-consistent way, taking into account this interaction of acceleration with the force of electromagnetic induction. To this end, we collect both terms containing d v dt <sup>D</sup> <sup>d</sup>2*<sup>z</sup>* dt <sup>2</sup> on the left side of the equality. As a result, after simple transformations, one obtains the motion equation for the TB in the presence of the MF as follows:

$$\begin{split} \frac{d^2 z}{d I^2} &= \begin{array}{c} \frac{1}{1 + \frac{A}{m\_0 \underline{d}\_r} B\_\varepsilon \frac{\partial B\_r}{\partial \varepsilon}}\\ \times \left\{ \mathbf{g}(z) - \frac{A \frac{d\underline{r}}{dI}}{m\_0} \left[ \left(\frac{\partial B\_r}{\partial z}\right)^2 + B\_z \frac{\partial^2 B\_r}{\partial z^2} \right] \right\}. \end{array} \end{split} \tag{8}$$

By comparing Eq. (8) and Eq. (6), we obtain the desired correction to FFA:

$$\begin{split} \Delta \mathbf{g}(\mathbf{z}) & \equiv \ \begin{array}{c} \mathbf{a}\_{\varepsilon.m.}(\mathbf{z}) = \ & \frac{1}{1 + \frac{A}{m\_{0}\frac{d\varepsilon}{dI}}B\_{\varepsilon}\frac{\partial B\_{\varepsilon}}{\partial \varepsilon}}\\ \times \left\{ \mathbf{g}(\mathbf{z}) - \ & \frac{A\frac{d\varepsilon}{dI}}{m\_{0}}\left[\left(\frac{\partial B\_{\varepsilon}}{\partial \varepsilon}\right)^{2} + B\_{\varepsilon}\frac{\partial^{2} B\_{\varepsilon}}{\partial \varepsilon^{2}}\right] \right\} - \mathbf{g}(\mathbf{z}). \end{split} \tag{9}$$

#### **5 Results**

The total weight of the TB, including the triple prism, is *P* D 58.8 g. The calculation shows that the combined coefficient *At* D 0.022 m<sup>4</sup> / Ohm. Let the stepper motor be in the *x, y* plane. It is also necessary to take into account the design features of the device, in particular, that its geometric center is offset from the origin by approximately 1.5 cm in the horizontal plane in both *x* and *y* directions. The origin is specified by the *z* axis, which is directed down along the trajectory of the stepper motor in the point. The end point of the trajectory is 17.2 cm below the stepper motor. In the last section of approximately 8 cm length, the trajectory is measured using a laser interferometer, and several hundred points *ti* (*zi*) are used for processing. The measurements have shown that the MF of the stepper motor is approximately a dipole having components of the magnetic moment:

$$M\_x = 0.16 \text{ Am}^2, \quad M\_y = 0.017 \text{ Am}^2, \quad M\_z = 0.014 \text{ Am}^2.$$

A numerical solution of Eq. (8) with typical initial conditions was obtained by the Runge-Kutta method. A representative trajectory *z(t)* is plotted in Fig. 2a. It is close to the conventional parabola. The numerical solution also gives us values of instant velocity v*(t)* along with acceleration d v.t / dt . Once these values are obtained, one can find instant values of the electromagnetic induction correction *ae*. *<sup>m</sup>*.(*z*), which is plotted in Fig. 2a. By fitting the shape of the obtained trajectory to the experimental points of the trajectory *ti (zi)*, one can find parameters v*0*, *g0*. The implementation of this algorithm depends on a particular device, and also on external factors such as seismic vibrations of the base, frequency modulation of the laser beam of the gravimeter, **Fig. 2** A plot of a part of the calculated trajectory of the falling TB *z(t)* (**a**), and the related correction to FFA (**b**), solid curve, the dashed curve being a linear fit

etc. It can be implemented in a number of ways. However, in order to evaluate the integral correction to the gravity acceleration *ae*. *<sup>m</sup>*.(*z*), one can proceed as follows.

Let us use the linear approximation of the curve for *ae*. *<sup>m</sup>*.(*z*), as shown in Fig. 2b:

$$
\Delta \mathbf{g}(z) = k + pz \equiv a + b \left( z - z\_m \right). \tag{10}
$$

In Eq. (10), *zm* is the point of measuring the *g0* value. A ¦<sup>2</sup> fit has resulted in the following values of the parameters: *k* D 1.48666 -Gal, and *p* D 0.194504 -Gal/cm. First, the latter values provide the following value of the desired correction in the lowest point of the trajectory: *a* D 0.011 -Gal. Bearing in mind that typical values of *ae*. *<sup>m</sup>*.(*z*) in the trajectory of Fig. 2 amount to tenth parts of microgal, one can assume that such a small final value of the correction might be due to a random compensation of individual terms. In any case, one can conclude that the desired correction is at the level of 0.1 -Gal. Second, using the *p* and *k* values obtained above, one can derive the value of *b* D 0.1945 -Gal/cm from Eq. (10) in a similar way. The value of *b* provides a correction to the gradient of FFA ". At the same time, the value of " can be measured independently by a relative gravimeter. In principle, this makes it possible to compare the present theory with the experiment.

#### **6 Discussion of the Results**

A consecutive self-consistent method of calculation of the desired correction to the FFA has been developed. The resulting Eq. (8) explicitly takes into account the fact that the correction is determined not only by the field itself, but also by its vertical gradient, as well as by the variation of the velocity along the trajectory. For this purpose, a consecutive consideration of all the three terms *F1, F2*, and *F3* in Eq. (7) is performed. The correction is obtained by solving Eq. (8) numerically for a representative trajectory of the TB. The need for self-consistency is due to the fact that the correction depends on the acceleration, on which it has a direct effect. It should be noted that in the previous papers (Absolutnye gravimetry 2017; Niebauer et al. 1995; Murata 1978, 1980) only the term *F3* was taken into account, and self-consistency was not considered. The obtained estimates indicate that the desired value of correction is within the tenth part of microgal.

The developed method also makes it possible to calculate the correction to the gravity gradient at the geographical point of measuring. In principle, this can serve as a test for comparing the theory with the experiment, after the gradient has been measured using a relative gravimeter.

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons. org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Evaluation of Systematic Errors in the Compact Absolute Gravimeter TAG-1 for Network Monitoring of Volcanic Activities**

Akito Araya, Keisuke Kasai, Masato Yoshida, Masataka Nakazawa, and Tsuneya Tsubokawa

### **Abstract**

Volcanic activities sometimes involve gravity changes, and this research is intended to establish an observation network surrounding an active volcano using compact absolute gravimeters. To simplify the configuration of absolute gravimeters, they are preferably operated with a light source distributed from a telecom band (wavelength of 1.5 m) laser through optical fibers. To evaluate the accuracy of the absolute gravimeter with the telecom band laser, we conducted observations using a prototype gravimeter (TAG-1) with frequency-stabilized lasers at both 1.5 m and 633 nm, and compared these results with the expected gravity at the site. Initially, both results showed offsets -187 -Gal and -9.6 -Gal for the 1.5-m laser and the 633-nm laser, respectively (1 Gal D 10-<sup>8</sup> m/s2). By correcting the systematic errors of the photo detectors measured by the synthetic chirp signal, the obtained absolute gravity was consistent with the expected value for both wavelengths; offsets from the expected gravity were reduced to 6.6 -Gal and 5.4 -Gal for 1.5 m and 633 nm, respectively. We also evaluated the errors associated with long-distance transmission of the 1.5-m laser using a reeled optical fiber (26 km) and an optical amplifier and found no degradation in the gravity data from the case of short transmission (10 m). These results allow networking of compact absolute gravimeters connected by telecom optical fibers that are operated using a common laser and expansion to volcanic areas to monitor the gravity change associated with volcanic activities.

#### **Keywords**

Absolute gravimeter Frequency stabilization Telecom band laser Volcanic observation

#### **1 Introduction**

Volcanic activities often involve earthquakes, crustal deformation, magnetic anomaly, and other geophysical phenomena. Gravity change is a direct probe of the mass movement

A. Araya (-)

K. Kasai · M. Yoshida · M. Nakazawa Research Institute of Electrical Communication, Tohoku University, Miyagi, Japan

T. Tsubokawa Shin-ei Keisoku, Iwate, Japan of magma that is monitored for the prediction of volcanic eruptions and to evaluate the transition of volcanic activities. Both types of gravimeters, relative gravimeters and absolute gravimeters, have been used, and the latter can measure longterm gravity changes without any instrumental drift with reference to an accurate wavelength of the frequency-stabilized laser. However, owing to the complex mechanism, large size, and high cost, absolute gravimeters have not commonly been used for volcanic observations. This research is intended to establish a network of compact absolute gravimeters for volcanic observations.

To construct this observation network, absolute gravimeters are preferably operated with telecom band (wavelength of 1.5 m) lasers distributed to each gravimeter via optical

Earthquake Research Institute, University of Tokyo, Tokyo, Japan e-mail: araya@eri.u-tokyo.ac.jp

J. T. Freymueller, L. Sánchez (eds.), *5th Symposium on Terrestrial Gravimetry: Static and Mobile Measurements (TG-SMM 2019)*, International Association of Geodesy Symposia 153, https://doi.org/10.1007/1345\_2020\_107

fibers because conventional lasers (wavelength of 633 nm) cannot be transmitted to distant sites because of the loss of optical fibers; which is typically 15–30 dB/km for 633 nm light, and 0.2–1 dB/km for 1.5-m light. To evaluate the accuracy of the absolute gravimeter with a telecom band laser, we conducted observations using a prototype gravimeter (TAG-1) with frequency-stabilized lasers at both 1.5-m and 633-nm wavelengths, and compared these results with the expected gravity of the site.

#### **2 Gravity Change Associated with Volcanic Activities**

Sakurajima is one of the most active volcanos in Japan. A devastating eruption occurred in 1914, and small eruptions still continue. It has been determined that the amount of magma in the magma chamber beneath the mountain is coming to that of 1914. Okubo et al. (2013) observed gravity changes associated with the volcanic eruption of Sakurajima using an absolute gravimeter. The gravity change was 10 -Gal (1 Gal D 10-<sup>8</sup> m/s2) at a distance of 2 km from the crater; this is only one factor larger than the background noise level due to local disturbances such as groundwater (Kazama and Okubo 2009). Therefore, observations near the crater and networking with a number of gravimeters surrounding the crater will significantly enhance the detectability of magma motion near the source by averaging the local disturbances using a number of sensors.

#### **3 TAG-1 Gravimeter**

Araya et al. (2014) developed a compact absolute gravimeter, TAG-1, and we used it for the evaluation of the availability of network monitoring of volcanic activities. To evaluate the accuracy of the absolute gravity with the telecom band laser, we conducted observations using TAG-1 with frequencystabilized lasers at wavelengths of 1.5 m and 633 nm, and compared these results with the expected gravity of the site. Figure 1 shows a picture and a schematic diagram of the TAG-1 gravimeter. It is comprised of a dropper for the free-fall mass and a built-in accelerometer for correction of seismic vibrations. By applying the built-in accelerometer and the small dropper, TAG-1 is compact and transportable for observations.

The laser light is introduced into the optical unit through the optical fiber and is incident to vacuum chambers confining the free-fall mass (Free-fall mirror) dropper and a reference pendulum (Reference mirror), both of which include retro reflective mirrors forming an interferometer. The interfered light is guided to photo detectors (PDs) through optical fibers. TAG-1 uses a quadrature interferometer for the displacement measurement of the free-fall mass, and the optical phase is calculated from the detected signals (Heydemann 1981; Svitlov and Araya 2014). From the quadratic dependence of the displacement with respect to time, the absolute gravity can be obtained. Effects of ground vibration acceleration on the gravity are corrected using data from the build-in accelerometer using the reference pendulum.

TAG-1 can be operated at both wavelengths of 1.5 m and 633 nm by using the PDs and the optical unit designed for each wavelength. InGaAs-type and Si-type PDs are used for the wavelengths of 1.5 m and 633 nm, respectively. Beam verticality can be adjusted so that the measuring laser beam and beam of the auto-collimator reflected on a reference alcohol surface are in parallel. For the vertical adjustment at the invisible 1.5-m wavelength, the measuring laser beam is monitored using an IR (infrared) viewer.

#### **4 Observation Using a Conventional 633-nm Laser and a Telecom Band (1.5 m) Laser**

We performed gravity measurements in a basement room of the main building of the Research Institute of Electrical Communication (RIEC), Tohoku University, using TAG-1 operated with both a conventional iodine-stabilized 633-nm He–Ne laser (wavelength of œ<sup>6</sup> D 632.99081163 nm), and a telecom band (1.5 m) laser (œ<sup>15</sup> D 1,538.803242 nm, Fig. 2) which was frequency-stabilized using the acetylene linear absorption spectrum with a linewidth of 500 MHz (Kasai et al. 2016) and whose frequency accuracy is estimated to be 10-<sup>9</sup> (Nakagawa and Onae 2004); the latter may realize the long-distance distribution of the light source and networking of gravimeters. The systematic errors in operation for both wavelengths were evaluated.

The free-fall mirror was dropped every 2 min. The obtained data were corrected for seismic noise measured by the build-in accelerometer. Figure 3 shows the measured gravity using the conventional 633-nm laser and the 1.5-m frequency-stabilized laser. The theoretical gravity variation is shown by the red line based on calculated tidal gravity and the absolute gravity measured by the relative measurement from a gravity reference point, as described in the following **Fig. 1** Picture (upper) and schematic diagram (lower) of the TAG-1 gravimeter. The optical unit shown in dashed red line in the lower figure can be replaced depending on the laser wavelengths

**Fig. 2** The 1.5-m frequency-stabilized laser used for the measurements

**Fig. 3** Gravity (blue dots) measured using the 1.5-m and 633-nm lasers

section. Slow tidal gravity variations were commonly observed for both cases, while offsets were apparent. The offsets averaged in the periods were -187 -Gal for the 1.5-m laser and -9.6 -Gal for the 633-nm laser. This may be due to the PDs because mechanical and optical configurations are essentially the same for both cases, and the laser wavelengths are well defined. We evaluated the systematic error of TAG-1 caused by the frequency response of the PDs.

#### **5 Systematic Error Evaluation of the PDs**

Because the frequency of the fringe signal is almost proportional to the velocity of the free-fall mass, the response of the PD causes a systematic error in the gravity measurement (Niebauer et al. 1995). To evaluate the error directly, synthetically modulated laser light that simulates the interferometer fringe was applied to the PD, and the difference in obtained gravity values from the measurement and from calculations could be regarded as estimates of the systematic error. The free-fall mass in gravity, *g,* generates a chirp interferometer signal with a frequency rate of d*f* /d*t* D 2 *g*/œ, where œ is the wavelength of the laser. Therefore, chirp frequency rates of 13.07 MHz/s and 30.96 MHz/s produce g D 9.8 m/s2 for œ D 1.5 m, and 633 nm, respectively. The laser intensity was modulated using an electro-optic amplitude modulator for the 1.5-m evaluation, while a laser diode was used for the 633 nm light source. In each case, the chirp signal was applied using a function generator whose clock and data sampling clock were both locked to the same Rb time base. For the measurement of the 1.5-m laser, the chirp frequency was set to change from 0.1 MHz to 2.6 MHz in 0.2 s, and it

**Fig. 4** The observed gravity (blue/light blue dots, left axis) and the estimated error of the PDs obtained from the synthetic chirp signal (green dots, right axis). Observed 1 and observed 2 are calculated from the 1.5-m datasets within 0:30–1:30 and 1:30–2:30, respectively, on 19 December, 2017, as shown in Fig. 3

gives 12.5 MHz/s and g15 D 9.6175202625 m/s2 for œ15. For the measurement of the 633-nm laser, the chirp frequency was set to change from 0.1 MHz to 6.3 MHz in 0.2 s, and it gives 31 MHz/s and g6 D 9.8113575803 m/s2 for œ6.

In the data processing of TAG-1, the measured displacement of the free-fall mirror from *t*<sup>0</sup> (time from the release) to *t*<sup>0</sup> C *t* was fitted using a quadratic function of time, and the gravity acceleration, *g*(*t*0), was obtained from the quadratic coefficient. For small *t*0, just after a short time from release, *g*(*t*0) showed disagreement with the theoretical dependence on *t*<sup>0</sup> because the release of the free-fall mirror induced a slight recoil vibration; although the total free-fall time was approximately 180 ms, the analysis interval was fixed at *t* D 80 ms, and *g*(*t*0) was calculated with changing *t*<sup>0</sup> to assess the effect of the recoil vibration. The same calculation was applied to the detection of the synthetic chirp signal and the systematic errors of the PDs were estimated by the difference in *g*(*t*0) from the calculated gravity (g15 or g6) obtained from the frequency rate of the chirp signal. Figure 4 shows the observed *g*(*t*0) (blue and light blue dots, left axis) and the estimated error of the PDs obtained from the synthetic chirp signal (green dots, right axis). The observed data for *t*<sup>0</sup> > 60 ms agreed with the error estimation of the PDs. The estimation shows small systematic errors for whole *t*0, and smaller *t*<sup>0</sup> gives smaller gravity acceleration; the decrease of observed gravity at 1.5 m in Fig. 3, calculated with *t*<sup>0</sup> D 10 ms and *t* D 150 ms, is consistent with this estimates of errors of the PDs.

**Fig. 5** The observed data corrected using the estimated errors for the 633-nm and 1.5-m lasers. Theoretical tides were removed and then averaged. The expected level based on the relative measurements were referenced to AOB-B, as shown by the red dashed line. To see the consistency of the corrected gravity, two data sets for each laser were calculated and showed similar results

The observed data, calculated with *t*<sup>0</sup> > 60 ms and *t* D 80 ms, were corrected using the estimated errors, as shown in Fig. 5. Theoretical tides were removed from the data and then averaged. In contrast to Fig. 3, observed data with the 633-nm laser and 1.5-m laser agreed well after the correction; moreover, they were consistent with the expected level (red dashed line) based on the relative measurements using a LaCoste gravimeter referenced to the Aobayama gravity reference point (AOB-B) where the absolute gravity has been determined. The AOB-B is located 2.3 km westward from RIEC (Fig. 6). The mean offsets shown in Fig. 5 were 6.6 -Gal for the 1.5-m laser and 5.4 -Gal for the 633-nm laser, which were compared with -187 -Gal and -9.6 -Gal, respectively, without the correction in Fig. 3.

We also evaluated the errors associated with long-distance transmission of the 1.5-m laser through the optical fiber. The laser light at 1.5 m was introduced through short (10 m) or long (26 km) optical fibers. In this experiment, we used a reeled optical fiber in the laboratory. As shown in Fig. 7, the measured absolute gravity did not change and showed no degradation even when the laser was provided through a 26-km-long optical fiber and an optical amplifier. Nevertheless, to estimate errors in a practical system in the field, environmental effects on the optical fibers, such as vibration and thermal disturbances, need to be measured.

2.3km

**Fig. 6** Location of the Aobayama gravity reference point, AOB-B

#### **6 Conclusions**

The compact absolute gravimeter, TAG-1, was successfully operated with both 633-nm and 1.5-m lasers. By correcting systematic errors of the PDs measured using a synthetic chirp signal, the obtained absolute gravity was consistent with the expected value for both wavelengths; the systematic error of 1.5-m PDs was estimated to be as much as -190 -Gal without the correction. These results can lead to networking of compact absolute gravimeters connected via telecom optical fibers operated using a common laser and can be expanded to volcanic areas to monitor the gravity change associated with volcanic activities.

**Fig. 7** Measured gravity (blue dots) with the 1.5-m laser through a short optical fiber (10 m, upper) and a long optical fiber (26 km) with an optical amplifier (lower)

**Acknowledgments** We appreciate Prof. Yoshiyuki Tanaka, University of Tokyo, and Prof. Satoshi Miura, Tohoku University, for support and arrange of the relative gravity measurements between AOB-B and RIEC. This study was partly supported by the Ministry of Education, Culture, Sports, Science, and Technology (MEXT) of Japan, under its Earthquake and Volcano Hazards Observation and Research Program.

Part of this work was carried out under the Cooperative Research Project Program of the Research Institute of Electrical Communication, Tohoku University.

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons. org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Design of New Launch and Interferometer Systems for the IMGC-02 Absolute Gravimeter**

Andrea Prato, Sergio Desogus, Claudio Origlia, Marco Bisi, and Alessandro Germak

### **Abstract**

For the measurement of the acceleration due to gravity, INRiM developed a transportable ballistic rise-and-fall absolute gravimeter, the IMGC-02. It uses laser interferometry to measure the symmetrical free rising and falling motion of a test mass in the gravity field. The launch system is composed of a moveable carriage fixed to two pairs of springs loaded by an electric stepper motor, which vertically throw up a corner cube retroreflector in vacuum. The interferometer system is a modified Mach-Zehnder interferometer where the launched corner cube acts as the reflector in one of the optical arms of the interferometer and the other retroreflector acts as inertial reference during the measurement. However, both systems entail some practical problems and uncertainty contributions that need to be reduced. In particular, the current launch system might cause beam shear and rotational effects due to unavoidable small different loadings of the springs, while the current interferometer system poses problems in the alignment of the mirrors, which is a highly time-consuming procedure and has to be performed before and, sometimes, during the measurement session. For this reason, a new launch system consisting of an electric linear motor which produces a linear force along its length, and a modified Jamin interferometer system entailing a simpler alignment and a better stability in time, have been designed. This works deals with the description of these new systems.

#### **Keywords**

Absolute gravimeter - Gravimetry - Interferometer -Rise-and-fall gravimeter

#### **1 Introduction**

Absolute measurements of the acceleration due to gravity, *g*, are performed by absolute gravimeters, traceable to the units of length and time through their primary standards. For this purpose, INRiM developed a transportable ballistic rise-andfall absolute gravimeter, the IMGC-02, which is recognized as the Italian primary standard because of the WGG/0441 resolution<sup>1</sup> and, as such, is listed in the BIPM KCDB.<sup>2</sup> The measurement of *g* is performed using the reconstructed rise and fall motion of a corner cube retroreflector, which moves vertically in vacuum. An interferometer system is implemented in order to obtain time and space coordinates of the trajectory using a visible laser beam. The interferometer measures the distance between a free-falling corner

A. Prato (-) · S. Desogus · C. Origlia · M. Bisi · A. Germak INRiM – Istituto Nazionale di Ricerca Metrologica, Torino, Italy e-mail: a.prato@inrim.it

<sup>1</sup>Resolution WGG/04-41 of the first joint meeting of the CCM WGG and SGCAG (26-27 May 2004, BIPM): *the first joint meeting of the CCM WGG and SGCAG recognized the absolute ballistic method of measurement of the acceleration due to gravity as a primary method*. 2https://www.bipm.org/kcdb/cmc/search?domain=PHYSICS&areaId= 4&keywords=italy&specificPart.branch=19&specificPart.service= 41&specificPart.subService=124&specificPart.individualService= 404&\_countries=1&publicDateFrom=&publicDateTo=&unit= &minValue=&maxValue=&minUncertainty=&maxUncertainty=.

J. T. Freymueller, L. Sánchez (eds.), *5th Symposium on Terrestrial Gravimetry: Static and Mobile Measurements (TG-SMM 2019)*, International Association of Geodesy Symposia 153, https://doi.org/10.1007/1345\_2020\_108

**Fig. 1** Schematic drawing and picture of the IMGC-02 absolute gravimeter

cube retroreflector and a second retroreflector mounted on the quasi-inertial mass of a vibration isolation system. The gravimeter is composed of five main parts: (1) a launch system in a vacuum cylinder, (2) an interferometer system connected to an anti-vibrating system, (3) a laser body, (4) a photodetector and (5) a supporting frame, as shown in Fig. 1. A detailed description of the gravimeter can be found in D'Agostino (2006) and D'Agostino et al. (2008). However, the current launch system, which is composed of a moveable carriage supported by two pairs of springs, and the interferometer system, which is a modified Mach-Zehnder interferometer, introduce uncertainty contributions and practical problems that need to be overcome. For this reason, new launch and interferometer systems have been designed. This work deals with the description of the current systems and their future improvements.

#### **2 The Launch System**

The current launch system is composed of the vacuum chamber, the test mass and the launching pad. The launch chamber can be approximated to a cylinder with a diameter of 14 cm and a height of 65 cm. The base of the vacuum chamber is made of stainless steel and is supported by three legs equipped with levelling screws, which allow the vertical alignment of the chamber. The main part of the vacuum chamber is a flanged glass pipe, whose bottom part is fitted to the base, whereas the upper part is sealed with an aluminium cover. A BK7 glass window (5 cm diameter, 1.27 cm thick, œ/10 flatness, parallelism better than 2 arcsec) is positioned at the centre of this cover and allows the laser beam to reach the test mass. Connections are sealed by O-rings and, inside the glass pipe, a Faraday cage shields the test mass from electrostatic charges. The base of the chamber has two arm pipes: the first one is connected to a low noise turbo pump, the second to an ionisation vacuum gauge. A rotary pump carries out the coarse evacuation. A pressure of 1 10-<sup>3</sup> Pa is reached after about 5 h.

The launching pad (Fig. 2) can be approximated to a <sup>5</sup> <sup>5</sup> 20 cm<sup>3</sup> parallelepiped and is tightly connected to the base of the vacuum chamber. It supports, throws up and catches the test mass at the end of its trajectory. The test mass is a corner cube retroreflector, which has the property to reflect an incident ray parallel to itself regardless of the angular orientation. The corner cube has a mass of 0.15 kg and lies horizontally on three supports, which are the vertices of an equilateral triangle. These supports constitute the upper ending part of the catcher, which is fixed to a moveable carriage (0.15 kg) on a vertical linear bearing rail, screwed to an aluminium rectangular support. The carriage, in turn, is fixed to two pairs of series springs and is retained by an electromagnet. The starting horizontal position of the test mass, resting on the catcher, corresponds to the best alignment of the interferometer, i.e. the best overlapping between the test and reference beams. The thrust of the system (around 27 N) is given by the two pairs of series springs (elastic modulus approximately 1.1 N mm-<sup>1</sup> each),

**Fig. 2** Schematic drawing (left) and picture (right) of the current launching pad

which are arranged in parallel and loaded through a screw gear driven by an electric stepper motor. The vertical motion of the test mass is triggered by cutting off the exciting current of the magnet. The impulse is transmitted through the three supporting points to the test mass. When the mechanical action ends, the body hovers and its trajectory is tracked. The resulting total stroke of the mechanical system, before the release of the test mass, is around 25 mm. When released, the test mass has a speed of about 2 m s-<sup>1</sup> and travels in the vacuum chamber a distance of about 20 cm. The launch system is automatically reloaded after the pad catches the free-falling test mass.

Ideally, the corner cube retroreflector realizes a point in space. Unfortunately, the corner cube's centre of mass is not coincident with its optical centre. For this reason, the corner cube is fitted in an aluminium frame, designed in such a way as to move the centre of mass of the assembly as close as possible to the optical centre. In this way, any rotation of the test mass around its centre of mass should affect neither the interferometer alignment nor the measured *g* value. Despite

the highest accuracy in the realization of the support, a small but unavoidable difference between the optical centre and the centre of mass is still present. This entails that beam shear effects and rotational effects might occur due to movements and rotations of the test mass on the horizontal plane, caused by (small) different loadings of the springs which generate, during the upward launch, lateral forces and moments on the moveable carriage. As a consequence, the two interfering beams can translate relative to each other; launches affected by a fringe visibility reduction above a fixed threshold (usually 10–20% of the maximum intensity) are rejected. The rotational effect, instead, introduces a centripetal acceleration, whose vertical component is added to the local gravity acceleration *g*. Experimental tests (D'Agostino 2006) show that the average rotation of the test mass during the trajectory is 10 mrad and the associated angular velocity is 25 mrad s-1, given a total flying time of around 400 ms. Considering the uncertainty in locating the optical centre and the centre of mass as a random variable with a uniform distribution within ˙0.1 mm, the expected uncertainty due to the rotation of the

**Fig. 3** Schematic drawing (left) and picture (right) of the new launching pad

corner cube is 3.6 -Gal and represents one of the largest contribution to the uncertainty budget (D'Agostino et al. 2008).

To overcome these issues, a new launch system has been designed (Fig. 3). It can be approximated to a 7 7 26 cm<sup>3</sup> parallelepiped and consists of an electric linear motor, which has its stator and rotor unrolled to produce a linear force along its length, to which the moveable carriage and the catcher are fixed. As in the current system, the moveable carriage slides on a linear bearing rail screwed to an aluminium rectangular support. In this way, once aligned, beam shear and rotational effects should be minimized. A 50% decrease in the number of rejected launches is foreseen, and the rotational effect uncertainty is expected at least to halve. The linear motor has a maximum stroke of 780 mm, a peak force of 67 N and a continuous maximum force of 14 N. The presence of motor flanges enables the easy mounting of the linear motor on the aluminium structure, while the clamping plate design enables quick assembly and disassembly of the linear motors without disassembling the flange. The rate of movement of the magnetic field is electronically controlled via software, which allows one to track the motion of the rotor. In this way, it is also possible to design the downward motion profile in order to collect the test mass at the end of its falling motion.

#### **3 The Interferometer System**

The current interferometer system is a modified Mach-Zehnder interferometer (Germak et al. 2002; D'Agostino 2006). A scheme is reported in Fig. 4. The light emitted by the He-Ne laser (œ 633 nm) passes through an optical Faraday isolator (FI) which avoids any backreflection into the laser, possibly changing its wavelength and thus affecting the measurement. After the Faraday isolator, the beam is enlarged to 2 mm spot size by a beam expander (BE) and proceeds vertically toward a beam splitter (BS). The reflected beam proceeds horizontally toward a small corner cube (SCC) used to observe the vertical orientation of the beam through a telescope (T). The transmitted beam, vertically directed, is deviated by a moveable mirror (M1) which can be rotated around two axes to adjust the direction of the shifting arm of the interferometer to the true vertical. These elements are needed to check the verticality of the laser beam. Afterwards the beam enters horizontally into the optical prism (OP) situated beneath the reference corner cube retroreflector (RR). The OP has been designed in such a way that one half is a reflecting mirror and the other half is a beam splitter. The incident light is divided into two beams by OP. One of the beams, which represents the fixed arm of the interferometer, proceeds horizontally and goes straight toward the detector (D); the second beam travels vertically down to the test mass's corner cube retroreflector (TR) and forms the shifting arm of the interferometer; it is reflected back to RR and strikes the movable mirror (M2) after being reflected by the OP. The adjustment of M2 allows the beams to recombine at the second beam splitting portion of OP. The interference fringe causes cyclic variations in intensity, which are due to the change in the difference between the two optical path lengths at every half-wavelength displacement of the test mass retroreflector TR. The detector converts the output light of the interferometer to an electric signal. The test mass trajectory is measured by timing this electric signal. Since the recombining beams have to be coaxial in order to avoid distortions on the laser interference fringes, the angular position of the mirror M2 has to be adjusted every time before the measurement and has to be monitored during the measurement. The alignment is achieved by rotating the mirror M2 around its two axes through a piezoelectric tilt actuator.

**Fig. 4** Scheme of the current modified Mach-Zehnder interferometer (left) compared to the new modified Jamin interferometer (right). Abbreviations: OP optical prism, BE beam expander, FI Faraday iso-

lator, BS beam splitter, D detector, T telescope, IP interference pattern control, M mirror, RR reference retroreflector corner-cube, TR test-mass retroreflector corner-cube, SCC small corner-cube

**Fig. 5** 3D-Scheme of the new modified Jamin interferometer on the horizontal plane

**Fig. 6** Bottom view of the 3D-Scheme of the new modified Jamin interferometer integrated in the IMGC-02 rise-and-fall ballistic gravimeter

Unfortunately, this operation entails practical problems, is highly time consuming and has to be performed before and during the measurement session. For this reason, a new modified Jamin interferometer (Shamir 1999) has been devised (Fig. 4). Such system is similar to the modified Mach-Zehnder interferometer except that the two beams directly recombine on the OP, thus the movable mirror M2

is removed. A 3-D detail of the main part of the new system is depicted in Fig. 5, while its assembly in the gravimeter is shown in Fig. 6. The alignment of the recombined beams is possible by just shifting the reference corner cube retroreflector RR along the horizontal plane. The main advantages are a simpler and faster alignment of the two beams and a better stability in time entailing an expected setting time reduction of around 25% and a further decrease of rejected throws. The new scheme has a potential Abbe error because of the removal of M2, which was used to correct the misalignment of the beams traveling in the two arms of the interferometer; using high optical quality corner cubes and prism (angular accuracy within 1 arcsec and a flatness within œ/10) the uncertainty contribution due to the Abbe error is minimized.

#### **4 Conclusions**

For the measurement of the acceleration due to gravity, INRiM developed a transportable ballistic rise-and-fall absolute gravimeter: the IMGC-02. Currently, the launch system is composed of a moveable carriage supported by two pairs of series springs in parallel, which are loaded by moving down the carriage through a screw gear driven by an electric stepper motor. Nevertheless, it is likely that lateral forces arising during the upward launch of the test mass, caused by unavoidable small different loadings of the two springs, make the test mass move on the horizontal plane and rotate. Therefore, the two interfering beams can translate relative to each other, so that beam shear and rotational effects might occur. To overcome these issues, a new launch system has been designed. It consists of an electric linear motor, which has its stator and rotor unrolled to produce a linear force along its length, fixed to the moveable carriage. The rate of movement of the magnetic field is electronically controlled to track the motion of the rotor. In this way, beam shear and rotational effects should be minimized. For what concern the interferometer system, at present a modified Mach-Zehnder interferometer is adopted. Unfortunately, the alignment operation entails practical problems, is highly time-consuming and has to be performed before and, sometimes, during the measurement session. For this reason, a new modified Jamin interferometer has been devised. This system is similar to the modified Mach-Zehnder interferometer except that the two beams directly recombine on the optical prism, thus the movable mirror is removed. The main advantages are a simpler alignment of the two beams and better stability in time. The new interferometer system should guarantee measurements that are more robust and a faster setup of the gravimeter.

#### **References**


Shamir J (1999) Optical systems and processes. SPIE, Bellingham

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons. org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Twelve Years of High Frequency Absolute Gravity Measurements at the UK's Space Geodesy Facility: Systematic Signals and Comparison with SLR Heights**

## Victoria Anne Smith, Graham Appleby, Marek Ziebart, and Jose Rodriguez

#### **Abstract**

Absolute gravity measurements taken on a near-weekly basis at a single location is a rarity. Twelve years of data at the UK's Space Geodesy Facility (SGF) provides evidence to show that the application of results from international comparisons of absolute gravimeters should be applied to data and are critical to the interpretation of theSGF gravity time series of data from 2007 to 2019. Though residual biases in the data are seen. The SGF time series comprises near weekly data, with exceptions for manufacturer services and participation in international instrument comparisons. Each data set comprises hourly data taken over 1 day, with between 100 and 200 drops per hour. Environmental modelling indicates that the annual groundwater variation at SGF of some 2 m influences the gravity data by 3.1 -Gal, based upon some measured and estimated soil parameters. The soil parameters were also used in the calculation of the effect of an additional telescope dome, built above the gravity laboratory, and have been shown to be realistic. Sited in close proximity to the long-established satellite laser ranging (SLR) system and the global navigation satellite systems (GNSS) the absolute gravimetry (AG) measurements provide a complimentary geodetic technique, which is non space-based. The SLR-derived height time series provides an independent measurement of vertical motion at the site which may be used to assess the AG results, which are impacted by ground motion as well as mass changes above and below the instruments.

#### **Keywords**

Absolute gravimeter - Gravimetry - International Terrestrial Reference Frame (ITRF) - Satellite Laser Ranging (SLR) -Time series

V. A. Smith (-)

NERC/BGS Space Geodesy Facility, Herstmonceux Castle, East Sussex, UK e-mail: vism@nerc.ac.uk

G. Appleby British Geological Survey Honorary Research Associate, Nottingham, UK e-mail: gapp@nerc.ac.uk

#### M. Ziebart

Department of Civil, Environmental and Geomatic Engineering, Faculty of Engineering Science, University College London, London, UK e-mail: m.ziebart@ucl.ac.uk

#### **1 Introduction**

Characterisation of key geodetic sites for the improvement of their products is important for the demands of geodesy and the various reference frames which are built upon the data from the worldwide geodetic observatories. The addition of absolute gravity (AG) measurement capabilities at the Space Geodesy Facility (SGF) was prompted by the growing demands of the Global Geodetic Observing System (GGOS)

J. Rodriguez

Centro Astronómico de Yebes, Guadalajara, Spain e-mail: jc.rodriguez@oan.es

and the European Combined Geodetic Network (ECGN). An early goal for the observatory was to determine if the geodetic products from the long established satellite laser ranging (SLR) and GNSS, could be improved through better characterisation of any un-modelled or mis-modelled signals by a non-space based AG technique. Determination of the glacial isostatic rate for the southeast of the UK was an additional long term goal.

In this paper, the results of 12 years of near-weekly AG data are presented. The environmental effects and calculations are briefly discussed. Emphasis is placed on the interpretation of the time series as a whole, where instrumentational inter-comparisons, data offsets and bias offsets are discussed and analysed. Finally, a comparison is made between the AG time series and station heights derived from SLR.

#### **2 Site Information**

The UK's Space Geodesy Facility is funded by the UK Natural Environment Research Council through the British Geological Survey. It is located in East Sussex, less than 5 km from the English Channel. The geodetic techniques of SLR, GNSS and absolute gravity measurements form the principle operations at the observatory, though numerous additional sensor data, including sun photometry and automated measurements of groundwater depth, atmospheric visibility, pressure, temperature and humidity, are recorded. The site is compact with each geodetic technique contained within 25 m. An additional GNSS antenna is located around 100 m distant. The gravimetry laboratory is located SE from the SLR, three meters below ground level as measured at the borehole, in a semi-sunken bunker that has approximately 1.2 m of soil above it. Unusually for geodetic sites, the subsurface is comprised of clay, with no bedrock beneath within 30 m. The local soil type is known to be Weald clay (British Geological Survey website1). Absolute gravity measurements have been taken at the observatory since 2006 on a near-weekly basis. The standard measurements comprise of hourly data sets with 100–200 drops per set taken for 25 h once a week.

#### **3 Hydrology**

The hydrology surrounding the gravity laboratory is slightly complex due both to its semi-sunken nature and the slope on which the SGF sits (Fig. 1). It has been well documented that water movement should be accounted for when analysing gravity data (Makinen and Tattari 1990; Harnisch

**Fig. 1** Schematic to show the semi-sunken nature of the gravity laboratory

and Harnisch 2006) groundwater, soil moisture and rainfall all change the mass around gravimeters and therefore, it should be understood in order to interpret the data correctly. Groundwater data have been recorded at Herstmonceux since the mid 1990s by an automatic data logging system in a borehole, measurements are recorded as depth below ground level. There is currently no capability for soil moisture to be recorded.

However, using estimates made by testing soil samples, approximations of the influence of these hydrologicallyinduced signals on the AG data have been made. The density, porosity and specific yield were needed to estimate the differences between wet and dry soils both above and below the gravimeter. The modified Bouguer plate corrections for soil moisture and ground water are:

$$
\delta \mathbf{g}\_{sm} = 2\pi G H \rho\_{\text{w}} \delta P = 4.2 H \delta P \tag{1}
$$

$$
\delta \mathbf{g}\_{\mathcal{g}\mathbf{w}} = 2\pi GP \rho\_{\mathbf{w}} \delta H = 4.2 P \delta H \tag{2}
$$

Where <sup>ı</sup>*gsm* indicates the effect on gravity given by a change in soil moisture, <sup>ı</sup>*ggw* indicates the effect on gravity given by a change in groundwater depth, G is the gravitational constant, *<sup>w</sup>* is the density of water, <sup>ı</sup>*<sup>P</sup>* is the change in the water filled pore spaces in the soil and ı*<sup>H</sup>* is effective change in depth to groundwater.

The effect on gravity, due to the variation of water content in the 1.2 m of soil above the gravimeter, was calculated to be less than 1 -Gal. The influence on the gravity data due to seasonalvariation in the groundwater height, of between10 and 12 m below the gravity laboratory, has been calculated

<sup>1</sup>http://mapapps.bgs.ac.uk/geologyofbritain/home.html?.

to be 3.1 -Gal at maximum. Application of the groundwater correction gives an almost imperceptible change to the appearance of the AG time series. Full details of these calculations can be found in the PhD thesis by Smith (2018).<sup>2</sup>

#### **4 Interpretation of the Time Series**

The weekly gravity data as presented in Fig. 2 clearly contain some interesting features. However, these features become of particular significance when important events such as changes to the instrumentation or environment are highlighted, as shown in Fig. 2. These events appear to cause clear offsets in the time series of AG data. If the data are split into slices, dictated by each event, then each resulting 'regime' of data, as shown in Fig. 3, can be analysed and the differences between them can be studied.

It is interesting to note that the data in each regime falls within a 10 -Gal range but are offset in mean value. Also, the standard deviation about the mean of each regime is similar: 1.17, 1.31, 1.39, 1.56 and 2.28 -Gal.

It should also be noted that the precision of the AG measurements has been decreasing for several years, which is thought to be driven by the environment, as comparison data has not indicated any problem with the instrument.

<sup>2</sup>Available from the British Library or UCL Discovery website.

If the data are to be interpreted as a whole, it is clear that the offset changes, coincident with instrumental visits to the manufacturers, need to be accounted for. Since the measurement of gravity is dependent on instrumentation, as well as spatial and temporal variables, no estimation of the accuracy of the instruments is possible. To quantify results of absolute gravimeters the best solution is to obtain relative offsets between instruments during international gravimeter comparisons (Francis et al. 2013). Relative offset values obtained from comparison results were obtained three times during this time series for the FG5 used (FG5-229). The comparisons results for 2007 and 2015 gave offset results for the SGF instrument of <sup>C</sup>1.5 ˙ 0.8and 0.08 ˙ 0.8 -Gal respectively (Francis and van Dam 2010; Pálinkáš et al. 2017). A basic mini-comparison was carried out in 2012, courtesy of O. Francis, at the Walferdange Underground Laboratory, with another FG5 (FG5-216) which itself was found to have an offset of <sup>C</sup>1.8 ˙ 3.1 -Gal during the international comparison of 2011 (Francis et al. 2013). Since FG5-229 was found to be in very close agreement (0.2 ˙ <sup>1</sup> -Gal) with FG5-216, the offset of C1.8 has been used for FG5- 229. Unfortunately, since there is no supporting comparison for the early data, before the SIM (system interface module, through which the majority of signals are passed and contains control electronics) repair in 2006, the data from the first regime shown in Figs. 2 and 3 have been discounted from further analysis at this time.

The offset corrections from the comparisons have been applied and the resulting time series is shown in Fig. 4. The data from 2007 to 2013 have been significantly smoothed and now appears to give consistent data. However, there are significant discontinuities remaining thereafter.

The offset in the data from 2016 onwards is of known origin; a small telescope dome was built above the gravity laboratory. An approximate calculation of mass change, based upon volume of clay extracted from above the gravity laboratory and the addition of concrete and brick, implies an offset in the gravity measurement of 3.6 -Gal. This is in reasonable agreement with the observed mean offset as discussed below. However, the data from 2013 are concerning. They have a large residual offset after the comparison results are applied. We have no evidence to support any local environment changes that could have accounted for this level of change over the 3 months that the instrument was away at the manufacturers. Although heavily researched, this bias remains of unknown origin.

Several methods could be employed to determine the biases in the 2013–2016 and 2016–2019 data sets. However, using the simplest method, calculating the mean values of each regime, proves to be an acceptable method in this case, since the difference between the mean values of each regime compares favourably with both the offsets obtained from the 2007 and 2011 comparisons and the estimated effect of mass changes due to the additional telescope dome. The comparisons in question gave a total offset of 3.5 -Gal, while the difference between the mean values is 3.6 -Gal. The calculated change due to the additional dome built above the laboratory, based upon soil samples, indicated an expected offset of 3.6 -Gal, whereas the difference between the mean values gives 4.8 -Gal. It should have noted that the dome calculations did not account for the additional masses from the telescope, telescope mount, pier or dome, which would all increase the expected offset. Furthermore, in earlier work, the magnitude of the bias in 2013 was taken to be 7.3 -Gal

**Fig. 6** Campaign simulation, uncorrected for bias, dome or comparisons

obtained by 'best-fit', whilst the difference between the mean values gives an offset close to this of 7.2 -Gal. When these offsets are applied, in addition to the offsets provided by comparisons of instruments, the results are shown in Fig. 5.

#### **5 Campaign Simulation**

To emulate low frequency FG5 data a 'campaign style' simulation was applied to the SGF time series; by taking one data point per year from the full SGF time series. Project data, comprising 25 h of data, were selected around the same epoch each year and plotted in Fig. 6. Interpretation of these results could be critical for a country-wide project: it could be tempting to apply an offset for the 2013 to 2017 data and estimate a glacial isostatic adjustment (GIA) rate from the residual. Such an interpretation would be critically flawed given our understanding of events in the whole SGF time series.

Bias and comparison corrected results for this campaign simulation are given in Fig. 7.

#### **5.1 AG vs SLR Heights**

The SGF is an International Laser Ranging Service (ILRS) Analysis Centre (AC) and is responsible for computing reference frame solutions that are subsequently made freely

**Fig. 8** AG data, converted to heights, plotted in red, with SLR derived station height changes plotted in green

available to the community via ILRS Data Centres.<sup>3</sup> The SGF AC contributed multi-year solutions towards the realisation of the International Terrestrial Reference Frame (e.g., ITRF2014, (Altamimi et al. 2016)) using global SLR measurements to the LAGEOS and Etalon satellites. In addition, research carried out by the AC, (Appleby et al. 2016; Rodríguez et al. 2019) after the work that contributed to ITRF2014, has identified subtle systematic effects in the range measurements at each station of the ILRS network that impact in particular the derived station heights at the few mm to 10 mm level. As a consequence of this work, new solutions for station coordinates that take account of these small systematics have been computed and may be considered close-to bias free. The resulting weekly-averaged heights of the Herstmonceux site thus provide a reference height series against which to measure the stability and potential systematics in the AG gravity observations; the AG measurements will be impacted by height changes, as tracked by the SLR solutions, as well as mass changes above and below the AG instrument, primarily hydrological changes as discussed in Sect. 3. Therefore, the data from the two techniques are not expected to match precisely. The SLR and AG time series (converted to height by the use of a multiplication factor of 5 mm/-Gal (Teferle 2009)) are shown in Fig. 8, where the data has been aligned in the vertical axis by matching the first data point of both the AG and SLR data. It is very interesting and promising to note that the AG series provides a smoother representation of the

<sup>3</sup>https://ilrs.cddis.eosdis.nasa.gov/data\_and\_products/data\_centers/ index.html.

height variation of the site than is determined by the SLR time series.

#### **6 Conclusion**

The high frequency, near weekly, AG data from SGF show clear evidence of the importance of testing the instruments after they have been serviced and after changes to the instrument have been carried out by the manufacturer. The implementation of the comparison results from international meetings is seen to be critical, though not the complete story for the bias offset seen in 2013. It is alarming that this bias coincides with the major instrument change of the upgrade of the FG5 to an FG5X. The results indicate that instrument comparisons should be carried out as often as is practical for all users of FG5s.

SLR-derived site height variations are used to validate the concept of applying the corrections made to the AG time series that results from inter-comparisons, bias estimation and new dome building. SGF is currently testing an additional site in the UK with an aim to be able to offer a minicomparison site for the community.

#### **References**


scale: estimation of systematic errors in LAGEOS observations 1993–2014. J Geod 90(12):1371–1388, 1*–*18. https://doi.org/10. 1007/s00190-016-0929-2


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/ licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

**Part III**

**Relative Gravimetry and Applications**

# **Studying the Evolution of Resolution Capabilities and Approximation Accuracy of Global Models by Spectral Characteristics**

Irina G. Ganagina, Alexander P. Karpik, Vadim F. Kanushin, Denis N. Goldobin, and Elena M. Mazurova

#### **Abstract**

The paper presents research of the evolution of resolution capabilities and approximation accuracy of the last decade global geopotential models based on space gravimetric measurements of CHAMP, GRACE and GOCE missions using their spectral characteristics.

The comparison between the model data with point measurements data on gravity anomalies and quasigeoid heights for the territory of Novosibirsk region is shown. Based on the research results the conclusion was drawn that accuracy characteristics of current global models under test built by the results of satellite gravimetry missions do not achieve the specified accuracy of 1 cm and 1 mGal on the territory under investigation. The research has made it possible to state that at the current technological and methodical level the potential has been reached as concerns EGF models resolution and accuracy enhancement.

#### **Keywords**

Current global geopotential model - Resolution - Approximation accuracy - Spectral characteristics - Degree dispersion -Territory of West Siberia

#### **1 Introduction**

The study of the resolution and approximation accuracy of the global gravity field (EGF) models is a vital task as regards fine structure of the gravitational field. Current models in the form of standardized coefficients of geopotential spherical harmonics make it possible to construct detailed high accuracy digital models of gravity fields characteristics (Koneshov et al. 2013; Aka Blush Ulfred 2019; Erol et al. 2020; Abd-Elmotaal et al. 2020). Reliable estimate of EGF models accuracy increases opportunity of their application in geodynamics, geophysics and navigation.

The object of this paper is the study of evolution of resolution and approximation accuracy for global models by their spectral characteristics.

Estimated approximation accuracy criterion for quasigeoid height and gravity anomalies implies the values of 1 cm and 1 mGal respectively, as specified by project GOCE developers.

#### **1.1 Theoretical Part**

When modelling the Earth's gravity field, the attraction potential expansion in a Fourier series is used by the system of spherical harmonics of geocentric coordinates—radius vector *r*, latitude -, and longitude in the form (Hoffman-Wellenhof and Moritz 2005)

$$V(\boldsymbol{\varphi}, \lambda, r) = \frac{fM}{r} \left[ 1 + \sum\_{n=2}^{\infty} \left( \frac{a\_{\ell}}{r} \right)^{n} \sum\_{m=0}^{n} \left( \overline{C}\_{nm} \cos m\lambda \right) \right] \tag{1}$$

$$+ \overline{S}\_{nm} \sin m\lambda \left| \overline{P}\_{nm} \left( \sin \boldsymbol{\varphi} \right) \right| \tag{1}$$

I. G. Ganagina (-) · A. P. Karpik · V. F. Kanushin · D. N. Goldobin Siberian State University of Geosystems and Technologies, Novosibirsk, Russia

e-mail: gam0209@yandex.ru; kaf.asnronomy@ssga.ru; kaf. asnronomy@ssga.ru; kaf.asnronomy@ssga.ru

E. M. Mazurova

Federal Scientific-Technical Center of Geodesy, Cartography and Spatial Data Infrastructure, Moscow, Russia e-mail: e\_mazurova@mail.ru

J. T. Freymueller, L. Sánchez (eds.), *5th Symposium on Terrestrial Gravimetry: Static and Mobile Measurements (TG-SMM 2019)*, International Association of Geodesy Symposia 153, https://doi.org/10.1007/1345\_2021\_137

where *fM*—geocentric gravitational constant;

*a <sup>e</sup>*—Earth's equatorial radius; Cnm Ë Snm—normalized dimensionless harmonic factors of geopotential;

Pnm .sin '/—normalized associated Legendre polynomials.

Summing-up in formula (1) is performed to infinity, and the Earth's gravity field models are limited by the maximum degree *N*0. Series (1), limited by maximum *N*<sup>0</sup> presents spectral expansion of the Earth's gravity field structure by waves lengths 360ı/2*N*<sup>0</sup> that corresponds to the spatial angular resolution.

$$
\Delta = \frac{180^{\circ}}{N\_0}.\tag{2}
$$

Simulation error of the global Earth's gravity field *V*(®*, , r*) by Fourier series *VN* (®*, , r*) in point P(®*, , r*)2¨ with approximate harmonic coefficients ˚ Cnm; Snm is written as

$$
\Delta\_N = \rho(V, V\_N) + \max\_{P \in \alpha \nu} |\: \delta N \: \mid, \tag{3}
$$

where

$$\rho\left(V, V\_N\right) = \max\_{P \in \alpha} |\left|\left.V(P) - V\_N(P)\right|\right|;\tag{4}$$

ı*N*—simulation error of field *V*(*P*) due to the error of harmonic factors fı*Cnm*, ı*Snm*g.

Geopotential approximation dispersion error *V*(®*, , r*), taking into account the error of harmonic factors, is estimated by formula (Kanushin et al. 2015)

$$D\left\{\Delta\_{N}\right\} = D - \sum\_{n=0}^{N} \left(D\_{n} - d\_{n}\right) \tag{5}$$

where *D* is a dispersion of the initial potential *V* with unlimited spectrum. The degree dispersion of geopotential is calculated by formula

$$D\_n = \left(\frac{a\_\epsilon}{r}\right)^{2n} \sum\_{m=0}^n \left(\overline{C}\_{nm}^2 + \overline{S}\_{nm}^2\right). \tag{6}$$

Formula (7) is used to calculate the degree dispersion of errors ıC<sup>2</sup> nm and ıS<sup>2</sup> nm of geopotential harmonic factors *V*

$$d\_n = \left(\frac{a\_e}{r}\right)^{2n} \sum\_{m=0}^{n} \left(\delta \overline{C}\_{nm}^2 + \delta \overline{S}\_{nm}^2\right). \tag{7}$$

Dependence of approximation error " of gravity anomalies and quasigeoid heights due to series (1) limit in the form of

$$\begin{Bmatrix} \varepsilon\_{\Delta g\_N}^2 \\ \varepsilon\_{\xi\_N}^2 \end{Bmatrix} = \begin{Bmatrix} D\_{\Delta g} \\ D\_{\xi} \end{Bmatrix} - \sum\_{n=2}^N \left( \begin{Bmatrix} D\_{\Delta g\_n} \\ D\_{\xi\_n} \end{Bmatrix} - \begin{Bmatrix} d\_{\Delta g\_n} \\ d\_{\xi\_n} \end{Bmatrix} \right), \tag{8}$$

spherical degree dispersions of gravity anomalies and quasigeoid heights expressed by factors Cnm and Snm

$$\left\{ \begin{array}{c} D\_{\Delta g\_{n}} \\ D\_{\zeta\_{n}} \end{array} \right\} = \left\{ \begin{array}{c} \chi^{2}(n-1)^{2} \\ R^{2} \left(\frac{\mu\_{\omega}}{r}\right)^{2n} \end{array} \right\} \sum\_{m=0}^{n} \left( \Delta \overline{C}\_{nm}^{2} + \overline{S}\_{nm}^{2} \right), \quad (9)$$

where Cnm <sup>D</sup> Cnm <sup>C</sup><sup>0</sup> nm—the difference between coefficients of normalized spherical functions of real and normal gravity fields;

C0 nm—coefficients of normal geopotential are referred to ellipsoid WGS- 84.

The degree dispersion of the errors of gravity anomalies and quasigeoid heights coefficients is written as

$$\left\{ \begin{array}{c} d\_{\Delta g\_{n}} \\ d\_{\zeta\_{n}} \end{array} \right\} = \left\{ \begin{array}{c} \nu^{2} (n-1)^{2} \\ R^{2} \left( \frac{a\_{e}}{r} \right)^{2n} \end{array} \right\} \sum\_{m=0}^{n} \left( \delta \overline{C}\_{nm}^{2} + \delta \overline{S}\_{nm}^{2} \right). \quad (10)$$

**Experimental.** The study of evolution of the resolution and accuracy for 70 models under test for global Earth's gravity field was conducted by the comparison results of geopotential harmonic coefficients degree dispersions and their errors as well as the computed (by them) degree dispersions of gravity anomalies and quasigeoid heights.

For satellite models under test their resolution and accuracy results according to their spectral characteristics are presented in the time range of 2008–2015 (Fig. 1). At the present moment the maximum expansion degree of the satellite models has achieved N0 D 300. The given expansion degree is limitary for satellite models built by current space gravimetry missions.

The evolution of the spatial resolution of the models built by the last decade satellite data shows that the gravity anomalies models resolution averages to 100 km, and for quasigeoid heights 200 km.

The results of global models resolution obtained by the combined data with error approximation 1 cm and 1 mGal are presented for the period of 1996–2019 (Fig. 2). The degree dispersions for the combined models under test practically coincide to the 240th degree. This testifies to the fact that the low-frequency harmonic segment of these geopotential models has been thoroughly studied.

The data on the evolution of resolution capability for the gravity anomalies and quasigeoid heights models obtained by the combined data shows that the resolution capability was changing insignificantly during the last 15 years, i.e.

**Fig. 1** Resolution capabilities for spectral characteristics obtained according to low-level models, when approximation errors of 1 cm and 1 mGal are achieved (period 2008–2015)

about 180 km for quasigeoid and about 70 km for gravity anomalies.

The results of spectral estimates of global geopotential ultrahigh-degree models show that their degree dispersions practically coincide over all the frequency range under study. This means that the same initial data were mostly used for constructing these models. Significant improvements for the resolution of gravity field ultrahigh-degree models (Fig. 3)

**Fig. 2** The results of calculating the resolution capabilities from the spectral characteristics obtained from the combined data, with approximation error of 1 cm and 1 mGal (period 1996–2019)

built in 2008–2019 have not been observed. For the gravity anomalies models the resolution is about 45 km and that for quasigeoid height—about 170 km.

Analysis of the obtained research results as regards evolution of gravity field models resolution capability (for combined ultrahigh-degree models made in 2008–2019) revealed that it has not practically changed. Resolution of satellite models (in the period of 2008–2019) changed for the gravity anomalies model from 200 km to 100 km; for the quasigeoid model it changed insignificantly, i.e. from 230 km to 200 km.

Analyzing the spectral estimator changes from one geopotentional model to another a clear view may be formed of the resolution capability of certain harmonics or groups of harmonics of these models, and of the uncertainty state

of each parameter under consideration due to the model determination features.

The comparison of model data with land point data on gravity anomalies and quasigeoid heights on the territory of Novosibirsk region is presented.

On the territory under study 17 second-order gravity base network stations were chosen with gravity determination accuracy being ˙ 0.05 mGal. The territory under investigation includes 208 stations where geometric leveling (first to fourth order) was conducted and the values of normal heights were obtained. At the same points, satellite coordinates were determined for developing active base stations networks. As a result of satellite network adjustment geodetic heights were obtained, with mean square errors ranging from 1.5 cm to

**Fig. 3** The results of calculating the resolution capabilities by the geopotential ultrahigh-degree models characteristics with approximation error of 1 cm and 1 mGal (period 2008–2019)

3.1 cm, 1.8 cm, on average (Karpik et al. 2010). On the territory under consideration gravity anomalies variations and elevation differences are insignificant.

The scheme of differences between the quasigeoid heights obtained from model EIGEN-6C4 and those of geoimetric leveling and GNSS measurements is given in Fig. 4. Applying model EIGEN-6C4 made it possible to achieve better results on the territory under consideration, with maximum difference being 22 cm and mean square error about 8 cm.

The scheme of differences between the gravity anomalies obtained by EIGEN-6C4 model and those based on the ground data is shown in Fig. 5.

Maximum difference amounted to 8mGal, with mean square error being about 4 mGal.

The comparison of model and ground data for the gravity and quasigeoid heights anomalies on the territory of Novosibirsk region (for 9 geopotential models) is given in Table 1.

#### **2 Conclusions. Interpretation**

Comparison of degree dispersions of geopotential harmonic coefficients and their errors is presented for 70 current global models built by the results of measurements of space gravimetric missions CHAMP, GRACE and GOCE.

**Fig. 4** Scheme of differences between the heights of the quasigeoid obtained from the EIGEN-6C4 model and those obtained from geometric leveling and GNSS measurements, m

The models under study revealed: the dependences of gravity anomalies and quasigeoid heights approximation errors on the degree of geopotential expansion into Fourier transform; ratio errors of geopotential harmonic coefficients determination; coefficients dispersions and those of gravity anomalies coefficient errors.

Satellite models of the Earth's gravity field built by only space gravimetry missions are not sufficient for obtaining quasigeoid heights with accuracy of 1–2 cm and spatial resolution less than 100 km. Satellite models are to be used in combination with land and topographic data to reestablish the global quasigeoid.

Models resolution at the satellite altitude level satisfies the requirements on the models accuracy stated by the developers. In case of necessity for recomputation of the gravity field characteristics for the Earth's surface the specified resolution capability is reduced.

The compared values of the model and ground data on the Novosibirsk region areas demonstrated that the least values of differences were obtained on the global geopotential models based on the GOCE mission data (Kanushin et al. 2015). The global models accuracy values for the territory of Novosibirsk region do not come up to 1 mGal and 1 cm, being in the range of 4–10 mGal and 7–20 cm, respectfully.

The obtained accuracy values of the gravity anomalies and quasigeoid heights in current global geopotential models under study may be considered as maximum for the newly constructed models by the data of GOCE project.

Space gravity mission GOCE made it possible to essentially improve harmonic coefficients for degrees 100–250.

Harmonic coefficients of geopotential obtained for current EGF models are in conformity with each other within mean square errors.

Further research of the global models is required using the results of space gravity mission GOCE as the most successful mission for investigating EGF.

Accuracy characteristics of the models under study approximate the specified ones in process of their construction, but fall short of them as concerns the required resolution capability.

Current global models of gravity field create new opportunities and essentially extend the range of challenges in studying the Earth's gravity field and its figure.

Comparing the evolution of resolution capability and approximation accuracy as regards current global models of the Earth's gravity field we may come to the conclusion that at the current stage of technological and methodical level, EGF model resolution capability and accuracy improvement has reached its maximum. For further improvement of EGF approximation accuracy new technological approaches and methodological principles should be developed. This implies that a new concept of the Earth's gravity field study should be developed.

**Fig. 5** Scheme of differences between the EIGEN-6C4 model reconstructed and ground values of gravity anomalies obtained for the Novosibirsk region in mGal



#### **References**

Abd-Elmotaal HA, Kühtreiber N, Seitz K, Heck B (2020) A precise geoid model for Africa: AFRgeo2019, International Association of Geodesy Symposia, Springer, Berlin, https://doi.org/10.1007/1345\_ 2020\_122


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/ licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **The Role of Non-tidal Atmospheric Loading in the Task of Gravity Field Estimation by Inter-Satellite Measurements**

### I. O. Skakun, V. V. Mitrikas, and V. V. Ianishevskii

#### **Abstract**

The paper reviews models of tidal and non-tidal variations of the Earth's gravitational field. Proposing an algorithm for the estimation of the Stokes coefficients based on inter-satellite measurements of low-orbit spacecrafts. By processing measurements of the GRACE mission, we obtained experimental estimates of gravity field monthly variations. The analysis of these values was carried out by calculating the change in the equivalent water height for a given area.

#### **Keywords**

GRACE - Gravity field -Atmospheric tides

#### **1 Introduction**

In 2002 satellites of the GRACE mission (Greicius 2013) were launched into near-earth orbit, allowing the mapping of the Earth's gravitational field at monthly intervals. The spatial resolution of derived models is several hundred kilometers, determining the height of a quasigeoid with an accuracy of a few centimeters.

The basis for calculating the monthly models of the gravity field is the joint processing of inter-satellite measurements and data on the position of spacecraft obtained from measurements of the navigation receiver (Beutler et al. 2010). The procedure consists of three stages:


As a result of the estimation, independent monthly models of geopotential are obtained.

As is known, in addition to the tidal and relativistic effects corresponding to the IERS 2010 convention (Petit and Luzum 2010), the processing of spacecraft in low orbits takes into account the influence of non-tidal atmospheric loading calculated according to the NWP (Numerical Weather Prediction model) of the European Center for Medium-Term weather forecasts (ECMWF – European Center for Medium-Range Weather Forecasts) at the German Research Center for Geoscience (GFZ German Research Center for Geoscience) (Dobslaw et al. 2016).

In current work we compare the first results of in-house estimation of monthly gravity field models based on satelliteto-satellite measurements with other solutions, and also study the contribution of non-tidal atmospheric loading using a specially developed experimental programming and mathematical software KAGOR (KArtography Geodesy and Orbit determination). KAGOR allows simulating orbits, kinematic trajectories and corresponding inter-satellite measurements with any given parameters, as well as to process real satellite-to-satellite measurements to estimate gravity field models.

e-mail: Ivan.skakun@glonass-iac.ru

© The Author(s) 2021

J. T. Freymueller, L. Sánchez (eds.), *5th Symposium on Terrestrial Gravimetry: Static and Mobile Measurements (TG-SMM 2019)*, International Association of Geodesy Symposia 153, https://doi.org/10.1007/1345\_2021\_134

I. O. Skakun (-) · V. V. Mitrikas · V. V. Ianishevskii JSC "TSNIImash", Korolyov, Russia

#### **2 Non-tidal Atmospheric Loading**

Based on the analysis and forecasting of the data from the operational high-resolution global model of the numerical weather forecast of the European Center for Medium-Range Weather Forecasts and bottom pressure data obtained by simulation without imposing additional restrictions on the global general circulation model MPIOM and ECMWF data, a product of non-tidal gravitational variations in the atmosphere and oceans is formed (Atmosphere and Ocean De-Aliasing Level-1B, AOD1B) which provides a priori information on temporal variations in the gravity field due to global mass redistribution in the atmosphere and the ocean (Dobslaw et al. 2017).

AOD1B was determined by the fully normalized Stokes coefficients of the anomalous external gravitational field of the Earth caused by the redistribution of masses predicted by the numerical models described above.

AOD1B consists of four sets of coefficients:


This product is set with a resolution of 3 h. Since January 1, 1976, atmospheric anomalies and bottom pressure anomalies have been arranged in a series of spherical functions up to the 100 d/o to ensure the processing of space geodetic missions equipped with laser retroreflectors, in particular, LAGEOS. Since January 1, 2000 AOD1B is set by decomposition to a degree and order of 180 and, in addition to atmospheric anomalies and anomalies in bottom pressure, takes into account density anomalies of the upper atmosphere, which allows to meet the high requirements of modern space geodetic missions aimed at the study of gravity field (CHAMP, GRACE, GOCE).

The influence of the atmosphere consists of two parts: non-tidal and tidal. Depending on the direct gravitational attraction of atmospheric masses by bodies generating tidal potential and, more significantly, on the periodic motion of the lower boundary of the atmosphere caused by the tide in the Earth's solid body and oceans and on the absorption of solar radiation by the atmosphere, which causes a signal with a period close to 24 h. The tidal part includes 12 waves (P1, S1, K1, N2, M2, L2, T2, S2, R2, T3, S3, R3) and is calculated by updating the coefficients of the corresponding waves of the tidal potential in the observation interval 2007– 2014. A non-tidal atmospheric load is formed by subtracting the obtained waves from the measurements, the study of the influence of which this work is devoted to.

An example of non-tidal atmospheric load for the 1st January, 2003, 00:00 in the form of an equivalent water height (EWH) Using a Gaussian filter (300 km) (Wahr 2007) is shown in Fig. 1.

#### **3 Estimation of Monthly Gravity Filed Models**

As previously noted, the gravity field monthly model estimation procedure consists of three stages:


Independent monthly models of geopotential are the result of an estimation.

Any processing is based on a measurement model. In general, it has the following form:

$$\mathbf{y}(t) = f\left(t, r\_0, v\_0, C\_{nm}, S\_{nm}, p, \dots\right) + \epsilon,\qquad(1)$$

where *y*(*t*) measurements, *t* time, *r*0, *v*<sup>0</sup> initial conditions in inertial reference frame, *Cnm*, *Snm* Stokes coefficients of degree and order less or equal then n and m, *p* instrumental parameters, – measurement error

At the first stage, the kinematic trajectory of the spacecraft (the position of the spacecraft for all navigation receiver measurements epochs on the processing interval calculated with PPP technique) was used as measurements. In the second, besides the trajectories of both spacecraft, the relative velocity was added.

At the first stage, we estimated the initial conditions of the spacecraft and the parameters of the accelerometer, while the Stokes coefficients was not estimated. For the calculation of the spacecraft coordinate system, star camera data was used. The orbit parametrization is generally similar to that used in the work devoted to the analysis of high-precision methods of ballistic-navigational support of space geodetic complexes using the Geo-IK2 spacecraft as an example (Zaliznyuk et al. 2019). At the second stage, Stokes coefficients up to the 80th degree and order and model empirical parameters of the measurement errors of the inter-satellite line were added

**Fig. 1** Non-tidal atmospheric load in the form of equivalent water level for the 1st January, 2003 00:00

to the refined parameters. At both stages, the relationship between the dependent variables and the independent ones is nonlinear, so linearization is performed in the vicinity of the a priori orbit. The a priori orbit was calculated by numerical integration of the system of equations of motion with the initial conditions, taking into account a number of perturbations in accordance with the IERS 2010 convention. Table 1 shows the perturbations taken into account. At both stages, the derivatives of the measurements with respect to the orbital parameters were calculated by numerical integration of the variations equations.

**Table 1** The perturbations taken into account when integrating the system of equations of motion


The contribution of various gravitational effects and a typical change in gravity filed over a monthly interval are shown in Fig. 2 in the values of the root of the degree variance in the logarithmic scales. Obviously, the static gravitational field has the greatest influence on the motion of spacecraft in the entire gravity spectrum. Then comes the ocean tide, whose contribution to the low-wave part of the spectrum is smaller than the changes in the gravity at monthly intervals. The atmospheric non-tidal load, for its part, has a greater effect on the gravity filed than its monthly variability up to 10–15 harmonics and less on higher harmonics.

Also, notice that when integrating the equations of motion, in addition to the conventional perturbations, measurements of non-conservative forces made by accelerometers installed on board both GRACE spacecraft were taken into account. Because of that, we need to add to the vector of the estimated parameters a systematic error for each of the three axes of the accelerometer and scaling factors in all of them as well.

At the third stage, from the matrices of conditional equations formed at the second stage corresponding to 3-h intervals, the matrix of normal equations was formed on a monthly interval. The vector of specified parameters of the monthly interval is a superposition of 3-h vectors in terms of orbital and instrumental parameters, and includes one set of

Stokes coefficients for the entire interval. On average, one monthly solution includes about 3,600,000 measurements, 10,000 estimated parameters, takes about 200 GB on the hard drive and lasts 3 days on PC.

#### **4 Results**

To assess the contribution of non-tidal atmospheric load to the estimation of monthly gravity field models, two series of solutions were carried out: with and without AOD1B. In Fig. 3 is a diagram of the solutions made to estimate the monthly gravity models. The missing of 2 months of both sets of solutions is due to the lack of accelerometer data at these intervals. This interval was chosen because of the relatively more successful functioning of the mission as a whole (a relatively small number of hardware failures). Considering the informativeness of the results of comparing the obtained time series of monthly gravity models and the high cost of computational resources for calculations, solutions without the use of non-tidal atmospheric loading for 2005 were not implemented.

To evaluate the characteristics of both obtained time series of the gravity models, the equivalent water height was calculated using a Gaussian filter (300 km) in the Greenland region, since for this region a pronounced annual signal of changes in the gravity field and a long-period trend are observed. A rectangle was used (60..85 deg. of latitude; 70..20 deg. of longitude). For the entire 3 year interval, the average value of gravitational potential was found, the pairwise difference of all monthly models with the obtained average was constructed, then for each point in the region with an increment of 1<sup>ı</sup> the equivalent water level of the difference of the monthly model with the average was considered and the average value for the entire region was considered. The values obtained for both series of gravity models and similar values calculated from models of other official GRACE measurement processing

**Fig. 3** Diagram of monthly gravity field solutions

**Fig. 2** The contribution of various gravitational effects

**Fig. 4** The equivalent water level in Greenland (60..85 deg. of latitude; -70..-20 deg. of longitude), calculated and averaged over the grid with a step of 1 degree using a Gaussian filter (300 km) over a 3-

centers (CSR, GFZ, JPL) are shown in Fig. 4. Considering that in the C20 series obtained from GRACE measurements, an unknown 161-day period is observed, the values of the C20 coefficient were taken from the results of processing five satellites (LAGEOS-1/2, Stella, Starlette, AJISAI) laser location data (Satellite Laser Ranging, SLR) (Loomis et al. 2019).

The results of the three official GRACE data centers are in good agreement with each other, due to the similarity of the models and processing algorithms used. The solutions we obtained taking into account the non-tidal atmospheric loading agree much better than the solutions obtained without taking into account AOD1B. It can be seen that ignoring the non-tidal loading in processing introduces a distortion in the estimation of the gravity field models the same size as the signal of its change. This is the first result for the KAGOR software and this software has potential capability to increase the accuracy of the solutions obtained, in particular, in outlier search algorithms and reducing the influence of anomalous measurements. Nevertheless, this SW already provides an account of all the necessary effects and allows to obtain an adequate estimate of the gravity field signals.

#### **5 Conclusions**

1. The KAGOR (cartography, geodesic research and orbit determination) experimental software has been developed, which allows to model orbits, kinematic trajectories

year interval according to monthly models of gravity filed by different analysis centers CSR, GFZ, JPL and IAC

and to take corresponding satellite-to-satellite measurements with any given parameters, as well as to process real satellite-to-satellite measurements to estimate gravity field models.


#### **References**


Loomis BD, Rachlin KE, Luthcke SB (2019) Improved earth oblateness rate reveals increased ice sheet losses and mass-driven sea level rise. Geophys Res Lett 46(12):6910–6917

Petit G, Luzum B (2010) IERS Technical Note; No. 36, vol. 36, p 179 Wahr JM (2007) 3.08 - time variable gravity from satellites. In: Schubert


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/ licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Formation of Self-consistent Navigational Gravity Maps of Local Areas and Joint Assessment of Their Navigation Quality**

## A.V. Sholokhov , S.B. Berkovich , N.I. Kotov , and M.G. Belonozhko

#### **Abstract**

The actual problem of creation of gravity maps which can be used as correctors in mapmatching- navigation systems (MMNS) is considered. It is assumed to use information from two sources: available gravity maps and measurement data of various parameters of the Earth's gravity field (EGF) (acceleration of gravity, deflection of the vertical, horizontal gradients of the potential).

Obtaining the required amount of measurement data of a specified accuracy is considered as a problem of planning primary measurements. At the first stage of its solution, the number of measurements necessary to achieve the required accuracy of the parameters of the gravitational map is minimized. It takes into account the accuracy and performance of measuring devices. The second stage assumes that the human operator specifies the position of the points at which to measure the parameters of the EGF with a known accuracy.

Joint processing of information is carried out on the basis of the collocation method. The estimation of EGF navigation quality is considered in a specified point and its nearest neighborhood. This approach is reasonable, if in the future it is required to optimally form the trajectory of the MMNS in a specified area or to a specified destination point. It is assumed rectilinear motion of the object from this point with a known initial level of coordinate errors. Position errors are defined as functions of distance from the initial point and direction of movement.

#### **Keywords**

Bouguer anomaly - Gravity field - Gravity map -Gravity survey

#### **1 Introduction**

Currently, there is a growing interest in the use of EGF for autonomous navigation of mobile objects. This is due to the following circumstances:


One of the actual problems of autonomous navigation with the use of EGF is the task of gravity maps formation (Beloglazov et al. 1985; Berkovich et al. 2018; Dzhandzhgava et al. 2013; Kotov et al. 2017). The basic requirements for such maps are determined by the level of errors, ensuring their effective use in the MMNS as correctors, as well as the initial alignment. Existing EGF maps and models do not fully

© The Author(s) 2022

J. T. Freymueller, L. Sánchez (eds.), *5th Symposium on Terrestrial Gravimetry: Static and Mobile Measurements (TG-SMM 2019)*, International Association of Geodesy Symposia 153, https://doi.org/10.1007/1345\_2021\_138

A.V. Sholokhov (-) · S.B. Berkovich · N.I. Kotov · M.G. Belonozhko Interregional Public Institution "Institute of Engineering Physics", Serpukhov, Moscow Region, Russia e-mail: naviserp5@iifmail.ru

meet the accuracy requirements. Requirements are increased strongly due to the improvement of EGF parameter sensors.

Precise measurement of EGF parameters (an order of magnitude more accurate than maps) at separately located points remains very costly at the present time. That is why it is important to improve methods of information processing, allowing to form (quickly update existing) gravity maps with the attraction of additional precise measurements.

The article deals with three relatively independent problems of ensuring the creation of gravitational maps for the required accuracy of MMNS. The first one is to ensure the specified accuracy of maps by performing minimum measurements of the gravitational parameters at sites on the Earth's surface. The problem of optimal estimation of map errors using additional measurements is the basis of the proposed method for creating updated maps. Its peculiarity is the relatively small number of measurements performed mainly at sites on highways, as well as the different accuracy of measurements and the correlation of their errors. The problem of evaluating the navigation quality (informativeness) of the navigation field is solved for the case when its results are used in forming the trajectories of the MMNS movement in a specified area or to a predetermined endpoint.

#### **2 Optimization of Projects for Primary Measurements in the Creation of Maps**

This optimization is due to the need to increase the accuracy of the parameters of the EGF with limitations on the possibility of carrying out measurements at points called basic (BP) on the ground surface. These capabilities are established by the number of BP and the accuracy of the EGF parameters. In the BP the composition of primary measurements of EGF parameters is: altitude, acceleration of gravity, and components of the plumb deviation, horizontal potential gradients. As a rule, in an arbitrary BP on the earth's surface, it is assumed to measure the height and gravity. In addition, deflection of the vertical and horizontal gradients can be measured. The set of the parameters and the requirements to the accuracy of their measurement in a particular BP are the result of solving the optimization problem under consideration. BP's position is a priori unknown. Therefore, an integrated index of planning quality is introduced, namely, the mean value of the root-mean square error (RMSE) of the EGF parameter in the local area of the terrain. It is used as an optimality criterion to be minimized.

Among the limitations in the problem are the edges of the area, the number of the same type of BP (such as where the same set of the measured parameters of the EGF). The problem is to find the coordinates of all BP, for which the integrated quality index is minimal. The main difficulty of finding its optimal solution is seen in a large number of optimized parameters – BP coordinates. It is also necessary to take into account that such solution is only the "first approximation" of the optimal solution of the problem. The final solution made by the human operator must take into account other significant factors that cannot or should not be included in the problem to be solved. Given this, it is permissible to search for a suboptimal solution. The proposed solution to the problem of BP placement, which provides a demanded accuracy of EGF parameters, includes two stages:


The minimum required number of BP is found through the density of the location of the BP on the ground surface, which is determined by the average distance between the BP. In General RMSE *<sup>p</sup>* of geodesic parameter *p* or integral index can be represented by the function

$$
\sigma\_p = f\_p\left(r, \sigma\_m, \mathbf{Q}\right), \tag{1}
$$

where *<sup>r</sup>* – grid step, *<sup>m</sup>* – RMSE of measurements in BP, **Q** – other cartographic source data BP coordinates, parameters of the covariance function of the gravity anomaly error. The function *fp* is formed using the RMSE *<sup>p</sup>* of the EGF parameters obtained by the collocation method for the values of its arguments specified in the grid nodes. Between nodes values *<sup>p</sup>* are found by interpolation.The relationship between the desired RMSE *<sup>p</sup>*, the number of BP, which is funded by the grid step *r*, and accuracy of additional measurements in BP allows to solve the inverse problem of finding the function's argument, i.e. when *<sup>p</sup>* is given. The program allows quickly solving the problem of planning the finding of these data, based on the specified levels of error (or the integral index of accuracy). The result of solving such problems provides the final specification of the BP position (their placement on the map) by the operator on the local terrain area.

#### **3 Methods for Solving the Problem of Creating Maps**

The creation of gravity navigation maps includes the solution of the following particular problems:


**Fig. 1** The example of gravimetric map (Bouguer anomaly) with additional measurement points


Joint data processing of available EGF maps and precise additional measurements is based on the collocation method (Berkovich et al. 2018; Sholokhov and Druzhinin 2011; Moritz 1980). It allows finding the optimal estimates of the EGF parameters at an arbitrary point, taking into account additional measurements made beforehand at BP within the local area. Characteristics of the accuracy of the obtained estimates can also be found.

One of the problematic issues of finding parameter estimates is the choice of a priori values of parameters characterizing the errors of EGF maps. It significantly affects the final estimates of the EGF parameters, since the amount of measurement data is small. To solve this problem, we use a special approach that allows us to take into account the optimal covariance of estimates of the desired navigation parameter at a given point, obtained with different combinations of values of a priori data on the errors of maps.

The solution to this problem is realized by specialized software. A sample main program window depicting Bouguer anomaly layer with points in which additional measurements are made (marked with the symbol ), shown in Fig. 1.

Program in conjunction with geographic information systems (Professional GIS "Panorama" 2019) provides a solution to the following applied tasks: import and visualization of gravimetric information, manual input of information from the various sensors, the joint processing of the measurement results, the formation of databases with the visualization of the results, estimation of accuracy of geodetic and gravimetric parameters, formation navigation maps layers of different gravimetric parameters (Bouguer and free-air gravity anomaly, the components of the deviation of the plumb line, the second derivative of the potential), adjustment of the maps according new measurement data, the reduction of gravimetric parameters at points with different heights to determine a gravimetric parameter of an arbitrary point of the route at the specified values of coordinates, evaluation of the navigation quality of the navigation field.

#### **4 Estimation of the Navigation Quality of the Navigation Field**

The informativeness of the EGF parameter is determined at a specified point, from which the rectilinear motion of the object in a known direction is assumed. This approach differs from the traditional one (Beloglazov et al. 1985; Koneshov et al. 2016; Peshekhonov et al. 2017). It allows you to prepare the initial data for the subsequent optimal formation of the trajectory of the MMNS in a specified area or to a specified endpoint. An indicator of navigation quality is the RMSE of the object's coordinates. Finding simple analytical expressions that allow to obtain achievable levels of RMSE of coordinates depending on these factors is difficult and in practice impractical.

Therefore, a numerical assessment of these levels is carried out, based on modeling the process of correction the coordinates of the MMNS. In it random values of the data

**Fig. 2** Example of navigation quality estimation for the specified point

are generated in accordance with the required RMSE. You need to find the following functions *fX*, *fY* for RMSE -*X*, -*Y* of coordinates. They are allow us to solve inverse problems – to find arguments of functions when levels of RMSE of coordinates are given. Such problems are of the principal interest, since in practice the levels of RMSE of required parameters – coordinates are set, and the conditions for achieving the desired accuracy, determined by the arguments of functions, are subject to definition. The solution of inverse problems is obtained by finding the roots of nonlinear equations

$$\begin{cases} \sigma\_X - f\_X\left(\sigma\_\emptyset, r\_\emptyset, \sigma\_m, L\_c\right)\big|\_{\sigma\_\emptyset, r\_\emptyset} = 0\\ \sigma\_Y - f\_Y\left(\sigma\_\emptyset, r\_\emptyset, \sigma\_m, L\_c\right)\big|\_{\sigma\_\emptyset, r\_\emptyset} = 0 \end{cases},\tag{2}$$

where *<sup>g</sup>*, *rg* – RMSE and spatial correlation radius of the EGF errors are specified, the accuracy *<sup>m</sup>* of the EGF parameters measurements is known, and the length *Lc* of the MMNS correction section is to be determined. The result of navigation quality estimation is presented in the form of a map of the minimum possible (achievable) location error for one point of a given object trajectory. An example is shown in Fig. 2. The contour lines show the achievable levels of distance RMSE - D q -2 X <sup>C</sup> -2 Y in the rectilinear motion of the object from the starting point of the trajectory (marked "C") in the appropriate direction. The error levels of the coordinates at the reference point, as well as other parameters of the object motion and the parameters of the MMNS are set by the operator.

#### **5 Conclusion**

The problems of creating gravity maps to specify the position of movable objects are considered. It is planned to use information from two sources: available gravity maps and measurement data of various parameters of the EGF. Joint data processing is performed on the basis of the standard collocation method. The peculiarity of the application of this method to find estimates of the EGF parameters is revealed. The estimation of EGF navigation quality is considered in relation to a given point and its nearest neighborhood. This approach is reasonable, if in the future it is required to optimally form the trajectory of the movable object in a specified area or to a given endpoint.

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/ licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Methods of Bistatic GNSS-Radio Altimetry for Determining Height Profile of the Ocean and Their Experimental Verification**

### Vladislav Lopatin and Vyacheslav Fateev

#### **Abstract**

The current state of satellite altimetry development allows the use of altimetry data in determining the detailed characteristics of the Earth's gravity field on the surface of the ocean in the form of models of geoid heights, deflection of vertical and gravity anomalies. The article presents the bistatic altimetry method based on the reflected signals of global navigation satellite systems (GNSS). In this method, measurement redundancy is necessary to solving the problem of determining the height of the geoid with high spatial resolution. The experiments showed the possibility of using the code and phase of the received signal. The experiments showed the possibility of using method of bistatic GNSS-altimetry to determine the height of the geoid.

#### **Keywords**

Altimetry - Geoid - GNSS-R -Measurement methods

#### **1 Introduction**

The definition of the geoid is one of the main tasks of geodesy, gravimetry and oceanography. Geoid is a surface of equal gravitational potential on the Earth, containing a point taken as the beginning of the height count. The surface of the geoid coincides with the surface of the World Oceans in the absence of disturbing forces such as wind, ocean currents, tides and conditionally continues under the continents. The height of the geoid is the elevation of the geoid above the ellipsoid.

Currently, the method of active satellite altimetry is used to measure the height of the geoid on the oceans. This method is based on the use of radio altimeters placed on board special spacecraft for geodetic purposes.

Satellite altimetry measurements with improved accuracy and spatial resolution can come closer in accuracy and resolution to the results of marine gravity surveys.

V. Lopatin (-) · V. Fateev

VNIIFTRI, Mendeleevo, Moscow Region, Russia e-mail: lopatin@vniiftri.ru;; fateev@vniiftri.ru

#### **2 Ways to Improve the Accuracy of Altimetry**

To date, over 10 national and international satellite altimetry projects have been implemented. Foreign altimetry missions are shown in Table 1.

To improve the spatial resolution of altimetry measurements in the medium term, it is necessary to use:


Satellite altimetry based on active monolocators has high measurement accuracy. However, the following disadvantages of active satellite altimetry can be identified:


It is suggested that GNSS-R be used to eliminate the disadvantages of active altimetry. GNSS-R is a remote sensing

© The Author(s) 2022

J. T. Freymueller, L. Sánchez (eds.), *5th Symposium on Terrestrial Gravimetry: Static and Mobile Measurements (TG-SMM 2019)*, International Association of Geodesy Symposia 153, https://doi.org/10.1007/1345\_2022\_139

**Table 1** Existing and future foreign altimetry missions


method of the Earth that uses GNSS signals reflected from the water surface, allowing simultaneous observations at several points in a wide band. This method is promising for mesoscale altimetry. GNSS altimetry advantages:


#### **3 Operating Principle**

The GNSS altimetry bistatic system is based on the reception of direct signals from navigation satellites by an antenna directed at the zenith. Signals reflected from the ocean surface in the specular zone are received by an antenna directed to the "nadir", The operating principle is presented in Fig. 1 (Zavorotny and Voronovich 2000; Gleason et al. 2005). To solve the problem of determining the height to the reflecting surface, it is necessary to determine the difference in the propagation time of the direct and reflected signals, the coordinates and speeds of the receiver and the satellite. The accuracy of determining the height depends on the noise of the equipment, the accuracy of accounting for the ionospheric delay, tropospheric delay, processing algorithms, etc. (Camps et al. 2017; Li et al. 2016).

#### **4 Methods of Bistatic GNSS-Altimetry**

Today we can talk about four methods:

(1) *Code method*. Code method can be called traditional. Receivers of bistatic altimetry system determine pseudorange and pseudodoppler frequency bias for direct and reflected signals, as well as commercial receivers. The received direct and reflected signals are correlated with a replica of the open code signal generated in the receiver. The method ignores the encrypted signal components, thereby reducing the achievable accuracy of bistatic GNSS altimetry (Zavorotny et al. 2014).


$$SNR(t) = \text{Accos}\left(\frac{4\pi}{\lambda}h\cos\left(\theta(t)\right) + \varphi\_o\right),\qquad(1)$$

where (*t*)—elevation angle of the satellite at time *t*, signal wavelength, *h*—height to the reflecting surface, ® *<sup>o</sup>* initial phase, A—amplitude of the total signal. Unknown parameters are *h*, ® *<sup>o</sup>*, A.

(4) *Phase method.* Phase measurements have higher accuracy than code measurements, but such measurements are ambiguous due to the uncertainty of the initial integer number of phase shift periods (Semmling et al. 2014).

**Fig. 1** Principle of operation of bistatic GNSS-R altimeter

**Fig. 2** Place of experiment p.1—place of measurements with the code method; p.2—place of measurements with the SNR method

However, if it is necessary to determine the change (increment) in height to the reflecting surface, the phase measurements can be used as more accurate on condition that there are no cycle slips. The change in height from time *t <sup>j</sup>* to *t <sup>i</sup>* is proportional to the change in the pseudophase difference:

$$
\Delta h = \frac{\Delta \varphi \left( t\_i \right)}{2 \sin \left( \theta \left( t\_i \right) \right)} - \frac{\Delta \varphi \left( t\_j \right)}{2 \sin \left( \theta \left( t\_j \right) \right)}, \tag{2}
$$

where ®—difference between the pseudophases of the direct and reflected signals.It is necessary that the measurements satisfy the Rayleigh criterion reference required. The Rayleigh criterion is fulfilled at low elevation angles of the spacecraft or at a small surface roughness. The accuracy of phase measurements can be centimeters.

#### **5 Experimental Verification of the Methods**

Experimental measurements were performed to test bistatic GNSS altimetry methods from a bridge over a reservoir. In the first case, the GNSS altimetry code method was tested. The place of measurements is shown in Fig. 2. The receiving antennas were located at point 1. The results of measurements of the height to the reflecting surface are preceded in Fig. 3. From the measurement results, we can conclude that the optimal elevation angles of satellites for this method are 30ı–90ı, and the measurement error may be less than 1 m.

The SNR method was tested in measurements at p. 2, the real and measured height of which to the reflecting surface

**Fig. 3** Height measurement results using the code method, true height from reflective surface 19.2 m

**Fig. 4** Results of height measurements using the method of SNR (blue—L1, green—L2)

was 17.1 m, Fig. 4. Such a noticeable difference between the heights of p.1 and p.2 is associated with the difference in the height of the bridge and different levels of the reservoir's water in different seasons in which measurements were taken.

An anechoic shielded chamber was chosen for testing the phase method. The measurement scheme is shown in Fig. 5. Providing the necessary roughness of the reflecting surface and the necessary elevation angles of the satellite is not always a feasible task for field measurements.

Semi-natural modeling using a GNSS simulator in anechoic chamber conditions allows (Frolov et al. 2017):


• perform measurements both on the code and phase delay of the signal.

Antennas receiving direct and reflected signals were located on a mast with a variable height. GNSS satellites signals were simulated by a simulator manufactured by Navis. The results of measuring the difference between the pseudophases of the direct and reflected signals with a decrease and a further increase in mast height are shown in Fig. 6. The standard deviation of measurements was no more than 0.2 cm.

The phase altimetry measurements with the CYGNSS bistatic GNSS system over two independent passages above Qinghai Lake are shown in (Li et al. 2018). As the water surface of the lake approaches the equipotential surface of gravity, altimetry measurements above the lake can give an independent estimate of the height of the geoid along the track. The measurement data are relative due to an unknown

**Fig. 5** The scheme of the experiment in an anechoic chamber

**Fig. 6** The scheme of the experiment in an anechoic chamber

integer number of phases. The measurements obtained after introducing all the necessary corrections were compared with the EGM2008 model along the track. The measurement results well repeat EGM2008, but anomalies of ˙ 10 cm are distinguished, which are repeated in time for different GNSS satellites.

#### **6 Summary and Conclusions**

Thus, the experiments show the consistency known methods of the bistatic GNSS-altimetry. An anechoic chamber and a GNSS signal simulator can be used to simulate the operation of a GNSS altimetry system using any of four known methods. The highest requirements for a GNSS signal simulator are imposed during phase measurements. The accuracy of determining the height can reach less than 0.1 m with space placement of a bistatic GNSS altimeter.

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/ licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Absolute and Relative Gravity Measurements at Volcanoes: Current State and New Developments Under the NEWTON-g Project**

## F. Greco, D. Carbone, F. Cannavò, A. A. Messina, and G. Siligato

#### **Abstract**

Gravity changes associated with volcanic processes occur over a wide range of time scales, from minutes to years and with magnitudes between a few and a few hundred microGal. High-precision instruments are needed to detect such small signals and both timelapse surveys along networks of stations, and continuous measurements at single points, are accomplished. Continuous volcano gravimetry is mostly carried out through relative gravimeters, either superconducting instruments, providing higher quality data, or the more widely used spring meters. On the other hand, time-lapse surveys can be carried out with relative (spring) gravimeters, that measure gravity differences between pairs of stations, or by absolute gravimeters, capable of measuring the absolute value of the gravitational acceleration at the observation point. Here we present the state-of-the-art of terrestrial gravity measurements to monitor and study active volcanoes and the possibilities of new gravimeters that are under development. In particular, we present data from a mini array of three iGrav superconducting gravimeters (SGs) at Mount Etna (the first network of SGs ever installed on an active volcano). A comparison between continuous gravity measurements recorded through the iGrav#016 superconducting gravimeter at Serra La Nave station (1730 m a.s.l.) and absolute gravity data collected with the Microg LaCoste FG5#238 gravimeter in the framework of repeated campaigns is also presented. Furthermore, we introduce the Horizon 2020 NEWTON-g project (New Tools for Terrain Gravimetry), funded under the FET-OPEN Research and Innovation Actions call, Work Programme 2016–2017 (Grant Agreement No 801221). In the framework of this project, we aim to develop a field-compatible *gravity imager*, including an array of low-costs Micro-Electro-Mechanical Systems (MEMS)-based relative gravimeters, anchored on an absolute quantum gravimeter. After the design and production phases, the *gravity imager* will be field-tested at Mt. Etna (Italy) during the last 2 years of the project.

**Keywords**

Gravimeters - Laser-cooled atom -MEMS

F. Greco (-) · D. Carbone · F. Cannavò · G. Siligato Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Catania - Osservatorio Etneo, Catania, Italy e-mail: filippo.greco@ingv.it

#### **1 Introduction**

High-precision gravity measurements provide a powerful tool for volcano monitoring, since they highlight processes that induce bulk mass/density changes. This technique thus provides an important contribution to the understanding of the processes driving magma ascent during the preparatory phases of eruptive events.

© The Author(s) 2020

J. T. Freymueller, L. Sánchez (eds.), *5th Symposium on Terrestrial Gravimetry: Static and Mobile Measurements (TG-SMM 2019)*, International Association of Geodesy Symposia 153, https://doi.org/10.1007/1345\_2020\_126

A. A. Messina Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Roma 2, Roma, Italy

In spite of its potential, volcano gravimetry is not widely adopted by research institutions and agencies in charge of monitoring active volcanoes. This depends on several factors, including: (1) the cost of instrumentation, (2) problems that are inherent with the use and deployment of instruments intended for laboratory conditions in the harsh environments that characterize the summit zones of most active volcanoes, and (3) difficulty in interpreting gravity changes that may result from volcanic processes, as well as hydrological effects (changes in the undergroung water mass), and instrumental artefacts. These difficulties are not insurmountable and efforts in a variety of volcanic settings highlight the value of time-variable gravimetry for understanding volcanic hazards, as well as revealing fundamental insight into the volcano dynamics (Carbone et al. 2017).

The Istituito Nazionale di Geofisica e Vulcanologia (INGV) has operated a relative gravity network for the monitoring of Etna volcano since 1986 (Budetta et al. 1989). The network has been developed and has evolved over the years and currently it consists of: (a) 71 benchmarks, covering an area of about 400 km2, for relative gravity campaigns (LaCoste & Romberg model D and Scintrex CG-3M/CG-5 gravimeters were used over time); (b) three continuously running gravity stations equipped with iGravs superconductive gravimeters by GWR Instruments, Inc.; (c) 14 stations for absolute gravity (AG) measurements using the Microg LaCoste FG5#238. In all stations for AG measurements and in some stations for relative measurements, the vertical gravity gradient is also measured.

The whole Etna gravity network is routinely occupied every summer (in the same season of the year to minimize seasonal variation effects). More frequent measurements are carried out through relative spring gravimeters in same elements of the network (almost monthly along the East-West profile; Fig. 1; Carbone et al. 2009; Greco et al. 2010; Del Negro et al. 2013; Bonforte et al. 2017) or during volcano activities. Since 2007, several absolute field surveys have been routinely carried out at Etna volcano to check the long term gravity variations and to confirm the gravity changes obtained by relative measurements (Pistorio et al. 2011; Greco et al. 2012, 2015). To collect the absolute gravity data, the IMGC-02 absolute gravimeter (D'Agostino et al. 2008), developed by Istituto Nazionale di Ricerca Metrologica (INRiM, Italy), was used in the period 2007– 2009. Since 2009 up to date the Micro-g FG5#238 absolute gravimeter is used.

In this paper, we present data collected with SGs and absolute gravimeters at Mt. Etna. Furthermore, we introduce the field-compatible *gravity imager*, that will be developed under the NEWTON-g project, including an array of lowcosts MEMS-based relative gravimeters, anchored on an absolute quantum gravimeter. This measuring system will provide continuous imaging of gravity changes associated with variations in subsurface fluid properties, with unparalleled spatio-temporal resolution and, consequently, with fundamental implications for risk management. We will deploy the new gravity imager at Etna volcano (Italy), where frequent gravity fluctuations, easy access to the active structures and the presence of a multiparameter monitoring system ensure an excellent natural laboratory for testing the new tools.

#### **2 The Mini-Array of Three iGrav SGs at Mount Etna**

Continuous gravity measurements have been successfully carried out at a number of volcanoes around the world using spring gravimeters. Nevertheless, these instruments do not provide reliable continuous data over intervals of weeks or more, because they are influenced by environmental factors and are subject to instrumental drift. Accordingly, most studies of continuous gravity at active volcanoes have focused on the analysis of changes over time-scales of minutes to a few days (Branca et al. 2003; Carbone et al. 2006, 2013, 2019).

An alternative to spring gravimeters for continuous measurements is given by superconducting gravimeters (SGs) that feature a much higher precision and stability than spring gravimeters. SGs are free from instrumental effects and thus allow to track even small gravity changes (1–2 - Gal) over a wide range of time scales (minutes to months). However, even the most portable SGs (e.g., the iGrav® by GWR) are not ideal for installation in the vicinity of active volcanic structures. Indeed, they require AC power at the installation site and some kind of hut or vault to house the instrumentation (Carbone et al. 2019).

At Mt. Etna, the installation of a mini-array of three SGs (distances of 3.5–15.5 km from the active craters; Fig. 1) began in September 2014, when the iGrav#16 was installed in the facilities of the Serra La Nave (SLN) Astrophysical Observatory (1,730 m above sea level (asl); 6.5 km from the summit craters; Fig. 1). In the summer of 2016, iGrav#20 and iGrav#25 were also installed at La Montagnola hut (MNT; 2,600 m asl; 3.5 km SE of the summit craters; Fig. 1) and in the facilities of INGV – OE located in the village of Nicolosi (NIC; 720 m asl; 15 km SE of the summit craters; Fig. 1), respectively.

To our knowledge, these are the first SGs ever installed on an active volcano.

Signals from these instruments have shown hydrologically-induced components superimposed on gravity changes that are related to volcanic processes (Carbone et al. 2019).

Figure 2 shows the signals acquired in the period 23 July– 05 September 2016 at NIC (top; iGrav#25), SLN (middle; iGrav#16) and MNT (bottom; iGrav#20) stations, after the corrections for the effect of Earth tides (Wenzel 1996), local atmospheric pressure variations (Merriam 1992), and polar motion (Hinderer et al. 2015), are accomplished.

**Fig. 1** Sketch map of Mt. Etna showing the mini-array of three SGs has been in continuous operation at Mount Etna since the summer of 2016. The stations belonging the East-West profile (blue circles) where

The reduced signals show common anomalies over intervals of 10–15 days. In particular, a gravity decrease is observed in the signals from SLN and MNT stations during August 8–18. The MNT/SLN amplitude ratio during this period points to the activation of a mass source located a few kilometers below sea level.

Figure 2 also shows positive variations on time scales of the order of a few hours on the signal acquired at MNT (red arrows). Since these variations (average amplitude equal to about 2 microGal) are not present on the SLN and NIC signals and occur during heavy rainfall, they are most probably due to local hydrological effects.

### **3 Comparison Between the Time Series from a SG and Absolute Gravity Data**

Figure 3 shows a 56-month long time series from iGrav#16 SG (at SLN station), after removing the effect of Earth Tides, atmospheric pressure and polar motion. Figure 3 also shows absolute gravity observations (red points) at the same

more frequent measurements are carried out through relative spring gravimeters are also indicated

station (the total error on AG measurements at this station is typically ˙ 3 microGal).

Once a linear fit, obtained through AG data, is removed, there is a good fit between the two data sets. Continuous gravity changes are within 1–2 microGal of gravity changes measured using AG at the same site. This comparison allows to estimate the long-term drift of the iGrav SG, which is of the order of 9 microGal/year (Fig. 3).

This example shows how repeated AG measurements and continuous gravity observations can be used together to get a fuller and more accurate picture of the temporal characteristics of the studied processes.

#### **4 New Tools for Terrain Gravimetry**

Gravimetry is a powerful geophysical tool that, through sensing changes in subsurface mass, can supply unique information on the dynamics of underground fluids, like water, magma, hydrocarbons, etc. This is critically important for both resource management and risk reduction. In spite

**Fig. 2** Gravity time series from NIC (top), SLN (middle) and MNT stations during the 23 July to 03 September 2016 interval, corrected for the effect of Earth tides and local atmospheric pressure. The time series are also high - pass filtered (cutoff of 0.01 Hz)

of this potential, gravimetry is currently underexploited. Indeed, continuous gravity measurements have been rarely performed, owing to the fact that the available instruments are not well suited for continuous measurements under severe field conditions. In addition, the high cost of available gravimeters limits deployment to, at most, a few sensors in a given area, implying that the achievable spatial resolution is always lower than needed.

To overcomethe current limits of terrain gravimetry, the NEWTON-g project, that received funding from the European Commission's Horizon 2020 programme, under the FET-OPEN 2016–2017 call, proposes the development of a new generation of gravimeters, based on MEMS and quantum technologies (www.newton-g.eu). MEMS technology will make possible the production of a device that is two orders of magnitude cheaper and lighter than currently available gravimeters. In the framework of NEWTON-g, one goal is focused on strategies to increase the stability of ordinary MEMS accelerometers, thus making them capable of measuring, beside inertial, also gravitational acceleration, with a sensitivity of a few tens of microGal/hz1/2 (Middlemiss et al. 2016).

The quantum gravimeter that will be developed under NEWTON-g will be the ideal complement to the MEMS devices. We expect that the new quantum gravimeter would: (a) measure the absolute value of the gravity acceleration (absolute instrument), with a sensitivity and a stability at the -Gal level; (b) perform continuous measurements at a high rate (2 Hz) over a long time interval (several years).

A network of low-cost MEMS gravimeters, anchored to the absolute quantum device, will form the so-called *gravity imager*, which will provide imaging of sub-surface mass changes with unparalleled spatio-temporal resolution (Fig. 4).

Once developed, the gravity imager will be field-tested at Etna volcano (Italy) during the last 2 years of the NEWTONg project (summer 2020 – summer 2022). Insights from the gravity imager will also be used for volcanic hazards analysis, to demonstrate the importance of using gravity to face problems of societal relevance. A successful implementation of NEWTON-g will open new doors and represents a fundamental step that can move gravimetry into a cornerstone resource for geophysical monitoring and research.

#### **5 Conclusions**

Many geophysical phenomena are associated with mass transport and may induce gravity changes measurable at the surface over a wide range of time scales. Gravimetry can thus supply unique information on mass transport within the Earth system.

**Fig. 3** (Top) time series from iGrav#16 at SLN station during the September 2014–April 2019 interval. Data are corrected for the effect of Earth tides, local atmospheric pressure and polar motion (grey signal). The time series is also high – pass filtered (cutoff of 0.01 Hz). The green

Through the experiments accomplished at Etna during the past few decades many steps forward have been made towards the regular acquisition of high-quality gravity data. The main issues with continuous gravity measurements at active volcanoes is to obtain and maintain data to a suitable standard (as for quality and continuity) against an adverse environment (high altitude, inaccessibility for several months, lack of mains electricity for power, variable temperature, pressure and humidity, seismicity, corrosive gases) and using instruments which are intended for use under laboratory conditions.

The combined use of discrete (absolute and relative) and continuous gravity measurements is a unique tool both for studying the internal dynamic of a volcano and for surveil-

curve is the same signal, after correction for a linear trend deduced from the AG data (red points). (Bottom) same signals as in the top panel, but during the January 2017–April 2019 interval

lance purposes. Such combined system has already allowed important conclusions to be drawn at Mt. Etna.

The state-of-the-art instruments are, in most cases, heavy, expensive, power-hungry and difficult to operate. This prevents the installation of dense networks of continuously running gravimeters, implying that suitable spatio-temporal resolution cannot be attained.

Less expensive and easier to deploy instruments will open new horizons for terrain gravimetry. In particular, the gravity imager that will be developed and field-tested under the NEWTON-g project will demonstrate the possibilities of a new generation of MEMS and quantum gravimeters, allowing a step that can turn gravimetry into a cornerstone resource for volcano monitoring and hazard assessment.

**Fig. 4** The NEWTON-g *gravity imager* (currently under development) consist of an array of low-cost MEMS relative gravimeters, anchored to an absolute quantum gravimeter

#### **References**


absolute gravimeter: measurement apparatus and applications in geophysics and volcanology. Ann Geophys 51(1):39–49


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/ licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **The Effect of Measurement Preprocessing in the Gravity-Aided Navigation**

Oleg Andreevich Stepanov and Aleksei Sergeevich Nosov

#### **Abstract**

The paper analyzes the effect of preliminary processing of gravity measurements on the accuracy of the marine gravity-aided navigation. The preliminary processing of the measurements is implemented in the filtering and smoothing modes. Obtained results are illustrated by a one-dimensional example of gravity-aided navigation problem.

#### **Keywords**

Accuracy analysis - Filtering and smoothing - Gravity-aided navigation - Preliminary measurement processing

#### **1 Introduction**

Map-aided navigation systems using geophysical fields for a position update are employed for a wide class of vehicles (Stepanov 1998; Bergman 1999). The vehicle position is updating through the comparison of the measured and reference samples (profiles) of the field along the vehicle trajectory. In marine navigation it is common to use the field of Earth gravity anomalies (GA), since it is stable in time and can be measured using well-developed inertial sensors: high-precision accelerometers and gravimeters (Bishop 2002; Sokolov et al. 2019).

However, the measurement of the gravitational field onboard the moving vehicles has its own peculiarities: the signal measured by the gravimeter consists of both GA and inertial accelerations of the object, considered as errors. There are two possible ways for solve the gravityaiding problem: use raw measurements and consider the detailed error model in the map-aiding algorithm or use the preprocessed measured samples in which the GA along the vehicle trajectory is estimated some way. In the second case, often used in practice, the measurements of the gravimeter are pre-processed (Wu et al. 2017). After that, the amplitude of the measurement errors of the GA significantly decreases, which allows us to simplify the map-aiding algorithm. At the same time, as it is shown in (Nosov and Stepanov 2018), the preliminary processing of measurements and the subsequent simplification of the algorithm can lead to an increase in the navigation error.

This paper analyzes the effect of preliminary processing of gravity measurements on the accuracy of the marine gravity-aided navigation. It continues the authors' study on a similar problem for underwater terrain-aided navigation in a presence of a white noise measurement error (Nosov and Stepanov 2018).

#### **2 Optimal and Suboptimal Solutions of the Gravity-Aided Navigation Problem**

Following (Nosov and Stepanov 2018), we consider a gravity-aided navigation problem in the simplest formulation, namely, we need to estimate a constant scalar random variable - (for example, a navigation system (NS) error from one of the coordinates) using a previously constructed GA map. This problem can be stated in Bayesian framework as follows (Stepanov 1998; Bergman 1999): to estimate the unknown parameter Eq. (1) using *p* scalar

O.A. Stepanov · A.S. Nosov (-)

CSRI Elektropribor, JSC, ITMO University, Saint Petersburg, Russia

measurements Eq. (2)

$$
\Delta\_i = \Delta\_{i-1} = \Delta,\tag{1}
$$

$$y\_i = \phi \left(\mathbf{x}\_i - \Delta\right) + \varepsilon\_i = \phi\_i \left(\Delta\right) + \varepsilon\_i,\tag{2}$$

where *xi* are coordinates, provided by the reference navigation system, e.g. inertial one, in a discrete set of points; "*<sup>i</sup>* are random measurement errors; i D 1:::p; (•) is a exactly known function (map) describing the dependence of the field on coordinate. For simplicity, we consider and "*<sup>i</sup>* as Gaussian with known stochastic properties.

To solve this problem in optimal way, the well-known numerical point-mass and sequential Monte-Carlo methods are used (Stepanov 1998; Bergman 1999; Nordlund 2002; Gustafsson et al. 2002). However complex behaviour of measurement errors "*<sup>i</sup>* which in turn requires high-dimension model makes the estimation algorithm implementation computationally expensive, because nonlinear Bayesian estimation methods are subject to "curse of dimensionality" (Daum and Huang 2003). In this case it is suitable to pre-process measurements to estimate the GA along the vehicle trajectory, i.e. the value of *i*(-). The measurements Eq. (2) are represented as a sum of some random process which describes the field values *gi* (*xi* -) and the error same as above *yi* D *gi* C "*i*. The problem of estimating *gi* is linear one, which allows the use of Kalman-type filtering and smoothing algorithms (Kalman 1960; Gelb et al. 1974; Peshekhonov and Stepanov 2017).

After the preliminary processing, the measurements used to solve the navigation problem can be written as follows:

$$
\hat{y}\_k = \phi \left( \mathbf{x}\_k - \Delta \right) + \xi\_k,\tag{3}
$$

where yO<sup>k</sup> is the estimate, provided by the preliminary filtering or smoothing algorithm; −*<sup>k</sup>* is the corresponding estimation error; k D 1:::n.

Note that preliminary processing in itself does not make the algorithm for estimating simpler. Moreover, if we use the whole set of measurements *Yp* after the pre-processing and take a proper account of the estimation errors, we can show that the estimation accuracy of will remain at the same level. However, since preliminary processing significantly reduces the measurement errors it becomes possible to reducethe number of measurements used to estimate - from *p* to *n*. This is achieved by increasing the interval *t*1 for measurements Eq. (2) compared with the interval *t*<sup>2</sup> for measurements Eq. (3). Furthermore if the interval between measurements Eq. (3) is chosen to exceed the correlation interval for the error −*k*, then its model can be approximated by discrete white noise, the level of which will correspond to the solution of a filtering problem or a smoothing one. This simplifies the model used in nonlinear algorithm, since there is no need to describe the complex behavior of measurement errors.

As indicated in the introduction, such a two-stage scheme for gravity-aided navigation algorithm can lead to an increase in positioning errors. To evaluate losses in accuracy, we will use the procedure based on comparing the unconditional calculated and real root-mean square errors (RMSE) for optimal and suboptimal (two-stage) algorithms (Nosov and Stepanov 2018).

#### **3 Accuracy Analysis Example**

Let us consider example of the gravity-aided navigation problem and compare the unconditional RMS errors for optimal and suboptimal (two-stage) algorithms. We have to specify stochastic models for random process which describes the GA values and measurement errors. For the GA profile g.t / Q along a rectilinear trajectory we use the Jordan model in the form of a stationary process (Jordan 1972). Its correlation function can be written as follows:

$$K\_{\tilde{\mathcal{g}}}\left(\rho\right) = \sigma\_{\tilde{\mathcal{g}}}^2 \left(1 + \alpha \rho - \frac{\left(\alpha \rho\right)^2}{2}\right) e^{-\alpha \rho},\tag{4}$$

where ¡ 0. The corresponding shaping filter is written as (Koshaev and Stepanov 2010; Peshekhonov and Stepanov 2017):

$$\begin{cases} \dot{\mathbf{x}}\_1 = -\beta \mathbf{x}\_1 + \mathbf{x}\_2, \\ \dot{\mathbf{x}}\_2 = -\beta \mathbf{x}\_2 + \mathbf{x}\_3, \\ \dot{\mathbf{x}}\_3 = -\beta \mathbf{x}\_3 + q\_{\bar{\mathbf{g}}} w\_{\bar{\mathbf{g}}}, \\ \tilde{\mathbf{g}} = -\beta \zeta \mathbf{x}\_1 + \mathbf{x}\_2. \end{cases} \tag{5}$$

In Eqs. (4) and (5) is a distance along the trajectory; ˛– inverse value of the correlation window; " D V ¢@g=@¡ <sup>Q</sup> = <sup>p</sup>2¢g<sup>Q</sup> ; *V* is the vessel speed; @g=@ <sup>Q</sup> is the parameter defining GA spatial variability; qgQ*w*<sup>g</sup><sup>Q</sup> is forcing white noise with power-spectrum density (PSD) q<sup>2</sup> <sup>g</sup> <sup>D</sup> 10ˇ32 <sup>g</sup><sup>Q</sup> ; ¢<sup>2</sup> <sup>g</sup><sup>Q</sup> is the GA variance; and D <sup>p</sup><sup>5</sup> <sup>1</sup> = <sup>p</sup><sup>5</sup> is the dimensionless coefficient.

To describe the errors of GA measurements on a sea vessel, we consider vertical movement due to heaving and whitenoise error. The model describing vertical accelerations *av* can be represented in the form:

$$\begin{cases} \dot{\mathbf{x}}\_4 = \mathbf{x}\_5, \\ \dot{\mathbf{x}}\_5 = \mathbf{x}\_6, \\ \dot{\mathbf{x}}\_6 = -b\_3 \mathbf{x}\_4 - b\_2 \mathbf{x}\_5 - b\_1 \mathbf{x}\_6 + \mathbf{w}\_v, \\ a\_v = \mathbf{x}\_6, \end{cases} \tag{6}$$

**Fig. 1** Example of the GA measurements and results of preliminary processing

where *b*<sup>3</sup> D (<sup>2</sup> C 2) ; *b*<sup>2</sup> D <sup>2</sup> C <sup>2</sup> C 2 ; *b*<sup>1</sup> D 2 C ; C D v p2b3 .b1b2 b3/ =b1; *wv* is white noise of unit PSD; *<sup>v</sup>* is the RMS value of vertical movement *x*4; is the prevailing heaving frequency; is the coefficient of heaving irregularity; and is the dimensionless coefficient.

In this case, the gravimeter measurements can be written as

$$\mathbf{y} = \tilde{\mathbf{g}} + a\_v + v\_{gr} = -\beta \zeta \mathbf{x}\_1 + \mathbf{x}\_2 + \mathbf{x}\_6 + v\_{gr}, \quad (7)$$

where *vgr* is the white-noise measurement error with a known PSD.

Typical parameters for marine gravimetric survey were chosen for modeling the problem according to the discrete representation of Eqs. (5)–(7): the RMS value of the GA derivative was set to be 3 mGal/km, the period of vertical displacements was 7 s, and their RMS value was 0.2 m. Gravimeter measurements in the form Eq. (7) were simulated on a fixed section of several 30 km length GA profiles with coordinates (5,000 20,000) m. With spatial interval of 1 m it gives *N* D 15,000 raw measurements presented in Fig. 1a. Note that although the RMS value of the sea vessel vertical displacements is only 0.2 m, the RMS error of field measurements exceeds 14,700 mGal. It is obvious that against the background of such errors, it is impossible to estimate of the GA without processing. Figure 1b shows the results of preliminary processing in filtering and smoothing modes, as well as 3 bounds corresponding to the estimates. Values were obtained using the Kalman filter and the Rauch– Tung–Striebel smoother based on models Eqs. (5)–(7).

Using the GA estimates presented in Fig. 1b we can simplify the map-aiding by replacing model Eq. (6) in nonlinear part of the algorithm with white noise error model. It is feasible by subsampling the estimates, corrupted by correlated pre-processing errors.

To select a subsampling interval, we use the correlation functions of pre-processing errors. The choice of subsampling interval based on value of 0.3 for the normalized correlation function. Under this condition, the spatial intervals for the pre-processed measurements were approximately 1,200 and 800 m for the filtering and smoothing modes, respectively. The variance of the residual discrete white noise error was selected corresponding to the variance of the GA estimation error calculated during the pre-processing.

The raw and pre-processed measurements, represented by Eqs. (2) and (3) were used in optimal and suboptimal algorithms, respectively, to simulate the solution of the gravity-aided navigation problem. In the optimal algorithm, a four-dimensional state vector including the and vector Eq. (6) was estimated. The estimate was calculated as the mathematical expectation of the conditional probability density function (p.d.f). In the same time suboptimal algorithms estimated only the using pre-processed measurements. For calculations in both cases, we used the point-mass method with the number of nodes *L* D 3000, which was supplemented with the Rao-Blackwellization procedure for the optimal algorithm (Stepanov and Toropov 2015). The a priori RMS position error was set at 700 m.

Figure 2 shows the examples of p.d.f. graphs for each algorithm. Orange lines indicate the true error of the NS,

**Fig. 2** Examples of p.d.f. for optimal and suboptimal algorithms

**Table 1** Real RMSEs and calculation time for optimal and two-stage algorithms


blue ones indicate its estimates obtained in the gravity-aiding algorithm, green ones indicate graphs of a posteriori p.d.f. depending on the distance covered.

Graphs below show calculated and real unconditional RMSEs for the algorithms under study. They were calculated using 250 Monte Carlo simulations. Table 1 contains results including the values of real RMSE obtained at the end of the observation for *n* measurements and mean calculation time on the reference computer.

From the Fig. 3 and table above, it is obvious that the twostage suboptimal algorithm with preliminary smoothing is comparable to the optimal algorithm in accuracy: when it was applied, the unconditional real RMSE of these algorithms were in 50–70 m range. In addition, calculated RMSE provided by the suboptimal algorithm with smoothing is close to the real values.

Algorithm with preliminary filtering of measurements underperforms both in accuracy and identity of real and calculated accuracy characteristics. Besides the smaller number of measurements used and greater variance of the estimation error, this can be explained by the presence of a phase delay in the field estimates produced by the filter. Although this algorithm is the fastest, low accuracy precludes its use.

As for the amount of calculations, based on time averaging of 250 runs of algorithms on the test computer, we determined that it took 30 s to process all measurements in the optimal algorithm and an order of magnitude less, that is, 3 s, for the two-stage algorithm which includes preliminary smoothing.

**Fig. 3** Real and calculated unconditional RMSEs of the optimal and suboptimal algorithms

#### **4 Conclusion**

The effect of preliminary gravity measurement processing on the accuracy of gravity-aided navigation has been analysed. It is based on comparison unconditional real RMS estimation errors for the two-stage suboptimal algorithms that use the measurement pre-processing with potential accuracy achieved by optimal Bayesian algorithm. The preliminary processing of the measurements is implemented in the filtering and smoothing modes.

The example of gravity-aided navigation problem has been considered. It has been shown that the accuracy of the suboptimal algorithm with preliminary smoothing is close to the potential one, and the amount of calculations has been reduced by an order of magnitude in comparison with optimal algorithm.

In the further research, it is supposed to consider map errors model and generalize the results to the case of update two-dimensional position of the vehicle.

**Acknowledgements** This work was supported by the Russian Science Foundation, project no. 18-19-00627.

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/ licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Sensitivity of Algorithms for Estimating the Gravity Disturbance Vector to Its Model Uncertainty**

O. A. Stepanov, D. A. Koshaev, O. M. Yashnikova, A. V Motorin, and L. P. Staroseltsev

#### **Abstract**

The work considers the results of filtering and smoothing of the gravity disturbance vector horizontal components and focuses on the sensitivity of these results to the model parameters in the case when the inertial-geodesic method is applied in the framework of a marine survey on a sea vessel.

#### **Keywords**

Filtering - Gravity disturbance - Sensitivity -Smoothing

#### **1 Introduction**

The complexity of the Earth's surface and its interior is the reason why the direction of the gravity vector !g (vertical line) does not coincide with the direction of the normal gravity vector !n at points on the Earth's physical surface. This difference is characterized by the gravity disturbance vector (GDV) (Peshekhonov et al. 2017; Torge 2001). The problem of determining all three components of this vector is a problem of vector gravimetry. Note that the vertical component of the GDV coincides with the free-air gravity anomaly up to the accuracy of correction for the height anomaly (Jekeli 1999, 2016). The horizontal components of the GDV are deviations of the vertical, which are determined by two angles of deviation: in the planes of the Meridian Ÿ and the first vertical ˜ (Peshekhonov et al. 2017; Torge 2001).

Among the methods used to determine the GDV horizontal components on a moving base, including marine vessels, the gravimetric method, based on the measurement of gravity anomalies, and the astronomical-geodetic method, based on the comparison of astronomical and geodetic coordinates, and its modifications are particularly widespread (Koneshov et al. 2016; Peshekhonov et al. 1995; Hirt et al. 2010). Also well-known are the method of gravitational gradiometry, based on the measurement of the second derivatives of the geopotential, and the method of satellite or aircraft altimetry, which makes use of measurements of trajectory altitude (Koneshov et al. 2016). We cannot but mention methods based on the construction of global models of the Earth's gravity field using the data on perturbations of satellite orbits (satellite-to-satellite tracking and other satellite missions). Combinations of these methods (for example, the astronomical-gravimetric method) are used (Peshekhonov et al. 2017; Koneshov et al. 2016). The inertial-geodetic method, which is one of the modifications of the astronomical-geodetic method, is used in the framework of the vector gravimetry problem solution. It is based on a comparison of astronomical latitude and longitude ®, œ, generated by the inertial navigation system (INS), and geodetic latitude and longitude *B, L*, obtained from the GNSS data (Koneshov et al. 2016; Emel'yantsev et al. 2015). Thus, the horizontal components of the GDV are defined as:

$$\begin{aligned} \varphi - B &= -\xi \\ (\lambda - L) &= \eta / \cos \ B. \end{aligned} \tag{1}$$

In contrast to the conventional astronomical-geodetic method, in the inertial-geodetic INS method, the INS generates both astronomical coordinates and their derivatives, which makes it possible to use complementary velocity measurements (Koneshov et al. 2016; Dmitriev 1991).

O. A. Stepanov · D. A. Koshaev · O. M. Yashnikova (-) ·

A. V. Motorin · L.P. Staroseltsev

CSRI Elektropribor, JSC, ITMO University, St. Petersburg, Russia

J. T. Freymueller, L. Sánchez (eds.), *5th Symposium on Terrestrial Gravimetry: Static and Mobile Measurements*

*<sup>(</sup>TG-SMM 2019)*, International Association of Geodesy Symposia 153, https://doi.org/10.1007/1345\_2021\_132

The idea of using precision INS to determine horizontal components is based on the use of the relation between the errors of the INS output navigation parameters and the GDV components (Dmitriev 1997; Schultz and Winokur 1969; Li and Jekeli 2008). The analysis of this error, which is actually a methodical one, opens up the possibility, in principle, of determining the GDV horizontal components directly relative to the reference value at a reference point by solving the estimation problem using differences in the measurements of INS and GNSS.

The effectiveness of solving such a problem in real time (in filtering mode) was analyzed earlier in (Anuchin 1992; Nesenyuk et al. 1980; Staroseltsev and Yashnikova 2016). At the same time, experience shows that application of the smoothing mode, which does not require obtaining estimates of the desired parameters in real time, significantly increases the accuracy of the problem solution. Note that the feature of the smoothing problem is that when obtaining estimates at a current time moment, we can use both past and future measurements (Stepanov and Koshaev 2010; Loparev and Yashnikova 2012). The specifics of the inertial-geodetic method make it possible to implement optimal Kalman filtering and smoothing algorithms. However, such algorithms are based on the error models of the GNSS, INS and its sensors, as well as GDV horizontal components themselves in the form of random processes. To derive an adequate model of GDV horizontal components, it is necessary to have a priori information about the nature of the field in the survey area, for example, such parameters as the GDV variance and the correlation radius. Incorrect assignment of parameters can significantly reduce the accuracy of the resulting estimate. The analysis of the algorithm sensitivity is aimed to study the decrease in the estimation accuracy in the case that stochastic models in an algorithm are defined incorrectly (Stepanov and Koshaev 2003, 2011). Sensitivity of algorithms for determining the vertical component of the GDV was discussed in (Stepanov and Motorin 2019; Stepanov et al. 2015a, b).

The paper considers the results of filtering and smoothing of the GDV horizontal components and focuses on the sensitivity of these results to the model parameters in the case when the inertial-geodetic method is applied in the framework of a marine survey on a sea vessel. The studies involve direct calculation of the covariance matrices of estimation errors (without using the Monte Carlo method).

#### **2 Description of the Error Model**

To solve the problems of filtering and smoothing in this study, we used the model of INS errors described in (Gusinsky et al. 1996, 1997). The state vector of this model includes errors in constructing the local vertical, errors in

generation of the East and North velocity components, errors of the sensing elements, and the model of GDV horizontal components. Loosely coupled INS-GNSS system is assumed with coordinate, velocity, orientation angle, and sensor errors feedback for the filtering mode during survey. The fixed interval smoothing for every tack was fully done after survey.

The approach to creation of an integrated system to measure GDV horizontal components implies very high requirements imposed on the instrumental errors of the INS. Gyroscope drift and measurement errors of linear acceleration must not exceed 0.0001 deg./h and 1 arcsec, respectively (Emel'yantsev et al. 2015; Nesenyuk et al. 1980). Therefore, the coefficients of the gyroscope error model were described by first-order Markov processes with the standard deviation (SD) at a level of 3 - 10-<sup>5</sup>ı/h and a correlation interval of 40 h; nonwhite-noise errors of the accelerometers were set by the stationary first-order Markov processes with an SD at a level of 1 arcsec and a correlation interval of 10 h; the additive white noise of the accelerometers was also set at a level of 1 arcsec/pHz. These processes define stability of inertial sensors. The initial root mean square errors (RMSE) of the sensors were set equal to the steady-state values of the corresponding Markov process SD. The RMSE of the INS initial alignment were set at a level of 1 arcsec for tilt angles and at a level of 0.1 m/s for the horizontal components of the velocity.

As noted above, for the solution of the estimation problem by the inertial-geodetic method in the formulation considered here, models of the GDV horizontal components must be given in the form of random processes describing the variability of the components along the vehicle trajectory (Staroseltsev and Yashnikova 2016; Nash and Jordan 1978). The expressions of spectral densities for the longitudinal *S*"*L*(!) and transverse *S*"*<sup>T</sup>* (!) horizontal components of the GDV in the rectilinear motion of the vehicle with a constant velocity *V* were assumed to be similar to those given in (Nesenyuk et al. 1980):

$$S\_{\varepsilon L} \left( \omega \right) = \frac{4 \alpha\_{\varepsilon} \sigma\_{\varepsilon}^{2} \omega^{2}}{\left( \omega^{2} + \alpha\_{\varepsilon}^{2} \right)^{2}}; \quad S\_{\varepsilon T} \left( \omega \right) = \frac{2 \alpha\_{\varepsilon} \sigma\_{\varepsilon}^{2}}{\omega^{2} + \alpha\_{\varepsilon}^{2}}; \quad \alpha\_{\varepsilon} = \frac{V}{2.3d}, \tag{2}$$

where ", ˛", *d* are SD, damping coefficient, and inverse correlation radius, respectively.

For the solution of the estimation problem, the model included measurements from the INS and the GNSS receiver, which are parts of the integrated system. It was assumed that the position measurements of the GNSS had white noise at a level of 3 cm/pHz and a component represented by a stationary first-order Markov process with a SD of 3 cm and a correlation interval of 1 h. Velocity measurements of the GNSS had only a white-noise error at a level of 3 cm/s/pHz.

### **3 Estimating the Accuracy of Determining the GDV Horizontal Components: The Results**

To simulate the solution of the smoothing problem, we modified software in the Matlab package (Stepanov and Koshaev 2003, 2011), which allowed us to calculate the RMSE of filtering and smoothing estimates and display their plots. Error-tolerant algorithms were used to calculate error covariance matrices based on UD decomposition. Owing to this, all the necessary calculations were performed without additional computational errors.

In the simulation, it was assumed that the vessel was moving North along the meridian at a speed of 10 knots; the initial inertial longitude was 30<sup>ı</sup> E; the initial latitude was 60<sup>ı</sup> N. The solution interval was 24 h, and the discreteness – 10 s.

We set different characteristics of the gravitational field variability: from weakly disturbed with a of 5 arcsec and a correlation radius of 20 nmi to the highly disturbed with an RMSE of 20 arcsec and a correlation radius of 5 nmi. The steady-state values of the RMSE in the estimation of the GDV horizontal components obtained in the simulation for filtering and smoothing algorithms are shown in Table 1.

Using the smoothing mode to calculate the GDV horizontal components provides a significant gain in accuracy compared to the filtering mode. As may be seen from the table, for a weakly disturbed field (SD of 5 arcsec) only filtering algorithms can be used because there is practically no gain from smoothing. This fact is easy to explain: it is well known that the estimation accuracy of constant values (random bias) in the filtering and smoothing modes are the same. As for a strongly disturbed field (SD of 20 arcsec), the use of smoothing algorithms can give almost a two-fold increase in the accuracy of the estimation problem solution.

#### **4 Estimating the Sensitivity of Algorithms for Determining the GDV Horizontal Components to the Model Parameters: The Results**

For sensitivity analysis, we used the same parameters of the GDV model as in Table 1. When calculating the sensitivity for each combination of real (true) and estimated values of a model parameter, we determined three different RMSE in the estimation of the GDV horizontal components.

The first one is the RMSE of the error in estimating the optimal algorithm, which corresponds to the estimate in the case that the model specified in the algorithm corresponds to the real one. This RMSE characterizes the maximum achievable estimation accuracy of the problem under given conditions.

The second one is the real RMSE that correspond to the real estimation accuracy of the algorithm tuned to an incorrect model of a disturbed field. Real RMSE calculated using method from (Stepanov and Koshaev 2003, 2011). The method is based on the solution of the covariance equation for the augmented vector, which includes both optimal and suboptimal estimates of state vectors corresponding to correct and incorrect models. This RMSE will always be higher than the RMSE for the optimal algorithm.

The third RMSE are those calculated using the data from the covariance equation of the algorithms that stand for the characteristics of the accuracy of the obtained solutions. If for different values of the parameters of the real models in a certain range the calculated RMSE is not lower than the real one, it may be stated that the algorithm has a guaranteeing property. We assume suitable to provide such RMSE along with others because it is one of the ways to judge about accuracy in real survey.

**Table 1** The steady-state values of the RMSE in the estimation of the GDV horizontal components obtained in the simulation for filtering and smoothing algorithms


**Fig. 1** RMSE of the GDV horizontal components for the case of motion in a strongly disturbed field (field SD – 20 arcsec, correlation radius – 5 nmi) with the filter and smoother tuned to a weakly disturbed field (field SD – 5 arcsec, correlation radius – 20 nmi)

An example of simulation for optimistic tuning of the algorithm, that is, such that the algorithm is tuned to a less variable field is shown in Fig. 1. It is a case of motion in a strongly disturbed field (SD of the field is 20 arcsec, correlation radius – 5 nmi) with a setting to a weakly disturbed field (SD of the field – 5 arcsec, correlation radius – 20 nmi).

Analyzing the curves in Fig. 1, we can draw the following conclusions: for the motion in a highly disturbed field with the filter, which use incorrect weakly disturbed model. The estimation accuracy will be approximately 1.5 times worse than the optimal one and 2–3 times worse than the calculated one, which will be shown by the covariance equation of the suboptimal filter. These conclusions are also valid for the case of smoothing algorithms. However, due to the fact that the smoothing mode is more accurate than the filtering mode, the difference in the absolute values of the RMSE in estimating the GDV orizontal components between the suboptimal and optimal algorithms is not so significant here. As can be seen from Fig. 1, for the preset parameters, the guaranteeing property of the filter does not hold either.

The situation that arises is not satisfactory for estimation of the GDV horizontal components, which is why an attempt was made to find such a filter setting that would ensure min-

**Fig. 2** RMSE of GDV horizontal components for the case of motion in a strongly disturbed field (field SD – 20 arcsec correlation radius – 5 nmi) with the filter and smoother tuned to a field with SD of 15 arcsec and correlation radius of 10 nmi

imal loss to the optimal algorithm both in weakly disturbed and strongly disturbed fields. Figure 2 presents similar plots for the motion in a strongly disturbed field with a SD of 20 arcsec and a correlation radius of 5 nmi with the filter setting for the field SD of 15 arcsec and a correlation radius of 10 nmi.

With this setting, the suboptimal filter insignificantly loses to the optimal one; however, the guaranteeing property of estimation does not hold either. The results of the simulation of motion in a weakly disturbed field (the field SD – 5 arcsec, the correlation radius – 20 nmi) with the same setting are presented in Fig. 3.

It is obvious that in this case, too, the real RMSE of the error in estimating the GDV horizontal components is close to the optimal one. However, the calculated RMSE significantly differs from both the optimal and real RMSE.

**Fig. 3** RMSE of GDV horizontal components for the case of motion in a weakly disturbed field (field SD – 5 arcsec, correlation radius – 20 nmi) with the filter and smoother tuned to a field with SD of 15 arcsec and the correlation radius of 10 nmi

#### **5 Conclusion**

The accuracy of the problem solution in determining GDV horizontal components by the inertial-geodetic method on a marine moving vessel in the filtering and smoothing modes has been studied. It is shown that the use of the smoothing mode for estimation of GDV horizontal components leads to a significant gain in accuracy in a strongly disturbed field as compared to the filtering mode. However, for a weakly disturbed field, there is practically no gain from smoothing.

Sensitivity of the geodetic method algorithms for estimating GDV horizontal components to the accuracy of setting the model of these components was analyzed. The analysis has shown that with optimistic filter tuning (tuning to a field that is less variable than in reality), suboptimal algorithms can significantly lose in accuracy to the optimal algorithm, and they do not provide guaranteeing estimation. At the same time, within the known limits of changes in the correlation radiuses and the RMSE of the field in the survey area, it is possible to choose the model parameters so that the suboptimal algorithm will not differ noticeably in accuracy from the optimal one. However, it is difficult to obtain an adequate calculated characteristic of accuracy. The solution to this problem might be in the development of adaptive algorithms.

**Acknowledgement** This work was supported by the Russian Science Foundation, project no. 18-19-00627.

#### **References**


sionnogo polya Zemli (up-to-date methods and tools for measuring the earth gravity field parameters). Concern CSRI Elektropribor, ITMO University, St. Petersburg. (in Russian)


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons. org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **List of Reviewers**

Akito Araya Andrey Golovan Andrei Motorin Alexander Sokolov Andrea Prato Aleksei Sholokhov Christoph Foerste Christian Voigt Derek Van Westrum Dmitry Bobrov Dmitry Koshaev Filippo Greco Hartmut Wziontek Irina G. Ganagina Natalya Kuz'mina Olga Yashnikova Roland Pail Roman Sermyagin Sergey Petrov Silvia Alicia Miranda Vladimir Vasiliev Vadim Vyazmin Valerii Granovskii Yury Bolotin

# **Author Index**

#### **A**

Appleby, G., 97–103 Araya, A., v, 81–87

#### **B**

Becker, D., 15–22, 51 Becker, M., 15–22 Belonozhko, M.G., 121–125 Berkovich, S.B., 121–125 Bisi, M., 89–95 Bolotin, Y.V., 23–30, 52

#### **C**

Cai, T., 40, 45–49 Cannavò, F., 133–138 Carbone, D., 133–138 Chistiakova, E., 67–71

#### **D**

Desogus, S., 89–95

#### **F**

Fateev, V., v, 59–65, 127–131

#### **G**

Ganagina, I.G., 107–113 Germak, A., 89–95 Gidaspov, D.D., 75–79 Goldobin, D.N., 107–113 Golovan, A.A., 23, 24, 27, 39–42, 54 Gorushkina, E.V., 39–42 Greco, F., v, 133–138

#### **I**

Ianishevskii, V.V., 115–119

#### **J**

Johann, F., 15–21

#### **K**

Kanushin, V.F., 107–113 Karpeshin, F.F., 75–79

Karpik, A.P., 107–113 Kasai, K., 81–87 Khasiev, I.S., 75–79 Kolpashchikova, T.N., 31–38 Koshaev, D.A., 4, 142, 147–152 Kotov, N.Y., 121–125 Krasnov, A.A., v, 3–12, 17 Kulinich, R.G., 31–38 Kuz'mina, N.V., 3–12

#### **L**

Lopatin, V., 127–131

#### **M**

Mazurova, E.M., 107–113 Messina, A.A., 133–138 Mitrikas, V.V., 115–119 Motorin, A.B., 147–152 Murzabekov, M., 59–65

## **N**

Nakazawa, M., 81–87 Nosov, A.S., 141–145

#### **O**

Origlia, C., 89–95 Orlov, O.A., 75–79

#### **P**

Papusha, I.A., 39–42 Prato, A., 89–95 Proshkina, Z.N., 31–38 Pruglo, A., 59–65

#### **R** Ravdin, S., 59–65 Rodriguez, J., 97–103

#### **S**

Sheremet, V.I., 75–79 Sholokhov, A.V., 121–125 Siligato, G., 133–138

© The Author(s) 2023 J. T. Freymueller, L. Sánchez (eds.), *5th Symposium on Terrestrial Gravimetry: Static and Mobile Measurements (TG-SMM 2019)*, International Association of Geodesy Symposia 153, https://doi.org/10.1007/978-3-031-25902-9

#### **T**

Tsubokawa, T., 81–87

**V** Valitov, M.G., 31–38 Vitushkin, L.F., v, 4, 75–79 Vyazmin, V.S., 23–30, 51–56 **Y** Yashnikova, O.M., 147–152 Yoshida, M., 81–87 Yu, M., 40, 46

#### **Z**

Ziebart, M., 97–103

# **Subject Index**

#### **A**

Absolute ballistic gravimeters (ABG), 5, 75–79 Absolute gravimeter, 4, 5, 7, 9, 10, 81–87, 89–95, 100, 134 Absolute gravity, 3–12, 18, 82, 85, 86, 97–103, 134, 135, 137 Accelerometers, 5, 7, 16, 17, 23, 24, 39–41, 45–49, 52, 54, 82, 116–118, 136, 141, 148 Accuracy analysis, 142–144 Airborne gravimetry, v, 17, 20, 23–30, 39, 51 Altimetry, 127–131, 147 Approximation accuracy, 107–113 Atmosphere loading estimates, 70 Atmospheric tides, 67

#### **B**

Bouguer anomaly, 123

#### **C**

Calibration, 16, 17, 20, 23–29, 40–42, 46, 56, 61, 68 Calibration coefficients, 60–64 Center gravitational gradients, 47–49 Current global geopotential model, 112

#### **D**

Deflection of vertical (DOV), 3–12, 59–65 Degree dispersion, 108, 109, 111 Digital zenith camera, 4

#### **E**

Earth's crust, 31, 35–38, 69 Earth tides, 67–69, 117, 134, 136, 137

#### **F**

Filtering, 4, 23–25, 51–54, 142, 144, 148–150, 152 Filtering and smoothing, 142, 143, 148, 149, 152 Free fall acceleration (FFA), 75–79 Frequency stabilization, 81, 82, 84

#### **G**

Geoid, 18, 107, 108, 127, 130 Geomagnetic field, 75–77 GNSS-R, 127, 129 GRACE, 111, 115–119

81–87, 89–95, 98, 100, 134–138, 141, 143 Gravimetry, v, 15–21, 23, 24, 26, 31, 32, 37, 51, 52, 56, 75, 79, 98, 108, 112, 127, 134–137, 147 Gravity-aided navigation, 141–145 Gravity disturbance, 16, 17, 19–21, 24–28, 55, 56 Gravity disturbance vector (GDV), 51–56, 147–152 Gravity field, 20, 45, 52, 54, 67, 69, 107–110, 112, 115–119, 147 Gravity gradiometer, 39–42, 45–49 Gravity map, 121–125 Gravity measurements, 3, 5, 7, 8, 10, 11, 18, 68, 69, 75–79, 82, 84, 97–103, 133–138, 141, 145 Gravity survey, 3–5, 127 Gravity variations, 17, 18, 67–71, 82, 84, 134 GT-2A gravimeter, 23–30

Gravimeters, v, 3–7, 9–11, 15–18, 21, 23–30, 32, 33, 67–71, 75–79,

#### **I**

Inertial measurement unit (IMU), 16–20 Instrument errors, 39–42 Integrated gravimetric system, 3 Interferometer, 77, 82, 84, 89–95 International Terrestrial Reference Frame (ITRF), 102

#### **K**

Kalman filter, 16, 23–25, 27, 28, 41, 42, 51–54, 143, 148

#### **M**

Magnetic induction, 77 Marine gravimetric studies, 31–34 Measurement methods, 128, 130 Micro-electro-mechanical systems (MEMS), 134, 136–138 Mohoroviciˇ c discontinuity, 31, 35, 37 ´ Moscow, 31, 67–71

#### **N**

Non-tidal gravity variations, 68, 70, 116 Numerical model, 45–49, 116

#### **O**

Observation technique, 59–65

#### **P**

Preliminary measurement processing, 141, 142, 145

© The Author(s) 2023 J. T. Freymueller, L. Sánchez (eds.), *5th Symposium on Terrestrial Gravimetry: Static and Mobile Measurements (TG-SMM 2019)*, International Association of Geodesy Symposia 153, https://doi.org/10.1007/978-3-031-25902-9

#### **R**

Repeat drape lines, 23, 24, 26–30 Resolution, 3, 15, 17, 19, 25, 28, 29, 45, 51–54, 56, 68, 69, 89, 105–113, 115, 116, 127, 134, 136, 137 Rise-and-fall gravimeter, 89, 94, 95 Rotating accelerometer gravity gradiometer, 45–49

#### **S**

Satellite laser ranging (SLR), 97–103, 119 Sea of Japan, 31–38 Sensitivity, 24, 39, 40, 64, 136, 145–152 Shipborne gravimetry, 15–21 Smoothing, 10, 11, 25, 142, 144, 145, 148–150, 152 Spectral characteristics, 107–113

Spherical scaling functions, 51–56 Splines, 23–30 Strapdown, 15–22, 51 Structural density modeling, 31, 35, 36 Systematic error, 17, 24, 25, 51, 52, 54, 81–87, 117

#### **T**

Telecom band laser, 81–84 Time series, 68, 98–103, 118, 135–137

#### **V**

Volcanic observation, 81