鲁达 《测绘学报》
构建与学术的桥梁 拉近与权威的距离
复制链接，关注《测绘学报》抖音！
【测绘学报的个人主页】长按复制此条消息，长按复制打开抖音查看TA的更多作品##7NsBSynuc88##[抖音口令]
本文内容来源于《测绘学报》2020年第8期，审图号GS（2020）4062号。
三维坐标转换公共点最优权值的单纯形搜索算法
郭迎钢 ,李宗春,何华,王志颖
信息工程大学, 河南 郑州 450001
基金项目：国家自然科学基金（41974216）
摘要：为提高三维坐标转换参数的求解质量，本文基于最优化算法提出了一种稳健的公共点加权坐标转换方法。以坐标转换后公共点的点位残差加权平方和最小为目标函数，利用Nelder-Mead单纯形直接搜索算法，寻找公共点坐标分量在解算坐标转换参数时的最优权重组合。以粒子加速器磁铁的准直安装为应用场景，利用模拟数据和实测数据对本文方法进行验证。结果表明：本文方法能够有效降低粗差观测值及质量不佳观测值的权重。与最小二乘、抗差估计等方法相比，本文方法解算结果的点位残差加权平方和更小，坐标转换参数质量更优。本文方法能提高三维坐标转换参数的求解质量，尤其适用于验前精度未知、观测数据质量不佳的情况。
关键词：公共点加权坐标转换，最优化算法，单纯形直接搜索算法，点位残差，抗差估计
引文格式：郭迎钢, 李宗春, 何华, 等. 三维坐标转换公共点最优权值的单纯形搜索算法. 测绘学报，2020，49(8)：1004-1013. DOI: 10.11947.
阅读全文：
全文概述
坐标系的转换与统一是空间测量数据处理最基本的问题之一，在大地测量、工程测量、地理信息系统等领域有着广泛的应用[1-3]。在精密工程和工业安装制造中，需要单一测量系统在不同测站测量或不同的测量系统相互配合来完成测量任务，如粒子加速器等大型科学实验装置的元器件准直安装[4-6]，航天器、飞机、船舶、水轮发电机组等大型构件的精密装配[7-9]，都涉及坐标系的转换和统一。
由于外界环境、通视条件、人为操作等因素，会使公共点坐标含有误差，从而影响坐标转换参数的求解精度。利用已有的公共点坐标求解更高精度的坐标转换参数成为研究重点[10-13]。权代表了平差过程中对观测值的信赖程度。有学者考虑了公共点坐标中的粗差，利用选权迭代法等抗差估计的方法来自动定位和改正粗差[14-16]。有学者同时考虑两套公共点的坐标误差，利用顾及系数矩阵和观测向量误差的整体最小二乘法进行坐标转换[17-20]，但存在求解参数与传统最小二乘估计结果并无二致的问题。文献[21]分析了不同权重对混合整体最小二乘平差结果的影响，但并未给出确定权重的方法。此外，有研究表明人工智能方法能够取得较高的坐标转换精度，不少学者致力于将神经网络、模糊逻辑、遗传编程等人工智能方法应用于二维和三维坐标转换中，研究更加自动化、智能、稳健的坐标转换方法[22-24]。实际应用中，有些公共点的坐标精度较差但并不是大粗差，此类低精度(尤其指观测误差处于2倍至3倍中误差之间)观测值会对坐标转换参数产生一定的影响，现有的抗差估计方法虽能一定程度上抑制此类观测值的影响，但难以达到转换坐标与目标坐标最佳契合的目的。为此，本文提出一种稳健的公共点加权三维坐标转换方法，以公共点的点位残差加权平方和最小为目标函数，以Nelder-Mead单纯形直接搜索算法为工具，实现每个公共点坐标分量在解算坐标转换参数时权重的自动搜索确定，以期提高坐标转换参数的求解质量。
1 三维七参数坐标转换模型
坐标转换模型是两个坐标系之间转换关系的描述，根据维度的不同，可以分为平面坐标转换模型和三维坐标转换模型等；根据坐标转换参数个数的不同，可以分为四参数、七参数等转换模型。此外还有计算机视觉中的仿射变换、透视变换等坐标转换模型。本文以三维七参数转换模型为例，展开公共点加权坐标转换方法的讨论。三维七参数转换模型一般应用于精密工程测量、三维点云数据处理等领域，设Q点在o-xyz坐标系中的坐标为(x,y,z)，在O-XYZ坐标系中的坐标为(X,Y,Z)。两个坐标系之间的转换关系为
(1)
式中
(Tx, Ty, Tz)为平移参数；(α,β,γ)为旋转参数；k为尺度参数。
由于测量误差的存在，Q点从o-xyz经坐标转换后在O-XYZ坐标系中的坐标(X,Y,Z)，与其目标坐标(X0,Y0,Z0)存在差异，则其坐标残差为
(2)
2 算法设计
2.1 定权依据
图 1所示为Q点坐标的估计值和观测值。M1、M2为两个具有代表性的观测值，M1代表精度较低的观测值，M2代表精度较高的观测值，d1、d2分别为M1、M2到估计值的距离，即观测值改正数或点位残差。在无粗差的情况下，点位残差应只包含随机误差的影响，服从正态分布。利用公共点求得坐标转换参数初值后，可以计算得到每个点坐标转换后的点位残差。以所有公共点的点位残差加权平方和最小为目标函数，采用最优化算法不断调节公共点的权重，给M2赋予较大的权重，给M1赋予较小的权重，预期可以获得更优的坐标转换结果。
图 1 控制点的观测值与估计值 Fig. 1 The observed and estimated coordinates of control points
图选项
2.2 权重优化算法
坐标转换残差反映了转换坐标与目标坐标的契合程度，能够用来衡量坐标转换的精度，点位残差加权平方和越小表明利用该组公共点求解的坐标转换参数质量更好。用f(P)表示点位残差加权平方和，则有
(3)
式中，P矩阵内为公共点各坐标分量的权重；V表示点位残差向量。
本文以n个公共点的坐标分量在计算坐标转换参数时的权重为决策变量，以转换后点位残差加权平方和最小为优化目标建立优化模型
(4)
式中，Ω=[0, 1]为P阵元素的取值范围。
为了求解此最优化问题，需要用到优化算法。传统的基于梯度下降、导数理论的优化算法等都需要计算目标函数的梯度，该过程不仅增加了计算工作量，而且还会影响优化算法的收敛性。随着计算机技术的发展，采用不依赖于导数的直接搜索优化算法得到了广泛的应用，如单纯形法、模拟退火和遗传算法等。其中，单纯形法是一种多维直接搜索算法，1947年由文献[25]提出，经由文献[26]加以改善，被应用于各个领域来提高工作效率、节约能量、减小损耗等，取得了很好的应用效果[27]。Nelder-Mead单纯形法通常用于多维空间寻找目标函数的最优值问题，若函数有n个变量，则单纯形为n+1个点集合A={p0,p1, …,pn}构成的几何对象，其中pi(i=0, …,n)为n维列向量。利用单纯形法求解的基本思想是：首先计算出目标函数在单纯形每个点集合位置处的对应函数值，然后对函数值排序，将函数值最大的元素不断迭代替换，同时更新中心点的值，直至单纯形收敛到函数最小值附近时停止迭代。算法流程为：
(1) 构建点集合A={p0,p1, …,pn}，然后计算这些点对应的目标函数值，并对目标函数值进行排序f(p0)≤f(p1)≤…≤f(pn)。
(2) 计算前n个点的中心点
。
(3) 通过映射操作(反射、扩展反射、外收缩、内收缩等)对集合A中的最大元素pn进行替换，并更新中心点p。
(4) 重复上述过程，直到单纯形收敛到函数最小值附近。
本文采用Nelder-Mead单纯形法来搜索使点位残差加权平方和最小的最优权重组合。以三维七参数坐标转换模型为基础，权重优化过程可按如下步骤进行：
(1) 权重初始化，取P=I3n，即所有公共点的坐标分量权重都赋值为1。
(2) 计算坐标转换参数的初值计算，计算点位残差加权平方和。
(3) 以点位残差加权平方和最小为目标函数，利用Nelder-Mead单纯形法搜索最优权重组合。
(4) 当迭代次数达到上限时，输出点位残差加权平方和随迭代次数的变化曲线，并输出最优的权重组合及坐标转换参数。
3 试验与分析
粒子加速器磁铁的准直安装是精密工程测量中坐标系转换与统一的典型应用场景[6]，如图 2所示。首先，需要在磁铁部件的特征位置安装靶标；然后，采用激光跟踪仪等测量系统获取测量坐标系下磁铁部件上特征点的坐标和安装控制网中部分控制点的坐标；最后，利用控制点在安装控制网坐标系和测量坐标系下的两组坐标值进行公共点转换，获得坐标转换参数，即可得到磁铁部件特征点在安装控制网坐标系下的坐标，以进一步指导其准直安装。为了验证本文方法应用于三维七参数坐标转换模型的正确性和合理性，以粒子加速器磁铁准直安装为例进行了仿真试验和实测试验。
图 2 粒子加速器磁铁准直安装 Fig. 2 Magnet alignment of particle accelerator
图选项
3.1 仿真试验
模拟了呈长方体状分布的12个控制点Q1~Q12，其坐标见表 1，空间分布如图 3所示。
表 1 模拟试验控制点坐标Tab. 1 Control points coordinates of simulation experiment
| 点名 | X | Y | Z |
| Q1 | 0.000 | 0.000 | 0.000 |
| Q2 | 3 000.000 | 0.000 | 0.000 |
| Q3 | 6 000.000 | 0.000 | 0.000 |
| Q4 | 0.000 | 4 000.000 | 0.000 |
| Q5 | 3 000.000 | 4 000.000 | 0.000 |
| Q6 | 6 000.000 | 4 000.000 | 0.000 |
| Q7 | 0.000 | 0.000 | 2 000.000 |
| Q8 | 3 000.000 | 0.000 | 2 000.000 |
| Q9 | 6 000.000 | 0.000 | 2 000.000 |
| Q10 | 0.000 | 4 000.000 | 2 000.000 |
| Q11 | 3 000.000 | 4 000.000 | 2 000.000 |
| Q12 | 6 000.000 | 4 000.000 | 2 000.000 |
表选项
图 3 模拟试验控制点空间分布 Fig. 3 Spatial distribution of control points in simulation experiment
图选项
激光跟踪仪采用绝对测距模式，全量程精度优于±10 μm，各测站之间的尺度参数差异非常小。因此，取尺度参数k=1。将这12个控制点绕X轴旋转30°、绕Y轴旋转60°、绕Z轴旋转90°，然后分别沿X、Y、Z方向平移1000 mm，得到Q1′~Q12′的坐标。激光跟踪仪点坐标测量的标称精度为±(15 μm+6 μm/m·D)，模拟数据的测量范围在10 m以内，其测量精度约为±0.075 mm。据此，按照如下3种方案生成模拟数据。
模拟数据1：在Q1′~Q12′的三维坐标分量上各加入[-0.075 mm, 0.075 mm]内的随机误差。
模拟数据2：在模拟数据1的基础上，在Q1′、Q2′的X坐标分量上加入0.500 mm的观测误差，将这两个公共点构造为粗差观测值。
模拟数据3：在模拟数据1的基础上，考虑到2σ< 0.200 mm < 3σ，在Q1′、Q2′的X坐标分量上加入0.200 mm的观测误差，将这两个公共点构造为低精度观测值。
将所有公共点等权处理按最小二乘法求解坐标转换参数，称其为传统方法；利用抗差估计方法(IGGIII)定位并剔除粗差点，称为抗差方法。用这两种方法与本文方法分别处理以上3组数据，用Σ表示点位残差加权平方和，则3种方法处理3组数据的结果见表 2。
表 2 点位残差加权平方和对比Tab. 2 Comparison of sum of weighted squared coordinate residual2
| 处理方法 | 模拟数据1 Σ | 模拟数据2 Σ | 模拟数据3 Σ |
| 传统方法 | 0.013 2 | 0.333 1 | 0.067 0 |
| 抗差方法 | 0.013 2 | 0.012 2 | 0.032 1 |
| 本文方法 | 0.004 0 | 0.012 0 | 0.006 3 |
表选项
由表 2可知，对于模拟数据1，在公共点数据仅包含随机误差时，传统方法与抗差方法的Σ一致，本文方法在每个公共点坐标分量残差平方和的基础上乘以小于1的权因子，因此Σ的数值小于传统方法和抗差方法。对于模拟数据2，当公共点数据包含大的粗差时，传统的最小二乘法受到较为严重的影响，Σ达到了0.333 mm2，与模拟数据1相比增大了25倍；抗差方法和本文方法的Σ值基本一致，且与模拟数据1的传统方法处理结果相近，表明两种方法在面对大的粗差时均表现出了良好的抗差性。对于模拟数据3，当公共点数据包含观测质量较差的观测值(不一定是粗差)时，传统方法受到了较大的影响，与模拟数据1相比Σ值增大了近5倍；抗差方法也受到了一定的影响，Σ值与模拟数据1相比增大了近2倍，且与模拟数据2的Σ值相比也有明显的增大，说明抗差方法在处理观测质量较差的观测值时，抗差性能较差。相比之下，本文方法的Σ值仍保持在较小的水平，且与模拟数据1的Σ值基本一致，表明本文方法在面临观测质量较差的观测值时，仍保持着优良的抗差性。
利用本文方法处理3组数据，迭代次数上限设置为2000次，Σ随迭代次数的变化如图 4所示。
图 4 Σ随迭代次数变化 Fig. 4 The variation ofΣwith iterations
图选项
从图 4可以看出，Nelder-Mead单纯形直接搜索算法保持着较好的搜索性能，随着迭代次数的增加，Σ的值越来越小；随着迭代次数的增加，Σ的变化趋于平缓，输出值逐渐逼近最优结果。各组数据的迭代次数均为1700次左右，计算耗时约17 s。
利用传统方法、抗差方法和本文方法求解的坐标转换参数与模拟参数真值的偏差见表 3。其中
表 3 解算的旋转、平移参数与真值的偏差Tab. 3 The deviation between the calculated value and the truth value of rotation and translation parameters
| 数据 | 方法 | dX/mm | dY/mm | dZ/mm | dQ/mm | dRX/(″) | dRY/(″) | dRZ/(″) |
| 模拟数据1 | 传统方法 | 0.008 1 | -0.003 2 | 0.008 3 | 0.012 0 | -0.36 | 0.72 | -0.36 |
| 抗差方法 | 0.008 1 | -0.003 2 | 0.008 3 | 0.012 0 | -0.36 | 0.72 | -0.36 | |
| 本文方法 | 0.004 7 | -0.008 8 | 0.004 0 | 0.010 7 | -0.36 | 0.36 | 0.00 | |
| 模拟数据2 | 传统方法 | 0.019 5 | 0.136 9 | -0.234 7 | 0.272 4 | -6.48 | -7.20 | -6.84 |
| 抗差方法 | 0.008 2 | -0.003 3 | 0.007 7 | 0.011 7 | -0.72 | 0.72 | -0.72 | |
| 本文方法 | 0.008 8 | 0.001 1 | 0.000 2 | 0.008 9 | -0.72 | 0.36 | -0.72 | |
| 模拟数据3 | 传统方法 | 0.013 0 | 0.052 5 | -0.088 5 | 0.103 7 | -2.88 | -2.52 | -2.88 |
| 抗差方法 | 0.010 5 | 0.020 2 | -0.033 6 | 0.040 6 | -1.80 | -0.72 | -1.80 | |
| 本文方法 | 0.012 0 | 0.003 5 | 0.009 6 | 0.015 8 | -1.80 | 0.72 | -1.80 |
表选项
由表 3可知，对于只含有随机误差的模拟数据1，3种方法求得的旋转平移参数基本一致，精度相当。对于含有粗差的模拟数据2，传统方法难以抵抗粗差的影响，解算的旋转平移参数与真值的偏差较大；抗差方法和本文方法均表现出了较好的抗差性，两种方法求解的旋转平移参数的精度均与仅含有随机误差的模拟数据1的解算精度相当。对于含有较差质量观测值的模拟数据3，传统方法和抗差方法均难以抵抗其影响，解算的旋转平移参数与利用模拟数据1的解算结果相比，质量有些变差，而本文方法受低质量观测值的影响较小，解算结果仍保持较好的质量，说明本文方法在处理有些含有低质量观测值时更能发挥其优势。
本文方法在处理3组数据时，给12个点各坐标分量的赋权情况见表 4。以坐标残差绝对值为横坐标、以权重为纵坐标，做散点图如图 5所示。
表 4 本文方法公共点赋权Tab. 4 Weights table of common points of the proposed method
| 点名 | 模拟数据1 | 模拟数据2 | 模拟数据3 | ||||||
| X | Y | Z | X | Y | Z | X | Y | Z | |
| Q1 | 0.774 | 0.172 | 0.989 | 0.001 | 0.982 | 1.000 | 0.001 | 0.838 | 0.986 |
| Q2 | 0.001 | 0.255 | 0.747 | 0.001 | 0.997 | 0.996 | 0.001 | 0.062 | 0.423 |
| Q3 | 0.001 | 0.683 | 0.973 | 0.478 | 0.966 | 0.973 | 0.999 | 0.016 | 0.450 |
| Q4 | 0.999 | 0.781 | 0.876 | 0.925 | 0.990 | 0.946 | 0.069 | 0.101 | 0.362 |
| Q5 | 0.238 | 0.961 | 0.004 | 0.995 | 0.972 | 0.951 | 0.857 | 0.550 | 0.993 |
| Q6 | 0.001 | 0.389 | 0.335 | 0.997 | 0.999 | 0.905 | 0.987 | 0.201 | 0.485 |
| Q7 | 0.920 | 0.154 | 0.056 | 0.002 | 0.978 | 1.000 | 0.480 | 0.228 | 0.993 |
| Q8 | 0.702 | 0.049 | 0.874 | 0.682 | 0.984 | 0.999 | 0.960 | 0.983 | 0.803 |
| Q9 | 0.658 | 0.340 | 0.071 | 0.991 | 0.961 | 0.962 | 0.894 | 0.897 | 0.001 |
| Q10 | 0.007 | 0.890 | 0.211 | 0.921 | 0.979 | 0.994 | 0.009 | 0.998 | 0.411 |
| Q11 | 0.990 | 0.131 | 0.959 | 0.919 | 0.989 | 0.942 | 0.655 | 0.508 | 0.002 |
| Q12 | 0.244 | 0.138 | 0.996 | 0.935 | 0.972 | 0.822 | 0.848 | 0.028 | 0.415 |
表选项
图 5 坐标分量残差与权重的散点图 Fig. 5 Scatter plot of coordinates residual and weight
图选项
观察表 4和图 5可知，对于模拟数据1，权重在[0, 1]上分布无明显的规律。对于模拟数据2，除了Q1、Q2的X坐标分量权重为0.001之外，其余绝大部分观测值的权重都分布于0.9~1之间，表明本文方法能够自动给粗差观测值赋予小的权值。对于模拟数据3，Q1、Q2的X坐标分量权重为0.001，其余观测值的权重在[0, 1]上分布无明显规律，表明本文方法能够自动给低质量观测值赋予小的权值。总的来看，本文方法以点位残差加权平方和最小为目标，通过迭代不断更新权重组合，最优结果对应的权重组合中，最小权重对应的观测值集合有效包含了粗差和低质量观测值，降低了这些观测值对坐标参数解算结果的影响，提高了解算结果的精度。
为了分析本文方法在面对不同数量、不同位置的粗差和低质量观测值时的表现，设计了以下几组模拟数据。
模拟数据4：在模拟数据1的基础上，在Q1′的X、Y、Z坐标分量上分别加入0.500 mm的观测误差，将该点构造为粗差点。
模拟数据5：在模拟数据1的基础上，在Q1′的X、Y、Z坐标分量上分别加入0.200 mm的观测误差，将该点构造为低质量观测点。
模拟数据6：在模拟数据1的基础上，在Q6′的X、Y、Z坐标分量上分别加入0.500 mm的观测误差，将该点构造为粗差点。
模拟数据7：在模拟数据1的基础上，在Q6′的X、Y、Z坐标分量上分别加入0.200 mm的观测误差，将该点构造为低质量观测点。
模拟数据8：在模拟数据1的基础上，在Q1′、Q2′的X坐标分量上分别加入0.500 mm的观测误差，在Q3′、Q4′的X坐标分量上分别加入0.200 mm的观测误差，构造一种粗差点与低质量观测点同时存在的情形。
模拟数据9：在模拟数据1的基础上，在Q1′、Q2′的X、Y、Z坐标分量上分别加入0.500 mm的观测误差，在Q3′、Q4′的X、Y、Z坐标分量上分别加入0.200 mm的观测误差，构造一种粗差点与低质量观测点同时存在的情形。
模拟数据10：在模拟数据一的基础上，在Q1′~Q5′的X、Y、Z坐标分量上分别加入0.500 mm的观测误差，构造一种近一半公共点为粗差点的情形。
模拟数据11：在模拟数据1的基础上，在Q1′~Q5′的X、Y、Z坐标分量上分别加入0.200 mm的观测误差，构造一种近一半公共点为低质量观测点的情形。
利用本文方法求解以上模拟数据对应的坐标转换参数，用解算的坐标转换参数与参数真值的偏差来衡量解算精度，见表 5。
表 5 模拟数据坐标转换参数的解算精度Tab. 5 The accuracy of coordinate transformation parameters of simulated data
| 数据 | 方法 | dX/mm | dY/mm | dZ/mm | dQ/mm | dRX/(″) | dRY/(″) | dRZ/(″) |
| 模拟数据1 | 本文方法 | 0.004 7 | -0.008 8 | 0.004 0 | 0.010 7 | -0.36 | 0.36 | 0.00 |
| 抗差方法 | 0.008 1 | -0.003 2 | 0.008 3 | 0.012 0 | -0.36 | 0.72 | -0.36 | |
| 模拟数据4 | 本文方法 | -0.007 8 | 0.001 3 | 0.022 5 | 0.023 8 | -0.38 | 1.03 | -0.55 |
| 抗差方法 | 0.001 3 | 0.005 9 | 0.008 1 | 0.010 1 | -0.20 | 0.45 | -0.32 | |
| 模拟数据5 | 本文方法 | -0.011 0 | 0.019 5 | -0.000 8 | 0.022 4 | 0.00 | -0.01 | -0.48 |
| 抗差方法 | 0.018 7 | 0.007 1 | -0.019 4 | 0.027 8 | -2.33 | 0.25 | -1.91 | |
| 模拟数据6 | 本文方法 | 0.008 5 | 0.000 0 | 0.002 5 | 0.008 8 | -0.99 | 0.51 | -0.86 |
| 抗差方法 | 0.012 4 | 0.000 5 | 0.001 0 | 0.012 4 | -1.04 | 0.28 | -0.79 | |
| 模拟数据7 | 本文方法 | 0.013 9 | 0.002 6 | -0.009 4 | 0.017 0 | -1.80 | 0.57 | -0.94 |
| 抗差方法 | 0.016 8 | -0.011 2 | 0.000 1 | 0.020 2 | -2.11 | -0.11 | -1.76 | |
| 模拟数据8 | 本文方法 | 0.009 0 | 0.001 9 | -0.002 3 | 0.009 5 | -0.72 | 0.33 | -0.56 |
| 抗差方法 | -0.000 8 | 0.004 3 | 0.013 9 | 0.014 6 | 4.01 | -1.81 | 4.45 | |
| 模拟数据9 | 本文方法 | -0.064 5 | 0.088 3 | -0.017 6 | 0.110 8 | -2.35 | 1.45 | -5.41 |
| 抗差方法 | -0.073 6 | 0.128 8 | -0.017 8 | 0.149 4 | 3.58 | -3.28 | 1.10 | |
| 模拟数据10 | 本文方法 | -0.527 1 | 0.392 9 | 0.067 5 | 0.660 8 | 38.38 | -22.98 | 26.04 |
| 抗差方法 | -0.420 2 | 0.388 2 | -0.165 7 | 0.595 6 | 16.66 | -17.60 | 6.02 | |
| 模拟数据11 | 本文方法 | -0.168 9 | 0.109 7 | 0.061 3 | 0.210 5 | 14.49 | -8.09 | 13.14 |
| 抗差方法 | -0.059 4 | 0.231 0 | -0.074 5 | 0.249 9 | 1.65 | -5.44 | -0.31 |
表选项
对比模拟数据1与模拟数据4、5、6、7的参数精度，可以看出不同位置的粗差或低精度观测值对于本文方法无显著的影响，本文方法解算参数的外符合精度保持在较高的水平。从模拟数据8的解算结果可以看出，当粗差观测值和低质量观测值共存但数量不多时，本文方法仍能表现出优良的抗差性，解算参数的精度较高。从模拟数据9的解算结果来看，当粗差观测值和低质量观测值总量达到总观测数的1/3时，本文方法与抗差方法的解算结果受到较为严重的影响，解算参数的精度较低。从模拟数据10和11的解算结果来看，当粗差观测值或低质量观测值总量接近总观测值的1/2时，本文方法与抗差估计方法均失去了抗差性。当然，这种测量场景是不会存在的。
3.2 实测数据试验
利用1台Leica AT402激光跟踪仪，在上海光源实验大厅，按自由设站法采集了两个测站的数据，共有10个公共点。测量场景如图 6所示。
图 6 实测试验场景 Fig. 6 Experimental scene
图选项
按传统方法进行公共点坐标转换，每个公共点的坐标残差统计情况见表 6，公共点的点位残差与点位残差中误差σ相比的大小统计情况见表 7。
表 6 公共点坐标残差Tab. 6 Coordinates residual of common points
| 点名 | dX | dY | dZ | 点位 |
| BD620 | -0.176 | -0.065 | -0.031 | 0.190 |
| BD703 | 0.075 | 0.162 | -0.015 | 0.179 |
| BD704 | 0.026 | -0.029 | -0.029 | 0.048 |
| BD705 | 0.112 | -0.085 | 0.061 | 0.154 |
| BD708 | -0.100 | -0.133 | 0.041 | 0.171 |
| BD709 | -0.237 | -0.316 | 0.025 | 0.396 |
| BD711 | -0.200 | -0.249 | 0.000 | 0.320 |
| BD712 | 0.071 | 0.322 | -0.010 | 0.330 |
| BD713 | 0.739 | 1.012 | -0.018 | 1.254 |
| BD720 | -0.311 | -0.618 | -0.024 | 0.692 |
表选项
表 7 公共点点位残差大小情况统计Tab. 7 Statistics table of common points coordinates residual
| ≤1σ | (1~2)σ | (2~3)σ | ＞3σ | σ/mm |
| 8 | 1 | 1 | 0 | 0.504 |
表选项
由表 6和表 7可知，用10个公共点进行公共点转换，σ=0.504 mm，BD713号点的点位残差大小在(2~3)σ之间，可以考虑舍弃该点。删除BD713后，重新进行公共点转换，转换后σ=0.235 mm，此时BD712为点位残差最大的点，点位残差达0.557 mm。同理，BD712的点位残差仍在(2~3)σ之间，继续删除BD712后，再次进行公共点转换，得到其余点的点位残差及其统计情况见表 8和表 9。
表 8 删除2点后的公共点残差Tab. 8 Coordinates residual of common points after 2 points are deleted
| 点名 | dX | dY | dZ | 点位 |
| BD620 | -0.243 | -0.161 | -0.008 | 0.291 |
| BD703 | 0.117 | 0.173 | -0.029 | 0.210 |
| BD704 | 0.056 | 0.024 | -0.022 | 0.064 |
| BD705 | 0.178 | 0.024 | 0.062 | 0.190 |
| BD708 | 0.068 | 0.062 | 0.003 | 0.092 |
| BD709 | -0.095 | -0.079 | 0.017 | 0.124 |
| BD711 | -0.043 | 0.022 | -0.007 | 0.048 |
| BD720 | -0.038 | -0.064 | -0.016 | 0.076 |
表选项
表 9 删除2点后的公共点点位残差大小情况统计Tab. 9 Statistics table of common points coordinates residual after 2 points are deleted
| ≤1σ | (1~2)σ | (2~3)σ | ＞3σ | σ/mm |
| 5 | 3 | 0 | 0 | 0.159 |
表选项
由表 6~表 8可以看出，依次删去BD713、BD712之后，σ由0.504 mm减小至0.235 mm，再减小至0.159 mm，对应的点位残差加权平方和Σ由2.609 mm2减小至2.502 mm2，再减小至0.212 mm2。表明通过人工筛选公共点、依次删去最大点位残差对应的公共点后，单位权中误差σ及点位残差平方和Σ均明显变小，提高了坐标系转换的契合程度。
利用本文方法对实测数据进行公共点坐标转换，得到Σ随迭代次数的变化如图 7所示。
图 7 实测试验Σ随迭代次数变化 Fig. 7 Variation ofΣof the measured data along with iterations
图选项
从图 7可以看出，本文方法在处理实测数据时与处理模拟数据时的表现一致，都保持着良好的搜索性能，随着迭代次数的增大，Σ的数值逐渐变小，逐渐逼近最优结果。
在利用抗差估计(IGGⅢ)的方法进行坐标转换时，验前单位权中误差σ0的取值要合理。σ0取值不同，解算结果及对应的Σ也会有不同。利用传统方法、抗差方法与本文方法处理实测数据，抗差方法分别以σ0=±1.0 mm、σ0=±0.1 mm两种情况进行计算，处理结果见表 10。
表 10 3种方法实测数据结果对比Tab. 10 Comparison of measured data adjustment results with three methods
| 数据处理策略 | Σ/mm2 | Tx/mm | Ty/mm | Tz/mm | RX/(″) | RY/(″) | (RZ-63.4970°)/(″) | |
| 传统方法 | 所有公共点 | 2.609 | -8 938.035 | 6 995.110 | 2.980 | 1.5 | 1.1 | 0.8 |
| 人工筛选公共点 | 0.212 | -8 937.915 | 6 995.204 | 2.952 | 1.3 | -0.5 | 2.5 | |
| 抗差方法 | σ0=±1.0 mm | 2.609 | -8 937.035 | 6 995.110 | 2.980 | 1.5 | 1.1 | 0.8 |
| σ0=±0.1 mm | 0.286 | -8 937.931 | 6 995.216 | 2.964 | 1.3 | 0.0 | 3.1 | |
| 本文方法 | 0.032 | -8 937.922 | 6 995.191 | 2.958 | 1.2 | -0.2 | 2.3 | |
表选项
由表 10可以看出，按传统方法将所有公共点等权处理，解算结果的Σ值达2.609 mm2，表明求解的坐标转换参数质量较差；人工删除点位残差较大的公共点后，Σ减小为0.212 mm2，获得了较高质量的坐标转换参数，但人工筛选公共点的效率较低，且容易误判。抗差估计方法对验前单位权中误差的取值要求较高，当σ0取值较大时(如σ0=±1.0 mm的情况)，起不到剔除质量不佳观测值的作用，此时的解算结果与传统方法采用所有公共点的解算结果一致，精度较差；当σ0取值合理时(如σ0=±0.1 mm的情况)，Σ的值明显变小，且解算的平移旋转参数也与人工筛选公共点的结果更为相近。相比之下，本文方法解算结果的Σ值最小，且解算的平移旋转参数与人工筛选公共点、σ0取值合理时抗差估计的解算结果都十分相近，充分验证了本文方法的正确性和有效性。
4 结 论
本文基于最优化算法提出了一种稳健的公共点加权三维坐标转换方法。以坐标转换后公共点的点位残差加权平方和最小为目标函数，利用Nelder-Mead单纯形直接搜索算法寻找公共点的最优权重组合，来求解更高质量的坐标转换参数。通过模拟试验和实测试验对本文方法进行了验证，结果表明：与传统的最小二乘法坐标转换方法相比，本文方法的点位残差加权平方和明显变小，具有抗差、自动化赋权等优点；与抗差估计坐标转换方法相比，本文方法不仅克服了验前单位权中误差的影响，还能够自动给低质量观测值赋予小的权值，使转换坐标与目标坐标更为契合。本文方法能够提高三维坐标转换参数的求解质量，尤其适用于验前精度未知、观测数据质量不佳的情况。同时，本文方法的思路也可推广应用于三维仿射变换、平面坐标转换等坐标转换模型中，可以提高转换坐标与目标坐标的契合程度。
不足之处是，当粗差或低质量观测值个数达总观测数的1/3时，本文方法受到较为严重的影响；当粗差或低质量观测值个数接近总观测数的1/2时，本文方法基本丧失抗差性。此外，由于解算过程中利用最优化算法迭代求解，因此该方法的计算效率较低。未来将对该方法关于转换矩阵的几何结构、协方差阵、误差的敏感程度和响应程度继续开展研究，同时也将采用多种公共点分布数据进行更加全面的测试。
作者简介
第一作者简介：郭迎钢(1992-), 男, 博士生, 研究方向为精密工程测量。E-mail:fariel_gyg@163.com
《测绘学报（英文版）》（JGGS）专刊征稿：LiDAR数据处理
论文推荐 | 王旭, 柴洪洲,王昶,种洋:优选小波函数的小波神经网络预报GPS卫星钟差
资讯 | 中国内地大学ESI排名出炉：313所高校上榜
人物 | 无止境的科学追求——记大地测量学家陈俊勇院士
权威 | 专业 | 学术 | 前沿
微信、抖音小视频投稿邮箱 | song_qi_fan@163.com
欢迎加入《测绘学报》作者QQ群： 751717395
进群请备注：姓名+单位+稿件编号