一种基于加速高斯混合模型的空间目标重建方法与流程

本发明涉及航天器近距离感知领域，特别涉及一种基于加速高斯混合模型的空间目标重建方法。

背景技术

点云可用于对目标外形进行直观描述，但对于进一步接近或抓捕等任务来说，这种形式的描述远不能满足任务需求，还需要在此基础之上进一步获取目标外形的解析描述。

高斯混合模型(gaussianmixturemodel,gmm)作为一种功能强大的参数化概率模型，核心思想在于将目标建模为多个高斯分布加权和的形式，已广泛应用于背景建模、图像分割、超像素生成、语音识别等各个领域，取得了丰硕的成果。

jian等最早在两点云配准过程中使用gmm对点云进行建模，最终将配准问题转化为两个gmm的对齐问题。在建模过程中，点云中的每个点使用一个高斯分量进行表示，各高斯分量的权重和协方差矩阵相同，均值为各点在3d空间中的坐标。eckart等和evangelidis等分别在两点云配准和多点云配准中使用gmm对点云进行建模。二者相同之处在于均使用少量gmm表示整个点云，而并非像jian等那样使用gmm表示每个点，这大大减少了所需gmm的数量；不同之处在于前者并不限定协方差矩阵的形式，而后者为简化计算过程使用各向同性的协方差矩阵。在evangelidis等基础之上，danelljan等还在对点云建模过程中引入除点云3d坐标之外的颜色信息，通过提高gmm维数以期提高点云配准的精度。目标点云中往往包含数量巨大的点，可能有几万、几十万、乃至上百万的量级，使用原始gmm对目标点云进行重建会消耗大量的时间，因此需要对其进行加速。

在上文提到的文献中，只有eckart等利用gmm算法的并行特性，利用gpu对建模过程进行加速，但文献中并未涉及加速过程的实现细节。此外，eckart等还在近期提出一种自顶向下的分层gmm模型，用于加速点云配准之前的建模过程。

使用原始高斯混合模型(gaussianmixturemodel，gmm)对目标点云进行重建主要有以下两个方面的缺点：目标点云中包含数量巨大的点，可能有几万、几十万乃至上百万的量级，重建会消耗大量的时间，需要对其进行加速；重建过程中的不确定性对最终结果影响很大，鲁棒性较差，如对参数初值的依赖性较强，初值选取不好容易导致参数估计无法收敛。因此，需要一种新的方法，解决目前方法所存在的计算量大、鲁棒性弱等问题。

技术实现要素：

本发明的目的是通过以下技术方案实现的。

本发明设计了一种基于加速高斯混合模型的空间目标重建方法。该重建方法针对航天器近距离感知中的目标三维结构重建问题进行研究，解决了目前方法所存在的计算量大、鲁棒性弱等问题。

一种基于加速高斯混合模型的空间目标重建方法，包括：

读入点云数据；

根据给定高斯混合模型分量的个数k，使用k-means++算法从点云中筛选出k个点作为初始聚类中心；

使用k-means算法将点云聚类为k个簇；

根据每个簇中所包含点的情况，分别对权重、均值、协方差矩阵初始化，作为高斯混合模型参数的初值；

使用高斯混合模型通过期望最大化算法对参数进行迭代更新，获取最终的估计值。

优选地，读入点云数据包括：

点云通过一系列离散点对目标进行描述，离散点数据采用结构体数组的形式存储；

点云读入过程中，将离散点数据以数组结构体形式存储在宿主内存；

gpu对全局内存进行合并访问；

将以数组结构体形式存储的数据上传至设备内存以供gpu访问。

优选地，k-means++算法包括：

使用k-means++算法从点云中筛选出k个点作为初始聚类中心；

k-means++算法的伪代码在cpu中进行计算；

将初始聚类中心以结构体数组形式存储于宿主内存并上传至设备内存以供gpu访问。

优选地，k-means算法包括：

计算各点云到各聚类中心的距离，将得到的距离矩阵存储于设备内存；

更新各点云分类标记，在gpu中进行计算，将所得标记存储于设备内存，将设备内存中的标记下载至宿主内存以供cpu访问；

更新各聚类中心，在cpu中进行计算，将更新后的聚类中心重新上传至设备内存；

收敛判断，满足任意一个条件即停止迭代。

优选地，收敛判断的两个条件分别为：

最大迭代次数超过阈值；

分类标记发生变化的点在点云中所占的比重低于阈值。

优选地，对均值、权重、协方差矩阵初始化包括：

均值利用k-means算法最后一次迭代过程中计算获得的聚类中心进行初始化；

使用与更新各聚类中心相同的方法对协方差矩阵初始化，在cpu中进行计算；

使用与更新各聚类中心相同的方法对权重初始化，cpu中进行计算；

将宿主内存中存储的初始化后的均值、权重、协方差矩阵上传至设备内存以供gpu访问。

优选地，对参数进行迭代更新，包括：

计算责任矩阵；

更新权重；

更新均值；

更新协方差矩阵；

收敛判断。

优选地，收敛判断使用两种准则进行收敛性判断，满足以下任意一个即停止迭代：

最大迭代次数超过阈值；

本次迭代中似然函数值的变化小于阈值。

优选地，期望最大化算法包括如下步骤：

使用参数计算期望；

rnk表示第n个样本由第k个高斯分量产生的概率，表示第k个高斯分量，其均值协方差矩阵表示第k个高斯分量的权重且满足约束样本由包含k个分量的高斯混合模型参数产生。

优选地，期望最大化算法包括如下步骤：

通过最大化期望算法的完备数据似然函数，对参数进行更新；

本发明的优点在于：

本发明提供的一种基于加速高斯混合模型的空间目标重建方法，针对航天器近距离感知中的目标三维重建问题进行研究，针对目标点云中包含数量巨大的点，可能有几万、几十万乃至上百万的量级，结合算法流程设计和cpu-gpu异构实现，计算速度快；重建过程中，鲁棒性强。

附图说明

通过阅读下文具体实施方式的详细描述，各种其他的优点和益处对于本领域普通技术人员将变得清楚明了。附图仅用于示出具体实施方式的目的，而并不认为是对本发明的限制。而且在整个附图中，用相同的参考符号表示相同的部件。在附图中：

附图1示出了根据本发明实施方式的加速高斯混合模型算法流程图；

附图2示出了根据本发明实施方式的加速高斯混合模型cuda实现示意图；

附图3示出了根据本发明实施方式的点云的结构体数组与数组结构体存储。

具体实施方式

下面将参照附图更详细地描述本公开的示例性实施方式。虽然附图中显示了本公开的示例性实施方式，然而应当理解，可以以各种形式实现本公开而不应被这里阐述的实施方式所限制。相反，提供这些实施方式是为了能够更透彻地理解本公开，并且能够将本公开的范围完整的传达给本领域的技术人员。

根据本发明的实施方式，提出一种基于加速高斯混合模型的空间目标重建方法。本发明技术方案的算法流程如图1所示，分为五个步骤：首先读入待重建的点云数据；根据给定高斯混合模型(gaussianmixturemodel,gmm)分量的个数k，使用k-means++算法从点云中筛选出k个点作为初始聚类中心；随后使用k-means算法将点云聚类为k个簇；然后根据每个簇中所包含点的情况，分别计算其均值μk、协方差矩阵σk和权重πk作为gmm参数θ的初值；最后使用gmm通过期望最大化算法(em算法，expectation-maximizationalgorithm)对参数θ进行迭代更新，获取最终的估计值。

假设d维空间中样本由包含k个分量的gmm产生，则其概率密度函数可以表示为

其中表示第k个高斯分量，其均值协方差矩阵表示第k个高斯分量的权重，且满足约束表示该模型所有参数的集合。

给定由该gmm产生的样本集直观上可以使用最大似然(maximumlikelihood，ml)来估计参数θ，即求取θ使得对数似然函数值最大

由公式(2)可以看出，k个高斯分量的加权和出现在对数函数内，导致最大似然无法求得闭合形式解，因此需要考虑使用其他方法。期望最大化(expectation-maximization，em)作为求解此种问题的常用方法，核心思想在于通过引入隐藏变量使得不完备观测数据集x变为完备观测数据集{x,c}，并求取θ使得新的对数似然函数值最大

其中cnk＝1表明第n个样本由第k个高斯分量产生，否则cnk＝0。

由公式(3)可以看出，与原似然函数相比，新的似然函数是易于求解的。考虑到隐藏变量无法进行直接观测，em算法在给定参数θ初始估计的情况下，迭代进行e步骤和m步骤直至收敛或达到迭代次数：

e步骤：使用当前参数θ计算cnk∈c的期望(expectation)，也称为责任(responsibility)，即第n个样本由第k个高斯分量产生的概率，记为rnk

m步骤：使用e步骤中计算的rnk替代公式(3)中的隐藏变量cnk，通过最

大化(maximization)期望的完备数据似然函数，对参数θ进行更新。

一方面，观察公式(4)至公式(7)可以看出，em算法的e步骤和m步骤都是易于实现并行计算的。另一方面，可以证明em算法仅能够确保参数θ收敛到局部极大值，需要提供一个较为合理的初值，因此要特别注意参数θ的初始化问题。

本发明技术方案的cuda(computeunifieddevicearchitecture，统一计算设备架构)实现如图2所示，其中实线箭头表示控制流，虚线箭头表示数据流；左侧部分代表cpu中执行的程序，右侧部分代表gpu中执行的程序，橙色矩形代表内核；宿主(host)内存中的变量以h_开头，设备(device)内存中的变量以d_开头。以下分别对五个步骤的实现进行介绍。

1.读入点云数据

点云通过一系列离散点(x,y,z)来对目标进行描述，在存储上采用图3(a)中结构体数组(arrayofstructure，aos)的形式。而gpu(图形处理器)中对全局内存进行合并访问(coalescedaccess)，为提高访问效率，需要访问连续、对齐的内存。因此在读入点云的过程中，将数据以图3(b)中数组结构体(structureofarray，soa)的形式存储在宿主内存，记为h_data；随后将h_data上传至设备内存以供gpu访问，记为d_data。

2.k-means++

为解决k-means算法依赖初始聚类中心、容易陷入局部极值的问题，在进行聚类之前，首先使用k-means++从点云中筛选出k个点作为初始聚类中心，其伪代码如算法1中所示。

从中可以看出，每个新聚类中心的选取都依赖于已选定的聚类中心，不易于并行处理，因此该模块在cpu中进行计算。在筛选每个新聚类中心的过程中都进行三次尝试，并取d最小的那个点作为结果。得到k个初始聚类中心后以结构体数组形式存储于宿主内存，记为h_center；随后将h_center上传至设备内存以供gpu访问，记为d_center。

3.k-means

本发明目的在于使用gmm对点云进行重建，侧重点在于求取gmm的参数θ，对点云的聚类只是附属品。因此使用计算更为简单的k-means对数据进行初步聚类，根据聚类结果计算各个簇的均值μk、协方差矩阵σk和权重πk，将其作为gmm参数θ较为合理的初值。这一方面有助于避免局部极值，加快gmm的收敛速度；另一方面巧妙地绕过了初期gmm计算量较大的m步骤，实现整体上的加速效果。

(1)计算各点到各聚类中心的距离

点云中的n个点到某聚类中心μk的距离可以并行计算，某点xn到k个聚类中心的距离也可以并行计算，因此该部分在gpu中进行计算，其伪代码如内核1中所示。

该内核计算得到存储于设备内存，记为d_dist，可视为k×n的二维矩阵，其第k行、第n列的元素d_dist[k*n+n]代表点xn到聚类中心μk的距离。若每个block中包含的thread数量设置为threadsperblock，则在调用该内核时block的维度设置为dimblock(threadsperblock,1,1)；考虑到n的值远大于k，grid的维度设置为其中符号表示向上取整。每个block中计算不同点到同一聚类中心μk的距离，因此将d_center中对应的数据载入block的共享内存中以提升访问速度。

(2)更新各点分类标记

在更新分类标记的过程中，每个点xn要比较自身到k个聚类中心的距离，而n个点的更新相互独立，可以并行进行，因此该部分在gpu中进行计算，其伪代码如内核2中所示。

该内核计算得到n个点的标记存储于设备内存，记为d_label；随后将d_label下载至宿主内存以供cpu访问，记为h_label，其第n个元素存储点xn的分类标记情况。在调用该内核时，block的维度设置为dimblock(threadsperblock,1,1)，grid的维度设置为

(3)更新各聚类中心

在更新各聚类中心的过程中，需要判断每个点的归属情况，进而对相应的中心进行更新，这会在程序中引入分支结构。gpu的基本执行单元warp必须运行相同程序，在遇到分支结构时会依次执行各个分支，并阻塞不满足分支条件的thread，这极大地降低了程序性能。

为解决这一问题，baih.t等首先将gpu中计算得到的各点分类情况下载至宿主内存，发挥cpu擅长分支预测的优势，在cpu中对数据重新进行排列，使得属于同一簇的点在内存上连续分布；然后将重排后的数据上传至设备内存后，在gpu中更新各聚类中心，每个block中的thread都只处理属于同一簇的点。采用这种策略虽然消除了分支结构的影响，但忽略了数据重排列和大量数据(3×k×n个数据)在宿主/设备内存中传输所带来的时间成本。

本发明更新聚类中心的部分全部在cpu中进行计算，每次只将更新后的聚类中心h_center(3×k个数据)重新上传至设备内存，其伪代码如算法2中所示。

(4)收敛判断

使用两种准则进行收敛性判断，满足任意一个即停止迭代：一是最大迭代次数超过阈值；二是分类标记发生变化的点在点云中所占的比重低于阈值。

4.权重、均值、协方差矩阵初始化

均值可以利用上一步k-means最后一次迭代过程中计算获得的聚类中心进行初始化，而协方差矩阵和权重的初始化与k-means更新聚类中心的过程类似，也在cpu中进行计算，其伪代码如算法3中所示。

初始化完成后将宿主内存中存储权重、均值和协方差矩阵的h_weight、h_mean和h_covar上传至设备内存以供gpu访问，分别记为d_weight、d_mean和d_covar。

5.高斯混合模型(gaussianmixturemodel,gmm)

k-means需要“硬性”地对各点进行分类，并根据各点的归属情况对相应聚类中心进行更新；gmm在e步骤中获取的责任矩阵则“软性”地衡量各点由某个簇产生的概率，在m步骤中每个点对于每个高斯分量权重、均值和协方差矩阵的更新都有所贡献，更加适合并行处理，因此e步骤和m步骤的计算都在gpu中进行。

(1)计算责任矩阵

首先计算各协方差矩阵的行列式其伪代码如内核3中所示。因为矩阵行列式的计算过程并不易于并行，在调用该内核时，使用k个block分别对每个高斯分量的协方差矩阵进行处理，每个block包含一个thread，即block的维度设置为dimblock(1,1,1)，grid的维度设置为dimgrid(k,1,1)。

随后计算责任矩阵，该过程与k-means中计算各点到各聚类中心的过程类似，其伪代码如内核4中所示。在调用该内核时，block的维度设置为dimblock(threadsperblock,1,1)，每个block中计算不同点由某个高斯分量产生的概率；grid的维度设置为

最后对内核4中获取的责任矩阵进行归一化，其伪代码如内核5中所示。在调用该内核时，block的维度设置为dimblock(threadsperblock,1,1)，每个thread计算一行元素之和，随后对该行元素进行初始化；grid的维度设置为作为中间变量的每行元素之和存储于设备内存中以备后面使用。

(2)更新权重

欲对gmm第k个高斯分量的权重πk进行更新，需要求得责任矩阵第k列中所有元素之和这是一个典型的规约(reduction)操作，非常适合于并行计算。nvidiacudasdk中给出了规约的多种实现方式，其中效率最高的一种采用了合并寻址、循环展开等策略。本文的实现即基于这种规约方式，为方便描述将循环展开部分记为函数reduce()，其功能在于将block共享内存中的数据进行加法规约，结果存储在共享内存的首地址中。

因为各block间无法通信，所以权重的更新过程分为两个步骤。第一步使用多个block并行进行规约，将当前block中的规约结果存储于全局内存中的相应位置，其伪代码如内核6中所示。

在调用该内核时，block的维度设置为dimblock(threadsperblock,1,1)；因为在将数据从全局内存读入共享内存的过程中进行了一次规约，并且责任矩阵每一列的规约可以并行，grid的维度设置为规约的初步结果sum可以看做的二维矩阵，其第j行、k列元素sumjk对应blockidx.x为j的block对责任矩阵中第k列2*threadsperblock个元素的规约结果。

第二步对sum中的每一列进行规约，获取责任矩阵每列元素之和以及相应的权重其伪代码如内核7中所示。在调用该内核时，block的维度设置为dimblock(threadsperblock,1,1)，grid的维度设置为dimgrid(k,1,1)。

(3)更新均值

欲对gmm第k个高斯分量的均值μk进行更新，分母项已经在更新权重的过程中获得，只需获取分子项分子项的每个元素rnkxn由标量rnk乘以三维向量xn获得，因此均值的更新过程与权重的更新过程类似，仍旧分两步，但需要在三个不同维度分别进行。其第一步的伪代码如内核8中所示，其中表示点xn的第i∈{1,2,3}维坐标。

在调用该内核时，block的维度设置为dimblock(threadsperblock,1,1)，grid的维度设置为规约的初步结果sum可以看做的三维矩阵，其第j行、k列、i页元素对应索引为(j,k,i)的block对2*threadsperblock个元素的规约结果。

第二步对sum不同页中的每一列进行规约，获取相应的均值其伪代码如内核9中所示。在调用该内核时，block的维度设置为dimblock(threadsperblock,1,1)，grid的维度设置为dimgrid(k,3,1)。

(4)更新协方差矩阵

欲对gmm第k个高斯分量的协方差矩阵σk进行更新，分母项同样已经在更新权重的过程中获得，只需获取分子项分子项的每个元素由标量rnk乘以3×3矩阵(xn-μk)(xn-μk)^t获得，因此协方差矩阵的更新过程亦与权重的更新过程类似，分两步进行。其第一步的伪代码如内核10中所示。

在调用该内核时，block的维度设置为dimblock(threadsperblock,1,1)，grid的维度设置为规约的初步结果sum可以看做的三维矩阵。

第二步对sum不同页中的每一列进行规约，获取相应的协方差矩阵其伪代码如内核11中所示。在调用该内核时，block的维度设置为dimblock(threadsperblock,1,1)，grid的维度设置为dimgrid(k,9,1)。

(5)收敛判断

与k-means类似，gmm也使用两种准则进行收敛性判断，满足任意一个即停止迭代：一是最大迭代次数超过阈值；二是本次迭代中log似然函数值的变化小于阈值。其中已经在归一化责任矩阵的过程中计算得到，只需要对其取对数后进行规约操作。计算log似然函数的伪代码仍旧分两步进行，如内核12、13所示。在调用内核12时，block的维度设置为dimblock(threadsperblock,1,1)，grid的维度设置为在调用内核13时，仅使用一个block，其维度设置为dimblock(threadsperblock,1,1)。

本发明提供的一种基于加速高斯混合模型的空间目标重建方法，针对航天器近距离感知中的目标三维结构重建问题进行研究，目标点云中包含数量巨大的点，可能有几万、几十万乃至上百万的量级，重建过程中计算速度快，鲁棒性强。

以上，仅为本发明示例性的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以权利要求的保护范围为准。