区域侦察是军用或民用无人机的主要功能之一。利用无人机进行区域侦察以获取地图信息并根据地图进行路径规划是当今一大研究趋势,随着侦察区域面积的增大,单架无人机有限的续航能力便成了束缚侦察能力的一大桎梏[1-2]。因此,多无人机协同侦察开始进入人们的视野。将区域进行分割处理,再将由分割得到的子区域预分配给各无人机,保证各无人机可以在自身的续航能力下协同完成区域的全覆盖侦察是当下热门的研究方案[3-4]。
子区域划分方法主要有Voronoi图法[5]、微元积分法[6]、二分位移迭代法[7]、中点二分迭代法[8],和多边形随机分割算法[9]等。其中,Voronoi图法的复杂度会随着样本点的增多而导致样本点定位时间变长,并且实现细节复杂[10];微元积分法虽然可用于栅格图形的连续性分割,但是算法精度较低[6],生成的子区域并不能保证较好的聚集性;二分位移和中点二分2种迭代算法分别采用不同的分割方法,前者算法为实现均分需要不断地调整分割线,并且当被分割区域呈凹多边形时,分割区域面积的计算成为一个难题;后者算法属于一种平行线分割算法[11],算法需要在初始化时确定平行线方向,不同方向产生的结果并不精确和稳定;多边形随机分割算法适用于对呈规则多边形的农作区域进行分割规划,算法中采用的图斑模型[12]会随着实际区域边界的变化发生改变,因此不适合应用于栅格化后区域边界相对固定的区域均分。
针对未知区域形状呈不规则多边形特点,综合考虑了区域的不规则性,以及特殊多边形可能存在的特殊情况(如凹多边形的凹点等),提出了适用于栅格地图区域渐进均分算法,增强算法对栅格地图的适用性,并结合实际无人机对未知区域进行侦察任务时的工作模式(边界→中心→边界),等分不规则多边形区域栅格地图。
综上,未知区域均分算法包括2个方面,首先建立栅格搜索的改进算法,确定全覆盖侦察区域,其次,建立从区域边界向内逐步划分的均分算法,获得相互邻接的子区域。仿真实验证明,对不同形状的侦察区域进行均分时,子区域划分数不同,均分算法始终保证子区域划分的均匀性、聚集性和连续性,实现未知区域的稳定均分。
考虑到无人机本身所携带设备的局限性,无人机在获取地面信息时,应做到按一定周期将信息发送回通信基站,这便需要无人机在任务过程中与通信基站的间距保持在一定范围内,如图1所示。然而,平行线等均分方法不能较好的保证上述条件。为此,提出了一种由外向内的区域均分方法,利用栅格法构建栅格地图,建立区域渐进均分算法,实现区域均分。
图1 多无人机同步区域侦察示意图
Fig.1 The multi-UAV synchronized area reconnaissance schematic
区域渐进均分算法包含4个阶段。阶段1,迭代法确定侦察区域的边界栅格,构造边界栅格中心数组。阶段2,采用扫描法和双特征标识改进的射线法,赋予栅格特征标识,实现对区域内部栅格的搜索,确定侦察区域内部栅格中心数组。阶段3,模仿水波扩散的子区域栅格的搜索算法,形成初始分割的子区域。阶段4,建立补偿规则,进行子区域优化。栅格地图划分后,形成所有子区域的栅格中心数组,供无人机侦察使用。算法整体流程如图2所示。
图2 区域渐进均分划分算法流程图
Fig.2 Flowchart of the progressive region equalization algorithm
图2中,扫描区域rS和侦察区域rR的关系如图3所示。图3中白色区域为rR,扫描区域rS为rR和黑色区域构成的整体区域。区域渐进算法涉及的数组和相关符号如表1所示。
表1 符号及其含义
Table 1 Symbol and its meaning
名称意义rS扫描区域rR侦察区域l无人机扫描宽度nA无人机扫描的完整栅格数y+表示沿y轴正方向Δk跳跃迭代法的迭代间隔Pb侦察区域边界栅格中心数组PR侦察区域栅格中心数组nR数组Pb包含的元素数Mav预均分栅格数数组Mnv初步优化得到的各子区域栅格数数组Mdv各子区域应调整的栅格数数组In_arrdv子区域栅格中心坐标数组,其中In_arrdv(i)表示第i个子区域包含的栅格中心坐标数组
图3 扫描区域和侦察区域关系示意
Fig.3 The relation between reconnaissance region and scanning region
跳跃迭代法的基本思想为:以间隔Δk在侦察区域的每条边上迭代取点,通过确定迭代点所属栅格以获得区域边界经过的栅格中心数组。当与上一迭代点所属相同栅格时,当前迭代点需要进行“跳跃”以避免后续同一栅格内迭代点的冗余迭代。迭代点跳跃原理如图4所示。
图4 迭代点跳跃原理图
Fig.4 The iteration point jumping schematic
首先,沿迭代终点向x和y方向做栅格平行线,记录2条线与栅格最近的交点;接着,确定区域边界上与两交点分别对应的参考点Pin_x和Pin_y;最后,对比参考点与迭代起点的水平间距,取间距较小的参考点作为跳跃目标点,从而完成跳跃过程。
侦察区域rR的边界线按照上述原理进行迭代和跳跃,不重复记录栅格中心,最终确定边界栅格中心数组为
Pb=[pb_1,pb_2,…,pb_i,…,pb_n]
边界栅格中心数组元素pb_i=(xb_i,yb_i,m,t)T为边界栅格数组中第i个栅格的标识向量,其中(xb_i,yb_i)为栅格中心坐标;m表示迭代过程中的第m条边,作为第一特征标识;t为第二特征标识。
rR内部栅格确定的相关算法主要有射线法[13]、面积法[14]以及夹角法[15]等。其中,射线法相较于其他算法,运算复杂度低并且结果稳定。因此,采用扫描法遍历扫描区域rS栅格的同时,使用射线法作为栅格位置属性的判定算法。
如图5所示,不同于凸多边形,射线法应用于凹多边形时,会因为凹点的存在而导致算法结果不准确。针对射线法的特殊情况(射线经过多边形顶点或与多边形一边重合[16],如图5(b)和图5(c)所示),结合栅格特点,提出了改进的射线法。
图5 射线法示意图
Fig.5 The schematic diagram of the ray method
2.2.1 数组预筛选
当多边形存在凹点时,凹点所对应的内角角度必然大于180°。通过比较两邻边向量的向量积可以初步确定区域多边形凹点数,如式(1)所示:
(1)
式(1)中,Pv_i(xv_i,yv_i)为区域rR按顺时针排列的顶点中第i个顶点。若φ>0,则顶点i为凹点,否则为凸点。改进射线法采用的射线方向为y+方向,所以需要从数组Pb中排除凹点,条件为
xv_i≤min{xv_i-1,xv_i+1}
(2)
或
xv_i≥max{xv_i-1,xv_i+1}
(3)
当侦察区域中一边角系数较大(或较小)时,边界栅格会存在与凹点同列(行)的相邻栅格,如图6所示。因此,在确定满足上述条件的凹点后,算法从Pb中筛选并排除与此凹点栅格同列(行)的相邻栅格(图6中蓝色栅格)。需要注意的是,当出现图中红线所示情况(蓝色栅格同时被另一条边经过)时,则保留此栅格,即保留栅格A。
图6 与凹点所在栅格相邻的栅格特征示意图
Fig.6 The schematic of the raster features adjacent to the raster where the concave point is located
2.2.2 确定区域顶点栅格标识
筛选后的数组Pb,需要对留下的顶点进行特征标识,如图7所示。图7(a)和图7(b)在被射线经过的顶点处有相同的边界栅格,其中以顶点A为端点的邻边l1和l2位于射线的两侧;以顶点B为端点的邻边l1和l2位于射线的一侧。对于与顶点B特征相同的顶点,赋值顶点所属栅格第二特征t标识为“1”,否则栅格第二特征标识t赋值为“0”。
图7 射线经过不同特征顶点示意图
Fig.7 The schematic diagram of the different characteristic vertices through which the rays pass
结合以上2点,侦察区域栅格搜索的改进射线算法流程如图8所示。
图8 改进射线法流程图
Fig.8 The flow chart of improved ray method
rS内的栅格有3种标识α、β和γ(对应栅格向量分别表示为(xS_m,yS_m,α/β/γ)T),分别表征当前栅格属于rR内部、rR边界(即栅格中心属于数组Pb)和rR外部。算法具体过程为:
1) 假设rS内所有栅格均属于rR,修改rS内所有栅格标识为α;
2) 根据数组Pb,将同属于Pb的栅格中心标识修改为β;
3) rS顶行栅格作为起始扫描线,修改起始扫描线中标识为α的栅格,修改为γ标识;
4) 开始扫描,通过比对扫描线和待扫描区域的标识来确定栅格的最终标识,标识比对和修改规则如表2所示。在比对过程中,当扫描线栅格和被扫描栅格标识分别为β和α时,对比结果会出现2种情况(如表2中第3行所示)时,需要进一步采用双特征标识射线法进行二次确定,即步骤5)。
表2 栅格标识比对和修改判定表
Table 2 Raster mark comparison and modification decision
扫描线栅格标识扫描前栅格标识扫描后栅格标识ααααβββαα或γβββγαγγββ
5) 二次确定栅格标识
双特征标识改进射线法通过识别由步骤4)给出的栅格标识来确定射线穿越边界的次数,算法中栅格的第二特征标识的识别优先级高于第一特征标识,若第二特征标识值为“1”,则跳过对当前连续栅格段第一标识的识别过程。
双特征标识改进射线算法的伪代码为:
Input: PN:the grid center point that will be scanned currently; Pcro:the boundary grid array in Pb that is penetrated by the ray emitted by PN
Output: PN(3):the marking value of the point to be scanned(PN(1)=xPN,PN(2)=yPN)
1: Pcro ← sortrows(Pcro,2,“descend” )
2: line_num ← 1
3: diff_arr ←[]
4: For i=1:length(:,1)-1 Do
5: If Pcro(i,2)-Pcro(i+1,2)≠1 Then
6: diff_arr{line_num,:} ← Pcro(i,:)
7: line_num ← line_num+1
8: End If
9: End For
10: diff_arr{line_num,:} ← Pcro(i+1,:)
11: cross_all ← 0
12: cross_num ← 1
13: For i=1:line_num Do
14: If length(find(diff_arr{line_num}(:,4)==1))≠0
Then
15: cross_all ← cross_all + 2
16: continue
17: End If
18: For j=1:length(diff_arr{i}(:,1))-1 Do
19: If diff_arr {i}(j,3) ≠ diff_arr{i}(j+1,3) Then
20: cross_num ← cross_num + 1
21: End If
22: End For
23: End For
24: cross_all ← cross_all + cross_num
25: If cross_all %2==1 Then
26: PN(3) ← α
27: Else
28: PN(3) ← γ
29: End If
1—9行计算射线经过的边界栅格分段数;10—16行判断每段边界栅格中每个栅格的第二特征标识值;17—22行统计每段边界栅格中第一特征标识值不同值的次数;23—28行通过射线穿过区域边界次数的奇偶性确定当前被扫描栅格标识值。
经过算法搜索,最终确定属于区域rR的栅格中心数组为
PR=[pR_1,pR_2,..pR_m,…,pR_nD]
其中,pR_m=(xR_m,yR_m, α/β)T。
区域渐进均分算法的后续阶段以此数组作为基础数据进行侦察区域栅格的均分。
首先,计算子区域划分数,即需要的无人机数量nsp由式(4)计算:
(4)
接着,随机选取rR边界栅格上等栅格数的nsp个栅格作为开始划分各子区域的起点栅格。结合栅格的规则排布特点和栅格中心等距的特点,设计了一种模仿水波的扩散算法。算法原理是模仿水波扩散,通过逐步增加扩散半径来实现子区域的划分。子区域的划分流程如图9所示。
图9 邻边扩散法流程图
Fig.9 The neighborhood diffusion method flowchart
统计每次扩散后各子区域已搜索的栅格数。根据各子区域栅格数,确定下次扩散的优先级:当前总栅格数越少的子区域在下次扩散时具有越高的扩散优先级。算法搜索效果如图10所示,算法随机等栅格数选取nsp(图中nsp=5)个栅格作为起点栅格(图10中红色点)。
图10 扩散算法搜索效果图
Fig.10 The effect map of diffusion algorithm search
在确定的子区域中可能存在栅格数超出nA的子区域,因此需要进一步优化。在算法的第4阶段采用邻边补偿法优化子区域栅格数,如2.4节所示。
设划分给nsp架无人机的均分栅格数的数组为
Mav=[mav_1,mav_2,…,mav_k,…,mav_nsp]
邻边扩散法确定的初始各子区域栅格数的数组为
Mnv=[mnv_1,mnv_2,…,mnv_k,…,mnv_nsp]
则确定各子区域应调整的栅格数的数组为
Mdv=[mdv_1,mdv_2,…,mdv_k,…,mdv_nsp]
其中,mdv_k=mav_k-mnv_k
根据数组Mdv和算法补偿规则表(表3)来确定当前子区域需要进行调整的边界栅格,从而实现当前子区域包含的栅格数与数组Mdv中对应的元素值相等,实现各自区域栅格数的均等化。
表3 算法补偿规则
Table 3 Algorithm compensation rule
第1影响因素第2影响因素第3影响因素需要调整的子区域边界mdv_i>0Δxst_i≥Δyst_iΔxst_i<Δyst_ixst_i>xst_i+1min(In_ardv(i)(:,1))xst_i
表3中,Δxst_i和Δyst_i分别为第i个子区域起点栅格中心的横坐标和纵坐标。
Δxst_i=xst_i-xst_i+1,Δyst_i=yst_i-yst_i+1
由表3可知,对于任一子区域,所要调整的邻边取决于3个影响因素,分别为数组Mdv元素的正负性、相邻子区域起点栅格横纵坐标的差值以及相邻子区域起点栅格的相对位置(注:若mdv_k=0,则跳过本次子区域的邻边补偿优化过程),通过以上影响因素可以确定两子区域之间的最长邻接边界,从而调整最长边界实现子区域间的栅格补偿优化。
为了证明算法的适用性,在随机生成区域中选择了3种不同类型子区域,如图11所示。通过图中3种子区域来验证算法的鲁棒性和普适性。其中,图11(a)采用常规多边形用于验证算法的普适性。图11(b)用于验证跳跃迭代法的适用性。由2.2.1可知,跳跃迭代法在应用过程中需要计算每条边的角系数,以此确定迭代间隔映射到水平方向上的取值,因此当角系数趋于无穷时,算法需要将其作为一种特殊情况进行处理。结合图5,射线法应用于凹多边形时易产生特殊的判定情况,图11(c)用于验证双特征标识改进射线法对凹多边形的适用性。
图11 3种实验多边形区域
Fig.11 Three polygonal regions used in the experiment
经算法优化后得到的数据如表4所示。表4中第4列表示扩散法得到的各子区域栅格数,第5列表示补偿法优化后,各子区域包含的栅格数。文中无人机的续航能力用栅格数表示,假设nA≤30。
表4 实验结果对比表
Table 4 The comparison of experimental results
nDnsp初步子区域栅格数最终子区域栅格数第1组188727,29,26,33,19,28,2627,27,27,27,27,27,26第2组112430,30,29,2328,28,28,28第3组156625,23,31,26,25,2626,26,26,26,26,26
从表4可以看到补偿算法的稳定性和均匀性。算法实际优化过程如图12所示。图12按4个阶段展示了区域渐进均分算法的工作过程。阶段1实现了区域栅格化(设无人机有效覆盖宽度为ls=0.5(单位长度),文中采用了文献[17]中的无人机模型,无人机飞行时仅对一侧进行有效地覆盖扫描,每次覆盖1/2栅格宽度,故栅格宽度为2ls×2ls=1×1,并且提取出(取Δk=0.01)区域的边界栅格(图12中阶段1蓝色点所示)。阶段2通过双特征标识改进射线法确定了待划分区域的栅格(图12中阶段2黄色点所示)。阶段3通过邻边扩散法确定了初步的子区域。经过邻边补偿法优化后,阶段4确定了最终分割的子区域。
图12 不同实验区域下的算法效果图
Fig.12 The effect of the algorithm under different experimental areas
注:子区域起点栅格采用随机策略选取,不同的选取方案可能会产生不同的优化结果,图12为其中一种优化结果,但不同优化结果均能保证以下特点:
1) 均匀性。经过补偿后的各子区域包含的栅格数始终可以对应等于均分数组Mav中的元素值。
2) 聚集性。各个子区域的形状呈球形聚集,保证子区域中栅格尽可能与同一子区域内的栅格邻接。
3) 连续性。各个子区域包含的每个栅格至少有一个同区域的栅格与其邻接。
图12中第2组实验和第3组实验的阶段一和阶段2证明了跳跃迭代法和双特征标识改进射线法在特殊情况下仍具有较好的适用性。通过图12中的阶段3和阶段4可以得知,不论是扩散法还是补偿法,算法计算得到的子区域保证了各子区域栅格之间的聚集性和连续性,并且经过补偿算法的优化,均分后的区域划分效果要优于初步区域划分效果。表4表明,图12中阶段3的分配方案(采用邻边扩散法初步确定的子区域栅格数)并不均匀,经过邻边补偿法的优化后,确定的各子区域包含栅格数等于均分数组Mav中各元素值。综上,经过整体算法的优化可以保证每架无人机可以在其续航能力内完成每块子区域的全覆盖侦察。
区域渐进均分算法与基于k-Means分组的Voronoi图[18-20]均分算法的对比图13。
图13 2种算法均分对比图
Fig.13 Comparison of the two algorithms on the average score
由图13可知:① 在无人机续航能力有限的条件下,当区域面积较大导致子区域划分数增大到一定程度时,基于k-Means分组的Voronoi图均分算法会在区域内部形成子区域,而在实际过程中无人机是由外向内进行区域的覆盖侦察,无人机起点应布置在区域边界上,所以Voronoi图法具有局限性;② Voronoi图法侧重的是多边形的面积均分,而区域渐进均分算法侧重的是栅格数的均分。
将Voronoi图法应用于栅格地图时,无法保证划分边界线所经过的栅格能较好的划分给各个子区域,因此算法在栅格地图的应用效果不尽人意。在同一情况下,区域渐近均分算法始终可以实现由外向内进行区域均分,因此适用于栅格地图的均分,对多无人机协同区域侦察也具有一定的参考意义。
针对单无人机无法覆盖面积较大的任意多边形侦察区域,提出了一种面向栅格地图的区域渐进均分算法。
1) 结合栅格化区域的特点,首先采用“跳跃迭代法”,以一定的迭代间隔确定区域边界栅格数组,接着通过对栅格进行凹点的筛选和栅格的双特征标识对射线法进行了改进,结合判定表确定了区域内部的栅格。
2) 建立了适用于栅格地图的区域均分算法。首先,模仿水波扩散从各子区域的起点栅格开始向内逐层扩散,并且在扩散过程中根据各子区域的规模对扩散优先级进行实时调整,避免了子区域扩散出现两极分化的情况。
3) 经过扩散法确定的各子区域栅格数可能超出无人机预分配栅格数甚至超过无人机最大可覆盖栅格数的情况,提出了邻边补偿法的优化方法,通过比较相邻子区域的特征,建立了补偿规则表,实现了子区域邻边补偿。
4) 区域渐进均分算法实验,并与应用较为广泛的Voronoi图均分算法进行了对比。当区域面积较大,导致子区域划分数较多时,相较于Voronoi图法,区域渐进均分算法始终保持从区域边界向内进行区域均分分割,对于无人机从区域外侧向内侧逐步进行全覆盖侦察具有更好的适用性和合理性。
需要注意的是,文中确定各子区域起点栅格中心时采用的是随机策略。由于侦察区域的不规则性,不同区域多边形对于不同起点栅格可能会产生不同的均分效果,但区域渐进均分算法均能保证划分的子区域满足连续性和均匀性,而聚集性会受到一定程度的影响,因此未来可考虑增强算法的稳定性和子区域聚集性。
[1]刘朝君,孙健,刘尚民.长航时无人机续航性能试飞方法优化研究[J].飞机设计,2019,39(1):5-10.LIU Chaojun,SUN Jian,LIU Shangmin.Optimization study of test flight method for long-endurance UAV endurance performance[J].Aircraft Design,2019,39(1):5-10.
[2]王湛,王江东,杨帆.关于无人机续航性能测试的研究[J].质量与认证,2019(2):73-75.WANG Zhan,WANG Jiangdong,YANG Fang.A study on the endurance performance testing of UAVs[J].Quality and Certification,2019(2):73-75.
[3]梁智韬,贾高伟,王建峰.无人机集群任务分配中群智能算法性能对比[J].火力与指挥控制,2022,47(6):28-37.LIANG Zhitao,JIA Gaowei,WANG Jianfeng.The comparison of swarm intelligence algorithm about the mission assignment of UAV swarm[J].Fire Control &Command Control,2022,47(6):28-37.
[4]袁学松.基于地理位置信息的无人机可靠路由算法[J].重庆工商大学学报(自然科学版),2021,38(1):50-56.YUAN Xuesong.Reliable routing algorithms for UAVs based on geographic location information[J].Journal of Chongqing Technology and Business University(Natural Science Edition),2021,38(1):50-56.
[5]YAN D M,WANG W P,LEVY B,et al.Efficient computation of clipped Voronoi diagram for mesh generation[J].Computer-Aided Design,2013,45(4):843-852.
[6]吴建华,张文朋,胡烈云,等.一种快速等面积分割平面简单多边形的算法[J].地理与地理信息科学,2020,36(1):1-6.WU Jianhua,ZHANG Wenpeng,HU Lieyun,et al.An algorithm for fast equal-area partitioning of planar simple polygons[J].Geography and Geographic Information Science,2020,36(1):1-6.
[7]刘苏,植江瑜.多边形分割算法在农村土地确权中的应用[J].城市勘测,2017(6):146-151.LIU Shu,ZHI Jiangyu.Application of polygon segmentation algorithm in rural land titling[J].Urban Survey,2017(6):146-151.
[8]赖靖敏.一种简单平面多边形的快速分割方法[J].北京测绘,2017(5):123-126.LAI Jingmin.A fast segmentation method for simple planar polygons[J].Beijing Surveying and Mapping,2017(5):123-126.
[9]谢小魁,王日明,冯国禄.限定面积的图斑随机分割算法设计和实现[J].测绘地理信息,2018,43(2):93-96.XIE Xiaokui,WANG Riming,FENG Guolu.Design and implementation of a random segmentation algorithm for patch with limited area[J].Mapping and Geographic Information,2018,43(2):93-96.
[10]董雪,刘润涛.基于Voronoi图的空间区域划分算法[J].哈尔滨商业大学学报(自然科学版),2011,27(6):867-869,880.DONG Xue,LIU Runtao.Spatial region segmentation algorithm based on Voronoi diagram[J].Journal of Harbin University of Commerce (Natural Science Edition),2011,27(6):867-869,880.
[11]YU X.Optimization approaches for a dubins vehicle in coverage planning problem and traveling salesman problems[D].Auburn University,Auburn,USA,2015.
[12]付婷,程雄,郝琨妍.基于邻近度的土地利用图斑综合方法[J].测绘地理息,2015,40(6):57-60.FU Ting,CHENG Xiong,HAO Kunyan.Proximity-based approach to land use patch synthesis[J].Mapping and Geographic Information,2015,40(6):57-60.
[13]GALETZKA M,GLAUNER P O.A simple and correct even-odd algorithm for the point-in-polygon problem for complex polygons[C]//Proceedings of the 12th International Joint Conference on Computer Vision,Imaging and Computer Graphics Theory and Applications (VISIGRAPP 2017).Porto,Portugal,2017:175-178.
[14]张明武,冷文韬,沈华.隐私保护的点与任意多边形位置关系判定[J].密码学报,2019,6:443-454.ZHANG Mingwu,LENG Wentao,SHEN Hua.Privacy-preserving point and arbitrary polygon position relationship determination[J].Journal of Cryptography,2019(6):443-454.
[15]HORMANN K,AGATHOS A.The point in polygon problem for arbitrary polygons[J].Computational Geometry:Theory and Applications,2001,20(3):131-144.
[16]郭飞,贾璐,吴佳静,等.基于PnPoly算法的杆塔所属风区研判及风速配置校核研究[J].宁夏电力,2021(3):35-40.GUO Fei,JIA Lu,WU Jiajing,et al.Research on wind zone research and wind speed configuration calibration based on PnPoly algorithm for wind zones belonging to towers[J].Ningxia Power,2021(3):35-40.
[17]LI Y,CHEN H,ER M J,et al.Coverage path planning for UAVs based on enhanced exact cellular decomposition method[J].Mechatronics,2011,21(5):876-885.
[18]LAZAR E A,LU J Y,RYCROFT C H.Voronoi cell analysis:The shapes of particle systems[J].American journal of physics,2022,90(3):469-480.
[19]SUHAILAN S,SAMAD S,BURHANUDDIN M A.et al.A guided ranking-based clustering using K-Means[C]//4th Mechanical Engineering Research Day (MERD).2017:255-256.
[20]YU Q,DAI B R.Accelerating k-means by grouping points automatically[C]//International Conference on Big Data Analytics and Knowledge Discovery,DAWAK 2017:199-213.