三维空间直角坐标转换中，布尔莎(Bursa)模型、莫洛金斯基(Molodensky)模型和武测模型等被广泛应用。常用的三维坐标转换模型采用最小二乘原理^[1]，将七参数模型误差方程线性化后，建立高斯-马尔可夫迭代模型求解坐标转换七参数，在小角度坐标转换情况下，误差可以忽略，因此这些坐标转换模型及其参数的求解方法适用于小角度的坐标旋转变换^[2-4]。为解决大角度的坐标旋转问题，国内外学者提出了多种三维空间直角坐标转换方法，其中最具代表性的包括方向余弦法^[5]、罗德里格矩阵法^[6]以及旋转矩阵的SVD估计法^[7]等，均取得了较好的效果。然而，上述方法与传统的欧拉角迭代方法一样，是基于坐标转换七参数模型的迭代解法^[8]，在解算过程中，计算繁琐且计算周期长。本文针对三维空间坐标转换模型的七参数求解问题，通过坐标重心化^[9]求解坐标转换的尺度参数，引入对偶四元数^[10]来构造坐标转换旋转矩阵和平移向量，推导基于对偶四元数的三维空间坐标转换模型的直接解法。最后，通过算例以小角度、中角度以及大角度等不同情况验证该算法的正确性和适用性。

1 对偶四元数及其旋转平移 1.1 对偶四元数

对偶四元数的数学形式为^[11]：

(1)

式中，ṗ为对偶四元数的实部，为对偶四元数的对偶部分，ε为对偶运算符且对偶运算符满足ε²=0，ε≠0。

四元数是一个超复数，假设任意一个四元数记为，该四元数由一个实部s₀和3个虚部s₁、s₂、s₃组成，其数学形式为：

(2)

式中，i、j、k为四元数的虚部单位，且满足i²= j²= k²=-1，jk =-kj = i，ik =-ki = j，ij =-ji = k。

1.2 旋转和平移的对偶四元数表达形式

利用对偶四元数描述坐标系间的旋转和平移变换，如图 1所示。为了简单和方便描述，这里把旋转轴和平移轴重合处理，假设向量l为旋转轴，角θ为旋转角，水平距离d表示沿着向量l方向平移的距离。坐标系间的旋转和平移变换(刚体的螺旋运动)可以用一个对偶四元数进行描述。假设该对偶四元数为，其一般形式为：

(3)

式中，四元数仅用于描述定点旋转，四元数，t表示为描述平移的平移轴矢，q表示四元数的虚部。

按照四元数格拉斯曼乘积以及四元数的矩阵形式^[11]，描述旋转的旋转矩阵R和描述平移的向量T可以同时用对偶四元数的实部和对偶部进行表示，其形式为：

(4)

式中，，四元数表示描述平移向量T的四元数，且，即四元数为实部为0的四元数。结合式(2)、(3)和(4)整理得到旋转矩阵R和平移向量T的具体形式为：

(5)

(6)

2 基于对偶四元数的三维空间坐标转换模型 2.1 三维空间坐标转换数学模型

假设三维空间两个坐标系分别为O-X₁Y₁Z₁和O′-X₂Y₂Z₂，空间一点A在坐标系O-X₁Y₁Z₁的坐标[X₁ Y₁ Z₁]^T通过绕Z₁轴旋转α角、绕X₁轴旋转β角以及绕Y₁旋转γ角作3次旋转，最终得到空间坐标系O′-X₂Y₂Z₂的坐标[X₂ Y₂ Z₂]^T。其相应的坐标转换公式可表示为：

(7)

式中，a₁、a₂、a₃分别表示X₁轴与X₂、Y₂、Z₂轴的方向余弦，b₁、b₂、b₃分别表示Y₁轴与X₂、Y₂、Z₂轴的方向余弦，c₁、c₂、c₃分别表示Z₁轴与X₂、Y₂、Z₂轴的方向余弦，λ为坐标转换的尺度参数，[X₀ Y₀ Z₀]^T表示坐标转换的平移向量。将式(7)写成矢量的形式为：

(8)

式中，向量r_T表示坐标系O′-X₂Y₂Z₂的坐标向量，向量r表示坐标系O-X₁Y₁Z₁的坐标向量，矩阵R为两个坐标系之间的旋转矩阵，向量r₀为两个坐标系之间的平移向量。

2.2 坐标转换参数的求解

假设两个不同的坐标系O-X₁Y₁Z₁和O′-X₂Y₂Z₂拥有n个控制点相对应的两套坐标系坐标，分别为[X_1i Y_1i Z_1i]^T和[X_2i Y_2i Z_2i]^T(i=1, …, n)，利用相对应的两套不同坐标系坐标求解两个坐标系间的坐标转换参数。求取两个坐标系间的尺度参数可参考文献[9]。首先分别将两个坐标系的n个控制点坐标重心化。其中，坐标系O-X₁Y₁Z₁的重心化坐标表示为[X_1i, Y_1i, Z_1i](i=1, …, n)，坐标系O′-X₂Y₂Z₂的重心化坐标表示为[X_2i, Y_2i, Z_2i](i=1, …, n)。建立误差方程之后，利用最小二乘原理计算求取坐标系间转换的尺度参数λ，即

(9)

求解出尺度参数之后，再求解姿态参数描述的旋转矩阵R。由式(8)可得，坐标转换的误差方程形式为：

(10)

将误差方程式(10)中n个控制点的平均误差最小，令最小误差函数为：

(11)

假设这里用于描述旋转和平移的对偶四元数为，则将式(8)用对偶四元数的形式表示为：

(12)

将式(12)中对偶四元数描述的坐标转换方程代入式(11)，整理得：

(13)

式中，

(14)

(15)

(16)

采用拉格朗日法对式(13)求极值，然后对四元数ṗ和分别求偏导数，整理得：

(17)

式中，η₁为矩阵N的特征值。要让函数F(λ, R, r₀)最小，则只需使特征值η₁最大，即求取矩阵N的最大特征值，矩阵N的最大特征值对应的特征向量即为四元数ṗ的矩阵形式，将其代入式(5)中即可求出旋转矩阵R；将四元数ṗ代入式(17)，即可求解出四元数；将四元数ṗ和代入式(6)中，即可求出平移向量T。最终实现无需迭代，直接解算两个不同坐标系间坐标转换七参数的解。

3 算例分析

为验证基于对偶四元数的三维空间坐标转换算法的正确性和适用性，设计了坐标转换小角度、中角度以及大角度等3种情况进行模拟实验。假设坐标系O-X₁Y₁Z₁为源坐标系，坐标系O′-X₂Y₂Z₂为目标坐标系，源坐标系中有12个控制点，相对应的目标坐标系中同时也有12个目标坐标系坐标，源坐标系按照绕旋转轴Z轴、X轴和Y轴的顺序依次逆时针旋转，相应的旋转角即为α、β、γ；然后按照[X₀ Y₀ Z₀]^T进行平移变换，同时两个坐标系之间存在一个尺度参数λ。源坐标系12个控制点坐标如表 1所示，不同角度模型的参数如表 2所示。结合表 1和2，通过式(7)计算得到的目标坐标系坐标如表 3所示。

以MATLAB为仿真平台，采用文中所提对偶四元数的坐标转换参数直接解法(算法1)和传统欧拉角迭代解法^[12](算法2)，分别对源坐标系控制点坐标和对应的目标坐标系控制点坐标进行坐标系间转换参数的求解，并从计算结果和计算时间两个方面进行比较(表 4)。

表 4的解算结果比较验证了文中所提基于对偶四元数的三维空间坐标转换直接解法是正确的，算法1的解算结果要优于算法2，且克服了算法2不适用于大旋转角度的劣势。在解算时间上，算法1的计算时间远小于算法2，如果用于区域网平差，效果会更加明显。基于对偶四元数的三维空间坐标转换直接解法无需计算初始值，可以直接计算，且无需线性化和迭代计算，对于大旋转角度的坐标系间的旋转依然能求解出较高精度的坐标转换参数。

4 结语

本文从对偶四元数的定义及其性质入手，描述了对偶四元数在表达坐标系间旋转和平移方面的便利表达形式。针对三维空间坐标转换模型的七参数求解问题，通过坐标重心化求解坐标转换的尺度参数，同时考虑到对偶四元数在描述旋转矩阵和平移向量中的优势，可以有效地克服旋转和平移的耦合误差，引入对偶四元数来构造坐标转换旋转矩阵和平移向量，推导了一种基于对偶四元数的三维空间坐标转换模型的直接解法。该算法与传统的欧拉角迭代算法相比，有效地克服了繁琐的三角运算，无需对旋转参数线性化、无需迭代而直接求解坐标转换的七参数，模型更加简单，计算简便。最后通过仿真实验表明，基于对偶四元数的三维空间坐标转换直接解法具有无需迭代、计算速度快、解算精度高等特点，且适用于大旋转角度的三维空间坐标系间的转换。