Methods of Modelling the Distortion Caused by Different Amount and Orientation of Coordinate Corrections between Two Coordinate Systems at Various Locations

The paper explains the ways of modelling a distortion as a phenomenon caused by different amount and orientation of the residuals-coordinate corrections between initial and target coordinate system at different locations and using the lattice grid models. In order to show more clearly the topographic terrain modelling, the application of triangular meshes i.e. TIN models is explained. The paper describes the collocation by least squares method that takes into account the impact of the distortions of neighbouring points in the observation point depending on the distance from the observation point. To model the final distortion computed by least squares collocation as smoothly as possible, the input values of distortion (positional and/or altitude residuals remaining after calculating transformation parameters) the concept of moving average is described, as well as the covariance function that is used to describe the influence of spatial distortion as function of distance. Using Briggs algorithm throughout GRID, the difference between two types of nodes is explained those that have perceived the distortion in their neighbourhood and those without close observation. Further, the deficiency of Briggs's method is processed that is primarily intended for highly polished so-called potential models. At the end, the formula for assessing the accuracy of the transformation GRID model is given.


DISTORTION MODELLING (INTRODUCTION)
Distortion -a model of alteration of the shape of transformed network -is the phenomenon caused by different amount and orientation of residuals -coordinate corrections between the initial and target coordinate system at various locations.It is locally to be expected that the amounts of distortion will be harmonized, i.e. that the residuals will mutually coincide in the direction and size of vectors.The most popular and intuitively the most explicit way of modelling the distortion, as well as of any potential physical surface, is the grid model, while topographic features of the terrain are for example more often modelled by means of TIN (Triangulated Irregular Network).In the first case, the requested mathematical process is the interpolation of discrete values of grid nodes, and in the second case, the connection of triangle networks according to the principle of Delaunay triangulation -which means that no point adjacent to the triangle should be inside of its circumcircle [1,2].

GRID-Lattice Model
Each observation area must be defined primarily by its boundaries in the plane (Y min , X min )−(Y max , X max ) or on ellipsoid (φ min , λ min )−( φ max , λ max ).The other important characteristic of GRID is its density, i.e. the number of columns and rows in the lattice, or the distance between its columns and rows.If the GRID is declared in the first mentioned way, then the distance between the columns and rows will be: and if the declaration refers to the given distances in both directions, then the number of rows and columns will be: ( The increase by '1' in the previous equation, i.e. the subtraction in the denominators (1), results from the fact that the number of 'wires' in the lattice is always larger by one than the number of rectangles between two neighbouring 'wires' along one of both axes.By analogy to this, the arbitrarily given point within the boundaries of the GRID area T(φ, λ) is positionally indexed (referenced) in the following way: where i is the number of the first, concretely more southern 'wire', i.e. row.It is analogously applied for the columns as well, i.e. j and the arbitrarily given point will be actually encircled by the rectangle of the GRID that has the value of the nodal point (i, j) on its southwestern part, (i + 1, j) northwest, (i + 1, j +1) northeast, and finally (i, j +1) southwest.
Such method of indexing is much more adequate for two-dimensional field, while the traditional method of indexing requires taking into consideration whether the discrete values of the GRID nodal point are written in rows and columns in external or internal loop.If the rows are in the external loop, i.e. if the columns are more frequently changed for each individual row, then the indexing is as follows: for one-dimensional array of the dimension n φ × n λ , where i and j are calculated in the same way as in (3).The interpolation of discrete values within an individual GRID cell is traditionally performed in 2 ways, i.e. by means of a simpler bilinear method [3] or by means of cubic spline method [4].Considering the level of complexity, the cubic spline interpolation has not shown better results than the bilinear transformation in densely enough approximated models [5].The unknown value in the observed point P is calculated from the known values in 4 closest neighbouring GRID points (Fig. 1).In the calculation of the values in the point P, the method of bilinear interpolation is used.The expressions ( 5) and ( 6) present the calculation of values δφ P in the point P: where: ) ( ) ) ( ) In the calculation of height distortion in the point P, the coefficients δφ or δλ in Eqs. ( 5) and ( 6) are replaced by the coefficients δh, hence, a new set of parameters c i is obtained for the height distortion.
Normally, any value of the quantity to be modelled is interpolated in the same way.

TIN-Triangular Model
About a hundred years ago, the first algorithms were developed for the modelling of the triangular network (TIN -engl.Triangular Irregular Network), so called Delaunay triangulation [2] or the Voroni polygons diagram [1].The relation between both models can be seen in Fig. 2. It is presumed that the terrain cannot have the points on the same x and y coordinates and different elevations.This presumption is necessary in order to reduce the problem of triangulation into the plane.The problem of triangulation is mathematically analogous to the problem of convex hull in the space, i.e. to the stretching of an elastic band through the given points.
The condition of triangulation is set: It is valid for each triangle ABC that all given points, except A, B and C are placed outside of the circumcircle of a triangle ABC (Fig. 3).This condition provides a unique solution of triangulation.
Figure 2 A point within the associated GRID cell [6] There are a few efficient algorithms for Delaunay triangulation, and one of the simplest and fastest is the socalled growing algorithm [7].The basic principle of the growing algorithm is adding the points into the existing network.The algorithm also requires that each point added to the network be placed within the convex polygon surrounding the existing network.Hence, the initial network is defined as a triangle with all given points included into it.The triangulation obtained in this way has 3 points more than the given points, and at the end, it is necessary to delete all triangles containing any of the 3 initial points.Figure 5 The point added to the network [6] Figure 6 The triangles to be deleted [6] Figure 7 Triangles have been deleted [6] The point is added (red in the Figure ) in the following way: First, it is checked which triangles need to be removed, and these are the ones with the circumcircle containing the point that needs to be added into the network (Fig. 6).After removing the triangle from the network, the point is placed within the polygon whose number of sides depends on the number and mutual position of removed triangles (Fig. 7).New triangles that are added to the network contain one side of the mentioned polygon each, and the point that needs to be added (Fig. 8).The procedure is repeated for all given points until obtaining the final network.
Figure 8 The point is added [6] For the purpose of comparing it with other models, TIN will finally be transformed into GRID again by means of the mentioned method Real Time Fineltra from chapter 2.3.[8], but also with certain modification considering not only the distortions in the vertices of the associated angle, but also of the neighbouring triangles.If we calculate the distortion as a weight mean at the midpoint of each triangle side from the distortion of 4 vertices defining two neighbouring triangles, then we actually have 6 distortion values in the associated triangle -3 in each vertex and 3 at the midpoints, so the final distortion of the point is within the triangle (7): (7) where the weights have been calculated as inverse values of the distance between the point and the vertices, i.e. midpoints of the sides raised to the selected power p i = 1/d n .Such approach reduced the effect of the "sharp" transition on the edges of the neighbouring triangles, i.e. sudden changes of distortion, as it is the case in Real Time Fineltra method.

Least-Squares Collocation
Least-squares collocation [9] is a technique that takes the influence of distortion of the neighbouring points into consideration in the observed point depending on the distance from the observed point.When calculating the distortion GRID, it is to use the accidentally arranged data in order to estimate the components of distortion (δφ, δλ) in each nodal point of the GRID from the known distortions of 7 points.The distance between each point and the nodal point of the GRID is a known quantity.

Moving Average
In order to make the final distortion model calculated by means of the least-squares collocation as smoothly as possible, the input values of distortion (positional and/or height residuals remaining after the calculation of transformation parameters) are reduced by the values of moving average.Moving average [10] is the method used to calculate general model trend at certain location (φ, λ) (in a certain radius s max , i.e.: where the weight is Finally, each individual quantity is additionally freed from the internal trend, i.e. the difference: because the final value of the distortion in each point by means of which the calculation is continued must be centred, i.e.: Detrending and centring of the distortion values by means of the moving average method is performed in order to avoid the singularity in solving the equations, i.e. of the cross-correlation matrix inversion, i.e. that Technical Gazette 25, 1(2018), 1-9 δ in order for the modelled quantity (distortion) to be freed from any systematic influence.

Empirical Covariance Function
Covariance function is used for describing the spatial influence of distortion as the function of distance (Fig. 9).The covariance function shows that the distortion of close points will be significantly correlated and that the distortion of the distant points will be non-correlated.The analytical model is then included into the empirical model.The analytical model is the tool for determining relevant characteristics of covariance between the input data and the observed point.Unlike the modelling of variogram being necessary for the application of Krigin method [11], the task of the covariogram modelling is to find optimal values of socalled correlation distance KU, i.e. of the quantity by means of which the exponent of the Gauss function will be attributed in the collocation itself.KU is the minimum value of the abscissa that reached at the ordinate value of the half of total covariance by the function interpolated on the basis of empirical values.It means that the abscissa presents the distances expressed in (m), (km) or some other units of distance, while the ordinate is the number line of covariance values -i.e. of the squared values of distances, i.e. (km 2 ), (m 2 ) or (cm 2 ).
The variance, being always a positive quantity, is the measure of probability distribution of some quantity, and practically it is the squared standard deviation of this quantity.In statistical terms, it is the expectation of random variable -e.g. the arithmetic mean value of the measured quantity.
On the other hand, two quantities, i.e. two random variables are observed in the same statistical model with the covariance, and it does not have to be positive, but can be positively and negatively correlated, i.e. dependent or generally non-correlated when two random variables are not dependent at all, which can be seen on covariogram when the graph asymptotically goes to zero.The input quantities for the modelling of covariogram are the range limit, the lat according to which the classes are formed.Maximum range is selected first of all intuitively, but also logically.It is the distance at which the measured data of one point should not have any real physical influence on the result of any other arbitrarily selected point, thus, the real physical influence can be negated by means of various mathematical approaches, or it can be reinforced also outside of its real range.
The variance of the whole set is calculated as: It is to be pointed out that C 0 is completely independent and exact value regardless of the values of input parameters in the modelling itself, i.e. of the maximum range and lat.The next step is the division of the area of maximum range into optimal number of classes or into optimal distance of each lat.In geostatistics, the size of the class lat of variogram should be determined.Similar is applied in covariograms.The decision is made on the basis of the following parameters: sample size (number of residuals, i.e. number of points) and the area of model capture.In [12], a concrete simple formula is given: where S is the area of the given territory, BT the number of points, and  the class lat in (m).For the area of about 3,026×10 11 (m 2 ) and for 1921 points, the optimal size of lat is 12,6 km.In all concrete examples given in the following chapter, the empirical variance function is modelled with the lat size of 10 km and maximum range of 100 km.
Then, the simplified covarinaces are calculated (because we deal here with correction/residuals, and not with mean errors of adjusted measurements) for all pairs of points (m, n), and the values are added into each class according to the obtained result of covariance.
The spherical distance d m,n is calculated further between the first pair of points (m, n), i.e. from the coordinates M(φ M , λ M ) and N(φ N , λ N ).According to the obtained result, it is then decided which class the pair belongs to, i.e. in which class the covariance of the pair shall be added that has not been calculated yet.The class index is determined as the integer value for each pair (m=1, BT−1;n=m, BT) where (m<n): where: d m,n is the spherical distance between the points M and N, and κ is the lat, i.e. the distance of each class.
The number of classes is: and the calculated covarinace is added to each class of index i that belongs to it according to the size of the distance between the points, d m,n : but the so called statistical frequency is also remembered, i.e. the number of pairs belonging to individual class on the basis of the classification according to mutual distance: because the final covariance of individual class of the indeks i is calculated afterwards as: and on the abscissa 2 ) 1 ( According to [12], the analytical function for the modelling of empirical covariance function is the following: where: d is the distance between the points, KU is the correlation distance, C 0 is the variance, the value of C(d) if d = 0, while [14] is simplified and modified: The correlation distance KU is determined numerically: where b κ and p κ are calculated as:  The adding of the numerator and denominator in Eq. ( 21) is performed only for those classes whose covariance value is placed within the interval between 30 % and 70 % of the variance of the whole set, and for those classes in which the number of the pairs of points (class frequency) is larger than the half of the first class frequency.It is interesting that a different selection of probability interval yields also completely different correlation distances for the same data set (Tab. 1).

Prediction by Least Squares
The prediction equation by least squares is [9]: The elements of the vector C l are calculated (Eq.( 25)) from the analytical covariance function depending on the distance between the known points and the observed GRID nodal point (d 1 , d 2 , d 3 , ..., d 7 in Fig. 10), the elements of the matrix C D are calculated (Eq.( 26)) from the analytical covariance function depending on the distances among all point combinations (d ij is the distance between the points i, j), and the Eq. ( 27) presents the vector l that contains centred distortions freed from the trend in all seven points.
The estimation of the expression presenting the prediction by least squares is a relatively simple process.The spatial distribution of data can cause the singularity of matrix covariance, which can result in the instability of the mathematical process.The consequence of the matrix covariance singularity can be the inability to predict the distortion in the GRID points or to predict the inaccurate values due to the rounding errors in the inverted matrix.
The other possibility is more dangerous due to the complex detecting of inaccurately predicted value.Since the matrix C D is always symmetrical, Cholesky algorithm is convenient for its inversion.
Finally, for the final value of distortion in any GRID nodal point, the trend of the moving average should be returned to the predicted δ ˆvalue, i.e.

Minimum Curvature Surface
The minimum curvatures surface of MCS was used for the positional transformation of coordinates between the North-American data NAD27 and NAD83 [15].The process of producing a surface can be considered as analogous to the process of mechanical bending of a thin metal elastic panel covering the observation area that tries to adapt to the points representing the values of distortion components at all locations under the influence of vertical forces and without deformations, and also without any influence of torsional, compression or extension mechanical forces accompanied by the smallest possible panel bending (minimum surface curvature).The practical implementation of this method using a computer program was done for the first time by [16].

Briggs Algorithm
In the whole GRID, we distinguish two types of nodes -those having the observed distortion in their environment and those without any close observations.For the GRID nodes observed in their vicinity (Fig. 11), the system of equations is set for the purpose of finding the value of surface curvature (29) that is minimized.First, the values of 5 coefficients b k are determined, and hence: , where b k are the coefficients that are calculated on the basis of 4 surrounding GRID nodes (the points A, B, C and D) and of one observation *E not focused on the GRID node.If the observation *E is not laced in the 1. quadrant looking from the node to be calculated, the points A, B, C and D are selected by rotating over 8 surrounding nodes.
Basically, the system of 5 equations and 5 unknowns is to be solved: where u i,j is the distortion in the node (i, j), u k is the distortion in the points A, B, C and D, and w n is the observation or given distortion in the node close to the point *E.
If there are more observations -given distortion placed near some node, it is suggested to approximate their values, where the differential coordinate differences (ξ 5 , η 5 ) should be approximated as well, apart from the distortion itself.
The nodes without the observations (given distortion) in their vicinity should be first classified according to their position in GRID.According to the position in GRID, the differential equations are the following [16].
The normal case when the node is further from the edge: For the node on the edge: ( ) One row or one column from the edge: In each of 4 GRID angles: Inside the angle: ( ) ( ) On the edge immediately to the angle: ( ) Out of the above mentioned 6 equations ( 33) -(38), the value u i,j is determined iteratively in each new iteration.For example, the value in a new iteration is determined from the Eq.(34) as: It is suggested for the initial values of the distortion in the GRID nodes to set some approximate values, e.g. the values of moving average.The iterative procedure is repeated, and the absolute value of maximum difference at any of the nodes between the two iteration circles does not go beyond the given value.
It should be pointed out related to the above mentioned equations that the indexing u i+1,j+1 actually presents the value of GRID nodal point that is moved from the nodal point presently calculated by one row and one column, i.e. if the nodal point in the 7 th row and the 3 rd column is calculated, u i+1,j+1 actually presents the value from the 8 th row and the 4 th column.

Splines with Tension
Two authors [17] have noticed the deficiencies of the Briggs method that is intended primarily for much smoothed so-called potential models.These deficiencies have been reflected first in the form of great harmonic oscillations in the areas of GRID with poorly or not at all present observations, i.e. first in the inflexion points.Therefore, the term of tension is introduced, and the original Briggs equations are modified.The GRID nodes with observation in their surroundings and those without observation are distinguished again.
For the GRID nodes having the observation in their vicinity, the coefficients b k are determined analogously to the original Briggs method, but the value of the GRID node is determined differently: The indexing z 11 is analogous to u i+1,j+1 from the original Briggs method, and z 0−2 would be written u i,j−2 according to Briggs.T I is the coefficient of internal tension, and  is the anisotropy of the GRID itself, i.e. the relation: The GRID nodes having no observation in the vicinity are determined in this case as follows (42): This equation can be graphically substantiated with the following presentation of weights in the surroundings of the observed node -black quadrant in Fig. 12.

Figure 12
The observed node and indexing of the neighbouring nodes [6] Instead of special conditions for boundaries and angles, the whole GRID is extended by 2 additional rows, i.e. columns on the left, up and down (Fig. 13).Technical Gazette 25, 1(2018), 1-9 Figure 13 The extension of GRID with additional columns and rows [6] The extension nodes are then calculated for the internal closer added row or column -in the concrete case, for the column of nodes added to the left (white triangles in Fig. 13): For the lower angle (crossed circle on the Fig. 13): For the external "most to the left" added column (white rectangles in Fig. 13 The boundary values are calculated analogously for the columns added to the right and for the rows added above and below the original GRID.In the Eqs.(43) to (45), T B is the coefficient of boundary tension.For relatively smoothed natural phenomena, e.g. for the gravity potential [17], the maximum amounts for internal and external tension are suggested: T I <0,25 and T B <0,1, while for the interpolation of topographic features, T I >0,35 and T B <0,2 are suggested.If T I =0 and T B =0 are selected for both tensions, the interpolated surface will be identical to the one according to Briggs method.
Since the Australian model [12] in which the method of multiple regression yielded significantly worse results than the minimum curvature surfaces or collocations, the research of the application of the method of multiple regression to be used in transformation was abandoned.

ACCURACY ESTIMATION BY MEANS OF GRID-MODEL TRANSFORMATION
The method of calculating the transformation accuracy is based on the following Ed.[12]: where: δ is the distortion component in the interdatum point interpolated by modelling from GRID, δ i is the known component of distortion in each point of input data, n is the number of input data, i.e. of interdatum points from which the GRID has been modelled.

CONCLUSION
It is to be pointed out in the Conclusion that the first comparison of stochastic and deterministic approaches, i.e. of GRID (collocations) and TIN (final elements) variants can be made only on the basis of equally arranged interdatum points on the entire state territory.
The accuracy of all three invariant surfaces of minimum curvature, the original Briggs method, as well as both spline methods with tension fall slightly behind the accuracies of the least-squares collocations (LSC solutions).
If the least-squares collocation should continue to be used for the modelling of distortion in the future, the following parameters should be optimized by all meansthe selection of optimal maximum distance for the correction of moving average, optimal maximum collocation distances, and research other analytical models of approximating the empirical covarinace function, because apart from various Gauss functions, there are very often logarithmic, exponential and spherical functions used today in geostatistics, as well as polynomials of various degrees.

Figure 4
Figure 4 Triangle with circumcircle [6] the distance between the observed location and each i th point within the radius s max , and i δ the real distortion of individual point.The rest of the point distortion is then calculated as:

Figure 9
Figure 9Principle of calculating the distortion components in the GRID nodal point[6]

Figure 10
Figure 10Covariance function depending on the distance[6] in the equations present; d κ the distance of the individual class (km), C κ the covariance of individual class, and C 0 the variance of the whole set.

Figure 11
Figure11 Calculation of a node when the observation *E is in the vicinity[6]