基于直接点云间距统计采样的三维点云形变追踪方法

目录

• 1 引言
• • 1.1 现有研究综述
  • 1.2 本文研究创新点
• 2 研究方法
• • 2.1 高密度点云生成
  • 2.2 点云预处理
  • 2.3 空间配准
  • 2.4 三维点云变形分析
  • • 2.4.1 直接点云对点云（C2C）距离估算
    • 2.4.2 M3C2距离估算
  • 2.5 统计采样
• 3 试验与结果
• 4 结果分析讨论
• 5 研究局限性与未来研究展望
• 6 结论

摘要：依托激光扫描仪或多视图摄影测量技术生成的基础设施高密度三维（3D）点云，具备作为无损检测信息数据源的巨大潜力。随着相关技术日趋成熟，现已能够重建出毫米级精度、具备照片级真实感的三维模型，完全满足结构检测与评估工作的精度要求。但人工解析点云数据往往耗时耗力，且无法直观反映结构受力性能与健康状态，因此亟需研发新型算法实现三维模型的高效分析。本文提出一种全新三维点云变化分析方法，可通过局部空间分析手段追踪结构随时间产生的微小位移。该方法融合逐点直接距离度量与统计采样策略，从中提取结构变形信息。通过识别并跟踪点云时序变化，能够量化结构力学变形，并同步计算对应的应变与应力；所得测量结果可用于评判结构使用工况及剩余承载能力。文中同时开展一系列室内试验，对所提方法进行验证。试验结果表明，该方法测量精度可达±0.2毫米（95%置信区间），适用于各类基础设施检测场景下高精度、自动化的几何分析与形变识别。课题组后续将开展两方面延伸研究：一是将该算法与有限元模型自动更新技术相结合；二是开展现场实测，验证该测量方法在工程实景下的适用性。

1 引言

  近年来，三维扫描技术已广泛应用于结构检测与状态评估领域，成为现代化检测手段。目前有多种方式可生成高密度三维扫描数据（即点云）[1]。激光雷达（LiDAR）是最主流的点云采集手段，而摄影测量是其替代方案——该技术通过图像测算完成三维点云重建。其中，密集运动恢复结构算法（DSfM）仅需低成本设备，即可生成精度与激光雷达扫描成果相当的点云数据[2,3]。结合计算机视觉分析技术，三维点云能够完善结构全生命周期状态评估体系，提升工程结构的抗灾耐久性能。   对工程人员而言，依托点云变化分析算法，能够识别并量化结构构件随时间产生的变形、损伤区域，还可基于监测数据更新竣工三维结构模型[4]。在土木工程领域，该类三维变化分析方法的典型应用场景为结构力学变形测量。通过测算变形量，可直接推导结构弹性、非弹性应变与累积损伤程度，以此校核结构正常使用极限状态性能[5]。

1.1 现有研究综述

  文献[6]与[7]的综述系统梳理了三维点云变化检测算法的理论基础。文献[6]指出，三维变化检测主要分为两大类方法：（1）几何对比法：通过测算三维几何差异实现变化识别；（2）几何-光谱联合分析法：融合几何信息与光谱信息开展变化分析。若场景三维几何形态复杂，且需要在单点尺度上对比两套三维模型的形变差异（即本文研究场景），几何对比法是量化单点形变的最优方案[7]。   现有相关研究中，已有学者利用激光雷达技术测量桥梁挠度[8]；也有算法通过向变形三维曲面拟合标准几何基元实现梁体变形建模[9]。多项研究验证了激光雷达等现代遥感技术在桥梁净空检测中的精度与效率：通过处理高密度三维点云数据，该类手段测量精度可媲美传统测绘设备（全球导航卫星系统GNSS、全站仪），且现场作业成本与工作量大幅降低[10-13]。文献[14]将实测激光雷达点云与平面、圆柱等标准几何模型开展逐点偏差分析，识别地震作用下建筑产生的变形。   但上述全部方法均存在局限：需要对曲面进行网格划分、多项式曲线拟合，或预先假定结构初始变形量。上述任一环节都会引入测量误差，且难以拓展至大型、复杂结构的检测场景[5,15]。

1.2 本文研究创新点

  本文提出并验证一种全新的点云变形测量方法。该方法结合逐点直接距离度量与统计采样技术，无需网格划分、曲线拟合，直接基于原始点云测算变形，最大限度保留原始数据完整性、降低测量误差。该测量方案的精度可接近传统专业监测仪器；同时测量以局部区域为单位，能够适配大尺度结构与任意不规则变形形态。   本文另一项研究贡献为：对比分析两种距离度量指标，并论证二者在本算法框架下的适配性优劣。

2 研究方法

图1展示了整套测量流程。首先分别生成未变形构件与变形构件的三维点云用于对比；若不存在结构未变形状态的点云基准数据，可由计算机辅助设计（CAD）图纸生成合成点云[16]。对点云执行预处理裁剪与去噪操作，所采用的去噪方法为文献[17]提出的统计离群点剔除算法。完成预处理后，将两组点云配准对齐，再计算两组点云间各点的距离，得到两组点云间逐点相对变形量，通常以热力图形式可视化展示（图2）。

随后对感兴趣区域（ROI）内的点测量值开展统计采样，消除距离计算过程中出现的离群值与错误测量带来的干扰。

2.1 高密度点云生成

  三维点云的生成方式多种多样[1]，激光扫描是目前最常用手段，且图1所示整套处理流程兼容任意点云生成方案。本文采用二维数字图像结合密集运动恢复结构（DSfM）技术生成点云[18]。DSfM是一种摄影测量技术，可生成分辨率高、密度大的三维点云，数据质量可与激光扫描结果持平。 传统软拷贝摄影测量需要预先已知相机三维坐标与姿态（即相机标定），而运动恢复结构（SfM）无需先验标定，依靠多幅图像间匹配显著特征点，通过冗余光束平差法同步求解相机参数（相机姿态）与场景三维几何形态[19]。该流程可精确解算图像集中所有相机姿态，同时生成稀疏（无照片真实感）三维点集。   重建流程首先提取图像特征描述子并完成图像间匹配，目的是在整套图像数据中识别同一实体特征。本文选用尺度不变特征变换（SIFT）特征[20]，该特征对图像缩放、仿射变换、噪声与光照变化均具备不变性。   获取相机参数后，执行第二步流程提升点云密度，还原具备照片真实感的三维模型。“密集”代表提取输入二维图像全部像素信息参与重建，与之相对，稀疏点云仅使用提取出的特征点进行三维重建。现有密集重建算法种类较多[18]，本文选用半全局匹配算法（SGM）[21]，该算法精度高、可拓展性强，能够生成超高密度点云，完整保留模型细节。

  该算法先为每一张视图重建深度（视差）图（图3），再采用文献[22]的分层策略融合所有深度图，形成全局统一的三维表征，高效消除多视图重叠带来的数据冗余。深度图指一对图像中同一三维点投影位置的差值。关于DSfM完整流程与可选方案的详细说明，可参考文献[2]、[18]。

2.2 点云预处理

  生成高密度点云后，需开展预处理工作。首先裁剪无关场景数据，降低计算开销；随后采用文献[17]的统计方法对点云平滑去噪。   计算每个点

p

i

p_i

pi​与其k个近邻点欧氏距离

d

i

d_i

di​的均值

μ

\\mu

μ与标准差

σ

\\sigma

σ；所有距离落在

μ

±

α

σ

\\mu \\pm \\alpha \\sigma

μ±ασ区间外的点判定为离群噪声点，从三维点云中剔除。本文按照文献[17]推荐参数设置：

k

=

30

k=30

k=30，

α

=

1

\\alpha=1

α=1，适配各处不均匀的点云密度。

2.3 空间配准

  采用迭代最近点算法（ICP）[24]对变形、未变形两组点云自动对齐配准。ICP算法核心目标：求解源点云

P

P

P中点

p

i

p_i

pi​与待配准点云

Q

Q

Q中点

q

i

q_i

qi​之间的刚性变换矩阵。通过近邻搜索与欧氏距离计算，匹配两组点云中距离最近的对应点。为提升搜索效率，采用k-d树空间划分算法[25]加速最近点检索，将单次搜索复杂度降至

O

(

log

⁡

n

)

O(\\log n)

O(logn)。   为求解

p

i

p_i

pi​、

q

i

q_i

qi​之间的旋转矩阵

R

R

R与平移向量

t

t

t，在全部对应点约束下最小化式(1)误差函数；

p

i

p_i

pi​、

q

i

q_i

qi​均以三维空间坐标

(

x

,

y

,

z

)

(x,y,z)

(x,y,z)数值表示。ICP算法整体复杂度为

O

[

w

×

n

×

l

o

g

(

n

)

]

O[w ×n ×log(n)]

O[w×n×log(n)]，

n

n

n为点数量，

w

w

w为点云数量。求解得到变换参数后，将两组点云统一映射至同一参考坐标系完成配准。

E

(

R

,

t

)

=

min

⁡

R

,

t

∑

i

∥

p

i

−

(

R

q

i

+

t

)

∥

2

E(R, t)=\\min _{R, t} \\sum_{i}\\left\\| p_{i}-\\left(R q_{i}+t\\right)\\right\\| ^{2}

E(R,t)=R,tmin​i∑​∥pi​−(Rqi​+t)∥2

  ICP本质是局部收敛算法，配准结果高度依赖初始位姿。为提升刚性变换精度，可先手动选取至少4组两组点云内同名点完成粗配准，根据对齐精度需求选择粗配或精配方案。本文结合手动粗配+ICP自动精配的组合方案优化测量精度；同时在刚性变换前完成点云去噪，降低离群点对点匹配与配准结果的干扰。

2.4 三维点云变形分析

  两组三维点云完成空间配准后，即可开展对比形变分析。本文采用逐点距离估算法，在单点尺度量化两组点云的位移差，最大限度保留原始数据信息。文中对比两种距离度量方法：基于最近点的直接点云对点云（C2C）对比法[26]、多尺度模型对点云对比法（M3C2）[27]。

2.4.1 直接点云对点云（C2C）距离估算

  该方法是快速、直接的三维对比算法，无需对数据集网格划分，也无需计算点法向量。C2C算法采用八叉树结构划分三维空间，基于八叉树检索已配准两组点云的最近点，再通过豪斯多夫距离计算两组点云间整体距离[28]。 设两个有限点集

A

=

{

a

1

,

.

.

.

,

a

p

}

A=\\{a_{1}, …, a_{p}\\}

A={a1​,…,ap​}、

B

=

{

b

1

,

.

.

.

,

b

q

}

B=\\{b_{1}, …, b_{q}\\}

B={b1​,…,bq​}，集合间豪斯多夫距离定义为式(2)：

H

(

A

,

B

)

=

max

⁡

(

h

(

A

,

B

)

,

h

(

B

,

A

)

)

(2)

H\\left( A, B\\right) =\\operatorname{max}\\left( h\\left( A, B\\right) , h\\left( B, A\\right) \\right) \\tag{2}

H(A,B)=max(h(A,B),h(B,A))(2) 式中

h

h

h为单向豪斯多夫距离：

h

(

A

,

B

)

=

max

⁡

a

∈

A

min

⁡

b

∈

B

∥

a

−

b

∥

h(A, B)=\\max _{a \\in A} \\min _{b \\in B}\\| a-b\\|

h(A,B)=a∈Amax​b∈Bmin​∥a−b∥

∥

⋅

∥

\\|\\cdot\\|

∥⋅∥代表向量范数。单向豪斯多夫距离

h

(

A

,

B

)

h(A,B)

h(A,B)不具备对称性（存在曲面满足

h

(

A

,

B

)

≠

h

(

B

,

A

)

h(A, B) ≠h(B, A)

h(A,B)=h(B,A)）；双向豪斯多夫距离

H

(

A

,

B

)

H(A,B)

H(A,B)取两组单向距离的最大值，具备对称性。 采用豪斯多夫距离做三维点云对比的优势：无需为集合A、B内的点建立显式一一对应关系。该距离仅表征纯几何差异，不预设点云密度均匀分布。

2.4.2 M3C2距离估算

  M3C2距离计算流程：首先从基准点云P中按指定最小三维间距下采样选取“核心点”，该步骤仅用于缩短计算耗时，为可选操作。本文未执行核心点下采样，本研究不考量计算成本。 随后基于局部协方差矩阵求解每个点的法向量，协方差矩阵为半正定对称矩阵，对应最小特征值的特征向量即为曲面法向量[29]：

C

3

×

3

=

1

k

∑

i

=

1

k

(

p

i

−

p

‾

)

(

p

i

−

p

‾

)

T

,

p

‾

=

1

k

∑

i

=

1

k

p

i

C_{3 × 3}=\\frac{1}{k} \\sum_{i=1}^{k}\\left(p_{i}-\\overline{p}\\right)\\left(p_{i}-\\overline{p}\\right)^{T},\\quad \\overline{p}=\\frac{1}{k} \\sum_{i=1}^{k} p_{i}

C3×3​=k1​i=1∑k​(pi​−p​)(pi​−p​)T,p​=k1​i=1∑k​pi​ 式中

C

C

C为协方差矩阵，

p

ˉ

\\bar{p}

pˉ​为点

p

i

p_i

pi​的k近邻局部几何中心。对协方差矩阵执行奇异值分解（SVD）[30]，得到特征向量

v

2

、

v

1

、

v

0

v_2、v_1、v_0

v2​、v1​、v0​与对应特征值

λ

2

≥

λ

1

≥

λ

0

\\lambda_2 ≥\\lambda_1 ≥\\lambda_0

λ2​≥λ1​≥λ0​；其中

v

0

v_0

v0​近似为点

p

i

p_i

pi​处曲面法向量

N

N

N，

v

2

、

v

1

v_2、v_1

v2​、v1​张成

p

i

p_i

pi​点切平面。该计算等价于最小二乘平面拟合。 特征值

λ

0

\\lambda_0

λ0​定量表征沿法向的数据离散程度，可作为局部噪声水平估计，代表k个邻域点相对局部拟合平面的偏离量；

λ

0

\\lambda_0

λ0​越小，平面拟合质量越高。   完成点法向量求解后，沿核心点

i

i

i的法向量

N

N

N构建圆柱形投影区域，圆柱直径

d

d

d为用户定义的投影尺度；限定圆柱最大长度以降低计算量。分别提取两组点云落在圆柱范围内的点子集，数量记为

n

1

、

n

2

n_1、n_2

n1​、n2​，形成两组空间位置分布。取分布均值得到两组点云沿法向的平均位置

i

1

、

i

2

i_1、i_2

i1​、i2​，两点间欧氏距离即为M3C2距离度量；两组分布标准差可作为沿法向点云粗糙度

σ

1

(

d

)

、

σ

2

(

d

)

\\sigma_1(d)、\\sigma_2(d)

σ1​(d)、σ2​(d)的局部估计。 均值与标准差易受离群值干扰，因此实际采用距离分布的中位数计算

i

1

、

i

2

i_1、i_2

i1​、i2​，采用四分位距表征粗糙度。

  图5为本文两种几何距离计算方法原理示意图。C2C方法主要缺陷：对噪声、离群点敏感，计算结果受点云空间采样率影响大；M3C2通过更精细的局部采样与均值平滑机制抑制点云噪声与离群点干扰，试验结果将直观验证该优势。

2.5 统计采样

  逐点测量方案存在固有缺陷：点云中每个点都会生成一组独立变形测量值，局部噪声、错误近邻匹配会扭曲测量结果，对于主变形方向垂直的曲面，该问题尤为突出；两组点云重叠缺失区域也会产生无效数据。为解决上述问题，本文采用统计采样方法，从局部区域海量单点测量值中提取可靠变形结果。   针对构件需要测算挠度的位置，在点云中划定感兴趣区域（ROI），汇总该区域全部挠度测量值，绘制数据频率分布直方图（图6），直观分析数据分布特性（正态分布、离群值、偏态等）；同时直方图可图形化归纳海量点云数据，是处理大规模点云数据集的重要工具。   如图6所示，受噪声与部分点错误近邻匹配影响，感兴趣区域内测量值分布通常不服从标准统计分布。因此本文基于感兴趣区域内数据的统计百分位数阈值筛选有效测量值；该阈值需通过试验标定，且C2C、M3C2两种算法的最优阈值差异显著。

3 试验与结果

  为评估本文所提测量方法（C2C、M3C2两种算法）的精度，本文开展一系列弯曲试验。试验采用三种不同厚度（3.18 mm、6.40 mm、12.7 mm）矩形截面铝制试件，一方面探究不同弹性挠度下该方法的测量精度，另一方面研究构件尺寸带来的影响。按照简支梁弹性梁理论，将试件分级加载至弯曲屈服承载力的75%，每级加载量为预估最大挠度的10%。所有弯曲试验均在Tinius-Olsen H50KT万能材料试验机（UTM）上完成，采用标准三点弯曲加载装置。   首先利用密集运动恢复结构（DSfM）技术生成未发生弯曲变形试件的点云；随后按预设10%分级逐级加载，每完成一级加载后，再次通过DSfM采集变形试件的点云。图像采集设备为尼康D800E相机搭配50 mm f/1.8定焦镜头，采用Agisoft Photoscan（1.2.4版本）软件完成DSfM重建全流程。   将变形与未变形状态的点云导入CloudCompare（2.6.1版本）完成配准与分析。三点弯曲试验的两处简支支座作为手动对齐基准点，辅助ICP算法完成后续配准，试验过程中支座固定无位移。分别计算未变形点云与各级挠度变形点云之间的C2C距离与M3C2距离。沿试件总跨长以1/8跨长为间隔划定多个感兴趣区域（ROI），单个ROI宽度6.35 mm，每个ROI平均包含44000个点。各ROI布设位置如图7所示。受三点弯曲加载装置遮挡，跨中位置未设置感兴趣区域。   将各ROI对应的点云与逐点距离计算结果导入MATLAB，开展去噪与统计分析。基于直方图分析得到的挠度估算值，与试验机实测挠度进行对比。受试验机硬件限制，仅能在试件跨中获取实测挠度，其余位置的理论挠度采用经典欧拉–伯努利梁理论计算得到。对于C2C算法，通过试验标定得到最优统计阈值为95百分位数；对于M3C2算法，则采用50百分位数（中位数）作为筛选阈值。

  图8展示C2C、M3C2两种方法测得挠度与理论预期挠度的对比曲线，3.18 mm厚试件完整测量误差见表1、表2。整体规律为：挠度数值越小，测量绝对误差越大。例如试件最大挠度6.3 mm工况下，C2C平均测量误差仅0.27 mm；当挠度降至3.78 mm时，C2C平均误差上升至0.38 mm。其余两种厚度试件的试验结果均呈现相同变化规律。

6.40 mm厚度试件的试验结果如图9所示，该组试验明确体现了本文测量方法在微小挠度（小于1 mm）工况下的精度下限。当跨中挠度仅0.3 mm时，受固有测量误差影响，算法无法准确还原试件挠曲线形态（图9a）；12.7 mm厚试件在小挠度工况下也出现相同现象。微小挠度工况下，测量绝对误差与大挠度工况基本持平，因此相对误差会显著增大。

3.18 mm试件全部测试结果的误差分布见图10，其余试件误差分布规律与之相近。两种测量方法的误差均值均小于1 mm，且所有试验工况下M3C2测量结果稳定性更强。综合全部试验数据，C2C算法整体平均误差0.57 mm，标准差0.54；M3C2算法整体平均误差0.40 mm，标准差0.47。

4 结果分析讨论

  图8展示了两种算法测得的挠度与试验机实测得到的理论挠度对比结果。3.18mm厚度试件完整的测量误差数据如表1、表2所示。总体规律为：挠度越小，测量绝对误差越大。例如，试件最大挠度6.3mm（本试验最大变形工况）时，C2C算法平均测量误差仅0.27mm；当跨中挠度降至3.78mm时，C2C平均误差上升至0.38mm。其余两种厚度试件的试验结果也呈现出相同规律。   6.40mm厚试件的试验结果如图9所示，该组试验揭示了本测量方法在挠度小于1mm工况下的精度极限。当跨中挠度仅0.3mm时，受测量固有误差影响，算法无法还原出准确的挠曲线形态（图9a）。12.7mm厚试件在微小变形工况下也观测到了相同现象。对于小挠度工况，测量绝对误差与大变形工况基本处于同一水平，这就导致相对误差大幅升高。   3.18mm试件全部试验的误差分布如图10所示，其余试件的误差分布规律与之类似。两种测量方法的误差均值均低于1mm；在全部试验工况中，M3C2的测量结果稳定性更好。整合全部试验数据后，C2C算法整体平均误差为0.57mm，标准差0.54；M3C2算法整体平均误差为0.40mm，标准差0.47。

5 研究局限性与未来研究展望

  挠度越小，测量绝对误差越大，造成该现象最主要的原因是变形点云与未变形基准点云在ICP配准过程中产生的配准误差。当变形量极小时，两组点云形态高度近似，难以完成高精度刚性配准，进而增大测量误差。该配准偏差带来的影响如图11所示：虽然测得的挠曲线整体形态准确，但试件整体围绕中点出现顺时针旋转，这正是点云对齐偏差造成的。

  其余误差来源还包括：密集运动恢复结构（DSfM）重建过程自身存在缺陷、C2C算法近邻点匹配出错、M3C2法中法向量估算存在偏差。总体而言，直接点云对比（C2C）对点云粗糙度（噪声）、离群点、配准不确定性以及点云采样间距均较为敏感，以上因素都会引发变形计算错误，这也凸显了两种距离计算方法的核心差异。   对于包含海量点云的三维模型，直接开展变化分析与实时可视化渲染会带来极大计算开销。计算机中央处理器（CPU）与图形处理器（GPU）硬件资源存在上限，海量三维点云极易超出常规计算机内存容量。此外，现有可视化方案通常会对三维模型进行抽稀或栅格化处理以缩减数据体量，该操作会损失测量精度，丢失细微结构特征信息[32]。   后续研究将重点实现核外空间数据结构算法，解除数据处理对数据集总规模的限制，突破内存容量瓶颈。同时采用优化并行计算框架与GPU加速方案，实现海量三维点云模型的高效分析。

6 结论

  综合试验结果表明，本文提出的测量方法具备逐点测算任意复杂结构变形的能力，且无需对曲面进行网格划分或曲线拟合。两种算法的平均测量误差均小于1毫米，测量结果离散程度相近。   对比C2C与M3C2两种算法可以发现：M3C2的测量精度更高、鲁棒性更强，原因是其在变形分析过程中内置了空间可变置信区间作为变化检测阈值。此外，M3C2在存在遮挡、数据缺失场景下表现更稳定——该算法仅在沿法向与对比三维模型存在相交区域时才计算形变，不会输出无依据的变化量。   无论采用何种距离度量指标，该方法达到的测量精度足以支撑遥感监测相关应用，尤其适用于需要精准捕捉微小局部变形、无法布设传统传感器的工程场景。例如构件局部翼缘屈曲、冲击损伤检测等场景，若采用全局曲面拟合类方法，会掩盖这类局部损伤特征。

下一阶段研究将在可控现场工况下开展实测，明确现场环境下该测量方法的实际测量精度。当前正在推进、后续计划开展的研究工作如下：