第11讲-点云数据处理20191111.pptx

第11讲-点云数据处理20191111.pptx

《第11讲-点云数据处理20191111.pptx》由会员分享，可在线阅读，更多相关《第11讲-点云数据处理20191111.pptx（92页珍藏版）》请在冰豆网上搜索。

,点云数据处理,中国科学院自动化研究所邹伟2018.12.03,人工智能学院课程：

智能传感与信息处理第11讲,1/92,点云：

在同一空间参考系下表达目标空间分布和目标表面特性的海量点集合。

在获取物体表面每个采样点的空间坐标后，得到的是一个点的集合，称之为点云（PointCloud）。

获取设备：

激光扫描仪、雷达、立体视觉、RGB-D相机等,2/92,根据激光测量原理得到的点云：

包括三维坐标（XYZ）和激光反射强度,点云来源,根据摄影测量原理得到的点云，包括三维坐标（XYZ）、颜色、纹理等信息。

3/92,点云滤波,点云精简,点云配准,点云分割,曲面建模,三维模型构建、环境理解等,位姿估计、空间感知融合与增强等,4/92,目的：

通过机载激光雷达快速获取高精度三维地理数据，对这些点云数据进行滤波过程就是将LIDAR点云数据中的地面点和非地面点进行分离的过程。

数字高程模型（DEM：

DigitalElevationModel）：

地面地形的数字化模拟，是用一组有序数值阵列形式表示地面高程的一种实体地面模型。

2.1面向地形构建的LIDAR点云数据滤波,5/92,机载Lidar+GPS+INS,2.1面向地形构建的LIDAR点云数据波,基于坡度理论的方法,2.1面向地形构建的LIDAR点云数据波,平面拟合的方法,曲面拟合的方法,聚类与分割的方法,7/92,2.1面向地形构建的LIDAR点云数据波,基本原理：

认为地形表面平缓光滑，局部区域内地形发生急剧变化的可能性较小，通过比较两点间的高差值是否满足高差函数来判断点是否为地面点。

=|:

滤波核函数（）：

非递减,=+.,：

最大地形坡度百分比；两点间距离；为标准偏差。

最大距离设为10m,2.1.1基于坡度理论的方法,8/92,=,2.1.1基于坡度理论的方法,2.1面向地形构建的LIDAR点云数据波,9/92,Step1:

采用较大移动窗口，在窗口内搜索最低点，由每次移动窗中找出的最低点的全体，拟合地面的平面方程（地面模型）。

Step2：

将所有点与地面模型比较，计算每点到平面的距离（高度差），凡高度差超过某个阈值，则认为是非地面点，将这些点滤除；Step3：

缩小窗口，对剩余地面点按照第1和第2步骤计算地面模型，改变阈值，进一步滤除与当前地面模型之间的高度差较大的点，重复数次后得到比较精确的地面点集合。

2.1.2平面拟合方法,2.1面向地形构建的LIDAR点云数据波,窗口小，有可能将一些大房屋顶点保留下来；窗口太大，则会将地面过于平滑，微小地形变化部分被滤除。

阈值太大，会将一些植被点作为地面点保留下来；阈值太小，可能将真实的较小的地形突变点去掉。

10/92,2.1.3移动曲面方法,2.1面向地形构建的LIDAR点云数据波,基本原理：

激光点云的空间关系反应了地形表面的空间变化，任何一个复杂的空间曲面，其局部面元可利用一个简单的二次曲面拟合：

=,=+,当局部面元小到一定程度，甚至可以将其表达为一个平面,=,=+,11/92,2.1.3移动曲面方法,2.1面向地形构建的LIDAR点云数据波,Step1:

选择局部最低的3点作为初始种子点,Step2:

初始平面拟合,Step3:

基于平面方程判别邻近激光点,Step4:

融入新的地面点进行平面拟合，当达到6点后进行曲面拟合,Step5：

基于二次曲面方程进行地面点的迭代判别，并不断更新地形曲面Step6：

移动窗口，重复复上述步骤,12/92,2.1.3移动曲面方法,2.1面向地形构建的LIDAR点云数据波,滤波前（11）,滤波后（11）,算法难点：

种子点选择以及滤波阈值的确定。

种子点选择不恰当会使得曲面迭代拟合结果陷入极值，无法得到正确结果；同时滤波阈值需要根据地形起伏自适应变化，否则难以取得较好的效果。

13/92,2.2双边滤波器,双边滤波器（Bilateralfilter）：

用在图像处理领域，可以保边去噪。

同时考虑中心像素点与邻域像素点的距离以及像素亮度差值确定权值进行平滑滤波,=,（,）,（,）,=,=+,（,）,滤波器公式：

14/92,滤波目标值：

深度信息（x,y方向通常决定视场）如何确定两部分信息？

（对应于图像中的像素点灰度差值和像素点距离）,2.2双边滤波器,利用中心扫描点的局部邻域点确定一局部平面S，作为点云表征曲面的局部逼近。

对于邻域区域内的任意一点Pi，Pi相对于平面S的符号距离可以看作是双边过滤时的像素点强度差值Pi与中心扫描点之间距离在S上的投影可认为是双边过滤的像素点距离,15/92,P,n,Pi,=,=,=（）,（）,=,滤波器：

距离计算：

权值计算：

2.2双边滤波器,16/92,知识点补充：

PCA算法,如何理解向量？

+,向量准确描述：

给定一组基，确定向量在基向量方向上的投影。

基变换：

基向量1：

基向量2：

行向量基向量,列向量原始数据向量,=,=,=,17/92,选择不同基可以对同一组数据给出不同表示，而且如果基的数量少于向量本身的维数，则可以达到降维的效果。

N维,K个,K维,关键问题：

如何选择基，才能最大程度上保留原有信息？

减掉均值,知识点补充：

PCA算法,=,18/92,协方差矩阵,协方差矩阵：

m个n维数据，将其按列排成维矩阵X，协方差矩阵确定为,=,维实对称矩阵,思想：

设基转换后的m个矢量组成矩阵Y，如果Y的协方差矩阵为对角阵（除对角线外的其它元素化为0），则表征该m个矢量可以将主要信息都较为独立的投影至各基矢量方向上，可以通过选用更少的基矢量表征这些矢量，达到降维效果。

知识点补充：

PCA算法,19/92,协方差矩阵对角化,由基矢量组成的矩阵P，为将原始协方差矩阵C进行对角化的相似变换矩阵。

任务描述：

寻找矩阵P，满足为对角阵，并且对角元素按从大到小排列，则P的前K行就是要寻找的基，用P的前K行组成的矩阵乘以X就能将X从N维降到K维，并满足优化条件。

根据矩阵理论，实对称矩阵满足不同特征值对应特征向量必然正交对于重数为r的特征值，必有r个线性无关向量对应于，可将其单位正交化处理,知识点补充：

PCA算法,20/92,=,=,满足：

=argmin=argmin,理论分析：

选择一个最佳向量,=,=,由1维恢复为n维,Frobenius范数,=argmin,=argmin,=（）=（）,=argmin+,=argmin,=,=argmin,=（）,=argmax,为最大特征值对应的特征向量,=argmax,知识点补充：

PCA算法,21/92,PCA（主成分分析：

principalcomponentsanalysis）基本步骤：

Step1：

将原始数据按列组成n行m列矩阵X；Step2：

将X的每一行进行零均值化，即减去这一行的均值Step3：

求协方差矩阵=Step4：

求协方差矩阵的特征值及其对应的归一化特征向量Step5：

将特征向量按对应特征值大小从上到下按行排列成矩阵，取前K行组成矩阵P；Step6：

Y=即为降维到K维后的数据,知识点补充：

PCA算法,22/92,经典应用：

人脸识别中的特征脸方法Eigenface,知识点补充：

PCA算法,23/92,=（）,（）,=,滤波器：

权值计算：

问题及关键：

平面S法线n？

-基于PCA的法线估计,=（）,=,协方差矩阵：

描述了邻域点云几何空间分布的统计特性,特征值分解：

C为对称半正定，特征值均为实数,=min,与之对应的特征矢量作为平面法线矢量n,2.2双边滤波器,24/92,原始点云,滤波器权值分布,滤波后点云,2.2双边滤波器,25/92,基本原理：

使用结构元素的窗口模板作为处理单元，利用形态学中的膨胀与腐蚀算法相组合，形成开、闭两种算子进行综合处理。

2.3基于数学形态学的滤波方法,结构元素：

二维结构元素由数值为0或1的矩阵组成，相当于滤波窗口,输入图像,输出图像,行i,列j,列j,行i,输出像素,逻辑运算,在结构元素每个像素位置上与二值图像对应的区域进行特定的逻辑运算,运算效果：

取决于结构元素的形状、大小、内容以及逻辑运算性质,26/92,=:

+,腐蚀运算：

2.3基于数学形态学的滤波方法,与模板B覆盖下的点进行逻辑或,被滤波像素点,图像中的目标区域,腐蚀过程：

拿模板B的原点和输入图像X上的点一个一个对比，如果B上的所有为1的元素所覆盖的点都在目标A的范围内，则B原点所对应的点判断为属于目标A，否则不属于。

A,X,B,27/92,2.3基于数学形态学的滤波方法,腐蚀运算：

=:

+,最终结果,可以看出：

结构元素的原点坐标也很重要！

28/92,=（）,膨胀运算：

2.3基于数学形态学的滤波方法,膨胀过程：

拿模板B的原点和输入图像X上的点一个一个对比，如果B上有一个为1的元素所对应的图像点在目标A的范围内，则B原点所对应的点判断为属于目标A，否则不属于。

边缘不处理,29/92,=（）,膨胀运算：

2.3基于数学形态学的滤波方法,（）,（）,30/92,白,绿,白,绿,2.3基于数学形态学的滤波方法,腐蚀,膨胀,31/92,开运算：

先腐蚀后膨胀,=,闭运算：

先膨胀后腐蚀,=,2.3基于数学形态学的滤波方法,32/92,2.3基于数学形态学的滤波方法,=,=,噪声滤除，尺寸减小,尺寸复原，有间断,尺寸扩大，间断修复,尺寸复原,开运算作用：

去除孤立小点、毛刺和小桥，总的位置和形状不变闭运算作用：

填平小孔，弥合小裂缝，总的位置和形状不变,33/92,点云数据的数学形态学运算：

假设三维点为,，对该点的膨胀腐蚀等操作就是对其值进行运算。

膨胀和腐蚀分别定义为：

=max（,）（）,=min（,）（）,W为邻域窗口，窗口的形状可以是一维、二维方形或其它形状。

膨胀操作是高度提升为邻域内的最高值；腐蚀将高度降低为邻域内的最低值。

膨胀,腐蚀,2.3基于数学形态学的滤波方法,34/92,2.3基于数学形态学的滤波方法,35/92,概念与目的：

原始数据量庞大，数据点密集，直接进行曲面重构或三角网格化，三角面片的数量巨大，增加计算时间；数据过细影响重构曲面的光滑性。

以较少数量的点来逼近原始点云模型的过程就称为点云精简。

算法设计约束与思路：

给定冗余误差阈值计算最小采样数分布；给定目标采样点数据搜索最小误差分布。

算法评估度量：

精度：

依据精简数据重建所得表面与真实表面差别要尽可能小，保持原特征。

简度：

精简后的三维点云数目尽可能少。

速度：

精简过程时间短、速度快，,36/92,根据点云的数据类型分类,散乱点云精简方法,三角网格点云精简方法,多边形点云精简方法,扫描线点云精简方法,包围盒法,均匀网格法,随机采样法,曲率采样法,聚类法,迭代法,边界保留法,等分布密度法,最小包围区域法,间距缩短法,等址缩减法,倍率缩减法,均匀弦长采样法,曲率累加采样法,弦高差重采样法,37/92,原理：

首先把所得到的数据点进行均匀网格划分，然后从每个网格中提取样本点，网格中的其余点将被去除掉。

步骤：

在垂直于扫描方向建立网格平面。

将所有点投影至网格平面上，每个网格与对应的数据点匹配。

每个网格中的点按照点到网格平面的距离远近排序，如果某个点位于各个点的中间，那么这个点就被选中保留。

3.1均匀精简法,38/92,3.2基于三角面片的散乱点云数据精简,基本思想：

将扫描获取的数据点云直接三角网格化，比较数据点所在三角面片的邻近三角面片法向量，根据法向量相似性在平面或近似平面的较平坦区域用大的三角面片取代小的三角面片删除多余点，从而实现数据精简。

39/92,Delaunay三角化准则：

空圆特性：

在Delaunay三角形网中任一三角形的外接圆范围内不会有其它点存在。

最大化最小角特性：

在散点集可能形成的三角化分中，Delaunay三角剖分所形成的三角形的最小角最大。

具体的说是指在两个相邻的三角形构成凸四边形的对角线，在相互交换后，六个内角的最小角不再增大。

3.2基于三角面片的散乱点云数据精简,40/92,三角形单位法矢量计算：

=,平坦曲面上的点构成的三角网包含的信息量最少，即这些三角形的法线向量的夹角越小。

如果以为顶点的三角形有n个，两两之间的法线向量的夹角有个，求这个夹角的最大值，如果该最大值小于设定阈值，则该点予以删除，否则保留。

3.2基于三角面片的散乱点云数据精简,41/92,3.3面向特征保持的点云数据精简,第1大步：

利用均值聚类法在空间域对点云全局聚类。

Step1：

选择初始化聚类中心；Step2：

根据点云到各类中心的距离，每个点选择距离最近类作为类别标号；Step3：

根据新的分类结果，计算同类中点云位置平均值作为新的聚类中心；Step4：

重复步骤2-3，直到类中心的位置不在发生变化。

基于K-D树的聚类中心初始化：

保证点云聚类快速收敛，最好能均匀分布在模型各处。

对点云数据构建K-D树，以K-D树的部分节点作为初始化聚类中心。

42/92,K-D树：

一种分割k维数据空间的数据结构，是二叉搜索树的扩展。

在每个内部节点中，用一个k-1维与坐标轴平行的超平面将节点所表示的k维空间分成两部分，存储在子树中的点大约一半落入一侧，另一半落入另一侧。

当一个节点中的点数少于给定的最大点数时，划分结束。

超平面在k个方向上交替出现，每个超平面中至少包括一个数据点。

3.3面向特征保持的点云数据精简,第1大步：

利用均值聚类法在空间域对点云全局聚类。

43/92,中位数划分,3.3面向特征保持的点云数据精简,第1大步：

利用均值聚类法在空间域对点云全局聚类。

44/92,点的最近邻回溯搜索,第2大步：

用主成分分析法估计点云法矢和候选特征点。

遍历每个聚类，若类中包含特征点则将该类细分为多个子类，细分时将聚类映射到高斯球。

=,=min,与之对应的特征矢量为平面法线矢量n,3.3面向特征保持的点云数据精简,主成分分析PCA：

45/92,最小特征值表示了点沿着曲面法矢的变化情况，可利用各点的最小特征值来估计高曲率区域。

=+,=时，为平面点；越大，曲面该点处越尖锐。

大于设定阈值，则设为候选特征点,3.3面向特征保持的点云数据精简,46/92,平面的高斯图为高斯球上的一点如果一连通曲面的高斯图为高斯球上的一点，则该曲面为平面。

估计的法矢总是存在误差，因此平面的高斯图为高斯球上的密集分布聚类,第3大步：

遍历每个聚类，若含候选特征点（需细分）则该聚类映射到高斯球。

曲面的高斯映射：

曲面每点的法向向量映射到一个单位球上，此球称为高斯球。

=（,）|+=,映射过程：

将法向矢量的起点平移至高斯参考球的球心，每个法向矢量会与高斯参考球的球面有一个交点，这个交点就称为法向矢量在高斯球上的高斯映像。

=，=，=,=,3.3面向特征保持的点云数据精简,47/92,第4大步：

用均值漂移算法对具有特征点的聚类的高斯图进行聚类细分。

高斯球上的各子聚类对应为空间点云的子聚类。

无法事先指定聚类细分数目，需要使用无参聚类算法对高斯图数据进行自动划分。

具有特征点的聚类包括多个平坦区域和特征点，多个平坦区域映射为高斯图上的多个聚类，特征点映射为独立的一个或者多个聚类。

均值漂移算法（Mean-Shift）,高斯图上各聚类中心所对应的点云为局部曲面细分结果,3.3面向特征保持的点云数据精简,48/92,点云,基于点云的三角化模型,最大误差：

原始曲面S上采样点到精简点云曲面S上投影点之间欧式距离的最大值平均误差：

原始曲面S上采样点到精简点云曲面S上投影点之间欧式距离的平均值,3.3面向特征保持的点云数据精简,49/92,50/92,基本原理：

Mean-Shift是一个迭代步骤，即先算出当前点的偏移均值，将该点移动到此偏移均值，然后以此为新的起始点，继续移动，直到满足最终条件。

51/92,是一个半径为的高维球区域（前图中的蓝圆形区域），定义为：

存在的问题：

在的区域内，每一个点对x的贡献是一样的。

实际上每一个样本的重要程度也是不一样的。

基本的MeanShift向量形式：

对于给定的d维空间Rd中的n个样本点xi，i=1,n，则对于x点，其MeanShift向量的基本形式为：

=,=,52/92,改进的MeanShift向量形式：

在基本的MeanShift向量形式中增加核函数和样本权重，得到如下的改进的MeanShift向量形式：

（）=（）=（）,（）=（）,=,（）是每一个样本的权重。

为单位核函数；H是一个正定的对称dd矩阵，称为带宽矩阵，是一个对角阵。

=,（）=（）=（）,MeanShift向量（）是归一化的概率密度梯度,53/92,概率密度梯度,54/92,（）,核密度估计（Parzen窗估计）：

概率论中用来估计未知的密度函数。

对一个概率密度函数f（x），已知d维空间中n个采样点,=,，f（x）的核函数估计为：

概率密度梯度,（）是一个赋给采样点的权重K（x）是一个核函数h为核函数窗口宽度h很大，变化缓慢，估计分辨率低。

h很小，变化迅速，估计统计变化大。

=（）=（）,55/92,概率密度梯度,=（）=（）=（）=（）,=（）+=（）,=（）=（）=（）=（）,56/92,=（）,（）=（）=（）,=（）=（）,=（）,计算（）令+=（）如果（+）（），结束循环，否则，重复上述步骤,MeanShift向量修正：

MeanShift步骤：

meanshift就是沿着密度上升的方向寻找同属一个簇的数据点,57/92,常用核函数：

均匀核：

=,=,Gauss核：

=,Epanechnikov核：

=（）,58/92,59/92,点云配准：

将两个或两个以上坐标系中的三维点云数据转换至统一坐标系统中的数学计算过程。

点云配准的实质就是空间坐标变换。

空间坐标变换可以由三类参数唯一确定：

尺度、旋转和平移。

60/92,配准方法分类,按照配准算法所利用的信息，可分为：

基于特征点的配准方法：

根据靶标或者手工选取特征点进行配准（六参数法）基于整体迭代优化的配准方法：

通过迭代计算两个点集最优匹配的几何变换（ICP算法）按自动化程度可以分为人工、半自动和全自动三种类型。

61/92,六参数法点云配准,点云配准实质是空间坐标变换。

空间坐标变换由三类参数确定：

尺度、旋转和平移。

点云配准前后相对大小没有发生改变，即不发生尺度变化，因此只需求解三个旋转参数和三个平移参数，所以称为六参数配准模型,设=,表示一个自由坐标系下的点云坐标，=,表示点在绝对坐标系下的坐标。

两个点云坐标系坐标轴之间的夹角为、和，平移参数为、和，则坐标变换关系为,=（,）+,至少需要3组点云数据，通过手动或者靶标提供。

求解方程组。

62/92,六参数法点云配准-靶标提供数据,63/92,六参数法点云配准-手动提供数据,64/92,ICP点云配准,迭代最近邻点配准算法（IterativeClosestPointAlgorithm，ICP）：

最为经典的数据配准算法。

其特征在于，通过求取源点云和目标点云之间的对应点对，基于对应点对采用最小二乘构造旋转平移矩阵，并利用所求矩阵，将源点云变换到目标点云的坐标系下，估计变换后源点云与目标点云的误差函数，若误差函数值大于阈值，则迭代进行上述运算直到满足给定的误差要求。

Step1：

寻找目标点云的特征点及其与特征点最邻近的源点云中的点Step2：

根据点之间的对应关系，利用最小二乘估计空间变换关系；Step3：

利用估计出来的空间变换关系对源点云进行作用和变换，然后利用变换后数据进行Step1Step3迭代，直至误差满足要求。

65/92,ICP匹配算法（IterativeClosestPoint）,目标点云,源点云,66/92,ICP匹配算法（IterativeClosestPoint）,67/92,ICP匹配算法（IterativeClosestPoint）,求取两组点云的重心：

=（+）,=,=,=,=,=+,点云旋转平移后二者的重心重合：

问题描述：

给定一组对应点和，求取旋转矩阵和平移矢量，使得下面的目标函数值最小,点云变换：

68/92,ICP匹配算法（IterativeClosestPoint）,=（+）,=+,=,=（）（）,=（+）,最小化函数推导：

minmax=,=,其中：

=max（=）,=max（）,69/92,ICP匹配算法（IterativeClosestPoint）,因此最终结果可确定为：

max（）,=,根据施瓦茨不等式，对于任意正交矩阵B，有：

H进行SVD分解：

=,=,令：

=,（）（）,=,=,70/92,ICP匹配算法（IterativeClosestPoint）,71/92,ICP匹配算法（IterativeClosestPoint）,72/92,ICP匹配算法（IterativeClosestPoint）,ICP优缺点能够实现点云自动化配准全局配准保证配准可靠性收敛效率及结果取决于初始迭代条件，即初始配准参数点在目标点中的最近邻点耗时，导致效率较低,ICP改进措施-大量工作可供参考点云重采样对应点搜索对应点权重选择误匹配点的去除目标函数的选择,ICP改进措施点云重采样：

降低样本数量，提高配准效率措施1：

对点云进行一致性采样：

格网采样措施2：

对点云进行随机采样措施3：

对点云进行特征采样：

对点云按法向量变化、曲率强度或颜色变化等信息进行特征采样，保留点云中的特征信息用于配准,73/92,ICP匹配算法（IterativeClosestPoint）,ICP改进措施对应点搜索：

提高对应点的搜索效率措施1：

建立空间索引提高搜索效率（Kd-Tree等）措施2：

在距离约束上，增加颜色、强度、法向量等约束，提高最近点为同名点的可靠性。

ICP改进措施对应点权重选择：

考虑每对对应点的可靠性和搜索准确性对转换参数估计的影响，提高参数估计的可靠性措施1：

常数权重值措施2：

分配较小的权重值给较大的点到点距离匹配点对措施3：

引入点对的法矢量信息来设置权重,74/92,ICP匹配算法（IterativeClosestPoint）,ICP改进措施误匹配点的去除：

提剔除一定数量的误匹配点对，降低弱匹配点对或误匹配点对对参数估计的影响措施1：

根据点对距离按比例或中误差阈值去除弱对应点对措施2：

去除那些与邻域点对不一致的点对，如两组对应点对之间的距离差大于一定阈值，则认为二者不一致措施3：

去除那些在点云边缘上的点对，边缘不稳定。

ICP改进措施目标函数的选择：

提高迭代收敛稳定性以及效率措施1：

对应点对直接的距离平方和最小措施2：

引入点对的颜色、强度或法向量差异作为约束措施3：

每个源数据点到以其对应点及对应点法向量建立的局部平面的距离平方和最小。

75/92,点云分割目的：

将点云分为多个同质区域，分割在同一区域中的点云具有相似的性质。

点云分割难点：

数据的冗余；点密度的不均匀性；在数据结构表示上缺乏明确统一的点云数据结构;,76/92,基于RANSAC的点云分割：

RANSAC,应用场合：

复杂场景中物体可归结于简单的几何形状，可以用方程来描述的，或者说，可以用有限的参数来描述复杂的物体，如桌面墙面多半是大平面，桌上的罐子应该是圆柱体，长方体的盒子可能是牛奶等等。

RanSaC算法可以很好的将此类物体分割出来。

77/92,基于RANSAC的点云分割：

RANSAC,RANSAC（RandomSampleConsensus，随机一致性采样）采用迭代的方式从一组包含离群（噪声）的被观测数据中估算出数学模型的参数。

输入：

样本数据：

既包含正确数据（inliers，可以被模型描述的数据），也包含异常数据（outliers，偏离正常范围很远、无法适应数学模型的数据），即噪声。

用于解释观测数据的参数化模型，并且在给定一组正确数据的情况下，存在符合这些数据的模型参数计算方法。

示例问题描述：

直线拟合问题。

噪音较大的数据集中提取直线模型（比方说只有20%的数据是符合模型的），最小二乘法就显得力不从心。

78/92,基于RANSAC的点云分割：

RANSAC,基本思想：

设样本集为P，模型阶次为n，集合P的样本数#（P）n。

从P中随机抽取包含n个样本的P的子集S，利用S中样本计