CN111369606B - 基于无控制扫描点云的文物对象高精度微变形监测方法 - Google Patents

Abstract

本发明公开了一种基于无控制扫描点云的文物对象高精度微变形监测方法，包括：步骤一、在一个时期内对文物对象连续进行两次数据采集，获得两个数据：应用ICP算法进行配准，得出文物对象表面未发生变化的情况下关节臂扫描点云的配准精度σ₀；随后求出变形度量D的中误差σ_d；步骤二、在另一个时期内对文物对象再进行数据采集，将两期第一次数据进行格网迭代法配准，使其配准精度小于或等于σ₀；步骤三、用两期的第一次数据进行变形度量D的计算，将D值变化超出置信区间(‑mσ_d～+mσ_d)的值视为文物对象的变化值。本发明在不布设标靶不损害本体的前提下，评价能够监测出的最小变化量及变化量分布，为文物保护提供科学的数据参考。

Description

基于无控制扫描点云的文物对象高精度微变形监测方法

技术领域

本发明属于建筑领域，涉及一种基于无控制扫描点云的文物对象高精度微变形监测方法。

背景技术

造像文物表面风化脱落是引起造像消亡的主要因素，由于文物表面变化往往是片状连续性的微小变化，基于监测点的传统变形监测方法不能捕捉其表面整体变形情况，应用三维激光扫描技术进行变形监测无疑是理想的方法。然而，基于标靶点的激光扫描点云监测也不适用于文物对象的微变形监测，因为其表面或者其赋存环境不允许布设标靶点；所以需要研究基于无控制激光扫描点云的文物微变形监测方法，建立基于激光三维激光扫描技术的监测评价体系。

文物表面风化脱落变形属于微小变形，并且变形区域可能随机分布在对象表面的任意位置，故扫描点云需要具有一定的密度和完整性。普通的地面激光雷达的点云分辨率、精度、完整性等还不能达到对上述微小随机变形进行监测的要求，相比之下，高分辨率高精度的关节臂扫描仪可以满足此类对象的监测需求，然而由于文物对象表面及其赋存环境不能安置标靶，所以研究无标靶关节臂点云的高精度配准方法是此类微小变形的关键问题。目前，无标靶激光扫描点云配准主要是提取对象点云的特征点、线、面，但该方法不适用于非同期有变形的扫描点云。看上去似乎无标靶有变形数据配准是个无解的问题。本发明因此而产生。

发明内容

本发明的一个目的是解决至少上述问题和/或缺陷，并提供至少后面将说明的优点。

本发明还有一个目的是提供一种基于无控制扫描点云的文物对象高精度微变形监测方法。

为此，本发明提供的技术方案为：

一种基于无控制扫描点云的文物对象高精度微变形监测方法，包括如下步骤：

步骤一、在一个时期内对文物对象进行连续两次数据采集，获得两个数据：第一期第一次数据和第一期第二次数据，应用ICP(Iterative Closest Point)算法对所述第一期第一次数据和第一期第二次数据进行配准，得出文物对象表面未发生变化的情况下关节臂扫描点云的配准精度σ₀，ICP算法精细配准后，基于点云数据构建不规则三角网，并逐点计算同期数据间的变形度量D，D表示第一期第二次数据的三角网顶点到第一期第一次数据的三角网面片的有向距离，并求出变形度量的中误差σ_d；

步骤二、在另一个时期内对所述文物对象再进行数据采集，获得第二期第一次数据，将所述第二期第一次数据和所述第一期第一次数据进行格网迭代法配准，使其配准精度小于或等于所述配准精度σ₀，其中，所述第一期第一次数据、第二次数据和所述第二期第一次数据均为基于关节臂扫描的精细点云数据；

步骤三、对满足配准精度的第二期第一次数据进行变形度量的计算：计算变形度量时，基于点云数据构建不规则三角网，以第二期第一次数据的三角网顶点到第一期第一次数据的三角网面片的距离作为不同期数据间的变形度量D，将所述不同期数据间的变形度量D值与所述同期数据间的变形度量D值相比，若其变化超出置信区间(-mσ_d～+mσ_d D_μ)的值视为所述文物对象的变化值，其中，m为不大于3的正整数。

优选的是，所述的基于无控制扫描点云的文物对象高精度微变形监测方法中，步骤二中，所述配准包括如下步骤：

2.1)首先应用点云同名特征计算所述第一期第一次数据和第二期第一次数据点云的粗配准参数，并运用空间相似变换模型把两期点云数据进行坐标统一；

2.2)以点云重心为原点，按照一定的步长对上述两期点云数据进行空间格网划分，并计算同名格网点云相似度，以同名格网点云相似度作为甄选无变形格网点云的基础；

其中，以D1形状函数距离直方图这一特征的相似度衡量所述同名格网点云相似度，D1为格网内各点到格网点云重心距离，格网点云相似度计算公式如(11)(12)和(13)所示：

H_k为参考点云第k个格网内点云D1形状函数距离直方图，H_k(i)表示H_k直方图第i个距离间隔所包含点数，

表示H_k直方图各距离间隔内的点数均值；H_k′为待检测点云第k′个格网内点云D1形状函数距离直方图，H_k′(i)表示H_k′直方图第i个距离间隔所包含点数，

表示H_k′直方图各距离间隔内的点数均值；

式中，d(H_k,H_k′)表示对应第k对格网点云的D1形状函数距离直方图相似度，n为直方图横向距离间隔数；

2.3)对上述两期的同名格网点云进行迭代配准，甄选出无变形格网点云，之后对上述两期点云数据进行高精度配准，具体步骤为：

2.3.1)针对所述第一期第一次数据和第二期第一数据每个象限内格网D1距离直方图相似度排序结果，选出相似性大的前n个格网数据作为配准基元，其中，n≥3；

2.3.2)对选出的8*n个格网点云数据进行ICP配准，如果配准精度没有达到σ₀阈值要求，去掉8*n个格网点云中相似度最低的格网数据，再进行所述ICP配准；

2.3.3)按照2.3.2)步骤的约束条件进行格网迭代配准，当某个象限中的格网数目变为零时仍然没有达到σ₀的配准精度，则缩小格网划分的粒度，再次进行2.3.1))-2.3.2)，直到迭代满足条件为止。

优选的是，所述的基于无控制扫描点云的文物对象高精度微变形监测方法中，步骤一中，ICP算法包括如下计算公式：

式中：式(1)所示为配准模型，公式中f(H)为目标函数，反映在变换矩阵H下，两个控制点集的差异程度；Q_i和P_i为待配准点云，R为旋转矩阵，T为平移矩阵，n为参与配准的点对数量，σ₀即为配准精度，表示配准后数据最邻近点对的欧氏距离中误差；在步骤一中，Q_i为第一期第一次点云数据，P_i为第一期第二次点云数据；ICP算法步骤如下：

1)在第一期第二次数据P中取点P_i∈P；

2)寻找第一期第一次数据中的对应点Q_i∈Q，要求Q_i与P_i欧氏距离最小；

3)计算旋转矩阵R和平移矩阵T，使得f(H)最小；

4)对P_i使用上一步求取的旋转矩阵R和平移矩阵T进行六参数变换，得到新的对应点集P′_i，用以代替P_i参与下一次迭代计算；

5)计算当前旋转矩阵中的角度变化值

Δω,Δκ；

6)角度变化值

Δω,Δκ均小于0.1秒时，则停止迭代计算，否则回到第一步，直至收敛为止；

7)根据上述公式(1)和(2)计算f(H)、σ₀。

优选的是，所述的基于无控制扫描点云的文物对象高精度微变形监测方法中，所述同期数据间的变形度量D和变形度量的中误差σ_d的计算步骤如下：

1)以第一期第一次数据T_ij为参考数据，以ICP配准后的第一期第二次数据T_i,j+1为待检测数据；

2)以T_i,j+1的三角网顶点K为出发点，并以该点的法线

作为方向构建射线方程R(t)：

式中，K为待检测数据T_i,j+1的三角网顶点坐标；t为系数，

为K点的单位法向量。

3)假定从K点出发的射线R(t)与参考数据T_ij交于点K′，且K′位于参考数据T_ij某个三角形ΔABC上，那么这个交点就可以用以下两个方程来表示：

K′(u,v)＝(1-u-v)*A+u*B+v*C

式中，K′为从K点出发的射线R(t)与参考数据T_ij的交点坐标；A，B，C分别代表三角形ΔABC的三个顶点坐标；(1-u-v)，u，v代表顶点A，B，C坐标的系数；

联立公式得：

求得参数u,v,t后，根据公式求出K′点的坐标值。

4)计算所述变形度量D：

|D|＝||K-K′||         (7)

式中，变形度量D为待检测三角网

的三角形顶点K到参考三角网

面片的有向距离。当解算的参数t值为正数时，有向距离为负，表示缩进；当解算的参数t值为负数时，有向距离为正，表示凸出。实际上，D的值与参数值-t相同。

计算出各顶点的D值后，根据公式并求出变形度量中误差σ_d；

鉴于D为离散型随机变量，且E(D)＝0则：

式中，σ_d为变形度量的中误差，是判断不同期次数据是否发生变化的指标；D_i为各顶点对应的变形度量值。

优选的是，所述的基于无控制扫描点云的文物对象高精度微变形监测方法中，步骤三中，计算变形度量值时，基于点云数据构建不规则三角网，以第二期第一次数据的三角网顶点到第一期第一次数据的三角网面片的有向距离作为变形度量，具体方法包括如下步骤：

4.1)以第一期第一次数据T_ij为参考数据，以ICP配准后的第二期第一次数据T_i+1,j为待检测数据；

4.2)以T_i+1,j的三角网顶点K为出发点，并以该点的法线

作为方向构建射线方程R(t)：

式中，K为待检测数据T_i+1,j的三角网顶点坐标；t为系数，

为K点的单位法向量。

4.3)假定从K点出发的射线R(t)与参考数据T_ij交于点K′，且K′位于参考数据T_ijT_ij某个三角形ΔABC上，那么这个交点就可以用以下两个方程来表示：

式中，K′为从K点出发的射线R(t)与参考数据T_ij的交点坐标；A，B，C分别代表三角形ΔABC的三个顶点坐标；(1-u-v)，u，v代表顶点A，B，C坐标的系数。

联立公式得：

求得参数u,v,t后，根据公式求出K′点的坐标值。

4.4)计算所述变形度量D：

|D|＝||K-K′||            (19)

式中，变形度量D为待检测三角网T_i+1,j的三角形顶点K到参考三角网T_ij面片的有向距离。当解算的参数t值为正数时，有向距离为负，表示缩进；当解算的参数t值为负数时，有向距离为正，表示凸出。实际上，D的值与参数值-t相同。

优选的是，所述的基于无控制扫描点云的文物对象高精度微变形监测方法中，m＝2。

优选的是，所述的基于无控制扫描点云的文物对象高精度微变形监测方法中，应用关节臂三维激光扫描系统进行数据的采集，设置扫描对象点云分辨率0.2-0.3mm，单次扫描点云个数约800万。

优选的是，所述的基于无控制扫描点云的文物对象高精度微变形监测方法中，所述一个时期为3个日历天。

本发明至少包括以下有益效果：

本方法能够在不布设标靶不损害本体的前提下，配准包含潜在变化量的邻期点云数据，基于配准结果计算变化量，并评价能够监测出的最小变化量及变化量及分布，为文物保护提供科学的数据参考。

本发明的其它优点、目标和特征将部分通过下面的说明体现，部分还将通过对本发明的研究和实践而为本领域的技术人员所理解。

附图说明

图1为本发明其中一个实施例中的点云数据的配准的流程示意图；

图2为本发明其中一个实施例中单个格网点云D1形状函数示意图；

图3为本发明其中一个实施例中微变形检测实验对象图，其中A为造像上半身，B为造像整体，C为造像下半身；

图4为本发明其中一个实施例中非同期两次关节臂扫描的点云数据图；

图5为本发明其中一个实施例中同期点云数据ICP配准结果图；

图6为本发明其中一个实施例中以5cm基准格网划分点云数据图之一；

图7为本发明其中一个实施例中以5cm基准格网划分的甄选结果图之一；

图8为本发明其中一个实施例中以5cm基准格网划分的ICP迭代计算的去除过程及结果图之一；

图9为本发明其中一个实施例中以2.5cm基准格网划分的甄选结果图之一；

图10为本发明其中一个实施例中以2.5cm基准格网划分的ICP迭代计算的去除过程及结果图之一；

图11本发明其中一个实施例中非同期数据微变形量检测结果图；

图12为本发明其中一个实施例中非同期数据微变形量检测局部结果图。

具体实施方式

下面结合附图对本发明做进一步的详细说明，以令本领域技术人员参照说明书文字能够据以实施。

应当理解，本文所使用的诸如“具有”、“包含”以及“包括”术语并不配出一个或多个其它元件或其组合的存在或添加。

申请人分析到，文物表面的风化脱落等变形是随着时间渐变式发展变化的，并且由于文物工艺的非标准性及文物赋存环境等因素的影响，文物表面变形是具有随机性的，再加之获取监测数据的时间间隔、根据其变化速度认为确定：在两期监测数据之间，一定会存在不变的区域，其位置和范围亦具有一定的随机性。鉴于此重要的发现，结合上述对变形监测中配准问题的一系列分析，本发明提出了一种基于空间格网迭代，甄选不变格网作为配准基元，同时顾及配准基元几何强度的有变形点云无标靶高精度监测方法。主要内容包括空间格网划分、配准格网选取策略设定以及迭代目标函数与特征度量选择等问题。该方法能够在不布设标靶不损害本体的前提下，配准包含潜在变化量的邻期点云数据，应用配准结果计算变化量，并评价能够监测出的最小变化量及变化量及分布，为文物保护提供科学的数据参考。

如图1所示，本发明提供一种基于无控制扫描点云的文物对象高精度微变形监测方法，其特征在于，包括如下步骤：

步骤一、在一个时期内对文物对象进行连续两次数据采集，获得两个数据：第一期第一次数据和第一期第二次数据，应用ICP(Iterative Closest Point)算法对所述第一期第一次数据和第一期第二次数据进行配准，得出文物对象表面未发生变化的情况下关节臂扫描点云的配准精度σ₀；同时逐点计算同期数据间的变形度量D，统计D值概率密度分布情况并求出变形度量中误差σ_d；作为优选，计算公式如下：

式中：配准模型如公式(1)所示，公式中f(H)为目标函数，反映在变换矩阵H下，两个控制点集的差异程度；Q_i为第一期第一次点云数据，P_i为第一期第二次点云数据，R为旋转矩阵，T为平移矩阵，n为参与配准的点对数量，σ₀即为配准精度，表示配准后数据最邻近点对的欧氏距离中误差。

ICP算法步骤如下：

1)在第一期第二次数据P中取点P_i∈P；

2)寻找第一期第一次数据中的对应点Q_i∈Q，要求Q_i与P_i欧氏距离最小；

3)计算旋转矩阵R和平移矩阵T，使得f(H)最小；

4)对P_i使用上一步求取的旋转矩阵R和平移矩阵T进行六参数变换，得到新的对应点集P′_i，用以代替P_i参与下一次迭代计算；

5)计算当前旋转矩阵中的角度变化值

Δω,Δκ；

6)角度变化值

Δω,Δκ均小于0.1秒时，则停止迭代计算，否则回到第一步，直至收敛为止；

7)根据上述公式(1)和(2)计算f(H)、σ₀。

ICP算法精细配准后，基于点云数据构建不规则三角网，并逐点计算同期数据间的变形度量D，D表示第一期第二次数据的三角网顶点到第一期第一次数据的三角网面片的有向距离。

计算步骤如下：

1)以第一期第一次数据T_ij为参考数据，以ICP配准后的第一期第二次数据T_i,j+1为待检测数据；

2)以T_i,j+1的三角网顶点K为出发点，并以该点的法线

作为方向构建射线方程R(t)：

式中，K为待检测数据T_i,j+1的三角网顶点坐标；t为系数，

为K点的单位法向量。

3)假定从K点出发的射线R(t)与参考数据T_ij交于点K′，且K′位于参考数据T_ij某个三角形ΔABC上，那么这个交点就可以用以下两个方程来表示：

式中，K′为从K点出发的射线R(t)与参考数据T_ij的交点坐标；A，B，C分别代表三角形ΔABC的三个顶点坐标；(1-u-v)，u，v代表顶点A，B，C坐标的系数；

联立公式得：

求得参数u,v,t后，根据公式求出K′点的坐标值。

4)计算所述变形度量D：

|D|＝||K-K′||                 (7)

式中，变形度量D为待检测三角网

的三角形顶点K到参考三角网

面片的有向距离。当解算的参数t值为正数时，有向距离为负，表示缩进；当解算的参数t值为负数时，有向距离为正，表示凸出。实际上，D的值与参数值-t相同。

计算出各顶点的D值后，根据公式并求出变形度量中误差σ_d；

鉴于D为离散型随机变量，且E(D)＝0则：

式中，σ_d为变形度量的中误差，是判断不同期次数据是否发生变化的指标；D_i为各顶点对应的变形度量值。

步骤二、在另一个时期内对所述文物对象再进行数据采集，获得第二期第一次数据，将所述第二期第一次数据和所述第一期第一次数据进行格网迭代配准，使其配准精度小于或等于所述配准精度σ₀，其中，所述第一期第一次数据、第二次数据和所述第二期第一次数据均为基于关节臂扫描的精细点云数据；

步骤三、对满足配准精度的第二期第一次数据进行变形度量的计算：计算变形度量时，基于点云数据构建不规则三角网，以第二期第一次数据的三角网顶点到第一期第一次数据的三角网面片的有向距离作为不同期数据间的变形度量D，将所述不同期数据间的变形度量D值与所述同期数据间的变形度量D值相比，若变化超出置信区间(-mσ_d～+mσ_dD_μ)的值视为所述文物对象的变化值，其中，m为不大于3的正整数。

在本发明的其中一个实施例中，作为优选，步骤二中，所述配准包括如下步骤：

2.1)首先应用点云同名特征计算所述第一期第一次数据和第二期第一次数据点云的粗配准参数，并运用空间相似变换模型把两期点云数据进行坐标统一；

2.2)以点云重心为原点，按照一定的步长对上述两期点云数据进行空间格网划分，并计算同名格网点云相似度，以同名格网点云相似度作为甄选无变形格网点云的基础；

其中，以D1形状函数距离直方图这一特征的相似度衡量所述同名格网点云相似度，计算公式如(11)(12)和(13)所示：

H_k为参考点云第k个格网内点云D1形状函数距离直方图，H_k(i)表示H_k直方图第i个距离间隔所包含点数，

表示H_k直方图各距离间隔内的点数均值；H_k′为待检测点云第k′个格网内点云D1形状函数距离直方图，H_k′(i)表示H_k′直方图第i个距离间隔所包含点数，

表示H_k′直方图各距离间隔内的点数均值。

式中，d(H_k,H_k′)表示对应第k对格网点云的D1形状函数距离直方图相似度，n为直方图横向距离间隔数。

2.3)对上述两期的同名格网点云进行迭代配准，甄选出无变形格网点云，之后对上述两期点云数据进行高精度配准，具体步骤为：

2.3.1)针对所述第一期第一次数据和第二期第一数据每个象限内格网D1距离直方图相似度排序结果，选出相似性大的前n个格网数据作为配准基元，其中，n≥3；

2.3.2)对选出的8*n个格网点云数据进行ICP配准，如果配准精度没有达到σ₀阈值要求，去掉8*n个格网点云中相似度最低的格网数据，再进行所述ICP配准；

2.3.3)按照2.3.2)步骤的约束条件进行格网迭代配准，当某个象限中的格网数目变为零时仍然没有达到σ₀的配准精度，则缩小格网划分的粒度，再次进行2.3.1)-2.3.2)，直到迭代满足条件为止。

在本发明的其中一个实施例中，作为优选，步骤三中，计算变形度量值时，基于点云数据构建不规则三角网，以第二期第一次数据的三角网顶点到第一期第一次数据的三角网面片的有向距离作为变形度量，具体方法包括如下步骤(与步骤一种变形量D的计算方式相同，此处以第二期第一次数据为待检测数据)：

4.1)以第一期第一次数据T_ij为参考数据，以ICP配准后的第二期第一次数据T_i+1,j为待检测数据；

4.2)以T_i+1,j的三角网顶点K为出发点，并以该点的法线

作为方向构建射线方程R(t)：

式中，K为待检测数据T_i+1,j的三角网顶点坐标；t为系数，

为K点的单位法向量。

4.3)假定从K点出发的射线R(t)与参考数据T_ij交于点K′，且K′位于参考数据T_ijT_ij某个三角形ΔABC上，那么这个交点就可以用以下两个方程来表示：

式中，K′为从K点出发的射线R(t)与参考数据T_ij的交点坐标；A，B，C分别代表三角形ΔABC的三个顶点坐标；(1-u-v)，u，v代表顶点A，B，C坐标的系数。

联立公式得：

求得参数u,v,t后，根据公式求出K′点的坐标值。

4.4)计算所述变形度量D：

|D|＝||K-K′||                      (19)

式中，变形度量D为待检测三角网T_i+1,j的三角形顶点K到参考三角网T_ij面片的有向距离。当解算的参数t值为正数时，有向距离为负，表示缩进；当解算的参数t值为负数时，有向距离为正，表示凸出。实际上，D的值与参数值-t相同。

在本发明的其中一个实施例中，作为优选，m＝2。

在本发明的其中一个实施例中，作为优选，应用关节臂三维激光扫描系统进行数据的采集，设置扫描对象点云分辨率0.2-0.3mm，单次扫描点云个数约800万。

在本发明的其中一个实施例中，作为优选，所述一个时期为3个日历天。

为使本领域技术更好地理解本发明的技术方案，现提供如下的实施例进行说明：

本文主要方法包括：格网点云相似度计算、格网迭代高精度配准、微变形检测与评价。在数据获取方面，由于造像的风化残损属于微小变形量，所以应用精度和密度较高的关节臂激光扫描仪获取造型表面点云数据。在同等观测条件下，获取造像同期两次关节臂点云数据，把同期两次数据的配准精度作为关节臂点云系统误差、环境误差、偶然误差等综合因素的体现，也作为两期点云数据配准及检测变形量的指标；在数据精确配准方面，对关节臂点云数据进行空间格网划分并计算格网点云相似度，在图形强度、相似性测度、及同期点云配准精度的约束下，甄选微变形数据中无变形格网点云为配准基元，实现无标靶微变形点云的高精度配准；在变形量评价方面，应用本文微变形度量方法计算微变形量，并以95％置信度、2倍同期数据配准标准差作为阈值进行变形评价分析，并应用色谱图对变形量进行表达。配准过程如下图所示。主要内容包括格网点云的相似度计算，格网迭代监测点云高精度配准，微变形监测与评价。

其中，ICP算法计算公式如下：

式中：配准模型如公式(1)所示，公式中f(H)为目标函数，反映在变换矩阵H下，两个控制点集的差异程度；Q_i和P_i为待配准点云，R为旋转矩阵，T为平移矩阵，n为参与配准的点对数量，σ₀即为配准精度，表示配准后数据最邻近点对的欧氏距离中误差。

ICP算法步骤如下：

1)在第一期第二次数据P中取点P_i∈P；

2)寻找第一期第一次数据中的对应点Q_i∈Q，要求Q_i与P_i欧氏距离最小；

3)计算旋转矩阵R和平移矩阵T，使得f(H)最小；

4)对P_i使用上一步求取的旋转矩阵R和平移矩阵T进行六参数变换，得到新的对应点集P_i′，用以代替P_i参与下一次迭代计算；

5)计算当前旋转矩阵中的角度变化值

Δω,Δκ；

6)角度变化值

Δω,Δκ均小于0.1秒时，则停止迭代计算，否则回到第一步，直至收敛为止；

7)根据上述公式(1)和(2)计算f(H)、σ₀。

ICP算法精细配准后，基于点云数据构建不规则三角网，并逐点计算同期数据间的变形度量D，D表示第一期第二次数据的三角网顶点到第一期第一次数据的三角网面片的有向距离。

计算步骤如下：

1)以第一期第一次数据T_ij为参考数据，以ICP配准后的第一期第二次数据T_i,j+1为待检测数据；

2)以T_i,j+1的三角网顶点K为出发点，并以该点的法线

作为方向构建射线方程R(t)：

式中，K为待检测数据T_i,j+1的三角网顶点坐标；t为系数，

为K点的单位法向量。

3)假定从K点出发的射线R(t)与参考数据T_ij交于点K′，且K′位于参考数据T_ij某个三角形ΔABC上，那么这个交点就可以用以下两个方程来表示：

式中，K′为从K点出发的射线R(t)与参考数据T_ij的交点坐标；A，B，C分别代表三角形ΔABC的三个顶点坐标；(1-u-v)，u，v代表顶点A，B，C坐标的系数；

联立公式得：

求得参数u,v,t后，根据公式求出K′点的坐标值。

4)计算所述变形度量D：

|D|＝||K-K′||            (7)

式中，变形度量D为待检测三角网

的三角形顶点K到参考三角网

面片的有向距离。当解算的参数t值为正数时，有向距离为负，表示缩进；当解算的参数t值为负数时，有向距离为正，表示凸出。实际上，D的值与参数值-t相同。

计算出各顶点的D值后，根据公式并求出变形度量中误差σ_d；

鉴于D为离散型随机变量，且E(D)＝0则：

式中，σ_d为变形度量的中误差，是判断不同期次数据是否发生变化的指标；D_i为各顶点对应的变形度量值。

1点云格网划分及相似度计算

首先应用点云同名特征计算两期点云的粗配准参数，运用空间相似变换模型把两期点云数据进行坐标统一；然后以点云重心为原点，按照一定的步长对两期监测点云进行空间格网划分；最后计算对应格网点云的相似度，并对其进行排序，以格网点云相似度作为甄选无变形格网点云的基础。每个网格相当于一个控制标靶，理论上通过四个及以上的控制格网点云即可实现监测数据的快速高精度配准，为了保证配准的图像强度，本文以八象限均有控制格网为最优配准方式。空间格网划分的步骤如下：

第一步，确定格网步长d。将T₁₁点云作为参考数据，基于粗配准的旋转平移参数，并以参考数据的重心作为坐标系原点，d_x、d_y、d_z，分别代表在X轴、Y轴和Z轴三个轴向上的步长。将每个网格点云视作一个控制标靶，网格步长可视为实际拼接用的标靶，比如分米立方体的边长。

第二步，计算网格个数。根据网格步长，计算出监测数据T₁₁在三个轴向上的网格数量。先比较T₁₁数据各点在X、Y、Z三个轴向上的最大值和最小值，分别计为X_max,X_min,Y_max,Y_min,Z_max,Z_min，然后根据已知的步长d_x、d_y、d_z，分别求出在X轴、Y轴和Z轴向上格网个数。划分格数为n_x＝X_max/d_x-X_min/d_x+2，同理，n_y＝Y_max/d_y-Y_min/d_y+2，n_z＝Z_max/d_z-Z_min/d_z+2。总的格网数为N_total＝n_x*n_y*n_z。

第三步，统计非零格网。监测数据仅有物体表面的三维几何信息，而格网是由大小相同且排列有序的网格构成的体状空间，所以必定存在不含任何监测数据的空格网，为加速网格迭代检索，需将这些格网予以剔除，计算T₁₁各格网的点数，标记并统计点数为零的格网，计其数量为N₀，最后统计出非零格网的数量N_valid＝N_total-N₀。

划分后各个对应格网点云可以看成是有变形或者无变形的同名格网点云。为了评价同名格网点云是否发生变形或者发生变形的程度，可以通过计算两期数据同名格网点云的相似度来实现。可描述点云形状的函数有D1,D2,D3,D4距离等，D1距离函数是指一个固定点到模型表面上任意随机点之间的距离，D2距离函数是指模型表面上任意两个随机点之间的距离，D3距离函数是指模型表面上任意三个随机点之间的三角形面积的平方根，D4距离函数是指模型表面上任意四个随机点之间构成的四面体的体积的三次方立方根，其中D1形状函数计算方法效率较高，距离函数是一个向量集，其原理如图5-3所示，本发明选择D1形状函数距离直方图为特征，应用该特征的相似度作为衡量两期数据表面几何变形的程度，相似度计算如公式(11)(12)(13)所示。

H_k为参考点云第k个格网内点云D1形状函数距离直方图，H_k(i)表示H_k直方图第i个距离间隔所包含点数，

表示H_k直方图各距离间隔内的点数均值；H_k′为待检测点云第k′个格网内点云D1形状函数距离直方图，H_k′(i)表示H_k直方图第i个距离间隔所包含点数，

表示H_k′直方图各距离间隔内的点数均值。

式中，d(H_k,H_k′)表示对应第k对格网点云的D1形状函数距离直方图相似度，n为直方图横向距离间隔数。

按上述方法，计算两期数据同名格网间直方图的相似度，并对各象限格网的D1距离相似度进行排序。从理论上分析，相似度越低，发生微变形的可能性越大，当然，该相似度与格网的粒度及微变形程度具有一定的相关性。

2格网点云迭代高精度配准

微变形量检测的时间间隔根据残损程度判断而定，为了通过检测多期变化量而进行更有价值的分析和检测，一般的检测时间段内并不会发生全面的大量的变化。各个象限内格网的D1距离相似度代表了格网内点云的变化程度，在检测时间段内，一定粒度格网情况下一定存在不变形的格网数据。根据上述观点，对格网点云进行迭代配准，甄选出不变形格网点云，以其配准结果作为计算和评价微变形量及微变形情况的指标，完成无标靶精细点云的高精度配准及变形量的检测。

为确保甄选的格网点云的配准精度，本文在待配准格网选择过程中遵循以下条件：

1)格网点云个数条件：计算中，点云数量小于格网点云个数均值1/3的格网不参与甄选和配准(点云数量除以有效格网数，即为格网点云个数均值，数量小于该均值1/3的格网由于其点数太少而不参与甄选)；

2)配准图像强度条件：要求空间每个象限至少有一个格网作为待配准格网；

3)迭代配准截止条件：以同期两次无变形点云配准精度为非同期微变形点云配准的阈值。

格网迭代高精度配准的步骤如下：

(1)对同期扫描数据T_i，1和T_i，2进行ICP配准，配准模型如公式(1)所示，配准精度为σ₀；

(2)针对T_i，1和T_i+1,1每个象限格网D1距离直方图相似度排序结果，选出相似度最大的前n个格网数据作为配准基元；

(3)对选出的8*n个格网点云数据进行ICP配准，如果配准精度没有达到σ₀阈值要求，去掉8*n个格网点云中相似度最低的格网数据，再进行ICP配准；

(4)按照(3)步的约束条件进行格网迭代配准，当某个象限中的格网数目变为零时仍然没有达到σ₀的配准精度，则缩小格网划分的粒度，再次进行上述(2)-(3)，直到迭代满足条件为止。

n的经验值为3，理论上n为1应该能够减少迭代的次数，但为了避免误差而导致的甄选格网错误造成的数据退化，本文以3为经验值。格网的粒度受对象的体量、残损微变形的程度及分布等因素的影响，初始值可根据经验而定。本文在实验部分中将对格网粒度的敏感度进行实验和分析。

3微变形检测与评价

3.1微变形检测

相比于点到点的配准特征度量，采用点到面的微变形度量更能反映物体表面真实的变化情况。本文将基于高精度关节臂点云数据构建不规则三角网，然后以前一期的三角网顶点到后一期三角网面片的距离作为变形度量。

具体算法如下：

1)以第一期第一次数据T_ij为参考数据，以ICP配准后的第二期第一次数据T_i+1,j为待检测数据；

2)以T_i+1,j的三角网顶点K为出发点，并以该点的法线

作为方向构建射线方程R(t)：

式中，K为待检测数据T_i+1,j的三角网顶点坐标；t为系数，

为K点的单位法向量。

3)假定从K点出发的射线R(t)与参考数据T_ij交于点K′，且K′位于参考数据T_ijT_ij某个三角形ΔABC上，那么这个交点就可以用以下两个方程来表示：

式中，K′为从K点出发的射线R(t)与参考数据T_ij的交点坐标；A，B，C分别代表三角形ΔABC的三个顶点坐标；(1-u-v)，u，v代表顶点A，B，C坐标的系数；

联立公式得：

求得参数u,v,t后，根据公式求出K′点的坐标值。

4)计算所述变形度量D：

|D|＝||K-K′||               (19)

式中，变形度量D为待检测三角网T_i+1,j的三角形顶点K到参考三角网T_ij面片的有向距离。当解算的参数t值为正数时，有向距离为负，表示缩进；当解算的参数t值为负数时，有向距离为正，表示凸出。实际上，D的值与参数值-t相同。

3.2微变形评价

引起不同期次检测数据比较后发生变化的原因有两种，第一是整体系统误差所致；第二是对象表面确实发生几何变化所致。变化评价的主要内容就是要建立一种微变形评价方法，用于评价对象表面是否确实发生变化以及变化量是多少。为此，本文将依据变化阈值来构建微变形评价方法。变化阈值是指对象表面实际发生变化的临界值，超出该临界值可认为相应区域已发生变化，而低于该临界值说明该区域未发生变化。

事实上，激光扫描数据的精度评价是一个非常复杂的问题，它取决于仪器设备，地理参考、对象材质、激光入射角度等。而利用激光扫描数据来检测对象表面几何变化，其监测精度评价过程就更加复杂，它不仅要考虑上述单次数据采集的影响因子，还要考虑不同期次数据的配准误差，很难通过构建严密的误差方程并按误差传播律求得监测精度值。

因此，本文将抛弃传统的误差评定方法，而是根据统计学原理，先对同期邻次监测数据ICP配准，然后逐点计算同期数据间的变形度量值D，最后统计D值概率密度分布情况并求出其中误差σ_d。由于同期数据T_Cj代表同一时段的数据，其表面实际未发生几何变化，所以基于本文方法对监测数据进行比较所产生的变形值是由于监测误差所致，其变形度值中误差σ_d即为本文微变形监测变形度量的精度指标。在变化评价上，采用置信概率为95％的2σ_d作为阈值，即变化超出置信区间(-2σ_d～+2D_μ)的值视为变化值。当然在实际应用，如果希望比对结果更加保守，亦可采用置信概率为99.7％的3σ_d作为变化评价阈值。

4实验与结论

为验证本文微变形检测方法，实验共采集3组试验数据，分别是X₁₁,X₁₂,X₂₁，其中，X₁₁,X₁₂为同期数据，通过X₁₁,X₁₂两次同期数据的比较，得出对象表面未发生变化的情况下关节臂扫描点云的配准精度σ₀，以此作为邻期格网迭代法的配准精度收敛条件。T₂₁为后一期监测数据，两期数据时间间隔为六个月。通过格网迭代法配准T₁₁,T₂₁两次邻期数据的配准，得出造像表面的实际变化情况。

试验数据均为基于关节臂扫描的精细点云数据及其三角网模型，其顶点点位精度优于0.1mm，高出相位式地面激光扫描仪2个数量级；三角网的平均边长优于0.3mm，高出相位式地面激光扫描仪1个数据级。

选择中国某石窟寺实际造像为实验对象，该造像材质为砂岩，受风化影响，表面已严重受损，整体外观面目全非，如图3所示，据石窟寺技术管理人员介绍，造像风化剥落现象还在持续进行，近些年有愈演愈烈的迹象。

获取实验对象表面高精密点云数据，要求各组实验数据间的范围、三角网边距基本一致。为了获取仪器误差及系统配准误差，应用关节臂数据进行同期次两次数据和非同期两次数据的采集，设置扫描对象点云分辨率0.2-0.3mm，单次扫描点云个数约800万。图4分别展示了同期两次数据T₁₁T₁₂及时间间隔6个月的非同期数据T₂₁。

两次扫描的各种条件及参数设置保持一致，应用ICP算法对两次数据T₁₁T₁₂进行配准，并计算其标准中误差D_σ(T₁₁T₁₂)，，配准结果如图5所示，从配准精度表1可以看出，同期数据的配准精度为0.1897mm，其配准误差主要体现在仪器误差等相关的系统误差。该成果用来作为非同期有变形数据高精度配准的标准。变形度量中误差σ_d为0.1973mm。

表1 同期数据ICP配准结果

Table 1 Registration results of data ICP in the same period

4.1微变形数据空间格网迭代高精度配准

通过上述人为变形数据的验证，对间隔六个月的实际风化微变形量数据T₁₁和T₂₁应用本文的方法进行格网迭代高精度配准。首次格网划分为5cm，划分结果如图6所示。

针对上述划分数据，应用粗配准结果计算其D1距离及直方图相似度，分别选出8个象限中D1距离直方图相似度最大的前3个共24个格网数据，结果如图7所示。

对上述数据进行顾及空间图形强度的ICP迭代计算，每次去除相似度最差格网点云，其迭代过程及配准精度如图8及表3所示。

表3 5cm格网迭代及精度变化

Table 3 5cm Grid iteration and the change of precision

从图8及表3结果中可以看出，5cm格网划分情况下，迭代到12次时，第8象限格网数变为0，而此时的配准精度为0.232458mm，还没有满足同期数据0.1897mm的配准精度，所以需要再次缩小格网粒度，以2.5cm为基准划分格网。针对2.5cm粒度划分点云数据，应用粗配准结果计算其D1距离及直方图相似度，分别选出8个象限中D1距离直方图相似度最大的前3个格网数据，结果如图9所示。

对上述数据进行顾及空间图形强度的ICP迭代计算，每次去除相似度最差格网点云，其迭代情况及配准精度如图10及表4所示

表4 2.5cm格网迭代及精度变化

Table 4 2.5cm Grid iteration and the change of precision

从图10及表4的结果中可以看出，2.5cm格网划分情况下，迭代到7次时，配准精度为0.178258mm，满足同期数据0.1897mm的配准精度，配准结束。

对两期非数据进行常规的ICP配准，其配准精度与本文方法的配准精度进行对比，如表5所示，从表中可以看出传统ICP算法配准精度不能达到要求，体现了配准精度方面的优越性。

表5 配准方法比较

Table 5 Comparison of registration methods

4.2实际微变形数据变形量距离计算及变形评价

对满足配准精度的微变形数据进行变形距离的计算。以T₁₁为基准，计算T₂₁到T₁₁的有向距离。循环遍历T₂₁每一个点，应用本文微变形度量计算方法计算其到T₁₁的距离，以同期两次数据T₁₁、T₁₂间变形度量值中误差σ_d为基准，变化阈值设为2σ_d。同期数据配准精度为0.1973mm，所以能检测出的最小变化量为0.3946mm，按照此阈值，对实际间隔6个月的两期数据进行微变形量统计和表达，变化范围为0.4-2mm，变化步长为0.2mm，T₁₁与T₂₁的变化色谱图如图11所示，此色谱图表达了实际微变形数据应用本文方法所能检测的风化微变形量的结果。图12为基座表面微变形情况，检测结果显示基座之上多为蓝色区域，代表该区域出现了较大范围的缩进现象，也就是说，基座表面的积土变薄了，这显然不符合常理，但经过与管理人员确认后方知，在检测期间，工作人员对基座进行了清扫，使得原有的积土变薄，并且基座变薄处还留有脚印，这恰恰体现了本文方法的正确性。

5结论

本章节以石窟寺造像风化微变形实际同期及非同期关节臂扫描点云数据为研究对象，应用本文提出的方法进行高精度配准及变形评价。为文物表面微变形量高精度检测提供了新的研究方法。针对于两期微变形数据的无标靶配准，应用同期两次扫描数据配准精度作为两期数据配准时的阈值，应用格网划分甄选无变形格网为高精度配准策略，并顾及图形强度进行格网粒度的迭代，以同期次配准精度为约束进行ICP迭代，甄选微变形数据中无变形量进行高精度配准。本文应用有风化造像实际间隔6个月的扫描数据，以本文的方法进行格网迭代，甄选出非变形区域点云，通过ICP迭代配准达到配准精度的要求，按照2倍同期数据变形度量值中误差为阈值，检测出了0.4mm以上的变化量，对变化区域及变化量给出了正确的科学的判断和表达。本文微变形量检测方法打破了现有ICP迭代配准的策略，以甄选微变形数据中无变形量区域为思想导向，以标靶配准思想用8象限数据非空顾及空间配准图形强度，以同期数据配准精度微变化标准进行高精度配准，为微变形检测提供了创新性的思路，并用实际数据证明了方法的正确性及其在文物表面微变形量检测中的极大优势，具有重要的理论和现实意义。本文的研究只是初步探索，后续会在配准空间图形强度的自适应迭代、多种格网相似性度量等多方面展开更加细致的研究。

尽管本发明的实施方案已公开如上，但其并不仅仅限于说明书和实施方式中所列运用，它完全可以被适用于各种适合本发明的领域，对于熟悉本领域的人员而言，可容易地实现另外的修改，因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节和这里示出与描述的图例。

Claims (8)

1.一种基于无控制扫描点云的文物对象高精度微变形监测方法，其特征在于，包括如下步骤：

步骤一、在一个时期内对文物对象进行连续两次数据采集，获得两个数据：第一期第一次数据和第一期第二次数据，应用ICP算法对所述第一期第一次数据和第一期第二次数据进行配准，得出文物对象表面未发生变化的情况下关节臂扫描点云的配准精度σ₀，ICP算法精细配准后，基于点云数据构建不规则三角网，并逐点计算同期数据间的变形度量D，D表示第一期第二次数据的三角网顶点到第一期第一次数据的三角网面片的有向距离，并求出变形度量的中误差σ_d，所述ICP算法为Iterative Closest Point算法；

步骤二、在另一个时期内对所述文物对象再进行数据采集，获得第二期第一次数据，将所述第二期第一次数据和所述第一期第一次数据进行格网迭代法配准，使其配准精度小于或等于所述配准精度σ₀，其中，所述第一期第一次数据、第二次数据和所述第二期第一次数据均为基于关节臂扫描的精细点云数据；

步骤三、对满足配准精度的第二期第一次数据进行变形度量的计算：计算变形度量时，基于点云数据构建不规则三角网，以第二期第一次数据的三角网顶点到第一期第一次数据的三角网面片的距离作为不同期数据间的变形度量D，将所述不同期数据间的变形度量D值与所述同期数据间的变形度量D值相比，若其变化超出置信区间(-mσ_d～+mσ_d)的值视为所述文物对象的变化值，其中，m为不大于3的正整数。

2.如权利要求1所述的基于无控制扫描点云的文物对象高精度微变形监测方法，其特征在于，步骤二中，所述配准包括如下步骤：

2.1)首先应用点云同名特征计算所述第一期第一次数据和第二期第一次数据点云的粗配准参数，并运用空间相似变换模型把两期点云数据进行坐标统一；

2.2)以点云重心为原点，按照一定的步长对上述两期点云数据进行空间格网划分，并计算同名格网点云相似度，以同名格网点云相似度作为甄选无变形格网点云的基础；

其中，以D1形状函数距离直方图这一特征的相似度衡量所述同名格网点云相似度，D1为格网内各点到格网点云重心距离，格网点云相似度计算公式如(11)(12)和(13)所示：

H_k为参考点云第k个格网内点云D1形状函数距离直方图，H_k(i)表示H_k直方图第i个距离间隔所包含点数，表示H_k直方图各距离间隔内的点数均值；H_k′为待检测点云第k′个格网内点云D1形状函数距离直方图，H_k′(i)表示H_k′直方图第i个距离间隔所包含点数，表示H_k′直方图各距离间隔内的点数均值；

式中，d(H_k，H_k′)表示对应第k对格网点云的D1形状函数距离直方图相似度，n为直方图横向距离间隔数；

2.3)对上述两期的同名格网点云进行迭代配准，甄选出无变形格网点云，之后对上述两期点云数据进行高精度配准，具体步骤为：

2.3.1)针对所述第一期第一次数据和第二期第一数据每个象限内格网D1距离直方图相似度排序结果，选出相似性大的前n个格网数据作为配准基元，其中，n≥3；

2.3.2)对选出的8*n个格网点云数据进行ICP配准，如果配准精度没有达到σ₀阈值要求，去掉8*n个格网点云中相似度最低的格网数据，再进行所述ICP配准；

2.3.3)按照2.3.2)步骤的约束条件进行格网迭代配准，当某个象限中的格网数目变为零时仍然没有达到σ₀的配准精度，则缩小格网划分的粒度，再次进行2.3.1))-2.3.2)，直到迭代满足条件为止。

3.如权利要求1或2任一所述的基于无控制扫描点云的文物对象高精度微变形监测方法，其特征在于，ICP算法包括如下计算公式：

式中：式(1)所示为配准模型，公式中f(H)为目标函数，反映在变换矩阵H下，两个控制点集的差异程度；Q_i和P_i为待配准点云，R为旋转矩阵，T为平移矩阵，n为参与配准的点对数量，σ₀即为配准精度，表示配准后数据最邻近点对的欧氏距离中误差；在步骤一中，Q_i为第一期第一次点云数据，P_i为第一期第二次点云数据；ICP算法步骤如下：

1)在第一期第二次数据P中取点P_i∈P；

2)寻找第一期第一次数据中的对应点Q_i∈Q，要求Q_i与P_i欧氏距离最小；

3)计算旋转矩阵R和平移矩阵T，使得f(H)最小；

4)对P_i使用上一步求取的旋转矩阵R和平移矩阵T进行六参数变换，得到新的对应点集P′_i，用以代替P_i参与下一次迭代计算；

5)计算当前旋转矩阵中的角度变化值Δω，Δκ；

6)角度变化值Δω，Δκ均小于0.1秒时，则停止迭代计算，否则回到第一步，直至收敛为止；

7)根据上述公式(1)和(2)计算f(H)、σ₀。

4.如权利要求1所述的基于无控制扫描点云的文物对象高精度微变形监测方法，其特征在于，所述同期数据间的变形度量D和变形度量的中误差σ_d的计算步骤如下：

1)以第一期第一次数据T_ij为参考数据，以ICP配准后的第一期第二次数据T_i，j+1为待检测数据；

2)以T_i，j+1的三角网顶点K为出发点，并以该点的法线作为方向构建射线方程R(t)：

式中，K为待检测数据T_i，j+1的三角网顶点坐标；t为系数，为K点的单位法向量；

3)假定从K点出发的射线R(t)与参考数据T_ij交于点K′，且K′位于参考数据T_ij某个三角形ΔABC上，那么这个交点就可以用以下两个方程来表示：

式中，K为从K点出发的射线R(t)与参考数据T_ij的交点坐标；A，B，C分别代表三角形ΔABC的三个顶点坐标；(1-u-v)，u，v代表顶点A，B，C坐标的系数；

联立公式得：

求得参数u，v，r后，根据公式求出K′点的坐标值；

4)计算所述变形度量D：

|D|＝||K-K′||     (7)

式中，变形度量D为待检测三角网T_i，j+1的三角形顶点K到参考三角网T_i，j面片的有向距离；当解算的参数t值为正数时，有向距离为负，表示缩进；当解算的参数t值为负数时，有向距离为正，表示凸出；实际上，D的值与参数值-t相同；

计算出各顶点的D值后，根据公式并求出变形度量中误差σ_d；

鉴于D为离散型随机变量，且E(D)＝0则：

式中，σ_d为变形度量的中误差，是判断不同期次数据是否发生变化的指标；D_i为各顶点对应的变形度量值。

5.如权利要求1所述的基于无控制扫描点云的文物对象高精度微变形监测方法，其特征在于，步骤三中，计算变形度量值时，基于点云数据构建不规则三角网，以第二期第一次数据的三角网顶点到第一期第一次数据的三角网面片的有向距离作为变形度量，具体方法包括如下步骤：

4.1)以第一期第一次数据T_ij为参考数据，以ICP配准后的第二期第一次数据T_i+1，j为待检测数据；

4.2)以T_i+1，j的三角网顶点K为出发点，并以该点的法线作为方向构建射线方程R(t)：

式中，K为待检测数据T_i+1，j的三角网顶点坐标；t为系数，为K点的单位法向量；

4.3)假定从K点出发的射线R(t)与参考数据T_ij交于点K′，且K′位于参考数据T_ij某个三角形ΔABC上，那么这个交点就可以用以下两个方程来表示：

式中，K′为从K点出发的射线R(t)与参考数据T_ij的交点坐标；A，B，C分别代表三角形ΔABC的三个顶点坐标；(1-u-v)，u，v代表顶点A，B，C坐标的系数；

联立公式得：

求得参数u，v，t后，根据公式求出K′点的坐标值；

4.4)计算所述变形度量D：

|D|＝||K-K′||     (19)

式中，变形度量D为待检测三角网T_i+1，j的三角形顶点K到参考三角网T_ij面片的有向距离；当解算的参数t值为正数时，有向距离为负，表示缩进；当解算的参数t值为负数时，有向距离为正，表示凸出；实际上，D的值与参数值-t相同。

6.如权利要求1所述的基于无控制扫描点云的文物对象高精度微变形监测方法，其特征在于，m＝2。

7.如权利要求1所述的基于无控制扫描点云的文物对象高精度微变形监测方法，其特征在于，应用关节臂三维激光扫描系统进行数据的采集，设置扫描对象点云分辨率0.2-0.3mm，单次扫描点云个数约800万。

8.如权利要求1所述的基于无控制扫描点云的文物对象高精度微变形监测方法，其特征在于，所述一个时期为3个日历天。

Images