On the determination of a locally optimized Ellipsoidal model of the Geoid surface in sea areas

Deriving a local geoid model has drawn much research interest in the last decade, in an endeavour to minimize the errors in orthometric heights calculations, inherited by the use of global geoid reference models. In most parts of the earth, the local geoid surface may be tens of meters away from the Global Reference biaxial Ellipsoid (WGS84), which create numerus problems in topographic, environmental and navigational applications. Several methods have been developed for optimizing the precision of the calculation of the geoid heights undulations and the accuracy of the corresponding orthometric heights calculations. The optimization refers either to the method used for data acquisition, or to the geometrical method used for the determination of the best fit local geoid model. In the present work, we focus on the reference ellipsoid used for the geometric and geoid heights determination and develop a method to provide the one that fits best to the local geoid surface. Moreover, we consider relatively small sea regions and near to coast areas, where the usual methods for data acquisition fail more or less, and we pay attention in two directions: To obtain accurate measured data and to have the best possible reference ellipsoid for the area at hand. In this due, we use the “GNSS-on-boat” methodology to obtain direct sea level data, which we induce in a Moore Penrose pseudoinverse procedure to calculate the best fit triaxial ellipsoid. This locally optimized reference ellipsoid minimizes the geometric heights in the region at hand. The method is applied in two closed sea areas in Greece, namely Corinthian and Patra’s gulf and also in four regions in the Ionian Sea, which exhibit significant geoid alterations. Taking into account all factors of uncertainty, the precision of the mean sea level surface, produced by the “GNSS on boat” methodology, had been estimated at 5.43 cm for the gulf of Patras, at 3.76 cm for the Corinthian gulf and at 3.31 for the Ionian and Adriatic Sea areas. The average difference of this surface and the local triaxial reference ellipsoid, calculated in this work, is found to be less than 15 cm, whereas the corresponding difference with respect to WGS84 is of the order of 30m.


Introduction
The geoid is an equipotential surface of the Earth's gravitational field which is approximated by the global Mean Sea Level (MSL), with high accuracy in sea areas if other influences (winds, tides) were absent and theoretical extension through the continents with hypothetical canals. The determination of the geoid surface is a fundamental problem in geodesy and due to its complexity and the lack of accurate data for the mass distribution of the earth, it is an open scientific challenge, the so-called geoid problem [1]. Modelling the geoid with an ideal geometric solid, suitable for geodetic use, has led to different geometrical modelling options, ranging from sphere to spheroids [2]. Τhe predominant geometric model of the geoid is the World Geodetic System 1984 (WGS84) [2,3], which is an Ellipsoid of Revolution, a global geocentric system on which the GPS navigation systems, used to geographically locate points on the Earth's surface, are based.
Since WGS84 is an ideal surface that minimizes the mean square geometric height of all the points on the Earth's surface, it inevitably introduces locally significant deviations in places where the Earth's relief has peculiarities. This applies to areas both on land and coastlines or areas of submarine faults where the MSL is quite different from the geoid. These deviations introduce errors in GPS navigation systems, showing parts of the Earth at different heights than they really are, which might have crucial impact in certain applications [3]. Therefore, the investigation for a high-precision local geoid model has drawn much scientific interest in high-reliability navigation systems technology, see for example [4,5] and references therein.
In this regard, most scientists focus on approximating the local geoid surface based on methods that exploit known geoid heights N, or geoid undulation differences ΔN, or orthometric and ellipsoidal heights, always with respect to the global reference ellipsoid WGS84. The optimized geoid surface is then calculated by an interpolation method and Least Square technique [4][5][6].
The present study focuses on locally optimizing the geometric model of the geoid, i.e. the reference surface with respect to which, the geoid heights and the geometric heights are measured. Moreover, we focus on closed sea areas or sea coast areas, where the geoid is well described by the sea surface while other methods (gravimetric, astrogeodetic or satellite techniques) used for geoid description in such regions exhibit low accuracy, due to large geoid alterations.
In particular, we seek for a triaxial ellipsoidal surface, denoted by Ɛ in this context, which minimizes the geometric heights in regions that exhibit special geodetic interest. This triaxial ellipsoid Ɛ is assumed to have the same center, equatorial plane and polar axis as WGS84, as these refer to the geocentric, with scientifically acceptable accuracy and with assumptions that are widely studied. The accuracy of the approximation depends on the accuracy of the data used for the geometric heights. In this due, the geographical coordinates of the points in the region of interest are provided by on-the-spot field measurements using a pioneering recording system, known as the GNSS-on-boat method [7][8][9][10][11], which provides very high accuracy especially in seaside or closed sea regions. With such accurate point measurements available, a Moore-Penrose analytical procedure is used to calculate the best fit ellipsoidal surface Ɛ and then, the corresponding geometric heights are analytically calculated with respect to Ɛ.
The method is applied to six areas with geodetic interest, two closed sea areas in Greece, namely Corinthian and Patra's gulf and also in four regions in the Ionian and Adriatic Sea, which exhibit significant geoid alterations. The results suggest an improvement in the local approximation of the geoid of the order of a few centimeters compared to the tens of meters provided by WGS84. The paper is structured as follows. In Section 2 the methodology of the work is discussed, with regard both to the GNSS-on-boat method and to the mathematical optimization calculations. Section 3 provides the implementation of the method in the six regions mentioned above, with the mathematical technicalities presented analytically for the first region of interest, namely the Patras gulf, while the corresponding results for the rest five regions are also included in Section 3. In Section 4 we discuss the results and conclusions of the work.

Α. The method of GNSS-on-boat
The GNSS-on-boat technique [12 -17] based on kinematic GNSS methodology [18] refers to the simultaneous satellite tracking by a number of moving (rover) GNSS receivers and by one or more stable (base) receivers on ground. The apparent changes of the base receiver position permit corrections to the computed coordinates of the rover receiver. In GNSS-on-boat methodology, the rover GNSS (GPS) receiver directly calculates the instantaneous elevation of the boat, relatively to the reference ellipsoid. Filtering of the short-term effects (sea-level fluctuations shown as two different levels) permits an accurate estimation of the sea surface topography (SST) along the boat tracks.
Other remote sensing techniques such as airborne altimetry (AA) use two or more sensors: a sensor (usually GNSS) to compute the geometric height of the plane, plus an altimeter to compute the elevation difference between the seas and the plane [19]. In those techniques the estimated accuracy is lower than the accuracy of GNSS-on-boat, since rover receiver records directly the sea surface heights. (Figure 1) In the present work, GNSS results from three research missions and six selected areas were used: in the Gulf of Patras, in the Gulf of Corinth and in the Ionian and Adriatic Sea.

Β. Determination of the best fit ellipsoidal in a region
We consider a geocentric reference coordinate system Oxyz with polar axis parallel to the axis of rotation of the Earth and with Earth's equatorial plane as the Oxy plane. Let a triaxial ellipsoid centered at the origin with semiaxes a, b, c defined by equation that is satisfied by the Cartesian coordinates, denoted by (x i , y i , z i ), i = 1, 2, …, k , of all measurement points in the area of interest. The following system is solved with respect to the unknown semiaxes Since (2) is overdetermined, while it is applied in a grid of more than 50 points, classical Moore -Penrose pseudoinverse procedure is applied to obtain the least square solution of (2), a method rather common in such problems [20,21], in the form Provided that the inverse matrix exists and introducing the suitable notation , the solution reads as Hence, equations (4)-(6) provide the triaxial ellipsoid Ɛ, centered and oriented as EGS84 with semiaxes adjusted so that it best fits to the region under investigation.
C. Calculation of the geometric height with respect to the best fit triaxial ellipsoid Next, we turn to calculate the distance of every given point i = (x i , y i , z i ), i = 1, 2, …, k, from the ellipsoid Ɛ, which is actually the geometric height h Ɛ,i , of the point i with respect to Ɛ. In particular, we calculate the minimum distance of each point ri from the tangent plane of the ellipsoid Ɛ, at the point d,i = (x d,i , y d,i , z d,i ), i = 1, 2, …, k, on Ɛ [22]. In order to calculate h Ɛ,i , we first find the corresponding d,i .
Since the distance function between two points in three-dimensional Euclidean space is: , for i = 1, 2, …, k, the point d,i is the one that minimizes the function F, under the condition that it belongs to the surface Ɛ, namely: Using the method of Lagrange Multipliers the Cartesian coordinates of d,i = (x d,i , y d,i , z d,i ) are given as The parameter λ i is calculated by solving the following system which produces a sixth degree equation with respect to λi: The real root λ i that corresponds to the minimum value of F is introduced in (9) to provide the point of tangency d,i , suitable for calculating the distance of point i from the ellipsoid Ɛ, i.e. the geometric height h Ɛ,i Then, the corresponding geometric height h Ɛ,i for every point ri , i = 1, 2, …, k is given in the form: Finally, the mean geometric height is calculated for all points of each region.

D. Comparison of Triaxial Ellipsoid E and Ellipsoid of Revolution WGS84
The final step is to estimate the contribution of the local triaxial ellipsoidal model calculated above, to the geometric description of the geoid. In this due, we compare the mean square heights resulting from the above calculations, with the corresponding high altitude area provided by WGS84, for the data collected by the GNSS-on-boat methodology, for the six regions mentioned and are discussed in detail in section 3. The range of the Earth's surface area in which the calculated triaxial ellipsoid is satisfactory as a local standard and, therefore, as an optimal reference surface for the orthometric height determination of the points of this region is estimated.

Implementation and Results
As mentioned already, the selected areas we used in the present work were part of the Gulf of Patras, part of the Corinthian Gulf as well as four areas of the Ionian Sea and the Adriatic Sea (A, B, C and D). These areas belong to three different research missions [8,9,11]. Figure 2 illustrates these areas on the location map.  Α. Gulf of Patras First we describe the highlights of the method for the case of the Gulf of Patras, where we had 50 data points in Cartesian coordinates, derived by the GNSS-on-boat methodology [7,8] and their geometric heights with respect to the WGS84. Their geographical coordinates range between φ = 38.29 , λ = 21.5 northwest and φ = 38.21, λ = 21.68 southeast of gulf of Patras, distributed on a grid with degree 2 cm [5]. Substituting these coordinates in (1) results to the algebraic system (3) with 50 equations and three unknowns. Applying the pseudoinverse of Moore -Penrose we obtain the three semiaxes of the best fit ellipsoid ƐPat for this region, via (4)-(6), a = 6.3779710 6 m, b = 6.3790710 6 m, c = 6.3568510 6 m.
The differences between semiaxes of the triaxial ellipsoid ƐPat and the Ellipsoid of Revolution WGS84 are: The mean geometric height of the considered points in Patras Gulf with respect to WGS84 is calculated from given data Also, the mean geometric heights of the same points with respect to the Ellipsoid ƐPat is calculated by means of equation (12) In the next subsection we obtain the corresponding results for the five more regions investigated.

B. Corinth gulf, Ionian and Adriatic Sea areas
For each of the other five regions located in Corinth gulf, the Ionian and the Adriatic Sea (see: Cor, A, B, C and D in Figure 2.) we used a number of points which we received from previous research with the GNSSοn-boat method [7,9,11] and we applied the method developed in the present work. The results are presented in the following tables. The three semiaxes of best fit Ellipsoid are presented in the table below. Αs we observe the differences in the semiaxes from the WGS84 model range from 20 meters to 22 kilometers. While, of course, the differences between them are much smaller as they are neighboring areas.
But, how well the best fit Ellipsoid describes the different sea areas? Table 2 presents the mean geometric heights of the mean sea level, with respect to each triaxial ellipsoid calculated above to best fit to each region. As we can see, the best fit Ellipsoid for each area limits the average geometric heights to 1 or maximum 15 cm. (central diagonal of the table). These values are vastly lower than the 25 to 37 meters of the WGS84. Thus, we have an excellent local description that is practically identical to the mean sea level. Of course, the best fit Ellipsoid of one area is not ideal for another area, no matter how close it is. For example, if we take as a reference the best fit Ellipsoid of area A, the average geometric heights of area D reach 11.3 meters. This is expected as the areas are miles away and the corresponding geoid endulation is significant also. The range of the validity of the local approximation with the best fit Ellipsoid calculated by means of the proposed method depends on the accuracy we want, as discussed in the next subsection. C. Combining data from different regions We turn now to apply our calculations for larger regions, by combining data from different regions. So first, based on the best fitting Ellipsoid of the Gulf of Patras, we used the data of the two bays together and then combinations from areas a, b, c and d. The summary results are presented in the following Table 3. In Figure 3, the mean geometric heights with respect to ƐPat for the above regions combinations are presented in a summary diagram. As expected, one can see the increase in the mean geometrical height as the data points depart from the area tangent to the best fit ellipsoid. For example, in the case of the combination of the two Gulfs where the maximum distance between them is 110 Km, the mean value of hƐ is limited to 1.07m, while the mean geometric height calculated for all regions, from Corinthos Gulf to area D in Adriatic Sea, with maximum distance 455 km, the mean geometric height hƐ increase to 10.62m. This diagram can be used to determine the range of the validity of the best fit ellipsoid. For example, if one defines the acceptable accuracy of the model up to 1m, then the validity of ƐPat as best reference ellipsoid should be restricted to the two Gulfs, while if 3m are also acceptable for geometric heights, then the validity of ƐPat can be extended to include Ionian regions A and B as well.

Conclusions
In the present study, we developed a method to determine the triaxial ellipsoidal model Ɛ of the geoid that best locally fits to sea areas, especially closed seas and near coast regions. The method was applied in specific areas, where accurate data measurements where available, using the GNSS-on-boat methodology. The results show that the geometric heights measured with respect to the best fit ellipsoid provided by the proposed method for each region are of the order of a few centimeters. Compared to the corresponding heights with respect to WGS84 which is of the order of 25 to 30m, these results suggest that the method provides good agreement of the geometric surface Ɛ and the sea level. The creation of a mathematical model, such as the Ɛ, which is locally identified with the geoid, is multifacetedly important as it enables the determination of real heights and elaboration of detailed geophysical studies, as well as knowledge of its micro-and macro-disturbances of the geoid depict geological changes in the earth's solid crust (such as underwater faults, deposits, etc.), which remain largely unknown.