Estimating the 7-Parameter Datum Transformation #

The general approach for converting geodetic data from one geodetic coordinate system to another is a transformation involving three translations, three rotations and a scale change. If we denote geocentric coordinates of a point in the source coordinate system as $(u,v,w)$ and the geocentric coordinates in the destination coordinate system as $(x,y,z)$ the following equations describe the seven parameters transformation:

$$ \begin{bmatrix}x\cr y\cr z\end{bmatrix} = (1+k)R(r_x,r_y,r_z)\begin{bmatrix}u\cr v\cr w\end{bmatrix} + T $$

where $k$ is the scale change, $T$ is a translation vector:

$$ T = \begin{bmatrix}t_x\cr t_y\cr t_z\end{bmatrix} $$

and $R$ is a rotation matrix:

$$ R(r_x, r_y,r_z) = R_3(r_z) \cdot R_2(r_y) \cdot R_1(r_x) $$

In the above formula $R_1,R_2$ and $R_3$ are rotation matrixes around each axis of the coordinate system:

$$ \eqalign{ R_1(\phi) =& \begin{bmatrix}1 & 0 & 0 \cr 0 & \cos\phi & \sin\phi \cr 0 & -\sin\phi & \cos\phi\end{bmatrix} \cr R_2(\phi) =& \begin{bmatrix}\cos\phi & 0 & -\sin\phi \cr 0 & 1 & 0 \cr \sin\phi & 0 & \cos\phi\end{bmatrix} \cr R_3(\phi) =& \begin{bmatrix}\cos\phi & \sin\phi & 0 \cr -\sin\phi & \cos\phi & 0 \cr 0 & 0 & 1\end{bmatrix} } $$

For nearly-aligned coordinate systems the rotation matrix can be approximated by:

$$ R(r_x, r_y, r_z) = \begin{bmatrix}1 & r_z & -r_y \cr -r_z & 1 & r_x \cr r_y & -r_z & 1\end{bmatrix} $$

and the equation for a 7-parameter transformation model becomes:

$$ \begin{bmatrix}x \cr y \cr z\end{bmatrix} = (1+k)\begin{bmatrix}1 & r_z & -r_y \cr -r_z & 1 & r_x \cr r_y & -r_z & 1\end{bmatrix}\begin{bmatrix}u \cr v \cr w\end{bmatrix}+\begin{bmatrix}t_x \cr t_y \cr t_z\end{bmatrix} $$

neglecting second order terms in $k, r_x, r_y, r_z$ and their products this further simplifies to:

$$ \begin{bmatrix}x \cr y \cr z\end{bmatrix} = \begin{bmatrix}1+k & r_z & -r_y \cr -r_z & 1+k & r_x \cr r_y & -r_z & 1+k\end{bmatrix}\begin{bmatrix}u \cr v \cr w\end{bmatrix}+\begin{bmatrix}t_x \cr t_y \cr t_z\end{bmatrix} \tag{ 1 } $$

In order to determine the 7 parameters, one needs to know geocentric coordinates in both systems for at least 3 stations. A least-squares approximation can then be obtained by minimizing the following function:

$$ \eqalign { S(t_x, t_y, t_z, r_x, r_y, r_z, k) = \sum_i[&(k\,u_i + r_zv_i - r_yw_i + t_x + u_i -x_i)^2 + \cr &(-r_zu_i +k\,v_i+r_xw_i+ t_y + v_i - y_i)^2 + \cr &(r_yu_i-r_xv_i+k\;w_i+t_z+w_i-z_i)^2]} $$

which is equivalent to solving the system:

$$ \eqalign{ \frac{\partial}{\partial t_x}S(t_x, t_y, t_z, r_x, r_y, r_z, k) &= 0 \cr \frac{\partial}{\partial t_y}S(t_x, t_y, t_z, r_x, r_y, r_z, k) &= 0 \cr \frac{\partial}{\partial t_z}S(t_x, t_y, t_z, r_x, r_y, r_z, k) &= 0 \cr \frac{\partial}{\partial r_x}S(t_x, t_y, t_z, r_x, r_y, r_z, k) &= 0 \cr \frac{\partial}{\partial r_y}S(t_x, t_y, t_z, r_x, r_y, r_z, k) &= 0 \cr \frac{\partial}{\partial r_z}S(t_x, t_y, t_z, r_x, r_y, r_z, k) &= 0 \cr \frac{\partial}{\partial k}S(t_x, t_y, t_z, r_x, r_y, r_z, k) &= 0 \cr } $$

That gives:

$$ \eqalign{ k\sum_i u_i + r_z\sum_iv_i-r_y\sum_iw_i + n\,t_x &= \sum_i(x_i-u_i) \cr k\sum_i v_i + r_x\sum_iw_i-r_z\sum_iu_i + n\,t_y &= \sum_i(y_i-v_i) \cr k\sum_i w_i - r_x\sum_iv_i+r_y\sum_iu_i + n\,t_z &= \sum_i(y_i-v_i) \cr r_x\sum_i(w_i^2+v_i^2) - r_y\sum_i(u_iv_i) - r_z\sum_i(u_iw_i) + t_y\sum_iw_i - t_z\sum_iv_i &= \sum_i(v_iz_i - w_iy_i)\cr -r_x\sum_i(u_iv_i) + r_y\sum_i(u_i^2+w_i^2) - r_z\sum(w_iv_i)-t_x\sum_iw_i+t_z\sum_iu_i &= \sum_i(u_iz_i-w_ix_i)\cr -r_x\sum(u_iw_i) -r_y\sum_i(v_iw_i) +r_z\sum_i(v_i^2+u_i^2)+t_x\sum_iv_i-t_y\sum_iu_i &= \sum_i(v_ix_i-u_iy_i)\cr k\sum_i(u_i^2+v_i^2+w_i^2)+t_x\sum_iu_i+t_y\sum_iv_i+t_z\sum_iw_i &= \sum_i(u_ix_i+v_iy_i+w_iz_i-u_i^2-v_i^2-w_i^2) } $$

Rewriting the system as:

$$ A\begin{bmatrix}t_x \cr t_y \cr t_z \cr r_x \cr r_y \cr r_z \cr k\end{bmatrix} = B $$

with:

$$ A = \begin{bmatrix} n & 0 & 0 & 0 & -\sum{w_i} & \sum{v_i} & \sum{u_i} \cr 0 & n & 0 & \sum{w_i} & 0 & -\sum{u_i} & \sum{v_i} \cr 0 & 0 & n & -\sum{v_i} & \sum{u_i} & 0 & \sum{w_i} \cr 0 & \sum{w_i} & -\sum{v_i} & \sum{w_i^2+v_i^2} & -\sum{u_iv_i} & -\sum{u_iw_i} & 0 \cr -\sum{w_i} & 0 & \sum{u_i} & -\sum{u_iv_i} & \sum{u_i^2+w_i^2} & -\sum{w_iv_i} & 0 \cr \sum{v_i} & -\sum{u_i} & 0 & -\sum{u_iw_i} & -\sum{v_iw_i} & \sum{v_i^2+u_i^2} & 0 \cr \sum{u_i} & \sum{v_i} & \sum{w_i} & 0 & 0 & 0 & \sum{u_i^2+v_i^2+w_i^2} \end{bmatrix} $$

and

$$ B = \begin{bmatrix}\sum{(x_i-u_i)} \cr \sum{(y_i-v_i)} \cr \sum{(z_i-w_i)} \cr \sum{(w_iy_i-v_iz_i)} \cr \sum{(u_iz_i-w_ix_i)} \cr \sum{(v_ix_i - u_iy_i)} \cr \sum{(u_ix_i + v_iy_i + w_iz_i -u_i^2 -v_i^2-w_i^2)} \end{bmatrix} $$

The solution will then be:

$$ \begin{bmatrix}t_x \cr t_y \cr t_z \cr r_x \cr r_y \cr r_z \cr k\end{bmatrix} = A^{-1}\, B \tag{ 2 } $$

Numerical Example #

(Data provided by Bureau Pierre Martin SA, Thierrens, Switzerland)

Geocentric coordinates on WGS84 ellipsoid are:

$$ u = \pmatrix{ 4331297.24\cr 4273147.84\cr 4253563.45\cr 4377795.40\cr 4390113.24} v = \pmatrix{567555.67\cr 575368.33\cr 733522.39\cr 468008.67\cr 696884.2} w = \pmatrix{4633133.80\cr 4684903.72\cr 4681452.19\cr 4601077.35\cr 4561133.03} $$

Geocentric coordinates on Bessel ellipsoid for the same stations are:

$$ x = \pmatrix{4330623.01\cr 4272474.04\cr 4252889.03\cr 4377121.28\cr 4389437.87} y = \pmatrix{567540.82\cr 575352.96\cr 733505.05\cr 467994.72\cr 696869.17} z = \pmatrix{4632728.32\cr 4684498.14\cr 4681047.3\cr 4600671.53\cr 4560728.3} $$

Using the method described above (equation 2), the datum transformation parameters are:

$$ \pmatrix{t_x\cr t_y\cr t_z\cr\ r_x\cr r_y\cr r_z\cr k} =\pmatrix{-651.287\cr -14.197\cr -362.266\cr -2.905\cdot cc\cr -1.698\cdot cc\cr -3.611\cdot cc\cr -7.399\cdot 10^{-6}} $$

All rotation angles are expressed in centesimal seconds. Conversion factor to radians is:

$$ cc = \pi/(200\cdot10^4) $$

Using these parameters and equation (1), we calculate transformed coordinates:

$$ x' = \pmatrix{4330623.04\cr 4272474.16\cr 4252889.01\cr 4377121.33\cr 4389437.68} y' = \pmatrix{567540.69\cr 575352.73\cr 733505.52\cr 467994.84\cr 696868.93} z' = \pmatrix{4632728.29\cr 4684498.02\cr 4681047.29\cr 4600671.50\cr 4560728.49} $$

Residual error is:

$$ r = |x-x'|^2 + |y-y'|^2+|z-z'|^2 = 0.474 $$

The results compare well against official datum transformation parameters published by Federal Office of Topography (Aug 1990):

$$ \pmatrix{t_x\cr t_y\cr t_z\cr\ r_x\cr r_y\cr r_z\cr k} =\pmatrix{-660.077\cr -13.551\cr -369.34\cr -2.484\cdot cc\cr -1.783\cdot cc\cr -2.939\cdot cc\cr -5.66\cdot 10^{-6}} $$

Using official datum transformation parameters, the transformed coordinates and residual error are:

$$ x' = \pmatrix{4330623.00\cr 4272474.04\cr 4252889.02\cr 4377121.27\cr 4389437.87} y' = \pmatrix{567540.82\cr 575352.97\cr 733506.06\cr 467994.73\cr 696869.18} z' = \pmatrix {4632728.32\cr 4684498.14\cr 4681047.30\cr 4600671.53\cr 4560728.30} $$

$$ r = |x-x'|^2 + |y-y'|^2+|z-z'|^2 = 1.016 $$