**Automatic Calibration of Process Noise Matrix and Measurement Noise Covariance for Multi-GNSS Precise Point Positioning**

**Xinggang Zhang 1,2,3, Pan Li 3,\*, Rui Tu 1,2, Xiaochun Lu 1,2, Maorong Ge 3,4 and Harald Schuh 3,4**


Received: 28 January 2020; Accepted: 12 March 2020; Published: 2 April 2020

**Abstract:** The Expectation-Maximization algorithm is adapted to the extended Kalman filter to multiple GNSS Precise Point Positioning (PPP), named EM-PPP. EM-PPP considers better the compatibility of multiple GNSS data processing and characteristics of receiver motion, targeting to calibrate the process noise matrix *Qt* and observation matrix *Rt*, having influence on PPP convergence time and precision, with other parameters. It is possibly a feasible way to estimate a large number of parameters to a certain extent for its simplicity and easy implementation. We also compare EM-algorithm with other methods like least-squares (co)variance component estimation (LS-VCE), maximum likelihood estimation (MLE), showing that EM-algorithm from restricted maximum likelihood (REML) will be identical to LS-VCE if certain weight matrix is chosen for LS-VCE. To assess the performance of the approach, daily observations from a network of 14 globally distributed International GNSS Service (IGS) multi-GNSS stations were processed using ionosphere-free combinations. The stations were assumed to be in kinematic motion with initial random walk noise of 1 mm every 30 s. The initial standard deviations for ionosphere-free code and carrier phase measurements are set to 3 m and 0.03 m, respectively, independent of the satellite elevation angle. It is shown that the calibrated *Rt* agrees well with observation residuals, reflecting effects of the accuracy of different satellite precise product and receiver-satellite geometry variations, and effectively resisting outliers. The calibrated *Qt* converges to its true value after about 50 iterations in our case. A kinematic test was also performed to derive 1 Hz GPS displacements, showing the RMSs and STDs w.r.t. real-time kinematic (RTK) are improved and the proper *Qt* is found out at the same time. According to our analysis despite the criticism that EM-PPP is very time-consuming because a large number of parameters are calculated and the first-order convergence of EM-algorithm, it is a numerically stable and simple approach to consider the temporal nature of state-space model of PPP, in particular when *Qt* and *Rt* are not known well, its performance without fixing ambiguities can even parallel to traditional PPP-RTK.

**Keywords:** EM-algorithm; multi-GNSS; PPP; process noise; observation covariance matrix; extended Kalman filter; machine learning

#### **1. Introduction**

Since Precise Point Positioning (PPP) emerged [1,2], people are primarily focusing on improving precise orbit and clock products, developing new algorithms to solve for ambiguities, to accelerate its

convergence and expand its applications such as PPP-real-time-kinematic (PPP-RTK), triple frequency PPP [3], ionosphere-constraint PPP and low-cost receiver PPP [4–10].

Generally, PPP can be realized by the least-squares method (including sequential least-squares) or extended Kalman filter. The least-squares method is for static state estimation and thus does not reflect varying user dynamics. To work the same as Kalman filter, the process noise matrix is added to the gain matrix of the sequential least-squares method to adjust receiver clock behavior and atmospheric activity and so on, which is named as a sequential filter [2]. Hence, in the following paper, the authors will only consider the Kalman filter for PPP data processing.

Although both the process noise *Qt* and observation covariance matrix *Rt* are the key to Kalman filter, limited attention is paid to the fundamental problem for multi-GNSS PPP. *Qt* and *Rt* must be consistent with state dynamics and measurement accuracy, respectively. For example, if the value of *Qt* is too small, the estimated state will lose its minimum mean squared error property, and if the value of *Qt* is too large with respect to the correct one, the estimated state will oscillate around the true value. Moreover, because of ground deformation and specific surroundings, *Qt* should not be kept fixed to calculate the optimal estimates. In other words, *Qt* should evolve with time and a proper *Qt* will shorten the PPP convergence time. If *Qt* is improper, it may damage PPP convergence sometimes.

As for GNSS, *Rt* not only depends on the measurement accuracy, elevation of the satellite, orbit and clock error, atmospheric delay error, multipath, missing data, etc. but possibly deteriorate after a lapse of a period. What is more, the assumption, frequently used in geodesy, that different types of measurements have a fixed error ratio is not always true, because the ratio is closely linked to receivers and antenna types, and to the performance of satellite system itself. For example, while fusing multiple GNSS, the weight of measurements of GPS is intended to be higher, other GNSS systems to be lower. It is not easy to give a prior accurate ratio.

Generally, four methods are often applicable to calibrate *Qt* and *Rt*. The first one is based on the innovation property of the Kalman filter, in which a moving-window recursive way is used to identify *Qt* and *Rt* [11–14]. However, none of them can maintain the positive semi-definiteness of the estimated covariances. To solve this problem, Odelson developed the autocovariance least-squares method for estimating covariances using a lagged autocovariance function [15]. This kind of least-squares method depends on the user-defined autocovariance function.

The second scheme to recognize *Qt* and *Rt* is the multiple model adaptive estimation (MMAE) [13,16]. MMAE runs a bank of Kalman filter in parallel, every one of them is driven by its pair of *Qt* and *Rt*. The final *Qt* and *Rt* are thought of as the weighted sum of the estimates of individual Kalman filter.

In the third scheme, M-estimator is introduced into an adaptive Kalman filter to increase its resistance to outliers, where an adaptive factor α to state error covariance matrix is constructed [17,18]. Yet, choosing a value for α is still very challenging. An improper α will result in biased results.

Another attractive scheme is the least-squares variance component estimation (LS-VCE) [19], which is based on least-squares principles. Similar to restricted maximum likelihood (REML), LS-VCE does not use observations directly but combine observations to exclude any fixed effects. However, LS-VCE needs to define the weight matrix on the user's own and increase its complexity.

In this contribution, a machine learning algorithm, the Expectation-Maximization (EM) algorithm, is developed to the extended Kalman filter to estimate PPP states, <sup>→</sup> *x t*, together with a large number of *Qt* and *Rt*. The EM-algorithm, which can be classified as the first scheme, works in an iterative procedure to locate maximum likelihood estimates of parameters. Its iteration consists of two steps: Expectation and Maximization. In the Expectation step, a function for the expectation of the log-likelihood is computed using the estimates of the current parameters. In the Maximization step, estimates of parameters are updated by maximizing the expected log likelihood function.

On the one hand, the main drawback of the EM-algorithm is that it converges slowly and needs heavy computation. Here the convergence refers to finding maximum likelihood estimator of parameters, not the PPP convergence time. However, for example, its convergence can be accelerated using the Aitken method or conjugate method [20].

On the other hand, it is fairly simple, and has robust convergence and deals conveniently with problems having a lot of parameters. For such problems, it is often the only algorithm to a large extent [19]. It is also capable of finding Kalman parameters even if we have missing data. In addition, it can detect outliers by introducing small weights for large outliers and can even estimate the outliers [21]. In contrast, the outliers are not removed but automatically downweighed in our article, since outliers sometimes take some useful information.

The paper is organized as follows, after clarifying the importance to calibrate *Qt* and *Rt*, the EM algorithm is introduced in the first section. Next, the state space model of PPP is reviewed from a point of machine learning and the methodology to adapt the EM-algorithm for extended Kalman filter for PPP is explained theoretically in detail. Thirdly, we compare the EM-algorithm with other methods, the static and kinematic instances are also given to demonstrate EM-PPP performance to improve the accuracy and reliability of PPP. Finally, the results are analyzed and the conclusion is drawn.

#### **2. State Space Model for PPP**

The state-space model allows to process GNSS data in a uniformed form. It is characterized as two equations: the state equation, which comprises a series of vector <sup>→</sup> *x <sup>t</sup>*,(1 ≤ *t* ≤ *N*, *N* is the number of epochs), and the observation equation. The state <sup>→</sup> *x t* cannot be observed directly, usually called hidden states, which is driven by hidden process noise. In this article the state-space model is described as the Kalman filter.

#### *2.1. State Equation*

The hidden state <sup>→</sup> *x <sup>t</sup>* of multi-GNSS PPP Kalman filter involves five types of parameters: three components of receiver coordinates, receiver clock error, system time difference w.r.t. GPS, troposphere zenith wet delay and ambiguities. Using subscript *t* to denote a specific time epoch, the state at time *t* evolves from the state at (*t* − 1) according to:

$$
\overrightarrow{\dot{\mathbf{x}}}\_{t} = \Phi\_{t}\overrightarrow{\mathbf{x}}\_{t-1} + \overrightarrow{\omega}\_{t} \tag{1}
$$

where Φ*<sup>t</sup>* is the transition matrix and <sup>→</sup> ω*<sup>t</sup>* the state process noise, which is assumed to be drawn from a zero-mean multivariate normal distribution, with covariance: <sup>→</sup> <sup>ω</sup>*<sup>t</sup>* ∼ N(0, *Qt*). Initial condition <sup>→</sup> *x* <sup>0</sup> is assumed to be a Gaussian vector with the a priori information *E*{ → *<sup>x</sup>* <sup>0</sup>} <sup>=</sup> <sup>→</sup> μ0,*Cov*( → *x* <sup>0</sup>) = *P*0.

The state transition matrix and the process noise matrix in static mode is defined for the position block:

$$
\Phi\_t = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}, \ Q\_t = \begin{bmatrix} Q\_x & 0 & 0 \\ 0 & Q\_y & 0 \\ 0 & 0 & Q\_z \end{bmatrix}\_t \tag{2}
$$

where *Qx*, *Qy* and *Qz* are the random process noise in *X*, *Y* and *Z* direction, respectively.

In addition to position parameters, the velocity parameters are also included in the state vector for our kinematic processing, whose system model for position and velocity block in the extended Kalman filter is given as:

$$\Phi\_{t} = \begin{bmatrix} 1 & \Delta t & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & \Delta t & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & \Delta t \\ 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix}; Q\_{t} = \begin{bmatrix} \frac{q\_{x}\Delta t^{3}}{3} & \frac{q\_{y}\Delta t^{2}}{2} & 0 & 0 & 0 & 0 \\ \frac{q\_{x}\Delta t^{2}}{2} & q\_{x}\Delta t & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{q\_{y}\Delta t^{3}}{3} & \frac{q\_{y}\Delta t^{2}}{2} & 0 & 0 \\ 0 & 0 & \frac{q\_{y}\Delta t^{2}}{2} & q\_{y}\Delta t & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{q\_{z}\Delta t^{3}}{3} & \frac{q\_{z}\Delta t^{2}}{2} \\ 0 & 0 & 0 & 0 & \frac{q\_{z}\Delta t^{2}}{2} & q\_{z}\Delta t \end{bmatrix} \tag{3}$$

The corresponding state vector is <sup>→</sup> *<sup>x</sup> <sup>t</sup>* <sup>=</sup> *x vx y vy z vz* for position and velocity block. The process noise matrix *Qt* is uniquely determined by *qx*, *qy* and *qz*, which are named as acceleration variance.

#### *2.2. Observation Equation*

The observation equation, that we used, is the double frequency ionosphere-free combination for multi-GNSS [22]:

$$PC\_r^{\iota S} = \rho\_r^i + \varepsilon dt\_r^G + T\_r^S + m\_r^i \cdot \varepsilon t d\_r + \varepsilon\_{r, \text{PC}}^{\iota, S} \tag{4}$$

$$\text{LC}\_{r}^{i,S} = \rho\_r^i + \text{cdt}\_r^G + T\_r^S + m\_r^i \cdot \text{ztd}\_r + \lambda\_{\text{LC}}^S \cdot \text{N}\_{r,\text{LC}}^S + \varepsilon\_{r,\text{LC}}^{i,S} \tag{5}$$

where subscript *r* indicates the receiver, superscript *i* represents the satellites, superscript *S* indicates GNSS constellation, following the convention of Rinex3.x (G: GPS, E: GALILEO, R: GLONASS and C: BEIDOU). *PCi*,*<sup>S</sup> <sup>r</sup>* and *LCi*,*<sup>S</sup> <sup>r</sup>* are the ionosphere-free combinations of code pseudo-range and phase observations (unit: m) respectively, which have already corrected satellite clock, the relativity effect, solid Earth tides, polar tides, ocean tides, phase wind up and a priori troposphere delay using troposphere model [23,24]. ρ*<sup>i</sup> <sup>r</sup>* is the geometry distance between receiver and satellite, *cdt<sup>G</sup> <sup>r</sup>* is receiver clock (unit: m), the superscript *G* of *cdtG <sup>r</sup>* implies that GPS time is selected as the reference time, *T<sup>S</sup> <sup>r</sup>* is the system time difference in meters of system *S* to GPS time. Specifically, for *S* = *G*, *T<sup>G</sup> <sup>r</sup>* is zero. *m<sup>i</sup> <sup>r</sup>* is troposphere mapping function and *ztdr* is troposphere zenith wet delay, λ*<sup>S</sup> LC* is the wavelength of LC combination corresponding to system *S* and *N<sup>S</sup> <sup>r</sup>*,*LC* is the LC ambiguity. <sup>ε</sup>*i*,*<sup>S</sup> <sup>r</sup>*,*PC* and <sup>ε</sup>*i*,*<sup>S</sup> <sup>r</sup>*,*LC* indicate other unmodeled errors or noise.

Equations (4) and (5) are nonlinear, the extended Kalman filter (EKF) can be used for nonlinear state estimation. For easy description, they are rewritten in a general form:

$$
\overrightarrow{\boldsymbol{y}}\_t = \boldsymbol{h}(\overrightarrow{\boldsymbol{x}}\_t) + \overrightarrow{\boldsymbol{v}}\_t \tag{6}
$$

where <sup>→</sup> *yt* <sup>=</sup> *y*<sup>1</sup> ··· *yj* ... *ykt* is the observation vector, *kt* is the number of observations at epoch *<sup>t</sup>*, *yj* ∈ {*PCi*,*<sup>S</sup> <sup>r</sup>* , *LCi*,*<sup>S</sup> <sup>r</sup>* }, → *v <sup>t</sup>* is the observation noise satisfying <sup>→</sup> *v <sup>t</sup>* ∼ N(0,*Rt*), *Rt* is the observation noise covariance matrix at epoch *t*.

#### *2.3. Kalman Filter*

Let *Ym* <sup>=</sup> <sup>→</sup> *y*1, ... , → *ym* denote all observations from epoch 1 to epoch *m*, and <sup>→</sup> *x <sup>t</sup>*|*<sup>m</sup>* represent the estimate of <sup>→</sup> *x <sup>t</sup>* given observations *Ym*, we have predicted state estimate and predicted covariance estimate:

$$
\overrightarrow{\dot{\boldsymbol{x}}}\_{t|t-1} = \Phi\_t \overrightarrow{\dot{\boldsymbol{x}}}\_{t-1|t-1} \tag{7}
$$

$$P\_{t|t-1} = \Phi\_t P\_{t-1|t-1} \Phi\_t' + Q\_t \tag{8}$$

After linearization of Equation (6) at predicted state <sup>→</sup> *x <sup>t</sup>*|*t*−1,

$$
\overrightarrow{\boldsymbol{\varepsilon}}\_{t|t-1} \approx H\_t \left( \overrightarrow{\mathbf{x}}\_t - \overrightarrow{\mathbf{x}}\_{t|t-1} \right) + \overrightarrow{\boldsymbol{\upsilon}}\_t \tag{9}
$$

where *Ht* <sup>=</sup> <sup>∂</sup> → *y* ∂ → *x* → *x <sup>t</sup>*|*t*−<sup>1</sup> , <sup>→</sup> *et*|*t*−<sup>1</sup> <sup>=</sup> <sup>→</sup> *yt* − *h* → *x <sup>t</sup>*|*t*−<sup>1</sup> . <sup>→</sup> *et*|*t*−<sup>1</sup> is called innovations or measurement residuals, then the Kalman filter is obtained:

$$K\_t = P\_{t|t-1} H\_t' \left( H\_t P\_{t|t-1} H\_t' + R\_l \right)^{-1} \tag{10}$$

$$
\overrightarrow{\dot{\boldsymbol{x}}}\_{t|t} = \overrightarrow{\dot{\boldsymbol{x}}}\_{t|t-1} + \boldsymbol{K}\_t \ \overrightarrow{\boldsymbol{e}}\_{t|t-1} \tag{11}
$$

$$P\_{t|t} = (I - K\_t H\_t) P\_{t|t-1} \tag{12}$$

where <sup>→</sup> *x <sup>t</sup>*|*<sup>t</sup>* and *Pt*|*<sup>t</sup>* are the updated Kalman estimate and the updated covariance estimate.

#### *2.4. Kalman Smoothing*

The Kalman smoother estimator could be obtained [25]:

$$
\overrightarrow{\mathbf{x}}\_{t-1|N} = \overrightarrow{\mathbf{x}}\_{t-1|t-1} + f\_{t-1} \left\{ \overrightarrow{\mathbf{x}}\_{t|N} - \overrightarrow{\mathbf{x}}\_{t|t-1} \right\} \tag{13}
$$

$$P\_{t-1|N} = P\_{t-1|t-1} + f\_{t-1}(P\_{t|N} - P\_{t|t-1})f\_{t-1}'\tag{14}$$

where *Jt*−<sup>1</sup> = *Pt*−1|*t*−1φ*<sup>t</sup>* [*Pt*|*t*−1] <sup>−</sup><sup>1</sup> , 1 <sup>≤</sup> *<sup>t</sup>* <sup>≤</sup> *<sup>N</sup>*, *<sup>N</sup>* is the number of epochs.

Kalman lag-one covariance holds with the starting condition

$$P\_{N,N-1\mid N} = \left(I - K\_N H\_N\right) \Phi\_N P\_{N-1\mid N-1} \tag{15}$$

for *t* = *N*, *N* − 1, ... , 2

$$P\_{t-1,t-2|N} = P\_{t-1|t-1} l\_{t-2}' + f\_{t-1} (P\_{t,t-1|N} - \Phi\_t P\_{t-1|t-1}) l\_{t-2}' \tag{16}$$

#### **3. EM-PPP**

The EM-algorithm is based on the innovation of the likelihood function to compute maximum likelihood estimation [25,26]. The likelihood function describes the probability of the observations, given a set of parameters. The parameters are found such that they maximize the likelihood function. The derivative of the likelihood function or log-likelihood is not always tractable. Therefore, iterative methods like Expectation-Maximization algorithms are very effective to find numerical solutions for the parameter estimates.

Denoting Θ = { → μ0, *P*0, *Qt*,*Rt <sup>t</sup>* <sup>=</sup> 1, ... , *<sup>N</sup>*}, *<sup>X</sup>* <sup>=</sup> { → *x* 0, → *x* 1, ... , → *x <sup>N</sup>*}, *Y* = { → *y*1, ... , → *yN*}. *Y* is thought of as incomplete data, and {*X*, *Y*} as complete data. Specifically for PPP, the log likelihood of the parameters of the state space model is approximately derived (ignoring constant):

$$\begin{split} 2\log L\_{X,Y}(\boldsymbol{\Theta}) &= -\log \Big| \Sigma\_{0} \Big| - \left( \overrightarrow{\mathbf{x}}\_{0} - \overrightarrow{\boldsymbol{\mu}}\_{0} \right)' \Sigma\_{0}^{-1} (\overrightarrow{\mathbf{x}}\_{0} - \overrightarrow{\boldsymbol{\mu}}\_{0}) \\ &- \sum\_{t=1}^{N} \log |Q\_{t}| - \sum\_{t=1}^{N} \left( \overrightarrow{\mathbf{x}}\_{t} - \boldsymbol{\Phi}\_{t} \overrightarrow{\mathbf{x}}\_{t} \right)' Q\_{t}^{-1} (\overrightarrow{\mathbf{x}}\_{t} - \boldsymbol{\Phi}\_{t} \overrightarrow{\mathbf{x}}\_{t}) \\ &- \sum\_{t=1}^{N} \log |R\_{t}| - \sum\_{t=1}^{N} \left( \overrightarrow{\mathbf{y}}\_{t} - h(\overrightarrow{\mathbf{x}}\_{t|N}) \right)' R\_{t}^{-1} (\overrightarrow{\mathbf{y}}\_{t} - h(\overrightarrow{\mathbf{x}}\_{t|N})) \end{split} \tag{17}$$

Since the hidden states <sup>→</sup> *x t* are unknown, only the expected value of the log likelihood conditioned on *Y* is accessible, as a result, the observation equation is expanded at smoother point <sup>→</sup> *x <sup>t</sup>*|*N*.

The Expectation (E-step) of EM algorithm for PPP requires computing the expected log-likelihood at the *jth* iteration:

$$\Omega(\Theta|\Theta^{(j-1)}) = E\langle 2\log L\_{X,Y}(\Theta)|Y, \Theta^{(j-1)}\rangle \tag{18}$$

then the parameters are recalculated at the Maximization step (M-step):

$$\Theta^{(j)} = \arg\max\_{\Theta} \Omega\Big(\Theta \big| \Theta^{(j-1)}\Big) \tag{19}$$

The two steps are repeated until the Θ(*j*) converges.

The EM-PPP is terminated when the following convergence criterion is reached:

$$R\text{-log} = \left| \frac{\ell^{(j)} - \ell^{(j-1)}}{\ell^{(j)}} \right| < \varepsilon \text{ or } j \ge \text{maximum number of iterations} \tag{20}$$

where ε is a small predefined amount and (*j*) is equal to

$$\mathcal{L}^{(j)} = \sum\_{t=1}^{N} \log \left| H\_t P\_{t|t-1} H\_t' + R\_t \right| + \sum\_{t=1}^{N} \left( \overrightarrow{\mathcal{c}}\_{t|t-1} \right)' \left( H\_t P\_{t|t-1} H\_t' + R\_t \right)^{-1} \left( \overrightarrow{\mathcal{c}}\_{t|t-1} \right) \tag{21}$$

A flowchart of our EM-PPP procedure is shown in Figure 1. In the initialization step, GNSS data preprocessing is performed including data integrity checking, cycle slips and outliers detection, phase center offset (PCO) and phase center variations (PCV) correction, synchronization of receiver clock using only code range measurements and initialization of parameters Θ. It also initializes the hidden state <sup>→</sup> *x* 0, sharing of the same processing noise *Qt* across the time step *t* = 1, ... , *N*.

**Figure 1.** Flowchart of the Expectation-Maximization (EM)-Precise Point Positioning (PPP) loop.

In the next step, the extended Kalman filter Equations (10)–(12) are implemented to compute a series of hidden states and their covariance. If the convergence condition Equation (20) is not satisfied, then Kalman smoothing is used to calculate smoothed state estimates and involved covariance matrix Equations (13)–(16), which is prepared for the E-step: calculation of the expected log likelihood function Equation (22). All parameters of Θ are updated during the M-step and prepared for the next iteration Equation (27).

#### *3.1. E-Step*

Taking the expectation upon Equation (17) over conditional distribution of the latent given observed data, we find immediately:

$$\begin{split} \Omega \left( \Theta \middle| \Theta^{(j-1)} \right) &= E \Big\{-2ln L\_{X,Y}(\Theta) \middle| Y \right\} \\ &= \ln[\Sigma\_{0}] + tr[\Sigma\_{0}^{-1}E\big( (\overrightarrow{\boldsymbol{x}\_{0}} - \overrightarrow{\boldsymbol{\mu}}\_{0})(\overrightarrow{\boldsymbol{x}\_{0}} - \overrightarrow{\boldsymbol{\mu}}\_{0})^{\prime} \big| Y \big] + \sum\_{t=1}^{N} \ln[Q\_{t}] \\ &+ \sum\_{t=1}^{N} tr[Q\_{t}^{-1}E\big( (\overrightarrow{\boldsymbol{x}\_{t}} - \Phi\_{t} \overrightarrow{\boldsymbol{x}\_{t-1}})(\overrightarrow{\boldsymbol{x}\_{t}} - \Phi\_{t} \overrightarrow{\boldsymbol{x}\_{t-1}})^{\prime} \big| Y \big] + \sum\_{t=1}^{N} \ln[R\_{t}] \\ &+ \sum\_{t=1}^{N} tr[R\_{t}^{-1}E\big( (\overrightarrow{\boldsymbol{y}}\_{t} - h(\overrightarrow{\boldsymbol{x}\_{t}}))(\overrightarrow{\boldsymbol{y}}\_{t} - h(\overrightarrow{\boldsymbol{x}\_{t}})) \big| \Big{)} \Big{]} \end{split} \tag{22}$$

where

$$E((\overrightarrow{\mathbf{x}}\_0 - \overrightarrow{\boldsymbol{\mu}}\_0)(\overrightarrow{\mathbf{x}}\_0 - \overrightarrow{\boldsymbol{\mu}}\_0)' | \mathbf{Y}) = P\_{0 \mid \mathbf{N}} + (\overrightarrow{\mathbf{x}}\_{0 \mid \mathbf{N}} - \overrightarrow{\boldsymbol{\mu}}\_0)(\overrightarrow{\mathbf{x}}\_{0 \mid \mathbf{N}} - \overrightarrow{\boldsymbol{\mu}}\_0)' \tag{23}$$

$$\begin{array}{ll} E\{\left(\overrightarrow{\mathbf{x}}\_{t} - \Phi\_{t}\overrightarrow{\mathbf{x}}\_{t-1}\right) \mid \left(\overrightarrow{\mathbf{x}}\_{t} - \Phi\_{t}\overrightarrow{\mathbf{x}}\_{t-1}\right)'\}\} \\ \qquad \to \begin{array}{ll} P\_{t|\mathcal{N}} + \overrightarrow{\mathbf{x}}\_{t|\mathcal{N}}\overrightarrow{\mathbf{x}}\_{t|\mathcal{N}} + \Phi\_{t}(P\_{t-1|\mathcal{N}} + \overrightarrow{\mathbf{x}}\_{t-1|\mathcal{N}}\overrightarrow{\mathbf{x}}\_{t-1|\mathcal{N}})\Phi\_{t}\mathsf{y} \\ \to \begin{array}{ll} \to & \overset{\sim}{\to} & \overset{\sim}{\to} & \overset{\sim}{\to} & \overset{\sim}{\to} \end{array} \end{array} \end{array} \tag{24}$$

Using Taylor series expression

$$
\overrightarrow{\boldsymbol{y}}\_{t} \approx h\left(\overrightarrow{\boldsymbol{x}}\_{t|\mathcal{N}}\right) + H\_{t|\mathcal{N}}\left(\overrightarrow{\boldsymbol{x}}\_{t} - \overrightarrow{\boldsymbol{x}}\_{t|\mathcal{N}}\right) + \overrightarrow{\boldsymbol{\upsilon}}\_{t\prime} \tag{25}
$$

where *Ht*|*<sup>N</sup>* <sup>=</sup> <sup>∂</sup> → *y* ∂ → *x* → *x <sup>t</sup>*|*<sup>N</sup>* , and let <sup>→</sup> *et*|*<sup>N</sup>* <sup>=</sup> <sup>→</sup> *yt* − *h* → *x <sup>t</sup>*|*<sup>N</sup>* , Δ<sup>→</sup> *<sup>x</sup> <sup>t</sup>*|*<sup>N</sup>* <sup>=</sup> <sup>→</sup> *<sup>x</sup> <sup>t</sup>* <sup>−</sup> <sup>→</sup> *x <sup>t</sup>*|*N*, we get

$$E\left\{ \left( \overrightarrow{\boldsymbol{y}}\_{t} - \boldsymbol{h}\left( \overrightarrow{\boldsymbol{x}}\_{t} \right) \right) \left( \overrightarrow{\boldsymbol{y}}\_{t} - \boldsymbol{h}\left( \overrightarrow{\boldsymbol{x}}\_{t} \right) \right)' \middle| \boldsymbol{Y} \right\} \approx \overrightarrow{\boldsymbol{\varepsilon}}\_{t \vert \mathbb{N}} \overrightarrow{\boldsymbol{\varepsilon}}\_{t \vert \mathbb{N}}^{\prime} + H\_{t \vert \mathbb{N}} P\_{t \vert \mathbb{N}} H\_{t \vert \mathbb{N}}^{\prime} \tag{26}$$

#### *3.2. M-Step*

Similar to complete-data weighted maximum likelihood estimation, from the first differential of Ω Θ Θ(*j*−1) , the maximum likelihood estimators are updated as follows:

$$\begin{array}{c} \stackrel{\rightarrow}{\mu}\_{0} = \stackrel{\rightarrow}{\mathbf{x}}\_{0|N} \\ \Sigma\_{0} = P\_{0|N} \end{array} \tag{27}$$

$$\begin{array}{c} Q\_{t} = E \{ (\stackrel{\rightarrow}{\mathbf{x}}\_{t} - \Phi\_{t} \stackrel{\rightarrow}{\mathbf{x}}\_{t-1}) (\stackrel{\rightarrow}{\mathbf{x}}\_{t} - \Phi\_{t} \stackrel{\rightarrow}{\mathbf{x}}\_{t-1}) \} \mid Y \} \\ Q\_{t} = E \{ (\stackrel{\rightarrow}{\mathbf{x}}\_{t} - \Phi\_{t} \stackrel{\rightarrow}{\mathbf{x}}\_{t-1}) (\stackrel{\rightarrow}{\mathbf{x}}\_{t} - \Phi\_{t} \stackrel{\rightarrow}{\mathbf{x}}\_{t-1}) \} \mid Y \} \end{array} \tag{27}$$

For simplicity, the initial covariance *P*0, and the measurement covariance *Rt* are assumed to be a diagonal matrix:

$$\begin{aligned} P\_0 &= \operatorname{diag} \left( q\_{01}, q\_{02}, \dots, q\_{0k\_0} \right), \\ R\_t &= \operatorname{diag} \left( r\_{t1}, r\_{t2}, \dots, r\_{tm\_t} \right) \end{aligned}$$

where *k*<sup>0</sup> indicates the dimension of the hidden state vector at initial epoch and *mt* is the dimension of the observation vector at epoch *t*.

#### **4. EM Compared to MLE and LS-VCE**

In literature [19], a comprehensive comparison is demonstrated between different estimation principles such as LS-VCE, best linear unbiased estimator (BLUE), best invariant quadratic unbiased estimator (BIQUE), minimum norm quadratic unbiased estimator (MINQUE) and restricted maximum likelihood estimator (REML). As shown previously in Section 3.1, the EM algorithm may be thought of as maximum likelihood estimation (MLE), but which finds the ML estimator in an iterative way. EM can be realized based on REML as well. Therefore, an additional comparison between EM variance estimation and LS-VCE and MLE is adequate.

To make theoretical analysis easy and consistent, in the following we first introduce how to covert Kalman filter to least-squares. Then we directly give different (co)variance estimators according to their distribution assumptions and the reason why the EM-algorithm is preferable in our solution.

#### *4.1. From Kalman Filter to Least Squares*

The linear (extended) Kalman equation can be transformed into the least-squares function model, which allows the following EM algorithm to be compared with LS-VCE on the same function model

and makes the theoretical analysis easy and convenient. To do so, the state Equation (1) and the observation Equation (6) are expanded at a priori value and organized as the function model:

$$\stackrel{\rightarrow}{y} = A\stackrel{\rightarrow}{x} + \stackrel{\rightarrow}{w}$$

$$E\{\stackrel{\rightarrow}{w}\} = 0, \; D\{\stackrel{\rightarrow}{w}\} = E\{\stackrel{\rightarrow}{w}\stackrel{\rightarrow}{w}'\} = Q = Q\_0 + \sum\_{k=1}^{p} \sigma\_k Q\_k \tag{28}$$

$$\begin{array}{c} \begin{array}{ccccc} \overrightarrow{y} = \begin{bmatrix} 0 & d\overrightarrow{y} \end{bmatrix}' & 0 & \cdots & 0 & d\overrightarrow{y} \end{bmatrix}' \\\\ A = \begin{bmatrix} \Phi\_1 & -1 & 0 & \cdots & 0 & 0\\ 0 & H\_1 & 0 & \cdots & 0 & 0\\ 0 & \Phi\_2 & -1 & \cdots & 0 & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & 0 & \cdots & \Phi\_N & -1\\ 0 & 0 & 0 & \cdots & 0 & H\_N \end{bmatrix} \\\\ \overrightarrow{x} = \begin{bmatrix} \overrightarrow{d}\overrightarrow{x}\_0' & \overrightarrow{d}\overrightarrow{x}\_1' & \overrightarrow{d}\overrightarrow{x}\_2' & \cdots & \overrightarrow{d}\overrightarrow{x}\_{N-1}' & \overrightarrow{d}\overrightarrow{x}\_N' & \cdots \end{bmatrix} \end{array}$$

where the a priori <sup>→</sup> *x* 0 *<sup>i</sup>* ,(*<sup>i</sup>* <sup>=</sup> 0, ... , *<sup>N</sup>*) value is subtracted from the original state vector <sup>→</sup> *x <sup>i</sup>*, leading to *d* → *<sup>x</sup> <sup>i</sup>* = <sup>→</sup> *<sup>x</sup> <sup>i</sup>* <sup>−</sup> <sup>→</sup> *x* 0 *<sup>i</sup>* and *d* → *yi* <sup>=</sup> <sup>→</sup> *yi* <sup>−</sup> <sup>→</sup> *y* 0 *<sup>i</sup>* , *N* is the number of epochs. The *m* × *n* matrix A is full column rank. The cofactor matrices *Qk*, are assumed to be known and their weighted sum *Q*<sup>0</sup> + *p k*=1 σ*kQk* is assumed to be positive definite and *Qk*,(*k* = 1, ... *p*) are linearly independent, and the (co)variance components σ*<sup>k</sup>* are unknown parameters. Matrix *Q*<sup>0</sup> is the known part of the variance matrix [19].

#### *4.2. Least-Squares EM*

Similar to Section 3, we can calculate the non-constant part of the full log-likelihood function and then conditional expectation on the observation <sup>→</sup> *y* and *Q*(*j*), given the data <sup>→</sup> *y* and the *jth* iteration estimates of (co)variance components <sup>σ</sup>(*j*) *<sup>k</sup>* or *<sup>Q</sup>*(*j*):

$$\begin{split} \Omega(\Omega|\overrightarrow{\boldsymbol{y}},\boldsymbol{Q}^{(j)}) &= \operatorname{E} \|\overrightarrow{\boldsymbol{\omega}}\|\_{\boldsymbol{\mathcal{Y}}} \boldsymbol{Q}^{(j)} \| = \log \boldsymbol{Q} + \operatorname{E} \|\boldsymbol{\pi}(\overrightarrow{\boldsymbol{\up{y}}} - \boldsymbol{A}\overrightarrow{\boldsymbol{\up{x}}})^{\prime} \boldsymbol{Q}^{-1} (\overrightarrow{\boldsymbol{y}} - \boldsymbol{A}\overrightarrow{\boldsymbol{\up{x}}})) \| \overrightarrow{\boldsymbol{y}} \| \\ &= \log \boldsymbol{Q} + \operatorname{tr}(\mathbf{Q}^{-1} \boldsymbol{E} ((\overrightarrow{\boldsymbol{y}} - \boldsymbol{A}\overrightarrow{\boldsymbol{\up{x}}}) (\overrightarrow{\boldsymbol{y}} - \boldsymbol{A}\overrightarrow{\boldsymbol{\up{x}}})^{\prime} \overrightarrow{\boldsymbol{y}})) \end{split} \tag{29}$$

where

$$\begin{aligned} E(\overrightarrow{\boldsymbol{y}} - \boldsymbol{A}\overrightarrow{\boldsymbol{x}})(\overrightarrow{\boldsymbol{y}} - \boldsymbol{A}\overrightarrow{\boldsymbol{x}})' \|\overrightarrow{\boldsymbol{y}}\| &= \left(\overrightarrow{\boldsymbol{y}} - \boldsymbol{A}\hat{\boldsymbol{x}}\right)(\overrightarrow{\boldsymbol{y}} - \boldsymbol{A}\hat{\boldsymbol{x}})' + \boldsymbol{A}Q\_{\boldsymbol{\lambda}}^{(j)}\boldsymbol{A}'\boldsymbol{Z}\_{\text{ML}}^{(j)}\\ \pounds &= \left(\mathbf{A}'(\mathbf{Q}^{(j)})\,^{-1}\boldsymbol{A}\right)^{-1}\boldsymbol{A}'(\mathbf{Q}^{(j)})^{-1}\overrightarrow{\boldsymbol{y}} \end{aligned} \tag{30}$$

M-step:

Maximizing the likelihood of the completed data based on Equation (29), the new estimates <sup>σ</sup>(*j*+1) *k* are calculated as

$$\begin{bmatrix} \sigma\_1^{(j+1)} \\ \vdots \\ \sigma\_p^{(j+1)} \end{bmatrix} = \left( A'\_{vEM} W\_{vEM} A\_{vEM} \right)^{-1} A'\_{vEM} W\_{vEM} \text{vec} \left( Z\_{EM}^{(j)} - Q\_0 \right) \tag{31}$$

where *AvEM* = " *vec <sup>Q</sup>*(*j*) 1 ··· *vec <sup>Q</sup>*(*j*) *p* #, *WvEM* = (Q(*j*) ) −1 <sup>⊗</sup> (Q(*j*) ) −1 , ⊗ is the Kronecker product and vec is *vec*-operator.

Equations (29) and (31) are the EM algorithm for ML estimation. If convergence is reached, set <sup>σ</sup>*<sup>k</sup>* <sup>=</sup> <sup>σ</sup>(*j*) *<sup>k</sup>* , otherwise increase *j* by one and return to E-step.

#### *4.3. MLE*

Once the general structure of probability density function is known, MLE can be simply realized and therefore used widely. If a multivariate normal distribution is given, the (co)variance components takes form:

$$
\begin{bmatrix}
\sigma\_1 \\ \vdots \\ \sigma\_p
\end{bmatrix} = \left(A'\_{\text{vMLE}} \mathcal{W}\_{\text{vMLE}} A\_{\text{vMLE}}\right)^{-1} A'\_{\text{vMLE}} \mathcal{W}\_{\text{vMLE}} \text{vec}(Z\_{\text{MLE}} - Q\_0) \tag{32}
$$

where *AvMLE* <sup>=</sup> *vecQ*<sup>1</sup> ··· *vecQp* , *WvMLE* = *<sup>Q</sup>*−<sup>1</sup> <sup>⊗</sup> *<sup>Q</sup>*<sup>−</sup>1.

From Equations (31) and (32), we know that least-squares EM and MLE estimators share the same design matrix and weight matrix. Their difference is mainly caused by the pseudo observation *vec <sup>Z</sup>*(*j*) *EM* − *Q*<sup>0</sup> and *vec*(*ZMLE* <sup>−</sup> *<sup>Q</sup>*0). *<sup>Z</sup>*(*j*) *EM* for EM-algorithm includes the effects of both observation post-residuals and the accuracy of the estimates *x*ˆ. In contrast, *ZMLE* for MLE consider only the observation post-residuals. Therefore, MLE estimator is probably over-optimistic to EM.

However, if the REML principle is used to derive the EM-algorithm, the effect of *x*ˆ is implicitly removed. Then, the EM-algorithm based on REML will be equivalent to the REML estimator.

#### *4.4. LS-VCE*

Another important problem is that the EM-algorithm and MLE do not take the loss of degrees of freedom from the estimation of <sup>→</sup> *x* into account. Borrowing the idea of REML, LS-VCE overcomes this problem based on (*n* − *p*) independently error contrasts. Specifically, let

$$
\hat{\mathfrak{f}} = B' \overrightarrow{\hat{\mathfrak{y}}}, \hat{\mathbb{Q}}\_{\mathbb{I}} = B' \textbf{Q} \textbf{B} \tag{33}
$$

where <sup>→</sup> *t* is misclosure vector, *B* is *m* × (*m* − *n*) matrix satisfying *B A* = 0, *rank*(*B*) = *m* − *n* = *b*. Then LS-VCE estimator is obtained:

$$\begin{bmatrix} \sigma\_1\\ \vdots\\ \sigma\_p \end{bmatrix} = \left( A'\_{\text{vLS-VCE}} W\_{\text{vLS-VCE}} A\_{\text{vLS-VCE}} \right)^{-1} A'\_{\text{vLS-VCE}} W\_{\text{vLS-VCE}} \text{vec}(Z\_{LS-VCE} - Q\_0) \tag{34}$$

with *AvLS*−*VCE* <sup>=</sup> *vec*(*B <sup>Q</sup>*1*B*) ··· *vec B QpB* , *ZLS*−*VCE* <sup>=</sup> <sup>ˆ</sup>*t*ˆ*<sup>t</sup>* , *WvLS*−*VCE* is a user-defined weight matrix. If *WvLS*−*VCE* is set to *Q*ˆ*<sup>t</sup>* ⊗ *Q*ˆ*t*, LS-VCE is the same as EM-algorithm based on REML.

LS-VCE is derived purely based on the least-squares method and we do not make any assumption on a probability density function (PDF). In contrast, EM-algorithm, MLE and REML are built upon a certain distribution, which explains why when applying LS-VCE, it is necessary for users to set weight matrix on their own.

#### *4.5. Preference for Recursive EM*

As discussed previously, the EM-algorithm can be implemented as either recursive form or batch form like MLE and LS-VCE. In our solution, we prefer the EM-algorithm based Kalman filter to other methods.

Recursive EM discriminates between the process noise and the observation noise. For GNSS, the process noise is usually different from observation noise. The process noise is directly connected to the geophysical phenomenon, which has not only linear but also non-linear variations, and suffers both time and spatial correlations [27–29]. As a result, it is relatively more difficult to estimate the process noise than the observation covariance matrix. Other batch methods mix the two types of stochastic processes with different behavior, which will bring us extra trouble.

Recursive EM can process data in not only post mode, but also in real-time mode. In both modes, data is processed epoch by epoch, allowing us to dynamically adjust the weight matrix, monitoring time-varying behavior and detecting abrupt motion.

#### **5. Validation**

#### *5.1. Static PPP Scheme*

A total of 14 IGS multi-GNSS stations are selected to assess the performance of EM-PPP (Figure 2). Those stations are evenly distributed on the Earth and track as many GNSS constellations as possible, covering not only the continent but also coastal and polar regions.

**Figure 2.** Distribution of selected IGS multi-GNSS stations.

Daily GNSS measurements from those IGS stations observed during DOY 119 in 2017 are used in this study. The true coordinate benchmarks are from IGS weekly solutions. The final GFZ Beidou multi-GNSS (GBM) products of the satellite orbits and clocks are applied (ftp://cddis.gsfc.nasa.gov/ pub/gps/products/mgex/). The used precise orbit correction has a sampling interval of 15 min while the precise clock has a sampling interval of 30 s, both generated at GFZ. For GPS, we aligned C1 to P1 using CODE differential code bias (DCB) product (ftp://ftp.aiub.unibe.ch/CODE/2017/). The frequency bands we used are L1 and L2 for GPS, G1 and G2 for Glonass, E1 and E5a for Galileo, B1 and B2 for BeiDou [9]. Receiver and satellite PCO and PCV were corrected using igs14.atx. solid Earth tides, pole tide and ocean tides are removed according to IERS Conventions 2010. For troposphere delay estimation, the zenith dry component of tropospheric delays was corrected with the a priori Niell model [24]. The zenith wet delay (ZWD) is estimated as an unknown parameter. Then, 24-h observation data sets with a sampling interval of 30 s were processed for all stations. The elevation cut-off was set to 6 degrees.

The initial guess of receiver coordinates is intentionally deviated by 100 m from IGS station solution. A priori standard deviation (std) of PC is set to 3 m, a priori std of LC to 0.03 m for pseudo-range combination (PC) and carrier phase combination (LC) combination, respectively.

The starting values for *Qt* are shown in Table 1. Random walk noise process with a spectral density equal to 1.0 mm/ √ 30 s is imposed on coordinates, which means a 1.0 mm disturbance very 30 s for IGS station in North, East and Up. It is not true in reality of course, but useful for test purposes. The receiver clock offset is supposed to be white noise. Zenith wet delay (ZWD) and inter-system bias (ISB) are modeled as random walk noise. Ambiguities can be considered as constant or random walk noise with very small spectral density.


**Table 1.** Initial guess *Qt* for static PPP: ZWD is the zenith wet delay of the troposphere, ISB is the system time difference with respect to GPS time.

If the maximum number of iteration is reached and EM-PPP does not converge, smaller spectral density should be assigned for *X*, *Y* and *Z* for the next cycle.

#### *5.2. Kinematic PPP Scheme*

Another kinematic dataset was used to further validate the performance of the EM-PPP. The data was collected at Wuhan, China, in November 14, 2013. The sampling interval is 1 Hz and the observed time span is about one hour. The final CODE precise satellite orbit and 5 s clock products are used to estimate the 1 Hz GPS displacements. The ambiguity-fixed double differenced real-time kinematic (RTK) solutions are adopted as the reference to assess the performance of kinematic EM-PPP solution.

The initial acceleration variance is assumed to be 10 m2 s3 for position and velocity states (Table 2), which can be used to calculate the process noise matrix for position and velocity. The initial *Qt* for receiver clock is modeled as white noise and estimated epoch-wisely. ZWD and ambiguities are also modeled as random walk processes with initial spectral density 0.01 cm/ √ *s* and 0.0 m/ √ *s*, respectively.

**Table 2.** Initial guess *Qt* for kinematic PPP: ZWD is zenith wet delay of troposphere.


#### **6. Results and Discussion**

#### *6.1. Static EM-PPP Solution*

It is found that EM-PPP usually converges after the iteration counter reaches 50. The positioning errors are shown in Table 3, including PPP results at 1st iteration with the biased stochastic model, and the results after 50 iterations of calibration to assess EM-PPP performance. EM-PPP convergence in our research means that the square root of 3D positioning errors of the last 20 consecutive epochs is less than 10 cm.

Table 3 indicates that when the biased *Qt* and *Rt* are fed in the beginning, PPP 3D errors are up to decimeters for a few stations. Horizontal errors are often greater than vertical errors, which is not consistent with the property of GNSS, because of inappropriate process noise matrix and measurement noise covariance matrix. After 50 iterations, the position errors are reduced to within 1 cm in North, East and Up direction on average using our EM-PPP algorithm. The mean 3D error is reduced to 1.77 cm without fixing ambiguities. The overall decrease percentage on average is 66.91%, 66.16%, 71.60% in North, East and Up direction, respectively, 69.95% for 3D errors.


**Table 3.** Statistics of EM-PPP absolute positioning errors (cm) with respect to IGS stations solution. All sites except YEL2 are processed using multiple GNSS observations. YEL2 is processed using only GPS data to verify the algorithm for single GNSS constellation. The ambiguities are not fixed.

To see what happened to the process noise *Qt* and the observation covariance matrix *Rt* before and after calibration, an example of JFNG station located in China is illustrated.

The residuals for PC and the corresponding formal errors are shown in Figure 3. The residuals for LC and the corresponding formal errors against satellite elevation angles after calibration are shown in Figure 4. To be clear, PC and LC residuals for BeiDou are plotted separately from those for GPS, Glonass and Galileo.

**Figure 3.** EM-PPP absolute pseudo-range combination (PC) residuals and formal errors of JFNG station at 50th iteration: Left two pictures are the PC residuals and their formal errors (square root of observation matrix *Rt*) of GPS, Glonass and Galileo, respectively. Similarly, the right two pictures show the PC residuals and formal errors of BeiDou.

**Figure 4.** EM-PPP carrier phase combination (LC) residuals and formal errors of JFNG station at 50th iteration: Left two pictures are the LC residuals and formal errors (square root of the observation covariance matrix) for GPS, Glonass and Galileo, respectively. Similarly, the right two pictures show the LC residuals and formal errors for BeiDou.

The results allow us to examine the relationship between *Rt* and the residuals. It is clearly shown that formal errors of *Rt* are changed from 3 m to vary about between 0 and 10 m for PC, and from 0.03 m to vary about between 0 and 0.6 m for LC. The big outliers of LC are due to the non-convergence at the initial epoch and EM-PPP does not fully converge at some epochs. Although *Rt* variations, similar to the residuals-varying pattern, are dependent on the satellite elevation angle, LC formal error of BeiDou *Rt* at the elevation of about 20 degrees is greater than the lower degree formal error (refer to the bottom right picture in Figure 4). Clearly, it is not advisable to choose observation weight only according to the satellite elevation angle. Another example is station ISTA, where Glonass *Rt* for PC at satellite elevation higher than 50 degrees is almost as great as lower degree errors (not plotted here). In addition, random or systematic outliers are downweighed accordingly for both PC and LC.

It can also be observed that BeiDou PC residuals peak are almost 6 m, worse than GPS, Glonass and Galileo at JFNG. In fact, the biggest error source comes from GEO C05 and IGSO C06, C07 and C08, probably due to their poor orbit accuracy and clock offset, because the residuals of those satellites show less independence of satellite elevation angle. *Rt* is different among GPS, Glonass and Galileo.

Thus *Rt* not only reflects the accuracy of measurement type itself, elevation-dependent, characteristics of GNSS constellation, but also the quality of GNSS satellite orbit and clock products, and should be adjusted dynamically. Consequently, EM-PPP is effective to calibrate *Rt* automatically and suppress outliers.

In general, *Qt* can be adjusted to its correct value, which is zero in our case, after less than ten iterations, as shown in Figure 5. The initial square root of *Qt* in North, East and Up directions are 0.001 m/ √ 30 s. *Qt* becomes zero at the sixth iteration in all three components and position solution at the last epoch varies little.

It has to be pointed out that EM-PPP suffers local extrema problems like other alternative methods. It is not an algorithm to locate the global maxima, therefore EM-PPP is sensitive to initial guesses. To escape from this local extrema trap, several sets of initialization schemes can be used, and select the best one selected as the final result.

**Figure 5.** Effects of EM-PPP iteration number on the position errors (top) and variations of the square root of *Qt* at last epoch on station JFNG.

#### *6.2. Kinematic EM-PPP Solution*

Figure 6 is the tracking route recovered by the EM-PPP at the 500th iteration. The relative log-likelihood versus iterations is plotted in Figure 7, where the results at the first iteration correspond to the solution to the traditional PPP (TPPP). The relative log likelihood decreases gradually from 1.0 to 8.665 <sup>×</sup> <sup>10</sup>−<sup>7</sup> and forms a concave curve. The less the relative log-likelihood, the less perturbative the EM-PPP solution becomes with respect to the RTK solution. As the iterative number grows to about 150, the EM-PPP solution converges nearly completely.

**Figure 6.** Driving route of moving vehicle in Wuhan, China on 14 November 2013.

**Figure 7.** Relative log-likelihood of two consecutive iterations.

Illustrated in Figure 8 are the positioning errors of the traditional PPP solution and EM-PPP solution at the 500th iteration with respect to the reference coordinates in the Up, East and North directions. It is observed that EM-PPP solutions are much more stable in all directions when compared to TPPP solutions. The TPPP solution in the East changes wobbly in comparison to the North and Up due to greater acceleration in the East (Figure 9), which leads to a larger bias in the East.

**Figure 8.** GPS-only kinematic positioning errors with respect to real-time kinematic (RTK): (**a**) traditional PPP errors, (**b**) EM-PPP positioning errors at the 500th iteration. Traditional PPP and EM-PPP share the same initial conditions.

D

E 

D

E 

**Figure 9.** Acceleration time series derived from second-order differencing RTK position time series in Up, East and North: (**a**) acceleration in Up. (**b**) Accelerations in the East and North.

In fact, the position and velocity (PV) model of our kinematic state equation is assumed to be a constant velocity model, which means the acceleration is zero. However, the realistic acceleration of the vehicle is no zero, which results in systematic bias and unstable solution.

Figures 10 and 11 describe the EM-PPP RMSs and STDs with respect to RTK against the number of iterations, respectively. Obviously, after about 40 iterations, RMSs in Up and North are decreasing and the solutions are improved in those two directions. In contrast to Up and North direction, in the East direction it presents a decrease from 9.7 cm falling to 7.6 cm and then increases up to 8.5 cm for the RMS. However, the combined effects of all three directions get decreasingly 3D RMS, proving that EM-PPP does improve the positioning accuracy in the kinematic mode in our case, and apparently converges with increasing iterations, consistent with the EM theory, though there is a little disturbance because of the existing system bias and outliers of pseudo-range. It can be imagined that a better result can be expected if acceleration observations are also observed.

**Figure 10.** EM-PPP RMSs w.r.t. RTK solution in Up, East and North (unit: m).

**Figure 11.** EM-PPP STDs w.r.t. RTK solution in Up, East and North (unit: m).

As for STDs with respect to RTK, they are consistently increasing in all three directions when EM-PPP continues its iteration. The STDs decrease from 4.0 cm, 8.4 cm and 2.5 cm to 3.5 cm, 1.6 cm and 1.6 cm in Up, East and North components, respectively. In other words, the STDs are improved by 12.5%, 80.9% and 36.0% in Up, East and North, accordingly.

D E  Given Figure 12 is the estimates of geodetic coordinates B, L and H, and their estimates of the square root of *Qt* after the 500th iteration. It is noted that the coordinates simultaneously stay stable or alternatively change sharply in all three directions, telling that the moving patterns of the vehicle switch between static and kinematic status from time to time. Theoretically, the process *Qt* should change between zero and positive values. Obviously, the estimates of the square root *Qt*, displayed in Figure 12b–f agree well with variations of coordinates. Time-varying *Qt* is identified, which takes the concrete dynamic mode into consideration.

**Figure 12.** EM-PPP estimates of the Geodetic coordinates B, L, H and their correspondingly calibrated *Qt* after the 500th iteration: (**a**) Latitude estimates B (degree). (**b**) The square root of *Qt* of B (m). (**c**) Longitude estimates L (degree). (**d**) The square root of *Qt* of L (m). (**e**) Height estimates H (m). (**f**) The square root of *Qt* of H (m).

If much smaller initial acceleration variances are given, for example 0.01 *<sup>m</sup>*<sup>2</sup> *<sup>s</sup>*<sup>3</sup> rather than <sup>10</sup> *<sup>m</sup>*<sup>2</sup> *<sup>s</sup>*<sup>3</sup> for position and velocity states, it can be shown that our algorithm can still recover the process noise matrix as well and makes no difference. As a result, the user can choose values for *Qt* randomly to a large extent.

#### **7. Conclusions**

A machine learning algorithm, EM algorithm, is adapted particularly to the extended Kalman filter to calibrate the process noise matrix and observation covariance noise matrix for PPP. The main advantage of EM-PPP is in fact that it is straightforward, simple, (locally) optimal and able to estimate large amounts of parameters and thus competent in calibrating the time-varying process noise and observation covariance for PPP state-space model, though its execution is time-consuming. The basic framework of EM-PPP is not limited to multi-PPP and can be applied to other fields of geodesy.

The whole procedure of EM-PPP is comprised of three parts: initialization, feedforward and backpropagation. In the beginning, the GNSS preprocess is performed to check the availability of the required data and mainly recognize cycle slips. Next, the whole process iterates between the estimation of hidden state and expectation and maximization.

The EM-algorithm is then compared with MLE and LS-VCE methods. We choose the recursive algorithm because it is superior to separate the process noise and observation variance, and to monitor time-varying behavior.

The approach was verified by selecting a global distribution of 14 IGS multi-GNSS station without fixing ambiguities. Based on the presented results, it concluded that EM-PPP is well suited for dynamically determining the time-varying process noise and observation noise. The calibrated observation variance matches the observation residuals from low satellite elevation angle to high satellite elevation angle. It resists orbit and clock errors and outliers through downweighing abnormal observations at different epochs, which is an alternative reasonable solution in contrast to the popular way that assigns weight according to the satellite elevation angle. People do not need to worry about separating observations into different categories carefully based on different GNSS constellations to estimate variance components like variance component estimation (VCE).

The spectral density of the assumed kinematic IGS station with 1 mm disturbance every 30 s in North, East and Up direction was estimated to be zero, implying that stations are static, which is consistent with reality.

An additional kinematic test was also implemented and reasonable values of *Qt* are found when biased initial *Qt* guess was given. The position errors are reduced in Up, East and North direction, respectively, w.r.t. RTK. In particular the STDs with respect to RTK are improved by 12.5%, 80.9% and 36.0% in Up, East and North, further showing that EM-PPP is also beneficial to kinematic PPP.

It has been confirmed that the EM-PPP is competitive for the calibration of the PPP stochastic model dynamically. The main drawback of this approach is that it converges slowly due to its first-order convergence. In the future, online EM-PPP may be derived to process GNSS data in real-time to overcome this problem if a large number of observations are available.

**Author Contributions:** X.Z. conceived and designed the algorithm; P.L. provided the kinematic GNSS data and helped validate the algorithm; M.G. supervised the whole procedure, and continuously discussed and analyzed the results and gave constructive suggestions; R.T. and X.L. participated in the experimental investigation; H.S. helped edit and revise the paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research is funded by the China Scholarship Council (CSC).

**Acknowledgments:** X.Z. is financially supported by the China Scholarship Council (CSC) for his study at the German Research Centre for Geosciences (GFZ). We thank the IGS and GFZ for providing GNSS observations, DCBs, precise orbit, and clock products. Our research is also partly supported by the Chinese Academy of Sciences (CAS) program of "Light of West China" (Grant No. 29Y607YR000103), Chinese Academy of Sciences, Russia and Ukraine and other countries of special funds for scientific and technological cooperation (Grant No. 2BY711HZ000101). This work is also partly supported by the National Natural Science Foundation of China (Grant No: 41674034, 41974032, 11903040), the Chinese Academy of Sciences (CAS) programs of "Pioneer Hundred

Talents" (Grant No: Y923YC1701), and the Chinese Academy of Sciences (CAS) program of "Western Youth Scholar" (Grant No: Y712YR4701, Y916YRa701), as well as "The Frontier Science Research Project" (Grant No: QYZDB-SSW-DQC028). We also thank the reviewers for their careful review.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Elementary Error Model Applied to Terrestrial Laser Scanning Measurements: Study Case Arch Dam Kops**

#### **Gabriel Kerekes \* and Volker Schwieger**

Institute of Engineering Geodesy, University of Stuttgart, Geschwister-Scholl-Str. 24D, 70174 Stuttgart, Germany; volker.schwieger@iigs.uni-stuttgart.de

**\*** Correspondence: gabriel.kerekes@iigs.uni-stuttgart.de

Received: 18 March 2020; Accepted: 10 April 2020; Published: 15 April 2020

**Abstract:** All measurements are affected by systematic and random deviations. A huge challenge is to correctly consider these effects on the results. Terrestrial laser scanners deliver point clouds that usually precede surface modeling. Therefore, stochastic information of the measured points directly influences the modeled surface quality. The elementary error model (EEM) is one method used to determine error sources impact on variances-covariance matrices (VCM). This approach assumes linear models and normal distributed deviations, despite the non-linear nature of the observations. It has been proven that in 90% of the cases, linearity can be assumed. In previous publications on the topic, EEM results were shown on simulated data sets while focusing on panorama laser scanners. Within this paper an application of the EEM is presented on a real object and a functional model is introduced for hybrid laser scanners. The focus is set on instrumental and atmospheric error sources. A different approach is used to classify the atmospheric parameters as stochastic correlating elementary errors, thus expanding the currently available EEM. Former approaches considered atmospheric parameters functional correlating elementary errors. Results highlight existing spatial correlations for varying scanner positions and different atmospheric conditions at the arch dam Kops in Austria.

**Keywords:** elementary error model; terrestrial laser scanning; variance-covariance matrix

#### **1. Introduction**

One of the main tasks in engineering geodesy is deformation and displacement monitoring of structures such as buildings, bridges, towers, dams, tunnels or other infrastructure works (cf. [1,2]). Independent of the measurement method, geodetic sensors are used to gather data either in a continuously manner or within different epochs. In both cases, these prerequisites are essential: a common geodetic reference system for all the epochs, knowledge about the deformation process and a stochastic model that describes the uncertainty of the measurements. Classical geodetic measurement methods like Global Navigation Satellite System (GNSS), total station, leveling, etc. have been used for decades in terrestrial point-wise monitoring and have well established and broadly accepted stochastic models [3]. Although highly reliable, point-wise acquisition methods have their limitations if objects with complex shapes like curved facades, high-rise buildings or arch dams, require deformation monitoring. Here is where area-wise deformation analysis covers the gap by implying measurement methods capable of remotely measuring a large area of the observed object [4]. To gain an impression of recent applications, the reader is referred to [5–7]. One recent method is Terrestrial Laser Scanning (TLS). Terrestrial Laser Scanners (TLSs) are active multi-sensor systems used to measure the three-dimensional geometry of a given surrounding within a certain range (cf. [8,9]). Laser scanners got more precise, compact and affordable in the past 20 years [10], but neither instrument manufacturers nor the scientific communities have reached common ground in what concerns all TLS influencing

error sources. This is commonly known as the TLS error budget or TLS stochastic model; which is currently still unsatisfactory [11]. Neuner et al. [12] give an overview of the available point cloud modeling methods used in engineering geodesy together with their stochastic models and state that none of them is established. Generally, a stochastic model is a mathematical model that describes real-life phenomena that are characterized by the presence of uncertainty [13]. In any case of direct and indirect measurements [14], the stochastic model can be expressed by a variance-covariance matrix (VCM) [15]. If knowledge about the existing correlations between all observations is missing, the VCM is reduced to a diagonal matrix poorly resembling the complex nature of all the error sources (cf. [11]). This consequently leads to possibly wrong decisions in the TLS deformation analysis [16] or inappropriate estimations of a specific surface. (cf. [7]).

To overcome this issue, the Elementary Error Model (EEM) can be used to define the stochastic model of TLS observations in form of a VCM that considers correlations. Previous work of Kauker and Schwieger [17] sets the foundation of applying the EEM on TLS measurements. To that point, the EEM model was applied on a TLS of panoramic type [9] and the atmospheric elementary errors were considered functional correlating. Continuing this line of work, the current contribution introduces a model for long-range hybrid type TLSs and classifies the atmospheric elementary errors as stochastic correlating for the first time. The latter is possible due to derived correlations between atmospheric parameters in the research area. Results are shown on airside point clouds of the Kops arch dam in Vorarlberg, Austria.

In the second section of this paper, the EEM theory is reviewed for comprehension. Section 3 describes the application of EEM on a Riegl VZ-2000 hybrid TLS (RIEGL Laser Measurement Systems GmbH, Horn, Austria) together with meteorological elementary errors and their influences on the distance measurements and vertical angles. The study case and outcomes are presented in Sections 4 and 5 concludes this contribution.

#### **2. Elementary Error Model Theory**

#### *2.1. General Remarks about Stochastic Models*

The purpose of a stochastic model is to describe the statistical properties of variables [18]. There are many possibilities of describing the propagation of uncertainty of these variables. Out of these, some are most commonly used in measurements; these are: Guide to the Expression of Uncertainty in Measurement (GUM) [19], Monte Carlo Method (MCM) [20] and the variance covariance propagation law (cf. [21]). Only the last two will be briefly discussed with regard to the assumed models. On one side, in the MCM *n* random variables are numerically processed without having any knowledge about neither the linear/non-linear nature of the random variables nor their statistical distribution. Based on the outcomes, the statistical distribution is derived with corresponding parameters such as expected value, standard deviation, skewness and kurtosis. One disadvantage is that the model is computed *n* times, which increases computation time drastically. For more details the reader is referred to [20,22]. On the other side, variance covariance propagation law assumes normal distributed random values and linear or linearized models. Outcomes are likewise normally distributed and the statistical parameters are completely described by the expected value and standard deviation. This is an advantageous method, since the linear or linearized functional model is only computed once [18], therefore reducing computation time. It is also the main reason of adopting it for the EEM of TLS measurements, where the observations number easily reaches a few hundred thousand or a few million. To support this hypothesis, Aichinger and Schwieger [23] proved after using MCM for TLS observations for different scanning configurations that in 90% of the cases, linear models can be assumed with a significance level of α = 0.003. Therefore, assuming a linear model for TLS observations is acceptable for most cases, even if observations have a non-linear nature. Regarding the numerical estimates introduced later, it is mentioned that no method of estimating the outcome's precision is currently used. This may be achieved in the future with the help of Variance Component Estimation

(VCE) based on [24,25], or a review of [26]. Our intention is to use sensitivity analysis (cf. [27,28]) and inspect how the input estimates influence the outcomes. All of these aspects will be prospectively presented in a different publication.

#### *2.2. Elementary Error Theory*

The general theory of the elementary error model was simultaneously defined by Hagen [29] and Bessel [30]. Later on, the model was elegantly presented by Pelzer [21] and extended by Schwieger [31]. Some of its applications can be found in exemplifying the error impact on several geodetic measurement methods like electronic distance measurement (EDM) instruments [32], GNSS observations [33] or recently TLS measurements [17].

According to the EEM theory, each realization of a measured random quantity differs from its expected value by a random deviation ε [31]. It is assumed that ε is comprised by the sum of countless, small elementary errors. Their absolute value is supposed to be equal and the probability of a positive and negative sign is likewise presumably equal [29]. The presumption of standard normal distribution of these errors is supported by an infinite number of elementary errors with infinitely small absolute values. Their impact on the observations can be modeled by using error vectors and influencing matrices. These matrices resemble the effect on the covariance matrix of observations. Three types of impacts are considered: non-correlating error vectors δ*k*, functional correlating error vector ξ and stochastic correlating error vectors γ*<sup>h</sup>* [31]. For each error type, corresponding influencing matrices are defined as follows: *p* matrices *Dk* for non-correlating errors, one matrix *F* for functional correlating errors and *q* matrices *Gh* for stochastic correlating errors. Therefore, the random deviation vector ε results as a sum of all elementary errors accordingly:

$$\mathfrak{e} = \sum\_{k=1}^{p} D\_k \cdot \delta\_k + F \cdot \xi + \sum\_{h=1}^{q} G\_h \cdot \mathcal{Y}\_h. \tag{1}$$

These influencing matrices have different structures depending on the elementary errors effects on the observations. Hereby, matrices *Dk* and *Gh* are symmetrical diagonal matrices because each elementary error of the non-correlating and stochastic correlating group influences exactly one measurement quantity functionally. The matrix *F* is fully populated because one functional correlating error may impact several measurement quantities [31]. Defining the functional relationships between observations *l***<sup>1</sup>** ... *ln* and the elementary errors δ, ξ and γ allows the calculation of the partial derivatives that populate the influencing matrixes as follows:

$$D\_{\mathbf{k}} = \begin{bmatrix} \frac{\partial l\_{1}}{\partial \delta\_{1k}} & 0 & \cdots & 0\\ 0 & \frac{\partial l\_{2}}{\partial \delta\_{2k}} & 0 & \vdots\\ \vdots & 0 & \ddots & \vdots\\ 0 & \cdots & \cdots & \frac{\partial l\_{n}}{\partial \delta\_{n}} \end{bmatrix}, \quad F = \begin{bmatrix} \frac{\partial l\_{1}}{\partial \delta\_{1}} & \frac{\partial l\_{1}}{\partial \delta\_{2}} & \cdots & \frac{\partial l\_{1}}{\partial \delta\_{n}}\\ \frac{\partial l\_{2}}{\partial \zeta\_{1}} & \frac{\partial l\_{2}}{\partial \zeta\_{2}} & \cdots & \frac{\partial l\_{2}}{\partial \zeta\_{n}}\\ \vdots & \vdots & \ddots & \vdots\\ \frac{\partial l\_{a}}{\partial \zeta\_{1}} & \frac{\partial l\_{a}}{\partial \zeta\_{2}} & \cdots & \frac{\partial l\_{a}}{\partial \zeta\_{n}} \end{bmatrix}, \quad G\_{\mathbf{k}} = \begin{bmatrix} \frac{\partial l\_{1}}{\partial \gamma\_{1k}} & 0 & \cdots & 0\\ 0 & \frac{\partial l\_{2}}{\partial \gamma\_{2k}} & 0 & \vdots\\ \vdots & 0 & \ddots & \vdots\\ 0 & \cdots & \cdots & \frac{\partial l\_{a}}{\partial \gamma\_{n}} \end{bmatrix}. \tag{2}$$

Applying the law of propagation of variance on Equation (1) yields the so called "synthetic covariance matrix" **Σ***ll* which by definition has the following form [33]:

$$
\Sigma\_{\rm II} = \sum\_{k=1}^{p} D\_k \cdot \Sigma\_{\delta\delta,k} \cdot D\_k^T + F \cdot \Sigma\_{\xi\xi} \cdot F^T + \sum\_{h=1}^{q} G\_h \cdot \Sigma\_{\mathcal{YY},h} \cdot G\_{h'}^T \tag{3}
$$

where each covariance matrix of the elementary errors is defined and structured as shown below. The covariance matrices for non-correlating errors **Σ**δδ,*<sup>k</sup>* and functional correlating errors **Σ**ξξ are diagonal matrices having the elementary error's variances on the main diagonal. As a result of

the possible covariances of the stochastic correlating errors, the corresponding matrix **Σ**γγ,*<sup>h</sup>* may be fully populated.

$$
\boldsymbol{\Sigma}\_{\mathsf{d}\mathsf{d},\mathsf{k}} = \begin{bmatrix}
\sigma\_{1\mathsf{k}}^{2} & 0 & \cdots & 0 \\
0 & \sigma\_{2\mathsf{k}}^{2} & 0 & \vdots \\
\vdots & 0 & \ddots & \vdots \\
0 & \cdots & \cdots & \sigma\_{n\mathsf{k}}^{2} \\
\end{bmatrix} \boldsymbol{\Sigma}\_{\mathsf{d}\mathsf{f}} = \begin{bmatrix}
\sigma\_{1}^{2} & 0 & \cdots & 0 \\
0 & \sigma\_{2}^{2} & 0 & \vdots \\
\vdots & 0 & \ddots & \vdots \\
0 & \cdots & \cdots & \sigma\_{m}^{2} \\
\end{bmatrix} \boldsymbol{\Sigma}\_{\mathsf{Y}\mathsf{Y}\mathsf{h}} = \begin{bmatrix}
\sigma\_{1\mathsf{h}}^{2} & \sigma\_{12\mathsf{h}} & \cdots & \sigma\_{1n\mathsf{h}} \\
\sigma\_{12\mathsf{h}} & \sigma\_{2\mathsf{h}}^{2} & \cdots & \sigma\_{2n\mathsf{h}} \\
\vdots & \vdots & \ddots & \vdots \\
\sigma\_{1n\mathsf{h}} & \cdots & \cdots & \sigma\_{n\mathsf{h}}^{2} \\
\end{bmatrix}. \tag{4}
$$

The challenging part is finding variances for all groups of errors and covariances for stochastic correlating errors. Correlations between the elementary errors are assumed to be zero. The variances may however be extracted from instrument manufacturers reports (cf. Section 3.2), empirical values (cf. Section 3.3) or an estimation based on the maximum error impact. In the last case, Pelzer [21] states that if the probability distribution is known, the standard deviation of an elementary error can be estimated with regard to its maximum error. Therefore, if a variable follows a rectangular distribution, the standard deviation is retrieved by multiplying the maximum error with 0.6. In case of a triangular distribution, multiplication is done by 0.4 and for normal distributions by 0.3. In what concerns the stochastic correlating group, values for the correlations must be supported by empirical values or literature. They represent stochastic relations for multi-dimensional normal distributed observations [17]. If the terms of Equation (3) are to be summed up, it can be seen that according to the matrices structures (see Equations (2) and (4)) the individual results are as follows: for non-correlating errors—a diagonal matrix, for functional and stochastic correlating errors—fully populated matrices. Thus, the synthetic variance-covariance matrix **Σ***ll* is also fully populated and illustrates the existing observation variances and covariances and indirectly their correlations.

#### **3. Elementary Error Theory for Terrestrial Laser Scanners**

#### *3.1. Error Soruces in Terrestrial Laser Scanning*

In order to apply the EEM on laser scanners, the error sources need to be identified and classified. As any measurement instrument, TLSs are realizations of an idealistic measurement system, therefore affected by physical manufacturing limitations. Even if the instrument itself would be hypothetically flawless, all measurements would be affected by the environment through which the electromagnetic waves are traveling (cf. [34]). Other error sources are related to the measured object properties such as surface material, roughness and color. These play an important role on the distance measurements and strongly depend on the used wavelength [35]. According to other authors (cf. [36]), scanning geometry is also considered an error source. Only instrumental and environmental error sources are treated in this contribution.

For a better understanding of how the TLS observations affect the coordinates (Figure 1), the mathematical relations between range (*R*), horizontal angle (λ), vertical angle (θ) and Cartesian coordinates (*X*, *Y*, *Z*) are described generically as follows:

$$\begin{aligned} X &= R \cdot \sin(\lambda) \cdot \cos(\theta), & R &= \sqrt{X^2 + Y^2 + Z^2}, \\ Y &= R \cdot \sin(\lambda) \cdot \sin(\theta), & \lambda &= \operatorname{atan}\left(\frac{X}{Y}\right), \\ Z &= R \cdot \cos(\theta), & \theta &= \arccos\left(\frac{Z}{R}\right). \end{aligned} \tag{5}$$

#### *3.2. Instrumental Elementary Errors*

In comparison to the panorama TLSs architecture, the hybrid scanner architecture is less present in commercially available TLSs. This may be a reason for the reduced amount of scientific publications on calibration models for hybrid scanners. Even though it measures basically the same type of polar coordinates, calibration parameters (CPs) are of a more complex nature (cf. [37]). The most common example is the rotating polygon mirror used for deflecting the laser beam. On one side, the distance varies at each mirror position and on the other side, the EDM source is usually mounted with an offset from the rotation axis, not to mention that it may be intentionally tilted. Further on, a classification of the instrumental errors is necessary. Firstly, an explanation is given about how the errors are considered and afterwards numerical values are given.

**Figure 1.** Example of Riegl Terrestrial Laser Scanner (TLS) and main axes.

To begin with, the non-correlating elementary errors are considered. These are measurement noise for angle and range measurements, which are not directly specified by the manufacturer. For range measurement, there is an entry for accuracy and one for precision. As defined by Riegl Laser Measurement Systems GmbH (Horn, Austria) [38], precision is the degree to which further measurements show the same result. If the definition of standard deviation is considered, it expresses how widely the random variable is spread out relative to the mean value of the sample [39]. Therefore, the given value for precision will be used as an indicator for instrument internal range noise at all measured ranges (see Table 1). For angle measurements, the data sheet of the instrument offers only "angle resolution" without further details. According to Wunderlich et al. [40], angular resolution can be interpreted as measurement precision (one sigma), therefore the same convention is used (see Table 1). The terms are generally presented in Equation (6) and their values are found in Table 1. Having this, the first term of Equation (2) and the first term of Equation (4) are now defined as follows:

$$
\Sigma\_{\delta\delta,k} = \begin{bmatrix}
\sigma\_A^2 & 0 & 0 \\
0 & \sigma\_\partial^2 & 0 \\
0 & 0 & \sigma\_R^2
\end{bmatrix}, D = I. \tag{6}
$$


**Table 1.** Classification of instrumental errors and their dimensions.

The influencing matrix *D* is the identity matrix, because no transformation from the coordinate space into observation space is needed at this moment. Only after the complete synthetic VCM is computed, a transformation based on Equation (5) is made from observation space to coordinate space.

Regarding the functional model of observations, a model defined by Lichti [41] and Lichti [42], later simplified by Schneider [43] is adopted. The latter applied it on a Riegl LMS-420i and could successfully improve the results after calibration. The simplification is mostly justified by the fact that not all of the CPs can be classified as significant after a calibration. Furthermore, if they are highly correlated they only reduce the validity of the model. Some of them are negligible, some are not determinable or separable, and therefore the used model is restricted to the minimum number of CPs identified as significant [37]. For more details about these parameters, the reader is advised to consult [37,43]. Following [43], the CPs for each observation can be defined as stated:

$$\begin{aligned} \Delta \mathcal{R} &= a\_0 + a\_1 \mathcal{R} + a\_2 \mathcal{R}^2, \\ \Delta \lambda &= b\_1 \sec(\theta) + b\_2 \tan(\theta) + b\_3 \sin(\lambda) + b\_4 \cos(\lambda) + \arcsin(\frac{b\_5}{\mathcal{R}}) + b\_6 \sin(2\lambda) + b\_7 \cos(2\lambda) + b\_8 \cos(3\lambda), \\ \Delta \theta &= c\_0 + c\_1 \sin(\theta) + c\_2 \cos(\theta) + \arcsin(c\_3/\mathcal{R}) + c\_4 \cos(3\lambda), \end{aligned} \tag{7}$$

where *a*<sup>0</sup> is the zero point error, *a*<sup>1</sup> scale error, *a*<sup>2</sup> quadric scale error, *b*<sup>1</sup> collimation axis error, *b*<sup>2</sup> horizontal axis error, *b*<sup>3</sup> and *b*<sup>4</sup> first and second horizontal circle eccentricity, *b*<sup>5</sup> eccentricity of the collimation axis with respect to the vertical axis, *b*<sup>6</sup> and *b*<sup>7</sup> non-orthogonality of the plane containing the horizontal angle encoder and the vertical axis, *b*<sup>8</sup> empirical parameter for compensation of remaining systematic effects, *c*<sup>0</sup> vertical circle index error, *c*<sup>1</sup> and *c*<sup>2</sup> first and second vertical circle eccentricity, *c*<sup>3</sup> eccentricity of the collimation axis with respect to the trunnion axis, *c*<sup>4</sup> empirical parameter to model a sinusoidal errors function of the horizontal direction with period of 120◦ (cosine term).

Out of all CPs, only some of them have numerical values and were determined as significant after calibration by Schneider [43]. For the EEM, variances of the CPs are introduced in the middle term of Equation (4) and the *F* matrix contains the partial derivatives of Equation (7). The values for the variances are presented in Table 1 with adopted dimensions for the Riegl VZ-2000 scanner. Further investigations on the hybrid scanner architecture are in progress based on the foundations set in [44].

#### *3.3. Meteorological Elementary Errors*

#### 3.3.1. Influences on the Distance Measurement

Similar to EDM of total stations, distance measurements in TLS, are influenced by air temperature and air pressure; in any case for long distances. Partial water vapor pressure is intentionally neglected due to its small influence. Most TLSs use near-infrared light for measuring distances. As known, the speed of light traveling through the atmosphere's different layers is diminished in comparison to the speed of light in vacuum. The atmospheric correction increases proportionally with the measured distance [34]. In case of ranges up to 200 m, these corrections may be neglected, but it cannot be neglected for long range scanners (e.g., Riegl VZ-2000) that measures up to 2050 m. According to manufacturer's specifications, the Riegl scanner has an atmospheric correction model implemented in the instrument, meaning that distances are corrected based on the introduced parameters for temperature, pressure and relative humidity. Information of how this happens can be taken from the RiSCAN Pro software documentation [45] and further inspected in the IAG 1999 resolutions [46]. The authors retain from explaining the whole process of retrieving the influencing coefficients for distance measurement and directly give the formula implemented in the EEM:

$$
\Delta n \cdot 10^{-6} = -0.93 \cdot \Delta t + 0.27 \cdot \Delta p\_\prime \\
\Delta R = -R \cdot \Delta n\_\prime \tag{8}
$$

where Δ*n* is the the change of the group refractive index of light, Δ*t* is change in temperature (◦C) and Δ*p* is change in pressure (hPa). Finally, the change in range Δ*R* is given. Note that these parameters are calculated for a mean atmosphere of 17 ◦C, 1000 hPa pressure and a wavelength of λ = 1550 nm. Interpreting this in terms of parts per million (ppm) depending on the two atmospheric parameters in standard conditions, a change in t of 1 ◦C affects the distance and refractive index by 0.93 ppm, a change in air pressure of 10 hPa yields a −2.7 ppm correction on the distance. For more details about this topic, the reader can consult [47] or [34].

#### 3.3.2. Influences on the Vertical Angle Measurement

In addition to the effects on distance measurements of any electro-optical measurement, atmospheric refraction also influences the vertical angle measurements. This effect causes image scintillation, often obvious in its extreme case when temperature gradients near the ground are high (e.g., in the desert or on a highway in hot summer days). This mostly affects angle measurements and is likewise important in geodesy, receiving much attention in transferring heights by trigonometric leveling. Nevertheless, this effect also occurs in TLS measurements and has been empirically studied by Friedli et al. [48]. The reader is advised to consult this work for understanding how refraction angles can be determined with the aid of reference values from total station measurements.

Figure 2 denotes the effects of atmospheric refraction out of which the refraction angle correction δ/2 is of further interest. This angle is given between the expected wave path and apparent line of sight also called tangent to refracted wave path. For more details about how δ/2 is deduced, refer to [49]. There are different ways of expressing the refraction angle correction, but only one has been chosen based on its simplicity and implemented in the EEM. The choice is not relevant in case of stochastic modeling. The corrected vertical angle can be therefore computed [49]:

$$
\theta = \theta' + \frac{\delta}{2}, \frac{\delta}{2} = \frac{R}{2 \cdot E\_R} \cdot k \cdot \rho,\tag{9}
$$

where θ is the corrected vertical angle, θ the measured vertical angle, *R* measured range, *ER* Earth's middle radius (6381 km), *k* refraction coefficient and ρ conversion constant between angle measurement units (degrees or grads) and radians. The coefficient of refraction *k* is usually needed to account for the curved light path from one point to another. It is defined as a ratio between the Earth radius and the radius of the line of sight which is mostly convex [50]. Very often, the Gaussian value of *k* = +0.13 is used by default as a setting for total station measurements, hoping that it holds true for most applications [51]. Nevertheless, *k* strongly varies throughout the day and is directly dependent on the temperature gradient ∂*T*/∂*Z* (K/m) (cf. [52]). If the refraction coefficient of a particular point is of interest, the local refraction coefficient *kloc* is given as a function dependent on temperature, pressure and the local temperature gradient (cf. [49,52]):

$$k\_{\rm loc} = 503 \cdot \frac{p}{T^2} \cdot \left(0.0343 + \frac{\partial T}{\partial Z}\right) \tag{10}$$

where *p* is pressure (hPa), *T* is temperature in (K) and ∂*T*/∂*Z* (K/m) is the temperature gradient at a certain point. The term *kloc* is used instead of an average *k* in Equation (9) for further purposes. As noticed in Equation (10), temperature gradient strongly determines the size of the local refraction coefficient; hence, its variation from ground level up to 100 m above the ground, as relevant for the later given examples, will be discussed. This is defined in meteorology or climate research under the name of micro- and local climate [53]. Hirt et al. [54] use the terms higher, intermediate and lower atmosphere to define the variation of the vertical temperature gradient (VTG) within a given range. By higher atmosphere, the layers from 100 m and above the ground surface are addressed. The VTG in this part of the troposphere has values around −0.006 K/m and is fairly independent of the Earth's surface temperature [54]. The next layer, the intermediate atmosphere between 20–30 m and 100 m is weakly influenced by the ground temperature and has an average value for the VTG of −0.01 K/m. This is where the refraction coefficient has an average value of +0.15 and it is also the layer to which the Gaussian value is most appropriate. Going a level lower, the first layer, considered the lower

atmosphere is where ground temperature reaches its maximum influence on the VTG. Several studies, summarized in [54], showed variations of the refraction coefficient between −3.5 and 3.5. Noteworthy are the empirical findings of Hennes [55] in which the local refraction coefficient reaches values of −2.9 (from a VTG of −0.5 K/m) leading to a concave curvature of the light path, contrasting the common belief about the chord being convex in almost all cases. Nevertheless, a less drastic value of −0.2 K/m is used in the current study, resembling an average value for this layer.

**Figure 2.** The effects of atmospheric refraction (after [3]).

Similar as in Section 3.3.1, the influencing coefficients are determined after computing the partial derivatives of Equations (9) and (10). Therefore, numeric values have been exemplary computed with the same conditions as stated before (*t* = 17 ◦C, *p* = 1000 hPa, VTG = −0.01 K/m) at a distance of 1000 m. The change in the measured vertical angle (in radians) is given by:

$$
\Delta\theta \cdot 10^{-6} = -0.08 \cdot \Delta t + 0.01 \cdot \Delta p + 468.17 \cdot \Delta VGT,\tag{11}
$$

where Δθ is the the change in the measured vertical angle, Δ*t* is change in temperature (◦C), Δ*p* is change in pressure (hPa) and Δ*VGT* is the change in vertical temperature gradient. In other words, a change in temperature of 10 ◦C affects the vertical angle by −0.8 μrad (−0.05 mgon), a change in air pressure of 10 hPa affects the vertical angle by 0.1 μrad (+0.006 mgon) and the most significant factor, a change in the VGT of 1 K/m results in a change of the angle with 468.17 μrad (29.8 mgon). This is not to be confused with the systematic effect of the refraction angle correction δ/2. For comprehension, for the above stated conditions and at 1000 m, δ/2 has a value of 0.7 mgon, which leads to a value of the linear error e = 11.4 mm. The intention is not to correct these systematic effects, but to show how varying temperature and pressure influence the error of position.

Although often not considered, air pressure also follows a gradient. This is less variable than the VGT and according to the Deutscher Wetterdienst Lexikon [56], the pressure gradient throughout the mentioned atmospheric layers is δ*p*/δ*z* = −0.125 hPa. This information will also be further used.

All these layer definitions and given values are adopted further in Section 3.4 to derive the necessary variances and covariances needed in the EEM.

#### *3.4. Atmospheric Errors as Stochastic Correlating Errors*

In most terrestrial precision measurements, if the atmospheric parameters temperature, pressure and relative humidity are needed, they are measured at the station point and in some cases near the observed object or second station point. For corrections, an averaged value of these parameters is used in most cases. This may hold true for airborne laser scanning, where the average between aircraft and ground temperature is sufficient [57], but in TLS, the situation changes. In addition to this, it was shown (cf. [48,54]) that even within a short time span, these may present strong variations and there is no straightforward method of correcting the measurement values dependent on these temporal variations. For this reason, modeling the impacts on the observations needs to be done stochastically. To do so, a VCM of the varying terms air temperature (*t*), air pressure (*p*) and VGT (*g*) is needed. The challenge is to fully populate the VCM **Σ**γγ,*<sup>h</sup>* from Equation (4) so that existing correlations between all elementary errors are known. This will have the general form in case of *t*, *p* and *g* (only upper diagonal presented) as follows:

**Σ**γγ,*atm* = ⎡ ⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ *s*2 *<sup>t</sup>*<sup>1</sup> *st*<sup>1</sup>*p*<sup>1</sup> *st*<sup>1</sup>*g*<sup>1</sup> *st*<sup>1</sup>*t*<sup>2</sup> *st*<sup>1</sup>*p*<sup>2</sup> *st*<sup>1</sup>*p*<sup>2</sup> ··· *st*<sup>1</sup>*tm st*<sup>1</sup>*pm st*<sup>1</sup>*gm s*2 *<sup>p</sup>*<sup>1</sup> *sp*<sup>1</sup>*g*<sup>1</sup> *sp*<sup>1</sup>*t*<sup>2</sup> *sp*<sup>1</sup>*p*<sup>2</sup> *sp*<sup>1</sup>*g*<sup>2</sup> ··· *sp*<sup>1</sup>*tm sp*<sup>1</sup>*pm sp*<sup>1</sup>*gm s*2 *<sup>g</sup>*<sup>1</sup> *sg*<sup>1</sup>*t*<sup>2</sup> *sg*<sup>1</sup>*p*<sup>2</sup> *sg*<sup>1</sup>*g*<sup>2</sup> ··· *sg*<sup>1</sup>*tm sg*<sup>1</sup>*pm sg*<sup>1</sup>*gm s*2 *<sup>t</sup>*<sup>2</sup> *st*<sup>2</sup>*p*<sup>2</sup> *st*<sup>2</sup>*p*<sup>2</sup> ··· *st*<sup>2</sup>*tm st*<sup>2</sup>*pm st*<sup>2</sup>*gm s*2 *<sup>p</sup>*<sup>2</sup> *sp*<sup>2</sup>*g*<sup>2</sup> ··· . . . . . . . . . *s*2 *<sup>g</sup>*<sup>2</sup> ··· ... . . . . . . . . . *s*2 *tm stmpm stmgm s*2 *tm spmgm s*2 *gm* ⎤ ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ (12)

The main diagonal is not difficult to fill in, according to what will be explained next, but the rest of the elements on all the upper and implicitly lower part of the same matrix are actually the challenge. To overcome this, correlation coefficients between *t*, *p* and *g* are computed for the given spatial distribution of all observations. This will be further on explained.

As in any terrestrial electro-optical measurement, in TLS observations light travels from the instrument to the measured object and back. Due to varying atmospheric conditions, it will be perturbed throughout the whole path and in order to evaluate how air temperature, air pressure and VTG vary along the path, pre-knowledge about these parameters (cf. Section 3.3.2) is used in a combined way with spatial information. Suppose a laser scanner is stationed at a certain distance near a tall object and observations are possible from the base to the tip of that object. If a rough digital terrain model (DTM) of the area is available, then the local topography is known, which further allows a classification of the VGT depending on how the topography varies. Simply explained, the limits of the gradient layers can be defined as surfaces with an offset from ground level according to how meteorologists have defined these limits (Figure 3 left). The yellow surface defines the separating layer at about 25 m between the lower and intermediate atmosphere; the red layer is the separation between intermediate and higher atmosphere at about 100 m above the ground. In order to have a better overview of the further steps, a vertical section is selected and exemplified in Figure 3 right.

**Figure 3.** (**Left**) Spatial separation of vertical temperature gradient (VGT) layers; (**right**) one section of the spatial VGT model.

It is necessary to roughly know the position of the scanner on the DTM. This is often referred to as georeferencing, but in this case the accuracy of the scanner's position is not of high importance, therefore an approximation suffices. In most cases, the air temperature and air pressure is measured near the laser scanner, usually at the instrument height. According to the situation depicted above, this is true for temperature and pressure only near the laser scanner, but the interest is in gaining information about the atmospheric parameters along the whole measurement path. Therefore, in the next step the observation lines are reconstructed in relation to the scanner position on the DTM. This directly shows which observation line passes through which atmospheric layer. Only two of them are depicted by the blue lines in Figure 3 right, but the same principle applies for all the rest. Further on, a series of points along these lines are selected; denoted by the yellow and pink circles. Point spacing along the observation line may be done subjectively; but a uniform distribution between scanner and object is suggested for representative results. Each of these points receive values for *t*, *p* and *g*, determined according to their position in space and in relation to the measured atmospheric parameters at the station point denoted by "TLS" in Figure 3 right. For example, along the lowest observation line, each yellow point will receive a value starting from *t*11, *p*11, *g*<sup>11</sup> to *t*1*n*, *p*1*n*, *g*1*n*. The individual values are extrapolated according to the individual position. If more measurements of the atmospheric parameters at other positions within the DTM are available, they can be considered within the extrapolation processes. This applies for all the other observation lines, until the highest one is reached. In this case, the pink points (Figure 3 right) receive *tm*1, *pm*1, *gm*<sup>1</sup> to *tmn*, *pmn*, *gmn*. Notice that the points are chosen on the same vertical line, fact that will be explained later. Having series of values for all three parameters along all observation lines allows the computation of variances and covariances along and between each observation line according to:

$$\begin{aligned} s\_a^2 &= \frac{1}{n-1} \cdot \sum\_{i=1}^n \left( a\_i - \overline{a\_i} \right)^2, \\ s\_{\overline{a\_i}, a\_j} &= \frac{1}{n-1} \cdot \sum\_{i,j=1}^n \left( a\_i - \overline{a\_i} \right) \cdot \left( a\_j - \overline{a\_j} \right), \end{aligned} \tag{13}$$

where *s*<sup>2</sup> *<sup>a</sup>* is the empirical standard deviation computed for each of the three atmospheric parameters. The value *ai* and *aj* are replaced consecutively by *ti*, *pi*, *gi* and *n* is the number of points along each observation line, *sai*,*aj* is the empirical covariance between pairs of the three parameters. To exemplify this, the general VCM **Σ**γγ,*atm* for one observation line has the following form:

$$
\Sigma\_{\mathsf{Y}\mathsf{Y}\mathsf{Z}\mathsf{atm}} = \begin{bmatrix}
\mathbf{s}\_{t}^{2} & \mathbf{s}\_{t\mathsf{p}} & \mathbf{s}\_{t\mathsf{\mathcal{K}}} \\
\mathbf{s}\_{t\mathsf{p}} & \mathbf{s}\_{p}^{2} & \mathbf{s}\_{p\mathsf{\mathcal{K}}} \\
\mathbf{s}\_{t\mathsf{\mathcal{K}}} & \mathbf{s}\_{p\mathsf{\mathcal{K}}} & \mathbf{s}\_{\mathsf{\mathcal{K}}}^{2}
\end{bmatrix}.
\tag{14}$$

Out of this VCM, a correlation matrix is computed, out of which correlation coefficients are extracted. For example, along each line, the correlation matrix *Ratm* is obtained and contains the following correlation coefficients *rtp*,*rtg*,*rpg*. This is valid for all lines up to the *n*-th observation.

$$\mathcal{R}\_{\text{atm}} = \frac{1}{\sqrt{\text{diag}\{\Sigma\_{\text{YY},\text{atm}}\}}} \cdot \Sigma\_{\text{YY},\text{atm}} \cdot \frac{1}{\sqrt{\text{diag}\{\Sigma\_{\text{YY},\text{atm}}\}}} = \begin{bmatrix} 1 & r\_{tp} & r\_{t\mathcal{S}} \\ r\_{tp} & 1 & r\_{p\mathcal{S}} \\ r\_{t\mathcal{g}} & r\_{p\mathcal{S}} & 1 \end{bmatrix}. \tag{15}$$

This is valid along observation lines and helps at filling in the block matrices on the main diagonal of VCM from Equation (12) with submatrices like in Equation (14). For all other elements of **Σ**γγ,*atm* the covariances are computed with the help of the correlation coefficients according to:

$$s\_{i\bar{j}} = r\_{i\bar{j}} \cdot s\_{i\bar{"\bar{s}}\text{'}\bar{s}}.\tag{16}$$

The values for *si* and *sj* are computed along each observation line with the help of Equation (13). Following Equation (15), a set of values for *rtp*,*rtg*,*rpg* is obtained. In addition to these, the correlation coefficients of the same parameters (*t*–*t*, *p*–*p* and *g*–*g*) between the observation lines have to be determined *rtt*,*rpp*,*rgg* and to accomplish this, VCMs and correlation matrices are computed between each parameter of the observation lines (e.g., *t*<sup>11</sup> and *tm*1). This is also the reason why it is relevant to have the yellow and pink points on the same vertical line. In this way values for each of the atmospheric parameters are treated as series of values for the same dimension, in this case vertical direction. A drawback of this proposal is that a number of *n* TLS observations leads to a number of *n*!/(*n* − 3)! permutations of correlation coefficients. This means that for e.g., 200 values taken three times (*t*, *p*, *g*) one would obtain a number of 7,880,400 correlation coefficients that need to be properly arranged in **Σ**γγ,*atm*. This is currently not achievable due to technical reasons for TLS observations where the number of observations easily reaches a few million. Therefore, one generic value is taken for each of the coefficients *rtp*,*rtg*,*rpg* and *rtt*,*rpp*,*rgg*. Numeric values are computed between the lowest and highest observation line taken from the vertical section as denoted in Figure 3 right. Finally, the individual values for the covariances are computed based on Equation (16) and then introduced in **Σ**γγ,*atm* like in Equation (12). Returning to the EEM, now that the matrix **Σ**γγ,*atm* is available, the influencing matrix *Gh* from Equation (2) must be properly filled. The complete matrix is a block matrix that has the partial derivatives of the observations with respect to *t*, *p* and *g*, as presented in Equation (17).

$$\mathbf{G}\_{\text{atm}} = \begin{bmatrix} G\_1 & 0 & \cdots & 0 \\ 0 & G\_2 & 0 & \vdots \\ \vdots & 0 & \ddots & \vdots \\ 0 & \cdots & \cdots & G\_n \end{bmatrix}, \mathbf{G}\_l = \begin{bmatrix} \frac{\partial \lambda}{\partial t} & \frac{\partial \lambda}{\partial p} & \frac{\partial \lambda}{\partial g} \\ \frac{\partial \Omega}{\partial t} & \frac{\partial \Omega}{\partial p} & \frac{\partial \Omega}{\partial g} \\ \frac{\partial \mathbf{R}}{\partial t} & \frac{\partial \mathbf{R}}{\partial p} & \frac{\partial \mathbf{R}}{\partial g} \end{bmatrix}, \tag{17}$$

where *i* = 1 ... *n* and *G***<sup>1</sup>** to *Gn* are block matrices that have the partial derivatives as shown above. Effects on the horizontal angles have not been discussed and are not considered in this model. Therefore, the first line of each block *g* will be filled with 0. The second line includes the coefficients presented in Section 3.3.2 in Equation (11) and this is the only line that has influencing values for all variations in the fully populated VCM (Equation (12)); the last line of *g* has the influencing values presented in Section 3.3.1, Equation (8) with the last element 0. The numerical values must be computed with regard to given atmospheric conditions and at a given range for each situation. Finally, the last term of the synthetic VCM can be computed and therefore the influences of instrumental and atmospheric parameters can be combined. In the upcoming section, the whole methodology presented until now will be applied for TLS point clouds of an arch dam.

#### **4. Study Case: Arch Dam Kops**

The Kops water dam is a storage concrete dam built between 1962 and 1969 in Vorarlberg, Austria. It is considered a hybrid dam made out of a gravity dam and an arch dam with artificial counterfort or abutment. It retains a volume of almost 43 million m<sup>3</sup> of water, thus creating the 1 km<sup>2</sup> "Kopssee" lake [58]. Only measurements of the downstream (airside) arch dam are considered. For this reason, its dimensions are mentioned to give a general impression. The crown spans over 400 m, its height is 122 m from foundation to crest and has a crest width of 6 m. Between 29 July and 2 August 2019, a first measurement campaign of the Kops dam took place and part of the results from another type of laser scanner are presented by Kerekes and Schwieger [59]. Further on, the EEM is applied on point clouds acquired with the Riegl VZ-2000 from varying positions.

In order to apply the EEM for meteorological elementary errors as described in Section 3.4, a DTM for the area of interest was cordially made available by the "Landesamt für Vermessung und Geoinformation" Land Vorarlberg, Austria. TLS Point clouds were acquired from four different station points. Figure 4 shows the distribution of these on the DTM together with an example of a vertical

section plane out of which temperature and pressure are extracted. Results will be presented for all four station points (S1–S4).

**Figure 4.** TLS station points and vertical section on digital terrain model (DTM) (Source for DTM: LiDAR Data obtained from Land Vorarlberg, Austria—data.vorarlberg.gv.at).

In order to have an overview about the varying scanning configurations and atmospheric conditions, Table 2 summarizes all the relevant parameters.


**Table 2.** Overview of scanning configurations and atmospheric conditions.

Considering the harsh local topography with steep slopes and vegetation, only the four station points depicted in Figure 4 were measurable within reasonable time, effort and coverages of the dam airside. With exception of S4, all other point clouds cover more than 80% of the airside surface. Weather recordings were made at each station point at the instrument height approx. 2 m above ground height with the Greisinger precision Thermo-, Barometer GTD1100 (Greisinger GmbH, Regenstauf, Germany). According to the technical specifications, air temperature is measured with an accuracy of +/−1% of the reading in the interval −10 ◦C to +50 ◦C and air pressure with +/−1.5 hPa in the interval of 750 hPa to 1100 hPa [60]. These accuracies are considered in the process of determining *t* and *p*. As regards the VGT, an empirical value for the uncertainty can be found in [55] for an Alpine region:

$$
\sigma\_{VGT} = \sqrt{\frac{1}{12} \cdot \left( VGT\_{\text{max}} - VGT\_{\text{min}} \right)^2} \tag{18}
$$

where a value for σ*VGT* = 0.25 K/m in case of the Alpine region in the lower atmosphere is given and *VGTmax* and *VGTmin* are the upper and lower numerical recorded values for the VGT [55].

In case of the other two layers, no empirical values for variances have been found to the best of the author's knowledge. Due to this, values are obtained by multiplying the VGT value with 0.3 following the explanation in Section 2.1; these variations are likewise considered in determining the values for *t*.

For all station points the same methodology is applied, but the vertical section is only visualized in the case of S4 (Figure 5). The authors consider this case the most interesting since the distance to the dam is the longest and observations pass through different atmospheric layers more than once.

**Figure 5.** Vertical section through the terrain and temperature gradient fictive separation.

Almost half of the observation lines in this profile pass through the lower layer twice, meaning that variances of temperature and VGT affect the lowest points in the point cloud more than the ones obtained from observations that travel through a more stable atmospheric layer. This is confirmed further when analyzing the error of position in the point cloud. To have an overview of all variances and spatial correlation coefficients, Table 3 presents them for all station points.

**Table 3.** Overview of variances and spatial correlation coefficients for atmospheric elementary errors.


The last pairs of correlation coefficients *rpt*,*rgt*,*rgp* are not in the table because they have the same value with *rtp*,*rtg*,*rpg*. Note that these values resemble the spatial correlations only. The subject of temporal correlations will be addressed in a future publication after a second measurement epoch is available. The variances and correlation coefficients from Table 3 are used to finally create a VCM **Σ**γγ,*atm* for the atmospheric elementary errors for each station point. In case of the instrumental errors, all the values are the same as stated in Table 1 since the same instrument was used in all station points.

The EEM is implemented in Matlab—MathWorks and currently limited to handling VCMs having sizes of up to 21,000 × 21,000 cells. More details about this can be found in [59]. Before applying the EEM on the point clouds, only points on the dam are selected and a subsampling is done. Consequently, the complete point cloud contains points on the dam airside with a distance of 1.5 m between them. Additional to this, vertical sections on the dam are analyzed since much attention was accorded to how atmospheric parameters vary along vertical profiles. The point spacing in the section is denser with an average distance of 15 cm between the points. This is done due to technical restriction mentioned above.

Coordinates (*X*, *Y*, *Z*) are considered instead of observations (*R*, λ, θ), because the VCM will be used for estimating the geometric primitives of a B-Spline surface in the future [59]. Therefore, the synthetic VCM is computed in observation space and then transformed with the help of equations 5 to coordinate space. Results are presented with regard to the error of position, spatial correlations along a vertical line chosen to be as long as possible along the scan (Figure 6) and contribution of all variances to the error budget in case of a single point.

**Figure 6.** Example of error of position distribution on point clouds, selected vertical sections (red lines) and selected points (red circles). Scanning station points are (**a**) S1, (**b**) S2, (**c**) S3, (**d**) S4.

#### *4.1. Error of Position*

The main diagonal of the synthetic VCM in coordinate space contains the variances for each point in the point cloud. The error of position is computed according to [21]:

$$
\sigma\_{xyz\text{\tiny\text{\tiny\text{\tiny\text{\}}}} = \sqrt{\sigma\_{xi}^{2} + \sigma\_{yi}^{2} + \sigma\_{zi}^{2}} \tag{19}
$$

Points are visualized with respect to their position in space (local coordinate system) relative to the laser scanner position (0, 0, 0) and the dimension of the error of position is given on a color scale. Figure 6 exemplary shows the results obtained from all four station points, the selected vertical section and a point for which the percentage contribution of the variances is given later on.

Here the coverage of the air side can be seen for the first time. As expected from Figure 4, S1 to S3 cover most of the dam's airside, whilst S4 (Figure 6d) has the smallest coverage due to height difference and vegetation. Outcomes led to average errors of position, as follows: S1 − σ*xyz* = 15.2 mm, S2 − σ*xyz* = 9.6 mm, S3 − σ*xyz* = 16.4 mm and S4 − σ*xyz* = 36.5 mm. A common color scale was chosen to maintain comparability; this is why the reader is asked to consult the digital version of the paper. At a first glance, it can be seen in which way the errors of position are distance dependent, with the smallest values for S2 (Figure 6b) which is the nearest point and the biggest value for S4 (Figure 6d) at a distance of 466 m. Both S1 and S3 (Figure 6a,c) present similar results due to the similar measurement configuration. At a closer inspection of the scan from S4, it can be seen how the error of position decreases with height, reaching a minimum at the crest (dark yellow). The lowest part has the highest values for the error of position (bright yellow). This resembles the smaller variation of the VGT (see Table 3) for TLS observations that pass through more stable atmospheric layers.

#### *4.2. Spatial Correlations Along Vertical Sections*

Existing correlations can be analyzed after computing a correlation matrix based on the VCM (Equation (15)). Generally, high correlations are an indication for high variances. Previous publications of the EEM for TLS topic (cf. [17]) treated atmospheric elementary errors as functional correlating. In this contribution, atmospheric elementary errors were treated stochastic correlating (cf. Section 3.4) for the first time. This leads to different spatial correlations than have not been discussed before. For this reason, special attention is offered to the stochastic correlating errors and presented in parallel with the ones obtained from the complete VCM where instrumental errors also influence the results. For each station point, one vertical section has been selected and spatial correlations are presented between the lowest point of the section and all other. The analysis is made only for the height coordinates Z. Results do not change if observations would be analyzed instead. The reason for choosing only one vertical section is that emphasis is put on how the stochastic correlating errors influence the correlations and height error of position in comparison with the complete VCM of the same station point. These cases are shown in parallel. In the first case, the EEM considers instrumental and atmospheric errors (Figure 6 left side) and in the second only atmospheric errors are considered in the EEM (Figure 6 right side).

Analyzing all correlations when instrumental and atmospheric elementary errors are considered, values are in almost all cases higher than 0.5 (Figure 7a,c,d) and present a linear decrease with increasing height. The same effect is noticed with the standard deviations of the heights that remain at approximately the same level. The exception to this is S2 (Figure 7b) where the correlations decrease with height and standard deviation increases. This is the station point with the smallest distance to the dam; therefore, an explanation for this effect is analyzed in the next section. Correlations in case of the atmospheric errors all show a linear behavior (Figure 7e–h), but at different levels. As maybe presumed, the standard deviations are very small at these distances and given level of variation for the atmospheric parameters, but the most interesting finding is at S4 (Figure 7h) where the decrease in standard deviation of height is obvious, explaining the presumption made in Section 4.1 that the upper observations travel through more stable layers of the atmosphere and are less affected by variations. This may lead to the thought that scans for e.g., objects over a valley are more reliable than those acquired from parallel to the ground level. This is partly true, since VGT may also be stable for a short period of time at ground level; therefore, this issue strongly depends on the topography and local conditions. If the conditions in Table 2 are reviewed, it can be seen that similar atmospheric conditions do not necessarily lead to a similar level of correlation, for example point clouds from S2 and S4 where acquired under differing conditions, but led to a similar level of correlation for the atmospheric elementary errors.

**Figure 7.** Left: spatial correlations, standard deviations (height)—complete synthetic VCM (**a**) S1, (**b**) S2, (**c**) S3, (**d**) S4; right: spatial correlations, standard deviations (height)—synthetic VCM with atmospheric elementary errors only (**e**) S1, (**f**) S2, (**g**) S3, (**h**) S4.

#### *4.3. Contribution of the Elementary Errors Variances to the Overall Error Budget*

The points depicted in Figure 6 (red circles) are chosen to analyze the contribution of variances to the whole error budget. Note that here, all dimensions (*X*, *Y*, *Z*) are considered and not only Z as in the section before. For the first three station points, they are approximately at the same level. In case of S4, this is not possible since that part of the dam cannot be scanned.

In all situations, instrumental errors make up the majority of the error budget (Figure 8). The parameters *b*4—horizontal circle eccentricity, *a*1—scale error and *c*1—horizontal circle eccentricity comprise in all cases more than 50% of the error budget. It can also be seen that some instrumental errors are negligible (e.g., *a*0) since they remain under the 1% quote. In the previous section, the standard deviations for Z are affected by instrumental errors of the vertical angle encoder (*c*0, *c*1, *c*4) and vertical angle noise σθ. It is also seen how this phenomenon is almost independent of distance and scanning configuration. Just to outline some of the instrumental errors, σθ is always between 9% and 11%, *c*<sup>1</sup> between 9% and 13% and *c*<sup>0</sup> 4% with exception to S4.

**Figure 8.** The contribution of variances for instrumental and atmospheric elementary errors in percent for (**a**) S1, (**b**) S2, (**c**) S3, (**d**) S4. For parameter names, see Sections 3.2 and 3.4.

The elementary errors for air temperature and air pressure at this level of variance (see Table 3) fall into the same category remaining under the quote of 1%, fact related to the small influencing coefficients in Equation (11). The contribution of the VGT reaches 3% of the error budget for S1 (Figure 8a) and S3 (Figure 8c). This is explainable by the fact that a large part of the observations are in the lower atmosphere layer, where the VGT instability is known to be higher (cf. [53]). A higher level of variance would generally lead to higher contributions to the error budget, but this is planned to be studied further for longer ranges and within different timespans.

#### **5. Conclusions**

Throughout this contribution, a method of defining a stochastic model for TLS observations was presented. At the time being, this model considers instrumental errors of long-range laser scanners and meteorological error sources.

The EEM is improved by directly introducing the existing spatial correlations into the synthetic VCM, without the need to compute them each time. This will confirm the EEM advantages when a second measurement epoch is available. It was also shown that some of the instrumental and atmospheric elementary errors can be neglected at the given level of variance, having a contribution of less than 1% to the complete error budget. Within this contribution, the line of previous EEM publications was continued with the introduction of atmospheric elementary errors treated as stochastic correlating and the introduction of a functional model for hybrid laser scanners. The newly achieved VCM plays an important role in surface estimations as already shown in [7]. Other possible applications that may benefit from this model can be encountered in landslide, glacier or rock cliff monitoring.

Resuming the newly discussed topics, it can be mentioned that:


As regards the topics under research, the EEM for TLS measurements is still lacking object related elementary errors. This implies using existing studies for different materials scanned from different positions at different ranges and include these errors into the stochastic correlating group. The authors intend to use intensity values of the reflected laser beam, as introduced by Wujanz et al. [61]. Having done this would make the EEM a powerful tool for generating a TLS stochastic model that considers all important error sources. Another topic under research is the sensitivity analysis of the input parameters. After completing the stochastic model, attention will be accorded to the impact of each individual input parameter on the outcomes. This allows an estimation of the optimal scanning configuration.

**Author Contributions:** Methodology, G.K. and V.S.; software, G.K.; validation, V.S. and G.K.; formal analysis, G.K.; resources, G.K. and V.S.; data acquisition, G.K.; writing—original draft preparation, G.K.; writing—review and editing, V.S.; supervision, V.S.; project administration, V.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by DFG (German Research Foundation), SCHW 838/7-3 within the project IMKAD II "Integrated space-time modeling based on correlated measurements for the determination of survey configurations and the description of deformation processes". The project is also associated with Cluster of Excellence Integrative Computational Design and Construction for Architecture (IntCDC) and supported by the DFG under Germany's Excellence Strategy—EXC 2120/1—390831618.

**Acknowledgments:** The authors also express their gratitude to Illwerke vkv AG, especially thanking Ralf Laufer for supporting the measurement campaign in 2019. Acknowledgment also goes to "Landesamt für Vermessung und Geoinformation" Land Vorarlberg for making LiDAR data available for the research area around the Kops lake.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
