空间直角坐标转换的统一模型

吴继忠 王安怡

1 南京工业大学测绘学院，南京市浦珠南路30号，211816

2 大连海洋大学海洋与土木工程学院，大连市黑石礁街52号，116023

根据空间维数的不同，空间直角坐标转换一般分为二维平面和三维空间坐标转换。二维平面坐标转换可直接变换成线性模型，其解算过程容易实现。三维空间直角坐标转换中，在旋转角度很小的情况下可以对旋转矩阵作一定的近似处理，从而转化为线性模型进行求解。当旋转角为大角度时，在小范围测区内7参数之间存在强相关性［1］，更难以处理的是旋转矩阵的各个参数关系复杂，通常进行线性化或迭代运算［2－3］进行求解，这些方法在实际应用中均有一定的条件限制或者使用复杂。本文提出了空间直角坐标转换的统一模型，结合空间直角坐标转换中旋转矩阵的正交特性，将模型求解转换为正交Procrustes问题，实现了二维和三维空间直角坐标转换的直接线性解算。

1 空间直角坐标转换的数学模型

二维平面坐标转换模型可以表示为：

式中，λ为尺度参数；为平移参数；旋转矩阵，φ是旋转角。显然，R2是正交矩阵，满足

三维空间直角坐标的转换模型表示为：

式中，λ为尺度参数；为平移参数；，其中为3个轴向旋转对应的旋转矩阵，其构成可参见相关文献，均为正交矩阵，因此R3也是正交矩阵，满足

对比式（1）和式（2）可知，二维平面坐标转换可以看成是三维空间直角坐标转换的一种特殊形式，二者性质完全相同，区别仅在于矩阵的维数不同。以三维空间直角坐标转换为例，设两套坐标系的公共点数目为m，各点在不同坐标系下的坐标记为：

B为公共点在目标坐标系统下的坐标，A为公共点在原坐标系统下的坐标。显然，每个公共点均满足式（2）的条件，将第i（i＝1，2，…，m）个点代入式（2）并将左右两边矩阵进行转置可得：

其中J＝［1，1，…，1］T为m×1的常数向量。由上述分析，空间直角坐标转换的统一模型可写为如下形式：

式中，k为常量，当k＝2或3时分别对应二维平面坐标转换和三维空间直角坐标转换，m为转换中公共点的数目，R∈Rk×k是转置后的旋转矩阵，是满足正交特性的方阵。式（6）即为空间直角坐标转换的统一模型。

2 正交Procrustes分析

利用公共点求解转换参数实际是要求出参数λ、R、T，使得λAR＋JT尽可能接近B，即求解如下问题：

其中‖·‖F表示Frobenius范数，R满足RTR＝I，式（7）本质上属于正交Procrustes问题。正交Procrustes问题可表述为：给定矩阵A∈Rn×p，B∈Rn×q，其中n＞q，求正交矩阵R∈Rp×q满足‖AR－B‖F＝min。

根据正交Procrustes分析的解算方法［4］，针对式（6）的求解过程如下：记ε＝λAR＋JT－B，根据Lagrange原理构造极值函数：

式中，L为乘数因子矩阵。将式（8）分别对λ、R、T求偏导数：

由式（11）可得：

由式（10）左乘RT得：

由于RTATAR和（L＋LT）／2均是对称矩阵，RTATB－RTATJTT也必须是对称矩阵，即

将式（12）代入式（14）得：

将上式展开，根据对称矩阵的特点同理可得：

记D＝AT（I－JJT／m）B，要使式（16）成立，必须

满足：

将式（17）分别左乘R、右乘RT，得：

将式（18）左右相乘得：

将D进行奇异值分解，svd（D）＝UΣVT，其中U和V均为正交矩阵，代入式（19）得：

故旋转矩阵R＝UVT，将式（12）代入式（9）得：

将R、λ代入式（12）即可得到平移参数T。根据旋转矩阵R的展开式，可以反算出不同轴向上的旋转角，其过程不再赘述。最后用单位权中误差来评定精度，单位权中误差计算公式为：

式中，ε＝λAR＋JT－B；f为自由度，在二维转换和三维转换的情况下分别等于2m－4和3m－7，m为公共点数量。

在计算过程中，首先将公共点坐标按式（6）构建矩阵A、B和常量矩阵J，再根据A、B、J构建矩阵D＝AT（I－JJT／m）B。对矩阵D进行奇异值分解，依次计算得到旋转矩阵R、尺度参数λ和平移参数T，最后进行精度评定。

3 算例分析

3.1 算例1

以文献［5］中的数据为例，5 个公共点在WGS－84和独立坐标系下的坐标如表1所示。为便于比较，取1、2、3号点求转换参数，再利用转换参数将4、5号点转换到地方坐标系，并与其已知值进行比较。

表1 平面坐标转换公共点坐标Tab.1 Cartesian coordinates of common points

利用本文模型和方法，首先计算出旋转矩阵，得到的旋转角为359°01′20.575 4″，进一步得到尺度因子为3.064 875×10－6，最后由旋转矩阵和尺度参数得到平移参数Δx＝－2 465 635.256 m，Δy＝－433 223.055m。利用上述转换参数将4、5号点转换到地方坐标系，结果见表2。

表2 平面坐标转换公共点坐标Tab.1 Cartesian coordinates of common points

由表2可知，本文计算结果与已知结果间仅存在微小的差异，最大坐标较差为3 mm。造成这一差异的原因是点位坐标是实测结果，均包含一定的观测误差。

3.2 算例2

以4个点的三维空间直角坐标为原始坐标，模拟10 套不同转换参数（表3），利用式（2）将4个控制点的原始坐标转换到目标坐标系下，最后利用本文算法计算出转换参数并与其模拟的真值进行比较，得出各个转换参数的计算误差及转换参数求解的中误差，结果列于表4。表4 最后一列为各个方案计算转换参数的中误差。

表3 模拟参数Tab.3 Simulation parameters and solutions

表4 模拟参数的求解误差Tab.4 Errors of simulation parameters

从表4可以看出，本文计算方法对旋转角、尺度参数和平移参数的数值范围没有限制，因而适合于任意条件下的三维空间直角坐标转换。其计算过程为矩阵运算，是一种纯粹的线性计算，突破了以往三维坐标转换方法中对旋转角的限制及计算方法的瓶颈问题，实现了三维坐标转换的直接解算。由于模拟数据的点位坐标不存在误差，在不同方案中，旋转角的计算误差小于10－13，尺度参数计算误差小于10－14，平移参数计算误差小于10－8，计算中误差的数量级也在10－8以下，这些计算误差是计算机运算过程中数据截断误差的累计所造成的，其大小对实际应用而言完全可以忽略不计。

4 结 语

本文提出了空间直角坐标转换的统一方法。首先建立二维平面坐标转换和三维空间直角坐标转换的统一模型，采用正交Procrustes分析对模型进行解算。解算过程不需要参数的先验近似值，也不需要线性化处理和迭代计算，对转换参数的数值范围没有限制，实现了空间直角坐标的转换在表示模型和解算方法上的统一。

［1］王解先.七参数转换中参数之间的相关性［J］.大地测量与地球动力学，2007，27（2）：43－46（Wang Jiexian.Correlations among Parameters in Seven－Parameter Transformation Model［J］.Journal of Geodesy and Geodynamics，2007，27（2）：43－46）

［2］曾文宪，陶本藻.三维坐标转换的非线性模型［J］.武汉大学学报：信息科学版，2003，28（5）：566－568（Zeng Wenxian，Tao Benzao.Non－Linear Adjustment Model of Three－Dimensional Coordinate Transformation［J］.Geomatics and Information Science of Wuhan University，2003，28（5）：566－568）

［3］姚宜斌，黄承猛，李程春，等.一种适用于大角度的三维坐标转换参数求解算法［J］.武汉大学学报：信息科学版，2012，37（3）：253－256（Yao Yibin，Huang Chengmeng，Li Chengchun et al.A New Algorithm for Solution of Transformation Parameters of Big Rotation Angle’s 3D Coordinate［J］.Geomatics and Information Science of Wuhan University，2012，37（3）：253－256）

［4］Crosilla F，Beinat A.Use of Generalised Procrustes Analysis for the Photogrammetric Block Adjustment by Independent Models［J］.ISPRS Journal of Photogrammetry ＆ Remote Sensing，2002，56（2）：195－209

［5］姚宜斌.平面坐标系统相互转换的一种简便算法［J］.测绘信息与工程，2001（1）：1－3（Yao Yibin.A Simple and Convenient Arithmetic for the Transformation of Plane Coordinate System［J］.Journal of Geomatics，2001（1）：1－3