DEVELOPMENT OF A LOCAL GEOID MODEL AT THE FEDERAL DISTRICT , BRAZIL , PATCH BY THE REMOVE-COMPUTE-RESTORE TECHNIQUE , FOLLOWING HELMERT ' S CONDENSATION METHOD Desenvolvimento

There are several techniques for determining geoid heights using ground gravity data, the geopotential models, the astro-geodetic components or a combination of them. Among the techniques used, the Remove-Compute-Restore (RCR) technique has been widely applied for the accurate determination of the geoid heights. This technique takes into account short, medium and long wavelength components derived from the elevation data obtained from Digital Terrain Models (DTM), ground gravity data and global geopotential models, respectively. This technique can be applied after adopting the procedures to compute gravity anomalies and, then, the geoid model, considering the integration of different wavelengths mentioned, and their compatibility with the vertical datum adopted. Thus, this paper presents the procedures, involving the RCR technique, following Helmert's condensation method, and its application to compute one local geoid model for the Federal District, Brazil. As a result, the local geoid model computed for the studied area was consistent with the available values of geoid heights derived from geometrical levelling technique supported by GNSS positioning.


Introduction
Height determination and vertical control with a precise geoid model constitutes one of the most challenging research subjects of geodesy, and it attracts more attention since 1980s, related with the wide spread and intensive use of GNSS techniques in surveying (Erol and Erol 2013).According to Sjoberg (2005) and Hirt (2011), many strategies used in gravity field modeling were developed at a time when the precision goal to determine the geoid height was 10 cm or less.Currently, according to Hirt (2011), to determine the geoid and quasi-geoid heights with an precision of centimeters or better, it is necessary to evaluate carefully, and, if necessary, correct the approaches that are inherent to the methods and the techniques used.
There are several methods for determining geoid heights using groung gravity data, the geopotential models, the astro-geodetic components or a combination of them.Among the techniques used to determine the geoid models using the gravity data at regional level, the bestknown approach in the literature is the RCR, according to Schwarz et al. (1990) and Abbak et al. (2012).This approach has been used in many parts of the world, and among them Canada, Turkey, Austria, United States, Australia and Brazil (Schwarz et al. 1990;Ayhan 1993;Zhang et al. 1998;Fotopoulos et al. 1999;Smith and Small 1999;Featherstone et al. 2004;Abbak et al. 2012;Blitzkow et al. 2012).
The RCR technique, according Sansò and Sideris (2013), takes into account the short, medium and long wavelength components that are derived from the elevation data obtained from Digital Terrain Models (DTM), ground gravity data and global geopotential models, respectively.This technique requires adopting procedures to compute gravity anomalies and then of the geoid model, considering the integration of the different wavelengths mentioned, and their compatibility to the vertical datum adopted.
Given the above, this work presents the procedures, involving the RCR technique, following Helmert's condensation method, and its application to compute one local geoid model for the Federal District, Brazil.The motivation for this work is due to cities development within the Federal District occur in flat areas with several infrastructure problems, such as water supply and drainage of rainwater and sewage, which demand accurate knowledge of orthometric height to solve them.

RCR Approach
The RCR technique for calculating the geoid model can be divided in three distinct stages.The first is the removal of the long wavelength component of the gravity anomaly generated by Helmert's second condensation method (  HEL g To develop the technique,  GM g and GM N can be estimated according to Smith (1998), using the geopotential coefficients adopted to a pre-set degree, to comprise only the long wavelength components.

     
To calculate  GM g and GM N , it is also necessary to subtract the fully normalized spherical harmonic coefficients of the gravitational potential of the coefficients implicit in the reference Bull.Geod.Sci, Articles section, Curitiba, v. 23, n°3, p.520-538, Jul -Sept, 2017.ellipsoid.This is done by the zonal spherical harmonic coefficients of the gravitational potential ( 2 ,0 n C ), according to Moritz (1984) and Smith (1998).Thus: where and 2 J is calculated as demonstrated by Cook (1959): GM ,  and f are the geocentric gravitational constant, angular velocity and the flattening of the reference ellipsoid, respectively.
For all other coefficients, it is assumed: , nm C and , nm S are the fully normalized spherical harmonic coefficients of the gravitational potential.
According to Blitzkow (1986), the equations 13, 17 and 18 represent, generically, the relationship between the coefficients linked to disturbing and gravitational potentials.In practice, the aforementioned equations show that the gravitational potential of the normal earth use only 0  m and n pair, and that does not contain terms which depend of Equations 3 and 4 do not consider the zero degree term in gravity anomaly ( 0 g ) and co-geoid ( 0 N ).Therefore, to compute  GM g and GM N considering a reference ellipsoid adopted, this term must be added on the equations 3 and 4, respectively.According to Kirby and Featherstone (1997), the degree zero term may be computed by: Development of... Bull.Geod.Sci, Articles section, Curitiba, v. 23, n°3, p.520-538, Jul -Sept, 2017.
W is the gravity potential on the surface of the geoid.
U is the normal gravity potential on the surface of the normal ellipsoid and may be computed by: tg e ab According Vaníček and Kleusberg (1987), where according Wong and Gore (1969): L is the maximum degree,   The discrete calculation of N presents a singularity when 0   .To work around this problem, Sideris and She (1995) proposed: Where 0 s is the radius of the next considered area.
Then, the calculation of the geoid model using the Stokes discrete formula is given by In the RCR technique, To calculate  RES g , it is necessary to find the gravity anomaly.The second condensation method of Helmert (  HEL g ) is the most often used because it produces the small indirect effect of topography (Heiskanen and Moritz 1985).
 FA g is the free-air anomaly, ATM C the atmospheric correction, T C is the terrain correction and  g is the indirect effect of topography (Heiskanen and Moritz, 1967), also known as the indirect secondary effect.
0.8658 9.727 10 3.482 10 obs g is the gravity observed on the physical surface of the Earth. FA g is calculated according to Featherstone and Dentith (1997), ATM C is calculated according to Kuroishi (1995), in mGal.x , y and z represent the planar coordinates and orthometric heights of the integration points and of the computation points ( p ).
The IE N is calculated also using the planar approach, according Wichiencharoen (1982).
g and estimate  HEL g , the height data are extracted from a previously defined DTM.This is necessary to eliminate the differences in the height values determined by different source data.
To estimate N , the  HEL g values were interpolated to generate a regular grid and enable the operations using the RCR technique.The inverse distance squared was used as the interpolation method.In general, first, the Bouguer correction (

Adjustment of the Gravimetric Geoid Model
The geoid height computed using gravity data can be evaluated by comparing To perform the evaluation, it is necessary to make the geoid height computed using gravity data compatible with the vertical reference datum location.As described by Sansò and Sideris (2013), the RCR technique refers to the geocentric reference system implicit in the geopotential model used.Also, the local levelling datum to which the orthometric heights refer will not likely correspond to the reference potential value of the geopotential model or the GPS reference system.
To solve this problem, it is necessary to combine the heterogeneous height data.
where, A is the design matrix; 0 X , vector of initial parameters; a X , vector of adjusted parameters; X , vector of corrections; b L , vector of observed values ( N ); and V , residue vector.
Among the functional models adopted, we have the classical four-parameter linear model presented by Sanso and Sideris (2013): N a bcos cos ccos sin dsin where b , c and d are the translation parameters; a is the change of the reference value of the potential; and  i and  i are Latitude and longitude of the GNSS/Levelling points.
After compute the parameters by LSM, the obtained geoid height is compatible with the local vertical datum adopted.
Besides making the vertical data compatible, it is correct to affirm that the LSM using the parametric model also takes into account: the random errors derived from N , h and H ; systematic effects and distortions of height data; theoretical assumptions and approximations made when processing the observed data; and the instability of the monument of the reference station over time, for example.

Evaluation of the local geoid model
The local geoid models computed can be evaluated on two ways, as presented by Tocho et al. (2013).
The first involves descriptive statistics of the absolute differences between the geoid heights ( N ) extracted from the computed geoid models ( N ) and from GNSS/levelling ( / GNSS Lev N ) points.Those differences can be expressed by: i is the point used in the evaluation.
The second involves descriptive statistics of the relative geoid heights differences (  rel N ) formed for the baselines computed from pairs of points, using the computed geoid model and the GNSS / levelling points.It can be shown by: Development of... Bull.Geod.Sci, Articles section, Curitiba, v. 23, n°3, p.520-538, Jul -Sept, 2017.
i and j are the points used to form the baseline in the evaluation.ij S is the baseline distance.
If the value of N is in mm and the value of S is in km,  rel N has the value in ppm.Bull.Geod.Sci, Articles section, Curitiba, v. 23, n°3, p.520-538, Jul -Sept, 2017.

A local geoid model at the Federal District, Brazil
The region of the Federal District of Brazil was chosen to compute the local geoid model (LGM).This region is located between 48.25ºW and 47.33ºW and 16.06ºS and 15.45ºS (Figure 2), with a slightly wavy relief, ranging from 600 to 1340 meters above sea level.
To compute a local geoid model at the study region, the following material was used: -GECO (Goce and Egm2008 COmbination) geopotential model (GGM), developed by Gilardoni et al. (2015) and made available by the ICGEM.GECO was chosen because it is the newest highest resolution geopotencial model based on the integration of the EGM2008 (Earth Gravitational Model 2008) and of the GOCE (Gravity field and steady-state Ocean Circulation Explorer) satellite tracking data (fifth release of the time-wise GOCE solution).
-DTM from the Shuttle Radar Topography Mission (SRTM), with spatial resolution of 90 m.The DTM are located between 49.75ºW and 45.83ºW and 17.56ºS and 13.95ºS (Figure 2).
-2312 ground gravity stations (GGS) provided by the Brazilian Institute of Geography and Statistics -IBGE, the National Petroleum Agency -ANP, including 323 new stations acquired by the authors.In addition, to complete the ground gravity data in regions without ground gravity stations, GECO up to degree and order 2190 was used.In this case, the ground gravity data were computed using the gravity anomaly and the height data extracted by ETOPO1 model, provided by National Oceanic and Atmospheric Administration -NOAA.All of the gravity data are located between 49.25ºW and 46.33ºW and 17.06ºS and 14.45ºS (Figure 2).
To analyze the results and to adjust the local geoid model to the local vertical datum (Imbituba vertical Datum), 24 points whose geoid heights were obtained by GNSS positioning and geometric levelling, provided by IBGE, were used.GECO has contribution of GOCE data up to the 250 degree.However, the choice of the degree up to 360 to compute  GM g and GM N on this was made because it presented less dispersion of the differences between GM N and / GNSS Lev N .The constants used in this work are shown in Table 1.These constants are presented by Moritz (1984) and IERS (International Earth Rotation and Reference Systems Service) Technical Note (Petit and Luzum, 2010).The local geoid model (Equation 2 and Figure 6) with a spatial resolution of 2.5' was obtained following the calculation of GM N (Equation 4), RES N (Equations 22 and 28), the zero degree term (Equations 19 and 20) and IE N (Equation 34).The zero degree term computed and added in  GM g and GM N was -0.152 mGal and -0.442 m, respectively.
The LSM (Equations 38 and 39) was used to adjust the local geoid model to the local vertical datum, using as reference 24 points whose geoid heights were obtained by GNSS and geometrical levelling (Figure 6).The altimetric precision of the points used as reference is approximately 0.073 m.As there are only a few points to apply the SLM, and that the lack of one of them can affect considerably the results of the adjustment, this study did not include part of these points as checkpoints.After adjusting the parameters (Equation 39), the local geoid model was estimated (Equation 40and Figure 8) free of systematic components - a N .The systematic components are presented on   To evaluate the results, it was analyzed the absolute and relative differences between the geoid heights extracted from Grav N and / GNSS Lev N .In both analyzes, the official geoid model adopted in Brazil (MAPGEO2015) was included to verify the performance of this work computed models.
Although this study did not include checkpoints to analyze the Final N , the residual value for the reference points extracted from the LSM was used too.
Table 2 and Figure 9 shown the descriptive statistics of the absolute differences between the geoid heights of the local geoid models and the geoid heights computed from the 24 GNSS/levelling points (Equation 41).It can be seen that the Quartile Coefficient of kurtosis is similar for all the models, but the Grav N (Figure 6) and second eccentricity of the reference ellipsoid.The RES N is calculated based on the principle of Stokes (Stokes 1849), which allows to estimate the values of the geoid height ( N ) using the gravity anomaly values ( g ) obtained on the physical surface of the Earth, considered as spherical.In the discrete form of the surface elements, N becomes (Sideris and She 1995): gravity anomaly of an area in a grid with n parallels and m meridians;   and   are the variations in geodetic coordinates, latitude and longitude, which comprises each area; '  and '  are the geodetic coordinates at the center of the area; 0  is the average normal gravity of each area;  is the spherical distance between two points; and    M S is the modified Stokes function, used to remove the low-degree terms of the Legendre polynomials from the    S (original Stokes function).

BC
) is added to each point ( p ) for which  HEL g has been calculated, followed by the interpolation of values for the points of the regular grid.Finally, B C is eliminated from the generated grid, thus yielding the Helmert anomaly estimated for the regular grid (  Grid HEL g ).The values of B C are used to smooth the values of  geometric altitudes ( h ), determined by GNSS positioning techniques and ( H ) orthometric heights, determined by geometric levelling taking as origin the local vertical reference datum.

Figure 1
Figure1shows the flowchart of computations used to estimate the local geoid model according to equations presented to implement the RCR technique and to adapt the geoid height.The flowchart includes a whole set of routines developed for reading the input data, for computation procedures and results generation.All routines, here called GRAVTool, were implemented based on the MATLAB® software.As shown in Figure1, the input data used for the calculation procedures include: global geopotential model (*.gfc), provided by the International Centre for Global Earth Models -ICGEM; DTM image (*.tiff and *.tfw); ground gravity data in ASCII (*.txt), containing the geodetic coordinates, orthometric height and observed gravity of used points; constants related to the reference ellipsoid and average density; and terrestrial data that originated from the GNSS positioning and geometric levelling, in ASCII (*.txt) containing geodetic coordinates, and the geometric and orthometric altitude of each point.All gravity anomaly and geoid height results, calculated using the equations shown in previous sections, are available as ASCII (*.txt).

Figure 1 :
Figure 1: Flowchart of the sequence of calculations to estimate the gravimetric geoid model.1) input data, and 2) sequence of calculations and output results.

Figure 2 :
Figure 2: Spatial distribution of the ground gravity stations (GGS) and ground gravity data computed from GECO geopotential model (GGM), tide free, up to degree and order 2190.Boundaries of the DTM, GGM, GGS, LGM and states are presented, too.

Figure
Figure 3:  GM g calculated for the study area, using Figure 4: GM N calculated for the study area, using

Figure
Figure 5:  HEL g calculated using the ground gravity data and DTM of the study area.The black polygon shows the geographical boundaries of the Federal District.

Figure 6 :
Figure 6: Calculation of the local geoid model ( Grav N ).The red points are geometrical levelling technique associated with GNSS positioning used to evaluate and to adjust Grav N to the local vertical datum in the study area.The black polygon shows the geographical boundaries of the Federal District.

Figure 7 :
Figure 7:  a N (systematic component) adjusted by the LSM, using as reference points whose geoid heights were obtained by GNSS positioning and geometric levelling for the study area.The black polygon shows the geographical boundaries of the Federal District.

Figure 8 :
Figure 8: N after applying  a N in the study area.The black polygon shows the geographical boundaries of the Federal District.
, the Grav N presented more symmetric than the other models and the discrepancies of the differences (  maximum minimumdifferences ) of the Grav N (0.254 m) and Final N (0.251 are similar. 

Table 1 :
Constants values used to compute the geoid model.

Table 2 :
Descriptive statistics of the differences between geoid heights from different models (