CN103978690B - 一种面向3d打印的物体内部结构优化方法 - Google Patents

Abstract

本发明公开了一种面向3D打印的物体内部结构优化方法，该方法基于Voronoi图和FEM(Finite？Element？Analysis？Method,有限元分析)，通过FEM计算模型在外力下的应力图，并通过Voronoi图结合应力图计算并产生类蜂窝状的内部结构；蜂窝状的结构可以在提供很好的结构强度的同时最小化耗材，不仅可以在3D打印的过程中节省材料，还可以赋予打印物体类蜂窝状的内部结构，承受给定的外力，使得物体更加坚固、结实。

Description

一种面向3D打印的物体内部结构优化方法

技术领域

本发明涉及本发明涉及一种面向3D打印的物体内部结构优化方法，具体涉及一种基于Voronoi图和有限元分析的优化物体强度体积比的方法。

背景技术

3D打印，即快速成型技术的一种，它是一种以数字模型文件为基础，运用粉末状金属或塑料等可粘合材料，通过逐层打印的方式来构造物体的技术。

最近几年，3D打印越来越多的被人们所提及，将数字模型文件打印成实物的能力，使得它在许多行业得到广泛应用。从技术实现上来看，3D打印通常是将打印材料“一层一层”地叠加起来，从而将数字蓝图变为实物。所以，模型的体积直接关系到打印过程中的耗材，大型复杂的模型将会耗费大量的材料。而打印过程的耗材直接决定了打印的成本。针对这个问题，Wang等人在2013年提出了用蒙皮框架结构代替内部材料的方法(WANG,W.,WANG,T.Y.,YANG,Z.,LIU,L.,TONG,X.,TONG,W.,DENG,J.,CHEN,F.,ANDLIU,X.2013.Cost-effectiveprintingof3Dobjectswithskin-framestructures.ACMTransactionsonGraphics(Proc.SIGGRAPHAsia)32,5)，为了减少内部材料，该方法将问题抽象为在物理和几何的条件约束下的框架结构(数目)的优化问题，但是该方法只考虑了自承重的情形，即只考虑所打印物体受重力的影响，并未考虑到物体受外力的情形，所以该方法的结果虽然可以很大程度的节省材料，但是所打印出来的物体不能保证承受足够的外力。另一方面，出于对打印物体的平衡性的考虑，Prevost等人在2013年提出的方法(PR′EVOST,R.,WHITING,E.,LEFEBVRE,S.,ANDSORKINEHORNUNG,O.2013.Makeitstand:balancingshapesfor3Dfabrication.ACMTrans.Graph.32,4(July),81:1–81:10.)，该方法可以使物体保持平衡，同时也能在一定程度上节省材料，该方法将这一问题形式化为包含平衡性与形变程度两项加权和的目标函数，并将问题定义为能量最小优化问题，来求解出使目标函数目标最优的模型，但是该方法仅仅是为了保持物体站立平衡。

上述的方法在减少打印耗材的同时，都在一定程度上改变了打印物体的物理强度，因为物体的用料与物理强度是紧密相连的。考虑到这一点，Zhou等人在2013年将这样一个结构问题看做带约束的优化问题(ZHOU,Q.,PANETTA,J.,ANDZORIN,D.2013.Worst-casestructuralanalysis.ACMTrans.Graph.32,4(July),137:1–137:12)，该文章为我们提供了物体结构分析的方法，该方法可以对于任意的模型，计算分析得出模型中最薄弱易碎的区域。另外，Stava等人在2012年提出了迭代优化的方法(STAVAO.VANEK,J.BENES,B.CARR,N.,ANDMˇE545CH,R.2012.Stressrelief:improvingstructuralstrengthof3Dprintableobjects.ACMTrans.Graph.31,4(July),48:1–48:11),通过在模型内部挖洞，对模型表面增厚和增加内部支撑结构的方式来增强物体的结构强度，同时也达到了节省材料的目的，问题在于该方法会改变物体的表面结构，对模型的外观产生影响。

发明内容

本发明为了解决上述问题，提出了一种面向3D打印的物体内部结构优化方法，该方法基于Voronoi图和FEM(FiniteElementAnalysisMethod,有限元分析)，通过FEM计算模型在外力下的应力图，并通过Voronoi图结合应力图计算并产生类蜂窝状的内部结构；蜂窝状的结构可以在提供很好的结构强度的同时最小化耗材，不仅可以在3D打印的过程中节省材料，还可以赋予打印物体类蜂窝状的内部结构，承受给定的外力，使得物体更加坚固、结实。

为了实现上述目的，本发明采用如下技术方案：

一种面向3D打印的物体内部结构优化方法，包括以下步骤：

(1)计算初始给定的实心模型S的应力分布SM(S，F)，确定目标优化函数的参数值，并将其初始化；

(2)根据计算得到的应力分布，进行重心Voronoi分割，得到α个Voronoi分割区域单元c＝{c₁,…c_i,…,c_α}，计算分割区域单元Harmonic场的等值曲面；

(3)将Voronoi分割区域单元中等值曲面内的部分挖空，建立新的实心模型S′，令S＝S′，重新计算其应力分布；

(4)根据新的应力分布，对比每个Voronoi分割区域单元，如果每个Voronoi单元内部都没有大于屈服点的应力点，转到(7)；如果Voronoi分割区域单元中有大于屈服点的应力点且该分割区域单元Harmonic场的等值面与每轮迭代等值面系数的减少步长的差值小于每轮迭代的等值面系数阈值，则转到步骤(5)；如果Voronoi分割区域单元中有大于屈服点的应力点且该分割区域单元Harmonic场的等值面与每轮迭代等值面系数的减少步长的差值不小于每轮迭代的等值面系数阈值，否则，转到步骤(6)；

(5)将S内的相邻且等值面参数一样的单元合并，产生新的Voronoi分割区域单元，新的单元的等值面系数不变，转到步骤(3)；

(6)减小该分割区域单元Harmonic场的等值面，增加该Voronoi分割区域单元的材料；转到步骤(3)；

(7)根据Voronoi分割区域单元，将各个单元内等值面内的部分挖空，得到新的模型，将该模型输出。

所述步骤(1)中，应力分布SM(S，F)的计算方法，具体包括以下步骤：

(1-1)将实心模型S四面体化；

(1-2)对四面体化后的模型进行有限元分析，得到模型内部的应力分布。

所述步骤(2)中，重心Voronoi分割的具体方法，包括以下步骤：

(2-1)在S内部使用误差扩散的方法(errordiffusion)来生成α个初始站点，这样在应力大的区域分布的站点数多于应力小的地方，使得站点的分布与应力分布相匹配，具体为：

(2-1-a)在S内部随机分布一个点p；

(2-1-b)定义点p关于p点处应力值的接受概率函数，该函数与p点处的应力值成正比，并通过随机数生成来判定点p是否被接受；

(2-1-c)重复(a)(b)过程，直至有α个点被接受，即生成了α个初始站点。

(2-2)根据公式 $\mathrm{Vor}\left(c_i\right)=\{x\∈S|\left|\right|x-c_i\left|\right|\≤\left|\right|x-c_j\left|\right|,\∀j\≠i\},$ 用Lloyd方法计算重心Voronoi分割，其中，i＝1,…,α；j＝1,…,α；χ为特定材料的屈服点。

所述步骤(2)中，计算分割区域单元Harmonic场的等值面的具体方法，包括以下步骤：

(2-a)对于每个Voronoi分割区域单元c_i∈c，计算c_i内部的Harmonic距离场，其中i＝1,…,α，；

(2-b)使用MarchingCubes方法构造出对应β_i的等值曲面，其中，β_i∈β；β＝{β₁,…,β_α}表示每个Voronoi分割区域单元中Harmonic场的等值曲面，产生三角形网格曲面；

(2-c)将新构造的α个三角形网格曲面与初始模型S合并，即得到每个Voronoi分割区域单元被挖空一部分的新的模型。

所述步骤(4)的具体方法为：对于每个Voronoi单元c_i∈c，如果每个单元内部都没有大于屈服点的应力点，则转到(7)；如果c_i中有大于屈服点的应力点且β_i-∈＜β_thre，则转到步骤(5)；如果c_i中有大于屈服点的应力点且β_i-∈≥β_thre，转到步骤(6)；其中，∈为每轮迭代等值面系数的减少步长；β_thre为每轮迭代的等值面系数阈值，β_i∈β；β＝{β₁,…,β_α}表示每个Voronoi分割区域单元中Harmonic场的等值曲面。

所述步骤(6)中，具体计算方法包括以下步骤：

(6-1)选出所有内部存在大于屈服点的应力点的Voronoi分割区域单元集合c_y，进行判断；

(6-2)对于所有的c_i∈c_y，将这些Voronoi分割区域单元的等值面减小为β_i-∈；

(6-3)根据新的β，将Voronoi分割区域单元中等值曲面内的部分挖空，建立新的实心模型S′，令S＝S′；转至步骤(3)。

所述步骤(6-1)的具体方法为：选出所有的含有大于屈服点的应力点的Voronoi分割区域单元集合c_y＝{c_i|SM(c_i,F)＞χ,i＝1,…,α}，如果集合c_y为空，则转到步骤(7)；如果集合c_y非空，若β_i-∈＜β_thre，则转到步骤(5)，否则，转到(6-2)；其中，∈为每轮迭代等值面系数的减少步长；β_thre为每轮迭代的等值面系数阈值，并且随着迭代次数的增加该阈值呈线性减小。

本发明的有益效果为：在解决挖空实心3D打印模型，节省材料的问题的同时，提出新的优化算法，在消除施加在物体上的应力的同时，最大化实体模型的物理强度体积比。

附图说明

图1为将模型四面体化后的切面图；

图2为对于模型施加外力后初始应力分布；

图3为包含50个初始站点的分布图；

图4为包含50个站点的重心Voronoi图；

图5(a)为一个单元内部的Harmonic距离场分布图；

图5(b)为根据Harmonic等值面距离场挖空一个单元内部的结果示意图；

图6(a)为优化迭代开始前初始化α＝50,β＝0.9的应力分布；

图6(b)为第一轮迭代结束后的结果示意图；

图6(c)为第二轮迭代开始前对相邻且等值面参数相同的单元进行合并之后的结果示意图；

图6(d)为最终优化后的结果示意图；

图7为本发明的流程框架图。

其中，201所指箭头表示施加在模型上的外力，202所指处的区域分布有大量大于屈服点的应力点，203所指箭头表示物体所受的重力。

具体实施方式：

下面结合附图与实施例对本发明作进一步说明。

如图7所示，一种基于Voronoi图和FEM的优化物体内部结构的3D打印技术，对于一个初始给定的实心模型S，不同的外力集合F，这里我们假设模型在实心的情况下是可以承受这样的外力的，即不会在实心模型内部出现应力大于屈服点的区域；否则的话，实心模型也无法承受这样大的外力，也就没有求解优化的必要性。

接下来，在保证模型可以承受外力集合F的情况下，优化模型的体积W_S，目标优化函数为

$\underset{\mathrm{\α},\mathrm{\β}}{\min }{W}_S(\mathrm{\α},\mathrm{\β})$ s.t.SM(S,F)＜χ

其中α为Voronoi图的站点数，β＝{β₁,…,β_α}表示每个Voronoi单元中Harmonic场的β_i,(i＝1,…,α)等值面，SM(S,F)为将F施加于S上后计算得到的应力分布，χ为特定材料的屈服点(例如，对PC/ABS等树脂材料，屈服点缺省值为4.1e⁷N/m²)。

屈服点即钢材或试样在拉伸时，当应力超过弹性极限，即使应力不再增加，而钢材或试样仍继续发生明显的塑性变形，称此现象为屈服，而产生屈服现象时的最小应力值即为屈服点。

有限元分析(FEM，FiniteElementAnalysisMethod)利用数学近似的方法对真实物理系统(几何和载荷工况)进行模拟。还利用简单而又相互作用的元素，即单元，就可以用有限数量的未知量去逼近无限未知量的真实系统。

Marchingcubes是Lorensen等人与1987年提出的一种体素级的等值面绘制方法(Lorensen,W.E.；Cline,HarveyE.(1987)."Marchingcubes:Ahighresolution3dsurfaceconstructionalgorithm".ACMComputerGraphics21(4):163–169)。在Marchingcubes方法中，对每个Voronoi单元进行体素化，再对每个体素进行处理，以三角面片来逼近其内部的等值面，在对每个体素的“扫描”过程中，构造出等值面。

Harmonic场即利用调和函数(HarmonicFunction)在每个Voronoi单元内的三维空间生成的距离场，距离场中的等值面隐式代表了孔洞的曲面；

本发明通过求解使得目标函数最优的α和β，确定最终的内部结构：

步骤一：计算初始的应力分布SM(S,F)；

步骤二：初始化α＝α⁰，β＝β⁰；

步骤三：根据计算得到的SM(S,F)，在S内部初始化α个种子站点，计算重心Voronoi分割，得到α个Voronoi分割区域单元c＝{c₁,…,c_α}；

步骤四：根据c_i∈c中的β_i∈β等值曲面，将c_i单元中β_i等值面内的部分挖空，得到新的模型S＝S’；

步骤五：重新计算应力分布SM(S,F)；

步骤六：根据新的应力分布，对比每个Voronoi分割区域单元，如果每个Voronoi单元内部都没有大于屈服点的应力点，转到(九)；如果Voronoi分割区域单元中有大于屈服点的应力点且该分割区域单元Harmonic场的等值面与每轮迭代等值面系数的减少步长的差值小于每轮迭代的等值面系数阈值，则转到步骤(七)；如果Voronoi分割区域单元中有大于屈服点的应力点且该分割区域单元Harmonic场的等值面与每轮迭代等值面系数的减少步长的差值不小于每轮迭代的等值面系数阈值，否则，转到步骤(八)；

步骤七：将S内的相邻且等值面参数一样的单元合并，产生新的Voronoi分割区域单元，新的单元的等值面系数不变，转到步骤(四)；

步骤八：减小该分割区域单元Harmonic场的等值面，增加该Voronoi分割区域单元的材料；转到步骤(四)；

步骤九：根据分割区域c＝{c₁,…,c_α},β＝{β₁,…,β_α}，将各个单元内等值面内的部分挖空，得到新的模型，将该模型输出；结束。

所述步骤一、五中涉及应力分布计算的具体步骤为：

(1-1)将输入模型S四面体化；

为了保证下一步有限元分析的准确性，一般情况我们控制四面体的数量在1×10⁴个左右。当然对于不同的模型，由于模型的复杂程度不同四面体的数量可以更多。图1是将模型四面体化后的结果；

(1-2)对四面体化后的模型进行有限元分析，输出应力分布；

图2是对于模型施加外力后初始应力分布。其中201所指箭头表示施加在模型上的外力，202所指处的区域分布有大量大于屈服点的应力点，203所指箭头表示物体所受的重力。

所述步骤三的具体步骤为：

(3-1)在S内部使用误差扩散的方法(errordiffusion)来生成α个初始站点，这样在应力大的区域分布的站点数多于应力小的地方，使得站点的分布与应力分布相匹配，具体做法如下：

(3-1-a)在S内部随机分布一个点p；

(3-1-b)定义点p关于p点处应力值的接受概率函数，该函数与p点处的应力值成正比，并通过随机数生成来判定点p是否被接受；

(3-1-c)重复(a)(b)过程，直至有α个点被接受，即生成了α个初始站点。

如图3，模型内部的点是根据图2所示的应力分布用error-diffusion的方法初始化50个点后的结果。显然在应力分布中图202所指的区域分布的点的密度最大，其余的区域则只有少量的点。这样就使得模型内部应力大的地方可以分布更多的Voronoi单元，将这些单元挖空之后将会产生更多的支撑结构来减轻应力；

(3-2)根据公式 $\mathrm{Vor}\left(c_i\right)=\{x\∈S|\left|\right|x-c_i\left|\right|\≤\left|\right|x-c_j\left|\right|,\∀j\≠i\},$ 用Lloyd方法(LLOYD,S.P.1982.LeastsquaresquantizationinPCM.IEEETransactionsonInformationTheory28,2,129–136.)计算重心Voronoi分割Lloyd。

Lloyd方法是一种经典的求解重心Voronoi图的方法，该方法是一个迭代的确定性算法，Lloyd方法在每次的迭代的过程对于每个站点c_i，计算其Voronoi单元的重心c′_i，令c′_i代替原来的c_i；如图4所示，应力大的区域分布着小且多的重心Voronoi单元。

所述步骤四的具体步骤为：

(4-1)对于每个Voronoi单元c_i∈c，i＝1,…,α，计算c_i内部的Harmonic距离场，距离场中的等值面隐式代表了孔洞的曲面；

图5(a)描述的是一个单元内部的Harmonic距离场的分布情况，由图5(b)挖空单元内部的结果来看，使用Harmonic距离场产生的结构类似于拱形结构，而拱形结构能将局部的应力均匀分布到全局表面上，从而减轻局部应力，并且这种拱形结构还有省材的特点。

(4-2)使用MarchingCubes方法构造出对应β_i的等值曲面，产生三角形网格曲面；

(4-3)将新构造的α个三角形网格曲面与初始模型S合并，即得到每个Voronoi单元被挖空一部分的新的模型。

所述步骤(六)的具体方法为：对于每个Voronoi单元c_i∈c，如果每个单元内部都没有大于屈服点的应力点，则转到(九)；如果c_i中有大于屈服点的应力点且β_i-∈＜β_thre，则转到步骤(七)；如果c_i中有大于屈服点的应力点且β_i-∈≥β_thre，转到步骤(八)；其中，∈为每轮迭代等值面系数的减少步长；β_thre为每轮迭代的等值面系数阈值，β_i∈β；β＝{β₁,…,β_α}表示每个Voronoi分割区域单元中Harmonic场的等值曲面。

所述步骤八具体步骤为：

(8-1)选出所有内部存在大于屈服点的应力点的Voronoi分割区域单元集合c_y，进行判断；

(8-2)对于所有的c_i∈c_y，将这些Voronoi分割区域单元的等值面减小为β_i-∈；

(8-3)根据新的β，将Voronoi分割区域单元中等值曲面内的部分挖空，建立新的实心模型S′，令S＝S′，；转至步骤(四)。

所述步骤(8-1)的具体方法为：选出所有的含有大于屈服点的应力点的Voronoi分割区域单元集合c_y＝{c_i|SM(c_i,F)＞χ,i＝1,…,α}，如果集合c_y为空，则转到步骤(九)；如果集合c_y非空，若β_i-∈＜β_thre，则转到步骤(七)，否则，转到(8-2)；其中，∈为每轮迭代等值面系数的减少步长；β_thre为每轮迭代的等值面系数阈值，并且随着迭代次数的增加该阈值呈线性减小。

将这些单元的等值面减小，就减少了这些单元内部被挖去的材料，也就是增加了单元内部的支撑材料，本方法就是以迭代的方式逐渐减轻应力至屈服点以下。当然，极有可能地，仅一次迭代增加的材料不能够使得该单元内部的应力点降到屈服点以下，所以下一次的迭代会再次选中该单元，并再次增加材料，直到该单元内部不存在大于屈服点的应力点。

同时，在求解过程中，存在β_thre这样一个阈值参数，目的在于控制每轮迭代的等值面的值，如果无法再阈值之内很快的减轻模型应力，则表示局部的改变对整体模型的应力改变影响很小，则转入步骤七，通过改变整体的模型构造，来快速的增强模型的强度，加快算法的求解过程。

图6中，经过两轮的迭代优化，(d)中即为最终优化的结果，可以看到模型应力大的区域大部分单元已经是被填充至实心，从应力分布上来看也已经不存在大于屈服点的应力点。至此，优化结束，将最终模型(d)输出。

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。

Claims (8)

1.一种面向3D打印的物体内部结构优化方法，其特征是：包括以下步骤：

(1)计算初始给定的实心模型S的应力分布SM(S，F)，确定目标优化函数的参数值，并将其初始化；

(2)根据计算得到的应力分布，进行重心Voronoi分割，得到α个Voronoi分割区域单元c＝{c₁,…c_i,…,c_α}，计算分割区域单元Harmonic场的等值曲面；

(3)将Voronoi分割区域单元中等值曲面内的部分挖空，建立新的实心模型S′，令S＝S′，重新计算其应力分布；

(4)根据新的应力分布，对比每个Voronoi分割区域单元，如果每个Voronoi单元内部都没有大于屈服点的应力点，转到步骤(7)；如果Voronoi分割区域单元中有大于屈服点的应力点且该分割区域单元Harmonic场的等值面与每轮迭代等值面系数的减少步长的差值小于每轮迭代的等值面系数阈值，则转到步骤(5)；如果Voronoi分割区域单元中有大于屈服点的应力点且该分割区域单元Harmonic场的等值面与每轮迭代等值面系数的减少步长的差值不小于每轮迭代的等值面系数阈值，则转到步骤(6)；

(5)将S内的相邻且等值面参数一样的单元合并，产生新的Voronoi分割区域单元，新的单元的等值面系数不变，转到步骤(3)；

(6)减小该分割区域单元Harmonic场的等值面，增加该Voronoi分割区域单元的材料，转到步骤(3)；

(7)根据Voronoi分割区域单元，将各个单元内等值面内的部分挖空，得到新的模型，将该模型输出，进行3D打印。

2.如权利要求1所述的一种面向3D打印的物体内部结构优化方法，其特征是：所述步骤(1)中，应力分布SM(S，F)的计算方法，具体包括以下步骤：

(1-1)将实心模型S四面体化；

(1-2)对四面体化后的模型进行有限元分析，得到模型内部的应力分布。

3.如权利要求1所述的一种面向3D打印的物体内部结构优化方法，其特征是：所述步骤(2)中，重心Voronoi分割的具体方法，包括以下步骤：

(2-1)在S内部使用误差扩散的方法(errordiffusion)来生成α个初始站点，这样在应力大的区域分布的站点数多于应力小的地方，使得站点的分布与应力分布相匹配；

(2-2)根据公式用Lloyd方法计算重心Voronoi分割，其中，i＝1,…,α；j＝1,…,α；χ为特定材料的屈服点。

4.如权利要求3所述的一种面向3D打印的物体内部结构优化方法，其特征是：所述步骤(2-1)具体包括以下步骤：

(2-1-a)在S内部随机分布一个点p；

(2-1-b)定义点p关于p点处应力值的接受概率函数，该函数与p点处的应力值成正比，并通过随机数生成来判定点p是否被接受；

(2-1-c)重复(a)(b)过程，直至有α个点被接受，即生成了α个初始站点。

5.如权利要求1所述的一种面向3D打印的物体内部结构优化方法，其特征是：所述步骤(2)中，计算分割区域单元Harmonic场的等值面的具体方法，包括以下步骤：

(2-a)对于每个Voronoi分割区域单元c_i∈c，计算c_i内部的Harmonic距离场，其中i＝1,…,α；

(2-b)使用MarchingCubes方法构造出对应β_i的等值曲面，其中，β_i∈β；β＝{β₁,…,β_α}表示每个Voronoi分割区域单元中Harmonic场的等值曲面，产生三角形网格曲面；

(2-c)将新构造的α个三角形网格曲面与初始模型S合并，即得到每个Voronoi分割区域单元被挖空一部分的新的模型。

6.如权利要求1所述的一种面向3D打印的物体内部结构优化方法，其特征是：所述步骤(4)的具体方法为：对于每个Voronoi单元c_i∈c，如果每个单元内部都没有大于屈服点的应力点，则转到步骤(7)；如果c_i中有大于屈服点的应力点且β_i-∈<β_thre，则转到步骤(5)；如果c_i中有大于屈服点的应力点且β_i-∈≥β_thre，转到步骤(6)，其中，∈为每轮迭代等值面系数的减少步长，β_thre为每轮迭代的等值面系数阈值，β_i∈β，β＝{β₁,…,β_α}表示每个Voronoi分割区域单元中Harmonic场的等值曲面。

7.如权利要求1所述的一种面向3D打印的物体内部结构优化方法，其特征是：所述步骤(6)中，具体计算方法包括以下步骤：

(6-1)选出所有内部存在大于屈服点的应力点的Voronoi分割区域单元集合c_y，进行判断；

(6-2)对于所有的c_i∈c_y，将这些Voronoi分割区域单元的等值面减小为β_i-∈；

(6-3)根据新的β，将Voronoi分割区域单元中等值曲面内的部分挖空，建立新的实心模型S′，令S＝S′，转至步骤(3)。

8.如权利要求7所述的一种面向3D打印的物体内部结构优化方法，其特征是：所述步骤(6-1)的具体方法为：选出所有的含有大于屈服点的应力点的Voronoi分割区域单元集合c_y＝{c_i|SM(c_i,F)>χ,i＝1,…,α}，如果集合c_y为空，则转到步骤(7)；如果集合c_y非空，若β_i-∈<β_thre，则转到步骤(5)，否则，转到(6-2)；其中，∈为每轮迭代等值面系数的减少步长；β_thre为每轮迭代的等值面系数阈值，并且随着迭代次数的增加该阈值呈线性减小。