本发明属于图像处理技术领域，更进一步涉及计算机视觉技术领域中的一种已知扫描点位置下的两点点云配准方法。本发明可以应用于3d重建等具体应用场景，针对扫描点空间位置确知的相邻点云，实现对两站三维数据同时旋转配准，从而构建完整的三维点云模型。
背景技术：
点云配准是将多视角扫描到的多站点云转换到统一的坐标系中，得到完整的三维模型。随机抽样一致性算法ransac在点云配准中是一种从样本中准确拟合出数学模型的方法。chen首次将ransac思想应用于点云配准中，之后aiger又提出了基于ransac的四点法。目前提出的基于ransac思想的点云配准方法至少需要三组对应点，需要的样本点越少，取得误匹配点的概率就越低。
北京建筑大学在其拥有的专利技术“一种基于球标靶探测的激光点云数据自动配准方法”(专利申请号：201310746219.1，授权公告号：103646156b)中公开了一种基于球标靶探测的激光点云数据自动配准方法，包括以下步骤：分别从相邻两个站点的激光点云数据中提取球标靶特征；通过寻找相邻点云中由任意三个非共线球标靶球心构成的同名三角形，计算待配准站和参考站的变换关系。该方法通过在场景中设置标靶计算相邻点云的变换坐标，破坏了场景信息。并且该方法通过球形拟合的方式获得球标靶的球心，因此该方法存在的不足之处是，具有一定误差，无法保证配准结果的精确度。
王新、刘永山、朱代春、刘新、张小潮在其发表的论文“三维点云配准算法的研究”(燕山大学学报，2016,40(6)：524-531)中提出了一种基于几何特征和ransac思想的粗配准算法以及基于点的领域几何特征的迭代配准方法。该方法在初始配准过程中采用的三点ransac，每次需要采样3组对应点作为样本点，通过判断由样本点构成的两个三角形是否满足空间几何中全等三角形的性质，将满足该性质的样本点视为合理的初始对应关系，并对合理的样本点进行刚体变换矩阵的求解。由于该方法采样需要的样本点越多，样本点之间的约束关系越复杂，存在的误差就越多，并且使得在采样过程中取得误匹配的概率变高，因此该方法存在的不足之处是，所用的三点ransac在采样过程中取得误匹配点的概率高，存在的误差多，配准效果难以达到理想目标，初始配准精度差。
技术实现要素：
本发明的目的在于克服上述已有技术的不足，提出了一种已知扫描点位置下的两点点云配准方法。
本发明实现的思路是，利用地面三维激光扫描仪获得待扫描物体的点云数据，并使用地面三维激光扫描仪配备的高精度全球定位系统gps设备确定对应站点的扫描点位置，取相邻两站点云作为待配准点云，通过局部特征描述匹配的方法获取待两站待配准点云的对应点集，通过不断在对应点集中随机采样两对对应点作为样本，将样本绕扫描点位置旋转配准，直到满足迭代条件为止，将获得的最优旋转矩阵作用于两站待配准点云，完成配准。
本发明实现的具体步骤包括如下：
(1)获取待扫描物体的点云数据：
用地面三维激光扫描仪获取待扫描物体的点云数据；
(2)预处理点云数据：
(2a)从点云数据中任意选取一个站点的点云数据，将选取的点云数据的原点平移至所选取站点的扫描点位置，将平移后的点云数据作为第一站待配准点云；
(2b)从第一站待配准点云平移前的点云数据附近任意选取一个站点的点云数据，将选取的点云数据的原点平移至所选取站点的扫描点位置，将平移后的点云数据作为第二站待配准点云；
(2c)利用局部特征匹配方法，提取第一站待配准点云中与第二站待配准点云中相匹配的点对，组成两站待配准点云的对应点集，对应点集中包含正确对应点和错误对应点；
(3)利用迭代次数计算函数，计算迭代次数，确定迭代次数的最大值和最小值；
(4)判断当前迭代次数是否小于最小迭代次数或小于上次最大迭代次数，若是，则执行步骤(5)；否则，执行步骤(12)；
(5)计算两站待配准点云的第一次旋转矩阵：
(5a)从两站待配准点云的对应点集中，任意选取一组对应点作为样本a；
(5b)以第一站待配准点云的扫描点位置为球心o1，第一站待配准点云的扫描点位置到样本a中第一个点的欧式距离为半径，做球面s1；
(5c)以第二站待配准点云的扫描点位置为球心o2，第二站待配准点云的扫描点位置到样本a中第二个点的欧式距离为半径，做球面s2；
(5d)判断两个球面的半径之和是否大于两个球心之间的欧式距离，若是，则执行步骤(5e)；否则，执行步骤(5a)；
(5e)在两个球面的相交圆环的边缘上任取一点作为配准点；
(5f)按照下式，计算样本a中的第一个点、球心o1、配准点三点所在平面的法向量：
其中，n1表示样本a中的第一个点、球心o1、配准点三点所在平面的法向量，
(5g)按照下式，计算球心o1和样本a中的第一个点确定的直线与球心o1和配准点确定的直线间的夹角：
其中，α1表示球心o1和样本a中的第一个点确定的直线与球心o1和配准点确定的直线间的夹角，arccos表示反余弦操作，r1表示球面s1的半径，d1表示样本a中第一个点到配准点的欧式距离；
(5h)将法向量n1和夹角α1代入绕任意轴旋转矩阵，得到第一站待配准点云的第一次旋转矩阵；
(5i)按照下式，计算样本a中的第二个点、球心o2、配准点三点所在平面的法向量：
其中，n2表示样本a中的第二个点、球心o2、配准点三点所在平面的法向量，
(5j)通过下式，计算球心o2和样本a中第二个点确定的直线与球心o2和配准点确定的直线间的夹角：
其中，α2表示球心o2和样本a中第二个点确定的直线与球心o2和配准点确定的直线间的夹角，r2表示球面s2的半径，d2表示样本a中第一个点到配准点的欧式距离；
(5k)将法向量n2和夹角α2代入绕任意轴旋转矩阵，得到第二站待配准点云的第一次旋转矩阵；
(6)计算两站待配准点云的第二次旋转矩阵：
(6a)从两站待配准点云的对应点集中，任意选取除样本a之外的一组对应点，作为样本b；
(6b)以球心o1到匹配点的向量为旋转轴，将样本b中的第一个点旋转一个
(6c)以球心o2到匹配点的向量为旋转轴，将样本b中的第二个点旋转
(6d)判断样本b中的第二个点的旋转角度是否大于360°，若是，则将样本b中第二个点的旋转角度置为0后执行步骤(6e)；否则，执行步骤(6c)；
(6e)判断样本b中的第一个点的旋转角度是否大于360°，若是，则执行步骤(6f)；否则，执行步骤(6b)；
(6f)在所有旋转记录中，寻找样本b中两个点的欧式距离最小时对应的样本b中第一个点的旋转角度β1和第二个点的旋转角度β2；
(6g)将球心o1到匹配点的向量和旋转角度β1，代入绕任意轴旋转矩阵中，得到第一站待配准点云的第二次旋转矩阵；
(6h)将球心o2到匹配点的向量和旋转角度β2，代入绕任意轴旋转矩阵中，得到第二站待配准点云的第二次旋转矩阵；
(7)计算当前的两站待配准点云的旋转矩阵：
(7a)将第一站待配准点云的第二次旋转矩阵与第一次旋转矩阵相乘，得到第一站待配准点云的旋转矩阵；
(7b)将第一站待配准点云的第二次旋转矩阵与第一次旋转矩阵相乘，得到第二站待配准点云的旋转矩阵；
(8)计算当前正确对应点个数：
(8a)对两站待配准点云的对应点集进行配准，得到配准后的对应点集：
(8b)将配准后的对应点集中所有满足欧式距离阈值的对应点数，作为当前正确对应点的个数；
(9)利用代价函数，计算当前代价评估值；
(10)更新参数：
(10a)判断当前正确对应点个数是否等于上次正确对应点个数，若是，执行步骤(10b)；否则，执行步骤(10c)；
(10b)判断当前代价评估值是否小于上次代价评估值，若是，用当前代价评估值更新上次代价评估值后执行步骤(10e)；否则，执行步骤(10c)；
(10c)判断当前正确对应点个数是否大于上次正确对应点个数，若是，执行步骤(10d)；否则，执行步骤(11)；
(10d)用当前正确对应点个数更新上次正确对应点个数，利用迭代次数计算函数，计算当前最大迭代次数，用当前最大迭代次数更新上次最大迭代次数；
(10e)用当前两站待配准点云的旋转矩阵更新上次两站待配准点云的旋转矩阵；
(11)判断当前正确对应点个数是否等于两站待配准点云的对应点集包含的点对数，若是，执行步骤(12)；否则，将当前迭代次数加1后执行步骤(4)；
(12)配准两站待配准点云：
按照下式，对两站待配准点云中的每个点进行配准，完成两站待配准点云的配准：
xi＝r1*(pi-o1)+o1
yi＝r2*(qi-o2)+o2
其中，xi表示对第一站待配准点云中第i个点进行配准后得到的点，r1表示第一站待配准点云的旋转矩阵，pi表示第一站待配准点云中第i个点，yi表示对第二站待配准点云中第i个点进行配准后得到的点，r2表示第二站待配准点云的旋转矩阵，qi表示第二站待配准点云中第i个点。
本发明与现有技术相比具有如下优点：
第一，由于本发明在获取待扫描物体的点云数据时，仅利用地面三维激光扫描仪进行采集，克服了现有技术点云配准时，通过在场景中设置的标靶计算相邻点云的变换坐标，破坏场景信息的缺点，使得本发明点云配准方法可以更好的保护场景信息，提高了配准结果真实性。
第二，由于本发明在计算两站待配准点云的旋转矩阵过程中，仅获取了两组样本，克服了现有技术点云配准时，采集的样本点较多，取得误匹配点的概率较大的缺点，使得本发明在点云配准过程中取得误匹配点的概率更小，提高了配准结果可靠性。
第三，由于本发明在计算待配准点云旋转矩阵过程中，利用高精度全球定位系统gps设备确定地面三维激光扫描仪扫描点中心位置作为约束条件，克服了现有技术点云配准时，约束关系复杂的缺点，使得本发明在点云配准过程中约束关系简单精确，提高了配准结果的精确性。
附图说明
图1是本发明的流程图；
图2是本发明两站待配准点云的第一次旋转矩阵的示意图；
图3是本发明两站待配准点云的第二次旋转矩阵的示意图；
图4是本发明计算两站待配准点云的第二次旋转矩阵步骤的流程图。
具体实施方式
下面结合附图对本发明作进一步的描述。
参照图1，本发明的具体实现步骤如下。
步骤1，获取待扫描物体的点云数据。
用地面三维激光扫描仪获取待扫描物体的点云数据。
步骤2，预处理点云数据。
从点云数据中任意选取一个站点的点云数据，将选取的点云数据的原点平移至所选取站点的扫描点位置，将平移后的点云数据作为第一站待配准点云，扫描点位置是指，利用地面三维激光扫描仪配备的高精度全球定位系统gps设备所确定的站点对应的地面三维激光扫描仪的扫描点位置。
从第一站待配准点云平移前的点云数据附近任意选取一个站点的点云数据，将选取的点云数据的原点平移至所选取站点的扫描点位置，将平移后的点云数据作为第二站待配准点云。
利用局部特征匹配方法，提取第一站待配准点云中与第二站待配准点云中相匹配的点对，组成两站待配准点云的对应点集，对应点集中包含正确对应点和错误对应点，局部特征匹配方法是指，使用局部特征描述子对两站待配准点云中的点做局部特征描述，用局部特征对两站待配准点云之间的点进行相似性匹配，将所有的匹配点的集合作为两站待配准点云的对应点集。
步骤3，利用下式，计算迭代次数，确定迭代次数的最大值和最小值：
其中，k表示迭代次数，
步骤4，判断当前迭代次数是否小于最小迭代次数或小于上次最大迭代次数，若是，则执行步骤5；否则，执行步骤12。
上次最大迭代次数是指，在第一次迭代时的最大迭代次数为106，其余迭代过程中的最大迭代次数为当前迭代之前更新的最大迭代次数。
步骤5，计算两站待配准点云的第一次旋转矩阵。
为了不破坏已知扫描点位置信息，在配准第一组对应点时，将对应点绕各自的扫描位置旋转进行配准。
参照图2，对本发明计算两站待配准点云的第一次旋转矩阵过程作进一步的描述。
(5.1)从两站待配准点云的对应点集中，任意选取一组对应点作为样本a。
(5.2)如图2所示，以第一站待配准点云的扫描点位置为球心o1，第一站待配准点云的扫描点位置到样本a中第一个点的欧式距离为半径，做球面s1。
(5.3)如图2所示，以第二站待配准点云的扫描点位置为球心o2，第二站待配准点云的扫描点位置到样本a中第二个点的欧式距离为半径，做球面s2。
(5.4)判断两个球面的半径之和是否大于两个球心之间的欧式距离，若是，则执行本步骤的(5.5)；否则，执行本步骤的(5.1)。
(5.5)当两个球面的半径之和大于两个球心之间的欧式距离时，如图2中所示，两个球面会相交于一个圆环，在两个球面的相交圆环的边缘上任取一点作为配准点。
(5.6)按照下式，计算样本a中的第一个点、球心o1、配准点三点所在平面的法向量：
其中，n1表示样本a中的第一个点、球心o1、配准点三点所在平面的法向量，
(5.7)按照下式，计算球心o1和样本a中的第一个点确定的直线与球心o1和配准点确定的直线间的夹角：
其中，α1表示球心o1和样本a中的第一个点确定的直线与球心o1和配准点确定的直线间的夹角，arccos表示反余弦操作，r1表示球面s1的半径，d1表示样本a中第一个点到配准点的欧式距离。
(5.8)将法向量n1和夹角α1代入绕任意轴旋转矩阵，得到第一站待配准点云的第一次旋转矩阵，任意旋转矩阵如下：
nx2+ny2+nz2＝1
其中，r表示绕任意轴旋转任意角度的旋转矩阵，nx表示旋转轴在x轴方向的投影，ny表示旋转轴在y轴方向的投影，nz表示旋转轴在z轴方向的投影，
(5.9)按照下式，计算样本a中的第二个点、球心o2、配准点三点所在平面的法向量：
其中，n2表示样本a中的第二个点、球心o2、配准点三点所在平面的法向量，
(5.10)通过下式，计算球心o2和样本a中第二个点确定的直线与球心o2和配准点确定的直线间的夹角：
其中，α2表示球心o2和样本a中第二个点确定的直线与球心o2和配准点确定的直线间的夹角，r2表示球面s2的半径，d2表示样本a中第一个点到配准点的欧式距离。
(5.11)将法向量n2和夹角α2代入绕任意轴旋转矩阵，得到第二站待配准点云的第一次旋转矩阵。
步骤6，计算两站待配准点云的第二次旋转矩阵。
为了不破坏已知扫描点位置信息和步骤5中已经配准好的第一组对应点，在配准第二组对应点时，以扫描点位置到匹配点的向量为旋转轴，将第二组对应点配准到一起。参照图3的示意图和图4的流程图，对本发明计算两站待配准点云的第二次旋转矩阵过程作进一步的描述。
(6.1)从两站待配准点云的对应点集中，任意选取除样本a之外的一组对应点，作为样本b。
(6.2)如图3所示，以球心o1到匹配点的向量为旋转轴，将样本b中的第一个点旋转一个
(6.3)如图3所示，以球心o2到匹配点的向量为旋转轴，将样本b中的第二个点旋转
(6.4)判断样本b中的第二个点的旋转角度是否大于360°，若是，则将样本b中第二个点的旋转角度置为0后执行本步骤的(6.5)；否则，执行本步骤的(6.3)。
(6.5)判断样本b中的第一个点的旋转角度是否大于360°，若是，则执行本步骤的(6.6)；否则，执行本步骤的(6.2)。
(6.6)在所有旋转记录中，寻找样本b中两个点的欧式距离最小时对应的样本b中第一个点的旋转角度β1和第二个点的旋转角度β2。
(6.7)将球心o1到匹配点的向量和旋转角度β1，代入绕任意轴旋转矩阵中，得到第一站待配准点云的第二次旋转矩阵。
(6.8)将球心o2到匹配点的向量和旋转角度β2，代入绕任意轴旋转矩阵中，得到第二站待配准点云的第二次旋转矩阵。
步骤7，计算当前的两站待配准点云的旋转矩阵。
将第一站待配准点云的第二次旋转矩阵与第一次旋转矩阵相乘，得到第一站待配准点云的旋转矩阵。将第一站待配准点云的第二次旋转矩阵与第一次旋转矩阵相乘，得到第二站待配准点云的旋转矩阵。
步骤8，计算当前正确对应点个数。
按照下式，对两站待配准点云的对应点集进行配准，得到配准后的对应点集：
m＝{(r1*(pi-o1)+o1,r2*(qi-o2)+o2)}
其中，m表示两站待配准点云的对应点集进行配准后的对应点集，{}表示集合符号，(·,·)表示对应点对符号，r1表示当前第一站待配准点云的旋转矩阵，pi表示在两站待配准点云的对应点集中属于第一站待配准点云的第i个点，r2表示当前第二站待配准点云的旋转矩阵，qi表示在两站待配准点云的对应点集中属于第二站待配准点云的第i个点。
将配准后的对应点集中所有满足欧式距离阈值的对应点数，作为当前正确对应点的个数。
步骤9，利用代价函数，计算当前代价评估值，代价函数如下：
其中，j表示两站待配准点云的完成一次配准的代价评估值，n表示两站待配准点云的对应点集包含的点对数，∑表示求和操作，ri在配准后的对应点集中第i个对应点的欧式距离大于对应点的欧式距离阈值时，ri表示对应点的欧式距离阈值，否则，ri表示配准后的对应点集中第i个对应点的欧式距离；ij在配准后的对应点集中第i个对应点的欧式距离大于对应点的欧式距离阈值时，ij表示0，否则，ij表示1。
步骤10，更新参数。
(10.1)判断当前正确对应点个数是否等于上次正确对应点个数，若是，执行本步骤的(10.2)；否则，执行本步骤的(10.3)。上次正确对应点个数是指，在第一次迭代时的正确对应点个数为2，其余迭代过程中的正确对应点个数为当前迭代之前更新的正确对应点个数。
(10.2)判断当前代价评估值是否小于上次代价评估值，若是，用当前代价评估值更新上次代价评估值后执行本步骤的(10.5)；否则，执行本步骤的(10.3)。上次代价评估值是指，第一次迭代时的代价评估值为106，其余迭代过程中的代价评估值为当前迭代之前更新的代价评估值。
(10.3)判断当前正确对应点个数是否大于上次正确对应点个数，若是，执行本步骤的(10.4)；否则，执行步骤11。
上次正确对应点个数是指，在第一次迭代时的正确对应点个数为2，其余迭代过程中的正确对应点个数为当前迭代之前更新的正确对应点个数。
(10.4)用当前正确对应点个数更新上次正确对应点个数，利用迭代次数计算函数，计算迭代次数，用计算得到的迭代次数更新上次最大迭代次数。
上次正确对应点个数是指，在第一次迭代时的正确对应点个数为2，其余迭代过程中的正确对应点个数为当前迭代之前更新的正确对应点个数。上次最大迭代次数是指，在第一次迭代时的最大迭代次数为106，其余迭代过程中的最大迭代次数为当前迭代之前更新的最大迭代次数。
迭代次数计算函数如下：
其中，k表示迭代次数，
(10.5)用当前两站待配准点云的旋转矩阵更新上次两站待配准点云的旋转矩阵。上次两站待配准点云的旋转矩阵是指，在第一次迭代时的两站待配准点云的旋转矩阵为106，其余迭代过程中的两站待配准点云的旋转矩阵为当前迭代之前更新的两站待配准点云的旋转矩阵。
步骤11，判断当前正确对应点个数是否等于两站待配准点云的对应点集包含的点对数，若是，执行步骤12；否则，将当前迭代次数加1后执行步骤4。
步骤12，配准两站待配准点云。
按照下式，对两站待配准点云中的每个点进行配准，完成两站待配准点云的配准：
xi＝r1*(pi-o1)+o1
yi＝r2*(qi-o2)+o2
其中，xi表示对第一站待配准点云中第i个点进行配准后得到的点，r1表示第一站待配准点云的旋转矩阵，pi表示第一站待配准点云中第i个点，yi表示对第二站待配准点云中第i个点进行配准后得到的点，r2表示第二站待配准点云的旋转矩阵，qi表示第二站待配准点云中第i个点。