一、CPD 非刚性配准算法基础
1.1 点云配准概述
点云配准作为计算机视觉和三维数据处理领域的关键技术,旨在寻求一种空间变换,将不同来源的点云数据在同一坐标系下精确对齐。其核心目标是通过计算合适的旋转矩阵 R R R和平移向量 t t t,使源点云 P = { p 1 , p 2 , ⋯ , p n } P = \{p_1, p_2, \cdots, p_n\} P={p1,p2,⋯,pn}经过变换 T ( p ) = R p + t T(p) = Rp + t T(p)=Rp+t后,与目标点云 Q = { q 1 , q 2 , ⋯ , q m } Q = \{q_1, q_2, \cdots, q_m\} Q={q1,q2,⋯,qm}达到最佳重合状态,通常以某种距离度量(如欧氏距离)来衡量重合程度 。例如,在三维重建中,从不同角度扫描物体得到的多个点云需要配准,才能拼接成完整的三维模型;在自动驾驶中,车辆通过激光雷达实时获取的环境点云,需要与预先构建的地图点云进行配准,以确定车辆的位置和姿态。
根据点云在配准过程中是否发生形状变化,点云配准可分为刚性配准和非刚性配准。刚性配准假设物体在变换过程中保持刚性,即仅发生平移和旋转,不改变物体的形状和大小 。刚性配准方法在处理刚体物体点云时表现出色,如机械零件的点云配准。常见的刚性配准算法是迭代最近点(ICP)算法,它通过不断迭代寻找最近点对,并根据这些点对计算最优的刚性变换,使源点云逐步逼近目标点云。然而,在现实世界中,许多物体存在复杂的非刚性形变,如人体的运动、生物组织的生长和变形等。非刚性配准则允许物体在配准过程中发生形状变化,能够处理包含反射、旋转、缩放、平移以及各种非线性变形的情况 。非刚性配准通过建立更复杂的变换模型,如薄板样条函数、径向基函数或基于物理模型的变形场,来准确描述物体的形变,从而实现点云的精确对齐。
点云配准技术在众多领域有着广泛的应用。在医学领域,点云配准用于医学图像分析,如配准不同模态(如 MRI、CT)的医学图像,结合不同成像方式的优势,为医生提供更全面的信息,辅助疾病诊断和治疗规划;在工业检测中,通过将测量得到的工业零部件点云与标准模型点云进行配准,可以检测出产品的尺寸偏差、表面缺陷等质量问题,提高产品质量控制水平;在机器人视觉中,点云配准帮助机器人识别和抓取复杂环境中的物体,通过配准实时获取的物体点云与预先存储的模型点云,机器人能够确定物体的姿态和位置,实现精确操作;在文化遗产保护领域,点云配准用于文物的三维数字化重建,将从不同角度扫描文物得到的点云数据进行配准和拼接,还原文物的完整形态,为文物的保护、修复和研究提供重要依据 。
1.2 CPD 算法核心思想
CPD 算法的核心思想是将点云配准问题巧妙地转化为概率密度估计问题,借助高斯混合模型(GMM)来对目标点云进行概率建模,并结合正则化形变场实现点云的非刚性变换,从而达成高精度的点云配准 。
在 CPD 算法中,目标点云被视作由高斯混合模型生成的观测数据。具体而言,假设源点云为 X = { x n } n = 1 N X = \{x_n\}_{n = 1}^N X={xn}n=1N,目标点云为 Y = { y m } m = 1 M Y = \{y_m\}_{m = 1}^M Y={ym}m=1M 。将源点云 X X X经过非刚性变换 T ( x n ; θ ) T(x_n;\theta) T(xn;θ)后的点作为高斯混合模型中各个高斯分布的中心,其中 θ \theta θ表示变换参数。每个高斯分布的概率密度函数为 p ( y ∣ n ) = N ( y ∣ T ( x n ; θ ) , σ 2 I ) p(y|n) = \mathcal{N}(y | T(x_n;\theta), \sigma^2 I) p(y∣n)=N(y∣T(xn;θ),σ2I),这里 N \mathcal{N} N表示正态分布, y y y是目标点云中的点, σ 2 \sigma^2 σ2为方差, I I I是单位矩阵 。同时,为了处理可能存在的离群点,引入一个额外的均匀分布 P ( N + 1 ) = ω 1 − ω P(N + 1) = \frac{\omega}{1 - \omega} P(N+1)=1−ωω,其中 ω \omega ω代表离群点在点云中所占的比例。这样,目标点云的概率密度函数可表示为 p ( y ) = ∑ n = 1 N + 1 P ( n ) p ( y ∣ n ) p(y) = \sum_{n = 1}^{N + 1} P(n) p(y|n) p(y)=∑n=1N+1P(n)p(y∣n),其中 P ( n ) P(n) P(n)为每个成分的权重,当 n ≤ N n \leq N n≤N时, P ( n ) = 1 N P(n) = \frac{1}{N} P(n)=N1,表示源点云经过变换后生成目标点云的概率;当 n = N + 1 n = N + 1 n=N+1时, P ( N + 1 ) P(N + 1) P(N+1)对应离群点的分布权重。通过这种方式,将点云配准问题转化为求解高斯混合模型参数 θ \theta θ和 σ 2 \sigma^2 σ2,使得目标点云的似然函数最大化的问题,即寻找最优的变换参数,使源点云经过变换后生成目标点云的概率最大 。
为了求解上述问题,CPD 算法采用期望最大化(EM)算法。在 E 步,根据当前的高斯混合模型参数,计算目标点云中每个点 y m y_m ym来自各个高斯分布(即与源点云经过变换后的各个点的对应关系)的概率 P ( n ∣ y m ) P(n|y_m) P(n∣ym) 。这个概率反映了目标点与源点云变换后各点的匹配程度,通过它可以确定点云之间的软对应关系,而非传统刚性配准中寻找精确的硬对应点对。在 M 步,利用 E 步计算得到的概率,更新高斯混合模型的参数,包括变换参数 θ \theta θ和方差 σ 2 \sigma^2 σ2,以最大化目标点云的似然函数。通过不断迭代 E 步和 M 步,高斯混合模型的参数逐渐收敛,从而得到最优的非刚性变换参数,实现源点云与目标点云的配准 。
为了确保非刚性变换的合理性和稳定性,CPD 算法引入了正则化形变场 。该形变场对源点云的变形进行约束,避免出现不合理的扭曲和变形。具体来说,通过定义一个正则化项,通常基于形变场的平滑性或其他几何约束条件,将其加入到目标函数中。在迭代优化过程中,不仅要最大化目标点云的似然函数,还要使正则化项的值保持在合理范围内,以平衡数据拟合和形变场的平滑性。例如,常见的正则化项可以基于薄板样条函数或其他径向基函数构建,通过调整正则化参数,控制形变场的平滑程度。这样,CPD 算法能够在考虑点云分布概率的同时,准确地描述物体的复杂形变,实现高质量的非刚性点云配准 。
1.3 相关理论基础
1.3.1 高斯混合模型(GMM)
高斯混合模型(Gaussian Mixture Model,GMM)是一种强大的概率模型,它通过多个高斯分布的线性组合来描述复杂的数据分布 。在 CPD 非刚性配准算法中,GMM 起着关键作用,用于对目标点云的概率分布进行建模。
从数学定义上看,GMM 的概率密度函数可表示为:
p ( x ) = ∑ k = 1 K π k N ( x ∣ μ k , Σ k ) p(x) = \sum_{k = 1}^{K} \pi_k \mathcal{N}(x | \mu_k, \Sigma_k) p(x)=∑k=1KπkN(x∣μk,Σk)
其中, K K K表示高斯分布的个数,即混合成分的数量; π k \pi_k πk是第 k k k个高斯分布的权重,满足 ∑ k = 1 K π k = 1 \sum_{k = 1}^{K} \pi_k = 1 ∑k=1Kπk=1且 0 ≤ π k ≤ 1 0 \leq \pi_k \leq 1 0≤πk≤1,它反映了每个高斯分布在混合模型中所占的比重; N ( x ∣ μ k , Σ k ) \mathcal{N}(x | \mu_k, \Sigma_k) N(x∣μk,Σk)是第 k k k个高斯分布的概率密度函数, μ k \mu_k μk是该高斯分布的均值向量,决定了分布的中心位置, Σ k \Sigma_k Σk是协方差矩阵,决定了分布的形状和尺度 。例如,当 x x x是二维数据时, μ k = ( μ k 1 , μ k 2 ) \mu_k = (\mu_{k1}, \mu_{k2}) μk=(μk1,μk2), Σ k \Sigma_k Σk是一个 2 × 2 2 \times 2 2×2的矩阵,其元素决定了数据在两个维度上的分布特性。
在 CPD 算法中,将目标点云视为由 GMM 生成的数据。假设源点云 X = { x n } n = 1 N X = \{x_n\}_{n = 1}^N X={xn}n=1N经过非刚性变换 T ( x n ; θ ) T(x_n;\theta) T(xn;θ)后得到的点作为 GMM 中各个高斯分布的中心,即 μ k = T ( x k ; θ ) \mu_k = T(x_k;\theta) μk=T(xk;θ) 。这里的非刚性变换 T ( x n ; θ ) T(x_n;\theta) T(xn;θ)是一个依赖于参数 θ \theta θ的函数,它可以描述源点云的复杂形变,使得源点云能够更好地与目标点云的分布相匹配 。例如,对于一个三维点云,非刚性变换可能包括平移、旋转、缩放以及更复杂的非线性变形,通过调整 θ \theta θ参数,可以实现源点云在空间中的各种形变,从而使以形变后的源点云为中心的高斯分布能够准确地拟合目标点云的分布 。
GMM 的优势在于其强大的拟合能力。理论上,只要混合成分 K K K足够大,GMM 可以逼近任何连续的概率分布 。这使得它在处理复杂的数据分布时表现出色,能够捕捉到数据中的各种特征和规律。在点云配准中,目标点云的分布往往是复杂多样的,可能包含噪声、离群点以及各种非线性形变带来的影响,GMM 能够有效地对这样的复杂分布进行建模,为后续的配准计算提供坚实的基础 。
1.3.2 期望最大化(EM)算法
期望最大化(Expectation-Maximization,EM)算法是一种用于求解含有隐变量的概率模型参数的迭代优化算法 。在 CPD 算法中,通过 EM 算法来估计高斯混合模型的参数,从而实现点云配准。
EM 算法的基本思想是通过迭代的方式,交替执行 E 步(期望步)和 M 步(最大化步),逐步逼近最优的模型参数 。具体步骤如下:
E 步(期望步):在这一步中,基于当前估计的模型参数,计算每个数据点属于各个高斯分布的概率,即后验概率 。对于 CPD 算法中的高斯混合模型,设当前估计的参数为 θ ( t ) \theta^{(t)} θ(t)( t t t表示迭代次数),对于目标点云中的每个点 y m y_m ym,计算它来自第 n n n个高斯分布的概率 P ( n ∣ y m , θ ( t ) ) P(n|y_m, \theta^{(t)}) P(n∣ym,θ(t)),根据贝叶斯公式,有:
P ( n ∣ y m , θ ( t ) ) = P ( n ) p ( y m ∣ n , θ ( t ) ) ∑ n = 1 N + 1 P ( n ) p ( y m ∣ n , θ ( t ) ) P(n|y_m, \theta^{(t)}) = \frac{P(n) p(y_m|n, \theta^{(t)})}{\sum_{n = 1}^{N + 1} P(n) p(y_m|n, \theta^{(t)})} P(n∣ym,θ(t))=∑n=1N+1P(n)p(ym∣n,θ(t))P(n)p(ym∣n,θ(t))
其中, P ( n ) P(n) P(n)是第 n n n个高斯分布的先验概率(在 CPD 算法中,当 n ≤ N n \leq N n≤N时, P ( n ) = 1 N P(n) = \frac{1}{N} P(n)=N1;当 n = N + 1 n = N + 1 n=N+1时, P ( N + 1 ) = ω 1 − ω P(N + 1) = \frac{\omega}{1 - \omega} P(N+1)=1−ωω, ω \omega ω为离群点比例), p ( y m ∣ n , θ ( t ) ) p(y_m|n, \theta^{(t)}) p(ym∣n,θ(t))是在参数 θ ( t ) \theta^{(t)} θ(t)下,点 y m y_m ym来自第 n n n个高斯分布的概率密度 。这个后验概率 P ( n ∣ y m , θ ( t ) ) P(n|y_m, \theta^{(t)}) P(n∣ym,θ(t))反映了点 y m y_m ym与源点云经过变换后的各个点(即高斯分布中心)的匹配程度,通过它可以建立点云之间的软对应关系,即每个目标点都以一定的概率与源点云变换后的多个点相对应,而不是像传统刚性配准那样寻找唯一的硬对应点对 。
M 步(最大化步):利用 E 步计算得到的后验概率,更新高斯混合模型的参数,以最大化目标点云的似然函数 。对于 CPD 算法,需要更新的参数包括非刚性变换参数 θ \theta θ和高斯分布的方差 σ 2 \sigma^2 σ2 。通过对似然函数关于这些参数求偏导并令其为零,可得到参数更新的公式 。例如,对于非刚性变换参数 θ \theta θ的更新,通常通过最小化一个包含数据拟合项和正则化项的目标函数来实现,以平衡数据匹配和形变场的合理性 。在更新参数后,得到新的参数估计 θ ( t + 1 ) \theta^{(t + 1)} θ(t+1) 。
迭代过程:重复执行 E 步和 M 步,直到模型参数收敛,即相邻两次迭代中参数的变化小于某个预设的阈值 。此时,得到的参数估计即为高斯混合模型的最优参数,对应的非刚性变换参数 θ \theta θ能够使源点云经过变换后与目标点云达到最佳的匹配状态,从而实现点云配准 。
EM 算法的收敛性基于其单调增加似然函数的特性 。在每次迭代中,M 步保证了似然函数在当前 E 步计算的后验概率下达到局部最大值,而 E 步则基于更新后的参数重新计算后验概率,为下一次 M 步提供更准确的估计 。通过不断迭代,似然函数逐渐增大,最终收敛到一个局部最优解 。虽然 EM 算法不能保证收敛到全局最优解,但在实际应用中,对于许多问题,其收敛到的局部最优解已经能够满足需求,并且通过合理选择初始参数和多次运行算法取最优结果等方式,可以在一定程度上提高得到全局最优解的概率 。
1.3.3 正则化形变场
在 CPD 非刚性配准算法中,正则化形变场是实现准确非刚性点云配准的关键要素之一,它主要用于约束源点云在非刚性变换过程中的形变,确保变换的合理性和稳定性 。
在点云配准中,非刚性变换可能会导致源点云产生不合理的扭曲和变形,使得配准结果失去物理意义或准确性 。为了避免这种情况,引入正则化形变场对源点云的变形进行约束 。正则化形变场通常基于某种几何或物理先验知识来构建,其核心思想是通过添加一个正则化项到目标函数中,来限制形变场的某些特性 。例如,常见的基于平滑性约束的正则化项,其目的是使形变场在空间上尽可能平滑,避免出现突变和不连续的变形 。从数学角度看,假设形变场可以表示为一个函数 u ( x ) u(x) u(x),它描述了源点云中点 x x x的位移,基于薄板样条(Thin Plate Spline,TPS)函数的正则化项可以定义为:
Ω ( u ) = ∫ R d ( ∥ ∇ 2 u ( x ) ∥ 2 ) d x \Omega(u) = \int_{\mathbb{R}^d} \left( \left\|\nabla^2 u(x) \right\|^2 \right) dx Ω(u)=∫Rd( ∇2u(x) 2)dx
其中, ∇ 2 \nabla^2 ∇2是拉普拉斯算子, d d d是点云的维度(通常为 3), ∥ ∇ 2 u ( x ) ∥ 2 \left\|\nabla^2 u(x) \right\|^2 ∇2u(x) 2衡量了形变场 u ( x ) u(x) u(x)的二阶导数的范数,积分表示对整个点云空间进行求和 。这个正则化项的作用是惩罚形变场的剧烈变化,使得形变场在空间上保持平滑 。当 Ω ( u ) \Omega(u) Ω(u)的值较小时,说明形变场变化较为平缓,源点云的变形更加合理;反之,若 Ω ( u ) \Omega(u) Ω(u)的值较大,则表示形变场存在较大的突变,可能导致源点云的变形不合理 。
在 CPD 算法的目标函数中,将正则化项与数据拟合项(即基于高斯混合模型的似然函数项)相结合,得到完整的目标函数:
L ( θ , σ 2 , u ) = ∑ m = 1 M log ( ∑ n = 1 N + 1 P ( n ) p ( y m ∣ n , θ , u ) ) − λ Ω ( u ) L(\theta, \sigma^2, u) = \sum_{m = 1}^{M} \log \left( \sum_{n = 1}^{N + 1} P(n) p(y_m|n, \theta, u) \right) - \lambda \Omega(u) L(θ,σ2,u)=∑m=1Mlog(∑n=1N+1P(n)p(ym∣n,θ,u))−λΩ(u)
其中, λ \lambda λ是正则化参数,用于平衡数据拟合和正则化约束的相对重要性 。当 λ \lambda λ取值较大时,正则化约束的作用更强,形变场更加平滑,但可能会牺牲一定的数据拟合精度;当 λ \lambda λ取值较小时,数据拟合的作用更强,可能会导致形变场不够平滑,出现不合理的变形 。因此,合理选择 λ \lambda λ的值对于获得准确且合理的配准结果至关重要 。在实际应用中,通常需要通过实验来确定最优的 λ \lambda λ值,以适应不同的点云数据和配准需求 。通过这种方式,正则化形变场能够在保证源点云与目标点云数据拟合的同时,有效约束源点云的非刚性变换,实现高质量的非刚性点云配准 。
二、CPD 算法原理深度剖析
2.1 概率建模
CPD 算法的基石是将点云配准问题巧妙地转化为概率密度估计问题,其核心在于运用高斯混合模型(GMM)对目标点云进行精确的概率建模 。
在这一建模过程中,假设源点云为 X = { x n } n = 1 N X = \{x_n\}_{n = 1}^N X={xn}n=1N,目标点云为 Y = { y m } m = 1 M Y = \{y_m\}_{m = 1}^M Y={ym}m=1M 。CPD 算法创新性地将目标点云视为由高斯混合模型生成的观测数据,其中每个高斯分布的中心对应源点云经过非刚性变换 T ( x n ; θ ) T(x_n;\theta) T(xn;θ)后的位置 。这里的非刚性变换 T ( x n ; θ ) T(x_n;\theta) T(xn;θ)是一个依赖于参数 θ \theta θ的函数,它能够描述源点云的复杂形变,使源点云通过形变更好地与目标点云的分布相匹配 。例如,在处理人体运动点云时,非刚性变换可以准确地模拟人体关节的弯曲、伸展等复杂动作引起的点云形状变化 。
高斯混合模型的概率密度函数表示为:
p ( y ) = ∑ n = 1 N + 1 P ( n ) p ( y ∣ n ) p(y) = \sum_{n = 1}^{N + 1} P(n) p(y|n) p(y)=∑n=1N+1P(n)p(y∣n)
其中, P ( n ) P(n) P(n)为均匀分布的权重 。当 n ≤ N n \leq N n≤N时, P ( n ) = 1 N P(n) = \frac{1}{N} P(n)=N1,这表示假设源点云经过变换后生成目标点云的各个点是等概率的 。 p ( y ∣ n ) p(y|n) p(y∣n)是第 n n n个高斯分布的概率密度,具体形式为 p ( y ∣ n ) = N ( y ∣ T ( x n ; θ ) , σ 2 I ) p(y|n) = \mathcal{N}(y | T(x_n;\theta), \sigma^2 I) p(y∣n)=N(y∣T(xn;θ),σ2I),这里 N \mathcal{N} N表示正态分布, y y y是目标点云中的点, T ( x n ; θ ) T(x_n;\theta) T(xn;θ)是源点云 x n x_n xn经过非刚性变换后的位置,作为高斯分布的均值, σ 2 \sigma^2 σ2为方差,用于控制高斯分布的宽度, I I I是单位矩阵 。方差 σ 2 \sigma^2 σ2的大小直接影响高斯分布的形状,较小的 σ 2 \sigma^2 σ2会使高斯分布更加集中,表明点云数据的分布更为紧凑;较大的 σ 2 \sigma^2 σ2则使高斯分布更加分散,适用于描述数据分布较为广泛的情况 。
在实际的点云数据中,不可避免地会存在离群点,这些离群点可能是由于数据采集过程中的噪声干扰、设备误差等原因产生的 。如果不加以处理,离群点会对配准结果产生严重的负面影响,导致配准精度下降甚至配准失败 。为了有效处理离群点,CPD 算法添加了一个特殊的均匀分布项 。设离群点在点云中所占的比例为 ω \omega ω,则添加的均匀分布项为 P ( N + 1 ) = ω 1 − ω P(N + 1) = \frac{\omega}{1 - \omega} P(N+1)=1−ωω 。这个均匀分布项的作用是为离群点提供一个合理的概率解释,使得离群点在概率模型中也能得到恰当的处理 。例如,当某个点被判定为离群点的概率较高时,它更有可能来自这个均匀分布,从而减少对整体配准过程的干扰 。通过这种方式,CPD 算法能够在复杂的点云数据中,准确地建模目标点云的分布,为后续基于期望最大化算法的参数估计和点云配准奠定坚实的基础 。
2.2 期望最大化(EM)算法在 CPD 中的应用
期望最大化(EM)算法是 CPD 非刚性配准算法中用于求解高斯混合模型参数的核心迭代优化方法,通过交替执行 E 步(期望步)和 M 步(最大化步),逐步逼近最优的模型参数,实现源点云与目标点云的精确配准 。
在 E 步中,主要任务是根据当前估计的高斯混合模型参数,计算目标点云中每个点 y m y_m ym来自各个高斯分布(即与源点云经过变换后的各个点的对应关系)的后验概率 P ( n ∣ y m , θ ( t ) ) P(n|y_m, \theta^{(t)}) P(n∣ym,θ(t)),其中 θ ( t ) \theta^{(t)} θ(t)表示第 t t t次迭代时的模型参数 。根据贝叶斯公式,后验概率的计算如下:
P ( n ∣ y m , θ ( t ) ) = P ( n ) p ( y m ∣ n , θ ( t ) ) ∑ n = 1 N + 1 P ( n ) p ( y m ∣ n , θ ( t ) ) P(n|y_m, \theta^{(t)}) = \frac{P(n) p(y_m|n, \theta^{(t)})}{\sum_{n = 1}^{N + 1} P(n) p(y_m|n, \theta^{(t)})} P(n∣ym,θ(t))=∑n=1N+1P(n)p(ym∣n,θ(t))P(n)p(ym∣n,θ(t))
在 CPD 算法的设定下,当 n ≤ N n \leq N n≤N时, P ( n ) = 1 N P(n) = \frac{1}{N} P(n)=N1,代表源点云经过变换后生成目标点云的点的先验概率;当 n = N + 1 n = N + 1 n=N+1时, P ( N + 1 ) = ω 1 − ω P(N + 1) = \frac{\omega}{1 - \omega} P(N+1)=1−ωω,用于处理离群点, ω \omega ω为离群点在点云中所占的比例 。 p ( y m ∣ n , θ ( t ) ) p(y_m|n, \theta^{(t)}) p(ym∣n,θ(t))是在参数 θ ( t ) \theta^{(t)} θ(t)下,点 y m y_m ym来自第 n n n个高斯分布的概率密度,其形式为 p ( y m ∣ n , θ ( t ) ) = N ( y m ∣ T ( x n ; θ ( t ) ) , σ 2 I ) p(y_m|n, \theta^{(t)}) = \mathcal{N}(y_m | T(x_n;\theta^{(t)}), \sigma^2 I) p(ym∣n,θ(t))=N(ym∣T(xn;θ(t)),σ2I) 。这个后验概率反映了目标点 y m y_m ym与源点云经过变换后的各个点(即高斯分布中心)的匹配程度,通过它建立了点云之间的软对应关系 。例如,在处理人体点云配准时,对于目标点云中的某个点,通过计算后验概率,可以确定它与源点云经过不同变形后的多个点的对应可能性,而不是像刚性配准那样寻找唯一的对应点,这种软对应关系能够更好地适应非刚性形变的复杂情况 。
进入 M 步,其核心是利用 E 步计算得到的后验概率,更新高斯混合模型的参数,以最大化目标点云的似然函数 。在 CPD 算法中,需要更新的参数主要包括非刚性变换参数 θ \theta θ和高斯分布的方差 σ 2 \sigma^2 σ2 。对于非刚性变换参数 θ \theta θ的更新,通常是通过最小化一个包含数据拟合项和正则化项的目标函数来实现 。数据拟合项基于目标点云与源点云经过变换后的匹配程度,即通过 E 步得到的后验概率来衡量;正则化项则用于约束形变场,确保源点云的变形合理 。例如,基于薄板样条函数的正则化项会惩罚形变场的剧烈变化,使得形变场更加平滑 。通过对目标函数关于 θ \theta θ求偏导并令其为零,可得到 θ \theta θ的更新公式 。对于方差 σ 2 \sigma^2 σ2的更新,也有相应的数学推导公式,其目的是调整高斯分布的宽度,以更好地拟合目标点云的分布 。例如,当点云数据的分布较为集中时,方差 σ 2 \sigma^2 σ2会相应减小,使高斯分布更加紧凑;当点云数据分布较为分散时,方差 σ 2 \sigma^2 σ2会增大,以适应数据的分布特性 。
在整个迭代过程中,不断重复执行 E 步和 M 步 。每一次迭代,E 步基于当前的模型参数计算后验概率,为 M 步提供更准确的点云对应关系估计;M 步则根据 E 步的结果更新模型参数,使得高斯混合模型能够更好地拟合目标点云 。随着迭代次数的增加,模型参数逐渐收敛,即相邻两次迭代中参数的变化小于某个预设的阈值 。当达到收敛条件时,得到的参数估计即为高斯混合模型的最优参数,对应的非刚性变换参数 θ \theta θ能够使源点云经过变换后与目标点云达到最佳的匹配状态,从而完成点云配准 。例如,在实际应用中,经过多次迭代后,源点云逐渐变形并精确地对齐到目标点云,实现了高质量的非刚性点云配准 。
2.3 刚性变换与非刚性变换模型参数求解
在点云配准中,刚性变换和非刚性变换是两种重要的变换类型,它们的模型参数求解方式存在显著差异,理解这些差异对于准确实现点云配准至关重要。
刚性变换模型假设物体在变换过程中保持刚性,即只发生平移和旋转,不发生形状改变。其变换模型通常可以表示为 T ( p ) = R p + t T(p) = Rp + t T(p)=Rp+t,其中 R R R是 3 × 3 3\times3 3×3的旋转矩阵, t t t是三维平移向量, p p p是点云中的点 。求解刚性变换模型参数的过程,通常是基于点云之间的对应关系,通过最小化某种距离度量来实现 。以迭代最近点(ICP)算法为例,该算法通过不断迭代寻找最近点对,并根据这些点对计算最优的刚性变换 。具体来说,在每次迭代中,首先为源点云中的每个点在目标点云中寻找最近点,形成点对;然后,基于这些点对,利用最小二乘法求解旋转矩阵 R R R和平移向量 t t t,使得源点云经过变换后与目标点云之间的均方误差最小 。通过不断重复这个过程,直到变换参数收敛,从而得到最优的刚性变换,实现点云配准 。例如,在对机械零件的点云进行配准时,由于零件在不同测量视角下主要发生平移和旋转,刚性变换模型能够很好地描述这种变化,通过 ICP 算法可以准确求解刚性变换参数,实现零件点云的精确对齐 。
非刚性变换模型则允许物体在配准过程中发生形状变化,能够处理包含反射、旋转、缩放、平移以及各种非线性变形的复杂情况 。在 CPD 算法中,非刚性变换采用薄板样条(Thin Plate Spline,TPS)函数来构建形变场,从而实现对源点云的复杂变形描述 。薄板样条函数是一种基于径向基函数的插值方法,其核心思想是通过在源点云的控制点上定义位移,然后利用薄板样条函数将这些位移插值到整个点云空间,从而得到源点云的非刚性变换 。具体来说,设源点云为 X = { x n } n = 1 N X = \{x_n\}_{n = 1}^N X={xn}n=1N,目标点云为 Y = { y m } m = 1 M Y = \{y_m\}_{m = 1}^M Y={ym}m=1M,薄板样条函数可以表示为 u ( x ) = s u m i = 1 N w i p h i ( l e f t ∣ x − x i r i g h t ∣ ) + a 0 + a 1 x 1 + a 2 x 2 + a 3 x 3 u(x) = \\sum_{i = 1}^{N} w_i \\phi(\\left\\|x - x_i\\right\\|) + a_0 + a_1 x_1 + a_2 x_2 + a_3 x_3 u(x)=sumi=1Nwiphi(left∣x−xiright∣)+a0+a1x1+a2x2+a3x3(对于三维点云),其中 u ( x ) u(x) u(x)表示点 x x x的位移, w i w_i wi是权重系数, p h i ( l e f t ∣ x − x i r i g h t ∣ ) \\phi(\\left\\|x - x_i\\right\\|) phi(left∣x−xiright∣)是径向基函数(在薄板样条中通常为 r 2 l o g r r^2\\log r r2logr, r = l e f t ∣ x − x i r i g h t ∣ r = \\left\\|x - x_i\\right\\| r=left∣x−xiright∣), a 0 , a 1 , a 2 , a 3 a_0, a_1, a_2, a_3 a0,a1,a2,a3是线性项系数 。在 CPD 算法中,通过最小化包含数据拟合项(基于高斯混合模型的似然函数)和正则化项(基于薄板样条函数的平滑约束)的目标函数来求解薄板样条函数的参数(即权重系数 w i w_i wi和线性项系数 a 0 , a 1 , a 2 , a 3 a_0, a_1, a_2, a_3 a0,a1,a2,a3) 。例如,在处理人体运动点云时,人体的关节弯曲、肌肉变形等复杂形变无法用刚性变换描述,而通过薄板样条构建的非刚性变换模型能够准确模拟这些形变,通过求解薄板样条函数的参数,使源点云能够根据人体的实际变形情况进行调整,从而与目标点云实现高精度配准 。
刚性变换模型参数求解主要基于点云之间的对应关系,通过最小化简单的距离度量(如均方误差)来实现,适用于物体形状不变的情况;而非刚性变换模型参数求解则通过构建复杂的形变场(如基于薄板样条函数),并最小化包含数据拟合和正则化约束的目标函数来实现,能够处理物体的复杂形变 。在实际应用中,应根据点云数据的特点和配准需求,选择合适的变换模型和参数求解方法,以实现准确的点云配准 。
2.4 算法流程与实现细节
CPD 非刚性配准算法的实现是一个复杂且精细的过程,其流程主要包括数据初始化、基于期望最大化(EM)算法的迭代优化以及迭代停止条件的判断等关键环节。
在数据初始化阶段,需要对源点云 X = { x n } n = 1 N X = \{x_n\}_{n = 1}^N X={xn}n=1N和目标点云 Y = { y m } m = 1 M Y = \{y_m\}_{m = 1}^M Y={ym}m=1M进行预处理 。首先,对原始点云数据进行去噪处理,去除由于数据采集设备噪声、环境干扰等因素引入的离群点和噪声点,以提高数据质量 。常见的去噪方法包括统计滤波、双边滤波等 。例如,统计滤波通过计算点云邻域内点的统计信息(如均值、方差),将偏离统计特征的点视为噪声点并去除 。其次,对去噪后的点云进行归一化处理,将点云数据的坐标范围统一到 [0, 1] 区间,这样可以消除点云数据在尺度上的差异,使算法在迭代过程中更加稳定 。归一化处理通常通过计算点云坐标的最大值和最小值,然后对每个点的坐标进行线性变换实现 。同时,需要初始化高斯混合模型(GMM)的参数,包括方差 σ 2 \sigma^2 σ2和离群点比例 ω \omega ω 。方差 σ 2 \sigma^2 σ2的初始值一般根据点云数据的分布特性进行设定,若点云数据分布较为集中, σ 2 \sigma^2 σ2可设置较小的值;若分布较为分散,则 σ 2 \sigma^2 σ2设置较大的值 。离群点比例 ω \omega ω的初始值可以根据经验或对数据的初步分析进行设定,一般在 0.05 - 0.2 之间 。例如,对于噪声较多的点云数据, ω \omega ω可适当取较大值 。
在迭代优化阶段,CPD 算法基于 EM 算法进行迭代 。在每次迭代中,E 步首先根据当前估计的高斯混合模型参数 θ ( t ) \theta^{(t)} θ(t)( t t t表示迭代次数),计算目标点云中每个点 y m y_m ym来自各个高斯分布(即与源点云经过变换后的各个点的对应关系)的后验概率 P ( n ∣ y m , θ ( t ) ) P(n|y_m, \theta^{(t)}) P(n∣ym,θ(t)) 。这个过程涉及到复杂的数学计算,需要根据贝叶斯公式,结合高斯分布的概率密度函数进行计算 。M 步则利用 E 步计算得到的后验概率,更新高斯混合模型的参数,包括非刚性变换参数 θ \theta θ和方差 σ 2 \sigma^2 σ2 。对于非刚性变换参数 θ \theta θ的更新,通常通过最小化一个包含数据拟合项和正则化项的目标函数来实现 。数据拟合项基于目标点云与源点云经过变换后的匹配程度,通过 E 步得到的后验概率来衡量;正则化项则用于约束形变场,确保源点云的变形合理 。在实现过程中,为了提高计算效率,可以采用一些优化技术,如矩阵运算的并行化处理 。利用多线程或 GPU 并行计算,加速矩阵乘法、加法等运算,从而缩短每次迭代的计算时间 。同时,合理的数据结构设计也能提高算法效率,例如使用 KD 树等数据结构来快速查找最近邻点,减少计算最近邻点对的时间开销 。
迭代优化的停止条件是判断算法是否收敛的关键依据 。常见的停止条件有两种 。一是基于参数变化的判断,当相邻两次迭代中高斯混合模型的参数变化小于某个预设的阈值时,认为算法收敛 。例如,对于非刚性变换参数 θ \theta θ,若 ∥ θ ( t + 1 ) − θ ( t ) ∥ < ϵ \left\|\theta^{(t + 1)} - \theta^{(t)}\right\| < \epsilon θ(t+1)−θ(t) <ϵ( ϵ \epsilon ϵ为预设的阈值,通常取值较小,如 1 0 − 5 10^{-5} 10−5),则满足停止条件 。二是基于目标函数值变化的判断,当相邻两次迭代中目标函数值的变化小于某个预设的阈值时,停止迭代 。即若 ∣ L ( θ ( t + 1 ) , σ 2 ( t + 1 ) , u ( t + 1 ) ) − L ( θ ( t ) , σ 2 ( t ) , u ( t ) ) ∣ < δ \left|L(\theta^{(t + 1)}, \sigma^{2(t + 1)}, u^{(t + 1)}) - L(\theta^{(t)}, \sigma^{2(t)}, u^{(t)})\right| < \delta L(θ(t+1),σ2(t+1),u(t+1))−L(θ(t),σ2(t),u(t)) <δ( δ \delta δ为预设的阈值,如 1 0 − 4 10^{-4} 10−4),则停止迭代 。在实际应用中,还可以结合最大迭代次数进行判断,若迭代次数达到预先设定的最大值(如 100 次),即使未满足上述收敛条件,也停止迭代,以避免算法陷入无限循环 。通过合理设置停止条件和实现细节优化,CPD 算法能够在保证配准精度的前提下,高效地完成非刚性点云配准任务 。