The Concept of Geodetic Analyses of the Measurement Results Obtained by Hydrostatic Leveling

: The article discusses the issue of hydrostatic leveling. Its application is presented in structural health monitoring systems in order to determine vertical displacements of controlled points. Moreover, the article includes a complete computation scheme that utilizes the estimation from observation di ﬀ erences, allowing the elimination of the inﬂuence of individual sensors’ systematic errors. The authors suggest two concepts of processing the measurement results depending on the sensors’ connection method. Additionally, the second concept is extended by the elements allowing the prediction of the displacements by means of Kalman ﬁltering.


Introduction
Hydrostatic leveling is mentioned as one of the oldest methods used to determine the height difference between controlled points. This method is mainly applied in measurements where highly accurate and highly reliable final results are necessary. However, it can be also used when geometric or trigonometric leveling is difficult to carry out. A relevant example is provided by the measurements made in 1939-hydrostatic leveling enabled establishing leveling connection between the Danish and Swedish height systems [1].
Hydrostatic leveling is often implemented in structural health monitoring (SHM) systems because of its high accuracy. The SHM systems with hydrostatic leveling were applied in objects such as dams [2], embankments [3], and metro tunnels; the latter systems are established during tunnel excavation [4]. Wei et al. [5] presented an application of hydrostatic leveling to determine vertical displacements of the SLS (Swiss Light Source) storage ring girders. Morishita and Ikegami [6] explained the slow-ground-motion monitoring system for the J-PARC (Japan Proton Accelerator Research Complex) linac, while Zhang et al. [7] focused on processing the power transmission tower vertical displacements. Other examples of hydrostatic leveling systems address the European Synchrotron Radiation Facility (ESRF) [8] and the TT1 tunnel at CERN (European Organization for Nuclear Research) [9]. The review above shows an enormous variety of objects involving hydrostatic leveling for displacement study.
The paper utilizes the adopted classification of hydrostatic leveling systems of Wilde et al. [10]. Thus, open and closed leveling systems are listed here. The open leveling systems are characterized by the simplest design in which the fluid freely flows between the sensors and is directly exposed to the environment. Contrary to open leveling systems, the closed leveling systems are characterized by absolute tightness of the sensors filled with working fluid, in order to reduce the effect of external factors on the hydrostatic leveling system performance.
In further analysis, the authors focus on introducing two concepts of processing the measurement results of the open leveling system. Concept I is based on the assumption that each sensor is connected directly to the expansion tank. Concept II, however, presupposes that sensors are mutually interconnected. Additionally, the second concept is supplemented by the Kalman filtering method that additionally enables the prediction of the identified displacements. The idea of hydrostatic leveling is based on the well-known Bernoulli equation in the following form [11]: where: ρ-fluid density, v-fluid velocity, p-hydrostatic pressure, g-standard acceleration of gravity, Z-fluid column height.
Assuming that for the majority of hydrostatic leveling systems the fluid is at rest (v = 0), Equation (1) takes the simplified form Besides the Bernoulli equation, this kind of leveling is also based on the principle of equal fluid-level in the communicating vessels. As a result, relation (2) can be used to determine the height difference between two sensors involved in the measurement of the height of fluid column Z. Therefore, the relation schematically illustrated in Figure 1 can be written as follows: where: subscript 1, 2-identification of first/second sensor, dZ 12 -height difference between the fulcrums of two sensors.
Geosciences 2019, 9, 406 2 of 12 In further analysis, the authors focus on introducing two concepts of processing the measurement results of the open leveling system. Concept I is based on the assumption that each sensor is connected directly to the expansion tank. Concept II, however, presupposes that sensors are mutually interconnected. Additionally, the second concept is supplemented by the Kalman filtering method that additionally enables the prediction of the identified displacements. The idea of hydrostatic leveling is based on the well-known Bernoulli equation in the following form [11]: where: ρ -fluid density, v -fluid velocity, p -hydrostatic pressure, g -standard acceleration of gravity, Z -fluid column height.
Assuming that for the majority of hydrostatic leveling systems the fluid is at rest ( v = 0 ), Equation (1) takes the simplified form Besides the Bernoulli equation, this kind of leveling is also based on the principle of equal fluid-level in the communicating vessels. As a result, relation (2) can be used to determine the height difference between two sensors involved in the measurement of the height of fluid column Z. Therefore, the relation schematically illustrated in Figure 1 can be written as follows: where: subscript , 1 2 -identification of first/second sensor, dZ 12 -height difference between the fulcrums of two sensors. After converting Equation (3), the result becomes Considering additional assumptions regarding equal hydrostatic pressure, fluid density, and standard gravity in both vessels ( , , p p g g 2 ), the height difference between two sensors may be derived from the relation This, however, is a simplified version of the hydrostatic leveling relation.
Observations collected through hydrostatic leveling may be affected by measurement errors, possibly caused by: Considering additional assumptions regarding equal hydrostatic pressure, fluid density, and standard gravity in both vessels (p 1 = p 2 , ρ 1 = ρ 2 , g 1 = g 2 ), the height difference between two sensors may be derived from the relation This, however, is a simplified version of the hydrostatic leveling relation.
Observations collected through hydrostatic leveling may be affected by measurement errors, possibly caused by:

•
Density variation of the fluid filling the sensor connecting pipe which is a function of temperature T and pressure p, ρ = f (T, p). It is difficult or even impossible to assess precisely the density variations; • Errors in determining the exact fluid column height in sensors; • The changes of atmospheric pressure; • Standard acceleration due to gravity changes; • A tidal phenomenon (low and high tides).
The errors mentioned above may trigger systematic errors s.

Materials and Methods
The aim of the paper is presenting the concept of geodetic processing of the hydrostatic leveling results. A number of relations are required to reach the final solution. It is assumed that Z (1) i is the fluid column height in sensor no. i (i = 1, . . . , n) in the initial epoch and Z (2) i is the fluid column height in the actual epoch.
Equation (5) for the initial epoch can be presented as follows: where: Z B -backward reading on sensor, Z F -forward reading on sensor, j-the number of height differences (j = 1, . . . , m).
Similarly, Equation (5) in the actual epoch takes the form where: Z B -backward reading on sensor, Z F -forward reading on sensor. When carrying out the hydrostatic leveling measurements, we can expect that the readings of the fluid column height may be affected by systematic errors. Consequently, the following can be noted: where s i -sensor no. i systematic error.
Thus: d Z A similar relation can be determined for readings in the actual epoch: where: Z In order to reject the impact of systematic errors on final results of the estimation, the method of observation differences may be applied as where: ds (2) = s (2) i − s (2) i+1 and ds (1) It is important that none of the analyzed papers describe the estimation based on observation differences to determine the sensor displacements.
Based on the assumption that systematic errors of individual sensors are constant over time (s Assuming that u B means the back sensor displacement and u F means the forward sensor displacement, the following formula can be applied: Converting relation (14), the equation of pseudo-observation corrections reads The equation above can be written in the matrix notation as follows: where In the present experimental course, the measurement results were obtained by means of the open leveling system. In both open and closed leveling systems, there were no redundant observations. Thus, the achieved results are unambiguous and there is no need to apply the optimization task written below: where: P = Q y −1 , Q y = P −1 (P-weight matrix, Q y -cofactor matrix of measurement results).
As no redundant observations appear, consequently the correction values equal zero (v = 0). Therefore, Equation (16) can be applied in the form of Au + y = 0, and so In order to assess the accuracy of computation results (in this case the displacement u i of each controlled point-sensor), the propagation law of cofactors can be used as follows: Assuming that D = A −1 and Q y = P −1 , the cofactor matrix Q u may be computed as The elements of the weight matrix P are determined further on. Applying the error propagation law for two sensors (notation: 1, 2) to relation (4), the following outcome is obtained: Geosciences 2019, 9, 406 5 of 12 or, in a simplified form, corresponding to relation (5): Based on relations (21) or (22), the mean square error values of height difference between two sensors may be obtained, both in the initial and actual epoch. The functional model dz j = dZ applied in the paper (described by Equation (13)) addresses the observation differences made in two epochs, therefore m 2 dz j can be specified as follows: where: -the mean square errors of height difference between two sensors computed from relation (21) or (22), in the initial and actual epochs, respectively. While the mean error values are known, the weight of individual pseudo-observations may be derived (p j = 1 m 2 dz j ), and the weight matrix is given as P = Diag(p 1 , . . . , p m ).
The performed analysis utilizes the Kalman filtering method. General assumptions regarding this method and its application in geodetic measurement result processing are discussed in [12][13][14][15][16]. In the paper, Kalman filtering is implemented to predict sensor displacements. In order to apply this methodology in the conducted analysis, it is essential to define a valid kinematic or dynamic model of the investigated object. Let us assume a general kinematic model including position, velocity, and acceleration in the form of where: u k -vector of displacements at epoch k, υ k -velocity vector at epoch k, a k -acceleration vector at epoch k, subscript k+1 -values denotes the prediction results, ∆t-time interval between epochs k and k+1 (∆t = t k+1 − t k ), and e u , e υ , e a -vectors of random disturbances, due to displacements, velocity, and acceleration at epoch k, respectively.
Assuming that X k+1 = [u k+1 , υ k+1 , a k+1 ] T is a prediction state vector, Φ k/k+1 is a transition matrix between states at epochs k and k+1, X k = [u k , υ k , a k ] T is a state vector at epoch k, and e = [e u , e υ , e a ] T is a vector of random disturbances, the relation (24) may be defined as follows (we assume that the expected value E(e) = 0): where the transition matrix takes the form In case of a computed prediction state vector X k+1 described by dependence (25), the covariance matrix C X k+1 is defined as follows: Kalman filtering enables simultaneous adjustment of the information obtained from prediction and observation on sensors. Thus, assuming that A k+1 is a known coefficient matrix at epoch k+1, v X k+1 is a vector of corrections for prediction at epoch k+1, and v k+1 is a vector of measurement residuals (called innovation) at epoch k+1, the system of observation equations can be presented in the form The weight matrix of innovation P v k+1 required for further consideration may be derived from the relation where: P y k+1 -weight matrix of observations at epoch k+1 and P −1 It allows the computation of the two dedicated gain matrices (gain matrix for observations K y k+1 and Kalman gain matrix K X k+1 ) according to the following formulas: Based on relation (29), there may be obtained the corrections for observations y k+1 made at epoch k+1: and consequently, the vector of adjusted observations in the following form: The corrections for prediction should be determined similarly: Then, the state vector at epoch k+1 is given by The covariance matrix corresponding to the state vector X k+1 described by the relation (34), can be written in the form where: I-unit matrix.

Results
The practical example introduces measurement results of the open levelling system without a gas reference system. An inquiry involving a five-sensor geodetic network is presented below. However, the computation scheme is similar for networks composed of a larger number of sensors. The sensors' location (for concept II) is shown schematically in Figure 2. In the presented leveling system, the expansion tank also acts as a reference sensor.

Results
The practical example introduces measurement results of the open levelling system without a gas reference system. An inquiry involving a five-sensor geodetic network is presented below. However, the computation scheme is similar for networks composed of a larger number of sensors. The sensors' location (for concept II) is shown schematically in Figure 2. In the presented leveling system, the expansion tank also acts as a reference sensor.  Table 1 includes the readings of the fluid column height obtained on each sensor in both measuring epochs. Further computational stages were carried out according to the two concepts of processing the measurement results, according to the sensors connection method: Concept I-each sensor is connected directly to the expansion tank, Concept II-sensors are connected in series and displacement prediction using Kalman filtering is implemented.
The displacement prediction by means of Kalman filtering may be implemented in both presented concepts. However, after performing analyses for the purposes of this paper it was decided to present the prediction on the example of concept II.

Concept I
Concept I is based on height differences between each sensor and the expansion tank (the reference sensor). The scheme of sensor location for concept I is shown in Figure 3. The height difference values in initial and actual epochs ( dZ (1) and dZ (2) ), as well as the difference of two height differences ( j dz ), are shown in Table 2. The reading on the reference sensor is adopted as the backward reading (according to relations (6) and (7)).  Table 1 includes the readings of the fluid column height obtained on each sensor in both measuring epochs. Further computational stages were carried out according to the two concepts of processing the measurement results, according to the sensors connection method: Concept I-each sensor is connected directly to the expansion tank, Concept II-sensors are connected in series and displacement prediction using Kalman filtering is implemented.
The displacement prediction by means of Kalman filtering may be implemented in both presented concepts. However, after performing analyses for the purposes of this paper it was decided to present the prediction on the example of concept II.

Concept I
Concept I is based on height differences between each sensor and the expansion tank (the reference sensor). The scheme of sensor location for concept I is shown in Figure 3. The height difference values in initial and actual epochs (dZ (1) and dZ (2) ), as well as the difference of two height differences (dz j ), are shown in Table 2. The reading on the reference sensor is adopted as the backward reading (according to relations (6) and (7)). Geosciences 2019, 9, 406 8 of 12

Concept II
Assuming that sensors are connected in series, as shown in Figure 2, concept II applies the height differences between individual sensors. Such a solution is often used in hydrostatic leveling systems, which is presented in [6,9,10]. If we consider the case with the controlled points equally spaced on the investigated object, then the implementation of either concept (I or II) is unrestricted. However, according to the paper's authors, in the case of extensive longitudinal networks (e.g., occurring in tunnels, large dams), duplicating pipes in order to connect each sensor directly to the expansion tank may be highly uneconomical and more difficult to perform.
The input data for concept II are: the height difference values in initial and actual epochs ( dZ (1) and dZ (2) ) as well as the difference of these height differences ( j dz ); they are presented in Table 3.
In order to assess the accuracy of the adjustment results by means of the cofactor matrix (relation (20)), it is necessary to determine the mean errors for each dz j . It is assumed that the reading mean error is identical for all sensors in both measuring epochs (m Z B = m Z F = 0.1 mm). As a consequence, the mean errors of displacements for each sensor m u i equals 0.2 mm. If |u 1 | ≥ 2m u i , the displacement is considered significant, i.e., all sensor displacements are significant.

Concept II
Assuming that sensors are connected in series, as shown in Figure 2, concept II applies the height differences between individual sensors. Such a solution is often used in hydrostatic leveling systems, which is presented in [6,9,10]. If we consider the case with the controlled points equally spaced on the investigated object, then the implementation of either concept (I or II) is unrestricted. However, according to the paper's authors, in the case of extensive longitudinal networks (e.g., occurring in tunnels, large dams), duplicating pipes in order to connect each sensor directly to the expansion tank may be highly uneconomical and more difficult to perform. The input data for concept II are: the height difference values in initial and actual epochs (dZ (1) and dZ (2) ) as well as the difference of these height differences (dz j ); they are presented in Table 3. Table 3. The input data for concept II.

Height Difference dZ j Initial Epoch [mm] Actual Epoch [mm] dz j [mm]
j = re f . It is assumed again that the reading mean error for each sensor equals 0.1 mm. Thus, the following mean errors of displacements were collected: m u 1 = 0.2 mm, m u 2 = 0.3 mm, m u 3 = 0.3 mm, m u 4 = 0.4 mm, and m u 5 = 0.4 mm. This case considers all sensor displacements significant.
Both presented solutions allow the derivation of identical results regarding sensor displacements. However, the mean errors of displacements due to concept II were generally 0.1-0.2 mm higher than in concept I.
The next computational stage was the displacement prediction using Kalman filtering. The displacement prediction is one of the aspects that the authors recognize to be crucial in the field of structural health monitoring systems. It is related to the fact that if the displacements approach their maximum values provided for the investigated object, it is possible to predict if the ultimate limit state will be exceeded soon and when exactly this may happen. In the paper, this method in particular is implemented because Kalman filtering enables simultaneous adjustment of the information obtained from prediction and observation on sensors. This allows, in the considered epoch, the determination of corrections for observations, as well as corrections for the prediction state vector. After including the obtained corrections, the updated kinematic model parameters make it possible to more effectively predict the displacements in the next epoch. To conclude, adopting simultaneous adjustment each time enables the improvement of the definition of kinematic model parameters, and therefore grants a more accurate displacement prediction. The significant argument for applying this approach is also the fact that during displacement prediction there may be measurement results used that have been obtained from other sensors.
A general model of the investigated object (specified by Equation (24)) in the conducted analysis is reduced by acceleration a. This is the effect of prior computations, hence the adoption of the model excluding acceleration does not cause significant result distortion. As a consequence, the kinematic model (24) can be written in a simplified form: Based on available data, the identical velocity corresponding to sensor nos. 2, 3, 4, and 5 is assumed, which equals υ k = −0.35 mm/h, while for sensor no. 1 it equals υ1 k = 0.4 mm/h. In this particular case, there is an issue of uneven sensor vertical displacements, thus the assumption of identical velocity for all sensors would lead to errors in final prediction.
For computation purposes, a two-hour prediction is carried out. The estimated state vector at epoch k is given as making it possible to determine the prediction state vector at epoch k+1 as follows: The next stage of the performed analysis is computing the prediction state vector covariance matrix at epoch k+1, which enables specifying the prediction accuracy. According to relation (26), the state vector covariance matrix at epoch k-C X k is required. The C X k matrix part referred to sensor displacements involving the cofactor matrix of displacements (20). Additionally, the mean errors of velocities are assumed as m υ1 k = 0.2 mm/h and m υ k = 0.35 mm/h. No more detailed information is available on the nature of the kinematic model random disturbances, hence the following is assumed C e = C X k . Based on the above, the mean error values of parameters obtained from prediction are assessed as follows:m u 1 k+1 = 0.5 mm, m u 2 k+1 = 0.8 mm, m u 3 k+1 = 0.9 mm, m u 4 k+1 = 0.9 mm, m u 5 k+1 = 0.9 mm, m υ1 k+1 = 0.3 mm/h, and m υ k+1 = 0.5 mm/h.
It is also noteworthy that nowadays, when utilizing hydrostatic leveling systems, there is continuous access to the current data from electronic sensors. Contrary to measurements made at selected moments of time, e.g., once a month or once a year, when the implementation of the total concept of Kalman filtration is justified, in the case of almost continuous measurement results updating, the further part of the analysis may seem to be of little use. However, for the purposes of this paper, the authors decided to present the complete results obtained by means of Kalman filtering.
In order to continue the computations, the readings of the fluid column height on each sensor at epoch k+1 were carried out. The following readings were obtained: Z no.5 = 85.5 mm. The height difference values in initial and k+1 epochs (dZ (1) and dZ (k+1) ), as well as the difference of these height differences (dz j ), were adopted as the input data for further analysis. They are presented in Table 4. As a consequence of the conducted analysis, corrections for observation and corrections for prediction at epoch k+1 were determined. The obtained values are presented in tabular form, as shown in Table 5. The corrections gained for prediction enable the determination of the state vector X k+1 at epoch k+1 in the form X k+1 = u 1 k+1 u 2 k+1 u 3 k+1 u 4 k+1 u 5 k+1 | υ1 k+1 υ k+1

Discussion
The article discusses the methods of processing the results of a hydrostatic leveling system application. A complete computation scheme made it possible to determine the sensor displacements employing estimation from observation differences. This approach has not been suggested in any of the available literature items concerning hydrostatic leveling. This approach allows the elimination of systematic errors possible in geodetic monitoring systems. In case of its rejection, it would be necessary to consider potential systematic errors.
Two concepts of measurement result processing were analyzed; both of them led to the same values of sensor displacements. The major difference addressed the mean errors of displacements, however, the lower mean error values were obtained in case of a direct connection between each sensor and the expansion tank (concept I). Concept II proposed the implementation of Kalman filtering in order to predict displacements, which the authors consider the crucial aspect related to structural health monitoring systems. Based on the correct displacement prediction, the decision-making process concerning an object's safety may be much improved. As previously mentioned, Kalman filtering may be also successfully applied in computation according to concept I, which the authors confirmed at the stage of preliminary analysis. Since the state vector in epoch k for concept I is similar to the one for concept II (the same displacement values were obtained for each sensor), the prediction state vector in epoch k+1 for concept I is given as X k+1 = u 1 k+1 u 2 k+1 u 3 k+1 u 4 k+1 u 5 k+1 | υ1 k+1 υ k+1 The achieved values were consistent with concept II sensor displacements. Since each sensor is connected directly to the expansion tank, in computations the mean errors of velocities are assumed as. m υ1 k = m υ k = 0.2 mm/h. By analogy with previous considerations, for concept I there were computed lower mean error values of parameters obtained from prediction (m u i for each sensor equals 0.5 mm, m υ1 k+1 = 0.3 mm/h and m υ k+1 = 0.3 mm/h).