Spatial Coordinate Transformations with Noisy Data

Abstract

The parametric transformation of spatial coordinates between different reference frames is traditionally dealt with a stepwise approach which imposes a suboptimal treatment in the presence of noisy data. The chapter explains briefly the drawbacks of this approach and then presents an alternative scheme for spatial coordinate transformations that improves the classic stepwise solution when using noisy coordinates of known stochastic structure. The proposed methodology is simple in principle, although its numerical implementation with nonlinear parametric models is a bit more involved and it relies on the joint least squares adjustment of the observed coordinates using their full stochastic model over all points of interest. The mathematical framework and the related properties of this “stacking” approach are presented in detail, along with a numerical example that demonstrates its feasibility for practical problems in geospatial applications.

1. Introduction

Spatial coordinate-based positions often need to be transformed from their own reference frame to another reference frame by an analytic parametric model. This is a standard problem in several fields of geosciences and engineering, including geodesy [11, 18, 19, 30, 33], land surveying and cadastral planning [2, 8, 16, 32], cartography and digital mapping [1, 5, 28, 38], photogrammetry and remote sensing [22, 25, 39], robotics and computer vision [7, 13, 31], among others.

Spatial coordinate transformations (SCTs) are utilized in practice either explicitly for determining the unknown coordinates of scattered points in a desired frame from their observed coordinates in a different frame or implicitly in the context of more composite procedures such as the self-calibration of terrestrial laser scanners [20, 23, 26], the conflation of digital maps and geographical databases [2, 4, 6, 27], the reconstruction of 3D models from multi-sensor data [9, 29, 35], and the integration of aerial or satellite images in ground-based systems of geographic coordinates [14, 15, 21]. Various technical terms have actually been used with regard to SCTs in practical problems, for instance, spatial adjustment, image registration, absolute orientation, geo-referencing, and frame transformation, to name a few. Despite their linkage to different application fields, all these terms refer, more or less, to the same archetypical problem, that is, the optimal fusion of partially overlapping configurations of spatial points using their coordinate-based representations in separate frames and an application-specific model to describe their systematic differences. Although this general viewpoint includes also cases with raster-type and/or vector-type data in multiple frames, the present chapter is confined to the vectorized version of SCT problems between two different frames.

2. Problem statement and research motivation

The transformation of spatial coordinates between different reference frames using a parametric model is traditionally resolved via a two-step procedure. Firstly, the parameters of the adopted model are estimated by least squares adjustment (LSA) techniques on the basis of control points with known coordinates in both frames of interest. The estimated parameters are then used to transform any set of known coordinates from the original frame (also called source or initial frame) to the desired target frame. The second step is applied not only to control points but also to additional points whose spatial coordinates are originally known only in the initial frame. Their transformed coordinates are often the primary objective in many practical applications, whereas the post-fit residuals at control points are commonly used as quality metrics of the transformation process.

Despite its rational character, the stepwise methodology imposes a suboptimal treatment in the presence of noisy data. The reason is that the accuracy of the original coordinates in the initial frame is ignored during their transformation to the target frame. Indeed, the second step entails only the forward implementation of the parametric model without attempting to minimize or, at least, reduce the propagated random errors of the original coordinates. As shown in [18], the stepwise procedure may actually enhance the data noise into the final results, in the sense that the transformed coordinates in the target frame could become less accurate than the original coordinates in the initial frame.

An additional weakness of the traditional stepwise procedure is that the control and non-control points (abbreviated hereafter as CPs and NPs) are handled independently throughout the transformation process. Yet, the initial coordinates of these points are usually acquired from the same pre-analysis phase or observational procedure (e.g., geodetic network adjustment, map digitization, etc.), and therefore they are affected by common error sources. Due to their separate treatment, the cross-correlated part of the initial coordinate errors at CPs and NPs is always ignored, a fact that may weaken the accuracy of the transformed coordinates in the target frame.

The focus of the present research is the formulation of a single-stage estimation scheme that can improve the classic stepwise solution in SCT problems. The proposed scheme is based on a properly weighted least squares adjustment of all observed coordinates, using their known variances and covariances (CV) in the entire group of transformation points. This stacking approach permits the rigorous treatment of intra-frame error correlations among CPs and NPs, and it generally leads to higher-accuracy results for the transformed coordinates. Our contribution provides easy-to-use optimal estimators for the transformed coordinates under any parametric model, regardless of the structure of the error covariance matrices of the input data. A numerical example is also given at the end of the chapter to demonstrate the feasibility of the proposed methodology for practical applications.

3. The stepwise approach in spatial coordinate transformations

3.1 Preamble

The mathematical setting is based on the general case of nonlinear transformation models, in accordance to the vectorized expression:

$$\mathit{X}=\mathit{f}\left({\mathit{X}}^{\prime },\mathit{\theta }\right)$$E1

where ${\mathit{X}}^{\prime }$ and $\mathit{X}$ contain the Cartesian coordinates for the same group of physical points (or homologous groups of points) in different reference frames. The spatial dimensionality does not need to be specified here, and it can refer to any case that occurs in practice. The vector $\mathit{\theta }$ represents the parameters of the transformation model which enable the coordinate mapping from the initial frame to the target frame.

In the following, Eq. (1) is considered as an exact formula for noise-free coordinates and provides the general framework for the LSA of observed coordinates in the involved frames. Simpler types of transformation models with joint or marginal linearity in ${\mathit{X}}^{\prime }$ and $\mathit{\theta }$ (e.g., errors-in-variables models, differential or close-to-identity models) can be also analyzed under the previous setting.

For the purpose of this contribution, the user’s data shall consist of:

1. the observed coordinates for CPs and NPs in the initial frame (denoted by ${\mathit{X}}^{\prime }$ and ${\mathit{Z}}^{\prime }$, respectively);

2. the observed coordinates for CPs in the target frame (denoted by $\mathit{X}$); and

3. the error CV matrices of the previous vectors (denoted by ${\mathit{\Sigma }}_{{\mathit{X}}^{\prime }}$, ${\mathit{\Sigma }}_{{\mathit{Z}}^{\prime }}$, ${\mathit{\Sigma }}_{\mathit{X}}$).

An additional matrix of special importance is the cross-CV matrix ${\mathit{\Sigma }}_{{\mathit{X}}^{\prime }{\mathit{Z}}^{\prime }}$ which reflects the intra-frame error correlation between CPs and NPs, and it is totally ignored in the traditional stepwise procedure.

3.2 Estimation of transformation parameters

The first step refers to the estimation of the transformation parameters using a sufficient number of known CPs. Following a statistical estimation perspective, the optimal parameter values are obtained by solving the nonlinear LSA problem

$$\underset{\mathit{\theta }}{min}\left({\mathit{v}}_{\mathit{X}}^T{\mathit{\Sigma }}_{\mathit{X}}^{-1}{\mathit{v}}_{\mathit{X}}+{\mathit{v}}_{{\mathit{X}}^{\prime }}^T{\mathit{\Sigma }}_{{\mathit{X}}^{\prime }}^{-1}{\mathit{v}}_{{\mathit{X}}^{\prime }}\right)$$E2

subject to

$$\mathit{X}+{\mathit{v}}_{\mathit{X}}=\mathit{f}\left({\mathit{X}}^{\prime }+{\mathit{v}}_{{\mathit{X}}^{\prime }},\mathit{\theta }\right)$$E3

where the vectors ${\mathit{v}}_{\mathit{X}}$ and ${\mathit{v}}_{{\mathit{X}}^{\prime }}$ represent the zero-mean random errors in the observed coordinates. After appropriate linearization, the above problem can be reduced to a linear LSA problem for the so-called Gauss-Helmert (GH) model [12, 17], and it leads to an iterative solution via successive refinements of the preliminary estimate:

$$\widehat{\mathit{\theta }}={\mathit{\theta }}_o+{\left({\mathit{J}}_{\mathit{\theta }}^T{\mathit{WJ}}_{\mathit{\theta }}\right)}^{-1}{\mathit{J}}_{\mathit{\theta }}^T\mathit{W}\left(\mathit{X}-\mathit{f}\left({\mathit{X}}^{\prime },{\mathit{\theta }}_o\right)\right)$$E4

$$\mathit{W}={\left({\mathit{\Sigma }}_{\mathit{X}}+{\mathit{J}}_{{\mathit{X}}^{\prime }}{\mathit{\Sigma }}_{{\mathit{X}}^{\prime }}{\mathit{J}}_{{\mathit{X}}^{\prime }}^T\right)}^{-1}$$E5

where ${\mathit{\theta }}_o$ contains approximate values for the transformation parameters. The recursive updating of the previous solution is performed by the Newton-Gauss iteration algorithm in accordance to a more complex expression that will be presented later in this chapter. The matrices ${\mathit{J}}_{{\mathit{X}}^{\prime }}$ and ${\mathit{J}}_{\mathit{\theta }}$ are the Jacobians with respect to the initial frame coordinates and the transformation parameters, that is,

$${\mathit{J}}_{{\mathit{X}}^{\prime }}=\frac{\partial \mathit{f}\left({\mathit{X}}^{\prime },\mathit{\theta }\right)}{\partial {\mathit{X}}^{\prime }},{\mathit{J}}_{\mathit{\theta }}=\frac{\partial \mathit{f}\left({\mathit{X}}^{\prime },\mathit{\theta }\right)}{\partial \mathit{\theta }}$$E6

and they need to be re-evaluated at each iteration step using the adjusted values from the previous step. For more details on nonlinear least squares adjustment and iterative computational algorithms, the reader should consult the excellent treatise in [24] (see also [3, 10, 34]).

3.3 Determination of transformed coordinates

After estimating the transformation parameters, an additional step is required to complete the solution of the problem at hand, that is, the computation of the transformed coordinates in the target frame. This is performed by a simple forward evaluation of the transformation model at the CPs and NPs, using the respective nonlinear formulae:

$$\widehat{\mathit{X}}=\mathit{f}\left({\mathit{X}}^{\prime },\widehat{\mathit{\theta }}\right)$$E7

$$\widehat{\mathit{Z}}=\mathit{f}\left({\mathit{Z}}^{\prime },\widehat{\mathit{\theta }}\right).$$E8

Note that $\widehat{\mathit{\theta }}$ corresponds to the estimated parameters from the first step, whereas ${\mathit{X}}^{\prime }$ and ${\mathit{Z}}^{\prime }$ refer to the observed coordinates in the initial frame. The following Jacobian matrices are also defined here (to be used later on):

$${\mathit{J}}_{{\mathit{Z}}^{\prime }}=\frac{\partial \mathit{f}\left({\mathit{Z}}^{\prime },\mathit{\theta }\right)}{\partial {\mathit{Z}}^{\prime }},{\tilde{\mathit{J}}}_{\mathit{\theta }}=\frac{\partial \mathit{f}\left({\mathit{Z}}^{\prime },\mathit{\theta }\right)}{\partial \mathit{\theta }}$$E9

which differ from their previous counterparts in Eq. (6) as they refer to a separate group of points (NPs).

3.4 Deficiency of the stepwise approach

A number of drawbacks exist in the stepwise approach for SCT problems with noisy data. More specifically, (a) the noise of the original coordinates remains unfiltered during their transformation to the target frame, (b) the correlated errors in the original coordinates between CPs and NPs are not taken into account, and (c) the accuracy of the transformed coordinates is not optimized under any statistical principle. All these drawbacks relate to the same modeling deficiency that is summarized as follows: the observed coordinates in the initial frame are contaminated by random errors which remain uncontrolled during the second step of the transformation process, and they are fully absorbed by the transformed coordinates of CPs and NPs.

The aforesaid deficiency is irrelevant for practical applications only in two cases:

• if the sole objective is to determine a set of transformation parameters between different frames, without the need to perform any coordinate transformation at specific points; or

• if spatial objects (e.g., point cloud, network, digitized map) need to be transferred from an initial frame to another frame, without any “quality improvement” of their transformed coordinates.

However, if the user’s goal is the optimal referencing of spatial objects with respect to a target frame, then the unfiltered data noise becomes a critical error source for SCT problems. This does not mean that the stepwise approach leads to wrong results, but it signifies that the composite estimators in Eqs. (7) and (8) do not provide an optimal solution of maximum accuracy for the transformed coordinates.

It is worth noting that the stepwise approach is not compelled to reproduce the prior reference coordinates of CPs in the target frame, that is, $\widehat{\mathit{X}}\ne \mathit{X}$, even if these coordinates are perfectly known without any errors!

3.5 Best-fitting transformation solutions

In some cases, the estimation of transformation parameters is performed via the alternative nonlinear least squares principle [36, 37]:

$$\underset{\mathit{\theta }}{min}{‖\mathit{X}-\mathit{f}\left({\mathit{X}}^{\prime },\mathit{\theta }\right)‖}^2$$E10

where $‖\cdot ‖$ denotes the standard form of the Euclidean vector norm. The rationale of the above principle is to bring in the best alignment two different coordinate sets over a group of CPs, and it does not lead to the same parameter estimates as the statistical least squares formulation of Section 3.2. Their formal equivalency occurs if the known coordinates in the initial frame are treated as noiseless quantities and the respective coordinates in the target frame are affected by uncorrelated random errors of equal variance. Nevertheless, Eq. (10) has a strong geometrical significance, and it is often used in practice regardless of the noise characteristics of the available data.

If the transformation parameters are obtained by the alternative principle of Eq. (10), then it obviously holds that

$${‖\mathit{X}-\widehat{\mathit{X}}‖}^2\to min$$E11

which implies that the transformed coordinates of CPs will be optimally fitted, in a least squares sense, to their prior known values in the target frame. This best-fitting property does not enforce statistical optimality to the accuracy of the transformed coordinates—the latter will still absorb the entire observation noise according to Eqs. (7) and (8). Therefore, the point to be stressed here is that a high-quality transformation solution should not just rely on the fitting performance at CPs, but it has to exploit in an optimal sense the stochastic error model of the observed coordinates over all points of interest.

4. The stacking approach in spatial coordinate transformations

4.1 Theoretical aspects

A unified optimal solution for SCT problems can be obtained in a single stage through the rigorous combination of all available data. This requires the joint LSA of the nonlinear transformation equations:

$$\mathit{X}=\mathit{f}\left({\mathit{X}}^{\prime },\mathit{\theta }\right)$$E12

$$\mathit{Z}=\mathit{f}\left({\mathit{Z}}^{\prime },\mathit{\theta }\right)$$E13

which should be performed in a linearized context via the Newton-Gauss iteration method [3, 24, 34]. The algebraic setup of this stacking adjustment and the basic properties of the resulting estimators for the transformed coordinates are presented in this section.

4.1.1 Linearization

At first, we need to approximate the nonlinear Eqs. (12) and (13) by the truncated multivariate Taylor’s series expansions:

$$\mathit{X}=\mathit{f}\left({\mathit{X}}_o^{\prime },{\mathit{\theta }}_o\right)+{\mathit{J}}_{{\mathit{X}}^{\prime }}\left({\mathit{X}}^{\prime }-{\mathit{X}}_o^{\prime }\right)+{\mathit{J}}_{\mathit{\theta }}\left(\mathit{\theta }-{\mathit{\theta }}_o\right)$$E14

$$\mathit{Z}=\mathit{f}\left({\mathit{Z}}_o^{\prime },{\mathit{\theta }}_o\right)+{\mathit{J}}_{{\mathit{Z}}^{\prime }}\left({\mathit{Z}}^{\prime }-{\mathit{Z}}_o^{\prime }\right)+{\tilde{\mathit{J}}}_{\mathit{\theta }}\left(\mathit{\theta }-{\mathit{\theta }}_o\right)$$E15

where ${\mathit{\theta }}_o$ is a vector of approximate values for the transformation parameters and ${\mathit{X}}_o^{\prime }$, ${\mathit{Z}}_o^{\prime }$ are vectors of approximate coordinates for the respective points in the initial frame. Taking into account that the observables correspond to the coordinate vectors $\mathit{X}$, ${\mathit{X}}^{\prime }$, and ${\mathit{Z}}^{\prime }$, the previous formulae should be further augmented as follows:

$$\mathit{X}+{\mathit{v}}_{\mathit{X}}=\mathit{f}\left({\mathit{X}}_o^{\prime },{\mathit{\theta }}_o\right)+{\mathit{J}}_{{\mathit{X}}^{\prime }}\left({\mathit{X}}^{\prime }+{\mathit{v}}_{{\mathit{X}}^{\prime }}-{\mathit{X}}_o^{\prime }\right)+{\mathit{J}}_{\mathit{\theta }}\left(\mathit{\theta }-{\mathit{\theta }}_o\right)$$E16

$$\mathit{Z}=\mathit{f}\left({\mathit{Z}}_o^{\prime },{\mathit{\theta }}_o\right)+{\mathit{J}}_{{\mathit{Z}}^{\prime }}\left({\mathit{Z}}^{\prime }+{\mathit{v}}_{{\mathit{Z}}^{\prime }}-{\mathit{Z}}_o^{\prime }\right)+{\tilde{\mathit{J}}}_{\mathit{\theta }}\left(\mathit{\theta }-{\mathit{\theta }}_o\right)$$E17

where the added vectors ${\mathit{v}}_{\mathit{X}}$, ${\mathit{v}}_{{\mathit{X}}^{\prime }}$, ${\mathit{v}}_{{\mathit{Z}}^{\prime }}$ denote the zero-mean random errors of the observed coordinates. The linearized expressions (16) and (17) can be equivalently written in the block-matrix form:

$$\left[\begin{array}{cc}-{\mathit{J}}_{\mathit{\theta }}& \mathbf{0}\\ -{{\displaystyle \tilde{\mathit{J}}}}_{\mathit{\theta }}& \mathit{I}\end{array}\right]\left[\begin{array}{c}\mathit{\theta }-{\mathit{\theta }}_o\\ \mathit{Z}\end{array}\right]+\left[\begin{array}{ccc}\mathit{I}& -{\mathit{J}}_{{\mathit{X}}^{\prime }}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& -{\mathit{J}}_{{\mathit{Z}}^{\prime }}\end{array}\right]\left[\begin{array}{c}{\mathit{v}}_{\mathit{X}}\\ {\mathit{v}}_{{\mathit{X}}^{\prime }}\\ {\mathit{v}}_{{\mathit{Z}}^{\prime }}\end{array}\right]+\left[\begin{array}{c}\mathit{X}-\mathit{f}\left({\mathit{X}}_o^{\prime },{\mathit{\theta }}_o\right)-{\mathit{J}}_{{\mathit{X}}^{\prime }}\left({\mathit{X}}^{\prime }-{\mathit{X}}_o^{\prime }\right)\\ -\mathit{f}\left({\mathit{Z}}_o^{\prime },{\mathit{\theta }}_o\right)-{\mathit{J}}_{{\mathit{Z}}^{\prime }}\left({\mathit{Z}}^{\prime }-{\mathit{Z}}_o^{\prime }\right)\end{array}\right]=\left[\begin{array}{c}\mathbf{0}\\ \mathbf{0}\end{array}\right]$$E18

which conforms to the usual structure of Gauss-Helmert linear models of statistical estimation theory [12, 17, 24]. Our objective here is to invert the above stacked system of the general form Ax + Bv + w = 0 using the general least squares principle ${\mathit{v}}^T\mathit{P}\mathit{v}=min$, in conjunction with the data weight matrix:

$$\mathit{P}={\mathit{\Sigma }}^{-1}={\left[\begin{array}{ccc}{\mathit{\Sigma }}_{\mathit{X}}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\mathit{\Sigma }}_{{\mathit{X}}^{\prime }}& {\mathit{\Sigma }}_{{\mathit{X}}^{\prime }{\mathit{Z}}^{\prime }}\\ \mathbf{0}& {\mathit{\Sigma }}_{{\mathit{Z}}^{\prime }{\mathit{X}}^{\prime }}& {\mathit{\Sigma }}_{{\mathit{Z}}^{\prime }}\end{array}\right]}^{-1}$$E19

which reflects the total statistical accuracy of the observables. Note that the inter-frame correlations of the observed coordinates are assumed to be zero, whereas the intra-frame correlations between CPs and NPs are taken into account by the cross-CV matrix ${\mathit{\Sigma }}_{{\mathit{Z}}^{\prime }{\mathit{X}}^{\prime }}\left(={\mathit{\Sigma }}_{{\mathit{X}}^{\prime }{\mathit{Z}}^{\prime }}^T\right)$.

If applied under a proper iterative setting, the LSA of Eq. (18) leads to the sought optimal solution of the problem at hand. Specifically, the transformation parameters and the coordinates of NPs in the target frame are both contained into the “parameter vector” of the stacked GH-type model, and they can be directly obtained via the respective least squares estimator (see next section). On the other hand, the estimated coordinates of CPs in the target frame shall be deduced in an implicit way by correcting the observed values $\mathit{X}$ for the effect of their random errors (${\mathit{v}}_{\mathit{X}}$) which are also estimable from the iterative least squares inversion of Eq. (18).

4.1.2 Optimal least squares estimators

By applying the general LSA solution of linear GH models (see [12, 17]) to the stacked system of Eq. (18) and after some extra lengthy derivations using analytic inversions of 2 × 2 block matrices, we obtain the explicit estimators for the transformation parameters:

$$\widehat{\mathit{\theta }}={\mathit{\theta }}_o+{\left({\mathit{J}}_{\mathit{\theta }}^T{\mathit{WJ}}_{\mathit{\theta }}\right)}^{-1}{\mathit{J}}_{\mathit{\theta }}^T\mathit{W}\left(\mathit{X}-\mathit{f}\left({\mathit{X}}_o^{\prime },{\mathit{\theta }}_o\right)-{\mathit{J}}_{{\mathit{X}}^{\prime }}\left({\mathit{X}}^{\prime }-{\mathit{X}}_o^{\prime }\right)\right)$$E20

and for the coordinates of NPs in the target frame

$$\begin{array}{c}\widehat{\mathit{Z}}=\mathit{f}\left({\mathit{Z}}_o^{\prime },{\mathit{\theta }}_o\right)+{\mathit{J}}_{{\mathit{Z}}^{\prime }}\left({\mathit{Z}}^{\prime }-{\mathit{Z}}_o^{\prime }\right)+{\tilde{\mathit{J}}}_{\mathit{\theta }}\left(\widehat{\mathit{\theta }}-{\mathit{\theta }}_o\right)\\ +{\mathit{J}}_{{\mathit{Z}}^{\prime }}{\mathit{\Sigma }}_{{\mathit{Z}}^{\prime }{\mathit{X}}^{\prime }}{\mathit{J}}_{{\mathit{X}}^{\prime }}^T\mathit{W}\left(\mathit{X}-\mathit{f}\left({\mathit{X}}_o^{\prime },{\mathit{\theta }}_o\right)-{\mathit{J}}_{{\mathit{X}}^{\prime }}\left({\mathit{X}}^{\prime }-{\mathit{X}}_o^{\prime }\right)-{\mathit{J}}_{\mathit{\theta }}\left(\widehat{\mathit{\theta }}-{\mathit{\theta }}_o\right)\right)\end{array}$$E21

whereas the estimated errors for each subset of observed coordinates are given by the equation

$$\left[\begin{array}{c}{\widehat{\mathit{v}}}_{\mathit{X}}\\ {\widehat{\mathit{v}}}_{{\mathit{X}}^{\prime }}\\ {\widehat{\mathit{v}}}_{{\mathit{Z}}^{\prime }}\end{array}\right]=\left[\begin{array}{c}-{\mathit{\Sigma }}_{\mathit{X}}\\ {\mathit{\Sigma }}_{{\mathit{X}}^{\prime }}{\mathit{J}}_{{\mathit{X}}^{\prime }}^T\\ {\mathit{\Sigma }}_{{\mathit{Z}}^{\prime }{\mathit{X}}^{\prime }}{\mathit{J}}_{{\mathit{X}}^{\prime }}^T\end{array}\right]\mathit{W}\left(\mathit{X}-\mathit{f}\left({\mathit{X}}_o^{\prime },{\mathit{\theta }}_o\right)-{\mathit{J}}_{{\mathit{X}}^{\prime }}\left({\mathit{X}}^{\prime }-{\mathit{X}}_o^{\prime }\right)-{\mathit{J}}_{\mathit{\theta }}\left(\widehat{\mathit{\theta }}-{\mathit{\theta }}_o\right)\right)$$E22

The auxiliary matrix $\mathit{W}$ that appears in the previous equations was defined earlier in Section 3.2. Finally, if we combine the first error component from Eq. (22) with the basic formula $\widehat{\mathit{X}}=\mathit{X}+{\widehat{\mathit{v}}}_{\mathit{X}}$, we get the estimated coordinates of CPs in the target frame:

$$\begin{array}{c}\widehat{\mathit{X}}=\mathit{f}\left({\mathit{X}}_o^{\prime },{\mathit{\theta }}_o\right)+{\mathit{J}}_{{\mathit{X}}^{\prime }}\left({\mathit{X}}^{\prime }-{\mathit{X}}_o^{\prime }\right)+{\mathit{J}}_{\mathit{\theta }}\left(\widehat{\mathit{\theta }}-{\mathit{\theta }}_o\right)\\ +{\mathit{J}}_{{\mathit{X}}^{\prime }}{\mathit{\Sigma }}_{{\mathit{X}}^{\prime }}{\mathit{J}}_{{\mathit{X}}^{\prime }}^T\mathit{W}\left(\mathit{X}-\mathit{f}\left({\mathit{X}}_o^{\prime },{\mathit{\theta }}_o\right)-{\mathit{J}}_{{\mathit{X}}^{\prime }}\left({\mathit{X}}^{\prime }-{\mathit{X}}_o^{\prime }\right)-{\mathit{J}}_{\mathit{\theta }}\left(\widehat{\mathit{\theta }}-{\mathit{\theta }}_o\right)\right)\end{array}$$E23

To facilitate a comprehensive analysis of the stacking approach, it is useful to rewrite Eqs. (21) and (23) in the combined Kalman-like form:

$$\left[\begin{array}{c}\widehat{\mathit{X}}\\ \widehat{\mathit{Z}}\end{array}\right]=\left[\begin{array}{c}\stackrel{⌢}{\mathit{X}}\\ \stackrel{⌢}{\mathit{Z}}\end{array}\right]+\left[\begin{array}{c}{\mathit{J}}_{{\mathit{X}}^{\prime }}{\mathit{\Sigma }}_{{\mathit{X}}^{\prime }}{\mathit{J}}_{{\mathit{X}}^{\prime }}^T\\ {\mathit{J}}_{{\mathit{Z}}^{\prime }}{\mathit{\Sigma }}_{{\mathit{Z}}^{\prime }{\mathit{X}}^{\prime }}{\mathit{J}}_{{\mathit{X}}^{\prime }}^T\end{array}\right]\underset{\mathit{W}}{\underbrace{{\left({\mathit{\Sigma }}_{\mathit{X}}+{\mathit{J}}_{{\mathit{X}}^{\prime }}{\mathit{\Sigma }}_{{\mathit{X}}^{\prime }}{\mathit{J}}_{{\mathit{X}}^{\prime }}^T\right)}^{-1}}}\left(\mathit{X}-\stackrel{⌢}{\mathit{X}}\right)$$E24

where the auxiliary terms $\stackrel{⌢}{\mathit{X}}$ and $\stackrel{⌢}{\mathit{Z}}$ are strictly given by the expressions

$$\left[\begin{array}{c}\stackrel{⌢}{\mathit{X}}\\ \stackrel{⌢}{\mathit{Z}}\end{array}\right]=\left[\begin{array}{c}\mathit{f}\left({\mathit{X}}_o^{\prime },{\mathit{\theta }}_o\right)+{\mathit{J}}_{{\mathit{X}}^{\prime }}\left({\mathit{X}}^{\prime }-{\mathit{X}}_o^{\prime }\right)+{\mathit{J}}_{\mathit{\theta }}\left(\widehat{\mathit{\theta }}-{\mathit{\theta }}_o\right)\\ \mathit{f}\left({\mathit{Z}}_o^{\prime },{\mathit{\theta }}_o\right)+{\mathit{J}}_{{\mathit{Z}}^{\prime }}\left({\mathit{Z}}^{\prime }-{\mathit{Z}}_o^{\prime }\right)+{\tilde{\mathit{J}}}_{\mathit{\theta }}\left(\widehat{\mathit{\theta }}-{\mathit{\theta }}_o\right)\end{array}\right]$$E25

which, to a first-order approximation, mimic the result of the traditional stepwise approach, that is,

$$\left[\begin{array}{c}\stackrel{⌢}{\mathit{X}}\\ \stackrel{⌢}{\mathit{Z}}\end{array}\right]\simeq \left[\begin{array}{c}\mathit{f}\left({\mathit{X}}^{\prime },\widehat{\mathit{\theta }}\right)\\ \mathit{f}\left({\mathit{Z}}^{\prime },\widehat{\mathit{\theta }}\right)\end{array}\right]$$E26

All previous estimators refer to a single execution of the weighted LSA in the linearized system of Eq. (18). Their use in practical applications with nonlinear transformation models requires a recursive algorithm, as explained in more detail in Section 4.2.

4.1.3 Basic features of the stacking approach

Compared to the traditional stepwise methodology, the stacking approach leads to the same least squares estimate for the transformation parameters but to different values for the estimated coordinates in the target frame. This partial equivalency is expected since the inclusion of NPs into the adjustment procedure does not contribute additional information for the transformation parameters. On the other hand, the estimated coordinates contain extra corrections which are derived from stochastic filtering of the coordinate residuals $\mathit{X}-\stackrel{⌢}{\mathit{X}}$ and kriging-like prediction over all points of interest [see Eq. (24)]. Loosely speaking, the effect of those corrections resembles a rubber-sheeting process in the sense of “stretching” the classic stepwise solution to counteract the propagated data noise in the entire set of transformed coordinates.

The stacking approach permits also the exact fit over all CPs regardless of the noise level in the initial frame. This essential property is easily verified by Eq. (24) which implies that

$${\mathit{\Sigma }}_{\mathit{X}}=\mathbf{0}\to \widehat{\mathit{X}}=\mathit{X}$$E27

or in a loosened version

$${\mathit{\Sigma }}_{\mathit{X}}\ll {\mathit{J}}_{{\mathit{X}}^{\prime }}{\mathit{\Sigma }}_{{\mathit{X}}^{\prime }}{\mathit{J}}_{{\mathit{X}}^{\prime }}^T\to \widehat{\mathit{X}}\simeq \mathit{X}$$E28

The first condition dictates that the transformed coordinates of CPs will match their prior values, if the latter are assumed to be of perfect quality. The second condition is also useful for practical applications, as it allows the users to improve the fitting performance of the transformation results via a simple tuning of the CV matrix ${\mathit{\Sigma }}_{\mathit{X}}$. This last option is essentially equivalent to stochastic constraining of the prior coordinates of CPs in the target frame.

As a final note, let us point out that both approaches give similar results in the presence of noiseless data in the initial frame. In such case the least squares estimators of the previous section admit the conditional behavior:

$${\mathit{\Sigma }}_{{\mathit{X}}^{\prime }}=\mathbf{0}\to {\widehat{\mathit{v}}}_{{\mathit{X}}^{\prime }}=\mathbf{0}\to \widehat{\mathit{X}}=\stackrel{⌢}{\mathit{X}}\simeq \mathit{f}\left({\mathit{X}}^{\prime },\widehat{\mathit{\theta }}\right)$$E29

$${\mathit{\Sigma }}_{{\mathit{Z}}^{\prime }{\mathit{X}}^{\prime }}=\mathbf{0}\to {\widehat{\mathit{v}}}_{{\mathit{Z}}^{\prime }}=\mathbf{0}\to \widehat{\mathit{Z}}=\stackrel{⌢}{\mathit{Z}}\simeq \mathit{f}\left({\mathit{Z}}^{\prime },\widehat{\mathit{\theta }}\right).$$E30

Interestingly, the CV matrix ${\mathit{\Sigma }}_{{\mathit{Z}}^{\prime }}$ does not play an active role within the stacking approach, in contrast to the cross-CV matrix ${\mathit{\Sigma }}_{{\mathit{Z}}^{\prime }{\mathit{X}}^{\prime }}$ which is of crucial importance for the optimal transformation at the NPs [see Eq. (22)]. In Table 1 all relevant cases that can appear in SCT problems are classified with regard to the stochastic model of the observed coordinates in the respective frames.

	CV matrices of observed coordinates			Does data noise filtering occur in the transformation process ?
	Initial frame		Target frame
Control points	${\mathit{\Sigma }}_{{\mathit{X}}^{\prime }}\ne \mathbf{0}$		${\mathit{\Sigma }}_{\mathit{X}}\ne \mathbf{0}$	Yes
	${\mathit{\Sigma }}_{{\mathit{X}}^{\prime }}\ne \mathbf{0}$		${\mathit{\Sigma }}_{\mathit{X}}=\mathbf{0}$	Yes—perfect fit to prior values
	${\mathit{\Sigma }}_{{\mathit{X}}^{\prime }}=\mathbf{0}$		${\mathit{\Sigma }}_{\mathit{X}}\ne \mathbf{0}$	No
Non-control points	${\mathit{\Sigma }}_{{\mathit{Z}}^{\prime }}\ne \mathbf{0}$	${\mathit{\Sigma }}_{{\mathit{Z}}^{\prime }{\mathit{X}}^{\prime }}\ne \mathbf{0}$		Yes
	${\mathit{\Sigma }}_{{\mathit{Z}}^{\prime }}\ne \mathbf{0}$	${\mathit{\Sigma }}_{{\mathit{Z}}^{\prime }{\mathit{X}}^{\prime }}=\mathbf{0}$		No
	${\mathit{\Sigma }}_{{\mathit{Z}}^{\prime }}=\mathbf{0}$	${\mathit{\Sigma }}_{{\mathit{Z}}^{\prime }{\mathit{X}}^{\prime }}=\mathbf{0}$		No

Table 1.

Different cases in the stacking approach with regard to the stochastic model of the observed coordinates.

4.2 Computational aspects

The numerical computation of the stacking solution in nonlinear SCT problems requires a recursive implementation of the least squares estimators given in Section 4.1.2. The Newton-Gauss iteration method is suitable for this purpose and entails the updating of the approximate vectors ${\mathit{X}}_o^{\prime }$, ${\mathit{Z}}_o^{\prime }$, ${\mathit{\theta }}_o$ at each step by their adjusted values from the previous step until sufficient convergence is achieved in all estimated quantities of interest [3, 24, 34].

The aforesaid procedure should be applied for computing both the transformation parameters and the coordinates of CPs/NPs in the target frame, based on the following algorithm:

$${\widehat{\mathit{\theta }}}^{\left(k\right)}={\mathit{\theta }}_o^{\left(k\right)}+{\left({\mathit{J}}_{\mathit{\theta }}^T{\mathit{WJ}}_{\mathit{\theta }}\right)}^{-1}{\mathit{J}}_{\mathit{\theta }}^T\mathit{W}\left(\mathit{X}-\mathit{f}\left({{\mathit{X}}^{\prime }}_o^{\left(k\right)},{\mathit{\theta }}_o^{\left(k\right)}\right)-{\mathit{J}}_{{\mathit{X}}^{\prime }}\left({\mathit{X}}^{\prime }-{{\mathit{X}}^{\prime }}_o^{\left(k\right)}\right)\right)$$E31

$$\left[\begin{array}{c}{\stackrel{⌢}{\mathit{X}}}^{\left(k\right)}\\ {\stackrel{⌢}{\mathit{Z}}}^{\left(k\right)}\end{array}\right]=\left[\begin{array}{c}\mathit{f}\left({{\mathit{X}}^{\prime }}_o^{\left(k\right)},{\mathit{\theta }}_o^{\left(k\right)}\right)+{\mathit{J}}_{{\mathit{X}}^{\prime }}\left({\mathit{X}}^{\prime }-{{\mathit{X}}^{\prime }}_o^{\left(k\right)}\right)+{\mathit{J}}_{\mathit{\theta }}\left({\widehat{\mathit{\theta }}}^{\left(k\right)}-{\mathit{\theta }}_o^{\left(k\right)}\right)\\ \mathit{f}\left({{\mathit{Z}}^{\prime }}_o^{\left(k\right)},{\mathit{\theta }}_o^{\left(k\right)}\right)+{\mathit{J}}_{{\mathit{Z}}^{\prime }}\left({\mathit{Z}}^{\prime }-{{\mathit{Z}}^{\prime }}_o^{\left(k\right)}\right)+{\tilde{\mathit{J}}}_{\mathit{\theta }}\left({\widehat{\mathit{\theta }}}^{\left(k\right)}-{\mathit{\theta }}_o^{\left(k\right)}\right)\end{array}\right]$$E32

$$\left[\begin{array}{c}{\widehat{\mathit{v}}}_{{\mathit{X}}^{\prime }}^{\left(k\right)}\\ {\widehat{\mathit{v}}}_{{\mathit{Z}}^{\prime }}^{\left(k\right)}\end{array}\right]=\left[\begin{array}{c}{\mathit{\Sigma }}_{{\mathit{X}}^{\prime }}{\mathit{J}}_{{\mathit{X}}^{\prime }}^T\\ {\mathit{\Sigma }}_{{\mathit{Z}}^{\prime }{\mathit{X}}^{\prime }}{\mathit{J}}_{{\mathit{X}}^{\prime }}^T\end{array}\right]{\left({\mathit{\Sigma }}_{\mathit{X}}+{\mathit{J}}_{{\mathit{X}}^{\prime }}{\mathit{\Sigma }}_{{\mathit{X}}^{\prime }}{\mathit{J}}_{{\mathit{X}}^{\prime }}^T\right)}^{-1}\left(\mathit{X}-{\stackrel{⌢}{\mathit{X}}}^{\left(k\right)}\right)$$E33

$$\left[\begin{array}{c}{\widehat{\mathit{X}}}^{\left(k\right)}\\ {\widehat{\mathit{Z}}}^{\left(k\right)}\end{array}\right]=\left[\begin{array}{c}{\stackrel{⌢}{\mathit{X}}}^{\left(k\right)}\\ {\stackrel{⌢}{\mathit{Z}}}^{\left(k\right)}\end{array}\right]+\left[\begin{array}{c}{\mathit{J}}_{{\mathit{X}}^{\prime }}{\widehat{\mathit{v}}}_{{\mathit{X}}^{\prime }}^{\left(k\right)}\\ {\mathit{J}}_{{\mathit{Z}}^{\prime }}{\widehat{\mathit{v}}}_{{\mathit{Z}}^{\prime }}^{\left(k\right)}\end{array}\right]$$E34

where the index $k=1,2,\dots$ denotes the LSA iteration step. All Jacobian matrices shown in these equations should be re-evaluated at each step as follows:

$${\mathit{J}}_{{\mathit{X}}^{\prime }}={\left.\frac{\partial \mathit{f}\left({\mathit{X}}^{\prime },\mathit{\theta }\right)}{\partial {\mathit{X}}^{\prime }}\right|}_{{\mathit{\theta }}_o^{\left(k\right)},{{\mathit{X}}^{\prime }}_o^{\left(k\right)}},{\mathit{J}}_{{\mathit{Z}}^{\prime }}={\left.\frac{\partial \mathit{f}\left({\mathit{Z}}^{\prime },\mathit{\theta }\right)}{\partial {\mathit{Z}}^{\prime }}\right|}_{{\mathit{\theta }}_o^{\left(k\right)},{{\mathit{Z}}^{\prime }}_o^{\left(k\right)}}$$E35

$${\mathit{J}}_{\mathit{\theta }}={\left.\frac{\partial \mathit{f}\left({\mathit{X}}^{\prime },\mathit{\theta }\right)}{\partial \mathit{\theta }}\right|}_{{\mathit{\theta }}_o^{\left(k\right)},{{\mathit{X}}^{\prime }}_o^{\left(k\right)}},{\tilde{\mathit{J}}}_{\mathit{\theta }}={\left.\frac{\partial \mathit{f}\left({\mathit{Z}}^{\prime },\mathit{\theta }\right)}{\partial \mathit{\theta }}\right|}_{{\mathit{\theta }}_o^{\left(k\right)},{{\mathit{Z}}^{\prime }}_o^{\left(k\right)}}.$$E36

Note that the auxiliary weight matrix $\mathit{W}$ that appears in Eq. (31) depends on ${\mathit{J}}_{{\mathit{X}}^{\prime }}$ [see Eq. (5)] and it is also required to be updated at each step.

To initialize the Newton-Gauss iteration process, a simple choice is to set the approximate coordinates equal to the observed values (${\mathit{X}}_o^{\prime }={\mathit{X}}^{\prime }$, ${\mathit{Z}}_o^{\prime }={\mathit{Z}}^{\prime }$), while the approximate transformation parameters are typically obtained via empirical procedures. The initial computation of $\widehat{\mathit{\theta }}$ is thus reduced to the simpler form given already in Section 3.2, whereas for subsequent iterations the rigorous expression of Eq. (31) should be used. The updating of all approximate vectors at each step should be performed by the following equations:

$${\mathit{\theta }}_o^{\left(k\right)}={\widehat{\mathit{\theta }}}^{\left(k-1\right)}$$E37

$${{\mathit{X}}^{\prime }}_o^{\left(k\right)}={\mathit{X}}^{\prime }+{\widehat{\mathit{v}}}_{{\mathit{X}}^{\prime }}^{\left(k-1\right)}$$E38

$${{\mathit{Z}}^{\prime }}_o^{\left(k\right)}={\mathit{Z}}^{\prime }+{\widehat{\mathit{v}}}_{{\mathit{Z}}^{\prime }}^{\left(k-1\right)}$$E39

Special cases with noise-free coordinates in the initial frame (${\mathit{\Sigma }}_{{\mathit{X}}^{\prime }}=\mathbf{0}$) and/or uncorrelated coordinates between CPs and NPs (${\mathit{\Sigma }}_{{\mathit{Z}}^{\prime }{\mathit{X}}^{\prime }}=\mathbf{0}$) can be easily treated under the previous framework, and they lead to identical results as the traditional stepwise approach.

4.3 Statistical accuracy assessment in SCT solutions

The error CV matrices of $\widehat{\mathit{\theta }}$, $\widehat{\mathit{X}}$, and $\widehat{\mathit{Z}}$ are the fundamental elements for the formal quality assessment in SCT solutions. Their rigorous expressions are obtained by covariance propagation to the respective estimators given in previous sections, and they are presented here without their full mathematical proofs.

Both the stepwise and the stacking approach lead to the same optimal estimate for the transformation parameters, whose error CV matrix is given by the formula:

$${\mathit{\Sigma }}_{\widehat{\mathit{\theta }}}={\left({\mathit{J}}_{\mathit{\theta }}^T{\left({\mathit{\Sigma }}_{\mathit{X}}+{\mathit{J}}_{{\mathit{X}}^{\prime }}{\mathit{\Sigma }}_{{\mathit{X}}^{\prime }}{\mathit{J}}_{{\mathit{X}}^{\prime }}^T\right)}^{-1}{\mathit{J}}_{\mathit{\theta }}\right)}^{-1}.$$E40

Regarding the accuracy assessment of the transformed coordinates by the stepwise approach, the following expressions should be used:

$$\begin{array}{c}{\overline{\mathit{\Sigma }}}_{\widehat{\mathit{X}}}={\mathit{J}}_{\mathit{\theta }}{\mathit{\Sigma }}_{\widehat{\mathit{\theta }}}{\mathit{J}}_{\mathit{\theta }}^T+{\mathit{J}}_{{\mathit{X}}^{\prime }}{\mathit{\Sigma }}_{{\mathit{X}}^{\prime }}{\mathit{J}}_{{\mathit{X}}^{\prime }}^T-\left({\mathit{J}}_{\mathit{\theta }}{\mathit{\Sigma }}_{\widehat{\mathit{\theta }}}{\mathit{J}}_{\mathit{\theta }}^T\right){\left({\mathit{\Sigma }}_{\mathit{X}}+{\mathit{J}}_{{\mathit{X}}^{\prime }}{\mathit{\Sigma }}_{{\mathit{X}}^{\prime }}{\mathit{J}}_{{\mathit{X}}^{\prime }}^T\right)}^{-1}\left({\mathit{J}}_{{\mathit{X}}^{\prime }}{\mathit{\Sigma }}_{{\mathit{X}}^{\prime }}{\mathit{J}}_{{\mathit{X}}^{\prime }}^T\right)\\ -\left({\mathit{J}}_{{\mathit{X}}^{\prime }}{\mathit{\Sigma }}_{{\mathit{X}}^{\prime }}{\mathit{J}}_{{\mathit{X}}^{\prime }}^T\right){\left({\mathit{\Sigma }}_{\mathit{X}}+{\mathit{J}}_{{\mathit{X}}^{\prime }}{\mathit{\Sigma }}_{{\mathit{X}}^{\prime }}{\mathit{J}}_{{\mathit{X}}^{\prime }}^T\right)}^{-1}\left({\mathit{J}}_{\mathit{\theta }}{\mathit{\Sigma }}_{\widehat{\mathit{\theta }}}{\mathit{J}}_{\mathit{\theta }}^T\right)\end{array}$$E41

$$\begin{array}{c}{\overline{\mathit{\Sigma }}}_{\widehat{\mathit{Z}}}={\tilde{\mathit{J}}}_{\mathit{\theta }}{\mathit{\Sigma }}_{\widehat{\mathit{\theta }}}{\tilde{\mathit{J}}}_{\mathit{\theta }}^T+{\mathit{J}}_{{\mathit{Z}}^{\prime }}{\mathit{\Sigma }}_{{\mathit{Z}}^{\prime }}{\mathit{J}}_{{\mathit{Z}}^{\prime }}^T-\left({\tilde{\mathit{J}}}_{\mathit{\theta }}{\mathit{\Sigma }}_{\widehat{\mathit{\theta }}}{\mathit{J}}_{\mathit{\theta }}^T\right){\left({\mathit{\Sigma }}_{\mathit{X}}+{\mathit{J}}_{{\mathit{X}}^{\prime }}{\mathit{\Sigma }}_{{\mathit{X}}^{\prime }}{\mathit{J}}_{{\mathit{X}}^{\prime }}^T\right)}^{-1}\left({\mathit{J}}_{{\mathit{X}}^{\prime }}{\mathit{\Sigma }}_{{\mathit{X}}^{\prime }{\mathit{Z}}^{\prime }}{\mathit{J}}_{{\mathit{Z}}^{\prime }}^T\right)\\ -\left({\mathit{J}}_{{\mathit{Z}}^{\prime }}{\mathit{\Sigma }}_{{\mathit{Z}}^{\prime }{\mathit{X}}^{\prime }}{\mathit{J}}_{{\mathit{X}}^{\prime }}^T\right){\left({\mathit{\Sigma }}_{\mathit{X}}+{\mathit{J}}_{{\mathit{X}}^{\prime }}{\mathit{\Sigma }}_{{\mathit{X}}^{\prime }}{\mathit{J}}_{{\mathit{X}}^{\prime }}^T\right)}^{-1}\left({\mathit{J}}_{\mathit{\theta }}{\mathit{\Sigma }}_{\widehat{\mathit{\theta }}}{\tilde{\mathit{J}}}_{\mathit{\theta }}^T\right)\end{array}$$E42

which refer to the CPs and NPs, respectively. The overbar symbol is used to distinguish the above error CV matrices from the respective expressions that apply in the stacking approach. The latter are given by the general formulae:

$${\mathit{\Sigma }}_{\widehat{\mathit{X}}}={\overline{\mathit{\Sigma }}}_{\widehat{\mathit{X}}}-{\mathit{K\Sigma }}_{\mathit{e}}{\mathit{K}}^T$$E43

$${\mathit{\Sigma }}_{\widehat{\mathit{Z}}}={\overline{\mathit{\Sigma }}}_{\widehat{\mathit{Z}}}-{\mathit{Q\Sigma }}_{\mathit{e}}{\mathit{Q}}^T$$E44

where the auxiliary matrices $\mathit{K}$ and $\mathit{Q}$ are defined as

$$\mathit{K}=\left({\mathit{J}}_{{\mathit{X}}^{\prime }}{\mathit{\Sigma }}_{{\mathit{X}}^{\prime }}{\mathit{J}}_{{\mathit{X}}^{\prime }}^T\right){\left({\mathit{\Sigma }}_{\mathit{X}}+{\mathit{J}}_{{\mathit{X}}^{\prime }}{\mathit{\Sigma }}_{{\mathit{X}}^{\prime }}{\mathit{J}}_{{\mathit{X}}^{\prime }}^T\right)}^{-1}$$E45

$$\mathit{Q}=\left({\mathit{J}}_{{\mathit{Z}}^{\prime }}{\mathit{\Sigma }}_{{\mathit{Z}}^{\prime }{\mathit{X}}^{\prime }}{\mathit{J}}_{{\mathit{X}}^{\prime }}^T\right){\left({\mathit{\Sigma }}_{\mathit{X}}+{\mathit{J}}_{{\mathit{X}}^{\prime }}{\mathit{\Sigma }}_{{\mathit{X}}^{\prime }}{\mathit{J}}_{{\mathit{X}}^{\prime }}^T\right)}^{-1}$$E46

and ${\mathit{\Sigma }}_{\mathit{e}}$ is the CV matrix of the coordinate residuals $\mathit{X}-\mathit{f}\left({\mathit{X}}^{\prime },\widehat{\mathit{\theta }}\right)$, that is,

$${\mathit{\Sigma }}_{\mathit{e}}={\mathit{\Sigma }}_{\mathit{X}}+{\mathit{J}}_{{\mathit{X}}^{\prime }}{\mathit{\Sigma }}_{{\mathit{X}}^{\prime }}{\mathit{J}}_{{\mathit{X}}^{\prime }}^T-{\mathit{J}}_{\mathit{\theta }}{\mathit{\Sigma }}_{\widehat{\mathit{\theta }}}{\mathit{J}}_{\mathit{\theta }}^T.$$E47

Equations (43) and (44) reveal the expected improvement of the statistical accuracy in the SCT solution by the stacking approach. The diagonal elements (i.e., coordinate error variances) of ${\mathit{\Sigma }}_{\widehat{\mathit{X}}}$ and ${\mathit{\Sigma }}_{\widehat{\mathit{Z}}}$ are always smaller than the respective elements of ${\overline{\mathit{\Sigma }}}_{\widehat{\mathit{X}}}$ and ${\overline{\mathit{\Sigma }}}_{\widehat{\mathit{Z}}}$, a fact that is attributed to the noise filtering of the observed coordinates during the transformation process.

5. Numerical example

To demonstrate the potential of the stacking approach in practical transformation problems, a simple example is given here for a simulated 2D network with seven CPs and four NPs. The true coordinates of all network points are listed in Table 2, and they are related by a second-order polynomial transformation:

	Initial frame		Target frame
	x′ (m)	y′ (m)	x (m)	y (m)
CP1	100.000	250.000	146.000	287.000
CP2	200.000	423.205	210.768	467.597
CP3	286.602	373.205	239.979	435.802
CP4	157.735	150.000	181.119	177.203
CP5	125.000	200.000	159.250	231.438
CP6	225.000	250.000	222.875	294.188
CP7	250.000	400.000	226.250	452.750
NP1	200.000	300.000	209.250	342.500
NP2	159.000	230.000	184.574	264.665
NP3	220.000	340.000	217.850	386.660
NP4	170.000	270.000	192.750	308.030

Table 2.

True coordinates of CPs and NPs in the simulated test network with respect to the initial and target frame.

$$x_i=a_o+a_1{x}_i^{\prime }+a_2{y}_i^{\prime }+a_3{x}_i^{\prime }{y}_i^{\prime }+a_4{x^{\prime }}_i^2+a_5{y^{\prime }}_i^2$$E48

$$y_i=b_o+b_1{x}_i^{\prime }+b_2{y}_i^{\prime }+b_3{x}_i^{\prime }{y}_i^{\prime }+b_4{x^{\prime }}_i^2+b_5{y^{\prime }}_i^2$$E49

whose associated parameters are provided in Table 3.

ao	a1	a2	a3	a4	a5
10.25	1.20	0.20	−0.0013	−0.0008	0.0001
bo	b1	b2	b3	b4	b5
18.50	−0.25	1.20	−0.0002	0.0011	−0.0002

Table 3.

True parameter values of the second-order polynomial transformation model.

The observed coordinates for our experiments stem by adding simulated Gaussian noise to the true values of Table 2. The known coordinates of NPs in the target frame are not included in the observables, but they were used only for cross-validation of the transformation results. The generated random errors at the CPs in the target frame are uncorrelated with a common standard deviation of 0.1 cm for the x and y coordinates. On the other hand, the generated random errors at the CPs/NPs in the initial frame are spatially correlated in terms of the simplified Gaussian-type covariance model:

$${\sigma }_{x_i^{\prime }{x}_k^{\prime }}={\sigma }_{y_i^{\prime }{y}_k^{\prime }}={\sigma }^2{e}^{-A{\left(x_i^{\prime }-x_k^{\prime }\right)}^2-B{\left(y_i^{\prime }-y_k^{\prime }\right)}^2}$$E50

$${\sigma }_{x_i^{\prime }{y}_k^{\prime }}=\rho {\sigma }^2{e}^{-A{\left(x_i^{\prime }-x_k^{\prime }\right)}^2-B{\left(y_i^{\prime }-y_k^{\prime }\right)}^2}$$E51

where $\sigma$ is the common error standard deviation for the x′ and y′ coordinates (set equal to 5 cm) and $\rho$ is their error correlation coefficient at each point (set equal to −0.2). The values of the auxiliary parameters A and B were fixed to 6 × 10⁻⁷ and 7 × 10⁻⁶, respectively, which ensure the positive definiteness of the resulting CV matrix for the observed coordinates in the initial frame.

Using a Monte Carlo sampling scheme and a Cholesky-based algorithm for the stochastic simulation of correlated random vectors, a total of 1000 noisy ensembles were produced for the triplet of coordinate vectors $\mathit{X}$, ${\mathit{X}}^{\prime }$, and ${\mathit{Z}}^{\prime }$. These synthetic datasets were used with the stepwise and stacking approach to determine the transformed coordinates and their associated accuracy, over all points of the simulated network.

The differences between the true and the transformed coordinates in the target frame, as obtained by all data ensembles under each approach, are shown in Figures 1 and 2. The cloud plots in these figures refer only to a subset of the CPs/NPs, yet similar results are acquired at all other network points. It is clear that the stacking approach yields significantly better results than the traditional stepwise approach, and it effectively filters the existing noise of the initial coordinates. The accuracy improvement ranges from 88 to 92% at the CPs, while it is a bit lower (63–78%) at the NPs (see detailed results in Table 4).

Figure 1.

Differences between the true and the transformed coordinates at two CPs over 1000 Monte Carlo sampling experiments in the simulated network. The point error ellipses (99% confidence level) by each transformation approach are also shown in red color. The scaling of the horizontal axes is in meters.

Figure 2.

Differences between the true and the transformed coordinates at two NPs over 1000 Monte Carlo sampling experiments in the simulated network. The point error ellipses (99% confidence level) by each transformation approach are also shown in red color. The scaling of the horizontal axes is in meters.

	Stepwise approach			Stacking approach			Accuracy improvement (%)
	${\sigma }_{\widehat{x}}$	${\sigma }_{\widehat{y}}$	$\sqrt{{\sigma }_{\widehat{x}}^2+{\sigma }_{\widehat{y}}^2}$	${\sigma }_{\widehat{x}}$	${\sigma }_{\widehat{y}}$	$\sqrt{{\sigma }_{\widehat{x}}^2+{\sigma }_{\widehat{y}}^2}$
CP1	0.6	0.7	0.9	0.1	0.1	0.1	89
CP2	0.3	0.9	0.9	0.1	0.1	0.1	89
CP3	0.3	0.7	0.7	0.1	0.1	0.1	86
CP4	0.6	1.0	1.2	0.1	0.1	0.1	92
CP5	0.6	0.8	1.0	0.1	0.1	0.1	90
CP6	0.4	0.6	0.8	0.1	0.1	0.1	88
CP7	0.3	0.8	0.8	0.1	0.1	0.1	88
NP1	0.5	0.7	0.8	0.2	0.2	0.3	63
NP2	0.5	0.7	0.9	0.1	0.1	0.2	78
NP3	0.4	0.7	0.8	0.2	0.3	0.3	63
NP4	0.5	0.7	0.8	0.2	0.2	0.3	63

Table 4.

Standard deviations of the transformed coordinates in the simulated network by different approaches.

All values given in cm.

It should be emphasized that the stochastic model of the observed coordinates plays a key role in the performance of the stacking approach. This means that the results shown here may exhibit different behavior—displaying either insignificant or even more profound accuracy improvement for the transformed coordinates—for varied choices of the CV matrices ${\mathit{\Sigma }}_{\mathit{X}}$, ${\mathit{\Sigma }}_{{\mathit{X}}^{\prime }}$, ${\mathit{\Sigma }}_{{\mathit{Z}}^{\prime }}$, and ${\mathit{\Sigma }}_{{\mathit{X}}^{\prime }{\mathit{Z}}^{\prime }}$.