CN110838115B - 运用轮廓线提取和四维曲面拟合的古文物三维模型变化检测方法 - Google Patents

Abstract

本发明公开了一种运用轮廓线提取和四维曲面拟合的古文物三维模型变化检测方法，包括以下步骤：对古文物两个不同时期的LiDAR点云数据集合进行配准；对LiDAR点云数据进行重分类，生成LiDAR点云模型；根据基于高斯球平行投影的轮廓点快速计算算法获取LiDAR点云模型的轮廓点；依次使用最小生成树算法、短边剪枝方法、光滑技术获得光滑的点云轮廓线；对提取到的点云轮廓线进行离散程度筛选，去除点云轮廓线中的噪声特征，得到两个时期的轮廓特征图；对筛选得到的轮廓特征图进行特征对比，得到对比结果；对得到的对比结果进行筛选，得到差别特征。本发明充分利用了LiDAR点云数据的优越性，具有周期短，成本低，精度高，更新时间快等特点，能够保证最终结果的精度。

Description

运用轮廓线提取和四维曲面拟合的古文物三维模型变化检测 方法

技术领域

本发明涉及测绘科学与技术领域，尤其涉及一种运用轮廓线提取和四维曲面拟合的古文物三维模型变化检测方法。

背景技术

变化检测的数据源有正射影像、遥感影像、LiDAR点云数据等。其中，LiDAR点云数据的获取不受季节天气光照等环境因素影响，同时获取数据的周期短，成本低，精度高，数据更新快，是进行古文物表面三维变换检测的良好数据源。但缺乏纹理信息是LiDAR点云数据的一个严重的缺点。因此通过基于高斯球平行投影的轮廓点快速计算等一系列算法获取物体表面的轮廓特征有利于对古文物进行高精度的三维变换检测。

现有的影像变化检测方法主要分为三个层次：(1)像素级变换检测；(2)特征级变换检测；(3)目标级变换检测。其中，特征级变换检测方法通常都是采用一定的算法先从原始图像中提取特征信息，然后对这些特征信息进行综合分析与对比，获取检测对象的变换信息。特征级的变化检测方法通过特征进行关联处理，将特征分成了有意义的组合，从而提高了特征属性判断的可信度和准确性。但是由于特征级变换检测过程相对复杂，目前还没有实用的成果出现。本发明使用微移动的方法进行特征的变换检测，巧妙的选取了每次移动的格网间隔，加快了计算机处理的速度。同时采用基于最小二乘法的多项式四维曲面拟合方法以及改进的随机游走算法，避免了局部最优解的出现，提高了对比检测结果的精度。

目前，市场上有很多使用二维影像的平面信息进行变换检测的方法，而使用物体表面的三维信息进行变换检测的实用算法较少。对于古文物来说，三维表面的变化情况是其变化的重要部分。因此该检测方法的研究对于进行古文物高精度，高效率的三维变化检测以及LiDAR点云数据轮廓线提取的处理和应用都有着十分重要的研究价值。

发明内容

本发明要解决的技术问题在于针对现有技术中基于特征的变换检测方法的缺陷，提供一种运用轮廓线提取和四维曲面拟合的古文物三维模型变化检测方法。

本发明解决其技术问题所采用的技术方案是：

本发明提供一种运用轮廓线提取和四维曲面拟合的古文物三维模型变化检测方法，该方法包括以下步骤：

步骤1、获取古文物在两个不同时期的LiDAR点云数据集合，对两个不同时期的LiDAR点云数据集合进行配准；

步骤2、对LiDAR点云数据进行重分类，生成LiDAR点云模型；

步骤3、根据基于高斯球平行投影的轮廓点快速计算算法获取LiDAR点云模型的轮廓点；

步骤4、根据获取的轮廓点，依次使用最小生成树算法、短边剪枝方法、光滑技术获得光滑的点云轮廓线；

步骤5、对提取到的点云轮廓线进行离散程度筛选，去除点云轮廓线中的噪声特征，得到两个时期的轮廓特征图；

步骤6、对筛选得到的轮廓特征图进行特征对比，得到对比结果；

步骤7、对得到的对比结果进行筛选，得到差别特征。

进一步地，本发明的步骤1中对两个不同时期的LiDAR点云数据集合进行配准的方法为：通过平移、旋转操作使两个LiDAR点云数据集合的位置重合。

进一步地，本发明的步骤2中进行重分类生成LiDAR点云模型的方法为：

记LiDAR点云数据的个数为M个，点云数据集合的长，宽，高之比为L：W：H；将点云数据集合按照

的格网进行划分；统计每个格网中点云数据的个数，记为m；若m≥1，则记该格网中心点为有效点云数据，否则该格网中无有效点云数据。

进一步地，本发明的步骤3中获取LiDAR点云模型轮廓点的方法为：

步骤31、通过分析LiDAR点云模型上各点的K近邻的协方差矩阵来计算该点的法向量；记LiDAR点云模型上某一点为P，找出点P邻域点的集合，记为N_p，计算出N_p的重心，以此建立起一个协方差矩阵C：

其中，

协方差矩阵C的最小特征值对应的特征向量就是该点的法向量；

步骤32、构建高斯球，利用四叉树的方法建立分层数据结构进行数据存储；将LiDAR点云模型中的点的法向映射到单位矢量球上；

设N＝(N_x，N_y，N_z)是点云模型上点P的正则化法向量，通过以下公式获得点P在高斯球上的映射位置：

建立以高斯球球心为坐标系原点的球坐标系，用经度1d和纬度ud来表示法向的球坐标；通过转换公式将高斯球坐标转换为以高斯球球心为坐标原点的直角坐标，转换公式如下：

x＝cos(ud)cos(ld),y＝cos(ud)sin(ld),z＝sin(ld)

利用四叉树模型建立分层数据结构组织点云数据；

步骤33、获取轮廓点：首先求出当前视向量V；随后从四叉树的根节点开始检测，如果节点的四个角点的法向与投影向量点积同号，说明它们的可见性状态一致，该节点内无轮廓点；丢弃该节点及其所包含的子节点，返回父节点；否则继续检测该节点的子节点，做同样的判断，直到检测到叶节点，记录该节点内点云数据轮廓点的索引；同样处理其它子节点，直到最后一个子节点；根据记录的节点内的索引计算出三维轮廓点集，记为Q。

进一步地，本发明的步骤32中利用四叉树模型建立分层数据结构组织点云数据的方法为：

对高斯球沿经纬线进行层次划分，得到四边区域，称为单元；利用四叉树结构存储这个层次划分，四叉树的根节点对应包含整个球的单元；对单元逐层划分，如果单元内的点小于一个设定的阈值或者单元的四个角点经纬度差别小于一个设定的阈值，则不再划分，这些不再做划分的单元就是四叉树的叶节点；对于四叉树的中间节点，只记录其四个儿子节点及其四个角点的信息，以及投射到此节点对应单元的模型上点的个数，而对于四叉树的叶节点，除了上述信息，还记录所有投射到此叶节点内的对应的模型上的点索引序号以及点的个数。

进一步地，本发明的步骤4中获取点云轮廓线的方法为：

步骤41、根据轮廓点的位置信息创建图G1；

在三维轮廓点集Q中确定一个起始点q₀；计算点集Q中的点与q₀的距离d，如果存在点q_i，使d<阈值d_max，d_max的值需根据三维模型点的相对位置情况进行确定，则连接q_i与q₀；再将q_i作为起始点重复上述步骤；当一次连接中所有的距离都不小于阈值D时，连接结束；再确定一个起始点，重复进行连接步骤，直到所有的点都遍历完毕；

步骤42、根据Kruskal最小生成树算法创建图G2；

把图G1中的顶点看成独立的树，以顶点连接的线的条数作为权值，按权值从小到大选择边，所选的边连接的两个顶点属于两颗不同的树，则成为最小生成树的一条边，并将这两颗树合并作为一颗树；重复上述步骤，直到所有顶点都在一颗树内或者有n-1条边为止；

步骤43、通过剪枝构造创建图G3；

检测图G2中所有点的分支数，记为m；若m>2,则检测各分支所占的格网个数；若至少有两个分支所占的格网个数≥阈值L，阈值L根据先验知识确定，则将所有小于阈值L的分支除去；若只有一个分支所占的格网个数≥L并且至少有两个分支所占的格网个数<L，则保留格网个数≥L的分支并在剩下的分支中选择最长的一条；若m≤2，则保留这些分支，不做任何处理；

步骤44、采用平面曲线演化法光滑处理轮廓线；

将图G3的点进行投影，得到由折线表示的初始轮廓线，计算折线顶点处的两折线的角度大小及角平分线；根据折线顶点处的角度移动该顶点的位置，顶点的角度越小，该处的局部曲率越大，顶点移动的距离越大；重复上述移动过程，直到折线上每个顶点的角度大于设定阀值，表示曲线达到了要求的光滑程度；上述过程的公式为：

其中，P_i ⁿ是点P_i第n次移动的结果；k_i为P_i处的角度大小，表示曲率：N_i为点P_i角平分线方向，代表法向。

进一步地，本发明的步骤5中进行离散程度筛选的方法为：

对提取到的轮廓线进行多项式拟合，多项式最高次项项数选取方法如下：

设特征点的集合为P_m，其对应的坐标集合为V(x，y，z)；分别计算x,y,z坐标的方差D(x)，D(y)和D(z)；根据先验知识以及古文物类型确定六个阈值，x坐标的阈值记为xmin，xmax，xmin<xmax，y坐标的阈值记为ymin,ymax，ymin<ymax，z坐标的阈值记为zmin，zmax，zmin<zmax；通过比较D(x)与xmin，xmax的大小，D(y)与ymin，ymax的大小以及d(z)与zmin，zmax的大小来确定多项式最高次项项数；

若D(z)<zmin，当D(x)<xmin时，x的最高次项为一次项；当xmin<D(x)<xmax，x的最高次项为二次项；当D(x)>xmax时，x的最高次项为三次项，y与x同理；

若zmin<D(z)<zmax，当D(x)<xmin时，x的最高次项为二次项；当xmin<D(x)<xmax，x的最高次项为三次项；当D(x)>xmax时，x的最高次项为四次项，y与x同理；

若D(z)>zmax，当D(x)<xmin时，x的最高次项为三次项；当xmin<D(x)<xmax，x的最高次项为四次项；当D(x)>xmax时，x的最高次项为五次项，y与x同理；

计算轮廓线与拟合的多项式的差值，若差值小于相应的阈值，则判断该特征为噪声，否则保留经过多项式拟合后的轮廓，相应阈值根据古文物的类型确定。

进一步地，本发明的步骤6中进行特征对比的方法为：

分别记两张正射影像的轮廓特征图为C1时期特征图和C2时期特征图；

采用分块对比方法，将点云数据分为4*4*4的模块，对应模块之间进行对比；

对比时采用基于特征的对比方法；先以C1时期的特征图为底图，C2时期的特征图为对比图进行对比；选取C2时期的特征图中的特征M，在C1时期的特征图中找到与特征M相对应的坐标区域，在该区域的Δ范围内进行微移动；移动方法如下：

设C2时期的特征图中特征M的边界坐标集合为K(X，Y，Z)；记C1时期的特征图中与K(X，Y，Z)相对应的坐标集合为K(X1，Y1，Z1)，记坐标集合K(X，Y，Z)中X坐标的最大值为maxX，Y坐标的最大值为maxY，Z坐标的最大值为maxZ；记坐标集合K(X1，Y1，Z1)中X1坐标的最大值为MaxX1，最小值为MinX1，Y1坐标的最大值为MaxY1，最小值为MinY1，Z1坐标的最大值为MaxZ1，最小值为MinZ1；则特征M在C1时期特征图中可进行微移动的立方体范围的角点坐标为：

(MaxX1+ΔX，MaxY1+ΔY，MaxZ1+ΔZ)，(MaxX1+ΔX，MinY1-ΔY，MaxZ1+ΔZ)，(MinX1-ΔX，MaxY1+ΔY，MaxZ1+ΔZ)，(MinX1-ΔX，MinY1-ΔY，MaxZ1+ΔZ)，(MaxX1+ΔX，MaxY1+ΔY，MinZ1-ΔZ)，(MaxX1+ΔX，MinY1-ΔY，MinZ1-ΔZ)，(MinX1-ΔX，MaxY1+ΔY，MinZ1-ΔZ)，(MinX1-ΔX，MinY1-ΔY，MinZ1-ΔZ)

分别记为A点，B点，C点，D点，E点，F点，G点，H点；

根据先验知识以及古文物的类别设定阈值ε，当特征M在两特征图中未重合的格网个数小于ε时，则认为特征M在两特征图中完全重合。

进一步地，本发明的步骤6中进行特征对比的具体步骤为：

步骤61、将C2时期特征图上的(maxX，maxY，maxZ)点与C1时期特征图上的(MaxX1+ΔX，MaxY1+ΔY，MaxZ1+ΔZ)点对齐，此时特征M位于ABCDEFGH点所构成的立方体中；

步骤62、在X方向上，特征M每隔

个格网移动一次，在Y方向上，特征M每隔

个格网移动一次，在Z方向上，特征M每隔

个格网移动一次；

每次移动后计算特征M与C1时期特征图中对应特征相差的格网个数，记为M_i；移动方向如下所示，先在x方向移动，当x方向无法继续移动时，在y方向移动一次，继续在x方向移动；当特征M在xy平面移动完毕后在z方向移动一次；

步骤63、基于最小二乘法进行四维曲面的多项式拟合，以坐标集合K(X1，Y1，Z1)中的X1坐标值为x坐标值，Y1坐标值为y坐标值，Z1坐标值为z坐标值，M_i值为t坐标值构造四维曲面T(x，y，z)；

步骤64、根据改进的随机游走算法求出曲面T中t坐标值的最小值，记为MinT；比较MinT与ε值的大小，若MinT＜ε，则特征M在两特征图中完全重合，否则特征M发生了变化；

步骤65、重复步骤61-步骤64，直到C2时期特征图中的所有特征都判断完毕，记作对比结果N1；

以C2时期的特征图为底图，C1时期的特征图为对比图按照步骤61-步骤65进行一次对比，得到对比结果N2。

进一步地，本发明的步骤64中改进的随机游走算法步骤为：

①给定初始迭代点(x₀，y₀，z₀)，初次行走步长λ，控制精度ε1；

②给定迭代控制次数I，K为当前迭代次数，置K＝1；

③当K＜I时，随机生成n个(-1，1)之间的三维向量u_i，u_i＝(x_i1′，y_i1′，z_i1′)(-1＜x_i1′＜1，-1＜y_i1′＜1，-1＜z_i1′＜1)；并将其标准化得到

令(x_i1，y_i1，z_i1)＝(x₀，y₀，z₀)+λu_i，取min{(x₁₁，y₁₁，z₁₁)，(x₂₁，y₂₁，z₂₁)...(x_n1，y_n1，z_n1)}，记为(x_m1，y_m1，z_m1)；完成第一步游走；

④计算函数值，如果T(x_m1，y_m1，z_m1)＜T(x₀，y₀，z₀)，即找到了一个比初始值好的点，那么K重新置为1，将(x_m1，y_m1，z_m1)作为初始值，回到第②步；否则K＝K+1，回到第③步；

⑤如果连续I次都找不到更优的值，则认为，最优解就在以当前最优解为中心，当前步长为半径的I维球内；此时，如果λ＜ε1，则结束算法；否则，令

回到第①步，开始新一轮游走。

本发明产生的有益效果是：本发明的运用轮廓线提取和四维曲面拟合的古文物三维模型变化检测方法，充分利用了LiDAR点云数据检测古文物表面的三维变化情况，数据更新快，成本低，精度高，获取的周期短，保证了最终获得的结果的精度。同时采用基于高斯球平行投影的轮廓点快速计算等一系列算法克服了LiDAR点云数据缺乏纹理信息的不足。进行轮廓对比时，选取了适当的间隔进行特征对比，使用基于最小二乘法的多项式四维曲面拟合方法和改进的随机游走算法，避免了局部最优解的出现，并缩短了计算机的处理时间。为当前文物保护领域提供了一种高精度，高效率的三维变化检测方法。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1是本发明的运用轮廓线提取和四维曲面拟合的古文物三维模型变化检测方法流程示意图。

图2是本发明的具体实施方案中进行轮廓点的筛选时对高斯球坐标和直角坐标进行转换的示意图。

图3是本发明的具体实施方案中特征微移动方向的示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

随着计算机技术的飞速发展，计算机在越来越多的行业有了越来越深度的利用。利用计算机技术有效地保护和研究古文物遗产，对推动重大文物抢救保护和考古工程具有重大意义。

本发明实施例的运用轮廓线提取和四维曲面拟合的古文物三维模型变化检测方法。输入不同时期的古文物LiDAR点云数据，通过计算机软件对该数据进行基于特征的变换检测。

具体实现方法包括以下步骤：

步骤1：对两个不同时期的LiDAR点云数据集合进行配准。通过平移，旋转等操作使两个LiDAR点云集合的位置大致重合。

步骤2：对LiDAR点云数据进行重分类，生成LiDAR点云模型。

记LiDAR点云数据的个数为M个，点云数据集合的长，宽，高之比为L：W：H。将点云数据集合按照

的格网进行划分。统计每个格网中点云数据的个数，记为m。若m≥1，则记该格网中心点为有效点云数据，否则该格网中无有效点云数据。

步骤3：根据基于高斯球平行投影的轮廓点快速计算算法获取LiDAR点云模型的轮廓点。

基于高斯球平行投影的轮廓点快速计算算法主要步骤如下：

(1)通过分析点云模型上点的K近邻的协方差矩阵来计算该点的法向量。

记点云模型上某一点为P。找出点P邻域点的集合，记为N_p，计算出N_p的重心，以此建立起一个协方差矩阵C：

其中，

协方差矩阵C的最小特征值对应的特征向量就是该点的法向量。

(2)构建高斯球，利用四叉树的方法建立分层数据结构以进行数据存储。

首先将点云模型中的点的法向映射到单位矢量球上。

设N＝(N_x，N_y，N_z)是点云模型上点P的正则化法向量，通过公式：

获得点P在高斯球上的映射位置。

由于点云模型上的每个点的法向都会和高斯球上的一个点相对应。因此可以建立以高斯球球心为坐标系原点的球坐标系。此时，用经度1d和纬度ud来表示法向的球坐标。

通过转换公式将高斯球坐标转换为以高斯球球心为坐标原点的直角坐标。转换公式如下：

x＝cos(ud)cos(ld)

y＝cos(ud)sin(ld)

z＝sin(ld)

利用四叉树模型建立分层数据结构组织点云数据。

对高斯球沿经纬线进行层次划分，得到四边区域，称为单元。利用四叉树结构存储这个层次划分，四叉树的根节点对应包含整个球的单元。对单元逐层划分，如果单元内的点小于一个设定的阈值或者单元的四个角点经纬度差别小于一个设定的阈值，则不再划分，这些不再做划分的单元就是四叉树的叶节点。对于四叉树的中间节点，只记录其四个儿子节点及其四个角点的信息，以及投射到此节点对应单元的模型上点的个数，而对于四叉树的叶节点，除了上述信息，还需要记录所有投射到此叶节点内的对应的模型上的点索引序号以及点的个数。

(3)获取轮廓点。

首先求出当前视向量V。随后从四叉树的根节点开始检测，如果节点的四个角点的法向与投影向量点积同号，说明它们的可见性状态一致，该节点内无轮廓点。丢弃该节点及其所包含的子节点，返回父节点。否则继续检测该节点的子节点，做同样的判断，直到检测到叶节点，记录该节点内点云数据轮廓点的索引。同样处理其它子节点，直到最后一个子节点。根据记录的节点内的索引计算出三维轮廓点集，记为Q。

步骤4：依次使用最小生成树算法，短边剪枝方法，光滑技术获得光滑的轮廓线。

(1)根据点的位置信息创建图G1。

在点集中确定一个起始点q₀。计算点集Q中的点与q₀的距离d，如果存在点q_i，使d＜阈值d_max(d_max的值需要根据三维模型点的相对位置情况进行确定)，则连接q_i与q₀。再将q_i作为起始点重复上述步骤。当一次连接中所有的距离都不小于阈值D时，连接结束。再确定一个起始点，重复进行连接步骤，直到所有的点都遍历完毕。

(2)根据Kruskal最小生成树算法创建图G2。

把G1图中的顶点看成独立的树，以顶点连接的线的条数作为权值，按权值从小到大选择边，所选的边连接的两个顶点应属于两颗不同的树，则成为最小生成树的一条边，并将这两颗树合并作为一颗树。重复上述步骤，直到所有顶点都在一颗树内或者有n-1条边为止。

(3)通过剪枝构造创建图G3。

检测图G2中所有点的分支数，记为m。若m＞2，则检测各分支所占的格网个数。若至少有两个分支所占的格网个数≥阈值L(阈值L根据先验知识确定。通过计算，发现L＝5时会得到比较好的检测结果)，则将所有小于阈值L的分支除去。若只有一个分支所占的格网个数≥L并且至少有两个分支所占的格网个数＜L，则保留格网个数≥L的分支并在剩下的分支中选择最长的一条。若m≤2，则保留这些分支，不做任何处理。

(4)采用平面曲线演化法光滑处理轮廓线。

将图G3的点进行投影，得到由折线表示的初始轮廓线，计算折线顶点处的两折线的角度大小及角平分线。根据折线顶点处的角度移动该顶点的位置，顶点的角度越小，该处的局部曲率越大，顶点移动的距离越大。重复上述移动过程，直到折线上每个顶点的角度大于某个阀值，表示曲线达到了要求的光滑程度。上述过程的公式为：

其中，P_i ⁿ是点P_i第n次移动的结果；k_i为P_i处的角度大小，表示曲率：N_i为点P_i角平分线方向，代表法向。

步骤5：对提取到的点云轮廓线进行离散程度筛选。

对提取到的轮廓线进行多项式拟合。多项式最高次项项数选取方法如下：

设特征点的集合为P_m，其对应的坐标集合为V(x，y，z)。分别计算x，y，z坐标的方差D(x)，D(y)和D(z)。根据先验知识以及古文物类型确定六个阈值，x坐标的阈值记为xmin，xmax(xmin＜xmax)，y坐标的阈值记为ymin，ymax(ymin＜ymax)，z坐标的阈值记为zmin，zmax(zmin＜zmax)。通过比较D(x)与xmin，xmax的大小，D(y)与ymin，ymax的大小以及D(z)与zmin，zmax的大小来确定多项式最高次项项数。

若D(z)＜zmin，当D(x)＜xmin时，x的最高次项为一次项；当xmin＜D(x)＜xmax，x的最高次项为二次项；当D(x)＞xmax时，x的最高次项为三次项，y与x同理。

若zmin＜D(z)＜zmax，当D(x)＜xmin时，x的最高次项为二次项；当xmin＜D(x)＜xmax，x的最高次项为三次项；当D(x)＞xmax时，x的最高次项为四次项，y与x同理。

若D(z)＞zmax，当D(x)＜xmin时，x的最高次项为三次项；当xmin＜D(x)＜xmax，x的最高次项为四次项；当D(x)＞xmax时，x的最高次项为五次项，y与x同理。

计算轮廓线与拟合的多项式的差值，若差值小于相应的阈值，则判断该特征为噪声，否则保留经过多项式拟合后的轮廓。相应阈值需要根据古文物的类型确定。

步骤6：对步骤5筛选出的轮廓特征进行对比。

分别记两张正射影像的轮廓特征图为C1时期特征图和C2时期特征图。

由于LiDAR点云数据量过大，因此采用分块对比方法。将点云数据分为4*4*4的模块，对应模块之间进行对比。

对比时采用基于特征的对比方法。先以C1时期的特征图为底图，C2时期的特征图为对比图进行对比。选取C2时期的特征图中的特征M。在C1时期的特征图中找到与特征M相对应的坐标区域，在该区域的Δ范围内进行微移动。移动方法如下：

设C2时期的特征图中特征M的边界坐标集合为K(X，Y，Z)。记C1时期的特征图中与K(X，Y，Z)相对应的坐标集合为K(X1，Y1，Z1)，记坐标集合K(X，Y，Z)中X坐标的最大值为maxX，Y坐标的最大值为maxY，Z坐标的最大值为maxZ。记坐标集合K(X1，Y1，Z1)中X1坐标的最大值为MaxX1，最小值为MinX1，Y1坐标的最大值为MaxY1，最小值为MinY1，Z1坐标的最大值为MaxZ1，最小值为MinZ1。

则特征M在C1时期特征图中可进行微移动的立方体范围的角点坐标为：(MaxX1+ΔX，MaxY1+ΔY，MaxZ1+Δz)，(MaxX1+ΔX，MinY1-ΔY，MaxZ1+ΔZ)，(MinX1-Δx，MaxY1+ΔY，MaxZ1+ΔZ)，(MinX1-ΔX，MinY1-ΔY，MaxZ1+Δz)，(MaxX1+ΔX，MaxY1+ΔY，MinZ1-ΔZ)，(MaxX1+ΔX，MinY1-ΔY，MinZ1-ΔZ)，(MinX1-ΔX，MaxY1+ΔY，MinZ1-ΔZ)，(MinX1-ΔX，MinY1-ΔY，MinZ1-ΔZ)

分别记为A点，B点，C点，D点，E点，F点，G点，H点。

根据先验知识以及古文物的类别设定一个趋近于0的阈值ε。当特征M在两特征图中未重合的格网个数小于ε时，则认为特征M在两特征图中完全重合。

主要步骤如下：

(1)将C2时期特征图上的(maxX，maxY，maxZ)点与C1时期特征图上的(MaxX1+ΔX，MaxY1+ΔY，MaxZ1+ΔZ)点对齐，此时特征M位于ABCDEFGH点所构成的立方体中。

(2)在X方向上，特征M每隔

个格网移动一次，在Y方向上，特征M每隔

个格网移动一次，在Z方向上，特征M每隔

个格网移动一次。每次移动后计算特征M与C1时期特征图中对应特征相差的格网个数，记为M_i。移动方向如下所示，先在x方向移动，当x方向无法继续移动时，在y方向移动一次，继续在x方向移动。当特征M在xy平面移动完毕后在z方向移动一次。

(3)基于最小二乘法进行四维曲面的多项式拟合，以坐标集合K(X1，Y1，Z1)中的X1坐标值为x坐标值，Y1坐标值为y坐标值，Z1坐标值为z坐标值，M_i值为t坐标值构造四维曲面T(x，y，z)。

(4)根据改进的随机游走算法求出曲面T中t坐标值的最小值，记为MinT。比较MinT与ε值的大小。若MinT＜ε，则特征M在两特征图中完全重合，否则特征M发生了变化。

改进的随机游走算法步骤如下：

①给定初始迭代点(x₀，y₀，z₀)，初次行走步长λ，控制精度ε1(ε1是一个非常小的正数，用于控制结束算法)。

②给定迭代控制次数I，K为当前迭代次数，置K＝1。

③当K＜I时，随机生成n个(-1，1)之间的三维向量：u_i，u_i＝(x_i1′，y_i1′，z_i1′)(-1＜x_i1′＜1，-1＜y_i1′＜1，-1＜z_i1′＜1)。并将其标准化得到

令(x_i1，y_i1，z_i1)＝(x₀，y₀，z₀)+λu_i，取min{(x₁₁，y₁₁，z₁₁)，(x₂₁，y₂₁，z₂₁)...(x_n1，y_n1，z_n1)}，记为(x_m1，y_m1，z_m1)。完成第一步游走。

④计算函数值，如果T(x_m1，y_m1，z_m1)＜T(x₀，y₀，z₀)，即找到了一个比初始值好的点，那么K重新置为1，将(x_m1，y_m1，z_m1)作为初始值，回到第②步；否则K＝K+1，回到第③步。

⑤如果连续I次都找不到更优的值，则认为，最优解就在以当前最优解为中心，当前步长为半径的I维球内(如果是三维，则刚好是空间中的球体)。此时，如果λ＜ε1，则结束算法；否则，令

回到第①步，开始新一轮游走。

(5)重复(1)-(4)的步骤，直到C2时期特征图中的所有特征都判断完毕，记作对比结果N1。

以C2时期的特征图为底图，C1时期的特征图为对比图按照步骤(1)-(5)进行一次对比。得到对比结果N2

步骤7：对步骤6得到的对比结果进行筛选。

将对比结果N1和对比结果N2进行相互对比。若特征M在对比结果N1和对比结果N2中都判断为发生了变化，则特征M为有差别特征，进行标记。若特征M在某一个对比结果中判断为发生了变化，在另一个对比结果中判断为未发生变化，则特征M为可疑特征，进行可疑特征标记。若特征M在对比结果N1和对比结果N2中都判断为未发生变化，则特征为无差别特征。

本发明的优点在于充分利用了LiDAR点云数据的优越性，LiDAR点云数据的获取具有周期短，成本低，精度高，更新时间快等特点，能够保证最终结果的精度。同时本发明克服了LiDAR点云数据缺乏纹理信息的不足，通过基于高斯球平行投影的轮廓点快速计算等一系列算法获取了物体的轮廓特征。除此之外，本发明在基于特征的变换检测方法中使用了基于最小二乘法的多项式四维曲面拟合和改进的随机游走算法，直接求出了全局最优解，避免了局部最优结果的影响。

应当理解的是，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，而所有这些改进和变换都应属于本发明所附权利要求的保护范围。

Claims (7)

1.一种运用轮廓线提取和四维曲面拟合的古文物三维模型变化检测方法，其特征在于，该方法包括以下步骤：

步骤1、获取古文物在两个不同时期的LiDAR点云数据集合，对两个不同时期的LiDAR点云数据集合进行配准；

步骤2、对LiDAR点云数据进行重分类，生成LiDAR点云模型；

步骤3、根据基于高斯球平行投影的轮廓点快速计算算法获取LiDAR点云模型的轮廓点；

步骤4、根据获取的轮廓点，依次使用最小生成树算法、短边剪枝方法、光滑技术获得光滑的点云轮廓线；

步骤5、对提取到的点云轮廓线进行离散程度筛选，去除点云轮廓线中的噪声特征，得到两个时期的轮廓特征图；

步骤6、对筛选得到的轮廓特征图进行特征对比，得到对比结果；

步骤7、对得到的对比结果进行筛选，得到差别特征；

步骤6中进行特征对比的方法为：

分别记两张正射影像的轮廓特征图为C1时期特征图和C2时期特征图；

采用分块对比方法，将点云数据分为4*4*4的模块，对应模块之间进行对比；

对比时采用基于特征的对比方法；先以C1时期的特征图为底图，C2时期的特征图为对比图进行对比；选取C2时期的特征图中的特征M，在C1时期的特征图中找到与特征M相对应的坐标区域，在该区域的Δ范围内进行微移动；移动方法如下：

设C2时期的特征图中特征M的边界坐标集合为K(X，Y，Z)；记C1时期的特征图中与K(X，Y，Z)相对应的坐标集合为K(X1，Y1，Z1)，记坐标集合K(X，Y，Z)中X坐标的最大值为maxX，Y坐标的最大值为maxY，Z坐标的最大值为maxZ；记坐标集合K(X1，Y1，Z1)中X1坐标的最大值为MaxX1，最小值为MinX1，Y1坐标的最大值为MaxY1，最小值为MinY1，Z1坐标的最大值为MaxZ1，最小值为MinZ1；则特征M在C1时期特征图中可进行微移动的立方体范围的角点坐标为：

(MaxX1+ΔX，MaxY1+ΔY，MaxZ1+ΔZ)，(MaxX1+ΔX，MinY1-ΔY，MaxZ1+ΔZ)，(MinX1-ΔX，MaxY1+ΔY，MaxZ1+ΔZ)，(MinX1-ΔX，MinY1-ΔY，MaxZ1+ΔZ)，(MaxX1+ΔX，MaxY1+ΔY，MinZ1-ΔZ)，(MaxX1+ΔX，MinY1-ΔY，MinZ1-ΔZ)，(MinX1-ΔX，MaxY1+ΔY，MinZ1-ΔZ)，(MinX1-ΔX，MinY1-ΔY，MinZ1-ΔZ)

分别记为A点，B点，C点，D点，E点，F点，G点，H点；

根据先验知识以及古文物的类别设定阈值ε，当特征M在两特征图中未重合的格网个数小于ε时，则认为特征M在两特征图中完全重合；

步骤6中进行特征对比的具体步骤为：

步骤61、将C2时期特征图上的(maxX，maxY，maxZ)点与C1时期特征图上的(MaxX1+ΔX，MaxY1+ΔY，MaxZ1+ΔZ)点对齐，此时特征M位于ABCDEFGH点所构成的立方体中；

步骤62、在X方向上，特征M每隔

个格网移动一次，在Y方向上，特征M每隔

个格网移动一次，在Z方向上，特征M每隔

个格网移动一次；

每次移动后计算特征M与C1时期特征图中对应特征相差的格网个数，记为M_i；移动方向如下所示，先在x方向移动，当x方向无法继续移动时，在y方向移动一次，继续在x方向移动；当特征M在xy平面移动完毕后在z方向移动一次；

步骤63、基于最小二乘法进行四维曲面的多项式拟合，以坐标集合K(X1，Y1，Z1)中的X1坐标值为x坐标值，Y1坐标值为y坐标值，Z1坐标值为z坐标值，M_i值为t坐标值构造四维曲面T(x，y，z)；

步骤64、根据改进的随机游走算法求出曲面T中t坐标值的最小值，记为MinT；比较MinT与ε值的大小，若MinT<ε，则特征M在两特征图中完全重合，否则特征M发生了变化；

步骤65、重复步骤61-步骤64，直到C2时期特征图中的所有特征都判断完毕，记作对比结果N1；

以C2时期的特征图为底图，C1时期的特征图为对比图按照步骤61-步骤65进行一次对比，得到对比结果N2；

步骤64中改进的随机游走算法步骤为：

①给定初始迭代点(x₀,y₀,z₀)，初次行走步长λ，控制精度ε1；

②给定迭代控制次数I，K_d为当前迭代次数，置K_d＝1；

③当K_d<I时，随机生成n个(-1,1)之间的三维向量u_i,u_i＝(x_i1′,y_i1′,z_i1′)，-1<x_i1′<1,-1<y_i1′<1,-1<z_i1′<1；并将其标准化得到

令(x_i1,y_i1,z_i1)＝(x₀,y₀,z₀)+λu_i,取min{(x₁₁,y₁₁,z₁₁),(x₂₁,y₂₁,z₂₁)…(x_n1,y_n1,z_n1)},记为(x_m1,y_m1,z_m1)；完成第一步游走；

④计算函数值，如果T(x_m1,y_m1,z_m1)<T(x₀,y₀,z₀)，即找到了一个比初始值好的点，那么K_d重新置为1，将(x_m1,y_m1,z_m1)作为初始值，回到第②步；否则K_d＝K_d+1，回到第③步；

⑤如果连续I次都找不到更优的值，则认为，最优解就在以当前最优解为中心，当前步长为半径的I维球内；此时，如果λ<ε1,则结束算法；否则，令

回到第①步，开始新一轮游走。

2.根据权利要求1所述的运用轮廓线提取和四维曲面拟合的古文物三维模型变化检测方法，其特征在于，步骤1中对两个不同时期的LiDAR点云数据集合进行配准的方法为：通过平移、旋转操作使两个LiDAR点云数据集合的位置重合。

3.根据权利要求1所述的运用轮廓线提取和四维曲面拟合的古文物三维模型变化检测方法，其特征在于，步骤2中进行重分类生成LiDAR点云模型的方法为：

记LiDAR点云数据的个数为M_p个，点云数据集合的长，宽，高之比为L：W：H；将点云数据集合按照

的格网进行划分；统计每个格网中点云数据的个数，记为m；若m≥1，则记该格网中心点为有效点云数据，否则该格网中无有效点云数据。

4.根据权利要求1所述的运用轮廓线提取和四维曲面拟合的古文物三维模型变化检测方法，其特征在于，步骤3中获取LiDAR点云模型轮廓点的方法为：

步骤31、通过分析LiDAR点云模型上各点的K近邻的协方差矩阵来计算该点的法向量；记LiDAR点云模型上某一点为P，找出点P邻域点的集合，记为N_p，计算出N_p的重心，以此建立起一个协方差矩阵C：

其中，

协方差矩阵C的最小特征值对应的特征向量就是该点的法向量；

步骤32、构建高斯球，利用四叉树的方法建立分层数据结构进行数据存储；将LiDAR点云模型中的点的法向映射到单位矢量球上；

设N＝(N_x，N_y，N_z)是点云模型上点P的正则化法向量，通过以下公式获得点P在高斯球上的映射位置：

建立以高斯球球心为坐标系原点的球坐标系，用经度1d和纬度ud来表示法向的球坐标；通过转换公式将高斯球坐标转换为以高斯球球心为坐标原点的直角坐标，转换公式如下：

x＝cos(ud)cos(ld)，y＝cos(ud)sin(ld)，z＝sin(ld)

利用四叉树模型建立分层数据结构组织点云数据；

步骤33、获取轮廓点：首先求出当前视向量V；随后从四叉树的根节点开始检测，如果节点的四个角点的法向与投影向量点积同号，说明它们的可见性状态一致，该节点内无轮廓点；丢弃该节点及其所包含的子节点，返回父节点；否则继续检测该节点的子节点，做同样的判断，直到检测到叶节点，记录该节点内点云数据轮廓点的索引；同样处理其它子节点，直到最后一个子节点；根据记录的节点内的索引计算出三维轮廓点集，记为Q。

5.根据权利要求4所述的运用轮廓线提取和四维曲面拟合的古文物三维模型变化检测方法，其特征在于，步骤32中利用四叉树模型建立分层数据结构组织点云数据的方法为：

对高斯球沿经纬线进行层次划分，得到四边区域，称为单元；利用四叉树结构存储这个层次划分，四叉树的根节点对应包含整个球的单元；对单元逐层划分，如果单元内的点小于一个设定的阈值或者单元的四个角点经纬度差别小于一个设定的阈值，则不再划分，这些不再做划分的单元就是四叉树的叶节点；对于四叉树的中间节点，只记录其四个儿子节点及其四个角点的信息，以及投射到此节点对应单元的模型上点的个数，而对于四叉树的叶节点，除了上述信息，还记录所有投射到此叶节点内的对应的模型上的点索引序号以及点的个数。

6.根据权利要求4所述的运用轮廓线提取和四维曲面拟合的古文物三维模型变化检测方法，其特征在于，步骤4中获取点云轮廓线的方法为：

步骤41、根据轮廓点的位置信息创建图G1；

在三维轮廓点集Q中确定一个起始点q₀；计算点集Q中的点与q₀的距离d，如果存在点q_i，使d<阈值d_max，d_max的值需根据三维模型点的相对位置情况进行确定，则连接q_i与q₀；再将q_i作为起始点重复上述步骤；当一次连接中所有的距离都不小于阈值D时，连接结束；再确定一个起始点，重复进行连接步骤，直到所有的点都遍历完毕；

步骤42、根据Kruskal最小生成树算法创建图G2；

把图G1中的顶点看成独立的树，以顶点连接的线的条数作为权值，按权值从小到大选择边，所选的边连接的两个顶点属于两颗不同的树，则成为最小生成树的一条边，并将这两颗树合并作为一颗树；重复上述步骤，直到所有顶点都在一颗树内或者有n_p-1条边为止；

步骤43、通过剪枝构造创建图G3；

检测图G2中所有点的分支数，记为m；若m>2,则检测各分支所占的格网个数；若至少有两个分支所占的格网个数≥阈值L，阈值L根据先验知识确定，则将所有小于阈值L的分支除去；若只有一个分支所占的格网个数≥L并且至少有两个分支所占的格网个数<L，则保留格网个数≥L的分支并在剩下的分支中选择最长的一条；若m≤2，则保留这些分支，不做任何处理；

步骤44、采用平面曲线演化法光滑处理轮廓线；

将图G3的点进行投影，得到由折线表示的初始轮廓线，计算折线顶点处的两折线的角度大小及角平分线；根据折线顶点处的角度移动该顶点的位置，顶点的角度越小，该处的局部曲率越大，顶点移动的距离越大；重复上述移动过程，直到折线上每个顶点的角度大于设定阀值，表示曲线达到了要求的光滑程度；上述过程的公式为：

其中，P_i ^ns是点P_i第ns次移动的结果；K_i为P_i处的角度大小，表示曲率：N_i为点P_i角平分线方向，代表法向。

7.根据权利要求6所述的运用轮廓线提取和四维曲面拟合的古文物三维模型变化检测方法，其特征在于，步骤5中进行离散程度筛选的方法为：

对提取到的轮廓线进行多项式拟合，多项式最高次项项数选取方法如下：

设特征点的集合为P_m，其对应的坐标集合为V_p(x，y，z)；分别计算x,y,z坐标的方差D(x)，D(y)和D(z)；根据先验知识以及古文物类型确定六个阈值，x坐标的阈值记为xmin，xmax，xmin<xmax，y坐标的阈值记为ymin,ymax，ymin<ymax，z坐标的阈值记为zmin，zmax，zmin<zmax；通过比较D(x)与xmin，xmax的大小，D(y)与ymin，ymax的大小以及d(z)与zmin，zmax的大小来确定多项式最高次项项数；

若D(z)<zmin，当D(x)<xmin时，x的最高次项为一次项；当xmin<D(x)<xmax，x的最高次项为二次项；当D(x)>xmax时，x的最高次项为三次项，y与x同理；

若zmin<D(z)<zmax，当D(x)<xmin时，x的最高次项为二次项；当xmin<D(x)<xmax，x的最高次项为三次项；当D(x)>xmax时，x的最高次项为四次项，y与x同理；

若D(z)>zmax，当D(x)<xmin时，x的最高次项为三次项；当xmin<D(x)<xmax，x的最高次项为四次项；当D(x)>xmax时，x的最高次项为五次项，y与x同理；

计算轮廓线与拟合的多项式的差值，若差值小于相应的阈值，则判断该特征为噪声，否则保留经过多项式拟合后的轮廓，相应阈值根据古文物的类型确定。

Images