**Preface to "Stochastic Models for Geodesy and Geoinformation Science"**

In geodesy and geoinformation science, as well as in many other technical disciplines, it is often not possible to directly determine the desired target quantities, for example, the 3D Cartesian coordinates of an object. Therefore, the unknown parameters must be linked with measured values, for instance, directions, angles, and distances, by a mathematical model. This consists of two fundamental components—the functional and the stochastic models. The functional model describes the geometrical–physical relationship between the measured values and the unknown parameters. This relationship is sufficiently well known for most applications.

With regard to stochastic models in geodesy and geoinformation science, two problem domains of fundamental importance arise:

1. How can stochastic models be set up as realistically as possible for the various geodetic observation methods and sensor systems, such as very-long-baseline interferometry (VLBI), global navigation satellite systems (GNSS), terrestrial laser scanners, and multisensor systems?

2. How can the stochastic information be adequately considered in appropriate least squares adjustment models?

Further questions include the interpretation of the stochastic properties of the computed target values for quality assessment of the results in terms of precision and reliability and for the detection of outliers in the input data (measurements).

In this Special Issue, current research results on these general questions are presented in ten peer-reviewed articles. The basic findings can be applied to all technical scientific fields where measurements are used for the determination of parameters to describe geometric or physical phenomena.

> **Frank Neitzel** *Editor*

## *Article* **Total Least-Squares Collocation: An Optimal Estimation Technique for the EIV-Model with Prior Information**

#### **Burkhard Scha**ff**rin**

Geodetic Science Program, School of Earth Sciences, The Ohio State University, Columbus, OH 43210, USA; schaffrin.1@osu.edu

Received: 3 May 2020; Accepted: 1 June 2020; Published: 13 June 2020

**Abstract:** In regression analysis, oftentimes a linear (or linearized) Gauss-Markov Model (GMM) is used to describe the relationship between certain unknown parameters and measurements taken to learn about them. As soon as there are more than enough data collected to determine a unique solution for the parameters, an estimation technique needs to be applied such as 'Least-Squares adjustment', for instance, which turns out to be optimal under a wide range of criteria. In this context, the matrix connecting the parameters with the observations is considered fully known, and the parameter vector is considered fully unknown. This, however, is not always the reality. Therefore, two modifications of the GMM have been considered, in particular. First, 'stochastic prior information' (p. i.) was added on the parameters, thereby creating the – still linear – Random Effects Model (REM) where the optimal determination of the parameters (random effects) is based on 'Least Squares collocation', showing higher precision as long as the p. i. was adequate (Wallace test). Secondly, the coefficient matrix was allowed to contain observed elements, thus leading to the – now nonlinear – Errors-In-Variables (EIV) Model. If not using iterative linearization, the optimal estimates for the parameters would be obtained by 'Total Least Squares adjustment' and with generally lower, but perhaps more realistic precision. Here the two concepts are combined, thus leading to the (nonlinear) 'EIV-Model with p. i.', where an optimal estimation (resp. prediction) technique is developed under the name of 'Total Least-Squares collocation'. At this stage, however, the covariance matrix of the data matrix – in vector form – is still being assumed to show a Kronecker product structure.

**Keywords:** Errors-In-Variables Model; Total Least-Squares; prior information; collocation vs. adjustment

#### **1. Introduction**

Over the last 50 years or so, the (linearized) Gauss-Markov Model (GMM) as standard model for the estimation of parameters from collected observation [1,2] has been refined in a number of ways. Two of these will be considered in more detail, namely


After a compact review and comparison of several key formulas (parameter estimates, residuals, variance component estimate, etc.) within the above three models, a new model will be introduced that allows the strengthening of the parameters and the weakening of the coefficient matrix at the same time. The corresponding estimation technique will be called "TLS collocation" and follows essentially the outline that had first been presented by this author in June 2009 at the Intl. Workshop on Matrices and Statistics in Smolenice Castle (Slovakia); cf. Schaffrin [9].

Since then, further computational progress has been made, e.g. Snow and Schaffrin [10]; but several open questions remain that need to be addressed elsewhere. These include issues related to a rigorous error propagation of the "TLS collocation" results; progress may be achieved here along similar lines as in Snow and Schaffrin [11], but is beyond the scope of the present paper.

For a general overview of the participating models, the following diagram (Figure 1) may be helpful.

**Figure 1.** Model Diagram.

The most informative model defines the top position, and the least informative model is at the bottom. The new model can be formed at the intermediate level (like the GMM), but belongs to the "nonlinear world" where *nonlinear* normal equations need to be formed and subsequently solved by iteration.

#### **2. A Compact Review of the "Linear World"**

#### *2.1. The (linearized) Gauss-Markov Model (GMM)*

Let the Gauss-Markov Model be defined by

$$\text{If } y = \,\_nA\,\,\xi + e\,,\,\,\eta := \text{rk}\,A \le \min\,\{m, n\}\,,\,e \sim \left(0, \sigma\_0^2 I\_{\mathbb{N}}\right)\,,\tag{1}$$

possibly after linearization and homogeneization; cf. Koch [2], or Grafarend and Schaffrin [12] among many others.

Here, *y* denotes the *n* × 1 vector of (incremental) observations,


Now, it is well known that the Least-Squares Solution (LESS) is based on the principle

$$
\epsilon^T \epsilon = \min \quad \text{s.t.} \; \epsilon = y - A \; \xi \tag{2}
$$

which leads to the "normal equations"

$$N\,\,\hat{\xi} = \mathfrak{c} \,\,\, \text{for} \quad [N,\,\,\mathfrak{c}] := A^{\mathrm{T}}[A,\,\,\mathfrak{y}].\tag{3}$$

Depending on rk *N*, the rank of the matrix *N*, the LESS may turn out *uniquely* as

$$
\mathcal{L}\_{\text{LESS}} = N^{-1} \mathcal{c} \quad \text{iff} \quad \text{rk} \; N = \text{rk} \; A = q = m; \tag{4}
$$

or it may belong to a *solution (hyper)space* that can be characterized by certain generalized inverses of *N*, namely

$$\begin{aligned} \left| \mathcal{L}\_{\text{LESS}} \right| &\in \{ \mathbf{N}^- \mathbf{c} \middle| \mathbf{N} \, \mathbf{N}^- \mathbf{N} = \mathbf{N} \} = \\ &= \left\{ \mathbf{N}\_{\text{rs}}^- \mathbf{c} \middle| \mathbf{N} \, \mathbf{N}\_{\text{rs}}^- \mathbf{N} = \mathbf{N} \, \right. \left. \begin{aligned} \left( \mathbf{N}\_{\text{rs}}^- \mathbf{N} \right)^- \mathbf{N} = \mathbf{N}\_{\text{rs}}^- \left( \mathbf{N}\_{\text{rs}}^- \right)^\mathrm{T} \right\}; \end{aligned} \tag{5}$$

where *N*− *rs* denotes an arbitrary *reflexive symmetric g-inverse* of *N* (including the "pseudo-inverse" *N*+),

$$\text{diff} \quad \text{rk} \, N = \text{rk} \, A = q < m. \tag{6}$$

In the latter case (6), any LESS will be *biased* with

$$\text{bias}\left(\mathcal{E}\_{\text{LESS}}\right) = E\left\{\mathcal{E}\_{\text{LESS}} - \underline{\boldsymbol{\xi}}\right\} \in \left\{ -\left(\boldsymbol{I}\_{\text{m}} - \boldsymbol{N}^{-}\boldsymbol{N}\right) \,\,\xi\right\} = \left\{ -\left(\boldsymbol{I}\_{\text{m}} - \boldsymbol{N}\_{\text{rs}}^{-}\boldsymbol{N}\right) \,\,\xi\right\},\tag{7}$$

its dispersion matrix will be

$$D\left\{\hat{\xi}\_{\text{LESS}}\right\} = \sigma\_0^2 \mathbf{N}^- \mathbf{N} (\mathbf{N}^-)^\mathbf{T} \in \left\{\sigma\_0^2 \mathbf{N}\_{\text{rs}}^- \middle| \mathbf{N}^\cdot \mathbf{N}\_{\text{rs}}^- \mathbf{N} = \mathbf{N}\_r \ \mathbf{N}\_{\text{rs}}^- \mathbf{N} \ \mathbf{N}\_{\text{rs}}^- = \mathbf{N}\_{\text{rs}}^- = \left(\mathbf{N}\_{\text{rs}}^-\right)^\mathbf{T}\right\} \tag{8}$$

and its Mean Squared Error (MSE) matrix, therefore,

$$\mathbf{MSE}\left\{\xi\_{\rm LESS}\right\} \in \left\{\sigma\_0^2 \middle| N\_{rs}^- + (I\_m - N\_{rs}^- \mathcal{N}) \left(\xi \left.\sigma\_0^{-2} \xi^T\right\|\left(I\_m - N\_{rs}^- \mathcal{N}\right)^T\right) \right\}.\tag{9}$$

In contrast to the choice of LESS in case (6), the *residual vector* will be *unique*; for *any* ξˆ LESS:

$$\overline{\sigma}\_{\text{LES}} = y - A \, \underline{\xi}\_{\text{LES}} \quad \sim \quad \text{(0, } \sigma\_0^2 (I\_n - A \, N\_{rs}^{-} A^T) = D\{y\} - D \{A \, \underline{\xi}\_{\text{LES}}\}). \tag{10}$$

Hence, the *optimal variance component estimate* will also be *unique*:

$$\boldsymbol{\hat{\sigma}}\_{0}^{2} = \left(\boldsymbol{n} - \boldsymbol{q}\right)^{-1} \boldsymbol{\overline{\varepsilon}}\_{\text{LESS}}^{\mathrm{T}} \cdot \mathbb{E}\_{\text{LESS}} = \left(\boldsymbol{n} - \boldsymbol{q}\right)^{-1} \left(\boldsymbol{y}^{\mathrm{T}} \boldsymbol{y} - \boldsymbol{c}^{\mathrm{T}} \boldsymbol{\xi}\_{\text{LESS}}\right),\tag{11}$$

$$E\left\|\mathfrak{H}\_{0}^{2}\right\| = \sigma\_{0}^{2}\,\,\,\,D\left\|\mathfrak{H}\_{0}^{2}\right\| \approx 2\left(\sigma\_{0}^{2}\right)^{2}/(n-q)\,\,\,\text{under quasi-normality.}\tag{12}$$

For the corresponding formulas in the full-rank case (4), the reflexive symmetric g-inverse *N*− *rs* simply needs to be replaced by the regular inverse *N*<sup>−</sup>1, thereby showing that the LESS turns into an unbiased estimate of ξ while its **MSE**-matrix coincides with its dispersion matrix (accuracy precision).

#### *2.2. The Random E*ff*ects Model (REM)*

Now, additional prior information (p. i.) is introduced to *strengthen* the parameter vector within the GMM (1). Based on Harville's [13] notation of a *m* × 1 vector 0 <sup>∼</sup> of so-called "*pseudo-observations*", the new *m* × 1 vector

$$\mathbf{x} := \underline{\mathbf{x}} - \underline{\mathbf{0}}\_{\sim} \text{ of \text{\textquotedblleft random effects\textquotedblright} (unknown) is formed with} \tag{13}$$

$$E\{\mathbf{x}\} = \beta\_0 \text{ as given } m \times 1 \text{ vector of p. i. and }\tag{14}$$

$$D|\mathbf{x}\rangle = \sigma\_0^2 Q\_0 \text{ as } n \times n \text{ positive-(semi)definite dispersion matrix of the p.i.d.} \tag{15}$$

Consequently, the Random Effects Model (REM) can be stated as

$$\underset{\sim}{y} := y - A \cdot 0 = \underset{n \times m}{A} \text{ x} + \text{e} \text{ , } q := \text{rk} \, A \le \underset{m \times 1}{\text{min}} \{ m, n \} \text{ , } \not p\_0 := \text{x} + \text{e}\_{0\prime} \tag{16}$$

$$\left( \begin{array}{cc} \varepsilon \\ \varepsilon 0 \end{array} \right) \sim \left( \begin{bmatrix} 0 \\ 0 \end{bmatrix}, \sigma\_0^2 \begin{bmatrix} I\_n & 0 \\ 0 & Q\_0 \end{bmatrix} \right), \ Q\_0 \text{ symmetric and } \text{nnd (non-negative definite)}.\tag{17}$$

If *P*<sup>0</sup> := *Q*−<sup>1</sup> <sup>0</sup> exists, the estimation/prediction of *x* may now be based on the principle

$$\left(\mathbf{e}^{\mathrm{T}}\mathbf{e} + \mathbf{e}\_{0}^{\mathrm{T}}\mathbf{P}\_{0}\,\mathbf{e}\_{0} = \min \quad \text{s.t.} \quad \mathbf{e} = \underset{\sim}{\mathbf{y}} - A\mathbf{x} \; \text{s.t.} \; \mathbf{e}\_{0} = \beta\_{0} - \mathbf{x} \tag{18}$$

which leads to the "normal equations" for the "*Least-Squares Collocation (LSC)*" solution

$$(N + P\_0) \cdot \mathfrak{x} = \underset{\sim}{\mathfrak{c}} + P\_0 \cdot \mathfrak{z}\_0 \tag{19}$$

which is *unique* even in the rank-deficient case (6).

If *<sup>Q</sup>*<sup>0</sup> is singular, but *<sup>e</sup>*<sup>0</sup> ∈ (*Q*0) with probability 1, then there must exist an *<sup>m</sup>* <sup>×</sup> 1 vector <sup>ν</sup><sup>0</sup> with

$$Q\_0 \,\nu^0 = -\epsilon\_0 = \text{x} - \beta\_0 \text{ with probability 1} \text{ (a.s.} \text{ \textdegree almost surely)}. \tag{20}$$

Thus, the principle (18) may be *equivalently* replaced by

$$\text{s.}\,^\text{T}\boldsymbol{\varepsilon} + \left(\boldsymbol{\nu}^0\right)^\text{T}\boldsymbol{Q}\_0\,\boldsymbol{\nu}^0 = \text{min} \quad \text{s.t.} \,\,\boldsymbol{\varepsilon} = \underset{\boldsymbol{\varepsilon}}{\boldsymbol{y}} - A\boldsymbol{x} \text{ and } (20), \tag{21}$$

which generates the *LSC solution* uniquely from the "modified normal equations"

$$\begin{array}{c} \hline \begin{array}{c} (I\_{\text{m}} + Q\_{0}N) \cdot \mathfrak{x}\_{\text{LSC}} = \beta\_{0} + Q\_{0}\mathfrak{c} \\ \hline \end{array} \end{array} \text{for } \mathfrak{c} := A^{\mathsf{T}} \underline{y} \tag{22}$$

or, alternatively, via the "*update formula*"

$$\text{Ar}\_{\text{LSC}} = \beta\_0 + Q\_0 \left( I\_{\text{w}} + N Q\_0 \right)^{-1} \left( \underline{c} - N \beta\_0 \right) \tag{23}$$

which exhibits the "*weak (local) unbiasedness*" of *x*˜LSC via

$$E\{\mathbf{\bar{x}}\_{\rm LSC}\} = \beta\_0 + Q\_0 \left(I\_{\rm ll} + NQ\_0\right)^{-1} \left(E\left\{\underline{c}\right\} - N\beta\_0\right) = \beta\_0 + 0 = E\{\mathbf{x}\}.\tag{24}$$

Consequently, the **MSE**-matrix of *x*˜LSC can be obtained from

$$\text{MSE } \{ \text{\"x}\_{\text{LSC}} \} = D \left\{ \text{\"x}\_{\text{LSC}} - \text{x\"} \right\} = \sigma\_0^2 Q\_0 \left( l\_m + N Q\_0 \right)^{-1} = \tag{25}$$

$$\rho = \sigma\_0^2 \left(N + P\_0\right)^{-1} \text{ if } P\_0 := Q\_0^{-1} \text{ exists},\tag{26}$$

whereas the dispersion matrix of *x*˜LSC itself is of *no relevance* here. It holds:

$$D\left[\mathbb{\ddot{x}}\_{\text{LSC}}\right] = D\left\{\mathbf{x}\right\} - D\left\{\mathbf{x} - \mathbb{\ddot{x}}\_{\text{LSC}}\right\} = \sigma\_0^2 Q\_0 - \sigma\_0^2 Q\_0 \left(I\_m + N Q\_0\right)^{-1}.\tag{27}$$

Again, *both residual vectors* are uniquely determined from

$$\left(\left(\overline{\mathbf{e}}\_{0}\right)\_{\text{LSC}} = \beta\_{0} - \overline{\mathbf{x}}\_{\text{LSC}} = -Q\_{0}\,\widehat{\nu}^{0}\_{\text{LSC}}\text{ for}\tag{28}$$

$$
\boldsymbol{\psi}\_{\text{LSC}}^{0} := \left(\boldsymbol{I}\_{\text{m}} + \boldsymbol{N}\boldsymbol{Q}\_{0}\right)^{-1} \left(\boldsymbol{c}\_{\text{v}} - \boldsymbol{N}\boldsymbol{\beta}\_{0}\right), \text{ and} \tag{29}
$$

$$\overline{\mathbf{x}}\_{\text{LSC}} = \underbrace{\mathcal{y}}\_{\sim} - A \,\, \overline{\mathbf{x}}\_{\text{LSC}} = \left[ I\_m - A \, Q\_0 \left( I\_m + N Q\_0 \right)^{-1} A^T \right] \left( \underline{\mathbf{y}} - A \beta\_0 \right) \tag{30}$$

with the control formula

$$\begin{array}{|c|}\hline \mathbf{A}^{\mathrm{T}}\vec{\boldsymbol{\varepsilon}}\_{\mathrm{LSC}} = \boldsymbol{\upvartheta}\_{\mathrm{LSC}}^{0} \\\\ \hline \end{array} \tag{31}$$

An optimal estimate of the *variance component* may now be obtained from

$$\begin{array}{rcl} \left(\boldsymbol{\delta}\_{0}^{2}\right)\_{\text{LSC}} &= \boldsymbol{n}^{-1} \cdot \left(\boldsymbol{\delta}\_{\text{LSC}}^{\text{T}} \cdot \mathbb{A}\_{\text{LSC}} + \left(\boldsymbol{\mathbb{t}}\_{\text{LSC}}^{0}\right)^{\text{T}} Q\_{0} \,\boldsymbol{\vartheta}\_{\text{LSC}}^{0}\right) = \\ &= \boldsymbol{n}^{-1} \cdot \left(\boldsymbol{\mathcal{y}}^{\text{T}} \boldsymbol{\mathcal{y}} - \boldsymbol{\mathcal{c}}\_{\text{LSC}}^{\text{T}} \mathbb{A}\_{\text{LSC}} - \boldsymbol{\beta}\_{0}^{\text{T}} \boldsymbol{\vartheta}\_{\text{LSC}}^{0}\right). \end{array} \tag{32}$$

#### **3. An Extension into the "Nonlinear World"**

#### *3.1. The Errors-In-Variables (EIV) Model*

In this scenario, the Gauss-Markov Model (GMM) is further *weakened* by allowing some or all of the entries in the coefficient matrix *A* to be observed. So, after introducing a corresponding *n* × *m* matrix of unknown random errors *EA*, the original GMM (1) turns into the EIV-Model

$$y = \underbrace{(A - E\_A) \cdot \xi}\_{n \times m} \cdot \xi + \varepsilon \; \; \; q := \text{rk} \; A = m < n \; \; \; \; \varepsilon \sim (0, \sigma\_0^2 I\_n) \; \; \; \tag{33}$$

with the vectorized form of *EA* being characterized through

$$\varepsilon\_{A} \, := \text{vec}\, E\_{A} \sim \left(0, \sigma\_{0}^{2} (I\_{m} \otimes I\_{n}) = \sigma\_{0}^{2} I\_{mn}\right), \text{ C} \{\varepsilon\_{i} e\_{A}\} = 0 \text{ (assumed)}.\tag{34}$$

Here, the vec operation transform a matrix into a vector by stacking all columns underneath each other while ⊗ denotes the *Kronecker-Zehfuss product* of matrices, defined by

$$\underset{p \times q}{\text{G}} \otimes \underset{r \times s}{\text{H}} = \left[ \underline{\text{g}}\_{ij} \cdot \text{H} \right]\_{pr \times qs} \text{ if } \text{G} = \left[ \underline{\text{g}}\_{ij} \right]. \tag{35}$$

In particular, the following key formula holds true:

$$\text{vec}\left(A\,B\,\mathbb{C}^{\text{I}}\right) = \left(\mathbb{C}\otimes A\right)\cdot \text{vec}\,B\tag{36}$$

for matrices of suitable size. Note that, in (34), the choice *QA* := *Imn* is a very special one. In general, *QA* may turn out singular whenever some parts of the matrix *A* remain unobserved (i.e., nonrandom).

In any case, thanks to the term

$$E\_A \cdot \xi = \left(\xi^\mathrm{T} \otimes I\_\mathrm{n}\right) \varepsilon\_{A\prime} \tag{37}$$

the model (33) needs to be treated in the "*nonlinear world*" even though the vector ξ may contain only incremental parameters. From now on, *A* is assumed to have *full column rank*, rk *A* =: *q* = *m*.

Following Schaffrin and Wieser [7], for instance, the "*Total-Least Squares Solution (TLSS)*" can be based on the principle

$$
\varepsilon^{\mathrm{T}}\varepsilon + \varepsilon\_A^{\mathrm{T}}\varepsilon\_A = \min \,\mathrm{s.t.}\,(\mathrm{33-34}),\tag{38}
$$

and the Lagrange function

$$\begin{array}{rcl}\Phi\left(\boldsymbol{\varepsilon},\boldsymbol{\varepsilon}\_{A},\boldsymbol{\xi},\boldsymbol{\lambda}\right):&=\boldsymbol{\varepsilon}^{\mathrm{T}}\boldsymbol{\varepsilon}+\boldsymbol{\varepsilon}\_{A}^{\mathrm{T}}\boldsymbol{\varepsilon}\_{A}+2\boldsymbol{\lambda}^{\mathrm{T}}\left(\boldsymbol{y}-\boldsymbol{A}\,\boldsymbol{\xi}-\boldsymbol{\varepsilon}+\boldsymbol{E}\_{A}\boldsymbol{\xi}\right)=\\&=\boldsymbol{\varepsilon}^{\mathrm{T}}\boldsymbol{\varepsilon}+\boldsymbol{\varepsilon}\_{A}^{\mathrm{T}}\boldsymbol{\varepsilon}\_{A}+2\boldsymbol{\lambda}^{\mathrm{T}}\left[\boldsymbol{y}-\boldsymbol{A}\,\boldsymbol{\xi}-\boldsymbol{\varepsilon}+\left(\boldsymbol{\xi}^{\mathrm{T}}\otimes\boldsymbol{I}\_{n}\right)\boldsymbol{e}\_{A}\right]\end{array}\tag{39}$$

which needs to be made *stationary*. The necessary Euler-Lagrange conditions then read:

$$\frac{1}{2}\frac{\partial}{\partial \varepsilon} \frac{\Phi}{\varepsilon} = \overline{\varepsilon} - \overline{\lambda} \doteq 0 \implies \overline{\lambda} = \overline{\varepsilon} \tag{40}$$

$$\frac{1}{2}\frac{\partial}{\partial \ e\_A} = \overline{e}\_A + \left(\xi\_{\text{TLS}} \otimes I\_n\right)\overline{\lambda} \doteq 0 \implies \mathcal{E}\_A = -\overline{\lambda}\,\mathcal{E}\_{\text{TLS}}^{\text{T}}\tag{41}$$

$$\frac{1}{2}\frac{\partial}{\partial \,\xi} = -(A - E\_A)^{\mathrm{T}}\lambda \neq 0 \implies A^{\mathrm{T}}\lambda = E\_A^{\mathrm{T}}\lambda = -\xi\_{\mathrm{TLS}} \cdot (\lambda^{\mathrm{T}}\lambda) \tag{42}$$

$$\frac{1}{2}\frac{\partial}{\partial \lambda} = y - A\,\xi\_{\rm TLS} - \overline{c} + \mathcal{E}\_A\,\xi\_{\rm TLS} \doteq 0 \Rightarrow y - A\,\xi\_{\rm TLS} = \lambda \,(1 + \xi\_{\rm TLS}^{\rm T}\,\xi\_{\rm TLS})\tag{43}$$

$$\Rightarrow \mathscr{c} - N\xi\_{\rm TLS}^{\xi} := A^{\rm T}(y - A \cdot \xi\_{\rm TLS}) = A^{\rm T}\lambda \, (1 + \xi\_{\rm TLS}^{\rm T} \cdot \xi\_{\rm TLS}) = -\xi\_{\rm TLS} \cdot \mathfrak{h}\_{\rm TLS} \tag{44}$$

for

$$\boldsymbol{\upmu\_{\rm TLS}} := (\boldsymbol{\hat{\lambda}^{\rm T}} \boldsymbol{\hat{\lambda}}) \cdot (\mathbf{1} + \boldsymbol{\xi}\_{\rm TLS}^{\rm T} \boldsymbol{\xi\_{\rm TLS}}) = \boldsymbol{\hat{\lambda}^{\rm T}} (\boldsymbol{y} - \boldsymbol{A} \, \boldsymbol{\xi\_{\rm TLS}}) = \boldsymbol{\tag{45}}$$

$$=\left(1+\hat{\xi}\_{\rm TLS}^{\rm T}\hat{\xi}\_{\rm TLS}\right)^{-1}\cdot\left(y-A\,\hat{\xi}\_{\rm TLS}\right)^{\rm T}\left(y-A\,\hat{\xi}\_{\rm TLS}\right)=\tag{46}$$

$$=\left(1+\mathcal{E}\_{\rm TLS}^{\rm T}\mathcal{E}\_{\rm TLS}\right)^{-1}\left[y^{\rm T}(y-A\,\mathcal{E}\_{\rm TLS})-\mathcal{E}\_{\rm TLS}^{\rm T}\left(c-N\,\mathcal{E}\_{\rm TLS}\right)\right] \quad \geq \quad \text{(47)}$$

$$\left(\begin{array}{ccccc} & & & \star & \star & \star\\ & \star & \star & \star & \star & \star \end{array}\right) \quad \geq \quad \text{(48)}$$

$$\mathbf{J} \Rightarrow (\mathbf{1} + \boldsymbol{\xi}\_{\mathrm{TLS}}^{\mathrm{T}} \boldsymbol{\xi}\_{\mathrm{TLS}}) \cdot \boldsymbol{\uprho}\_{\mathrm{TLS}} = \mathbf{y}^{\mathrm{T}} \mathbf{y} - \mathbf{c}^{\mathrm{T}} \boldsymbol{\upxi}\_{\mathrm{TLS}} + (\boldsymbol{\xi}\_{\mathrm{TLS}}^{\mathrm{T}} \boldsymbol{\upxi}\_{\mathrm{TLS}}) \cdot \boldsymbol{\uprho}\_{\mathrm{TLS}}$$

$$\Rightarrow \boxed{\mathring{\nu}\_{\rm TLS} = \y^{\rm T} y - c^{\rm T} \xi\_{\rm TLS}} \tag{48}$$

which needs to be solved in connection with the "modified normal equations" from (44), namely

$$\boxed{\left(N-\mathbb{M}\_{\text{TLS}}I\_{\text{m}}\right)\not{\mathcal{E}\_{\text{TLS}}}=\mathcal{C}}.\tag{49}$$

Due to the nonlinear nature of ξˆ TLS, it is not so easy to determine if it is an unbiased estimate, or how its MSE-matrix may exactly look like. First attempts of a rigorous error propagation have recently been undertaken by Amiri-Simkooei et al. [14] and by Schaffrin and Snow [15], but are beyond the scope of this paper.

Instead, both the *optimal residual vector <sup>e</sup>*˜TLS and the optimal residual matrix (*E*˜*A*)TLS are readily available through (40) and (43) as

$$
\hat{\sigma}\_{\rm TLS} = \hat{\lambda} = \left( y - A \, \hat{\xi}\_{\rm TLS} \right) \cdot \left( 1 + \hat{\xi}\_{\rm TLS}^{\rm T} \, \hat{\xi}\_{\rm TLS} \right)^{-1} \, \tag{50}
$$

and through (41) as

$$(E\_A)\_{\rm TLS} = -\mathbb{E}\_{\rm TLS} \cdot \xi\_{\rm TLS}^{\rm T} = -\left(y - A \cdot \xi\_{\rm TLS}\right) \cdot \left(1 + \xi\_{\rm TLS}^{\rm T} \cdot \xi\_{\rm TLS}\right)^{-1} \cdot \xi\_{\rm TLS}.\tag{51}$$

As optimal variance component estimate, it is now proposed to use the formula

$$(\boldsymbol{\delta}\_0^2)\_{\rm TLS} = (\boldsymbol{n} - \boldsymbol{m})^{-1} \cdot \left[ \boldsymbol{\mathbb{\tilde{e}}}\_{\rm TLS}^T \boldsymbol{\mathfrak{e}}\_{\rm TLS} + (\boldsymbol{\tilde{e}}\_A)\_{\rm TLS}^T (\boldsymbol{\tilde{e}}\_A)\_{\rm TLS} \right] = \boldsymbol{0}$$

$$=(n-m)^{-1} \cdot \mathbb{J}^{\text{T}} \mathbb{\hat{\lambda}}(1 + \xi\_{\text{TLS}}^{\text{T}} \xi\_{\text{TLS}}) = \mathbb{\hat{\nu}\_{\text{TLS}}}/(n-m),\tag{52}$$

in analogy to the previous estimates (11) and (32).

#### *3.2. A New Model: The EIV-Model with Random E*ff*ects (EIV-REM)*

In the following, the above EIV-Model (33-34) is *strengthened* by introducing stochastic prior information (p. i.) on the parameters which thereby change their character and become "random effects" as in (13-15). The EIV-REM can, therefore, be stated as

$$\underbrace{\partial y}\_{\sim} = \underbrace{(A - E\_A)}\_{n \times m} \cdot \mathbf{x} + \mathbf{c} \; , \; q := \text{rk} \; A \le \min\{m, n\} \; , \; \beta\_0 = \mathbf{x} + \mathbf{c}\_0 \; (\text{given}) , \tag{53}$$

with

$$\left[\boldsymbol{\varepsilon}\_{A} = \operatorname{vec} \boldsymbol{E}\_{A}\right] \sim \left( \begin{bmatrix} 0\\0\\0 \end{bmatrix}, \sigma\_{0}^{2} \begin{bmatrix} I\_{\text{lt}} & 0 & 0\\0 & I\_{\text{mn}} & 0\\0 & 0 & Q\_{0} \end{bmatrix} \right), \text{ Q}\_{0} \text{ symmetric and nnd.} \tag{54}$$

The first set of formulas will be derived by assuming that the weight matrix *P*<sup>0</sup> := *Q*−<sup>1</sup> <sup>0</sup> exists uniquely for the p. i. Then, the "*TLS collocation (TLSC)*" may be based on the principle

$$
\boldsymbol{e}^{\mathrm{T}}\boldsymbol{e} + \boldsymbol{e}\_{A}^{\mathrm{T}}\boldsymbol{e}\_{A} + \boldsymbol{e}\_{0}^{\mathrm{T}}\boldsymbol{P}\_{0}\boldsymbol{e}\_{0} = \min \,\mathrm{s.t.}\,(53\text{-}54),\tag{55}
$$

resp. on the *Lagrange function*

⎡ ⎢⎢⎢⎢⎢⎢⎢⎢⎣

$$\Phi\left(\boldsymbol{e},\boldsymbol{e}\_{A},\boldsymbol{e}\_{0},\lambda\right) := \boldsymbol{e}^{\mathrm{T}}\boldsymbol{e} + \boldsymbol{e}\_{A}^{\mathrm{T}}\boldsymbol{e}\_{A} + \boldsymbol{e}\_{0}^{\mathrm{T}}\boldsymbol{P}\_{0}\,\boldsymbol{e}\_{0} + 2\lambda^{\mathrm{T}}[\underbrace{\left(\boldsymbol{y} - \boldsymbol{A}\,\boldsymbol{\beta}\_{0} - \boldsymbol{e} + \left(\boldsymbol{\beta}\_{0}^{\mathrm{T}}\otimes\boldsymbol{I}\_{n}\right)\,\boldsymbol{e}\_{A} + \boldsymbol{A}\,\boldsymbol{e}\_{0} \quad \underbrace{-\mathrm{E}\_{A}\,\boldsymbol{e}\_{0}}\_{=-\left(\boldsymbol{e}\_{0}^{\mathrm{T}}\otimes\boldsymbol{I}\_{n}\right)\,\boldsymbol{e}\_{A}}]\tag{56}$$

which needs to be made stationary. The necessary Euler-Lagrange conditions then read:

$$\frac{1}{2}\frac{\partial}{\partial \varepsilon} \frac{\Phi}{\varepsilon} = \tilde{\varepsilon} - \hat{\lambda} \doteq 0 \implies \hat{\lambda} = \tilde{\varepsilon} \tag{57}$$

$$\frac{1}{2}\frac{\partial \Phi}{\partial \,\varepsilon\_{A}} = \tilde{\varepsilon}\_{A} + \left[\left(\beta \mathbf{o} - \tilde{\varepsilon}\mathbf{o}\right) \otimes I\_{\mathrm{I}}\right] \hat{\lambda} \triangleq 0 \implies \tilde{E}\_{A} = -\hat{\lambda} \cdot \left(\beta \mathbf{o} - \tilde{\varepsilon}\mathbf{o}\right)^{\mathrm{T}} =: -\hat{\lambda} \cdot \left.\tilde{\mathbf{x}}\_{\mathrm{T}\mathrm{LSC}}^{\mathrm{T}}\right| \tag{58}$$

$$\frac{1}{2}\frac{\partial}{\partial \varepsilon\_0} = P\_0 \ \overline{\varepsilon}\_0 + \left(A - \overline{E}\_A\right)^T \hat{\lambda} \ \triangleq 0 \ \Rightarrow \ A^T \hat{\lambda} = \bar{E}\_A^T \hat{\lambda} - P\_0 \ \overline{\varepsilon}\_0 \ \Rightarrow \tag{59}$$

$$\Rightarrow A^{\mathsf{T}}\hat{\lambda} = -\mathsf{\tilde{x}}\_{\text{TLSC}} \cdot (\hat{\lambda}^{\mathsf{T}}\hat{\lambda}) + \mathsf{\hat{v}}\_{\text{TLSC}}^{0} \text{ for } \mathsf{\hat{v}}\_{\text{TLSC}}^{0} := -P\_{0}\ \mathsf{\tilde{x}}\_{0} = P\_{0}(\boldsymbol{\beta}\_{0} - \mathsf{\tilde{x}}\_{\text{TLSC}}) \tag{60}$$

$$\frac{1}{2}\frac{\partial\Phi}{\partial\lambda} = \underset{\sim}{y} -A\left(\left\{\mathbf{0} - \tilde{\mathbf{e}}\_{0}\right\} - \tilde{\mathbf{e}} + \tilde{E}\_{A}\left(\left\{\mathbf{0} - \tilde{\mathbf{e}}\_{0}\right\} \right) \triangleq \mathbf{0} \implies \underset{\sim}{y} -A\,\tilde{\mathbf{x}}\_{\text{TLSC}} = \hat{\lambda}\left(1 + \hat{\mathbf{x}}\_{\text{TLSC}}^{\text{T}}\hat{\mathbf{x}}\_{\text{TLSC}}\right) \Rightarrow \tag{61}$$

$$\hat{\lambda} \Rightarrow \hat{\lambda} = \underbrace{(y - A \,\, \text{\"x}\_{\text{TLSC}}) \cdot \left(1 + \text{\"x}\_{\text{TLSC}}^{\text{T}} \,\, \text{\"x}\_{\text{TLSC}}\right)^{-1}}\_{\sim \text{(62)}}.\tag{62}$$

Combining (60) with (62) results in

(*c* <sup>∼</sup> <sup>−</sup> *<sup>N</sup> <sup>x</sup>*˜TLSC) · (<sup>1</sup> <sup>+</sup> *<sup>x</sup>*˜ T TLSC *x*˜TLSC) <sup>−</sup><sup>1</sup> <sup>=</sup> *<sup>A</sup>*Tλ<sup>ˆ</sup> <sup>=</sup> <sup>−</sup>*x*˜TLSC(<sup>1</sup> <sup>+</sup> *<sup>x</sup>*˜ T TLSC *x*˜TLSC) −2 (*y* ∼ − *A x*˜TLSC) <sup>T</sup>(*y* ∼ − *A x*˜TLSC) + νˆ 0 TLSC, (63)

and finally in

$$\left(\text{N} + (1 + \hat{\mathfrak{x}}\_{\text{TLSC}}^{\text{T}} \hat{\mathfrak{x}}\_{\text{TLSC}})P\_0 - \hat{\mathfrak{y}}\_{\text{TLSC}}I\_m\right)\hat{\mathfrak{x}}\_{\text{TLSC}} = c\_s + P\_0 \left. \beta\_0 \cdot \left(1 + \hat{\mathfrak{x}}\_{\text{TLSC}}^{\text{T}} \hat{\mathfrak{x}}\_{\text{TLSC}}\right)\right|\tag{64}$$

where

$$\begin{array}{lcl}\hline \boldsymbol{\upmu}\_{\text{TLSC}} := (1 + \boldsymbol{\upmu}\_{\text{TLSC}}^{\text{T}} \boldsymbol{\upmu}\_{\text{TLSC}})^{-1} (\boldsymbol{y} - A \boldsymbol{\upmu}\_{\text{TLSC}})^{\text{T}} (\boldsymbol{y} - A \boldsymbol{\upmu}\_{\text{TLSC}}), \text{ and} \\\boldsymbol{\upmu}\_{\text{TLSC}}^{0} := -P\_{0} (\boldsymbol{\upmu}\_{0} - \boldsymbol{\upmu}\_{\text{TLSC}}) = -P\_{0} \boldsymbol{\upmu}\_{0}, \text{ provided that } P\_{0} \text{ exists.} \\\hline \end{array} \tag{65}$$

In the more general case of a *singular* matrix *Q*0, an approach similar to (20) can be followed, leading to the equation system

$$\mathbb{E}\left[\left(1+\mathbb{\tilde{x}}\_{\text{TLSC}}^{\text{T}}\mathbb{\tilde{x}}\_{\text{TLSC}}\right)\cdot I\_{\text{m}}+Q\_{0}N-\mathbb{\hat{\nu}}\_{\text{TLSC}}\cdot Q\_{0}\right]\mathbb{\tilde{x}}\_{\text{TLSC}}=\beta\_{0}\cdot\left(1+\mathbb{\tilde{x}}\_{\text{TLSC}}^{\text{T}}\mathbb{\tilde{x}}\_{\text{TLSC}}\right)+Q\_{0}\mathbb{\tilde{x}}\_{\text{n}}\Big|\tag{66}$$

that needs to be solved in connection with (65). Obviously,

$$\begin{array}{ccccc}\mathfrak{x}\_{\text{TLSC}} & \longmapsto & \beta\_0 \qquad \text{if} \quad Q\_0 & \longmapsto & 0, \quad \text{and} \\ \mathfrak{x}\_{\text{TLSC}} & \longmapsto & \mathfrak{x}\_{\text{TLS}} \quad \text{if} \quad P\_0 & \longmapsto & 0. \end{array} \tag{67}$$

Again, it is still unclear if *x*˜TLSC represents an unbiased prediction of the vector *x* of random effects. Also, very little (if anything) is known about the corresponding MSE-matrix of *x*˜TLSC. The answers to these open problems will be left for a future contribution. It is, however, possible to find the respective *residual vectors*/*matrices* represented as follows:

$$\tilde{\alpha}\_{\rm TLSC} = \hat{\lambda} = \left(\underline{y} - A \,\, \tilde{\mathbf{x}}\_{\rm TLSC}\right) \cdot \left(1 + \tilde{\mathbf{x}}\_{\rm TLSC}^{\rm T} \,\, \tilde{\mathbf{x}}\_{\rm TLSC}\right)^{-1},\tag{68}$$

$$(\mathbf{E}\_A)\_{\rm TLS} = -\boldsymbol{\mathbb{\hat{A}}} \cdot \hat{\mathbf{x}}\_{\rm TLS}^{\rm T} = -(y - A \,\hat{\mathbf{x}}\_{\rm TLS}) \cdot \left(1 + \hat{\mathbf{x}}\_{\rm TLS}^{\rm T} \,\hat{\mathbf{x}}\_{\rm TLS}\right)^{-1} \cdot \hat{\mathbf{x}}\_{\rm TLS'}^{\rm T} \tag{69}$$

$$(\mathfrak{F}\_0)\_{\rm TLSC} = -\mathrm{Q}\_0 \cdot \mathfrak{h}\_{\rm TLSC}^0 = \mathfrak{f}\_0 - \mathfrak{F}\_{\rm TLSC} \tag{70}$$

while a suitable formula for the variance component is suggested as

$$\left(\left(\delta\_{0}^{2}\right)\_{\text{TLSC}} = n^{-1}\cdot\;\left[\;\left(\mathbb{M}\_{\text{TLSC}} + \left(\mathbb{M}\_{\text{TLSC}}^{0}\right)^{\text{T}}Q\_{0}\left(\mathbb{M}\_{\text{TLSC}}^{0}\right)\right]\right.\tag{71}$$

#### **4. Conclusions and Outlook**

Key formulas have been developed successfully to optimally determine the parameters and residuals within the new 'EIV-Model with p. i.' (or EIV-REM) which turns out to be more general than the other three models considered here (GMM, REM, EIV-Model). In particular, it is quite obvious that


Hence the new EIV-REM can indeed serve as a universal representative of the whole class of models presented here.

Therefore, in a follow-up paper, it is planned to also cover more general dispersion matrices for *e* and *eA* in (54), similarly to the work by Schaffrin et al. [16] for the EIV-Model with singular dispersion matrices for *eA*.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


© 2020 by the author. 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* **Weighted Total Least Squares (WTLS) Solutions for Straight Line Fitting to 3D Point Data**

#### **Georgios Malissiovas 1,\*, Frank Neitzel 1, Sven Weisbrich <sup>1</sup> and Svetozar Petrovic 1,2**


Received: 14 July 2020; Accepted: 25 August 2020; Published: 29 August 2020

**Abstract:** In this contribution the fitting of a straight line to 3D point data is considered, with Cartesian coordinates *xi*, *yi*, *zi* as observations subject to random errors. A direct solution for the case of equally weighted and uncorrelated coordinate components was already presented almost forty years ago. For more general weighting cases, iterative algorithms, e.g., by means of an iteratively linearized Gauss–Helmert (GH) model, have been proposed in the literature. In this investigation, a new direct solution for the case of pointwise weights is derived. In the terminology of total least squares (TLS), this solution is a direct weighted total least squares (WTLS) approach. For the most general weighting case, considering a full dispersion matrix of the observations that can even be singular to some extent, a new iterative solution based on the ordinary iteration method is developed. The latter is a new iterative WTLS algorithm, since no linearization of the problem by Taylor series is performed at any step. Using a numerical example it is demonstrated how the newly developed WTLS approaches can be applied for 3D straight line fitting considering different weighting cases. The solutions are compared with results from the literature and with those obtained from an iteratively linearized GH model.

**Keywords:** 3D straight line fitting; total least squares (TLS); weighted total least squares (WTLS); nonlinear least squares adjustment; direct solution; singular dispersion matrix; laser scanning data

#### **1. Introduction**

Modern geodetic instruments, such as terrestrial laser scanners, provide the user directly with 3D coordinates in a Cartesian coordinate system. However, in most cases these 3D point data are not the final result. For an analysis of the recorded data or for a representation using computer-aided design (CAD), a line, curve or surface approximation with a continuous mathematical function is required.

In this contribution the fitting of a spatial straight line is discussed considering the coordinate components *xi*, *yi*, *zi* of each point *Pi* as observations subject to random errors, which results in a nonlinear adjustment problem. An elegant direct least squares solution for the case of equally weighted and uncorrelated observations has already been presented in 1982 by Joviˇci´c et al. [1]. Unfortunately, this article was very rarely cited in subsequent publications, which is probably due to the fact that it was written in Croatian language. Similar least squares solutions, direct as well, have been published by Kahn [2] and Drixler ([3], pp. 46–47) some years later. In these contributions, it was shown that the problem of fitting a straight line to 3D point data can be transformed into an eigenvalue problem. An iterative least squares solution for fitting a straight line to equally weighted

and uncorrelated 3D points has been presented by Späth [4], by minimizing the sum of squared orthogonal distances of the observed points to the requested straight line.

For more general weighting schemes iterative least squares solutions have been presented by Kupferer [5] and Snow and Schaffrin [6]. In both contributions various nonlinear functional models were introduced and tested and the Gauss-Newton approach has been employed for an iterative least squares solution, which involves the linearization of the functional model to solve the adjustment problem within the Gauss–Markov (GM) or the Gauss–Helmert (GH) models. The linearization of originally nonlinear functional models is a very popular approach in adjustment calculation, where the solution is then determined by iteratively linearized GM or GH models, see e.g., the textbooks by Ghilani [7], Niemeier [8] or Perovic [9]. Pitfalls to be avoided in the iterative adjustment of nonlinear problems have been pointed out by Pope [10] and Lenzmann and Lenzmann [11].

Fitting a straight line to 3D point data can also be considered as an adjustment problem of type total least squares (TLS) for an errors-in-variables (EIV) model, as already pointed out by Snow and Schaffrin [6]. For some weighting cases, problems expressed within the EIV model can have a direct solution using singular value decomposition (SVD). This solution is known under the name Total Least Squares (TLS) and has been firstly proposed by Golub and Van Loan [12], Van Huffel and Vandewalle [13] or Van Huffel and Vandewalle ([14], p. 33 ff.). The TLS solution is obtained by computing the roots of a polynomial, i.e., by solving the characteristic equation of the eigenvalues, which is identical in the case of fitting a straight line to 3D point data to the solutions of Joviˇci´c et al. [1], Kahn [2] and Drixler ([3], pp. 46–47). This is something that has been observed by Malissiovas et al. [15], where a relationship has been presented between TLS and the direct least squares solution of the same adjustment problem, while postulating uncorrelated and equally weighted observations.

To involve more general weight matrices in the adjustment procedure, iterative algorithms have been presented in the TLS literature without linearizing the underlying problem by Taylor series at any step of the solution process. These are algorithmic approaches known as weighted total least squares (WTLS), presented e.g., by Schaffrin and Wieser [16], Shen et al. [17] or Amiri and Jazaeri [18]. A good overview of such algorithms, as well as alternative solution strategies can be found in the dissertation theses of Snow [19] and Malissiovas [20]. An attempt to find a WTLS solution for straight line fitting to 3D point data was made by Guo et al. [21].

To avoid confusion, it is to clarify that the terms TLS and WTLS refer to algorithmic approaches for obtaining a least squares solution, which is either direct or iterative but without linearizing the problem by Taylor series at any step. This follows the statement of Snow ([19], p. 7), that "the terms total least-squares (TLS) and TLS solution [...] will mean the least-squares solution within the EIV model without linearization". Of course, a solution within the GH model is more general in the sense that it can be utilized to solve any nonlinear adjustment problem, while TLS and WTLS algorithms can treat only a certain class of nonlinear adjustment problems. This has been firstly discussed by Neitzel and Petrovic [22] and Neitzel [23], who showed that the TLS estimate within an EIV model can be identified as a special case of the method of least squares within an iteratively linearized GH model.

To the extent of our knowledge, a WTLS algorithm for fitting a straight line to 3D point data has not been presented yet. Therefore, in this study we derive two novel WTLS algorithms for the discussed adjustment problem considering two different cases of stochastic models:


The adjustment problem resulting from case (i) can still be solved directly, i.e., a direct WTLS solution, presented in Section 2.1 of this paper. For case (ii) an iterative solution without linearizing the problem by Taylor series is derived in Section 2.2, i.e., an iterative WTLS solution. Both solutions are based on the work of Malissiovas [20], where similar algorithms have been presented for the solution of other typical geodetic tasks, such as straight line fitting to 2D point data, plane fitting to 3D point data and 2D similarity coordinate transformation. The WTLS solution for straight line fitting to 3D point data will be derived from a geodetic point of view by means of introducing residuals for all observations, formulating appropriate condition and constraint equations, setting up and solving the resulting normal equation systems.

#### **2. Straight Line Fitting to 3D Point Data**

A straight line in 3D space can be expressed by

$$\frac{y - y\_0}{a} = \frac{\mathbf{x} - \mathbf{x}\_0}{b} = \frac{z - z\_0}{c},\tag{1}$$

as explained e.g., in the handbook of mathematics by Bronhstein et al. ([24], p. 217). This equation defines a line that passes through a point with coordinates *x*0, *y*<sup>0</sup> and *z*<sup>0</sup> and is parallel to a direction vector with components *a*, *b* and *c*. Since the number of unknown parameters is six, two additional constraints between the unknown parameters have to be taken into account for a solution, as Snow and Schaffrin [6] have clearly explained. Proper constraints will be selected at a later point of the derivations.

Considering an overdetermined configuration for which we want to obtain a least squares solution, we introduce for all points *Pi* residuals *vxi*, *vyi* , *vzi* for the observations *xi*, *yi*, *zi* assuming that they are normally distributed with **<sup>v</sup>** ∼ (**0**, *<sup>σ</sup>*<sup>2</sup> <sup>0</sup> **Q**LL). Here the 1 × 3*n* vector **v** contains the residuals *vxi*, *vyi* , *vzi*, **<sup>Q</sup>**LL is the 3*<sup>n</sup>* × <sup>3</sup>*<sup>n</sup>* corresponding cofactor matrix and *<sup>σ</sup>*<sup>2</sup> <sup>0</sup> the theoretical variance factor. Based on (1), the nonlinear conditions equations

$$\begin{aligned} a(\mathbf{x}\_i + \mathbf{v}\_{\mathbf{x}\_i} - \mathbf{x}\_0) - b(\mathbf{y}\_i + \mathbf{v}\_{\mathbf{y}\_i} - \mathbf{y}\_0) &= \mathbf{0}, \\\\ b(z\_i + \mathbf{v}\_{z\_i} - z\_0) - c(\mathbf{x}\_i + \mathbf{v}\_{\mathbf{x}\_i} - \mathbf{x}\_0) &= \mathbf{0}, \\\\ c(y\_i + \mathbf{v}\_{\mathbf{y}\_i} - \mathbf{y}\_0) - a(z\_i + \mathbf{v}\_{z\_i} - z\_0) &= \mathbf{0}, \end{aligned} \tag{2}$$

can be formulated for each point *Pi*, with *i* = 1, ... , *n* and *n* being the number of observed points. A functional model for this problem can be expressed by two of these three nonlinear condition equations per observed point. Using all three condition equations for solving the problem would lead to a singular normal matrix due to linearly dependent normal equations.

#### *2.1. Direct Total Least Squares Solution for Equally Weighted and Uncorrelated Observations*

Considering all coordinate components *xi*, *yi*, *zi* for all points *Pi* as equally weighted and uncorrelated observations with

$$p\_{y\_i} = p\_{x\_i} = p\_{z\_i} = 1 \,\,\forall i,\tag{3}$$

a least squares solution can be derived by minimizing the objective function

$$
\Omega = \sum\_{i=1}^{n} v\_{x\_i}^2 + v\_{y\_i}^2 + v\_{z\_i}^2 \to \min. \tag{4}
$$

A direct least squares solution of this problem, respectively a TLS solution, has been presented by Malissiovas et al. [15] and Malissiovas [20]. According to these investigations, equally weighted residuals correspond to the normal distances

$$D\_i^2 = v\_{x\_i}^2 + v\_{y\_i}^2 + v\_{z\_i \prime}^2 \tag{5}$$

as deviation measures between the observed 3D point cloud and the straight line to be computed. Thus, the objective function can be written equivalently as

$$
\Omega = \sum\_{i=1}^{n} v\_{x\_i}^2 + v\_{y\_i}^2 + v\_{z\_i}^2 = \sum\_{i=1}^{n} D\_i^2 \to \min. \tag{6}
$$

An expression for the squared normal distances

$$D^2 = \frac{\left[a(\mathbf{x} - \mathbf{x}\_0) - b(y - y\_0)\right]^2 + \left[b(z - z\_0) - c(\mathbf{x} - \mathbf{x}\_0)\right]^2 + \left[c(y - y\_0) - a(z - z\_0)\right]^2}{a^2 + b^2 + c^2},\tag{7}$$

can be found in the handbook of Bronhstein et al. ([24], p. 218).

An appropriate parameterization of the problem involves the substitution of the unknown parameters *y*0, *x*<sup>0</sup> and *z*<sup>0</sup> with the coordinates of the center of mass of the observed 3D point data. A proof for this parameter replacement has been given by Joviˇci´c et al. [1] and concerns only the case of equally weighted and uncorrelated observations. A solution for the line parameters can be computed by minimizing the objective function (6), under the constraint

$$a^2 + b^2 + c^2 = 1,\tag{8}$$

or by searching for stationary points of the Lagrange function

$$K = \sum\_{i=1}^{n} D\_i^2 - k(a^2 + b^2 + c^2 - 1),\tag{9}$$

with *k* denoting the Lagrange multiplier. The line parameters can be computed directly from the normal equations, either by solving a characteristic cubic equation or an eigenvalue problem. A detailed explanation of this approach was given by Malissiovas ([20], p. 74 ff.). In the following subsections we will give up the restriction of equally weighted and uncorrelated coordinate components and derive least squares solutions considering more general weighting schemes for the observations.

#### *2.2. Direct Weighted Total Least Squares Solution*

In this case we consider the coordinate components *xi*, *yi*, *zi* of each point *Pi* to be uncorrelated and of equal precision with

$$
\sigma\_{y\_i} = \sigma\_{xi} = \sigma\_{zi} = \sigma\_i \quad \forall i,\tag{10}
$$

yielding the corresponding (pointwise) weights

$$p\_i = \frac{1}{\sigma\_i^2} \quad \forall i. \tag{11}$$

A least squares solution of the problem can be found by minimizing the objective function

$$\begin{split} \Omega &= \sum\_{i=1}^{n} p\_{X\_i} v\_{x\_i}^2 + p\_{Y\_i} v\_{y\_i}^2 + p\_{z\_i} v\_{z\_i}^2 \to \min \\ &= \sum\_{i=1}^{n} p\_i \left( v\_{x\_i}^2 + v\_{y\_i}^2 + v\_{z\_i}^2 \right) \to \min. \end{split} \tag{12}$$

Taking into account the relation of the squared residuals with the normal distances of Equation (5), it is possible to express the objective function as

$$
\Omega = \sum\_{i=1}^{n} p\_i D\_i^2 \to \min. \tag{13}
$$

Using the constraint (8), the expression of the normal distances can be written as

$$D\_i^2 = \left[ a(\mathbf{x}\_i - \mathbf{x}\_0) - b(y\_i - y\_0) \right]^2 + \left[ b(z\_i - z\_0) - c(\mathbf{x}\_i - \mathbf{x}\_0) \right]^2 + \left[ c(y\_i - y\_0) - a(z\_i - z\_0) \right]^2. \tag{14}$$

A further simplification of the problem is possible, by replacing the unknown parameters *y*0, *x*<sup>0</sup> and *z*<sup>0</sup> with the coordinates of the weighted center of mass of the 3D point data

$$y\_0 = \frac{\sum\_{i=1}^n p\_i y\_i}{\sum\_{i=1}^n p\_i}, \mathbf{x}\_0 = \frac{\sum\_{i=1}^n p\_i \mathbf{x}\_i}{\sum\_{i=1}^n p\_i}, z\_0 = \frac{\sum\_{i=1}^n p\_i z\_i}{\sum\_{i=1}^n p\_i}. \tag{15}$$

It can be proven that this parameterization is allowed and possible, following the same line of thinking as in the study of Joviˇci´c et al. [1]. Therefore, for this weighted case we can write the normal distances as

$$D\_i^2 = \left(a\mathbf{x}\_i' - by\_i'\right)^2 + \left(bz\_i' - c\mathbf{x}\_i'\right)^2 + \left(cy\_i' - az\_i'\right)^2,\tag{16}$$

with the coordinates reduced to the weighted center of mass of the observed 3D point data

$$y\_i' = y\_i - \frac{\sum\_{i=1}^n p\_i y\_i}{\sum\_{i=1}^n p\_i}, \ x\_i' = x\_i - \frac{\sum\_{i=1}^n p\_i x\_i}{\sum\_{i=1}^n p\_i}, \ z\_i' = z\_i - \frac{\sum\_{i=1}^n p\_i z\_i}{\sum\_{i=1}^n p\_i}. \tag{17}$$

Searching for stationary points of the Lagrange function

$$\begin{aligned} K &= \sum\_{i=1}^{n} D\_i^2 - k(a^2 + b^2 + c^2 - 1) \\ &= \sum\_{i=1}^{n} \left( a\mathbf{x}\_i' - b\mathbf{y}\_i' \right)^2 + \left( b\mathbf{z}\_i' - c\mathbf{x}\_i' \right)^2 + \left( c\mathbf{y}\_i' - a\mathbf{z}\_i' \right)^2 - k(a^2 + b^2 + c^2 - 1) \end{aligned} \tag{18}$$

leads to the system of normal equations

$$\frac{\partial K}{\partial a} = 2a \left( \sum\_{i=1}^{n} p\_i \mathbf{x}\_i'^2 + \sum\_{i=1}^{n} p\_i \mathbf{z}\_i'^2 - k \right) - 2b \sum\_{i=1}^{n} p\_i y\_i' \mathbf{x}\_i' - 2c \sum\_{i=1}^{n} p\_i y\_i' \mathbf{z}\_i' = 0,\tag{19}$$

$$\frac{\partial K}{\partial b} = 2b\left(\sum\_{i=1}^{n} p\_i y\_i'^2 + \sum\_{i=1}^{n} p\_i z\_i'^2 - k\right) - 2a\sum\_{i=1}^{n} p\_i y\_i' \mathbf{x}\_i' - 2c\sum\_{i=1}^{n} p\_i \mathbf{x}\_i' \mathbf{z}\_i' = 0,\tag{20}$$

$$\frac{\partial K}{\partial c} = 2c \left( \sum\_{i=1}^{n} p\_i y\_i'^2 + \sum\_{i=1}^{n} p\_i \mathbf{x}\_i'^2 - k \right) - 2a \sum\_{i=1}^{n} p\_i y\_i' z\_i' - 2b \sum\_{i=1}^{n} p\_i \mathbf{x}\_i' z\_i' = 0,\tag{21}$$

and

$$\frac{\partial K}{\partial k} = -\left(a^2 + b^2 + c^2 - 1\right) = 0.\tag{22}$$

Equations (19) to (21) are a homogeneous system of equations which are linear in the unknown line parameters *a*, *b* and *c*, and has a nontrivial solution if

$$
\begin{vmatrix}
(d\_1 - k) & d\_4 & d\_5 \\
 & d\_4 & (d\_2 - k) & d\_6 \\
 & d\_5 & d\_6 & (d\_3 - k)
\end{vmatrix} = 0,\tag{23}
$$

with the elements of the determinant

$$\begin{aligned} d\_1 &= \sum\_{i=1}^n p\_i \mathbf{x}\_i'^2 + \sum\_{i=1}^n p\_i \mathbf{z}\_i'^2, \ d\_2 = \sum\_{i=1}^n p\_i \mathbf{y}\_i'^2 + \sum\_{i=1}^n p\_i \mathbf{z}\_i'^2, \ d\_3 = \sum\_{i=1}^n p\_i \mathbf{y}\_i'^2 + \sum\_{i=1}^n p\_i \mathbf{x}\_i'^2, \\ d\_4 &= -\sum\_{i=1}^n p\_i \mathbf{y}\_i' \mathbf{x}\_i', \ d\_5 = -\sum\_{i=1}^n p\_i \mathbf{y}\_i' \mathbf{z}\_i', \ d\_6 = -\sum\_{i=1}^n p\_i \mathbf{x}\_i' \mathbf{z}\_i'. \end{aligned} \tag{24}$$

Equation (23) is a cubic characteristic equation with the unknown parameter *k* and a solution is given by Bronhstein et al. ([24], p. 63). The unknown line parameters *a*, *b* and *c* can be computed either by substituting *k*min into Equations (19)–(21) under the constraint (22) or by transforming the equation system and solving an eigenvalue problem, i.e., the direct WTLS solution.

#### *2.3. Iterative Weighted Total Least Squares Solution*

In this case, we consider the fact that the coordinate components *xi*, *yi*, *zi* of each point *Pi* are not original observations. They are either

(i) derived from single point determinations, e.g., using polar elementary measurements of slope distance *ρi*, direction *φ<sup>i</sup>* and tilt angle *θi*, e.g., from measurements with a terrestrial laser scanner, so that the coordinates result from the functional relationship

$$\begin{aligned} x\_i &= \rho\_i \cos \phi\_i \sin \theta\_{i\prime} \\ y\_i &= \rho\_i \sin \phi\_i \cos \theta\_{i\prime} \\ z\_i &= \rho\_i \cos \theta\_i. \end{aligned} \tag{25}$$

Using the standard deviations *σρ<sup>i</sup>* , *σφ<sup>i</sup>* and *σθ <sup>i</sup>* of the elementary measurements, a 3 × 3 variance-covariance matrix can be provided for the coordinate components *xi*, *yi*, *zi* of each point *Pi* by variance-covariance propagation based on (25). Covariances among the *n* different points do not occur in this case.

or

(ii) derived from a least squares adjustment of an overdetermined geodetic network. From the solution within the GM or GH model, the 3*n* × 3*n* variance-covariance matrix of the coordinate components *xi*, *yi*, *zi* of all points *Pi* can be obtained. This matrix is a full matrix, considering also covariances among the *n* different points. It may even have a rank deficiency in case the coordinates were determined by a free network adjustment.

The variances and covariances from (i) or (ii) can then be assembled in a dispersion matrix

$$\mathbf{Q}\_{\rm LL} = \begin{bmatrix} \mathbf{Q}\_{\rm xx} & \mathbf{Q}\_{\rm xy} & \mathbf{Q}\_{\rm xz} \\ \mathbf{Q}\_{\rm yx} & \mathbf{Q}\_{\rm yy} & \mathbf{Q}\_{\rm yz} \\ \mathbf{Q}\_{\rm xx} & \mathbf{Q}\_{\rm xy} & \mathbf{Q}\_{\rm xz} \end{bmatrix}, \tag{26}$$

for the coordinate components as "derived observations". In case of a non-singular dispersion matrix, it is possible to compute the respective weight matrix

$$\mathbf{P} = \begin{bmatrix} \mathbf{P}\_{\text{xx}} & \mathbf{P}\_{\text{xy}} & \mathbf{P}\_{\text{xz}} \\ \mathbf{P}\_{\text{yx}} & \mathbf{P}\_{\text{yy}} & \mathbf{P}\_{\text{yz}} \\ \mathbf{P}\_{\text{xx}} & \mathbf{P}\_{\text{xy}} & \mathbf{P}\_{\text{zz}} \end{bmatrix} = \mathbf{Q}\_{\text{LL}}^{-1}. \tag{27}$$

Starting with the functional model (2), we choose to work with the nonlinear condition equations

$$\begin{aligned} \mathcal{L}(y\_i + v\_{y\_i} - y\_0) - a(z\_i + v\_{z\_i} - z\_0) &= 0, \\\\ \mathcal{L}(\mathbf{x}\_i + v\_{x\_i} - \mathbf{x}\_0) - b(z\_i + v\_{z\_i} - z\_0) &= 0. \end{aligned} \tag{28}$$

Two additional constraints must be taken into account in this case. For example, Snow and Schaffrin [6] proposed the constraints

$$\begin{aligned} a^2 + b^2 + c^2 &= 1, \\\\ a y\_0 + b x\_0 + c z\_0 &= 0. \end{aligned} \tag{29}$$

A further development of an algorithmic solution using these constraints is possible. However, we choose to parameterize the functional model so that an additional constraint is unnecessary and the equations become simpler. Therefore, we choose to fix one coordinate of the point on the line

$$z\_0 = \frac{\sum\_{i=1}^{n} z\_i}{n} = \overline{z}\_\prime \tag{30}$$

as well as one of the unknown components of the direction vector

$$
\mathfrak{c} = 1.\tag{31}
$$

This simplification of the functional model has been used by Boroviˇcka et al. [25] for a similar practical example, this of estimating the trajectory of a meteor. As already mentioned in that work, the resulting direction vector of the straight line can be transformed to a unit vector afterwards, i.e., a solution using the constraint *a*<sup>2</sup> + *b*<sup>2</sup> + *c*<sup>2</sup> = 1. The simplified functional model can be written in vector notation as

$$\begin{aligned} \left( (\mathbf{y}\_{\mathbf{c}} + \mathbf{v}\_{\mathbf{y}} - \mathbf{e}y\_0) - a(\mathbf{z}\_{\mathbf{c}} + \mathbf{v}\_{\mathbf{z}} - \mathbf{e}\overline{z}) = \mathbf{0}, \\\\ \left( \mathbf{x}\_{\mathbf{c}} + \mathbf{v}\_{\mathbf{x}} - \mathbf{e}x\_0 \right) - b(\mathbf{z}\_{\mathbf{c}} + \mathbf{v}\_{\mathbf{z}} - \mathbf{e}\overline{z}) = \mathbf{0}, \end{aligned} \tag{32}$$

with the vectors **x**c, **y**<sup>c</sup> and **z**<sup>c</sup> containing the coordinates of the observed points and vectors **v**x, **v**<sup>y</sup> and **v**<sup>z</sup> the respective residuals. Vector **e** is a vector of ones, with length equal to the number of observed points *n*.

A weighted least squares solution of this problem can be derived by minimizing the objective function

$$\boldsymbol{\Omega} = \mathbf{v}\_{\mathbf{x}}^{\mathrm{T}} \mathbf{P}\_{\mathrm{xx}} \mathbf{v}\_{\mathbf{x}} + \mathbf{v}\_{\mathbf{y}}^{\mathrm{T}} \mathbf{P}\_{\mathrm{YY}} \mathbf{v}\_{\mathbf{y}} + \mathbf{v}\_{\mathbf{z}}^{\mathrm{T}} \mathbf{P}\_{\mathrm{zz}} \mathbf{v}\_{\mathbf{z}} + 2 \left. \mathbf{v}\_{\mathbf{x}}^{\mathrm{T}} \mathbf{P}\_{\mathrm{XY}} \mathbf{v}\_{\mathbf{y}} + 2 \left. \mathbf{v}\_{\mathbf{x}}^{\mathrm{T}} \mathbf{P}\_{\mathrm{xx}} \mathbf{v}\_{\mathbf{z}} + 2 \left. \mathbf{v}\_{\mathbf{y}}^{\mathrm{T}} \mathbf{P}\_{\mathrm{YY}} \mathbf{v}\_{\mathbf{z}} \right| \right. \tag{33}$$

or by searching for stationary points of the Lagrange function

$$\begin{aligned} K &= \Omega - 2\mathbf{k}\_1^T \left[ (\mathbf{y}\_\mathbf{c} + \mathbf{v}\_\mathbf{y} - \mathbf{e}y\_0) - a(\mathbf{z}\_\mathbf{c} + \mathbf{v}\_\mathbf{z} - \mathbf{e}\mathbf{\tilde{z}}) \right] \\ &- 2\mathbf{k}\_2^T \left[ (\mathbf{x}\_\mathbf{c} + \mathbf{v}\_\mathbf{x} - \mathbf{e}x\_0) - b(\mathbf{z}\_\mathbf{c} + \mathbf{v}\_\mathbf{z} - \mathbf{e}\mathbf{\tilde{z}}) \right] \end{aligned} \tag{34}$$

with two distinct vectors **k**<sup>1</sup> and **k**<sup>2</sup> of Lagrange multipliers. Taking the partial derivatives with respect to all unknowns and setting them to zero results in the system of normal equations

$$\frac{1}{2}\frac{\partial K}{\partial \mathbf{v}\_{\mathbf{x}}^{\rm T}} = \mathbf{P}\_{\mathbf{x}\mathbf{x}}\mathbf{v}\_{\mathbf{x}} + \mathbf{P}\_{\mathbf{x}\mathbf{y}}\mathbf{v}\_{\mathbf{y}} + \mathbf{P}\_{\mathbf{x}\mathbf{z}}\mathbf{v}\_{\mathbf{z}} - \mathbf{k}\_{2} = \mathbf{0},\tag{35}$$

$$\frac{1}{2}\frac{\partial \bar{K}}{\partial \mathbf{v}\_{\mathbf{y}}^{\rm T}} = \mathbf{P}\_{\mathbf{yy}}\mathbf{v}\_{\mathbf{y}} + \mathbf{P}\_{\mathbf{yx}}\mathbf{v}\_{\mathbf{x}} + \mathbf{P}\_{\mathbf{yx}}\mathbf{v}\_{\mathbf{z}} - \mathbf{k}\_{1} = \mathbf{0},\tag{36}$$

$$\frac{1}{2}\frac{\partial \mathcal{K}}{\partial \mathbf{v}\_{\mathbf{z}}^{\mathrm{T}}} = \mathbf{P}\_{\mathrm{YY}}\mathbf{v}\_{\mathbf{y}} + \mathbf{P}\_{\mathrm{xx}}\mathbf{v}\_{\mathbf{x}} + \mathbf{P}\_{\mathrm{yz}}\mathbf{v}\_{\mathbf{y}} + a\mathbf{k}\_1 + b\mathbf{k}\_2 = \mathbf{0},\tag{37}$$

$$\frac{1}{2}\frac{\partial K}{\partial \mathbf{k}\_1^T} = (\mathbf{y}\_c + \mathbf{v}\_y - \mathbf{e}y\_0) - a(\mathbf{z}\_c + \mathbf{v}\_x - \mathbf{e}\mathbf{Z}) = \mathbf{0},\tag{38}$$

$$\frac{1}{2}\frac{\partial K}{\partial \mathbf{k}\_2^T} = (\mathbf{x}\_{\mathbf{c}} + \mathbf{v}\_{\mathbf{x}} - \mathbf{e}\mathbf{x}\_0) - a(\mathbf{z}\_{\mathbf{c}} + \mathbf{v}\_{\mathbf{z}} - \mathbf{e}\overline{z}) = \mathbf{0},\tag{39}$$

$$\frac{1}{2}\frac{\partial K}{\partial a} = \mathbf{k}\_1^\mathrm{T} \left(\mathbf{z}\_\mathrm{c} + \mathbf{v}\_\mathrm{z} - \mathbf{e}\Xi\right) = 0,\tag{40}$$

$$\frac{1}{2}\frac{\partial K}{\partial b} = \mathbf{k}\_2^T \left(\mathbf{z}\_c + \mathbf{v}\_x - \mathbf{e}\mathbf{\overline{z}}\right) = 0,\tag{41}$$

$$\frac{1}{2}\frac{\partial K}{\partial y\_0} = \mathbf{e}^{\mathsf{T}} \mathbf{k}\_{\mathsf{I}} = 0,\tag{42}$$

$$\frac{1}{2}\frac{\partial \mathbf{k}}{\partial \mathbf{x}\_0} = \mathbf{e}^{\mathsf{T}} \mathbf{k}\_2 = 0. \tag{43}$$

To derive explicit expressions for the residual vectors, we write Equations (35)–(37) as

$$
\begin{bmatrix}
\mathbf{P\_{xx}} & \mathbf{P\_{xy}} & \mathbf{P\_{xz}} \\
\mathbf{P\_{yx}} & \mathbf{P\_{yy}} & \mathbf{P\_{yz}} \\
\mathbf{P\_{zx}} & \mathbf{P\_{zy}} & \mathbf{P\_{zz}}
\end{bmatrix}
\begin{bmatrix}
\mathbf{v\_x} \\
\mathbf{v\_y} \\
\mathbf{v\_z}
\end{bmatrix} = 
\begin{bmatrix}
\mathbf{k\_2} \\
\mathbf{k\_1} \\
\end{bmatrix},
\tag{44}
$$

which leads to the solution for the residual vectors

$$
\begin{aligned}
\begin{bmatrix}
\mathbf{v}\_{\mathbf{x}} \\
\mathbf{v}\_{\mathbf{y}} \\
\mathbf{v}\_{\mathbf{z}}
\end{bmatrix} &= \begin{bmatrix}
\mathbf{P}\_{\mathbf{x}\mathbf{x}} & \mathbf{P}\_{\mathbf{xy}} & \mathbf{P}\_{\mathbf{x}\mathbf{z}} \\
\mathbf{P}\_{\mathbf{y}\mathbf{x}} & \mathbf{P}\_{\mathbf{yy}} & \mathbf{P}\_{\mathbf{yz}} \\
\mathbf{P}\_{\mathbf{z}\mathbf{x}} & \mathbf{P}\_{\mathbf{zy}} & \mathbf{P}\_{\mathbf{z}\mathbf{z}}
\end{bmatrix}^{-1} \begin{bmatrix}
\mathbf{k}\_{2} \\
\mathbf{k}\_{1} \\
\end{bmatrix} \\
&= \begin{bmatrix}
\mathbf{Q}\_{\mathbf{x}\mathbf{x}} & \mathbf{Q}\_{\mathbf{xy}} & \mathbf{Q}\_{\mathbf{xz}} \\
\mathbf{Q}\_{\mathbf{yx}} & \mathbf{Q}\_{\mathbf{yy}} & \mathbf{Q}\_{\mathbf{yz}} \\
\mathbf{Q}\_{\mathbf{xx}} & \mathbf{Q}\_{\mathbf{xy}} & \mathbf{Q}\_{\mathbf{zz}}
\end{bmatrix} \begin{bmatrix}
\mathbf{k}\_{2} \\
\mathbf{k}\_{1} \\
\end{bmatrix},
\end{aligned}
\tag{45}$$

or explicitly

$$\mathbf{v}\_{\mathbf{x}} = \left(\mathbf{Q}\_{\mathbf{x}\mathbf{y}} - a\mathbf{Q}\_{\mathbf{x}\mathbf{z}}\right)\mathbf{k}\_1 + \left(\mathbf{Q}\_{\mathbf{x}\mathbf{x}} - b\mathbf{Q}\_{\mathbf{x}\mathbf{z}}\right)\mathbf{k}\_{\mathbf{2}}.\tag{46}$$

$$\mathbf{v}\_{\rm Y} = \left(\mathbf{Q}\_{\rm yy} - a\mathbf{Q}\_{\rm yz}\right)\mathbf{k}\_1 + \left(\mathbf{Q}\_{\rm xy} - b\mathbf{Q}\_{\rm yz}\right)\mathbf{k}\_{2\prime} \tag{47}$$

$$\mathbf{v}\_{\mathbf{z}} = \left(\mathbf{Q}\_{\text{yz}} - a\mathbf{Q}\_{\text{z}}\right)\mathbf{k}\_1 + \left(\mathbf{Q}\_{\text{xz}} - b\mathbf{Q}\_{\text{xz}}\right)\mathbf{k}\_2. \tag{48}$$

Inserting these expressions for the residual vectors into Equations (38) and (39) yields

$$\mathbf{W}\_1 \mathbf{k}\_1 + \mathbf{W}\_2 \mathbf{k}\_2 + \mathbf{y}\_\mathbf{c} - \mathbf{e} y\_0 - a \mathbf{z}\_\mathbf{c} = \mathbf{0},\tag{49}$$

$$\mathbf{W}\_2 \mathbf{k}\_1 + \mathbf{W}\_3 \mathbf{k}\_2 + \mathbf{x}\_\mathbf{c} - \mathbf{e} \mathbf{x}\_0 - b \mathbf{z}\_\mathbf{c} = \mathbf{0},\tag{50}$$

with the auxiliary matrices **W**1, **W**<sup>2</sup> and **W**<sup>3</sup> as

$$\mathbf{W}\_1 = \mathbf{Q}\_\mathbf{Y} - 2a\mathbf{Q}\_\mathbf{Y} + a^2 \mathbf{Q}\_\mathbf{Z} \tag{51}$$

$$\mathbf{W}\_2 = \mathbf{Q}\_{\mathbf{Y}\mathbf{y}} - a\mathbf{Q}\_{\mathbf{X}\mathbf{z}} - b\mathbf{Q}\_{\mathbf{Y}\mathbf{z}} + ab\mathbf{Q}\_{\mathbf{z}\mathbf{z}} \tag{52}$$

$$\mathbf{W}\mathbf{\hat{s}} = \mathbf{Q}\_{\text{X}} - 2b\mathbf{Q}\_{\text{XZ}} + b^2 \mathbf{Q}\_{\text{Z}} \tag{53}$$

Using Equations (49), (50) and (40)–(43), the nonlinear equation system

$$\begin{aligned} \mathbf{W}\_1 \mathbf{k}\_1 + \mathbf{W}\_2 \mathbf{k}\_2 - \mathbf{e}\_3 \mathbf{e}\_0 - a \mathbf{z}\_c &= -\mathbf{y}\_c, \\\\ \mathbf{W}\_2 \mathbf{k}\_1 + \mathbf{W}\_3 \mathbf{k}\_2 - \mathbf{e} \mathbf{x}\_0 - b \mathbf{z}\_c &= -\mathbf{x}\_c, \\\\ -\left(\mathbf{z}\_c + \mathbf{v}\_z - \mathbf{e} \overline{\mathbf{z}}\right)^\mathrm{T} \mathbf{k}\_1 &= 0, \\\\ -\left(\mathbf{z}\_c + \mathbf{v}\_z - \mathbf{e} \overline{\mathbf{z}}\right)^\mathrm{T} \mathbf{k}\_2 &= 0, \\\\ -\mathbf{e}^\mathrm{T} \mathbf{k}\_1 &= 0, \\\\ -\mathbf{e}^\mathrm{T} \mathbf{k}\_2 &= 0 \end{aligned} \tag{54}$$

can be set up. To express this system of equations into a block matrix form and to obtain a symmetrical matrix of normal equations in the following, it is advantageous to add the terms −*a* (**v**<sup>z</sup> − **e***z*) and −*b* (**v**<sup>z</sup> − **e***z*) to both sides of the first two equations, leading to the system of equations

$$\begin{aligned} \mathbf{W\_1k\_1} + \mathbf{W\_2k\_2} - \mathbf{e}\_y \boldsymbol{\rho}\_0 - a \left( \mathbf{z\_c} + \mathbf{v\_z} - \mathbf{e}\widetilde{\boldsymbol{\Xi}} \right) &= -\mathbf{y\_c} - a \left( \mathbf{v\_z} - \mathbf{e}\widetilde{\boldsymbol{\Xi}} \right), \\\\ \mathbf{W\_2k\_1} + \mathbf{W\_3k\_2} - \mathbf{e}\mathbf{x\_0} - b \left( \mathbf{z\_c} + \mathbf{v\_z} - \mathbf{e}\widetilde{\boldsymbol{\Xi}} \right) &= -\mathbf{x\_c} - b \left( \mathbf{v\_z} - \mathbf{e}\widetilde{\boldsymbol{\Xi}} \right), \\\\ &- \left( \mathbf{z\_c} + \mathbf{v\_z} - \mathbf{e}\widetilde{\boldsymbol{\Xi}} \right)^T \mathbf{k\_1} = 0, \\\\ &- \left( \mathbf{z\_c} + \mathbf{v\_z} - \mathbf{e}\widetilde{\boldsymbol{\Xi}} \right)^T \mathbf{k\_2} = 0, \\\\ &- \mathbf{e}^T \mathbf{k\_1} = 0, \\\\ &- \mathbf{e}^T \mathbf{k\_2} = 0. \end{aligned} \tag{55}$$

Arranging the unknowns in a vector

$$\mathbf{X} = \begin{bmatrix} \mathbf{k}\_1 \\ \mathbf{k}\_2 \\ a \\ b \\ y\_0 \\ \mathbf{x}\_0 \end{bmatrix},\tag{56}$$

the equation system (55) can be written as

$$\mathbf{N}\mathbf{X} = \mathbf{n}.\tag{57}$$

with the matrix of normal equations

$$\mathbf{N} = \begin{bmatrix} \mathbf{W} & \mathbf{A} \\ \mathbf{A}^T & \mathbf{0} \end{bmatrix},\tag{58}$$

constructed using the matrices

$$\mathbf{W} = \begin{bmatrix} \mathbf{W}\_1 & \mathbf{W}\_2 \\ \mathbf{W}\_2 & \mathbf{W}\_3 \end{bmatrix} \tag{59}$$

and

$$\mathbf{A}^{\mathsf{T}} = \begin{bmatrix} -\left(\mathbf{z}\_{\mathsf{c}} + \mathbf{v}\_{\mathsf{z}} - \mathbf{e}\,\overline{\mathbf{z}}\right) & \mathbf{0} \\ \mathbf{0} & -\left(\mathbf{z}\_{\mathsf{c}} + \mathbf{v}\_{\mathsf{z}} - \mathbf{e}\,\overline{\mathbf{z}}\right) \\ -\mathbf{e} & \mathbf{0} \\ \mathbf{0} & -\mathbf{e} \end{bmatrix}.\tag{60}$$

The right-hand side of the equation system (57) reads

$$\mathbf{n} = \begin{bmatrix} -\mathbf{y}\_{\text{c}} - a\left(\mathbf{z}\_{\text{c}} + \mathbf{v}\_{\text{z}} - \mathbf{e}\tilde{\mathbf{z}}\right) \\ -\mathbf{x}\_{\text{c}} - b\left(\mathbf{z}\_{\text{c}} + \mathbf{v}\_{\text{z}} - \mathbf{e}\tilde{\mathbf{z}}\right) \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} \tag{61}$$

It is important to point out that matrix **N** and vector **n** contain the unknown parameters *a*, *b* and **v**z in their entries. Therefore, to express and solve these normal equations that live in the "nonlinear world" with the help of vectors and matrices (that only exist in the "linear world"), appropriate approximate values *a*0, *b*<sup>0</sup> and **v**<sup>0</sup> <sup>z</sup> have to be introduced for all auxiliary matrices required for the setup of matrix **N** and vector **n**. The solution for the vector of unknowns can be computed by

$$
\hat{\mathbf{X}} = \mathbf{N}^{-1} \mathbf{n},
\tag{62}
$$

without the need of a linearization by Taylor series at any step of the calculation process. The WTLS solution for the line parameters can be computed following the *ordinary iteration method*, as it is explained for example by Bronhstein ([24], p. 896). Therefore, the solutions *a*ˆ, ˆ *b* and **v**˜ *z*, after stripping them of their random character, are to be substituted as new approximate values as long as necessary until a sensibly chosen break-off condition is met. For example, the maximum absolute difference between two consecutive solutions must be smaller than a predefined threshold , which in this problem can be formulated as

$$\max \left| \left\lceil \begin{array}{c} a^0 \\ b^0 \end{array} \right\rceil - \left\lceil \begin{array}{c} \mathfrak{a} \\ \mathfrak{b} \end{array} \right\rceil \right| \leq \mathfrak{c}\_{\prime} \tag{63}$$

with |·| denoting the absolute value. The predicted residual vectors **v**˜ x, **v**˜ y, **v**˜ <sup>z</sup> can be computed from Equations (46)–(48). The iterative procedure for the presented WTLS solution can be found in Algorithm 1.

#### **Algorithm 1** Iterative WTLS solution

Choose approximate values for *a*0, *b*<sup>0</sup> and **v**<sup>0</sup> z . Define parameter *c* = 1. Set threshold for the break-off condition of the iteration process. Set parameter *da* = *db* = ∞, for entering the iteration process. **while** *da* > and *db* > **do** Compute matrices **W** and **A**. Build matrix **N** and vector **n**. Estimate the vector of unknowns **X**ˆ . Compute the residual vector ˜**v**z . Compute parameters *da* <sup>=</sup> <sup>|</sup>*a*<sup>ˆ</sup> <sup>−</sup> *<sup>a</sup>*0<sup>|</sup> and *db* <sup>=</sup> <sup>|</sup><sup>ˆ</sup> *<sup>b</sup>* − *<sup>b</sup>*0|. Update the approximate values with the estimated ones, with *a*<sup>0</sup> = *a*ˆ, *b*<sup>0</sup> = ˆ *b* and **v**<sup>0</sup> <sup>z</sup> = **v**˜ z. **end while return** *a*ˆ and ˆ *b*, with *c* = 1.

After computing the line parameters *a*ˆ and ˆ *b* and putting them into a vector, we can easily scale it into a unit vector by dividing each component with the length of the vector

$$
\begin{bmatrix} a\_{\mathbf{n}} \\ b\_{\mathbf{n}} \\ c\_{\mathbf{n}} \end{bmatrix} = \frac{\begin{bmatrix} \hat{a} \\ \hat{b} \\ 1 \end{bmatrix}}{\begin{bmatrix} 1 \\ \hat{a} \\ \hat{b} \\ 1 \end{bmatrix}} \tag{64}
$$

with ||·|| denoting the Euclidean norm. The derived parameters *a*n, *b*<sup>n</sup> and *c*<sup>n</sup> refer to the normalized components of a unit vector, that is parallel to the requested line, with

$$a\_n^2 + b\_n^2 + c\_n^2 = 1\tag{65}$$

#### *2.4. WTLS Solution with Singular Dispersion Matrices*

The algorithmic approach presented in Section 2.3 can also cover cases when dispersion matrices are singular. Such a solution depends on the inversion of matrix **N** (58), which depends on the rank deficiency of matrix **W** (59). Following the argumentation of Malissiovas ([20], p. 112), a criterion that ensures a unique solution of the problem can be described in this case by

$$\text{rank}\left(\left[\mathbf{W} \mid \mathbf{A}\right]\right) = 2n,\tag{66}$$

with



In cases of singular dispersion matrices, the rank of matrix **W** will be smaller than 2*n*. A unique solution will exist if the rank of the augmented matrix [**W** | **A**] is equal to the number of condition equations 2*n* of the problem. It is important to mention that the developed idea is based on the Neitzel–Schaffrin criterion, which has been firstly proposed by Neitzel and Schaffrin [26,27], particularly for a solution of an adjustment problem within the GH model when singular dispersion matrices must be employed.

#### *2.5. A Posteriori Error Estimation*

In this section we want to determine the variance-covariance matrix of the estimated parameters. The following derivations can be used for computing the a posteriori stochastic information for all weighted cases discussed in this investigation, i.e., the direct and the iterative WTLS solutions. Therefore, we will employ the fundamental idea of variance-covariance propagation. This is a standard procedure explained by many authors, like for example in the textbooks of Wells and Krakiwsky ([28], p. 20), Mikhail ([29], p. 76 ff.) or Niemeier ([8], p. 51 ff.) and has been further employed in the GM and the GH model, so that the a posteriori stochastic results can be computed using directly the developed matrices from each model. A detailed explanation is given by Malissiovas ([20], p. 31 ff.).

As we have already mentioned in the introduction of this article, a TLS, respectively a WTLS solution, can be regarded as a special case of a least squares solution within the GH model. From the presented WTLS algorithm we observe that the derived matrix of normal Equation (58) is equal to the matrix if it was computed within the GH model. Therefore, it is possible to compute

$$\mathbf{N}^{-1} = \begin{bmatrix} \mathbf{Q}\_{11} & \mathbf{Q}\_{12} \\ \mathbf{Q}\_{21} & \mathbf{Q}\_{22} \end{bmatrix} \tag{67}$$

and extract the dispersion matrix for the unknown parameters from

$$
\mathbf{Q}\_{\\$\&} = -\mathbf{Q}\_{22}.\tag{68}
$$

The a posteriori variance factor is

$$
\hat{\sigma}\_0^2 = \frac{\mathbf{v}^T \mathbf{P} \mathbf{v}}{r} \tag{69}
$$

with vector **v** holding all residuals and *r* denoting the redundancy of the adjustment problem. In case of a singular dispersion matrix, it is not possible to compute the weight matrix **P**, as in Equation (27). Therefore, we can make use of the solution for the residual vectors from Equation (45) and insert them in Equation (69) to obtain

$$\mathbf{p}\_0^2 = \frac{\mathbf{v}\_\mathbf{x}^\mathrm{T}\mathbf{k}\_2 + \mathbf{v}\_\mathbf{y}^\mathrm{T}\mathbf{k}\_1 - a\mathbf{v}\_\mathbf{z}^\mathrm{T}\mathbf{k}\_1 - b\mathbf{v}\_\mathbf{z}^\mathrm{T}\mathbf{k}\_2}{r}.\tag{70}$$

The variance-covariance matrix of the unknown parameters can be then derived by

$$
\Sigma\_{\mathbb{k}\mathbb{k}} = \hat{\sigma}\_0^2 \mathbf{Q}\_{\mathbb{k}\mathbb{k}} = \begin{bmatrix}
\sigma\_a^2 & \sigma\_{ab} & \sigma\_{ay\_0} & \sigma\_{ax\_0} \\
\sigma\_{ba} & \sigma\_b^2 & \sigma\_{by\_0} & \sigma\_{bx\_0} \\
\sigma\_{ya} & \sigma\_{yb}b & \sigma\_{yb}^2 & \sigma\_{y\_0x\_0} \\
\sigma\_{x0a} & \sigma\_{x0b} & \sigma\_{x0y\_0} & \sigma\_{x0}^2
\end{bmatrix} . \tag{71}$$

To derive the variance-covariance matrix of the normalized vector components (64), we can explicitly write the equations

$$a\_{\rm n} = \frac{\hat{a}}{\sqrt{\hat{a}^2 + \hat{b}^2 + c^2}},$$

$$b\_{\rm n} = \frac{\hat{b}}{\sqrt{\hat{a}^2 + \hat{b}^2 + c^2}},\tag{72}$$

$$c\_{\rm n} = \frac{c}{\sqrt{\hat{a}^2 + \hat{b}^2 + c^2}}.$$

with *c* = 1. Following the standard procedure of the variance-covariance propagation in nonlinear cases, we can write the Jacobian matrix

$$\mathbf{F} = \begin{bmatrix} \frac{\partial a\_{\mathbf{n}}}{\partial a} & \frac{\partial a\_{\mathbf{n}}}{\partial b} \\\\ \frac{\partial b\_{\mathbf{n}}}{\partial a} & \frac{\partial b\_{\mathbf{n}}}{\partial b} \\\\ \frac{\partial c\_{\mathbf{n}}}{\partial a} & \frac{\partial c\_{\mathbf{n}}}{\partial b} \end{bmatrix}. \tag{73}$$

Taking into account the variances and covariances of the line parameters *a*ˆ and ˆ *b* from (71)

$$
\Sigma\_{\hat{\mathbb{A}}\hat{\mathbb{b}}} = \Sigma\_{\mathbb{A}\mathbb{R}}(1:2, 1:2) = \left[ \begin{array}{cc} \sigma\_{a}^{2} & \sigma\_{ab} \\ \sigma\_{ba} & \sigma\_{b}^{2} \end{array} \right], \tag{74}
$$

we can compute the variance-covariance matrix of the normalized components

$$\boldsymbol{\Sigma}\_{\mathbf{a}\_{\mathsf{h}}\mathbf{b}\_{\mathsf{h}}\mathbf{c}\_{\mathsf{h}}} = \mathbf{F}\boldsymbol{\Sigma}\_{\mathbf{A}\mathbf{\hat{b}}}\mathbf{F}^{\mathrm{T}} = \begin{bmatrix} \sigma\_{a\_{n}}^{2} & \sigma\_{a\_{n}b\_{n}} & \sigma\_{a\_{n}c\_{n}}\\ \sigma\_{b\_{n}a\_{n}} & \sigma\_{b\_{n}}^{2} & \sigma\_{b\_{n}c\_{n}}\\ \sigma\_{c\_{n}a\_{n}} & \sigma\_{c\_{n}b\_{n}} & \sigma\_{c\_{n}}^{2} \end{bmatrix} . \tag{75}$$

#### **3. Numerical Examples**

In this section we present the solutions for fitting a straight line to 3D point data using the TLS approach from Section 2.1 and the newly developed WTLS approaches from Sections 2.2 and 2.3. The utilized dataset for this investigaiton consists of *n* = 50 points and originates from the work of Petras and Podlubny [30]. It has been utilized also by Snow and Schaffrin, see Table A1 in [6], for a solution within the GH model, which will be used here for validating the results of the presented WTLS solutions. Three different stochastic models will be imposed in the following:


#### *3.1. Equal Weights*

For the first case under investigation, we consider all coordinate components *xi*, *yi*, *zi* as equally weighted and uncorrelated observations, yielding the weights as shown in (3). The least squares solution of this problem within the GH model, presented by Snow and Schaffrin [6], is listed in Table 1.


**Table 1.** Least squares solution within the Gauss–Helmert (GH) model of Snow and Schaffrin [6].

A direct TLS solution for this problem can be derived using the approach presented in Section 2.1. The results are shown in Table 2. Numerically equal results have been derived by using the direct WTLS approach of Section 2.2 by setting all weights equal to one.


**Table 2.** Direct TLS solution from Section 2.1.

Comparing the solution with the one presented in Table 1, it can be concluded that the numerical results for the parameters coincide within the specified decimal places. Regarding the small difference in the numerical value for the variance factor, it is to be noted that the value in Table 2 was confirmed by two independent computations.

Furthermore, a point on the line can be easily computed using the equations of the functional model (2), as long as the direction vector parallel to the requested line is known. Alternatively, all the adjusted points will lie on the requested straight line, which can be simply computed by adding the computed residuals to the measured coordinates.

#### *3.2. Pointwise Weights*

For the second weighted case under investigation, we consider the coordinate components *xi*, *yi*, *zi* of each point *Pi* to be uncorrelated and of equal precision. From the standard deviations listed in Table 3, the corresponding pointwise weights can be obtained from (11).


**Table 3.** Pointwise precision *σxi* = *σyi* = *σzi* = *σ<sup>i</sup>* for each point *Pi*.

A direct WTLS solution is derived, following the approach presented in Section 2.2. The determinant (23)

$$
\begin{vmatrix}
(d\_1 - k) & d\_4 & d\_5 \\
 & d\_4 & (d\_2 - k) & d\_6 \\
 & d\_5 & d\_6 & (d\_3 - k)
\end{vmatrix} = 0,
$$

can be built with the components

$$\begin{aligned} d\_1 &= 213.505250528675, \\ d\_2 &= 206.458905097029, \\ d\_3 &= 198.273122927545, \\ \text{and} \\ d\_4 &= -20.9837621443375, \\ d\_5 &= -12.1835697465792, \\ d\_6 &= -81.6787394185243. \end{aligned} \tag{76}$$

which leads to the solutions for the Lagrange multiplier

$$\hat{k} = \begin{cases} 115.0596477492 & (\hat{k}\_{\text{min}}) \\ 218.3490470615 & \\ 284.8285837425 & \end{cases} \tag{77}$$

The direct WTLS solution for the line orientation components is shown in Table 4.


**Table 4.** Direct WTLS solution from Section 2.2.

The presented results are numerically equal to the iterative WTLS solution using the algorithmic approach of Section 2.3, as well as the solution within the GH model.

#### *3.3. General Weights*

For the last weighted case in this investigation, we impose the most general case, i.e., correlated coordinate components with individual precision resulting in a singular dispersion matrix. To obtain such a matrix for our numerical investigations, we firstly solved the adjustment problem within the GH model with an identity matrix as dispersion matrix. From the resulting 150 × 150 dispersion matrix of the residuals

$$\mathbf{Q}\_{\rm{VV}} = \begin{bmatrix} \mathbf{Q}\_{\rm{V}\_{\rm{X}}\rm{V}\_{\rm{X}}} & \mathbf{Q}\_{\rm{V}\_{\rm{X}}\rm{V}\_{\rm{Y}}} & \mathbf{Q}\_{\rm{V}\_{\rm{X}}\rm{V}\_{\rm{X}}} \\ \mathbf{Q}\_{\rm{V}\_{\rm{Y}}\rm{V}\_{\rm{X}}} & \mathbf{Q}\_{\rm{V}\_{\rm{Y}}\rm{V}\_{\rm{Y}}} & \mathbf{Q}\_{\rm{V}\_{\rm{Y}}\rm{V}\_{\rm{X}}} \\ \mathbf{Q}\_{\rm{V}\_{\rm{X}}\rm{V}\_{\rm{X}}} & \mathbf{Q}\_{\rm{V}\_{\rm{X}}\rm{V}\_{\rm{Y}}} & \mathbf{Q}\_{\rm{V}\_{\rm{X}}\rm{V}\_{\rm{X}}} \end{bmatrix},\tag{78}$$

computed as e.g., presented by Malissiovas ([20], p. 46), we take the variances and covariances between the individual point coordinates, but not among the points, i.e., the diagonal elements of each sub-matrix in (78) and arrange them in a new 150 × 150 matrix

$$\mathbf{Q}\_{\rm LL} = \begin{bmatrix} \text{diag}\left(\mathbf{Q}\_{\rm v\_{\rm x}\rm v\_{\rm x}}\right) & \text{diag}\left(\mathbf{Q}\_{\rm v\_{\rm x}\rm v\_{\rm y}}\right) & \text{diag}\left(\mathbf{Q}\_{\rm v\_{\rm x}\rm v\_{\rm x}}\right) \\ \text{diag}\left(\mathbf{Q}\_{\rm v\_{\rm y}\rm v\_{\rm x}}\right) & \text{diag}\left(\mathbf{Q}\_{\rm v\_{\rm y}\rm v\_{\rm y}}\right) & \text{diag}\left(\mathbf{Q}\_{\rm v\_{\rm y}\rm v\_{\rm x}}\right) \\ \text{diag}\left(\mathbf{Q}\_{\rm v\_{\rm x}\rm v\_{\rm x}}\right) & \text{diag}\left(\mathbf{Q}\_{\rm v\_{\rm x}\rm v\_{\rm y}}\right) & \text{diag}\left(\mathbf{Q}\_{\rm v\_{\rm x}\rm v\_{\rm x}}\right) \end{bmatrix},\tag{79}$$

with "diag()" denoting that only the diagonal elements are taken into account. This is an example of pointwise variances and covariances, described as case (i) in Section 2.3, but now yielding a singular dispersion matrix for the observations with

$$\text{rank}\,(\mathbf{Q}\_{\text{LL}}) = 100\text{.}$$

Before deriving an iterative WTLS solution for this weighted case, we must check if the criterion (66) for a unique solution of the adjustment problem is fulfilled. Therefore, we computed the 100 × 100 matrix **W** with

$$\text{rank}\,(\mathsf{W}) = 100\_{\mathsf{M}}$$

and the 100 × 4 matrix **A** with

rank (**A**) = 4.

The criterion ensures that a unique solution exists when using the presented singular dispersion matrix, while

$$\text{rank}\left(\left[\mathbf{W} \mid \mathbf{A}\right]\right) = 100,\tag{80}$$

since *n* = 50 observed points are used in this example, cf. (66). As for all iterative procedures, appropriate starting values for the unknowns must be provided. However, they can be obtained easily by first generating a direct solution with a simplified stochastic model. The iterative WTLS solution for the direction vector of the requested straight line is presented in Table 5.


**Table 5.** Iterative WTLS solution from Section 2.3.

The presented WTLS solution has been found to be numerically equal to the least squares solution within the GH model. Detailed numerical investigations of the convergence behaviour, e.g., in comparison to an adjustment within the GH model, are beyond the scope of this article. However, in many numerical examples it could be observed that the iterative WTLS approach showed a faster convergence rate compared to an adjustment within an iteratively linearized GH model.

#### **4. Conclusions**

For the problem of straight line fitting to 3D point data, two novel WTLS algorithms for two individual weighting schemes have been presented in this study:


Both approaches are based on the work of Malissiovas [20], where similar algorithms have been presented for adjustment problems that belong to the same class, i.e., nonlinear adjustments that can be expressed within the EIV model. The approach presented in Section 2.1 provides a direct TLS solution assuming equally weighted and uncorrelated coordinate components. The fact that this assumption is inappropriate, e.g. for the evaluation of laser scanning data, has often been accepted in the past to provide a direct solution for large data sets. With the newly developed approach in Section 2.2 it is now possible to compute a direct WTLS solution at least for a more realistic stochastic model, namely pointwise weighting schemes.

If more general weight matrices must be taken into account in the stochastic model, including correlations or singular dispersion matrices, the presented algorithm of Section 2.3 can be utilized for an iterative solution without linearizing the problem by Taylor series at any step, following the algorithmic idea of WTLS. A criterion that ensures a unique solution of the problem when employing singular dispersion matrices has also been presented, which is based on the original ideas of Neitzel and Schaffrin [26,27], for a solution within the GH model.

Numerical examples have been presented in Section 3 for testing the presented WTLS algorithms. The utilized dataset of the observed 3D point data has also been employed by Snow and Schaffrin [6] for a solution of the problem within the GH model and originates from the study of Petras and Podlubny [30]. The results of the presented algorithms have been compared in all cases with existing solutions or the solutions coming from existing algorithms and have been found to be numerically equal.

**Author Contributions:** Conceptualization, G.M., F.N. and S.P.; methodology, G.M. and S.P.; software, G.M. and S.W.; validation, F.N. and S.P.; formal analysis, G.M.; investigation, G.M. and S.W.; data curation, G.M. and S.W.; writing—original draft preparation, G.M. and F.N.; writing—review and editing, G.M., F.N., S.P. and S.W.; supervision, F.N. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors acknowledge support by the German Research Foundation and the Open Access Publication Fund of TU Berlin.

**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/).
