空间坐标与投影系统系列(四)：七参数转换实例

2018-08-17

上一篇介绍了坐标转换。本篇文章我们将以实例来介绍七参数的转换过程。

七参数转换法

七参数法（包括布尔莎模型，一步法模型，海尔曼特等）是解决此问题的比较严密和通用的方法¹。其涉及到的七个参数为：X平移，Y平移，Z平移，X旋转，Y旋转，Z旋转，尺度变化m。此七个参数可以通过在需要转化的区域里选取3个以上的转换控制点对而获取。
如果区域范围不大，最远点间的距离不大于30Km，可以用三参数（莫洛登斯基模型），即X平移，Y平移，Z平移，而此时将X旋转，Y旋转，Z旋转，尺度变化m视为0。所以三参数只是七参数的一种特例。三参数只需通过1个控制点对就能获取。

七参数公式如下:

$$
\begin{bmatrix} {X}\\ {Y}\\ {Z} \end{bmatrix}=\begin{bmatrix} {\Delta X}\\ {\Delta Y}\\ {\Delta Z} \end{bmatrix}+(1+m)\begin{bmatrix} {1}&{\theta_z}&{-\theta_y} \\ {-\theta_z}&{1}&{\theta_x}\\ {\theta_y}&{-\theta_x}&{1} \end{bmatrix}\begin{bmatrix} {X_0}\\ {Y_0}\\ {Z_0} \end{bmatrix}
$$

式中, $\Delta X$, $\Delta Y$, $\Delta Z$为三个平移参数; $θ_x$, $θ_y$, $θ_z$为3个旋转参数，m为尺度参数。将七参数作为未知数，并且按参数线性化，则可以转换为：

$$
\begin{bmatrix}
{X - X(X_0)}\\ {Y - Y(Y_0)}\\ {Z - Z(Z_0)}
\end{bmatrix}

\begin{bmatrix}
{1}&{0}&{0}&{0}&{-(1+m)Z_0}&{(1+m)Y_0}&{X_0+\theta_zY_0-\theta_yZ_0} \\ {0}&{1}&{0}&{(1+m)Z_0}&{0}&{-(1+m)X_0}&{-\theta_zX_0 +Y_0+\theta_xZ_0} \\ {0}&{0}&{1}&{-(1+m)Y_0}&{(1+m)X_0}&{0}&{\theta_yX_0-\theta_xY_0+Z_0}
\end{bmatrix}
\begin{bmatrix} {\Delta X}\\ {\Delta Y}\\ {\Delta Z}\\ {\theta_x}\\ {\theta_y}\\ {\theta_z}\\ {m}\end{bmatrix}
$$

其中，X(X0),Y(Y0),Z(Z0)由七参数公式给出，利用最小二乘法法即可对该线性化的方程进行迭代求解。

七参数转换法实例

下面以python代码为例²，利用已有四个控制点坐标，根据最小二乘来求解七参数模型的参数，并利用一个检查点来验证其转换精度。

导入数据

使用4个同名点的原坐标(x, y, z)和目标坐标（x’, y’, z’)来求解七参数

import numpy as np

source_points = [
    (3381400.980, 395422.030, 32.956),
    (3381404.344, 395844.239, 32.207),
    (3382149.810, 396003.592, 33.290),
    (3382537.793, 395985.359, 33.412)
]
target_points = [
    (3380968.194, 539468.888, 13.875),
    (3380977.154, 539890.934, 13.179),
    (3381724.612, 540040.47,  14.273),
    (3381724.636, 540040.485, 14.282)
]

#归一化处理便于模型的快速迭代
ave_src = np.mean(np.array(source_points), axis=0)
ave_tar = np.mean(np.array(target_points), axis=0)
source_points -= ave_src
target_points -= ave_tar

正向计算

根据坐标转换公式来列立系数方程定义函数返回系数矩阵。使用4个同名点的坐标转换公式作为误差方程进行最小二乘平差，根据坐标转换公式来列立系数方程 $W x = b$ :

#  定义函数返回系数矩阵 B， l
#  定义函数： point2matrix，
#     通过给定的同名点坐标列立误差方程B系数阵的部分
#     x, y, z： 原坐标值
#    args: 七参数误差值[Delta_X, Delta_Y, Delta_Z, theta_x, theta_y, theta_z, m]
#   返回值： W系数阵的部分
       
def point2matrix(x, y, z, args):
    array = [
        [1, 0, 0,              0, -(1+args[6])*z,  (1+args[6])*y,          x + args[5]*y - args[4]*z],
        [0, 1, 0,  (1+args[6])*z,              0, -(1+args[6])*x, -args[5]*x +         y + args[3]*z],
        [0, 0, 1, -(1+args[6])*y,  (1+args[6])*x,              0,  args[4]*x - args[3]*y +         z]
    ]
    return np.array(array)

# 定义函数： points2W
#      通过同名点序列列立误差方程B系数阵的整体
#       x, y, z： 同名点序列
#       args: 七参数误差值[Delta_X, Delta_Y, Delta_Z, theta_x, theta_y, theta_z, m]
#       返回值： W系数阵

def points2W(points, args):   
    big_mat = None 
    for (x, y, z) in points:
        mat = point2matrix(x, y, z, args)
        if big_mat is None:
            big_mat = mat
        else:
            big_mat = np.concatenate((big_mat, mat))

    return big_mat

# 定义函数： points2b
#       通过同名点坐标转换关系列立误差方程B系数阵的整体
#       x, y, z： 同名点的原坐标和目标坐标对组成的序列
#       args: 七参数误差值[Delta_X, Delta_Y, Delta_Z, theta_x, theta_y, theta_z, m]
#       返回值： b系数阵

def points2b(source, target, args):   
    big_mat = [0] * len(source) * 3

    for i, ((x1, y1, z1), (x2, y2, z2)) in enumerate(zip(source, target)):
        (x0, y0, z0) = ordinationConvert(x1, y1, z1, args)
        big_mat[3*i + 0] = x2 - x0
        big_mat[3*i + 1] = y2 - y0
        big_mat[3*i + 2] = z2 - z0

    return np.array(big_mat).transpose()

def ordinationConvert(x1, y1, z1, args):
    x2 = args[0] + (1+args[6])*(         x1 + args[5]*y1 - args[4]*z1)
    y2 = args[1] + (1+args[6])*(-args[5]*x1 +         y1 + args[3]*z1)
    z2 = args[2] + (1+args[6])*( args[4]*x1 - args[3]*y1 +         z1)
    return (x2, y2, z2)

系数求解

组成法方程 (W’W) x = (W’b) 并利用最小二乘求解x

Args = np.array([0, 0, 0, 0, 0, 0, 0], dtype='float64')
parameters = np.array([1, 1, 1, 1, 1, 1, 1])

# 当七参数的误差值之和大于1e-10时，迭代运算得到更精确的结果
while np.fabs(np.array(parameters)).sum() > 1e-10 :
    W = points2W(source_points, Args)
    b = points2b(source_points, target_points, Args)
    qxx = np.linalg.inv(np.dot(W.transpose(), W))
    parameters = np.dot(np.dot(qxx, W.transpose()), b)
    Args += parameters
	
# 打印七参数
print("Args:")
print(np.round(Args, 3))

七参数输出结果为
$$
[-0. -0. -0. -0. 0. -0.072 -0.215]
$$

模型精度评定和结果验证

评定最小二乘模型的精度结果，并使用预留的已知坐标同名点来验证七参数模型的效果

# 检查点坐标
source_test_points = [
    (3381402.058, 395657.940, 32.728)
]

target_test_points = [
    (3380972.424, 539704.811, 13.665)
]

#归一化处理
source_test_points = np.array(source_test_points - ave_src)

# 单位权标准差，即所得模型的坐标精度
sigma0 = np.sqrt((b*b).sum() / 2)
# 参数标准差，即所得模型的七参数精度
sigmax = sigma0 * np.sqrt(np.diag(qxx))
print('单位权中误差: %.3f' % (sigma0))
print('参数中误差:')
print(np.round(sigmax,3))
(x2, y2, z2) = ordinationConvert(source_test_points[0, 0], source_test_points[0, 1], source_test_points[0, 2], Args) + ave_tar
print('模型预测结果: ')
print('[(%.3f, %.3f, %.3f)]'%(x2, y2, z2))
print('真实结果: ')
print(target_test_points)

精度检查
目标点与预测值对比

单位权中误差: 162.711
参数中误差:
[ 81.356  81.356  81.356   0.65    0.311   0.191   0.15 ]
模型预测结果:
[(3380987.508, 539711.311, 13.654)]
真实结果:
[(3380972.424, 539704.811, 13.665)]

结束语

坐标系系列的文章在此就告一段落了。坐标系以及与其相关的信息处理一直与我们测绘工作者息息相关。因此我们有必要了解相关概念。祝大家工作顺利！

1: 七参数转换是针对将地理坐标系所对应的空间直角坐标转换为另一坐标系的空间直接坐标。空间直角坐标系的原点位于地球参考椭球的中心，Z轴与地球自转轴平行并指向参考椭球的北极，X轴指向参考椭球的本初子午线，Y轴与X轴和Z轴相互垂直最终构成一个右手系。大地坐标系是以大地基准为基础建立起来的，大地基准又以参考椭球为基础，由此大地坐标系又被称为椭球坐标系。

2: 代码存在中文注释，在Python3中可正常运行，若使用Python2，则需要去掉中文注释来正常运行。

(作者：廖溥昀，孙玉峰，陈关州)