Some mathematical assumptions for accurate transformation parameters between WGS84 and Nord Sahara geodetic systems

1 Introduction

3D conformal transformations, also known as similarity transformations (since conformal transformations preserve shape and angles between vectors in space remain unchanged), are commonly used in surveying, photogrammetry, and geodesy. These transformations are necessary for the purposes of cartography, geodesy, and other topographical and cadastral work when the coordinates are generally expressed in the system in effect in a country. The present work consists of developing a strategy of observation, processing double geodetic networks with a view to the comparison between the different transformation models (Bursa–Wolf, Molodensky–Badekas, 2D-multiple regression, etc.), and the realization of a geographical coordinates planimetric conversion grid between the WGS84 system and other national systems like the Algerian North Sahara Datum (Kheloufi et al., 2009).

3D transformations are usually applied to convert coordinates related to one geodetic datum to another; this operation is commonly known as a datum transformation. In such applications, the rotations between the two 3D coordinate axes are small (usually less than 1 s of arc), and certain approximations are used to simplify rotation matrices; these simplified matrices are a common feature of the Bursa–Wolf and the Molodensky–Badekas transformations.

The transformation of GPS data into the geodetic local system, which is the basis of the mapping system in force in every country and in which the results of geodetic work must be expressed, requires knowledge of the passage parameters with the use of a suitable model. The choice of transformation model requires the availability of a dataset known in each of the two systems.

The transformation is presented either in similarity form nor under geographic formalism known as the Molodensky model (Figure 1) or a bi-dimensional 2D affine regression.

Figure 1

Translations and rotations between two systems.

2 Earlier development of transformation formalisms

The names of the two transformations are an acknowledgment to the authors Bursa (1962), Wolf (1963), Molodensky et al. (1960), and Badekas (1969) of technical papers and reports on transformation methods related to the orientation of reference ellipsoids and 3D geodetic networks. 3D conformal transformations for single points are often given in the vectorial form (Collilieux et al., 2016).

The two most commonly used 3D models (Molodensky and Bursa) have the following global simplified formalism:

Molodensky–Badekas:

Bursa–Wolf:

where X → is the coordinate vector in the target system, x → is the coordinate vector in the origin system, x → 0 is the coordinate vector of barycenter in the origin system, T → is the translation vector, ΔR is the variation of the two systems rotation, I is the identity matrix, and Δk is the scale factor variation between metrics of the two systems, generally equal to zero, except in the case of projection from spherical curved coordinates to planar topographical surface.

The mathematical relationship can be written in such a way as to be able to consider the transformation parameters as unknown and thus be able to estimate them. Such a form is that we can estimate the transformation parameters knowing the coordinates in both systems by the following equation (Boateng, 2016):

These two equations could not be obtained without supposing some assumptions for a geodetic mathematical hypothesis in order to lead to these canonical linear formalisms (equation (7)).

In the Molodensky formalism Molodensky (1960), the computation of transformation parameters goes in general through the barycenter coordinates written in both systems (Kheloufi, 2012):

As a result, in matricial form, the Molodensky formalism for a set of n-points is as follows:

where x ¯ i , y ¯ i , and z ¯ i are the coordinates of the isobarycenters in both systems (Novikova, 2020). This model (Molodensky–Badekas) is known as a Cartesian formalism and is obtained by the introduction of a vector and an initial point x ₀ around which (E _x , E _y , E _z ) rotations are performed. This point is generally considered an isobarycenter in both coordinate systems Collier (2020). This model is often used for local, wide-scale transformations (land survey, civil engineering, etc.; equation (7)).

3 Some problems in Molodensky–Badekas and Bursa–Wolf transformations

In reality, the Molodensky transformation is described as shown in Figure 2, where the centroid G is displaced from O2 by translations t _X , t _Y , and t _Z measured in the directions of the X, Y, and Z axes of system 2 and t 2 = t X t Y t Z 2 T is the position of the vector of the centroid. The mathematical relationship between the two systems via this transformation is, in reality, more complex than its final known form and is written as follows (Zavoti and Kalmár, 2016):

Figure 2

Illustration of Molodensky–Badekas barycenter and Bursa–Wolf transformation.

The Figure 2 illustrates graphically the complexity and differences in the implementation of the two 3D transformations relevant to the type of data (point coordinates) and the expected accuracy as a function of the purpose of the transformation.

4 Observations quality analysis in the least squares adjustment

5 Statistical study and results validation

For the validation purpose of this study, a set of 12 benchmarks (points whose coordinates are known in both systems) from the zero-order North Algerian network (Figure 3) are used. Indeed, the results obtained by implementing the 3D Molodensky–Badekas and Bursa–Wolf models show how much the two methods are almost similar but differ slightly in the triangular superior elements of variance–covariance matrix (VCV) (Table 2).

Figure 3

Distribution of zero-order North Algerian network of used benchmarks set for algorithm validation.

The two global methods lead practically to the closed results on residues except that the Molodensky–Badekas formalism gives non-correlated parameters which means that minimal variance is well reached in the case of Molodensky–Badekas than Bursa–Wolf. Indeed the reduction of benchmark coordinates according to the centroid (barycenter) ones mitigates the effect of rotation between the two systems and thus decreases the errors over the geodetic points positions (Hashemi et al., 2013) and even their root mean square (RMS) (Table 1). The VCV matrix in Tables 2 and 3 clearly indicates the difference in precision between Molodensky and Buras–Wolf as displayed in bold, higher variance values indicate less precise results. The normal character denotes less correlation over the transformation parameter errors, thus more precision in the case of Molodensky–Badekas. In order to increase the reliability of transformation parameters, another 2D formalism should be adopted and compared to the 3D ones to carry out the effect of the third dimension on the accuracy of these parameters (Zavoti and Kalmár, 2016).

Among these most commonly used formalisms, we can note multiple regression MREs, affine transformation, and others. For this purpose, we have converted geographic coordinates in both systems to Cartesian ones (Ligas and Banasik, 2011).

Institute of Cartography and Remote Sensing (INCT) undertook the computation of transformation parameters over the same GNSS observation points (12 points in North of Algeria) and found the values reported in Table 1.

The global RMS of these computations, which is approximatively equal to 0.9288 m for translations and 0.143 s for rotations for all components, is considered acceptable for cartographic work. In comparison to the standard deviation of INCT transformation parameters reported in the last column of Table 1, the accuracy reached in our case seems to be satisfactory (the maximum value was 0.03 m for Bursa–Wolf in Z-translation). In the histogram of Figure 5, RMS values of the determined parameters are far less than the values of the adopted parameters (INCT) because, wherever the method is used, it does not exceed 0.25 m. The geographic transformation as “VEIS” formalism is converted to a Cartesian one by the conversion of geographic coordinates (Featherstone, 2008), (George, 1999).

In Table 2, it is well demonstrated that the number of correlation coefficient values greater than half is 19 for the Bursa Wolf (BW) model and six for the Molodensky Badekas (MB) model. It denotes that Molodensky–Badekas formalism Aström (2001), by introducing the barycenter, decreases the correlations between the double points and thus attenuates the propagation of errors on the transformation parameters. If a point is suspect or badly defined, the reduction of its coordinates to those of the barycenter is supposed to decorrelate the errors and thus reduce their propagation.

5.1 Comparison between linear (general least square, GLS) and nonlinear adjustment (LM)

The comparison is not only limited to the use of transformation formalism but also concerns the adjustment technique used in relation to the type of problem to be resolved. Indeed, we are sometimes faced with particular cases of nonlinear systems that need special fitting techniques like “Levenberg–Marquardt,” “Steepest Decent,” “Newton algorithm” and many others (Ruffhead, 2020, 2021).

In this work, I have opted for LM to analyze and compare its purpose with that of the linear method “GLS” as it occurs. As mentioned in the third assumption above, this method does not take into consideration the linearization problem and thus avoids truncation of the equation using the Taylor development until certain order and no information about the system would be lost, thus preserving the accuracy of the results shown in Tables 3 and 4. Considering the Bursa–Wolf transformation as an example, the equation can be written as follows:

(25) I 2 = ( 1 + d s ) ( I + U ) I 1 + t 2 ,

where I 1 = [ X Y Z ] 1 T and I 2 = [ X Y Z ] 2 T ,

Table 4

Accuracies obtained by GLS resolution for two global models with assumptions compared to MRE’s

Model	Component	Minimum (m)	Maximum (m)	Mean (m)	RMS (m)
Bursa–Wolf	X	−1.426	2.686	1.293	0.238
	Y	−0.301	1.164	0.697	0.434
	Z	−1.491	1.307	0.461	1.018
Molodensky–Badekas model	X	−1.417	0.407	−0.841	0.801
	Y	0.238	0.541	0.404	0.354
	Z	−2.033	2.019	−1.308	1.158
Multiple regression equations	Λ	−0.063	0.319	0.110	0.164
	Φ	−0.407	0.634	−0.164	0.366

ds is a very small quantity, and t 2 = [ t X t Y t Z ] 2 T is a vector of translation between the two systems (Deakin, 2006).

The above equation is considered simultaneously for a nonlinear process using the LM algorithm and a linear one illustrated by GLS. In the following, we present the comparison of results between the two fitting methods for both transformations.

The Table 4 present the results of statistical analysis for the global models and the multiple regression MRE’s algorithm. Table 3 concerns the generalized least squares, and Table 4 shows the nonlinear fitting of the global (3D) models. Globally, the RMS turn around the meters in the worst case (>2.058 for Molodensky–Badekas) and less than 20 cm for the MRE’s model. This denotes the possibility of improving the precision of position in the optimum case; geodetic positioning sometimes needs sub-centimeter (geodynamics and technical auscultation) as accuracy to be satisfied. In Table 5, was reported the same models processed by LM algorithm with resolving the equation (12), δ _X is thus a solution of the nonlinear system applied for all models. For this purpose, all models and both linear and nonlinear algorithms have been implemented in the « Delphi10 » application (Figure 4). A GLS solution has been introduced to determine a unique set of geodetic coordinates, with accompanying accuracy prediction coordinates from different observing scenarios (Soler et al., 2012).

Figure 4

Screenshot of the application interface including both techniques.

According to the Table 5, the LM algorithm improves the accuracy of the global models if we take into consideration different assumptions (Acar, 2006). In this case, there is no need to appeal the multiple regression method because, for certain applications (geodetic works), we need to transform 3D cartesian components (X, Y, Z).

Nevertheless, the MRE remains interesting for limited area transformation in the case of cadastral surveys. In Figure 4, the choice of the degree of MRE is decisive for better accuracy; in our example, the choice of degree two in the interface (red framed degrees choices in Figure 4) is more suitable to obtain less error in transformation calculus (Figure 5), whereas degree three is not numerically stable because of the minus sign in the regression factors.

Figure 5

Variation of parameters RMS in meters for the used techniques (global models adjusted with GLS) and MREs.

We could improve the quality of transformation parameters issued by GLS by Zoning fitting and considering a zonage for a wide area, i.e., processing the method by set of benchmarks belonging to several restricted areas. In this case, we have partially resolved the problem because the obtained parameters are regional and not global, which is useful for a restricted area (Hok Sum, 2003).

In Tables 3–5, we present the results of data processing by both 3D and 2D models (Bursa–Wolf, Molodensky–Badekas, and MRE) with and without the assumptions explained above for carrying out the contribution of each model and also the contribution of assumptions taken into account in the improvement of precision in the processing of geodesic networks.

Table 5

Results and analysis of the nonlinear (LM) process over the two global models

Model	Component	Minimum (m)	Maximum (m)	Mean (m)	RMS (m)
Bursa–Wolf	X	0.126	0.686	1.293	0.213
	Y	−0.101	0.164	0.697	0.762
	Z	0.572	1.102	0.461	1.015
Molodensky–Badekas model	X	−1.223	0.311	−0.841	0.906
	Y	0.103	0.264	0.404	0.325
	Z	−1.695	1.019	−1.308	1.065

For the global models, it is clear from tales that the three considered assumptions have brought improvements in accuracy (Table 4) better than the improvements brought about by the implementation of nonlinear least squares analysis “LM”. Indeed, the precisions revolve around 35 cm with Molodensky–Badekas and 40 cm with Bursa–Wolf once the assumption is considered (Table 4) compared to the results of Table 3, where these hypotheses have been omitted (Yang, 2014).

The results of Tables 3 and 4 show that 2D model MRE minimizes errors, and even its statistical analysis is very good considering the tests undertaken on it.

6 Discussion

After considering these assumptions, we can say that we are in conformity with the directives and resolutions of the “UN-GGIM”¹ concerning the use of the international geodetic frames and attachment of our local systems to the latest versions of ITRF to minimize position errors during compensations because the difference between the ITRF2018 and the WGS84 is some centimeters. The global 3D models processed with OLS seem to be the most used and famous for this purpose. Nevertheless, in different fields of engineering and earth science, certain cases need more accuracy; the ordinary linear least squares prove to be limited.

Thus, we have to use new numerical methods of resolution that can provide great efficiency in polynomial modelization. Therefore, instead of being limited to linear models, we have to apply the nonlinear least squares resolution to resolve the transformation problem between geodetic systems.

This need appears especially in the case of the Nord Sahara datum (Algeria), on which the linear models are not very appropriate because of the lack of information about the geoid’s undulation (the third component is not accurately determined). In this article, our main aim is to carry out the importance of using the nonlinear least squares method (LM) applied over two 3D (global transformation) models and a 2D one (MRE) on a huge area benchmarks network (North Algeria).

Coordinate transformation parameters and their RMS have been computed by both the OLS and LM algorithms compared with multiple regression (Kheloufi et al., 2009) to compare, on the one hand, the linear adjustment within its two variants (local and global) with the MRE’s and, on the other hand, the linear adjustment and the nonlinear one. In this context, a set of 12 benchmarks localized in North Algeria have been integrated to compute the transformation parameters (3D and 2D). Tables 3, 4, and 5 show the difference in accuracies between different models (BW, MB, and MRE’s) processed by the same method (GLS) and also the difference between the methods themselves, i.e., GLS vs. LM. For example, the Molodensky formalism (Deakin, 2006) generally gave better precision when processed by GLS for the X and Y components (1.014 m against 1.838 m) except for the Y component, where Bursa–Wolf is near the meter when Molodensky–Badekas is about 2 m. Multiple regression, which is a 2D formalism, ensures sub-meter accuracy and is adapted for wide-area transformation processing like in North Algeria. In our case, for a near 500 km baseline length, the method turned out to be very efficient because it did not exceed a meter. For the global models, we could reach these performances with precision only because we have taken into consideration the assumptions cited above (assumptions 1, 2, and 3).

Previous studies that have been performed without including a hypothesis have given error-blemished results with less accuracy (Table 5). For both global models, we could not reach the submeter for several components, where even the X and Y components outperformed the geodetic expected errors.

The purpose of this work was to develop a solution to the crucial problem of transformations between geodetic reference frames, particularly between the system in force in Algeria, the North Sahara – Clarke 1880, and the system in which the coordinates resulting from the GPS campaigns are expressed, i.e., the WGS84. Through this multicriteria analysis, we could prove that the coordinate transformation problem could be solved to carry out reliable transformations within acceptable precisions for different applications (geodetic, geodynamic, and land survey). The calculation of these parameters required an investigation into the existence of double-point files of impeccable quality and homogeneous and dense exhaustiveness. This led us to use some files from previous undertaken by different institutions (centre des techniques spatiales, INCT, entreprise nationale de geophysique) and in different projects: PUBLIC WORKS DIVISION-DTP-Oran network, Reseau geodesique de Sonatrach-Sonatrach project, North Algeria densification campaign, etc., where the set of points are determined with a priori centimeter precision standard deviations.^[1]

As the determination of the transformation parameters is also based on the numerical models and mathematical formalisms, our work focused on the treatment of double points in our disposition by different methods:

• 3D model of Molodensky–Badekas and Bursa–Wolf,

• 2D model of multiple regression equations,

• Nonlinear LM algorithm confronted with LLS.

The results presented in this study highlight the strengths and weaknesses of each model on the basis of the parameter standard deviation and even the position RMS, as well as the residuals after transformations, as illustrated in tables and figures. Indeed, according to the obtained results, the comparison between the 3D models gives an advantage for the Molodensky–Badekas formalism compared to the Bursa–Wolf formalism and the advantage for the nonlinear adjustment at the expense of the linear one. The analysis made on the estimated parameters also shows a difference in precision between the 3D models and the multiple regression method, which gave low RMS (as shown in Figure 5, 2 cm for MREs and 4.5 cm for Bursa–Wolf).

As a perspective, we recommend the use of planimetric grid method with interpolation methods (plus proche voisin, bilinear, and kriging) for more accuracy improvement.