Adjustment of an Integrated Geodetic Network Composed of GNSS Vectors and Classical Terrestrial Linear Pseudo-Observations

: The paper proposes a new method for adjusting classical terrestrial observations (total station) together with satellite (GNSS-Global Navigation Satellite Systems) vectors. All the observations are adjusted in a single common three-dimensional system of reference. The proposed method does not require the observations to be projected onto an ellipsoid or converted between reference systems. The adjustment process follows the transformation of a classical geodetic network (distances and horizontal and vertical angles) into a spatial linear (distance) network. This step facilitates easy integration with GNSS vectors when results are numerically processed. The paper offers detailed formulas for calculating pseudo-observations (spatial distances) from input terrestrial observations (horizontal and vertical angles, horizontal distances, height of instrument and height of target). The next stage was to set observation equations and transform them into a linear form (functional adjustment model of geodetic observations). A method was provided as well for determining the mean errors of the pseudo-observations, necessary to assess the accuracy of the values following the adjustment (point coordinates). The proposed algorithm was veriﬁed in practice whereby an integrated network made up of a GNSS vector network and a classical linear-angular network was adjusted.


Introduction
Integrated measurement methods are usually employed for various surveying engineering jobs, such as monitoring land surface displacements or structure deformation [1][2][3][4][5][6][7]. Classical (terrestrial) surveying techniques are usually based on control networks referred to as a local (national) system of coordinates [8,9]. Survey results are usually processed by a simultaneous adjustment of classical observations (angles and distances) and GNSS vectors in a common mathematical space [5]. Integrated networks may be adjusted on the GRS'80 (Geodetic Reference System '80) reference ellipsoid surface or a horizontal projection plane in a local system. It is necessary to pre-process the observations in both cases. This process can include the projection of GNSS vectors onto an ellipsoid (calculating the length of a geodetic line and its original azimuth), projection of classical observations (horizontal distances) onto the surface of an ellipsoid (calculating the projection corrections), or transformation of GNSS vectors (∆X, ∆Y, ∆Z) onto a horizontal plane [10][11][12][13]. The determination of the height of the GNSS network points (e.g., calculation of ellipsoidal heights and their conversion into values referenced to the local model of geoid) is a separate computational stage [14][15][16][17]. All the pre-processing is rather labour-intensive and requires practical experience and knowledge [13,18,19]. What is more, one cannot avoid errors resulting from the transforming of the original observations into pseudo-observations on a common plane of reference (such as projection errors or errors of the geoid model).
Integrated networks are proposed to be used if satellite signal exposure is insufficient (in forests or difficult topography). It is then that classical observations can be used to improve the GNSS vector network. Detailed investigations into integrated geodetic networks 2 of 13 can be found in many available publications [6,14,[20][21][22]. For example, Kutoglu [8] showed a method for adjusting a GNSS network as a linear (trilateration) network. He proved that slope distances measured with any surveying method can be adjusted in any reference system (cf. [23]). Gargula [24] proposed an alternative adjustment method whereby both distances and angles between GNSS vectors are calculated in a geocentric spatial system (XYZ). A resulting set of linear-angular non-reduced pseudo-observations can be adjusted in reference to a local reference system.
Land-surveying measurements are adjusted because they are burdened with random errors (cf. [8,25]). The errors are treated like normally distributed random variables (Gauss distribution, Figure 1). observations on a common plane of reference (such as projection errors or errors of the geoid model).
Integrated networks are proposed to be used if satellite signal exposure is insufficient (in forests or difficult topography). It is then that classical observations can be used to improve the GNSS vector network. Detailed investigations into integrated geodetic networks can be found in many available publications [6,14,[20][21][22]. For example, Kutoglu [8] showed a method for adjusting a GNSS network as a linear (trilateration) network. He proved that slope distances measured with any surveying method can be adjusted in any reference system (cf. [23]). Gargula [24] proposed an alternative adjustment method whereby both distances and angles between GNSS vectors are calculated in a geocentric spatial system (XYZ). A resulting set of linear-angular non-reduced pseudo-observations can be adjusted in reference to a local reference system.
Land-surveying measurements are adjusted because they are burdened with random errors (cf. [8,25]). The errors are treated like normally distributed random variables (Gauss distribution, Figure 1). The error of measurement (ε) exhibits normal distribution if the error's density function f(ε) can be expressed as: where: e-the base of a natural logarithm. The probability that the measurement error ε (as a random variable) will take a value from the interval (ε1; ε2) can be expressed in the form ( Figure 1): where: dε-the differential of ε.
As the actual values of errors of measurement ε are unknown, they are replaced with observation corrections v, to be determined during adjustment. The goal of adjusting observations with the least-squares method is to select such correction values v that the sum of their squares multiplied by the weights (p) is the smallest.
The principal idea behind the adjustment method proposed in this paper involves the transformation of a classical geodetic network (distances and horizontal and vertical angles) into a spatial linear (distance) network independent of the local reference system. The objective is to adjust a classical network and a GNSS vector network together in a common, geocentric XYZ coordinate system (referenced to the GRS'80 ellipsoid). The error of measurement (ε) exhibits normal distribution if the error's density function f (ε) can be expressed as: where: e-the base of a natural logarithm. The probability that the measurement error ε (as a random variable) will take a value from the interval (ε 1 ; ε 2 ) can be expressed in the form ( Figure 1): where: dε-the differential of ε.
As the actual values of errors of measurement ε are unknown, they are replaced with observation corrections v, to be determined during adjustment. The goal of adjusting observations with the least-squares method is to select such correction values v that the sum of their squares multiplied by the weights (p) is the smallest.
The principal idea behind the adjustment method proposed in this paper involves the transformation of a classical geodetic network (distances and horizontal and vertical angles) into a spatial linear (distance) network independent of the local reference system. The objective is to adjust a classical network and a GNSS vector network together in a common, geocentric XYZ coordinate system (referenced to the GRS'80 ellipsoid).

Creating Linear Pseudo-Observations from the Classical Terrestrial Measurements
Mathematical equations needed for the task are developed using the principle of indirect levelling ( Figure 2) used in the traditional topographic survey (total station). This way one can calculate (using the Pythagorean theorem) the spatial (actual) distance between two points of a geodetic network (j-instrument station, k-target position): Appl. Sci. 2021, 11, x FOR PEER REVIEW 3 of 14

Creating Linear Pseudo-Observations from the Classical Terrestrial Measurements
Mathematical equations needed for the task are developed using the principle of indirect levelling ( Figure 2) used in the traditional topographic survey (total station). This way one can calculate (using the Pythagorean theorem) the spatial (actual) distance between two points of a geodetic network (j-instrument station, k-target position): ( ) Next, the formulation to calculate the height difference h (see Figure 2) is substituted into Equation (3).
This yields a general equation for the spatial distance between point j (station) and the measured point k as a function of the initial observations (d̅ , α, i, s): Apart from this type of distance (5), the adjustment of the spatial network will require the length of the section between two points (L-left target and R-right target), measured from station S ( Figure 3).  Next, the formulation to calculate the height difference h (see Figure 2) is substituted into Equation (3).
This yields a general equation for the spatial distance between point j (station) and the measured point k as a function of the initial observations (d, α, i, s): Apart from this type of distance (5), the adjustment of the spatial network will require the length of the section between two points (L-left target and R-right target), measured from station S ( Figure 3).

Creating Linear Pseudo-Observations from the Classical Terrestrial Measurements
Mathematical equations needed for the task are developed using the principle of indirect levelling ( Figure 2) used in the traditional topographic survey (total station). This way one can calculate (using the Pythagorean theorem) the spatial (actual) distance between two points of a geodetic network (j-instrument station, k-target position): ( ) Next, the formulation to calculate the height difference h (see Figure 2) is substituted into Equation (3).
This yields a general equation for the spatial distance between point j (station) and the measured point k as a function of the initial observations (d̅ , α, i, s): Apart from this type of distance (5), the adjustment of the spatial network will require the length of the section between two points (L-left target and R-right target), measured from station S ( Figure 3).  First, the horizontal distance d LR is determined with the law of cosines: The relationship between the horizontal distance d LR and the spatial distance d LR is shown in Figure 4 and Equation (7).
Appl. Sci. 2021, 11, x FOR PEER REVIEW 4 of 14 First, the horizontal distance d̅ LR is determined with the law of cosines: The relationship between the horizontal distance d̅ LR and the spatial distance dLR is shown in Figure 4 and Equation (7).
Substitution of Equation (6) to Equation (7) and employment of the general Equation for height difference (4) yields an equation for distance dLR as a function of results of a classical field survey (horizontal distances d̅ SL, d̅ SR; horizontal angle βLSR; vertical angles αSL, αSR and signal heights sL, sR): The distances calculated with (5) and (8) will be considered linear pseudo-observations.

Stochastic Adjustment Model
Proper numerical processing of geodetic survey data involves the adjustment of observations (according to the method of least squares) and assessment of the accuracy of the results. To this end, it is necessary to transform the mean errors of the original observations (d̅ , α, β, I, s-see Figures 2 and 3) into mean errors of the pseudo-observations. To do this, one can employ the propagation of mean error [7,25].
The mean error of the pseudo-observations djk (5), which are distances between the instrument station j and signal k (target), is expressed as: Partial derivatives (∂) of each variable (survey result) are determined as follows: Figure 4. Determination of the height difference between two ground points L and R.
Substitution of Equation (6) to Equation (7) and employment of the general Equation for height difference (4) yields an equation for distance d LR as a function of results of a classical field survey (horizontal distances d SL , d SR ; horizontal angle β LSR ; vertical angles α SL , α SR and signal heights s L , s R ): The distances calculated with (5) and (8) will be considered linear pseudo-observations.

Stochastic Adjustment Model
Proper numerical processing of geodetic survey data involves the adjustment of observations (according to the method of least squares) and assessment of the accuracy of the results. To this end, it is necessary to transform the mean errors of the original observations (d, α, β, I, s-see Figures 2 and 3) into mean errors of the pseudo-observations. To do this, one can employ the propagation of mean error [7,25].
The mean error of the pseudo-observations d jk (5), which are distances between the instrument station j and signal k (target), is expressed as: Partial derivatives (∂) of each variable (survey result) are determined as follows: The mean error of the pseudo-observations d LR (4), which are distances between the left (L) and right (R) target (Figure 3), is determined as shown below (also using the propagation of mean error): Partial derivatives of each variable (survey results) are expressed with the following equations: The formulas for partial derivatives appearing in Equations (9) and (10) were verified also for units.
The law of propagation of variance and covariance can be used to determine the mean errors instead of the propagation of mean error (e.g., [8,25]) because consecutive pseudo-observations can depend on the same angles and distances. Nevertheless, based on previous tests (cf. [24]), the actual effect of correlation of angles and distances on the values of calculated pseudo-observations is negligible.
The stochastic model for the integrated network is complemented with information on a priori mean errors (m (∆x) ; m (∆y) ; m (∆z) ) of components of the GNSS vector (∆x jk ; ∆y jk ; ∆z jk ), which can be obtained in post-processing [26].

Functional Adjustment Model
The creation of the functional model of the adjustment of a geodetic network involves the listing of observation equations and transforming them into linear equations of correction. The formulas below are general equations for three types of observations (in an integrated spatial geodetic network): (1) the station-target distance; (2) the left target-right target distance; (3) components of the GNSS vector.

The Procedure for Adjusting the Integrated Network
The adjustment of the integrated geodetic network (using the method of least squares) will be based on an overdetermined system of equations of correction type (12), (15) and (18), expressed as the following matrix form: where: T -the vector of corrections type (12), (15) and (18) to be determined (curly brackets { . . . } stand for all elements of a type); A-the matrix of coefficients of the unknowns (partial derivatives) in Equations (12), (15) and (18); X = (δx L ; δy L ; δz L ; δx R ; δy R ; δz R ); δx j ; δy j ; δz j ; δx k ; δy k ; δz k T -the vector of the unknowns-increments to approximate coordinates; T -the vector of absolute terms type (13), (16) and (19).
The estimated vector of the unknownsX is calculated with the method known from the adjustment calculus [13], which stems from the imposition of the least square condition (V T ·P·V = min.) on the system (21): where: P-the matrix of weights set up from mean errors of the linear pseudo-observations (9), (10) and mean errors of GNSS vector measurements (m (∆x) ; m (∆y) ; m (∆z) ): The next step is to substitute the vector of unknowns X (21) with the calculated vector X (22) and calculate the vector of observation corrections V, which are used to adjust the observations-the left sides of the observation Equations (11), (14) and (17).
Information on mean errors of the adjusted coordinates (m x , m y , m z ) can be found on the diagonal of the covariance matrix (Q x ) of the vector X: where: m 0 -the standard error of unit weight; r-the number of redundant observations.
Next, a single parameter characterising the point's accuracy is calculated for each point-the error of position (m P ) in a three-dimensional Cartesian system:

Results and Discussion (Numerical Example)
The proposed method for adjusting an integrated network was verified with a simple practical example (Figure 5b). Calculations were performed for a GNSS vector network (without classical linear-angular observations, Figure 5a) as well, to compare the results. The test calculations involved actual survey results for a section of a control network for monitoring ground displacement in an active mining area.  (Table 1) were obtained from post-processing of GNSS data. Approximate coordinates of the points to be determined ( Table 2) were calculated from measured (non-adjusted) vectors. For the purpose of the comparative analyses, coordinates and distances were recorded down to 0.0001 m.  Vectors between vertices of the network were measured with a static GNSS method (two receivers, survey duration 60 min). GNSS vectors were used in both test variants: GNSS network (Figure 5a) and integrated network (Figure 5b).
Classical terrestrial measurements (horizontal and vertical angles, horizontal distances) were performed with a precise total station. The surveying instrument and the signal (reflector) were placed on carefully levelled and centred tripods with tribrachs. The height of both the instrument (i) and the signal (s) was measured repeatedly with a special device. The classical observations were used in the integrated network (Figure 5b).
For calculation purposes, points 2 and 6 were assumed to be fixed (reference points), while the other points (3, 4 and 5) were to be determined-in both cases (Figure 5a,b). Tables 1 and 2 show the input data necessary to adjust the vector network (Figure 5a). Values of components of the GNSS vectors (∆X, ∆Y and ∆Z) and their mean errors (Table 1) were obtained from post-processing of GNSS data. Approximate coordinates of the points to be determined (Table 2) were calculated from measured (non-adjusted) vectors. For the purpose of the comparative analyses, coordinates and distances were recorded down to 0.0001 m. Linear pseudo-observations, that is, actual spatial distances, and their mean errors were calculated from the classical observations (Tables 3 and 4). Note the values of errors m jk (d) . They are identical to the mean error of the measured horizontal distance (Table 3) due to small differences between horizontal distances d ( Table 3) and spatial distances d jk (Table 4).  Table 5 shows part of a table with a matrix of coefficients A (20), built as the integrated network is being adjusted. The remaining part of matrix A is filled with coefficients from equations of classical GNSS vectors (18). The absolute terms (the last table column) in the equations of corrections (12) and (15) are necessary to create the matrix L (20). Table 5. Coefficients at the unknowns (partial derivatives) and the absolute terms (13), (16) in the equations of corrections for the pseudo-observations (12), (15).  Tables 6 and 7 (for the two test variants, respectively). A comparison of results of the adjustment of the GNSS vector network and the integrated network (Table 8) reveals only minor differences in the coordinates X, Y, Z (up to about 5 mm). The difference in the point location is illustrated by the resultant linear discrepancy (δ XYZ ), which assumes values from about 3 mm to about 5 mm. The values demonstrate the impact of the classical measurements on the adjustment of the GNSS vector network (this order of differences in coordinates can be significant when the absolute ground displacements are measured, for example). Note also the mean errors of the coordinates (m X , m Y , m Z ) and the error of position (m P ). They are almost identical for both variants. This might be indicative of the correctness of equations of the pseudoobservation mean errors (9) and (10), always assuming that both the procedures employed similar observation weighting principles. For ground displacement monitoring, the relative positions of points (so-called relative displacement) are measured as well. Therefore, the adjusted coordinates (Tables 6 and 7) were used to calculate spatial distances between the points (Table 9). Despite differences of several millimetres in point coordinates (Table 8), the computed values are practically identical for both variants (the differences are below 1 mm in most cases). This fact may indicate cumulative GNSS survey errors (caused by antenna phase centre variations, for example). The calculated distances (Table 9) were juxtaposed with reference distances (d CL ), which were obtained from the additional precise (repeated several times) classical surveys: the horizontal distances were measured with a precise total station, and the height differences were measured using the precise geometric levelling method. Allowing for the hypothetical assumption that spatial distances d CL are free of errors, one can note a positive impact of the additional classical measurements (Difference 3) on most distances measured with GNSS (Difference 2).

Summary and Conclusions
The paper proposes a new method for adjusting an integrated network made up of GNSS vectors and classical terrestrial observations. The first computing stage involves a list of linear pseudo-observations that are original linear-angular observations converted into spatial distances. The next step is to transform (a priori) the mean errors of the classical measurements into mean errors of the pseudo-observations (the stochastic model). The functional model of the adjustment is made up of a set of observation equations for the GNSS vectors and for the new pseudo-observations (expressed as a function of the original linear-angular measurements). The objective of the pre-processing is the concurrent adjustment of the pseudo-observations and GNSS vectors in a common mathematical space XYZ.
The second part of the paper presented a practical application of the method to adjust an integrated geodetic network. The results (coordinates of points, their mean errors and spatial distances between the points) were juxtaposed with results of the adjustment of the GNSS vector network.
The comparative analysis demonstrated that the new method for adjusting integrated networks yields similar results (coordinates) as an adjustment of a vector network. It has further been demonstrated that an integrated network provides more accurate (close to real) spatial distances between points (in relation to reference distances obtained from precise classical surveys).
The primary advantage of the proposed adjustment method is the ease of integration of GNSS vectors with linear pseudo-observations obtained from classical measurements. The preparation of pseudo-observations is much easier than for other available methods for adjusting integrated networks (where it is necessary to determine the lengths of geodetic lines and their azimuths on the ellipsoid, project classical observations onto the ellipsoid, or convert ellipsoidal heights into orthometric values, etc.). The proposed method can be employed to adjust periodic measurements of control networks for ground displacement monitoring. However, the use of this calculation method is limited to short lines. If distances between the points are longer than 200-300 m, the effect of refraction and the curvature of the earth should be considered [27].
The research reported here will be continued. A detailed computing algorithm for the method and its implementation as a computer application is planned. Furthermore, an attempt will be made to test the new adjustment method on a network with much longer GNSS vectors and pseudo-observations (about a few hundred metres).